      SUBROUTINE FITELL
C ESTIMATE THE POSITION ANGLES AT THE OBSERVED TIMES BASED ON THE 
C APPARENT ELLIPSE AND THE EQUAL-AREAS LAW.
C VARIABLES WITHIN THIS SUBROUTINE ARE:
C     AREA(N)      AREA OF LITTLE SLICES FROM ONE CALCULATED POINT 
C                  OF THE ELLIPSE TO THE NEXT (mas**2)
C     TIME(N)      TIME TO SWEEP OUT AREA(N) (days)
C     AREA_ELL     AREA OF THE WHOLE ELLIPSE (mas**2)
C     AREA_TIME    AREA SWEPT OUT PER UNIT TIME (mas**2/day)
C     SUM_AREA     TOTAL AREA SWEPT OUT FROM REFERENCE TO CURRENT POS.
C     SUM_TIME     TOTAL TIME REQUIRED TO SWEEP OUT SUM_AREA
C     THETA_PT(N)  ANGLE BTW ORIGIN (i.e., PRIMARY STAR) AND 
C                  POINT ALONG CALCULATED ELLIPSE (radians)
C     TOBSMOD(N)   OBSERVATION TIMES MODULO THE PERIOD, SHIFTED TO
C                  BRING TOBSMOD(IREF) TO ZERO (days)
C     IREF         INDEX OF THE DATA POINT FROM WHICH TO START COUNTING
C                  TIME
C     INEAR        INDEX OF THE POINT ALONG THE CALCULATED ELLIPSE
C                  NEAREST TO DATA POINT 'IREF'
C     ELLORDER(N)  LIST OF INDEXES OF OBSERVATION POINTS, IN ORDER AS
C                  THEY APPEAR AROUND THE ELLIPSE STARTING FROM IREF
C     RA_REF       RA OF ESTIMATE OF REFERENCE DATA POINT POSITION
C     DEC_REF      DEC OF ESTIMATE OF REFERENCE DATA POINT POSITION
C     TOBSREF      TIME OF REFERENCE DATA POINT OBSERVATION
C
C     STIMEINT(3)  } ARRAYS OF CURRENT THREE POINTS OF SUM_TIME, 
C     TH_PTINT(3)  } THETA_PT, RA_ELL, DEC_ELL FROM WHICH TO INTERPOLATE
C     RAELLINT(3)  } VALUES FOR ESTIMATES OF WHERE THE DATA POINT IS
C     DECELLINT(3) } PREDICTED TO LIE GIVEN THE MODEL ELLIPSE
C
C     DIST         } FOR SEARCHING FOR CALCULATED ELLIPSE POINT NEAREST
C     DISTANCE     } THE REFERENCE POSITION
C
C     CHISQ        CHI SQUARED OF (ESTIMATED - OBSERVED)

      IMPLICIT UNDEFINED (A-Z)
      INCLUDE 'PARAMETR.INC'
      INCLUDE 'ELLIPSE.INC'

      CHARACTER*5  FITSTRING
      INTEGER*4    I, J, K, L, M, LL, MM, IREF, INEAR, NEXTPT, DUMMY
      INTEGER*4    ELLORDER(MAXPTS)

      REAL*8       XA(3), YA(3), WA(3), ZA(3)
      REAL*8       STIMEINT(3), TH_PTINT(3), RAELLINT(3), DECELLINT(3)
      REAL*8       AREA_ELL, SUM_AREA, SUM_TIME, ERRSIZE, ERRS
      REAL*8       RA_REF, DEC_REF
      REAL*8       DISTANCE, DIST, CHISQ



c-------- b e g i n n i n g    o f  p r o g r a m    e l l f i t ------

c Before running fitell, let's do a grid search to improve
c the apparent ellipse you fit to the data. This is accomplished
c by some simple subroutines, which don't share the COMMON blocks
c and get their data with the subroutine call. (C.A. Hummel, 5/16/1991)
c ellfit uses the external subroutines and functions fchisq, gridls
c and fixell. The latter calculates with the current ellipse
c parameters the angles of the data points and the radii of
c corresponding ellipse points. The radial differences are used to
c calculate a weighted chisqr (mode=1), using weights which are
c calculated as the projection of the major axis (or minor, whichever 
c is larger) onto the radius each iteration.
c           
      real*8 radius(maxpts),phi(maxpts),sigmar(maxpts)
      real*8 radfit(maxpts)
      real*8 chiold,chisqre,chifr
      real*8 angle
      real*8 ellpar(5),dellpar(5),sellpar(5)
      real*8 fchisq
      integer*4 mode,nterms,nfree,iter
      character*80 text
      mode=1
      nterms=5
      chifr=0.001
      ellpar(1)=dble(-cra)
      ellpar(2)=dble(cdec)
      ellpar(3)=a_angle+pi/2.
      ellpar(4)=a_app
      ellpar(5)=b_app
      dellpar(1)=0.5
      dellpar(2)=0.5
      dellpar(3)=0.1
      dellpar(4)=0.5
      dellpar(5)=0.5
      call fixell(phi,radius,npts,ellpar,radfit,ra,dec,maxpts)
      do i=1,npts
       angle=phi(i)+ellpar(3)+pi/2.-phi_err(i)
       sigmar(i)=max(abs(dra(i)*cos(angle)),
     &                abs(ddec(i)*sin(angle)))
      enddo
      nfree=2*npts-nterms
      chiold=fchisq(radius,sigmar,npts,nfree,mode,radfit,maxpts)
      do iter=1,100
       call gridls(phi,radius,sigmar,npts,nterms,mode,
     &      ellpar,dellpar,sellpar,radfit,chisqre,ra,dec,maxpts)
       if(1.-chisqre/chiold.lt.chifr)goto 1111
       chiold=chisqre
       do i=1,npts
       angle=phi(i)+ellpar(3)+pi/2.-phi_err(i)
       sigmar(i)=max(abs(dra(i)*cos(angle)),
     &                abs(ddec(i)*sin(angle)))
       enddo
      enddo
 1111 write(text,1112)-ellpar(1),ellpar(2),ellpar(3)*180./pi-pi/2.,
     & ellpar(4),ellpar(5),chisqre
 1112 format('cra= ',f5.1,' cdec= ',f5.1,' pa= ',f5.1,' a= ',f4.1,
     & ' b= ',f4.1,' chisqre= ',f6.3)
      call locate(0,2)
      call write_string(text)
      cra=sngl(-ellpar(1))
      cdec=sngl(ellpar(2))
      a_angle=ellpar(3)-pi/2.
      a_app=ellpar(4)
      b_app=ellpar(5)
      bra(1)=cra-sngl(b_app*cos(a_angle))
      bdec(1)=cdec+sngl(b_app*sin(a_angle))
      bra(2)=2.*cra-bra(1)
      bdec(2)=2.*cdec-bdec(1)
      ara(1)=cra-sngl(a_app*sin(a_angle))
      adec(1)=cdec-sngl(a_app*cos(a_angle))
      ara(2)=2.*cra-ara(1)
      adec(2)=2.*cdec-adec(1)
      e_app=sqrt(a_app**2-b_app**2)/a_app
 
c-------------- e n d   o f   p r o g r a m    e l l f i t ----------

D     WRITE(PR,*) 'FITELL:  NPTS = ', NPTS

      DISTANCE    = 1.E9
      IREF        = 1

      IF (.NOT.ELLFIT) THEN

C---------------------------------------------------------------
C FIND THE POINT WITH THE SMALLEST UNCERTAINTY
C---------------------------------------------------------------


          DO I = 1, NPTS
              DIST    = DRA(I)**2 + DDEC(I)**2
              IF (DIST.LE.DISTANCE.AND..NOT.DELETE(I)) THEN
                DISTANCE= DIST
                IREF  = I                  ! "IREF" LABELS REFERENCE
              END IF                       ! DATA POINT
          END DO

      ELSE

C---------------------------------------------------------------
C GET THE DATA POINT TO BE USED AS THE REFERENCE FROM WHICH TO
C START COUNTING AREAS
C POSSIBLE VARIATION:  TAKE FIRST POINT AS REFERENCE, THEN 
C VARY ITS ANGLE TO GET MINIMUM CHI SQUARED.
C---------------------------------------------------------------

          CALL LOCATE (0,1)
          CALL WRITE_STRING('Place cursor at            ')
          CALL LOCATE (0,2)  
          CALL WRITE_STRING('reference data point       ')
          CALL LOCATE (0,3)  
          CALL WRITE_STRING('and hit any key           ')           
          CALL LOCATE (0,4)
          CALL WRITE_STRING('                  ')
          CALL PGCURSE(XR,YR,COMND)

          DO I = 1, NPTS
              DIST        = (RA(I) - XR)**2 + (DEC(I) - YR)**2
              IF (DIST.LT.DISTANCE.AND..NOT.DELETE(I)) THEN
                DISTANCE  = DIST
                IREF      = I                  ! "IREF" LABELS REFERENCE
              END IF                           ! DATA POINT
          END DO

          XR              = RA(IREF)
          YR              = DEC(IREF)

      END IF

C---------------------------------------------------------------
C GET THE ELLIPSE
C---------------------------------------------------------------

      CALL CALCELL

C---------------------------------------------------------------
C FIND ELLIPSE POINT NEAREST IN ANGLE TO REFERENCE OBSERVATION 
C---------------------------------------------------------------

      DISTANCE        = 1000.

      DO I = 1, ELLPTS
          DIST        = ABS(THETA_PT(I) - THETAOBS(IREF))
          IF (DIST.LT.DISTANCE) THEN
            DISTANCE  = DIST
            INEAR     = I
          END IF
      END DO                 

C---------------------------------------------------------------
C INTERPOLATE TO FIND POSITION ON ELLIPSE AT ANGLE OF IREF
C---------------------------------------------------------------

C      WRITE(PR,*) ' FITELL INIT FOR IREF SEARCH'

      DO I = 1, 3
          J           = INEAR - SENSE*(I-1)
          IF (J.LE.0)      J = J + (ELLPTS-5)
          IF (J.GT.ELLPTS) J = J - (ELLPTS-5)
          TH_PTINT(I) = THETA_PT(J)
          RAELLINT(I) = RA_ELL(J)
          DECELLINT(I)= DEC_ELL(J)
      END DO
C      WRITE(PR,'(A,3(F6.2))')
C     1   ' TH_PTINT = ',TH_PTINT(1),TH_PTINT(2),TH_PTINT(3)

C CHECK TO SEE IF TH_PTINT ARRAY STRADDLES ZERO; IF SO, FIX IT
      IF(SENSE.EQ.1.AND.TH_PTINT(3).GT.TH_PTINT(2)) THEN
          TH_PTINT(1) = TH_PTINT(1) + 2.*PI
          TH_PTINT(2) = TH_PTINT(2) + 2.*PI
      END IF
      IF(SENSE.EQ.-1.AND.TH_PTINT(3).LT.TH_PTINT(2)) THEN
          TH_PTINT(1) = TH_PTINT(1) - 2.*PI
          TH_PTINT(2) = TH_PTINT(2) - 2.*PI
      END IF
C      WRITE(PR,'(1X,A,3(F6.2,1X),1X,A,I2,A,F5.2)')
C     1  'FITELL: TH_PTINT = ',TH_PTINT(1),TH_PTINT(2),TH_PTINT(3),
C     2  '; THETAOBS(',IREF,')=',THETAOBS(IREF)

C CHECK TO SEE IF TH_PTINT HAS BEEN FIXED SO IT'S IN THE NEIGHBORHOOD
C OF THETAOBS(IREF)
      IF(TH_PTINT(3)-THETAOBS(IREF).GT.PI) THEN
          DO I = 1, 3
              TH_PTINT(I)=TH_PTINT(I) - 2.*PI
          END DO
      END IF
      IF(TH_PTINT(3)-THETAOBS(IREF).LT.-PI) THEN
          DO I = 1, 3
              TH_PTINT(I)=TH_PTINT(I) + 2.*PI
          END DO
      END IF

      CALL POLINT( TH_PTINT, RAELLINT,  THETAOBS(IREF), RA_REF  )
      CALL POLINT( TH_PTINT, DECELLINT, THETAOBS(IREF), DEC_REF )

      RA_EST(IREF)    = RA_REF
      DEC_EST(IREF)   = DEC_REF
      R_EST(IREF)     = SQRT( RA_REF**2 + DEC_REF**2 )
      THETA_EST(IREF) = THETAOBS(IREF)

C---------------------------------------------------------------
C CALCULATE OBSERVATION TIMES MODULO PERIOD REFERRED TO IREF
C---------------------------------------------------------------

      DO I = 1, NPTS
          TOBSMOD(I)  = MOD(TOBS(I),PERIOD)
      END DO

      TOBSREF         = TOBSMOD(IREF)

      DO I = 1, NPTS
          TOBSMOD(I)  = TOBSMOD(I) - TOBSREF
          IF (TOBSMOD(I).LT.0.) TOBSMOD(I) = TOBSMOD(I) + PERIOD
      END DO

      TOBSMOD(IREF)   = 0.

C      DO I = 1,NPTS
C          WRITE(PR,'(1X,A,I2,A,F8.3)')
C     1        'FITELL: TOBSMOD(',I,') = ',TOBSMOD(I)
C      END DO

C---------------------------------------------------------------
C FIGURE OUT THE ORDER IN WHICH THE OBSERVATIONS FALL AROUND
C THE ELLIPSE
C---------------------------------------------------------------

      CALL INDEXX(NPTS,TOBSMOD,ELLORDER)

C      WRITE(PR,*) 'FITELL: ORDER OF POINTS:'
C      WRITE(PR,'(1X,I2)') (ELLORDER(I),I=1,NPTS)

C---------------------------------------------------------------
C CALCULATE AREA AND TIME FROM REFERENCE POSITION TO NEXT ELLIPSE 
C POINT
C---------------------------------------------------------------

      SUM_AREA        = 0.5 * SENSE * ( RA_ELL(INEAR)*DEC_REF 
     1                                  - DEC_ELL(INEAR)*RA_REF )
      SUM_TIME        = SUM_AREA / AREA_TIME

C---------------------------------------------------------------
C INITIALIZE INTERPOLATION ARRAYS
C---------------------------------------------------------------

C      WRITE(PR,*) ' FITELL INIT FOR PREDICTED POSITION SEARCH'

C      WRITE(PR,'(A,I3,A,I2)') 'FITELL: INEAR=',INEAR,';  SENSE=',SENSE
C      WRITE(PR,'(1X,A,F7.3,A,F7.3,A,F10.5)') 'RA_REF=',RA_REF,
C     1               '   DEC_REF=',DEC_REF,
C     2               '   AREA_TIME=',AREA_TIME
      DO I = 1, 3
          J           = INEAR + SENSE * (I-3)
          IF(J.LT.1)      J = J + (ELLPTS-5)
          IF(J.GT.ELLPTS) J = J - (ELLPTS-5)
          STIMEINT(I) = 0.5 * SENSE * ( DEC_REF*RA_ELL(J)
     1                      - RA_REF*DEC_ELL(J) ) / AREA_TIME
          TH_PTINT(I) = THETA_PT(J)
          RAELLINT(I) = RA_ELL  (J)
          DECELLINT(I)= DEC_ELL (J)
C          WRITE(PR,'(1X,I3,A,F7.3,A,F7.3)')
C     1        J, '  RA=',RA_ELL(J), '  DEC=',DEC_ELL(J)
      END DO

C          WRITE(PR,'(A,3(F6.2))')
C     1       ' TH_PTINT = ',TH_PTINT(1),TH_PTINT(2),TH_PTINT(3)

C CHECK TO SEE IF TH_PTINT ARRAY STRADDLES ZERO; IF SO, FIX IT
          IF(SENSE.EQ.1.AND.TH_PTINT(3).LT.TH_PTINT(2)) THEN
              TH_PTINT(1) = TH_PTINT(1) - 2.*PI
              TH_PTINT(2) = TH_PTINT(2) - 2.*PI
C          WRITE(PR,'(A,3(F6.2),A)')
C     1       ' TH_PTINT = ',TH_PTINT(1),TH_PTINT(2),TH_PTINT(3),
C     2       ' AFTER FIXING'
          END IF
          IF(SENSE.EQ.-1.AND.TH_PTINT(3).GT.TH_PTINT(2)) THEN
              TH_PTINT(1) = TH_PTINT(1) + 2.*PI
              TH_PTINT(2) = TH_PTINT(2) + 2.*PI
C          WRITE(PR,'(A,3(F6.2),A)')
C     1       ' TH_PTINT = ',TH_PTINT(1),TH_PTINT(2),TH_PTINT(3),
C     2       ' AFTER FIXING'
          END IF

C---------------------------------------------------------------
C LOOP FOR FINDING PREDICTED POSITIONS AT OBSERVED TIMES
C---------------------------------------------------------------

      K               = 1
      NEXTPT          = ELLORDER(K)
C      WRITE(PR,1000)
 1000 FORMAT('  K NXT TOBSMOD ______STIMEINT______  ____TH_PTINT_____',
     1  '  ____RAELLINT_____  ____DECELLINT____',
     2  ' RAEST  DECEST REST  THEST')
C      WRITE(PR,1001) K,NEXTPT,TOBSMOD(NEXTPT),
C     1               STIMEINT(1),STIMEINT(2),STIMEINT(3),
C     2               TH_PTINT(1),TH_PTINT(2),TH_PTINT(3),
C     3               RAELLINT(1),RAELLINT(2),RAELLINT(3),
C     4               DECELLINT(1),DECELLINT(2),DECELLINT(3),
C     5               RA_EST(NEXTPT),DEC_EST(NEXTPT),R_EST(NEXTPT),
C     6               THETA_EST(NEXTPT)
 1001 FORMAT(1X,I2,1X,I2,2X,F6.2,2X,F6.2,1X,F6.2,1X,F6.2,2X,
     1       F5.2,1X,F5.2,1X,F5.2,2X,F5.1,1X,F5.1,1X,F5.1,2X,
     2       F5.1,1X,F5.1,1X,F5.1,2X,F5.1,2X,F5.1,2X,F4.1,2X,F5.2)

      K               = 2
      NEXTPT          = ELLORDER(K)

C      WRITE(PR,*) ' FITELL: LOOP FOR FINDING PREDICTED POSITIONS'

      DO I = 1, ELLPTS

          J           = INEAR + I * SENSE        ! SET INDEX
          IF (J.GT.ELLPTS) J = J - (ELLPTS-5)
          IF (J.LT.   1  ) J = J + (ELLPTS-5)

          SUM_AREA    = SUM_AREA + AREA(J)       ! INCREMENT AREA
          SUM_TIME    = SUM_TIME + TIME(J)       ! AND TIME

C          WRITE(PR,'(1X,A,I4,1X,I4)') 
C     1       ' FITELL: PERIASTR LOOP POINT 2. I,J = ',I,J
C          WRITE(PR,'(1X,A,E11.4,2X,E11.4)')
C     1       ' FITELL: RA_ELL, DEC_ELL = ',RA_ELL(J), DEC_ELL(J)

          CALL PGSCI(9)
          CALL PGPOINT(1,RA_ELL4(J),DEC_ELL4(J),3)
          CALL EGA_RESTORE_DEFAULT

          STIMEINT(1) = STIMEINT(2)              ! UPDATE INTERPOLATION
          STIMEINT(2) = STIMEINT(3)              ! ARRAYS
          STIMEINT(3) = SUM_TIME

          TH_PTINT(1) = TH_PTINT(2)
          TH_PTINT(2) = TH_PTINT(3)
          TH_PTINT(3) = THETA_PT(J)

C          WRITE(PR,'(A,I3,A,3(F6.2))')
C     1       'I= ',I,' TH_PTINT = ',TH_PTINT(1),TH_PTINT(2),TH_PTINT(3)

C CHECK TO SEE IF TH_PTINT ARRAY STRADDLES ZERO; IF SO, FIX IT
          IF(SENSE.EQ.1.AND.TH_PTINT(3).LT.TH_PTINT(2)) THEN
              TH_PTINT(1) = TH_PTINT(1) - 2.*PI
              TH_PTINT(2) = TH_PTINT(2) - 2.*PI
C          WRITE(PR,'(A,I3,A,3(F6.2),A)')
C     1       'I= ',I,' TH_PTINT = ',TH_PTINT(1),TH_PTINT(2),TH_PTINT(3),
C     2       ' AFTER FIXING'
          END IF
          IF(SENSE.EQ.-1.AND.TH_PTINT(3).GT.TH_PTINT(2)) THEN
              TH_PTINT(1) = TH_PTINT(1) + 2.*PI
              TH_PTINT(2) = TH_PTINT(2) + 2.*PI
C          WRITE(PR,'(A,I3,A,3(F6.2),A)')
C     1       'I= ',I,' TH_PTINT = ',TH_PTINT(1),TH_PTINT(2),TH_PTINT(3),
C     2       ' AFTER FIXING'
          END IF

          RAELLINT(1) = RAELLINT(2)
          RAELLINT(2) = RAELLINT(3)
          RAELLINT(3) = RA_ELL(J)

          DECELLINT(1)= DECELLINT(2)
          DECELLINT(2)= DECELLINT(3)
          DECELLINT(3)= DEC_ELL(J)

   50     IF (SUM_TIME.GT.TOBSMOD(NEXTPT)) THEN  ! INTERPOLATE CALC-
                                                 ! ULATED POSITION

              CALL POLINT(STIMEINT, TH_PTINT,  TOBSMOD(NEXTPT),   
     1                    THETA_EST(NEXTPT))
              CALL POLINT(TH_PTINT, RAELLINT,  THETA_EST(NEXTPT), 
     1                    RA_EST(NEXTPT)   )
              CALL POLINT(TH_PTINT, DECELLINT, THETA_EST(NEXTPT), 
     1                    DEC_EST(NEXTPT)  )
              IF(THETA_EST(NEXTPT).LT.0.)
     1              THETA_EST(NEXTPT)=THETA_EST(NEXTPT)+2.*PI
              IF(THETA_EST(NEXTPT).GT.2.*PI)
     1              THETA_EST(NEXTPT)=THETA_EST(NEXTPT)-2.*PI

              R_EST(NEXTPT) = SQRT( RA_EST(NEXTPT)**2 
     1                              + DEC_EST(NEXTPT)**2 )
C      WRITE(PR,1001) K,NEXTPT,TOBSMOD(NEXTPT),
C     1               STIMEINT(1),STIMEINT(2),STIMEINT(3),
C     2               TH_PTINT(1),TH_PTINT(2),TH_PTINT(3),
C     3               RAELLINT(1),RAELLINT(2),RAELLINT(3),
C     4               DECELLINT(1),DECELLINT(2),DECELLINT(3),
C     5               RA_EST(NEXTPT),DEC_EST(NEXTPT),R_EST(NEXTPT),
C     6               THETA_EST(NEXTPT)

              K       = K + 1
              IF (K.GT.NPTS) THEN
                  CALL LOCATE(0,24)
                  WRITE(PR,'(A,I3,A)') 'FITELL: ',K-1,' POINTS FOUND'
                  GO TO 100
              ELSE
                  NEXTPT  = ELLORDER(K)
              END IF

              GO TO 50                   ! CHECK TO SEE IF NEXT POINT
                                         ! HAS BEEN PASSED TOO
          END IF

      END DO

  100 CONTINUE         ! JUMP TO HERE AFTER ALL POINTS FOUND

      IF (K.LE.NPTS) THEN
          CALL LOCATE(0,24)
          CALL WRITE_STRING ('NOT ALL POINTS WERE FOUND')
          WRITE(PR,*) ' FITELL: NOT ALL POINTS FOUND'
      END IF

      CALL CHISQR(R_EST,ROBS,THETA_EST,THETAOBS,DRA,DDEC,PHI_ERR,
     1            NPTS,DELETE,CHISQ,DELRA,DELDEC)

      WRITE(TITLE(6),1002) CHISQ, NPTS-DEL_COUNT
 1002 FORMAT(' CHISQ = ',F10.3,' FOR ',I3,' POINTS')

      ELLFIT          = .TRUE.

      RETURN
      END


