      SUBROUTINE CALCELL
C CALCULATES AN ELLIPSE FROM THE MAJOR AXIS, POSITION ANGLE,
C AND ELLIPTICITY, CONVERTS IT TO RECTANGULAR COORDINATES, 
C THEN TRANSLATES IT TO THE DESIRED CENTER.  ALSO CALCULATES 
C THE POSITION ANGLE OF EACH CALCULATED POINT ALONG THE ELLIPSE.

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

      INTEGER*4  I, J

      THETA_ELL(1)   = 0.

      DO I = 1,ELLPTS-5
C CALCULATE ELLIPSE; ORIGIN IS AT ONE FOCUS
          R_ELL(I)   = A_APP * (1. - E_APP**2)
     1                           / ( 1. + E_APP * COS(THETA_ELL(I)) )
C ROTATE AND CONVERT TO RECTANGULAR COORDINATES
          RA_ELL(I)  = R_ELL(I) * SIN( THETA_ELL(I) + A_ANGLE ) 
          DEC_ELL(I) = R_ELL(I) * COS( THETA_ELL(I) + A_ANGLE ) 
C TRANSLATE SO CENTER IS AT ORIGIN
          RA_ELL(I)  = RA_ELL(I)  
     1                     + SQRT(A_APP**2 - B_APP**2) * SIN(A_ANGLE)
          DEC_ELL(I) = DEC_ELL(I) 
     1                     + SQRT(A_APP**2 - B_APP**2) * COS(A_ANGLE)
C TRANSLATE TO NEW CENTER
          RA_ELL(I)  = RA_ELL(I)  + CRA
          DEC_ELL(I) = DEC_ELL(I) + CDEC
C
          THETA_ELL(I+1) = THETA_ELL(I) + 2.*PI / (ELLPTS - 5)
      END DO
C FILL IN THE REST OF THE ARRAY
      DO I = 1, 5
          RA_ELL(ELLPTS+I-5) = RA_ELL(I)
          DEC_ELL(ELLPTS+I-5)= DEC_ELL(I)
      END DO

C CALCULATE POSITION ANGLE FOR EACH CALCULATED POINT
      DO I = 1, ELLPTS
          THETA_PT(I)= ATAN2( RA_ELL(I) , DEC_ELL(I) )
          IF(THETA_PT(I).LT.0.) THETA_PT(I) = THETA_PT(I)+2.*PI
      END DO

C CALCULATE AREAS AND TIMES
      AREA_ELL       = PI * A_APP * B_APP
      AREA_TIME      = AREA_ELL / PERIOD

C      CUM_TIME(1)    = TOBSREF
C      CUM_AREA(1)    = 0.

      DO I = 1, ELLPTS
          J          = I - SENSE
          IF (J.GT.ELLPTS) J = J - (ELLPTS-5)
          IF (J.LT.1     ) J = J + (ELLPTS-5)
          AREA(I)    = 0.5 * SENSE * ( RA_ELL(I)*DEC_ELL(J)
     1                                 - DEC_ELL(I)*RA_ELL(J) )
          TIME(I)    = AREA(I) / AREA_TIME

      END DO

C MAKE A SINGLE-PRECISION COPY OF RA_ELL, DEC_ELL FOR PLOTTING
      DO I = 1, ELLPTS
          RA_ELL4(I) = RA_ELL(I)
          DEC_ELL4(I)= DEC_ELL(I)
      END DO

C      WRITE(PR,'(1X,A)') 'CALCELL:'
C      WRITE(PR,'(1X,A)')
C     1 'PT    THETA   TIME     PT    THETA   TIME     '//
C     2 'PT    THETA   TIME     PT    THETA   TIME'
      DO I = 1, ELLPTS/4+1
C          WRITE (PR,'(1X,4(I3,1X,F8.5,1X,F6.3,4X))')
C     1     I,THETA_PT(I),TIME(I),
C     2     I+ELLPTS/4+1,THETA_PT(I+ELLPTS/4+1),TIME(I+ELLPTS/4+1),
C     3     I+ELLPTS/4+1+ELLPTS/4+1,THETA_PT(I+ELLPTS/4+1+ELLPTS/4+1),
C     4                      TIME(I+ELLPTS/4+1+ELLPTS/4+1),
C     5     I+ELLPTS/4+1+ELLPTS/4+1+ELLPTS/4+1,
C     6              THETA_PT(I+ELLPTS/4+1+ELLPTS/4+1+ELLPTS/4+1),
C     7              TIME(I+ELLPTS/4+1+ELLPTS/4+1+ELLPTS/4+1)
      END DO
C      WRITE(PR,'(1X,A)') 'CALCELL:'
C      WRITE(PR,'(1X,A)')
C     1 'PT     R.A.    TIME     PT     R.A.    TIME     '//
C     2 'PT     R.A.    TIME     PT     R.A.    TIME'
      DO I = 1, ELLPTS/4+1
C          WRITE (PR,'(1X,4(I3,1X,F9.5,1X,F6.3,4X))')
C     1     I,RA_ELL(I),TIME(I),
C     2     I+ELLPTS/4+1,RA_ELL(I+ELLPTS/4+1),TIME(I+ELLPTS/4+1),
C     3     I+ELLPTS/4+1+ELLPTS/4+1,RA_ELL(I+ELLPTS/4+1+ELLPTS/4+1),
C     4                      TIME(I+ELLPTS/4+1+ELLPTS/4+1),
C     5     I+ELLPTS/4+1+ELLPTS/4+1+ELLPTS/4+1,
C     6              RA_ELL(I+ELLPTS/4+1+ELLPTS/4+1+ELLPTS/4+1),
C     7              TIME(I+ELLPTS/4+1+ELLPTS/4+1+ELLPTS/4+1)
      END DO
C      WRITE(PR,'(1X,A)') 'CALCELL:'
C      WRITE(PR,'(1X,A)')
C     1 'PT      DEC    TIME     PT      DEC    TIME     '//
C     2 'PT      DEC    TIME     PT      DEC    TIME'
      DO I = 1, ELLPTS/4+1
C          WRITE (PR,'(1X,4(I3,1X,F9.5,1X,F6.3,4X))')
C     1     I,DEC_ELL(I),TIME(I),
C     2     I+ELLPTS/4+1,DEC_ELL(I+ELLPTS/4+1),TIME(I+ELLPTS/4+1),
C     3     I+ELLPTS/4+1+ELLPTS/4+1,DEC_ELL(I+ELLPTS/4+1+ELLPTS/4+1),
C     4                      TIME(I+ELLPTS/4+1+ELLPTS/4+1),
C     5     I+ELLPTS/4+1+ELLPTS/4+1+ELLPTS/4+1,
C     6              DEC_ELL(I+ELLPTS/4+1+ELLPTS/4+1+ELLPTS/4+1),
C     7              TIME(I+ELLPTS/4+1+ELLPTS/4+1+ELLPTS/4+1)
      END DO

      RETURN
      END

