	PROGRAM PRP1
C
C
C**********************************************************************P
C  A. WOOTTEN JULY 1980
C	PRP PERFORMS PRECESSION AND RADIAL VELOCITY CALCULATIONS
C	AND PLACES THE RESULTS IN THE APPROPRIATE LINE OF FILE
C	RUN.DAT
C	SUBROUTINE SRCPAN PERFORMS PRECESSION AND ALLIED CALCULATIONS
C	ACCURACY IS .015 ARCSEC. ENTER A PROPER MOTION AS PER 
C      100 YEARS.
C	SUBROUTINE RADVEL CALCULATES RADIAL VELOCITIES TOWARD THE 
C	POSITION OF THE DESIRED SOURCE USING P. STUMPFF'S SUBROUTINE
C	BARVEL TO CALCULATE THE BARYCENTRIC VELOCITY COMPONENTS OF
C	THE EARTH.  ITS REPORTED ACCURACY IS 42 CM/S.
C
C**********************************************************************
	DIMENSION S(3),NAME(8),NTYPE(12),ANTLOC(3),IT(12)
	DIMENSION RR(18),IOUT(12)
	 REAL*8 DJDAY,VTOT,VE,V,CONST,TP,PMRA,PMDEC,TIMINT,RR
	REAL*8 PMA,PMD,S,DJAN0,DU,TU,UTDA,SMD,T,START,C1,GST,WLONG
	INTEGER*2 IOUT
	REAL*8	XLST,RA1,RA,DEC,RAP,DCP,ASP,DSP,VSUN,ANTLOC,VT
	DATA SIGNT/'-'/
	DATA SIGN0/' '/
	TP=3.14159265358979D0*2.D0
	C=2.997935E+05
	CONST=TP/360.D0
	CALL ASSIGN(3,'RUN.DAT')
	READ (3,1) (NTYPE(I),I=1,12)
	READ (3,2) NAME
	READ (3,*) (IOUT(K),K=1,9)
        L=1
        DO 25 I=1,9,3
        ANTLOC(L)=(100000000.*IOUT(I)+10000.*IOUT(1+I)+IOUT(2+I)*1.)
        L=L+1
25      CONTINUE
C ANTLOC CONTAINS (E,N,U) COORDINATES IN MICROMETERS
	READ (3,*) IT(1)
	READ (3,*) IDAY,IYR
	READ (3,*) IUTH,IUTM,UTS,TIMINT,ICOHIT
	CALL JULDAT(IYR,IDAY,IUTH,IUTM,UTS,DJDAY,DJAN0)
C NOW CALCULATE THE MEAN SIDEREAL TIME AS IN DOP
C DJDAY IS EPHEMERIS JULIAN DATE,DJAN0 IS GMT FROM JAN 0.0 TO NOW (DAYS)
	DU=DJAN0-2415020.D0
	TU=DU/36525.D0
	UTDA=DJDAY-DJAN0
	SMD=DU+UTDA
	T=SMD/36525.D0
	START=45.836D0+129.1794D0+8640184.542D0*(TU-.7D0)+.0929D0*TU**2
	START=(6.D0+38.D0/60.D0+START/3600.D0)/24.D0
	C1=.997269566414D0
	GST=START+UTDA/C1
C C1 IS FOR OVRO, C2 FOR MMT, C3 FOR IRTF
C1	WLONG=.3285618827D0
        WLONG=.3080122686D0
C3      WLONG=.4318657D0  
  	XLST=GST-WLONG
	XLST=XLST-IDINT(XLST)
        TYPE 4, IUTH,IUTM,UTS,UTDA,'UT'
        TYPE 36, XLST
	NTYP1=1HC 
	NTYP2=1HP 
	NTYP3=1HN
C XLST IS LOCAL MEAN SIDEREAL TIME
	IF(NTYPE(1).EQ.NTYP1) GO TO 15
10	IF(NTYPE(1).EQ.NTYP3) GO TO 20 
	IF(NTYPE(1).NE.NTYP2) GO TO 202
	CALL PLNPOS(NAME,IYR,DJDAY,XLST,UTDA,TIMINT,RAP,DCP,VT,S,RR,IFL)
	IF (IFL.GT.0) GO TO 45
	GO TO 21
45	TYPE 46,IFL
46	FORMAT(X,'OHMYGOD, ERROR # ',X,I4,X, 'HAS BEEN COMMITTED!!')
	CALL EXIT(1)
202	TYPE 203
203	FORMAT(X,'UNKNOWN SOURCE TYPE SPECIFIED')
	CALL EXIT(1)
20	CONTINUE
C	THIS ASSIGN AND READ IS NECESSARY TO OPEN A FILE SO THAT CLOSE ON
C	EXIT DOESN'T BOMB
	CALL ASSIGN(4,'MMCAT.DAT')
	READ (4,5) (NTYPE(1))
22	FORMAT(12A1)
C READ IN SOURCE COORDINATES
	READ (3,*)  RARE,DECRE
        IRAH=IFIX(RARE/1.E4)
        IRAM=IFIX(RARE/1.E2)-IRAH*100
        RAS=RARE-IRAH*10000.-IRAM*100.
        IF (DECRE.LT.0.) SIGN=SIGNT
        DECRE=ABS(DECRE)
        IDD=IFIX(DECRE/1.E4)
        IDM=IFIX(DECRE/1.E2)-IDD*100
        DS=DECRE-IDD*10000.-IDM*100.
        PMRA=0.D0
        PMDEC=0.D0
	GO TO 19
15      CALL ASSIGN(4,'MMCAT.DAT')
16	READ (4,5) (NTYPE(I),I=4,9),IRAH,IRAM,RAS,SIGN,IDD,IDM,DS,PMRA,PMDEC
	DO 18 I=1,6
	K=3+I
C	IF SOURCE NAME IS EN...... (END) NAME IS READ FROM RUN.DAT
	IF ((NTYPE(4).EQ.1HE).And. (NTYPE(5).EQ.1HN)) GO TO 20
	IF(NTYPE(K).NE.NAME(I)) GO TO 16
18	CONTINUE
	READ (3,*) DUMMY
19      ISIGN=1
	IF(SIGN.EQ.SIGNT) ISIGN=-1
	RA1=(IRAH+IRAM/60.D0+RAS/3600.D0)/24.D0
C	TYPE 6, NAME,IRAH,IRAM,RAS,SIGN,IDD,IDM,DS,PMRA,PMDEC
	RA=TP*RA1
	DEC=ISIGN*TP*(IDD+IDM/60.D0+DS/3600.D0)/360.D0
         PMRA=PMRA*15.D0*CONST/360000.D0
        PMDEC=PMDEC*CONST/360000.D0
	CALL SRCPAN(DJDAY,RA,DEC,RAP,DCP,PMRA,PMDEC)
	CALL SIDMAT(RR,TP,C1)
	GO TO 23
C	THE  OLD POSITION, SOURCE VECTOR AND ROTATION MATRICES ARE READ
21	READ (3,*) DUMMY
23	READ (3,*) DUMMY
	DO 37 I=1,7
37	READ (3,*) IT(1)
  	RAP=RAP+PMA
	DCP=DCP+PMD
	IF(RAP.GE.TP) RAP=RAP-(1.D0/TP)
	IF(RAP.LT.0.0D0) RAP=RAP+(1.D0/TP)
	ASP=RAP/CONST/15.D0
	DSP=DCP/CONST
	IH=ASP
	M=(ASP-IH)*60.
	SS=(ASP-IH-M/60.)*3600.
	ID=DSP
	IDDM=(DSP-ID)*60.
	DDS=(DSP-ID-IDDM/60.)*3600.
	ID=IABS(ID)
	 IDDM=IABS(IDDM)
	DDS=ABS(DDS)
	RARE=1.E4*IRAH+1.E2*IRAM+RAS
	RAREP=1.E4*IH+1.E2*M+SS
	DECRE=ISIGN*(1.E4*IDD+1.E2*IDM+DS)
	DECREP=1.E4*ID+1.E2*IDDM+DDS
        SIGN2=SIGN0
        SIGN=SIGN0
        IF (DEC.LT.0.) SIGN=SIGNT
        IF (DSP.LT.0.) SIGN2=SIGNT
	IF(SIGN2.EQ.SIGNT) DECREP=-1.*DECREP
	CALL RADVEL(DJDAY,XLST,RAP,DCP,VSUN,VTOT,VE,S)
	READ (3,*) (IOUT(K),K=1,6)
        FREQIF=IOUT(1)*10000000.+IOUT(2)*10000.+IOUT(3)*1.
        FREQ=IOUT(4)*100000000.+IOUT(5)*10000.+IOUT(6)*1.
        FREQIF=FREQIF/1.E3
        FREQ=FREQ/1.E3
	READ (3,*) IHARM
C	DETERMINE IF THE VELOCITY IS LSR ,HEL(IOCENTRIC) OR PLA(NETARY)
	READ (3,*) V
C        TYPE 7, V,NTYPE(10)
	IF (NTYPE(10).EQ.1HH) GO TO 30
	IF (NTYPE(10).EQ.1HP) GO TO 31
	IF (NTYPE(10).EQ.1HL) GO TO 32
	TYPE 102,(NTYPE(I),I=10,12)
102	FORMAT(X,'UNKNOWN VELOCITY TYPE',X,3A1)
30	VTOT=V+VTOT-VSUN
	GO TO 33
31      VTOT=VT
	GO TO 33
32	VTOT=VTOT+V
C	CALCULATE THE RELATIVISTIC DOPPLER FACTOR DLAMBDA/LAMBDA
33      DOPFAC=((1.+(VTOT/C))/(SQRT(1.-(VTOT/C)**2)))-1.
  	SYNTH=(FREQ*(1.-(DOPFAC))-FREQIF)/IHARM
	rewind 3
	WRITE (3,1) (NTYPE(I),I=1,8 ),SIGN0,(NTYPE(K),K=10,12) 
	WRITE (3,2) NAME
	CALL PRUIP(ANTLOC,0,IOUT)
	WRITE (3,13) (IOUT(I),I=2,4),(IOUT(J),J=6,8),(IOUT(K),K=10,12) 
C NOW TRANSFORM ANTLOC (ENU COORDINATES) INTO EQUATORIAL (XYZ) COORDINAT
C T AND TU ARE TEMPORARY STORAGE, LATITUDE IS 37.2340556 NORTH DEGREES
	TU=(37.23405556D0*TP)/360.D0
	T=-ANTLOC(2)*DSIN(TU) + ANTLOC(3)*DCOS(TU)
	ANTLOC(3)=T+2.D0*ANTLOC(2)*DSIN(TU)
	ANTLOC(2)=ANTLOC(1)
	ANTLOC(1)=T
	CALL PRUIP(ANTLOC,0,IOUT)
	WRITE (3,13) (IOUT(I),I=2,4),(IOUT(J),J=6,8),(IOUT(K),K=10,12) 
	WRITE (3,*) IDAY,IYR
        DUMMY=SNGL(TIMINT)
	WRITE (3,*) IUTH,IUTM,UTS,DUMMY,ICOHIT
	WRITE (3,*) RARE,DECRE
	WRITE (3,*) RAREP,DECREP
C	PRUIP PRINTS NOS IN A FORMAT CRAVED BY UIP
	CALL PRUIP(S,1,IOUT)
	WRITE (3,13) (IOUT(I),I=2,4),(IOUT(J),J=6,8),(IOUT(K),K=10,12) 
	S(1)=RR(1)
	S(2)=RR(2)
	S(3)=RR(3)
	CALL PRUIP(S,2,IOUT)
	WRITE (3,13) IOUT
        S(1)=RR(4)
        S(2)=RR(5)
        S(3)=RR(6)
        CALL PRUIP(S,2,IOUT)
        WRITE (3,13) IOUT
        S(1)=RR(7)
        S(2)=RR(8)
        S(3)=RR(9)
        CALL PRUIP(S,2,IOUT)
        WRITE (3,13) IOUT
        S(1)=RR(10)
        S(2)=RR(11)
        S(3)=RR(12)
        CALL PRUIP(S,2,IOUT)
        WRITE (3,13) IOUT
        S(1)=RR(13)
        S(2)=RR(14)
        S(3)=RR(15)
        CALL PRUIP(S,2,IOUT)
        WRITE (3,13) IOUT
        S(1)=RR(16)
        S(2)=RR(17)
        S(3)= RR(18)
        CALL PRUIP(S,2,IOUT)
        WRITE (3,13) IOUT
	S(1)=FREQIF*1000.
	S(2)=FREQ*1000.
	S(3)=0.
	CALL PRUIP(S,0,IOUT)
	WRITE (3,*) (IOUT(I),I=2,4),(IOUT(J),J=6,8) 
	WRITE (3,*) IHARM
        DUMMY=SNGL(V)
        DUMMY2=SNGL(VTOT)
	WRITE (3,*) DUMMY,DUMMY2

	S(1)=SYNTH*1.D6
	S(2)=0.
	CALL PRUIP(S,0,IOUT)
	WRITE (3,*) (IOUT(I),I=2,4)
	WRITE (3,11)
C	INDICATE SUCCESSFUL EXIT
	CALL EXIT(0)
1	FORMAT(12A1,X,A2)
2	FORMAT(8A1,X,A2)
3	FORMAT(X,2(I8,X),X,A2)
4	FORMAT(X,I2,X,I2,X,F6.3,X,F9.3,X,A2)
5       FORMAT(2X,6A1,3X,2I2,F6.3,1X,A1,2I2,F5.2,F7.3,F6.2,X,A2)
6       FORMAT(2X,8A1,X,2I2,F6.3,1X,A1,2I2,F5.2,F7.3,F6.2)
7	FORMAT(X,F8.3,X,A2) 
8	FORMAT(X,F10.3,X,A2)
9	FORMAT(X)
11	FORMAT(XX)
12	FORMAT(X,F10.5,X,3A1,X,A2)
13      FORMAT(X,12(I5,X))
14	FORMAT(X,F10.0,X,F10.0,X,F10.0,X,A2)
35	FORMAT(I5,X,A2)
36      FORMAT(X,E16.8)
	END



	SUBROUTINE SIDMAT(RR,TP,C1)
C	THIS SUBROUTINE CALCULATES THE SIDEREAL ROTATION MATRICES
C	FOR .1 SEC AND 30 SEC
	DIMENSION RR(18)
	REAL*8 RR,T,TH1,C1,TP
C	R IS THE .1 SECOND ROTATION MATRIX
	T=.1D0
13	TH1= (TP*T*C1)/(3600.D0*24.D0)
	RR(1)=DCOS(TH1)
	RR(2)=DSIN(TH1)
	RR(3)=0.D0
	RR(6)=RR(2)
	RR(5)=RR(1)
	RR(4)=-RR(6)
	RR(6)=0.D0
	RR(7)  =0.D0
	RR(8)=0.D0
	RR(9)=1.D0
C	R1 IS THE 30 SECOND ROTATION MATRIX
	TH1= TP*30.D0*C1/(24.D0*3600.D0)
	RR(10)=DCOS(TH1)
	RR(11)=DSIN(TH1)
	RR(12)=0.D0
	RR(15)=RR(11)
	RR(14)=RR(10)
	RR(13)=-RR(15)
	RR(15)=0.D0
	RR(16)  =0.D0
	RR(17)=0.D0
	RR(18)=1.D0

	RETURN
	END
	SUBROUTINE PRUIP(X,I,IOUT)
	DIMENSION X(3),IOUT(12)
	INTEGER*2 IOUT


	REAL*8 X,BASE,A,B,C,D
C IF I IS 0 THE BASE IS 1. (MICRON)
	IF (I.EQ.0) BASE=1.D0
C	IF I IS 1 THE BASE IS 2 EXP 31 = 1
	IF(I.EQ.1) BASE=2147483648.D0
C	IF I IS 2 THE BASE IS 2 EXP 47 = 1
	IF(I.EQ.2) BASE=140737488355328.D0
C	DO THE NORMALIZATION
	X(1)=X(1)*BASE
	X(2)=X(2)*BASE
	X(3)=X(3)*BASE
C IF X EQUALS 1, AN OVERFLOW WILL OCCUR UNLESS WE PREVENT IT THUS:
	DO 15 II=1,3
	IF (X(II).EQ.BASE) X(II)=X(II)-1.D0
15	CONTINUE
C	MAKE EACH VECTOR COMPONENT 4 FOUR-DIGIT INTEGERS
	K=1
	DO 10 II=1,3
	A=AINT(X(II)/1.D12)
	B= AINT((X(II)-A*1.D12)/1.D8)
	D=(X(II)-(A*1.D12))
        D=(D -(B*1.D8))
        C =AINT(D/1.D4)
	D=(D-(C*1.D4))
	IOUT(K)=IFIX(SNGL(A))
	IOUT(K+1)=IFIX(SNGL(B))
	IOUT(K+2)=IFIX(SNGL(C))
	IOUT(K+3)=IFIX(SNGL(D))
	K=K+4
10	CONTINUE
	IF(I.LT.2) GO TO 11
	RETURN
1	FORMAT(X,3(X,I3,X,A1,X,3(I4,X,A1,X)),X,A2)
11	CONTINUE

2	FORMAT(X,3(6X,I3,X,A1,X,2(I4,X,A1,X)),X,A2)
	RETURN
	END

      SUBROUTINE SRCPAN(DJ,AS,DS,RA,DC,PMRA,PMDEC)
        COMMON/DATES/DJ1950
      REAL*8 DS,AS,DC,RA,VECA(3),VECB(3)
      REAL*8 NORM,DELNG,CDEC,DJ,DELOBL
      REAL*8 NUTD(4)
      VECA(1)=DCOS(AS)*DCOS(DS)
      VECA(2)=DSIN(AS)*DCOS(DS)
      VECA(3)=DSIN(DS)
        YR=(DJ-DJ1950)/365.25D0
      CALL PRECES(DJ,VECA,VECB)
      CALL ERTVEL(DJ,VX,VY)
      DELOBL=0.D0
      CALL MECMEQ(DJ,VECB,VECA,+1,DELOBL,MOBL)
      VECA(1)=VECA(1)+VX
      VECA(2)=VECA(2)+VY
      NORM=0.D0
      DO 356 I=1,3
      NORM=NORM+VECA(I)*VECA(I)
356   CONTINUE
      NORM=DSQRT(NORM)
      DO 357 I=1,3
      VECA(I)=VECA(I)/NORM
357   CONTINUE
      CALL NUTATE(DJ,NUTD)
      DELNG=NUTD(1)
      DELOBL=NUTD(2)
      VECB(1)=DCOS(DELNG)*VECA(1)-DSIN(DELNG)*VECA(2)
      VECB(2)=DSIN(DELNG)*VECA(1)+DCOS(DELNG)*VECA(2)
      VECB(3)=VECA(3)
      CALL MECMEQ(DJ,VECB,VECA,-1,DELOBL,MOBL)
      DO 1001 I=1,3
      VECB(I)=VECA(I)
1001   CONTINUE
      RA=DATAN2(VECB(2),VECB(1))
        RA=RA + PMRA*YR
      IF (RA.LT.0.0D0)  RA = RA + 6.283185307179586D0
      CDEC=DSQRT(VECB(1)*VECB(1)+VECB(2)*VECB(2))
      DC=DATAN2(VECB(3),CDEC)
        DC=DC + PMDEC*YR
      RETURN
      END
C=======================================================================
      SUBROUTINE PRECES(DJ,VECIN,VECOUT)
C PRECESS FROM MEAN EQ. COORD. OF 1950.0 TO MEAN EQ. OF DATE
C  T=TROPICAL CENTURIES SINCE 1950.0
      DOUBLE PRECISION DJ50PT,DJ,DPERY
      DOUBLE PRECISION P(3,3),VECIN(3),VECOUT(3)
      DOUBLE PRECISION E,Z,TH,T
      DOUBLE PRECISION T2,T3,DEGRAD
      DOUBLE PRECISION CE,SE,CTH,STH,CZ,SZ
      DOUBLE PRECISION ASECRD,PI
	DATA ITIMES /0/
	ITIMES=ITIMES+1
	IF (ITIMES.GT.1) GO TO 12
C	PREVENT AN ATTEMPT TO OPEN PRECES.DAT MORE THAN ONCE.
	CALL ASSIGN(2,'PRECES.DAT')
12	CONTINUE
	READ (2,100,ERR=101) DJO,P
	IF ((DJ-DJO) .LT. 1.) GO TO 5
100	FORMAT(X,D16.9,/,3(1X,3(D16.9,1X)/))
101	CONTINUE
	REWIND	2
      PI=3.14159265358979D0
      DEGRAD=180.D0/PI
      ASECRD=3600.D0*DEGRAD
	DJ50PT=2433282.314D0
      DPERY=365.24219879D0
      T=(DJ-DJ50PT)/(DPERY*100.D0)
      T2=T*T
      T3=T*T*T
	E=2304.948*T+0.3022D0*T2+0.0183D0*T3
	Z=E+0.791D0*T2
      E=E/ASECRD
      Z=Z/ASECRD
	TH=2004.2555D0*T-0.4268D0*T2-0.0418*T3
      TH=TH/ASECRD
      CE=DCOS(E)
      SE=DSIN(E)
      CTH=DCOS(TH)
      STH=DSIN(TH)
      CZ=DCOS(Z)
      SZ=DSIN(Z)
      P(1,1)=CE*CTH*CZ-SE*SZ
      P(1,2)=-SE*CTH*CZ-CE*SZ
      P(1,3)=-STH*CZ
      P(2,1)=CE*CTH*SZ+SE*CZ
      P(2,2)=-SE*CTH*SZ+CE*CZ
      P(2,3)=-STH*SZ
      P(3,1)=CE*STH
      P(3,2)=-SE*STH
      P(3,3)=CTH
	DJO=DJ

	WRITE (2,100) DJO,P
5	CONTINUE
      DO 10 I=1,3
      VECOUT(I)=0.D0
      DO 11 J=1,3
      VECOUT(I)=VECOUT(I)+P(I,J)*VECIN(J)
11    CONTINUE
10    CONTINUE
	REWIND 2
      RETURN
      END
C=======================================================================
      SUBROUTINE ERTVEL(DJ,VX,VY)
C   EARTH VELOCITY INMEAN ECLIPTIC COORD. OF DATE
C   DJ JULIAN DAYS E.T.     VX,VY IN UNITS OF C
C  NEGLECTS 'ELLIPTIC ABBERATION' COMPONENT OF VELOCITY
C  SUNLNG=MEAN LONG OF SUN IN MEAN ECLIPTIC COORD OF DATE
C    DELLNG=CORRECTION TO LONGITUDE FOR ELLIPTICITY TO FIRST ORDER IN
C ECCENTRICITY;ACCURACY APPROX. .05 DEG.
C APPROX. .02 ASEC ABBER. CORRECTION ERROR
C ERTLNG=EARTH LONG MEAN ECLIPTIC COORD OF DATE
C PERLNG=MEAN LONG OF PERIGEE IN MEAN ECLIPTIC COORD OF DATE
      DOUBLE PRECISION DJ,SUNLNG,PI,ERTLNG
      DOUBLE PRECISION DJ1900,DT,D
      DOUBLE PRECISION PERLNG
      PI=3.14159265358979D0
      DJ1900=2415020.D0
      D=DJ-DJ1900
      DT=D/10000.D0
      ERTSPD=29.78
      C=299792.5
      SUNLNG=279.696678D0+0.9856473354D0*D+0.00002267D0*DT*DT
      SUNLNG=DMOD(SUNLNG,360.D0)
      SUNLNG=SUNLNG*PI/180.D0
      ERTLNG=SUNLNG+PI
      ECC=.01672
      PERLNG=281.220833D0+0.0000470684D0*D+0.0000339D0*DT*DT+
     + 0.00000007D0*DT*DT*DT
      PERLNG=PERLNG*PI/180.D0+PI
      DELLNG=ECC*DSIN(ERTLNG-PERLNG)*2
      ERTLNG=ERTLNG+DELLNG
      VX=-ERTSPD*DSIN(ERTLNG)/C
      VY=ERTSPD*DCOS(ERTLNG)/C
      RETURN
      END
C=======================================================================
      SUBROUTINE MECMEQ(DJ,VECIN,VECOUT,ISGN,DELOBL,MOBL)
C  TRANSFORMS BETWEEN MEAN EQ. &MEAN ECLIPTIC OF DATE (DELOBL=0.0)
C   ISGN=+1 IF MEQ TO MECL; ISGN=-1 IF MECL TO MEQ
C TO TRANSFORM FROM MEAN ECLIP. TO TRUE  EQ. OF DATE SET DELOBL=NUTATION
C  IN OBLIQ. & ISGN=-1
C   T=TROPICAL YRS SINCE 1950.0
      DOUBLE PRECISION DJ50PT,DJ,DPERY
      DOUBLE  PRECISION DELOBL,OBL
      DOUBLE PRECISION VECIN(3),VECOUT(3),R(3,3)
      REAL*8 MOBL,DMOBL,T
      DOUBLE PRECISION T2,T3,DEGRAD
      DOUBLE PRECISION ASECRD,PI
      DPERY=365.24219879D0
      DJ50PT=2433282.423D0
      T=(DJ-DJ50PT)/(DPERY*100.D0)
      T2=T*T
      T3=T*T*T
      PI=3.14159265358979D0
      DEGRAD=180.D0/PI
      ASECRD=3600.D0*DEGRAD
      DMOBL=-46.850D0*T-0.0034D0*T2+0.0018D0*T3
      DMOBL=DMOBL/ASECRD
      MOBL=(23.D0+26.D0/60.D0+44.84D0/3600.D0)/DEGRAD+DMOBL
      VECOUT(1)=VECIN(1)
      OBL=MOBL+DELOBL
      R(2,2)=DCOS(OBL)
      R(3,3)=R(2,2)
      R(2,3)=ISGN*DSIN(OBL)
      R(3,2)=-R(2,3)
      DO 10 I=2,3
      VECOUT(I)=0.D0
      DO 11 J=2,3
      VECOUT(I)=VECOUT(I)+R(I,J)*VECIN(J)
11    CONTINUE
10    CONTINUE
      RETURN
      END
C=======================================================================
      SUBROUTINE JULDAT( YRA,  DAY, HR, MIN, SEC, DJ, DJAN0)
C     PROGRAM TO CHANGE YEAR/ MONTH/ DAY INTO JULIAN DATE
C     DJ1950 IS THE JULIAN DATE OF JAN 0 , 1950.
        COMMON/DATES/DJ1950
      INTEGER YR, MO, YRA,      DAY , HR, MIN , YRDY
      DOUBLE PRECISION DJ,       DJAN0, DJ1950, DFRAC
      DJ1950 = 2433281.5D0
	YR=YRA-1900
      YRDY =   DAY
      IYRPQ = (YR- 50)*365 + (YR - 49)/4
      DJAN0 = DJ1950 + IYRPQ
      DFRAC=HR/24.D0+MIN/1440.D0+SEC/86400.D0
      DJ = DJAN0 + YRDY + DFRAC
      RETURN
      END
      SUBROUTINE NUTATE(TJD,NUT)                                               1
C     THIS IS SUBROUTINE OF JAY LIESKE, MODIFIED TO GIVE ONLY NUT(1), AND      2
C     NUT(2).  TIME RATE OF CHANGE TERMS NOT GIVEN.                            3
C     ALSO MODIFIED TO REMEMBER LAST TIME AT WHICH NUTATION OBTAINED           4
C     AND TO UTILIZE THIS IF POSSIBLE RATHER THAN RECALCULATE.                 5
      DOUBLE PRECISION TJDSAV,NUT1SV,NUT2SV                                    6
      DOUBLE PRECISION TJD                                              NU     7
C--SUBROUTINE TO COMPUTE NUTATION ANGLES AND RATES FROM EXPRESSIONS     NU     8
C--    GIVEN BY WOOLARD TABLE 26 OF APAE VOL. 15 PART 1 (1953)          NU     9
C--    WHICH IS EQUIVALENT TO EXPLANATORY SUPPLEMENT PAGES 44-45.       NU    10
C--  NOTE   RATES FOR TIME-VARIATION OF COEFFICIENTS WERE ORIGINALLY    NU    11
C--    GIVEN BY WOOLARD  AS RATES PER TROPICAL CENTURY.   HERE ALL RATESNU    12
C--    ARE EVALUATED IN TERMS OF JULIAN CENTURIES.  DIFFERENCE IS ZILCH.NU    13
C--    CODED BY JAY LIESKE 11/69                                        NU    14
C--NUT(1) CONTAINS NUTATION IN LONGITUDE (RADIANS)                      NU    15
C--NUT(2) CONTAINS NUTATION IN OBLIQUITY(RADIANS)                       NU    16
C--NUT(3) CONTAINS NUTATION RATE IN LONGITUDE (RADIANS PER DAY)         NU    17
C--NUT(4) CONTAINS NUTATION RATE IN OBLIQUITY (RADIANS PER DAY)         NU    18
      DOUBLE PRECISION ANGL(5),ANGRAT(5),NUT(4)                         NU    19
      DOUBLE PRECISION TWOPI,TB,T,THETA,THDOT,SINT,COST                 NU    20
      DOUBLE PRECISION ARG(4,5)                                               21
C--ARGUMENT L  = L0 + LAT + LAAT**2 + LAAAT**3                          NU    22
      DATA ARG              /0.822512800925925923D0,                    NU    23
     A1325.55235863425925D0, 2.55324074074074074D-5,                    NU    24
     B3.99691358024691357D-8,0.995766203703703701D0,
     C99.9973604166666660D0,-4.16666666666666667D-7,                    NU    28
     D-9.25925925925925926D-9,3.12524691358024691D-2, 
     E1342.22784763888888D0, -8.91975308641975310D-6,                   NU    32
     F-9.25925925925925926D-10,0.974270794753086417D0,
     G1236.85309504629629D0, -3.98919753086419752D-6,                   NU    36
     H5.24691358024691355D-9,0.719953541666666667D0,
     I-5.37261668981481480D0, 5.77160493827160492D-6,                   NU    40
     J6.17283950617283948D-9  /                                         NU    41
      DIMENSION IJKLM(5,69)                                                   42
      DATA IJKLM/                                                             43
     A4*0,1,  2*0,2,-2,2,  4*0,2,  2*0,2,0,2,  0,1,3*0,                 NU    44
     B1,4*0,  0,1,2,-2,2,  2*0,2,0,1,  1,0,2,0,2,  0,-1,2,-2,2,         NU    45
     C2*0,2,-2,1,  0,2,3*0,  0,2*2,-2,2  ,                              NU    46
     D -1,0,2,0,2,  1,3*0,1,  -1,3*0,1,  -1,0,3*2,  -2,0,2,0,1,         NU    47
     E 1,0,2,0,1,  2*0,3*2,  1,0,2,-2,2,  2,0,2,0,2,  -1,0,2,0,1,       NU    48
     F 0,1,2*0,1,  -1,2*0,2,1,  1,2*0,-2,1,  0,-1,2*0,1,  -1,0,2*2,1,   NU    49
     G 0,1,2,0,2,  0,-1,2,0,2,  1,0,3*2,  2,0,2,-2,2,  3*0,2,1,         NU    50
     H -2,2*0,2,1,  0,-1,2,-2,1,  1,0,2,-2,1,  3*0,-2,1,  2*0,2*2,1,    NU    51
     I 0,-2,2,-2,1,  2,2*0,-2,1,  2,0,2,0,1,  0,1,2,-2,1,  -2,0,2,0,2 , NU    52
     J 1,2*0,-2,0,  3*0,2,0,  2,2*0,-2,0,  2,4*0, 2*0,2,2*0,            NU    53
     K 2*0,2,-2,0,  2,0,-2,2*0,  2*1,0,-2,0,  1,2*0,2,0,  1,-1,3*0,     NU    54
     L 0,1,0,-2,0,  3*0,1,0,  1,0,-2,2*0,  1,0,2,2*0,  2*1,3*0,         NU    55
     M 1,2*0,-1,0,  1,-1,2,0,2,  1,-1,0,-1,0,  -2,3*0,1,  -1,0,2,-2,1,  NU    56
     N 2,3*0,1,  2*-1,3*2,  0,-1,3*2,  1,3*0,2,  2*1,2,0,2,             NU    57
     O 3,0,2,0,2  /                                                     NU    58
      DOUBLE PRECISION CT(2,13),ST(2,13),C(30),S(30),SL(26)             NU    59
C--OBLIQUITY.  COSINE TERMS                                             NU    60
      DATA CT                      /9.21D0,91.D-5,.5522D0,-29.D-5,      NU    61
     A-904.D-4,4.D-5, 884.D-4,-5.D-5, 4*0.D0, 216.D-4,-6.D-5, 183.D-4,  NU    62
     B0.D0,113.D-4,-1.D-5, -93.D-4,3.D-5, -66.D-4,3*0.D0, 7.D-4,0.D0/   NU    63
      DATA C/ -50.D-4,-31.D-4,30.D-4,22.D-4,-24.D-4,23.D-4,14.D-4,      NU    64
     C-11.D-4,11.D-4,-10.D-4,8.D-4,-7.D-4,7.D-4,2*5.D-4,-3.D-4,2*3.D-4, NU    65
     D-2.D-4,3*3.D-4,-3.D-4,2*3.D-4,2.D-4,-2.D-4,2.D-4,-2.D-4,2.D-4/    NU    66
C--LONGITUDE.  SINE TERMS                                               NU    67
      DATA ST                      / -17.2327D0,-1737.D-5,-1.2729D0,    NU    68
     E-13.D-5,2088.D-4,2.D-5, -2037.D-4,-2.D-5, 1261.D-4,-31.D-5,       NU    69
     F675.D-4,1.D-5, -497.D-4,12.D-5, -342.D-4,-4.D-5, -261.D-4,0.D0,   NU    70
     G214.D-4,-5.D-5, 124.D-4,1.D-5, 16.D-4,-1.D-5, -15.D-4,1.D-5/      NU    71
      DATA S/ 114.D-4,58.D-4,-57.D-4,-52.D-4,45.D-4,-44.D-4,-32.D-4,    NU    72
     H26.D-4,-26.D-4,19.D-4,-15.D-4,14.D-4,-13.D-4,-10.D-4,-9.D-4,7.D-4,NU    73
     I 2*-6.D-4,6.D-4,-6.D-4,2*-5.D-4,5.D-4,2*-5.D-4,-4.D-4,4.D-4,      NU    74
     J-4.D-4,3.D-4,-3.D-4/                                              NU    75
      DATA SL/ -149.D-4,60.D-4,45.D-4,28.D-4,25.D-4,-21.D-4,10.D-4,     NU    76
     K-7.D-4,6.D-4,4.D-4,2*-4.D-4,4.D-4,3.D-4,3*-3.D-4,3*-2.D-4,2.D-4,  NU    77
     L3*-2.D-4,2.D-4,-2.D-4  /                                          NU    78
      DATA TWOPI/6.28318530717958648D0/                                 NU    79
      DATA TB/2415020.0D0/                                              NU    80
      DATA TJDSAV/-1.D37/                                                    81.
	DATA ITIMES /0/
      IF (TJD.EQ.TJDSAV)  GO TO 201                                           82
	ITIMES=ITIMES+1
	IF(ITIMES.GT.1) GO TO 200
	CALL ASSIGN(1,'NUTATE.DAT')
200	CONTINUE
	READ (1,100,ERR=101) DJO,NUT
C	ITIMES COUNTS CALLS, ONLY ALLOWING NUTATE.DAT TO BE OPENED
C	WHEN ITIMES=1
100	FORMAT(X,D16.9,/,1X,4(D16.9,1X))
	IF((TJD-DJO).LT.1.) GO TO 15
101	CONTINUE
	REWIND 1
      T=(TJD-TB)/36525.D0                                               NU    83
C     DO 1 J=1,4                                                        NU    84
C   1 NUT(J)=0.D0                                                       NU    85
      NUT(1)=0.D0                                                             86
      NUT(2)=0.D0                                                             87
      DO 2 I=1,5                                                        NU    88
      ANGL(I)=DMOD(ARG(1,I)+T*(ARG(2,I)+T*(ARG(3,I)+T*ARG(4,I))),1.D0)  NU    89
     1*TWOPI                                                            NU    90
C   2 ANGRAT(I)=(ARG(2,I)+T*(2.D0*ARG(3,I)+3.D0*T*ARG(4,I)))*TWOPI/     NU    91
C    236525.D0                                                          NU    92
    2 CONTINUE                                                                93
      DO 6 J=1,69                                                       NU    94
      THETA=0.D0                                                        NU    95
C     THDOT=0.D0                                                        NU    96
      DO 3 K=1,5                                                        NU    97
      IF(IJKLM(K,J).EQ.0) GO TO 3                                       NU    98
      THETA=THETA+DBLE(FLOAT(IJKLM(K,J)))*ANGL(K)                       NU    99
C     THDOT=THDOT+DBLE(FLOAT(IJKLM(K,J)))*ANGRAT(K)                     NU   100
    3 CONTINUE                                                          NU   101
      SINT=DSIN(THETA)                                                  NU   102
      COST=DCOS(THETA)                                                  NU   103
      IF(J.GT.13) GO TO 4                                               NU   104
C--  J.LE.13                                                            NU   105
      NUT(1)=NUT(1)+(ST(1,J)+T*ST(2,J))*SINT                            NU   106
      NUT(2)=NUT(2)+(CT(1,J)+T*CT(2,J))*COST                            NU   107
C     NUT(3)=NUT(3)+(ST(1,J)+T*ST(2,J))*COST*THDOT+ST(2,J)*SINT/36525.D0NU   108
C     NUT(4)=NUT(4)-(CT(1,J)+T*CT(2,J))*SINT*THDOT+CT(2,J)*COST/36525.D0NU   109
      GO TO 6                                                           NU   110
    4 CONTINUE                                                          NU   111
      IF(J.GT.43) GO TO 5                                               NU   112
C--  14 .GE. J . LE. 43                                                 NU   113
      NUT(1)=NUT(1)+S(J-13)*SINT                                        NU   114
C     NUT(3)=NUT(3)+S(J-13)*COST*THDOT                                  NU   115
      NUT(2)=NUT(2)+C(J-13)*COST                                        NU   116
C     NUT(4)=NUT(4)-C(J-13)*SINT*THDOT                                  NU   117
      GO TO 6                                                           NU   118
    5 CONTINUE                                                          NU   119
C--  44 .GE. J .LE. 69                                                  NU   120
      NUT(1)=NUT(1)+SL(J-43)*SINT                                       NU   121
C     NUT(3)=NUT(3)+SL(J-43)*COST*THDOT                                 NU   122
    6 CONTINUE                                                          NU   123
C     DO 7 I=1,4                                                        NU   124
C   7 NUT(I)=NUT(I)/206264.8062470964D0                                      125
C*******************         CONVERT FROM ARC SECONDS TO RADIANS.            126
      NUT(1)=NUT(1)/206264.8062470964D0                                      127
      NUT(2)=NUT(2)/206264.8062470964D0                                      128
      TJDSAV=TJD                                                             129
      NUT1SV=NUT(1)                                                          130
      NUT2SV=NUT(2)                                                          131

	DJO=TJD
	WRITE (1,100) DJO,NUT
	rewind 1
      RETURN                                                            NU   132
  201 NUT(1)=NUT1SV                                                          133
      NUT(2)=NUT2SV                                                          134
15	CONTINUE
	REWIND 1
      RETURN                                                                 135
      END                                                               NU   136
	SUBROUTINE RADVEL(DJDAY,XLST,ADATE,DDATE,VSUN,VTOT,VOBS,S)
			         DIMENSION DVEL(3),S(3)
          DOUBLE PRECISION DVEL,VOBS,VSUN,S,XLST
        DOUBLE PRECISION DXDT,DYDT,DZDT,VBARVL,VTOT,SALPH,SDEC,DJDAY,HA
	DOUBLE PRECISION DDATE,ADATE,CAT,ALAT,OLONG,WLONG,DLAT,RHO,VRHO
	DAUKM=1.495978D+08
	TP=6.2831853071796D0
C	PRECESSION OF THE SOLAR MOTION
	RAS=18.D0*15.D0*TP/360.D0
	DCS=30.D0*TP/360.D0
	CALL SRCPAN(DJDAY,RAS,DCS,SALPH,SDEC,0.,0.)
C	TRANSFORMATION OF SOLAR MOTION TO CARTESIAN COORDINATES
	ZS=20.D0*DSIN(SDEC)
        XS=20.D0*DCOS(SALPH)*DCOS(SDEC)
        YS=20.D0*DSIN(SALPH)*DCOS(SDEC)
C	CALCULATE HOUR ANGLE
	HA=TP*XLST-ADATE
5       FORMAT(X,3E16.8)
C	PROJECT SOLAR MOTION ONTO STELLAR LINE OF SIGHT
	S(1)=DCOS(DDATE)*DCOS(ADATE)
	S(2)=DCOS(DDATE)*DSIN(ADATE)
	S(3)=DSIN(DDATE)
	VSUN=-XS*S(1)-YS*S(2)-ZS*S(3)
C1RADIUS VECTOR TO #1 10.4 M TELESCOPE
C2RADIUS VECTOR TO MMT TELESCOPE
C1	(LONG=118 16' 56".2, LAT 37 14' 02".6, EL 1222M)
C2      (LONG+110 53' 03".9, LAT 31 41' 19".6, EL 2608M)
C3      (LONG=155 28'.3, LAT 19 49'.6, EL 4000M????)
C	REDUCING GEODETIC TO GEOCENTRIC LATITUDE AND PROJECTING
C
C	ROTATIONAL VELOCITY ONTO THE STELLAR LINE OF SIGHT.
C1	ALAT=37.23405556D0
C2      ALAT=31.68877778D0
        ALAT=19.826667D0 
C1	OLONG=118.2822778D0
C2      OLONG=110.8844167D0
        OLONG=155.47167D0 
	CAT=ALAT*TP/360.D0
	WLONG=OLONG/360.D0
C ORIGINAL CODE HAD 1920.D0 M HERE, WHICH MAY NOT BE RIGHT 9/28/83
C1	EL=1920.D0
C2      EL=2608.D0
        EL=4000.D0
	RHO=6378160.D0*(.998327073D0+.001676438D0*DCOS(2.D0*CAT)-
     A   3.51D-6*DCOS(4.D0*CAT)+.000000008D0*DCOS(6.D0*CAT))+EL
	VRHO=TP*RHO
	C1=.997269566414D0
        VRHO=VRHO/24.0D3/3600.D0/C1
	DLAT=-(11.D0*60.D0+32.7430D0)*DSIN(2.D0*CAT)+1.1633D0*DSIN(4.D0*
     A	CAT)-.0026D0*DSIN(6.D0*CAT)
	CAT=CAT+DLAT*TP/3600.D0/360.D0
	VOBS=VRHO*DCOS(CAT)*DCOS(DDATE)*DSIN(HA)
C	NOW CALCULATE MOTION OF CENTER OF EARTH RELATIVE TO CENTER
C	OF SOLAR SYSTEM
	DEQ=0.0
	CALL BARVEL(DJDAY,DEQ,DVEL)
	DXDT=DVEL(1)*DAUKM
	DYDT=DVEL(2)*DAUKM
	DZDT=DVEL(3)*DAUKM
C	PROJECT THESE ON THE STELLAR LINE OF SIGHT
	VBARVL=DXDT*S(1)+DYDT*S(2)+DZDT*S(3)
C	NOW PUT INTO S THE SOURCE VECTOR IN K Y LO'S XYZ SYSTEM 
	S(1)=DCOS(DDATE)*DCOS(HA)
	S(2)=-DCOS(DDATE)*DSIN(HA)

	VBARVL=-VBARVL
C        SUM VELOCITIES
	VTOT=VBARVL+VOBS+VSUN 
	RETURN





	END
      SUBROUTINE BARVEL(DJE,DEQ,DVEL)                                   BARVL001


C                                                                       BARVL002
C   IBM-VERSION                                                         BARVL003
C                                                                       BARVL004
C***********************************************************************BARVL005
C   P.STUMPFF   3.10.1978                                               BARVL006
C   CALCULATION OF BARYCENTRIC VELOCITY COMPONENTS OF THE EARTH.        BARVL007
C   LARGEST DEVIATIONS FROM JPL-DE96 ARE ABOUT 42 CM/S  WITHIN          BARVL008
C   THE TIME INTERVAL 1945-2000 .                                       BARVL009
C                                                                       BARVL010
C   GIVEN ARE   DJE=JULIAN EPHEMERIS DATE                               BARVL011
C               DEQ=EPOCH OF MEAN EQUATOR AND EQUINOX OF DVEL(K)        BARVL012
C                   IF DEQ=0.0, DVEL(K) IS REFERRED TO MEAN EQUINOX     BARVL013
C                   AND MEAN EQUATOR OF DJE                             BARVL014
C                                                                       BARVL015
C   RESULT  DVEL(K)=VELOCITY COMPONENTS    K=1,2,3---DX/DT,DY/DT,DZ/DT  BARVL016
C              UNIT=A.U./SEC.                                           BARVL017
C***********************************************************************BARVL018
C                                                                       BARVL019
      IMPLICIT REAL*8 (D)                                               BARVL020
      DIMENSION DVEL(3),SN(4),PL(4),PO(4),PN(4),PE(4),PI(4),            BARVL021
     *          FORBEL(7),SORBEL(17),DPREMA(3,3)                        BARVL022
      DIMENSION DCFEL(3,8),DCEPS(3),CCSEL(3,17),DCARGS(2,15),           BARVL023
     *        CCAMPS(5,15),CCSEC(3,4),DCARGM(2,3),CCAMPM(4,3),CCPAMV(4) BARVL024
      EQUIVALENCE (SORBEL(1),E),(FORBEL(1),G),                          BARVL025
     *            (FORBEL(4),PL(1)),(SORBEL(2),PO(1)),(SORBEL(6),PN(1)),BARVL026
     *            (SORBEL(10),PE(1)),(SORBEL(14),PI(1))                 BARVL027
C                                                                       BARVL028
      DATA DC2PI/6.2831853071796D0/, CC2PI/6.283185/                    BARVL029
C                                                                       BARVL030
C   CONSTANTS OF FAST CHANGING ELEMENTS                                 BARVL031
      DATA DCFEL/                                                       BARVL032
     *  1.7400353D+00, 6.2833195099091D+02, 5.2796D-06,                 BARVL033
     *  6.2565836D+00, 6.2830194572674D+02,-2.6180D-06,                 BARVL034
     *  4.7199666D+00, 8.3997091449254D+03,-1.9780D-05,                 BARVL035
     *  1.9636505D-01, 8.4334662911720D+03,-5.6044D-05,                 BARVL036
     *  4.1547339D+00, 5.2993466764997D+01, 5.8845D-06,                 BARVL037
     *  4.6524223D+00, 2.1354275911213D+01, 5.6797D-06,                 BARVL038
     *  4.2620486D+00, 7.5025342197656D+00, 5.5317D-06,                 BARVL039
     *  1.4740694D+00, 3.8377331909193D+00, 5.6093D-06/                 BARVL040
C                                                                       BARVL041
C   CONSTANTS OF SLOWLY CHANGING ELEMENTS                               BARVL042
      DATA DCEPS/  4.093198D-01,-2.271110D-04,-2.860401D-08 /           BARVL043
      DATA CCSEL/                                                       BARVL044
     *  1.675104E-02,-4.179579E-05,-1.260516E-07,                       BARVL045
     *  2.220221E-01, 2.809917E-02, 1.852532E-05,                       BARVL046
     *  1.589963E+00, 3.418075E-02, 1.430200E-05,                       BARVL047
     *  2.994089E+00, 2.590824E-02, 4.155840E-06,                       BARVL048
     *  8.155457E-01, 2.486352E-02, 6.836840E-06,                       BARVL049
     *  1.735614E+00, 1.763719E-02, 6.370440E-06,                       BARVL050
     *  1.968564E+00, 1.524020E-02,-2.517152E-06,                       BARVL051
     *  1.282417E+00, 8.703393E-03, 2.289292E-05,                       BARVL052
     *  2.280820E+00, 1.918010E-02, 4.484520E-06,                       BARVL053
     *  4.833473E-02, 1.641773E-04,-4.654200E-07,                       BARVL054
     *  5.589232E-02,-3.455092E-04,-7.388560E-07,                       BARVL055
     *  4.634443E-02,-2.658234E-05, 7.757000E-08,                       BARVL056
     *  8.997041E-03, 6.329728E-06,-1.939256E-09,                       BARVL057
     *  2.284178E-02,-9.941590E-05, 6.787400E-08,                       BARVL058
     *  4.350267E-02,-6.839749E-05,-2.714956E-07,                       BARVL059
     *  1.348204E-02, 1.091504E-05, 6.903760E-07,                       BARVL060
     *  3.106570E-02,-1.665665E-04,-1.590188E-07/                       BARVL061
C                                                                       BARVL062
C   CONSTANTS OF THE ARGUMENTS OF THE SHORT-PERIOD PERTURBATIONS        BARVL063
C   OF THE EMB BY THE PLANETS                                           BARVL064
      DATA DCARGS/                                                      BARVL065
     *  5.0974222D+00,-7.8604195454652D+02,                             BARVL066
     *  3.9584962D+00,-5.7533848094674D+02,                             BARVL067
     *  1.6338070D+00,-1.1506769618935D+03,                             BARVL068
     *  2.5487111D+00,-3.9302097727326D+02,                             BARVL069
     *  4.9255514D+00,-5.8849265665348D+02,                             BARVL070
     *  1.3363463D+00,-5.5076098609303D+02,                             BARVL071
     *  1.6072053D+00,-5.2237501616674D+02,                             BARVL072
     *  1.3629480D+00,-1.1790629318198D+03,                             BARVL073
     *  5.5657014D+00,-1.0977134971135D+03,                             BARVL074
     *  5.0708205D+00,-1.5774000881978D+02,                             BARVL075
     *  3.9318944D+00, 5.2963464780000D+01,                             BARVL076
     *  4.8989497D+00, 3.9809289073258D+01,                             BARVL077
     *  1.3097446D+00, 7.7540959633708D+01,                             BARVL078
     *  3.5147141D+00, 7.9618578146517D+01,                             BARVL079
     *  3.5413158D+00,-5.4868336758022D+02/                             BARVL080
C                                                                       BARVL081
C   AMPLITUDES OF THE SHORT-PERIOD PERTURBATIONS                        BARVL082
      DATA CCAMPS/                                                      BARVL083
     * -2.279594E-5, 1.407414E-5, 8.273188E-6, 1.340565E-5,-2.490817E-7,BARVL084
     * -3.494537E-5, 2.860401E-7, 1.289448E-7, 1.627237E-5,-1.823138E-7,BARVL085
     *  6.593466E-7, 1.322572E-5, 9.258695E-6,-4.674248E-7,-3.646275E-7,BARVL086
     *  1.140767E-5,-2.049792E-5,-4.747930E-6,-2.638763E-6,-1.245408E-7,BARVL087
     *  9.516893E-6,-2.748894E-6,-1.319381E-6,-4.549908E-6,-1.864821E-7,BARVL088
     *  7.310990E-6,-1.924710E-6,-8.772849E-7,-3.334143E-6,-1.745256E-7,BARVL089
     * -2.603449E-6, 7.359472E-6, 3.168357E-6, 1.119056E-6,-1.655307E-7,BARVL090
     * -3.228859E-6, 1.308997E-7, 1.013137E-7, 2.403899E-6,-3.736225E-7,BARVL091
     *  3.442177E-7, 2.671323E-6, 1.832858E-6,-2.394688E-7,-3.478444E-7,BARVL092
     *  8.702406E-6,-8.421214E-6,-1.372341E-6,-1.455234E-6,-4.998479E-8,BARVL093
     * -1.488378E-6,-1.251789E-5, 5.226868E-7,-2.049301E-7, 1.678311E-8,BARVL094
     * -8.043059E-6,-2.991300E-6, 1.473654E-7,-3.154542E-7, 1.261480E-8,BARVL095
     *  3.699128E-6,-3.316126E-6, 2.901257E-7, 3.407826E-7, 2.457125E-8,BARVL096
     *  2.550120E-6,-1.241123E-6, 9.901116E-8, 2.210482E-7, 2.522960E-8,BARVL097
     * -6.351059E-7, 2.341650E-6, 1.061492E-6, 2.878231E-7,-1.738673E-7/BARVL098
C                                                                       BARVL099
C   CONSTANTS OF THE SECULAR PERTURBATIONS IN LONGITUDE                 BARVL100
      DATA CCSEC1/-7.757020E-08/                                        BARVL101
      DATA CCSEC/                                                       BARVL102
     *  9.124190E-06, 9.990265E-01, 2.622706E+00,                       BARVL103
     *  3.102810E-05, 4.035027E+00, 3.525565E-01,                       BARVL104
     *  1.289600E-06, 5.550147E-01, 2.076942E+00,                       BARVL105
     *  9.793240E-07, 5.508259E+00, 1.559103E+01/                       BARVL106
C                                                                       BARVL107
C   SIDEREAL RATE IN LONGITUDE, RATE IN MEAN ANOMALY                    BARVL108
      DATA DCSLD/ 1.990987D-07/,  CCSGD/ 1.990969E-07/                  BARVL109
C                                                                       BARVL110
C   SOME CONSTANTS USED IN THE CALCULATION OF THE LUNAR CONTRIBUTION    BARVL111
      DATA CCKM/3.122140E-05/, CCMLD/2.661699E-06/, CCFDI/2.399485E-07/ BARVL112
C                                                                       BARVL113
C   CONSTANTS OF THE ARGUMENTS OF THE PERTURBATIONS                     BARVL114
C   OF THE MOTION OF THE MOON                                           BARVL115
      DATA DCARGM/                                                      BARVL116
     *  5.1679830D+00, 8.3286911095275D+03,                             BARVL117
     *  5.4913150D+00,-7.2140632838100D+03,                             BARVL118
     *  5.9598530D+00, 1.5542754389685D+04/                             BARVL119
C                                                                       BARVL120
C   AMPLITUDES OF THE PERTURBATIONS OF THE MOON                         BARVL121
      DATA CCAMPM/                                                      BARVL122
     *  1.097594E-01, 2.896773E-07, 5.450474E-02, 1.438491E-07,         BARVL123
     * -2.223581E-02, 5.083103E-08, 1.002548E-02,-2.291823E-08,         BARVL124
     *  1.148966E-02, 5.658888E-08, 8.249439E-03, 4.063015E-08/         BARVL125
C                                                                       BARVL126
C   PRODUCT  (SEMIMAJOR AXIS * MASS * SIDEREAL RATE)                    BARVL127
C   FOR THE PLANETS                                                     BARVL128
      DATA CCPAMV/8.326827E-11,1.843484E-11,1.988712E-12,1.881276E-12/  BARVL129
C                                                                       BARVL130
C                                                                       BARVL131
C   EXECUTION                                                           BARVL132
C                                                                       BARVL133
C   TIME-ARGUMENTS                                                      BARVL134
      DT=(DJE-2415020.0D0)/36525.0D0                                    BARVL135
      T=DT                                                              BARVL136
      DTSQ=DT*DT                                                        BARVL137
      TSQ=DTSQ                                                          BARVL138
C                                                                       BARVL139
C   VALUES OF ALL ELEMENTS FOR THE INSTANT DJE                          BARVL140
      DO 100 I=1,8                                                      BARVL141
      DLOCAL=DMOD(DCFEL(1,I)+DT*DCFEL(2,I)+DTSQ*DCFEL(3,I), DC2PI)      BARVL142
      IF(I.EQ.1)  DML=DLOCAL                                            BARVL143
      IF(I.NE.1)  FORBEL(I-1)=DLOCAL                                    BARVL144
  100 CONTINUE                                                          BARVL145
      DEPS=DMOD(DCEPS(1)+DT*DCEPS(2)+DTSQ*DCEPS(3), DC2PI)              BARVL146
      DO 200 I=1,17                                                     BARVL147
      SORBEL(I)=AMOD(CCSEL(1,I)+T*CCSEL(2,I)+TSQ*CCSEL(3,I), CC2PI)     BARVL148
  200 CONTINUE                                                          BARVL149
C                                                                       BARVL150
C   SECULAR PERTURBATIONS IN LONGITUDE                                  BARVL151
      DO 300 I=1,4                                                      BARVL152
      A=AMOD(CCSEC(2,I)+T*CCSEC(3,I), CC2PI)                            BARVL153
      SN(I)=SIN(A)                                                      BARVL154
  300 CONTINUE                                                          BARVL155
C                                                                       BARVL156
C   PERIODIC PERTURBATIONS OF THE EMB (EARTH-MOON BARYCENTER)           BARVL157
      PERTL =(CCSEC(1,1)+T*CCSEC1)*SN(1) +CCSEC(1,2)*SN(2)              BARVL158
     *                 +CCSEC(1,3)*SN(3) +CCSEC(1,4)*SN(4)              BARVL159
      PERTLD=0.0                                                        BARVL160
      PERTR =0.0                                                        BARVL161
      PERTRD=0.0                                                        BARVL162
      DO 400 I=1,15                                                     BARVL163
      A=DMOD(DCARGS(1,I)+DT*DCARGS(2,I), DC2PI)                         BARVL164
      COSA=COS(A)                                                       BARVL165
      SINA=SIN(A)                                                       BARVL166
      PERTL =PERTL + CCAMPS(1,I)*COSA+CCAMPS(2,I)*SINA                  BARVL167
      PERTR =PERTR + CCAMPS(3,I)*COSA+CCAMPS(4,I)*SINA                  BARVL168
      IF(I.GE.11)   GO TO 400                                           BARVL169
      PERTLD=PERTLD+(CCAMPS(2,I)*COSA-CCAMPS(1,I)*SINA)*CCAMPS(5,I)     BARVL170
      PERTRD=PERTRD+(CCAMPS(4,I)*COSA-CCAMPS(3,I)*SINA)*CCAMPS(5,I)     BARVL171
  400 CONTINUE                                                          BARVL172
C                                                                       BARVL173
C   ELLIPTIC PART OF THE MOTION OF THE EMB                              BARVL174
      ESQ=E*E                                                           BARVL175
      DPARAM=1.0D0-ESQ                                                  BARVL176
      PARAM=DPARAM                                                      BARVL177
      TWOE=E+E                                                          BARVL178
      TWOG=G+G                                                          BARVL179
      PHI=TWOE*((1.0-ESQ*0.125)*SIN(G)+E*0.625*SIN(TWOG)                BARVL180
     *          +ESQ*0.5416667*SIN(G+TWOG) )                            BARVL181
      V=G+PHI                                                           BARVL182
      SINV=SIN(V)                                                       BARVL183
      COSV=COS(V)                                                       BARVL184
      DPSI=DPARAM/(1.0D0+E*COSV)                                        BARVL185
      PHID=TWOE*CCSGD*((1.0+ESQ*1.5)*COSV+E*(1.25-SINV*SINV*0.5))       BARVL186
      PSID=CCSGD*E*SINV/SQRT(PARAM)                                     BARVL187
C                                                                       BARVL188
C   PERTURBED HELIOCENTRIC ECLIPTICAL VELOCITY OF THE EMB, INCLUDING    BARVL189
C   A FACTOR = 1 - MASS(EARTH+MOON) = 0.99999696                        BARVL190
      D1PDRO=(1.0D0 + PERTR)*0.99999696D0                               BARVL191
      DRD=D1PDRO*(PSID+DPSI*PERTRD)                                     BARVL192
      DRLD=D1PDRO*DPSI*(DCSLD+PHID+PERTLD)                              BARVL193
      DTL=DMOD(DML+PHI+PERTL, DC2PI)                                    BARVL194
      DSINL=DSIN(DTL)                                                   BARVL195
      DCOSL=DCOS(DTL)                                                   BARVL196
      DXD= DRD*DCOSL-DRLD*DSINL                                         BARVL197
      DYD= DRD*DSINL+DRLD*DCOSL                                         BARVL198
C                                                                       BARVL199
C   INFLUENCE OF ECCENTRICITY, EVECTION AND VARIATION ON THE MOON"S     BARVL200
C   GEOCENTRIC VELOCITY                                                 BARVL201
      PERTL =0.0                                                        BARVL202
      PERTLD=0.0                                                        BARVL203
      PERTP =0.0                                                        BARVL204
      PERTPD=0.0                                                        BARVL205
      DO 500 I=1,3                                                      BARVL206
      A=DMOD(DCARGM(1,I)+DT*DCARGM(2,I), DC2PI)                         BARVL207
      SINA=SIN(A)                                                       BARVL208
      COSA=COS(A)                                                       BARVL209
      PERTL =PERTL+CCAMPM(1,I)*SINA                                     BARVL210
      PERTLD=PERTLD+CCAMPM(2,I)*COSA                                    BARVL211
      PERTP =PERTP+CCAMPM(3,I)*COSA                                     BARVL212
      PERTPD=PERTPD-CCAMPM(4,I)*SINA                                    BARVL213
  500 CONTINUE                                                          BARVL214
C                                                                       BARVL215
C   HELIOCENTRIC ECLIPTICAL VELOCITY OF THE EARTH                       BARVL216
      TL=FORBEL(2)+PERTL                                                BARVL217
      SINL=SIN(TL)                                                      BARVL218
      COSL=COS(TL)                                                      BARVL219
      SIGMA=CCKM/(1.0+PERTP)                                            BARVL220
      A=SIGMA*(CCMLD+PERTLD)                                            BARVL221
      B=SIGMA*PERTPD                                                    BARVL222
      DXD=DXD+A*SINL+B*COSL                                             BARVL223
      DYD=DYD-A*COSL+B*SINL                                             BARVL224
      DZD=-SIGMA*CCFDI*COS(FORBEL(3))                                   BARVL225
C                                                                       BARVL226
C   INFLUENCE OF THE PLANETS AND CALCULATION OF THE BARYCENTRIC         BARVL227
C   ECLIPTICAL VELOCITY OF THE EARTH                                    BARVL228
      DO 600 I=1,4                                                      BARVL229
      PLON=PL(I)                                                        BARVL230
      POMG=PO(I)                                                        BARVL231
      PECC=PE(I)                                                        BARVL232
      PAMV=CCPAMV(I)                                                    BARVL233
      TL=AMOD(PLON+2.0*PECC*SIN(PLON-POMG), CC2PI)                      BARVL234
      DXD=DXD+PAMV*(SIN(TL)+PECC*SIN(POMG))                             BARVL235
      DYD=DYD-PAMV*(COS(TL)+PECC*COS(POMG))                             BARVL236
      DZD=DZD-PAMV*PI(I)*COS(PLON-PN(I))                                BARVL237
  600 CONTINUE                                                          BARVL238
C                                                                       BARVL239
C   TRANSFORMATION TO MEAN EQUATOR OF DATE                              BARVL240
      DCOSEP=DCOS(DEPS)                                                 BARVL241
      DSINEP=DSIN(DEPS)                                                 BARVL242
      DYAD=DCOSEP*DYD-DSINEP*DZD                                        BARVL243
      DZAD=DSINEP*DYD+DCOSEP*DZD                                        BARVL244
C                                                                       BARVL245
      DVEL(1)=DXD                                                       BARVL247
      DVEL(2)=DYAD                                                      BARVL248
      DVEL(3)=DZAD                                                      BARVL249
      RETURN                                                            BARVL250
C                                                                       BARVL251
C   PRECESSION FROM DATE TO EQ                                          BARVL252
      END                                                               BARVL259
C
C
       SUBROUTINE PLNPOS(NPLAN,NYR,DJDAY,SLST,SUT,UTINT,ALPHAS,
     1   DELTAS,VRAD,UVEC,RMAT,IFLAG)
C
C PLANET POSITION SUBROUTINE TO BE USED BY PRP.  (GLB---JULY,1980)
C
       REAL*8 DJDAY,SUT,DELUT,SLST,VROT,VRAD,UVEC(3),RMAT(18),UTINT,
     1   ALPH(3),DELT(3),DIST(3),ETMUT,SEPT,EPOS(3,3),POS(3,2),EFREQ,
     2   SIDT(2),A,B,C,ALPHA(2),DELTA(2),T(2),HA(2),CH,SH,CD,SD,
     3   HORP(2),DELH,DELD,CDH,SDH,CDD,SDD,CONST(4),VE1,VE2,    
     4   S(3),ALPHAS,DELTAS,VSUN,VTOT
       DIMENSION NPLAN(8)
       CONST(1) = 2.062648062470964D5
       CONST(2) = 1.7314560D3
       CONST(3) = 1.002737811906D0
       CONST(4) = 2.0D0 * 3.14159265358979D0
C ET - UT IS ESTIMATED TO BE 51.4 SECONDS ON 08/01/80:
       ETMUT = 5.14D1 / 8.64D4
       SEPT = SUT + ETMUT
C CONVERT UTINT (MINUTES) TO DELUT (FRACTIONS OF A DAY):
       DELUT = UTINT / 1.44D3
       CALL PLNFND(NPLAN,NYR,SEPT,ALPH,DELT,DIST,EFREQ,IFLAG)
       IF (IFLAG .NE. 0) GO TO 900
       DO 50 J = 1,3
       EPOS(1,J) = ALPH(J)
       EPOS(2,J) = DELT(J)
       EPOS(3,J) = DIST(J)
   50  CONTINUE
       DO 200 I = 1,3
       A = 5.0D-1 * (EPOS(I,3) - 2.0D0 * EPOS(I,2) + EPOS(I,1))
       B = 5.0D-1 * (4.0D0 * EPOS(I,2) - 3.0D0 * EPOS(I,1) - EPOS(I,3))
       C = EPOS(I,1)
       T(1) = DMOD((EFREQ * DMOD(SEPT,1.0D0)),1.0D0)
       T(2) = T(1) + EFREQ * DELUT
       DO 100 L = 1,2
       POS(I,L) = A * T(L) * T(L) + B * T(L) + C
  100  CONTINUE
  200  CONTINUE
       VRAD = CONST(2) * (POS(3,2) - POS(3,1)) / DELUT
       SIDT(2) = SLST + CONST(3) * DELUT
       SIDT(1) = SLST
       ALPHA(1) = POS(1,1)
       ALPHA(2) = POS(1,2)
       DELTA(1) = POS(2,1)
       DELTA(2) = POS(2,2)
       HORP(1) = (8.794D0 / POS(3,1)) / CONST(1)
       HORP(2) = (8.794D0 / POS(3,2)) / CONST(1)
C THE FOLLOWING CONSTANTS ARE FROM THE AENA AND PRESUMABLY REFER TO
C THE OLD CENTER PAD:
C   RHO * COS(LAT:PRIME) = 0.79733
C         TAN(LAT:PRIME) = 0.75482
       DO 400 L = 1,2
       DO 300 J = 1,4
       HA(L) = CONST(4) * SIDT(L) - ALPHA(L)
       CH = DCOS(HA(L))
       SH = DSIN(HA(L))
       CD = DCOS(DELTA(L))
       SD = DSIN(DELTA(L))
       ALPHA(L) = POS(1,L) - 7.9733D-1 * HORP(L) * SH / CD
       DELTA(L) = POS(2,L) - 7.9733D-1 * HORP(L) * (7.5482D-1 * CD
     1   - CH * SD)
  300  CONTINUE
  400  CONTINUE
       ALPHAS = ALPHA(1)
       DELTAS = DELTA(1)
       CALL RADVEL(DJDAY,SIDT(1),ALPHA(1),DELTA(1),VSUN,VTOT,VE1,UVEC)
       DJEND = DJDAY + DELUT
       CALL RADVEL(DJEND,SIDT(2),ALPHA(2),DELTA(2),VSUN,VTOT,VE2,S)
       VRAD = VRAD + (VE1 + VE2) / 2.0D0
       CH = DCOS((HA(1) + HA(2)) / 2.0D0)
       SH = DSIN((HA(1) + HA(2)) / 2.0D0)
       DELH = HA(2) - HA(1)
       DELD = DELTA(2) - DELTA(1)
       DELH = DELH / (6.0D2 * UTINT)
       DELD = DELD / (6.0D2 * UTINT)
       DO 500 I = 1,2
       N = 9 * (I - 1)
       IF (I .EQ. 1) GO TO 450
       DELH = 3.0D2 * DELH
       DELD = 3.0D2 * DELD
  450  CDH = DCOS(DELH)
       SDH = DSIN(DELH)
       CDD = DCOS(DELD)
       SDD = DSIN(DELD)
C IN THE FOLLOWING MATRIX ELEMENTS THE COS AND SIN OF HA AND
C HA + DELH ARE APPROXIMATED BY THE SIN AND COS OF THE MEAN HA FOR
C THE SCAN.
       RMAT(N + 1) = CDH * CDD
       RMAT(N + 2) = SDH * CDD
       RMAT(N + 3) = -CH * SDD
       RMAT(N + 4) = -SDH * CDD
       RMAT(N + 5) = CDH * CDD
       RMAT(N + 6) = SH * SDD
       RMAT(N + 7) = CH * SDD
       RMAT(N + 8) = -SH * SDD
       RMAT(N + 9) = CDD
  500  CONTINUE
  900  RETURN
       END
C
C
       SUBROUTINE PLNFND(NPLAN,NYR,SEPT,ALPHA,DELTA,DIS,EFREQ,
     1   IFLAG)
C
C EPHEMERIS HUNTING SUBROUTINE USED BY PLNPOS.  (GLB---JULY,1980)
C
       REAL*8 SEPT,ALPHA(3),DELTA(3),DIS(3),PT(3),T,EFREQ
       DIMENSION IRA1(5),IDEC1(6),IRA2(5),IDEC2(6),IRA3(5),IDEC3(6)
       DIMENSION IYR(3),IDAY(3),IHR(3),NPLAN(8)
 1000  FORMAT(11X,5I1,F5.3,1X,A1,5I1,F4.2,1X,F12.9,1X,I4,1X,I3,1X,I2)
       IF ((NPLAN(1) .EQ. 1HM) .AND. (NPLAN(2) .EQ. 1HE)) GO TO 21
       IF ((NPLAN(1) .EQ. 1HV) .AND. (NPLAN(2) .EQ. 1HE)) GO TO 22
       IF ((NPLAN(1) .EQ. 1HM) .AND. (NPLAN(2) .EQ. 1HA)) GO TO 23
       IF ((NPLAN(1) .EQ. 1HJ) .AND. (NPLAN(2) .EQ. 1HU)) GO TO 24
       IF ((NPLAN(1) .EQ. 1HS) .AND. (NPLAN(2) .EQ. 1HA)) GO TO 25
       IF ((NPLAN(1) .EQ. 1HU) .AND. (NPLAN(2) .EQ. 1HR)) GO TO 26
       IF ((NPLAN(1) .EQ. 1HN) .AND. (NPLAN(2) .EQ. 1HE)) GO TO 27
       IF ((NPLAN(1) .EQ. 1HS) .AND. (NPLAN(2) .EQ. 1HU)) GO TO 28
       IF ((NPLAN(1) .EQ. 1HM) .AND. (NPLAN(2) .EQ. 1HO)) GO TO 29
       IF ((NPLAN(1) .EQ. 1HI) .AND. (NPLAN(2) .EQ. 1HO)) GO TO 30
       IF ((NPLAN(1) .EQ. 1HE) .AND. (NPLAN(2) .EQ. 1HU)) GO TO 31
       IF ((NPLAN(1) .EQ. 1HG) .AND. (NPLAN(2) .EQ. 1HA)) GO TO 32
       IF ((NPLAN(1) .EQ. 1HC) .AND. (NPLAN(2) .EQ. 1HA)) GO TO 33
       IF ((NPLAN(1) .EQ. 1HT) .AND. (NPLAN(2) .EQ. 1HI)) GO TO 34
       IF ((NPLAN(1) .EQ. 1HP) .AND. (NPLAN(2) .EQ. 1HL)) GO TO 35
       IF ((NPLAN(1) .EQ. 1HC) .AND. (NPLAN(2) .EQ. 1HO)) GO TO 36
       IF ((NPLAN(1) .EQ. 1HA) .AND. (NPLAN(2) .EQ. 1HS)) GO TO 37
       IF ((NPLAN(1) .EQ. 1HO) .AND. (NPLAN(2) .EQ. 1HT)) GO TO 38
       GO TO 520
   21  CALL ASSIGN(4,' MERCURY.DAT')
       GO TO 40
   22  CALL ASSIGN(4,' VENUS.DAT')
       GO TO 40
   23  CALL ASSIGN(4,' MARS.DAT')
       GO TO 40
   24  CALL ASSIGN(4,' JUPITER.DAT')
       GO TO 40
   25  CALL ASSIGN(4,' SATURN.DAT')
       GO TO 40
   26  CALL ASSIGN(4,' URANUS.DAT')
       GO TO 40
   27  CALL ASSIGN(4,' NEPTUNE.DAT')
       GO TO 40
   28  CALL ASSIGN(4,' SUN.DAT')
       GO TO 40
   29  CALL ASSIGN(4,' MOON.DAT')
       GO TO 40
   30  CALL ASSIGN(4,' IO.DAT')
       GO TO 40
   31  CALL ASSIGN(4,' EUROPA.DAT')
       GO TO 40
   32  CALL ASSIGN(4,' GANYMEDE.DAT')
       GO TO 40
   33  CALL ASSIGN(4,' CALLISTO.DAT')
       GO TO 40
   34  CALL ASSIGN(4,' TITAN.DAT')
       GO TO 40
   35  CALL ASSIGN(4,' PLANET.DAT')
       GO TO 40
   36  CALL ASSIGN(4,' COMET.DAT')
       GO TO 40
   37  CALL ASSIGN(4,' ASTEROID.DAT')
       GO TO 40
   38  CALL ASSIGN(4,' OTHER.DAT')
   40  DO 200 I = 1,365
       T = SEPT
       READ(4,1000,END=500,ERR=510) IRA1,AS1,IDEC1,DS1,DIS(1),IYR(1),
     1   IDAY(1),IHR(1)
   45  IF(IYR(1) - NYR) 200,50,490
   50  PT(1) = 1.D0*(IDAY(1)) + 1.D0*(IHR(1)) / 2.4D1
       IF (PT(1) - T) 60,60,490
   60  READ(4,1000,END=500) IRA2,AS2,IDEC2,DS2,DIS(2),IYR(2),IDAY(2),
     1   IHR(2)
       PT(2) = 1.D0*(IDAY(2)) + 1.D0*(IHR(2)) / 2.4D1
       IF (PT(2) - T) 70,70,80
   70  DO 75 IP = 1,6
       IDEC1(IP) = IDEC2(IP)
       IF (IP .EQ. 6) GO TO 75
       IRA1(IP) = IRA2(IP)
   75  CONTINUE
       AS1 = AS2
       DS1 = DS2
       DIS(1) = DIS(2)
       IYR(1) = IYR(2)
       IDAY(1) = IDAY(2)
       IHR(1) = IHR(2)
       GO TO 45
   80  READ(4,1000,END=500) IRA3,AS3,IDEC3,DS3,DIS(3),IYR(3),IDAY(3),
     1   IHR(3)
       PT(3) = 1.D0*(IDAY(3)) + 1.D0*(IHR(3)) / 2.4D1
       IF (ABS(PT(3) - PT(2) - PT(2) + PT(1)) .GT. 1.0D-4) GO TO 480
       GO TO 210
  200  CONTINUE
  210  CONTINUE
       EFREQ = 1.0D0 / (PT(3) - PT(2))
       CALL ANGLA(0,IRA1,AS1,IDEC1,DS1,ALPHA(1),DELTA(1))
       CALL ANGLA(0,IRA2,AS2,IDEC2,DS2,ALPHA(2),DELTA(2))
       CALL ANGLA(0,IRA3,AS3,IDEC3,DS3,ALPHA(3),DELTA(3))
       IFLAG = 0
       GO TO 900
  480  IFLAG = 1
C ENTRIES NOT AT EQUAL INTERVALS.
       GO TO 900
  490  IFLAG = 2
C PREMATURE EPHEMERIS.
       GO TO 900
  500  IFLAG = 3
C OUT-OF-DATE EPHEMERIS.
       GO TO 900
  510  IFLAG = 4
C NO FILE WITH CORRECT NAME.
       GO TO 900
  520  IFLAG = 5
C PROGRAM DOESN'T RECOGNIZE THE OBJECT NAME.
  900  RETURN
       END
C
C
       SUBROUTINE ANGLA(IFL,IRA,AS,IDEC,DS,ALP,DELT)
C
C ANGLE CONVERSION SUBROUTINE.  (GLB---JULY,1980)
C
       REAL*8 ALP,DELT,CONST
       DIMENSION IRA(5),IDEC(6)
       DATA INEG/'-'/
       CONST = 206264.8062470964D0
       IF (IFL .EQ. 1) GO TO 200
       ALP = 1.5D1 * (3.6D4 * IRA(1) + 3.6D3 * IRA(2) + 6.0D2 * IRA(3)
     1   + 6.0D1 * IRA(4) + 1.0D1 * IRA(5) + AS) / CONST
       ISIGN = 1
       IF (IDEC(1) .EQ. INEG) ISIGN = -1
       DELT = ISIGN * (3.6D4 * IDEC(2) + 3.6D3 * IDEC(3) + 6.0D2 *
     2   IDEC(4) + 6.0D1 * IDEC(5) + 1.0D1 * IDEC(6) + DS) / CONST
       GO TO 900
  200  CONTINUE
  900  RETURN
       END

	SUBROUTINE EXIT(I)
	CLOSE(UNIT=1)
	CLOSE(UNIT=2)
	CLOSE(UNIT=3)
	CLOSE(UNIT=4)
C       CALL EXST(I)
        
	RETURN
	END
