	parameter (n2048=2048*3)
        PARAMETER (MAXEXP=10, MAXPSF=207, MAXPAR=6)

        real par(MAXPAR)
   
	real         a_b(n2048*n2048)
	character*60  f1,f2,f3
	character*80  head(36)
        INTEGER NCOL, NROW
	logical logi
        if(iargc().lt.1)then
        write(*,*)
        write(*,*)
        write(*,*)'     **********   DAO_SUBSTAR  2004.5 ***********'
        write(*,*)
        write(*,*)
	write(*,*)'     Usage:  dao_substar 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 !'
 
        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)//'lst'			! 3 read  *.lst
        open(3,file=f2,status='old')

	call readfits(f1,a_b,ncol,nrow)
C
	wtach=-2.
        call SUBSTR(PAR,MAXPAR,MAXPSF,MAXEXP,a_b,NCOL,NROW,WATCH)
C
	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)
	call swap4(a_b,ncol*nrow*4)
	call writefits(f3,a_b,ncol*nrow*4,0)
	write(*,*)'  ok!'
	
      END!
C


      SUBROUTINE  SUBSTR  (PAR, MAXPAR, MAXPSF, MAXEXP, 
     .     F, NCOL, NROW, WATCH)
C
C=======================================================================
C
C This subroutine scales and shifts the point-spread function according
C to each star's magnitude and centroid, and subtracts the resulting
C profile from a copy of the original picture.
C
C             OFFICIAL DAO VERSION:  1991 April 18
C
C Arguments
C
C  WATCH (INPUT) governs whether information relating to the progress 
C        of the reductions is to be typed on the terminal screen
C        during execution.
C
C=======================================================================
C
      IMPLICIT NONE
      INTEGER MAXBOX, MAXPSF, MAXPAR, MAXEXC, MAXEXP
      PARAMETER  (MAXBOX=69, MAXEXC=300)
C
C Parameters
C
C MAXBOX is the side of the square subarray containing the largest
C        (circular) PSF that can be subtracted from the picture.
C
C MAXPSF is the largest permissible number of elements on a side of the
C        (square) look-up table for the point-spread function.
C
C        MAXBOX = (MAXPSF-7)/2.
C
      INTEGER NCOL, NROW
      REAL PAR(MAXPAR)
      REAL F(NCOL,NROW), PSF(MAXPSF,MAXPSF,MAXEXP)
      INTEGER IEXC(MAXEXC)
C
      REAL USEPSF
      INTEGER RDPSF, MAX0, MIN0, INT
C
      CHARACTER*80 LINE
      CHARACTER*30 SUBPIC
      CHARACTER*30 COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL, SWITCH
      CHARACTER CASE*5, ANSWER*1
      REAL LOBAD, X, Y, STRMAG, SKY, DX, DY, DYSQ, SCALE, COL, ROW
      REAL DELTAX, DELTAY, DVDXC, DVDYC, XPSF, YPSF, PSFRAD, DUM
      REAL HIBAD, THRESH, AP1, PHPADU, READNS, FRAD, WATCH, PSFMAG
      REAL BRIGHT, PSFRSQ
      INTEGER I, J, L, NL, IDUM, NEXC, LX, LY, NX, NY, ISTAR, ID
      INTEGER ISTAT, IPSTYP, NPSF, NPAR, NEXP, NFRAC
C
      COMMON /FILNAM/ COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL
C
C-----------------------------------------------------------------------
C
C SECTION 1
C
C Get file names and set up the needed numerical constants.
C
c      CALL TBLANK                                   ! Type a blank line
c  950 CALL GETNAM ('File with the PSF:', PSFFIL)
c      IF ((PSFFIL .EQ. 'END OF FILE') .OR.
c     .     (PSFFIL .EQ. 'GIVE UP')) THEN
c         PSFFIL = ' '
c         RETURN
c      END IF
      ISTAT = RDPSF (1, IPSTYP, PAR, MAXPAR, NPAR,
     .     PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, 
     .     PSFMAG, BRIGHT, XPSF, YPSF)
c      IF (ISTAT .LT. 0) THEN
c         PSFFIL = 'GIVE UP'
c         GO TO 950
c      END IF
C
      PSFRAD=(REAL(NPSF-1)/2. - 1.)/2.
      PSFRSQ=PSFRAD**2
      COL=FLOAT(NCOL)
      ROW=FLOAT(NROW)
C
c      CALL CLFILE (2)
	close(11)
C
c 1015 CALL GETNAM ('File with photometry:', PROFIL)
c      IF ((PROFIL .EQ. 'END OF FILE') .OR.
c     .     (PROFIL .EQ. 'GIVE UP')) THEN
c         PROFIL = ' '
c         RETURN
c      END IF
c      CALL INFILE (2, PROFIL, ISTAT)
c      IF (ISTAT .NE. 0) THEN
c         CALL STUPID ('Error opening input file '//PROFIL)
c         PROFIL = 'GIVE UP'
c         GO TO 1015
c      END IF
C
c      CALL GETYN ('Do you have stars to leave in?', ANSWER)
c      IF (ANSWER .EQ. 'Y') THEN
c         SUBPIC = SWITCH(PROFIL, CASE('.lst'))
c 1025    CALL GETNAM ('File with star list:', SUBPIC)
c         IF ((SUBPIC .EQ. 'END OF FILE') .OR. (SUBPIC .EQ.
c     .        'GIVE UP')) THEN
c            CALL CLFILE (2)
c            RETURN
c         END IF
c         CALL INFILE (3, SUBPIC, ISTAT)
c         IF (ISTAT .NE. 0) THEN
c            CALL STUPID ('Error opening file '//SUBPIC)
c            SUBPIC = 'GIVE UP'
c            GO TO 1025
c         END IF
c
c         CALL CHECK (3, NL)
         ANSWER = 'Y'
         NEXC = 0
 1035    NEXC = NEXC+1
	 if(NEXC.GT.MAXEXC)goto 1036
	 read(3,*,end=1036)IEXC(NEXC)
         goto 1035
 1036    close(3)
	 NEXC=NEXC-1
	 write(*,*)NEXC

c         NEXC = 0
c 1035    NEXC = NEXC+1
c 1036    CALL RDSTAR (3, NL, IEXC(NEXC), DUM, DUM, DUM, DUM)
c         IF (IEXC(NEXC) .EQ. 0) GO TO 1036
c         IF (IEXC(NEXC) .LT. 0) THEN
c            NEXC = NEXC-1
c         ELSE
c            IF (NEXC .LT. MAXEXC) GO TO 1035
c         END IF
c	 close(3)
c                write(*,*)NEXC,MAXEXC
c         IF (NEXC .LE. 0) ANSWER = 'N'

c         CALL CLFILE (3)
c      END IF
C
c      CALL RDHEAD (2, NL, IDUM, IDUM, LOBAD, HIBAD, THRESH, AP1, 
c     .     PHPADU, READNS, FRAD)
	read(2,'(a)')i
        read(2,*)NL,IDUM,IDUM,LOBAD,HIBAD,THRESH,AP1,PHPADU,READNS,FRAD
C
c      SUBPIC=SWITCH(PROFIL, CASE('s'))
c      CALL GETNAM ('Name for subtracted image:', SUBPIC)
c      IF (SUBPIC .EQ. 'END OF FILE') GO TO 9010    ! CTRL-Z was entered
C
C Copy the input picture verbatim into the output picture.
C
C     CALL COPPIC (SUBPIC, F, NCOL, NROW, ISTAT)
c      CALL COPPIC (SUBPIC, ISTAT)
c      IF (ISTAT .NE. 0) THEN
c         CALL STUPID ('Error creating output image.')
c         CALL CLPIC ('COPY')
c         CALL CLFILE (2)
c         RETURN
c      END IF
C
      IF (WATCH .GT. 0.5) THEN
c         CALL TBLANK
	  write(*,*)
         CALL OVRWRT (' Star', 1)
      END IF
C
C Read the entire image in from the copy.
C
c      LX = 1
c      LY = 1
c      NX = NCOL
c      NY = NROW
c      CALL RDARAY ('COPY', LX, LY, NX, NY, NCOL, F, ISTAT)
C
C-----------------------------------------------------------------------
C
C SECTION 2
C
C Loop over stars.
C
      ISTAR=0
 2000 ISTAR=ISTAR+1
 2010 CALL RDSTAR (2, NL, ID, X, Y, STRMAG, SKY)
      IF (ID .LT. 0) GO TO 9000                ! End-of-file encountered
      IF (ID .EQ. 0) GO TO 2010                ! Ignore a blank line
      IF ((1.-X .GT. PSFRAD) .OR. (X-COL .GT. PSFRAD) .OR.
     .    (1.-Y .GT. PSFRAD) .OR. (Y-ROW .GT. PSFRAD)) GO TO 2000
      IF (ANSWER .EQ. 'Y') THEN
         L = NEXC
         DO I=1,L
            IF (ID .EQ. IEXC(I)) THEN
               IEXC(I) = IEXC(NEXC)
               NEXC = NEXC-1
               IF (NEXC .EQ. 0) ANSWER = 'N'
               GO TO 2010
            END IF
         END DO
      END IF
      IF (STRMAG .GE. 99.) GO TO 2000          ! Ignore a bad star
      IF ((WATCH .GT. 0.5) .AND. (MOD(ISTAR,20) .EQ. 0)) THEN
         WRITE (LINE,620) ISTAR
  620    FORMAT (I5)
         CALL OVRWRT (LINE(1:5), 2)
      END IF
      DELTAX=(X-1.)/XPSF-1.
      DELTAY=(Y-1.)/YPSF-1.
      LX=MAX0(1, INT(X-PSFRAD)+1)
      LY=MAX0(1, INT(Y-PSFRAD)+1)
      NX=MIN0(NCOL, INT(X+PSFRAD))
      NY=MIN0(NROW, INT(Y+PSFRAD))
      SCALE=10.**(0.4*(PSFMAG-STRMAG))
C
C Subtract the shifted scaled PSF
C
      DO 2030 J=LY,NY
         DY=FLOAT(J)-Y
         DYSQ=DY**2
         DO 2020 I=LX,NX
            IF (F(I,J) .GT. HIBAD) GO TO 2020
            DX=FLOAT(I)-X
            IF (DX**2+DYSQ .GE. PSFRSQ) THEN
               IF (DX .GT. 0.) GO TO 2030
            ELSE
               F(I,J)=F(I,J)-SCALE*USEPSF(IPSTYP, DX, DY, BRIGHT, PAR, 
     .              PSF, NPSF, NPAR, NEXP, NFRAC, DELTAX, DELTAY, 
     .              DVDXC, DVDYC)
            END IF
 2020    CONTINUE
 2030 CONTINUE
      GO TO 2000                              ! Go to next star
C
C-----------------------------------------------------------------------
C
C Normal return.
C
 9000 CONTINUE
      IF (WATCH .GT. 0.5) THEN
         WRITE (LINE,620) ISTAR-1
         CALL OVRWRT (LINE(1:5), 2)
      END IF
      IF (WATCH .GT. 0.5) CALL OVRWRT (' ', 4)
C
C Write the modified image back into the copy.
C
c      LX = 1
c      LY = 1
c      NX = NCOL
c      NY = NROW
c      CALL WRARAY ('COPY', LX, LY, NX, NY, NCOL, F, ISTAT)
c      CALL CLPIC ('COPY')
c      IF (WATCH .GT. -1.5) write(*,*)' done. '
 9010 close(2)
      RETURN
C
      END!
C
C This file contains subroutines that are not I/O related, but rather
C involve arithmetic or character operations that my be somewhat 
C machine-specific-- it may be necessary or desirable to make some
C changes to optimize the code to run on your computer.
C
C              OFFICIAL DAO VERSION:  1991 April 18
C
C***********************************************************************
C
C Current contents:
C
C   INVERS  inverts a square matrix.
C
C     VMUL  multiplies a square matrix (on the left) by a column vector
C           (on the right) yielding a new column vector.
C
C    QUICK  does a quicksort on a vector of data.
C
C     SHOW  produces what is effectively an eleven-level gray scale plot
C           of a rectangular data array on the terminal.
C
C***********************************************************************
C
      SUBROUTINE  INVERS (A, MAX, N, IFLAG)
C
C Although it seems counter-intuitive, the tests that I have run
C so far suggest that the 180 x 180 matrices that NSTAR needs can
C be inverted with sufficient accuracy if the elements are REAL
C rather than DOUBLE PRECISION
C
C Arguments
C 
C     A (INPUT/OUTPUT) is a square matrix of dimension N.  The inverse 
C       of the input matrix A is returned in A.
C
C   MAX (INPUT) is the size assigned to the matrix A in the calling 
C       routine.  It's needed for the dimension statement below.
C
C IFLAG (OUTPUT) is an error flag.  IFLAG = 1 if the matrix could not
C       be inverted; IFLAG = 0 if it could.
C
      IMPLICIT NONE
      INTEGER MAX
      REAL A(MAX,MAX)
C
      INTEGER N, IFLAG, I, J, K
C
C-----------------------------------------------------------------------
C
      IFLAG=0
      I=1
  300 IF(A(I,I).EQ.0.0E0)GO TO 9100
      A(I,I)=1.0E0/A(I,I)
      J=1
  301 IF(J.EQ.I)GO TO 304
      A(J,I)=-A(J,I)*A(I,I)
      K=1
  302 IF(K.EQ.I)GO TO 303
      A(J,K)=A(J,K)+A(J,I)*A(I,K)
  303 IF(K.EQ.N)GO TO 304
      K=K+1
      GO TO 302
  304 IF(J.EQ.N)GO TO 305
      J=J+1
      GO TO 301
  305 K=1
  306 IF(K.EQ.I)GO TO 307
      A(I,K)=A(I,K)*A(I,I)
  307 IF(K.EQ.N)GO TO 308
      K=K+1
      GO TO 306
  308 IF(I.EQ.N)RETURN                                   ! Normal return
      I=I+1
      GO TO 300
C
C-----------------------------------------------------------------------
C
C Error:  zero on the diagonal.
C
 9100 IFLAG=1
      RETURN
C
      END!
C
C#######################################################################
C
      SUBROUTINE  VMUL (A, MAX, N, V, X)
C
C Multiply a matrix by a vector:
C
C                    A * V = X
C
C Arguments
C
C    A(column,row)  (INPUT) is a square matrix of dimension N.
C
C              MAX  (INPUT) is the size assigned to the array in the 
C                           calling routine.
C
C           V(row)  (INPUT) is a column vector of dimension N.
C
C           X(row) (OUTPUT) is a column vector of dimension N.
C
      IMPLICIT NONE
      INTEGER MAX
      REAL A(MAX,MAX), V(MAX)
      REAL X(MAX)
C
      DOUBLE PRECISION DBLE
      REAL SNGL
C
      DOUBLE PRECISION SUM
      INTEGER N, I, J
C
C-----------------------------------------------------------------------
C
      I=1
  200 SUM=0.0D0
      J=1
  201 SUM=SUM+DBLE(A(J,I))*DBLE(V(J))
      IF (J .EQ. N) GO TO 203
      J=J+1
      GO TO 201
  203 X(I)=SNGL(SUM)
      IF (I .EQ. N) RETURN                               ! Normal return
      I=I+1
      GO TO 200
      END!
C
C#######################################################################
C
      SUBROUTINE  DINVRS (A, MAX, N, IFLAG)
      IMPLICIT NONE
      INTEGER MAX
      DOUBLE PRECISION A(MAX,MAX)
C
      INTEGER I, J, K, N, IFLAG
C
C-----------------------------------------------------------------------
C
      IFLAG=0
      I=1
  300 IF(A(I,I).EQ.0.0D0)GO TO 9100
      A(I,I)=1.0D0/A(I,I)
      J=1
  301 IF(J.EQ.I)GO TO 304
      A(J,I)=-A(J,I)*A(I,I)
      K=1
  302 IF(K.EQ.I)GO TO 303
      A(J,K)=A(J,K)+A(J,I)*A(I,K)
  303 IF(K.EQ.N)GO TO 304
      K=K+1
      GO TO 302
  304 IF(J.EQ.N)GO TO 305
      J=J+1
      GO TO 301
  305 K=1
  306 IF(K.EQ.I)GO TO 307
      A(I,K)=A(I,K)*A(I,I)
  307 IF(K.EQ.N)GO TO 308
      K=K+1
      GO TO 306
  308 IF(I.EQ.N)RETURN                                   ! Normal return
      I=I+1
      GO TO 300
C
C-----------------------------------------------------------------------
C
C Error:  zero on the diagonal.
C
 9100 IFLAG=1
      RETURN
C
      END!
C
C#######################################################################
C
      SUBROUTINE  DVMUL (A, MAX, N, V, X)
      IMPLICIT NONE
      INTEGER MAX
      DOUBLE PRECISION A(MAX,MAX), V(MAX), X(MAX)
C
      DOUBLE PRECISION SUM
      INTEGER I, J, N
C
C-----------------------------------------------------------------------
C
      I=1
  200 SUM=0.0D0
      J=1
  201 SUM=SUM+A(J,I)*V(J)
      IF (J .EQ. N) GO TO 203
      J=J+1
      GO TO 201
  203 X(I)=SUM
      IF (I .EQ. N) RETURN                               ! Normal return
      I=I+1
      GO TO 200
      END!
C
C#######################################################################
C
      SUBROUTINE  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 (ID, NSTAR, INDEX, HOLD)
      IMPLICIT NONE
C
      REAL  HOLD(*)
      INTEGER INDEX(*),ID(*)
C
      INTEGER I, NSTAR
C
      DO I=1,NSTAR
         HOLD(I)=ID(I)
      END DO
      DO I=1,NSTAR
         ID(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

