c allstar 96,12
c itea 50-->200 97,4,29
        PARAMETER (NOPT=10)
	parameter (n2048=2048*3)

	real         a_b(n2048*n2048)
	real	     subt(n2048*n2048)
	real	     sigma(n2048*n2048)
	
	character*60  f1,f2,f3
	character*80  head(36),hhead(72)
      REAL OPT(NOPT)
C
      REAL PERERR, PROERR
      INTEGER NCOL, NROW, ISTAT
      LOGICAL CENTER,logi

C
        if(iargc().lt.1)then
        write(*,*)
        write(*,*)
        write(*,*)'     **********   DAO_STAR  2004.5 ***********'
        write(*,*)
        write(*,*)
	write(*,*)'     Usage:  dao_star file.fit'
	write(*,*)
	write(*,*)
	stop
	endif
	call getarg(1,f1)

        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='pip2.par'
        inquire(file=f3,exist=logi)
        if(.not.logi)f3='/vega2/qso5bin/pip2.par'
	open(1,file=f3,status='old')
21	read(1,'(a)')f3	
	if(f3(1:2).eq."**")goto 21
22	read(1,'(a)')f3	
	if(f3(1:2).ne."**")goto 22
23	read(1,'(a)')f3	
	if(f3(1:2).ne."**")goto 23
	do 30 i=1,nopt
31	read(1,'(a)')f3
	if(f3(1:2).eq."**")goto 31
30	read(f3(1:),'(29x,f9.2)')opt(i)
	close(1)

        k=index(f1,'.')
        f2=f1(1:k)//'psf'			! 1 read  *.psf
        open(1,file=f2,status='old')
        f2=f1(1:k)//'ap'			! 2 read  *.ap
        open(2,file=f2,status='old')
        f2=f1(1:k)//'als'			! 3 write *.als
        open(3,file=f2,status='UNKNOWN')
	call readfits(f1,a_b,ncol,nrow)
          if(ncol.eq.2048)is6=1
	  if(ncol.eq.6144)is6=3
	  opt(1)=opt(1)*is6
	  opt(9)=opt(9)*is6
	  opt(10)=opt(10)*is6
	write(*,'(7f10.2)')opt
C
      IF (OPT(3) .GE. 0.5) THEN
         CENTER=.TRUE.
      ELSE 
         CENTER=.FALSE.
      END IF
      PERERR = 0.01*OPT(7)
      PROERR = 0.01*OPT(8)
C
      CALL ALLSTR (a_b, NCOL, NROW, SUBT, SIGMA, OPT(1), 
     .     OPT(5), OPT(4), NINT(OPT(2)), CENTER, NINT(OPT(6)), 
     .     PERERR, PROERR, OPT(9), OPT(10), ISTAT)
C
	if(istat.eq.9000)then
	f3=f1(1:lnblnk(f1))//'s'//char(0)
	open(1,file=f1,status='UNKNOWN',access='direct',recl=2880)
	i=0
502	i=i+1
	read(1,rec=i)head
	call writefits(f3,head,2880,i)
	if(indexpos(head,"END     ",36).gt.36)goto 502
	close(1)
c	call swap4(subt,ncol*nrow*4)
c	call writefits(f3,subt,ncol*nrow*4,0)
	call swap4(a_b,ncol*nrow*4)
	call writefits(f3,a_b,ncol*nrow*4,0)

	open(1,file=f1,status='UNKNOWN',access='direct',recl=5760)
	read(1,rec=1)hhead
	close(1)
	call writemes(hhead,f1)
	write(*,*)'  ok!'
	endif
	
      END!

	subroutine writemes(head,f3)
	character*80 head(72)
	character*60 f3,ff
	logical logi
	
        k=index(f3,'.')
        f3=f3(1:k)//'mes'
	call system('date >'//f3)
	
	open(1,file=f3,status='old')
	read(1,'(a)')ff
	close(1)
	open(1,file=f3,status='old')
	write(1,'(a)')'********'//ff
        call wmes(head,'NAXIS   ')
        call wmes(head,'NAXIS1  ')
        call wmes(head,'NAXIS2  ')
        call wmes(head,'IMNAME  ')
        call wmes(head,'EXPOSURE')
        call wmes(head,'RA      ')
        call wmes(head,'DEC     ')
        call wmes(head,'VOLT2   ')
        call wmes(head,'VOLT3   ')
	call wmes(head,'OBJECT  ')
        call wmes(head,'OBSERVER')
        call wmes(head,'A81     ')
        call wmes(head,'A82     ')
        call wmes(head,'A83     ')
        call wmes(head,'A84     ')
        call wmes(head,'A85     ')
        call wmes(head,'A86     ')
        call wmes(head,'A87     ')
        call wmes(head,'A88     ')
        call wmes(head,'RA2000  ')
        call wmes(head,'DEC2000 ')
	call wmes(head,'GALLONG ')
	call wmes(head,'GALLATI ')
	call wmes(head,'EXTINCT ')
	close (1)
        inquire(file='pip2.par',exist=logi)
        if(logi)call system('cat pip2.par >>'//f3)
        if(.not.logi)call system('cat /vega2/qso5bin/pip2.par >>'//f3)
	end

	subroutine wmes(head,a)
	character*80 head(72)
	character*8 a
	ip=indexpos(head,a,72)
	if(ip.lt.72)write(1,'(a)')head(ip)
	end
C
C#######################################################################
C
      SUBROUTINE  ALLSTR  (DATA, NCOL, NROW, SUBT, SIGMA, FITRAD, 
     .     WATCH, HALF, IEXP, CENTER, MAXGRP, 
     .     PERERR, PROERR, SKYIN, SKYOUT, ISTAT)
C
C=======================================================================
C
C Photometry for many stars by simultaneous multiple PSF fits iterating
C the entire star list simultaneously.
C 
C              OFFICIAL DAO VERSION:  1987 June 30
C
C Input arguments:
C
C       DATA is the input image, stored as a two-dimensional REAL 
C            array.
C
C NCOL, NROW are the dimensions of the input image.
C
C       SUBT is the address of some working space where ALLSTAR can
C            keep a scratch copy of the image.
C
C     SIGMA is the address of some working space where ALLSTAR can
C            keep track of the standard errors of the individual pixels
C
C     FITRAD is the user-definable parameter "Fitting radius."  It
C            governs how many pixels out from the centroid of the star 
C            will actually be considered in computing the least-
C            squares profile fits.
C
C      WATCH is the 'watch progress' parameter specified by the user.  
C            If WATCH > 0, information relating to the progress of
C            the reductions will be typed on the terminal during 
C            execution.
C
C HALF, IEXP are the "Clipping radius" and "Clipping exponent" 
C            parameters which specify the degree to which discordant
C            pixels are to be ignored.
C
C     CENTER is a logical specifying whether the user wants the stellar
C            centroids to be redetermined.
C
C     MAXGRP is the largest group size for which the solution is to be
C            performed.
C
C=======================================================================
C
      IMPLICIT NONE
      INTEGER MAXSTR, MAXPSF, MAXIT, MAXMAX, MAXPAR, MAXEXP
      PARAMETER (MAXSTR=50000, MAXPSF=207, MAXIT=200, MAXMAX=100,
     .     MAXPAR=6, MAXEXP=10)
C
C Parameters:
C
C MAXSTR The maximum number of stars in a frame; also the maximum
C        number of pixels in the sky annulus.  
C
C MAXPSF the largest PSF look-up table that can be accomodated.  If
C        PSFMAX is the largest acceptable PSF radius (see main-line
C        program above), then MAXPSF = 2*[2*(PSFMAX+1)]+1.
C
C MAXMAX is the largest group for which a solution will ever be
C        attempted = maximum permissible value of MAXGRP.
C
      INTEGER NCOL, NROW
      CHARACTER LINE*80
      REAL C(3*MAXMAX+1,3*MAXMAX+1), V(3*MAXMAX+1), X(3*MAXMAX+1)
      REAL DATA(NCOL,NROW), SUBT(NCOL,NROW), SIGMA(NCOL,NROW)
      REAL PSF(MAXPSF,MAXPSF,MAXEXP), PAR(MAXPAR)
      REAL XC(MAXSTR+1), YC(MAXSTR+1), MAG(MAXSTR+1), SKY(MAXSTR+1)
      REAL RPIXSQ(MAXSTR), CHI(MAXSTR), SUMWT(MAXSTR)
      REAL NUMER(MAXSTR), DENOM(MAXSTR)
      REAL MAGERR(MAXSTR), DXOLD(MAXSTR), DYOLD(MAXSTR)
      REAL XCLAMP(MAXSTR), YCLAMP(MAXSTR)
      INTEGER ID(MAXSTR+1)
      INTEGER NPIX(MAXSTR)
      LOGICAL SKIP(MAXSTR), LAST(MAXSTR)
C
      REAL AMIN1, AMAX1, ABS, SQRT, USEPSF, PCTILE
      INTEGER MIN0, MAX0, INT, RDPSF
C
      REAL LOBAD, HIBAD, ERR, SEP, THRESH, AP1, PHPADU, RONOIS, DUM
      REAL FAINT, CLMPMX, WCRIT, RADIUS, XMIN, XMAX, YMIN, YMAX
      REAL PERERR, PROERR, SKYIN, SKYOUT, SKYISQ, SKYOSQ, RELERR
      REAL PSFMAG, BRIGHT, XPSF, YPSF, SEPCRT, SEPMIN, PSFRAD, PEAK
      REAL FITRAD, WATCH, HALF, PKERR, PSFRSQ, RADSQ, RSQ, DX, DYSQ
      REAL SIGSQ, WT, VAL, RHOSQ, DWT, SKYBAR, CHIGRP, SUMRES, GRPWT
      REAL DPOS, DFDSIG, XKWT, DF, DIFF, SHARP, DELTAX, DELTAY, D, DY
      REAL DVDXC, DVDYC, Y
      INTEGER I, J, K, L, IFAINT, IDUM, IY, IX, NSTAR, NCONV, NITER
      INTEGER MAXGRP, MINSKY, MAXUNK, LX, LY, ISTAT, INTRVL
      INTEGER IPSTYP, NPSF, NPAR, NEXP, NFRAC, NL, ITSKY
      INTEGER IEXP, IXMIN, IXMAX, IYMIN, IYMAX, MX, MY, NDISAP
      INTEGER I3, J3, LSTAR, NSTR, NTERM, I3M2, ISTAR, N
      LOGICAL CENTER, CLIP, OMIT, REDO
C
      DATA SUMWT /MAXSTR*1./
C
C-----------------------------------------------------------------------
C
C SECTION 1
C
C Get ready, get set, . . .
C
      SKYISQ = SKYIN**2
      SKYOSQ = AMIN1(SKYOUT**2, SKYISQ + (MAXSTR-1)/3.1415927)
      SKYOUT = AMIN1(SKYOUT, SQRT(SKYOSQ))
      MINSKY = 100.
      MAXUNK=MAXMAX*3+1            ! Largest possible number of unknowns
C
C Read the point-spread function into memory.
C
      ISTAT = RDPSF (1, IPSTYP, PAR, MAXPAR, NPAR,
     .     PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, 
     .     PSFMAG, BRIGHT, XPSF, YPSF)
	close(1)				! *.psf
      PEAK = USEPSF(IPSTYP, 0., 0., BRIGHT, PAR, PSF, NPSF, 
     .     NPAR, NEXP, NFRAC, 0., 0., DVDXC, DVDYC)
C
C Stars will be checked for merger if they are separated by less than
C 1 FWHM of the image core.
C
C     Crit. sep. = 2.355*sigma, where
C          sigma = SQRT [ (sigma(X)**2 + sigma(Y)**2)/2 ]
C
      SEPCRT=2.*(PAR(1)**2+PAR(2)**2)
C
C Stars will be considered unconditionally merged if they are separated
C by less than about 0.375 * FWHM.
C
      SEPMIN=AMIN1(1.,0.14*SEPCRT)
C
C SEPCRT contains the square of the critical separation.
C SEPMIN  contains the square of the minimum separation.
C
      PKERR=PROERR/(PAR(1)*PAR(2)) ! **2      ! See fitting errors below
C
      PSFRAD = (REAL(NPSF-1)/2.-1.)/2.
      PSFRSQ = PSFRAD**2
C
C Ascertain the name of the file with the input photometry, and open it.
C
	read(2,'(a)')idum			! *.ap
      read(2,*) NL,IDUM,IDUM,LOBAD,HIBAD,THRESH,AP1,PHPADU,RONOIS,DUM
C
C Inquire the name of the output file, and open it.
C
      CALL WRHEAD (3, 1, NCOL, NROW, 7, LOBAD, HIBAD, THRESH, AP1, 
     .     PHPADU, RONOIS, FITRAD)
C
C We need to keep track of the anticipated standard error of the 
C brightness value in a given pixel as the reductions proceed.
C
      RONOIS=RONOIS**2
C
      LX = 1
      L = 1
      DO IY=1,NROW
         DO IX=1,NCOL
            IF ((DATA(IX,IY) .GT. HIBAD) .OR. 
     .           (DATA(IX,IY) .LT. LOBAD)) THEN
               SIGMA(IX,IY) = -1.1E38
            ELSE
               SIGMA(IX,IY) = RONOIS
            END IF
         END DO
      END DO
C
      RADSQ=FITRAD**2
      NSTAR=0
C
C Read in all the stars.
C
      FAINT=PSFMAG+12.5
      CLMPMX=0.25*FITRAD
      I=0
 1110 I=I+1
      IF (I .GT. MAXSTR) THEN
         write(*,*)'Too many stars.'
         GO TO 2000
      END IF
 1120 CALL RDSTAR (2, NL, ID(I), XC(I), YC(I), MAG(I), SKY(I))
      IF (ID(I) .LT. 0) GO TO 2000        ! End-of-file was encountered
      IF (ID(I) .EQ. 0) GO TO 1120        ! A blank line was encountered
      IF (MAG(I) .LT. FAINT) THEN
         MAG(I)=10.**(0.4*(PSFMAG-MAG(I)))
      ELSE
         MAG(I)=0.003
      END IF
      DXOLD(I)=0.
      DYOLD(I)=0.
      XCLAMP(I)=CLMPMX
      YCLAMP(I)=CLMPMX
      GO TO 1110
C
C-----------------------------------------------------------------------
C
C SECTION 2
C
C GO.
C
 2000 NSTAR=I-1                                      ! Number of stars
	close(2)
      WRITE (6,5552) NSTAR
 5552 FORMAT (I6, ' stars.  <<')

c      WRITE (6,610)
c  610 FORMAT (//' I = iteration number', 
c     .        //' R = number of stars that remain'
c     .        //' D = number of stars that disappeared'
c     .        //' C = number of stars that converged'
c     .       //)
c      CALL OVRWRT ('     I     R     D     C', 1)
C
C Initialize accumulators and constraints on parameter corrections.
C
      NCONV=0
      NITER=0
      CALL STRIP (ID, XC, YC, MAG, SKY, SKIP, MAXSTR,
     .     NSTAR, NDISAP, SQRT(SEPMIN), NUMER, DENOM)
      IF (NDISAP .GT. 0) THEN
         WRITE (LINE,682) NITER, NSTAR, NDISAP, NCONV
         CALL OVRWRT (LINE(1:28), 3)
      END IF
C
      ITSKY = 1
 2100 NITER = NITER+1
      CLIP = (IEXP .NE. 0) .AND. (NITER .GE. 4)
C
C Set up critical errors for star rejection.
C
      WCRIT = 0.
      IF (NITER .GE. 4) WCRIT = AMIN1(2.,(NITER+5.)/10.)
C
C Sort stars by y-coordinate, for minimum paging.
C
      CALL QUICK (YC, NSTAR, NUMER)
      CALL RECTFYI (ID, NSTAR, NUMER, DENOM)
      CALL RECTFY (XC, NSTAR, NUMER, DENOM)
      CALL RECTFY (MAG, NSTAR, NUMER, DENOM)
      CALL RECTFY (SKY, NSTAR, NUMER, DENOM)
      CALL RECTFY (SUMWT, NSTAR, NUMER, DENOM)
      CALL RECTFY (DXOLD, NSTAR, NUMER, DENOM)
      CALL RECTFY (DYOLD, NSTAR, NUMER, DENOM)
      CALL RECTFY (XCLAMP, NSTAR, NUMER, DENOM)
      CALL RECTFY (YCLAMP, NSTAR, NUMER, DENOM)
C
C Determine how much of the image we need to copy into the working area.
C
C First determine the big rectangular area that contains all the pixels
C within one working radius of any star.  The "working radius" is
C either the fitting radius or the outer sky radius (if this is an
C iteration during which the sky is to be determined), whichever is
C larger.  Set their SIGMAs negative.
C
      IF ((NITER .EQ. ITSKY) .AND. (SKYOUT .GT. FITRAD))
     .     THEN
         RADIUS = SKYOUT
      ELSE
         RADIUS = FITRAD
      END IF
C
      XMIN=NCOL+RADIUS
      XMAX=-RADIUS
      DO I=1,NSTAR
         XMIN=AMIN1(XMIN, XC(I))
         XMAX=AMAX1(XMAX, XC(I))
      END DO
C
      IXMIN=MAX0(1, MIN0(NCOL, INT(XMIN-RADIUS)+1))
      IXMAX=MAX0(1, MIN0(NCOL, INT(XMAX+RADIUS)))
      IYMIN=MAX0(1, MIN0(NROW, INT(YC(1)-RADIUS)+1))
      IYMAX=MAX0(1, MIN0(NROW, INT(YC(NSTAR)+RADIUS)))
C
      IF (WATCH .GT. 0.5) THEN
C
C Set up reporting intervals for row counter.
C
         D=AMAX1(0.1, ALOG10(FLOAT(IYMAX-IYMIN+1)/10.))
         INTRVL=INT(D)
         D=D-FLOAT(INTRVL)
         IF (D .GT. 0.8) THEN
            INTRVL=10**(INTRVL+1)
         ELSE IF (D .GT. 0.5) THEN
            INTRVL=5*10**INTRVL
         ELSE IF (D .GT. 0.2) THEN
            INTRVL=2*10**INTRVL
         ELSE
            INTRVL=10**INTRVL
         END IF
      END IF
C
C Beginning of big loop over y-coordinates.
C
      RSQ = RADIUS**2
      N = 0
      L = 0
      DO 2190 IY=IYMIN,IYMAX
         Y = REAL(IY)
         DO IX=IXMIN,IXMAX
            SIGMA(IX,IY) = -ABS(SIGMA(IX,IY))
         END DO
C
C OK.  Now find all the points within one working radius of each 
C star.  Set the sigmas of these pixels positive again and copy 
C them from DATA to SUBT.
C 
C N points at the last star that has such a low x-coordinate
C that it can be omitted.
C
         DO 2120 I=N+1,NSTAR
            DY = Y - YC(I)
C
            IF (DY .GT. RADIUS) THEN
               N = I
               GO TO 2120
            ELSE IF (DY .LT. -RADIUS) THEN
               GO TO 2150
            ELSE
               DYSQ = (REAL(IY) - YC(I))**2
               LX=MAX0(1, MIN0(NCOL, INT(XC(I)-RADIUS)+1))
               MX=MAX0(1, MIN0(NCOL, INT(XC(I)+RADIUS)))
               DO IX=LX,MX
                  DX = REAL(IX) - XC(I)
                  IF (DX**2+DYSQ .LE. RSQ) THEN
                     IF (SIGMA(IX,IY) .GT. -1.E38) THEN
                        SIGMA(IX,IY) = ABS(SIGMA(IX,IY))
                        SUBT(IX,IY) = DATA(IX,IY)
                     END IF
                  ELSE
                     IF (DX .GE. 0.) GO TO 2120
                  END IF
               END DO
            END IF
 2120    CONTINUE
C
C Subtract from the working copy of the image (SUBT) each star which 
C has not yet converged by subtracting the shifted, scaled PSF 
C according to the current best parameter estimates; do this only
C for pixels with non-negative SIGMAs, indicating that
C they will actually be needed for this iteration.
C
C L points at the last star that has such a low y-coordinate
C that it can be omitted.
C
 2150    CONTINUE
         DO 2170 I=L+1,NSTAR
            DY = Y - YC(I)
            IF (DY .GT. PSFRAD) THEN
               L = I
               GO TO 2170
            ELSE IF (DY .LT. -PSFRAD) THEN
               GO TO 2180
            ELSE
               LX=MAX0(IXMIN, MIN0(IXMAX, INT(XC(I)-PSFRAD)+1))
               MX=MAX0(IXMIN, MIN0(IXMAX, INT(XC(I)+PSFRAD)))
               DELTAX=(XC(I)-1.)/XPSF-1.
               DELTAY=(YC(I)-1.)/YPSF-1.
               DYSQ=DY**2
               DO 2160 IX=LX,MX
                  IF (SIGMA(IX,IY) .GE. 0.) THEN
                     DX=FLOAT(IX)-XC(I)
                     IF (DX**2+DYSQ .LT. PSFRSQ) THEN
                        SUBT(IX,IY)=SUBT(IX,IY)-
     .                       MAG(I)*USEPSF(IPSTYP, DX, DY, BRIGHT, 
     .                       PAR, PSF, NPSF, NPAR, NEXP, NFRAC, 
     .                       DELTAX, DELTAY, DVDXC, DVDYC)
                     ELSE
                        IF (DX .GT. 0.) GO TO 2170
                     END IF
                  END IF
 2160          CONTINUE
            END IF
 2170    CONTINUE
 2180    IF ((WATCH .GT. 0.5) .AND. (NSTAR .GT. 100) .AND.
     .        (INTRVL*(IY/INTRVL) .EQ. IY)) THEN
            WRITE (LINE, 5551) IYMIN, IY, IYMAX
 5551       FORMAT (3I5)
            CALL OVRWRT (LINE(1:15), 2)
         END IF
 2190 CONTINUE
C
C If sky values are to be redetermined this iteration, do it now.
C
      IF ((SKYOUT .GT. 0.5) .AND. (NITER .EQ. ITSKY))
     .     THEN
         ITSKY = ITSKY + 1 + NITER/5
         DO 2198 ISTAR=1,NSTAR
            LX=MAX0(IXMIN, MIN0(IXMAX, INT(XC(ISTAR)-SKYOUT)+1))
            MX=MAX0(IXMIN, MIN0(IXMAX, INT(XC(ISTAR)+SKYOUT)))
            LY=MAX0(IYMIN, MIN0(IYMAX, INT(YC(ISTAR)-SKYOUT)+1))
            MY=MAX0(IYMIN, MIN0(IYMAX, INT(YC(ISTAR)+SKYOUT)))
            N = 0
c           nbad = 0
            DO 2196 J=LY,MY
               DYSQ = (REAL(J) - YC(ISTAR))**2
               DO 2194 I = LX,MX
                  DX = REAL(I) - XC(ISTAR)
                  RSQ = DX**2 + DYSQ
                  IF (RSQ .GT. SKYOSQ) THEN
                     IF (DX .GT. 0.) THEN
                        GO TO 2196
                     ELSE
                        GO TO 2194
                     END IF
                  END IF
C
                  IF (RSQ .GE. SKYISQ) THEN
                     IF (SIGMA(I,J) .GT. -1.E38) THEN
                        N = N+1
                        NUMER(N) = SUBT(I,J)
                        IF (N .GE. MAXSTR) GO TO 2197
c                       else
c                          nbad = nbad+1
                     END IF
                  END IF
C
 2194          CONTINUE
 2196       CONTINUE
C
 2197       IF (N .GT. MINSKY) THEN
               D = PCTILE (NUMER, N, (N+1)/2)
               DF = 4.*(RONOIS + AMAX1(0.,D)/PHPADU)
               DX = 0.
               DY = 0.
               DO I=1,N
                  DPOS = NUMER(I) - D
                  WT = 1./(1.+(DPOS**2/DF)**2)
                  DX = DX + WT*DPOS
                  DY = DY + WT
               END DO
               SKY(ISTAR) = D + DX/DY
            END IF
 2198    CONTINUE
      END IF
C
C Create provisional star groups.
C
      CALL REGRP (ID, XC, YC, MAG, SKY, SUMWT, DXOLD, DYOLD, 
     .     XCLAMP, YCLAMP, MAXSTR, NSTAR, FITRAD, LAST, NUMER, 
     .     DENOM)
C
C Now get ready to do the next iteration, group by group.
C
      IF (WATCH .GT. 0.5) THEN
C
C Set up reporting intervals for star counter.
C
         D=AMAX1(0.1, ALOG10(FLOAT(NSTAR)/10.))
         INTRVL=INT(D)
         D=D-FLOAT(INTRVL)
         IF (D .GT. 0.8) THEN
            INTRVL=10**(INTRVL+1)
         ELSE IF (D .GT. 0.5) THEN
            INTRVL=5*10**INTRVL
         ELSE IF (D .GT. 0.2) THEN
            INTRVL=2*10**INTRVL
         ELSE
            INTRVL=10**INTRVL
         END IF
      END IF
      RADIUS=FITRAD
      ISTAR=1
 2200 CONTINUE
C
C Find the last star in the current group.
C
      DO LSTAR = ISTAR,NSTAR
         SKIP(LSTAR) = .FALSE.
         IF (LAST(LSTAR)) GO TO 2210
      END DO
      LSTAR=NSTAR
 2210 NSTR=LSTAR-ISTAR+1
C
C Start crunching on this group.
C
      IF (NSTR .GT. MAXGRP) THEN
         IF (WATCH .GT. 1.5) THEN
            WRITE (LINE,61) NSTR, RADIUS
   61       FORMAT ('Group too large:', I6, '  (', F5.2, ')')
            CALL OVRWRT (LINE(1:30), 3)
         END IF
         RADIUS=0.95*RADIUS
         IF (CENTER .AND. (NITER .GE. 2)) THEN
            IF (RADIUS .LT. 1.2) THEN
               IF (WATCH .GT. 1.5) THEN
                  WRITE (LINE,62)
   62             FORMAT ('Group too dense to reduce.')
                  CALL OVRWRT (LINE(1:26), 3)
               END IF
C
C Mark the faintest star in this group for deletion.
C
               FAINT=100.
               DO I=ISTAR,LSTAR
                  IF (MAG(I) .LT. FAINT) THEN
                     FAINT=MAG(I)
                     L=I
                  END IF
               END DO
               IF (WATCH .GT. 1.5) THEN
                  WRITE (6,63) ID(L), XC(L), YC(L), 
     .                 PSFMAG-1.085736*ALOG(MAG(L)), ' group too dense'
   63             FORMAT (I6, 2F9.3, F9.3, 9X, A, F6.3)
               END IF
               IF (.NOT. SKIP(L)) NDISAP=NDISAP+1
               SKIP(L)=.TRUE.                        ! Flag for deletion
               GO TO 3000
            END IF
         ELSE 
            IF (RADIUS .LT. 0.8) THEN
               IF (WATCH .GT. 1.5) THEN
                  WRITE (LINE,62)
                  CALL OVRWRT (LINE(1:26), 3)
               END IF
C
C Mark the faintest star in this group for deletion.
C
               FAINT=100.
               DO I=ISTAR,LSTAR
                  IF (MAG(I) .LT. FAINT) THEN
                     FAINT=MAG(I)
                     L=I
                  END IF
               END DO
               IF (WATCH .GT. 1.5) WRITE (6,63) ID(L), XC(L), YC(L), 
     .              PSFMAG-1.085736*ALOG(MAG(L)), ' group too dense'
               IF (.NOT. SKIP(L)) NDISAP=NDISAP+1
               SKIP(L)=.TRUE.                        ! Flag for deletion
               GO TO 3000
            END IF
         END IF
         L = NSTR
         CALL REGRP (ID(ISTAR), XC(ISTAR), YC(ISTAR), MAG(ISTAR), 
     .        SKY(ISTAR), SUMWT(ISTAR), DXOLD(ISTAR), DYOLD(ISTAR), 
     .        XCLAMP(ISTAR), YCLAMP(ISTAR), L, NSTR, RADIUS, 
     .        LAST(ISTAR), NUMER(ISTAR), DENOM(ISTAR))
         GO TO 2200
      END IF
C
      RADIUS=FITRAD
      XMIN=NCOL+FITRAD
      XMAX=-FITRAD
      YMIN=NROW+FITRAD
      YMAX=-FITRAD
C
C Get the extent of the rectangular area of the frame relevant to this
C group.  At the same time, determine the average SKY and the average
C CHI value for the group, and zero the various statistical 
C accumulators.  What we would REALLY like to do is to compute the
C anticipated error for each pixel the same way as PEAK and NSTAR do.
C We can't really do that, because stellar groups aren't defined the
C same way as in NSTAR: here they change as stars move from one group
C to another, or converge and are subtracted from the frame.  Therefore,
C there is no consistent way to define a group's CHI value the same as
C was done in NSTAR.  I will take the cheap way out:  using the
C individual smoothed CHI values for the stars which were obtained
C during the previous iteration (which should still be stored in the
C array SUMWT) I will simply average those values and call the mean the
C group's effective smoothed CHI.  This way, at least, the reductions
C will be the same as before for isolated stars.  (Note that,
C for the first iteration, the SUMWT's have been set to 1.0 by a DATA
C statement above.)
C
      SKYBAR = 0.
      CHIGRP = 0.
      DO I=ISTAR,LSTAR
         XMIN = AMIN1(XMIN, XC(I))
         XMAX = AMAX1(XMAX, XC(I))
         YMIN = AMIN1(YMIN, YC(I))
         YMAX = AMAX1(YMAX, YC(I))
         SKYBAR = SKYBAR+SKY(I)
         CHIGRP = CHIGRP+SUMWT(I)
         CHI(I) = 0.
         SUMWT(I) = 0.
         NUMER(I) = 0.
         DENOM(I) = 0.
         NPIX(I) = 0
      END DO
      SKYBAR = SKYBAR/FLOAT(NSTR)
      CHIGRP = CHIGRP/FLOAT(NSTR)
      IF (CENTER .AND. (NITER .GE. 2)) THEN
         NTERM=3*NSTR
      ELSE
         NTERM=NSTR
      END IF
C
C If sky is to be determined: NTERM=NTERM+1
C
C Now... on with the iteration.
C
      IXMIN=MIN0(NCOL, MAX0(1, INT(XMIN-FITRAD)+1))
      IXMAX=MIN0(NCOL, MAX0(1, INT(XMAX+FITRAD)))
      IYMIN=MIN0(NROW, MAX0(1, INT(YMIN-FITRAD)+1))
      IYMAX=MIN0(NROW, MAX0(1, INT(YMAX+FITRAD)))
C
C IXMIN, IXMAX, IYMIN, and IYMAX are now the limits of a rectangular 
C array containing all pixels within one fitting radius of any star in 
C the group.
C
C Zero the normal matrix and the vector of residuals.
C
      DO 2270 J=1,NTERM
      V(J)=0.0
      DO 2270 I=J,NTERM
 2270 C(I,J)=0.0
      SUMRES=0.
      GRPWT=0.
C
C Now deal with the pixels one by one.
C
      DO 2390 IY=IYMIN,IYMAX
      DO 2380 IX=IXMIN,IXMAX
      IF (SIGMA(IX,IY) .LT. 0.) GO TO 2380
C
C If this pixel is within one fitting radius of at least one star in
C the current group, include it in the calculation.  Otherwise, skip 
C it.  While figuring this out, compute the squared distance of this 
C pixel from the centroid of each star in the group.  
C
      OMIT=.TRUE.
      DO I=ISTAR,LSTAR
         RPIXSQ(I)=(FLOAT(IX)-XC(I))**2+(FLOAT(IY)-YC(I))**2
         IF (RPIXSQ(I) .LT. RADSQ) THEN
            OMIT=.FALSE.
         END IF
      END DO
      IF (OMIT) GO TO 2380                      ! Do not need this pixel
C
      DO I=ISTAR,LSTAR
         IF (RPIXSQ(I) .LT. RADSQ) THEN
            SKIP(I)=.FALSE.
         ELSE
            SKIP(I)=.TRUE.
         END IF
      END DO
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 0.75% 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 
C proportional to H/sigma**4 (sigma is the Gaussian 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 (see definition of PKERR above.)
C
      D=SUBT(IX,IY)-SKYBAR                      ! Residual of this pixel
      DPOS=AMAX1(0., DATA(IX,IY)-D)
      IF ((DPOS .GT. HIBAD) .AND. (NITER .GE. 4)) THEN
         GO TO 2380
      END IF
C
C    DPOS = raw data minus residual
C         = model-predicted brightness in this pixel, consisting of sky
C           plus all stellar profiles, which presumably is non-negative.
C
C The four error sources in our noise model are:
C     (1) Readout noise
C     (2) Poisson noise
C     (3) Flat-field errors
C     (4) Errors in the PSF
C
C Numerically, the squares of these quantities are
C     (1) RONOIS = SIGMA(IX,IY)   (initially, at least)
C     (2) DPOS/PHPADU where DPOS = sum of the stellar profiles plus sky
C     (3) constant x DPOS**2
C     (4) constant x the sum of (stellar profile)**2
C
C Here's the thing:  at this point we don't have the 
C sum of (stellar profile)**2 available to us.  All we have is 
C (sum of stellar profile)**2, which isn't the same thing at all.
C Nevertheless, we have to assume that 
C
C     sum of (stellar profile)**2 = (sum of stellar profile)**2;
C
C otherwise I'd have to use a whole new REAL array as big as the
C original CCD image, and I don't want to do that.  This should be
C good enough for pixels where either (a) the other three error sources
C dominate, .OR. (b) one of the stars putting light into that pixel is
C much brighter than all the rest.  Thus, the approximation should be
C good UNLESS the pixel is near the center of two or more bright
C stars --- in which case things will be hairy anyway.
C
      SIGSQ = SIGMA(IX,IY) + DPOS/PHPADU + (PERERR*DPOS)**2 +
     .     (PKERR*(DPOS-SKYBAR))**2
      RELERR=ABS(D)/SQRT(SIGSQ)
      IF (CLIP .AND. (RELERR .GT. 100.)) GO TO 2380
      WT=0.
C
C Now include this pixel in the fitting equations for the group.
C
      DO 2320 I=ISTAR,LSTAR
      IF (SKIP(I)) GO TO 2320
      RSQ=RPIXSQ(I)/RADSQ
      IF (RSQ .GE. 0.999999) GO TO 2320  ! Safety check vs divide by 0
      WT=AMAX1(WT, 5./(5.+RSQ/(1.-RSQ)))
C
C The condition equation for pixel (IX,IY) is of the form
C
C data(IX,IY)-summation{scale*psf(IX-Xcenter,IY-Ycenter)}-sky=residual
C
C Then we will jigger the scale's, Xcenter's, and Ycenter's such that
C
C                Summation{weight * residual**2}
C
C is minimized.  'weight' will be a function (1) of the distance of this
C pixel from the center of the nearest star, (2) of the model-predicted
C brightness of the pixel (taking into consideration the readout noise, 
C the photons/ADU, and the interpolation error of the PSF), and (3) of 
C the size of the residual itself.  (1) is necessary to prevent the
C non-linear least-squares solution from oscillating:  oft-times it will
C come to pass that if you include a pixel in the solution, then the
C predicted shift of the centroid will cause that pixel to be excluded 
C in the next iteration, and the new predicted shift of the centroid
C will cause that pixel to be included again.  This could go on ad
C infinitum.  The cure is to have the weight of a pixel go 
C continuously to zero as its distance from the stellar centroid
C approaches the fitting radius.  In a case like that just described,
C the solution can then find a real minimum of the sum of the
C weighted squared residuals with that pixel at some low-weight position
C just inside the fitting radius.  (2) is just sensible weighting.
C (3) is just a crude attempt at making the solution more robust against
C bad pixels.
C
      VAL = USEPSF(IPSTYP, FLOAT(IX)-XC(I), FLOAT(IY)-YC(I), BRIGHT, 
     .     PAR, PSF, NPSF, NPAR, NEXP, NFRAC, (XC(I)-1.)/XPSF-1., 
     .     (YC(I)-1.)/YPSF-1., DVDXC, DVDYC)
      IF (NTERM .GT. NSTR) THEN
         I3=(I-ISTAR+1)*3
         K=I3-2
         X(K)=-VAL
         K=I3-1
         X(K)=-MAG(I)*DVDXC
         X(I3)=-MAG(I)*DVDYC
      ELSE
         K=I-ISTAR+1
         X(K)=-VAL
      END IF
      RHOSQ=((XC(I)-FLOAT(IX))/PAR(1))**2+
     .     ((YC(I)-FLOAT(IY))/PAR(2))**2
      IF (RHOSQ .LE. 36.) THEN
         RHOSQ = 0.6931472*RHOSQ
         DFDSIG = EXP(-RHOSQ)*(RHOSQ-1.)
         NUMER(I) = NUMER(I)+DFDSIG*D/SIGSQ
         DENOM(I) = DENOM(I)+DFDSIG**2/SIGSQ
      END IF
 2320 CONTINUE
C
C At this point, the vector X contains the first derivative of
C the condition equation for pixel (IX,IY) with respect to each of
C the fitting parameters for all of the stars. Now these derivatives 
C will be added into the normal matrix and the vector of residuals.
C
C Add this residual into the weighted sum of the absolute relative 
C residuals.
C
      DWT=WT*RELERR
      SUMRES=SUMRES+DWT
      GRPWT=GRPWT+WT
C
C SUMRES is the weighted sum of [ABS(residual)/sigma] for all the 
C pixels in the group.  Now also add the weighted value of 
C [ABS(residual)/sigma] into the accumulating sum for each of the
C stars.
C
      DO 2330 I=ISTAR,LSTAR
      IF (SKIP(I)) GO TO 2330
      CHI(I)=CHI(I)+DWT
      SUMWT(I)=SUMWT(I)+WT
      NPIX(I)=NPIX(I)+1
 2330 CONTINUE
C
C Up until now, WT represents only the radial weighting profile.  Now
C figure in the anticipated standard error of the pixel.  Reject any
C pixel with a 100-sigma residual after iteration 2.
C
      WT=WT/SIGSQ
      IF (CLIP) WT=WT/(1.+(RELERR/(CHIGRP*HALF))**IEXP)
C
C Now work this pixel into the normal matrix.
C
      DWT=D*WT
C     If sky is to be determined: C(NTERM,NTERM)=C(NTERM,NTERM)+WT
C     If sky is to be determined: V(NTERM)=V(NTERM)-DWT
      DO 2370 I=ISTAR,LSTAR
      IF (SKIP(I)) GO TO 2370
      IF (NTERM .GT. NSTR) THEN
         I3=(I-ISTAR+1)*3
         I3M2=I3-2
         DO 2340 K=I3M2,I3
C        If sky is to be determined: C(NTERM,K)=C(NTERM,K)-X(K)*WT
 2340    V(K)=V(K)+X(K)*DWT
         DO 2360 J=ISTAR,I
         IF (SKIP(J)) GO TO 2360
         J3=(J-ISTAR+1)*3
         DO 2350 K=I3M2,I3
         XKWT=X(K)*WT
         DO 2350 L=J3-2,MIN0(K, J3)
 2350    C(K,L)=C(K,L)+X(L)*XKWT
 2360    CONTINUE
      ELSE
         K=I-ISTAR+1
         V(K)=V(K)+X(K)*DWT
         XKWT=X(K)*WT
C        If sky is to be determined: C(NTERM,K)=C(NTERM,K)-XKWT
         DO 2365 J=ISTAR,I
         IF (SKIP(J)) GO TO 2365
         L=J-ISTAR+1
         C(K,L)=C(K,L)+X(L)*XKWT
 2365    CONTINUE
      END IF
 2370 CONTINUE
C
 2380 CONTINUE
 2390 CONTINUE
C
C Reflect the normal matrix across the diagonal.
C
      IF (NTERM .GT. 1) THEN
         DO 2410 L=2,NTERM
         DO 2410 K=1,L-1
 2410    C(K,L)=C(L,K)
      END IF
C
C Compute the estimate of the standard deviation of the residuals for 
C the group as a whole, and for each star.  This estimate starts out as
C SQRT(PI/2)*{SUM[weight*ABS(residual/sigma)]/SUM(weight)} and then gets
C corrected for bias by SQRT(no. of pixels/(no. of pixels - degrees of
C freedom)).
C
      IF (GRPWT .GT. NTERM) THEN
         CHIGRP=1.2533141*SUMRES/SQRT(GRPWT*(GRPWT-NTERM))
C
C But then I drive the value toward unity, depending on exactly how
C many pixels were involved:  if CHI is based on exactly a total 
C weight of 3, then it is extremely poorly determined, and we just
C want to keep CHI = 1.  The larger GRPWT is, the better determined
C CHI is, and the less we want to force it toward unity.  So,
C just take the weighted average of CHI and unity, with weights
C GRPWT-3 and 1, respectively.
C
         CHIGRP=((GRPWT-3.)*CHIGRP+3.)/GRPWT
      ELSE
         CHIGRP=1.
      END IF
C
C CHIGRP has been pulled toward its expected value of unity to keep the 
C statistics of a small number of pixels from compeletely dominating 
C the error analysis.  Similarly, the photometric errors for the 
C individual stars will be pulled toward unity now.  Later on, if the
C number of stars in the group is greater than one, CHI(I) will be 
C nudged toward the group average.  In order to work optimally, of 
C course, this requires that PHPADU, RONOIS, and the other noise 
C contributors which I have postulated properly represent the true 
C errors expected in each pixel.
C
C At the same time, be sure that every star in the group contains at 
C least 3 valid pixels if recentroiding is being performed, 1 valid
C pixel if not.  If any star in the group fails to meet this criterion
C mark that star for deletion, set REDO to true, and skip ahead to the 
C next group.
C
      REDO=.FALSE.
      DO I=ISTAR,LSTAR
         SKIP(I)=.FALSE.                         ! Flag for retention
         IF (CENTER .AND. (NITER .GE. 2)) THEN
            IF (NPIX(I) .LT. 3) THEN
               REDO=.TRUE.
               IF (WATCH .GT. 1.5) WRITE (6,63) ID(I), XC(I), YC(I), 
     .              PSFMAG-1.085736*ALOG(MAG(I)), 
     .              ' too few valid pixels'
               IF (.NOT. SKIP(I)) NDISAP=NDISAP+1
               SKIP(I)=.TRUE.                    ! Flag for deletion
            ELSE
               I3=(I-ISTAR+1)*3-2
               IF (SUMWT(I) .GT. 3.) THEN
                  CHI(I)=1.2533141*CHI(I)/SQRT(SUMWT(I)*(SUMWT(I)-3.))
C
C Store a smoothed CHI value for the star in SUMWT.
C
                  SUMWT(I) = ((SUMWT(I)-3.)*CHI(I) + 3.)/SUMWT(I)
               ELSE
                  CHI(I)=CHIGRP
               END IF
            END IF
         ELSE
            IF (NPIX(I) .LT. 1) THEN
               REDO=.TRUE.
               IF (.NOT. SKIP(I)) NDISAP=NDISAP+1
               SKIP(I)=.TRUE.                    ! Flag for deletion
               IF (WATCH .GT. 1.5) WRITE (6,63) ID(I), XC(I), YC(I), 
     .              PSFMAG-1.085736*ALOG(MAG(I)), 
     .              ' too few valid pixels'
            ELSE
               I3=I-ISTAR+1
               IF (SUMWT(I) .GT. 1.) THEN
                  CHI(I)=1.2533141*CHI(I)/SQRT(SUMWT(I)*(SUMWT(I)-1.))
C
C Store a smoothed CHI value for the star in SUMWT.
C
                  SUMWT(I) = ((SUMWT(I)-3.)*CHI(I) + 3.)/SUMWT(I)
               ELSE
                  CHI(I)=CHIGRP
               END IF
            END IF
         END IF
      END DO
      IF (REDO) THEN
         GO TO 3000
      END IF
      CALL INVERS (C, MAXUNK, NTERM, ISTAT)
      DO J=1,NTERM
         IF (C(J,J) .LE. 0.) THEN
C
C Uh-oh.  Zero on the diagonal. Booma Booma!
C
            IF (CENTER .AND. (NITER .GE. 2)) THEN
               I = (J+2)/3 + ISTAR - 1
            ELSE
               I = J + ISTAR - 1
            END IF
C
C I now points at the (first, at least) star that caused the problem.
C
            IF (.NOT. SKIP(I)) NDISAP=NDISAP+1
            SKIP(I) = .TRUE.
            IF (WATCH .GT. 1.5) WRITE (6,63) ID(I), XC(I), YC(I), 
     .           PSFMAG-1.085736*ALOG(MAG(I)), ' singular matrix'
            GO TO 3000
         END IF
      END DO
      CALL VMUL (C, MAXUNK, NTERM, V, X)
C If sky is to be determined: SKYBAR=SKYBAR-X(NTERM)
C If sky is to be determined: IF(ABS(X(NTERM)).GT.0.01)REDO=.TRUE.
C
C In the beginning, the brightness of each star will be permitted to
C change by no more than two magnitudes per iteration, and the x,y 
C coordinates of each centroid will be permitted to change by no more 
C than 0.4 pixel per iteration.  Any time that the parameter
C correction changes sign from one iteration to the next, the maximum 
C permissible change will be reduced by a factor of two.  These
C clamps are released any time a star in the group disappears.
C
      DO 2520 I=ISTAR,LSTAR
         IF (CENTER .AND. (NITER .GE. 2)) THEN
C
C Correct both magnitude and position.
C
            L=3*(I-ISTAR+1)
            K=L-1
            J=L-2
C
C Note that the sign of the correction is such that it must be 
C SUBTRACTED from the current value of the parameter to get the 
C improved parameter value.  This being the case, if the correction
C to the brightness is negative (the least-squares thinks that the
C star should be brighter) a change of 1 magnitude is a change of a
C factor of 2.5; if the brightness correction is positive (the star
C should be fainter) a change of 1 magnitude is a change of 60%.
C
            DWT=DXOLD(I)*X(K)
            IF (DWT .LT. 0.) THEN
               XCLAMP(I)=AMAX1(0.0001,0.5*XCLAMP(I))
            ELSE
               XCLAMP(I)=AMIN1(CLMPMX,1.04*XCLAMP(I))
            END IF
            XC(I)=XC(I)-X(K)/(1.+ABS(X(K)/XCLAMP(I)))
            DXOLD(I)=X(K)
            DWT=DYOLD(I)*X(L)
            IF (DWT .LT. 0.) THEN
               YCLAMP(I)=AMAX1(0.0001,0.5*YCLAMP(I))
            ELSE
               YCLAMP(I)=AMIN1(CLMPMX,1.04*YCLAMP(I))
            END IF
            YC(I)=YC(I)-X(L)/(1.+ABS(X(L)/YCLAMP(I)))
            DYOLD(I)=X(L)
            MAG(I)=MAG(I)-X(J)/
     .           (1.+AMAX1( X(J)/(0.84*MAG(I)) , -X(J)/(5.25*MAG(I)) ))
            MAGERR(I) = SUMWT(I)*SQRT(C(J,J))
            IF (NITER .GE. 4) THEN
               REDO = .FALSE.
               IF (ABS(X(J)) .GT.
     .              AMAX1( 0.1*MAGERR(I) , 0.0005*MAG(I) )) THEN
                  REDO=.TRUE.
               ELSE
                  DF = (0.1*SUMWT(I))**2
                  IF (X(K)**2 .GT. AMAX1(DF*C(K,K), 4.E-6)) THEN
                     REDO=.TRUE.
                  ELSE IF (X(L)**2 .GT. AMAX1(DF*C(L,L), 4.E-6)) THEN
                     REDO=.TRUE.
                  END IF
               END IF
            ELSE
               REDO = .TRUE.
            END IF
         ELSE
            J=I-ISTAR+1
C
C Correct magnitude only. Since this is an easy, linear problem,
C elaborate clamps are not needed.  A simple one will do to keep
C the brightness from ever going negative.
C
            MAG(I)=MAG(I)-X(J)/(1. + 1.2*ABS( X(J)/MAG(I) ))
            MAGERR(I) = SUMWT(I)*SQRT(C(J,J))
            IF (NITER .GE. 2) THEN
               REDO=.FALSE.
               IF (ABS(X(J)) .GT.
     .              AMAX1( 0.1*MAGERR(I), 0.0005*MAG(I) )) REDO = .TRUE.
            ELSE
               REDO=.TRUE.
            END IF
         END IF
C
C No star will be allowed to converge with a signal-to-noise ratio
C worse than 2.0.  Such stars must be retained until either they get
C better than S/N = 2, or they get eliminated some time after
C iteration 5 (if S/N < 1.0), 10 (if S/N < 1.5), or 15 (if S/N < 2.0).
C However, anything which is lucky enough to survive until the last
C iteration will be written out regardless.
C
         IF (MAG(I) .LT. 2.0*MAGERR(I)) REDO = .TRUE.
         IF (NITER .GE. MAXIT) REDO=.FALSE.
C
C If this star converged, write out the results for it,
C flag it for deletion from the star list (SKIP = .TRUE.), and
C subtract it from the reference copy of the image (DATA).
C
         IF (.NOT. REDO) THEN
            NCONV=NCONV+1
C
C Subtract from the reference copy of the image (DATA).
C
            LX=MAX0(0, MIN0(NCOL, INT(XC(I)-PSFRAD)))+1
            MX=MAX0(1, MIN0(NCOL, INT(XC(I)+PSFRAD)))
            LY=MAX0(0, MIN0(NROW, INT(YC(I)-PSFRAD)))+1
            MY=MAX0(1, MIN0(NROW, INT(YC(I)+PSFRAD)))
            DELTAX=(XC(I)-1.)/XPSF-1.
            DELTAY=(YC(I)-1.)/YPSF-1.
            DO 2195 J=LY,MY
               DY=FLOAT(J)-YC(I)
               DYSQ=DY**2
               DO 2195 K=LX,MX
                  IF (SIGMA(K,J) .GT. -1.E38) THEN
                     DX=FLOAT(K)-XC(I)
                     IF (DX**2+DYSQ .GE. PSFRSQ) GO TO 2195
                     DIFF=MAG(I)*USEPSF(IPSTYP, DX, DY, BRIGHT, PAR, 
     .                    PSF, NPSF, NPAR, NEXP, NFRAC, DELTAX, DELTAY, 
     .                    DVDXC, DVDYC)
                     DATA(K,J)=DATA(K,J)-DIFF
C 
C We must add this star's contribution to the pixel's noise into the
C noise map SIGMA (which, of course, actually contains sigma**2).
C Please note that of the four contributors to the variance
C
C     (1) RONOIS = SIGMA(IX,IY)   (initially, at least)
C     (2) DPOS/PHPADU where DPOS = sum of the stellar profiles plus sky
C     (3) constant x DPOS**2
C     (4) constant x the sum of (the stellar profiles)**2
C
C (1) is already included in SIGMA, and we can calculate this star's
C contribution to (2) and (4) directly.  (3) is tricker.  (3) is
C (sum of sky plus all stellar profiles)**2, and we need to determine
C the contribution of (this stellar profile).  We need to use algebra:
C
C (sum of sky and all OTHER stellar profiles + this stellar profile)**2
C
C = (sum of sky plus all OTHER stellar profiles)**2 +
C  2 x (sum of sky plus OTHER stellar profiles) x (this stellar profile)
C   + (this stellar profile)**2.  
C
C So we must add
C
C 2 x (sum of sky plus OTHER stellar profiles) x (this stellar profile)
C + (this stellar profile)**2
C
C into SIGMA.  That way, later on when the solution has forgotten that
C this star ever existed, so it only adds
C
C (sum of sky plus all OTHER stellar profiles)**2
C
C into the noise expected for the pixel, it will get the noise right.
C Since we have already subtracted this star's contribution to the image
C (array DATA) above, I will use the current DATA as my working
C approximation for (sum of sky plus all OTHER stellar profiles).
C I really shouldn't do this; I should really be calculating each star's
C model-predicted contribution to this pixel, but that would take too
C long.
C
                     IF (DIFF .GT. 0.) THEN
                       DIFF = DIFF/PHPADU +
     .                      PERERR*2.*AMAX1(0., DATA(K,J))*DIFF
     .                      + (PKERR*DIFF)**2
                       SIGMA(K,J)=SIGN(ABS(SIGMA(K,J))+DIFF, SIGMA(K,J))
                     END IF
                  END IF
 2195       CONTINUE                          ! End of loops over pixels
            SHARP=1.4427*PAR(1)*PAR(2)*NUMER(I)/(MAG(I)*PEAK*DENOM(I))
            SHARP=AMIN1( 99.999, AMAX1( SHARP, -99.999 ))
            ERR=1.085736*MAGERR(I)/MAG(I)
            MAG(I)=PSFMAG-1.085736*ALOG(MAG(I))
            IF ((SKY(I) .LT. 9999.99) .AND. (SKY(I) .GT. -999.99)) THEN
               WRITE (3,321) ID(I), XC(I), YC(I), MAG(I), ERR, SKY(I), 
     .                                    FLOAT(NITER), CHI(I), SHARP
  321          FORMAT (I6, 3F9.3, F9.4, F9.3, F9.0, 2F9.3)
            ELSE
               WRITE (3,322) ID(I), XC(I), YC(I), MAG(I), ERR, SKY(I), 
     .                                    FLOAT(NITER), CHI(I), SHARP
  322          FORMAT (I6, 3F9.3, F9.4, F9.2, F9.0, 2F9.3)
            END IF
            SKIP(I)=.TRUE.                    ! Remove from star list
         END IF
 2520 CONTINUE
C
C If there is more than one star remaining in this group, check to see 
C whether any two of them have merged.  This means, find the closest
C pair and see whether they are TOO close.
C
      K = 0
      RSQ = SEPCRT
      IF (NSTR .GT. 1) THEN
         DO 2230 I=ISTAR+1,LSTAR
         IF (SKIP(I)) GO TO 2230
C
         DO 2220 J=ISTAR,I-1
            IF (SKIP(J)) GO TO 2220
            SEP=(XC(I)-XC(J))**2+(YC(I)-YC(J))**2
            IF (SEP .GE. RSQ) GO TO 2220
            RSQ = SEP
            IF (MAGERR(I) .GT. MAGERR(J)) THEN
               K = I
               L = J
            ELSE
               K = J
               L = I
            END IF
 2220    CONTINUE
C
 2230    CONTINUE
C
C No two stars have merged if K still equals zero.
C
         IF (K .LE. 0) GO TO 2260
C
C The K-th star is now the one to eliminate, the J-th the one to
C keep.
C
C Now eliminate the K-th.
C
         IF (MAG(K)*(RSQ-SEPMIN)/(SEPCRT-SEPMIN) .GT. 
     .        WCRIT*MAGERR(K) ) GO TO 2260
C
C Now replace the centroid of the L-th star with the weighted mean of
C the most recent estimates of the centroids of the L-th and K-th
C stars, and the brightness of the L-th with the sum of the brightnesses
C of the L-th and K-th.
C
         ID(L)=MIN0(ID(L),ID(K))
         XC(L)=XC(L)*MAG(L)+XC(K)*MAG(K)
         YC(L)=YC(L)*MAG(L)+YC(K)*MAG(K)
         MAG(L)=MAG(L)+MAG(K)
         XC(L)=XC(L)/MAG(L)
         YC(L)=YC(L)/MAG(L)
C
C Remove the K-th star from the group.
C
         IF (.NOT. SKIP(K)) NDISAP=NDISAP+1
         SKIP(K)=.TRUE.                       ! Remove from star list
         IF (WATCH .GT. 1.5) THEN
            WRITE (6,63) ID(K), XC(K), YC(K), 
     .           PSFMAG-1.085736*ALOG(MAG(K)), ' blended with'
            WRITE (6,63) ID(L), XC(L), YC(L), 
     .           PSFMAG-1.085736*ALOG(MAG(L))
         END IF
C
C Loosen the clamps of every star in the group.
C
         DO I=ISTAR,LSTAR
            DXOLD(I)=0.
            DYOLD(I)=0.
            XCLAMP(I)=AMIN1(0.5*CLMPMX, 2.*XCLAMP(I))
            YCLAMP(I)=AMIN1(0.5*CLMPMX, 2.*YCLAMP(I))
         END DO
      END IF
C
C If the number of iterations completed is less than or equal to 3, 
C perform another iteration no questions asked.
C
 2260 IF (NITER .LE. 3) GO TO 3000
C
C > IF NO STAR HAS BEEN REMOVED FROM THIS GROUP DURING THIS ITERATION <
C
C Check whether any of the stars is too faint (more than 12.5
C magnitudes fainter than the PSF star).  If several stars are
C too faint, delete the faintest one, and set the brightnesses of
C the other faint ones exactly to 12.5 mag below the PSF star.
C That way on the next iteration we will see whether these stars
C want to grow or to disappear.
C
      FAINT=1.0E-5
      IFAINT=0
C
      DO I=ISTAR,LSTAR
         IF (SKIP(I)) GO TO 3000
         IF (MAG(I) .LT. FAINT) THEN
            FAINT=MAG(I)
            IFAINT=I
         END IF
         IF (MAG(I) .LT. 1.E-5) MAG(I)=1.E-5
      END DO
C
C If at least one star is more than 12.5 mag. fainter than the
C PSF, then  IFAINT  is the index of the faintest of them, and FAINT
C is the relative brightness of the faintest of them.
C
      IF (IFAINT .GT. 0) THEN
         IF (.NOT. SKIP(IFAINT)) NDISAP=NDISAP+1
         SKIP(IFAINT)=.TRUE.                     ! Remove from star list
         IF (WATCH .GT. 1.5) WRITE (6,63) ID(IFAINT), XC(IFAINT), 
     .           YC(IFAINT), PSFMAG-1.085736*ALOG(MAG(IFAINT)), 
     .           ' too faint'
C
C Loosen the clamps of every star in the group.
C
         DO I=ISTAR,LSTAR
            DXOLD(I)=0.
            DYOLD(I)=0.
            XCLAMP(I)=AMIN1(0.5*CLMPMX, 2.*XCLAMP(I))
            YCLAMP(I)=AMIN1(0.5*CLMPMX, 2.*YCLAMP(I))
         END DO
      ELSE IF (NITER .GE. 4) THEN
C
C If no star in this group is more than 12.5 mag. fainter than the PSF, 
C then after the fifth iteration delete the least certain star if it 
C is less than a one-sigma detection; after the tenth iteration delete 
C the least certain star if it is less than a 1.5-sigma detection; 
C after the fifteenth iteration delete the least certain star if it is 
C less than a two-sigma detection.
C
         FAINT=1.E38
         IFAINT=0
C
         DO 2550 I=ISTAR,LSTAR
            WT=MAG(I)/MAGERR(I)
            IF (WT .LT. FAINT) THEN
               FAINT=WT
               IFAINT=I
            END IF
 2550    CONTINUE
C
         IF (FAINT .LT. WCRIT) THEN
            IF (.NOT. SKIP(IFAINT)) NDISAP=NDISAP+1
            SKIP(IFAINT)=.TRUE.                  ! Remove from star list
            IF (WATCH .GT. 1.5) WRITE (6,63) ID(IFAINT), XC(IFAINT), 
     .              YC(IFAINT), PSFMAG-1.085736*ALOG(MAG(IFAINT)),
     .              ' error too big: ',
     .              1.085736*SQRT(MAGERR(IFAINT))/MAG(IFAINT)
C
C Loosen the clamps of every star in the group.
C
            DO I=ISTAR,LSTAR
               DXOLD(I)=0.
               DYOLD(I)=0.
               XCLAMP(I)=AMIN1(0.5*CLMPMX, 2.*XCLAMP(I))
               YCLAMP(I)=AMIN1(0.5*CLMPMX, 2.*YCLAMP(I))
            END DO
         END IF
      END IF
 3000 CONTINUE
      IF ((WATCH .GT. 0.5) .AND. (NSTAR .GT. 100)) THEN
         IF (LSTAR/INTRVL .GT. (ISTAR-1)/INTRVL) THEN
            WRITE (LINE,682) NITER, INTRVL*(LSTAR/INTRVL), 
     .           NDISAP, NCONV
            CALL OVRWRT (LINE(1:24), 2)
         END IF
      END IF
      ISTAR=LSTAR+1
      IF (ISTAR .LE. NSTAR) GO TO 2200           ! Go on to next group
C
C We've gone through the entire star list for this iteration.
C
C Find the last star in the list that still needs more work 
C (SKIP=.FALSE.).  If there aren't any, then we're done!  If there is
C one, point NSTAR at it.
C
      ISTAR=0
 3005 CONTINUE
      IF (SKIP(NSTAR)) THEN
         NSTAR=NSTAR-1
         IF (NSTAR .GT. 0) THEN
            GO TO 3005
         ELSE
            WRITE (LINE,682) NITER, NSTAR, NDISAP, NCONV
            CALL OVRWRT (LINE(1:28), 3)
            GO TO 9000                       ! None left!! We can go now
         END IF
      END IF
C
C Well, at this point we've established that there are at least some
C stars in the list that still need more work.  Remove any stars with
C which we are finished (SKIP = .TRUE.) by overwriting the first such
C star with the star at position NSTAR in the list (we have just
C confirmed that this is the last star in the list that still needs
C work).  Then decrement NSTAR and go back up to 3005 to find the new
C last star in the list that still needs work.
C
 3010 ISTAR=ISTAR+1
      IF (ISTAR .GE. NSTAR) THEN
c         WRITE (LINE,682) NITER, NSTAR, NDISAP, NCONV
  682    FORMAT (4I6, '  <<')
c         CALL OVRWRT (LINE(1:28), 3)
         GO TO 2100        ! The list is ready.  Start another iteration
      END IF
C
      IF (SKIP(ISTAR)) THEN
         ID(ISTAR)=ID(NSTAR)
         XC(ISTAR)=XC(NSTAR)
         YC(ISTAR)=YC(NSTAR)
         MAG(ISTAR)=MAG(NSTAR)
         SKY(ISTAR)=SKY(NSTAR)
         SUMWT(ISTAR)=SUMWT(NSTAR)
         DXOLD(ISTAR)=DXOLD(NSTAR)
         DYOLD(ISTAR)=DYOLD(NSTAR)
         XCLAMP(ISTAR)=XCLAMP(NSTAR)
         YCLAMP(ISTAR)=YCLAMP(NSTAR)
         SKIP(ISTAR)=.FALSE.
         NSTAR=NSTAR-1
         GO TO 3005
      ELSE
         GO TO 3010
      END IF
C
C-----------------------------------------------------------------------
C
C Normal return.
C
 9000 CONTINUE
	close(1)
	istat=9000
      RETURN
      END!

c *********************************************************************
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  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
         write(*,*)'Stack overflow in QUICK.  Increase MAXSTK.'
      END IF
      STKLO(NSTAK)=LIMLO
      STKHI(NSTAK)=LO-1
      LIMLO=LO+1
      GO TO 100
  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
         write(*,*)'Stack overflow in QUICK.  Increase MAXSTK.'
      END IF
      STKLO(NSTAK)=LO+1
      STKHI(NSTAK)=LIMHI
      LIMHI=LO-1
      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
      SUBROUTINE RECTFY (X, NSTAR, INDEX, HOLD)
      IMPLICIT NONE
C
      REAL X(*), HOLD(*)
      INTEGER INDEX(*)
C
      INTEGER I, NSTAR
C
      DO I=1,NSTAR
         HOLD(I)=X(I)
      END DO
      DO I=1,NSTAR
         X(I)=HOLD(INDEX(I))
      END DO
      RETURN
      END!
C
      SUBROUTINE RECTFYI (X, NSTAR, INDEX, HOLD)
      IMPLICIT NONE
C
      REAL HOLD(*)
      INTEGER INDEX(*),X(*)
C
      INTEGER I, NSTAR
C
      DO I=1,NSTAR
         HOLD(I)=X(I)
      END DO
      DO I=1,NSTAR
         X(I)=HOLD(INDEX(I))
      END DO
      RETURN
      END!
C
C#######################################################################
C
      REAL FUNCTION PCTILE(DATUM,N,NPCT)
C
C=======================================================================
C
C This is a modification of a quick-sorting algorithm, which is intended
C to take in a vector of numbers, and return the value of the PCT-th
C percentile in that vector:
C
C    DATUM (input real vector)     containing real data.
C    N     (input integer)         number of elements in DATUM.
C    NPCT  (input integer)         element of sorted vector whose value 
C                                  is desired.
C    PCTILE (output real)          the value of the NPCT-th element
C                                  in the sorted vector DATUM.
C
C-----------------------------------------------------------------------
C
C The quick-sorting algorithm was suggested by the discussion on pages 
C 114-119 of THE ART OF COMPUTER PROGRAMMING, Vol. 3, SORTING AND 
C SEARCHING, by D.E. Knuth, which was referenced in Don Wells' 
C subroutine QUIK.  This is my own attempt at encoding a quicksort-- 
C                                                             PBS.
C
C The array DATUM contains randomly ordered data. 
C
      IMPLICIT NONE
      REAL DATUM(*)
C
      INTEGER MIN0, MAX0
C
      REAL DKEY
      INTEGER LO, HI, N, NPCT, LIMLO, LIMHI, NTARG
C
C Which element of the sorted array will we be interested in?
C
      NTARG=MAX0(1,MIN0(N,NPCT))
C
C Initialize the pointers.
C
      LIMLO=1
      LIMHI=N
C
  100 DKEY=DATUM(LIMLO)
C     TYPE *,'LOW=',LIMLO,' HIGH=',LIMHI,' KEY=',DKEY
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
C If LO equals HI, we have tested all the elements in the current search
C interval.
C
      IF(LO.EQ.HI)GO TO 200
      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)
      LO=LO+1
  110 CONTINUE
      IF(LO.EQ.HI)GO TO 200
      IF(DATUM(LO).GE.DKEY)GO TO 119
      LO=LO+1
      GO TO 110
C
  119 DATUM(HI)=DATUM(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
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).  If LO = NTARG, then DATUM(LO) is
C the value we are looking for.  If NTARG < LO, then we want to sort the
C values of DATUM from LIMLO to LO-1, inclusive, whereas if NTARG > LO,
C then we want to sort the values of DATUM from LO+1 to LIMHI, 
C inclusive.
C
C     TYPE *,'NTARG=',NTARG,' LO=',LO
      IF(NTARG-LO)300,900,400
  300 LIMHI=LO-1
      GO TO 100
  400 LIMLO=LO+1
      GO TO 100
  900 PCTILE=DATUM(LO)
      RETURN
      END!
C
C#######################################################################
C
      REAL FUNCTION SMLLST (DATUM, N, M)
      IMPLICIT NONE
      REAL DATUM(*)
      INTEGER I, J, N, M, LEAST
      DO J=1,M
         SMLLST = DATUM(J)
         LEAST = J
         DO I=J,N
            IF (DATUM(I) .LT. SMLLST) THEN
               SMLLST = DATUM(I)
               LEAST = I
            END IF
         END DO
         DATUM(LEAST) = DATUM(J)
         DATUM(J) = SMLLST
      END DO
      RETURN
      END!
C
C#######################################################################
C
      REAL FUNCTION BIGGST (DATUM, N, M)
      IMPLICIT NONE
      REAL DATUM(*)
      INTEGER I, J, N, M, MOST
      DO J=N,N-M+1,-1
         MOST = J
         BIGGST = DATUM(MOST)
         DO I=1,J
            IF (DATUM(I) .GT. BIGGST) THEN
               BIGGST = DATUM(I)
               MOST = I
            END IF
         END DO
         DATUM(MOST) = DATUM(J)
         DATUM(J) = BIGGST
      END DO
      RETURN
      END!
C
C#######################################################################
C
      INTEGER FUNCTION  RDPSF  (id, IPSTYP, PAR, MAXPAR, NPAR,
     .     PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, 
     .     PSFMAG, BRIGHT, XPSF, YPSF)
C
C Read in the point-spread function
C
      IMPLICIT NONE
      integer id
      INTEGER MAXPSF, MAXTYP, MAXPAR, MAXEXP
      PARAMETER (MAXTYP=6)
C
      REAL PAR(MAXPAR), PSF(MAXPSF,MAXPSF,MAXEXP)
C
      CHARACTER*8 LABEL, CHECK
      REAL PSFMAG, BRIGHT, XPSF, YPSF
      INTEGER I, J, K, IPSTYP, NPSF, NPAR, NEXP, NFRAC, ISTAT
      INTEGER NTERM, NPARAM
C
      READ(id,302,IOSTAT=ISTAT)LABEL,NPSF,NPAR,NEXP,NFRAC,PSFMAG, 
     .     BRIGHT, XPSF, YPSF
  302 FORMAT (1X, A8, 4I5, F9.3, F15.3, 2F9.1)
      IF (ISTAT .NE. 0) THEN
        write(*,*)'Error reading PSF.'
	close(id)
         RDPSF = -1
         RETURN
      END IF
      DO IPSTYP=1,MAXTYP
         I = NPARAM(IPSTYP, 1., CHECK, PAR, MAXPAR)
         IF ((LABEL .EQ. CHECK) .AND. (I .EQ. NPAR)) GO TO 1100
      END DO
      write(*,*)'Inappropriate PSF: '//LABEL
      close(id)
      RDPSF = -1
      RETURN
C
 1100 READ (id,301,IOSTAT=ISTAT) (PAR(I), I=1,NPAR)
  301 FORMAT (1X, 6E13.6)
      IF (ISTAT .NE. 0) THEN
         write(*,*)'Error reading PSF.'
         RDPSF = -1
         close(id)
         RETURN
      END IF
      NTERM = NEXP+NFRAC
      IF (NTERM .GE. 1) THEN
         DO K=1,NTERM
            READ(id,311,IOSTAT=ISTAT)((PSF(I,J,K),I=1,NPSF),J=1,NPSF)
  311       FORMAT (1X, 6E13.6)
            IF (ISTAT .NE. 0) THEN
               close(id)
               RDPSF = -1
               RETURN
            END IF
         END DO
      END IF
      close(id)
      RDPSF = 0
      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
      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
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 DAORAN (IDUM)
C
C RAN2 from Numerical Recipes
C
      INTEGER IR(97)
      DATA IFF /0/
      DATA M/714025/, IA/1366/, IC/150889/, RM/1.400511E-6/
 1000 CONTINUE
      IF ((IDUM .LT. 0) .OR. (IFF .EQ. 0)) THEN
         IFF = 1
         IDUM = MOD(IABS(IC-IDUM), M)
         DO J=1,97
            IDUM = MOD(IA*IDUM+IC,M)
            IR(J) = IDUM
         END DO
         IDUM = MOD(IA*IDUM+IC,M)
         IY = IDUM
      END IF
      J = 1+(97*IY)/M
      IF ((J .GT. 97) .OR. (J .LT. 1)) PAUSE
      IY = IR(J)
      DAORAN = IY*RM
      IDUM = MOD(IA*IDUM+IC,M)
      IR(J) = IDUM
      IF (DAORAN .LE. 0.) GO TO 1000              ! Stetson's modification
C
      RETURN
      END!
C
C######################################################################
C
      REAL  FUNCTION  NRML  (RANNUM)
C
C Convert a uniform probability distribution to a Gaussian distribution
C with mean zero and standard deviation unity.
C
      IMPLICIT NONE
C
      REAL SQRT, ALOG
C
      REAL P, RANNUM, SIGN, T
 1000 P=RANNUM
      SIGN=-1.
      IF (P .GT. 0.5) THEN
         P=P-0.5
         SIGN=1.
      ELSE IF (P .LE. 0.) THEN
         NRML = -1.E20
         RETURN
      END IF
      T=SQRT(ALOG(1/P**2))
      NRML=SIGN*(T- (2.30753+.27061*T) / (1.+T*(.99229+T*.04481)) )
      RETURN
      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#######################################################################
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
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*************************************
      SUBROUTINE STRIP (ID, X, Y, MAG, SKY, SKIP, MAXSTR,
     .     NSTAR, NDISAP, RADIUS, INDEX, HOLD)
      REAL X(MAXSTR), Y(MAXSTR), MAG(MAXSTR), SKY(MAXSTR)
      REAL HOLD(MAXSTR)
      INTEGER ID(MAXSTR), INDEX(MAXSTR)
      LOGICAL SKIP(MAXSTR)
C
      NDISAP = 0
      RADSQ = RADIUS**2
      IF (NSTAR .LE. 1) RETURN
      DO I=1,NSTAR
         SKIP(I) = .FALSE.
      END DO
C
      CALL QUICK (Y, NSTAR, INDEX)
      CALL RECTFYI (ID, NSTAR, INDEX, HOLD)
      CALL RECTFY (X, NSTAR, INDEX, HOLD)
      CALL RECTFY (MAG, NSTAR, INDEX, HOLD)
      CALL RECTFY (SKY, NSTAR, INDEX, HOLD)
C
      DO 1200 I=1,NSTAR-1
         IF (SKIP(I)) GO TO 1200
         DO 1100 J=I+1,NSTAR
            IF (SKIP(J)) GO TO 1100
            DY = Y(J)-Y(I)
            IF (DY .GT. RADIUS) GO TO 1200
            DX = X(J)-X(I)
            IF (ABS(DX) .GT. RADIUS) GO TO 1100
            IF (DX**2+DY**2 .GT. RADSQ) GO TO 1100
            IF (MAG(J) .LE. MAG(I)) THEN
               SKIP(J) = .TRUE.
               GO TO 1100
            ELSE
               SKIP(I) = .TRUE.
               GO TO 1200
            END IF
 1100    CONTINUE
 1200 CONTINUE
C
      ISTAR = 0
 2000 CONTINUE
      IF (SKIP(NSTAR)) THEN
         NSTAR = NSTAR-1
         NDISAP = NDISAP+1
         GO TO 2000
      END IF
 2100 ISTAR = ISTAR+1
      IF (ISTAR .GE. NSTAR) RETURN
      IF (SKIP(ISTAR)) THEN
         ID(ISTAR) = ID(NSTAR)
         X(ISTAR) = X(NSTAR)
         Y(ISTAR) = Y(NSTAR)
         MAG(ISTAR) = MAG(NSTAR)
         SKY(ISTAR) = SKY(NSTAR)
         SKIP(ISTAR) = .FALSE.
         NSTAR = NSTAR-1
         NDISAP = NDISAP+1
         GO TO 2000
      ELSE
         GO TO 2100
      END IF
      END
C
C#######################################################################
C
      SUBROUTINE  REGRP (ID, X, Y, MAG, SKY, CHI, DXOLD, DYOLD, 
     .     XCLAMP, YCLAMP, MAXSTR, NSTAR, FITRAD, LAST, INDEX, HOLD)
C
C=======================================================================
C
C This subroutine accepts a list of stellar coordinates, and
C associates the stars into natural groups based on a critical 
C separation:  stars within one critical separation of each other are 
C put into the same group; no star is within one critical separation of 
C any star outside its group.  
C
C             OFFICIAL DAO VERSION:  1985 August 15
C
C======================================================================
C
      REAL X(MAXSTR), Y(MAXSTR), MAG(MAXSTR), SKY(MAXSTR)
      REAL CHI(MAXSTR), DXOLD(MAXSTR), DYOLD(MAXSTR)
      REAL XCLAMP(MAXSTR), YCLAMP(MAXSTR), HOLD(MAXSTR)
      INTEGER ID(MAXSTR), INDEX(MAXSTR)
      LOGICAL LAST(MAXSTR)
C
C-----------------------------------------------------------------------
C
C SECTION 1
C
C Get set up.
C
      IF (NSTAR .LE. 1) RETURN
      CRIT=2.*FITRAD
      CRITSQ=CRIT**2
C
C Check that the stars are sorted by y-coordinate on input.
C
      INDEX(1) = 1
      DO I=2,NSTAR
         IF (Y(I) .GE. Y(I-1)) THEN
            INDEX(I) = I
         ELSE
            CALL QUICK (Y, NSTAR, INDEX)
            CALL RECTFY (X, NSTAR, INDEX, HOLD)
            GO TO 900
         END IF
      END DO
C
  900 ITEST=0
      ITOP=2
C
C The stars are currently in a stack NSTAR stars long.  The variable 
C ITEST will point to the position in the stack occupied by the star 
C which is currently the center of a circle of the critical radius, 
C within which we are looking for other stars; this also starts out 
C with a value of 1.  ITOP points to the top position in the stack of 
C the stars which have not yet been assigned to groups; this starts 
C out with the value 2.  Each time through, the 
C program goes down through the stack from ITOP and looks for stars 
C within the critical distance from the star at stack position 
C ITEST.  When such a star is found, it changes places in the stack 
C with the star at ITOP and ITOP is incremented by one.  When the 
C search reaches a star J such that Y(J)-Y(ITEST) > CRIT it is known
C that no further stars will be found within the critical distance, the
C pointer ITEST is incremented by one, and the search proceeds again 
C from the new value of ITOP.  If the pointer ITEST catches up 
C with the pointer ITOP, that means that the group currently being 
C built up is complete.  Then a new group is started beginning with 
C the star at the current position ITEST, ( = the instantaneous 
C value of ITOP), ITOP is incremented by 1, and the next group is built
C up as before.
C
 2100 ITEST=ITEST+1
      LAST(ITEST)=.FALSE.
      IF (ITEST .EQ. ITOP) THEN
C
C ITEST has reached ITOP; no other unassigned stars are within a 
C critical separation of any member of the current group.  The group is 
C therefore complete.  Signify this by setting LAST(ITEST-1)=.TRUE.
C (ITEST = the current value of ITOP), and then increment the value of 
C ITOP by one.
C
         J=ITEST-1
         IF (J .GT. 0) LAST(J)=.TRUE.
         ITOP=ITOP+1                              ! Increment ITOP
C
C If ITOP is greater than NSTAR at this point, then we are finished 
C (the last group contains one star).  Otherwise, on with the search.
C
         IF (ITOP .GT. NSTAR) THEN
            LAST(ITEST)=.TRUE.
            GO TO 3000
         END IF
      END IF
C
C Now go through the list of unassigned stars, occupying positions ITOP
C through NSTAR in the stack, to look for stars within the critical 
C distance of the star at position ITEST in the stack.  If one is found,
C move it up to stack position ITOP and increment ITOP by one.
C 
      XTEST=X(ITEST)
      YTEST=Y(ITEST)
      J=ITOP
      DO 2120 I=J,NSTAR
         DY=Y(I)-YTEST
         IF (DY .GT. CRIT) GO TO 2100
         DX=X(I)-XTEST
         IF (ABS(DX) .GT. CRIT) GO TO 2120
         IF (DX**2+DY**2 .GT. CRITSQ) GO TO 2120
C
C This star is within the critical distance of the star at stack 
C position ITEST.  Therefore it should be added to the current group by
C moving it up to position ITOP in the stack, where the pointer ITEST 
C may eventually reach it.
C
         CALL ASWAP (MAXSTR, ITOP, I, X, Y, INDEX)
C
C Now increment ITOP by 1 to point at the topmost unassigned star in the
C stack.
C
         ITOP=ITOP+1
C
C If ITOP is greater than NSTAR, then all stars have been assigned to
C groups, and we are finished.  
C
         IF (ITOP .GT. NSTAR) THEN
            DO K=ITEST,NSTAR-1
               LAST(K)=.FALSE.
            END DO
            LAST(NSTAR)=.TRUE.
            GO TO 3000
         END IF
 2120 CONTINUE
      GO TO 2100
C
 3000 CONTINUE
C
C Rectify the remaining quantities.
C
      CALL RECTFYI (ID, NSTAR, INDEX, HOLD)
      CALL RECTFY (MAG, NSTAR, INDEX, HOLD)
      CALL RECTFY (SKY, NSTAR, INDEX, HOLD)
      CALL RECTFY (CHI, NSTAR, INDEX, HOLD)
      CALL RECTFY (DXOLD, NSTAR, INDEX, HOLD)
      CALL RECTFY (DYOLD, NSTAR, INDEX, HOLD)
      CALL RECTFY (XCLAMP, NSTAR, INDEX, HOLD)
      CALL RECTFY (YCLAMP, NSTAR, INDEX, HOLD)
C
      RETURN
      END
C
C#######################################################################
C
      SUBROUTINE  ASWAP (MAXSTR, I, J, X, Y, INDEX)
C
C=======================================================================
C
C Make the I-th and J-th stars in the stack change places (J > I),
C without otherwise altering the order of the stars.  The other 
C arguments are self-evident.
C
C=======================================================================
C
      IMPLICIT NONE
      INTEGER MAXSTR
C
      REAL X(MAXSTR), Y(MAXSTR)
      INTEGER INDEX(MAXSTR)
C
      REAL XHOLD, YHOLD
      INTEGER I, J, K, L, IHOLD
C
C-----------------------------------------------------------------------
C
c     call ovrwrt ('ASWAP', 2)
      XHOLD=X(J)
      YHOLD=Y(J)
      IHOLD=INDEX(J)
      DO K=J,I+1,-1
         L=K-1
         X(K)=X(L)
         Y(K)=Y(L)
         INDEX(K)=INDEX(L)
      END DO
      X(I)=XHOLD
      Y(I)=YHOLD
      INDEX(I)=IHOLD
      RETURN
      END
 
        function indexpos(head,f1,n)
        character*80 head(*),f1*8
        do 10 indexpos=1,n
10      if(head(indexpos)(1:8).eq.f1)return
        end
