+++ /dev/null
-* Culled from 970528-1.f in Burley's g77 test suite. Copyright
-* status not clear. Feel free to chop down if the bug is still
-* reproducible (see end of test case for how bug shows up in gdb
-* run of f771). No particular reason it should be a noncompile
-* case, other than that I didn't want to spend time "fixing" it
-* to compile cleanly (with -O0, which works) while making sure the
-* ICE remained reproducible. -- burley 1999-08-26
-
-* Date: Mon, 26 May 1997 13:00:19 +0200 (GMT+0200)
-* From: "D. O'Donoghue" <dod@da.saao.ac.za>
-* To: Craig Burley <burley@gnu.ai.mit.edu>
-* Cc: fortran@gnu.ai.mit.edu
-* Subject: Re: g77 problems
-
- program dophot
- parameter (napple = 4)
- common /window/nwindo,ixwin(50),iywin(50),iboxwin(50),itype(50)
- common/io/luout,ludebg
- common/search/nstot,thresh
- common /fitparms / acc(npmax),alim(npmax),mit,mpar,mfit1,
- + mfit2,ind(npmax)
- common /starlist/ starpar(npmax,nsmax), imtype(nsmax),
- 1shadow(npmax,nsmax),shaderr(npmax,nsmax),idstr(nsmax)
- common /aperlist/ apple(napple ,nsmax)
- common /parpred / ava(npmax)
- common /unitize / ufactor
- common /undergnd/ nfast, nslow
- common/bzero/ scale,zero
- common /ctimes / chiimp, apertime, filltime, addtime
- common / drfake / needit
- common /mfit/ psfpar(npmax),starx(nfmax),stary(nfmax),xlim,ylim
- common /vers/ version
- logical needit,screen,isub,loop,comd,burn,wrtres,fixedxy
- logical fixed,piped,debug,ex,clinfo
- character header*5760,rhead*2880
- character yn*1,version*40,ccd*4,infile*20
- character*30 numf,odir,record*80
- integer*2 instr(8)
- character*800 line
- external pseud0d, pseud2d, pseud4d, pseudmd, shape
-C
-C Initialization
- data burn, fixedxy,fixed, piped
- + /.false.,.false.,.false.,.false./
- data needit,screen,comd,isub
- + /.true.,.false.,.true.,.false. /
- data acc / .01, -.03, -.03, .01, .03, .1, .03 /
- data alim / -1.0e8, 2*-1.0e3, -1.0e8, 3*-1.0e3 /
-C
- version = 'DoPHOT Version 1.0 LINUX May 97 '
- debug=.false.
- clinfo=.false.
- line(1:800) = ' '
- odir = ' '
-C
-C
-C Read default tuneable parameters
- call tuneup ( nccd, ccd, piped, debug )
- version(33:36) = ccd(1:4)
-C
-
- ludebg=6
- if(piped)then
- yn='n'
- else
- write(*,'(''****************************************'')')
- write(*,1000) version
- write(*,'(''****************************************''//)')
-C
- write(*,'(''Screen output (y/[n])? '',$)')
- read(*,1000) yn
- end if
- if(yn.eq.'y'.or.yn.eq.'Y') then
- screen=.true.
- luout=6
- else
- luout=2
- end if
-C
- if(piped)then
- yn='y'
- else
- write(*,'(''Batch mode ([y]/n)? '',$)')
- read(*,1000) yn
- end if
- if(yn.eq.'n'.or.yn.eq.'N') comd = .false.
-C
- if(.not.comd) then
- write(*,
- * '(''Do you want windowing ([y]/n)? '',$)')
- read(*,1000)yn
- iwindo=1
- if(yn.eq.'n'.or.yn.eq.'N')then
- nwindo=0
- iwindo=0
- end if
-C
- write(*,
- * '(''Star classification info (y/[n]) ?'',$)')
- read(*,1000)yn
- clinfo=.false.
- if(yn.eq.'y'.or.yn.eq.'Y')clinfo=.true.
-C
- write(*,
- * '(''Create a star-subtracted frame (y/[n])? '',$)')
- read(*,1000) yn
- if(yn.eq.'y'.or.yn.eq.'Y') isub = .true.
-C
- write(*,'(''Apply after-burner (y/[n])? '',$)')
- read(*,1000) yn
- if ( yn.eq.'y'.or.yn.eq.'Y' ) burn = .true.
- wrtres = burn
-C
- write(*,'(''Read from fixed (X,Y) list (y/[n])? '',$)')
- read(*,1000) yn
- if ( yn.eq.'y'.or.yn.eq.'Y' ) then
- fixedxy = .true.
- fixed = .true.
- burn = .true.
- wrtres = .true.
- endif
- endif
- iopen=0
-C
-C This is the start of the loop over the input files
-c
- iframe=0
- open(10,file='timing',status='unknown',access='append')
-
-1 ifit = 0
- iapr = 0
- itmn = 0
- model = 1
- xc = 0.0
- yc = 0.0
- rc = 0.0
- ibr = 0
- ixy = 0
-C
- iframe=iframe+1
- tgetpar=0.0
- tsearch=0.0
- tshape=0.0
- timprove=0.0
-C
-C Batch mode ...
-
- if ( comd ) then
- if(iopen.eq.0)then
- iopen=1
- open(11,file='dophot.bat',status='old',err=995)
- end if
- read(11,1000,end=999)infile
-c now read in the parameter instructions. these are:
-c instr(1) : if 1, specifies uncrowded field, otherwise crowded
-c instr(2) : if 1, specifies sequential frames of same field
-c with a window around the stars of interest -
-c all other objects are ignored
-c instr(3) : if 0, takes cmin from dophot.inp (via tuneup)
-c if>0, sets cmin=instr(3)
-c instr(4) : if 0, does nothing
-c if 1, then opens a file called classifications
-c sets clinfo to .true. and writes out the star
-c typing info to this file
-c instr(5) : Delete the shd.nnnnnnn file
-c instr(6) : Delete the out.nnnnnnn file
-c instr(7) : Delete the input frame
-c instr(8) : Create a star-subtracted frame
- read(11,*)instr
- read(11,*)ifit,iapr,tmn,model,xc,yc,rc,ibr,ixy
- nocrwd = instr(1)
- iwindo=instr(2)
- if(iwindo.eq.0)nwindo=0
- itmn=tmn
- if ( instr(3).gt.0 ) cmin=instr(3)
- clinfo=.false.
- if ( instr(4).gt.0 )then
- clinfo=.true.
- open(12,file='classifications',status='unknown')
- ludebg=12
- end if
- if ( instr(8).ne.0 ) then
- isub = .true.
- else
- isub = .false.
- endif
-C
- if(ibr.ne.0) burn = .true.
- if(ixy.ne.0) then
- fixedxy = .true.
- fixed = .true.
- burn = .true.
- goto 20
- endif
- if(iwindo.eq.0)then
- write(6,10)iframe,infile(1:15)
- 10 format(' ***** DoPHOT-ing frame ',i4,': ',a)
- if(ludebg.eq.12)write(ludebg,11)iframe,infile(1:15)
- 11 format(////' ',62('*')/
- * ' * DoPHOT-ing frame ',i4,': ',a,
- * ' *'/' ',62('*'))
- end if
- if(iwindo.eq.1)then
- write(6,12)iframe,infile(1:15)
- 12 format(' ***** DoPHOT-ing frame ',i4,': ',a,
- * ' - Windowed *****')
- if(ludebg.eq.12)write(ludebg,13)iframe,infile(1:15)
- 13 format(////' ',62('*')/
- * ' * DoPHOT-ing frame ',i4,': ',a,
- * ' - Windowed *'/2x,62('*'))
- end if
-C
-C Interactive...
- else
- write(*,'(''Image name: '',$)')
- read(*,1000) infile
- if(infile(1:1).eq.' ') goto 999
-1000 format(a)
- write(*,'(''Crowded field mode ([y]/n) ? '',$)')
- read(*,1000)yn
- nocrwd=0
- if(yn.eq.'n'.or.yn.eq.'N')nocrwd=1
- if(.not.fixed) then
- write(*,1001)
-1001 format('Sky model ([1]=Plane, 2=Power, 3=Hubble)? ',$)
- read(*,1000)record
- if(record.ne.' ')then
- read(record,*) model
- else
- model=1
- end if
- else
- burn=.true.
- goto 20
- endif
- endif
-C
-C if windowing, open the file and read the window
- if(iwindo.eq.1)then
- inquire(file='windows',exist=ex)
- if(.not.ex)go to 997
- if(iframe.eq.1)open(9,file='windows',status='old')
- nwindo=0
- 2 read(9,*,end=3)intype,inx,iny,inbox
- nwindo=nwindo+1
- if(nwindo.gt.50)then
- print *,'too many windows - max = 50'
- stop
- end if
- ixwin(nwindo)=inx
- iywin(nwindo)=iny
- iboxwin(nwindo)=inbox
- itype(nwindo)=intype
- go to 2
-
- 3 rewind 9
- if(screen)print 4,(itype(j),ixwin(j),iywin(j),iboxwin(j),
- * j=1,nwindo)
- 4 format(' Windows: Type X Y Size'/
- * (I13,i6,i5,i5))
- end if
-
- t1 = cputime(0.0)
-C
-C Read FITS frame.
- call getfits(1,infile,header,nhead,nfast,nslow,numf,nc,line,ccd)
-C
-C Ignore frame if not the correct chip
- if(nc.lt.0) goto 900
-C
-C Estimate starting PSF parameters.
- 15 call getparams(nfast,nslow,gxwid,gywid,skyval,tmin,tmax,
- * iframe)
- tgetpar = cputime(t1) + tgetpar
- if(debug)write(ludebg,16)iframe,skyval,gxwid,gywid,tmin,tmax
- 16 format(' Getparams on frame ',i4,' sky ',f6.1,' gxwid ',f5.1,
- * ' gywid ',f5.1,' tmin ',f5.1,' tmax ',f5.1)
-C
-C Initialize
- do j=1,nsmax
- imtype(j) = 0
- do i=1,npmax
- shadow(i,j)=0.
- shaderr(i,j)=0.
- enddo
- enddo
-C
- skyguess=skyval
- tfac = 1.0
-C Use 4.5 X SD as fitting width
- fitr=fitfac*(gxwid*asprat*gywid)**0.25 + 0.5
- i=fitr
- irect(1)=i
- irect(2)=fitr/asprat
-C Use 4/3 X FitFac X SD as aperture width
- gmax = asprat*gywid
- if(gxwid.gt.gmax) gmax=gxwid
- aprw = 1.33*fitfac*sqrt(gmax) + 0.5
- i = aprw
- arect(1) = i
- i = aprw/asprat + 0.1
- arect(2) = i
-C
- if(irect(1).gt.50) irect(1)=50
- if(irect(2).gt.50) irect(2)=50
- if(arect(1).gt.45.) arect(1)=45.
- if(arect(2).gt.45.) arect(2)=45.
-C
- if (screen) call htype(line,skyval,.false.,fitr,ngr,ncon)
-C
-C Prompt for further information
- if ( .not.comd ) then
- write(*,1002)
- 1002 format(/'The above are the inital parameters DoPHOT'/
- * 'has found. You can change them now or accept'/
- * 'the values in [ ] by pressing enter'/)
-
- write(*,1004)tmin
- 1004 format('Enter Tmin: threshold for star detection',
- * ' [',f5.1,'] ',$)
- read(*,1000)record
- if(record.ne.' ')read(record,*)tmin
-
- write(*,1005)cmin
- 1005 format('Enter Cmin: threshold for PSF stars',
- * ' [',f5.1,'] ',$)
- read(*,1000)record
- if(record.ne.' ')read(record,*)cmin
-
- write(*,1006)
- 1006 format('Do you want to fix the aperture mag size ?',
- * ' (y/[n]) ')
- read(*,1000)record
- if(record.eq.'y'.or.record.eq.'Y')then
- write(*,1007)
- 1007 format('Enter the size in pixels: ',$)
- read(*,*)iapr
- if(iapr.gt.0) then
- arect(1)=iapr
- i = iapr/asprat + 0.1
- arect(2)=i
- end if
- endif
-C
- write(*,1008)
- 1008 format('Satisfied with other input parameters ? ([y]/n)?',$)
- read(*,1000) yn
- if(yn.eq.'n'.or.yn.eq.'N')then
- yn='n'
- else
- yn='y'
- end if
- if(.not.(yn.eq.'y'.or.yn.eq.'Y') ) call input
- else
- if ( ifit.ne.0 ) then
- irect(1)=ifit
- irect(2)=(ifit/asprat + 0.1)
- endif
- if ( iapr.ne.0 ) then
- arect(1)=iapr
- i = iapr/asprat + 0.1
- arect(2)=i
- endif
- if ( itmn.ne.0 ) tmin = itmn
- if ( .not.(xc.eq.0.0.and.yc.eq.0.0) ) then
- xcen = xc
- ycen = yc
- endif
- endif
-C
-C--------------------------------
-C
-C
- call setup ( numf,nc,screen,line,skyval,fitr,ngr,ncon,
- +nfast, nslow )
-C
-C if the uncrowded field option has been chosen, jump
-C straight to the minimum threshold
-C
- if(nocrwd.eq.1)tmax=tmin
-C
-C Adjust tfac so that thresh ends precisely on Tmin.
- if(tmin/tmax .gt. 0.999) then
- thresh = tmin
- tfac = 1.
- else
- thresh = tmax
- xnum = alog10(tmax/tmin)/alog10(2.**tfac)
- if(xnum.gt.1.5) then
- xnum = float(nint(xnum))
- else if(xnum.ge.1) then
- xnum = 2.0
- else
- xnum = 1.0
- endif
- tfac = alog10(tmax/tmin)/alog10(2.)/xnum
- endif
-C
-C------------------------------------------------------------------------
-C
-C This is the BIG LOOP which searches the frame for stars
-C with intensities > thresh.
-C
-C-----------------------------------------------------------------------
-C
- loop = .true.
- nstot = 0
- do while ( loop )
- loop = thresh/tmin .ge. 1.01
- write(luout,1050) thresh
-1050 format(/20('-')/'THRESHOLD: ', f10.3)
- if(ludebg.eq.12)write(ludebg,1050) thresh
-C
-C Fit given model to sky values.
-C
- call varipar(nstot, nfast, nslow )
- t1 = cputime(0.0)
-C
-C Identifies potential objects in cleaned array IMG
- nstar = isearch( pseud2d, nfast, nslow , clinfo)
- tsearch = cputime(t1) + tsearch
-C
- if ( (nstar .ne. 0).or.(xnum.lt.1.5) ) then
-C
-C Performs 7-parameter PSF fit and determines nature of object.
- t1 = cputime(0.0)
- call shape(pseud2d,pseud4d,nfast,nslow,clinfo)
- tshape = cputime(t1) + tshape
-C
-C Computes average sky values etc from star list
- call paravg
- t1 = cputime(0.0)
-C
-C Computes 4-parameter fits for all stellar objects using
-C new average shape parameters.
- call improve(pseud2d,nfast,nslow,clinfo)
- timprove = cputime(t1) + timprove
- end if
-C
-C Calculate aperture photometry on last pass.
- if(.not.loop) call aper ( pseud2d, nstot, nfast, nslow )
-C
- totaltime = (tgetpar+tsearch+tshape+timprove)
- write(3,1060) totaltime
- write(4,1060) totaltime
- write(luout,1060) totaltime
-1060 format('Total CPU time consumed:',F10.2,' seconds.')
- write(10,1070)infile,tgetpar,tsearch,tshape,timprove,
- * totaltime
-1070 format(a20,' T(getp/f)',f5.1,' T(search)',f5.1,
- * ' T(shape)',f5.1,' T(improve)',f5.1,
- * ' Total',f6.1)
- call title (line,skyval,.false.,fitr,ngr,ncon,strint,ztot,nums)
- rewind(2)
- rewind(3)
- rewind(4)
-C
- call output ( line )
-C
-C Now reduce the threshold and loop back
-C
- thresh = thresh/2.**tfac
- end do
-C
-C--------- END OF BIG LOOP ---------------------------------------
-C
-C If after-burner required, residuals from analytic PSF are computed
-C and stored in RES.
-C
-20 if ( burn ) then
-C
-C If using a fixed (X,Y) coordinate list, read it.
- if (fixed) then
-C Read the image frame
- call getfits(1,infile,header,nhead,nfast,nslow,numf,nc,line)
-C
-C Initialize arrays, open files etc.
- call setup ( numf,nc,screen,line,skyval,fitr,ngr,ncon,
- +nfast, nslow )
-C
-C Read the XY list
- write(luout,'(''Reading XY list ...'')')
- call xylist(numf, nc, ios )
- if(ios.ne.0) then
- fixed = .false.
- write(luout,'(''SXY file absent or incorrect...'')')
- goto 15
- endif
-C
- call htype(line,skyval,.false.,fitr,ngr,ncon)
-C
-C Remove good stars
- write(luout,'(''Cleaning frame of stars: '',i8)') nstot
- call clean ( pseud2d, nstot, nfast, nslow, -1)
-C
-C Calculate aperture photometry
-C call aper ( pseud2d, nstot, nfast, nslow )
- else
- rewind(3)
- rewind(4)
- endif
-C
-C-----------------------
-C Flag all stars close together in groups. Keep making the distance
-C criterion FITR smaller until the maximum number in a group is less
-C than NFMAX
-C
- fitr = amax1(arect(1),arect(2))
- fitr = fitr + 2.0
- nmax = 10000
- write(*,'(''Regrouping ...'')')
-C
- do while ( nmax.gt.nfmax )
- fitr = fitr - 1.0
- write(luout,'(''Min distance ='',f8.1)') fitr
- call regroup( fitr, ngr, nmax )
- enddo
-C
- xlim = irect(1)/2
- ylim = irect(2)/2
-C
-C Calculate normalized PSF residual from PSEUD2D
- call getres (pseud0d,pseud2d,strint,rmn,rmx,nfast,nslow,irect,
- +arect,ztot,nums)
- if(nums.eq.0) then
- write(luout,'(''No suitable PSF stars!'')')
- goto 30
- endif
-C
- write(luout,'(/''AFTERBURNER tuned ON!'')')
-C
-C Fit multiple stars in a group with enhanced PSF using box size IRECT.
- call mulfit( pseud2d,pseudmd,ngr,ncon,nfast,nslow,irect )
-C
-C Re-calculate aperture photometry
- call aperm ( pseudmd, nstot, nfast, nslow )
-C
- call skyadj ( nstot )
-C
- call title (line,skyval,.true.,fitr,ngr,ncon,strint,ztot,nums)
- call output ( line )
- endif
-C---------------------
-C
-C----- This section skipped if PSF residual not written out ------
-C
-30 if( isub ) then
-C
-C Write final Cleaned array.
- infile = 'x'//numf(1:nc)//'.fits'
- call putfits(2,infile,header,nhead,nfast,nslow)
- close(2)
-C
-C If afterburner used, then residual array also written out.
-C Find suitable scale for writing residual PSF to FITS "R" file.
-C
- if ( wrtres ) then
- scale=20000.0/(rmx-rmn)
- zero=-scale*rmn
- do j=-nres,nres
- jj=nres+j+1
- do i=-nres,nres
- ii=nres+i+1
- big(ii,jj)=scale*res(i,j)+zero
- enddo
- enddo
- nx=2*nres+1
-C
- infile = 'r'//numf(1:nc)//'.fits'
- zer=-zero/scale
- scl=1.0/scale
-C
-C Create a FITS header for the normalized PSF residual image
- call sethead(rhead,numf,nx,nx,zer,scl)
- scale=1.0
- zero=0.0
-C Write the normalized PSF residual image
- call putfits(2,infile,rhead,1,nx,nx)
- close(2)
- endif
-C
- end if
-C
-C
-900 close(1)
- close(3)
- close(4)
- if ( .not.screen ) close(luout)
- if(comd) then
- if(instr(5).eq.1)call system('rm shd.'//numf(1:nc))
- if(instr(6).eq.1)call system('rm out.'//numf(1:nc))
- n=1
- do while(infile(n:n).ne.' ')
- n=n+1
- end do
- if(instr(7).eq.1)call system('rm '//infile(1:n-1))
- end if
- fixed = fixedxy
- goto 1
-C
-995 print 996
-996 format(/'*** Fatal error ***'/
- * 'You asked for batch processing but'/
- * 'I cant open the "dophot.bat" file.'/
- * 'Please make one (using batchdophot)'/
- * 'and restart DoPHOT'/)
- go to 999
-
-C
-997 print 998
-998 format(/'*** Fatal error ***'/
- * 'You asked for "windowed" processing'/
- * 'but I cant open the "windows" file.'/
- * 'Please make one and restart DoPHOT'/)
-
-999 call exit(0)
- end
-
-* (gdb) r
-* Starting program: /home3/craig/gnu/f77-e/gcc/f771 -quiet < ../../play/19990826-4.f -O
-* [...]
-* Breakpoint 2, fancy_abort (
-* file=0x8285220 "../../g77-e/gcc/config/i386/i386.c", line=4399,
-* function=0x82860df "output_fp_cc0_set") at ../../g77-e/gcc/rtl.c:1010
-* (gdb) up
-* #1 0x8222fab in output_fp_cc0_set (insn=0x8382324)
-* at ../../g77-e/gcc/config/i386/i386.c:4399
-* (gdb) p insn
-* $1 = 0x3a
-* (gdb) up
-* #2 0x8222b81 in output_float_compare (insn=0x8382324, operands=0x82acc60)
-* at ../../g77-e/gcc/config/i386/i386.c:4205
-* (gdb) p insn
-* $2 = 0x8382324
-* (gdb) whatis insn
-* type = rtx
-* (gdb) pr
-* (insn 2181 2180 2191 (parallel[
-* (set (cc0)
-* (compare (reg:SF 8 %st(0))
-* (mem:SF (plus:SI (reg:SI 6 %ebp)
-* (const_int -9948 [0xffffd924])) 0)))
-* (clobber (reg:HI 0 %ax))
-* ] ) 29 {*cmpsf_cc_1} (insn_list 2173 (insn_list 2173 (nil)))
-* (expr_list:REG_DEAD (reg:DF 8 %st(0))
-* (expr_list:REG_UNUSED (reg:HI 0 %ax)
-* (nil))))
-* (gdb)