C @(#)peak.for 17.1.1.1 (ES0-DMD) 01/25/02 17:14:09 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C SUBROUTINE DAOPK (PAR, MAXPAR, PSF, MAXPSF, MAXEXP, F, . WATCH, FITRAD, PERERR, PROERR) C C======================================================================= C C Single-star profile-fitting routine. C C OFFICIAL DAO VERSION: 1991 April 18 C C This subroutine reads in an array around each star in an input C photometry file and fits the point-spread function to the observed C stellar profile by least-squares. The three parameters solved for C are the brightness ratio between the PSF and the program star, and C the coordinates of the centroid of the program star. C C Arguments C C WATCH (INPUT) governs whether information relating to the progress C of the reductions is to be typed on the terminal screen C during execution. C C FITRAD (INPUT) is the fitting radius. Only pixels within FITRAD of C the current estimate of a star's centroid will be included in C the least-squares determination of its centroid and magnitude. C C All are user-definable optional parameters. C C======================================================================= C IMPLICIT NONE INTEGER MAXPSF, MAXEXP, MAXBOX, MAXPAR PARAMETER (MAXBOX=69) C C Parameters C C MAXPSF is the length of the side of the largest PSF look-up table C allowed. (Note: half-pixel grid size) C C MAXBOX is the length of the side of the largest box within which the C PSF can be evaluated. (MAXBOX = (MAXPSF-7)/2 ) C CHARACTER*30 COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL, SWITCH CHARACTER LINE*80, EXTEND*30, CASE*3 REAL F(MAXBOX,MAXBOX) REAL PSF(MAXPSF,MAXPSF,MAXEXP), PAR(MAXPAR) INTEGER NINT, MIN0 C REAL LOBAD, AMAG, SHARP, CHI, ERRMAG, DELTAX, DELTAY REAL APMAG, RADIUS, FRAD, READNS, AP1, THRESH, HIBAD REAL PROERR, PERERR, FITRAD, WATCH, PHPADU, RONOIS REAL PERR, PKERR REAL X, Y, PSFMAG, BRIGHT, XPSF, YPSF, SCALE, SKY INTEGER RDPSF, NITER, IST, LX, LY, NX, NY INTEGER NFRAC, IBEG, ISTAR, NEXP, NPAR, NPSF, IPSTYP INTEGER IDUM, NL, ISTAT, NCOL, NROW, NTERM C COMMON /SIZE/ NCOL, NROW COMMON /FILNAM/ COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL COMMON /ERROR/ PHPADU, RONOIS, PERR, PKERR C C----------------------------------------------------------------------- C C SECTION 1 C C Ascertain the name of the file containing coordinates and magnitude C estimates for the program stars. Open it and read the header. C CALL TBLANK 950 CALL GETNAM ('File with aperture results:', MAGFIL) IF ((MAGFIL .EQ. 'END OF FILE') .OR. . (MAGFIL .EQ. 'GIVE UP')) THEN MAGFIL = ' ' RETURN END IF CALL INFILE (2, MAGFIL, ISTAT) IF (ISTAT .NE. 0) THEN MAGFIL = 'GIVE UP' GO TO 950 END IF C CALL RDHEAD (2, NL, IDUM, IDUM, LOBAD, HIBAD, THRESH, AP1, . PHPADU, READNS, FRAD) RONOIS=READNS**2 C C Ascertain the name of the PSF file, open it, and read the PSF. C 900 CALL GETNAM ('File with the PSF:', PSFFIL) IF ((PSFFIL .EQ. 'END OF FILE') .OR. . (PSFFIL .EQ. 'GIVE UP')) THEN PSFFIL = ' ' CALL CLFILE (2) RETURN END IF C ISTAT = RDPSF (PSFFIL, IPSTYP, PAR, MAXPAR, NPAR, . PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, . PSFMAG, BRIGHT, XPSF, YPSF) IF (ISTAT .NE. 0) THEN PSFFIL = 'GIVE UP' GO TO 900 END IF C NTERM = NEXP+NFRAC PERR = 0.01*PERERR PKERR = 0.01*PROERR/(PAR(1)*PAR(2)) ! **2 C C Get ready to do the PEAK fitting: get the name of the output file. C Open it and write the header. C PROFIL=SWITCH(MAGFIL, CASE('.pk')) 960 CALL GETNAM ('File for PEAK results:', PROFIL) IF ((PROFIL .EQ. 'END OF FILE') .OR. . (PROFIL .EQ. 'GIVE UP')) THEN CALL CLFILE (2) PROFIL = ' ' RETURN END IF PROFIL = EXTEND(PROFIL, CASE('pk')) CALL OUTFIL (3, PROFIL, ISTAT) IF (ISTAT .NE. 0) THEN CALL STUPID ('Error opening output file '//PROFIL) PROFIL = 'GIVE UP' GO TO 960 END IF C RADIUS = AMIN1( FITRAD, (REAL(NPSF-1)/2. - 1.)/2. ) CALL WRHEAD (3, 1, NCOL, NROW, 7, LOBAD, HIBAD, THRESH, AP1, . PHPADU, READNS, RADIUS) IF (WATCH .GT. 0.5) CALL TBLANK ! Type a blank line IBEG=1 C C----------------------------------------------------------------------- C C SECTION 2 C C Reduce the stars, one by one. C C Read data for the next star and initialize things. C 2000 CALL RDSTAR (2, NL, ISTAR, X, Y, APMAG, SKY) IF (ISTAR .LT. 0) GO TO 9000 ! End-of-file encountered IF (ISTAR .EQ. 0) GO TO 2000 ! Blank line encountered IF (APMAG .GE. 99.) APMAG=PSFMAG+5. ! The ol' college try DELTAX=(X-1.)/XPSF-1. DELTAY=(Y-1.)/YPSF-1. LX=MAX0(1, INT(X-RADIUS)+1) LY=MAX0(1, INT(Y-RADIUS)+1) NX=MIN0(NCOL, INT(X+RADIUS)) - LX + 1 NY=MIN0(NROW, INT(Y+RADIUS)) - LY + 1 CALL RDARAY ('DATA', LX, LY, NX, NY, MAXBOX, F, IST) C C At this point LX, LY are the coordinates in the big picture of the C first pixel in the subarray. The dimensions of the subarray are C NX, NY. C X=X-LX+1 Y=Y-LY+1 C C X,Y are now the coordinates of the star's centroid in the subframe C (where (x,y)=(1.0,1.0) are the coordinates of the center of the first C pixel and (x,y)=(NX,NY) are the coordinates of the center of the last C pixel in the subarray). C C Display the subarray on the terminal, if desired, and type out headers C if appropriate. C IF (WATCH .GT. 1.5) CALL SHOW (F, F(NINT(X),NINT(Y)), 0.9*SKY, . NX, NY, MAXBOX) IF ((WATCH .GT. 1.5) .OR. ((WATCH .GT. 0.5) .AND. (IBEG .EQ. 1))) . WRITE (6,620) 620 FORMAT (9X, 'STAR', 6X, 'X', 8X, 'Y', 8X, 'MAGNITUDE', 7X, . 'CHI SHARP IT') IBEG=2 SCALE=10.**(-0.4*(APMAG-PSFMAG)) ! Starting value CALL PKFIT (F, NX, NY, MAXBOX, X, Y, SCALE, SKY, RADIUS, LOBAD, . HIBAD, BRIGHT, IPSTYP, PAR, MAXPAR, NPAR, . PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, . DELTAX, DELTAY, ERRMAG, CHI, SHARP, NITER) IF (NITER .GT. 0) GO TO 2010 C C A singular matrix occurred during the least-squares solution. C WRITE (LINE,621) ISTAR 621 FORMAT (8X, I5, ' had a singular matrix.') CALL STUPID (LINE) GO TO 2000 C C Everything went fine. C 2010 X=X+LX-1. Y=Y+LY-1. AMAG=PSFMAG-2.5*ALOG10(SCALE) ERRMAG=AMIN1(2.0, 1.086*ERRMAG/SCALE) IF (WATCH .GT. 0.5) WRITE (6,622) ISTAR, X, Y, AMAG, ERRMAG, . CHI, SHARP, NITER 622 FORMAT (8X, I5, 2F9.2, F9.3, ' +- ', F5.3, 2F7.2, I5) WRITE (3,320) ISTAR, X, Y, AMAG, ERRMAG, SKY, FLOAT(NITER), CHI, . SHARP 320 FORMAT (I6, 3F9.3, F9.4, F9.3, F9.0, F9.2, F9.3) GO TO 2000 C C----------------------------------------------------------------------- C C Normal return. C 9000 CALL CLFILE (3) CALL CLFILE (2) CALL STUPID (' Done.') RETURN C END! C C####################################################################### C SUBROUTINE PKFIT (F, NX, NY, MAXBOX, X, Y, SCALE, SKY, RADIUS, . LOBAD, HIBAD, BRIGHT, IPSTYP, PAR, MAXPAR, NPAR, . PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, . DELTAX, DELTAY, ERRMAG, CHI, SHARP, NITER) C C======================================================================= C C This is the subroutine which does the actual one-star least-squares C profile fit for PEAK. C C OFFICIAL DAO VERSION: 1991 April 6 C C Arguments C F (INPUT) is an NX by NY array containing actual picture data. C C MAXBOX (INPUT) is the maximum value allowable for either NX or NY, C needed for the dimension statements below. PEAK and PSF will C provide different values of MAXBOX. C C SCALE (INPUT/OUTPUT) is the initial estimate of the brightness of C the star, expressed as a fraction of the brightness of the C PSF. Upon return, the final computed value of SCALE will C be passed back to the calling routine. C C X, Y (INPUT/OUTPUT) are the initial estimates of the centroid of C the star relative to the corner (1,1) of the subarray. Upon C return, the final computed values of X and Y will be passed C back to the calling routine. C C SKY (INPUT) is the local sky brightness value, carried on from C PHOTOMETRY via the data files. C C RADIUS (INPUT) is the fitting radius-- only pixels within RADIUS of C the instantaneous estimate of the star's centroid will be C included in the fit. C C LOBAD and HIBAD (INPUT) are bad pixel limits-- any pixel whose C brightness value falls outside this range will be presumed to C be bad, and will be ignored. C C BRIGHT (INPUT) contains the brightness normalization of the model C analytic point spread function C C PAR (INPUT) contains the values of the remaining NPARAM parameters C defining the analytic function which approximates the core of C the PSF. C C PSF (INPUT) is an NPSF by NPSF by NEXP+NFRAC look-up table C containing corrections from the analytic approximation of the C PSF to the true PSF. C C ERRMAG (OUTPUT) is the estimated standard error of the value of SCALE C returned by this routine. C C CHI (OUTPUT) is the estimated goodness-of-fit statistic: the ratio C of the observed pixel-to-pixel mean absolute deviation from C the profile fit, to the value expected on the basis of the C read-out noise and the photons/ADU (which are brought in C through COMMON block /ERROR/). C C SHARP (OUTPUT) is a goodness-of-fit statistic describing how much C broader the actual profile of the object appears than the C profile of the PSF. C C NITER (OUTPUT) is the number of iterations the solution required to C achieve convergence. If NITER = 50, the solution did not C converge. If for some reason a singular matrix occurs during C the least-squares solution, this will be flagged by setting C NITER = -1. C C======================================================================= C IMPLICIT NONE INTEGER MAXPSF, MAXEXP, MAXBOX, MAXPAR C C Parameter C C MAXPSF is the length of the side of the largest PSF look-up table C allowed. (Note: half-pixel grid spacing) C REAL C(3,3), V(3), CLAMP(3), DTOLD(3) REAL F(MAXBOX,MAXBOX), T(3), DT(3), NUMER REAL PSF(MAXPSF,MAXPSF,MAXEXP), PAR(MAXPAR) REAL USEPSF, AMAX1, AMIN1, ABS INTEGER MIN0, MAX0 C REAL LOBAD, DF, WT, DFDSIG, RHOSQ, RELERR, SIG, SIGSQ REAL FPOS, DVDXC, DVDYC, RSQ, DX, DXSQ, DY, DYSQ, DATUM REAL DENOM, SUMWT, CHIOLD, RADSQ, SHARP, CHI, ERRMAG REAL DELTAX, DELTAY, X, Y, PHPADU, RONOIS, PERR, PKERR REAL SCALE, SKY, RADIUS, HIBAD, BRIGHT INTEGER I, J, IX, IY, IXLO, IXHI, IYLO, IYHI, ISTAT INTEGER NPIX, NFRAC, NEXP, NITER, NPAR, NPSF, IPSTYP INTEGER NX, NY, NCOL, NROW LOGICAL CLIP, REDO C COMMON /SIZE/ NCOL, NROW COMMON /ERROR/ PHPADU, RONOIS, PERR, PKERR C C----------------------------------------------------------------------- C C Initialize a few things for the solution. C RADSQ=RADIUS**2 DO 1010 I=1,3 CLAMP(I)=1. 1010 DTOLD(I)=0.0 CHIOLD=1. NITER=0 SHARP=0. CLIP = .FALSE. C C----------------------------------------------------------------------- C C Here begins the big least-squares loop. C 2000 NITER=NITER+1 C C Initialize things for this iteration. CHI and SHARP will be C goodness-of-fit indices. CHI will also be used in determining the C weights of the individual pixels. As the solution iterates, the new C value of CHI will be built up in the variable CHI, while a smoothed C value of CHI computed from the previous iteration will be carried C along in CHIOLD. C CHI=0.0 SUMWT=0.0 NUMER=0.0 DENOM=0.0 DO 2010 I=1,3 V(I)=0.0 ! Zero the vector of residuals DO 2010 J=1,3 2010 C(I,J)=0.0 ! Zero the normal matrix C C Choose the little box containing points inside the fitting radius. C IXLO=MAX0(1, INT(X-RADIUS)) IYLO=MAX0(1, INT(Y-RADIUS)) IXHI=MIN0(NX, INT(X+RADIUS)+1) IYHI=MIN0(NY, INT(Y+RADIUS)+1) C C Now build up the normal matrix and vector of residuals. C NPIX=0 DO 2090 IY=IYLO,IYHI DY=FLOAT(IY)-Y DYSQ=DY**2 DO 2090 IX=IXLO,IXHI DATUM=F(IX,IY) IF ((DATUM .LT. LOBAD) .OR. (DATUM .GT. HIBAD)) GO TO 2090 DX=FLOAT(IX)-X DXSQ=DX**2 C C DX and DY are the distance of this pixel from the centroid of the C star. Is this pixel truly inside the fitting radius? C RSQ=(DXSQ+DYSQ)/RADSQ IF (1.-RSQ .LE. 2.E-6) GO TO 2090 ! Prevents floating underflows C C The fitting equation is of the form C C Observed brightness = C [SCALE + delta(SCALE)] * [PSF + delta(Xcen)*d(PSF)/d(Xcen) + C delta(Ycen)*d(PSF)/d(Ycen) ] C C and is solved for the unknowns delta(SCALE) ( = the correction to C the brightness ratio between the program star and the PSF) and C delta(Xcen) and delta(Ycen) ( = corrections to the program star's C centroid). C C The point-spread function is equal to the sum of the integral under C a two-dimensional analytic profile plus values interpolated from C a look-up table. C T(1)=USEPSF(IPSTYP, DX, DY, BRIGHT, PAR, PSF, NPSF, NPAR, NEXP, . NFRAC, DELTAX, DELTAY, DVDXC, DVDYC) IF ((SCALE*T(1)+SKY .GT. HIBAD) .AND. (NITER.GE.3)) GO TO 2090 T(2)=SCALE*DVDXC T(3)=SCALE*DVDYC DF=F(IX,IY)-SKY-SCALE*T(1) C C DF is the residual of the brightness in this pixel from the PSF fit. C C The expected random error in the pixel is the quadratic sum of C the Poisson statistics, plus the readout noise, plus an estimated C error of PERERR% of the total brightness for the difficulty of flat- C fielding and bias-correcting the chip, plus an estimated error of C some fraction of the fourth derivative at the peak of the profile, C to account for the difficulty of accurately interpolating within the C point-spread function. The fourth derivative of the PSF is roughly C proportional to H/sigma**4 (sigma is the width parameter for C the stellar core); using the geometric mean of sigma(x) and sigma(y), C this becomes H/[sigma(x)*sigma(y)]**2. The ratio of the fitting C error to this quantity is estimated from a good-seeing CTIO frame to C be approximately 0.027. NOWADAYS, PAR(1) and PAR(2) are the HWHM, C not the sigma, so the coefficient is now of order 0.05 = 5.0%. C FPOS=AMAX1(0., F(IX,IY)-DF) C C FPOS = raw data minus residual = model-predicted value of the C intensity at this point (which presumably is non-negative). C SIGSQ=FPOS/PHPADU+RONOIS+(PERR*FPOS)**2+(PKERR*(FPOS-SKY))**2 SIG=SQRT(SIGSQ) RELERR=ABS(DF/SIG) C C SIG is the anticipated standard error of the intensity in this pixel, C including readout noise, Poisson photon statistics, and an estimate C of the standard error of interpolating within the PSF. C WT=5./(5.+RSQ/(1.-RSQ)) ! Weight as function of radius C C Now add this pixel into the quantities which go to make up the SHARP C index. C RHOSQ=DXSQ/PAR(1)**2+DYSQ/PAR(2)**2 C C Include in the sharpness index only those pixels within six C HWHMs of the centroid of the star. (This saves time and C floating underflows by excluding pixels which contribute less than C about one part in a million to the index.) C IF (RHOSQ .LE. 36.) THEN RHOSQ=0.6931472*RHOSQ DFDSIG=EXP(-RHOSQ)*(RHOSQ-1.) FPOS=AMAX1(0., F(IX,IY)-SKY)+SKY NUMER=NUMER+DFDSIG*DF/SIGSQ DENOM=DENOM+DFDSIG**2/SIGSQ END IF C C Derive the weight of this pixel. First of all, the weight depends C upon the distance of the pixel from the centroid of the star-- it C is determined from a function which is very nearly unity for radii C much smaller than the fitting radius, and which goes to zero for C radii very near the fitting radius. Then reject any pixels with C 10-sigma residuals (after the first iteration). C CHI=CHI+WT*RELERR SUMWT=SUMWT+WT C C Now the weight is scaled to the inverse square of the expected mean C error. C WT=WT/SIGSQ C C Reduce the weight of a bad pixel. A pixel having a residual of 2.5 C sigma gets reduced to half weight; a pixel having a residual of 5. C sigma gets weight 1/257. C IF (CLIP) WT=WT/(1.+(0.4*RELERR/CHIOLD)**8) C C Now add the pixel into the vector of residuals and the normal matrix. C DO 2030 I=1,3 V(I)=V(I)+DF*T(I)*WT DO 2030 J=1,3 2030 C(I,J)=C(I,J)+T(I)*T(J)*WT C NPIX=NPIX+1 2090 CONTINUE ! End of loop over pixels C C Compute the (robust) goodness-of-fit index CHI. C IF (SUMWT .GT. 3.) THEN CHI=1.2533141*CHI/SQRT(SUMWT*(SUMWT-3.)) C C CHI is pulled toward its expected value of unity before being stored C in CHIOLD to keep the statistics of a small number of pixels from C completely dominating the error analysis. C CHIOLD=((SUMWT-3.)*CHI+3.)/SUMWT ELSE CHI = 1. CHIOLD = 1. END IF C C Compute the parameter corrections and check for convergence. C IF (NPIX .LT. 3) THEN NITER=-1 RETURN END IF CALL INVERS (C, 3, 3, ISTAT) IF (ISTAT .EQ. 0) GO TO 2100 NITER=-1 RETURN C C Everything OK so far. C 2100 CALL VMUL (C, 3, 3, V, DT) C C In the beginning, the brightness of the star will not be permitted C to change by more than two magnitudes per iteration (that is to say, C if the estimate is getting brighter, it may not get brighter by C more than 525% per iteration, and if it is getting fainter, it may C not get fainter by more than 84% per iteration). The x and y C coordinates of the centroid will be allowed to change by no more C than one-half pixel per iteration. Any time that a parameter C correction changes sign, the maximum permissible change in that C parameter will be reduced by a factor of 2. C DO 2110 I=1,3 IF (DTOLD(I)*DT(I) .LT. 0.) THEN CLAMP(I)=0.5*CLAMP(I) ELSE CLAMP(I)=MIN(1., 1.1*CLAMP(I)) END IF 2110 DTOLD(I)=DT(I) C SCALE=SCALE+DT(1)/ . (1.+AMAX1(DT(1)/(5.25*SCALE),-DT(1)/(0.84*SCALE))/CLAMP(1)) X=X+DT(2)/(1.+ABS(DT(2))/(0.4*CLAMP(2))) Y=Y+DT(3)/(1.+ABS(DT(3))/(0.4*CLAMP(3))) IF (NITER .LE. 1) GO TO 2000 REDO=.FALSE. C C Convergence criteria: if the most recent computed correction to the C brightness is larger than 0.01% or than 0.05 * sigma(brightness), C whichever is larger, OR if the absolute change in X or Y is C greater than 0.001 pixels, convergence has not been achieved. C ERRMAG=CHIOLD*SQRT(C(1,1)) IF (CLIP) THEN IF (ABS(DT(1)) .GT. . AMAX1( 0.05*ERRMAG, 0.0001*SCALE )) THEN REDO=.TRUE. ELSE IF (AMAX1(ABS(DT(2)),ABS(DT(3))) .GT. 0.001) THEN REDO=.TRUE. END IF ELSE IF (ABS(DT(1)) .GT. . AMAX1( ERRMAG, 0.002*SCALE )) THEN REDO = .TRUE. ELSE IF (AMAX1(ABS(DT(2)),ABS(DT(3))) .GT. 0.02) THEN REDO = .TRUE. END IF END IF C IF (NITER .GE. 50) REDO = .FALSE. IF (REDO) GO TO 2000 IF ((NITER .LT. 50) .AND. (.NOT. CLIP)) THEN CLIP = .TRUE. DO I=1,3 DTOLD(I) = 0. CLAMP(I) = AMAX1(CLAMP(I), 0.25) END DO GO TO 2000 ELSE SHARP=1.4427*PAR(1)*PAR(2)*NUMER/(BRIGHT*SCALE*DENOM) SHARP=AMIN1(99.999,AMAX1(SHARP,-99.999)) RETURN END IF C END!