C psf 96,12
        PARAMETER (MAXEXP=10, MAXPAR=6)
        PARAMETER (NOPT=20, MAXPSF=207)

        real opt(nopt)
	real psf(MAXPSF,MAXPSF,MAXEXP)
	real par(MAXPAR)
	integer ifail,iflag(MAXPSF)
	
	real         a_b(4096*4032)
	
	character*60  f1,f2,f3,ch50
	character*80  head(72)
	logical logi
	
        if(iargc().lt.1)then
        write(*,*)
        write(*,*)
        write(*,*)'     **********   DAO_PSF  2010.3 ***********'
        write(*,*)
        write(*,*)
	write(*,*)'     Usage:  u_psf4k file.fit [-1,0,1,2] [!]'
	write(*,*)
	write(*,*)
        write(*,*)'     if add !, u_psf4k DONOT delete any star'
	stop
	endif
	call getarg(1,f1)

10      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 !'
 
	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)
c	write(*,'(7f10.2)')opt
	close(1)

        open(1,file=f1,status='old',access='direct',recl=2880*2)
        read(1,rec=1)head
        ipos=indexpos(head,'SEEING  ')            ! SEEING
        if(ipos.lt.72)then
	   read(head(ipos)(25:),*)x               ! 2001,10
           if(x.ne.0)opt(5)=x/0.45*2
	   write(*,*)'FWHM(seeing in pixel)=',opt(5)
	endif
	close(1) 

	call readfits(f1,a_b,ncol,nrow)

        k=index(f1,'.')
        f2=f1(1:k)//'lst'
        open(1,file=f2,status='old')               ! 1 -- *.lst
        f2=f1(1:k)//'ap'
        open(2,file=f2,status='old')               ! 2 -- *.ap
        f2=f1(1:k)//'psf'
        open(3,file=f2,status='UNKNOWN')           ! 3 -- *.psf
        f2=f1(1:k)//'nei'
        open(4,file=f2,status='UNKNOWN')           ! 4 -- *.nei

	if(iargc().gt.1)then
	  call getarg(2,f2)
	  read(f2(1:),*)opt(14)
	endif
	write(*,*)'variable PSF: ',opt(14)

	ifail=0
	if(iargc().gt.2)then
	  call getarg(3,f2)
	  if(f2(1:1).eq.'!')ifail=-1
	endif

	do i=1,MAXPSF
	iflag(i)=0
	enddo	

        call  GETPSF (a_b,NCOL,NROW,PAR,PSF,OPT,NOPT,ifail,iflag)
	if(ifail.eq.1)then
	  write(*,*)' *********** delete Bad PSF_star, redo !'
          k=index(f1,'.')
          f2=f1(1:k)//'lst'
          open(1,file=f2,status='old')               ! 1 -- *.lst
          f3=f1(1:k)//'tmp'
          open(2,file=f3,status='UNKNOWN')
	  i=0             
40	  i=i+1
	  read(1,'(a)',end=50)ch50	  
	  if(iflag(i).eq.0) write(2,'(a)')ch50
	  goto 40
50	  close(1)
	  close(2)
	  call unlink(f2)
          call system("mv "//f3//" "//f2) 
	  goto 10
	endif
 
	end	

      SUBROUTINE  GETPSF (PIC,NCOL,NROW,PAR,PSF,OPT,NOPT,ifail,iflag)
      IMPLICIT NONE
C
      INTEGER MAXBOX, MAXPSF, MAXSTR, MAXN, MAXPAR, MAXEXP, NCOL, NROW
      INTEGER NOPT
      PARAMETER  (MAXPAR=6, MAXPSF=207, MAXEXP=10) 
      PARAMETER  (MAXBOX=69, MAXSTR=50000, MAXN=200) 

	integer ifail,iflag(MAXPSF)
C
      CHARACTER LINE*80, LABEL*8
      CHARACTER  FLAG(MAXN)*1, ANSWER*1
      DOUBLE PRECISION CC(2,2), VV(2), AA(2)
      REAL C(MAXEXP,MAXEXP), V(MAXEXP), A(MAXEXP)
      REAL CON(MAXPSF,MAXPSF), TERM(MAXEXP)
      REAL PIC(NCOL,NROW)
      REAL PAR(MAXPAR)
      REAL PSF(MAXPSF,MAXPSF,MAXEXP)
      REAL XCEN(MAXSTR), YCEN(MAXSTR), APMAG(MAXSTR), SKY(MAXSTR)
      REAL HPSF(MAXN), WEIGHT(MAXN), OPT(NOPT)
      REAL HJNK(MAXN), XJNK(MAXN), YJNK(MAXN)
      REAL PROFIL, BICUBC, SQRT
      INTEGER ID(MAXSTR), NTAB(-1:3), NPARAM
      LOGICAL IN(MAXPSF,MAXPSF), EDGE(MAXPSF,MAXPSF), SATR8D(MAXN)
C
      DOUBLE PRECISION PSFMAG, DLOG10, DBLE
      REAL LOBAD, HIBAD, THRESH, AP1, PHPADU, READNS, RONOIS,
     .     FWHM, WATCH, FITRAD, RADSQ, PSFRAD, PSFRSQ, FMAX, 
     .     DFDX, DFDY, SIG, DX, DY, RDX, RDY, XMID, YMID, DP, OLD, W, 
     .     WT, RADIUS, SCALE, SUMSQ, SUMN, DATUM, RSQ, DX2, DY2
      REAL ERRMAG, CHI, SHARP, PERR, PKERR
      INTEGER I, J, K, L, N, LX, LY, MX, MY, NTOT, NSTAR, ISTAR, NL
      INTEGER IEXPAND, IFRAC, IPSTYP, NPSF, NPAR, ISTAT, IRADSQ,
     .     NEXP, MIDDLE, MIDSQ, IDX, JDY, JDYSQ, IX, JY, ITER, NPASS
      INTEGER II, JJ, IR, JRADSQ, NPSFSQ
      LOGICAL STAROK, SATUR8, REDO
      logical logi
C
      COMMON /ERROR/ PHPADU, RONOIS, PERR, PKERR
C
      DATA NTAB / 0, 1, 3, 6, 10 /
C
C-----------------------------------------------------------------------
C
C SECTION 1
C
C Set up the necessary variables, open the necessary files, read in the
C relevant data for all stars.
C
      FWHM = OPT(5)
      WATCH = OPT(11)
      FITRAD = OPT(12)
      RADSQ = FITRAD**2
      PSFRAD = OPT(13)
      PSFRSQ = (PSFRAD+2.)**2
      IEXPAND = NINT(OPT(14))
      IFRAC = NINT(OPT(15))
      NPASS = NINT(OPT(17))
C
      NPSF = 2*(NINT(2.*PSFRAD)+1)+1
      PSFRAD = (REAL(NPSF-1)/2. - 1.)/2.
      IRADSQ = NINT(2.*PSFRAD)**2
      JRADSQ = NINT(1.4*PSFRAD)**2
C
C Ascertain the name of the input (aperture, usually) photometry file, 
C and read in the relevant data for all stars.
C
	read(2,'(a)')i			! *.ap
      read(2,*) NL,I,I,LOBAD,HIBAD,THRESH,AP1,PHPADU,READNS,DX
      RONOIS = READNS**2
C
      I = 0
 1010 I = I+1
      IF (I .GT. MAXSTR) GO TO 1100
 1020 CALL RDSTAR (2, NL, ID(I), XCEN(I), YCEN(I), APMAG(I), SKY(I))
      IF (ID(I) .LT. 0) GO TO 1100             ! End-of-file encountered
      IF (ID(I) .EQ. 0) GO TO 1020             ! Blank line encountered
      GO TO 1010
C
 1100 NTOT = I-1
	close(2)
C
C-----------------------------------------------------------------------
C
C SECTIONS 2
C
C Learn the name of each PSF star and find it in the star list.  Display
C a subarray around it if desired, and check it for invalid pixels.
C
      NSTAR = 0
      N = NPARAM(1, FWHM, LABEL, PAR, MAXPAR)
C
C NSTAR will point at the last PSF star in the stack after the PSF
C stars have been brought to the top.
C
c	read(1,'(a)')i
c.	read(1,'(a)')i
c.	read(1,'(a)')i
 2000 CONTINUE
      read(1,*,end=3000,err=2000) I, DX, DY, RDX, RDY    ! read *.lst
      IF (I .EQ. 0) GO TO 2000
C
      DO 2010 ISTAR=1,NTOT
      IF (ID(ISTAR) .EQ. I) GO TO 2020
 2010 CONTINUE
      WRITE (6,701) I, ' was not found.'
  701 FORMAT (1X, I5, A)
      GO TO 2000
C
C If a given star appears in the .LST file more than once, employ it
C only once.
C
 2020 IF (ISTAR .LE. NSTAR) GO TO 2000
      IF ( ( INT(XCEN(ISTAR)-FITRAD) .LT. 0 ) .OR.
     .     ( INT(XCEN(ISTAR)+FITRAD) .GT. NCOL ) .OR.
     .     ( INT(YCEN(ISTAR)-FITRAD) .LT. 0 ) .OR.
     .     ( INT(YCEN(ISTAR)+FITRAD) .GT. NROW ) ) THEN
         WRITE (6,701) ID(ISTAR), ' is too near edge of frame.'
         GO TO 2000
      END IF
C
C Define the subarray containing the star and search it for "bad" 
C pixels.  This will include a three-pixel border around the PSF box 
C itself, to guarantee valid interpolations to within the PSF box, and
C to allow for a possible one-pixel shift in the star's centroid.
C
      LX = MAX0( 1, INT(XCEN(ISTAR)-PSFRAD)-2 )
      LY = MAX0( 1, INT(YCEN(ISTAR)-PSFRAD)-2 )
      MX = MIN0(NCOL, INT(XCEN(ISTAR)+PSFRAD)+3 )
      MY = MIN0(NROW, INT(YCEN(ISTAR)+PSFRAD)+3 )
C
      FMAX=-32768.
      STAROK=.TRUE.
      SATUR8=.FALSE.
C
 2040 REDO = .FALSE.
      DO 2070 J=LY,MY
         DP = (REAL(J)-YCEN(ISTAR))**2
         DO 2050 I=LX,MX
            DX = REAL(I)-XCEN(ISTAR)
            RSQ = DX**2 + DP
            IF (RSQ .GT. PSFRSQ) THEN
               IF (DX .GE. 0.) GO TO 2070
               GO TO 2050
            END IF
            DATUM = PIC(I,J)
            IF ((DATUM .LT. LOBAD) .OR. (DATUM .GT. HIBAD)) THEN
C
C There is a defective pixel.
C
               IF (RSQ .LE. RADSQ) THEN
C
C It is inside the fitting radius.
C
                 IF (DATUM .LT. LOBAD) THEN
C
C The pixel is bad low, so reject the star.
C
                   WRITE (6,701) ID(ISTAR), 
     .                  ' is not a good star.' 
	write(*,*)datum,lobad
                   GO TO 2000
                 ELSE IF (.NOT. SATUR8) THEN
C
C The pixel is bad high, so presume the star is saturated.
C
                   WRITE (6,701) ID(ISTAR),
     .                  ' is saturated.'
                   IF (WATCH .GT. -0.5) THEN
	             write(*,'(a,$)')'Use it anyway? '
	             read(*,'(a)')ANSWER
                     IF (ANSWER .EQ. 'Y') THEN
                       SATUR8 = .TRUE.
                     ELSE
                       GO TO 2000
                     END IF
                   ELSE
                     SATUR8 = .TRUE.
                   END IF
                 END IF
                 GO TO 2050
               END IF
C
C The bad pixel isn't inside the fitting radius.
C
               IF (STAROK) THEN
                 STAROK = .FALSE.
                 IF (WATCH .GT. -0.5) THEN
                   WRITE (6,628) ID(ISTAR), DATUM, 7, I, J
 628               FORMAT (1X, I5, ' has a bad pixel:',  F9.1, 
     .                  ' at position', A1, 2I5)
                 END IF
               ELSE
                 IF (WATCH .GT. -0.5) 
     .                WRITE (6,628) ID(ISTAR), DATUM, 0, I, J
               END IF
C
C Replace the defective pixel with an average of surrounding
C pixels.
C
               DATUM = 0.
               SUMN = 0.
               DO L=MAX0(1,J-1),MIN0(NROW,J+1)
                  DO K=MAX0(1,I-1),MIN0(NCOL,I+1)
                     IF ((PIC(K,L) .GE. LOBAD) .AND.
     .                   (PIC(K,L) .LE. HIBAD)) THEN
                        DATUM = DATUM + PIC(K,L)
                        SUMN = SUMN+1.
                     END IF
                  END DO
               END DO
C
C If there are fewer than three valid neighboring pixels, do not
C fudge the pixel yet.  Instead, we will go through again (and again
C if necessary) and allow the fudge to eat into the defect from the
C outside.
C
               IF (SUMN .LT. 2.5) THEN
                  REDO = .TRUE.
               ELSE
                  PIC(I,J) = DATUM/SUMN
               END IF
            ELSE
C
C The pixel wasn't bad.
C
               IF (DATUM .GT. FMAX) FMAX=DATUM
            END IF
 2050    CONTINUE
 2070 CONTINUE
      IF (REDO) GO TO 2040
C
      IF (STAROK) THEN
         IF ((.NOT. SATUR8) .AND. (WATCH .LT. -0.5) .AND.
     .        (WATCH .GT. -1.5)) THEN
C           WRITE (6,701) ID(ISTAR), ' seems fine.'
         END IF
      ELSE IF (.NOT. SATUR8) THEN
         IF ((WATCH .LT. -0.5) .AND. (WATCH .GT. -1.5)) THEN
            WRITE (6,701) ID(ISTAR), ' has bad pixels.'
         ELSE
            IF (WATCH .GT. -0.5) THEN
	       write(*,'(a,$)')'TRy this on anyway? '
	       read(*,'(a)')ANSWER
               IF (ANSWER .EQ. 'E') RETURN
               IF (ANSWER .NE. 'Y') GO TO 2000
            END IF
         END IF
      END IF
C
      IF (WATCH .GT. 0.5) THEN
         CALL SHOW (PIC(LX,LY), FMAX, SKY(ISTAR), 
     .        MX-LX+1, MY-LY+1, NCOL)
         WRITE (6,622) ID(ISTAR), XCEN(ISTAR), YCEN(ISTAR),
     .        APMAG(ISTAR), NINT(FMAX)
  622    FORMAT (1X, I5, 2F9.2, F9.3, 2X, 'Brightest pixel:', I7)
	       write(*,'(a,$)')'Use this one? '
	       read(*,'(a)')ANSWER
         IF (ANSWER .EQ. 'E') RETURN
         IF (ANSWER .EQ. 'N') GO TO 2000
      END IF
C
C Switch this star (at position ISTAR) with the one at position NSTAR+1
C (which is the highest star in the stack not already a PSF star).
C Estimate the height of the best-fitting profile of type 1 (probably
C a Gaussian) on the basis of the central pixel.
C
      NSTAR = NSTAR+1
      CALL SWAPP (ID, XCEN, YCEN, APMAG, SKY, ISTAR, NSTAR)
      HPSF(NSTAR) = (FMAX-SKY(NSTAR))/
     .     PROFIL(1, 0., 0., PAR, DFDX, DFDY, TERM, 0)
      SATR8D(NSTAR) = SATUR8
      IF (NSTAR .LT. MAXN) GO TO 2000
C
 3000 close (1)					! read *.lst end
      NEXP = NTAB(IEXPAND) + 2*IFRAC
      IF (NSTAR .LT. NEXP) THEN
         write(*,*)
     .     'There aren''t enough PSF stars for a PSF this variable.'
         WRITE (6,*) ' Please change something.'
         RETURN
      END IF
C
C If the first star is saturated, exchange it for the first 
C unsaturated one.
C
      IF (SATR8D(1)) THEN
         DO ISTAR=2,NSTAR
            IF (.NOT. SATR8D(ISTAR)) GO TO 3005
         END DO
         write(*,*)'Every single PSF star is saturated.'
         RETURN
C
 3005    CONTINUE
         CALL SWAPP (ID, XCEN, YCEN, APMAG, SKY, 1, ISTAR)
         SATR8D(ISTAR) = .TRUE.
         SATR8D(1) = .FALSE.
      END IF
C
      CALL OVRWRT ('    Chi     Parameters...', 1)
C
C If OPT(16) (Analytic model PSF) is negative, then try all PSF types
C from 1 to | OPT(16) |, inclusive.
C
      K = NINT(OPT(16))
      J = MAX0(1,K)
      K = IABS(K)
      OLD = 1.E38
      IPSTYP = 0
      DO 3050 I=J,K
C
         DO L=1,NSTAR
            HJNK(L) = HPSF(L)
            XJNK(L) = XCEN(L)
            YJNK(L) = YCEN(L)
         END DO
C
         N = NPARAM(I, FWHM, LABEL, PAR, MAXPAR)
         CALL  FITANA  (PIC, NCOL, NROW, HJNK, XJNK, YJNK, SKY, 
     .       SATR8D, NSTAR, FITRAD, WATCH, I, PAR, N, SIG)
C
C SIG is the root-mean-square scatter about the best-fitting analytic
C function averaged over the central disk of radius FITRAD, expressed
C as a fraction of the peak amplitude of the analytic model.
C
         IF (PAR(1) .LT. 0.) GO TO 3050
C
C If this latest fit is better than any previous one, store the
C parameters of the stars and the model for later reference.
C
         IF (SIG .LT. OLD) THEN
            IPSTYP = I
            WRITE (LINE,201) (PAR(L), L=1,N)
  201       FORMAT (1X, 1P, 6E13.6)
            DO L=1,NSTAR
               HPSF(L) = HJNK(L)
               XCEN(L) = XJNK(L)
               YCEN(L) = YJNK(L)
            END DO
            OLD = SIG
         END IF
 3050 CONTINUE
      IF (IPSTYP .LE. 0) RETURN
      SIG = OLD
      XMID = REAL(NCOL-1)/2.
      YMID = REAL(NROW-1)/2.

        xmid=0.
        ymid=0.
        do 111 i=1,nstar
        xmid=xmid+xcen(i)
111     ymid=ymid+ycen(i)
        xmid=xmid/nstar
        ymid=ymid/nstar
 

      NPAR = NPARAM(IPSTYP, FWHM, LABEL, PAR, MAXPAR)
      READ (LINE,201) (PAR(L), L=1,NPAR)
C
C This last bit ensures that the subsequent computations will be 
C performed with the parameters rounded off exactly the same as they
C appear in the output .PSF file.
C
C=======================================================================
C
C At this point, we have values for the parameters of the best-fitting
C analytic function.  Now we may want to generate the look-up table of
C corrections from the best-fitting analytic function to the actual
C data.  This will be generated within a square box of size NPSF x NPSF
C with half-pixel spacing, centered on the centroid of the star.
C
C First, subtract the analytic function from all the PSF stars in the
C original image.
C
 3300 DP = 0.
      SATUR8 = .FALSE.
      DO 3400 ISTAR=1,NSTAR
         RDX = AMAX1(PROFIL(IPSTYP,0.,0.5,PAR,DFDX,DFDY,TERM,0),
     .               PROFIL(IPSTYP,0.,-0.5,PAR,DFDX,DFDY,TERM,0),
     .               PROFIL(IPSTYP,0.5,0.,PAR,DFDX,DFDY,TERM,0),
     .               PROFIL(IPSTYP,-0.5,0.,PAR,DFDX,DFDY,TERM,0))
         IF (HPSF(ISTAR)*RDX+SKY(ISTAR) .GT. HIBAD) SATR8D(ISTAR)=.TRUE.
C
C The centroids and peak heights are not yet known for saturated
C stars.
C
      IF (SATR8D(ISTAR)) THEN
         SATUR8 = .TRUE.
         GO TO 3400
      END IF
C
      LX = MAX0( 1, INT(XCEN(ISTAR)-PSFRAD)-1 )
      LY = MAX0( 1, INT(YCEN(ISTAR)-PSFRAD)-1 )
      MX = MIN0(NCOL, INT(XCEN(ISTAR)+PSFRAD)+2 )
      MY = MIN0(NROW, INT(YCEN(ISTAR)+PSFRAD)+2 )
      RDX = 0.0
      RDY = 0.0
      DO 3350 J=LY,MY
         DY = REAL(J) - YCEN(ISTAR)
         W = DY**2
         DO I=LX,MX
            DX = REAL(I) - XCEN(ISTAR)
            PIC(I,J) = PIC(I,J) - HPSF(ISTAR) *
     .           PROFIL(IPSTYP, DX, DY, PAR, DFDX, DFDY, TERM, 0)
            IF (DX**2 + W .LT. RADSQ) THEN
               RDX = RDX + (PIC(I,J)-SKY(ISTAR))**2
               RDY = RDY + 1.
            END IF
         END DO
 3350 CONTINUE
      DP = DP + RDY                      ! Total number of pixels
      RDX = SQRT(RDX/RDY)                ! Scatter inside fit radius
      PSF(ISTAR,1,1) = RDX/(HPSF(ISTAR)*
     .     PROFIL(IPSTYP, 0., 0., PAR, DFDX, DFDY, TERM, 0))
 3400 CONTINUE
C
      DP = SQRT(DP/(DP-REAL(NEXP+3*NSTAR))) ! Degrees of freedom
      DO ISTAR=1,NSTAR
         IF (.NOT. SATR8D(ISTAR)) THEN
            PSF(ISTAR,1,1) = DP * PSF(ISTAR,1,1)
            FLAG(ISTAR) = ' '
            IF (SIG .GT. 0.) THEN
               RDX = PSF(ISTAR,1,1)/SIG
               IF (HPSF(ISTAR) .LE. 0) THEN
                  WEIGHT(ISTAR) = 0.
               ELSE
                  WEIGHT(ISTAR) = 1. /(1. + (RDX/2.)**2)
               END IF
               IF (RDX .GE. 3.) THEN
                  FLAG(ISTAR) = '*'
                  if(ifail.eq.0)ifail=1
	          iflag(istar)=1
               ELSE IF (RDX .GE. 2.) THEN
                  FLAG(ISTAR) = '?'
               END IF
            END IF
         END IF
      END DO
      K = (NSTAR-1)/5 + 1                   ! Number of lines
      WRITE (6,67)
   67 FORMAT (/' Profile errors:'/)
      DO I=1,K
         LINE = ' '
         DO ISTAR=I,NSTAR,K
            J = 16*(ISTAR-I)/K + 1
            IF (SATR8D(ISTAR)) THEN
               WRITE (LINE(J:J+15),68) ID(ISTAR)
   68          FORMAT (1X, I5, ' saturated')
            ELSE
               WRITE (LINE(J:J+15),69) ID(ISTAR), 
     .              PSF(ISTAR,1,1), FLAG(ISTAR)
   69          FORMAT (1X, I5, F7.3, 1X, A1, 1X)
            END IF
         END DO
         WRITE (6,*) LINE(1:79)
      END DO
        if(ifail.eq.1)return
      SCALE = HPSF(1)
      DO ISTAR=NSTAR,1,-1
         HPSF(ISTAR) = HPSF(ISTAR)/HPSF(1)
      END DO
	 
C
C Tabulate the constant part of the PSF.
C
      MIDDLE = (NPSF+1)/2
      MIDSQ = MIDDLE**2
      NPSFSQ = ((NPSF-1)/2)**2
      DO J=1,NPSF
         JDY = J-MIDDLE
         JDYSQ = JDY**2
         DY = REAL(JDY)/2.
         DO I=1,NPSF
            IDX = I-MIDDLE
            K = IDX**2 + JDYSQ
            IF (K .LE. MIDSQ) THEN
               IN(I,J) = .TRUE.
               DX = REAL(IDX)/2.
               CON(I,J) = 
     .              PROFIL(IPSTYP, DX, DY, PAR, DFDX, DFDY, V, 0)
               IF (K .GE. NPSFSQ) EDGE(I,J) = .TRUE.
            ELSE
               IN(I,J) = .FALSE.
               DO K=1,NEXP
                  PSF(I,J,K) = 0.
               END DO
            END IF
         END DO
      END DO
      IF (IEXPAND+IFRAC .LT. 0) THEN
         PSFMAG = 0.0D0
         DO J=1,NPSF
            JDY = J-MIDDLE
            JDYSQ = JDY**2
            DO I=1,NPSF
               IF (IN(I,J)) THEN
                  IDX = I-MIDDLE
                  IF (IDX**2+JDYSQ .LE. IRADSQ) PSFMAG = PSFMAG + 
     .                 CON(I,J)
               END IF
            END DO
         END DO
         PSFMAG = 25.D0-2.5D0*DLOG10(SCALE*PSFMAG/4.D0)
         GO TO 2900
      END IF
C
C=======================================================================
C
C Now compute look-up tables.
C
 3500 CONTINUE
C
C MIDDLE is the center of the look-up table, which will correspond to
C the centroid of the analytic PSF.  MIDSQ is the square of the
C radius of the PSF table, plus a pixel's worth of slack just to
C make sure you'll always have a 4x4 array suitable for interpolation.
C
      DO K=1,NEXP
         DO J=1,NPSF
            DO I=1,NPSF
               PSF(I,J,K) = 0.0
            END DO
         END DO
      END DO
C
      NPASS = AMIN0(NSTAR-1,NPASS)
      DO 4950 JY=1,NPSF
         IF (WATCH .GT. -1.5) THEN
            WRITE (LINE,*) '  Computed', JY, '  rows of', NPSF, 
     .           '  in the PSF.'
            CALL OVRWRT (LINE(1:60), 2)
         END IF
         RDY = REAL(JY-MIDDLE)/2.
C
         DO 4940 IX=1,NPSF
C
C Don't waste time with pixels outside a radius = (MIDDLE-1).
C
            IF (.NOT. IN(IX,JY)) GO TO 4940
            RDX = REAL(IX-MIDDLE)/2.
            DO 4935 ITER=0,NPASS
C
C Initialize accumulators.
C
            TERM(1) = 1.
            SUMSQ = 0.
            SUMN = 0.
            DO L=1,NEXP
               V(L) = 0.
               DO K=1,NEXP
                  C(K,L) = 0.
               END DO
            END DO
C
C Now, for this point in the set of lookup table(s) [(IX,JY), where
C (MIDDLE,MIDDLE) corresponds to the centroid of the PSF], consult
C each of the PSF stars to determine the value(s) to put into the
C table(s).
C
            DO 4900 ISTAR = 1,NSTAR
               IF (SATR8D(ISTAR)) GO TO 4900
C
C What pixels in the original image constitute a 4x4 box surrounding
C this (IX,JY) in the PSF, as referred to this star's centroid?
C
               DX = XCEN(ISTAR) + RDX
               DY = YCEN(ISTAR) + RDY
               I = INT(DX)
               J = INT(DY)
C
C The 4x4 box is given by  I-1 <= x <= I+2, J-1 <= y <= J+2.
C
               IF ((I .LT. 2) .OR. (J .LT. 2) .OR. 
     .              (I+2 .GT. NCOL) .OR.
     .              (J+2 .GT. NROW)) GO TO 4900           ! Next star
C
               DO L=J-1,J+2
                  DO K=I-1,I+2
                     IF ((PIC(K,L) .GT. HIBAD) .OR.
     .                    (PIC(K,L) .LT. -HIBAD)) GO TO 4900
                  END DO
               END DO
C
               DX = DX - I
               DY = DY - J
C
C The point which corresponds PRECISELY to the offset (RDX,RDY)
C from the star's centroid, lies a distance (DX,DY) from pixel
C (I,J).
C
C Use bicubic interpolation to evaluate the residual PSF amplitude
C at this point.  Scale the residual up to match the first PSF star.
C
               DP = BICUBC(PIC(I-1,J-1), NCOL, DX, DY, DFDX, DFDY)
     .              - SKY(ISTAR)
               DP = DP/HPSF(ISTAR)
               SUMSQ = SUMSQ + ABS(DP)
               SUMN = SUMN + 1.
               IF (IEXPAND .GE. 1) THEN
                  TERM(2) = (XCEN(ISTAR)-1.)/XMID-1.
                  TERM(3) = (YCEN(ISTAR)-1.)/YMID-1.
                  IF (IEXPAND .GE. 2) THEN
                     TERM(4) = 1.5*TERM(2)**2-0.5
                     TERM(5) = TERM(2)*TERM(3)
                     TERM(6) = 1.5*TERM(3)**2-0.5
                     IF (IEXPAND .GE. 3) THEN
                        TERM(7) = TERM(2)*(5.*TERM(4)-2.)/3.
                        TERM(8) = TERM(4)*TERM(3)
                        TERM(9) = TERM(2)*TERM(6)
                        TERM(10) = TERM(3)*(5.*TERM(6)-2.)/3.
                     END IF
                  END IF
               END IF
C
C               IF (IFRAC .GE. 1) THEN
C
C INSERT CODE HERE
C
               IF (ITER .GT. 0) THEN
                  OLD = 0.
                  DO K=1,NEXP
                     OLD = OLD + PSF(IX,JY,K) * TERM(K)
                  END DO
C
                  IF (ITER .LE. MAX0(3,NPASS/2)) THEN
                     W = HPSF(ISTAR)/(1. + (ABS(DP - OLD)/SIG))
                  ELSE 
                     W = HPSF(ISTAR)/(1. + ((DP - OLD)/SIG)**2)
                  END IF
               ELSE
                  W = HPSF(ISTAR)
               END IF
C
               W = W*WEIGHT(ISTAR)
               DO K=1,NEXP
                  WT = W*TERM(K)
                  V(K) = V(K) + WT*DP
                  DO L=1,NEXP
                     C(K,L) = C(K,L) + WT*TERM(L)
                  END DO
               END DO
 4900       CONTINUE
C
            IF (SUMN .LT. NEXP) THEN
               write(*,*)'Not enough PSF stars.  Please start over.'
               close(3)
               RETURN
            END IF
            CALL INVERS (C, MAXEXP, NEXP, ISTAT)
            CALL VMUL (C, MAXEXP, NEXP, V, A)
            DO K=1,NEXP
               PSF(IX,JY,K) = A(K)
            END DO
            IF (SUMN .LE. NEXP) GO TO 4940
            SIG = 1.2533*SUMSQ/SQRT(SUMN*(SUMN - NEXP))
 4935       CONTINUE
 4940    CONTINUE
 4950 CONTINUE
C
C Make the average value outside the PSF radius identically zero.
C
      DX = 0.
      K = 0
      DO J=1,NPSF
         DO I=1,NPSF
            IF (EDGE(I,J)) THEN
               DX = DX+SCALE*CON(I,J)+PSF(I,J,1)
               K = K+1
            END IF
         END DO
      END DO
      DX = DX/REAL(K)
      DO J=1,NPSF
         DO I=1,NPSF
            IF (IN(I,J)) PSF(I,J,1) = PSF(I,J,1)-DX
         END DO
      END DO
      IF (NEXP .LT. 2) GO TO 2800
C
C At this point, we must be sure that any higher order terms
C in the PSF [i.e. those that go as powers of (XCEN-XMID) and 
C (YCEN-YMID)] contain zero volume, so that the total volume
C of the PSF is independent of position.  This will be done
C by looking at the total flux contained in the look-up
C tables of index 2 and higher.  The analytic function
C will be scaled to the same total flux and subtracted from
C each of the higher lookup tables and added into the
C lookup table of index 1.  I scale the analytic function
C instead of simply transferring the net flux from the 
C higher tables to table 1, because it is poor fits of
C the analytic profile to the image data which has caused
C the net flux in the lookup tables to depend upon position.
C
C Compute the brightness weighted average value of each of 
C the terms in the polynomial expansion.
C
      DO K=1,NEXP
         TERM(K) = 0.
      END DO
C
      DO ISTAR=NSTAR,1,-1
         DX = (XCEN(ISTAR)-1.)/XMID-1.
         DY = (YCEN(ISTAR)-1.)/YMID-1.
         W = WEIGHT(ISTAR)*HPSF(ISTAR)
         TERM(1) = TERM(1) + W
         TERM(2) = TERM(2) + W*DX
         TERM(3) = TERM(3) + W*DY
         IF (IEXPAND .GE. 2) THEN
            DX2 = DX**2
            DY2 = DY**2
            TERM(4) = TERM(4) + W*(1.5*DX2-0.5)
            TERM(5) = TERM(5) + W*(DX*DY)
            TERM(6) = TERM(6) + W*(1.5*DY2-0.5)
            IF (IEXPAND .GE. 3) THEN
               TERM(7) = TERM(7) + DX*(5.*DX2-3.)/2.
               TERM(8) = TERM(8) + DY*(1.5*DX2-0.5)
               TERM(9) = TERM(9) + DX*(1.5*DY2-0.5)
               TERM(10) = TERM(10) + DY*(5.*DY2-3.)/2.
            END IF
         END IF
C
C            IF (IFRAC .GE. 1) THEN
C
C INSERT CODE HERE
C
      END DO
C
      DO K=NEXP,1,-1
         TERM(K) = TERM(K)/TERM(1)
      END DO
C
C Okey doke.  Now if there is any net volume contained in any of the
C higher order (variable) terms of the PSF, we will remove it in two
C components.  First of all, there is the possibility that there were some
C errors in the sky estimates which was correlated with position.
C Then there is the possibility that there was some systematic error
C in estimating HPSF across the frame.  We correct these by removing
C any net volume in the higer-order look-up tables and by removing
C any correlation with the analytic profile.  We want to make sure
C we do this only for the area within one PSF radius.
C
C
      DO K=2,NEXP
         DO JJ=0,1
            DO II=0,1
               DO J=1,2
                  VV(J) = 0.D0
                  DO I=1,2
                     CC(I,J) = 0.D0
                  END DO
               END DO
               DO J=1+JJ,NPSF-JJ,2
                  JDY = J-MIDDLE
                  JDYSQ = JDY**2
                  DO I=1+II,NPSF-II,2
                     IF (IN(I,J)) THEN
                        IDX = I-MIDDLE
                        IR = IDX**2+JDYSQ
                        IF (IR .LE. IRADSQ) THEN
                           VV(1) = VV(1)+PSF(I,J,K)
cxyz                       VV(2) = VV(2)+PSF(I,J,K)*CON(I,J)
                           CC(1,1) = CC(1,1)+1.
                           CC(2,1) = CC(2,1)+CON(I,J)
cxyz                       CC(2,2) = CC(2,2)+CON(I,J)**2
                        END IF
                        IF (IR .GE. JRADSQ) THEN
                           VV(2) = VV(2)+PSF(I,J,K)
                           CC(1,2) = CC(1,2)+1.
                           CC(2,2) = CC(2,2)+CON(I,J)
                        END IF
                     END IF
                  END DO
               END DO
cxyz           CC(1,2) = CC(2,1)
               CALL DINVRS (CC, 2, 2, ISTAT)
               CALL DVMUL (CC, 2, 2, VV, AA)
               SCALE = SCALE + TERM(K)*AA(2)
               DX = TERM(K)*AA(1)
               DO J=1+JJ,NPSF-JJ,2
                  JDY = J-MIDDLE
                  JDYSQ = JDY**2
                  DO I=1+II,NPSF-II,2
                     IF (IN(I,J)) THEN
                        IDX = I-MIDDLE
                        PSF(I,J,K) = PSF(I,J,K) -
     .                       AA(1)-CON(I,J)*AA(2)
                        PSF(I,J,1) = PSF(I,J,1) + DX
                     END IF
                  END DO
               END DO
            END DO
         END DO
      END DO
C
 2799 continue
C
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
 2800 CONTINUE
      PSFMAG = 0.0D0
      DO J=1,NPSF
         JDY = J-MIDDLE
         JDYSQ = JDY**2
         DO I=1,NPSF
            IF (IN(I,J)) THEN
               IDX = I-MIDDLE
               IF (IDX**2+JDYSQ .LE. IRADSQ) THEN
                  PSFMAG = PSFMAG + SCALE*CON(I,J) + DBLE(PSF(I,J,1))
               END IF
            END IF
         END DO
      END DO
      if(PSFMAG.lt.0)then
        inquire(file='psf.mag',exist=logi)
        if(.not.logi)then
	  PSFMAG=12.5
	else
	  open(55,file='psf.mag',status='old')
	  read(55,*)PSFMAG
	  PSFMAG=PSFMAG+0.5
	  close(55)
	endif
      else
        PSFMAG = 25.D0-2.5D0*DLOG10(PSFMAG/4.D0)
        open(55,file='psf.mag',status='unknown')
	write(55,*)PSFMAG
	close(55)
      endif	
C
C If there are saturated stars, we must now guesstimate their positions
C and peak brightnesses.  Determine these from single-profile fits a la
C the PEAK routine, with a fitting radius equal to half the PSF radius.
C Subtract the analytic part of their profile from the image, and
C go back and redetermine the lookup tables.
C
      IF (SATUR8 .AND. (OPT(18) .GT. 0.5)) THEN
         PERR = 0.01*OPT(19)
         PKERR = 0.01*OPT(20)/(PAR(1)*PAR(2)) ! **2
         DO ISTAR=1,NSTAR
            IF (SATR8D(ISTAR)) THEN
               LX = MAX0( 1, INT(XCEN(ISTAR)-PSFRAD)-1 )
               LY = MAX0( 1, INT(YCEN(ISTAR)-PSFRAD)-1 )
               MX = MIN0(NCOL, INT(XCEN(ISTAR)+PSFRAD)+2 )
               MY = MIN0(NROW, INT(YCEN(ISTAR)+PSFRAD)+2 )
               K = MX-LX+1
               L = MY-LY+1
               DX = (XCEN(ISTAR)-1.)/XMID-1.
               DY = (YCEN(ISTAR)-1.)/YMID-1.
               HPSF(ISTAR) = 3.*HPSF(1)            ! Lousy starting guess
               XCEN(ISTAR) = XCEN(ISTAR)-LX+1.
               YCEN(ISTAR) = YCEN(ISTAR)-LY+1.
               CALL PKFIT (PIC(LX,LY), K, L, NCOL, XCEN(ISTAR), 
     .              YCEN(ISTAR), HPSF(ISTAR), SKY(ISTAR), 0.5*PSFRAD,
     .              LOBAD, HIBAD, SCALE, IPSTYP, PAR, MAXPAR, NPAR,
     .              PSF, MAXPSF, MAXEXP, NPSF, NEXP, IFRAC, DX,
     .              DY, ERRMAG, CHI, SHARP, N, NCOL, NROW)
               XCEN(ISTAR) = XCEN(ISTAR)+LX-1.
               YCEN(ISTAR) = YCEN(ISTAR)+LY-1.
               WRITE (6,666) ID(ISTAR), XCEN(ISTAR), YCEN(ISTAR), 
     .              PSFMAG-2.5*ALOG10(HPSF(ISTAR)), 
     .              AMIN1(2.0, 1.086*ERRMAG/HPSF(ISTAR)), 
     .              REAL(N), CHI, SHARP
  666               FORMAT (1X, I5, 4F9.3, F9.0, F9.2, F9.3)
               LX = MAX0( 1, INT(XCEN(ISTAR)-PSFRAD)-1 )
               LY = MAX0( 1, INT(YCEN(ISTAR)-PSFRAD)-1 )
               MX = MIN0(NCOL, INT(XCEN(ISTAR)+PSFRAD)+2 )
               MY = MIN0(NROW, INT(YCEN(ISTAR)+PSFRAD)+2 )
               RDX = SCALE*HPSF(ISTAR)
               DO J=LY,MY
                  DY = REAL(J) - YCEN(ISTAR)
                  W = DY**2
                  DO I=LX,MX
                     IF (PIC(I,J) .LE. HIBAD) THEN
                        DX = REAL(I) - XCEN(ISTAR)
                        RDY = RDX * PROFIL(IPSTYP, 
     .                       DX, DY, PAR, DFDX, DFDY, TERM, 0)
                        IF (RDY+SKY(ISTAR) .LE. HIBAD) THEN
                           PIC(I,J) = PIC(I,J) - RDY
                        ELSE
                           PIC(I,J) = RDY+SKY(ISTAR)
C
C Ensuring it is greater than HIBAD, and hence will be ignored.
C
                        END IF
                     END IF
                  END DO
               END DO
               SATR8D(ISTAR) = .FALSE.
               WEIGHT(ISTAR) = 0.5
            END IF
         END DO
         SATUR8 = .FALSE.
         GO TO 3500
      END IF
C
C---------------------------------------------------------------------
C
 2900 CONTINUE
C
C Write PSF file header.
C
      WRITE (3,202) LABEL, NPSF, NPAR, NTAB(IEXPAND), 5*IFRAC,
     .     PSFMAG, SCALE, XMID, YMID
      
	      WRITE (*,202) LABEL, NPSF, NPAR, NTAB(IEXPAND), 5*IFRAC,
     .     PSFMAG, SCALE, XMID, YMID
  202 FORMAT (1X, A8, 4I5, F9.3, F15.3, 2F9.1)
      WRITE (3,201) (PAR(L), L=1,NPAR)
C
      IF (NEXP .GT. 0) THEN
         DO K=1,NEXP
            WRITE (3,204) ((PSF(I,J,K), I=1,NPSF), J=1,NPSF)
  204       FORMAT (1X, 1P, 6E13.6)
         END DO
      END IF
      close(3)
C
C Now generate the file of neighbors.
C
C First, write out the PSF stars.
C
      CALL WRHEAD (4, 3, NCOL, NROW, 7, LOBAD, HIBAD, THRESH, AP1, 
     .     PHPADU, READNS, RADIUS)
      DO I=1,NSTAR
         WRITE (4,280) ID(I), XCEN(I), YCEN(I), APMAG(I), SKY(I)
 280     FORMAT (1X, I5, 3F9.3, f11.3)
      END DO
      IF (NSTAR .GE. NTOT) GO TO 9000
C
C Now look for all stars within a radius equal to 
C 1.5*PSFRAD+2.*FITRAD+1 of any of the PSF stars.
C
      RSQ = (1.5*PSFRAD+2.*FITRAD+1.)**2

      L = NSTAR+1
C
C L will point at the first irrelevant star.
C
      DO I=1,NSTAR
         J = L
 5020    IF ((XCEN(J)-XCEN(I))**2 + (YCEN(J)-YCEN(I))**2 .LE. RSQ) THEN
            CALL SWAPP (ID, XCEN, YCEN, APMAG, SKY, J, L)
            WRITE (4,280) ID(L), XCEN(L), YCEN(L), APMAG(L), SKY(L)
            L = L+1
         END IF
         IF (J .LT. NTOT) THEN
            J = J+1
            GO TO 5020
         END IF
      END DO
C
C Now all stars in the stack between positions NSTAR+1 and L-1,
C inclusive, lie within the specified radius of the PSF stars and have
C been written out.  Finally, look for any stars within a radius 
C equal to 2*FITRAD+1 of any of THESE stars, and of each other.
C
      IF (L .GT. NTOT) GO TO 9000
      RSQ = (2.*FITRAD+1)**2
      I = NSTAR+1
 5110 J = L
 5120 IF ((XCEN(J)-XCEN(I))**2 + (YCEN(J)-YCEN(I))**2 .LE. RSQ) THEN
         CALL SWAPP (ID, XCEN, YCEN, APMAG, SKY, J, L)
         WRITE (4,280) ID(L), XCEN(L), YCEN(L), APMAG(L), SKY(L)
         L = L+1
      END IF
      IF (J .LT. NTOT) THEN
         J = J+1
         GO TO 5120
      END IF
C
      I = I+1
      IF (I .LT. L) GO TO 5110
 9000 close(4)
      RETURN
      END!
C
C#######################################################################
C
      SUBROUTINE  FITANA  (PIC, NCOL, NROW, H, XCEN, YCEN, SKY, 
     .     SATR8D, NSTAR, 
     .    FITRAD, WATCH, IPSTYP, PAR, NPAR, CHI)
C
C This subroutine fits the ANALYTIC profile to the selected stars.
C
C  OFFICIAL DAO VERSION:    1991 May 10
C
      IMPLICIT NONE
      INTEGER MAXPAR, MAXBOX, MAXN, NCOL, NROW
      PARAMETER  (MAXPAR=6, MAXBOX=69, MAXN=200)
      CHARACTER*80 LINE
      REAL PIC(NCOL,NROW)
      REAL H(*), XCEN(*), YCEN(*), SKY(*), CLAMP(MAXPAR), OLD(MAXPAR)
      REAL C(MAXPAR,MAXPAR), V(MAXPAR), PAR(MAXPAR), T(MAXPAR), 
     .     Z(MAXPAR)
      REAL XCLAMP(MAXN), YCLAMP(MAXN), XOLD(MAXN), YOLD(MAXN)
      REAL ABS, PROFIL
      LOGICAL SATR8D(*)
C
      REAL DHN, DHD, DXN, DXD, DYN, DYD, DX, DY, DYSQ
      REAL PEAK, DP, PROD, DHDXC, DHDYC, P, WT
      REAL RSQ, SUMWT, OLDCHI, CHI, WATCH, FITRAD
      INTEGER I, J, K, L, LX, LY, MX, MY, ISTAT
      INTEGER ISTAR, MPAR, NITER, NPAR, IPSTYP, NSTAR
      LOGICAL FULL
C
C-----------------------------------------------------------------------
C
      FULL = .FALSE.
      DO I=1,NPAR
         Z(I) = 0.
         OLD(I) = 0.
         CLAMP(I) = 0.5
      END DO
      CLAMP(1) = 2.
      CLAMP(2) = 2.
      NITER = 0
      OLDCHI = 0.
      DO I=1,NSTAR
         XOLD(I) = 0.
         YOLD(I) = 0.
         XCLAMP(I) = 1.
         YCLAMP(I) = 1.
      END DO
C
      MPAR = 2
C
C-----------------------------------------------------------------------
C
C SECTION 1
C
C Now we will fit an integrated analytic function to the central part 
C of the stellar profile.  For each star we will solve for three 
C parameters: (1) H, the central height of the model profile (above 
C sky); (2) XCEN, C the centroid of the star in x; and (3) YCEN, 
C likewise for y.  In addition, from ALL STARS CONSIDERED TOGETHER
C we will determine any other parameters, PAR(i), required to describe 
C the profile.  We will use a circle of radius FITRAD centered on the 
C position of each PSF star. NOTE THAT we do not fit the data to an
C actual analytic profile, but rather the function is numerically
C integrated over the area of each pixel, and the observed data are fit
C to these integrals.
C
      RSQ = FITRAD**2
 1000 NITER = NITER+1
      IF (NITER .GT. 600) THEN
         write(*,*)'Failed to converge.'
         PAR(1) = -1.
         RETURN
      END IF
C
C Initialize the big accumulators.
C
      DO I=1,NPAR
         V(I)=0.0                        ! Zero the vector of residuals
         DO J=1,NPAR
            C(J,I)=0.0                   ! Zero the normal matrix
         END DO
      END DO
      CHI = 0.0
      SUMWT = 0.0
C
C Using the analytic model PSF defined by the current set of parameters,
C compute corrections to the brightnesses and centroids of all the PSF
C stars.  MEANWHILE, accumulate the corrections to the model parameters.
C
      DO 1400 ISTAR=1,NSTAR
      IF (SATR8D(ISTAR)) GO TO 1400
      LX = INT(XCEN(ISTAR)-FITRAD) + 1
      LY = INT(YCEN(ISTAR)-FITRAD) + 1
      MX = INT(XCEN(ISTAR)+FITRAD)
      MY = INT(YCEN(ISTAR)+FITRAD)
C
      DHN = 0.0D0
      DHD = 0.0D0
      DXN = 0.0D0
      DXD = 0.0D0
      DYN = 0.0D0
      DYD = 0.0D0
      DO J=LY,MY
         DY = REAL(J) - YCEN(ISTAR)
         DYSQ = DY**2
         DO I=LX,MX
            DX = REAL(I) - XCEN(ISTAR)
            WT = (DX**2+DYSQ)/RSQ
            IF (WT .LT. 1.) THEN
               P = PROFIL(IPSTYP, DX, DY, PAR, DHDXC, DHDYC, T, 0)
               DP = PIC(I,J) - H(ISTAR)*P - SKY(ISTAR)
c              if (istar .le. 2) type 6661, i, j, pic(i,j),
c    .              h(istar)*p, sky(istar), dp
c6661          format (2i6, 4f12.3)
               DHDXC = H(ISTAR)*DHDXC
               DHDYC = H(ISTAR)*DHDYC
               WT = 5./(5.+WT/(1.-WT))
               PROD = WT*P
               DHN = DHN + PROD*DP
               DHD = DHD + PROD*P
               PROD = WT*DHDXC
               DXN = DXN + PROD*DP
               DXD = DXD + PROD*DHDXC
               PROD = WT*DHDYC
               DYN = DYN + PROD*DP
               DYD = DYD + PROD*DHDYC
            END IF
         END DO
      END DO
C
c     if (istar .le. 2) type 6662, h(istar), dhn, dhd, 
c    .     h(istar)+dhn/dhd
c6662 format (4f12.3)
c     if (istar .eq. 2) type *
c     if (istar .eq. 2) accept *
      H(ISTAR) = H(ISTAR) + DHN/DHD
C
      DXN = DXN/DXD
      IF (XOLD(ISTAR)*DXN .LT. 0.) XCLAMP(ISTAR) = 0.5*XCLAMP(ISTAR)
      XOLD(ISTAR) = DXN
      XCEN(ISTAR) = XCEN(ISTAR)+DXN/(1.+ABS(DXN)/XCLAMP(ISTAR))
C
      DYN = DYN/DYD
      IF (YOLD(ISTAR)*DYN .LT. 0.) YCLAMP(ISTAR) = 0.5*YCLAMP(ISTAR)
      YOLD(ISTAR) = DYN
      YCEN(ISTAR) = YCEN(ISTAR)+DYN/(1.+ABS(DYN)/YCLAMP(ISTAR))
C
      PEAK = H(ISTAR) * PROFIL(IPSTYP,0.,0., PAR, DHDXC, DHDYC, T, 0)
      DO J=LY,MY
         DY = REAL(J)-YCEN(ISTAR)
         DYSQ = DY**2
         DO I=LX,MX
            DX = REAL(I)-XCEN(ISTAR)
            WT = (DX**2+DYSQ)/RSQ
            IF (WT .LT. 1.) THEN
               P = PROFIL(IPSTYP, DX, DY, PAR, DHDXC, DHDYC, T, 1)
               DP = PIC(I,J) - H(ISTAR)*P - SKY(ISTAR)
               DO K=1,MPAR
                  T(K) = H(ISTAR)*T(K)
               END DO
               CHI = CHI + ABS(DP/PEAK)
               SUMWT = SUMWT + 1.
               WT = 5./(5.+WT/(1.-WT))
               IF (NITER .GE. 4) WT = WT/(1.+ABS(20.*DP/PEAK))
               DO K=1,MPAR
                  V(K) = V(K) + WT*DP*T(K)
                  DO L=1,MPAR
                     C(L,K) = C(L,K) + WT*T(L)*T(K)
                  END DO
               END DO
            END IF
         END DO
      END DO
 1400 CONTINUE
C
C Correct the fitting parameters.
C
      CALL INVERS (C, MAXPAR, MPAR, ISTAT)
      CALL VMUL (C, MAXPAR, MPAR, V, Z)
      DO I=1,MPAR
         IF (Z(I)*OLD(I) .LT. 0.) THEN
            CLAMP(I) = 0.5*CLAMP(I)
         ELSE
            CLAMP(I) = 1.1*CLAMP(I)
         END IF
         OLD(I) = Z(I)
         Z(I) = CLAMP(I)*Z(I)
      END DO
      Z(1) = AMAX1(-0.1*PAR(1), AMIN1(0.1*PAR(1), Z(1)))
      Z(2) = AMAX1(-0.1*PAR(2), AMIN1(0.1*PAR(2), Z(2)))
C     Z(3) = Z(3)/(1.+ABS(Z(3))/(AMIN1(0.1,1.-ABS(PAR(3)))))
      DO I=1,MPAR
         PAR(I) = PAR(I)+Z(I)
      END DO
C
      SUMWT = SUMWT - REAL(MPAR + 3*NSTAR)
      IF (SUMWT .GT. 0) THEN
         CHI = 1.2533*CHI/SUMWT
      ELSE
         CHI = 9.9999
      END IF
C
      WRITE (LINE,661) CHI, (PAR(I), I=1,MPAR)
  661 FORMAT (1X, F7.4, 6F10.5)
      I=10*MPAR+13
      IF (WATCH .GT. -1.5) CALL OVRWRT (LINE(1:I), 2)
      IF (MPAR .EQ. NPAR) THEN
         IF (ABS(OLDCHI/CHI-1.) .LT. 1.E-5) THEN
            CALL OVRWRT (LINE(1:I), 3)
            RETURN
         END IF
      ELSE
         IF (ABS(OLDCHI/CHI-1.) .LT. 1.E-3) THEN
            MPAR = MPAR+1
            OLDCHI = 0.
            GO TO 1000
         END IF
      END IF
C
      OLDCHI = CHI
      GO TO 1000
      END!
C
C#######################################################################
C
C Exchange two stars in the star list.
C
C   OFFICIAL DAO VERSION:       1991 May 10
C
      SUBROUTINE SWAPP (ID, X, Y, A, S, I, J)
      IMPLICIT NONE
      REAL X(*), Y(*), A(*), S(*)
      INTEGER ID(*)
C
      REAL HOLD
      INTEGER I, J, IHOLD
C
      IHOLD = ID(I)
      ID(I) = ID(J)
      ID(J) = IHOLD
      HOLD = X(I)
      X(I) = X(J)
      X(J) = HOLD
      HOLD = Y(I)
      Y(I) = Y(J)
      Y(J) = HOLD
      HOLD = A(I)
      A(I) = A(J)
      A(J) = HOLD
      HOLD = S(I)
      S(I) = S(J)
      S(J) = HOLD
      RETURN
      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!
C
C This file contains subroutines that are not I/O related, but rather
C involve arithmetic or character operations that my be somewhat 
C machine-specific-- it may be necessary or desirable to make some
C changes to optimize the code to run on your computer.
C
C              OFFICIAL DAO VERSION:  1991 April 18
C
C***********************************************************************
C
C Current contents:
C
C   INVERS  inverts a square matrix.
C
C     VMUL  multiplies a square matrix (on the left) by a column vector
C           (on the right) yielding a new column vector.
C
C    QUICK  does a quicksort on a vector of data.
C
C     SHOW  produces what is effectively an eleven-level gray scale plot
C           of a rectangular data array on the terminal.
C
C***********************************************************************
C
      SUBROUTINE  INVERS (A, MAX, N, IFLAG)
C
C Although it seems counter-intuitive, the tests that I have run
C so far suggest that the 180 x 180 matrices that NSTAR needs can
C be inverted with sufficient accuracy if the elements are REAL
C rather than DOUBLE PRECISION
C
C Arguments
C 
C     A (INPUT/OUTPUT) is a square matrix of dimension N.  The inverse 
C       of the input matrix A is returned in A.
C
C   MAX (INPUT) is the size assigned to the matrix A in the calling 
C       routine.  It's needed for the dimension statement below.
C
C IFLAG (OUTPUT) is an error flag.  IFLAG = 1 if the matrix could not
C       be inverted; IFLAG = 0 if it could.
C
      IMPLICIT NONE
      INTEGER MAX
      REAL A(MAX,MAX)
C
      INTEGER N, IFLAG, I, J, K
C
C-----------------------------------------------------------------------
C
      IFLAG=0
      I=1
  300 IF(A(I,I).EQ.0.0E0)GO TO 9100
      A(I,I)=1.0E0/A(I,I)
      J=1
  301 IF(J.EQ.I)GO TO 304
      A(J,I)=-A(J,I)*A(I,I)
      K=1
  302 IF(K.EQ.I)GO TO 303
      A(J,K)=A(J,K)+A(J,I)*A(I,K)
  303 IF(K.EQ.N)GO TO 304
      K=K+1
      GO TO 302
  304 IF(J.EQ.N)GO TO 305
      J=J+1
      GO TO 301
  305 K=1
  306 IF(K.EQ.I)GO TO 307
      A(I,K)=A(I,K)*A(I,I)
  307 IF(K.EQ.N)GO TO 308
      K=K+1
      GO TO 306
  308 IF(I.EQ.N)RETURN                                   ! Normal return
      I=I+1
      GO TO 300
C
C-----------------------------------------------------------------------
C
C Error:  zero on the diagonal.
C
 9100 IFLAG=1
      RETURN
C
      END!
C
C#######################################################################
C
      SUBROUTINE  VMUL (A, MAX, N, V, X)
C
C Multiply a matrix by a vector:
C
C                    A * V = X
C
C Arguments
C
C    A(column,row)  (INPUT) is a square matrix of dimension N.
C
C              MAX  (INPUT) is the size assigned to the array in the 
C                           calling routine.
C
C           V(row)  (INPUT) is a column vector of dimension N.
C
C           X(row) (OUTPUT) is a column vector of dimension N.
C
      IMPLICIT NONE
      INTEGER MAX
      REAL A(MAX,MAX), V(MAX)
      REAL X(MAX)
C
      DOUBLE PRECISION DBLE
      REAL SNGL
C
      DOUBLE PRECISION SUM
      INTEGER N, I, J
C
C-----------------------------------------------------------------------
C
      I=1
  200 SUM=0.0D0
      J=1
  201 SUM=SUM+DBLE(A(J,I))*DBLE(V(J))
      IF (J .EQ. N) GO TO 203
      J=J+1
      GO TO 201
  203 X(I)=SNGL(SUM)
      IF (I .EQ. N) RETURN                               ! Normal return
      I=I+1
      GO TO 200
      END!
C
C#######################################################################
C
      SUBROUTINE  DINVRS (A, MAX, N, IFLAG)
      IMPLICIT NONE
      INTEGER MAX
      DOUBLE PRECISION A(MAX,MAX)
C
      INTEGER I, J, K, N, IFLAG
C
C-----------------------------------------------------------------------
C
      IFLAG=0
      I=1
  300 IF(A(I,I).EQ.0.0D0)GO TO 9100
      A(I,I)=1.0D0/A(I,I)
      J=1
  301 IF(J.EQ.I)GO TO 304
      A(J,I)=-A(J,I)*A(I,I)
      K=1
  302 IF(K.EQ.I)GO TO 303
      A(J,K)=A(J,K)+A(J,I)*A(I,K)
  303 IF(K.EQ.N)GO TO 304
      K=K+1
      GO TO 302
  304 IF(J.EQ.N)GO TO 305
      J=J+1
      GO TO 301
  305 K=1
  306 IF(K.EQ.I)GO TO 307
      A(I,K)=A(I,K)*A(I,I)
  307 IF(K.EQ.N)GO TO 308
      K=K+1
      GO TO 306
  308 IF(I.EQ.N)RETURN                                   ! Normal return
      I=I+1
      GO TO 300
C
C-----------------------------------------------------------------------
C
C Error:  zero on the diagonal.
C
 9100 IFLAG=1
      RETURN
C
      END!
C
C#######################################################################
C
      SUBROUTINE  DVMUL (A, MAX, N, V, X)
      IMPLICIT NONE
      INTEGER MAX
      DOUBLE PRECISION A(MAX,MAX), V(MAX), X(MAX)
C
      DOUBLE PRECISION SUM
      INTEGER I, J, N
C
C-----------------------------------------------------------------------
C
      I=1
  200 SUM=0.0D0
      J=1
  201 SUM=SUM+A(J,I)*V(J)
      IF (J .EQ. N) GO TO 203
      J=J+1
      GO TO 201
  203 X(I)=SUM
      IF (I .EQ. N) RETURN                               ! Normal return
      I=I+1
      GO TO 200
      END!
C
C#######################################################################
C
 	SUBROUTINE  SHOW (F, FMAX, FZERO, NX, NY, NCOL)
C
C=======================================================================
C
C A simple subroutine to use alphanumeric characters to produce a 2-D
C gray-scale plot on an alphanumeric terminal.
C
C=======================================================================
C
      IMPLICIT NONE
      INTEGER NCHAR, MAXPLT, NCOL
      PARAMETER  (NCHAR=11, MAXPLT=36)
C
C Parameters
C
C   NCHAR is the number of discrete gray levels we wish to produce.
C
C MAXPLT is the widest plot that can be produced on the terminal
C         screen.  Since two characters will be typed out per pixel (to
C         make the overall plot more nearly square) MAXPLT should be 
C         equal to (N-2)/2 where N is the number of character positions 
C         per line on the screen; the extra two characters will be used
C         for vertical bars ('|') to delimit the picture.  Arrays
C         which are more than MAXPLT pixels on a side will be 
C         rebinned before display.
C
      REAL F(NCOL,*), FF(MAXPLT)
C
      REAL SQRT
      INTEGER MIN0
C
      CHARACTER BLANK*80, DASH*2
      CHARACTER*2 CHAR(NCHAR)
      REAL SUM, PIXELS, FMAX, FZERO, S
      INTEGER I, ISTEP, NX, MX, IX, IY, KX, JX, JY, NY, LOW
      DATA BLANK /' '/, DASH/'--'/
      DATA CHAR / '  ', '- ', '--', '::', '==', 'll', 
     .     'II', '%%', '00', 'HH', '##' /
C
C-----------------------------------------------------------------------
C
      ISTEP=((NX-1)/MAXPLT)+1
C
C ISTEP is the number of pixels in each row of the input array which
C will have to be averaged for each pixel of the display.
C
      MX=(NX+ISTEP-1)/ISTEP
C
C MX is the number of pixels per row which will be produced on the 
C output display.
C
      LOW=MAXPLT-MX+1
      S=SQRT(AMAX1(FLOAT(NCHAR), FMAX-FZERO))
      WRITE (6,610) BLANK(1:LOW), (DASH, I=1,MX), '+ '
  610 FORMAT (A, '+', 80A2)
C
      DO 1010 IY=1,NY,ISTEP
      KX=0
      DO 1007 IX=1,NX,ISTEP
      KX=KX+1
      PIXELS=0.0
      SUM=0.0
      DO 1005 JY=IY,MIN0(NY,IY+ISTEP-1)
      DO 1003 JX=IX,MIN0(NX,IX+ISTEP-1)
      PIXELS=PIXELS+1.0
 1003 SUM=SUM+F(JX,JY)
 1005 CONTINUE
 1007 FF(KX)=SUM/PIXELS
 1010 WRITE (6,602) BLANK(1:LOW), '|', (CHAR(MIN0(
     .     NCHAR,
     .     IFIX( NCHAR*SQRT(AMAX1(0., FF(IX)-FZERO)  )/S )+1)), 
     .     IX=1,MX), '|'
  602 FORMAT ( 80A )
      WRITE (6,610) BLANK(1:LOW), (DASH, I=1,MX), '+ '
      RETURN                                             ! Normal return
C
      END!
C
C#######################################################################
C
      REAL  FUNCTION  BICUBC  (F, NBOX, DX, DY, DFDX, DFDY)
C
C Perform a type of bicubic interpolation in a grid of values.
C For a point located DX, DY distant from the corner of the grid
C (defined to be 1,1), return both the interpolated value and
C its first derivatives with respect to x and y.
C
      IMPLICIT NONE
      INTEGER NBOX
      REAL F(NBOX,NBOX), TEMP(4), DFDXT(4)
C
      REAL DX, DY, DFDX, DFDY, C1, C2, C3, C4
      INTEGER JY
C
C By construction, the point at which we want to estimate the function
C will lie between the second and third columns, and between the second
C and third rows of F, at a distance of (DX,DY) from the (2,2) element
C of F.
C
      DO JY=1,4
         C1 = 0.5*(F(3,JY)-F(1,JY))
         C4 = F(3,JY) - F(2,JY) - C1
         C2 = 3.*C4 - 0.5*(F(4,JY)-F(2,JY)) + C1
         C3 = C4 - C2
         C4 = DX*C3
         TEMP(JY) = DX*(DX*(C4+C2)+C1)+F(2,JY)
         DFDXT(JY)= DX*(C4*3.+2.*C2)+C1
      END DO
      C1 = 0.5*(TEMP(3)-TEMP(1))
      C4 = TEMP(3) - TEMP(2) - C1
      C2 = 3.*C4 - 0.5*(TEMP(4)-TEMP(2)) + C1
      C3 = C4 - C2
      C4 = DY*C3
      BICUBC = DY*(DY*(C4+C2)+C1)+TEMP(2)
      DFDY = DY*(C4*3.+2.*C2)+C1
      C1 = 0.5*(DFDXT(3)-DFDXT(1))
      C4 = DFDXT(3) - DFDXT(2) - C1
      C2 = 3.*C4 - 0.5*(DFDXT(4)-DFDXT(2)) + C1
      C3 = C4 - C2
      DFDX = DY*(DY*(DY*C3+C2)+C1)+DFDXT(2)
      RETURN
      END!
C
C#######################################################################
C
      REAL  FUNCTION  PROFIL  (IPSTYP, DX, DY, PAR, DHDXC, DHDYC, 
     .     TERM, IDERIV)
C
C Compute the value of an ANALYTIC prfile for a point DX,DY distant
C from the centroid.  Return both the computed value and its
C first derivatives with respect to x and y.  If IDERIV .NE. 0,
C return also the first derivatives with respect to all the parameters
C defining the profile.
C
      IMPLICIT NONE
      INTEGER MAXPAR, MAXPT
      PARAMETER (MAXPAR=6, MAXPT=4)
C
      REAL PAR(MAXPAR), TERM(MAXPAR)
      REAL D(MAXPT,MAXPT), W(MAXPT,MAXPT)
      REAL X(MAXPT), XSQ(MAXPT), P1XSQ(MAXPT)
C
      REAL EXP, DAOERF
C
      REAL DX, DY, DHDXC, DHDYC, WFSQ, Y, WT, WF, ONEMP3
      REAL RSQ, E, TALPHA, P1SQ, P2SQ, XY, DENOM
      REAL FUNC, YSQ, WP4FOD, P4FOD, F, P1P2, ERFX, DHDSX, ERFY
      REAL DEBY, DFBY, DBYX0, DBYY0
      REAL DHDSY, ALPHA, P2YSQ
      INTEGER I, IPSTYP, IDERIV, IX, IY, NPT
C
      DATA D / 0.00000000,  0.0,        0.0       , 0.0       ,
     .        -0.28867513,  0.28867513, 0.0       , 0.0       ,
     .        -0.38729833,  0.00000000, 0.38729833, 0.0       ,
     .        -0.43056816, -0.16999052, 0.16999052, 0.43056816/
      DATA W / 1.00000000,  0.0       , 0.0       , 0.0       ,
     .         0.50000000,  0.50000000, 0.0       , 0.0       ,
     .         0.27777778,  0.44444444, 0.27777778, 0.0       ,
     .         0.17392742,  0.32607258, 0.32607258, 0.17392742/
C
      PROFIL = 0.
      DHDXC = 0.
      DHDYC = 0.
C
      IF (IDERIV .GT. 0) THEN
         DO I=1,MAXPAR
            TERM(I) = 0.
         END DO
      END IF
C
      IF (IPSTYP .EQ. 1) THEN
C
C GAUSSIAN
C
C     F = ERFX * ERFY / (PAR(1) * PAR(2))
C
C PAR(1) is the HWHM in X; sigma(x) = 0.8493218 * HWHM
C PAR(2) is the HWHM in Y; ditto
C
         P1P2 = PAR(1)*PAR(2)
         ERFX = DAOERF(DX, 0., PAR(1), DHDXC, DHDSX)
         ERFY = DAOERF(DY, 0., PAR(2), DHDYC, DHDSY)
         PROFIL = ERFX*ERFY/P1P2
         DHDXC = DHDXC*ERFY/P1P2
         DHDYC = DHDYC*ERFX/P1P2
         IF (IDERIV .GT. 0) THEN
            TERM(1) = (DHDSX-ERFX/PAR(1))*ERFY/P1P2
            TERM(2) = (DHDSY-ERFY/PAR(2))*ERFX/P1P2
         END IF
      ELSE IF (IPSTYP .EQ. 3) THEN
C
C MOFFAT FUNCTION  BETA = 2.5
C                            BETA-1
C F = --------------------------------------------------------
C      Ax * Ay * [1 + (X/Ax)**2 + (Y/Ay)**2 + (XY*Axy)]**BETA
C
C PAR(1) is the HWHM in x at y = 0: 
C
C             1/2 = 1/[1 + (PAR(1)/Ax)**2]**BETA
C so
C             2**(1/BETA) - 1 = (PAR(1)/Ax)**2
C
C             Ax**2 = PAR(1)**2/[2**(1/BETA) - 1]
C
C When BETA = 2.5, Ax**2 = 3.129813 * PAR(1)**2
C
C Hence, let us use
C
C                                  1
C F = ---------------------------------------------------------------
C     P(1)*P(2)*{1+0.3195079*[(X/P(1))**2+(Y/P(2))**2+(XY*P(3))]**2.5
C 
C neglecting a constant of proportionality.
C
         ALPHA = 0.3195079
         TALPHA = 0.6390158                 ! 2.*ALPHA
         P1SQ = PAR(1)**2
         P2SQ = PAR(2)**2
         P1P2 = PAR(1)*PAR(2)
         XY = DX*DY
C
         DENOM = 1. + ALPHA*(DX**2/P1SQ + DY**2/P2SQ + XY*PAR(3))
         IF (DENOM .GT. 1.E4) RETURN
         FUNC = 1. / (P1P2 * DENOM**PAR(4))
         IF (FUNC .GE. 0.046) THEN
            NPT = 4
         ELSE IF (FUNC .GE. 0.0022) THEN
            NPT = 3
         ELSE IF (FUNC .GE. 0.0001) THEN
            NPT = 2
         ELSE IF (FUNC .GE. 1.E-10) THEN
            PROFIL = (PAR(4) - 1.) * FUNC
            P4FOD = PAR(4)*ALPHA*PROFIL/DENOM
            DHDXC = P4FOD*(2.*DX/P1SQ + DY*PAR(3))
            DHDYC = P4FOD*(2.*DY/P2SQ + DX*PAR(3))
            IF (IDERIV .GT. 0) THEN
               TERM(1) = (2.*P4FOD*DX**2/P1SQ-PROFIL)/PAR(1)
               TERM(2) = (2.*P4FOD*DY**2/P2SQ-PROFIL)/PAR(2)
               TERM(3) = - P4FOD*XY
C              TERM(4) = PROFIL*(1./(PAR(4)-1.)-ALOG(DENOM))
            END IF
            RETURN
         ELSE
            RETURN
         END IF
C
         DO IX=1,NPT
            X(IX) = DX+D(IX,NPT)
            XSQ(IX) = X(IX)**2
            P1XSQ(IX) = XSQ(IX)/P1SQ
         END DO
C
         DO IY=1,NPT
            Y = DY+D(IY,NPT)
            YSQ = Y**2
            P2YSQ = YSQ/P2SQ
            DO IX=1,NPT
               WT = W(IY,NPT)*W(IX,NPT)
               XY=X(IX)*Y
               DENOM = 1. + ALPHA*(P1XSQ(IX) + P2YSQ + XY*PAR(3))
               FUNC = (PAR(4) - 1.) / (P1P2 * DENOM**PAR(4))
               P4FOD = PAR(4)*ALPHA*FUNC/DENOM
               WP4FOD = WT*P4FOD
               WF = WT*FUNC
               PROFIL = PROFIL + WF
               DHDXC = DHDXC + WP4FOD*(2.*X(IX)/P1SQ + Y*PAR(3))
               DHDYC = DHDYC + WP4FOD*(2.*Y/P2SQ + X(IX)*PAR(3))
               IF (IDERIV .GT. 0) THEN
                  TERM(1) = TERM(1) + 
     .                     (2.*WP4FOD*P1XSQ(IX)-WF)/PAR(1)
                  TERM(2) = TERM(2) + 
     .                     (2.*WP4FOD*P2YSQ-WF)/PAR(2)
                  TERM(3) = TERM(3) - WP4FOD*XY
C                 TERM(4) = TERM(4) + WF*(1./(PAR(4)-1.)-ALOG(DENOM))
               END IF
            END DO
         END DO
      ELSE IF (IPSTYP .EQ. 5) THEN
C
C Penny function --- Gaussian core plus Lorentzian wings.  The Lorentzian 
C is elongated along the x or y axis, the Gaussian may be tilted.
C
         P1SQ = PAR(1)**2
         P2SQ = PAR(2)**2
         ONEMP3 = 1.-PAR(3)
         XY = DX*DY
C
         RSQ = DX**2/P1SQ + DY**2/P2SQ
         IF (RSQ .GT. 1.E10) RETURN
C
         F = 1./(1.+RSQ)
         RSQ = RSQ + XY*PAR(4)
         IF (RSQ .LT. 34.) THEN
            E = EXP(-0.6931472*RSQ)
            FUNC = PAR(3)*E + ONEMP3*F
         ELSE
            E = 0.
            FUNC = ONEMP3*F
         END IF
C
         IF (FUNC .GE. 0.046) THEN
            NPT = 4
         ELSE IF (FUNC .GE. 0.0022) THEN
            NPT = 3
         ELSE IF (FUNC .GE. 0.0001) THEN
            NPT = 2
         ELSE IF (FUNC .GE. 1.E-10) THEN
            PROFIL = FUNC
            DFBY = ONEMP3*F**2
            DEBY = 0.6931472*PAR(3)*E
            DBYX0 = 2.*DX/P1SQ
            DBYY0 = 2.*DY/P2SQ
            DHDXC = DEBY*(DBYX0 + DY*PAR(4)) + DFBY*DBYX0
            DHDYC = DEBY*(DBYY0 + DX*PAR(4)) + DFBY*DBYY0
            IF (IDERIV .GT. 0) THEN
               DBYX0 = DBYX0*DX/PAR(1)
               DBYY0 = DBYY0*DY/PAR(2)
               DFBY = DFBY + DEBY
               TERM(1) = DFBY * DBYX0
               TERM(2) = DFBY * DBYY0
               TERM(3) = E - F
               TERM(4) = - DEBY * XY
     .              / (0.5 - ABS(PAR(4)))
            END IF
            RETURN
         ELSE
            RETURN
         END IF
C
         DO IX=1,NPT
            X(IX) = DX+D(IX,NPT)
            P1XSQ(IX) = X(IX)/P1SQ
         END DO
C
         DO IY=1,NPT
            Y = DY+D(IY,NPT)
            P2YSQ = Y/P2SQ
            DO IX=1,NPT
               WT = W(IY,NPT)*W(IX,NPT)
               XY = X(IX)*Y
               RSQ = P1XSQ(IX)*X(IX) + P2YSQ*Y
               F = 1./(1.+RSQ)
               RSQ = RSQ + XY*PAR(4)
               IF (RSQ .LT. 34.) THEN
                  E = EXP(-0.6931472*RSQ)
                  FUNC = PAR(3)*E + ONEMP3*F
                  DEBY = 0.6931472*WT*PAR(3)*E
               ELSE
                  E = 0.
                  FUNC = ONEMP3*F
                  DEBY = 0.
               END IF
               PROFIL = PROFIL + WT*FUNC
               DFBY = WT*ONEMP3*F**2
               DBYX0 = 2.*P1XSQ(IX)
               DBYY0 = 2.*P2YSQ
               DHDXC = DHDXC + 
     .              DEBY*(DBYX0 + DY*PAR(4)) + DFBY*DBYX0
               DHDYC = DHDYC +
     .              DEBY*(DBYY0 + DX*PAR(4)) + DFBY*DBYY0
               IF (IDERIV .GT. 0) THEN
                  DBYX0 = DBYX0*DX/PAR(1)
                  DBYY0 = DBYY0*DY/PAR(2)
                  TERM(1) = TERM(1) + (DFBY+DEBY)*DBYX0
                  TERM(2) = TERM(2) + (DFBY+DEBY)*DBYY0
                  TERM(3) = TERM(3) + WT*(E-F)
                  TERM(4) = TERM(4) - DEBY * XY
               END IF
            END DO
         END DO
      ELSE IF (IPSTYP .EQ. 6) THEN
C
C Penny function --- Gaussian core plus Lorentzian wings.
C The Lorentzian and Gaussian may be tilted in different
C directions.
C
         P1SQ = PAR(1)**2
         P2SQ = PAR(2)**2
         ONEMP3 = 1.-PAR(3)
         XY = DX*DY
C
         RSQ = DX**2/P1SQ + DY**2/P2SQ
         DFBY = RSQ + PAR(5)*XY
         IF (DFBY .GT. 1.E10) RETURN
         F = 1./(1.+DFBY)
C
         DEBY = RSQ + PAR(4)*XY
         IF (DEBY .LT. 34.) THEN
            E = EXP(-0.6931472*DEBY)
         ELSE
            E = 0.
         END IF
C
         FUNC = PAR(3)*E + ONEMP3*F
         IF (FUNC .GE. 0.046) THEN
            NPT = 4
         ELSE IF (FUNC .GE. 0.0022) THEN
            NPT = 3
         ELSE IF (FUNC .GE. 0.0001) THEN
            NPT = 2
         ELSE IF (FUNC .GE. 1.E-10) THEN
            PROFIL = FUNC
            DFBY = ONEMP3*F**2
            DEBY = 0.6931472*PAR(3)*E
            DBYX0 = 2.*DX/P1SQ
            DBYY0 = 2.*DY/P2SQ
            DHDXC = DEBY*(DBYX0 + DY*PAR(4)) + 
     .              DFBY*(DBYX0 + DY*PAR(5))
            DHDYC = DEBY*(DBYY0 + DX*PAR(4)) + 
     .              DFBY*(DBYY0 + DX*PAR(5))
            IF (IDERIV .GT. 0) THEN
               DBYX0 = DBYX0*DX/PAR(1)
               DBYY0 = DBYY0*DY/PAR(2)
               TERM(5) = -DFBY * XY
               DFBY = DFBY + DEBY
               TERM(1) = DFBY * DBYX0
               TERM(2) = DFBY * DBYY0
               TERM(3) = E - F
               TERM(4) = - DEBY * XY
            END IF
            RETURN
         ELSE
            RETURN
         END IF
C
         DO IX=1,NPT
            X(IX) = DX+D(IX,NPT)
            P1XSQ(IX) = X(IX)/P1SQ
         END DO
C
         DO IY=1,NPT
            Y = DY+D(IY,NPT)
            P2YSQ = Y/P2SQ
            DO IX=1,NPT
               WT = W(IY,NPT)*W(IX,NPT)
               XY = X(IX)*Y
               RSQ = P1XSQ(IX)*X(IX) + P2YSQ*Y
               F = 1./(1.+RSQ+PAR(5)*XY)
               DEBY = RSQ + PAR(4)*XY
               IF (DEBY .LT. 34.) THEN
                  E = EXP(-0.6931472*DEBY)
                  FUNC = PAR(3)*E + ONEMP3*F
                  DEBY = 0.6931472*WT*PAR(3)*E
               ELSE
                  E = 0.
                  FUNC = ONEMP3*F
                  DEBY = 0.
               END IF
               PROFIL = PROFIL + WT*FUNC
               DFBY = WT*ONEMP3*F**2
               DBYX0 = 2.*P1XSQ(IX)
               DBYY0 = 2.*P2YSQ
               DHDXC = DHDXC + 
     .              DEBY*(DBYX0 + DY*PAR(4)) + 
     .              DFBY*(DBYX0 + DY*PAR(5))
               DHDYC = DHDYC +
     .              DEBY*(DBYY0 + DX*PAR(4)) + 
     .              DFBY*(DBYY0 + DX*PAR(5))
               IF (IDERIV .GT. 0) THEN
                  DBYX0 = DBYX0*DX/PAR(1)
                  DBYY0 = DBYY0*DY/PAR(2)
                  TERM(1) = TERM(1) + (DFBY+DEBY)*DBYX0
                  TERM(2) = TERM(2) + (DFBY+DEBY)*DBYY0
                  TERM(3) = TERM(3) + WT*(E-F)
                  TERM(4) = TERM(4) - DEBY * XY
                  TERM(5) = TERM(5) - DFBY * XY
               END IF
            END DO
         END DO
      ELSE IF (IPSTYP .EQ. 2) THEN
C
C MOFFAT FUNCTION   BETA = 1.5
C
C                            BETA-1
C F = --------------------------------------------------------
C      Ax * Ay * [1 + (X/Ax)**2 + (Y/Ay)**2 + (XY*Axy)]**BETA
C
C PAR(1) is the HWHM in x at y = 0: 
C
C             1/2 = 1/[1 + (PAR(1)/Ax)**2]**BETA
C so
C             2**(1/BETA) - 1 = (PAR(1)/Ax)**2
C
C             Ax**2 = PAR(1)**2/[2**(1/BETA) - 1]
C
C When BETA = 1.5, Ax**2 = 1.7024144 * PAR(1)**2
C
C Hence, let us use
C
C                                  1
C F = ---------------------------------------------------------------
C     P(1)*P(2)*{1+0.5874011*[(X/P(1))**2+(Y/P(2))**2+(XY*P(3))]**1.5
C 
C neglecting a constant of proportionality.
C
         ALPHA = 0.5874011
         TALPHA = 1.1748021  ! 2.*ALPHA
         P1SQ = PAR(1)**2
         P2SQ = PAR(2)**2
         P1P2 = PAR(1)*PAR(2)
         XY = DX*DY
C
         DENOM = 1. + ALPHA*(DX**2/P1SQ + DY**2/P2SQ + XY*PAR(3))
         IF (DENOM .GT. 5.E6) RETURN
         FUNC = 1. / (P1P2 * DENOM**PAR(4))
         IF (FUNC .GE. 0.046) THEN
            NPT = 4
         ELSE IF (FUNC .GE. 0.0022) THEN
            NPT = 3
         ELSE IF (FUNC .GE. 0.0001) THEN
            NPT = 2
         ELSE IF (FUNC .GE. 1.E-10) THEN
            PROFIL = (PAR(4) - 1.) * FUNC
            P4FOD = PAR(4)*ALPHA*PROFIL/DENOM
            DHDXC = P4FOD*(2.*DX/P1SQ + DY*PAR(3))
            DHDYC = P4FOD*(2.*DY/P2SQ + DX*PAR(3))
            IF (IDERIV .GT. 0) THEN
               TERM(1) = (2.*P4FOD*DX**2/P1SQ-PROFIL)/PAR(1)
               TERM(2) = (2.*P4FOD*DY**2/P2SQ-PROFIL)/PAR(2)
               TERM(3) = - P4FOD*XY
C              TERM(4) = PROFIL*(1./(PAR(4)-1.)-ALOG(DENOM))
            END IF
            RETURN
         ELSE
            RETURN
         END IF
C
         DO IX=1,NPT
            X(IX) = DX+D(IX,NPT)
            XSQ(IX) = X(IX)**2
            P1XSQ(IX) = XSQ(IX)/P1SQ
         END DO
C
         DO IY=1,NPT
            Y = DY+D(IY,NPT)
            YSQ = Y**2
            P2YSQ = YSQ/P2SQ
            DO IX=1,NPT
               WT = W(IY,NPT)*W(IX,NPT)
               XY=X(IX)*Y
               DENOM = 1. + ALPHA*(P1XSQ(IX) + P2YSQ + XY*PAR(3))
               FUNC = (PAR(4) - 1.)/ (P1P2 * DENOM**PAR(4))
               P4FOD = PAR(4)*ALPHA*FUNC/DENOM
               WP4FOD = WT*P4FOD
               WF = WT*FUNC
               PROFIL = PROFIL + WF
               DHDXC = DHDXC + WP4FOD*(2.*X(IX)/P1SQ + Y*PAR(3))
               DHDYC = DHDYC + WP4FOD*(2.*Y/P2SQ + X(IX)*PAR(3))
               IF (IDERIV .GT. 0) THEN
                  TERM(1) = TERM(1) + 
     .                     (2.*WP4FOD*P1XSQ(IX)-WF)/PAR(1)
                  TERM(2) = TERM(2) + 
     .                     (2.*WP4FOD*P2YSQ-WF)/PAR(2)
                  TERM(3) = TERM(3) - WP4FOD*XY
C                 TERM(4) = TERM(4) + WF*(1./(PAR(4)-1.)-ALOG(DENOM))
               END IF
            END DO
         END DO
      ELSE IF (IPSTYP .EQ. 4) THEN
C
C LORENTZ FUNCTION
C                      1
C F = --------------------------------------
C     [1 + (X/Ax)**2 + (Y/Ay)**2 + (XY*Axy)]
C
C PAR(1) is the HWHM in x at y = 0.
C
         P1SQ = PAR(1)**2
         P2SQ = PAR(2)**2
         P1P2 = PAR(1)*PAR(2)
         XY = DX*DY
C
         DENOM = 1. + DX**2/P1SQ + DY**2/P2SQ + XY*PAR(3)
         IF (DENOM .GT. 1.E10) RETURN
         FUNC = 1. / DENOM
         IF (FUNC .GE. 0.046) THEN
            NPT = 4
         ELSE IF (FUNC .GE. 0.0022) THEN
            NPT = 3
         ELSE IF (FUNC .GE. 0.0001) THEN
            NPT = 2
         ELSE IF (FUNC .GE. 1.E-10) THEN
            PROFIL = FUNC
            WFSQ = FUNC**2
            DHDXC = WFSQ*(2.*DX/P1SQ + DY*PAR(3))
            DHDYC = WFSQ*(2.*DY/P2SQ + DX*PAR(3))
            IF (IDERIV .GT. 0) THEN
               TERM(1) = WFSQ*(2.*DX**2/P1SQ)/PAR(1)
               TERM(2) = WFSQ*(2.*DY**2/P2SQ)/PAR(2)
               TERM(3) = - WFSQ*XY
            END IF
            RETURN
         ELSE
            RETURN
         END IF
C
         DO IX=1,NPT
            X(IX) = DX+D(IX,NPT)
            XSQ(IX) = X(IX)**2
            P1XSQ(IX) = XSQ(IX)/P1SQ
         END DO
C
         DO IY=1,NPT
            Y = DY+D(IY,NPT)
            YSQ = Y**2
            P2YSQ = YSQ/P2SQ
            DO IX=1,NPT
               WT = W(IY,NPT)*W(IX,NPT)
               XY=X(IX)*Y
               DENOM = 1. + P1XSQ(IX) + P2YSQ + XY*PAR(3)
               FUNC = 1. / DENOM
               WF = WT*FUNC
               WFSQ = WF*FUNC
               PROFIL = PROFIL + WF
               DHDXC = DHDXC + WFSQ*(2.*X(IX)/P1SQ + Y*PAR(3))
               DHDYC = DHDYC + WFSQ*(2.*Y/P2SQ + X(IX)*PAR(3))
               IF (IDERIV .GT. 0) THEN
                  TERM(1) = TERM(1) + WFSQ*(2.*P1XSQ(IX))/PAR(1)
                  TERM(2) = TERM(2) + WFSQ*(2.*P2YSQ)/PAR(2)
                  TERM(3) = TERM(3) - WFSQ*XY
               END IF
            END DO
         END DO
      ELSE
         write(*,*)'Invalid PSF type.'
	 write(*,*)'Bye'
	 stop
      END IF
      RETURN
      END!
C
C=======================================================================
C
      FUNCTION NPARAM  (IPSTYP, FWHM, LABEL, PAR, MAXPAR)
      CHARACTER*8 LABEL
      INTEGER MAXPAR
      REAL PAR(MAXPAR)
      PAR(1) = FWHM/2.
      PAR(2) = PAR(1)
      IF (IPSTYP .EQ. 1) THEN
         NPARAM = 2
         LABEL = 'GAUSSIAN'
      ELSE IF (IPSTYP .EQ. 3) THEN
         NPARAM = 3
         PAR(3) = 0.
         PAR(4) = 2.5
         LABEL = 'MOFFAT25'
      ELSE IF (IPSTYP .EQ. 5) THEN
         NPARAM = 4
         PAR(3) = 0.75
         PAR(4) = 0.0
         LABEL = 'PENNY1  '
      ELSE IF (IPSTYP .EQ. 6) THEN
         NPARAM = 5
         PAR(3) = 0.75
         PAR(4) = 0.0
         PAR(5) = 0.0
         LABEL = 'PENNY2  '
      ELSE IF (IPSTYP .EQ. 2) THEN
         NPARAM = 3
         PAR(3) = 0.
         PAR(4) = 1.5
         LABEL = 'MOFFAT15'
      ELSE IF (IPSTYP .EQ. 4) THEN
         NPARAM = 3
         PAR(3) = 0.
         LABEL = 'LORENTZ '
      ELSE
         write(*,*)'Invalid PSF type: '//CHAR(IPSTYP+48)
      END IF
      RETURN
      END!
C
C#######################################################################
C
      REAL  FUNCTION  DAOERF (XIN, XO, BETA, DFDXO, DFDBET)
C
C Numerically integrate a Gaussian function 
C
C          F = EXP {-0.5*[(x-XO)/SIGMA]**2 },
C
C from XIN-0.5 to XIN+0.5 using Gauss-Legendre integration.  BETA
C is the half-width at half-maximum, which is equal to 1.17741 * SIGMA.
C Thus,
C
C          F = EXP {-0.6931472*[(x-XO)/BETA]**2 }.
C
C Also: provide the first derivative of the integral with respect to 
C Xo and BETA.  Use Gauss-Legendre integration.
C
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
      INTEGER MAXPT
      PARAMETER (MAXPT=4)
C
      REAL DX(MAXPT,MAXPT), WT(MAXPT,MAXPT)
C
      REAL EXP
C
      REAL X, XSQ
      REAL XIN, XO, BETA, DFDXO, DFDBET, BETASQ, DELTAX, F, WF
      INTEGER NPT, I
C
      DATA DX / 0.00000000,  0.0,        0.0       , 0.0       ,
     .         -0.28867513,  0.28867513, 0.0       , 0.0       ,
     .         -0.38729833,  0.00000000, 0.38729833, 0.0       ,
     .         -0.43056816, -0.16999052, 0.16999052, 0.43056816/
      DATA WT / 1.00000000,  0.0       , 0.0       , 0.0       ,
     .          0.50000000,  0.50000000, 0.0       , 0.0       ,
     .          0.27777778,  0.44444444, 0.27777778, 0.0       ,
     .          0.17392742,  0.32607258, 0.32607258, 0.17392742/
      DAOERF = 0.
      DFDXO = 0.
      DFDBET = 0.
      BETASQ=BETA**2
      DELTAX = XIN-XO
C
      XSQ = DELTAX**2
      F = XSQ/BETASQ
      IF (F .GT. 34.) RETURN
      F = EXP(-0.6931472*F)
      IF (F .GE. 0.046) THEN
         NPT = 4
      ELSE IF (F .GE. 0.0022) THEN
         NPT = 3
      ELSE IF (F .GE. 0.0001) THEN
         NPT = 2
      ELSE IF (F .GE. 1.E-10) THEN
         DAOERF = F
         DFDXO = 1.3862944 * DELTAX * F / BETASQ
         DFDBET = 1.3862944 * XSQ * F / (BETASQ*BETA)
         RETURN
      ELSE
         RETURN
      END IF
C
      DO I=1,NPT
         X = DELTAX + DX(I,NPT)
         XSQ = X**2
         F = EXP(-0.6931472*XSQ/BETASQ)
         WF = WT(I,NPT)*F
         DAOERF = DAOERF+WF
         DFDXO = DFDXO + X*WF
         DFDBET = DFDBET + XSQ*WF
      END DO
      DFDXO = 1.3862944*DFDXO/BETASQ
      DFDBET = 1.3862944*DFDBET/(BETASQ*BETA)
C
      RETURN
      END!
C
C#######################################################################
C
      SUBROUTINE  PKFIT (F, NX, NY, MAXBOX, X, Y, SCALE, SKY, RADIUS,
     .     LOBAD, HIBAD, BRIGHT, IPSTYP, PAR, MAXPAR, NPAR,
     .     PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, 
     .     DELTAX, DELTAY, ERRMAG, CHI, SHARP, NITER, NCOL, NROW)
C
C=======================================================================
C
C This is the subroutine which does the actual one-star least-squares
C profile fit for PEAK. 
C
C           OFFICIAL DAO VERSION:  1991 April 6
C
C Arguments
C      F (INPUT) is an NX by NY array containing actual picture data.
C
C MAXBOX (INPUT) is the maximum value allowable for either NX or NY,
C        needed for the dimension statements below.  PEAK and PSF will
C        provide different values of MAXBOX.
C
C  SCALE (INPUT/OUTPUT) is the initial estimate of the brightness of
C        the star, expressed as a fraction of the brightness of the
C        PSF.  Upon return, the final computed value of SCALE will
C        be passed back to the calling routine.
C
C   X, Y (INPUT/OUTPUT) are the initial estimates of the centroid of 
C        the star relative to the corner (1,1) of the subarray.  Upon
C        return, the final computed values of X and Y will be passed 
C        back to the calling routine.
C
C    SKY (INPUT) is the local sky brightness value, carried on from
C        PHOTOMETRY via the data files.
C
C RADIUS (INPUT) is the fitting radius-- only pixels within RADIUS of
C        the instantaneous estimate of the star's centroid will be
C        included in the fit.
C
C LOBAD and HIBAD (INPUT) are bad pixel limits-- any pixel whose 
C        brightness value falls outside this range will be presumed to 
C        be bad, and will be ignored.
C
C BRIGHT (INPUT) contains the brightness normalization of the model
C        analytic point spread function
C
C    PAR (INPUT) contains the values of the remaining NPARAM parameters 
C        defining the analytic function which approximates the core of 
C        the PSF.
C
C    PSF (INPUT) is an NPSF by NPSF by NEXP+NFRAC look-up table 
C        containing corrections from the analytic approximation of the 
C        PSF to the true PSF.
C
C ERRMAG (OUTPUT) is the estimated standard error of the value of SCALE
C        returned by this routine.
C
C    CHI (OUTPUT) is the estimated goodness-of-fit statistic:  the ratio
C        of the observed pixel-to-pixel mean absolute deviation from
C        the profile fit, to the value expected on the basis of the
C        read-out noise and the photons/ADU (which are brought in
C        through COMMON block /ERROR/).
C
C  SHARP (OUTPUT) is a goodness-of-fit statistic describing how much
C        broader the actual profile of the object appears than the
C        profile of the PSF.
C
C  NITER (OUTPUT) is the number of iterations the solution required to
C        achieve convergence.  If NITER = 50, the solution did not 
C        converge.  If for some reason a singular matrix occurs during
C        the least-squares solution, this will be flagged by setting
C        NITER = -1.
C
C=======================================================================
C        
      IMPLICIT NONE
      INTEGER MAXPSF, MAXEXP, MAXBOX, MAXPAR
C
C Parameter
C
C MAXPSF is the length of the side of the largest PSF look-up table
C        allowed.  (Note:  half-pixel grid spacing)
C
      REAL C(3,3), V(3), CLAMP(3), DTOLD(3)
      REAL F(MAXBOX,MAXBOX), T(3), DT(3), NUMER
      REAL PSF(MAXPSF,MAXPSF,MAXEXP), PAR(MAXPAR)
      REAL USEPSF, AMAX1, AMIN1, ABS
      INTEGER MIN0, MAX0
C
      REAL LOBAD, DF, WT, DFDSIG, RHOSQ, RELERR, SIG, SIGSQ
      REAL FPOS, DVDXC, DVDYC, RSQ, DX, DXSQ, DY, DYSQ, DATUM
      REAL DENOM, SUMWT, CHIOLD, RADSQ, SHARP, CHI, ERRMAG
      REAL DELTAX, DELTAY, X, Y, PHPADU, RONOIS, PERR, PKERR
      REAL SCALE, SKY, RADIUS, HIBAD, BRIGHT
      INTEGER I, J, IX, IY, IXLO, IXHI, IYLO, IYHI, ISTAT
      INTEGER NPIX, NFRAC, NEXP, NITER, NPAR, NPSF, IPSTYP
      INTEGER NX, NY, NCOL, NROW
      LOGICAL CLIP, REDO
C
      COMMON /ERROR/ PHPADU, RONOIS, PERR, PKERR
C
C-----------------------------------------------------------------------
C
C Initialize a few things for the solution.
C
      RADSQ=RADIUS**2
      DO 1010 I=1,3
      CLAMP(I)=1.
 1010 DTOLD(I)=0.0
      CHIOLD=1.
      NITER=0
      SHARP=0.
      CLIP = .FALSE.
C
C-----------------------------------------------------------------------
C
C Here begins the big least-squares loop.
C
 2000 NITER=NITER+1
C
C Initialize things for this iteration.  CHI and SHARP will be 
C goodness-of-fit indices.  CHI will also be used in determining the 
C weights of the individual pixels.  As the solution iterates, the new 
C value of CHI will be built up in the variable CHI, while a smoothed
C value of CHI computed from the previous iteration will be carried 
C along in CHIOLD.
C
      CHI=0.0
      SUMWT=0.0
      NUMER=0.0
      DENOM=0.0
      DO 2010 I=1,3
      V(I)=0.0                          ! Zero the vector of residuals
      DO 2010 J=1,3
 2010 C(I,J)=0.0                        ! Zero the normal matrix
C
C Choose the little box containing points inside the fitting radius.
C
      IXLO=MAX0(1, INT(X-RADIUS))
      IYLO=MAX0(1, INT(Y-RADIUS))
      IXHI=MIN0(NX, INT(X+RADIUS)+1)
      IYHI=MIN0(NY, INT(Y+RADIUS)+1)
C
C Now build up the normal matrix and vector of residuals.
C
      NPIX=0
      DO 2090 IY=IYLO,IYHI
      DY=FLOAT(IY)-Y
      DYSQ=DY**2
      DO 2090 IX=IXLO,IXHI
      DATUM=F(IX,IY)
      IF ((DATUM .LT. LOBAD) .OR. (DATUM .GT. HIBAD)) GO TO 2090
      DX=FLOAT(IX)-X
      DXSQ=DX**2
C
C DX and DY are the distance of this pixel from the centroid of the 
C star.  Is this pixel truly inside the fitting radius?
C
      RSQ=(DXSQ+DYSQ)/RADSQ
      IF (1.-RSQ .LE. 2.E-6) GO TO 2090   ! Prevents floating underflows
C
C The fitting equation is of the form
C 
C Observed brightness = 
C     [SCALE + delta(SCALE)] * [PSF + delta(Xcen)*d(PSF)/d(Xcen) +
C                                           delta(Ycen)*d(PSF)/d(Ycen) ]
C
C and is solved for the unknowns delta(SCALE) ( = the correction to 
C the brightness ratio between the program star and the PSF) and 
C delta(Xcen) and delta(Ycen) ( = corrections to the program star's 
C centroid).
C
C The point-spread function is equal to the sum of the integral under
C a two-dimensional analytic profile plus values interpolated from
C a look-up table.
C
      T(1)=USEPSF(IPSTYP, DX, DY, BRIGHT, PAR, PSF, NPSF, NPAR, NEXP, 
     .     NFRAC, DELTAX, DELTAY, DVDXC, DVDYC)
      IF ((SCALE*T(1)+SKY .GT. HIBAD) .AND. (NITER.GE.3)) GO TO 2090
      T(2)=SCALE*DVDXC
      T(3)=SCALE*DVDYC
      DF=F(IX,IY)-SKY-SCALE*T(1)
C
C DF is the residual of the brightness in this pixel from the PSF fit.
C
C The expected random error in the pixel is the quadratic sum of
C the Poisson statistics, plus the readout noise, plus an estimated
C error of PERERR% of the total brightness for the difficulty of flat-
C fielding and bias-correcting the chip, plus an estimated error of 
C some fraction of the fourth derivative at the peak of the profile,
C to account for the difficulty of accurately interpolating within the 
C point-spread function.  The fourth derivative of the PSF is roughly
C proportional to H/sigma**4 (sigma is the width parameter for
C the stellar core); using the geometric mean of sigma(x) and sigma(y), 
C this becomes H/[sigma(x)*sigma(y)]**2.  The ratio of the fitting 
C error to this quantity is estimated from a good-seeing CTIO frame to 
C be approximately 0.027.  NOWADAYS, PAR(1) and PAR(2) are the HWHM,
C not the sigma, so the coefficient is now of order 0.05 = 5.0%.
C
      FPOS=AMAX1(0., F(IX,IY)-DF)
C
C FPOS = raw data minus residual = model-predicted value of the 
C intensity at this point (which presumably is non-negative).
C
      SIGSQ=FPOS/PHPADU+RONOIS+(PERR*FPOS)**2+(PKERR*(FPOS-SKY))**2
      SIG=SQRT(SIGSQ)
      RELERR=ABS(DF/SIG)
C
C SIG is the anticipated standard error of the intensity in this pixel,
C including readout noise, Poisson photon statistics, and an estimate
C of the standard error of interpolating within the PSF.  
C
      WT=5./(5.+RSQ/(1.-RSQ))             ! Weight as function of radius
C
C Now add this pixel into the quantities which go to make up the SHARP
C index.
C
      RHOSQ=DXSQ/PAR(1)**2+DYSQ/PAR(2)**2
C
C Include in the sharpness index only those pixels within six
C HWHMs of the centroid of the star.  (This saves time and
C floating underflows by excluding pixels which contribute less than
C about one part in a million to the index.)
C
      IF (RHOSQ .LE. 36.) THEN
         RHOSQ=0.6931472*RHOSQ
         DFDSIG=EXP(-RHOSQ)*(RHOSQ-1.)
         FPOS=AMAX1(0., F(IX,IY)-SKY)+SKY
         NUMER=NUMER+DFDSIG*DF/SIGSQ
         DENOM=DENOM+DFDSIG**2/SIGSQ
      END IF
C
C Derive the weight of this pixel.  First of all, the weight depends
C upon the distance of the pixel from the centroid of the star-- it
C is determined from a function which is very nearly unity for radii
C much smaller than the fitting radius, and which goes to zero for
C radii very near the fitting radius.  Then reject any pixels with 
C 10-sigma residuals (after the first iteration).
C
      CHI=CHI+WT*RELERR
      SUMWT=SUMWT+WT
C
C Now the weight is scaled to the inverse square of the expected mean 
C error.
C
      WT=WT/SIGSQ
C
C Reduce the weight of a bad pixel.  A pixel having a residual of 2.5
C sigma gets reduced to half weight; a pixel having a residual of 5.
C sigma gets weight 1/257.
C
      IF (CLIP) WT=WT/(1.+(0.4*RELERR/CHIOLD)**8)
C
C Now add the pixel into the vector of residuals and the normal matrix.
C
      DO 2030 I=1,3
      V(I)=V(I)+DF*T(I)*WT
      DO 2030 J=1,3
 2030 C(I,J)=C(I,J)+T(I)*T(J)*WT
C
      NPIX=NPIX+1
 2090 CONTINUE                                 ! End of loop over pixels
C
C Compute the (robust) goodness-of-fit index CHI.
C
      IF (SUMWT .GT. 3.) THEN
         CHI=1.2533141*CHI/SQRT(SUMWT*(SUMWT-3.))
C
C CHI is pulled toward its expected value of unity before being stored
C in CHIOLD to keep the statistics of a small number of pixels from 
C completely dominating the error analysis.  
C
         CHIOLD=((SUMWT-3.)*CHI+3.)/SUMWT
      ELSE
         CHI = 1.
         CHIOLD = 1.
      END IF
C
C Compute the parameter corrections and check for convergence.
C
      IF (NPIX .LT. 3) THEN
         NITER=-1
         RETURN
      END IF
      CALL INVERS (C, 3, 3, ISTAT)
      IF (ISTAT .EQ. 0) GO TO 2100
      NITER=-1
      RETURN
C
C Everything OK so far.
C
 2100 CALL VMUL (C, 3, 3, V, DT)
C 
C In the beginning, the brightness of the star will not be permitted
C to change by more than two magnitudes per iteration (that is to say, 
C if the estimate is getting brighter, it may not get brighter by
C more than 525% per iteration, and if it is getting fainter, it may
C not get fainter by more than 84% per iteration).  The x and y
C coordinates of the centroid will be allowed to change by no more
C than one-half pixel per iteration.  Any time that a parameter
C correction changes sign, the maximum permissible change in that
C parameter will be reduced by a factor of 2.
C
      DO 2110 I=1,3
      IF (DTOLD(I)*DT(I) .LT. 0.) THEN
         CLAMP(I)=0.5*CLAMP(I)
      ELSE
         CLAMP(I)=MIN(1., 1.1*CLAMP(I))
      END IF
 2110 DTOLD(I)=DT(I)
C
      SCALE=SCALE+DT(1)/
     .  (1.+AMAX1(DT(1)/(5.25*SCALE),-DT(1)/(0.84*SCALE))/CLAMP(1))
      X=X+DT(2)/(1.+ABS(DT(2))/(0.4*CLAMP(2)))
      Y=Y+DT(3)/(1.+ABS(DT(3))/(0.4*CLAMP(3)))
      IF (NITER .LE. 1) GO TO 2000
      REDO=.FALSE.
C
C Convergence criteria:  if the most recent computed correction to the 
C brightness is larger than 0.01% or than 0.05 * sigma(brightness),
C whichever is larger, OR if the absolute change in X or Y is
C greater than 0.001 pixels, convergence has not been achieved.
C
      ERRMAG=CHIOLD*SQRT(C(1,1))
      IF (CLIP) THEN
         IF (ABS(DT(1)) .GT. 
     .        AMAX1( 0.05*ERRMAG, 0.0001*SCALE )) THEN
            REDO=.TRUE.
         ELSE IF (AMAX1(ABS(DT(2)),ABS(DT(3))) .GT. 0.001) THEN
            REDO=.TRUE.
         END IF
      ELSE
         IF (ABS(DT(1)) .GT. 
     .        AMAX1( ERRMAG, 0.002*SCALE )) THEN
            REDO = .TRUE.
         ELSE IF (AMAX1(ABS(DT(2)),ABS(DT(3))) .GT. 0.02) THEN
            REDO = .TRUE.
         END IF
      END IF
C
      IF (NITER .GE. 50) REDO = .FALSE.
      IF (REDO) GO TO 2000
      IF ((NITER .LT. 50) .AND. (.NOT. CLIP)) THEN
         CLIP = .TRUE.
         DO I=1,3
            DTOLD(I) = 0.
            CLAMP(I) = AMAX1(CLAMP(I), 0.25)
         END DO
         GO TO 2000
      ELSE
         SHARP=1.4427*PAR(1)*PAR(2)*NUMER/(BRIGHT*SCALE*DENOM)
         SHARP=AMIN1(99.999,AMAX1(SHARP,-99.999))
         RETURN
      END IF
C
      END!
c
      subroutine ovrwrt (line, iwhich)
      character*(*) line
      if (iwhich .eq. 1) then
         write (6,1) line
    1    format (/a)
      else if (iwhich .eq. 2) then
         write (6,2) line, char(13)
    2    format (a, a1, $)
      else if (iwhich .eq. 3) then
         write (6,3) line
    3    format (a)
      else
         write (6,4) line, char(13)
    4    format (/a, a1, $)
      end if
      call flush (6)
      return
      end
C
C#######################################################################
C
      REAL  FUNCTION  USEPSF (IPSTYP, DX, DY, BRIGHT, PAR, PSF, 
     .     NPSF, NPAR, NEXP, NFRAC, DELTAX, DELTAY, DVDXC, DVDYC)
C
C Evaluate the PSF for a point distant DX, DY from the center of a
C star located at relative frame coordinates DELTAX, DELTAY.
C
      IMPLICIT NONE
      INTEGER MAXPSF, MAXPAR, MAXEXP
      PARAMETER (MAXPSF=207, MAXPAR=6, MAXEXP=10)
C
      REAL PAR(*), PSF(MAXPSF,MAXPSF,*), JUNK(MAXEXP)
C
      REAL PROFIL, BICUBC
C
      REAL MIDDLE, BRIGHT, DVDXC, DVDYC, DELTAX, DELTAY, XX, YY
      REAL DX, DY, CORR, DFDX, DFDY
      INTEGER K, LX, LY, IPSTYP
      INTEGER NFRAC, NTERM, NPSF, NEXP, NPAR
C
      NTERM = NEXP + NFRAC
      USEPSF = BRIGHT*PROFIL(IPSTYP, DX, DY, PAR, DVDXC, DVDYC, 
     .     JUNK, 0)
      DVDXC = BRIGHT*DVDXC
      DVDYC = BRIGHT*DVDYC
      IF (NTERM .LT. 0) RETURN
      MIDDLE = (NPSF+1)/2
C
C The PSF look-up tables are centered at (MIDDLE, MIDDLE).
C
      IF (NEXP .GE. 0) THEN
         JUNK(1) = 1.
         IF (NEXP .GE. 2) THEN
            JUNK(2) = DELTAX
            JUNK(3) = DELTAY
            IF (NEXP .GE. 4) THEN
               JUNK(4) = 1.5*DELTAX**2-0.5
               JUNK(5) = DELTAX*DELTAY
               JUNK(6) = 1.5*DELTAY**2-0.5
               IF (NEXP .GE. 7) THEN
                  JUNK(7) = DELTAX*(5.*JUNK(4)-2.)/3.
                  JUNK(8) = JUNK(4)*DELTAY
                  JUNK(9) = DELTAX*JUNK(6)
                  JUNK(10) = DELTAY*(5.*JUNK(6)-2.)/3.
               END IF
            END IF
         END IF
      END IF
C
C     IF (NFRAC .GT. 0) THEN
C        J = NEXP+1
C        JUNK(J) = -2.*(DX - REAL(NINT(DX)))
C        J = J+1
C        JUNK(J) = -2.*(DY - REAL(NINT(DY)))
C        J = J+1
C        JUNK(J) = 1.5*JUNK(J-2)**2 - 0.5
C        J = J+1
C        JUNK(J) = JUNK(J-3)*JUNK(J-2)
C        J = J+1
C        JUNK(J) = 1.5*JUNK(J-3)**2 - 0.5
C     END IF
      XX = 2.*DX+MIDDLE
      LX = INT(XX)
      YY = 2.*DY+MIDDLE
      LY = INT(YY)
C
C This point in the stellar profile lies between columns LX and LX+1,
C and between rows LY and LY+1 in the look-up tables.
C
      DO K=1,NTERM
         CORR = BICUBC(PSF(LX-1,LY-1,K), MAXPSF, 
     .        XX-REAL(LX), YY-REAL(LY), DFDX, DFDY)
         USEPSF = USEPSF + JUNK(K)*CORR
         DVDXC = DVDXC-JUNK(K)*DFDX
         DVDYC = DVDYC-JUNK(K)*DFDY
      END DO
      RETURN
      END!
C
        function indexpos(head,f1)
        character*80 head(72),f1*8
        do 10 indexpos=1,72
10      if(head(indexpos)(1:8).eq.f1)return
        end


