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


