c photometry parameter (n1=4096) parameter (n2=4032) PARAMETER (MAXEXP=10, NOPT=20, MAXPSF=207) parameter (MAX=MAXPSF*MAXPSF*MAXEXP/2) parameter (MAXAP=12) real par(MAXAP+2) real w1(MAX) integer w2(MAX) real opt(nopt) real a_b(n1*n2) logical logi character*60 f1,f2,f3 character*80 head(36) common /SIZE/ ncol,nrow,corner,corner1 k=iargc() if(k.lt.1)then write(*,*) write(*,*) write(*,*)' ********** DAO_PHOT 2006.3 ***********' write(*,*) write(*,*) write(*,*)' Usage: u_phot4k file.fit [3]' write(*,*) write(*,*)' 3 edge_lines' write(*,*) stop endif call getarg(1,f1) corner=3. if(k.gt.1)then call getarg(2,f2) read(f2(1:),*)corner endif write(*,*)'set edge line: ~',corner k=index(f1,'.') if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit' write(*,'(a)')f1 inquire(file=f1,exist=logi) if(.not.logi)stop ' file not found !' open(1,file=f1,status='old',access='direct',recl=2880) read(1,rec=1)head ipos=indexpos(head,'NAXIS ') if(ipos.lt.36)read(head(ipos)(65:77),'(i5,i8)')n12,n49 if(n12.ne.0)then b_high=490000. if(n49.ne.0)b_high=n49 else b_high=60000. endif close(1) k=index(f1,'.') ! 2 for read f2=f1(1:k)//'coo' ! 3 for write open(2,file=f2,status='UNKNOWN') f2=f1(1:k)//'ap' open(3,file=f2,status='UNKNOWN') f3='bok.par' inquire(file=f3,exist=logi) if(.not.logi)f3='/vega2/rhbin/bok.par' open(1,file=f3,status='old') do 20 i=1,nopt 21 read(1,'(a)')f3 if(f3(1:2).eq."**")goto 21 20 read(f3(1:),'(29x,f9.2)')opt(i) if(opt(4).gt.b_high)b_high=opt(4) opt(4)=b_high c write(*,'(7f10.2)')opt do 30 i=1,MAXAP+2 31 read(1,'(a)')f3 if(f3(1:2).eq."**")goto 31 read(f3(1:),'(29x,f9.2)')par(i) 30 par(i)=par(i) write(*,'(7f10.2)')par close(1) call readfits(f1,a_b,ncol,nrow) corner1=ncol-corner CALL PHOTSB (a_b,w1,w2, MAX, OPT(4), OPT(11), par) end SUBROUTINE PHOTSB (F, SKY, INDEX, MAXSKY, INBAD, WATCH, par) PARAMETER (MINSKY=20, MAXAP=12) c IMPLICIT NONE INTEGER MINSKY, MAXSKY, MAXAP, NCOL, NROW C CHARACTER*80 LINE DOUBLE PRECISION MAGLIM, MAGSQ, WT, SUMWT, DLOG10, DBLE, DMAX1 DOUBLE PRECISION APMAG(MAXAP), AREA(MAXAP) REAL SKY(MAXSKY), MAGERR(MAXAP) REAL F(NCOL,NROW), ERROR(3) REAL PAR(MAXAP+2) REAL AMIN1, AMAX1, SQRT INTEGER INDEX(MAXSKY) INTEGER MIN0, MAX0 C REAL LOBAD, SKYMOD, SKYSIG, SKYSKW, SIGSQ, SKYVAR REAL DATUM, R, RSQ, FRACTN, EDGE, HIBAD, INBAD, THRESH, DUM REAL PHPADU, READNS, XC, YC, DMAG, WATCH, APMXSQ REAL RINSQ, ROUT, ROUTSQ, DYSQ INTEGER I, J, K, NAPER, IDUM, NMAG, LX, LY INTEGER MX, MY, NSKY, ISTAR, NL common /SIZE/ ncol,nrow,corner,corner1 C write(*,*)'HIGH_GOOD_DATUM= ',INBAD NAPER=MAXAP APMXSQ=-1. kcorner=0 DO 1010 I=1,MAXAP IF (PAR(I) .LE. 0.0) GO TO 1020 1010 APMXSQ=AMAX1(APMXSQ, (PAR(I)+0.5)**2) GO TO 1030 C 1020 NAPER=I-1 1030 CONTINUE C C NAPER is the number of apertures, whose radii are stored in C elements 1 through NAPER of the array PAR. C C APMXSQ is the outermost edge of the largest aperture-- if the C distance squared of the center of a pixel from the centroid of C the star is greater than APMXSQ, then we know that no part C of the pixel is to be included in any aperture. C C Now define the other variables whose values are in the table. C RINSQ=AMAX1(PAR(MAXAP+1), 0.)**2 ! Inner sky radius squared ROUTSQ = REAL(MAXSKY)/3.142 + RINSQ DUM = PAR(MAXAP+2)**2 IF (DUM .GT. ROUTSQ) THEN write(*,*)' *** You have specified too big a sky annulus. ***' WRITE (6,6) SQRT(ROUTSQ) 6 FORMAT (F10.2, ' pixels is the largest outer sky radius ', . 'currently permitted.') RETURN ELSE IF (DUM .LE. RINSQ) THEN write(*,*)'outer sky radius is no bigger than the inner radius.' return ELSE ROUT = PAR(MAXAP+2) ROUTSQ = DUM END IF C C Inquire the name of the input data file with the stellar positions, C and open it. C NL = -1 LOBAD = -10. HIBAD = INBAD write(*,*)hibad c CALL RDHEAD(2, NL, IDUM, IDUM, LOBAD, HIBAD, THRESH, DUM, c . PHPADU, READNS, R) READ(2,'(a)')dum READ(2,*)NL,IDUM,IDUM,LOBAD,HIBAD,THRESH,DUM,PHPADU,READNS,R READ(2,'(a)')dum IF (NL .LT. 1) NL=1 C C Inquire file name for output aperture photometry results, and open C the new file. C CALL WRHEAD (3, 2, NCOL, NROW, 7, LOBAD, HIBAD, THRESH, . PAR(1), PHPADU, READNS, R) C C If progress is being monitored, type out column headers. C IF (WATCH .GT. 0.5) WRITE (6,610) 610 FORMAT (/13X, 'STAR', 5X, 'X', 7X, 'Y', 9X, 'MAG.(1)', 8X, 'SKY') C C Initialize variables for the computation of the magnitude limit. C MAGLIM=0.0D0 MAGSQ=0.0D0 SUMWT=0.0D0 NMAG=0 C C----------------------------------------------------------------------- C C SECTION 2 C C Derive aperture photometry object by object. C C Get the coordinates of next object to be measured. C c LX = 1 c LY = 1 c NX = NCOL c NY = NROW c CALL RDARAY ('DATA', LX, LY, NX, NY, NCOL, F, ISTAT) 2000 CALL RDSTAR (2, NL, ISTAR, XC, YC, DMAG, DUM) IF (ISTAR .LT. 0) GO TO 9000 IF (ISTAR .EQ. 0) GO TO 2000 if(xc .lt. corner .or. xc .gt. corner1)goto 2000 if(yc .lt. corner .or. yc .gt. corner1)goto 2000 C C Compute the limits of the submatrix. C LX = MAX0(1, INT(XC-ROUT)+1) MX = MIN0(NCOL, INT(XC+ROUT)) LY = MAX0(1, INT(YC-ROUT)+1) MY = MIN0(NROW, INT(YC+ROUT)) EDGE=AMIN1(XC-0.5, (NCOL+0.5)-XC, YC-0.5, (NROW+0.5)-YC) C C EDGE is the distance of the star's centroid from the outermost C extremum of the array. C C Initialize star counts and aperture area. C DO 2010 I=1,NAPER APMAG(I) = 0.D0 C C If this star aperture extends outside the array, the magnitude C in this aperture will be no good. C IF (EDGE .LT. PAR(I)) APMAG(I)=-1.1D30 ! Null magnitude 2010 AREA(I)=0.0D0 C C Now read through the submatrix, picking out the data we want. C NSKY=0 C DO 2130 J=LY,MY DYSQ=(J-YC)**2 C DO 2125 I=LX,MX RSQ=DYSQ+(I-XC)**2 DATUM=F(I,J) C C Is this pixel within the sky annulus? C IF ((RSQ .LT. RINSQ) .OR. (RSQ .GT. ROUTSQ) .OR. . (NSKY .GT. MAXSKY) .OR. (DATUM .LT. LOBAD) .OR. . (DATUM .GT. HIBAD)) GO TO 2110 NSKY=NSKY+1 SKY(NSKY)=DATUM C C The inclusion of partial pixels inside the aperture is done as C follows: if the distance of the center of the current pixel from the C centroid of the star [radius vector r(i,j)] is exactly equal to the C radius of the aperture [R(k)], then one-half of the counts in the C pixel are included. If r(i,j) < R(k)-0.5, then the entire pixel is C included, while if r(i,j) > R(k)+0.5, the pixel is wholly excluded. C In between, viz. for R(k)-0.5 < r(i,j) < R(k)+0.5, the fraction of C the counts included varies linearly. Therefore a circular aperture C is approximated by an irregular (not even convex) polygon. C C If this pixel falls completely outside the LARGEST aperture, go on C to the next pixel. Notice that APMXSQ has actually been defined C as (R(k)+0.5)**2 for the largest value of R(k), in accordance with C the formula used for the partial pixels. C 2110 CONTINUE IF (RSQ .GT. APMXSQ) GO TO 2125 R=SQRT(RSQ)-0.5 C DO 2120 K=1,NAPER C C If this pixel falls completely outside THIS aperture, go on to the C next aperture. C IF (R .GT. PAR(K)) GO TO 2120 FRACTN=AMAX1(0.0, AMIN1(1.0,PAR(K)-R)) C C FRACTN is the fraction of the pixel that falls inside the C (irregular) aperture. C C If the pixel is bad, set the total counts in this aperture to a number C so negative that it will never be positive again. C IF ((DATUM .GE. LOBAD) .AND. (DATUM .LE. HIBAD) .AND. . (APMAG(K) .GT. -1.0D20)) THEN APMAG(K) = APMAG(K)+DBLE(FRACTN*DATUM) AREA(K) = AREA(K)+DBLE(FRACTN) ELSE IF (DATUM .GT. HIBAD) THEN APMAG(K) = -1.1D38 ELSE APMAG(K) = -1.0D30 END IF 2120 CONTINUE 2125 CONTINUE C 2130 CONTINUE C C We have accumulated the brightnesses of individual sky pixels in the C one-dimensional array SKY. Pixels falling above or below the BAD C limits have already been eliminated. Now sort SKY to place the C pixels in order of increasing brightness. C IF (NSKY .LT. MINSKY) THEN write(*,*)istar, * 'There aren''t enough pixels in the sky annulus.' do i=1,naper APMAG(I) = 99.999D0 MAGERR(I) = 9.999 enddo goto 2221 END IF CALL QUICK (SKY, NSKY, INDEX) C C Obtain the mode, standard deviation, and skewness of the peak in the C sky histogram. C CALL MMM (SKY, NSKY, HIBAD, READNS, . DUM, DATUM, SKYMOD, SKYSIG, SKYSKW) SKYVAR=SKYSIG**2 SIGSQ=SKYVAR/FLOAT(NSKY) C C SKYMOD has units of (ADU/pixel), and SKYSIG is the pixel-to-pixel C scatter of SKYMOD, in units of (ADU/pixel). SKYVAR is the C variance (square of the standard deviation) of the sky brightness, C (ADU/pixel)**2, and SIGSQ is the square of the standard error of the C mean sky brightness. C C Subtract the sky from the integrated brightnesses in the apertures, C convert the results to magnitudes, and compute standard errors. C DO 2220 I=1,NAPER IF (APMAG(I) .LT. -1.D31) THEN APMAG(I) = 99.999D0 MAGERR(I) = 9.999 GO TO 2220 END IF C C If the modal sky value could not be determined, set the magnitude C to 94.999. C IF (SKYSIG .LT. -0.5) GO TO 2210 APMAG(I)=APMAG(I)-DBLE(SKYMOD)*AREA(I) C C If the star + sky is fainter than the sky, or if the star aperture C extends beyond the limits of the picture, or if there is a bad pixel C in the star aperture, set the magnitude to 94.999. C IF (APMAG(I) .LE. 0.0D0) GO TO 2210 ERROR(1)=SNGL(AREA(I))*SKYVAR ERROR(2)=SNGL(APMAG(I))/PHPADU ERROR(3)=SIGSQ*SNGL(AREA(I))**2 C C These variables ERRORn are the respective variances (squares of the C mean errors) for: (1) random noise inside the star aperture, including C readout noise and the degree of contamination by other stars in the C neighborhood, as estimated by the scatter in the sky values (this C standard error increases as the square root of the area of the C aperture); (2) the Poisson statistics of the observed star brightness; C (3) the uncertainty of the mean sky brightness (this standard error C increases directly with the area of the aperture). C MAGERR(I)=AMIN1(9.999, . 1.0857*SQRT(ERROR(1)+ERROR(2)+ERROR(3))/SNGL(APMAG(I))) APMAG(I)=25.D0-2.5D0*DLOG10(APMAG(I)) IF (APMAG(I) .GT. 94.999D0) GO TO 2210 GO TO 2220 2210 APMAG(I)=94.999D0 MAGERR(I)=9.999 2220 CONTINUE 2221 continue C C Write out the answers. C IF (WATCH .GT. 0.5) WRITE (6,620)ISTAR, XC, YC, APMAG(1), . MAGERR(1), SKYMOD 620 FORMAT (12X, I5, 2F8.2, F9.3, ' +-', F6.3, F8.1) WRITE (3,320) ISTAR, XC, YC, (APMAG(I), I=1,NAPER) 320 FORMAT (/I6, 14F9.3) WRITE (3,321) SKYMOD, AMIN1(999.99,SKYSIG), . AMIN1(999.99, AMAX1(-99.99,SKYSKW)), . (MAGERR(I), I=1,NAPER) c 321 FORMAT (4X, F9.3, 2F6.2, F8.4, 11F9.4) 321 FORMAT (4X, F9.2,1x, 2F6.2, F8.4, 11F9.4) kcorner=kcorner+1 C IF (APMAG(1) .GT. 99.D0) GO TO 2000 WT=(2./(2.-DMAG))*(100./MAGERR(1))**2 MAGLIM=MAGLIM+WT*(SNGL(APMAG(1))-DMAG) MAGSQ=MAGSQ+WT*(SNGL(APMAG(1))-DMAG)**2 SUMWT=SUMWT+WT NMAG=NMAG+1 GO TO 2000 C C----------------------------------------------------------------------- C C Normal return. C C Estimate magnitude limit, close files, and return. C 9000 CONTINUE write(*,*)'left object: ',kcorner close(2) close(3) IF (SUMWT .LE. 0.0) RETURN MAGLIM=MAGLIM/SUMWT MAGSQ=MAGSQ/SUMWT-MAGLIM**2 MAGSQ=DSQRT(DMAX1(0.D0, MAGSQ)) WRITE (LINE,630) MAGLIM, MAGSQ 630 FORMAT (' Estimated magnitude limit (Aperture 1): ', F4.1, . ' +-', F4.1, ' per star.') RETURN C C----------------------------------------------------------------------- C END! C C####################################################################### C SUBROUTINE QUICK (DATUM, N, INDEX) C C======================================================================= C C A quick-sorting algorithm suggested by the discussion on pages 114-119 C of THE ART OF COMPUTER PROGRAMMING, Vol. 3, SORTING AND SEARCHING, by C D.E. Knuth, which was referenced in Don Wells' subroutine QUIK. This C is my own attempt at encoding a quicksort-- PBS. C C Arguments C C DATUM (INPUT/OUTPUT) is a vector of dimension N containing randomly C ordered real data upon input. Upon output the elements of C DATUM will be in order of increasing value. C C C INDEX (OUTPUT) is an integer vector of dimension N. Upon return to C the calling program the i-th element of INDEX will tell where C the i-th element of the sorted vector DATUM had been BEFORE C DATUM was sorted. C C======================================================================= C IMPLICIT NONE INTEGER MAXSTK, N PARAMETER (MAXSTK=28) C C Parameter C C MAXSTK is the maximum number of entries the stack can contain. C A limiting stack length of 14 restricts this quicksort C subroutine to vectors of maximum length of order 32,768 C (= 2**15). REAL DATUM(N) INTEGER INDEX(N), STKLO(MAXSTK), STKHI(MAXSTK) C REAL DKEY INTEGER I, HI, LO, NSTAK, LIMLO, LIMHI, IKEY C C Initialize INDEX. C DO I=1,N INDEX(I)=I END DO C C Initialize the pointers. C NSTAK=0 LIMLO=1 LIMHI=N C 100 DKEY=DATUM(LIMLO) IKEY=INDEX(LIMLO) C TYPE *, 'LO =', LIMLO, ' HI =', LIMHI C C Compare all elements in the sub-vector between LIMLO and LIMHI with C the current key datum. C LO=LIMLO HI=LIMHI 101 CONTINUE C IF (LO .EQ. HI)GO TO 200 C IF (DATUM(HI) .LE. DKEY) GO TO 109 HI=HI-1 C C The pointer HI is to be left pointing at a datum SMALLER than the C key, which is intended to be overwritten. C GO TO 101 C 109 DATUM(LO)=DATUM(HI) INDEX(LO)=INDEX(HI) LO=LO+1 110 CONTINUE C IF (LO .EQ. HI) GO TO 200 C IF (DATUM(LO) .GE. DKEY) GO TO 119 C LO=LO+1 GO TO 110 C 119 DATUM(HI)=DATUM(LO) INDEX(HI)=INDEX(LO) HI=HI-1 C C The pointer LO is to be left pointing at a datum LARGER than the C key, which is intended to be overwritten. C GO TO 101 C 200 CONTINUE C C LO and HI are equal, and point at a value which is intended to C be overwritten. Since all values below this point are less than C the key and all values above this point are greater than the key, C this is where we stick the key back into the vector. C DATUM(LO)=DKEY INDEX(LO)=IKEY C DO 1666 I=LIMLO,LO-1 C1666 TYPE *, DATUM(I) C TYPE *, DATUM(LO), ' KEY' C DO 2666 I=LO+1,LIMHI C2666 TYPE *, DATUM(I) C C At this point in the subroutine, all data between LIMLO and LO-1, C inclusive, are less than DATUM(LO), and all data between LO+1 and C LIMHI are larger than DATUM(LO). C C If both subarrays contain no more than one element, then take the most C recent interval from the stack (if the stack is empty, we're done). C If the larger of the two subarrays contains more than one element, and C if the shorter subarray contains one or no elements, then forget the C shorter one and reduce the other subarray. If the shorter subarray C contains two or more elements, then place the larger subarray on the C stack and process the subarray. C IF (LIMHI-LO .GT. LO-LIMLO) GO TO 300 C C Case 1: the lower subarray is longer. If it contains one or no C elements then take the most recent interval from the stack and go C back and operate on it. C IF (LO-LIMLO .LE. 1) GO TO 400 C C If the upper (shorter) subinterval contains one or no elements, then C process the lower (longer) one, but if the upper subinterval contains C more than one element, then place the lower (longer) subinterval on C the stack and process the upper one. C IF (LIMHI-LO .GE. 2) GO TO 250 C C Case 1a: the upper (shorter) subinterval contains no or one elements, C so we go back and operate on the lower (longer) subinterval. C LIMHI=LO-1 GO TO 100 C 250 CONTINUE C C Case 1b: the upper (shorter) subinterval contains at least two C elements, so we place the lower (longer) subinterval on the stack and C then go back and operate on the upper subinterval. C NSTAK=NSTAK+1 IF (NSTAK .GT. MAXSTK) THEN c CALL STUPID ('Stack overflow in QUICK. Increase MAXSTK.') c CALL OOPS END IF STKLO(NSTAK)=LIMLO STKHI(NSTAK)=LO-1 LIMLO=LO+1 C DO 3666 I=1,NSTAK C3666 TYPE *, 'STACK: ', I, STKLO(I), STKHI(I) GO TO 100 C 300 CONTINUE C C Case 2: the upper subarray is longer. If it contains one or no C elements then take the most recent interval from the stack and C operate on it. C IF (LIMHI-LO .LE. 1) GO TO 400 C C If the lower (shorter) subinterval contains one or no elements, then C process the upper (longer) one, but if the lower subinterval contains C more than one element, then place the upper (longer) subinterval on C the stack and process the lower one. C IF (LO-LIMLO .GE. 2) GO TO 350 C C Case 2a: the lower (shorter) subinterval contains no or one elements, C so we go back and operate on the upper (longer) subinterval. C LIMLO=LO+1 GO TO 100 C 350 CONTINUE C C Case 2b: the lower (shorter) subinterval contains at least two C elements, so we place the upper (longer) subinterval on the stack and C then go back and operate on the lower subinterval. C NSTAK=NSTAK+1 IF (NSTAK .GT. MAXSTK) THEN c CALL STUPID ('Stack overflow in QUICK. Increase MAXSTK.') c CALL OOPS END IF STKLO(NSTAK)=LO+1 STKHI(NSTAK)=LIMHI LIMHI=LO-1 C DO 4666 I=1,NSTAK C4666 TYPE *, 'STACK: ', I, STKLO(I), STKHI(I) GO TO 100 C 400 CONTINUE C C Take the most recent interval from the stack. If the stack happens C to be empty, we are done. C IF (NSTAK .LE. 0) THEN RETURN ! Normal return END IF C TYPE *, 'POP: ', NSTAK, STKLO(NSTAK), STKHI(NSTAK) LIMLO=STKLO(NSTAK) LIMHI=STKHI(NSTAK) NSTAK=NSTAK-1 GO TO 100 C END! C C C####################################################################### C SUBROUTINE MMM (SKY, NSKY, HIBAD, READNS, SKYMN, SKYMED, . SKYMOD, SIGMA, SKEW) C C======================================================================= C C Official DAO version: 1988 July 1 C C This version of MMM (modified by PBS 1984.IV.10ff) assumes that C the sky brightnesses in the one-dimensional array SKY are already C sorted on entering this routine, and that pixels outside the "bad" C limits have already been eliminated. C C This particular version of MMM also takes cognizance of the fact that, C pixels falling below the LOBAD threshold already having been C eliminated, the contaminated sky pixels values overwhelmingly display C POSITIVE departures from the true value. C C If for some reason it is impossible to obtain the mode of the sky C distribution, this will be flagged by setting SIGMA = -1.0. C C Arguments C C SKY (INPUT) is a real vector containing actual sorted sky values. C NSKY (INPUT) is the number of defined elements in SKY. C SKYMOD (OUTPUT) is the estimated mode of the sky values. C SIGMA (OUTPUT) is the computed standard deviation of the peak in C the sky histogram. C SKEW (OUTPUT) is the computed skewness of the peak in the sky C histogram. C C======================================================================= C IMPLICIT NONE INTEGER NSKY REAL SKY(NSKY) C DOUBLE PRECISION DSQRT, DBLE REAL ALOG10, AMIN1, AMAX1 C DOUBLE PRECISION SUM,SUMSQ REAL CUT, CUT1, CUT2, DELTA, SKYMID, SKYMED, SKYMN, SKYMOD REAL SIGMA, SKEW, R, SIGN, HIBAD, CENTER, SIDE, READNS REAL DMOD, OLD, CLAMP INTEGER I, J, K, L, M INTEGER MINIMM, MAXIMM, NITER, ISTEP, MAXIT, MINSKY, JSTEP LOGICAL REDO DATA MAXIT / 30 /, MINSKY / 20 / C C----------------------------------------------------------------------- C C SECTION 1 C IF (NSKY .LE. 0) THEN GO TO 9900 END IF SKYMID=0.5*(SKY((NSKY+1)/2)+SKY(NSKY/2+1)) C C SKYMID is the median value for the whole ensemble of sky pixels. C Notice that if NSKY is odd, then (NSKY+1)/2 and (NSKY/2)+1 will be the C same number, whereas if NSKY is even, they will be different numbers. C This same trick will be used again later. C C Initialize the variables for accumulating the mean and standard C deviation, and initialize the rejection limits. C SUM=0.D0 SUMSQ=0.D0 CUT1=AMIN1(SKYMID-SKY(1), SKY(NSKY)-SKYMID, HIBAD-SKYMID) C C For the first pass we will consider only pixels in a symmetric C interval of brightness values about the median value. This exploits C the assumption that all the bad pixels are already rejected from the C lower end of the brightness range. C CUT2=SKYMID + CUT1 CUT1=SKYMID - CUT1 C MINIMM=0 DO 1010 I=1,NSKY IF (SKY(I) .LT. CUT1) THEN MINIMM=I GO TO 1010 END IF IF (SKY(I) .GT. CUT2) GO TO 1020 DELTA=SKY(I)-SKYMID SUM=SUM+DELTA SUMSQ=SUMSQ+DELTA**2 MAXIMM=I 1010 CONTINUE C C Henceforth in this subroutine, MINIMM will point to the highest value C rejected at the lower end of the vector, and MAXIMM will point to the C highest value accepted at the upper end of the vector. C MAXIMM-MINIMM is the number of pixels within the acceptance range. C C Compute mean and sigma (from the first pass). C 1020 CONTINUE SKYMED=0.5*(SKY((MINIMM+MAXIMM+1)/2)+SKY((MINIMM+MAXIMM)/2+1)) SKYMN=SUM/DBLE(MAXIMM-MINIMM) SIGMA=DSQRT(SUMSQ/DBLE(MAXIMM-MINIMM)-SKYMN**2) SKYMN=SKYMN+SKYMID C C The middle sky value, SKYMID, was subtracted off up above and added C back in down here to reduce the truncation error in the computation C of SIGMA. C Note that this definition of SIGMA is incorrect by a factor of C SQRT [NSKY/(NSKY-1.)], but for all but pathological cases (where none C of this can be trusted anyway), it's close enough. C SKYMOD=SKYMN IF (SKYMED .LT. SKYMN) SKYMOD=3.*SKYMED-2.*SKYMN C C If the mean is less than the mode, that means the contamination is C slight, and the mean value is what we really want. Note that this C introduces a slight bias toward underestimating the sky when C the scatter in the sky is caused by random fluctuations rather than C by contamination, but I think this bias is negligible compared to the C problem of contamination. C C----------------------------------------------------------------------- C C SECTION 2 C C Rejection and recomputation loop: C NITER=0 OLD = 0. CLAMP = 1. 2000 NITER=NITER+1 IF ((NITER .GT. MAXIT) .OR. (MAXIMM-MINIMM .LT. MINSKY)) THEN GO TO 9900 END IF C C Compute Chauvenet rejection criterion. C R=ALOG10(FLOAT(MAXIMM-MINIMM)) R=AMAX1(2., (-.1042*R+1.1695)*R+.8895) C C Compute rejection limits (symmetric about the current mode). C CUT=R*SIGMA+0.5*ABS(SKYMN-SKYMOD) CUT=AMAX1(1.5,CUT) CUT1=SKYMOD-CUT CUT2=SKYMOD+CUT C C Recompute mean and sigma by adding and/or subtracting sky values C at both ends of the interval of acceptable values. C C At each end of the interval, ISTEP will show the direction we have to C step through the vector to go from the old partition to the new one. C Pixels are added or subtracted depending upon whether the limit is C moving toward or away from the mode. C REDO=.FALSE. C C Is CUT1 above or below the minimum currently-accepted value? C ISTEP=INT(SIGN(1.0001, CUT1-SKY(MINIMM+1))) JSTEP=(ISTEP+1)/2 C C If ISTEP = +1, JSTEP = 1. If ISTEP = -1, JSTEP=0. If ISTEP = +1, C then we know that at least one pixel must be deleted at the low end. C IF (ISTEP .GT. 0) GO TO 2120 2100 IF ((ISTEP .LT. 0) .AND. (MINIMM .LE. 0)) GO TO 2150 C C Quit when SKY(MINIMM) < CUT1 <= SKY(MINIMM+1) C IF ((SKY(MINIMM) .LE. CUT1) .AND. (SKY(MINIMM+1) .GE. CUT1)) . GO TO 2150 C C If ISTEP is positive, subtract out the sky value at MINIMM+1; if C ISTEP is negative, add in the sky value at MINIMM. C 2120 CONTINUE DELTA=SKY(MINIMM+JSTEP)-SKYMID SUM=SUM-REAL(ISTEP)*DELTA SUMSQ=SUMSQ-REAL(ISTEP)*DELTA**2 MINIMM=MINIMM+ISTEP REDO=.TRUE. ! A change has occured GO TO 2100 C 2150 CONTINUE C C Is CUT2 above or below the current maximum? C ISTEP=INT(SIGN(1.0001, CUT2-SKY(MAXIMM))) JSTEP=(ISTEP+1)/2 C C If ISTEP = +1, JSTEP = 1. If ISTEP = -1, JSTEP=0. If ISTEP = -1, C then we know that we must subtract at least one pixel from the high C end. C IF (ISTEP .LT. 0) GO TO 2220 2200 IF ((ISTEP .GT. 0) .AND. (MAXIMM .GE. NSKY)) GO TO 2250 C C Quit when SKY(MAXIMM) <= CUT2 < SKY(MAXIMM+1) C IF ((SKY(MAXIMM) .LE. CUT2) .AND. (SKY(MAXIMM+1) .GE. CUT2)) . GO TO 2250 C C If ISTEP is positive, add in the sky value at MAXIMM+1; if ISTEP is C negative, subtract off the sky value at MAXIMM. C 2220 DELTA=SKY(MAXIMM+JSTEP)-SKYMID SUM=SUM+REAL(ISTEP)*DELTA SUMSQ=SUMSQ+REAL(ISTEP)*DELTA**2 MAXIMM=MAXIMM+ISTEP REDO=.TRUE. ! A change has occured GO TO 2200 C 2250 CONTINUE C C Compute mean and sigma (from this pass). C SKYMN=SUM/DBLE(MAXIMM-MINIMM) SIGMA=DSQRT(SUMSQ/DBLE(MAXIMM-MINIMM)-SKYMN**2) SKYMN=SKYMN+SKYMID C C Obtain the median. To first approximation, the median would be the C value of the sky in the middle pixel in the sorted data (if the C total number is odd) or the mean of the two pixels straddling C the middle (if the total number of pixels is even). C C SKYMED=0.5*(SKY((MINIMM+MAXIMM+1)/2)+SKY((MINIMM+MAXIMM)/2+1)) C C However, this is not good enough. If you look at the estimator for C the mode, you will note that a tiny change in the list of sky pixels, C just sufficient to alter the median value of the sky brightness by C one unit, will change the estimator of the mode by three units. We C really want something more robust than this. As a first attempt C at a more robust median estimator, I estimated the median C of the distribution by the mean of the central twenty percent of sky C values. This involved considerable care to make sure you get C a perfectly symmetric sample of pixels about the median, whether C there is an even or an odd number of pixels within the acceptance C interval. However, even this is not good enough when you have C quantized intensities and detectors with low enough read noise that C a large fraction of the sky pixels have the same (quantized) intensity: C a given fixed percentage on either side of the median does not C necessarily give a good estimator of how the values differ to the two C sides. So I add another criterion: if either limiting value is not C sufficiently different from the median, expand the range by one pixel C on either side until we reach a significantly different intensity value. C In a normal distribution, the central 20% corresponds to +/- 0.25 sigma, C so expand the range until both limiting values differ from the central C value by at least 0.25 times the read noise. C SKYMED=0.0 CENTER = REAL(MINIMM+1 + MAXIMM)/2. SIDE = REAL(NINT(0.2*REAL(MAXIMM-MINIMM)))/2. + 0.25 J = NINT(CENTER-SIDE) K = NINT(CENTER+SIDE) L = NINT(CENTER-0.25) M = NINT(CENTER+0.25) R = 0.25*READNS 2305 CONTINUE IF ((J .GT. 1) .AND. (K .LT. NSKY) .AND. ( . (SKY(L)-SKY(J) .LT. R) .OR. (SKY(K)-SKY(M) .LT. R) )) THEN J = J-1 K = K+1 GO TO 2305 END IF C DO 2310 I=J,K 2310 SKYMED=SKYMED+SKY(I) C SKYMED=SKYMED/REAL(K-J+1) IF (SKYMED .LT. SKYMN) THEN DMOD=3.*SKYMED-2.*SKYMN - SKYMOD ELSE DMOD=SKYMN - SKYMOD END IF C C If the mean is less than the mode, that means the contamination is C slight, and the mean value is what we really want. Note that this C introduces a slight bias toward underestimating the sky when C the scatter in the sky is caused by random fluctuations rather than C by contamination, but I think this bias is negligible compared to the C problem of contamination. C C If the limits have not yet stopped moving, try again. C IF (DMOD*OLD .LT. 0.) CLAMP = 0.5*CLAMP SKYMOD = SKYMOD + CLAMP*DMOD OLD = DMOD IF (REDO) GO TO 2000 C C----------------------------------------------------------------------- C C Normal return. C SKEW=(SKYMN-SKYMOD)/AMAX1(1., SIGMA) NSKY=MAXIMM-MINIMM RETURN C C----------------------------------------------------------------------- C C An error condition has been detected. C 9900 SIGMA=-1.0 SKEW=0.0 RETURN C END! C C####################################################################### C SUBROUTINE WRHEAD (LUN, NL, NCOL, NROW, ITEMS, LOBAD, HIBAD, . THRESH, AP1, PHPADU, READNS, FRAD) C C======================================================================= C C Subroutine to write a standard header into an output sequential C data file on the disk. Same as RDHEAD, except that all the output C arguments are now input arguments. ITEMS tells how many of the C individual arguments are to be written into the header. C C Latest DAOPHOT II version: 1991 January 15 C======================================================================= C IMPLICIT NONE CHARACTER*8 HEAD(7) REAL A(7) C REAL LOBAD, HIBAD, THRESH, AP1, PHPADU, READNS, FRAD INTEGER LUN, NL, NCOL, NROW, ITEMS, I DATA HEAD /' LOWBAD', ' HIGHBAD', ' THRESH', ' AP1', . ' PH/ADU', ' RNOISE', ' FRAD'/ C C----------------------------------------------------------------------- C A(1)=LOBAD A(2)=HIBAD A(3)=THRESH A(4)=AP1 A(5)=PHPADU A(6)=READNS A(7)=FRAD WRITE (LUN,900) (HEAD(I), I=1,ITEMS) 900 FORMAT (' NL NX NY', 8A8) IF (A(2) .LT. 99999.9) THEN WRITE (LUN,901) NL, NCOL, NROW, (A(I), I=1,ITEMS) 901 FORMAT (1X, I2, 2I5, 2F8.1, 5F8.2) ELSE WRITE (LUN,902) NL, NCOL, NROW, (A(I), I=1,ITEMS) 902 FORMAT (1x, I2, 2I5, F8.0, F9.0, 5F8.2) END IF WRITE (LUN,901) ! Write a blank line RETURN ! Normal return C END! C C####################################################################### C SUBROUTINE RDSTAR (LU, NL, ID, X, Y, AMAG, SKY) C C======================================================================= C C Read in an ID number, x and y coordinates, magnitude, and sky value C for the next star in the input file, whatever the file type (i.e., C NL = 1, 2, or 3). If an end of file is encountered, set the star ID C to a negative number and return; if a blank line is encountered, the C star ID will be zero. LU is the logical unit number to be used; the C other arguments are obvious. C C======================================================================= C IMPLICIT NONE CHARACTER*133 LINE REAL X, Y, AMAG, SKY, DUMMY INTEGER LU, N, ISTAT, ID, NL C 900 CALL RDCHAR (LU, LINE, N, ISTAT) IF (ISTAT .GT. 0) THEN ID = -1 ! ID = -1 means END OF FILE RETURN ELSE IF (ISTAT .LT. 0) THEN write(*,*)'Unable to read line from input file' END IF C IF (N .LE. 1) THEN ID = 0 ! Blank line RETURN END IF C N = N+1 LINE(N:N) = '/' IF (NL .EQ. 1) THEN READ (LINE(1:N),*,ERR=1000) ID, X, Y, AMAG, DUMMY, SKY ELSE IF (NL .EQ. 2) THEN READ (LINE(1:N),*,ERR=1000,END=2000) ID, X, Y, AMAG IF (ID .NE. 0) READ(LU,321,ERR=1000) SKY 321 FORMAT (4X, F9.3) ELSE IF (NL .EQ. 3) THEN READ (LINE(1:N),*,ERR=1000) ID, X, Y, AMAG, SKY END IF RETURN ! Normal return 1000 write(*,*)'WARNING: Corrupt star data encountered in input:' GO TO 900 2000 ID=-1 RETURN END! C C####################################################################### C SUBROUTINE RDCHAR (LUN, LINE, N, ISTAT) IMPLICIT NONE CHARACTER*(*) LINE INTEGER LUN, N, ISTAT, I, J C READ (LUN,1,END=8000,ERR=9000) LINE 1 FORMAT (A) N = 0 DO I=1,LEN(LINE) J = ICHAR(LINE(I:I)) IF ((J .GE. 33) .AND. (J .LE. 126)) N = I END DO ISTAT = 0 RETURN 8000 LINE = ' ' N = 0 ISTAT = 1 ! ISTAT = 1 means END OF FILE RETURN 9000 LINE = ' ' N = 0 ISTAT = -1 ! ISTAT = -1 means ERROR RETURN END! function indexpos(head,f1) character*80 head(36),f1*8 do 10 indexpos=1,36 10 if(head(indexpos)(1:8).eq.f1)return end