1
C===========================================================================
2
C Copyright (C) 1995-2011 European Southern Observatory (ESO)
4
C This program is free software; you can redistribute it and/or
5
C modify it under the terms of the GNU General Public License as
6
C published by the Free Software Foundation; either version 2 of
7
C the License, or (at your option) any later version.
9
C This program is distributed in the hope that it will be useful,
10
C but WITHOUT ANY WARRANTY; without even the implied warranty of
11
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12
C GNU General Public License for more details.
14
C You should have received a copy of the GNU General Public
15
C License along with this program; if not, write to the Free
16
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
19
C Correspondence concerning ESO-MIDAS should be addressed as follows:
20
C Internet e-mail: midas@eso.org
21
C Postal address: European Southern Observatory
22
C Data Management Division
23
C Karl-Schwarzschild-Strasse 2
24
C D 85748 Garching bei Muenchen
28
C===========================================================================
33
INTEGER MAXFRM, PSFMAX, NOPT, MAPSIZ
34
PARAMETER (MAXFRM=105000000, PSFMAX=35, NOPT=10)
35
C MAXFRM=105000000 => 10000 x 10500
37
C PARAMETER (MAXFRM=655360, PSFMAX=35, NOPT=10) => 640 x 1024
38
C MAXFRM=16810000 => 4100 x 4100
40
C MAXFRM represents the maximum image size the program is designed to
41
C tolerate. 655360 is large enough for 640 x 1024, and is also
42
C more than adequate for 800 x 800.
43
C PN 3/99 this is not true anymore, therefore MAXFRM increased
45
C PSFMAX is the maximum PSF radius. A value of 35 is consistent with
48
C NOPT is just the number of user-definable options (see below).
50
CHARACTER*26 LBL(NOPT)
51
REAL OPT(NOPT), OMIN(NOPT), OMAX(NOPT)
52
C REAL DATA(MAXFRM), SUBT(MAXFRM), SIGMA(MAXFRM)
53
C SAVE DATA, SUBT, SIGMA
55
CHARACTER*30 INPICT, OPTFIL, CASE
58
INTEGER NCOL, NROW, ISTAT
59
INTEGER*8 PNTRA, PNTRB, PNTRC
60
INTEGER MADRID(1) ! MIDAS
63
INCLUDE 'MID_INCLUDE:ST_DEF.INC' ! MIDAS
64
COMMON /VMR/ MADRID ! MIDAS
65
COMMON /SIZE/ NCOL, NROW
66
INCLUDE 'MID_INCLUDE:ST_DAT.INC' ! MIDAS
69
DATA LBL/' FITTING RADIUS', ! 1
70
. ' CE (CLIPPING EXPONENT)', ! 2
71
. ' REDETERMINE CENTROIDS', ! 3
72
. ' CR (CLIPPING RANGE)', ! 4
73
. ' WATCH PROGRESS', ! 5
74
. ' MAXIMUM GROUP SIZE', ! 6
75
. ' PERCENT ERROR (in %)', ! 7
76
. ' PROFILE ERROR (in %)', ! 8
77
. ' IS (INNER SKY RADIUS)', ! 9
78
. ' OS (OUTER SKY RADIUS)'/ ! 10
79
DATA OPT / 2.5, 6., 1., 2.5, 1., 50., 0.75, 5., 0., 0./
80
DATA OMIN / 1.6, 0.0, 0.0, 0., 0.0, 1.0, 0., 0., 0., 0./
81
DATA OMAX / 10., 8.0, 1.0, 10., 2.0, 100., 100., 100.,
84
C Set up the values of the optional parameters.
86
C Call OPTION, first with OPTFIL = 'allstar.opt' to set initial
87
C values for the optional parameters. If the file isn't there, the
88
C routine will check that the default values (specified in the data
89
C statement above) are valid, and return here with those values intact.
91
CALL STSPRO ('allstar') ! MIDAS
92
CALL STECNT('PUT', 1, 0, 0) ! MIDAS
96
CALL STFXMP(MAPSIZ,D_R4_FORMAT,PNTRA,ISTAT)
97
CALL STFXMP(MAPSIZ,D_R4_FORMAT,PNTRB,ISTAT)
98
CALL STFXMP(MAPSIZ,D_R4_FORMAT,PNTRC,ISTAT)
99
INPICT = ' ' ! avoid strange output
101
C This section show at start up a short message to the user informing him
102
C about the maximum size of the frames,
105
2'--------------------------------------------------------------',
110
CALL STTPUT(' ',ISTAT)
112
2'This version of ALLSTAR can be used for frames with a maximum ',
114
WRITE(TXTBUF,12345) MAPSIZ
115
12345 FORMAT('size of 10000 * 10500 (= ',I9,') pixels.')
116
CALL STTPUT(' '//TXTBUF,ISTAT)
118
2'In case your frame is larger than this value,', ISTAT)
120
2'ALLSTAR will not be able to process it.',ISTAT)
122
2'In such case, please contact your local MIDAS support person.',
125
2'--------------------------------------------------------------',
129
OPTFIL=CASE('allstar.opt')
130
CALL OPTION (OPTFIL, NOPT, LBL, OPT, OMIN, OMAX, 'OPT>', ISTAT)
132
IF (OPT(3) .GE. 0.5) THEN
140
1900 CALL GETNAM ('Input image name:', INPICT)
141
IF ((INPICT .EQ. 'END OF FILE') .OR. (INPICT .EQ. 'EXIT'))
143
CALL ATTACH (INPICT, OPEN)
148
IF (NCOL*NROW .GT. MAXFRM) THEN
149
CALL STUPID ('Picture is too large!')
151
69 FORMAT (/' Maximum allowed number of pixels =', I8)
155
CALL ALLSTR (MADRID(PNTRA), NCOL, NROW, MADRID(PNTRB),
157
. OPT(1), OPT(5), OPT(4), NINT(OPT(2)), CENTER,
158
. NINT(OPT(6)), PERERR, PROERR, OPT(9), OPT(10))
164
C#######################################################################
166
SUBROUTINE ALLSTR (DATA, NCOL, NROW, SUBT, SIGMA,
167
. FITRAD, WATCH, HALF, IEXP, CENTER,
168
. MAXGRP, PERERR, PROERR, SKYIN, SKYOUT)
170
C=======================================================================
172
C Photometry for many stars by simultaneous multiple PSF fits iterating
173
C the entire star list simultaneously.
175
C OFFICIAL DAO VERSION: 1987 June 30
179
C DATA is the input image, stored as a two-dimensional REAL
182
C NCOL, NROW are the dimensions of the input image.
184
C SUBT is the address of some working space where ALLSTAR can
185
C keep a scratch copy of the image.
187
C SIGMA is the address of some working space where ALLSTAR can
188
C keep track of the standard errors of the individual pixels
190
C FITRAD is the user-definable parameter "Fitting radius." It
191
C governs how many pixels out from the centroid of the star
192
C will actually be considered in computing the least-
193
C squares profile fits.
195
C WATCH is the 'watch progress' parameter specified by the user.
196
C If WATCH > 0, information relating to the progress of
197
C the reductions will be typed on the terminal during
200
C HALF, IEXP are the "Clipping radius" and "Clipping exponent"
201
C parameters which specify the degree to which discordant
202
C pixels are to be ignored.
204
C CENTER is a logical specifying whether the user wants the stellar
205
C centroids to be redetermined.
207
C MAXGRP is the largest group size for which the solution is to be
210
C=======================================================================
213
INTEGER MAXSTR, MAXPSF, MAXIT, MAXMAX, MAXPAR, MAXEXP
214
PARAMETER (MAXSTR=50000, MAXPSF=145, MAXIT=200, MAXMAX=100,
215
. MAXPAR=6, MAXEXP=6)
219
C MAXSTR The maximum number of stars in a frame.
221
C MAXPSF the largest PSF look-up table that can be accomodated. If
222
C PSFMAX is the largest acceptable PSF radius (see main-line
223
C program above), then MAXPSF = 2*[2*(PSFMAX+1)]+7.
225
C MAXMAX is the largest group for which a solution will ever be
226
C attempted = maximum permissible value of MAXGRP.
229
CHARACTER*30 COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL, SWITCH
230
CHARACTER LINE*80, SUBPIC*30, EXTEND*30, CASE*5
231
REAL C(3*MAXMAX+1,3*MAXMAX+1), V(3*MAXMAX+1), X(3*MAXMAX+1)
232
REAL DATA(NCOL,NROW), SUBT(NCOL,NROW), SIGMA(NCOL,NROW)
233
REAL PSF(MAXPSF,MAXPSF,MAXEXP), PAR(MAXPAR)
235
REAL XC(MAXSTR+1), YC(MAXSTR+1), MAG(MAXSTR+1), SKY(MAXSTR+1)
236
REAL RPIXSQ(MAXSTR), CHI(MAXSTR), SUMWT(MAXSTR)
237
REAL NUMER(MAXSTR), DENOM(MAXSTR)
238
REAL MAGERR(MAXSTR), DXOLD(MAXSTR), DYOLD(MAXSTR)
239
REAL XCLAMP(MAXSTR), YCLAMP(MAXSTR)
241
SAVE XC, YC, MAG, SKY
242
SAVE RPIXSQ, CHI, SUMWT
244
SAVE MAGERR, DXOLD, DYOLD
249
LOGICAL SKIP(MAXSTR), LAST(MAXSTR)
251
REAL AMIN1, AMAX1, ABS, SQRT, USEPSF
252
INTEGER MIN0, MAX0, INT, RDPSF
254
REAL LOBAD, HIBAD, ERR, SEP, THRESH, AP1, PHPADU, RONOIS, DUM
255
REAL FAINT, CLMPMX, WCRIT, RADIUS, XMIN, XMAX, YMIN, YMAX
256
REAL PERERR, PROERR, SKYIN, SKYOUT, SKYISQ, SKYOSQ, RELERR
257
REAL PSFMAG, BRIGHT, XPSF, YPSF, SEPCRT, SEPMIN, PSFRAD, PEAK
258
REAL FITRAD, WATCH, HALF, PKERR, PSFRSQ, RADSQ, RSQ, DX, DYSQ
259
REAL SIGSQ, WT, VAL, RHOSQ, DWT, SKYBAR, CHIGRP, SUMRES, GRPWT
260
REAL DPOS, DFDSIG, XKWT, DF, DIFF, SHARP, DELTAX, DELTAY, D, DY
262
INTEGER I, J, K, L, IFAINT, IDUM, IY, IX, NSTAR, NCONV, NITER
263
INTEGER MAXGRP, ITSKY, MINSKY, MAXUNK, LX, LY, ISTAT, INTRVL
264
INTEGER IPSTYP, NPSF, NPAR, NEXP, NFRAC, NL
265
INTEGER IEXP, IXMIN, IXMAX, IYMIN, IYMAX, MX, MY, NDISAP
266
INTEGER I3, J3, LSTAR, NSTR, NTERM, I3M2, ISTAR, N
267
LOGICAL CENTER, CLIP, OMIT, REDO
269
COMMON /FILNAM/ COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL
271
DATA SUMWT /MAXSTR*1./
273
C-----------------------------------------------------------------------
277
C Get ready, get set, . . .
280
SKYOSQ = AMIN1(SKYOUT**2, SKYISQ + (MAXSTR-1)/3.1415927)
281
SKYOUT = AMIN1(SKYOUT, SQRT(SKYOSQ))
284
MAXUNK=MAXMAX*3+1 ! Largest possible number of unknowns
285
CALL TBLANK ! Type a blank line
287
C Read the point-spread function into memory.
289
950 CALL GETNAM ('File with the PSF:', PSFFIL)
290
IF ((PSFFIL .EQ. 'END OF FILE') .OR.
291
. (PSFFIL .EQ. 'GIVE UP')) RETURN
293
ISTAT = RDPSF (PSFFIL, IPSTYP, PAR, MAXPAR, NPAR,
294
. PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC,
295
. PSFMAG, BRIGHT, XPSF, YPSF)
296
IF (ISTAT .NE. 0) THEN
300
PEAK = USEPSF(IPSTYP, 0., 0., BRIGHT, PAR, PSF, NPSF,
301
. NPAR, NEXP, NFRAC, 0., 0., DVDXC, DVDYC)
303
C Stars will be checked for merger if they are separated by less than
304
C 1 FWHM of the image core.
306
C Crit. sep. = 2.355*sigma, where
307
C sigma = SQRT [ (sigma(X)**2 + sigma(Y)**2)/2 ]
309
SEPCRT=2.*(PAR(1)**2+PAR(2)**2)
311
C Stars will be considered unconditionally merged if they are separated
312
C by less than about 0.375 * FWHM.
314
SEPMIN=AMIN1(1.,0.14*SEPCRT)
316
C SEPCRT contains the square of the critical separation.
317
C SEPMIN contains the square of the minimum separation.
319
PKERR=PROERR/(PAR(1)*PAR(2))**2 ! See fitting errors below
321
PSFRAD = (REAL(NPSF-1)/2.-1.)/2.
324
C Ascertain the name of the file with the input photometry, and open it.
326
960 CALL GETNAM ('Input file:', MAGFIL)
327
IF ((MAGFIL .EQ. 'END OF FILE') .OR.
328
. (MAGFIL .EQ. 'GIVE UP')) RETURN
330
CALL INFILE (2, MAGFIL, ISTAT)
331
IF (ISTAT .NE. 0) THEN
332
CALL STUPID ('Error opening input file '//MAGFIL)
337
CALL RDHEAD (2, NL, IDUM, IDUM, LOBAD, HIBAD, THRESH, AP1,
338
. PHPADU, RONOIS, DUM)
339
IF ((NL .LT. 1) .OR. (NL .GT. 3)) THEN
340
CALL STUPID ('Not a valid input file.')
346
C Inquire the name of the output file, and open it.
348
PROFIL=SWITCH(MAGFIL, CASE('.als'))
349
970 CALL GETNAM ('File for results:', PROFIL)
350
IF ((PROFIL .EQ. 'END OF FILE') .OR.
351
. (PROFIL .EQ. 'GIVE UP')) RETURN
353
PROFIL = EXTEND(PROFIL, CASE('als'))
354
CALL OUTFIL (1, PROFIL, ISTAT)
355
IF (ISTAT .NE. 0) THEN
356
CALL STUPID ('Error opening output file '//PROFIL)
360
CALL WRHEAD (1, 1, NCOL, NROW, 7, LOBAD, HIBAD, THRESH, AP1,
361
. PHPADU, RONOIS, FITRAD)
363
C Name for the output image.
365
SUBPIC=SWITCH(PROFIL, CASE('s'))
366
CALL GETNAM ('Name for subtracted image:', SUBPIC)
368
C We need to keep track of the anticipated standard error of the
369
C brightness value in a given pixel as the reductions proceed.
377
CALL RDARAY ('DATA', LX, LY, NCOL, L, NCOL, DATA(1,IY), ISTAT)
378
IF (ISTAT .NE. 0) THEN
379
CALL STUPID ('Error reading picture.')
384
IF ((DATA(IX,IY) .GT. HIBAD) .OR.
385
. (DATA(IX,IY) .LT. LOBAD)) THEN
386
SIGMA(IX,IY) = -1.1E38
388
SIGMA(IX,IY) = RONOIS
396
C Read in all the stars.
402
IF (I .GT. MAXSTR) THEN
403
CALL STUPID ('Too many stars.')
407
1120 CALL RDSTAR (2, NL, ID(I), XC(I), YC(I), MAG(I), SKY(I))
408
IF (ID(I) .LT. 0) GO TO 2000 ! End-of-file was encountered
409
IF (ID(I) .EQ. 0) GO TO 1120 ! A blank line was encountered
410
IF (MAG(I) .LT. FAINT) THEN
411
MAG(I)=10.**(0.4*(PSFMAG-MAG(I)))
421
C-----------------------------------------------------------------------
427
2000 NSTAR=I-1 ! Number of stars
430
5552 FORMAT (/1X, I5, ' stars. <<')
432
610 FORMAT (//' I = iteration number',
433
. //' R = number of stars that remain'
434
. //' D = number of stars that disappeared'
435
. //' C = number of stars that converged'
437
CALL OVRWRT (' I R D C', 1)
439
C Initialize accumulators and constraints on parameter corrections.
443
CALL STRIP (ID, XC, YC, MAG, SKY, SKIP, MAXSTR,
444
. NSTAR, NDISAP, SQRT(SEPMIN), NUMER, DENOM)
445
IF (NDISAP .GT. 0) THEN
446
WRITE (LINE,682) NITER, NSTAR, NDISAP, NCONV
447
CALL OVRWRT (LINE(1:28), 3)
451
CLIP = (IEXP .NE. 0) .AND. (NITER .GE. 4)
453
C Set up critical errors for star rejection.
456
IF (NITER .GE. 5) WCRIT = 1.0
457
IF (NITER .GE. 10) WCRIT = 1.5
458
IF (NITER .GE. 15) WCRIT = 2.0
460
C Sort stars by y-coordinate, for minimum paging.
462
CALL QUICK (YC, NSTAR, NUMER)
463
CALL IRECTFY (ID, NSTAR, NUMER, DENOM)
464
CALL RECTFY (XC, NSTAR, NUMER, DENOM)
465
CALL RECTFY (MAG, NSTAR, NUMER, DENOM)
466
CALL RECTFY (SKY, NSTAR, NUMER, DENOM)
467
CALL RECTFY (SUMWT, NSTAR, NUMER, DENOM)
468
CALL RECTFY (DXOLD, NSTAR, NUMER, DENOM)
469
CALL RECTFY (DYOLD, NSTAR, NUMER, DENOM)
470
CALL RECTFY (XCLAMP, NSTAR, NUMER, DENOM)
471
CALL RECTFY (YCLAMP, NSTAR, NUMER, DENOM)
473
C Determine how much of the image we need to copy into the working area.
475
C First determine the big rectangular area that contains all the pixels
476
C within one working radius of any star. The "working radius" is
477
C either the fitting radius or the outer sky radius (if this is an
478
C iteration during which the sky is to be determined), whichever is
479
C larger. Set their SIGMAs negative.
481
IF ((MOD(NITER, ITSKY) .EQ. 0) .AND. (SKYOUT .GT. FITRAD))
491
XMIN=AMIN1(XMIN, XC(I))
492
XMAX=AMAX1(XMAX, XC(I))
495
IXMIN=MAX0(1, MIN0(NCOL, INT(XMIN-RADIUS)+1))
496
IXMAX=MAX0(1, MIN0(NCOL, INT(XMAX+RADIUS)))
497
IYMIN=MAX0(1, MIN0(NROW, INT(YC(1)-RADIUS)+1))
498
IYMAX=MAX0(1, MIN0(NROW, INT(YC(NSTAR)+RADIUS)))
500
IF (WATCH .GT. 0.5) THEN
502
C Set up reporting intervals for row counter.
504
D=AMAX1(0.1, ALOG10(FLOAT(IYMAX-IYMIN+1)/10.))
508
INTRVL=10**(INTRVL+1)
509
ELSE IF (D .GT. 0.5) THEN
511
ELSE IF (D .GT. 0.2) THEN
518
C Beginning of big loop over y-coordinates.
523
DO 2190 IY=IYMIN,IYMAX
526
SIGMA(IX,IY) = -ABS(SIGMA(IX,IY))
529
C OK. Now find all the points within one working radius of each
530
C star. Set the sigmas of these pixels positive again and copy
531
C them from DATA to SUBT.
533
C N points at the last star that has such a low x-coordinate
534
C that it can be omitted.
539
IF (DY .GT. RADIUS) THEN
542
ELSE IF (DY .LT. -RADIUS) THEN
545
DYSQ = (REAL(IY) - YC(I))**2
546
LX=MAX0(1, MIN0(NCOL, INT(XC(I)-RADIUS)+1))
547
MX=MAX0(1, MIN0(NCOL, INT(XC(I)+RADIUS)))
549
DX = REAL(IX) - XC(I)
550
IF (DX**2+DYSQ .LE. RSQ) THEN
551
IF (SIGMA(IX,IY) .GT. -1.E38) THEN
552
SIGMA(IX,IY) = ABS(SIGMA(IX,IY))
553
SUBT(IX,IY) = DATA(IX,IY)
556
IF (DX .GE. 0.) GO TO 2120
562
C Subtract from the working copy of the image (SUBT) each star which
563
C has not yet converged by subtracting the shifted, scaled PSF
564
C according to the current best parameter estimates; do this only
565
C for pixels with positive SIGMAs, indicating that
566
C they will actually be needed for this iteration.
568
C L points at the last star that has such a low y-coordinate
569
C that it can be omitted.
574
IF (DY .GT. PSFRAD) THEN
577
ELSE IF (DY .LT. -PSFRAD) THEN
580
LX=MAX0(IXMIN, MIN0(IXMAX, INT(XC(I)-PSFRAD)+1))
581
MX=MAX0(IXMIN, MIN0(IXMAX, INT(XC(I)+PSFRAD)))
582
DELTAX=(XC(I)-1.)/XPSF-1.
583
DELTAY=(YC(I)-1.)/YPSF-1.
586
IF (SIGMA(IX,IY) .GT. 0.) THEN
588
IF (DX**2+DYSQ .LT. PSFRSQ) THEN
589
SUBT(IX,IY)=SUBT(IX,IY)-
590
. MAG(I)*USEPSF(IPSTYP, DX, DY, BRIGHT,
591
. PAR, PSF, NPSF, NPAR, NEXP, NFRAC,
592
. DELTAX, DELTAY, DVDXC, DVDYC)
594
IF (DX .GT. 0.) GO TO 2170
600
2180 IF ((WATCH .GT. 0.5) .AND.
601
. (INTRVL*(IY/INTRVL) .EQ. IY)) THEN
602
WRITE (LINE, 5551) IYMIN, IY, IYMAX
604
CALL OVRWRT (LINE(1:15), 2)
608
C If sky values are to be redetermined this iteration, do it now.
610
IF ((SKYOUT .GT. 0.5) .AND. (MOD(NITER, ITSKY) .EQ. 0))
612
DO 2198 ISTAR=1,NSTAR
613
LX=MAX0(IXMIN, MIN0(IXMAX, INT(XC(ISTAR)-SKYOUT)+1))
614
MX=MAX0(IXMIN, MIN0(IXMAX, INT(XC(ISTAR)+SKYOUT)))
615
LY=MAX0(IYMIN, MIN0(IYMAX, INT(YC(ISTAR)-SKYOUT)+1))
616
MY=MAX0(IYMIN, MIN0(IYMAX, INT(YC(ISTAR)+SKYOUT)))
620
DYSQ = (REAL(J) - YC(ISTAR))**2
622
DX = REAL(I) - XC(ISTAR)
624
IF (RSQ .GT. SKYOSQ) THEN
632
IF (RSQ .GE. SKYISQ) THEN
633
IF (SIGMA(I,J) .GT. -1.E38) THEN
636
IF (N .GE. MAXSTR) GO TO 2197
645
2197 IF (N .GT. MINSKY) THEN
646
CALL QUICK (NUMER, N, DENOM)
649
DO I=(N+1)/2-J,(N/2)+1+J
652
SKY(ISTAR) = DX/REAL((N/2)+2*J+2-(N+1)/2)
657
C Create provisional star groups.
659
CALL REGRP (ID, XC, YC, MAG, SKY, SUMWT, DXOLD, DYOLD,
660
. XCLAMP, YCLAMP, MAXSTR, NSTAR, FITRAD, LAST, NUMER,
663
C Now get ready to do the next iteration, group by group.
665
IF (WATCH .GT. 0.5) THEN
667
C Set up reporting intervals for star counter.
669
D=AMAX1(0.1, ALOG10(FLOAT(NSTAR)/10.))
673
INTRVL=10**(INTRVL+1)
674
ELSE IF (D .GT. 0.5) THEN
676
ELSE IF (D .GT. 0.2) THEN
686
C Find the last star in the current group.
688
DO LSTAR = ISTAR,NSTAR
689
IF (LAST(LSTAR)) GO TO 2210
692
2210 NSTR=LSTAR-ISTAR+1
694
C Start crunching on this group.
696
IF (NSTR .GT. MAXGRP) THEN
697
IF (WATCH .GT. 1.5) THEN
698
WRITE (LINE,61) NSTR, RADIUS
699
61 FORMAT ('Group too large:', I5, ' (', F5.2, ')')
700
CALL OVRWRT (LINE(1:30), 3)
703
IF (CENTER .AND. (NITER .GE. 2)) THEN
704
IF (RADIUS .LT. 1.2) THEN
705
IF (WATCH .GT. 1.5) THEN
707
62 FORMAT ('Group too dense to reduce.')
708
CALL OVRWRT (LINE(1:26), 3)
711
C Mark the faintest star in this group for deletion.
715
IF (MAG(I) .LT. FAINT) THEN
720
IF (WATCH .GT. 1.5) THEN
721
WRITE (6,63) ID(L), XC(L), YC(L),
722
. PSFMAG-1.085736*ALOG(MAG(L)), ' group too dense'
723
63 FORMAT (1X, I5, 2F9.3, F9.3, 9X, A, F6.3)
725
SKIP(L)=.TRUE. ! Flag for deletion
730
IF (RADIUS .LT. 0.8) THEN
731
IF (WATCH .GT. 1.5) THEN
733
CALL OVRWRT (LINE(1:26), 3)
736
C Mark the faintest star in this group for deletion.
740
IF (MAG(I) .LT. FAINT) THEN
745
IF (WATCH .GT. 1.5) WRITE (6,63) ID(L), XC(L), YC(L),
746
. PSFMAG-1.085736*ALOG(MAG(L)), ' group too dense'
747
SKIP(L)=.TRUE. ! Flag for deletion
753
CALL REGRP (ID(ISTAR), XC(ISTAR), YC(ISTAR), MAG(ISTAR),
754
. SKY(ISTAR), SUMWT(ISTAR), DXOLD(ISTAR), DYOLD(ISTAR),
755
. XCLAMP(ISTAR), YCLAMP(ISTAR), L, NSTR, RADIUS,
756
. LAST(ISTAR), NUMER(ISTAR), DENOM(ISTAR))
766
C Get the extent of the rectangular area of the frame relevant to this
767
C group. At the same time, determine the average SKY and the average
768
C CHI value for the group, and zero the various statistical
769
C accumulators. What we would REALLY like to do is to compute the
770
C anticipated error for each pixel the same way as PEAK and NSTAR do.
771
C We can't really do that, because stellar groups aren't defined the
772
C same way as in NSTAR: here they change as stars move from one group
773
C to another, or converge and are subtracted from the frame. Therefore,
774
C there is no consistent way to define a group's CHI value the same as
775
C was done in NSTAR. I will take the cheap way out: using the
776
C individual smoothed CHI values for the stars which were obtained
777
C during the previous iteration (which should still be stored in the
778
C array SUMWT) I will simply average those values and call the mean the
779
C group's effective smoothed CHI. This way, at least, the reductions
780
C will be the same as before for isolated stars. (Note that,
781
C for the first iteration, the SUMWT's have been set to 1.0 by a DATA
787
XMIN = AMIN1(XMIN, XC(I))
788
XMAX = AMAX1(XMAX, XC(I))
789
YMIN = AMIN1(YMIN, YC(I))
790
YMAX = AMAX1(YMAX, YC(I))
791
SKYBAR = SKYBAR+SKY(I)
792
CHIGRP = CHIGRP+SUMWT(I)
799
SKYBAR = SKYBAR/FLOAT(NSTR)
800
CHIGRP = CHIGRP/FLOAT(NSTR)
801
IF (CENTER .AND. (NITER .GE. 2)) THEN
807
C If sky is to be determined: NTERM=NTERM+1
809
C Now... on with the iteration.
811
IXMIN=MIN0(NCOL, MAX0(1, INT(XMIN-FITRAD)+1))
812
IXMAX=MIN0(NCOL, MAX0(1, INT(XMAX+FITRAD)))
813
IYMIN=MIN0(NROW, MAX0(1, INT(YMIN-FITRAD)+1))
814
IYMAX=MIN0(NROW, MAX0(1, INT(YMAX+FITRAD)))
816
C IXMIN, IXMAX, IYMIN, and IYMAX are now the limits of a rectangular
817
C array containing all pixels within one fitting radius of any star in
820
C Zero the normal matrix and the vector of residuals.
829
C Now deal with the pixels one by one.
831
DO 2390 IY=IYMIN,IYMAX
832
DO 2380 IX=IXMIN,IXMAX
833
IF (SIGMA(IX,IY) .LE. 0.) GO TO 2380
835
C If this pixel is within one fitting radius of at least one star in
836
C the current group, include it in the calculation. Otherwise, skip
837
C it. While figuring this out, compute the squared distance of this
838
C pixel from the centroid of each star in the group.
842
RPIXSQ(I)=(FLOAT(IX)-XC(I))**2+(FLOAT(IY)-YC(I))**2
843
IF (RPIXSQ(I) .LT. RADSQ) THEN
847
IF (OMIT) GO TO 2380 ! Do not need this pixel
850
IF (RPIXSQ(I) .LT. RADSQ) THEN
857
C The expected random error in the pixel is the quadratic sum of
858
C the Poisson statistics, plus the readout noise, plus an estimated
859
C error of 0.75% of the total brightness for the difficulty of flat-
860
C fielding and bias-correcting the chip, plus an estimated error of
861
C some fraction of the fourth derivative at the peak of the profile,
862
C to account for the difficulty of accurately interpolating within the
863
C point-spread function. The fourth derivative of the PSF is
864
C proportional to H/sigma**4 (sigma is the Gaussian width parameter for
865
C the stellar core); using the geometric mean of sigma(x) and sigma(y),
866
C this becomes H/[sigma(x)*sigma(y)]**2. The ratio of the fitting
867
C error to this quantity is estimated from a good-seeing CTIO frame to
868
C be approximately 0.027 (see definition of PKERR above.)
870
D=SUBT(IX,IY)-SKYBAR ! Residual of this pixel
871
DPOS=AMAX1(0., DATA(IX,IY)-D)
873
C DPOS = raw data minus residual
874
C = model-predicted brightness in this pixel, consisting of sky
875
C plus all stellar profiles, which presumably is non-negative.
877
C The four error sources in our noise model are:
880
C (3) Flat-field errors
881
C (4) Errors in the PSF
883
C Numerically, the squares of these quantities are
884
C (1) RONOIS = SIGMA(IX,IY) (initially, at least)
885
C (2) DPOS/PHPADU where DPOS = sum of the stellar profiles plus sky
886
C (3) constant x DPOS**2
887
C (4) constant x the sum of (stellar profile)**2
889
C Here's the thing: at this point we don't have the
890
C sum of (stellar profile)**2 available to us. All we have is
891
C (sum of stellar profile)**2, which isn't the same thing at all.
892
C Nevertheless, we have to assume that
894
C sum of (stellar profile)**2 = (sum of stellar profile)**2;
896
C otherwise I'd have to use a whole new REAL array as big as the
897
C original CCD image, and I don't want to do that. This should be
898
C good enough for pixels where either (a) the other three error sources
899
C dominate, .OR. (b) one of the stars putting light into that pixel is
900
C much brighter than all the rest. Thus, the approximation should be
901
C good UNLESS the pixel is near the center of two or more bright
902
C stars --- in which case things will be hairy anyway.
904
SIGSQ = SIGMA(IX,IY) + DPOS/PHPADU + (PERERR*DPOS)**2 +
905
. (PKERR*(DPOS-SKYBAR))**2
906
RELERR=ABS(D)/SQRT(SIGSQ)
907
IF (CLIP .AND. (RELERR .GT. 100.)) GO TO 2380
910
C Now include this pixel in the fitting equations for the group.
912
DO 2320 I=ISTAR,LSTAR
913
IF (SKIP(I)) GO TO 2320
915
IF (RSQ .GE. 0.999999) GO TO 2320 ! Safety check vs divide by 0
916
WT=AMAX1(WT, 5./(5.+RSQ/(1.-RSQ)))
918
C The condition equation for pixel (IX,IY) is of the form
920
C data(IX,IY)-summation{scale*psf(IX-Xcenter,IY-Ycenter)}-sky=residual
922
C Then we will jigger the scale's, Xcenter's, and Ycenter's such that
924
C Summation{weight * residual**2}
926
C is minimized. 'weight' will be a function (1) of the distance of this
927
C pixel from the center of the nearest star, (2) of the model-predicted
928
C brightness of the pixel (taking into consideration the readout noise,
929
C the photons/ADU, and the interpolation error of the PSF), and (3) of
930
C the size of the residual itself. (1) is necessary to prevent the
931
C non-linear least-squares solution from oscillating: oft-times it will
932
C come to pass that if you include a pixel in the solution, then the
933
C predicted shift of the centroid will cause that pixel to be excluded
934
C in the next iteration, and the new predicted shift of the centroid
935
C will cause that pixel to be included again. This could go on ad
936
C infinitum. The cure is to have the weight of a pixel go
937
C continuously to zero as its distance from the stellar centroid
938
C approaches the fitting radius. In a case like that just described,
939
C the solution can then find a real minimum of the sum of the
940
C weighted squared residuals with that pixel at some low-weight position
941
C just inside the fitting radius. (2) is just sensible weighting.
942
C (3) is just a crude attempt at making the solution more robust against
945
VAL = USEPSF(IPSTYP, FLOAT(IX)-XC(I), FLOAT(IY)-YC(I), BRIGHT,
946
. PAR, PSF, NPSF, NPAR, NEXP, NFRAC, (XC(I)-1.)/XPSF-1.,
947
. (YC(I)-1.)/YPSF-1., DVDXC, DVDYC)
948
IF (NTERM .GT. NSTR) THEN
959
RHOSQ=((XC(I)-FLOAT(IX))/PAR(1))**2+
960
. ((YC(I)-FLOAT(IY))/PAR(2))**2
961
IF (RHOSQ .LE. 36.) THEN
962
RHOSQ = 0.6931472*RHOSQ
963
DFDSIG = EXP(-RHOSQ)*(RHOSQ-1.)
964
NUMER(I) = NUMER(I)+DFDSIG*D/SIGSQ
965
DENOM(I) = DENOM(I)+DFDSIG**2/SIGSQ
969
C At this point, the vector X contains the first derivative of
970
C the condition equation for pixel (IX,IY) with respect to each of
971
C the fitting parameters for all of the stars. Now these derivatives
972
C will be added into the normal matrix and the vector of residuals.
974
C Add this residual into the weighted sum of the absolute relative
981
C SUMRES is the weighted sum of [ABS(residual)/sigma] for all the
982
C pixels in the group. Now also add the weighted value of
983
C [ABS(residual)/sigma] into the accumulating sum for each of the
986
DO 2330 I=ISTAR,LSTAR
987
IF (SKIP(I)) GO TO 2330
993
C Up until now, WT represents only the radial weighting profile. Now
994
C figure in the anticipated standard error of the pixel. Reject any
995
C pixel with a 100-sigma residual after iteration 2.
998
IF (CLIP) WT=WT/(1.+(RELERR/(CHIGRP*HALF))**IEXP)
1000
C Now work this pixel into the normal matrix.
1003
C If sky is to be determined: C(NTERM,NTERM)=C(NTERM,NTERM)+WT
1004
C If sky is to be determined: V(NTERM)=V(NTERM)-DWT
1005
DO 2370 I=ISTAR,LSTAR
1006
IF (SKIP(I)) GO TO 2370
1007
IF (NTERM .GT. NSTR) THEN
1011
C If sky is to be determined: C(NTERM,K)=C(NTERM,K)-X(K)*WT
1012
2340 V(K)=V(K)+X(K)*DWT
1014
IF (SKIP(J)) GO TO 2360
1018
DO 2350 L=J3-2,MIN0(K, J3)
1019
2350 C(K,L)=C(K,L)+X(L)*XKWT
1025
C If sky is to be determined: C(NTERM,K)=C(NTERM,K)-XKWT
1027
IF (SKIP(J)) GO TO 2365
1029
C(K,L)=C(K,L)+X(L)*XKWT
1037
C Reflect the normal matrix across the diagonal.
1039
IF (NTERM .GT. 1) THEN
1045
C Compute the estimate of the standard deviation of the residuals for
1046
C the group as a whole, and for each star. This estimate starts out as
1047
C SQRT(PI/2)*{SUM[weight*ABS(residual/sigma)]/SUM(weight)} and then gets
1048
C corrected for bias by SQRT(no. of pixels/(no. of pixels - degrees of
1051
IF (GRPWT .GT. NTERM) THEN
1052
CHIGRP=1.2533141*SUMRES/SQRT(GRPWT*(GRPWT-NTERM))
1054
C But then I drive the value toward unity, depending on exactly how
1055
C many pixels were involved: if CHI is based on exactly a total
1056
C weight of 3, then it is extremely poorly determined, and we just
1057
C want to keep CHI = 1. The larger GRPWT is, the better determined
1058
C CHI is, and the less we want to force it toward unity. So,
1059
C just take the weighted average of CHI and unity, with weights
1060
C GRPWT-3 and 1, respectively.
1062
CHIGRP=((GRPWT-3.)*CHIGRP+3.)/GRPWT
1067
C CHIGRP has been pulled toward its expected value of unity to keep the
1068
C statistics of a small number of pixels from compeletely dominating
1069
C the error analysis. Similarly, the photometric errors for the
1070
C individual stars will be pulled toward unity now. Later on, if the
1071
C number of stars in the group is greater than one, CHI(I) will be
1072
C nudged toward the group average. In order to work optimally, of
1073
C course, this requires that PHPADU, RONOIS, and the other noise
1074
C contributors which I have postulated properly represent the true
1075
C errors expected in each pixel.
1077
C At the same time, be sure that every star in the group contains at
1078
C least 3 valid pixels if recentroiding is being performed, 1 valid
1079
C pixel if not. If any star in the group fails to meet this criterion
1080
C mark that star for deletion, set REDO to true, and skip ahead to the
1085
IF (CENTER .AND. (NITER .GE. 2)) THEN
1086
IF (NPIX(I) .LT. 3) THEN
1088
IF (WATCH .GT. 1.5) WRITE (6,63) ID(I), XC(I), YC(I),
1089
. PSFMAG-1.085736*ALOG(MAG(I)),
1090
. ' too few valid pixels'
1091
SKIP(I)=.TRUE. ! Flag for deletion
1094
SKIP(I)=.FALSE. ! Flag for retention
1096
IF (SUMWT(I) .GT. 3.) THEN
1097
CHI(I)=1.2533141*CHI(I)/SQRT(SUMWT(I)*(SUMWT(I)-3.))
1099
C Store a smoothed CHI value for the star in SUMWT.
1101
SUMWT(I) = ((SUMWT(I)-3.)*CHI(I) + 3.)/SUMWT(I)
1107
IF (NPIX(I) .LT. 1) THEN
1109
SKIP(I)=.TRUE. ! Flag for deletion
1111
IF (WATCH .GT. 1.5) WRITE (6,63) ID(I), XC(I), YC(I),
1112
. PSFMAG-1.085736*ALOG(MAG(I)),
1113
. ' too few valid pixels'
1115
SKIP(I)=.FALSE. ! Flag for retention
1117
IF (SUMWT(I) .GT. 1.) THEN
1118
CHI(I)=1.2533141*CHI(I)/SQRT(SUMWT(I)*(SUMWT(I)-1.))
1120
C Store a smoothed CHI value for the star in SUMWT.
1122
SUMWT(I) = ((SUMWT(I)-3.)*CHI(I) + 3.)/SUMWT(I)
1132
CALL INVERS (C, MAXUNK, NTERM, ISTAT)
1134
IF (C(J,J) .LE. 0.) THEN
1136
C Uh-oh. Zero on the diagonal. Booma Booma!
1138
IF (CENTER .AND. (NITER .GE. 2)) THEN
1144
C I now points at the (first, at least) star that caused the problem.
1148
IF (WATCH .GT. 1.5) WRITE (6,63) ID(I), XC(I), YC(I),
1149
. PSFMAG-1.085736*ALOG(MAG(I)), ' singular matrix'
1153
CALL VMUL (C, MAXUNK, NTERM, V, X)
1154
C If sky is to be determined: SKYBAR=SKYBAR-X(NTERM)
1155
C If sky is to be determined: IF(ABS(X(NTERM)).GT.0.01)REDO=.TRUE.
1157
C In the beginning, the brightness of each star will be permitted to
1158
C change by no more than two magnitudes per iteration, and the x,y
1159
C coordinates of each centroid will be permitted to change by no more
1160
C than 0.4 pixel per iteration. Any time that the parameter
1161
C correction changes sign from one iteration to the next, the maximum
1162
C permissible change will be reduced by a factor of two. These
1163
C clamps are released any time a star in the group disappears.
1165
DO 2520 I=ISTAR,LSTAR
1166
IF (CENTER .AND. (NITER .GE. 2)) THEN
1168
C Correct both magnitude and position.
1174
C Note that the sign of the correction is such that it must be
1175
C SUBTRACTED from the current value of the parameter to get the
1176
C improved parameter value. This being the case, if the correction
1177
C to the brightness is negative (the least-squares thinks that the
1178
C star should be brighter) a change of 1 magnitude is a change of a
1179
C factor of 2.5; if the brightness correction is positive (the star
1180
C should be fainter) a change of 1 magnitude is a change of 60%.
1183
IF (DWT .LT. 0.) THEN
1184
XCLAMP(I)=AMAX1(0.001,0.5*XCLAMP(I))
1186
XCLAMP(I)=AMIN1(CLMPMX,1.2*XCLAMP(I))
1188
XC(I)=XC(I)-X(K)/(1.+ABS(X(K)/XCLAMP(I)))
1191
IF (DWT .LT. 0.) THEN
1192
YCLAMP(I)=AMAX1(0.001,0.5*YCLAMP(I))
1194
YCLAMP(I)=AMIN1(CLMPMX,1.2*YCLAMP(I))
1196
YC(I)=YC(I)-X(L)/(1.+ABS(X(L)/YCLAMP(I)))
1199
. (1.+AMAX1( X(J)/(0.84*MAG(I)) , -X(J)/(5.25*MAG(I)) ))
1200
MAGERR(I) = SUMWT(I)*SQRT(C(J,J))
1201
IF (NITER .GE. 4) THEN
1204
. AMAX1( 0.1*MAGERR(I) , 0.0005*MAG(I) )) THEN
1207
DF = (0.1*SUMWT(I))**2
1208
IF (X(K)**2 .GT. AMAX1(DF*C(K,K), 4.E-6)) THEN
1210
ELSE IF (X(L)**2 .GT. AMAX1(DF*C(L,L), 4.E-6)) THEN
1220
C Correct magnitude only. Since this is an easy, linear problem,
1221
C elaborate clamps are not needed. A simple one will do to keep
1222
C the brightness from ever going negative.
1224
MAG(I)=MAG(I)-X(J)/(1. + 1.2*ABS( X(J)/MAG(I) ))
1225
MAGERR(I) = SUMWT(I)*SQRT(C(J,J))
1226
IF (NITER .GE. 2) THEN
1229
. AMAX1( 0.1*MAGERR(I), 0.0005*MAG(I) )) REDO = .TRUE.
1235
C No star will be allowed to converge with a signal-to-noise ratio
1236
C worse than 2.0. Such stars must be retained until either they get
1237
C better than S/N = 2, or they get eliminated some time after
1238
C iteration 5 (if S/N < 1.0), 10 (if S/N < 1.5), or 15 (if S/N < 2.0).
1239
C However, anything which is lucky enough to survive until the last
1240
C iteration will be written out regardless.
1242
IF (MAG(I) .LT. 2.0*MAGERR(I)) REDO = .TRUE.
1243
IF (NITER .GE. MAXIT) REDO=.FALSE.
1245
C If this star converged, write out the results for it,
1246
C flag it for deletion from the star list (SKIP = .TRUE.), and
1247
C subtract it from the reference copy of the image (DATA).
1249
IF (.NOT. REDO) THEN
1252
C Subtract from the reference copy of the image (DATA).
1254
LX=MAX0(0, MIN0(NCOL, INT(XC(I)-PSFRAD)))+1
1255
MX=MAX0(1, MIN0(NCOL, INT(XC(I)+PSFRAD)))
1256
LY=MAX0(0, MIN0(NROW, INT(YC(I)-PSFRAD)))+1
1257
MY=MAX0(1, MIN0(NROW, INT(YC(I)+PSFRAD)))
1258
DELTAX=(XC(I)-1.)/XPSF-1.
1259
DELTAY=(YC(I)-1.)/YPSF-1.
1264
IF (SIGMA(K,J) .GT. -1.E38) THEN
1266
IF (DX**2+DYSQ .GE. PSFRSQ) GO TO 2195
1267
DIFF=MAG(I)*USEPSF(IPSTYP, DX, DY, BRIGHT, PAR,
1268
. PSF, NPSF, NPAR, NEXP, NFRAC, DELTAX, DELTAY,
1270
DATA(K,J)=DATA(K,J)-DIFF
1272
C We must add this star's contribution to the pixel's noise into the
1273
C noise map SIGMA (which, of course, actually contains sigma**2).
1274
C Please note that of the four contributors to the variance
1276
C (1) RONOIS = SIGMA(IX,IY) (initially, at least)
1277
C (2) DPOS/PHPADU where DPOS = sum of the stellar profiles plus sky
1278
C (3) constant x DPOS**2
1279
C (4) constant x the sum of (the stellar profiles)**2
1281
C (1) is already included in SIGMA, and we can calculate this star's
1282
C contribution to (2) and (4) directly. (3) is tricker. (3) is
1283
C (sum of sky plus all stellar profiles)**2, and we need to determine
1284
C the contribution of (this stellar profile). We need to use algebra:
1286
C (sum of sky and all OTHER stellar profiles + this stellar profile)**2
1288
C = (sum of sky plus all OTHER stellar profiles)**2 +
1289
C 2 x (sum of sky plus OTHER stellar profiles) x (this stellar profile)
1290
C + (this stellar profile)**2.
1294
C 2 x (sum of sky plus OTHER stellar profiles) x (this stellar profile)
1295
C + (this stellar profile)**2
1297
C into SIGMA. That way, later on when the solution has forgotten that
1298
C this star ever existed, so it only adds
1300
C (sum of sky plus all OTHER stellar profiles)**2
1302
C into the noise expected for the pixel, it will get the noise right.
1303
C Since we have already subtracted this star's contribution to the image
1304
C (array DATA) above, I will use the current DATA as my working
1305
C approximation for (sum of sky plus all OTHER stellar profiles).
1306
C I really shouldn't do this; I should really be calculating each star's
1307
C model-predicted contribution to this pixel, but that would take too
1310
IF (DIFF .GT. 0.) THEN
1311
DIFF = DIFF/PHPADU +
1312
. PERERR*2.*AMAX1(0., DATA(K,J))*DIFF
1314
SIGMA(K,J)=SIGN(ABS(SIGMA(K,J))+DIFF, SIGMA(K,J))
1317
2195 CONTINUE ! End of loops over pixels
1318
SHARP=1.4427*PAR(1)*PAR(2)*NUMER(I)/(MAG(I)*PEAK*DENOM(I))
1319
SHARP=AMIN1( 99.999, AMAX1( SHARP, -99.999 ))
1320
ERR=1.085736*MAGERR(I)/MAG(I)
1321
MAG(I)=PSFMAG-1.085736*ALOG(MAG(I))
1322
WRITE (1,321) ID(I), XC(I), YC(I), MAG(I), ERR, SKY(I),
1323
. FLOAT(NITER), CHI(I), SHARP
1324
321 FORMAT (I6, 5F9.3, F9.0, F9.2, F9.3)
1325
SKIP(I)=.TRUE. ! Remove from star list
1329
C If there is more than one star remaining in this group, check to see
1330
C whether any two of them have merged. This means, find the closest
1331
C pair and see whether they are TOO close.
1335
IF (NSTR .GT. 1) THEN
1336
DO 2230 I=ISTAR+1,LSTAR
1337
IF (SKIP(I)) GO TO 2230
1340
IF (SKIP(J)) GO TO 2220
1341
SEP=(XC(I)-XC(J))**2+(YC(I)-YC(J))**2
1342
IF (SEP .GE. RSQ) GO TO 2220
1344
C Two stars are overlapping. Identify the fainter of the two.
1347
IF (MAG(I) .LT. MAG(J)) THEN
1358
C No two stars have merged if K still equals zero.
1360
IF (K .LE. 0) GO TO 2260
1362
C The K-th star is now the fainter of the two, the L-th, the brighter.
1364
C Now eliminate the fainter of the two if they are TOO close.
1366
IF ((RSQ .GT. SEPMIN) .AND.
1367
. (MAG(K) .GT. WCRIT*MAGERR(K))) GO TO 2260
1369
C Now replace the centroid of the L-th star with the weighted mean of
1370
C the most recent estimates of the centroids of the L-th and K-th
1371
C stars, and the brightness of the L-th with the sum of the brightnesses
1372
C of the L-th and K-th.
1374
XC(L)=XC(L)*MAG(L)+XC(K)*MAG(K)
1375
YC(L)=YC(L)*MAG(L)+YC(K)*MAG(K)
1376
MAG(L)=MAG(L)+MAG(K)
1380
C Remove the K-th star from the group.
1382
SKIP(K)=.TRUE. ! Remove from star list
1383
IF (WATCH .GT. 1.5) THEN
1384
WRITE (6,63) ID(K), XC(K), YC(K),
1385
. PSFMAG-1.085736*ALOG(MAG(K)), ' blended with'
1386
WRITE (6,63) ID(L), XC(L), YC(L),
1387
. PSFMAG-1.085736*ALOG(MAG(L))
1391
C Loosen the clamps of every star in the group.
1396
XCLAMP(I)=AMAX1(0.5*CLMPMX, XCLAMP(I))
1397
YCLAMP(I)=AMAX1(0.5*CLMPMX, YCLAMP(I))
1401
C If the number of iterations completed is less than or equal to 3,
1402
C perform another iteration no questions asked.
1404
2260 IF (NITER .LE. 3) GO TO 3000
1406
C > IF NO STAR HAS BEEN REMOVED FROM THIS GROUP DURING THIS ITERATION <
1408
C Check whether any of the stars is too faint (more than 12.5
1409
C magnitudes fainter than the PSF star). If several stars are
1410
C too faint, delete the faintest one, and set the brightnesses of
1411
C the other faint ones exactly to 12.5 mag below the PSF star.
1412
C That way on the next iteration we will see whether these stars
1413
C want to grow or to disappear.
1419
IF (SKIP(I)) GO TO 3000
1420
IF (MAG(I) .LT. FAINT) THEN
1424
IF (MAG(I) .LT. 1.E-5) MAG(I)=1.E-5
1427
C If at least one star is more than 12.5 mag. fainter than the
1428
C PSF, then IFAINT is the index of the faintest of them, and FAINT
1429
C is the relative brightness of the faintest of them.
1431
IF (IFAINT .GT. 0) THEN
1432
SKIP(IFAINT)=.TRUE. ! Remove from star list
1434
IF (WATCH .GT. 1.5) WRITE (6,63) ID(IFAINT), XC(IFAINT),
1435
. YC(IFAINT), PSFMAG-1.085736*ALOG(MAG(IFAINT)),
1438
C Loosen the clamps of every star in the group.
1443
XCLAMP(I)=AMAX1(0.5*CLMPMX, XCLAMP(I))
1444
YCLAMP(I)=AMAX1(0.5*CLMPMX, YCLAMP(I))
1446
ELSE IF (NITER .GE. 5) THEN
1448
C If no star in this group is more than 12.5 mag. fainter than the PSF,
1449
C then after the fifth iteration delete the least certain star if it
1450
C is less than a one-sigma detection; after the tenth iteration delete
1451
C the least certain star if it is less than a 1.5-sigma detection;
1452
C after the fifteenth iteration delete the least certain star if it is
1453
C less than a two-sigma detection.
1458
DO 2550 I=ISTAR,LSTAR
1460
IF (WT .LT. FAINT) THEN
1466
IF (FAINT .LT. WCRIT) THEN
1467
SKIP(IFAINT)=.TRUE. ! Remove from star list
1468
IF (WATCH .GT. 1.5) WRITE (6,63) ID(IFAINT), XC(IFAINT),
1469
. YC(IFAINT), PSFMAG-1.085736*ALOG(MAG(IFAINT)),
1470
. ' error too big: ',
1471
. 1.085736*MAGERR(IFAINT)/MAG(IFAINT)
1474
C Loosen the clamps of every star in the group.
1479
XCLAMP(I)=AMAX1(0.5*CLMPMX, XCLAMP(I))
1480
YCLAMP(I)=AMAX1(0.5*CLMPMX, YCLAMP(I))
1485
IF (WATCH .GT. 0.5) THEN
1486
IF (LSTAR/INTRVL .GT. (ISTAR-1)/INTRVL) THEN
1487
WRITE (LINE,682) NITER, INTRVL*(LSTAR/INTRVL),
1489
CALL OVRWRT (LINE(1:24), 2)
1493
IF (ISTAR .LE. NSTAR) GO TO 2200 ! Go on to next group
1495
C We've gone through the entire star list for this iteration.
1497
C Find the last star in the list that still needs more work
1498
C (SKIP=.FALSE.). If there aren't any, then we're done! If there is
1499
C one, point NSTAR at it.
1503
IF (SKIP(NSTAR)) THEN
1505
IF (NSTAR .GT. 0) THEN
1508
WRITE (LINE,682) NITER, NSTAR, NDISAP, NCONV
1509
CALL OVRWRT (LINE(1:28), 3)
1510
GO TO 9000 ! None left!! We can go now
1514
C Well, at this point we've established that there are at least some
1515
C stars in the list that still need more work. Remove any stars with
1516
C which we are finished (SKIP = .TRUE.) by overwriting the first such
1517
C star with the star at position NSTAR in the list (we have just
1518
C confirmed that this is the last star in the list that still needs
1519
C work). Then decrement NSTAR and go back up to 3005 to find the new
1520
C last star in the list that still needs work.
1523
IF (ISTAR .GE. NSTAR) THEN
1524
WRITE (LINE,682) NITER, NSTAR, NDISAP, NCONV
1525
682 FORMAT (4I6, ' <<')
1526
CALL OVRWRT (LINE(1:28), 3)
1527
GO TO 2100 ! The list is ready. Start another iteration
1530
IF (SKIP(ISTAR)) THEN
1534
MAG(ISTAR)=MAG(NSTAR)
1535
SKY(ISTAR)=SKY(NSTAR)
1536
SUMWT(ISTAR)=SUMWT(NSTAR)
1537
DXOLD(ISTAR)=DXOLD(NSTAR)
1538
DYOLD(ISTAR)=DYOLD(NSTAR)
1539
XCLAMP(ISTAR)=XCLAMP(NSTAR)
1540
YCLAMP(ISTAR)=YCLAMP(NSTAR)
1548
C-----------------------------------------------------------------------
1555
C Copy the input picture verbatim into the output picture.
1557
IF (SUBPIC(1:11) .NE. 'END OF FILE') THEN
1558
CALL COPPIC (SUBPIC, SUBT, NCOL, NROW, ISTAT)
1559
IF (ISTAT .NE. 0) THEN
1560
CALL STUPID ('Error opening output picture.')
1565
CALL WRARAY ('COPY', LX, LY, NCOL, NROW, NCOL, DATA, ISTAT)
1566
IF (ISTAT .NE. 0) THEN
1567
CALL STUPID ('Error writing picture.')
1572
CALL STUPID (' Done. ')
added
removed
contrib/daophot/src/allstar.for
