ACEPH.FTN
FTN77 
$FILES 0,1
      PROGRAM ACEPH (6,200), CALLABLE COMET EPHEMERIS <850822.1153> 
C     ************************************************************* 
C   *****  TWO BODY COMET EPHEMERIS GENERATOR  *****
C 
C   THE TWO BODY COMET EPHEMERIS GENERATOR (ACEPH) IS A CALLABLE
C   PROGRAM DESIGNED TO PRODUCE PREDICTED OBERVABLES OF A 
C   COMET FROM  OSCULATING ELEMENTS THAT ARE INPUT BY THE USER. 
C   IT ACCOMPLISHES THIS THROUGH THE COMPUTATION OF A TWO-BODY ORBIT
C   OF THE COMET FROM THE OSCULATING ELEMENTS AND THE COMPUTATION OF
C   SOLAR COORDINATES FROM NEWCOMB'S THEORY. THUS THE PROGRAM AVOIDS
C   THE NECESSITY OF NUMERICALLY INTEGRATING THE COMET ORBIT, AS WELL 
C   AS READING AN EXTERNAL PLANETARY EPHEMERIS. 
C 
C   THIS PROGRAM WAS WRITTEN AT THE JET PROPULSION LABORATORY BY
C   DR. M.S.W. KEESEY. SUBROUTINE HELXYZ WAS WRITTEN BY DR. D.K. YEOMANS. 
C   SUBROUTINE SUN WAS DEVELOPED AT THE U.S. NAVAL OBSERVATORY. 
C 
C   THIS PROGRAM DOES NOT USE MORE THAN ONE SET OF OSCULATING ELEMENTS FOR A
C   GIVEN EPHEMERIS INTERVAL; THAT IS, EACH EPHEMERIS INTERVAL MUST 
C   CORRESPOND TO ONE AND ONLY ONE SET OF OSCULATING ELEMENTS.
C 
C   INITIALIZING PARAMETERS READ IN ARE DEFINED AS FOLLOWS... 
c   NOTE THAT THIS PROGRAM READS FROM THE SAME SOURCE FILES AS 'TBCEG'. 
C   This version of the program is hard wired to run option 2 only. 
C 
C   input file line 1 : 
C      'NAME'   COMET NAME 80 characters maximum
C 
C   input file line 2 : 
C      'NAPPS'  IS THE NUMBER OF SETS OF ELEMENTS TO BE READ IN.
C               ie the number of repeats of line 4, max = 10 sets 
C 
C   input file line 3 : (this is read for option 1 which does not use it) 
C      'STACOR' IS THE OBSERVATORY COORDINATES ARRAY of 3 elements :
C               east longitude in degrees,
C               rho*cos (geocentric latitude),
C               rho*sin (geocentric latitude).
C               where rho is the earth radius at the station divided by 
C               the equatorial radius of 6378160 meters.
C      'ETMUTC' IS THE DIFFERENCE BETWEEN EPHEMERIS TIME AND UNIVERSAL
C                      TIME COORDINATED in seconds. 
C 
C   input file line 4 : 
C      'ELMTS'  IS THE ARRAY OF 6 OSCULATING ELEMENTS : 
C               Julian date of perihelion passage (ET), 
C               perihelion distance in a.u.,
C               eccentricity, 
C               argument of perihelion (degrees), 
C               longitude of ascending node (degrees),
C               inclination (degrees).
C 
C   input file line 5 : 
C      'NEL'    IS THE KEY FOR THE LOCATION OF THE CORRECT SET OF OSCULATING
C                      ELEMENTS FOR THE CURRENT EPHEMERIS INTERVAL. 
C      'TSTRT'  IS THE BEGINNING TIME OF THE CURRENT EPHEMERIS INTERVAL,
C               in julian days. 
C      'TSTOP'  IS THE ENDING TIME OF THE CURRENT EPHEMERIS INTERVAL, 
C               in julian days. 
C      'H'      IS THE STEP SIZE OF THE EPHEMERIS IN JULIAN DAYS. 
C      AND SEVEN PARAMETERS FOR COMPUTING THE MAGNITUDE ESTIMATES.
C      These are not necessary and the line is terminated with a / .
C 
C      SAMPLE INPUT PARAMETER FILE FOR HARTRAO: 
C   HALLEY'S COMET
C   1 
C   27.684722 0.900407079 -0.434097561 55.2 
C   2446470.95542 0.5871548 0.9672617 111.84343 58.14254 162.24011
C   1 2446290.5 2446470.5 1 / 
C 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      implicit integer*4 (i-n)
      integer*2 ipram(5), inbuf(20), ioutbuf(10), 
     *          lu, iprint
      character name*80, infile*20, csign*4, csign2*4 
      real*4    outra50, outdec50, outra, outdec, outvel
      DIMENSION ELMTS(10,6),II(10),JJ(10),KK(10),FHR(10),STACOR(3)
      DIMENSION RSG(6),RCH(3),RC(6),PUT(4),ROG(6) 
      DIMENSION RSGA(6,5),RCHA(3,5),RCA(6,5),RHOV(5),RCOMV(5) 
      DIMENSION TV(5),RUV(5),PUTA(4,5)
* 
      equivalence (inbuf,infile), (csign,isign), (csign2,isign2), 
     *            (ioutbuf(1), outra50), (ioutbuf(3), outdec50),
     *            (ioutbuf(5), outra),   (ioutbuf(7), outdec),
     *            (ioutbuf(9), outvel)
* 
      COMMON /ELM/  TT,Q,ECR,PX,PY,PZ,QX,QY,QZ
      COMMON /TOPO/ IOPT,AU,REM,STACOR,ETMUTC 
* 
      DATA TWOPI,RAD,CSE,SNE/6.283185307179586D0,0.0174532925199D0, 
     *     0.9174369524418883D0,0.3978811861525251D0/ 
      DATA DPR/57.2957795130823252D0/ 
C 
C     pick up parameters from the calling program 
      call rmpar (ipram)
C     pick up comet ephemeris file name from inbuf
      call exec (14,1,inbuf,20) 
C     i/o lu, print info option from input parameters 
      lu = ipram(1) 
      iprint = ipram(2) 
      iopt = 2
C 
C     open the file with the comet ephemeris parameters 
      open (5, file=infile) 
*     read the Comet name (80 characters) 
      read (5,'(a)') name 
      READ (5,*) NAPPS
C     READ AND WRITE THE OBSERVATORY COORDINATES AND THE DIFFERENCE 
C     BETWEEN ET AND UTC
      READ  (5,*) (STACOR(I),I=1,3),ETMUTC
      if (iprint .eq. 1) write (lu,601) 
  601 FORMAT(' OBSERVATORY COORDINATES AND E.T.-U.T.C. '/ 
     *' EAST LONGITUDE ',3X,' RHO*COS(GEOCENTRIC LAT.) ',3X,
     *' RHO*SIN(GEOCENTRIC LAT.) ',3X,' ET MINUS UTC ') 
      if (iprint .eq. 1) write (lu,602) (STACOR(I),I=1,3),ETMUTC
  602 FORMAT(F13.8,7X,F14.10,15X,F14.10,16X,F9.5) 
      STACOR(1)=STACOR(1)/DPR 
      ETMUTC=ETMUTC/86400.D0
C 
C     READ AND WRITE THE TABLE OF ORBITAL ELEMENTS. 
  700 if (iprint .eq. 1) write (lu,701) 
  701 FORMAT(' MATRIX OF ORBITAL ELEMENTS',/,2X,' YR  MN DY  HR', 
     *7X,'J.D.',11X,'Q',11X,'E',7X,'SOMEGA',5X,'LOMEGA',8X,'I') 
      DO I=1,NAPPS
          READ  (5,*) (ELMTS(I,J),J=1,6)
C         COMPUTE THE CALENDAR DATE OF THE INSTANT OF PERIHELION PASSAGE. 
          CALL DAYCNV(ELMTS(I,1),II(I),JJ(I),KK(I),FHR(I))
          if (iprint .eq. 1)
     *      write (lu,704) II(I),JJ(I),KK(I),FHR(I),(ELMTS(I,J),J=1,6)
      end do
  704     FORMAT(I5,2I3,F5.1,2X,F13.5,2F11.7,3F11.5)
C 
C     BEGIN MAIN EPHEMERIS COMPUTATION. 
c     Enter Key to elements in array
  709 read  (5,*) NEL 
      close (5) 
C 
C     WRITE THE PAGE HEAD FOR THE EPHEMERIS OUTPUT. 
C     IF(IOPT.EQ.1) write(6,711) name(1:80) 
C 711 FORMAT(//27X,'TWO BODY EPHEMERIS FOR ',A,//,1X,' YR',2X,'MN', 
C    *1X,'DY',2X,'HR',3X,'J.D.(ET)',4X,'R.A.',1X,'1950.0',1X,'DEC.',6X, 
C    *'R.A.',2X,'DATE',2X,'DEC.',4X,'DELTA',4X,'R',5X,'TMAG',2X,'NMAG', 
C    *2X,'THETA',3X,'BETA',2X,'HLONG',1X,'HLAT')
* 
      IF (IOPT .EQ. 2) write(lu,712) name(1:80) 
  712 FORMAT(/'TWO BODY EPHEMERIS FOR ',A,//,1X,' YR',2X,'MN',
     *1X,'DY',2X,'HR',3X,'J.D.(UT)',4X,'R.A.',1X,'1950.0',1X,'DEC.',6X, 
     *'R.A.',2X,'DATE',2X,'DEC.',4X,'DELTA',2X,'DELTAD',4X,'R',4X,
     *'RDOT',2X,'THETA',3X,'BETA')
* 
C     IF(IOPT.EQ.3) write(6,713) name(1:80) 
C 713 FORMAT(//27X,'TWO BODY EPHEMERIS FOR ',A,//,1X,' YR',2X,'MN', 
C    *1X,'DY',2X,'HR',3X,'J.D.(UT)',4X,'R.A.',1X,'1950.0',1X,'DEC.',6X, 
C    *'R.A.',2X,'DATA',2X,'DEC.',4X,'DELTA',5X,'R',4X,'PSANG',3X, 
C    *'AZ',4X,'ALT',2X,'TLONG',2X,'TLAT',2X,'THETA',3X,'BETA')
C 
C   REARRANGE STOR. OF ORB. ELEMENTS, AND CONV. ANGULAR QUANTITIES TO RADIANS.
      TT =ELMTS(NEL,1)
      Q  =ELMTS(NEL,2)
      ECR=ELMTS(NEL,3)
      WR =ELMTS(NEL,4)*RAD
      OMR=ELMTS(NEL,5)*RAD
      XIR=ELMTS(NEL,6)*RAD
C 
C   COMPUTE THE P AND Q ROTATION FACTORS. 
      PP=DSIN(WR)*DCOS(XIR)*DCOS(OMR)+DCOS(WR)*DSIN(OMR)
      QQ=DCOS(WR)*DCOS(XIR)*DCOS(OMR)-DSIN(WR)*DSIN(OMR)
      PX=-(DSIN(WR)*DCOS(XIR)*DSIN(OMR)-DCOS(WR)*DCOS(OMR)) 
      PY=PP*CSE-SNE*DSIN(XIR)*DSIN(WR)
      PZ=CSE*DSIN(XIR)*DSIN(WR)+PP*SNE
      QX=-(DCOS(WR)*DCOS(XIR)*DSIN(OMR)+DSIN(WR)*DCOS(OMR)) 
      QY=-(SNE*DSIN(XIR)*DCOS(WR)-QQ*CSE) 
      QZ=QQ*SNE+CSE*DSIN(XIR)*DCOS(WR)
C 
C   NTRIG=0 SECTION - 
C   COMPUTE THE HELIOCENTRIC COORDINATES OF THE COMET, THE GEOCENTRIC SOLAR 
C   COORDINATES, AND THE HELIOCENTRIC  COMET DISTANCE.
C   IF IOPT=1, ALSE COMPUTE THE GEOCENTRIC COMET DISTANCE, AND APPLY
C              THE LIGHT TIME CORRECTION. 
C   IF IOPT=2, COMPUTE THE TOPOCENTRIC COMET DISTANCE, AND APPLY THE
C                LIGHT TIME CORRECTION AND THE ELLIPTIC PORTION OF THE
C               ANNUAL PART OF THE STELLAR ABERRATION.
C     find current UT in Julian days
      call jd(tstrt)
C     observation time = start = end time, interval = 0.1 days
      tobs  = tstrt 
      tstop = tstrt 
      h = 0.1 
  714 NTRIG=0 
      IF(IOPT.GE.2) GO TO 716 
C 715 CALL APPAD1(NTRIG,TOBS,PUT,RSG,RCH,RC)
C     GO TO 718 
  716 CALL APPAD2(NTRIG,TOBS,RU,ECLOBL,PUT,RSG,RCH,RC,ROG,POSANG) 
C 
C   COMPUTE THE EPHEMERIS ALPHA DELTA, CH TO HRS, MIN 
718   RHO=DSQRT(RC(1)**2+RC(2)**2+RC(3)**2) 
      DEC= DASIN(RC(3)/RHO) 
      RA=DATAN2(RC(2),RC(1))
      IF(RA.LT.0.D0) RA=RA+TWOPI
      CALL ARGCHG(RA,DEC,IRAHR,RAMIN,ISIGN,IDCDG,DECMIN)
C 
C   COMPUTE EARTH-COMET-SUN ANGLE BETA, AND THE TOTAL AND NUCLEAR MAGNITUDES. 
      GRHO=DSQRT((RSG(1)+RCH(1))**2+(RSG(2)+RCH(2))**2+ 
     *(RSG(3)+RCH(3))**2) 
      RCOM=DSQRT(RCH(1)*RCH(1)+RCH(2)*RCH(2)+RCH(3)*RCH(3)) 
      RSUNSQ=RSG(1)*RSG(1)+RSG(2)*RSG(2)+RSG(3)*RSG(3)
      BETA= DACOS((RCOM*RCOM+GRHO*GRHO-RSUNSQ)/(2.D0*RCOM*GRHO))/RAD
      XMAG =XMT+XTNE*DLOG10(RHO)+XTNS*DLOG10(RCOM)
      XMAG2=XMN+XNNE*DLOG10(RHO)+XNNS*DLOG10(RCOM)+P9*BETA
C 
C   COMPUTE THE SUN-EARTH-COMET ANGLE.
      RSUN=DSQRT(RSUNSQ)
      THETA= DACOS((GRHO*GRHO+RSUNSQ-RCOM*RCOM)/(2.D0*GRHO*RSUN))/RAD 
      IF(IOPT.NE.1) GO TO 719 
C 
C   COMPUTE THE HELIOCENTRIC ECLIPTIC LONGITUDE AND LATITUDE. 
C     Y1= RCH(2)*CSE+RCH(3)*SNE 
C     Z1=-RCH(2)*SNE+RCH(3)*CSE 
C     RP=DSQRT(RCH(1)*RCH(1)+Y1*Y1) 
C     RE=DSQRT(RCH(1)*RCH(1)+Y1*Y1+Z1*Z1) 
C     BBETA=DASIN(Z1/RE)/RAD
C     XLAMDA=DATAN2(Y1/RP,RCH(1)/RP)/RAD
C     IF(XLAMDA.LT.0.D0)XLAMDA=XLAMDA+360.D0
  719 IF(IOPT.LE.2) GO TO 720 
C 
C   COMPUTE THE TOPOCENTRIC ECLIPTIC LONGITUDE AND LATITUDE.
C     YY1=RC(2)*CSE+RC(3)*SNE 
C     ZZ1=-RC(2)*SNE+RC(3)*CSE
C     RRP=DSQRT(RC(1)*RC(1)+YY1*YY1)
C     RRE=DSQRT(RC(1)*RC(1)+YY1*YY1+ZZ1*ZZ1)
C     TOPLAT=DASIN(ZZ1/RRE)/RAD 
C     TOPLNG=DATAN2(YY1/RRP,RC(1)/RRP)/RAD
C     IF(TOPLNG.LT.0.D0) TOPLNG=TOPLNG+360.D0 
C 
C   NTRIG SET TO 1 IN ORDER TO APPLY PRECESSION, NUTATION, AND THE
C   CIRCULAR PORTION OF THE ANNUAL PART AND THE DIURNAL PART OF THE 
C   STELLAR ABERRATION TO THE TOPOCENTRIC COMET COORDINATES.
  720 NTRIG=1 
      IF(IOPT.GE.2) GO TO 722 
C     CALL APPAD1(NTRIG,TOBS,PUT,RSG,RCH,RC)
C     GO TO 723 
  722 CALL APPAD2(NTRIG,TOBS,RU,ECLOBL,PUT,RSG,RCH,RC,ROG,POSANG) 
C 
C   COMPUTE APPARENT RIGHT ASCENSION AND DECLINATION. 
  723 RHO=DSQRT(RC(1)*RC(1)+RC(2)*RC(2)+RC(3)*RC(3))
      DEC2=DASIN(RC(3)/RHO) 
      RA2=DATAN2(RC(2),RC(1)) 
      IF(RA2.LT.0.D0) RA2=RA2+TWOPI 
      CALL ARGCHG(RA2,DEC2,IRAHR2,RAMIN2,ISIGN2,IDCD2,DECMN2) 
C 
C   CONVERT THE JULIAN DATE INTO A CALENDAR DATE. 
      CALL DAYCNV(TOBS,IYR,JMON,KDAY,FQR) 
      IF(IOPT.EQ.2) GO TO 725 
c     IF(IOPT.EQ.3) GO TO 760 
C 
C   PRINT OPTION 1 RESULTS. 
c     write(6,724)IYR,JMON,KDAY,FQR,TOBS,IRAHR,RAMIN,ISIGN,IDCDG,DECMIN,
c    XIRAHR2,RAMIN2,ISIGN2,IDCD2,DECMN2,RHO,RCOM,XMAG,XMAG2,THETA,
c    XBETA,XLAMDA,BBETA 
c 724 FORMAT(I5,2I3,F5.1,F11.1,I3,F7.3,2X,A1,I2,F6.2,I4,F7.3,2X,A1,I2,
c    XF6.2,F8.4,F7.4,2F6.2,1X,F6.2,1X,F6.2,F7.1,F6.1) 
c     GO TO 780 
C 
C   COMPUTE THE COMET-OBSERVER AND COMET-SUN RELATIVE VELOCITIES. 
  725 TM=TOBS-2.D0*H
      TLAST=TOBS+2.D0*H 
      IF(TM.GE.TSTRT) GO TO 730 
      DO 729 K=1,5
          CALL APPAD2(0,TM,RU,ECLOBL,PUT,RSG,RCH,RC,ROG,POSANG) 
          CALL APPAD2(1,TM,RU,ECLOBL,PUT,RSG,RCH,RC,ROG,POSANG) 
          TV(K)=TM
          RUV(K)=RU 
          DO 726 I=1,4
  726         PUTA(I,K)=PUT(I)
          RHOV(K)=DSQRT(RC(1)*RC(1)+RC(2)*RC(2)+RC(3)*RC(3))
          RCOMV(K)=DSQRT(RCH(1)*RCH(1)+RCH(2)*RCH(2)+RCH(3)*RCH(3)) 
          DO 727 I=1,6
              RSGA(I,K)=RSG(I)
  727         RCA(I,K)=RC(I)
          DO 728 I=1,3
  728         RCHA(I,K)=RCH(I)
          TM=TM+H 
  729     CONTINUE
      GO TO 739 
  730 DO 735 K=1,4
          TV(K)=TV(K+1) 
          RHOV(K)=RHOV(K+1) 
          RCOMV(K)=RCOMV(K+1) 
          DO 733 I=1,3
  733         RCHA(I,K)=RCHA(I,K+1) 
          DO 734 I=1,4
  734         PUTA(I,K)=PUTA(I,K+1) 
          DO 735 I=1,6
              RSGA(I,K)=RSGA(I,K+1) 
  735         RCA(I,K)=RCA(I,K+1) 
      CALL APPAD2(0,TLAST,RU,ECLOBL,PUT,RSG,RCH,RC,ROG,POSANG)
      CALL APPAD2(1,TLAST,RU,ECLOBL,PUT,RSG,RCH,RC,ROG,POSANG)
      TV(5)=TLAST 
      RUV(5)=RU 
      DO 736 I=1,4
  736     PUTA(I,5)=PUT(I)
      DO 737 I=1,6
          RSGA(I,5)=RSG(I)
  737     RCA(I,5)=RC(I)
      DO 738 I=1,3
  738     RCHA(I,5)=RCH(I)
      RHOV(5)=DSQRT(RC(1)*RC(1)+RC(2)*RC(2)+RC(3)*RC(3))
      RCOMV(5)=DSQRT(RCH(1)*RCH(1)+RCH(2)*RCH(2)+RCH(3)*RCH(3)) 
  739 CALL VELCTY(TOBS,H,RHOV,RCD)
      CALL VELCTY(TOBS,H,RCOMV,RCOMD) 
      RCD=RCD*AU/86400.D0 
      RCOMD=RCOMD*AU/86400.D0 
C 
C   PRINT OPTION 2 RESULTS. 
      write (lu,750) IYR,JMON,KDAY,FQR,TOBS,IRAHR,RAMIN,ISIGN,IDCDG,
     *decmin,IRAHR2,RAMIN2,ISIGN2,IDCD2,DECMN2,RHO,RCD,RCOM,RCOMD,
     *THETA,BETA
  750 FORMAT(I5,2I3,F5.1,F11.1,I3,F7.3,2X,A1,I2,F6.2,I4,F7.3,2X,A1,I2,
     *F6.2,F8.4,1X,F6.2,F7.4,F7.2,1X,F6.2,1X,F6.2)
      GO TO 780 
C 
C   COMPUTE TOPOCENTRIC ALTITUDE AND AZIMUTH OF THE COMET.
c 760 CALL AZALT(RU,RA2,DEC2,AZ,ALT)
C 
C   PRINT OPTION 3 RESULTS. 
c     POSANG=POSANG*DPR 
c     ALT=ALT*DPR 
c     AZ=AZ*DPR 
c     write(6,775)IYR,JMON,KDAY,FQR,TOBS,IRAHR,RAMIN,ISIGN,IDCDG,DECMIN,
c    *IRAHR2,RAMIN2,ISIGN2,IDCD2,DECMN2,RHO,RCOM,POSANG,AZ,ALT, 
c    *TOPLNG,TOPLAT,THETA,BETA
c 775 FORMAT(I5,2I3,F5.1,F11.1,I3,F7.3,2X,A1,I2,F6.2,I4,F7.3,2X,A1,I2,
c    *F6.2,F8.4,1X,F7.4,1X,F5.1,1X,2(1X,F5.1,F6.1),2(1X,F6.1))
C 
  780 TOBS=TOBS+H 
      IF(TOBS.LE.TSTOP) GO TO 714 
c 
C     store output data in ioutbuf for return to caller 
      outra50  = (ramin/60d0 + irahr) * 15d0
      outdec50 = decmin / 60d0 + idcdg
      if (csign(1:1) .eq. '-') outdec50 = -outdec50 
      outra  = (ramin2/60d0 + irahr2) * 15d0
      outdec = decmn2/60d0 + idcd2
      if (csign2(1:1) .eq. '-' ) outdec = -outdec 
      outvel = rcd
c 
c     write(6,800)II(NEL),JJ(NEL),KK(NEL),FHR(NEL), 
c    *            (ELMTS(NEL,J),J=1,6)
c 800 FORMAT(///,22X,' THE ABOVE EPHEMERIS IS BASED ON THE FOLLOWING ECL
c    *IPTIC (1950.0) ORBITAL ELEMENTS.',/,23X,' YR  MN DY  HR',7X,'J.D.'
c    *,11X,'Q',10X,'E',8X,'SOMEGA',5X,'LOMEGA',8X,'I',/ 
c    *21X,I5,2I3,F5.1,2X,F13.5,2F11.7,3F11.5) 
c     IF(IOPT.EQ.1) write(6,820) XMT,XTNE,XTNS,XMN,XNNE,XNNS,P9 
c 820 FORMAT(//'EXPLANATION OF SYMBOLS'//' J.D.= JULIAN DATE(EPHEMERIS T
c    *IME)',/' R.A. AND DEC. - 1950.0 = GEOCENTRIC RIGHT ASCENSION AND
c    *DECLINATION REFERRED TO THE MEAN EQUATOR AND EQUINOX OF 1950.0',/ 
c    *' R.A. AND DEC. - DATE = GEOCENTRIC RIGHT ASCENSION AND DECLINATIO
c    *N REFERRED TO THE APPARENT EQUATOR AND EQUINOX OF DATE,',/
c    *24X,'R.A. IN HOURS AND TIME MINUTES, DEC. IN DEGREES AND ARC MINUT
c    *ES',/' DELTA= GEOCENTRIC DISTANCE OF OBJECT IN A.U.',/
c    *' R= HELIOCENTRIC DISTANCE OF OBJECT IN A.U.',/ 
c    *' TMAG= TOTAL   MAGNITUDE = ',F4.1,' + ',F5.2,'*DLOG10(DELTA) + ',
c    *F5.2,'*DLOG10(R)',/ 
c    *' NMAG= NUCLEAR MAGNITUDE = ',F4.1,' + ',F5.2,'*DLOG10(DELTA) + ',
c    *F5.2,'*DLOG10(R) + ',F4.2,'*BETA',/ 
c    *' THETA= SUN-EARTH-OBJECT ANGLE IN DEGREES',/ 
c    *' BETA= SUN-OBJECT-EARTH ANGLE IN DEGREES',/
c    *' HLONG, HLAT = HELIOCENTRIC ECLIPTIC LONGITUDE AND LATITUDE (1950
c    *.0), IN DEGREES') 
* 
      IF (iprint .eq. 1 .and. IOPT .EQ. 2) write(LU,825)
  825 FORMAT(/'EXPLANATION OF SYMBOLS'/ 
     *' J.D.= JULIAN DATE (UNIVERSAL TIME)'/
     *' R.A. AND DEC. : 1950.0 = TOPOCENTRIC RIGHT ASCENSION AND '
     *    'DECLINATION REFERRED '/
     *26X,'TO THE MEAN EQUATOR AND EQUINOX OF 1950.0',/ 
     *' R.A. AND DEC. : DATE   = TOPOCENTRIC RIGHT ASCENSION AND '
     *    'DECLINATION REFERRED'/ 
     *26X,'TO THE APPARENT EQUATOR AND EQUINOX OF DATE,'/ 
     *26X,'R.A. IN HOURS AND TIME MINUTES, DEC. IN DEGREES '
     *    'AND ARC MINUTES'/
     *' DELTA= TOPOCENTRIC DISTANCE OF OBJECT IN A.U.',/
     *' DELTAD = D(DELTA)/DT, IN KM/SEC',/
     *' R= HELIOCENTRIC DISTANCE OF OBJECT IN A.U.',/ 
     *' RDOT = D(R)/DT, IN KM/SEC',/
     *' THETA= SUN-EARTH-OBJECT ANGLE IN DEGREES',/ 
     *' BETA= SUN-OBJECT-EARTH ANGLE IN DEGREES') 
* 
c     IF(IOPT.EQ.3) write(6,830)
c 830 FORMAT(/'EXPLANATION OF SYMBOLS',/' J.D.= JULIAN DATE(UNIVERSAL T 
c    *IME)',/' R.A. AND DEC. - 1950.0 = TOPOCENTRIC RIGHT ASCENSION AND 
c    *DECLINATION REFERRED TO THE MEAN EQUATOR AND EQUINOX OF 1950.0',/ 
c    *' R.A. AND DEC. - DATE = TOPOCENTRIC RIGHT ASCENSION AND DECLINATI
c    *ON REFERRED TO THE APPARENT EQUATOR AND EQUINOX OF DATE,',/ 
c    *24X,'R.A. IN HOURS AND TIME MINUTES, DEC. IN DEGREES AND ARC MINUT
c    *ES',/' DELTA= TOPOCENTRIC DISTANCE OF OBJECT IN A.U.',/ 
c    *' R= HELIOCENTRIC DISTANCE OF OBJECT IN A.U.',/ 
c    *' PSANG = POSITION ANGLE (EAST OF NORTH) OF EXTENDED HELIOCENTRIC 
c    *RADIUS VECTOR, IN DEGREES',/
c    *' AZ, ALT = AZIMUTH AND ALTITUDE OF OBJECT, IN DEGREES',/ 
c    *' TLONG, TLAT = TOPOCENTRIC ECLIPTIC LONGITUDE AND LATITUDE OF OBJ
c    *ECT (1950.0), IN DEGREES',/ 
c    *' THETA= SUN-EARTH-OBJECT ANGLE IN DEGREES',/ 
c    *' BETA= SUN-OBJECT-EARTH ANGLE IN DEGREES') 
* 
      call exec (14,2,ioutbuf,10) 
  200 STOP
      END 
*********************************************************************** 
      block data
      implicit double precision (A-H,O-Z) 
      implicit integer*4 (i-n)
      double precision stacor(3)
      common /ELM/ tt,q,ecr,px,py,pz,qx,qy,qz 
      common /TOPO/ iopt, au, rem, stacor, etmutc 
      data   au/1.4959787066d8/, rem/6378.140d0/
      end 
*********************************************************************** 
      SUBROUTINE ANABR3(DJ,X) 
C 
C   THIS SUBROUTINE APPLIES THE CIRCULAR PORTION OF THE ANNUAL PART OF
C   THE STELLAR ABERRATION TO INPUT COORDINATES.
C 
C   REF. - MEEUS,J., ASTRONOMICAL FORMULAE FOR CALCULATORS, MONOGRAFIEEN
C   OVER ASTRONOMIE EN ASTROFYSICA, UITGEGEVEN DOOR VOLKSSTERRENWACHT 
C   URANIA V.Z.W. EN VERENIGING VOOR STERRENKUNDE V.Z.W. VOL. 4, 1979.
C 
C   NOTE - THE INPUT RECTANGULAR COORDINATES ARE REFERRED TO THE MEAN 
C   EQUATOR AND EQUINOX OF DATE.
C 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION X(6)
      DATA DPR/57.29577951D0/ 
      DATA CONSAB/-9.933832328D-5/
      TJC=(DJ-2415020.0D0)/36525.D0 
      CMASUN=(358.47583D0+(35999.04975D0-(1.5D-4+3.3D-6*TJC)*TJC)*TJC)
     */DPR
      CC=(1.919460D0-(4.789D-3+1.4D-5*TJC)*TJC)*DSIN(CMASUN)
     *+(2.0094D-2-1.00D-4*TJC)*DSIN(2.D0*CMASUN)
     *+2.93D-4*DSIN(3.D0*CMASUN)
      CMLSUN=279.69668D0+(36000.76892D0+3.025D-4*TJC)*TJC 
      TLSUN=(DMOD(CMLSUN+CC,360.D0))/DPR
      ECLOBL=(23.452294D0-(1.30125D-2+(1.64D-6-5.03D-7*TJC)*TJC)*TJC) 
     */DPR
      SR=DSQRT(X(1)*X(1)+X(2)*X(2)+X(3)*X(3)) 
      RA=DATAN2(X(2),X(1))
      DEC=DASIN(X(3)/SR)
      SINRA=DSIN(RA)
      COSRA=DCOS(RA)
      SINDEC=DSIN(DEC)
      COSDEC=DCOS(DEC)
      SINTL=DSIN(TLSUN) 
      COSTL=DCOS(TLSUN) 
      SINEPS=DSIN(ECLOBL) 
      COSEPS=DCOS(ECLOBL) 
      TANEPS=SINEPS/COSEPS
      DELRA=CONSAB*(COSRA*COSTL*COSEPS+SINRA*SINTL)/COSDEC
      DELDEC=CONSAB*(COSTL*COSEPS*(TANEPS*COSDEC-SINRA*SINDEC)
     *+COSRA*SINDEC*SINTL)
      RA=RA+DELRA 
      DEC=DEC+DELDEC
      COSDEC=DCOS(DEC)
      X(1)=SR*COSDEC*DCOS(RA) 
      X(2)=SR*COSDEC*DSIN(RA) 
      X(3)=SR*DSIN(DEC) 
      RETURN
      END 
C******************************************************************** 
c     SUBROUTINE APPAD1(NTRIG,ETOBS,PUT,RSG,RCH,RCG)
C 
C   THIS SUBROUTINE COMPUTES THE GEOCENTRIC COMET COORDINATES.
C 
C   IF NTRIG=0, THE LIGHT TIME CORRECTION IS APPLIED TO THE COMET 
C   COORDINATES. IF NTRIG=1, PRECESSION, NUTATION, AND THE CIRCULAR 
C   PORTION OF THE ANNUAL ABERRATION ARE ALSO APPLIED.
C 
C   INPUT - 
C      NTRIG - 0 OR 1.
C      ETOBS - EPHEMERIS TIME OF OBSERVATION, IN J.D. 
C      PUT   - AN ARRAY CONTAINING DELTA PSI AND DELTA EPSILON, USED
C               IN THE NUTATION SUBROUTINE, IN RADIANS. 
C   OUTPUT- 
C      RSG - GEOCENTRIC SOLAR COORDINATES, WRT THE MEAN EQUATOR AND 
C             EQUINOX OF 1950.0, IN AU'S. 
C      RCH - HELIOCENTRIC COMET COORDINATES, WRT THE MEAN EQUATOR 
C             AND EQUINOX OF 1950.0, IN AU'S. 
C      RCG - GEOCENTRIC COMET COORDINATES, IN AU'S. 
C            IF NTRIG=0, WRT MEAN EQUATOR AND EQUINOX OF 1950.0.
C            IF NTRIG=1, WRT APPARENT EQUATOR AND EQUINOX OF DATE.
C 
c     IMPLICIT DOUBLE PRECISION(A-H,O-Z)
c     implicit integer*4 (i-n)
c     DIMENSION EL(21),C(16),RSG(6),RCH(3),RCG(6) 
c     DIMENSION PUT(4),OBM(2) 
c     COMMON/ELM/TT,Q,ECR,PX,PY,PZ,QX,QY,QZ 
c     DATA CAUPD/173.145D0/ 
c     DATA B1950/2433282.423D0/ 
c     DATA TWOPI/6.283185307179586D0/ 
C 
c     IF(NTRIG.EQ.1) GO TO 50 
c     DO 5 I=1,4
c   5 PUT(I)=0.D0 
C 
C   COMPUTE HELIOCENTRIC COMET COORDINATES. 
c     CALL HELXYZ(ETOBS,TT,Q,ECR,PX,PY,PZ,QX,QY,QZ,RCH(1),RCH(2),RCH(3))
C 
C   COMPUTE GEOCENTRIC SOLAR COORDINATES, AND PRECESS THEM TO THE 
C   MEAN EQUATOR AND EQUINOX OF 1950.0. 
c     CALL SUN(ETOBS,EL,C,RSG(1),RSG(2),RSG(3)) 
c     CALL PRECES(B1950,ETOBS,.FALSE.,RSG)
C 
C   START LIGHT TIME ITERATION. 
c  11 DO 12 I=1,3 
c  12 RCG(I)=RCH(I)+RSG(I)
c     SR=DSQRT(RCG(1)**2+RCG(2)**2+RCG(3)**2) 
c     DT=SR/CAUPD 
c     IF(ABS(DT-DTZ).LT.1.D-7) GO TO 15 
c     ETINC=ETINC+DTZ-DT
c     ETT2=ETOBS+ETINC
c     DTZ=DT
c     CALL HELXYZ(ETT2,TT,Q,ECR,PX,PY,PZ,QX,QY,QZ,RCH(1),RCH(2),RCH(3)) 
c     GO TO 11
C 
c  15 CONTINUE
C 
C   APPLY ETERMS. 
c     DEC=DASIN(RCG(3)/SR)
c     RA=DATAN2(RCG(2),RCG(1))
c     IF(RA.LT.0.D0) RA=RA+TWOPI
c     CALL ETERMS(RA,DEC) 
c     FACT=SR*DCOS(DEC) 
c     RCG(1)=FACT*DCOS(RA)
c     RCG(2)=FACT*DSIN(RA)
c     RCG(3)=SR*DSIN(DEC) 
c     GO TO 80
C 
C   PRECESS THE COMET GEOCENTRIC COORDINATES TO THE MEAN EQUATOR AND
C   EQUINOX OF DATE, AND APPLY THE CIRCULAR PORTION OF THE ANNUAL PART
C   OF THE STELLAR ABERRATION.
c  50 CALL PRECES(B1950,ETOBS,.TRUE.,RCG) 
c     CALL ANABR3(ETOBS,RCG)
C 
C   NUTATE THE COORDINATES TO THE APPARENT EQUATOR AND EQUINOX OF DATE. 
c     CALL NUTANG(ETOBS,DPSI,DEPS)
c     PUT(1)=DPSI 
c     PUT(2)=DEPS 
c     CALL NUTATE(ETOBS,PUT,OBM,.TRUE.,RCG) 
c  80 CONTINUE
c     RETURN
c     END 
C******************************************************************** 
      SUBROUTINE APPAD2(NTRIG,UTOBS,RU,ECLOBL,PUT,RSG,RCH,RCT,ROG,PSANG)
C 
C   THIS SUBROUTINE COMPUTES THE TOPOCENTRIC COMET COORDINATES. 
C 
C   IF NTRIG=0, THE LIGHT TIME CORRECTION AND THE ELLIPTIC PORTION
C   OF THE ANNUAL ABERRATION ARE APPLIED TO THE COMET COORDINATES.
C   IF NTRIG=1, PRECESSION, NUTATION, THE CIRCULAR PORTION OF THE 
C   ANNUAL ABERRATION, AND THE DIURNAL ABERRATION ARE ALSO APPLIED. 
C   THE POSITION ANGLE OF THE EXTENDED RADIUS VECTOR IS ALSO COMPUTED.
C 
C   INPUT - 
C      NTRIG  - 0 OR 1. 
C      UTOBS  - UNIVERSAL TIME OF OBSERVATION, IN J.D.
C   OUTPUT- 
C       RU     - LOCAL SIDEREAL TIME, IN RADIANS. 
C       ECLOBL - OBLIQUITY OF THE ECLIPTIC, IN RADIANS. 
C       PUT    - ARRAY CONTAINING DELTA PSI AND DELTA EPSILON, USED IN
C                 THE NUTATION SUBROUTINES, IN RADIANS. 
C       RSG    - GEOCENTRIC SOLAR COORDINATES, WRT THE MEAN EQUATOR 
C                 AND EQUINOX OF 1950.0, IN AU'S. 
C       RCH    - HELIOCENTRIC COMET COORDINATES, WRT THE MEAN EQUATOR 
C                 AND EQUINOX OF 1950.0, IN AU'S. 
C       RCT    - TOPOCENTRIC COMET COORDINATES, IN AU'S.
C                 IF NTRIG=0, WRT MEAN EQUATOR AND EQUINOX OF 1950.0. 
C                 IF NTRIG=1, WRT APPARENT EQUATOR AND EQUINOX OF DATE. 
C       ROG    - GEOCENTRIC OBSERVATORY COORDINATES, WRT THE MEAN 
C                 EQUATOR AND EQUINOX OF 1950.0, IN AU'S. 
C 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION EL(21),C(16),RSG(6),RCH(3),RCT(6),ROG(6)
      DIMENSION RCTTMP(6),RCHTMP(3),RCTAIL(6) 
      DIMENSION PUT(4),OBM(2) 
      DIMENSION STACOR(3) 
      LOGICAL FIRST 
      COMMON/ELM/TT,Q,ECR,PX,PY,PZ,QX,QY,QZ 
      COMMON/TOPO/IOPT,AU,REM,STACOR,ETMUTC 
      DATA TWOPI/6.283185307179586D0/ 
      DATA CAUPD/173.145D0/ 
      DATA DJ1900/2415020.0D0/
      DATA B1950/2433282.423D0/ 
      DATA DPR/57.2957795130823252D0/ 
      DATA RPTS/7.27220521664303933D-5/ 
C 
      FIRST=.TRUE.
      ETOBS=ETMUTC+UTOBS
      AUPER=REM/AU
      IF(NTRIG.EQ.1) GO TO 50 
      ETINC=0.D0
      DTZ=0.D0
      DO 5 I=1,4
    5 PUT(I)=0.D0 
C 
C   COMPUTE HELIOCENTRIC COMET COORDINATES. 
      CALL HELXYZ(ETOBS,TT,Q,ECR,PX,PY,PZ,QX,QY,QZ,RCH(1),RCH(2),RCH(3))
C 
C   COMPUTE GEOCENTRIC SOLAR COORDINATES AND PRECES THEM TO MEAN
C   EQUATOR AND EQUINOX OF 1950.0.
      CALL SUN(ETOBS,EL,C,RSG(1),RSG(2),RSG(3)) 
      CALL PRECES(B1950,ETOBS,.FALSE.,RSG)
C 
C   CONVERT OBSERVATORY COORDINATES FROM CYLINDRICAL TO RECTANGULAR,
C   AND ROTATE THEM TO MEAN EQUATOR AND EQUINOX OF 1950.0.
      ROG(1)=AUPER*STACOR(2)
      ROG(2)=0.D0 
      ROG(3)=AUPER*STACOR(3)
      TU=(UTOBS-DJ1900)/36525.D0
      RU=23925.836D0+(8640184.54D0+0.0929D0*TU)*TU
      J=UTOBS 
      DJ5=J 
      DJ5=DJ5-0.5D0 
      TSEC=(UTOBS-DJ5)*TWOPI
      IF(TSEC.GE.TWOPI) TSEC=TSEC-TWOPI 
      CALL NUTANG(ETOBS,DPSI,DEPS)
      PUT(1)=DPSI 
      PUT(2)=DEPS 
      TE=(ETOBS-DJ1900)/36525.D0
      ECLOBL=(23.45229444D0-(1.30125D-2+(1.64D-6-5.03D-7*TE)*TE)*TE)/DPR
      RU=DMOD((RU*RPTS+DPSI*DCOS(ECLOBL)+TSEC+STACOR(1)),TWOPI) 
      CALL ROTAT6(ROG,RU,-3,ROG)
      CALL NUTATE(ETOBS,PUT,OBM,.FALSE.,ROG)
      CALL PRECES(B1950,ETOBS,.FALSE.,ROG)
C 
C   COMPUTE TOPOCENTRIC COMET COORDINATES AND PERFORM THE LIGHT TIME
C   ITERATION.
   11 DO 12 I=1,3 
   12 RCT(I)=RCH(I)+RSG(I)-ROG(I) 
      SR=DSQRT(RCT(1)*RCT(1)+RCT(2)*RCT(2)+RCT(3)*RCT(3)) 
      DT=SR/CAUPD 
      IF(ABS(DT-DTZ).LT.1.D-7) GO TO 15 
      ETINC=ETINC+DTZ-DT
      ETT2=ETOBS+ETINC
      DTZ=DT
      CALL HELXYZ(ETT2,TT,Q,ECR,PX,PY,PZ,QX,QY,QZ,RCH(1),RCH(2),RCH(3)) 
      GO TO 11
C 
C   HERE AFTER CONVERGENCE
   15 CONTINUE
C 
C   APPLY ETERM CORRECTIONS.
      DEC=DASIN(RCT(3)/SR)
      RA=DATAN2(RCT(2),RCT(1))
      IF(RA.LT.0.D0) RA=RA+TWOPI
      CALL ETERMS(RA,DEC) 
      FACT=SR*DCOS(DEC) 
      RCT(1)=FACT*DCOS(RA)
      RCT(2)=FACT*DSIN(RA)
      RCT(3)=SR*DSIN(DEC) 
C 
      IF(IOPT.EQ.2) GO TO 80
      IF(NTRIG.EQ.0) GO TO 80 
C 
C   PRECESS TOPOCENTRIC COMET COORDINATES TO MEAN EQUATOR AND EQUINOX 
C   OF DATE, AND APPLY CIRCULAR PORTION OF ANNUAL ABERRATION. 
   50 CALL PRECES(B1950,ETOBS,.TRUE.,RCT) 
      CALL ANABR3(ETOBS,RCT)
C 
C   NUTATE COORDINATES TO APPARENT EQUATOR AND EQUINOX OF DATE, AND APPLY 
C   DIURNAL PORTION OF STELLAR ABERRATION.
      CALL NUTATE(ETOBS,PUT,OBM,.TRUE.,RCT) 
      SR=DSQRT(RCT(1)*RCT(1)+RCT(2)*RCT(2)+RCT(3)*RCT(3)) 
      DEC=DASIN(RCT(3)/SR)
      RA=DATAN2(RCT(2),RCT(1))
      IF(RA.LT.0.D0) RA=RA+TWOPI
      CALL DIABER(STACOR,RU,RA,DEC) 
      FACT=SR*DCOS(DEC) 
      RCT(1)=FACT*DCOS(RA)
      RCT(2)=FACT*DSIN(RA)
      RCT(3)=SR*DSIN(DEC) 
C 
   80 IF(IOPT.EQ.2) GO TO 99
      IF(NTRIG.EQ.0) GO TO 99 
      IF(.NOT.FIRST) GO TO 90 
C 
C   COMPUTE POSITION ANGLE OF COMET TAIL. 
      FIRST=.FALSE. 
      DO 82 I=1,6 
   82 RCTTMP(I)=RCT(I)
      DO 84 I=1,3 
      RCHTMP(I)=RCH(I)
      RCH(I)=1.1D0*RCH(I) 
   84 RCT(I)=RCH(I)+RSG(I)-ROG(I) 
      SR=DSQRT(RCT(1)*RCT(1)+RCT(2)*RCT(2)+RCT(3)*RCT(3)) 
      GO TO 15
   90 CONTINUE
      DO 92 I=1,3 
   92 RCH(I)=RCHTMP(I)
      DO 94 I=1,6 
      RCTAIL(I)=RCT(I)
   94 RCT(I)=RCTTMP(I)
      RHOT=DSQRT(RCTAIL(1)**2+RCTAIL(2)**2+RCTAIL(3)**2)
      DECT=DASIN(RCTAIL(3)/RHOT)
      RAT=DATAN2(RCTAIL(2),RCTAIL(1)) 
      IF(RAT.LT.0.D0) RAT=RAT+TWOPI 
      SR=DSQRT(RCT(1)*RCT(1)+RCT(2)*RCT(2)+RCT(3)*RCT(3)) 
      DEC=DASIN(RCT(3)/SR)
      RA=DATAN2(RCT(2),RCT(1))
      IF(RA.LT.0.D0) RA=RA+TWOPI
      DRAT=RAT-RA 
      DDECT=DECT-DEC
      SINP=DSIN(DRAT) 
      COSP=DTAN(DECT)*DCOS(DEC)-DSIN(DEC)*DCOS(DRAT)
      PSANG=DATAN2(SINP,COSP) 
      IF(PSANG.LT.0.D0) PSANG=PSANG+TWOPI 
C 
   99 CONTINUE
      RETURN
      END 
C******************************************************************** 
      SUBROUTINE ARGCHG(RA,DEC,IH,AM, ISGN,ID,AMIN) 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      implicit integer*4 (i-n)
      DATA IPLUS,IMINUS /'+','-'/ 
      DATA RAD/0.0174532925199D0/ 
C 
C**** KEY .NE. -1 ---- HOURS AND DEGREES TO RADIANS 
C**** KEY .EQ. -1 ---- RADIANS TO HOURS AND DEGREES 
C 
      KEY=-1
      IF(KEY.EQ.-1) GO TO 10
      H=IH
      D=ID
      RA=15.D0*(H+AM/6.D1)* RAD 
      DEC=(D+AMIN/6.D1)*RAD 
      IF(ISGN.EQ.1) DEC=-DEC
      RETURN
   10 HOURS=RA/RAD/15.D0
       HOURS=DMOD(HOURS,24.D0)
      IH=HOURS
      AM=(HOURS-IH)*6.D1
      AM=AM+1.D-6 
      IF(AM.LT.6.D1) GO TO 12 
      AM=AM-6.D1
      IH=IH+1 
   12 IF(IH.GE.24) IH=IH-24 
      IF(DEC) 15,20,20
   15 ISGN=IMINUS 
      GO TO 25
   20 ISGN=IPLUS
   25 DEG=DABS(DEC)/RAD 
      ID=DEG
      AMIN=(DEG-ID)*6.D1
      RETURN
      END 
C******************************************************************** 
c     SUBROUTINE AZALT(RU,RA,DEC,AZ,ALT)
C 
C   THIS SUBROUTINE CONVERTS TOPOCENTRIC RIGHT ASCENSION AND DECLINATION
C   TO AZIMUTH AND ALTITUDE, USING AN ALGORITHM PRESENTED ON P.402 OF 
C   ESCOBAL,P., METHODS OF ORBIT DETERMINATION. 
C 
c     IMPLICIT DOUBLE PRECISION(A-H,O-Z)
c     implicit integer*4 (i-n)
c     DIMENSION STACOR(3),CL(3),CS(3),CE(3),CZ(3),CLH(3)
c     COMMON/TOPO/IOPT,AU,REM,STACOR,ETMUTC 
c     DATA FLAT/3.35233D-3/ 
c     DATA TWOPI/6.283185307179586D0/ 
c     FLDFSQ=(1.D0-FLAT)**2 
c     TANPHI=STACOR(3)/STACOR(2)/FLDFSQ 
c     GEDLAT=DATAN(TANPHI)
c     COSGDL=DCOS(GEDLAT) 
c     SINGDL=DSIN(GEDLAT) 
c     COSRU=DCOS(RU)
c     SINRU=DSIN(RU)
c     CL(1)=DCOS(DEC)*DCOS(RA)
c     CL(2)=DCOS(DEC)*DSIN(RA)
c     CL(3)=DSIN(DEC) 
c     CS(1)=SINGDL*COSRU
c     CS(2)=SINGDL*SINRU
c     CS(3)=-COSGDL 
c     CE(1)=-SINRU
c     CE(2)=COSRU 
c     CE(3)=0.D0
c     CZ(1)=COSRU*COSGDL
c     CZ(2)=SINRU*COSGDL
c     CZ(3)=SINGDL
c     DO 20 I=1,3 
c  20 CLH(I)=0.D0 
c     DO 30 I=1,3 
c     CLH(1)=CLH(1)+CS(I)*CL(I) 
c     CLH(2)=CLH(2)+CE(I)*CL(I) 
c  30 CLH(3)=CLH(3)+CZ(I)*CL(I) 
c     COSH=DSQRT(1.D0-CLH(3)**2)
c     ALT=DATAN2(CLH(3),COSH) 
c     COSA=-CLH(1)/COSH 
c     SINA=CLH(2)/COSH
c     AZ=DATAN2(SINA,COSA)
c     IF(AZ.LT.0) AZ=AZ+TWOPI 
c     RETURN
c     END 
C******************************************************************** 
      SUBROUTINE DAYCNV(XJD,I,J,K,F)
C 
C**** CONVERTS JULIAN DATES TO GREGORIAN CLAENDAR DATES 
C 
C**** JD = JULIAN DATE
C**** I = YEAR
C**** J = MONTH 
C**** K = DAY 
C 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      implicit integer*4 (i-n)
      JD1=XJD 
      TEST=XJD-JD1
      JD=JD1+1
      IF(TEST.LT. 0.5) JD=JD1 
      FRAC=XJD-0.5-JD1
      IF(TEST.LT. 0.5) FRAC=XJD+0.5-JD1 
      F=FRAC*24.0 
      L=JD+68569
      N=4*L/146097
      L=L-(146097*N+3)/4
      I=4000*(L+1)/1461001
      L=L-1461*I/4+31 
      J=80*L/2447 
      K=L-2447*J/80 
      L=J/11
      J=J+2-12*L
      I=100*(N-49)+I+L
      RETURN
      END 
C******************************************************************** 
      SUBROUTINE DIABER(STACOR,RU,RA,DEC) 
C 
C   THIS SUBROUTINE APPLIES THE DIURNAL PORTION OF THE STELLAR
C   ABERRATION TO THE RIGHT ASCENSION AND DECLINATION.
C   REFERENCE - P.49 OF THE EXPLANATORY SUPPLEMENT TO THE AENA. 
C 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION STACOR(3) 
      DATA TWOPI/6.283185307179586D0/ 
      HAL=RU-RA 
      IF(HAL.LT.0.D0) HAL=HAL+TWOPI 
      DRA=1.5490D-6*STACOR(2)*DCOS(HAL)/DCOS(DEC) 
      DDEC=1.5514D-6*STACOR(2)*DSIN(HAL)*DSIN(DEC)
      RA=RA+DRA 
      DEC=DEC+DDEC
      RETURN
      END 
C******************************************************************** 
      SUBROUTINE ETERMS(RA,DEC) 
C 
C   THIS SUBROUTINE APPLIES THE ELLIPTIC PORTION (E-TERMS) OF THE ANNUAL
C   PART OF STELLAR ABERRATION TO THE RIGHT ASCENSION AND DECLINATION.
C   REFERENCE - P.144 OF THE EXPLANATORY SUPPLEMENT TO THE AENA.
C 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      implicit integer*4 (i-n)
      DRA=-1.653215D-6*DSIN(2.945243113D0+RA)/DCOS(DEC) 
      DDEC=-1.653215D-6*DCOS(2.945243113D0+RA)*DSIN(DEC)
     *-1.40596D-7*DCOS(DEC) 
      RA=RA+DRA 
      DEC=DEC+DDEC
      RETURN
      END 
C******************************************************************** 
      SUBROUTINE HELXYZ(TOBS,TT,Q,ECR,PX,PY,PZ,QX,QY,QZ,XG,YG,ZG) 
C 
C   THIS SUBROUTINE COMPUTES HELIOCENTRIC EQUATORIAL RECTANGULAR COORD- 
C   INATES FOR THE EQUATOR AND EQUINOX OF 1950.0 FROM COMETARY TWO-BODY 
C   ELEMENTS AND THE P AND Q ROTATION FACTORS.
C 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      implicit integer*4 (i-n)
      DATA XK/0.01720209895D0/
      EXCR=ECR-1.D0 
      IF(EXCR) 85,950,951 
  950 CONTINUE
C 
C   PARABOLIC CASE
      R2=2.0**(1.0/2.0) 
      Q32=Q**(3.0/2.0)
      XN=XK*(TOBS-TT)/(R2*Q32)
      KOUT=1
      SIG=1.0 
   70 SIG2=SIG**2 
      SIG3=SIG*SIG2 
      SIGMA=(XN+(2.0/3.0)*SIG3)/(1.0+SIG2)
      DEL=DABS(SIGMA)-DABS(SIG) 
      KOUT=KOUT+1 
      IF(KOUT.GT.10) GO TO 80 
      SIG=SIGMA 
      GO TO 70
   80 CONTINUE
      S=SIGMA 
      S2=SIGMA**2 
      XG=Q*PX*(1.0-S2)+2.0*Q*QX*S 
      YG=Q*PY*(1.0-S2)+2.0*Q*QY*S 
      ZG=Q*PZ*(1.0-S2)+2.0*Q*QZ*S 
      RETURN
  951 CONTINUE
C 
C   HYPERBOLIC CASE 
      AN1=Q/(1.D0-ECR)
      AN2=-AN1
      AN2=AN2**(3.D0/2.D0)
      AN3=XK/AN2
      ANON=AN3*(TOBS-TT)
      FHV=ANON
  953 FHVN=FHV-(ECR*DSINH(FHV)-FHV-ANON)/(ECR*DCOSH(FHV)-1.D0)
      TST=FHVN-FHV
      IF(DABS(TST).LT.1.D-13) GO TO 952 
      FHV=FHVN
      GO TO 953 
  952 CONTINUE
      XOB=AN1*(DCOSH(FHV)-ECR)
      AN4=ECR**2-1.D0 
      YOB=-AN1*DSQRT(AN4)*DSINH(FHV)
      XG=XOB*PX+YOB*QX
      YG=XOB*PY+YOB*QY
      ZG=XOB*PZ+YOB*QZ
      RETURN
C 
C   ELLIPTIC CASE 
   85 A=Q/(1.0-ECR) 
      ART=DSQRT(A)
      AR=XK*(TOBS-TT)/(ART**3)
      I=1 
      EE1=AR
      EE=EE1
   86 DE=EE1-EE+ECR*DSIN(EE)
      IF(I.GE. 10) GO TO 87 
      I=I+1 
      EE=EE+DE/(1.0D0-ECR*DCOS(EE +0.5D0*DE)) 
      GO TO 86
   87 CONTINUE
      B=A*DSQRT(1.0D0-ECR**2) 
      FAT1=DCOS(EE)-ECR 
      FAT2=DSIN(EE) 
      XG=A*PX*FAT1+B*QX*FAT2
      YG=A*PY*FAT1+B*QY*FAT2
      ZG=A*PZ*FAT1+B*QZ*FAT2
      RETURN
      END 
C******************************************************************** 
      SUBROUTINE NUTANG(TJD,DPSI,DEPS)
C 
C   THIS SUBROUTINE COMPUTES NUTATION ANGLES. 
C 
C   INPUT - 
C      TJD - EPHEMERIS TIME OF OBSERVATION (JED). 
C 
C   OUTPUT -
C      DPSI - DELTA PSI (RADIANS).
C      DEPS - DELTA EPSILON (RADIANS).
C 
C   REF. - J. MEEUS, ASTRONOMICAL FORMULAE FOR CALCULATORS, MONOGRAFIEEN
C   OVER ASTRONOMIE EN ASTROFYSICA, UITGEGEVEN DOOR VOLKSSTERRENWACHT 
C   URANIA V.Z.W. EN VERENIGING VOOR STERRENKUNDE V.Z.W. VOL. 4, 1979.
C 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      implicit integer*4 (i-n)
      DATA ASPR/206264.806247096355D0/
      DATA DPR/57.2957795130823252D0/ 
      TJC=(TJD-2415020.0D0)/36525.D0
      CMLSUN=(279.6967D0+(36000.7689D0+3.03D-4*TJC)*TJC)/DPR
      CMLMN=(270.4342D0+(481267.8831D0-1.133D-3*TJC)*TJC)/DPR 
      CMASUN=(358.4758D0+(35999.0498D0-1.5D-4*TJC)*TJC)/DPR 
      CMAMN=(296.1046D0+(477198.8491D0+9.192D-3*TJC)*TJC)/DPR 
      CLANMN=(259.1833D0-(1934.1420D0-2.078D-3*TJC)*TJC)/DPR
      DPSI=-(17.2327D0+1.737D-2*TJC)*DSIN(CLANMN) 
     *-(1.2729D0+1.3D-4*TJC)*DSIN(2.D0*CMLSUN)
     *+2.088D-1*DSIN(2.D0*CLANMN)-2.037D-1*DSIN(2.D0*CMLMN) 
     *+(1.261D-1-3.1D-4*TJC)*DSIN(CMASUN)+6.75D-2*DSIN(CMAMN) 
     *-(4.97D-2-1.2D-4*TJC)*DSIN(2.D0*CMLSUN+CMASUN)
     *-3.42D-2*DSIN(2.D0*CMLMN-CLANMN)-2.61D-2*DSIN(2.D0*CMLMN+CMAMN) 
     *+2.14D-2*DSIN(2.D0*CMLSUN-CMASUN) 
     *-1.49D-2*DSIN(2.D0*CMLSUN-2.D0*CMLMN+CMAMN) 
     *+1.24D-2*DSIN(2.D0*CMLSUN-CLANMN)+1.14D-2*DSIN(2.D0*CMLMN-CMAMN)
      DEPS=(9.21D0+9.1D-4*TJC)*DCOS(CLANMN) 
     *+(5.522D-1-2.9D-4*TJC)*DCOS(2.D0*CMLSUN)-9.04D-2*DCOS(2.D0*CLANMN)
     *+8.84D-2*DCOS(2.D0*CMLMN)+2.16D-2*DCOS(2.D0*CMLSUN+CMASUN)
     *+1.83D-2*DCOS(2.D0*CMLMN-CLANMN)+1.13D-2*DCOS(2.D0*CMLMN+CMAMN) 
     *-9.3D-3*DCOS(2.D0*CMLSUN-CMASUN)-6.6D-3*DCOS(2.D0*CMLSUN-CLANMN)
      DPSI=DPSI/ASPR
      DEPS=DEPS/ASPR
      RETURN
      END 
C******************************************************************** 
      SUBROUTINE NUTATE(DJ,PUT,OBM,FRWD,X)
C 
C     FRWD = TRUE IMPLIES NUTATION FROM 1950,  = FALSE IMPLIES TO 1950
C 
C     REF. - MELBOURNE, W.G., ET.AL., CONSTANTS AND RELATED INFORMATION 
C     FOR ASTRODYNAMIC CALCULATIONS, 1968, JPL TR 32-1306, 1968, P. 9.
C 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      implicit integer*4 (i-n)
      DIMENSION OBM(2),OBT(2),DPSI(2),X(6),PUT(4) 
      LOGICAL FRWD
      DATA SPR/206264.8062D0/ 
      T=(DJ-2433282.5D0)/36525.D0 
C----  CONSTANTS DOCUMENT   TR32-1306  1968  PG. 9
      OBM(1)=(84404.84D0+T*(-46.850D0+T*(-0.0034D0+T*0.0018D0)))/SPR
      OBM(2)=(-46.850D0-0.0068D0*T+0.0054D0*T*T)/(SPR*36525.D0) 
      DPSI(1)=PUT(1)
      DPSI(2)=PUT(3)
      OBT(1)=OBM(1)+PUT(2)
      OBT(2)=OBM(2)+PUT(4)
      IF(.NOT. FRWD) GO TO 1
      CALL ROTAT6(X,OBM,+1,X) 
      CALL ROTAT6(X,DPSI,-3,X)
      CALL ROTAT6(X,OBT,-1,X) 
      RETURN
   1  CALL ROTAT6(X,OBT,+1,X) 
      CALL ROTAT6(X,DPSI,+3,X)
      CALL ROTAT6(X,OBM,-1,X) 
      RETURN
      END 
C******************************************************************** 
      SUBROUTINE PRECES(TZ,TF,FRWD,X) 
C ----  REFERENCE  -  LIESKE  JPL TECH REPORT  TR32-1044
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      implicit integer*4 (i-n)
      DIMENSION ZE(2),TH(2),Z(2),X(3) 
      DOUBLE PRECISION JD1900 
      LOGICAL FRWD
      DATA JD1900/2415020.313516463D0/
      DATA TCENT/36524.219879D0/
      DATA SPR/206264.8062D0/ 
      SCALE=SPR*TCENT 
      IF(TZ .EQ. TZZ) GO TO 2 
      T1=(TZ-JD1900)/TCENT
      TZZ=TZ
      P1=2304.253D0+T1*(1.3972D0+T1*0.125D-3) 
      P2=0.3023D0-0.211D-3*T1 
      P3=0.0180D0 
      Q1=P1 
      Q2=1.0949D0+0.460D-3*T1 
      Q3=0.0183D0 
      R1=2004.684D0+T1*(-.8532D0-T1*0.317D-3) 
      R2=-0.4266D0-0.320D-3*T1
      R3=-0.0418D0
   2  T=(TF-TZ)/TCENT 
      ZE(1)=T*(P1+T*(P2+T*P3))/SPR
      ZE(2)=(P1+T*(2.D0*P2+T*3.D0*P3))/SCALE
      Z (1)=T*(Q1+T*(Q2+T*Q3))/SPR
      Z (2)=(Q1+T*(2.D0*Q2+T*3.D0*Q3))/SCALE
      TH(1)=T*(R1+T*(R2+T*R3))/SPR
      TH(2)=(R1+T*(2.D0*R2+T*3.D0*R3))/SCALE
      IF(FRWD) GO TO 1
C ----  HERE FOR PRECESSING FROM TF TO TZ 
      CALL ROTAT6(X,Z ,+3,X)
      CALL ROTAT6(X,TH,-2,X)
      CALL ROTAT6(X,ZE,+3,X)
      RETURN
C ----  HERE FOR PRECESSING FROM TZ TO TF 
   1  CALL ROTAT6(X,ZE,-3,X)
      CALL ROTAT6(X,TH,+2,X)
      CALL ROTAT6(X,Z ,-3,X)
      RETURN
      END 
C******************************************************************** 
      SUBROUTINE ROTAT6(X,A,L,Y)
      DOUBLE PRECISION X(6),A(2),Y(6),C,S,T,SIGN
      K=IABS(L) 
      I=MOD(K,3)+1
      J=6-I-K 
      SIGN=1.D0 
      IF(L .LT. 0) SIGN=-1.D0 
      C=DCOS(A(1))
      S=DSIN(A(1)*SIGN) 
      T=X(I+3)*C+X(J+3)*S+SIGN*A(2)*(X(J)*C-X(I)*S) 
      Y(J+3)=X(J+3)*C-X(I+3)*S-SIGN*A(2)*(X(I)*C+X(J)*S)
      Y(I+3)=T
      Y(K+3)=X(K+3) 
      T=X(I)*C+X(J)*S 
      Y(J)=X(J)*C-X(I)*S
      Y(I)=T
      Y(K)=X(K) 
      RETURN
      END 
C******************************************************************** 
      SUBROUTINE SUN (DJ,EL,C,XSUN,YSUN,ZSUN) 
C 
C                     NEWCOMB'S THEORY OF THE EARTH 
C 
C         FROM THE 'TABLES OF THE SUN' (A.P.A.E. VI, PART 1), NEWCOMB'S 
C     THEORY OF THE MOTION OF THE OBJECT IS EVALUATED, GIVING GEOMETRIC 
C     COORDINATES REFERRED TO THE MEAN EQUINOX AND EQUATOR OF DATE AS A 
C     FUNCTION OF EPHEMERIS TIME EXPRESSED IN JULIAN CENTURIES OF 36525 
C     DAYS. 
C 
C     THIS PROGRAM ORIGINATED AT THE U.S. NAVAL OBSERVATORY 
C******************************************************************** 
C 
C     DEFINITIONS 
C 
C     DJ    = INPUT ARGUMENT TIME IN JULIAN EPHEMERIS DAYS
C 
C     EL( 2)= ORBITAL ECCENTRICITY
C     EL( 5)= LONGITUDE OF PERIGEE
C     EL( 8)= MEAN SIDEREAL MOTION                      NOT GIVEN 
C     EL( 9)= UNPERTURBED MEAN LONGITUDE, RADIANS 
C     EL(10)= MEAN ANOMALY, AFFECTED BY LONG PERIOD PERTURBATIONS,
C             RADIANS 
C     EL(11)= UNPERTURBED RADIUS, ASTRONOMICAL UNITS
C     EL(12)= EQUATION OF THE CENTER, RADIANS 
C     EL(14)= MEAN LONGITUDE OF MOON
C     EL(15)= MEAN ANOMALY OF MOON
C     EL(16)= LUNAR MEAN ARGUMENT OF LATITUDE 
C     EL(17)= MEAN LONGITUDE OF LUNAR ASCENDING NODE
C     EL(21)= MEAN LONGITUDE OF MOON'S PERIGEE
C     REMAINING ELEMENTS OF EL ARRAY NOT USED 
C 
C     C( 1) = PERTURBED RADIUS VECTOR, ASTRONOMICAL UNITS 
C     C( 2) = TRUE LONGITUDE, DEGREES 
C     C( 3) = LATITUDE, DEGREES 
C     C( 4) = ECLIPTIC LONGITUDE MEAN EQUINOX OF DATE, RADIANS
C     C( 5) = ECLIPTIC LATITUDE MEAN EQUINOX AND EQUATOR OF DATE, RADIANS 
C     REMAINING ELEMENTS OF C ARRAY NOT USED
C     EL( 1)= SEMI-MAJOR AXIS 
C 
C 
C******************************************************************** 
C 
C     SOME NECESSARY PRELIMINARIES
C 
      DOUBLE PRECISION DJ,T,TP,RO,GME,GV,GMA,GJ,GS,DL,DR,DBLARG,DB,DG,
     1D,T18,EL(21),C(16),TWOPI,STR,RTD
      DOUBLE PRECISION S,EA,EB,EQ,ED,EPS,XSUNE,YSUNE,ZSUNE, 
     *XSUN,YSUN,ZSUN,FOUR,CONV11,CONV12 
      DIMENSION X(8,121),X1(160),X2(160),X3(160),X4(160),X5(160),X6(8)
      DIMENSION X0(160) 
      EQUIVALENCE(X(1,  1),X0(1)) 
      EQUIVALENCE(X(1, 21),X1(1)) 
      EQUIVALENCE(X(1, 41),X2(1)) 
      EQUIVALENCE(X(1, 61),X3(1)) 
      EQUIVALENCE(X(1, 81),X4(1)) 
      EQUIVALENCE(X(1,101),X5(1)) 
      EQUIVALENCE(X(1,121),X6(1)) 
      DATA TWOPI/6.283185307179586D0/ 
      DATA STR/206264806.2470964D0/ 
      DATA RTD/57.295779513082321D0/
C 
C 
C******************************************************************** 
C 
C     PART I    TABLES OF THE PERTURBATIONS 
C 
C     PERTURBATIONS BY MERCURY
C 
C                  J    I     VC      VS    RHOC    RHOS      BC     BS 
      DATAX0/    - 1.,+ 1., -   6.,-  11., +  26.,-  12.,     0.,    0.,
     1           - 1.,+ 2., -   3.,-   3., -   4.,+   5.,     0.,    0.,
     2           - 1.,+ 3., +  15.,-   1., -   1.,-  18.,     0.,    0.,
     3           - 1.,+ 4., +  19.,-  13., -   2.,-   4.,     0.,    0.,
C 
C 
C     PERTURBATIONS BY VENUS
C 
C                  J    I     VC      VS    RHOC    RHOS      BC     BS 
     4           - 1.,  0., +  33.,-  67., -  85.,-  39., +  24.,-  17.,
     5           - 1.,+ 1., +2353.,-4228., -2062.,-1146., -   4.,+   3.,
     6           - 1.,+ 2., -  65.,-  34., +  68.,-  14., +   6.,-  92.,
     7           - 1.,+ 3., -   3.,-   8., +  14.,-   8., +   1.,+   7.,
     8           - 2.,  0.,     0.,    0.,     0.,+   4.,     0.,    0.,
     9           - 2.,+ 1., -  99.,+  60., +  84.,+ 136., +  23.,-   3.,
     1           - 2.,+ 2., -4702.,+2903., +3593.,+5822., +  10.,-   6.,
     2           - 2.,+ 3., +1795.,-1737., - 596.,- 632., +  37.,-  56.,
     3           - 2.,+ 4., +  30.,-  33., +  40.,+  33., +   5.,-  13.,
     4           - 3.,+ 2., -  13.,    0.,     0.,+  21., +  13.,+   5.,
     5           - 3.,+ 3., - 666.,+  27., +  44.,+1044., +   8.,+   1.,
     6           - 3.,+ 4., +1508.,- 397., - 381.,-1448., + 185.,- 100.,
     7           - 3.,+ 5., + 763.,- 684., + 126.,+ 148., +   6.,-   3.,
     8           - 3.,+ 6., +  12.,-  12., +  14.,+  13., -   2.,+   4.,
     9           - 4.,+ 3.,     0.,    0.,     0.,+   6., +   4.,+   5.,
     1           - 4.,+ 4., - 188.,-  93., - 166.,+ 337.,     0.,    0./
      DATAX1/    - 4.,+ 5., - 139.,-  38., -  51.,+ 189., -  31.,-   1.,
     1           - 4.,+ 6., + 146.,-  42., -  25.,-  91., +  12.,    0.,
     2           - 4.,+ 7., +   5.,-   4., +   3.,+   5.,     0.,    0.,
     3           - 5.,+ 5., -  47.,-  69., - 134.,+  93.,     0.,    0.,
     4           - 5.,+ 6., -  28.,-  25., -  39.,+  43., -   8.,-   4.,
     5           - 5.,+ 7., - 119.,-  33., -  37.,+ 136., -  18.,-   6.,
     6           - 5.,+ 8., + 154.,    0.,     0.,-  26.,     0.,    0.,
     7           - 6.,+ 5.,     0.,    0.,     0.,    0., -   2.,+   6.,
     8           - 6.,+ 6., -   4.,-  38., -  80.,+   8.,     0.,    0.,
     9           - 6.,+ 7., -   4.,-  13., -  24.,+   7., -   2.,-   3.,
     1           - 6.,+ 8., -   6.,-   7., -  10.,+  10., -   2.,-   3.,
     2           - 6.,+ 9., +  14.,+   3., +   3.,-  12.,     0.,    0.,
     3           - 7.,+ 7., +   8.,-  18., -  38.,-  17.,     0.,    0.,
     4           - 7.,+ 8., +   1.,-   6., -  12.,-   3.,     0.,    0.,
     5           - 7.,+ 9., +   1.,-   3., -   4.,+   1.,     0.,    0.,
     6           - 7.,+10., -   2.,-   2., -   3.,+   3.,     0.,    0.,
     7           - 8.,+ 8., +   9.,-   7., -  14.,-  19.,     0.,    0.,
     8           - 8.,+ 9., +   2.,-   3., -   5.,-   4.,     0.,    0.,
     9           - 8.,+12., -  33.,-  54., -  43.,+   8., -   5.,-   9.,
     1           - 8.,+13.,     0.,    0., -   9.,-   8.,     0.,    0./
      DATAX2/    - 8.,+14.,     0.,    0., -  25.,+  22.,     0.,    0.,
     1           - 9.,+ 9., +   6.,-   1., -   2.,-  13.,     0.,    0.,
     2           - 9.,+10.,     0.,    0., -   1.,-   4.,     0.,    0.,
     3           -10.,+10.,     0.,    0., +   3.,-   7.,     0.,    0.,
C 
C 
C     PERTURBATIONS BY MARS 
C 
C                  J    I     VC      VS    RHOC    RHOS      BC     BS 
     4           + 1.,- 2., -   5.,-   4., -   5.,+   6.,     0.,    0.,
     5           + 1.,- 1., - 216.,- 167., -  92.,+ 119.,     0.,    0.,
     6           + 1.,  0., -   8.,-  47., +  27.,-   6.,     0.,    0.,
     7           + 2.,- 3., +  40.,-  10., -  13.,-  50.,     0.,    0.,
     8           + 2.,- 2., +1963.,- 567., - 573.,-1976.,     0.,-   8.,
     9           + 2.,- 1., -1659.,- 617., +  64.,- 137.,     0.,    0.,
     1           + 2.,  0., -  24.,+  15., -  18.,-  25., -   8.,+   2.,
     2           + 3.,- 4., +   1.,-   4., -   6.,    0.,     0.,    0.,
     3           + 3.,- 3., +  53.,- 118., - 154.,-  67.,     0.,    0.,
     4           + 3.,- 2., + 396.,- 153., -  77.,- 201.,     0.,    0.,
     5           + 3.,- 1., +   8.,+   1.,     0.,+   6.,     0.,    0.,
     6           + 4.,- 4., +  11.,+  32., +  46.,-  17.,     0.,    0.,
     7           + 4.,- 3., - 131.,+ 483., + 461.,+ 125., +   7.,+   1.,
     8           + 4.,- 2., + 526.,- 256., +  43.,+  96.,     0.,    0.,
     9           + 4.,- 1., +   7.,-   5., +   6.,+   8.,     0.,    0.,
     1           + 5.,- 5., -   7.,+   1.,     0.,+  12.,     0.,    0./
      DATAX3/    + 5.,- 4., +  49.,+  69., +  87.,-  62.,     0.,    0.,
     1           + 5.,- 3., -  38.,+ 200., +  87.,+  17.,     0.,    0.,
     2           + 5.,- 2., +   3.,+   1., -   1.,+   3.,     0.,    0.,
     3           + 6.,- 6.,     0.,    0., -   4.,-   3.,     0.,    0.,
     4           + 6.,- 5., -  20.,-   2., -   3.,+  30.,     0.,    0.,
     5           + 6.,- 4., - 104.,- 113., - 102.,+  94.,     0.,    0.,
     6           + 6.,- 3., -  11.,+ 100., -  27.,-   4.,     0.,    0.,
     7           + 7.,- 6., +   3.,-   5., -   9.,-   5.,     0.,    0.,
     8           + 7.,- 5., -  49.,+   3., +   4.,+  60.,     0.,    0.,
     9           + 7.,- 4., -  78.,-  72., -  26.,+  28.,     0.,    0.,
     1           + 8.,- 7., +   1.,+   3., +   5.,-   1.,     0.,    0.,
     2           + 8.,- 6., +   6.,-   8., -  12.,-   9.,     0.,    0.,
     3           + 8.,- 5., +  51.,-  10., -   8.,-  44.,     0.,    0.,
     4           + 8.,- 4., -  17.,-  12., +   5.,-   6.,     0.,    0.,
     5           + 9.,- 7., +   2.,+   3., +   5.,-   3.,     0.,    0.,
     6           + 9.,- 6., +  13.,-  25., -  30.,-  16.,     0.,    0.,
     7           + 9.,- 5., +  60.,-  15., -   4.,-  17.,     0.,    0.,
     8           +10.,- 7., +   2.,+   5., +   7.,-   3.,     0.,    0.,
     9           +10.,- 6., -   7.,+  18., +  14.,+   6.,     0.,    0.,
     1           +10.,- 5., +   5.,-   2.,     0.,    0.,     0.,    0./
      DATAX4/    +11.,- 7., +   9.,+  15., +  17.,-  10.,     0.,    0.,
     1           +11.,- 6., -  12.,+  42., +   8.,+   3.,     0.,    0.,
     2           +12.,- 7., -   4.,-   5., -   4.,+   3.,     0.,    0.,
     3           +13.,- 8., -  13.,-   1., -   1.,+  15.,     0.,    0.,
     4           +13.,- 7., -  30.,-  33., -   4.,+   3.,     0.,    0.,
     5           +15.,- 9., +  13.,-  16., -  17.,-  14.,     0.,    0.,
     6           +15.,- 8., + 200.,-  30., -   1.,-   6.,     0.,    0.,
     7           +17.,-10., -   2.,-   4., -   4.,+   2.,     0.,    0.,
     8           +17.,- 9., -  10.,+  24.,     0.,    0.,     0.,    0.,
C 
C 
C     PERTURBATIONS BY JUPITER
C 
C                  J    I     VC      VS    RHOC    RHOS      BC     BS 
     9           + 1.,- 3., -   3.,-   1., -   2.,+   5.,     0.,    0.,
     1           + 1.,- 2., - 155.,-  52., -  78.,+ 193., +   7.,    0.,
     2           + 1.,- 1., -7208.,+  59., +  56.,+7067., -   1.,+  17.,
     3           + 1.,  0., - 307.,-2582., + 227.,-  89., +  16.,    0.,
     4           + 1.,+ 1., +   8.,-  73., +  79.,+   9., +   1.,+  23.,
     5           + 2.,- 3., +  11.,+  68., + 102.,-  17.,     0.,    0.,
     6           + 2.,- 2., + 136.,+2728., +4021.,- 203.,     0.,    0.,
     7           + 2.,- 1., - 537.,+1518., +1376.,+ 486., +  13.,+ 166.,
     8           + 2.,  0., -  22.,-  70., -   1.,-   8.,     0.,    0.,
     9           + 3.,- 4., -   5.,+   2., +   3.,+   8.,     0.,    0.,
     1           + 3.,- 3., - 162.,+  27., +  43.,+ 278.,     0.,    0./
      DATAX5/    + 3.,- 2., +  71.,+ 551., + 796.,- 104., +   6.,-   1.,
     1           + 3.,- 1., -  31.,+ 208., + 172.,+  26., +   1.,+  18.,
     2           + 4.,- 4., -   3.,-  16., -  29.,+   5.,     0.,    0.,
     3           + 4.,- 3., -  43.,+   9., +  13.,+  73.,     0.,    0.,
     4           + 4.,- 2., +  17.,+  78., + 110.,-  24.,     0.,    0.,
     5           + 4.,- 1., -   1.,+  23., +  17.,+   1.,     0.,    0.,
     6           + 5.,- 5.,     0.,    0., -   1.,-   3.,     0.,    0.,
     7           + 5.,- 4., -   1.,-   5., -  10.,+   2.,     0.,    0.,
     8           + 5.,- 3., -   7.,+   2., +   3.,+  12.,     0.,    0.,
     9           + 5.,- 2., +   3.,+   9., +  13.,-   4.,     0.,    0.,
C 
C 
C     PERTURBATIONS BY SATURN 
C 
C                  J    I     VC      VS    RHOC    RHOS      BC     BS 
     1           + 1.,- 2., -   3.,+  11., +  15.,+   3.,     0.,    0.,
     2           + 1.,- 1., -  77.,+ 412., + 422.,+  79., +   1.,+   6.,
     3           + 1.,  0., -   3.,- 320., +   8.,-   1.,     0.,    0.,
     4           + 1.,+ 1.,     0.,-   8., +   8.,    0., -   1.,+   6.,
     5           + 2.,- 3.,     0.,    0., -   3.,-   1.,     0.,    0.,
     6           + 2.,- 2., +  38.,- 101., - 152.,-  57.,     0.,    0.,
     7           + 2.,- 1., +  45.,- 103., - 103.,-  44.,     0.,    0.,
     8           + 2.,  0., +   2.,-  17.,     0.,    0.,     0.,    0.,
     9           + 3.,- 2., +   7.,-  20., -  30.,-  11.,     0.,    0.,
     1           + 3.,- 1., +   6.,-  16., -  16.,-   6.,     0.,    0./
      DATAX6/    + 4.,- 2., +   1.,-   3., -   4.,-   1.,     0.,    0./
C 
C 
C******************************************************************** 
C 
C     PART II   REDUCTION OF ARGUMENT TIME
C 
C     T  = TIME IN JULIAN CENTURIES FROM 1900 JANUARY 0, GREENWICH MEAN NOON
      T  = (DJ - 2415020.D0)/36525.D0 
C 
C     TP = TIME IN JULIAN YEARS     FROM 1850 JANUARY 0, GREENWICH MEAN NOON
      TP = (DJ - 2396758.D0)/365.25D0 
C 
C     T18= TIME IN JULIAN CENTURIES FROM 1800 JANUARY 0, GREENWICH MEAN NOON
      T18= (DJ - 2378496.D0)/36525.D0 
C 
C 
C******************************************************************** 
C 
C     PART III  COMPUTATION OF ELLIPTIC ELEMENTS AND SECULAR TERMS
C 
      FOUR = 0.6283019457267408D+03 + 4.78D0/STR
      EL( 2) =  0.1675104D-01 - T*(0.4180D-04 + 0.126D-06*T)
      EL( 5) =  0.4908229466868880D+01 + (0.3000526416797348D-01
     1        +(0.7902463002085429D-05 +  0.5817764173314427D-07*T)*T)*T
      EL( 9) =  0.4881627934111871D+01 + (0.6283319509909086D+03
     1        + 0.5279620987282842D-05*T)*T 
      EL(10) =  0.6256583580497099D+01 + (  FOUR
     1        -(0.2617993877991491D-05 +  0.5817764173314427D-07*T)*T)*T
C 
C     REDUCTION TO THE FIRST REVOLUTION 
      EL( 9) = DMOD (EL( 9),TWOPI)
      EL(10) = DMOD (EL(10),TWOPI)
C 
C 
C 
C******************************************************************** 
C 
C     PART IV   LUNAR TERMS 
C 
C     FUNDAMENTAL ARGUMENTS OF HANSEN WITH NEWCOMB CORRECTIONS
      EL(14) =  0.5859485105530519D+01 + (0.8399709200339593D+04
     1        +(0.4619304753611657D-04 +  0.6544984694978728D-07*T18) 
     2        *T18)*T18 
      EL(17) =  0.5807638839584459D+00 - (0.3375728238223387D+02
     1        -(0.3964806284113784D-04 +  0.3470781143063166D-07*T18) 
     2        *T18)*T18 
      EL(21) =  0.3933938487925745D+01 + (0.7101833330913285D+02
     1        -(0.1752456013106637D-03 -  0.1773109075921903D-06*T18) 
     2        *T18)*T18 
      EL(14) = DMOD (EL(14),TWOPI)
      EL(17) = DMOD (EL(17),TWOPI)
      EL(21) = DMOD (EL(21),TWOPI)
C 
C     DERIVED ARGUMENTS 
      EL(15) = EL(14) - EL(21)
      EL(16) = EL(14) - EL(17)
      EL(15) = DMOD (EL(15),TWOPI)
      EL(16) = DMOD (EL(16),TWOPI)
C     MEAN ELONGATION 
      D      = EL(14) - EL(9) 
C 
C     COMBINATIONS OF ARGUMENTS AND THE PERTURBATIONS 
      D = DMOD (D,TWOPI)
      ARG = D 
      DL =    + 6454.*SIN(ARG) + 13.*SIN(3.*ARG)
      DR =    + 13360.*COS(ARG) +  30.*COS(3.*ARG)
C 
      DBLARG = D + EL(15) 
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DL = DL +  177.*SIN(ARG)
      DR = DR +  370.*COS(ARG)
C 
      DBLARG = D - EL(15) 
      DBLARG = DMOD(DBLARG,TWOPI) 
      ARG = DBLARG
      DL = DL -  424.*SIN(ARG)
      DR = DR - 1330.*COS(ARG)
C 
      DBLARG = 3.D0*D - EL(15)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DL = DL +   39.*SIN(ARG)
      DR = DR +   80.*COS(ARG)
C 
      DBLARG = D + EL(10) 
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DL = DL -   64.*SIN(ARG)
      DR = DR -  140.*COS(ARG)
C 
      DBLARG = D - EL(10) 
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DL = DL +  172.*SIN(ARG)
      DR = DR + 360.*COS(ARG) 
C 
      EL(16) = DMOD (EL(16),TWOPI)
      ARG = EL(16)
      DB =    + 576.*SIN(ARG) 
C 
      DBLARG = EL(16) + EL(15)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DB = DB +  16.*SIN(ARG) 
C 
      DBLARG = EL(16) - EL(15)
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DB = DB -  47.*SIN(ARG) 
C 
      DBLARG = EL(16) - 2.D0*(EL(9) - EL(17)) 
      DBLARG = DMOD (DBLARG,TWOPI)
      ARG = DBLARG
      DB = DB +  21.*SIN(ARG) 
C 
C 
C******************************************************************** 
C 
C     PART V    COMPUTATION OF PERIODIC PERTURBATIONS 
C 
C     THE PERTURBING MEAN ANOMALIES 
C 
      GME = 0.43296383D+01 + .2608784648D+02*TP 
      GV  = 0.19984020D+01 + .1021322923D+02*TP 
      GMA = 0.19173489D+01 + .3340556174D+01*TP 
      GJ  = 0.25836283D+01 + .5296346478D+00*TP 
      GS  = 0.49692316D+01 + .2132432808D+00*TP 
C 
C     REDUCTION TO FIRST REVOLUTION 
      GME = DMOD (GME,TWOPI)
      GV  = DMOD (GV ,TWOPI)
      GMA = DMOD (GMA,TWOPI)
      GJ  = DMOD (GJ ,TWOPI)
      GS  = DMOD (GS ,TWOPI)
C 
C 
C     MODIFICATION OF FUNDAMENTAL ARGUMENTS 
C 
C     APPLICATION OF THE JUPITER-SATURN GREAT INEQUALITY
C     TO JUPITER'S MEAN ANOMALY 
C 
      GJ = GJ + 0.579904067D-02 * DSIN(5.D0*GS - 2.D0*GJ + 1.1719644977 
     1     D0 - 0.397401726D-03 * TP) 
C 
C     LONG PERIOD PERTURBATIONS OF MEAN ANOMALY 
C 
C     TIME REPRESENTED IN SINGLE PRECISION BY 'ST'
      ST = T
C     ARGUMENT ( 4MARS - 7EARTH + 3VENUS ) OF SECOND ORDER AS TO MASSES 
      DG = 266.*SIN(.555015 + 2.076942*ST)
C 
C     ARGUMENT ( 3JUPITER - 8MARS + 4EARTH ) OF 2ND ORDER AS TO MASSES
     1    +6400.*SIN(4.035027 + .3525565*ST)
C 
C     ARGUMENT ( 13EARTH - 8VENUS ) FIRST ORDER IN MASSES, 5TH IN ECCENTRICITIES
C                                                          5TH IN INCLINATIONS
     2    +(1882.-16.*T)*SIN(.9990265 + 2.622706*ST)
C 
      GJ = DMOD (GJ,TWOPI)
C 
C 
C     COMPUTATION OF THE EQUATION OF THE CENTER 
C 
C     FORM PERTURBED MEAN ANOMALY 
      EL(10) = DG/STR + EL(10)
      EL(10) = DMOD(EL(10),TWOPI) 
      EL(12) =   DSIN(     EL(10)) * (6910057. - (17240. + 52.*T)*T)
     1         + DSIN(2.D0*EL(10)) * (  72338. -    361.*T) 
     2         + DSIN(3.D0*EL(10)) * (   1054. -      1.*T) 
     3         + DSIN(4.D0*EL(10)) *       18.
C 
C     THE UNPERTURBED RADIUS VECTOR 
      RO     =                          30570. -    150.*T
     1         - DCOS(     EL(10)) * (7274120. - (18140. + 50.*T)*T)
     2         - DCOS(2.D0*EL(10)) * (  91380. -    460.*T) 
     3         - DCOS(3.D0*EL(10)) * (   1450. -     10.*T) 
     4         - DCOS(4.D0*EL(10)) *       20.
      EL(11) = 10.D0**(RO*1.D-09) 
C 
C     PERTURBATIONS BY MERCURY IN LONGITUDE AND RADIUS VECTOR 
      DO 10 K=1,4 
C     ARGUMENT  J*MERCURY +   I * EARTH 
      DBLARG = X(1,K)*GME + X(2,K)*EL(10) 
      DBLARG = DMOD(DBLARG,TWOPI) 
      ARG = DBLARG
      CS  = COS(ARG)
      SS  = SIN(ARG)
      DL  =(X(3,K)*CS  + X(4,K)*SS )+ DL
      DR  =(X(5,K)*CS  + X(6,K)*SS )+ DR
   10 CONTINUE
C 
C     PERTURBATIONS BY VENUS
      DO 20 K=5,44
C     ARGUMENT J * VENUS +   I * EARTH
      DBLARG = X(1,K)*GV + X(2,K)*EL(10)
      DBLARG = DMOD(DBLARG,TWOPI) 
      ARG = DBLARG
      CS  = COS(ARG)
      SS  = SIN(ARG)
      DL  =(X(3,K)*CS  + X(4,K)*SS )+ DL
      DR  =(X(5,K)*CS  + X(6,K)*SS )+ DR
      DB  =(X(7,K)*CS  + X(8,K)*SS )+ DB
   20 CONTINUE
C 
C     PERTURBATIONS BY MARS 
      DO 30 K=45,89 
C     ARGUMENT  J * MARS  +   I * EARTH 
      DBLARG = X(1,K)*GMA + X(2,K)*EL(10) 
      DBLARG = DMOD(DBLARG,TWOPI) 
      ARG = DBLARG
      CS  = COS(ARG)
      SS  = SIN(ARG)
      DL  =(X(3,K)*CS  + X(4,K)*SS )+ DL
      DR  =(X(5,K)*CS  + X(6,K)*SS )+ DR
      DB  =(X(7,K)*CS  + X(8,K)*SS )+ DB
   30 CONTINUE
C 
C     PERTURBATIONS BY JUPITER
      DO 40 K=90,110
C     ARGUMENT J*JUPITER +   I * EARTH
      DBLARG = X(1,K)*GJ + X(2,K)*EL(10)
      DBLARG = DMOD(DBLARG,TWOPI) 
      ARG = DBLARG
      CS  = COS(ARG)
      SS  = SIN(ARG)
      DL  =(X(3,K)*CS  + X(4,K)*SS )+ DL
      DR  =(X(5,K)*CS  + X(6,K)*SS )+ DR
      DB  =(X(7,K)*CS  + X(8,K)*SS) +DB 
   40 CONTINUE
C 
C     PERTURBATIONS BY SATURN 
      DO 50 K=111,121 
C     ARGUMENT J*SATURN  +   I * EARTH
      DBLARG = X(1,K)*GS + X(2,K)*EL(10)
      DBLARG = DMOD(DBLARG,TWOPI) 
      ARG = DBLARG
      CS  = COS(ARG)
      SS  = SIN(ARG)
      DL  =(X(3,K)*CS  + X(4,K)*SS )+ DL
      DR  =(X(5,K)*CS  + X(6,K)*SS )+ DR
      DB  =(X(7,K)*CS  + X(8,K)*SS )+ DB
   50 CONTINUE
C 
C 
C******************************************************************** 
C 
C     PART VI   COMPUTATION OF ECLIPTIC COORDINATES 
C 
      C(1) = EL(11)*10.D0**(DR*1.D-09)
      C(4) = (DL + DG + EL(12))/STR + EL(9) 
      C(4) = DMOD(C(4),TWOPI) 
      C(5) = DB/STR 
      C(2) = C(4)*RTD 
      C(3) = C(5)*RTD 
C 
C 
C******************************************************************** 
C 
C     THE FOLLOWING COMPUTES THE RECTANGULAR COORDINATES
C     REFERRED TO THE MEAN EQUATOR AND EQUINOX OF DATE. 
C 
      CONV11=(TWOPI/360.D0)/3600.D0 
      CONV12=36524.219879D0-.00000614D0*T 
      S=(DJ-2415020.3135D0)/CONV12
      EA=(23.D0*3600.D0+27.D0*60.D0+8.26D0)*CONV11
      EB=(-46.845D0*S)*CONV11 
      EQ=(-.0059D0*S*S)*CONV11
      ED=(.00181D0*S*S*S)*CONV11
      EPS=EA+EB+EQ+ED 
      XSUNE=C(1)*DCOS(C(4))*DCOS(C(5))
      YSUNE=C(1)*DSIN(C(4))*DCOS(C(5))
      ZSUNE=C(1)*DSIN(C(5)) 
      XSUN=XSUNE
      YSUN=YSUNE*DCOS(EPS)-ZSUNE*DSIN(EPS)
      ZSUN=YSUNE*DSIN(EPS)+ZSUNE*DCOS(EPS)
      RETURN
      END 
C******************************************************************** 
      SUBROUTINE VELCTY(TOBS,H,RV,RCD)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      implicit integer*4 (i-n)
      DIMENSION T(5),RV(5),R(3),RCDV(2),TI(3),COEF(5) 
      DIMENSION TDI(5),TD(5,5)
      DO 20 I=1,5 
   20 T(I)=TOBS+(I-3)*H 
      DT=0.25D0 
      IF(DT.GT.H) DT=H
      TI(1)=TOBS-DT 
      TI(2)=TOBS
      TI(3)=TOBS+DT 
      DO 70 I=1,3 
      DO 35 K=1,5 
      TDI(K)=TI(I)-T(K) 
      DO 35 J=1,5 
   35 TD(J,K)=T(J)-T(K) 
      DO 60 L=1,5 
      PRODN=1.D0
      PRODD=1.D0
      DO 40 M=1,5 
      IF(M.EQ.L) GO TO 40 
      PRODN=PRODN*TDI(M)
   40 CONTINUE
      DO 50 M=1,5 
      IF(M.EQ.L) GO TO 50 
      PRODD=PRODD*TD(L,M) 
   50 CONTINUE
      COEF(L)=PRODN/PRODD 
   60 CONTINUE
      R(I)=0.D0 
      DO 65 N=1,5 
   65 R(I)=R(I)+COEF(N)*RV(N) 
   70 CONTINUE
      RCDV(1)=(R(2)-R(1))/DT
      RCDV(2)=(R(3)-R(2))/DT
      RCD=(RCDV(1)+RCDV(2))/2 
      RETURN
      END 
* 
      subroutine jd (djulda)
*     ======================
*     converts current UT form computer clock to Julian days
      implicit integer*4 (i-n)
      integer*2 it(5), iy 
      real*8   djulda 
* 
      call exec (11,it,iy)
      nyrm1 = iy - 1
      ic = nyrm1 / 100
      julda = 1721425 + 365*nyrm1 + nyrm1/4 - ic + ic/4 
      djulda = julda + (((dble(it(1))/100D0 + dble(it(2)))/60D0 
     *         + dble(it(3)))/60D0 + dble(it(4)) - 12D0)/24D0 
     *         + dble(it(5))
      return
      end 
************************************************************************
