	parameter (n512=512)
	parameter (npos=5)

        PARAMETER (MAXEXP=10, MAXPAR=6, MAXBOX=13, MAXSKY=10000)
        PARAMETER (NOPT=20, NCMD=23, MAXPSF=207)
        parameter (MAXPIC=3*n512*MAXBOX)

	real w1(MAXPIC)
	real w2(MAXBOX*MAXBOX*2) 
        real opt(nopt)
	
	real         a_b(n512*n512)
	real	     jnk(n512*n512)
	character*80 head(36)
	logical logi	
	character*60  f1,f2,f3
	common /SIZE/ ncol,nrow
	common /EDGE/ n12,ix1,ix2,iy1,iy2,seeing
	seeing=0.
	
        if(iargc().lt.1)then
        write(*,*)
        write(*,*)
        write(*,*)'     **********   DAO_FIND  2000.7 ***********'
        write(*,*)
        write(*,*)
	write(*,*)'     Usage:  dao_find file'
	write(*,*)
	write(*,*)
	stop
	endif
	call getarg(1,f1)
        k=index(f1,'.')
        if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit'
        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,'BITPIX  ')
	if(ipos.lt.36)read(head(ipos)(26:),*)n32
	ipos=indexpos(head,'NAXIS   ')
	if(ipos.lt.36)read(head(ipos)(65:77),'(i5,i8)')n12,n49
	if(n12.ne.0)then
	ipos=indexpos(head,'NAXIS1  ')
	read(head(ipos)(65:74),'(2i5)')ix1,ix2
	ipos=indexpos(head,'NAXIS2  ')
	read(head(ipos)(65:74),'(2i5)')iy1,iy2
          b_high=490000.
          if(n49.ne.0)b_high=n49
	else
	ipos=indexpos(head,'TTIME   ')
	if(head(ipos)(25:25).eq.">")then
	  write(*,*)" use cob?? on linux again please"
	  stop
	endif	
          b_high=28000.
	endif	
	ipos=indexpos(head,'VOLT3   ')            ! SEEING
	if(ipos.lt.36)read(head(ipos)(25:),*)seeingo
          seeing=seeingo/1.67
	close(1)

        k=index(f1,'.')
        f2=f1(1:k)//'coo'
        open(3,file=f2,status='UNKNOWN')

	f3='pip2.par'
        inquire(file=f3,exist=logi)
        if(.not.logi)f3='/vega2/qso5bin/pip2.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 2000 3 8 for our ccd only
          if(seeingo.gt.5..and.opt(5).lt.5.)then
	    opt(7)=opt(7)-(seeingo-5.)*0.1
	write(*,*)'larger seeing, change Low sharpness cutoff',opt(7)
	  endif
	close(1)	

	call readfits(f1,a_b,ncol,nrow)
        if(n32.eq.16)call int2real(a_b,a_b,ncol,nrow)
	if(n12.eq.0)then
	n12=1
	ix1=1
	iy1=1
	ix2=ncol
	iy2=nrow
	endif
	ix1=ix1+npos
	iy1=iy1+npos
	ix2=ix2-npos
	iy2=iy2-npos

        MAX = MAXPIC/3
        MAXCOL = MAX/MAXBOX
        CALL FIND (w1,w1(MAX+1),w1(2*MAX+1),w2,w2(MAXBOX*MAXBOX+1),
     *         MAX, MAXBOX, MAXCOL, MAXSKY, OPT, NOPT,a_b,jnk,NSTAR)
 	write(*,*)f1,NSTAR," stars"
	end	

	subroutine int2real(a,b,n1,n2)
	real a(1)
	integer*2 b(1),c(512*512)
	k=n1*n2
	do 10 i=1,k
10	c(i)=b(i)
	do 20 i=1,k
20	a(i)=c(i)
	end

      SUBROUTINE  FIND (D, H, JCYLN, G, SKIP, MAX, MAXBOX, MAXCOL, 
     .     MAXSKY, OPT, NOPT, a_b,jnk,NSTAR)
c      IMPLICIT NONE
      INTEGER MAX, MAXBOX, MAXCOL, MAXSKY, NOPT, n12,ix1,ix2,iy1,iy2
C
      REAL D(MAXCOL,MAXBOX), H(MAXCOL,MAXBOX), DATA(2), OPT(NOPT)
      REAL G(MAXBOX,MAXBOX), AMAX1
      INTEGER JCYLN(MAX) 
      LOGICAL SKIP(MAXBOX,MAXBOX)
C
C Variables
C
      CHARACTER LINE*5
      REAL PIXELS, RADIUS, FWHM, SIGSQ, RSQ, RELERR, SKYLVL, TEMP
      REAL HMIN, LOBAD, HIBAD, WATCH, P, DATUM, HEIGHT, DENOM, SGOP
      REAL SHARP, ROUND, SHRPLO, SHRPHI, RNDLO, RNDHI
      REAL SUMG, SUMGSQ, SUMGD, SUMD, SG, SGSQ, SGD, SD, WT, HX, HY
      REAL DGDX, SDGDX, SDGDXS, SDDGDX, SGDGDX, SUM2, SUM4
      REAL XCEN, YCEN, DX, DY, PHPADU, READNS, SKYMOD
      INTEGER NHALF, NBOX, MIDDLE, LASTCL, LASTRO, NCOL, NROW, JSQ
      INTEGER ISTAT, NROWS, NSTAR
      INTEGER I, J, K, N, IX, IY, JX, JY, KX, LX, LY
      LOGICAL SATUR8
C
C Common
C
	real a_b(512*512) 
	real jnk(512*512)
        common /SIZE/ ncol,nrow
	common /EDGE/ n12,ix1,ix2,iy1,iy2,seeing

      HIBAD = OPT(4)
      FWHM  = OPT(5)
      if(seeing.ne.0.)then
	FWHM=seeing
      endif
      SHRPLO= OPT(7)
      SHRPHI= OPT(8)
      RNDLO = OPT(9)
      RNDHI = OPT(10)
      WATCH = OPT(11)
C
      RADIUS=AMAX1(2.001, 0.637*FWHM)
      NHALF=MIN0((MAXBOX-1)/2, INT(RADIUS))
      NBOX=2*NHALF+1                ! Length of the side of the subarray
      MIDDLE=NHALF+1
C
      LASTRO=NROW-NHALF
      LASTCL=NCOL-NHALF
C
      SIGSQ=(FWHM/2.35482)**2
      RADIUS=RADIUS**2
C
      SUMG=0.0
      SUMGSQ=0.0
      PIXELS=0.0
      DO J=1,NBOX
         JSQ=(J-MIDDLE)**2
C
         DO I=1,NBOX
            RSQ=FLOAT((I-MIDDLE)**2+JSQ)
            G(I,J)=EXP(-0.5*RSQ/SIGSQ)
            IF (RSQ .LE. RADIUS) THEN
               SKIP(I,J)=.FALSE.
               SUMG=SUMG+G(I,J)
               SUMGSQ=SUMGSQ+G(I,J)**2
               PIXELS=PIXELS+1.0
            ELSE
               SKIP(I,J)=.TRUE.
            END IF
         END DO
      END DO
      DENOM=SUMGSQ-(SUMG**2)/PIXELS
      SGOP=SUMG/PIXELS
C
      RELERR=1.0/DENOM
      RELERR=SQRT(RELERR)
      CALL SKY (D, H, JCYLN, MIN0(MAX,MAXSKY), OPT(1), HIBAD, READNS,
     .     PHPADU, SKYMOD, IX, a_b)
C
c      write(*,'(a,$)')' Number of frames averaged, summed: '
c      read(*,*)data(1),data(2)
	data(1)=1
	data(2)=n12
      IF ((DATA(1) .LT. 0.5) .OR. (DATA(2) .LT. 0.5)) RETURN
      READNS = OPT(1)**2*DATA(2)/DATA(1)
      PHPADU = OPT(2)*DATA(1)
      HMIN = SQRT(READNS + AMAX1(0.,SKYMOD)/PHPADU)
      LOBAD = 0.1*NINT(10.*(SKYMOD-OPT(3)*HMIN))
      HMIN = 0.01*NINT(100.*OPT(6)*RELERR*HMIN)
      READNS = SQRT(READNS)
C
C
      DO JY=1,MIDDLE-1
         DO IX=1,NCOL
            D(IX,JY) = -1.1E38
         END DO
      END DO
C
C Now read the first NHALF rows into the cylinder buffer, putting
C row 1 of the image into row MIDDLE of the buffer.
C
      LX = 1
      LY = 1
      NROWS = NHALF
      CALL RDARAY (a_b,LX,LY,NCOL,NROWS,MAXCOL,D(1,MIDDLE),ISTAT)
      NROWS=1
C
C The cylinder buffer D now contains the actual image data for the
C first NHALF rows of the picture.  We will soon create the file
C containing the derived values of  h  (see above) one row at a time.
C
C Now we will step through the picture row by row. JY remembers which
C row in the big picture we are working on.  For each row JY, the 
C convolved data will be accumulated in the vector H(i,2), and then 
C written into the JY-th row of the scratch picture.
C
      JY=0
 2020 JY=JY+1                              ! Increment image-row pointer
      IF (JY .GT. NROW) GO TO 2100         ! Have we reached the bottom?
C
C Note that at any given time we have only NBOX rows of the original 
C image in memory, contained in the cylinder buffer D(i,j), 
C j = 1, ..., NBOX, but not necessarily in that order.  For instance, 
C if NBOX = 5, when JY = 1,
C 
C      row:     1    2    3    4    5     of G is to be fitted to
C
C      row:     *    *    1    2    3     of the original picture which
C                                         is contained in
C      row:     1    2    3    4    5     of D.
C 
C When row 1 of the picture is done, JY is set to 2, and row 4 
C of the original picture is read into row 1 of the cylinder buffer,
C D, overwriting the null values which we put there before.
C Hence, when JY = 2,
C 
C      row:     1    2    3    4    5     of G is to be fitted to
C
C      row:     *    1    2    3    4     of the original picture which
C                                         is contained in
C      row:     2    3    4    5    1     of D.
C
C As a final example, consider the situation for JY = 7:
C
C      row:     1    2    3    4    5     of G is to be fitted to
C
C      row:     5    6    7    8    9     of the original picture which
C                                         is contained in
C      row:     2    3    4    5    1     of D.
C
C In other words:
C
C      row:     1    2    3    4    5     variable J
C
C      row:     5    6    7    8    9     variable JY
C                                         
C      row:     2    3    4    5    1     vector JCYLN(J)
C
C The cylinder buffer, D, just rolls down through the picture like a
C caterpillar tread, dropping off rows of data when they are no longer
C necessary and picking up new ones in their place.  The data are
C handled in this way (a) to minimize the amount of memory required,
C by storing only those rows that are immediately wanted, consistent
C with (b) minimizing the number of data transfers.  Now, for the 
C CURRENT value of JY, which row of the cylinder buffer is to be fitted 
C to each row of G?  The answers will be contained in the vector 
C JCYLN.
C
C JCYLN(MIDDLE) is the row in the cylinder buffer where we will find 
C the data for row JY of the big picture, which is to be fitted to row
C MIDDLE of G.  Similarly, JCYLN returns the position in the cylinder 
C buffer of the row to be fitted to the J-th row of G 
C (J = 1, ..., NBOX). 
C
C Now that this is all straight, read in the data for row JY+NHALF 
C (overwriting the data for row JY-NHALF-1, which is no longer needed).
C
      DO J=1,NBOX
         IY = JY + (J - MIDDLE)
C
C IY is that row of the big picture which is to be matched up against
C row J of the Gaussian function, during the convolution of this
C row JY of the big picture.
C
C Which row of the cylinder buffer contains row IY of the big picture?
C
         I = IY + NHALF
C
C I now represents the position that row IY of the big picture would
C have had in the cylinder buffer if the cylinder buffer were
C arbitrarily long, i.e. row 1 of the image in row 3 of D, row 2
C in row 4, row 3 in row 5, row 4 in row 6, ... in the examples
C above.  Now we wrap this around.
C
         JCYLN(J) = MOD(I-1,NBOX) + 1
      END DO
C
      LY = JY+NHALF
      IF (LY .LE. NROW) THEN
         CALL RDARAY (a_b, LX, LY, NCOL, NROWS, MAXCOL, 
     .        D(1,JCYLN(NBOX)), ISTAT)
      ELSE
         K = JCYLN(NBOX)
         DO IX=1,NCOL
            D(IX,K) = -1.1E38
         END DO
      END IF
C
C Compute the local brightness enhancement for each pixel in the row,
C The enhancement is computed from a circular region contained
C within an NBOX x NBOX array centered on the current pixel, using the 
C array, G(I,J), and the constants SUMG, SUMGSQ, and PIXELS computed 
C above.  (These constants will need to be modified if the circular 
C region used in the calculation contains any bad pixels; we will use 
C the variables SG, SGSQ, and P for temporary storage of these 
C constants, and SGD and SD for the accumulation of [G*D] and [D] which
C are also needed.)
C
      DO 2050 JX=1,NCOL
C
      SGD=0.
      SD=0.
      SGSQ=SUMGSQ
      SG=SUMG
      P=PIXELS
C
      DO 2040 IX = JX-NHALF,JX+NHALF
        I = MIDDLE + (IX-JX)
        DO 2040 J=1,NBOX
          K = JCYLN(J)
          IF (SKIP(I,J)) GO TO 2040
            IF ((IX .GE. 1) .AND. (IX .LE. NCOL)) THEN
              DATUM = D(IX,K)
              IF ((DATUM .GE. LOBAD).AND.(DATUM .LE. HIBAD)) THEN
                SGD = SGD+G(I,J)*DATUM
                SD = SD+DATUM
                GO TO 2040
              END IF
            END IF
            SGSQ = SGSQ-G(I,J)**2
            SG = SG-G(I,J)
            P = P-1.
 2040 CONTINUE
C
C Compute the central height of the best fitting Gaussian function,
C temporarily storing it in the variable, then putting it into array
C element H(JX, 2).
C
      IF (P .GT. 1.5) THEN
         IF (P .LT. PIXELS) THEN
            SGSQ = SGSQ-(SG**2)/P
            IF (SGSQ .NE. 0.) THEN
               SGD = (SGD-SG*SD/P)/SGSQ
            ELSE
               SGD = 0.
            END IF
         ELSE
            SGD = (SGD-SGOP*SD)/DENOM
         END IF
      ELSE
         SGD = 0.
      END IF
      H(JX,2) = SGD
C
 2050 CONTINUE
C
C Write this newly-computed row of brightness enhancements to the
C scratch output picture.
C
      IF ((WATCH .GT. 0.5) .AND. (MOD(JY,20) .EQ. 0)) THEN
         WRITE (LINE,620) JY
  620    FORMAT(I5)
      END IF
      CALL WRARAY (jnk, LX, JY, NCOL, NROWS, MAXCOL, H(1,2), ISTAT)
      GO TO 2020
C
 2100 CONTINUE
C
C Later on, when we try to decide whether a local maximum represents
C a stellar profile or a delta function ( = bright bad pixel), we will 
C compare the brightness of the central pixel to the average of the 
C surrounding pixels.  To be ready for that, we here modify SKIP to 
C skip over the central pixel, and set PIXELS equal to the number of 
C pixels in the circular area not counting the central pixel.
C
      SKIP(MIDDLE,MIDDLE) = .TRUE.
c     PIXELS = PIXELS-1.0
C
C-----------------------------------------------------------------------
C
C SECTION 3
C
C Read in both the convolved data from the scratch disk file and the raw
C data from the original picture.  Search for local maxima in the 
C convolved brightness data.  When these are found, compute image-shape
C statistics from the raw data to eliminate non-stellar brightness
C enhancements (as well as possible) and estimate the position of the 
C centroid of the brightness enhancement.
C
 3000 CONTINUE
C
C Now the star search may begin.  The original image data will be read
C into the cylinder buffer D again, just as before.  At the same time, 
C the brightness enhancements will be read from the scratch disk file 
C into another cylinder buffer, H.  The brightness enhancements will 
C then be searched for local maxima.  When these are found, functions 
C of the original image data will be used to derive shape parameters 
C designed to identify bad pixels and bad columns or rows.
C
c----------------------------------------------------------
      CALL WRHEAD (3, 1, NCOL, NROW, 7, LOBAD, HIBAD, HMIN, 
     .     0., PHPADU, READNS, OPT(12))
      IF (WATCH .GT. 0.5) THEN
         WRITE (6,630)
  630    FORMAT(6X, '                              MAGS'/
     .     6X, '                              FROM '/
     .     6X, '       STAR     X      Y     LIMIT    SHARP    ROUND')
      END IF
C
C Read the first NHALF rows into each of the cylinder buffers,
C putting row 1 of the image into row MIDDLE of the buffer.
C
      DO JY=1,MIDDLE-1
         DO IX=1,NCOL
            D(IX,JY) = -1.1E38
            H(IX,JY) = 0.
         END DO
      END DO
C
      LX = 1
      LY = 1
      NROWS = NHALF
      CALL RDARAY (a_b, LX, LY, NCOL, NROWS, MAXCOL, 
     .     D(1,MIDDLE), ISTAT)
      CALL RDARAY (jnk, LX, LY, NCOL, NROWS, MAXCOL, 
     .     H(1,MIDDLE), ISTAT)
      NROWS = 1
C
C Now step through the picture row by row.  Again JY is the image-row 
C counter.
C
      NSTAR = 0
      JY = 0
 3020 JY = JY+1
C
C Have we just finished reducing the last row? If not, work on the 
C next row.  If so, go on to Section 4.
C
      IF (JY .GT. NROW) GO TO 4000
C
C Determine the position in the cylinder buffers of all the rows 
C contained in the box.
C
      DO J=1,NBOX
         IY = JY+(J-MIDDLE)+NHALF
         JCYLN(J) = MOD(IY-1,NBOX) + 1
      END DO
C
C Read in the data for row JY+NHALF.
C
      LY = JY+NHALF
      IF (LY .LE. NROW) THEN
         CALL RDARAY (a_b, LX, LY, NCOL, NROWS, MAXCOL, 
     .        D(1,JCYLN(NBOX)), ISTAT)
         CALL RDARAY (jnk, LX, LY, NCOL, NROWS, MAXCOL, 
     .        H(1,JCYLN(NBOX)), ISTAT)
      ELSE
         K = JCYLN(NBOX)
         DO IX=1,NCOL
            D(IX,K) = -1.1E38
            H(IX,K) = 0.0
         END DO
      END IF
C
C Now step across the row, pixel by pixel.
C
      JX = 1
 3040 HEIGHT = H(JX,JCYLN(MIDDLE))
C
C Sieve to locate a local maximum in the brightness enhancement.  To 
C be a local maximum, the brightness enhancement in a given pixel must
C be above the threshold, and it must also be greater than the 
C brightness enhancement of any pixel within a radius equal to 
C 1.5 sigma.
C
      IF (HEIGHT .LT. HMIN) GO TO 3300
      DO 3051 IX=JX-NHALF,JX+NHALF
         IF ((IX .LT. 1) .OR. (IX .GT. NCOL)) GO TO 3051
         I = MIDDLE + (IX-JX)
         DO 3050 J=1,NBOX
            K = JCYLN(J)
            IF (SKIP(I,J)) GO TO 3050
            IF (HEIGHT .LT. H(IX,K)) GO TO 3300
 3050    CONTINUE
 3051 CONTINUE
C
C The brightness enhancement of this pixel is now confirmed to be above
C the threshold, and to be larger than in any other pixel within a 
C radius of 1.5 sigma.
C
C Now we derive the shape indices.  First, is the object much more
C sharply peaked than the input FWHM?  Compare the central pixel to
C the mean of the surrounding (non-bad) pixels.  If this difference is
C greater than the originally estimated height of the Gaussian or less
C than two-tenths the height of the Gaussian, reject the star
C (assuming SHRPLO and SHRPHI have the default values of 0.2 and
C 1.0; otherwise, muta mutandis.)
C
      SHARP=0.
      ROUND=0.
      DATUM=D(JX,JCYLN(MIDDLE))
      IF ((DATUM .LT. LOBAD) .OR. (DATUM .GT. HIBAD)) GO TO 3068
      P=0.
      DO 3061 IX=JX-NHALF,JX+NHALF
      IF ((IX .LT. 1) .OR. (IX .GT. NCOL)) GO TO 3061
         I = MIDDLE + (JX-IX)
         TEMP = 0.0
         DO 3060 J=1,NBOX
            K=JCYLN(J)
            IF ((IX.EQ.JX).AND.(K.EQ.MIDDLE)) GO TO 3060
            DATUM=D(IX,K)
            IF ((DATUM .GE. LOBAD) .AND. (DATUM .LE. HIBAD)) THEN
               TEMP = TEMP+(DATUM-SKYMOD)
               P = P + 1.
            END IF
 3060    CONTINUE
         SHARP = SHARP+TEMP
 3061 CONTINUE
C
      SHARP=(D(JX,JCYLN(MIDDLE))-SKYMOD-SHARP/P)/HEIGHT
      IF ((SHARP .LT. SHRPLO) .OR. (SHARP .GT. SHRPHI)) GO TO 3200
 3068 CONTINUE
C
C Now check to see whether the object is strongly elongated in some
C arbitrary direction.  Compare a measure of the bilateral symmetry
C with a measure of the fourfold symmetry.
C
C SKIP ALL THIS IF THE STAR IS TOO NEAR THE EDGE OF THE FRAME!
C
      IF ((JX .LT. MIDDLE) .OR. (JX .GT. LASTCL) .OR.
     .     (JY .LT. MIDDLE) .OR. (JY .GT. LASTRO)) THEN
         XCEN = REAL(JX)
         YCEN = REAL(JY)
         GO TO 3190
      END IF
C start
      SUM2 = 0.
      SUM4 = 0.
      DO I=0,NHALF
         DO J=1,NHALF
c           if (jy .ge. 77) type 6667, h(jx-i,jcyln(middle-j)), 
c    .           h(jx+i,jcyln(middle+j)),
c    .           h(jx-j,jcyln(middle+i)), h(jx+j,jcyln(middle-i))
c6667       format (4f6.0)
            SUM2 = SUM2 + 
     .           H(JX-I,JCYLN(MIDDLE-J)) + H(JX+I,JCYLN(MIDDLE+J)) -
     .           H(JX-J,JCYLN(MIDDLE+I)) - H(JX+J,JCYLN(MIDDLE-I))
            SUM4 = SUM4 + 
     .           ABS(H(JX-I,JCYLN(MIDDLE-J))) + 
     .           ABS(H(JX+I,JCYLN(MIDDLE+J))) +
     .           ABS(H(JX-J,JCYLN(MIDDLE+I))) + 
     .           ABS(H(JX+J,JCYLN(MIDDLE-I)))
         END DO
      END DO
      ROUND = 2.*SUM2/SUM4
c     if (jy .ge. 77) type *, sum2, sum4, 4.*sum2/sum4
c     if (jy .ge. 77) accept *
      IF ((ROUND .LT. RNDLO) .OR. (ROUND .GT. RNDHI)) GO TO 3200
      IX = JX-MIDDLE
C
      SUMGD=0.0
      SUMGSQ=0.0
      SUMG=0.0
      SUMD=0.0
      SDGDX=0.0
      SDGDXS=0.0
      SDDGDX=0.0
      SGDGDX=0.0
      P=0.
      N=0
      DO 3073 I=1,NBOX
         SG=0.
         SD=0.
         KX = IX+I
         DO 3070 J=1,NBOX
            WT=FLOAT(MIDDLE-ABS(J-MIDDLE))
            K=JCYLN(J)
            DATUM=D(KX,K)
            IF ((DATUM .GE. LOBAD) .AND. (DATUM .LE. HIBAD)) THEN
               SD=SD+(DATUM-SKYMOD)*WT
               SG=SG+G(I,J)*WT
            END IF
 3070    CONTINUE
         IF (SG .GT. 0.0) THEN
            WT=FLOAT(MIDDLE-ABS(I-MIDDLE))
            SUMGD=SUMGD+WT*SG*SD
            SUMGSQ=SUMGSQ+WT*SG**2
            SUMG=SUMG+WT*SG
            SUMD=SUMD+WT*SD
            P=P+WT
            N=N+1
            DGDX=SG*(MIDDLE-I)
            SDGDXS=SDGDXS+WT*DGDX**2
            SDGDX=SDGDX+WT*DGDX
            SDDGDX=SDDGDX+WT*SD*DGDX
            SGDGDX=SGDGDX+WT*SG*DGDX
         END IF
 3073 CONTINUE
C
C We need at least three points to estimate the height and position
C of the star, and the local sky brightness.
C
      IF (N .LE. 2) GO TO 3200
      HX=(SUMGD-SUMG*SUMD/P)/(SUMGSQ-(SUMG**2)/P)
C
C DX is the height of the best-fitting marginal Gaussian.  If this is
C non-positive, this is not an acceptable star.
C
      IF (HX .LE. 0.) GO TO 3200
C
C Compute the first-order correction to the x-centroid of the star.
C Note that a factor of HX/SIGSQ is missing from SDGDX, SDDGDX, and 
C SGDGDX, and a factor of (HX/SIGSQ)**2 is missing from SDGDXS.
C
      SKYLVL=(SUMD-HX*SUMG)/P
      DX=(SGDGDX-(SDDGDX-SDGDX*(HX*SUMG+SKYLVL*P)))/(HX*SDGDXS/SIGSQ)
C
ccccc      XCEN=JX+DX/(1.+ABS(DX))               90427
       XCEN=JX+AMAX1(-1., AMIN1(1., DX))

C If the best estimate of the star's center falls outside the image,
C reject it.
C
      IF ((XCEN .LT. 0.5) .OR. (XCEN .GT. NCOL+0.5)) GO TO 3200
C
C Compute the height of the y-marginal Gaussian distribution.
C
      SUMGD=0.
      SUMGSQ=0.
      SUMG=0.
      SUMD=0.
      SDGDX=0.
      SDGDXS=0.
      SDDGDX=0.
      SGDGDX=0.
      P=0.
      N=0
      SATUR8 = .FALSE.
      DO 3078 J=1,NBOX
         K=JCYLN(J)
         SG=0.
         SD=0.
         DO 3076 I=1,NBOX
            WT=FLOAT(MIDDLE-ABS(I-MIDDLE))
            KX=IX+I
            DATUM=D(KX,K)
            IF (DATUM .LE. HIBAD) THEN
               IF (DATUM .GE. LOBAD) THEN
                  SD=SD+(DATUM-SKYMOD)*WT
                  SG=SG+G(I,J)*WT
               END IF
            ELSE
               IF (.NOT. SKIP(I,J)) SATUR8 = .TRUE.
            END IF
 3076    CONTINUE
C
         IF (SG .GT. 0.0) THEN
            WT=FLOAT(MIDDLE-ABS(J-MIDDLE))
            SUMGD=SUMGD+WT*SG*SD
            SUMGSQ=SUMGSQ+WT*SG**2
            SUMG=SUMG+WT*SG
            SUMD=SUMD+WT*SD
            P=P+WT
            DGDX=SG*(MIDDLE-J)
            SDGDX=SDGDX+WT*DGDX
            SDGDXS=SDGDXS+WT*DGDX**2
            SDDGDX=SDDGDX+WT*SD*DGDX
            SGDGDX=SGDGDX+WT*SG*DGDX
            N=N+1
         END IF
C
 3078 CONTINUE
C
      IF (N .LE. 2) GO TO 3200
      HY=(SUMGD-SUMG*SUMD/P)/(SUMGSQ-(SUMG**2)/P)
      IF (HY .LE. 0.0) GO TO 3200
      SKYLVL=(SUMD-HY*SUMG)/P
      DY=(SGDGDX-(SDDGDX-SDGDX*(HY*SUMG+SKYLVL*P)))/(HY*SDGDXS/SIGSQ)

cccc      YCEN=JY+DY/(1.+ABS(DY))                 90427
       YCEN=JY+AMAX1(-1., AMIN1(1., DY))

      IF ((YCEN .LT. 0.5) .OR. (YCEN .GT. NROW+0.5)) GO TO 3200
      DY = 2.*(HX-HY)/(HX+HY)
C
C Do not apply this roundness criterion if an off-central pixel is saturated.
C
      IF (.NOT. SATUR8) THEN
         IF ((DY .LT. RNDLO) .OR. (DY .GT. RNDHI)) GO TO 3200
      END IF
C
C
C The fully verified and located star may now be dignified with its own
C ID number.
C
 3190 	if(xcen.lt.ix1)goto 3200
      	if(xcen.gt.ix2)goto 3200
      	if(ycen.lt.iy1)goto 3200
      	if(ycen.gt.iy2)goto 3200
     
      NSTAR=NSTAR+1
      HEIGHT=-2.5*ALOG10(HEIGHT/HMIN)
      IF (WATCH .GT. 0.5) WRITE (6,631) NSTAR, XCEN, YCEN, HEIGHT, 
     .     SHARP, ROUND
  631 FORMAT(12X, I5, 2F7.1, F9.1, 2F9.2)
	if(hx.ge.10000. .or. hy.ge.10000.)then
      WRITE (3,331)NSTAR,XCEN,YCEN,HEIGHT,SHARP,ROUND,DY,hx,hy
	else
      WRITE (3,330)NSTAR,XCEN,YCEN,HEIGHT,SHARP,ROUND,DY,hx,hy
	endif
  330 FORMAT(I6, 8F9.3)
  331 FORMAT(I6, 6F9.3, 2f9.0)

 3200 CONTINUE
C
C If the sieve above (between statements 3040 and 3050) has detected a
C local maximum in the brightness enhancement, whether this enhancement
C was subsequently confirmed to be a star or not, then there is no need
C to check the other pixels in this row between JX+1 and JX+NHALF,
C inclusive, since we know there can't be a local maximum there.
C
      JX = JX+NHALF
 3300 JX = JX+1
C
C Have we passed the last pixel in the row?  If not, work on this 
C pixel.  If so, go to next row.
C
      IF (JX .LE. NCOL) GO TO 3040
      GO TO 3020
C
C-----------------------------------------------------------------------
C
C SECTION 4
C
C Find out whether the user is happy.  If so, delete the scratch picture
C and close up shop.  If not, return to the beginning of Section 3.
C
4000	return
C Normal return.
C
      END!
C
C#######################################################################
C
      SUBROUTINE  SKY  (D, S, INDEX, MAX, READNS, HIBAD, SKYMN, SKYMED,
     .     SKYMOD, N, a_b)
C
C=======================================================================
C
C This subroutine estimates an average sky value for a picture by taking
C individual pixels scattered over the picture.  The brightness values 
C are sorted, and the modal value is estimated using the MMM subroutine.
C
C               OFFICIAL DAO VERSION:  1991 April 18
C
C=======================================================================
C
      IMPLICIT NONE
      INTEGER MAX
C
C MAX    is the maximum number of sky pixels we can deal with,
C        given the limited amount of working space.
C
      REAL S(MAX), D(MAX)
      INTEGER INDEX(MAX)
C
      REAL READNS, HIBAD, SKYMN, SKYMED, SKYMOD, SKYSIG, SKYSKW
      INTEGER NCOL, NROW, ISTEP, LX, LY, NX, NY, IROW, I, N
      INTEGER ISTAT, IFIRST
	real	     a_b(512*512)
	common /SIZE/ ncol,nrow
C
C-----------------------------------------------------------------------
C
C The spacing between pixels that will be included in the sample is
C estimated by the ratio of the total number of pixels in the picture to
C the maximum number of pixels that can be accomodated in the vector S.
C
      ISTEP = NCOL*NROW/MAX+1
C
C Go through the disk file reading a row at a time and extracting every 
C ISTEP-th pixel.  If ISTEP is not equal to 1, make sure that the
C starting pixel for each row is staggered.
C

      LX = 1
      NX = NCOL
      NY = 1
      IFIRST = 0
      N = 0
      DO IROW=1,NROW
         LY = IROW
         CALL RDARAY (a_b, LX, LY, NX, NY, MAX, D, ISTAT)
         IF (ISTAT .NE. 0) RETURN
         IFIRST = IFIRST + 1
         IF (IFIRST .GT. ISTEP) IFIRST = IFIRST - ISTEP
         I = IFIRST
 1010    IF (ABS(D(I)) .LE. HIBAD) THEN
            N = N+1
            S(N) = D(I)
            IF (N .EQ. MAX) GO TO 1100
            I = I + ISTEP
         ELSE
            I = I+1
         END IF
         IF (I .LE. NCOL) GO TO 1010
      END DO
C
C Sort these values, then estimate the mode.
C
 1100 CONTINUE
      CALL QUICK (S, N, INDEX)
      CALL MMM (S, N, HIBAD, READNS, SKYMN, SKYMED, SKYMOD, SKYSIG, 
     .     SKYSKW)
C
C Normal return.
C
      RETURN
      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.1, F8.0, 5F8.2)
      END IF
      WRITE (LUN,901)                             ! Write a blank line
      RETURN                                        ! Normal return
C
      END!
c
        subroutine rdaray (id, lx, ly, mx, my, nx, func, ier)
        real func(nx,*),id(512*512)
        common /SIZE/ ncol,nrow
        mx = lx+mx-1
        my = ly+my-1
        lx = max(1,lx)
        ly = max(1,ly)
        mx = min(ncol,mx)
        my = min(nrow,my)
        my = my-ly+1
        do 10 j=1,my
        jy = ly+j-1
10      call imgs2r (id, func(1,j), lx, mx, jy, jy, ier, ncol)
        mx = mx-lx+1
        end

        subroutine wraray (id, lx, ly, mx, my, maxx, func, ier)
	real id(512*512)
        real func(maxx,*)
        common /size/ ncol, nrow
        mx = lx+mx-1
        my = ly+my-1
        lx = max(1,lx)
        ly = max(1,ly)
        mx = min(ncol, mx)
        my = min(nrow, my)
        nx = mx-lx+1
        ny = my-ly+1
        do 10 j=1,ny
        jy = ly+j-1
10      call imps2r (id, func(1,j), lx, mx, jy, jy, ier, ncol)
        mx = nx
        my = ny
        end
	
	subroutine imps2r(id,buf,i1,i2,j1,j2,ier,ncol)
	real id(ncol,*),buf(1)
 	k=0
	do 10 j=j1,j2
	do 10 i=i1,i2
	k=k+1
10	id(i,j)=buf(k)
	ier=0
	end

	subroutine imgs2r(id,buf,i1,i2,j1,j2,ier,ncol)
	real id(ncol,*),buf(1)
	k=0
	do 10 j=j1,j2
	do 10 i=i1,i2
	k=k+1
10	buf(k)=id(i,j)
	ier=0
	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
 

