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