1
C @(#)peak.for 19.1 (ES0-DMD) 02/25/03 13:23:50
2
C===========================================================================
3
C Copyright (C) 1995 European Southern Observatory (ESO)
5
C This program is free software; you can redistribute it and/or
6
C modify it under the terms of the GNU General Public License as
7
C published by the Free Software Foundation; either version 2 of
8
C the License, or (at your option) any later version.
10
C This program is distributed in the hope that it will be useful,
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
C GNU General Public License for more details.
15
C You should have received a copy of the GNU General Public
16
C License along with this program; if not, write to the Free
17
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge,
20
C Corresponding concerning ESO-MIDAS should be addressed as follows:
21
C Internet e-mail: midas@eso.org
22
C Postal address: European Southern Observatory
23
C Data Management Division
24
C Karl-Schwarzschild-Strasse 2
25
C D 85748 Garching bei Muenchen
27
C===========================================================================
29
SUBROUTINE DAOPK (PAR, MAXPAR, PSF, MAXPSF, MAXEXP, F,
30
. WATCH, FITRAD, PERERR, PROERR)
32
C=======================================================================
34
C Single-star profile-fitting routine.
36
C OFFICIAL DAO VERSION: 1991 April 18
38
C This subroutine reads in an array around each star in an input
39
C photometry file and fits the point-spread function to the observed
40
C stellar profile by least-squares. The three parameters solved for
41
C are the brightness ratio between the PSF and the program star, and
42
C the coordinates of the centroid of the program star.
46
C WATCH (INPUT) governs whether information relating to the progress
47
C of the reductions is to be typed on the terminal screen
50
C FITRAD (INPUT) is the fitting radius. Only pixels within FITRAD of
51
C the current estimate of a star's centroid will be included in
52
C the least-squares determination of its centroid and magnitude.
54
C All are user-definable optional parameters.
56
C=======================================================================
59
INTEGER MAXPSF, MAXEXP, MAXBOX, MAXPAR
64
C MAXPSF is the length of the side of the largest PSF look-up table
65
C allowed. (Note: half-pixel grid size)
67
C MAXBOX is the length of the side of the largest box within which the
68
C PSF can be evaluated. (MAXBOX = (MAXPSF-7)/2 )
70
CHARACTER*30 COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL, SWITCH
71
CHARACTER LINE*80, EXTEND*30, CASE*3
73
REAL PSF(MAXPSF,MAXPSF,MAXEXP), PAR(MAXPAR)
76
REAL LOBAD, AMAG, SHARP, CHI, ERRMAG, DELTAX, DELTAY
77
REAL APMAG, RADIUS, FRAD, READNS, AP1, THRESH, HIBAD
78
REAL PROERR, PERERR, FITRAD, WATCH, PHPADU, RONOIS
80
REAL X, Y, PSFMAG, BRIGHT, XPSF, YPSF, SCALE, SKY
81
INTEGER RDPSF, NITER, IST, LX, LY, NX, NY
82
INTEGER NFRAC, IBEG, ISTAR, NEXP, NPAR, NPSF, IPSTYP
83
INTEGER IDUM, NL, ISTAT, NCOL, NROW, NTERM
85
COMMON /SIZE/ NCOL, NROW
86
COMMON /FILNAM/ COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL
87
COMMON /ERROR/ PHPADU, RONOIS, PERR, PKERR
89
C-----------------------------------------------------------------------
93
C Ascertain the name of the file containing coordinates and magnitude
94
C estimates for the program stars. Open it and read the header.
97
950 CALL GETNAM ('File with aperture results:', MAGFIL)
98
IF ((MAGFIL .EQ. 'END OF FILE') .OR.
99
. (MAGFIL .EQ. 'GIVE UP')) THEN
103
CALL INFILE (2, MAGFIL, ISTAT)
104
IF (ISTAT .NE. 0) THEN
109
CALL RDHEAD (2, NL, IDUM, IDUM, LOBAD, HIBAD, THRESH, AP1,
110
. PHPADU, READNS, FRAD)
113
C Ascertain the name of the PSF file, open it, and read the PSF.
115
900 CALL GETNAM ('File with the PSF:', PSFFIL)
116
IF ((PSFFIL .EQ. 'END OF FILE') .OR.
117
. (PSFFIL .EQ. 'GIVE UP')) THEN
123
ISTAT = RDPSF (PSFFIL, IPSTYP, PAR, MAXPAR, NPAR,
124
. PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC,
125
. PSFMAG, BRIGHT, XPSF, YPSF)
126
IF (ISTAT .NE. 0) THEN
133
PKERR = 0.01*PROERR/(PAR(1)*PAR(2)) ! **2
135
C Get ready to do the PEAK fitting: get the name of the output file.
136
C Open it and write the header.
138
PROFIL=SWITCH(MAGFIL, CASE('.pk'))
139
960 CALL GETNAM ('File for PEAK results:', PROFIL)
140
IF ((PROFIL .EQ. 'END OF FILE') .OR.
141
. (PROFIL .EQ. 'GIVE UP')) THEN
146
PROFIL = EXTEND(PROFIL, CASE('pk'))
147
CALL OUTFIL (3, PROFIL, ISTAT)
148
IF (ISTAT .NE. 0) THEN
149
CALL STUPID ('Error opening output file '//PROFIL)
154
RADIUS = AMIN1( FITRAD, (REAL(NPSF-1)/2. - 1.)/2. )
155
CALL WRHEAD (3, 1, NCOL, NROW, 7, LOBAD, HIBAD, THRESH, AP1,
156
. PHPADU, READNS, RADIUS)
157
IF (WATCH .GT. 0.5) CALL TBLANK ! Type a blank line
160
C-----------------------------------------------------------------------
164
C Reduce the stars, one by one.
166
C Read data for the next star and initialize things.
168
2000 CALL RDSTAR (2, NL, ISTAR, X, Y, APMAG, SKY)
169
IF (ISTAR .LT. 0) GO TO 9000 ! End-of-file encountered
170
IF (ISTAR .EQ. 0) GO TO 2000 ! Blank line encountered
171
IF (APMAG .GE. 99.) APMAG=PSFMAG+5. ! The ol' college try
172
DELTAX=(X-1.)/XPSF-1.
173
DELTAY=(Y-1.)/YPSF-1.
174
LX=MAX0(1, INT(X-RADIUS)+1)
175
LY=MAX0(1, INT(Y-RADIUS)+1)
176
NX=MIN0(NCOL, INT(X+RADIUS)) - LX + 1
177
NY=MIN0(NROW, INT(Y+RADIUS)) - LY + 1
178
CALL RDARAY ('DATA', LX, LY, NX, NY, MAXBOX, F, IST)
180
C At this point LX, LY are the coordinates in the big picture of the
181
C first pixel in the subarray. The dimensions of the subarray are
187
C X,Y are now the coordinates of the star's centroid in the subframe
188
C (where (x,y)=(1.0,1.0) are the coordinates of the center of the first
189
C pixel and (x,y)=(NX,NY) are the coordinates of the center of the last
190
C pixel in the subarray).
192
C Display the subarray on the terminal, if desired, and type out headers
195
IF (WATCH .GT. 1.5) CALL SHOW (F, F(NINT(X),NINT(Y)), 0.9*SKY,
197
IF ((WATCH .GT. 1.5) .OR. ((WATCH .GT. 0.5) .AND. (IBEG .EQ. 1)))
199
620 FORMAT (9X, 'STAR', 6X, 'X', 8X, 'Y', 8X, 'MAGNITUDE', 7X,
202
SCALE=10.**(-0.4*(APMAG-PSFMAG)) ! Starting value
203
CALL PKFIT (F, NX, NY, MAXBOX, X, Y, SCALE, SKY, RADIUS, LOBAD,
204
. HIBAD, BRIGHT, IPSTYP, PAR, MAXPAR, NPAR,
205
. PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC,
206
. DELTAX, DELTAY, ERRMAG, CHI, SHARP, NITER)
207
IF (NITER .GT. 0) GO TO 2010
209
C A singular matrix occurred during the least-squares solution.
211
WRITE (LINE,621) ISTAR
212
621 FORMAT (8X, I5, ' had a singular matrix.')
216
C Everything went fine.
220
AMAG=PSFMAG-2.5*ALOG10(SCALE)
221
ERRMAG=AMIN1(2.0, 1.086*ERRMAG/SCALE)
222
IF (WATCH .GT. 0.5) WRITE (6,622) ISTAR, X, Y, AMAG, ERRMAG,
224
622 FORMAT (8X, I5, 2F9.2, F9.3, ' +- ', F5.3, 2F7.2, I5)
225
WRITE (3,320) ISTAR, X, Y, AMAG, ERRMAG, SKY, FLOAT(NITER), CHI,
227
320 FORMAT (I6, 3F9.3, F9.4, F9.3, F9.0, F9.2, F9.3)
230
C-----------------------------------------------------------------------
236
CALL STUPID (' Done.')
241
C#######################################################################
243
SUBROUTINE PKFIT (F, NX, NY, MAXBOX, X, Y, SCALE, SKY, RADIUS,
244
. LOBAD, HIBAD, BRIGHT, IPSTYP, PAR, MAXPAR, NPAR,
245
. PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC,
246
. DELTAX, DELTAY, ERRMAG, CHI, SHARP, NITER)
248
C=======================================================================
250
C This is the subroutine which does the actual one-star least-squares
251
C profile fit for PEAK.
253
C OFFICIAL DAO VERSION: 1991 April 6
256
C F (INPUT) is an NX by NY array containing actual picture data.
258
C MAXBOX (INPUT) is the maximum value allowable for either NX or NY,
259
C needed for the dimension statements below. PEAK and PSF will
260
C provide different values of MAXBOX.
262
C SCALE (INPUT/OUTPUT) is the initial estimate of the brightness of
263
C the star, expressed as a fraction of the brightness of the
264
C PSF. Upon return, the final computed value of SCALE will
265
C be passed back to the calling routine.
267
C X, Y (INPUT/OUTPUT) are the initial estimates of the centroid of
268
C the star relative to the corner (1,1) of the subarray. Upon
269
C return, the final computed values of X and Y will be passed
270
C back to the calling routine.
272
C SKY (INPUT) is the local sky brightness value, carried on from
273
C PHOTOMETRY via the data files.
275
C RADIUS (INPUT) is the fitting radius-- only pixels within RADIUS of
276
C the instantaneous estimate of the star's centroid will be
277
C included in the fit.
279
C LOBAD and HIBAD (INPUT) are bad pixel limits-- any pixel whose
280
C brightness value falls outside this range will be presumed to
281
C be bad, and will be ignored.
283
C BRIGHT (INPUT) contains the brightness normalization of the model
284
C analytic point spread function
286
C PAR (INPUT) contains the values of the remaining NPARAM parameters
287
C defining the analytic function which approximates the core of
290
C PSF (INPUT) is an NPSF by NPSF by NEXP+NFRAC look-up table
291
C containing corrections from the analytic approximation of the
292
C PSF to the true PSF.
294
C ERRMAG (OUTPUT) is the estimated standard error of the value of SCALE
295
C returned by this routine.
297
C CHI (OUTPUT) is the estimated goodness-of-fit statistic: the ratio
298
C of the observed pixel-to-pixel mean absolute deviation from
299
C the profile fit, to the value expected on the basis of the
300
C read-out noise and the photons/ADU (which are brought in
301
C through COMMON block /ERROR/).
303
C SHARP (OUTPUT) is a goodness-of-fit statistic describing how much
304
C broader the actual profile of the object appears than the
305
C profile of the PSF.
307
C NITER (OUTPUT) is the number of iterations the solution required to
308
C achieve convergence. If NITER = 50, the solution did not
309
C converge. If for some reason a singular matrix occurs during
310
C the least-squares solution, this will be flagged by setting
313
C=======================================================================
316
INTEGER MAXPSF, MAXEXP, MAXBOX, MAXPAR
320
C MAXPSF is the length of the side of the largest PSF look-up table
321
C allowed. (Note: half-pixel grid spacing)
323
REAL C(3,3), V(3), CLAMP(3), DTOLD(3)
324
REAL F(MAXBOX,MAXBOX), T(3), DT(3), NUMER
325
REAL PSF(MAXPSF,MAXPSF,MAXEXP), PAR(MAXPAR)
326
REAL USEPSF, AMAX1, AMIN1, ABS
329
REAL LOBAD, DF, WT, DFDSIG, RHOSQ, RELERR, SIG, SIGSQ
330
REAL FPOS, DVDXC, DVDYC, RSQ, DX, DXSQ, DY, DYSQ, DATUM
331
REAL DENOM, SUMWT, CHIOLD, RADSQ, SHARP, CHI, ERRMAG
332
REAL DELTAX, DELTAY, X, Y, PHPADU, RONOIS, PERR, PKERR
333
REAL SCALE, SKY, RADIUS, HIBAD, BRIGHT
334
INTEGER I, J, IX, IY, IXLO, IXHI, IYLO, IYHI, ISTAT
335
INTEGER NPIX, NFRAC, NEXP, NITER, NPAR, NPSF, IPSTYP
336
INTEGER NX, NY, NCOL, NROW
339
COMMON /SIZE/ NCOL, NROW
340
COMMON /ERROR/ PHPADU, RONOIS, PERR, PKERR
342
C-----------------------------------------------------------------------
344
C Initialize a few things for the solution.
355
C-----------------------------------------------------------------------
357
C Here begins the big least-squares loop.
361
C Initialize things for this iteration. CHI and SHARP will be
362
C goodness-of-fit indices. CHI will also be used in determining the
363
C weights of the individual pixels. As the solution iterates, the new
364
C value of CHI will be built up in the variable CHI, while a smoothed
365
C value of CHI computed from the previous iteration will be carried
373
V(I)=0.0 ! Zero the vector of residuals
375
2010 C(I,J)=0.0 ! Zero the normal matrix
377
C Choose the little box containing points inside the fitting radius.
379
IXLO=MAX0(1, INT(X-RADIUS))
380
IYLO=MAX0(1, INT(Y-RADIUS))
381
IXHI=MIN0(NX, INT(X+RADIUS)+1)
382
IYHI=MIN0(NY, INT(Y+RADIUS)+1)
384
C Now build up the normal matrix and vector of residuals.
392
IF ((DATUM .LT. LOBAD) .OR. (DATUM .GT. HIBAD)) GO TO 2090
396
C DX and DY are the distance of this pixel from the centroid of the
397
C star. Is this pixel truly inside the fitting radius?
399
RSQ=(DXSQ+DYSQ)/RADSQ
400
IF (1.-RSQ .LE. 2.E-6) GO TO 2090 ! Prevents floating underflows
402
C The fitting equation is of the form
404
C Observed brightness =
405
C [SCALE + delta(SCALE)] * [PSF + delta(Xcen)*d(PSF)/d(Xcen) +
406
C delta(Ycen)*d(PSF)/d(Ycen) ]
408
C and is solved for the unknowns delta(SCALE) ( = the correction to
409
C the brightness ratio between the program star and the PSF) and
410
C delta(Xcen) and delta(Ycen) ( = corrections to the program star's
413
C The point-spread function is equal to the sum of the integral under
414
C a two-dimensional analytic profile plus values interpolated from
417
T(1)=USEPSF(IPSTYP, DX, DY, BRIGHT, PAR, PSF, NPSF, NPAR, NEXP,
418
. NFRAC, DELTAX, DELTAY, DVDXC, DVDYC)
419
IF ((SCALE*T(1)+SKY .GT. HIBAD) .AND. (NITER.GE.3)) GO TO 2090
422
DF=F(IX,IY)-SKY-SCALE*T(1)
424
C DF is the residual of the brightness in this pixel from the PSF fit.
426
C The expected random error in the pixel is the quadratic sum of
427
C the Poisson statistics, plus the readout noise, plus an estimated
428
C error of PERERR% of the total brightness for the difficulty of flat-
429
C fielding and bias-correcting the chip, plus an estimated error of
430
C some fraction of the fourth derivative at the peak of the profile,
431
C to account for the difficulty of accurately interpolating within the
432
C point-spread function. The fourth derivative of the PSF is roughly
433
C proportional to H/sigma**4 (sigma is the width parameter for
434
C the stellar core); using the geometric mean of sigma(x) and sigma(y),
435
C this becomes H/[sigma(x)*sigma(y)]**2. The ratio of the fitting
436
C error to this quantity is estimated from a good-seeing CTIO frame to
437
C be approximately 0.027. NOWADAYS, PAR(1) and PAR(2) are the HWHM,
438
C not the sigma, so the coefficient is now of order 0.05 = 5.0%.
440
FPOS=AMAX1(0., F(IX,IY)-DF)
442
C FPOS = raw data minus residual = model-predicted value of the
443
C intensity at this point (which presumably is non-negative).
445
SIGSQ=FPOS/PHPADU+RONOIS+(PERR*FPOS)**2+(PKERR*(FPOS-SKY))**2
449
C SIG is the anticipated standard error of the intensity in this pixel,
450
C including readout noise, Poisson photon statistics, and an estimate
451
C of the standard error of interpolating within the PSF.
453
WT=5./(5.+RSQ/(1.-RSQ)) ! Weight as function of radius
455
C Now add this pixel into the quantities which go to make up the SHARP
458
RHOSQ=DXSQ/PAR(1)**2+DYSQ/PAR(2)**2
460
C Include in the sharpness index only those pixels within six
461
C HWHMs of the centroid of the star. (This saves time and
462
C floating underflows by excluding pixels which contribute less than
463
C about one part in a million to the index.)
465
IF (RHOSQ .LE. 36.) THEN
466
RHOSQ=0.6931472*RHOSQ
467
DFDSIG=EXP(-RHOSQ)*(RHOSQ-1.)
468
FPOS=AMAX1(0., F(IX,IY)-SKY)+SKY
469
NUMER=NUMER+DFDSIG*DF/SIGSQ
470
DENOM=DENOM+DFDSIG**2/SIGSQ
473
C Derive the weight of this pixel. First of all, the weight depends
474
C upon the distance of the pixel from the centroid of the star-- it
475
C is determined from a function which is very nearly unity for radii
476
C much smaller than the fitting radius, and which goes to zero for
477
C radii very near the fitting radius. Then reject any pixels with
478
C 10-sigma residuals (after the first iteration).
483
C Now the weight is scaled to the inverse square of the expected mean
488
C Reduce the weight of a bad pixel. A pixel having a residual of 2.5
489
C sigma gets reduced to half weight; a pixel having a residual of 5.
490
C sigma gets weight 1/257.
492
IF (CLIP) WT=WT/(1.+(0.4*RELERR/CHIOLD)**8)
494
C Now add the pixel into the vector of residuals and the normal matrix.
499
2030 C(I,J)=C(I,J)+T(I)*T(J)*WT
502
2090 CONTINUE ! End of loop over pixels
504
C Compute the (robust) goodness-of-fit index CHI.
506
IF (SUMWT .GT. 3.) THEN
507
CHI=1.2533141*CHI/SQRT(SUMWT*(SUMWT-3.))
509
C CHI is pulled toward its expected value of unity before being stored
510
C in CHIOLD to keep the statistics of a small number of pixels from
511
C completely dominating the error analysis.
513
CHIOLD=((SUMWT-3.)*CHI+3.)/SUMWT
519
C Compute the parameter corrections and check for convergence.
521
IF (NPIX .LT. 3) THEN
525
CALL INVERS (C, 3, 3, ISTAT)
526
IF (ISTAT .EQ. 0) GO TO 2100
530
C Everything OK so far.
532
2100 CALL VMUL (C, 3, 3, V, DT)
534
C In the beginning, the brightness of the star will not be permitted
535
C to change by more than two magnitudes per iteration (that is to say,
536
C if the estimate is getting brighter, it may not get brighter by
537
C more than 525% per iteration, and if it is getting fainter, it may
538
C not get fainter by more than 84% per iteration). The x and y
539
C coordinates of the centroid will be allowed to change by no more
540
C than one-half pixel per iteration. Any time that a parameter
541
C correction changes sign, the maximum permissible change in that
542
C parameter will be reduced by a factor of 2.
545
IF (DTOLD(I)*DT(I) .LT. 0.) THEN
546
CLAMP(I)=0.5*CLAMP(I)
548
CLAMP(I)=MIN(1., 1.1*CLAMP(I))
553
. (1.+AMAX1(DT(1)/(5.25*SCALE),-DT(1)/(0.84*SCALE))/CLAMP(1))
554
X=X+DT(2)/(1.+ABS(DT(2))/(0.4*CLAMP(2)))
555
Y=Y+DT(3)/(1.+ABS(DT(3))/(0.4*CLAMP(3)))
556
IF (NITER .LE. 1) GO TO 2000
559
C Convergence criteria: if the most recent computed correction to the
560
C brightness is larger than 0.01% or than 0.05 * sigma(brightness),
561
C whichever is larger, OR if the absolute change in X or Y is
562
C greater than 0.001 pixels, convergence has not been achieved.
564
ERRMAG=CHIOLD*SQRT(C(1,1))
567
. AMAX1( 0.05*ERRMAG, 0.0001*SCALE )) THEN
569
ELSE IF (AMAX1(ABS(DT(2)),ABS(DT(3))) .GT. 0.001) THEN
574
. AMAX1( ERRMAG, 0.002*SCALE )) THEN
576
ELSE IF (AMAX1(ABS(DT(2)),ABS(DT(3))) .GT. 0.02) THEN
581
IF (NITER .GE. 50) REDO = .FALSE.
583
IF ((NITER .LT. 50) .AND. (.NOT. CLIP)) THEN
587
CLAMP(I) = AMAX1(CLAMP(I), 0.25)
591
SHARP=1.4427*PAR(1)*PAR(2)*NUMER/(BRIGHT*SCALE*DENOM)
592
SHARP=AMIN1(99.999,AMAX1(SHARP,-99.999))
added
removed
contrib/daophot/libsrc/peak.for
