1* Culled from 970528-1.f in Burley's g77 test suite.  Copyright
2* status not clear.  Feel free to chop down if the bug is still
3* reproducible (see end of test case for how bug shows up in gdb
4* run of f771).  No particular reason it should be a noncompile
5* case, other than that I didn't want to spend time "fixing" it
6* to compile cleanly (with -O0, which works) while making sure the
7* ICE remained reproducible.  -- burley 1999-08-26
8
9* Date: Mon, 26 May 1997 13:00:19 +0200 (GMT+0200)
10* From: "D. O'Donoghue" <dod@da.saao.ac.za>
11* To: Craig Burley <burley@gnu.ai.mit.edu>
12* Cc: fortran@gnu.ai.mit.edu
13* Subject: Re: g77 problems
14
15	program dophot
16	parameter (napple = 4)
17        common /window/nwindo,ixwin(50),iywin(50),iboxwin(50),itype(50)
18        common/io/luout,ludebg
19	common/search/nstot,thresh
20	common /fitparms / acc(npmax),alim(npmax),mit,mpar,mfit1,
21     +                     mfit2,ind(npmax)
22	common /starlist/ starpar(npmax,nsmax), imtype(nsmax),
23	1shadow(npmax,nsmax),shaderr(npmax,nsmax),idstr(nsmax)
24	common /aperlist/ apple(napple ,nsmax)
25	common /parpred / ava(npmax)
26	common /unitize / ufactor
27	common /undergnd/ nfast, nslow
28	common/bzero/ scale,zero
29	common /ctimes / chiimp, apertime, filltime, addtime
30	common / drfake / needit
31	common /mfit/ psfpar(npmax),starx(nfmax),stary(nfmax),xlim,ylim
32	common /vers/ version
33 	logical needit,screen,isub,loop,comd,burn,wrtres,fixedxy
34	logical fixed,piped,debug,ex,clinfo
35	character header*5760,rhead*2880
36	character yn*1,version*40,ccd*4,infile*20
37	character*30 numf,odir,record*80
38	integer*2 instr(8)
39	character*800 line
40	external pseud0d, pseud2d, pseud4d, pseudmd, shape
41C
42C	Initialization
43	data burn,   fixedxy,fixed,  piped
44     +     /.false.,.false.,.false.,.false./
45	data needit,screen,comd,isub
46     + /.true.,.false.,.true.,.false. /
47	data acc / .01, -.03, -.03, .01, .03, .1, .03 /
48	data alim / -1.0e8, 2*-1.0e3, -1.0e8, 3*-1.0e3 /
49C
50	version = 'DoPHOT Version 1.0 LINUX May 97 '
51        debug=.false.
52        clinfo=.false.
53	line(1:800) = ' '
54	odir = ' '
55C
56C
57C	Read default tuneable parameters
58	call tuneup ( nccd, ccd, piped, debug )
59	version(33:36) = ccd(1:4)
60C
61
62        ludebg=6
63        if(piped)then
64          yn='n'
65        else
66	  write(*,'(''****************************************'')')
67	  write(*,1000) version
68	  write(*,'(''****************************************''//)')
69C
70          write(*,'(''Screen output (y/[n])? '',$)')
71	  read(*,1000) yn
72        end if
73	if(yn.eq.'y'.or.yn.eq.'Y') then
74          screen=.true.
75          luout=6
76        else
77          luout=2
78        end if
79C
80        if(piped)then
81          yn='y'
82        else
83          write(*,'(''Batch mode ([y]/n)? '',$)')
84          read(*,1000) yn
85        end if
86	if(yn.eq.'n'.or.yn.eq.'N') comd = .false.
87C
88	if(.not.comd) then
89          write(*,
90     *         '(''Do you want windowing ([y]/n)? '',$)')
91          read(*,1000)yn
92          iwindo=1
93          if(yn.eq.'n'.or.yn.eq.'N')then
94            nwindo=0
95            iwindo=0
96          end if
97C
98          write(*,
99     *       '(''Star classification info (y/[n]) ?'',$)')
100          read(*,1000)yn
101          clinfo=.false.
102          if(yn.eq.'y'.or.yn.eq.'Y')clinfo=.true.
103C
104	  write(*,
105     *        '(''Create a star-subtracted frame (y/[n])? '',$)')
106	  read(*,1000) yn
107 	  if(yn.eq.'y'.or.yn.eq.'Y') isub = .true.
108C
109	  write(*,'(''Apply after-burner (y/[n])? '',$)')
110	  read(*,1000) yn
111	  if ( yn.eq.'y'.or.yn.eq.'Y' ) burn = .true.
112	  wrtres = burn
113C
114	  write(*,'(''Read from fixed (X,Y) list (y/[n])? '',$)')
115	  read(*,1000) yn
116	  if ( yn.eq.'y'.or.yn.eq.'Y' ) then
117	    fixedxy = .true.
118	    fixed = .true.
119	    burn = .true.
120	    wrtres = .true.
121	  endif
122	endif
123        iopen=0
124C
125C       This is the start of the loop over the input files
126c
127        iframe=0
128        open(10,file='timing',status='unknown',access='append')
129
1301	ifit = 0
131	iapr = 0
132	itmn = 0
133	model = 1
134	xc = 0.0
135	yc = 0.0
136	rc = 0.0
137	ibr = 0
138	ixy = 0
139C
140        iframe=iframe+1
141        tgetpar=0.0
142        tsearch=0.0
143        tshape=0.0
144        timprove=0.0
145C
146C	Batch mode ...
147
148	if ( comd ) then
149          if(iopen.eq.0)then
150            iopen=1
151            open(11,file='dophot.bat',status='old',err=995)
152          end if
153          read(11,1000,end=999)infile
154c         now read in the parameter instructions. these are:
155c         instr(1) : if 1, specifies uncrowded field, otherwise crowded
156c         instr(2) : if 1, specifies sequential frames of same field
157c                          with a window around the stars of interest -
158c                          all other objects are ignored
159c         instr(3) : if 0, takes cmin from dophot.inp (via tuneup)
160c                    if>0, sets cmin=instr(3)
161c         instr(4) : if 0, does nothing
162c                    if 1, then opens a file called classifications
163c                    sets clinfo to .true. and writes out the star
164c                    typing info to this file
165c         instr(5) : Delete the shd.nnnnnnn file
166c         instr(6) : Delete the out.nnnnnnn file
167c         instr(7) : Delete the input frame
168c         instr(8) : Create a star-subtracted frame
169          read(11,*)instr
170          read(11,*)ifit,iapr,tmn,model,xc,yc,rc,ibr,ixy
171          nocrwd = instr(1)
172          iwindo=instr(2)
173          if(iwindo.eq.0)nwindo=0
174          itmn=tmn
175          if ( instr(3).gt.0 ) cmin=instr(3)
176          clinfo=.false.
177          if ( instr(4).gt.0 )then
178            clinfo=.true.
179            open(12,file='classifications',status='unknown')
180            ludebg=12
181          end if
182 	  if ( instr(8).ne.0 ) then
183	    isub = .true.
184	  else
185	    isub = .false.
186	  endif
187C
188	  if(ibr.ne.0) burn = .true.
189	  if(ixy.ne.0) then
190	    fixedxy = .true.
191	    fixed = .true.
192	    burn = .true.
193	    goto 20
194          endif
195          if(iwindo.eq.0)then
196            write(6,10)iframe,infile(1:15)
197   10       format('  ***** DoPHOT-ing frame ',i4,': ',a)
198            if(ludebg.eq.12)write(ludebg,11)iframe,infile(1:15)
199   11       format(////'  ',62('*')/
200     *                 '  *     DoPHOT-ing frame ',i4,': ',a,
201     *                 '                 *'/'  ',62('*'))
202          end if
203          if(iwindo.eq.1)then
204            write(6,12)iframe,infile(1:15)
205   12       format('  ***** DoPHOT-ing frame ',i4,': ',a,
206     *             '   - Windowed *****')
207            if(ludebg.eq.12)write(ludebg,13)iframe,infile(1:15)
208   13       format(////'  ',62('*')/
209     *                 '  *     DoPHOT-ing frame ',i4,': ',a,
210     *                 '   - Windowed    *'/2x,62('*'))
211          end if
212C
213C	Interactive...
214	else
215	  write(*,'(''Image name: '',$)')
216	  read(*,1000) infile
217	  if(infile(1:1).eq.' ') goto 999
2181000	  format(a)
219          write(*,'(''Crowded field mode ([y]/n) ? '',$)')
220          read(*,1000)yn
221          nocrwd=0
222          if(yn.eq.'n'.or.yn.eq.'N')nocrwd=1
223	  if(.not.fixed) then
224	    write(*,1001)
2251001        format('Sky model ([1]=Plane, 2=Power, 3=Hubble)? ',$)
226            read(*,1000)record
227            if(record.ne.' ')then
228	      read(record,*) model
229            else
230              model=1
231            end if
232	  else
233	    burn=.true.
234	    goto 20
235	  endif
236	endif
237C
238C       if windowing, open the file and read the window
239        if(iwindo.eq.1)then
240          inquire(file='windows',exist=ex)
241          if(.not.ex)go to 997
242          if(iframe.eq.1)open(9,file='windows',status='old')
243          nwindo=0
244    2     read(9,*,end=3)intype,inx,iny,inbox
245          nwindo=nwindo+1
246          if(nwindo.gt.50)then
247            print *,'too many windows - max = 50'
248            stop
249          end if
250          ixwin(nwindo)=inx
251          iywin(nwindo)=iny
252          iboxwin(nwindo)=inbox
253          itype(nwindo)=intype
254          go to 2
255
256    3     rewind 9
257          if(screen)print 4,(itype(j),ixwin(j),iywin(j),iboxwin(j),
258     *                       j=1,nwindo)
259    4     format(' Windows: Type   X    Y   Size'/
260     *           (I13,i6,i5,i5))
261        end if
262
263	t1 = cputime(0.0)
264C
265C	Read FITS frame.
266	call getfits(1,infile,header,nhead,nfast,nslow,numf,nc,line,ccd)
267C
268C	Ignore frame if not the correct chip
269	if(nc.lt.0) goto 900
270C
271C	Estimate starting PSF parameters.
272   15   call getparams(nfast,nslow,gxwid,gywid,skyval,tmin,tmax,
273     *                 iframe)
274        tgetpar = cputime(t1) + tgetpar
275        if(debug)write(ludebg,16)iframe,skyval,gxwid,gywid,tmin,tmax
276   16   format(' Getparams on frame ',i4,'  sky ',f6.1,'  gxwid ',f5.1,
277     *         '  gywid ',f5.1,'  tmin ',f5.1,'  tmax ',f5.1)
278C
279C	Initialize
280	do j=1,nsmax
281	  imtype(j) = 0
282	  do i=1,npmax
283	    shadow(i,j)=0.
284	    shaderr(i,j)=0.
285	  enddo
286	enddo
287C
288	skyguess=skyval
289	tfac = 1.0
290C	Use 4.5 X SD as fitting width
291	fitr=fitfac*(gxwid*asprat*gywid)**0.25 + 0.5
292	i=fitr
293	irect(1)=i
294	irect(2)=fitr/asprat
295C	Use 4/3 X FitFac X SD as aperture width
296	gmax = asprat*gywid
297 	if(gxwid.gt.gmax) gmax=gxwid
298	aprw = 1.33*fitfac*sqrt(gmax) + 0.5
299	i = aprw
300	arect(1) = i
301	i = aprw/asprat + 0.1
302	arect(2) = i
303C
304	if(irect(1).gt.50) irect(1)=50
305	if(irect(2).gt.50) irect(2)=50
306	if(arect(1).gt.45.) arect(1)=45.
307	if(arect(2).gt.45.) arect(2)=45.
308C
309	if (screen) call htype(line,skyval,.false.,fitr,ngr,ncon)
310C
311C       Prompt for further information
312	if ( .not.comd ) then
313          write(*,1002)
314 1002     format(/'The above are the inital parameters DoPHOT'/
315     *            'has found. You can change them now or accept'/
316     *            'the values in [ ] by pressing enter'/)
317
318          write(*,1004)tmin
319 1004     format('Enter Tmin: threshold for star detection',
320     *           ' [',f5.1,']  ',$)
321          read(*,1000)record
322          if(record.ne.' ')read(record,*)tmin
323
324          write(*,1005)cmin
325 1005     format('Enter Cmin: threshold for PSF stars',
326     *           '      [',f5.1,']  ',$)
327          read(*,1000)record
328          if(record.ne.' ')read(record,*)cmin
329
330          write(*,1006)
331 1006     format('Do you want to fix the aperture mag size ?',
332     *           ' (y/[n]) ')
333          read(*,1000)record
334          if(record.eq.'y'.or.record.eq.'Y')then
335            write(*,1007)
336 1007       format('Enter the size in pixels: ',$)
337            read(*,*)iapr
338 	    if(iapr.gt.0) then
339              arect(1)=iapr
340              i = iapr/asprat + 0.1
341              arect(2)=i
342            end if
343	  endif
344C
345	  write(*,1008)
346 1008     format('Satisfied with other input parameters ? ([y]/n)?',$)
347	  read(*,1000) yn
348          if(yn.eq.'n'.or.yn.eq.'N')then
349            yn='n'
350          else
351            yn='y'
352          end if
353	  if(.not.(yn.eq.'y'.or.yn.eq.'Y') ) call input
354	else
355	  if ( ifit.ne.0 ) then
356	    irect(1)=ifit
357	    irect(2)=(ifit/asprat + 0.1)
358	  endif
359	  if ( iapr.ne.0 ) then
360	    arect(1)=iapr
361	    i = iapr/asprat + 0.1
362	    arect(2)=i
363	  endif
364	  if ( itmn.ne.0 ) tmin = itmn
365 	  if ( .not.(xc.eq.0.0.and.yc.eq.0.0) ) then
366	    xcen = xc
367 	    ycen = yc
368          endif
369	endif
370C
371C--------------------------------
372C
373C
374	call setup ( numf,nc,screen,line,skyval,fitr,ngr,ncon,
375     +nfast, nslow )
376C
377C       if the uncrowded field option has been chosen, jump
378C       straight to the minimum threshold
379C
380        if(nocrwd.eq.1)tmax=tmin
381C
382C	Adjust tfac so that thresh ends precisely on Tmin.
383	if(tmin/tmax .gt. 0.999) then
384	  thresh = tmin
385	  tfac = 1.
386	else
387	  thresh = tmax
388	  xnum = alog10(tmax/tmin)/alog10(2.**tfac)
389	  if(xnum.gt.1.5) then
390	    xnum = float(nint(xnum))
391	  else if(xnum.ge.1) then
392	    xnum = 2.0
393	  else
394	    xnum = 1.0
395	  endif
396	  tfac = alog10(tmax/tmin)/alog10(2.)/xnum
397	endif
398C
399C------------------------------------------------------------------------
400C
401C         This is the BIG LOOP which searches the frame for stars
402C               with intensities > thresh.
403C
404C-----------------------------------------------------------------------
405C
406	loop = .true.
407	nstot = 0
408	do while ( loop )
409	  loop = thresh/tmin .ge. 1.01
410	  write(luout,1050) thresh
4111050	  format(/20('-')/'THRESHOLD: ', f10.3)
412	  if(ludebg.eq.12)write(ludebg,1050) thresh
413C
414C         Fit given model to sky values.
415C
416          call varipar(nstot, nfast, nslow )
417	  t1 = cputime(0.0)
418C
419C         Identifies potential objects in cleaned array IMG
420 	  nstar = isearch( pseud2d, nfast, nslow , clinfo)
421	  tsearch = cputime(t1) + tsearch
422C
423    	  if ( (nstar .ne. 0).or.(xnum.lt.1.5) ) then
424C
425C           Performs 7-parameter PSF fit and determines nature of object.
426	    t1 = cputime(0.0)
427	    call shape(pseud2d,pseud4d,nfast,nslow,clinfo)
428	    tshape = cputime(t1) + tshape
429C
430C           Computes average sky values etc from star list
431 	    call paravg
432  	    t1 = cputime(0.0)
433C
434C           Computes 4-parameter fits for all stellar objects using
435C           new average shape parameters.
436  	    call improve(pseud2d,nfast,nslow,clinfo)
437	    timprove = cputime(t1) + timprove
438	  end if
439C
440C         Calculate aperture photometry on last pass.
441	  if(.not.loop) call aper ( pseud2d, nstot, nfast, nslow )
442C
443 	  totaltime = (tgetpar+tsearch+tshape+timprove)
444	  write(3,1060) totaltime
445 	  write(4,1060) totaltime
446	  write(luout,1060) totaltime
4471060	  format('Total CPU time consumed:',F10.2,' seconds.')
448          write(10,1070)infile,tgetpar,tsearch,tshape,timprove,
449     *                  totaltime
4501070      format(a20,'   T(getp/f)',f5.1,'  T(search)',f5.1,
451     *               '  T(shape)',f5.1,'  T(improve)',f5.1,
452     *               '  Total',f6.1)
453	  call title (line,skyval,.false.,fitr,ngr,ncon,strint,ztot,nums)
454	  rewind(2)
455 	  rewind(3)
456  	  rewind(4)
457C
458	  call output ( line )
459C
460C         Now reduce the threshold and loop back
461C
462	  thresh = thresh/2.**tfac
463  	end do
464C
465C--------- END OF BIG LOOP ---------------------------------------
466C
467C	If after-burner required, residuals from analytic PSF are computed
468C	and stored in RES.
469C
47020	if ( burn ) then
471C
472C	If using a fixed (X,Y) coordinate list, read it.
473	 if (fixed) then
474C	 Read the image frame
475 	  call getfits(1,infile,header,nhead,nfast,nslow,numf,nc,line)
476C
477C	 Initialize arrays, open files etc.
478	  call setup ( numf,nc,screen,line,skyval,fitr,ngr,ncon,
479     +nfast, nslow )
480C
481C	 Read the XY list
482	  write(luout,'(''Reading XY list ...'')')
483	  call xylist(numf, nc, ios )
484	  if(ios.ne.0) then
485	   fixed = .false.
486	   write(luout,'(''SXY file absent or incorrect...'')')
487	   goto 15
488	  endif
489C
490	  call htype(line,skyval,.false.,fitr,ngr,ncon)
491C
492C	 Remove good stars
493	  write(luout,'(''Cleaning frame of stars: '',i8)') nstot
494	  call clean ( pseud2d, nstot, nfast, nslow, -1)
495C
496C	Calculate aperture photometry
497C	  call aper ( pseud2d, nstot, nfast, nslow )
498	 else
499	   rewind(3)
500	   rewind(4)
501	 endif
502C
503C-----------------------
504C	Flag all stars close together in groups.  Keep making the distance
505C	criterion FITR smaller until the maximum number in a group is less
506C	than NFMAX
507C
508	 fitr = amax1(arect(1),arect(2))
509	 fitr = fitr + 2.0
510	 nmax = 10000
511	 write(*,'(''Regrouping ...'')')
512C
513	 do while ( nmax.gt.nfmax )
514  	  fitr = fitr - 1.0
515	  write(luout,'(''Min distance ='',f8.1)') fitr
516	  call regroup( fitr, ngr, nmax )
517	 enddo
518C
519	 xlim = irect(1)/2
520	 ylim = irect(2)/2
521C
522C	Calculate normalized PSF residual from PSEUD2D
523	 call getres (pseud0d,pseud2d,strint,rmn,rmx,nfast,nslow,irect,
524     +arect,ztot,nums)
525	 if(nums.eq.0) then
526	  write(luout,'(''No suitable PSF stars!'')')
527	  goto 30
528	 endif
529C
530	 write(luout,'(/''AFTERBURNER tuned ON!'')')
531C
532C	Fit multiple stars in a group with enhanced PSF using box size IRECT.
533	 call mulfit( pseud2d,pseudmd,ngr,ncon,nfast,nslow,irect )
534C
535C	Re-calculate aperture photometry
536	 call aperm ( pseudmd, nstot, nfast, nslow )
537C
538	 call skyadj ( nstot )
539C
540 	 call title (line,skyval,.true.,fitr,ngr,ncon,strint,ztot,nums)
541	 call output ( line )
542	endif
543C---------------------
544C
545C-----  This section skipped if PSF residual not written out ------
546C
54730	if( isub ) then
548C
549C	Write final Cleaned array.
550 	 infile = 'x'//numf(1:nc)//'.fits'
551	 call putfits(2,infile,header,nhead,nfast,nslow)
552	 close(2)
553C
554C	If afterburner used, then residual array also written out.
555C	Find suitable scale for writing residual PSF to FITS "R" file.
556C
557 	 if ( wrtres ) then
558	  scale=20000.0/(rmx-rmn)
559	  zero=-scale*rmn
560	  do j=-nres,nres
561  	   jj=nres+j+1
562	   do i=-nres,nres
563  	    ii=nres+i+1
564	    big(ii,jj)=scale*res(i,j)+zero
565	   enddo
566	  enddo
567	  nx=2*nres+1
568C
569	  infile = 'r'//numf(1:nc)//'.fits'
570  	  zer=-zero/scale
571 	  scl=1.0/scale
572C
573C	Create a FITS header for the normalized PSF residual image
574	  call sethead(rhead,numf,nx,nx,zer,scl)
575 	  scale=1.0
576	  zero=0.0
577C	Write the normalized PSF residual image
578	  call putfits(2,infile,rhead,1,nx,nx)
579	  close(2)
580	 endif
581C
582	end if
583C
584C
585900	close(1)
586	close(3)
587 	close(4)
588	if ( .not.screen ) close(luout)
589	if(comd) then
590          if(instr(5).eq.1)call system('rm shd.'//numf(1:nc))
591          if(instr(6).eq.1)call system('rm out.'//numf(1:nc))
592          n=1
593          do while(infile(n:n).ne.' ')
594            n=n+1
595          end do
596          if(instr(7).eq.1)call system('rm '//infile(1:n-1))
597        end if
598	fixed = fixedxy
599	goto 1
600C
601995     print 996
602996     format(/'*** Fatal error ***'/
603     *          'You asked for batch processing but'/
604     *          'I cant open the "dophot.bat" file.'/
605     *          'Please make one (using batchdophot)'/
606     *          'and restart DoPHOT'/)
607        go to 999
608
609C
610997     print 998
611998     format(/'*** Fatal error ***'/
612     *          'You asked for "windowed" processing'/
613     *          'but I cant open the "windows" file.'/
614     *          'Please make one and restart DoPHOT'/)
615
616999	call exit(0)
617	end
618
619* (gdb) r
620* Starting program: /home3/craig/gnu/f77-e/gcc/f771 -quiet < ../../play/19990826-4.f -O
621* [...]
622* Breakpoint 2, fancy_abort (
623*     file=0x8285220 "../../g77-e/gcc/config/i386/i386.c", line=4399,
624*     function=0x82860df "output_fp_cc0_set") at ../../g77-e/gcc/rtl.c:1010
625* (gdb) up
626* #1  0x8222fab in output_fp_cc0_set (insn=0x8382324)
627*     at ../../g77-e/gcc/config/i386/i386.c:4399
628* (gdb) p insn
629* $1 = 0x3a
630* (gdb) up
631* #2  0x8222b81 in output_float_compare (insn=0x8382324, operands=0x82acc60)
632*     at ../../g77-e/gcc/config/i386/i386.c:4205
633* (gdb) p insn
634* $2 = 0x8382324
635* (gdb) whatis insn
636* type = rtx
637* (gdb) pr
638* (insn 2181 2180 2191 (parallel[
639*             (set (cc0)
640*                 (compare (reg:SF 8 %st(0))
641*                     (mem:SF (plus:SI (reg:SI 6 %ebp)
642*                             (const_int -9948 [0xffffd924])) 0)))
643*             (clobber (reg:HI 0 %ax))
644*         ] ) 29 {*cmpsf_cc_1} (insn_list 2173 (insn_list 2173 (nil)))
645*     (expr_list:REG_DEAD (reg:DF 8 %st(0))
646*         (expr_list:REG_UNUSED (reg:HI 0 %ax)
647*             (nil))))
648* (gdb)
649