      SUBROUTINE CALCTI 
C CALCULATE THIELE-INNES CONSTANTS: 
      IMPLICIT UNDEFINED (A-Z)
      INCLUDE 'PARAMETR.INC'
      INCLUDE 'ELLIPSE.INC'

      INTEGER*4  I, J, K, L, INEAR, INTERP_CT, COUNT

      REAL*4     RA_TI1_4, DEC_TI1_4, RA_TI2_4, DEC_TI2_4

      REAL*8     TH_PTINT(3), RAELLINT(3),DECELLINT(3),RMOD,THMOD
      REAL*8     RAMOD, DECMOD
      REAL*8     THETAJ, THETAK, THETAL
      REAL*8     X_TI1, X_TI2, Y_TI1, Y_TI2
      REAL*8     RA_TI1, DEC_TI1, RA_TI2, DEC_TI2, THETA_90
      REAL*8     DISTI, DISTJ, TANSQI, SUM, DIFF
      REAL*8     DUM1, DUM2, DUM3, DUM4
      REAL*8     SUM_TIME, TIMEL, STIMEINT(3), T_90, TORB
      REAL*8     I_PRINT, AN_PRINT, AP_PRINT, A(7), DYDA(7,2)

      CALL EGA_RESTORE_DEFAULT

C--------------------------------------------------------------------
C DETERMINE LOCATION AND TIME OF PERIASTRON 
C--------------------------------------------------------------------

      THETA_LINE      = ATAN2(CRA,CDEC) + PI
                                ! ANGLE OF LINE FROM CENTER TO ORIGIN;
                                ! ADD PI SO WE GET POINT CLOSER TO
                                ! ORIGIN

      IF (THETA_LINE.GE.2*PI) THETA_LINE = THETA_LINE - 2.*PI
      IF (THETA_LINE.LT.0.  ) THETA_LINE = THETA_LINE + 2.*PI
                                ! PUT THETA_LINE BETWEEN ZERO AND 2*PI

C--------------------------------------------------------------------
C FIRST, INITIALIZE INTERPOLATION ARRAYS
C--------------------------------------------------------------------

      TH_PTINT(1)     = THETA_PT(ELLPTS-5)
      TH_PTINT(2)     = THETA_PT(1)
      TH_PTINT(3)     = THETA_PT(2)

      RAELLINT(1)     = RA_ELL  (ELLPTS-5)
      RAELLINT(2)     = RA_ELL  (1)
      RAELLINT(3)     = RA_ELL  (2)

      DECELLINT(1)    = DEC_ELL (ELLPTS-5)
      DECELLINT(2)    = DEC_ELL (1)
      DECELLINT(3)    = DEC_ELL (2)

C      WRITE(PR,'(1X,A,F8.5,1X,F8.5,1X,F8.5)')
C     1   'CALCTI INIT TO FIND INEAR: TH_PTINT(1,2,3)= ',
C     2   TH_PTINT(1),TH_PTINT(2),TH_PTINT(3)

C--------------------------------------------------------------------
C SECOND, FIND POINT ALONG ELLIPSE WITH POSITION ANGLE EQUAL TO
C THETA_LINE; FIND ELLIPSE POINT NEAREST IN ANGLE TO PERIASTRON
C POSITION
C--------------------------------------------------------------------

      DO I = 1, ELLPTS-4
          J               = I + 1
          IF (J.GT.ELLPTS) J = J - (ELLPTS-5)
          K               = I + 2
          IF (K.GT.ELLPTS) K = K - (ELLPTS-5)
          IF (THETA_PT(I).LT.THETA_LINE  .AND.
     1                     THETA_PT(J).GE.THETA_LINE) THEN

              CALL POLINT( TH_PTINT, RAELLINT,  THETA_LINE, RA_TI1  ) 
              CALL POLINT( TH_PTINT, DECELLINT, THETA_LINE, DEC_TI1 )

              RA_TI1_4    = RA_TI1
              DEC_TI1_4   = DEC_TI1
 
              DISTI       = ABS(THETA_PT(I) - THETA_LINE)
              DISTJ       = ABS(THETA_PT(J) - THETA_LINE)

              IF (DISTI.LT.DISTJ) THEN
                  INEAR   = I
              ELSE
                  INEAR   = J
              END IF

              WRITE(PR,'(1X,A,F8.3,1X,F8.3)')
     1           ' CALCTI: RA_TI1, DEC_TI1 = ',RA_TI1, DEC_TI1
C              WRITE(PR,'(1X,A,I3,A,I2)') ' CALCTI: INEAR= ',INEAR,
C     1                 '; SENSE= ',SENSE

              GO TO 100

          ELSE
              TH_PTINT(1) = TH_PTINT(2)
              TH_PTINT(2) = TH_PTINT(3)
              TH_PTINT(3) = THETA_PT(K)

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

              DECELLINT(1)= DECELLINT(2)
              DECELLINT(2)= DECELLINT(3)
              DECELLINT(3)= DEC_ELL(K)
          END IF
      END DO

  100 CONTINUE              ! COME TO HERE WHEN RA_TI1, DEC_TI1 FOUND

C--------------------------------------------------------------------
C AND PLOT THE PERIASTRON ON THE PLOT
C--------------------------------------------------------------------

      CALL PGSCI(2)                              ! 2 = RED
      CALL PGPOINT(1, RA_TI1, DEC_TI1, 8)        ! 8 = CIRCLE W/ CROSS
      CALL PGPTEXT(RA_TI1, DEC_TI1, 0., 0., ' PERIASTR')
      CALL PGSCI(1)

C--------------------------------------------------------------------
C THEN CALCULATE ECCENTRICITY OF TRUE ORBIT FROM PERIASTRON POSITION
C--------------------------------------------------------------------

      E_TRUE              = SQRT((CRA**2 + CDEC**2)/
     1                           ((CRA-RA_TI1)**2 + (CDEC-DEC_TI1)**2))

C--------------------------------------------------------------------
C COUNT AHEAD TO POSITION OF TRUE ANOMALY = 90 DEG
C FIRST, CALCULATE TIME TO TRUE ANOMALY = 90 DEG
C--------------------------------------------------------------------

      T_90                = ( PERIOD/(2*PI) ) 
     1                        * ( ACOS(E_TRUE) 
     2                            - E_TRUE * SQRT((1-E_TRUE**2)) )

C--------------------------------------------------------------------
C SECOND, CALCULATE TIME FROM PERIASTRON TO NEXT ELLIPSE POINT
C--------------------------------------------------------------------

C      WRITE(OUTC,*) ' RA_ELL, INEAR, DEC_TI1,DEC_ELL, RA_TI1, AREA_TIME'
C      WRITE(OUTC,*) RA_ELL(INEAR)
C      WRITE(OUTC,*) INEAR
C      WRITE(OUTC,*) DEC_TI1
C      WRITE(OUTC,*) DEC_ELL(INEAR)
C      WRITE(OUTC,*) RA_TI1
C      WRITE(OUTC,*) AREA_TIME

      SUM_TIME            = 0.5 * SENSE * ( RA_ELL(INEAR)*DEC_TI1
     1                           - DEC_ELL(INEAR)*RA_TI1 )/ AREA_TIME

C--------------------------------------------------------------------
C THIRD, INITIALIZE INTERPOLATION ARRAYS
C--------------------------------------------------------------------

      DO I = 1, 3
          J               = INEAR+SENSE*(I-3)
          IF(J.LT.1)      J = J + (ELLPTS-1)
          IF(J.GT.ELLPTS) J = J - (ELLPTS-1)
          STIMEINT(I)     = 0.5 * SENSE * ( DEC_TI1*RA_ELL(J)
     1                         - RA_TI1*DEC_ELL(J) ) /AREA_TIME
          TH_PTINT(I)     = THETA_PT(J)
          RAELLINT(I)     = RA_ELL  (J)
          DECELLINT(I)    = DEC_ELL (J)
      END DO

C      WRITE(PR,'(1X,A,F8.5,1X,F8.5,1X,F8.5)')
C     1   'CALCTI INIT FOR TRUE ANOM. =90: TH_PTINT(1,2,3)= ',
C     2   TH_PTINT(1),TH_PTINT(2),TH_PTINT(3)
C      WRITE(PR,'(1X,A,F8.5,1X,F8.5,1X,F8.5)')
C     1   '                                STIMEINT(1,2,3)= ',
C     2   STIMEINT(1),STIMEINT(2),STIMEINT(3)

C--------------------------------------------------------------------
C FOURTH, ADD UP TIMES UNTIL TRUE ANOMALY = 90 DEG IS REACHED
C--------------------------------------------------------------------

      DO I = 1, ELLPTS
          J               = INEAR + I * SENSE
          IF (J.GT.ELLPTS) J = J - (ELLPTS-5)
          IF (J.LT.1     ) J = J + (ELLPTS-5)

          SUM_TIME        = SUM_TIME + TIME(J)

          STIMEINT(1)     = STIMEINT(2)
          STIMEINT(2)     = STIMEINT(3)
          STIMEINT(3)     = SUM_TIME   

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

          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)

                                            
          IF (SUM_TIME.GT.T_90) THEN
              CALL POLINT( STIMEINT, TH_PTINT,  T_90,     THETA_90 )
              CALL POLINT( TH_PTINT, RAELLINT,  THETA_90, RA_TI2   )
              CALL POLINT( TH_PTINT, DECELLINT, THETA_90, DEC_TI2  )
              RA_TI2_4    = RA_TI2
              DEC_TI2_4   = DEC_TI2
              GO TO 200
          END IF
      END DO

  200 CONTINUE                      ! JUMP TO HERE AFTER T_90 REACHED

      CALL PGSCI(7)                              ! 7 = YELLOW
      CALL PGPOINT(1, RA_TI2, DEC_TI2, 8)        ! 8 = CIRCLE W/ CROSS
      CALL PGPTEXT(RA_TI2, DEC_TI2, 0., 0., ' TA90')
      CALL PGSCI(1)

C--------------------------------------------------------------------
C DETERMINE TIME AT PERIASTRON BY COUNTING BACKWARDS FROM THE FIRST 
C OBSERVATION DATE.  
C FIRST, FIND THE ELLIPSE POINT NEAREST THE ESTIMATED POSITION
C OF THE FIRST OBSERVATION AND INITIALIZE INTERPOLATION ARRAYS.
C--------------------------------------------------------------------

      DISTI           = 1000.
      DO I = 1, ELLPTS-5
          DISTJ       = ABS(THETA_PT(I) - THETA_EST(1))
          IF (DISTJ.LT.DISTI) THEN
              DISTI   = DISTJ
              INEAR   = I
          END IF
      END DO

      DO I = 1, 3
          J           = INEAR - SENSE * (I-1)
          IF (J.GT.ELLPTS) J = J - (ELLPTS-5)
          IF (J.LT.1     ) J = J + (ELLPTS-5)
          TIMEL       = 0.5 * SENSE * ( RA_ELL(J)*DEC_EST(1)
     1                                    - DEC_ELL(J)*RA_EST(1) )
     2                    / AREA_TIME
          STIMEINT(I) = JULDAY(1) - TIMEL
          TH_PTINT(I) = THETA_PT(J)
C          WRITE(PR,'(1X,A,I3,1X,I3,1X,F9.5)') 
C     1         'CALCTI INIT FOR PERIASTR: I,J, TIMEL = ',
C     2                                 I,J,TIMEL
      END DO

C      WRITE(PR,'(1X,A,F8.5,1X,F8.5,1X,F8.5)')
C     1   'CALCTI INIT TO FIND PERIASTRON: TH_PTINT(1,2,3)= ',
C     2   TH_PTINT(1),TH_PTINT(2),TH_PTINT(3)
C      WRITE(PR,'(1X,A,F9.5,1X,F9.5,1X,F9.5)')
C     1   '                     STIMEINT(1,2,3) - 2447000 = ',
C     2   STIMEINT(1) - 2447000. ,STIMEINT(2) - 2447000. ,
C     3   STIMEINT(3) - 2447000. 

C--------------------------------------------------------------------
C SECOND, SEARCH BACKWARD ALONG THE ORBIT UNTIL PERIASTRON IS REACHED,
C THEN INTERPOLATE TO FIND THE EPOCH
C--------------------------------------------------------------------

C      WRITE(PR,'(1X,A)') 'CALCTI:'
C      WRITE(PR,'(1X,A)') '.I. .J. .K. .L.  '
C     1     //'.........TH_PTINT.........'

      DO I = 1, ELLPTS
                                              
          J           = INEAR - SENSE * (I-1) 
          K           = INEAR - SENSE * I     
          L           = INEAR - SENSE * (I+1) 
          IF (J.GT.ELLPTS) J = J - (ELLPTS-5)
          IF (J.LT.1     ) J = J + (ELLPTS-5)
          IF (K.GT.ELLPTS) K = K - (ELLPTS-5)
          IF (K.LT.1     ) K = K + (ELLPTS-5)
          IF (L.GT.ELLPTS) L = L - (ELLPTS-5)
          IF (L.LT.1     ) L = L + (ELLPTS-5)

C          THETAJ      = THETA_PT(J)
C          THETAK      = THETA_PT(K)
C          THETAL      = THETA_PT(L)

C DIRECT ORBITS
          IF(SENSE.EQ.1) THEN
            IF (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

C RETROGRADE ORBITS
          ELSE
            IF (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
          END IF

C          WRITE(PR,'(1X,I3,1X,I3,1X,I3,1X,I3,2X,
C     1           F8.5,1X,F8.5,1X,F8.5)')
C     2        I,J,K,L,
C     3        TH_PTINT(1),TH_PTINT(2),TH_PTINT(3)

C DIRECT ORBITS
          IF (       SENSE.EQ.1
     1                  .AND.
     2        TH_PTINT(3).LT.THETA_LINE
     3                  .AND.
     4        TH_PTINT(2).GE.THETA_LINE) THEN
            CALL POLINT(TH_PTINT,STIMEINT,THETA_LINE,EPOCH)
            GO TO 300
          END IF

C RETROGRADE ORBITS
          IF (       SENSE.EQ.-1
     1                  .AND.
     2        TH_PTINT(3).GT.THETA_LINE
     3                  .AND.
     4        TH_PTINT(2).LE.THETA_LINE) THEN
            CALL POLINT(TH_PTINT,STIMEINT,THETA_LINE,EPOCH)
            GO TO 300
          END IF

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

          STIMEINT(1) = STIMEINT(2)
          STIMEINT(2) = STIMEINT(3)
          STIMEINT(3) = STIMEINT(3) -TIME(L)

      END DO

  300 CONTINUE                    ! COME HERE AFTER EPOCH FOUND

C--------------------------------------------------------------------
C CALCULATE THIELE-INNES CONSTANTS
C--------------------------------------------------------------------

      X_TI1               = RA_TI1  / (1. - E_TRUE)
      Y_TI1               = DEC_TI1 / (1. - E_TRUE)
      X_TI2               = RA_TI2  / (1. - E_TRUE**2)
      Y_TI2               = DEC_TI2 / (1. - E_TRUE**2)

C      WRITE(PR,1000) X_TI1, Y_TI1, X_TI2, Y_TI2
 1000 FORMAT(1X,'X1 = ',F6.2,'  Y1 = ',F6.2,'   X2 = ',F6.2,
     1        '  Y2 = ',F6.2)

C--------------------------------------------------------------------
C CALCULATE ORBITAL ELEMENTS
C--------------------------------------------------------------------

      SUM                 = DATAN2( (X_TI1 - Y_TI2), (X_TI2 + Y_TI1) )
      DIFF                = DATAN2(-(X_TI1 + Y_TI2), -(X_TI2 - Y_TI1) )

C      WRITE(PR,*) 'SUM, DIFF = ',SUM,DIFF

      ASC_NODE            = 0.5 * (SUM - DIFF)
      IF (ASC_NODE.LT.0.) ASC_NODE = ASC_NODE + PI
      IF (ASC_NODE.GE.PI) ASC_NODE = ASC_NODE - PI

      ARG_PERI            = SUM - ASC_NODE
      IF (ARG_PERI.LT.0.) ARG_PERI = ARG_PERI+2.*PI

      TANSQI              = -(X_TI1 + Y_TI2) * SIN(SUM)
     1                         / ( (X_TI1 - Y_TI2) * SIN(DIFF) )
      INCLIN              = 2. * DATAN( SQRT(TANSQI) )

      A_TRUE              = (X_TI2 + Y_TI1) 
     1                         / ( COS(SUM) * (1.+COS(INCLIN)) )

      AN_PRINT            = ASC_NODE / RPDEG    ! CONVERT ANGLES TO
      AP_PRINT            = ARG_PERI / RPDEG    ! DEGREES
      I_PRINT             = INCLIN   / RPDEG

      WRITE(PR,1006) 
      WRITE(PR,1007) A_TRUE, E_TRUE, PERIOD, EPOCH, I_PRINT, AN_PRINT,
     1               AP_PRINT
 1006 FORMAT(' ORBITAL ELEMENTS FROM THIELE-INNES METHOD:')
 1007 FORMAT(' SEMIMAJOR AXIS:   ',4X,F6.2,' mas',  /,
     1       ' ELLIPTICITY:      ',6X,F5.3,          /,
     2       ' PERIOD:           ',3X,F10.5,' days', /,
     3       ' EPOCH:         JD ',F11.3,         /,
     4       ' INCLINATION:      ',3X,F6.1, ' deg',  /,
     5       ' ASCENDING NODE:   ',3X,F6.1, ' deg',  /,
     6       ' ARG. PERIASTRON:  ',3X,F6.1, ' deg'     )

C--------------------------------------------------------------------
C PLOT THE RESULTING TRUE ORBIT ON THE SCREEN.  IT SHOULD OVERLAP
C THE APPARENT ORBIT
C--------------------------------------------------------------------

      TORB             = EPOCH

      A(1)             = A_TRUE
      A(2)             = E_TRUE
      A(3)             = PERIOD
      A(4)             = EPOCH
      A(5)             = INCLIN
      A(6)             = ASC_NODE
      A(7)             = ARG_PERI

C      WRITE(PR,*) ' CALCTI:'
C      WRITE(PR,'(A)')
C     1  ' I    TORB-EPOCH     YMOD(1,2)   RA_MOD  DEC_MOD'

      IF(ORBPTS.GT.ELLPTS) ORBPTS = ELLPTS      ! Arrays are dimensioned
                                                ! with  ELLPTS
      DO I = 1, ORBPTS
          TORB         = TORB + PERIOD/(ORBPTS)
          CALL TRUE2APP(TORB,A,RMOD,THMOD,RAMOD,DECMOD,DYDA,PR)
          RA_MOD(I)    = RAMOD
          DEC_MOD(I)   = DECMOD
C          WRITE(PR,'(1X,I3,2X,F11.3,2X,F6.2,2X,F6.2,2X,F6.2,2X,F6.2)')
C     1      I,TORB,RMOD,THMOD,RA_MOD(I),DEC_MOD(I)
      END DO

      TRUFIT           = .TRUE.

C--------------------------------------------------------------------
C WRITE MODEL FILE WITH ELLIPSE ELEMENTS AND THIELE-INNES CONSTANTS
C     PARAMETERS TO WRITE ARE: MAJOR AXIS, ECCENTRICITY, P.A. OF 
C         MAJOR AXIS, COORDINATES OF CENTER, FOUR THIELE-INNES 
C         CONSTANTS, TIME OF PERIASTRON, PERIOD.
C     OTHER STUFF TO WRITE:  STAR NAME AND NUMBER, INPUT DATA FILE
C         NAME, DATE OF FIT, UNCERTAINTIES
C--------------------------------------------------------------------

      RETURN
      END

