diff --git a/CMakeLists.txt b/CMakeLists.txt index 7b0950f3a..f6496d893 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -268,11 +268,8 @@ set (wsjt_FSRCS lib/chkss2.f90 lib/coord.f90 lib/db.f90 - lib/dcoord.f90 lib/decode65a.f90 lib/decode65b.f90 - lib/fftw3mod.f90 - lib/jt9fano.f90 lib/decode4.f90 lib/decoder.f90 lib/decjt9.f90 @@ -280,18 +277,17 @@ set (wsjt_FSRCS lib/deg2grid.f90 lib/demod64a.f90 lib/determ.f90 - lib/dot.f90 lib/downsam9.f90 lib/encode232.f90 lib/encode4.f90 lib/entail.f90 + lib/ephem.f90 lib/extract.F90 lib/extract4.f90 - lib/geocentric.f90 - lib/getmet4.f90 lib/fano232.f90 lib/fchisq.f90 lib/fchisq65.f90 + lib/fftw3mod.f90 lib/fil3.f90 lib/fil3c.f90 lib/fil4.f90 @@ -311,6 +307,7 @@ set (wsjt_FSRCS lib/genwspr.f90 lib/geodist.f90 lib/getlags.f90 + lib/getmet4.f90 lib/graycode.f90 lib/graycode65.f90 lib/grayline.f90 @@ -323,13 +320,15 @@ set (wsjt_FSRCS lib/interleave63.f90 lib/interleave9.f90 lib/inter_wspr.f90 + lib/jplsubs.f lib/jt4.f90 lib/jt4a.f90 lib/jt65a.f90 + lib/jt9fano.f90 + lib/libration.f90 lib/lpf1.f90 lib/mixlpf.f90 - lib/moon2.f90 - lib/moondop.f90 + lib/moondopjpl.f90 lib/morse.f90 lib/move.f90 lib/options.f90 @@ -343,6 +342,7 @@ set (wsjt_FSRCS lib/savec2.f90 lib/sec_midn.f90 lib/setup65.f90 + lib/slasubs.f lib/sleep_msec.f90 lib/slope.f90 lib/smo.f90 @@ -358,8 +358,6 @@ set (wsjt_FSRCS lib/sync9.f90 lib/timer.f90 lib/timf2.f90 - lib/tm2.f90 - lib/toxyz.f90 lib/twkfreq.f90 lib/twkfreq65.f90 lib/wav11.f90 @@ -1002,6 +1000,12 @@ install (FILES #COMPONENT runtime ) +install (FILES + contrib/Ephemeris/JPLEPH + DESTINATION ${WSJT_SHARE_DESTINATION}/${WSJT_DATA_DESTINATION} + #COMPONENT runtime + ) + # # Mac installer files # diff --git a/contrib/Ephemeris/JPLEPH b/contrib/Ephemeris/JPLEPH new file mode 100644 index 000000000..8eea4aeda Binary files /dev/null and b/contrib/Ephemeris/JPLEPH differ diff --git a/lib/astro.f90 b/lib/astro.f90 index a0787998f..f80beb91e 100644 --- a/lib/astro.f90 +++ b/lib/astro.f90 @@ -42,11 +42,7 @@ subroutine astro(nyear,month,nday,uth,freq8,Mygrid, & call sun(nyear,month,nday,uth,lon,lat,RASun,DecSun,LST, & AzSun,ElSun,mjd,day) -! freq=nfreq*1.e6 -! if(nfreq.eq.2) freq=1.8e6 -! if(nfreq.eq.4) freq=3.5e6 - - call MoonDop(nyear,month,nday,uth,lon,lat,RAMoon,DecMoon, & + call MoonDopJPL(nyear,month,nday,uth,lon,lat,RAMoon,DecMoon, & LST,HA,AzMoon,ElMoon,vr,dist) ! Compute spatial polarization offset diff --git a/lib/astro0.f90 b/lib/astro0.f90 index 0075e1ffd..f9ce47bfb 100644 --- a/lib/astro0.f90 +++ b/lib/astro0.f90 @@ -1,14 +1,16 @@ subroutine astro0(nyear,month,nday,uth8,freq8,mygrid,hisgrid, & AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00, & dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,sd8,poloffset8,xnr8,dfdt,dfdt0, & - width1,width2,w501,w502,xlst8,techo8) + width1,width2,xlst8,techo8) parameter (DEGS=57.2957795130823d0) character*6 mygrid,hisgrid real*8 AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8 real*8 dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,xnr8,dfdt,dfdt0,dt - real*8 sd8,poloffset8,day8,width1,width2,w501,w502,xlst8 + real*8 sd8,poloffset8,day8,width1,width2,xlst8 real*8 uth8,techo8,freq8 + real*8 xl,b + common/librcom/xl(2),b(2) data uth8z/0.d0/ save @@ -19,18 +21,19 @@ subroutine astro0(nyear,month,nday,uth8,freq8,mygrid,hisgrid, & day,xlon2,xlat2,xlst,techo) AzMoonB8=AzMoon ElMoonB8=ElMoon + xl2=xl(1) + xl2a=xl(2) + b2=b(1) + b2a=b(2) call astro(nyear,month,nday,uth,freq8,mygrid,1,1, & AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler, & dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr, & day,xlon1,xlat1,xlst,techo) + xl1=xl(1) + xl1a=xl(2) + b1=b(1) + b1a=b(2) - day8=day - xlst8=xlst - techo8=techo - call tm2(day8,xlat1,xlon1,xl1,b1) - call tm2(day8,xlat2,xlon2,xl2,b2) - call tm2(day8+1.d0/1440.0,xlat1,xlon1,xl1a,b1a) - call tm2(day8+1.d0/1440.0,xlat2,xlon2,xl2a,b2a) fghz=1.d-9*freq8 dldt1=DEGS*(xl1a-xl1) dbdt1=DEGS*(b1a-b1) @@ -41,14 +44,6 @@ subroutine astro0(nyear,month,nday,uth8,freq8,mygrid,hisgrid, & rate2=sqrt((dldt1+dldt2)**2 + (dbdt1+dbdt2)**2) width2=0.5*6741*fghz*rate2 - fbend=0.7 - a2=0.0045*log(fghz/fbend)/log(1.05) - if(fghz.lt.fbend) a2=0.0 - f50=0.19 * (fghz/fbend)**a2 - if(f50.gt.1.0) f50=1.0 - w501=f50*width1 - w502=f50*width2 - AzSun8=AzSun ElSun8=ElSun AzMoon8=AzMoon diff --git a/lib/astrosub.f90 b/lib/astrosub.f90 index 10e717e1d..85eb86244 100644 --- a/lib/astrosub.f90 +++ b/lib/astrosub.f90 @@ -6,12 +6,16 @@ subroutine astrosub(nyear,month,nday,uth8,freq8,mygrid,hisgrid, & implicit real*8 (a-h,o-z) character*6 mygrid,hisgrid,c1*1 character*6 AzElFileName*(*),jpleph*(*) + character*80 jpleph_file_name logical*1 bTx + common/jplcom/jpleph_file_name + + jpleph_file_name=jpleph call astro0(nyear,month,nday,uth8,freq8,mygrid,hisgrid, & AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00, & dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,sd8,poloffset8,xnr8,dfdt,dfdt0, & - width1,width2,w501,w502,xlst8,techo8) + width1,width2,xlst8,techo8) imin=60*uth8 isec=3600*uth8 diff --git a/lib/dcoord.f90 b/lib/dcoord.f90 deleted file mode 100644 index 683be0a0a..000000000 --- a/lib/dcoord.f90 +++ /dev/null @@ -1,40 +0,0 @@ -SUBROUTINE DCOORD(A0,B0,AP,BP,A1,B1,A2,B2) - - implicit real*8 (a-h,o-z) -! Examples: -! 1. From ha,dec to az,el: -! call coord(pi,pio2-lat,0.,lat,ha,dec,az,el) -! 2. From az,el to ha,dec: -! call coord(pi,pio2-lat,0.,lat,az,el,ha,dec) -! 3. From ra,dec to l,b -! call coord(4.635594495,-0.504691042,3.355395488,0.478220215, -! ra,dec,l,b) -! 4. From l,b to ra,dec -! call coord(1.705981071d0,-1.050357016d0,2.146800277d0, -! 0.478220215d0,l,b,ra,dec) -! 5. From ecliptic latitude (eb) and longitude (el) to ra, dec: -! call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb) - - SB0=sin(B0) - CB0=cos(B0) - SBP=sin(BP) - CBP=cos(BP) - SB1=sin(B1) - CB1=cos(B1) - SB2=SBP*SB1 + CBP*CB1*cos(AP-A1) - CB2=SQRT(1.D0-SB2**2) - B2=atan(SB2/CB2) - SAA=sin(AP-A1)*CB1/CB2 - CAA=(SB1-SB2*SBP)/(CB2*CBP) - CBB=SB0/CBP - SBB=sin(AP-A0)*CB0 - SA2=SAA*CBB-CAA*SBB - CA2=CAA*CBB+SAA*SBB - TA2O2=0.0 !Shut up compiler warnings. -db - IF(CA2.LE.0.D0) TA2O2=(1.D0-CA2)/SA2 - IF(CA2.GT.0.D0) TA2O2=SA2/(1.D0+CA2) - A2=2.D0*atan(TA2O2) - IF(A2.LT.0.D0) A2=A2+6.2831853071795864D0 - - RETURN -END SUBROUTINE DCOORD diff --git a/lib/dot.f90 b/lib/dot.f90 deleted file mode 100644 index a68c5bd13..000000000 --- a/lib/dot.f90 +++ /dev/null @@ -1,11 +0,0 @@ -real*8 function dot(x,y) - - real*8 x(3),y(3) - - dot=0.d0 - do i=1,3 - dot=dot+x(i)*y(i) - enddo - - return -end function dot diff --git a/lib/ephem.f90 b/lib/ephem.f90 new file mode 100644 index 000000000..75efef533 --- /dev/null +++ b/lib/ephem.f90 @@ -0,0 +1,87 @@ +subroutine ephem(mjd0,dut,east_long,geodetic_lat,height,nspecial, & + RA,Dec,Az,El,techo,dop,fspread_1GHz,vr) + + implicit real*8 (a-h,o-z) + real*8 jd !Time of observationa as a Julian Date + real*8 mjd,mjd0 !Modified Julian Date + real*8 prec(3,3) !Precession matrix, J2000 to Date + real*8 rmatn(3,3) !Nutation matrix + real*8 rme2000(6) !Vector from Earth center to Moon, JD2000 + real*8 rmeDate(6) !Vector from Earth center to Moon at Date + real*8 rmeTrue(6) !Include nutation + real*8 raeTrue(6) !Vector from Earth center to Obs at Date + real*8 rmaTrue(6) !Vector from Obs to Moon at Date + logical km,bary,jplok !Set km=.true. to get km, km/s from ephemeris + common/stcomx/km,bary,pvsun(6) !Common used in JPL subroutines + common/librcom/xl(2),b(2) + + twopi=8.d0*atan(1.d0) !Define some constants + rad=360.d0/twopi + clight=2.99792458d5 + au2km=0.1495978706910000d9 + pi=0.5d0*twopi + pio2=0.5d0*pi + km=.true. + freq=1000.0d6 + + do jj=1,2 + mjd=mjd0 + if(jj.eq.1) mjd=mjd - 1.d0/1440.d0 + djutc=mjd + jd=2400000.5d0 + mjd + djtt=mjd + sla_DTT(jd)/86400.d0 + ttjd=jd + sla_DTT(jd)/86400.d0 + +! inquire(file='JPLEPH',exist=jplok) +! if(jplok) then + if(nspecial.ne.8) then + call pleph(ttjd,10,3,rme2000) !RME (J2000) from JPL ephemeris + + year=2000.d0 + (jd-2451545.d0)/365.25d0 + call sla_PREC (2000.0d0, year, prec) !Get precession matrix + rmeDate(1:3)=matmul(prec,rme2000(1:3)) !Moon geocentric xyz at Date + rmeDate(4:6)=matmul(prec,rme2000(4:6)) !Moon geocentric vel at Date + else + call sla_DMOON(djtt,rmeDate) !No JPL ephemeris, use DMOON + rmeDate=rmeDate*au2km + endif + + if(nspecial.eq.7) then + rmeTrue=rmeDate + else +!Nutation to true equinox of Date + call sla_NUT(djtt,rmatn) + call sla_DMXV(rmatn,rmeDate,rmeTrue) + call sla_DMXV(rmatn,rmeDate(4),rmeTrue(4)) + endif + +! Local Apparent Sidereal Time: + djut1=djutc + dut/86400.d0 + if(nspecial.eq.6) djut1=djutc + xlast=sla_DRANRM(sla_GMST(djut1) + sla_EQEQX(djtt) + east_long) + call sla_PVOBS(geodetic_lat,height,xlast,raeTrue) + rmaTrue=rmeTrue - raeTrue*au2km + + if(nspecial.ne.2) then +! Allow for planetary aberration + tl=499.004782D0*SQRT(rmaTrue(1)**2 + rmaTrue(2)**2 + rmaTrue(3)**2) + rmaTrue(1:3)=rmaTrue(1:3)-tl*rmaTrue(4:6)/au2km + endif + +!Topocentric RA, Dec, dist, velocity + call sla_DC62S(rmaTrue,RA,Dec,dist,RAdot,DECdot,vr) + dop=-2.d0 * freq * vr/clight !EME doppler shift + techo=2.d0*dist/clight !Echo delay time (s) + call libration(jd,RA,Dec,xl(jj),b(jj)) + enddo + + fspread_1GHz=0.0d0 + dldt=57.2957795131*(xl(2)-xl(1)) + dbdt=57.2957795131*(b(2)-b(1)) + rate=sqrt((2*dldt)**2 + (2*dbdt)**2) + fspread_1GHz=0.5*6741*rate + + call sla_DE2H(xlast-RA,Dec,geodetic_lat,Az,El) + + return +end subroutine ephem diff --git a/lib/geocentric.f90 b/lib/geocentric.f90 deleted file mode 100644 index 11ea7cb60..000000000 --- a/lib/geocentric.f90 +++ /dev/null @@ -1,17 +0,0 @@ -subroutine geocentric(alat,elev,hlt,erad) - - implicit real*8 (a-h,o-z) - -! IAU 1976 flattening f, equatorial radius a - f = 1.d0/298.257d0 - a = 6378140.d0 - c = 1.d0/sqrt(1.d0 + (-2.d0 + f)*f*sin(alat)*sin(alat)) - arcf = (a*c + elev)*cos(alat) - arsf = (a*(1.d0 - f)*(1.d0 - f)*c + elev)*sin(alat) - hlt = datan2(arsf,arcf) - erad = sqrt(arcf*arcf + arsf*arsf) - erad = 0.001d0*erad - - return -end subroutine geocentric - diff --git a/lib/jplsubs.f b/lib/jplsubs.f new file mode 100644 index 000000000..673690ec5 --- /dev/null +++ b/lib/jplsubs.f @@ -0,0 +1,899 @@ +C++++++++++++++++++++++++ +C + SUBROUTINE FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL) +C +C++++++++++++++++++++++++ +C +C Version 1.0 uses the INQUIRE statement to find out the the record length +C of the direct access file before opening it. This procedure is non-standard, +C but seems to work for VAX machines. +C +C THE SUBROUTINE ALSO SETS THE VALUES OF NRECL, NRFILE, AND NAMFIL. + +C ***************************************************************** +C ***************************************************************** +C +C THE PARAMETERS NAMFIL, NRECL, AND NRFILE ARE TO BE SET BY THE USER +C +C ***************************************************************** + +C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE + + CHARACTER*80 NAMFIL + +c NAMFIL='JPLEPH' + +C ***************************************************************** + +C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS +C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES +C (for a VAX, it is probably 1) +C + NRECL=4 + +C ***************************************************************** + +C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE + +c NRFILE=12 + +C ***************************************************************** + +C FIND THE RECORD SIZE USING THE INQUIRE STATEMENT + + +c IRECSZ=0 + + INQUIRE(FILE=NAMFIL,RECL=IRECSZ) + +C IF 'INQUIRE' DOES NOT WORK, USUALLY IRECSZ WILL BE LEFT AT 0 + + IF(IRECSZ .LE. 0) write(*,*) + . ' INQUIRE STATEMENT PROBABLY DID NOT WORK' + + KSIZE=IRECSZ/NRECL + + RETURN + + END +C++++++++++++++++++++++++ +C + SUBROUTINE FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL) +C +C++++++++++++++++++++++++ +C THIS SUBROUTINE OPENS THE FILE, 'NAMFIL', WITH A PHONY RECORD LENGTH, READS +C THE FIRST RECORD, AND USES THE INFO TO COMPUTE KSIZE, THE NUMBER OF SINGLE +C PRECISION WORDS IN A RECORD. +C +C THE SUBROUTINE ALSO SETS THE VALUES OF NRECL, NRFILE, AND NAMFIL. + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + + SAVE + + INTEGER OLDMAX + PARAMETER (OLDMAX = 400) + INTEGER NMAX + PARAMETER (NMAX = 1000) + CHARACTER*6 TTL(14,3),CNAM(NMAX) + CHARACTER*80 NAMFIL,jpleph_file_name + DIMENSION SS(3) + INTEGER IPT(3,13) + common/jplcom/jpleph_file_name + +C ***************************************************************** +C ***************************************************************** +C +C THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER +C +C ***************************************************************** + +C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS +C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES +C (for UNIX, it is probably 4) +C + NRECL=4 + +C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE + + NRFILE=12 + +C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE + +! NAMFIL='JPLEPH' + NAMFIL=jpleph_file_name + +C ***************************************************************** +C ***************************************************************** + +C ** OPEN THE DIRECT-ACCESS FILE AND GET THE POINTERS IN ORDER TO +C ** DETERMINE THE SIZE OF THE EPHEMERIS RECORD + + MRECL=NRECL*1000 + + OPEN(NRFILE, + * FILE=NAMFIL, + * ACCESS='DIRECT', + * FORM='UNFORMATTED', + * RECL=MRECL, + * STATUS='OLD') + + READ(NRFILE,REC=1)TTL,(CNAM(K),K=1,OLDMAX),SS,NCON,AU,EMRAT, + & ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3) + + CLOSE(NRFILE) + +C FIND THE NUMBER OF EPHEMERIS COEFFICIENTS FROM THE POINTERS + + KMX = 0 + KHI = 0 + + DO I = 1,13 + IF (IPT(1,I) .GT. KMX) THEN + KMX = IPT(1,I) + KHI = I + ENDIF + ENDDO + + ND = 3 + IF (KHI .EQ. 12) ND=2 + + KSIZE = 2*(IPT(1,KHI)+ND*IPT(2,KHI)*IPT(3,KHI)-1) + + RETURN + + END +C++++++++++++++++++++++++ +C + SUBROUTINE FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL) +C +C++++++++++++++++++++++++ +C +C THE SUBROUTINE SETS THE VALUES OF NRECL, KSIZE, NRFILE, AND NAMFIL. + + SAVE + CHARACTER*80 NAMFIL,jpleph_file_name + common/jplcom/jpleph_file_name + +C ***************************************************************** +C ***************************************************************** +C +C THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER + +C ***************************************************************** + +C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS +C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES + + NRECL=4 + +C ***************************************************************** + +C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE (DEFAULT: 12) + + NRFILE=12 + +C ***************************************************************** + +C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE + +! NAMFIL='JPLEPH' + NAMFIL=jpleph_file_name + +C ***************************************************************** + +C KSIZE must be set by the user according to the ephemeris to be read + +C For de200, set KSIZE to 1652 +C For de405, set KSIZE to 2036 +C For de406, set KSIZE to 1456 +C For de414, set KSIZE to 2036 +C For de418, set KSIZE to 2036 +C For de421, set KSIZE to 2036 +C For de422, set KSIZE to 2036 +C For de423, set KSIZE to 2036 +C For de424, set KSIZE to 2036 +C For de430, set KSIZE to 2036 + + KSIZE = 2036 + +C ******************************************************************* + + RETURN + + END +C++++++++++++++++++++++++++ +C + SUBROUTINE PLEPH ( ET, NTARG, NCENT, RRD ) +C +C++++++++++++++++++++++++++ +C NOTE : Over the years, different versions of PLEPH have had a fifth argument: +C sometimes, an error return statement number; sometimes, a logical denoting +C whether or not the requested date is covered by the ephemeris. We apologize +C for this inconsistency; in this present version, we use only the four necessary +C arguments and do the testing outside of the subroutine. +C +C THIS SUBROUTINE READS THE JPL PLANETARY EPHEMERIS +C AND GIVES THE POSITION AND VELOCITY OF THE POINT 'NTARG' +C WITH RESPECT TO 'NCENT'. +C +C CALLING SEQUENCE PARAMETERS: +C +C ET = D.P. JULIAN EPHEMERIS DATE AT WHICH INTERPOLATION +C IS WANTED. +C +C ** NOTE THE ENTRY DPLEPH FOR A DOUBLY-DIMENSIONED TIME ** +C THE REASON FOR THIS OPTION IS DISCUSSED IN THE +C SUBROUTINE STATE +C +C NTARG = INTEGER NUMBER OF 'TARGET' POINT. +C +C NCENT = INTEGER NUMBER OF CENTER POINT. +C +C THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS: +C +C 1 = MERCURY 8 = NEPTUNE +C 2 = VENUS 9 = PLUTO +C 3 = EARTH 10 = MOON +C 4 = MARS 11 = SUN +C 5 = JUPITER 12 = SOLAR-SYSTEM BARYCENTER +C 6 = SATURN 13 = EARTH-MOON BARYCENTER +C 7 = URANUS 14 = NUTATIONS (LONGITUDE AND OBLIQ) +C 15 = LIBRATIONS, IF ON EPH FILE +C +C (IF NUTATIONS ARE WANTED, SET NTARG = 14. FOR LIBRATIONS, +C SET NTARG = 15. SET NCENT=0.) +C +C RRD = OUTPUT 6-WORD D.P. ARRAY CONTAINING POSITION AND VELOCITY +C OF POINT 'NTARG' RELATIVE TO 'NCENT'. THE UNITS ARE AU AND +C AU/DAY. FOR LIBRATIONS THE UNITS ARE RADIANS AND RADIANS +C PER DAY. IN THE CASE OF NUTATIONS THE FIRST FOUR WORDS OF +C RRD WILL BE SET TO NUTATIONS AND RATES, HAVING UNITS OF +C RADIANS AND RADIANS/DAY. +C +C The option is available to have the units in km and km/sec. +C For this, set km=.true. in the STCOMX common block. +C + + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + INTEGER NMAX + PARAMETER (NMAX = 1000) + + DIMENSION RRD(6),ET2Z(2),ET2(2),PV(6,13) + DIMENSION PVST(6,11),PNUT(4) + DIMENSION SS(3),CVAL(NMAX),PVSUN(6),ZIPS(2) + DATA ZIPS/2*0.d0/ + + LOGICAL BSAVE,KM,BARY + LOGICAL FIRST + DATA FIRST/.TRUE./ + + INTEGER LIST(12),IPT(39),DENUM + COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT + COMMON/STCOMX/KM,BARY,PVSUN + common/jplcom/jpleph_file_name + +C INITIALIZE ET2 FOR 'STATE' AND SET UP COMPONENT COUNT +C + ET2(1)=ET + ET2(2)=0.D0 + GO TO 11 + +C ENTRY POINT 'DPLEPH' FOR DOUBLY-DIMENSIONED TIME ARGUMENT +C (SEE THE DISCUSSION IN THE SUBROUTINE STATE) + + ENTRY DPLEPH(ET2Z,NTARG,NCENT,RRD) + + ET2(1)=ET2Z(1) + ET2(2)=ET2Z(2) + + 11 IF(FIRST) CALL STATE(ZIPS,LIST,PVST,PNUT) + FIRST=.FALSE. + + 96 IF(NTARG .EQ. NCENT) RETURN + + DO I=1,12 + LIST(I)=0 + ENDDO + +C CHECK FOR NUTATION CALL + + IF(NTARG.NE.14) GO TO 97 + IF(IPT(35).GT.0) THEN + LIST(11)=2 + CALL STATE(ET2,LIST,PVST,PNUT) + DO I=1,4 + RRD(I)=PNUT(I) + ENDDO + RRD(5) = 0.d0 + RRD(6) = 0.d0 + RETURN + ELSE + DO I=1,4 + RRD(I)=0.d0 + ENDDO + WRITE(6,297) + 297 FORMAT(' ***** NO NUTATIONS ON THE EPHEMERIS FILE *****') + STOP + ENDIF + +C CHECK FOR LIBRATIONS + + 97 CONTINUE + DO I=1,6 + RRD(I)=0.d0 + ENDDO + + IF(NTARG.NE.15) GO TO 98 + IF(IPT(38).GT.0) THEN + LIST(12)=2 + CALL STATE(ET2,LIST,PVST,PNUT) + DO I=1,6 + RRD(I)=PVST(I,11) + ENDDO + RETURN + ELSE + WRITE(6,298) + 298 FORMAT(' ***** NO LIBRATIONS ON THE EPHEMERIS FILE *****') + STOP + ENDIF + +C FORCE BARYCENTRIC OUTPUT BY 'STATE' + + 98 BSAVE=BARY + BARY=.TRUE. + +C SET UP PROPER ENTRIES IN 'LIST' ARRAY FOR STATE CALL + + DO I=1,2 + K=NTARG + IF(I .EQ. 2) K=NCENT + IF(K .LE. 10) LIST(K)=2 + IF(K .EQ. 10) LIST(3)=2 + IF(K .EQ. 3) LIST(10)=2 + IF(K .EQ. 13) LIST(3)=2 + ENDDO + +C MAKE CALL TO STATE + + CALL STATE(ET2,LIST,PVST,PNUT) + + DO I=1,10 + DO J = 1,6 + PV(J,I) = PVST(J,I) + ENDDO + ENDDO + + IF(NTARG .EQ. 11 .OR. NCENT .EQ. 11) THEN + DO I=1,6 + PV(I,11)=PVSUN(I) + ENDDO + ENDIF + + IF(NTARG .EQ. 12 .OR. NCENT .EQ. 12) THEN + DO I=1,6 + PV(I,12)=0.D0 + ENDDO + ENDIF + + IF(NTARG .EQ. 13 .OR. NCENT .EQ. 13) THEN + DO I=1,6 + PV(I,13) = PVST(I,3) + ENDDO + ENDIF + + IF(NTARG*NCENT .EQ. 30 .AND. NTARG+NCENT .EQ. 13) THEN + DO I=1,6 + PV(I,3)=0.D0 + ENDDO + GO TO 99 + ENDIF + + IF(LIST(3) .EQ. 2) THEN + DO I=1,6 + PV(I,3)=PVST(I,3)-PVST(I,10)/(1.D0+EMRAT) + ENDDO + ENDIF + + IF(LIST(10) .EQ. 2) THEN + DO I=1,6 + PV(I,10) = PV(I,3)+PVST(I,10) + ENDDO + ENDIF + + 99 DO I=1,6 + RRD(I)=PV(I,NTARG)-PV(I,NCENT) + ENDDO + + BARY=BSAVE + + RETURN + END +C+++++++++++++++++++++++++++++++++ +C + SUBROUTINE INTERP(BUF,T,NCF,NCM,NA,IFL,PV) +C +C+++++++++++++++++++++++++++++++++ +C +C THIS SUBROUTINE DIFFERENTIATES AND INTERPOLATES A +C SET OF CHEBYSHEV COEFFICIENTS TO GIVE POSITION AND VELOCITY +C +C CALLING SEQUENCE PARAMETERS: +C +C INPUT: +C +C BUF 1ST LOCATION OF ARRAY OF D.P. CHEBYSHEV COEFFICIENTS OF POSITION +C +C T T(1) IS DP FRACTIONAL TIME IN INTERVAL COVERED BY +C COEFFICIENTS AT WHICH INTERPOLATION IS WANTED +C (0 .LE. T(1) .LE. 1). T(2) IS DP LENGTH OF WHOLE +C INTERVAL IN INPUT TIME UNITS. +C +C NCF # OF COEFFICIENTS PER COMPONENT +C +C NCM # OF COMPONENTS PER SET OF COEFFICIENTS +C +C NA # OF SETS OF COEFFICIENTS IN FULL ARRAY +C (I.E., # OF SUB-INTERVALS IN FULL INTERVAL) +C +C IFL INTEGER FLAG: =1 FOR POSITIONS ONLY +C =2 FOR POS AND VEL +C +C +C OUTPUT: +C +C PV INTERPOLATED QUANTITIES REQUESTED. DIMENSION +C EXPECTED IS PV(NCM,IFL), DP. +C +C + IMPLICIT DOUBLE PRECISION (A-H,O-Z) +C + SAVE +C + DOUBLE PRECISION BUF(NCF,NCM,*),T(2),PV(NCM,*),PC(18),VC(18) + +C + DATA NP/2/ + DATA NV/3/ + DATA TWOT/0.D0/ + DATA PC(1),PC(2)/1.D0,0.D0/ + DATA VC(2)/1.D0/ +C +C ENTRY POINT. GET CORRECT SUB-INTERVAL NUMBER FOR THIS SET +C OF COEFFICIENTS AND THEN GET NORMALIZED CHEBYSHEV TIME +C WITHIN THAT SUBINTERVAL. +C + DNA=DBLE(NA) + DT1=DINT(T(1)) + TEMP=DNA*T(1) + L=IDINT(TEMP-DT1)+1 + +C TC IS THE NORMALIZED CHEBYSHEV TIME (-1 .LE. TC .LE. 1) + + TC=2.D0*(DMOD(TEMP,1.D0)+DT1)-1.D0 + +C CHECK TO SEE WHETHER CHEBYSHEV TIME HAS CHANGED, +C AND COMPUTE NEW POLYNOMIAL VALUES IF IT HAS. +C (THE ELEMENT PC(2) IS THE VALUE OF T1(TC) AND HENCE +C CONTAINS THE VALUE OF TC ON THE PREVIOUS CALL.) + + IF(TC.NE.PC(2)) THEN + NP=2 + NV=3 + PC(2)=TC + TWOT=TC+TC + ENDIF +C +C BE SURE THAT AT LEAST 'NCF' POLYNOMIALS HAVE BEEN EVALUATED +C AND ARE STORED IN THE ARRAY 'PC'. +C + IF(NP.LT.NCF) THEN + DO 1 I=NP+1,NCF + PC(I)=TWOT*PC(I-1)-PC(I-2) + 1 CONTINUE + NP=NCF + ENDIF +C +C INTERPOLATE TO GET POSITION FOR EACH COMPONENT +C + DO 2 I=1,NCM + PV(I,1)=0.D0 + DO 3 J=NCF,1,-1 + PV(I,1)=PV(I,1)+PC(J)*BUF(J,I,L) + 3 CONTINUE + 2 CONTINUE + IF(IFL.LE.1) RETURN +C +C IF VELOCITY INTERPOLATION IS WANTED, BE SURE ENOUGH +C DERIVATIVE POLYNOMIALS HAVE BEEN GENERATED AND STORED. +C + VFAC=(DNA+DNA)/T(2) + VC(3)=TWOT+TWOT + IF(NV.LT.NCF) THEN + DO 4 I=NV+1,NCF + VC(I)=TWOT*VC(I-1)+PC(I-1)+PC(I-1)-VC(I-2) + 4 CONTINUE + NV=NCF + ENDIF +C +C INTERPOLATE TO GET VELOCITY FOR EACH COMPONENT +C + DO 5 I=1,NCM + PV(I,2)=0.D0 + DO 6 J=NCF,2,-1 + PV(I,2)=PV(I,2)+VC(J)*BUF(J,I,L) + 6 CONTINUE + PV(I,2)=PV(I,2)*VFAC + 5 CONTINUE +C + RETURN +C + END + +C+++++++++++++++++++++++++ +C + SUBROUTINE SPLIT(TT,FR) +C +C+++++++++++++++++++++++++ +C +C THIS SUBROUTINE BREAKS A D.P. NUMBER INTO A D.P. INTEGER +C AND A D.P. FRACTIONAL PART. +C +C CALLING SEQUENCE PARAMETERS: +C +C TT = D.P. INPUT NUMBER +C +C FR = D.P. 2-WORD OUTPUT ARRAY. +C FR(1) CONTAINS INTEGER PART +C FR(2) CONTAINS FRACTIONAL PART +C +C FOR NEGATIVE INPUT NUMBERS, FR(1) CONTAINS THE NEXT +C MORE NEGATIVE INTEGER; FR(2) CONTAINS A POSITIVE FRACTION. +C +C CALLING SEQUENCE DECLARATIONS +C + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + DIMENSION FR(2) + +C MAIN ENTRY -- GET INTEGER AND FRACTIONAL PARTS + + FR(1)=DINT(TT) + FR(2)=TT-FR(1) + + IF(TT.GE.0.D0 .OR. FR(2).EQ.0.D0) RETURN + +C MAKE ADJUSTMENTS FOR NEGATIVE INPUT NUMBER + + FR(1)=FR(1)-1.D0 + FR(2)=FR(2)+1.D0 + + RETURN + + END + + +C++++++++++++++++++++++++++++++++ +C + SUBROUTINE STATE(ET2,LIST,PV,PNUT) +C +C++++++++++++++++++++++++++++++++ +C +C THIS SUBROUTINE READS AND INTERPOLATES THE JPL PLANETARY EPHEMERIS FILE +C +C CALLING SEQUENCE PARAMETERS: +C +C INPUT: +C +C ET2 DP 2-WORD JULIAN EPHEMERIS EPOCH AT WHICH INTERPOLATION +C IS WANTED. ANY COMBINATION OF ET2(1)+ET2(2) WHICH FALLS +C WITHIN THE TIME SPAN ON THE FILE IS A PERMISSIBLE EPOCH. +C +C A. FOR EASE IN PROGRAMMING, THE USER MAY PUT THE +C ENTIRE EPOCH IN ET2(1) AND SET ET2(2)=0. +C +C B. FOR MAXIMUM INTERPOLATION ACCURACY, SET ET2(1) = +C THE MOST RECENT MIDNIGHT AT OR BEFORE INTERPOLATION +C EPOCH AND SET ET2(2) = FRACTIONAL PART OF A DAY +C ELAPSED BETWEEN ET2(1) AND EPOCH. +C +C C. AS AN ALTERNATIVE, IT MAY PROVE CONVENIENT TO SET +C ET2(1) = SOME FIXED EPOCH, SUCH AS START OF INTEGRATION, +C AND ET2(2) = ELAPSED INTERVAL BETWEEN THEN AND EPOCH. +C +C LIST 12-WORD INTEGER ARRAY SPECIFYING WHAT INTERPOLATION +C IS WANTED FOR EACH OF THE BODIES ON THE FILE. +C +C LIST(I)=0, NO INTERPOLATION FOR BODY I +C =1, POSITION ONLY +C =2, POSITION AND VELOCITY +C +C THE DESIGNATION OF THE ASTRONOMICAL BODIES BY I IS: +C +C I = 1: MERCURY +C = 2: VENUS +C = 3: EARTH-MOON BARYCENTER +C = 4: MARS +C = 5: JUPITER +C = 6: SATURN +C = 7: URANUS +C = 8: NEPTUNE +C = 9: PLUTO +C =10: GEOCENTRIC MOON +C =11: NUTATIONS IN LONGITUDE AND OBLIQUITY +C =12: LUNAR LIBRATIONS (IF ON FILE) +C +C OUTPUT: +C +C PV DP 6 X 11 ARRAY THAT WILL CONTAIN REQUESTED INTERPOLATED +C QUANTITIES (OTHER THAN NUTATION, STOERD IN PNUT). +C THE BODY SPECIFIED BY LIST(I) WILL HAVE ITS +C STATE IN THE ARRAY STARTING AT PV(1,I). +C (ON ANY GIVEN CALL, ONLY THOSE WORDS IN 'PV' WHICH ARE +C AFFECTED BY THE FIRST 10 'LIST' ENTRIES, AND BY LIST(12) +C IF LIBRATIONS ARE ON THE FILE, ARE SET. +C THE REST OF THE 'PV' ARRAYIS UNTOUCHED.) +C THE ORDER OF COMPONENTS STARTING IN PV(1,I) IS: X,Y,Z,DX,DY,DZ. +C +C ALL OUTPUT VECTORS ARE REFERENCED TO THE EARTH MEAN +C EQUATOR AND EQUINOX OF J2000 IF THE DE NUMBER IS 200 OR +C GREATER; OF B1950 IF THE DE NUMBER IS LESS THAN 200. +C +C THE MOON STATE IS ALWAYS GEOCENTRIC; THE OTHER NINE STATES +C ARE EITHER HELIOCENTRIC OR SOLAR-SYSTEM BARYCENTRIC, +C DEPENDING ON THE SETTING OF COMMON FLAGS (SEE BELOW). +C +C LUNAR LIBRATIONS, IF ON FILE, ARE PUT INTO PV(K,11) IF +C LIST(12) IS 1 OR 2. +C +C NUT DP 4-WORD ARRAY THAT WILL CONTAIN NUTATIONS AND RATES, +C DEPENDING ON THE SETTING OF LIST(11). THE ORDER OF +C QUANTITIES IN NUT IS: +C +C D PSI (NUTATION IN LONGITUDE) +C D EPSILON (NUTATION IN OBLIQUITY) +C D PSI DOT +C D EPSILON DOT +C +C * STATEMENT # FOR ERROR RETURN, IN CASE OF EPOCH OUT OF +C RANGE OR I/O ERRORS. +C +C COMMON AREA STCOMX: +C +C KM LOGICAL FLAG DEFINING PHYSICAL UNITS OF THE OUTPUT +C STATES. KM = .TRUE., KM AND KM/SEC +C = .FALSE., AU AND AU/DAY +C DEFAULT VALUE = .FALSE. (KM DETERMINES TIME UNIT +C FOR NUTATIONS AND LIBRATIONS. ANGLE UNIT IS ALWAYS RADIANS.) +C +C BARY LOGICAL FLAG DEFINING OUTPUT CENTER. +C ONLY THE 9 PLANETS ARE AFFECTED. +C BARY = .TRUE. =\ CENTER IS SOLAR-SYSTEM BARYCENTER +C = .FALSE. =\ CENTER IS SUN +C DEFAULT VALUE = .FALSE. +C +C PVSUN DP 6-WORD ARRAY CONTAINING THE BARYCENTRIC POSITION AND +C VELOCITY OF THE SUN. +C +C + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + SAVE + + INTEGER OLDMAX + PARAMETER ( OLDMAX = 400) + INTEGER NMAX + PARAMETER ( NMAX = 1000) + + DIMENSION ET2(2),PV(6,11),PNUT(4),T(2),PJD(4),BUF(1500), + . SS(3),CVAL(NMAX),PVSUN(6) + + INTEGER LIST(12),IPT(3,13) + + LOGICAL FIRST + DATA FIRST/.TRUE./ + + CHARACTER*6 TTL(14,3),CNAM(NMAX) + CHARACTER*80 NAMFIL + + LOGICAL KM,BARY + + COMMON/EPHHDR/CVAL,SS,AU,EMRAT,NUMDE,NCON,IPT + COMMON/CHRHDR/CNAM,TTL + COMMON/STCOMX/KM,BARY,PVSUN + +C +C ENTRY POINT - 1ST TIME IN, GET POINTER DATA, ETC., FROM EPH FILE +C + IF(FIRST) THEN + FIRST=.FALSE. + +C ************************************************************************ +C ************************************************************************ + +C THE USER MUST SELECT ONE OF THE FOLLOWING BY DELETING THE 'C' IN COLUMN 1 + +C ************************************************************************ + +C CALL FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL) +C CALL FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL) + CALL FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL) + + IF(NRECL .EQ. 0) WRITE(*,*)' ***** FSIZER IS NOT WORKING *****' + +C ************************************************************************ +C ************************************************************************ + + IRECSZ=NRECL*KSIZE + NCOEFFS=KSIZE/2 + + OPEN(NRFILE, + * FILE=NAMFIL, + * ACCESS='DIRECT', + * FORM='UNFORMATTED', + * RECL=IRECSZ, + * STATUS='OLD') + + READ(NRFILE,REC=1)TTL,(CNAM(K),K=1,OLDMAX),SS,NCON,AU,EMRAT, + & ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3) + & ,(CNAM(L),L=OLDMAX+1,NCON) + + IF(NCON .LE. OLDMAX)THEN + READ(NRFILE,REC=2)(CVAL(I),I=1,OLDMAX) + ELSE + READ(NRFILE,REC=2)(CVAL(I),I=1,NCON) + ENDIF + + NRL=0 + + ENDIF + +C ********** MAIN ENTRY POINT ********** + + IF(ET2(1) .EQ. 0.D0) RETURN + + S=ET2(1)-.5D0 + CALL SPLIT(S,PJD(1)) + CALL SPLIT(ET2(2),PJD(3)) + PJD(1)=PJD(1)+PJD(3)+.5D0 + PJD(2)=PJD(2)+PJD(4) + CALL SPLIT(PJD(2),PJD(3)) + PJD(1)=PJD(1)+PJD(3) + +C ERROR RETURN FOR EPOCH OUT OF RANGE + + IF(PJD(1)+PJD(4).LT.SS(1) .OR. PJD(1)+PJD(4).GT.SS(2)) GO TO 98 + +C CALCULATE RECORD # AND RELATIVE TIME IN INTERVAL + + NR=IDINT((PJD(1)-SS(1))/SS(3))+3 + IF(PJD(1).EQ.SS(2)) NR=NR-1 + + tmp1 = DBLE(NR-3)*SS(3) + SS(1) + tmp2 = PJD(1) - tmp1 + T(1) = (tmp2 + PJD(4))/SS(3) + +C READ CORRECT RECORD IF NOT IN CORE + + IF(NR.NE.NRL) THEN + NRL=NR + READ(NRFILE,REC=NR,ERR=99)(BUF(K),K=1,NCOEFFS) + ENDIF + + IF(KM) THEN + T(2)=SS(3)*86400.D0 + AUFAC=1.D0 + ELSE + T(2)=SS(3) + AUFAC=1.D0/AU + ENDIF + +C INTERPOLATE SSBARY SUN + + CALL INTERP(BUF(IPT(1,11)),T,IPT(2,11),3,IPT(3,11),2,PVSUN) + + DO I=1,6 + PVSUN(I)=PVSUN(I)*AUFAC + ENDDO + +C CHECK AND INTERPOLATE WHICHEVER BODIES ARE REQUESTED + + DO 4 I=1,10 + IF(LIST(I).EQ.0) GO TO 4 + + CALL INTERP(BUF(IPT(1,I)),T,IPT(2,I),3,IPT(3,I), + & LIST(I),PV(1,I)) + + DO J=1,6 + IF(I.LE.9 .AND. .NOT.BARY) THEN + PV(J,I)=PV(J,I)*AUFAC-PVSUN(J) + ELSE + PV(J,I)=PV(J,I)*AUFAC + ENDIF + ENDDO + + 4 CONTINUE + +C DO NUTATIONS IF REQUESTED (AND IF ON FILE) + + IF(LIST(11).GT.0 .AND. IPT(2,12).GT.0) + * CALL INTERP(BUF(IPT(1,12)),T,IPT(2,12),2,IPT(3,12), + * LIST(11),PNUT) + +C GET LIBRATIONS IF REQUESTED (AND IF ON FILE) + + IF(LIST(12).GT.0 .AND. IPT(2,13).GT.0) + * CALL INTERP(BUF(IPT(1,13)),T,IPT(2,13),3,IPT(3,13), + * LIST(12),PV(1,11)) + + RETURN + + 98 WRITE(*,198)ET2(1)+ET2(2),SS(1),SS(2) + 198 FORMAT(' *** Requested JED,',f12.2, + * ' not within ephemeris limits,',2f12.2,' ***') + + STOP + + 99 WRITE(*,'(2F12.2,A80)')ET2,'ERROR RETURN IN STATE' + + STOP + + END +C+++++++++++++++++++++++++++++ +C + SUBROUTINE CONST(NAM,VAL,SSS,N) +C +C+++++++++++++++++++++++++++++ +C +C THIS ENTRY OBTAINS THE CONSTANTS FROM THE EPHEMERIS FILE +C +C CALLING SEQEUNCE PARAMETERS (ALL OUTPUT): +C +C NAM = CHARACTER*6 ARRAY OF CONSTANT NAMES +C +C VAL = D.P. ARRAY OF VALUES OF CONSTANTS +C +C SSS = D.P. JD START, JD STOP, STEP OF EPHEMERIS +C +C N = INTEGER NUMBER OF ENTRIES IN 'NAM' AND 'VAL' ARRAYS +C + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + SAVE + + INTEGER NMAX + PARAMETER (NMAX = 1000) + + CHARACTER*6 NAM(*),TTL(14,3),CNAM(NMAX) + + DOUBLE PRECISION VAL(*),SSS(3),SS(3),CVAL(NMAX),ZIPS(2) + DOUBLE PRECISION PVST(6,11),PNUT(4) + DATA ZIPS/2*0.d0/ + + INTEGER IPT(3,13),DENUM,LIST(12) + logical first + data first/.true./ + + COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT + COMMON/CHRHDR/CNAM,TTL + +C CALL STATE TO INITIALIZE THE EPHEMERIS AND READ IN THE CONSTANTS + + IF(FIRST) CALL STATE(ZIPS,LIST,PVST,PNUT) + first=.false. + + N=NCON + + DO I=1,3 + SSS(I)=SS(I) + ENDDO + + DO I=1,N + NAM(I)=CNAM(I) + VAL(I)=CVAL(I) + ENDDO + + RETURN + + END diff --git a/lib/libration.f90 b/lib/libration.f90 new file mode 100644 index 000000000..7330d0497 --- /dev/null +++ b/lib/libration.f90 @@ -0,0 +1,38 @@ +subroutine libration(jd,RA,Dec,xl,b) + +! Compute optical libration of moon at jd: that is, the sub-observer +! point (xl,b) in selenographic coordinates. RA and Dec are +! topocentric values. + + implicit real*8 (a-h,o-z) + parameter (RADS=0.0174532925199433d0) + parameter (TWOPI=6.28318530717959d0) + real*8 jd,j2000,mjd,lambda + + j2000=2451545.0d0 + RA2000=RA + Dec2000=Dec + year=2000.0d0 + (jd-j2000)/365.25d0 + mjd=jd-2400000.d0 + call sla_PRECES('FK5',year,2000.d0,RA2000,Dec2000) + call sla_EQECL(RA2000,Dec2000,mjd,lambda,beta) + day=jd - j2000 + t = day / 36525.d0 + xi = 1.54242 * RADS + ft = 93.2720993 + 483202.0175273 * t - .0034029 * t * t + b= ft / 360 + a = 360 * (b - floor(b)) + if (a.lt.0.d0) a = 360 + a; + f=a/57.2957795131d0 + omega=sla_dranrm(2.182439196d0 - t*33.7570446085d0 + t*t*3.6236526d-5) + w = lambda - omega + y = sin(w) * cos(beta) * cos(xi) - sin(beta) * sin(xi) + x = cos(w) * cos(beta) + a = sla_dranrm(atan2(y, x)) + xl = a - f + if(xl.lt.-0.25*TWOPI) xl=xl+TWOPI !Fix 'round the back' angles + if(xl.gt.0.25*TWOPI) xl=xl-TWOPI + b = asin(-sin(w) * cos(beta) * sin(xi) - sin(beta) * cos(xi)) + + return +end subroutine libration diff --git a/lib/moon2.f90 b/lib/moon2.f90 deleted file mode 100644 index 8498267e5..000000000 --- a/lib/moon2.f90 +++ /dev/null @@ -1,166 +0,0 @@ -subroutine moon2(y,m,Day,UT,lon,lat,RA,Dec,topRA,topDec, & - LST,HA,Az,El,dist) - - implicit none - - integer y !Year - integer m !Month - integer Day !Day - real*8 UT !UTC in hours - real*8 RA,Dec !RA and Dec of moon - -! NB: Double caps are single caps in the writeup. - - real*8 NN !Longitude of ascending node - real*8 i !Inclination to the ecliptic - real*8 w !Argument of perigee - real*8 a !Semi-major axis - real*8 e !Eccentricity - real*8 MM !Mean anomaly - - real*8 v !True anomaly - real*8 EE !Eccentric anomaly - real*8 ecl !Obliquity of the ecliptic - - real*8 d !Ephemeris time argument in days - real*8 r !Distance to sun, AU - real*8 xv,yv !x and y coords in ecliptic - real*8 lonecl,latecl !Ecliptic long and lat of moon - real*8 xg,yg,zg !Ecliptic rectangular coords - real*8 Ms !Mean anomaly of sun - real*8 ws !Argument of perihelion of sun - real*8 Ls !Mean longitude of sun (Ns=0) - real*8 Lm !Mean longitude of moon - real*8 DD !Mean elongation of moon - real*8 FF !Argument of latitude for moon - real*8 xe,ye,ze !Equatorial geocentric coords of moon - real*8 mpar !Parallax of moon (r_E / d) - real*8 lat,lon !Station coordinates on earth - real*8 gclat !Geocentric latitude - real*8 rho !Earth radius factor - real*8 GMST0,LST,HA - real*8 g - real*8 topRA,topDec !Topocentric coordinates of Moon - real*8 Az,El - real*8 dist - - real*8 rad,twopi,pi,pio2 - data rad/57.2957795131d0/,twopi/6.283185307d0/ - - d=367*y - 7*(y+(m+9)/12)/4 + 275*m/9 + Day - 730530 + UT/24.d0 - ecl = 23.4393d0 - 3.563d-7 * d - -! Orbital elements for Moon: - NN = 125.1228d0 - 0.0529538083d0 * d - i = 5.1454d0 - w = mod(318.0634d0 + 0.1643573223d0 * d + 360000.d0,360.d0) - a = 60.2666d0 - e = 0.054900d0 - MM = mod(115.3654d0 + 13.0649929509d0 * d + 360000.d0,360.d0) - - EE = MM + e*rad*sin(MM/rad) * (1.d0 + e*cos(MM/rad)) - EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.d0 - e*cos(EE/rad)) - EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.d0 - e*cos(EE/rad)) - - xv = a * (cos(EE/rad) - e) - yv = a * (sqrt(1.d0-e*e) * sin(EE/rad)) - - v = mod(rad*atan2(yv,xv)+720.d0,360.d0) - r = sqrt(xv*xv + yv*yv) - -! Get geocentric position in ecliptic rectangular coordinates: - - xg = r * (cos(NN/rad)*cos((v+w)/rad) - & - sin(NN/rad)*sin((v+w)/rad)*cos(i/rad)) - yg = r * (sin(NN/rad)*cos((v+w)/rad) + & - cos(NN/rad)*sin((v+w)/rad)*cos(i/rad)) - zg = r * (sin((v+w)/rad)*sin(i/rad)) - -! Ecliptic longitude and latitude of moon: - lonecl = mod(rad*atan2(yg/rad,xg/rad)+720.d0,360.d0) - latecl = rad*atan2(zg/rad,sqrt(xg*xg + yg*yg)/rad) - -! Now include orbital perturbations: - Ms = mod(356.0470d0 + 0.9856002585d0 * d + 3600000.d0,360.d0) - ws = 282.9404d0 + 4.70935d-5*d - Ls = mod(Ms + ws + 720.d0,360.d0) - Lm = mod(MM + w + NN+720.d0,360.d0) - DD = mod(Lm - Ls + 360.d0,360.d0) - FF = mod(Lm - NN + 360.d0,360.d0) - - lonecl = lonecl & - -1.274d0 * sin((MM-2.d0*DD)/rad) & - +0.658d0 * sin(2.d0*DD/rad) & - -0.186d0 * sin(Ms/rad) & - -0.059d0 * sin((2.d0*MM-2.d0*DD)/rad) & - -0.057d0 * sin((MM-2.d0*DD+Ms)/rad) & - +0.053d0 * sin((MM+2.d0*DD)/rad) & - +0.046d0 * sin((2.d0*DD-Ms)/rad) & - +0.041d0 * sin((MM-Ms)/rad) & - -0.035d0 * sin(DD/rad) & - -0.031d0 * sin((MM+Ms)/rad) & - -0.015d0 * sin((2.d0*FF-2.d0*DD)/rad) & - +0.011d0 * sin((MM-4.d0*DD)/rad) - - latecl = latecl & - -0.173d0 * sin((FF-2.d0*DD)/rad) & - -0.055d0 * sin((MM-FF-2.d0*DD)/rad) & - -0.046d0 * sin((MM+FF-2.d0*DD)/rad) & - +0.033d0 * sin((FF+2.d0*DD)/rad) & - +0.017d0 * sin((2.d0*MM+FF)/rad) - - r = 60.36298d0 & - - 3.27746d0*cos(MM/rad) & - - 0.57994d0*cos((MM-2.d0*DD)/rad) & - - 0.46357d0*cos(2.d0*DD/rad) & - - 0.08904d0*cos(2.d0*MM/rad) & - + 0.03865d0*cos((2.d0*MM-2.d0*DD)/rad) & - - 0.03237d0*cos((2.d0*DD-Ms)/rad) & - - 0.02688d0*cos((MM+2.d0*DD)/rad) & - - 0.02358d0*cos((MM-2.d0*DD+Ms)/rad) & - - 0.02030d0*cos((MM-Ms)/rad) & - + 0.01719d0*cos(DD/rad) & - + 0.01671d0*cos((MM+Ms)/rad) - - dist=r*6378.140d0 - -! Geocentric coordinates: -! Rectangular ecliptic coordinates of the moon: - - xg = r * cos(lonecl/rad)*cos(latecl/rad) - yg = r * sin(lonecl/rad)*cos(latecl/rad) - zg = r * sin(latecl/rad) - -! Rectangular equatorial coordinates of the moon: - xe = xg - ye = yg*cos(ecl/rad) - zg*sin(ecl/rad) - ze = yg*sin(ecl/rad) + zg*cos(ecl/rad) - -! Right Ascension, Declination: - RA = mod(rad*atan2(ye,xe)+360.d0,360.d0) - Dec = rad*atan2(ze,sqrt(xe*xe + ye*ye)) - -! Now convert to topocentric system: - mpar=rad*asin(1.d0/r) -! alt_topoc = alt_geoc - mpar*cos(alt_geoc) - gclat = lat - 0.1924d0*sin(2.d0*lat/rad) - rho = 0.99883d0 + 0.00167d0*cos(2.d0*lat/rad) - GMST0 = (Ls + 180.d0)/15.d0 - LST = mod(GMST0+UT+lon/15.d0+48.d0,24.d0) !LST in hours - HA = 15.d0*LST - RA !HA in degrees - g = rad*atan(tan(gclat/rad)/cos(HA/rad)) - topRA = RA - mpar*rho*cos(gclat/rad)*sin(HA/rad)/cos(Dec/rad) - topDec = Dec - mpar*rho*sin(gclat/rad)*sin((g-Dec)/rad)/sin(g/rad) - - HA = 15.d0*LST - topRA !HA in degrees - if(HA.gt.180.d0) HA=HA-360.d0 - if(HA.lt.-180.d0) HA=HA+360.d0 - - pi=0.5d0*twopi - pio2=0.5d0*pi - call dcoord(pi,pio2-lat/rad,0.d0,lat/rad,ha*twopi/360,topDec/rad,az,el) - Az=az*rad - El=El*rad - - return -end subroutine moon2 diff --git a/lib/moondop.f90 b/lib/moondop.f90 deleted file mode 100644 index 210565866..000000000 --- a/lib/moondop.f90 +++ /dev/null @@ -1,72 +0,0 @@ -subroutine MoonDop(nyear,month,nday,uth4,lon4,lat4,RAMoon4, & - DecMoon4,LST4,HA4,AzMoon4,ElMoon4,vr4,dist4) - - implicit real*8 (a-h,o-z) - real*4 uth4 !UT in hours - real*4 lon4 !West longitude, degrees - real*4 lat4 !Latitude, degrees - real*4 RAMoon4 !Topocentric RA of moon, hours - real*4 DecMoon4 !Topocentric Dec of Moon, degrees - real*4 LST4 !Locat sidereal time, hours - real*4 HA4 !Local Hour angle, degrees - real*4 AzMoon4 !Topocentric Azimuth of moon, degrees - real*4 ElMoon4 !Topocentric Elevation of moon, degrees - real*4 vr4 !Radial velocity of moon wrt obs, km/s - real*4 dist4 !Echo time, seconds - - real*8 LST - real*8 RME(6) !Vector from Earth center to Moon - real*8 RAE(6) !Vector from Earth center to Obs - real*8 RMA(6) !Vector from Obs to Moon - real*8 rme0(6) - logical km - - data rad/57.2957795130823d0/,twopi/6.28310530717959d0/ - - km=.true. - dlat=lat4/rad - dlong1=lon4/rad - elev1=200.d0 - call geocentric(dlat,elev1,dlat1,erad1) - - dt=100.d0 !For numerical derivative, in seconds - UT=uth4 - -! NB: geodetic latitude used here, but geocentric latitude used when -! determining Earth-rotation contribution to Doppler. - - call moon2(nyear,month,nDay,UT-dt/3600.d0,dlong1*rad,dlat*rad, & - RA,Dec,topRA,topDec,LST,HA,Az0,El0,dist) - call toxyz(RA/rad,Dec/rad,dist,rme0) !Convert to rectangular coords - - call moon2(nyear,month,nDay,UT,dlong1*rad,dlat*rad, & - RA,Dec,topRA,topDec,LST,HA,Az,El,dist) - call toxyz(RA/rad,Dec/rad,dist,rme) !Convert to rectangular coords - - phi=LST*twopi/24.d0 - call toxyz(phi,dlat1,erad1,rae) !Gencentric numbers used here! - radps=twopi/(86400.d0/1.002737909d0) - rae(4)=-rae(2)*radps !Vel of Obs wrt Earth center - rae(5)=rae(1)*radps - rae(6)=0.d0 - - do i=1,3 - rme(i+3)=(rme(i)-rme0(i))/dt - rma(i)=rme(i)-rae(i) - rma(i+3)=rme(i+3)-rae(i+3) - enddo - - call fromxyz(rma,alpha1,delta1,dtopo0) !Get topocentric coords - vr=dot(rma(4),rma)/dtopo0 - - RAMoon4=topRA - DecMoon4=topDec - LST4=LST - HA4=HA - AzMoon4=Az - ElMoon4=El - vr4=vr - dist4=dist - - return -end subroutine MoonDop diff --git a/lib/moondopjpl.f90 b/lib/moondopjpl.f90 new file mode 100644 index 000000000..c97b8ff3e --- /dev/null +++ b/lib/moondopjpl.f90 @@ -0,0 +1,43 @@ +subroutine MoonDopJPL(nyear,month,nday,uth4,lon4,lat4,RAMoon4, & + DecMoon4,LST4,HA4,AzMoon4,ElMoon4,vr4,dist4) + + implicit real*8 (a-h,o-z) + real*4 uth4 !UT in hours + real*4 lon4 !East longitude, degrees + real*4 lat4 !Latitude, degrees + real*4 RAMoon4 !Topocentric RA of moon, hours + real*4 DecMoon4 !Topocentric Dec of Moon, degrees + real*4 LST4 !Locat sidereal time, hours + real*4 HA4 !Local Hour angle, degrees + real*4 AzMoon4 !Topocentric Azimuth of moon, degrees + real*4 ElMoon4 !Topocentric Elevation of moon, degrees + real*4 vr4 !Radial velocity of moon wrt obs, km/s + real*4 dist4 !Echo time, seconds + + twopi=8.d0*atan(1.d0) !Define some constants + rad=360.d0/twopi + clight=2.99792458d5 + + call sla_CLDJ(nyear,month,nday,djutc,j) + djutc=djutc + uth4/24.d0 + dut=-0.460d0 + + east_long=lon4/rad + geodetic_lat=lat4/rad + height=40. + nspecial=0 + + call ephem(djutc,dut,east_long,geodetic_lat,height,nspecial, & + RA,Dec,Az,El,techo,dop,fspread_1GHz,vr) + + RAMoon4=RA + DecMoon4=Dec + LST4=LST + HA4=HA + AzMoon4=Az*rad + ElMoon4=El*rad + vr4=vr + dist4=techo + + return +end subroutine MoonDopJPL diff --git a/lib/slasubs.f b/lib/slasubs.f new file mode 100644 index 000000000..16081f28b --- /dev/null +++ b/lib/slasubs.f @@ -0,0 +1,3396 @@ + SUBROUTINE sla_CLDJ (IY, IM, ID, DJM, J) +*+ +* - - - - - +* C L D J +* - - - - - +* +* Gregorian Calendar to Modified Julian Date +* +* Given: +* IY,IM,ID int year, month, day in Gregorian calendar +* +* Returned: +* DJM dp modified Julian Date (JD-2400000.5) for 0 hrs +* J int status: +* 0 = OK +* 1 = bad year (MJD not computed) +* 2 = bad month (MJD not computed) +* 3 = bad day (MJD computed) +* +* The year must be -4699 (i.e. 4700BC) or later. +* +* The algorithm is adapted from Hatcher 1984 (QJRAS 25, 53-55). +* +* Last revision: 27 July 2004 +* +* Copyright P.T.Wallace. All rights reserved. +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + INTEGER IY,IM,ID + DOUBLE PRECISION DJM + INTEGER J + +* Month lengths in days + INTEGER MTAB(12) + DATA MTAB / 31,28,31,30,31,30,31,31,30,31,30,31 / + + + +* Preset status. + J = 0 + +* Validate year. + IF ( IY .LT. -4699 ) THEN + J = 1 + ELSE + +* Validate month. + IF ( IM.GE.1 .AND. IM.LE.12 ) THEN + +* Allow for leap year. + IF ( MOD(IY,4) .EQ. 0 ) THEN + MTAB(2) = 29 + ELSE + MTAB(2) = 28 + END IF + IF ( MOD(IY,100).EQ.0 .AND. MOD(IY,400).NE.0 ) + : MTAB(2) = 28 + +* Validate day. + IF ( ID.LT.1 .OR. ID.GT.MTAB(IM) ) J=3 + +* Modified Julian Date. + DJM = DBLE ( ( 1461 * ( IY - (12-IM)/10 + 4712 ) ) / 4 + : + ( 306 * MOD ( IM+9, 12 ) + 5 ) / 10 + : - ( 3 * ( ( IY - (12-IM)/10 + 4900 ) / 100 ) ) / 4 + : + ID - 2399904 ) + +* Bad month. + ELSE + J=2 + END IF + + END IF + + END + DOUBLE PRECISION FUNCTION sla_DAT (UTC) +*+ +* - - - - +* D A T +* - - - - +* +* Increment to be applied to Coordinated Universal Time UTC to give +* International Atomic Time TAI (double precision) +* +* Given: +* UTC d UTC date as a modified JD (JD-2400000.5) +* +* Result: TAI-UTC in seconds +* +* Notes: +* +* 1 The UTC is specified to be a date rather than a time to indicate +* that care needs to be taken not to specify an instant which lies +* within a leap second. Though in most cases UTC can include the +* fractional part, correct behaviour on the day of a leap second +* can only be guaranteed up to the end of the second 23:59:59. +* +* 2 For epochs from 1961 January 1 onwards, the expressions from the +* file ftp://maia.usno.navy.mil/ser7/tai-utc.dat are used. +* +* 3 The 5ms time step at 1961 January 1 is taken from 2.58.1 (p87) of +* the 1992 Explanatory Supplement. +* +* 4 UTC began at 1960 January 1.0 (JD 2436934.5) and it is improper +* to call the routine with an earlier epoch. However, if this +* is attempted, the TAI-UTC expression for the year 1960 is used. +* +* +* :-----------------------------------------: +* : : +* : IMPORTANT : +* : : +* : This routine must be updated on each : +* : occasion that a leap second is : +* : announced : +* : : +* : Latest leap second: 2015 July 1 : +* : : +* :-----------------------------------------: +* +* Last revision: 5 July 2008 +* +* Copyright P.T.Wallace. All rights reserved. +*- + + IMPLICIT NONE + + DOUBLE PRECISION UTC + + DOUBLE PRECISION DT + + + + IF (.FALSE.) THEN + +* - - - - - - - - - - - - - - - - - - - - - - * +* Add new code here on each occasion that a * +* leap second is announced, and update the * +* preamble comments appropriately. * +* - - - - - - - - - - - - - - - - - - - - - - * + +* 2015 July 1 + ELSE IF (UTC.GE.57204D0) THEN + DT=36D0 + +* 2012 July 1 + ELSE IF (UTC.GE.56109D0) THEN + DT=35D0 + +* 2009 January 1 + ELSE IF (UTC.GE.54832D0) THEN + DT=34D0 + +* 2006 January 1 + ELSE IF (UTC.GE.53736D0) THEN + DT=33D0 + +* 1999 January 1 + ELSE IF (UTC.GE.51179D0) THEN + DT=32D0 + +* 1997 July 1 + ELSE IF (UTC.GE.50630D0) THEN + DT=31D0 + +* 1996 January 1 + ELSE IF (UTC.GE.50083D0) THEN + DT=30D0 + +* 1994 July 1 + ELSE IF (UTC.GE.49534D0) THEN + DT=29D0 + +* 1993 July 1 + ELSE IF (UTC.GE.49169D0) THEN + DT=28D0 + +* 1992 July 1 + ELSE IF (UTC.GE.48804D0) THEN + DT=27D0 + +* 1991 January 1 + ELSE IF (UTC.GE.48257D0) THEN + DT=26D0 + +* 1990 January 1 + ELSE IF (UTC.GE.47892D0) THEN + DT=25D0 + +* 1988 January 1 + ELSE IF (UTC.GE.47161D0) THEN + DT=24D0 + +* 1985 July 1 + ELSE IF (UTC.GE.46247D0) THEN + DT=23D0 + +* 1983 July 1 + ELSE IF (UTC.GE.45516D0) THEN + DT=22D0 + +* 1982 July 1 + ELSE IF (UTC.GE.45151D0) THEN + DT=21D0 + +* 1981 July 1 + ELSE IF (UTC.GE.44786D0) THEN + DT=20D0 + +* 1980 January 1 + ELSE IF (UTC.GE.44239D0) THEN + DT=19D0 + +* 1979 January 1 + ELSE IF (UTC.GE.43874D0) THEN + DT=18D0 + +* 1978 January 1 + ELSE IF (UTC.GE.43509D0) THEN + DT=17D0 + +* 1977 January 1 + ELSE IF (UTC.GE.43144D0) THEN + DT=16D0 + +* 1976 January 1 + ELSE IF (UTC.GE.42778D0) THEN + DT=15D0 + +* 1975 January 1 + ELSE IF (UTC.GE.42413D0) THEN + DT=14D0 + +* 1974 January 1 + ELSE IF (UTC.GE.42048D0) THEN + DT=13D0 + +* 1973 January 1 + ELSE IF (UTC.GE.41683D0) THEN + DT=12D0 + +* 1972 July 1 + ELSE IF (UTC.GE.41499D0) THEN + DT=11D0 + +* 1972 January 1 + ELSE IF (UTC.GE.41317D0) THEN + DT=10D0 + +* 1968 February 1 + ELSE IF (UTC.GE.39887D0) THEN + DT=4.2131700D0+(UTC-39126D0)*0.002592D0 + +* 1966 January 1 + ELSE IF (UTC.GE.39126D0) THEN + DT=4.3131700D0+(UTC-39126D0)*0.002592D0 + +* 1965 September 1 + ELSE IF (UTC.GE.39004D0) THEN + DT=3.8401300D0+(UTC-38761D0)*0.001296D0 + +* 1965 July 1 + ELSE IF (UTC.GE.38942D0) THEN + DT=3.7401300D0+(UTC-38761D0)*0.001296D0 + +* 1965 March 1 + ELSE IF (UTC.GE.38820D0) THEN + DT=3.6401300D0+(UTC-38761D0)*0.001296D0 + +* 1965 January 1 + ELSE IF (UTC.GE.38761D0) THEN + DT=3.5401300D0+(UTC-38761D0)*0.001296D0 + +* 1964 September 1 + ELSE IF (UTC.GE.38639D0) THEN + DT=3.4401300D0+(UTC-38761D0)*0.001296D0 + +* 1964 April 1 + ELSE IF (UTC.GE.38486D0) THEN + DT=3.3401300D0+(UTC-38761D0)*0.001296D0 + +* 1964 January 1 + ELSE IF (UTC.GE.38395D0) THEN + DT=3.2401300D0+(UTC-38761D0)*0.001296D0 + +* 1963 November 1 + ELSE IF (UTC.GE.38334D0) THEN + DT=1.9458580D0+(UTC-37665D0)*0.0011232D0 + +* 1962 January 1 + ELSE IF (UTC.GE.37665D0) THEN + DT=1.8458580D0+(UTC-37665D0)*0.0011232D0 + +* 1961 August 1 + ELSE IF (UTC.GE.37512D0) THEN + DT=1.3728180D0+(UTC-37300D0)*0.001296D0 + +* 1961 January 1 + ELSE IF (UTC.GE.37300D0) THEN + DT=1.4228180D0+(UTC-37300D0)*0.001296D0 + +* Before that + ELSE + DT=1.4178180D0+(UTC-37300D0)*0.001296D0 + + END IF + + sla_DAT=DT + + END + SUBROUTINE sla_DC62S (V, A, B, R, AD, BD, RD) +*+ +* - - - - - - +* D C 6 2 S +* - - - - - - +* +* Conversion of position & velocity in Cartesian coordinates +* to spherical coordinates (double precision) +* +* Given: +* V d(6) Cartesian position & velocity vector +* +* Returned: +* A d longitude (radians) +* B d latitude (radians) +* R d radial coordinate +* AD d longitude derivative (radians per unit time) +* BD d latitude derivative (radians per unit time) +* RD d radial derivative +* +* P.T.Wallace Starlink 28 April 1996 +* +* Copyright (C) 1996 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION V(6),A,B,R,AD,BD,RD + + DOUBLE PRECISION X,Y,Z,XD,YD,ZD,RXY2,RXY,R2,XYP + + + +* Components of position/velocity vector + X=V(1) + Y=V(2) + Z=V(3) + XD=V(4) + YD=V(5) + ZD=V(6) + +* Component of R in XY plane squared + RXY2=X*X+Y*Y + +* Modulus squared + R2=RXY2+Z*Z + +* Protection against null vector + IF (R2.EQ.0D0) THEN + X=XD + Y=YD + Z=ZD + RXY2=X*X+Y*Y + R2=RXY2+Z*Z + END IF + +* Position and velocity in spherical coordinates + RXY=SQRT(RXY2) + XYP=X*XD+Y*YD + IF (RXY2.NE.0D0) THEN + A=ATAN2(Y,X) + B=ATAN2(Z,RXY) + AD=(X*YD-Y*XD)/RXY2 + BD=(ZD*RXY2-Z*XYP)/(R2*RXY) + ELSE + A=0D0 + IF (Z.NE.0D0) THEN + B=ATAN2(Z,RXY) + ELSE + B=0D0 + END IF + AD=0D0 + BD=0D0 + END IF + R=SQRT(R2) + IF (R.NE.0D0) THEN + RD=(XYP+Z*ZD)/R + ELSE + RD=0D0 + END IF + + END + SUBROUTINE sla_DCC2S (V, A, B) +*+ +* - - - - - - +* D C C 2 S +* - - - - - - +* +* Cartesian to spherical coordinates (double precision) +* +* Given: +* V d(3) x,y,z vector +* +* Returned: +* A,B d spherical coordinates in radians +* +* The spherical coordinates are longitude (+ve anticlockwise looking +* from the +ve latitude pole) and latitude. The Cartesian coordinates +* are right handed, with the x axis at zero longitude and latitude, and +* the z axis at the +ve latitude pole. +* +* If V is null, zero A and B are returned. At either pole, zero A is +* returned. +* +* Last revision: 22 July 2004 +* +* Copyright P.T.Wallace. All rights reserved. +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION V(3),A,B + + DOUBLE PRECISION X,Y,Z,R + + + X = V(1) + Y = V(2) + Z = V(3) + R = SQRT(X*X+Y*Y) + + IF (R.EQ.0D0) THEN + A = 0D0 + ELSE + A = ATAN2(Y,X) + END IF + + IF (Z.EQ.0D0) THEN + B = 0D0 + ELSE + B = ATAN2(Z,R) + END IF + + END + SUBROUTINE sla_DCS2C (A, B, V) +*+ +* - - - - - - +* D C S 2 C +* - - - - - - +* +* Spherical coordinates to direction cosines (double precision) +* +* Given: +* A,B d spherical coordinates in radians +* (RA,Dec), (long,lat) etc. +* +* Returned: +* V d(3) x,y,z unit vector +* +* The spherical coordinates are longitude (+ve anticlockwise looking +* from the +ve latitude pole) and latitude. The Cartesian coordinates +* are right handed, with the x axis at zero longitude and latitude, and +* the z axis at the +ve latitude pole. +* +* Last revision: 26 December 2004 +* +* Copyright P.T.Wallace. All rights reserved. +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION A,B,V(3) + + DOUBLE PRECISION COSB + + + COSB = COS(B) + + V(1) = COS(A)*COSB + V(2) = SIN(A)*COSB + V(3) = SIN(B) + + END + SUBROUTINE sla_DE2H (HA, DEC, PHI, AZ, EL) +*+ +* - - - - - +* D E 2 H +* - - - - - +* +* Equatorial to horizon coordinates: HA,Dec to Az,El +* +* (double precision) +* +* Given: +* HA d hour angle +* DEC d declination +* PHI d observatory latitude +* +* Returned: +* AZ d azimuth +* EL d elevation +* +* Notes: +* +* 1) All the arguments are angles in radians. +* +* 2) Azimuth is returned in the range 0-2pi; north is zero, +* and east is +pi/2. Elevation is returned in the range +* +/-pi/2. +* +* 3) The latitude must be geodetic. In critical applications, +* corrections for polar motion should be applied. +* +* 4) In some applications it will be important to specify the +* correct type of hour angle and declination in order to +* produce the required type of azimuth and elevation. In +* particular, it may be important to distinguish between +* elevation as affected by refraction, which would +* require the "observed" HA,Dec, and the elevation +* in vacuo, which would require the "topocentric" HA,Dec. +* If the effects of diurnal aberration can be neglected, the +* "apparent" HA,Dec may be used instead of the topocentric +* HA,Dec. +* +* 5) No range checking of arguments is carried out. +* +* 6) In applications which involve many such calculations, rather +* than calling the present routine it will be more efficient to +* use inline code, having previously computed fixed terms such +* as sine and cosine of latitude, and (for tracking a star) +* sine and cosine of declination. +* +* P.T.Wallace Starlink 9 July 1994 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION HA,DEC,PHI,AZ,EL + + DOUBLE PRECISION D2PI + PARAMETER (D2PI=6.283185307179586476925286766559D0) + + DOUBLE PRECISION SH,CH,SD,CD,SP,CP,X,Y,Z,R,A + + +* Useful trig functions + SH=SIN(HA) + CH=COS(HA) + SD=SIN(DEC) + CD=COS(DEC) + SP=SIN(PHI) + CP=COS(PHI) + +* Az,El as x,y,z + X=-CH*CD*SP+SD*CP + Y=-SH*CD + Z=CH*CD*CP+SD*SP + +* To spherical + R=SQRT(X*X+Y*Y) + IF (R.EQ.0D0) THEN + A=0D0 + ELSE + A=ATAN2(Y,X) + END IF + IF (A.LT.0D0) A=A+D2PI + AZ=A + EL=ATAN2(Z,R) + + END + SUBROUTINE sla_DEULER (ORDER, PHI, THETA, PSI, RMAT) +*+ +* - - - - - - - +* D E U L E R +* - - - - - - - +* +* Form a rotation matrix from the Euler angles - three successive +* rotations about specified Cartesian axes (double precision) +* +* Given: +* ORDER c*(*) specifies about which axes the rotations occur +* PHI d 1st rotation (radians) +* THETA d 2nd rotation ( " ) +* PSI d 3rd rotation ( " ) +* +* Returned: +* RMAT d(3,3) rotation matrix +* +* A rotation is positive when the reference frame rotates +* anticlockwise as seen looking towards the origin from the +* positive region of the specified axis. +* +* The characters of ORDER define which axes the three successive +* rotations are about. A typical value is 'ZXZ', indicating that +* RMAT is to become the direction cosine matrix corresponding to +* rotations of the reference frame through PHI radians about the +* old Z-axis, followed by THETA radians about the resulting X-axis, +* then PSI radians about the resulting Z-axis. +* +* The axis names can be any of the following, in any order or +* combination: X, Y, Z, uppercase or lowercase, 1, 2, 3. Normal +* axis labelling/numbering conventions apply; the xyz (=123) +* triad is right-handed. Thus, the 'ZXZ' example given above +* could be written 'zxz' or '313' (or even 'ZxZ' or '3xZ'). ORDER +* is terminated by length or by the first unrecognized character. +* +* Fewer than three rotations are acceptable, in which case the later +* angle arguments are ignored. If all rotations are zero, the +* identity matrix is produced. +* +* P.T.Wallace Starlink 23 May 1997 +* +* Copyright (C) 1997 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + CHARACTER*(*) ORDER + DOUBLE PRECISION PHI,THETA,PSI,RMAT(3,3) + + INTEGER J,I,L,N,K + DOUBLE PRECISION RESULT(3,3),ROTN(3,3),ANGLE,S,C,W,WM(3,3) + CHARACTER AXIS + + + +* Initialize result matrix + DO J=1,3 + DO I=1,3 + IF (I.NE.J) THEN + RESULT(I,J) = 0D0 + ELSE + RESULT(I,J) = 1D0 + END IF + END DO + END DO + +* Establish length of axis string + L = LEN(ORDER) + +* Look at each character of axis string until finished + DO N=1,3 + IF (N.LE.L) THEN + +* Initialize rotation matrix for the current rotation + DO J=1,3 + DO I=1,3 + IF (I.NE.J) THEN + ROTN(I,J) = 0D0 + ELSE + ROTN(I,J) = 1D0 + END IF + END DO + END DO + +* Pick up the appropriate Euler angle and take sine & cosine + IF (N.EQ.1) THEN + ANGLE = PHI + ELSE IF (N.EQ.2) THEN + ANGLE = THETA + ELSE + ANGLE = PSI + END IF + S = SIN(ANGLE) + C = COS(ANGLE) + +* Identify the axis + AXIS = ORDER(N:N) + IF (AXIS.EQ.'X'.OR. + : AXIS.EQ.'x'.OR. + : AXIS.EQ.'1') THEN + +* Matrix for x-rotation + ROTN(2,2) = C + ROTN(2,3) = S + ROTN(3,2) = -S + ROTN(3,3) = C + + ELSE IF (AXIS.EQ.'Y'.OR. + : AXIS.EQ.'y'.OR. + : AXIS.EQ.'2') THEN + +* Matrix for y-rotation + ROTN(1,1) = C + ROTN(1,3) = -S + ROTN(3,1) = S + ROTN(3,3) = C + + ELSE IF (AXIS.EQ.'Z'.OR. + : AXIS.EQ.'z'.OR. + : AXIS.EQ.'3') THEN + +* Matrix for z-rotation + ROTN(1,1) = C + ROTN(1,2) = S + ROTN(2,1) = -S + ROTN(2,2) = C + + ELSE + +* Unrecognized character - fake end of string + L = 0 + + END IF + +* Apply the current rotation (matrix ROTN x matrix RESULT) + DO I=1,3 + DO J=1,3 + W = 0D0 + DO K=1,3 + W = W+ROTN(I,K)*RESULT(K,J) + END DO + WM(I,J) = W + END DO + END DO + DO J=1,3 + DO I=1,3 + RESULT(I,J) = WM(I,J) + END DO + END DO + + END IF + + END DO + +* Copy the result + DO J=1,3 + DO I=1,3 + RMAT(I,J) = RESULT(I,J) + END DO + END DO + + END + SUBROUTINE sla_DMOON (DATE, PV) +*+ +* - - - - - - +* D M O O N +* - - - - - - +* +* Approximate geocentric position and velocity of the Moon +* (double precision) +* +* Given: +* DATE D TDB (loosely ET) as a Modified Julian Date +* (JD-2400000.5) +* +* Returned: +* PV D(6) Moon x,y,z,xdot,ydot,zdot, mean equator and +* equinox of date (AU, AU/s) +* +* Notes: +* +* 1 This routine is a full implementation of the algorithm +* published by Meeus (see reference). +* +* 2 Meeus quotes accuracies of 10 arcsec in longitude, 3 arcsec in +* latitude and 0.2 arcsec in HP (equivalent to about 20 km in +* distance). Comparison with JPL DE200 over the interval +* 1960-2025 gives RMS errors of 3.7 arcsec and 83 mas/hour in +* longitude, 2.3 arcsec and 48 mas/hour in latitude, 11 km +* and 81 mm/s in distance. The maximum errors over the same +* interval are 18 arcsec and 0.50 arcsec/hour in longitude, +* 11 arcsec and 0.24 arcsec/hour in latitude, 40 km and 0.29 m/s +* in distance. +* +* 3 The original algorithm is expressed in terms of the obsolete +* timescale Ephemeris Time. Either TDB or TT can be used, but +* not UT without incurring significant errors (30 arcsec at +* the present time) due to the Moon's 0.5 arcsec/sec movement. +* +* 4 The algorithm is based on pre IAU 1976 standards. However, +* the result has been moved onto the new (FK5) equinox, an +* adjustment which is in any case much smaller than the +* intrinsic accuracy of the procedure. +* +* 5 Velocity is obtained by a complete analytical differentiation +* of the Meeus model. +* +* Reference: +* Meeus, l'Astronomie, June 1984, p348. +* +* P.T.Wallace Starlink 22 January 1998 +* +* Copyright (C) 1998 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE,PV(6) + +* Degrees, arcseconds and seconds of time to radians + DOUBLE PRECISION D2R,DAS2R,DS2R + PARAMETER (D2R=0.0174532925199432957692369D0, + : DAS2R=4.848136811095359935899141D-6, + : DS2R=7.272205216643039903848712D-5) + +* Seconds per Julian century (86400*36525) + DOUBLE PRECISION CJ + PARAMETER (CJ=3155760000D0) + +* Julian epoch of B1950 + DOUBLE PRECISION B1950 + PARAMETER (B1950=1949.9997904423D0) + +* Earth equatorial radius in AU ( = 6378.137 / 149597870 ) + DOUBLE PRECISION ERADAU + PARAMETER (ERADAU=4.2635212653763D-5) + + DOUBLE PRECISION T,THETA,SINOM,COSOM,DOMCOM,WA,DWA,WB,DWB,WOM, + : DWOM,SINWOM,COSWOM,V,DV,COEFF,EMN,EMPN,DN,FN,EN, + : DEN,DTHETA,FTHETA,EL,DEL,B,DB,BF,DBF,P,DP,SP,R, + : DR,X,Y,Z,XD,YD,ZD,SEL,CEL,SB,CB,RCB,RBD,W,EPJ, + : EQCOR,EPS,SINEPS,COSEPS,ES,EC + INTEGER N,I + +* +* Coefficients for fundamental arguments +* +* at J1900: T**0, T**1, T**2, T**3 +* at epoch: T**0, T**1 +* +* Units are degrees for position and Julian centuries for time +* + +* Moon's mean longitude + DOUBLE PRECISION ELP0,ELP1,ELP2,ELP3,ELP,DELP + PARAMETER (ELP0=270.434164D0, + : ELP1=481267.8831D0, + : ELP2=-0.001133D0, + : ELP3=0.0000019D0) + +* Sun's mean anomaly + DOUBLE PRECISION EM0,EM1,EM2,EM3,EM,DEM + PARAMETER (EM0=358.475833D0, + : EM1=35999.0498D0, + : EM2=-0.000150D0, + : EM3=-0.0000033D0) + +* Moon's mean anomaly + DOUBLE PRECISION EMP0,EMP1,EMP2,EMP3,EMP,DEMP + PARAMETER (EMP0=296.104608D0, + : EMP1=477198.8491D0, + : EMP2=0.009192D0, + : EMP3=0.0000144D0) + +* Moon's mean elongation + DOUBLE PRECISION D0,D1,D2,D3,D,DD + PARAMETER (D0=350.737486D0, + : D1=445267.1142D0, + : D2=-0.001436D0, + : D3=0.0000019D0) + +* Mean distance of the Moon from its ascending node + DOUBLE PRECISION F0,F1,F2,F3,F,DF + PARAMETER (F0=11.250889D0, + : F1=483202.0251D0, + : F2=-0.003211D0, + : F3=-0.0000003D0) + +* Longitude of the Moon's ascending node + DOUBLE PRECISION OM0,OM1,OM2,OM3,OM,DOM + PARAMETER (OM0=259.183275D0, + : OM1=-1934.1420D0, + : OM2=0.002078D0, + : OM3=0.0000022D0) + +* Coefficients for (dimensionless) E factor + DOUBLE PRECISION E1,E2,E,DE,ESQ,DESQ + PARAMETER (E1=-0.002495D0,E2=-0.00000752D0) + +* Coefficients for periodic variations etc + DOUBLE PRECISION PAC,PA0,PA1 + PARAMETER (PAC=0.000233D0,PA0=51.2D0,PA1=20.2D0) + DOUBLE PRECISION PBC + PARAMETER (PBC=-0.001778D0) + DOUBLE PRECISION PCC + PARAMETER (PCC=0.000817D0) + DOUBLE PRECISION PDC + PARAMETER (PDC=0.002011D0) + DOUBLE PRECISION PEC,PE0,PE1,PE2 + PARAMETER (PEC=0.003964D0, + : PE0=346.560D0,PE1=132.870D0,PE2=-0.0091731D0) + DOUBLE PRECISION PFC + PARAMETER (PFC=0.001964D0) + DOUBLE PRECISION PGC + PARAMETER (PGC=0.002541D0) + DOUBLE PRECISION PHC + PARAMETER (PHC=0.001964D0) + DOUBLE PRECISION PIC + PARAMETER (PIC=-0.024691D0) + DOUBLE PRECISION PJC,PJ0,PJ1 + PARAMETER (PJC=-0.004328D0,PJ0=275.05D0,PJ1=-2.30D0) + DOUBLE PRECISION CW1 + PARAMETER (CW1=0.0004664D0) + DOUBLE PRECISION CW2 + PARAMETER (CW2=0.0000754D0) + +* +* Coefficients for Moon position +* +* Tx(N) = coefficient of L, B or P term (deg) +* ITx(N,1-5) = coefficients of M, M', D, F, E**n in argument +* + INTEGER NL,NB,NP + PARAMETER (NL=50,NB=45,NP=31) + DOUBLE PRECISION TL(NL),TB(NB),TP(NP) + INTEGER ITL(5,NL),ITB(5,NB),ITP(5,NP) +* +* Longitude +* M M' D F n + DATA TL( 1)/ +6.288750D0 /, + : (ITL(I, 1),I=1,5)/ +0, +1, +0, +0, 0 / + DATA TL( 2)/ +1.274018D0 /, + : (ITL(I, 2),I=1,5)/ +0, -1, +2, +0, 0 / + DATA TL( 3)/ +0.658309D0 /, + : (ITL(I, 3),I=1,5)/ +0, +0, +2, +0, 0 / + DATA TL( 4)/ +0.213616D0 /, + : (ITL(I, 4),I=1,5)/ +0, +2, +0, +0, 0 / + DATA TL( 5)/ -0.185596D0 /, + : (ITL(I, 5),I=1,5)/ +1, +0, +0, +0, 1 / + DATA TL( 6)/ -0.114336D0 /, + : (ITL(I, 6),I=1,5)/ +0, +0, +0, +2, 0 / + DATA TL( 7)/ +0.058793D0 /, + : (ITL(I, 7),I=1,5)/ +0, -2, +2, +0, 0 / + DATA TL( 8)/ +0.057212D0 /, + : (ITL(I, 8),I=1,5)/ -1, -1, +2, +0, 1 / + DATA TL( 9)/ +0.053320D0 /, + : (ITL(I, 9),I=1,5)/ +0, +1, +2, +0, 0 / + DATA TL(10)/ +0.045874D0 /, + : (ITL(I,10),I=1,5)/ -1, +0, +2, +0, 1 / + DATA TL(11)/ +0.041024D0 /, + : (ITL(I,11),I=1,5)/ -1, +1, +0, +0, 1 / + DATA TL(12)/ -0.034718D0 /, + : (ITL(I,12),I=1,5)/ +0, +0, +1, +0, 0 / + DATA TL(13)/ -0.030465D0 /, + : (ITL(I,13),I=1,5)/ +1, +1, +0, +0, 1 / + DATA TL(14)/ +0.015326D0 /, + : (ITL(I,14),I=1,5)/ +0, +0, +2, -2, 0 / + DATA TL(15)/ -0.012528D0 /, + : (ITL(I,15),I=1,5)/ +0, +1, +0, +2, 0 / + DATA TL(16)/ -0.010980D0 /, + : (ITL(I,16),I=1,5)/ +0, -1, +0, +2, 0 / + DATA TL(17)/ +0.010674D0 /, + : (ITL(I,17),I=1,5)/ +0, -1, +4, +0, 0 / + DATA TL(18)/ +0.010034D0 /, + : (ITL(I,18),I=1,5)/ +0, +3, +0, +0, 0 / + DATA TL(19)/ +0.008548D0 /, + : (ITL(I,19),I=1,5)/ +0, -2, +4, +0, 0 / + DATA TL(20)/ -0.007910D0 /, + : (ITL(I,20),I=1,5)/ +1, -1, +2, +0, 1 / + DATA TL(21)/ -0.006783D0 /, + : (ITL(I,21),I=1,5)/ +1, +0, +2, +0, 1 / + DATA TL(22)/ +0.005162D0 /, + : (ITL(I,22),I=1,5)/ +0, +1, -1, +0, 0 / + DATA TL(23)/ +0.005000D0 /, + : (ITL(I,23),I=1,5)/ +1, +0, +1, +0, 1 / + DATA TL(24)/ +0.004049D0 /, + : (ITL(I,24),I=1,5)/ -1, +1, +2, +0, 1 / + DATA TL(25)/ +0.003996D0 /, + : (ITL(I,25),I=1,5)/ +0, +2, +2, +0, 0 / + DATA TL(26)/ +0.003862D0 /, + : (ITL(I,26),I=1,5)/ +0, +0, +4, +0, 0 / + DATA TL(27)/ +0.003665D0 /, + : (ITL(I,27),I=1,5)/ +0, -3, +2, +0, 0 / + DATA TL(28)/ +0.002695D0 /, + : (ITL(I,28),I=1,5)/ -1, +2, +0, +0, 1 / + DATA TL(29)/ +0.002602D0 /, + : (ITL(I,29),I=1,5)/ +0, +1, -2, -2, 0 / + DATA TL(30)/ +0.002396D0 /, + : (ITL(I,30),I=1,5)/ -1, -2, +2, +0, 1 / + DATA TL(31)/ -0.002349D0 /, + : (ITL(I,31),I=1,5)/ +0, +1, +1, +0, 0 / + DATA TL(32)/ +0.002249D0 /, + : (ITL(I,32),I=1,5)/ -2, +0, +2, +0, 2 / + DATA TL(33)/ -0.002125D0 /, + : (ITL(I,33),I=1,5)/ +1, +2, +0, +0, 1 / + DATA TL(34)/ -0.002079D0 /, + : (ITL(I,34),I=1,5)/ +2, +0, +0, +0, 2 / + DATA TL(35)/ +0.002059D0 /, + : (ITL(I,35),I=1,5)/ -2, -1, +2, +0, 2 / + DATA TL(36)/ -0.001773D0 /, + : (ITL(I,36),I=1,5)/ +0, +1, +2, -2, 0 / + DATA TL(37)/ -0.001595D0 /, + : (ITL(I,37),I=1,5)/ +0, +0, +2, +2, 0 / + DATA TL(38)/ +0.001220D0 /, + : (ITL(I,38),I=1,5)/ -1, -1, +4, +0, 1 / + DATA TL(39)/ -0.001110D0 /, + : (ITL(I,39),I=1,5)/ +0, +2, +0, +2, 0 / + DATA TL(40)/ +0.000892D0 /, + : (ITL(I,40),I=1,5)/ +0, +1, -3, +0, 0 / + DATA TL(41)/ -0.000811D0 /, + : (ITL(I,41),I=1,5)/ +1, +1, +2, +0, 1 / + DATA TL(42)/ +0.000761D0 /, + : (ITL(I,42),I=1,5)/ -1, -2, +4, +0, 1 / + DATA TL(43)/ +0.000717D0 /, + : (ITL(I,43),I=1,5)/ -2, +1, +0, +0, 2 / + DATA TL(44)/ +0.000704D0 /, + : (ITL(I,44),I=1,5)/ -2, +1, -2, +0, 2 / + DATA TL(45)/ +0.000693D0 /, + : (ITL(I,45),I=1,5)/ +1, -2, +2, +0, 1 / + DATA TL(46)/ +0.000598D0 /, + : (ITL(I,46),I=1,5)/ -1, +0, +2, -2, 1 / + DATA TL(47)/ +0.000550D0 /, + : (ITL(I,47),I=1,5)/ +0, +1, +4, +0, 0 / + DATA TL(48)/ +0.000538D0 /, + : (ITL(I,48),I=1,5)/ +0, +4, +0, +0, 0 / + DATA TL(49)/ +0.000521D0 /, + : (ITL(I,49),I=1,5)/ -1, +0, +4, +0, 1 / + DATA TL(50)/ +0.000486D0 /, + : (ITL(I,50),I=1,5)/ +0, +2, -1, +0, 0 / +* +* Latitude +* M M' D F n + DATA TB( 1)/ +5.128189D0 /, + : (ITB(I, 1),I=1,5)/ +0, +0, +0, +1, 0 / + DATA TB( 2)/ +0.280606D0 /, + : (ITB(I, 2),I=1,5)/ +0, +1, +0, +1, 0 / + DATA TB( 3)/ +0.277693D0 /, + : (ITB(I, 3),I=1,5)/ +0, +1, +0, -1, 0 / + DATA TB( 4)/ +0.173238D0 /, + : (ITB(I, 4),I=1,5)/ +0, +0, +2, -1, 0 / + DATA TB( 5)/ +0.055413D0 /, + : (ITB(I, 5),I=1,5)/ +0, -1, +2, +1, 0 / + DATA TB( 6)/ +0.046272D0 /, + : (ITB(I, 6),I=1,5)/ +0, -1, +2, -1, 0 / + DATA TB( 7)/ +0.032573D0 /, + : (ITB(I, 7),I=1,5)/ +0, +0, +2, +1, 0 / + DATA TB( 8)/ +0.017198D0 /, + : (ITB(I, 8),I=1,5)/ +0, +2, +0, +1, 0 / + DATA TB( 9)/ +0.009267D0 /, + : (ITB(I, 9),I=1,5)/ +0, +1, +2, -1, 0 / + DATA TB(10)/ +0.008823D0 /, + : (ITB(I,10),I=1,5)/ +0, +2, +0, -1, 0 / + DATA TB(11)/ +0.008247D0 /, + : (ITB(I,11),I=1,5)/ -1, +0, +2, -1, 1 / + DATA TB(12)/ +0.004323D0 /, + : (ITB(I,12),I=1,5)/ +0, -2, +2, -1, 0 / + DATA TB(13)/ +0.004200D0 /, + : (ITB(I,13),I=1,5)/ +0, +1, +2, +1, 0 / + DATA TB(14)/ +0.003372D0 /, + : (ITB(I,14),I=1,5)/ -1, +0, -2, +1, 1 / + DATA TB(15)/ +0.002472D0 /, + : (ITB(I,15),I=1,5)/ -1, -1, +2, +1, 1 / + DATA TB(16)/ +0.002222D0 /, + : (ITB(I,16),I=1,5)/ -1, +0, +2, +1, 1 / + DATA TB(17)/ +0.002072D0 /, + : (ITB(I,17),I=1,5)/ -1, -1, +2, -1, 1 / + DATA TB(18)/ +0.001877D0 /, + : (ITB(I,18),I=1,5)/ -1, +1, +0, +1, 1 / + DATA TB(19)/ +0.001828D0 /, + : (ITB(I,19),I=1,5)/ +0, -1, +4, -1, 0 / + DATA TB(20)/ -0.001803D0 /, + : (ITB(I,20),I=1,5)/ +1, +0, +0, +1, 1 / + DATA TB(21)/ -0.001750D0 /, + : (ITB(I,21),I=1,5)/ +0, +0, +0, +3, 0 / + DATA TB(22)/ +0.001570D0 /, + : (ITB(I,22),I=1,5)/ -1, +1, +0, -1, 1 / + DATA TB(23)/ -0.001487D0 /, + : (ITB(I,23),I=1,5)/ +0, +0, +1, +1, 0 / + DATA TB(24)/ -0.001481D0 /, + : (ITB(I,24),I=1,5)/ +1, +1, +0, +1, 1 / + DATA TB(25)/ +0.001417D0 /, + : (ITB(I,25),I=1,5)/ -1, -1, +0, +1, 1 / + DATA TB(26)/ +0.001350D0 /, + : (ITB(I,26),I=1,5)/ -1, +0, +0, +1, 1 / + DATA TB(27)/ +0.001330D0 /, + : (ITB(I,27),I=1,5)/ +0, +0, -1, +1, 0 / + DATA TB(28)/ +0.001106D0 /, + : (ITB(I,28),I=1,5)/ +0, +3, +0, +1, 0 / + DATA TB(29)/ +0.001020D0 /, + : (ITB(I,29),I=1,5)/ +0, +0, +4, -1, 0 / + DATA TB(30)/ +0.000833D0 /, + : (ITB(I,30),I=1,5)/ +0, -1, +4, +1, 0 / + DATA TB(31)/ +0.000781D0 /, + : (ITB(I,31),I=1,5)/ +0, +1, +0, -3, 0 / + DATA TB(32)/ +0.000670D0 /, + : (ITB(I,32),I=1,5)/ +0, -2, +4, +1, 0 / + DATA TB(33)/ +0.000606D0 /, + : (ITB(I,33),I=1,5)/ +0, +0, +2, -3, 0 / + DATA TB(34)/ +0.000597D0 /, + : (ITB(I,34),I=1,5)/ +0, +2, +2, -1, 0 / + DATA TB(35)/ +0.000492D0 /, + : (ITB(I,35),I=1,5)/ -1, +1, +2, -1, 1 / + DATA TB(36)/ +0.000450D0 /, + : (ITB(I,36),I=1,5)/ +0, +2, -2, -1, 0 / + DATA TB(37)/ +0.000439D0 /, + : (ITB(I,37),I=1,5)/ +0, +3, +0, -1, 0 / + DATA TB(38)/ +0.000423D0 /, + : (ITB(I,38),I=1,5)/ +0, +2, +2, +1, 0 / + DATA TB(39)/ +0.000422D0 /, + : (ITB(I,39),I=1,5)/ +0, -3, +2, -1, 0 / + DATA TB(40)/ -0.000367D0 /, + : (ITB(I,40),I=1,5)/ +1, -1, +2, +1, 1 / + DATA TB(41)/ -0.000353D0 /, + : (ITB(I,41),I=1,5)/ +1, +0, +2, +1, 1 / + DATA TB(42)/ +0.000331D0 /, + : (ITB(I,42),I=1,5)/ +0, +0, +4, +1, 0 / + DATA TB(43)/ +0.000317D0 /, + : (ITB(I,43),I=1,5)/ -1, +1, +2, +1, 1 / + DATA TB(44)/ +0.000306D0 /, + : (ITB(I,44),I=1,5)/ -2, +0, +2, -1, 2 / + DATA TB(45)/ -0.000283D0 /, + : (ITB(I,45),I=1,5)/ +0, +1, +0, +3, 0 / +* +* Parallax +* M M' D F n + DATA TP( 1)/ +0.950724D0 /, + : (ITP(I, 1),I=1,5)/ +0, +0, +0, +0, 0 / + DATA TP( 2)/ +0.051818D0 /, + : (ITP(I, 2),I=1,5)/ +0, +1, +0, +0, 0 / + DATA TP( 3)/ +0.009531D0 /, + : (ITP(I, 3),I=1,5)/ +0, -1, +2, +0, 0 / + DATA TP( 4)/ +0.007843D0 /, + : (ITP(I, 4),I=1,5)/ +0, +0, +2, +0, 0 / + DATA TP( 5)/ +0.002824D0 /, + : (ITP(I, 5),I=1,5)/ +0, +2, +0, +0, 0 / + DATA TP( 6)/ +0.000857D0 /, + : (ITP(I, 6),I=1,5)/ +0, +1, +2, +0, 0 / + DATA TP( 7)/ +0.000533D0 /, + : (ITP(I, 7),I=1,5)/ -1, +0, +2, +0, 1 / + DATA TP( 8)/ +0.000401D0 /, + : (ITP(I, 8),I=1,5)/ -1, -1, +2, +0, 1 / + DATA TP( 9)/ +0.000320D0 /, + : (ITP(I, 9),I=1,5)/ -1, +1, +0, +0, 1 / + DATA TP(10)/ -0.000271D0 /, + : (ITP(I,10),I=1,5)/ +0, +0, +1, +0, 0 / + DATA TP(11)/ -0.000264D0 /, + : (ITP(I,11),I=1,5)/ +1, +1, +0, +0, 1 / + DATA TP(12)/ -0.000198D0 /, + : (ITP(I,12),I=1,5)/ +0, -1, +0, +2, 0 / + DATA TP(13)/ +0.000173D0 /, + : (ITP(I,13),I=1,5)/ +0, +3, +0, +0, 0 / + DATA TP(14)/ +0.000167D0 /, + : (ITP(I,14),I=1,5)/ +0, -1, +4, +0, 0 / + DATA TP(15)/ -0.000111D0 /, + : (ITP(I,15),I=1,5)/ +1, +0, +0, +0, 1 / + DATA TP(16)/ +0.000103D0 /, + : (ITP(I,16),I=1,5)/ +0, -2, +4, +0, 0 / + DATA TP(17)/ -0.000084D0 /, + : (ITP(I,17),I=1,5)/ +0, +2, -2, +0, 0 / + DATA TP(18)/ -0.000083D0 /, + : (ITP(I,18),I=1,5)/ +1, +0, +2, +0, 1 / + DATA TP(19)/ +0.000079D0 /, + : (ITP(I,19),I=1,5)/ +0, +2, +2, +0, 0 / + DATA TP(20)/ +0.000072D0 /, + : (ITP(I,20),I=1,5)/ +0, +0, +4, +0, 0 / + DATA TP(21)/ +0.000064D0 /, + : (ITP(I,21),I=1,5)/ -1, +1, +2, +0, 1 / + DATA TP(22)/ -0.000063D0 /, + : (ITP(I,22),I=1,5)/ +1, -1, +2, +0, 1 / + DATA TP(23)/ +0.000041D0 /, + : (ITP(I,23),I=1,5)/ +1, +0, +1, +0, 1 / + DATA TP(24)/ +0.000035D0 /, + : (ITP(I,24),I=1,5)/ -1, +2, +0, +0, 1 / + DATA TP(25)/ -0.000033D0 /, + : (ITP(I,25),I=1,5)/ +0, +3, -2, +0, 0 / + DATA TP(26)/ -0.000030D0 /, + : (ITP(I,26),I=1,5)/ +0, +1, +1, +0, 0 / + DATA TP(27)/ -0.000029D0 /, + : (ITP(I,27),I=1,5)/ +0, +0, -2, +2, 0 / + DATA TP(28)/ -0.000029D0 /, + : (ITP(I,28),I=1,5)/ +1, +2, +0, +0, 1 / + DATA TP(29)/ +0.000026D0 /, + : (ITP(I,29),I=1,5)/ -2, +0, +2, +0, 2 / + DATA TP(30)/ -0.000023D0 /, + : (ITP(I,30),I=1,5)/ +0, +1, -2, +2, 0 / + DATA TP(31)/ +0.000019D0 /, + : (ITP(I,31),I=1,5)/ -1, -1, +4, +0, 1 / + + + +* Centuries since J1900 + T=(DATE-15019.5D0)/36525D0 + +* +* Fundamental arguments (radians) and derivatives (radians per +* Julian century) for the current epoch +* + +* Moon's mean longitude + ELP=D2R*MOD(ELP0+(ELP1+(ELP2+ELP3*T)*T)*T,360D0) + DELP=D2R*(ELP1+(2D0*ELP2+3D0*ELP3*T)*T) + +* Sun's mean anomaly + EM=D2R*MOD(EM0+(EM1+(EM2+EM3*T)*T)*T,360D0) + DEM=D2R*(EM1+(2D0*EM2+3D0*EM3*T)*T) + +* Moon's mean anomaly + EMP=D2R*MOD(EMP0+(EMP1+(EMP2+EMP3*T)*T)*T,360D0) + DEMP=D2R*(EMP1+(2D0*EMP2+3D0*EMP3*T)*T) + +* Moon's mean elongation + D=D2R*MOD(D0+(D1+(D2+D3*T)*T)*T,360D0) + DD=D2R*(D1+(2D0*D2+3D0*D3*T)*T) + +* Mean distance of the Moon from its ascending node + F=D2R*MOD(F0+(F1+(F2+F3*T)*T)*T,360D0) + DF=D2R*(F1+(2D0*F2+3D0*F3*T)*T) + +* Longitude of the Moon's ascending node + OM=D2R*MOD(OM0+(OM1+(OM2+OM3*T)*T)*T,360D0) + DOM=D2R*(OM1+(2D0*OM2+3D0*OM3*T)*T) + SINOM=SIN(OM) + COSOM=COS(OM) + DOMCOM=DOM*COSOM + +* Add the periodic variations + THETA=D2R*(PA0+PA1*T) + WA=SIN(THETA) + DWA=D2R*PA1*COS(THETA) + THETA=D2R*(PE0+(PE1+PE2*T)*T) + WB=PEC*SIN(THETA) + DWB=D2R*PEC*(PE1+2D0*PE2*T)*COS(THETA) + ELP=ELP+D2R*(PAC*WA+WB+PFC*SINOM) + DELP=DELP+D2R*(PAC*DWA+DWB+PFC*DOMCOM) + EM=EM+D2R*PBC*WA + DEM=DEM+D2R*PBC*DWA + EMP=EMP+D2R*(PCC*WA+WB+PGC*SINOM) + DEMP=DEMP+D2R*(PCC*DWA+DWB+PGC*DOMCOM) + D=D+D2R*(PDC*WA+WB+PHC*SINOM) + DD=DD+D2R*(PDC*DWA+DWB+PHC*DOMCOM) + WOM=OM+D2R*(PJ0+PJ1*T) + DWOM=DOM+D2R*PJ1 + SINWOM=SIN(WOM) + COSWOM=COS(WOM) + F=F+D2R*(WB+PIC*SINOM+PJC*SINWOM) + DF=DF+D2R*(DWB+PIC*DOMCOM+PJC*DWOM*COSWOM) + +* E-factor, and square + E=1D0+(E1+E2*T)*T + DE=E1+2D0*E2*T + ESQ=E*E + DESQ=2D0*E*DE + +* +* Series expansions +* + +* Longitude + V=0D0 + DV=0D0 + DO N=NL,1,-1 + COEFF=TL(N) + EMN=DBLE(ITL(1,N)) + EMPN=DBLE(ITL(2,N)) + DN=DBLE(ITL(3,N)) + FN=DBLE(ITL(4,N)) + I=ITL(5,N) + IF (I.EQ.0) THEN + EN=1D0 + DEN=0D0 + ELSE IF (I.EQ.1) THEN + EN=E + DEN=DE + ELSE + EN=ESQ + DEN=DESQ + END IF + THETA=EMN*EM+EMPN*EMP+DN*D+FN*F + DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF + FTHETA=SIN(THETA) + V=V+COEFF*FTHETA*EN + DV=DV+COEFF*(COS(THETA)*DTHETA*EN+FTHETA*DEN) + END DO + EL=ELP+D2R*V + DEL=(DELP+D2R*DV)/CJ + +* Latitude + V=0D0 + DV=0D0 + DO N=NB,1,-1 + COEFF=TB(N) + EMN=DBLE(ITB(1,N)) + EMPN=DBLE(ITB(2,N)) + DN=DBLE(ITB(3,N)) + FN=DBLE(ITB(4,N)) + I=ITB(5,N) + IF (I.EQ.0) THEN + EN=1D0 + DEN=0D0 + ELSE IF (I.EQ.1) THEN + EN=E + DEN=DE + ELSE + EN=ESQ + DEN=DESQ + END IF + THETA=EMN*EM+EMPN*EMP+DN*D+FN*F + DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF + FTHETA=SIN(THETA) + V=V+COEFF*FTHETA*EN + DV=DV+COEFF*(COS(THETA)*DTHETA*EN+FTHETA*DEN) + END DO + BF=1D0-CW1*COSOM-CW2*COSWOM + DBF=CW1*DOM*SINOM+CW2*DWOM*SINWOM + B=D2R*V*BF + DB=D2R*(DV*BF+V*DBF)/CJ + +* Parallax + V=0D0 + DV=0D0 + DO N=NP,1,-1 + COEFF=TP(N) + EMN=DBLE(ITP(1,N)) + EMPN=DBLE(ITP(2,N)) + DN=DBLE(ITP(3,N)) + FN=DBLE(ITP(4,N)) + I=ITP(5,N) + IF (I.EQ.0) THEN + EN=1D0 + DEN=0D0 + ELSE IF (I.EQ.1) THEN + EN=E + DEN=DE + ELSE + EN=ESQ + DEN=DESQ + END IF + THETA=EMN*EM+EMPN*EMP+DN*D+FN*F + DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF + FTHETA=COS(THETA) + V=V+COEFF*FTHETA*EN + DV=DV+COEFF*(-SIN(THETA)*DTHETA*EN+FTHETA*DEN) + END DO + P=D2R*V + DP=D2R*DV/CJ + +* +* Transformation into final form +* + +* Parallax to distance (AU, AU/sec) + SP=SIN(P) + R=ERADAU/SP + DR=-R*DP*COS(P)/SP + +* Longitude, latitude to x,y,z (AU) + SEL=SIN(EL) + CEL=COS(EL) + SB=SIN(B) + CB=COS(B) + RCB=R*CB + RBD=R*DB + W=RBD*SB-CB*DR + X=RCB*CEL + Y=RCB*SEL + Z=R*SB + XD=-Y*DEL-W*CEL + YD=X*DEL-W*SEL + ZD=RBD*CB+SB*DR + +* Julian centuries since J2000 + T=(DATE-51544.5D0)/36525D0 + +* Fricke equinox correction + EPJ=2000D0+T*100D0 + EQCOR=DS2R*(0.035D0+0.00085D0*(EPJ-B1950)) + +* Mean obliquity (IAU 1976) + EPS=DAS2R*(84381.448D0+(-46.8150D0+(-0.00059D0+0.001813D0*T)*T)*T) + +* To the equatorial system, mean of date, FK5 system + SINEPS=SIN(EPS) + COSEPS=COS(EPS) + ES=EQCOR*SINEPS + EC=EQCOR*COSEPS + PV(1)=X-EC*Y+ES*Z + PV(2)=EQCOR*X+Y*COSEPS-Z*SINEPS + PV(3)=Y*SINEPS+Z*COSEPS + PV(4)=XD-EC*YD+ES*ZD + PV(5)=EQCOR*XD+YD*COSEPS-ZD*SINEPS + PV(6)=YD*SINEPS+ZD*COSEPS + + END + SUBROUTINE sla_DMXV (DM, VA, VB) +*+ +* - - - - - +* D M X V +* - - - - - +* +* Performs the 3-D forward unitary transformation: +* +* vector VB = matrix DM * vector VA +* +* (double precision) +* +* Given: +* DM dp(3,3) matrix +* VA dp(3) vector +* +* Returned: +* VB dp(3) result vector +* +* To comply with the ANSI Fortran 77 standard, VA and VB must be +* different arrays. However, the routine is coded so as to work +* properly on many platforms even if this rule is violated. +* +* Last revision: 26 December 2004 +* +* Copyright P.T.Wallace. All rights reserved. +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DM(3,3),VA(3),VB(3) + + INTEGER I,J + DOUBLE PRECISION W,VW(3) + + +* Matrix DM * vector VA -> vector VW + DO J=1,3 + W=0D0 + DO I=1,3 + W=W+DM(J,I)*VA(I) + END DO + VW(J)=W + END DO + +* Vector VW -> vector VB + DO J=1,3 + VB(J)=VW(J) + END DO + + END + DOUBLE PRECISION FUNCTION sla_DRANGE (ANGLE) +*+ +* - - - - - - - +* D R A N G E +* - - - - - - - +* +* Normalize angle into range +/- pi (double precision) +* +* Given: +* ANGLE dp the angle in radians +* +* The result (double precision) is ANGLE expressed in the range +/- pi. +* +* P.T.Wallace Starlink 23 November 1995 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION ANGLE + + DOUBLE PRECISION DPI,D2PI + PARAMETER (DPI=3.141592653589793238462643D0) + PARAMETER (D2PI=6.283185307179586476925287D0) + + + sla_DRANGE=MOD(ANGLE,D2PI) + IF (ABS(sla_DRANGE).GE.DPI) + : sla_DRANGE=sla_DRANGE-SIGN(D2PI,ANGLE) + + END + DOUBLE PRECISION FUNCTION sla_DRANRM (ANGLE) +*+ +* - - - - - - - +* D R A N R M +* - - - - - - - +* +* Normalize angle into range 0-2 pi (double precision) +* +* Given: +* ANGLE dp the angle in radians +* +* The result is ANGLE expressed in the range 0-2 pi. +* +* Last revision: 22 July 2004 +* +* Copyright P.T.Wallace. All rights reserved. +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION ANGLE + + DOUBLE PRECISION D2PI + PARAMETER (D2PI=6.283185307179586476925286766559D0) + + + sla_DRANRM = MOD(ANGLE,D2PI) + IF (sla_DRANRM.LT.0D0) sla_DRANRM = sla_DRANRM+D2PI + + END + DOUBLE PRECISION FUNCTION sla_DTT (UTC) +*+ +* - - - - +* D T T +* - - - - +* +* Increment to be applied to Coordinated Universal Time UTC to give +* Terrestrial Time TT (formerly Ephemeris Time ET) +* +* (double precision) +* +* Given: +* UTC d UTC date as a modified JD (JD-2400000.5) +* +* Result: TT-UTC in seconds +* +* Notes: +* +* 1 The UTC is specified to be a date rather than a time to indicate +* that care needs to be taken not to specify an instant which lies +* within a leap second. Though in most cases UTC can include the +* fractional part, correct behaviour on the day of a leap second +* can only be guaranteed up to the end of the second 23:59:59. +* +* 2 Pre 1972 January 1 a fixed value of 10 + ET-TAI is returned. +* +* 3 See also the routine sla_DT, which roughly estimates ET-UT for +* historical epochs. +* +* Called: sla_DAT +* +* P.T.Wallace Starlink 6 December 1994 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION UTC + + DOUBLE PRECISION sla_DAT + + + sla_DTT=32.184D0+sla_DAT(UTC) + + END + SUBROUTINE sla_ECMAT (DATE, RMAT) +*+ +* - - - - - - +* E C M A T +* - - - - - - +* +* Form the equatorial to ecliptic rotation matrix - IAU 1980 theory +* (double precision) +* +* Given: +* DATE dp TDB (loosely ET) as Modified Julian Date +* (JD-2400000.5) +* Returned: +* RMAT dp(3,3) matrix +* +* Reference: +* Murray,C.A., Vectorial Astrometry, section 4.3. +* +* Note: +* The matrix is in the sense V(ecl) = RMAT * V(equ); the +* equator, equinox and ecliptic are mean of date. +* +* Called: sla_DEULER +* +* P.T.Wallace Starlink 23 August 1996 +* +* Copyright (C) 1996 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE,RMAT(3,3) + +* Arc seconds to radians + DOUBLE PRECISION AS2R + PARAMETER (AS2R=0.484813681109535994D-5) + + DOUBLE PRECISION T,EPS0 + + + +* Interval between basic epoch J2000.0 and current epoch (JC) + T = (DATE-51544.5D0)/36525D0 + +* Mean obliquity + EPS0 = AS2R* + : (84381.448D0+(-46.8150D0+(-0.00059D0+0.001813D0*T)*T)*T) + +* Matrix + CALL sla_DEULER('X',EPS0,0D0,0D0,RMAT) + + END + DOUBLE PRECISION FUNCTION sla_EPJ (DATE) +*+ +* - - - - +* E P J +* - - - - +* +* Conversion of Modified Julian Date to Julian Epoch (double precision) +* +* Given: +* DATE dp Modified Julian Date (JD - 2400000.5) +* +* The result is the Julian Epoch. +* +* Reference: +* Lieske,J.H., 1979. Astron.Astrophys.,73,282. +* +* P.T.Wallace Starlink February 1984 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE + + + sla_EPJ = 2000D0 + (DATE-51544.5D0)/365.25D0 + + END + SUBROUTINE sla_EQECL (DR, DD, DATE, DL, DB) +*+ +* - - - - - - +* E Q E C L +* - - - - - - +* +* Transformation from J2000.0 equatorial coordinates to +* ecliptic coordinates (double precision) +* +* Given: +* DR,DD dp J2000.0 mean RA,Dec (radians) +* DATE dp TDB (loosely ET) as Modified Julian Date +* (JD-2400000.5) +* Returned: +* DL,DB dp ecliptic longitude and latitude +* (mean of date, IAU 1980 theory, radians) +* +* Called: +* sla_DCS2C, sla_PREC, sla_EPJ, sla_DMXV, sla_ECMAT, sla_DCC2S, +* sla_DRANRM, sla_DRANGE +* +* P.T.Wallace Starlink March 1986 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DR,DD,DATE,DL,DB + + DOUBLE PRECISION sla_EPJ,sla_DRANRM,sla_DRANGE + + DOUBLE PRECISION RMAT(3,3),V1(3),V2(3) + + + +* Spherical to Cartesian + CALL sla_DCS2C(DR,DD,V1) + +* Mean J2000 to mean of date + CALL sla_PREC(2000D0,sla_EPJ(DATE),RMAT) + CALL sla_DMXV(RMAT,V1,V2) + +* Equatorial to ecliptic + CALL sla_ECMAT(DATE,RMAT) + CALL sla_DMXV(RMAT,V2,V1) + +* Cartesian to spherical + CALL sla_DCC2S(V1,DL,DB) + +* Express in conventional ranges + DL=sla_DRANRM(DL) + DB=sla_DRANGE(DB) + + END + DOUBLE PRECISION FUNCTION sla_EQEQX (DATE) +*+ +* - - - - - - +* E Q E Q X +* - - - - - - +* +* Equation of the equinoxes (IAU 1994, double precision) +* +* Given: +* DATE dp TDB (loosely ET) as Modified Julian Date +* (JD-2400000.5) +* +* The result is the equation of the equinoxes (double precision) +* in radians: +* +* Greenwich apparent ST = GMST + sla_EQEQX +* +* References: IAU Resolution C7, Recommendation 3 (1994) +* Capitaine, N. & Gontier, A.-M., Astron. Astrophys., +* 275, 645-650 (1993) +* +* Called: sla_NUTC +* +* Patrick Wallace Starlink 23 August 1996 +* +* Copyright (C) 1996 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE + +* Turns to arc seconds and arc seconds to radians + DOUBLE PRECISION T2AS,AS2R + PARAMETER (T2AS=1296000D0, + : AS2R=0.484813681109535994D-5) + + DOUBLE PRECISION T,OM,DPSI,DEPS,EPS0 + + + +* Interval between basic epoch J2000.0 and current epoch (JC) + T=(DATE-51544.5D0)/36525D0 + +* Longitude of the mean ascending node of the lunar orbit on the +* ecliptic, measured from the mean equinox of date + OM=AS2R*(450160.280D0+(-5D0*T2AS-482890.539D0 + : +(7.455D0+0.008D0*T)*T)*T) + +* Nutation + CALL sla_NUTC(DATE,DPSI,DEPS,EPS0) + +* Equation of the equinoxes + sla_EQEQX=DPSI*COS(EPS0)+AS2R*(0.00264D0*SIN(OM)+ + : 0.000063D0*SIN(OM+OM)) + + END + SUBROUTINE sla_GEOC (P, H, R, Z) +*+ +* - - - - - +* G E O C +* - - - - - +* +* Convert geodetic position to geocentric (double precision) +* +* Given: +* P dp latitude (geodetic, radians) +* H dp height above reference spheroid (geodetic, metres) +* +* Returned: +* R dp distance from Earth axis (AU) +* Z dp distance from plane of Earth equator (AU) +* +* Notes: +* +* 1 Geocentric latitude can be obtained by evaluating ATAN2(Z,R). +* +* 2 IAU 1976 constants are used. +* +* Reference: +* +* Green,R.M., Spherical Astronomy, CUP 1985, p98. +* +* Last revision: 22 July 2004 +* +* Copyright P.T.Wallace. All rights reserved. +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION P,H,R,Z + +* Earth equatorial radius (metres) + DOUBLE PRECISION A0 + PARAMETER (A0=6378140D0) + +* Reference spheroid flattening factor and useful function + DOUBLE PRECISION F,B + PARAMETER (F=1D0/298.257D0,B=(1D0-F)**2) + +* Astronomical unit in metres + DOUBLE PRECISION AU + PARAMETER (AU=1.49597870D11) + + DOUBLE PRECISION SP,CP,C,S + + + +* Geodetic to geocentric conversion + SP = SIN(P) + CP = COS(P) + C = 1D0/SQRT(CP*CP+B*SP*SP) + S = B*C + R = (A0*C+H)*CP/AU + Z = (A0*S+H)*SP/AU + + END + DOUBLE PRECISION FUNCTION sla_GMST (UT1) +*+ +* - - - - - +* G M S T +* - - - - - +* +* Conversion from universal time to sidereal time (double precision) +* +* Given: +* UT1 dp universal time (strictly UT1) expressed as +* modified Julian Date (JD-2400000.5) +* +* The result is the Greenwich mean sidereal time (double +* precision, radians). +* +* The IAU 1982 expression (see page S15 of 1984 Astronomical Almanac) +* is used, but rearranged to reduce rounding errors. This expression +* is always described as giving the GMST at 0 hours UT. In fact, it +* gives the difference between the GMST and the UT, which happens to +* equal the GMST (modulo 24 hours) at 0 hours UT each day. In this +* routine, the entire UT is used directly as the argument for the +* standard formula, and the fractional part of the UT is added +* separately. Note that the factor 1.0027379... does not appear in the +* IAU 1982 expression explicitly but in the form of the coefficient +* 8640184.812866, which is 86400x36525x0.0027379... +* +* See also the routine sla_GMSTA, which delivers better numerical +* precision by accepting the UT date and time as separate arguments. +* +* Called: sla_DRANRM +* +* P.T.Wallace Starlink 14 October 2001 +* +* Copyright (C) 2001 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION UT1 + + DOUBLE PRECISION sla_DRANRM + + DOUBLE PRECISION D2PI,S2R + PARAMETER (D2PI=6.283185307179586476925286766559D0, + : S2R=7.272205216643039903848711535369D-5) + + DOUBLE PRECISION TU + + + +* Julian centuries from fundamental epoch J2000 to this UT + TU=(UT1-51544.5D0)/36525D0 + +* GMST at this UT + sla_GMST=sla_DRANRM(MOD(UT1,1D0)*D2PI+ + : (24110.54841D0+ + : (8640184.812866D0+ + : (0.093104D0-6.2D-6*TU)*TU)*TU)*S2R) + + END + SUBROUTINE sla_NUTC (DATE, DPSI, DEPS, EPS0) +*+ +* - - - - - +* N U T C +* - - - - - +* +* Nutation: longitude & obliquity components and mean obliquity, +* using the Shirai & Fukushima (2001) theory. +* +* Given: +* DATE d TDB (loosely ET) as Modified Julian Date +* (JD-2400000.5) +* Returned: +* DPSI,DEPS d nutation in longitude,obliquity +* EPS0 d mean obliquity +* +* Notes: +* +* 1 The routine predicts forced nutation (but not free core nutation) +* plus corrections to the IAU 1976 precession model. +* +* 2 Earth attitude predictions made by combining the present nutation +* model with IAU 1976 precession are accurate to 1 mas (with respect +* to the ICRF) for a few decades around 2000. +* +* 3 The sla_NUTC80 routine is the equivalent of the present routine +* but using the IAU 1980 nutation theory. The older theory is less +* accurate, leading to errors as large as 350 mas over the interval +* 1900-2100, mainly because of the error in the IAU 1976 precession. +* +* References: +* +* Shirai, T. & Fukushima, T., Astron.J. 121, 3270-3283 (2001). +* +* Fukushima, T., Astron.Astrophys. 244, L11 (1991). +* +* Simon, J. L., Bretagnon, P., Chapront, J., Chapront-Touze, M., +* Francou, G. & Laskar, J., Astron.Astrophys. 282, 663 (1994). +* +* This revision: 24 November 2005 +* +* Copyright P.T.Wallace. All rights reserved. +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE,DPSI,DEPS,EPS0 + +* Degrees to radians + DOUBLE PRECISION DD2R + PARAMETER (DD2R=1.745329251994329576923691D-2) + +* Arc seconds to radians + DOUBLE PRECISION DAS2R + PARAMETER (DAS2R=4.848136811095359935899141D-6) + +* Arc seconds in a full circle + DOUBLE PRECISION TURNAS + PARAMETER (TURNAS=1296000D0) + +* Reference epoch (J2000), MJD + DOUBLE PRECISION DJM0 + PARAMETER (DJM0=51544.5D0 ) + +* Days per Julian century + DOUBLE PRECISION DJC + PARAMETER (DJC=36525D0) + + INTEGER I,J + DOUBLE PRECISION T,EL,ELP,F,D,OM,VE,MA,JU,SA,THETA,C,S,DP,DE + +* Number of terms in the nutation model + INTEGER NTERMS + PARAMETER (NTERMS=194) + +* The SF2001 forced nutation model + INTEGER NA(9,NTERMS) + DOUBLE PRECISION PSI(4,NTERMS), EPS(4,NTERMS) + +* Coefficients of fundamental angles + DATA ( ( NA(I,J), I=1,9 ), J=1,10 ) / + : 0, 0, 0, 0, -1, 0, 0, 0, 0, + : 0, 0, 2, -2, 2, 0, 0, 0, 0, + : 0, 0, 2, 0, 2, 0, 0, 0, 0, + : 0, 0, 0, 0, -2, 0, 0, 0, 0, + : 0, 1, 0, 0, 0, 0, 0, 0, 0, + : 0, 1, 2, -2, 2, 0, 0, 0, 0, + : 1, 0, 0, 0, 0, 0, 0, 0, 0, + : 0, 0, 2, 0, 1, 0, 0, 0, 0, + : 1, 0, 2, 0, 2, 0, 0, 0, 0, + : 0, -1, 2, -2, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=11,20 ) / + : 0, 0, 2, -2, 1, 0, 0, 0, 0, + : -1, 0, 2, 0, 2, 0, 0, 0, 0, + : -1, 0, 0, 2, 0, 0, 0, 0, 0, + : 1, 0, 0, 0, 1, 0, 0, 0, 0, + : 1, 0, 0, 0, -1, 0, 0, 0, 0, + : -1, 0, 2, 2, 2, 0, 0, 0, 0, + : 1, 0, 2, 0, 1, 0, 0, 0, 0, + : -2, 0, 2, 0, 1, 0, 0, 0, 0, + : 0, 0, 0, 2, 0, 0, 0, 0, 0, + : 0, 0, 2, 2, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=21,30 ) / + : 2, 0, 0, -2, 0, 0, 0, 0, 0, + : 2, 0, 2, 0, 2, 0, 0, 0, 0, + : 1, 0, 2, -2, 2, 0, 0, 0, 0, + : -1, 0, 2, 0, 1, 0, 0, 0, 0, + : 2, 0, 0, 0, 0, 0, 0, 0, 0, + : 0, 0, 2, 0, 0, 0, 0, 0, 0, + : 0, 1, 0, 0, 1, 0, 0, 0, 0, + : -1, 0, 0, 2, 1, 0, 0, 0, 0, + : 0, 2, 2, -2, 2, 0, 0, 0, 0, + : 0, 0, 2, -2, 0, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=31,40 ) / + : -1, 0, 0, 2, -1, 0, 0, 0, 0, + : 0, 1, 0, 0, -1, 0, 0, 0, 0, + : 0, 2, 0, 0, 0, 0, 0, 0, 0, + : -1, 0, 2, 2, 1, 0, 0, 0, 0, + : 1, 0, 2, 2, 2, 0, 0, 0, 0, + : 0, 1, 2, 0, 2, 0, 0, 0, 0, + : -2, 0, 2, 0, 0, 0, 0, 0, 0, + : 0, 0, 2, 2, 1, 0, 0, 0, 0, + : 0, -1, 2, 0, 2, 0, 0, 0, 0, + : 0, 0, 0, 2, 1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=41,50 ) / + : 1, 0, 2, -2, 1, 0, 0, 0, 0, + : 2, 0, 0, -2, -1, 0, 0, 0, 0, + : 2, 0, 2, -2, 2, 0, 0, 0, 0, + : 2, 0, 2, 0, 1, 0, 0, 0, 0, + : 0, 0, 0, 2, -1, 0, 0, 0, 0, + : 0, -1, 2, -2, 1, 0, 0, 0, 0, + : -1, -1, 0, 2, 0, 0, 0, 0, 0, + : 2, 0, 0, -2, 1, 0, 0, 0, 0, + : 1, 0, 0, 2, 0, 0, 0, 0, 0, + : 0, 1, 2, -2, 1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=51,60 ) / + : 1, -1, 0, 0, 0, 0, 0, 0, 0, + : -2, 0, 2, 0, 2, 0, 0, 0, 0, + : 0, -1, 0, 2, 0, 0, 0, 0, 0, + : 3, 0, 2, 0, 2, 0, 0, 0, 0, + : 0, 0, 0, 1, 0, 0, 0, 0, 0, + : 1, -1, 2, 0, 2, 0, 0, 0, 0, + : 1, 0, 0, -1, 0, 0, 0, 0, 0, + : -1, -1, 2, 2, 2, 0, 0, 0, 0, + : -1, 0, 2, 0, 0, 0, 0, 0, 0, + : 2, 0, 0, 0, -1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=61,70 ) / + : 0, -1, 2, 2, 2, 0, 0, 0, 0, + : 1, 1, 2, 0, 2, 0, 0, 0, 0, + : 2, 0, 0, 0, 1, 0, 0, 0, 0, + : 1, 1, 0, 0, 0, 0, 0, 0, 0, + : 1, 0, -2, 2, -1, 0, 0, 0, 0, + : 1, 0, 2, 0, 0, 0, 0, 0, 0, + : -1, 1, 0, 1, 0, 0, 0, 0, 0, + : 1, 0, 0, 0, 2, 0, 0, 0, 0, + : -1, 0, 1, 0, 1, 0, 0, 0, 0, + : 0, 0, 2, 1, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=71,80 ) / + : -1, 1, 0, 1, 1, 0, 0, 0, 0, + : -1, 0, 2, 4, 2, 0, 0, 0, 0, + : 0, -2, 2, -2, 1, 0, 0, 0, 0, + : 1, 0, 2, 2, 1, 0, 0, 0, 0, + : 1, 0, 0, 0, -2, 0, 0, 0, 0, + : -2, 0, 2, 2, 2, 0, 0, 0, 0, + : 1, 1, 2, -2, 2, 0, 0, 0, 0, + : -2, 0, 2, 4, 2, 0, 0, 0, 0, + : -1, 0, 4, 0, 2, 0, 0, 0, 0, + : 2, 0, 2, -2, 1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=81,90 ) / + : 1, 0, 0, -1, -1, 0, 0, 0, 0, + : 2, 0, 2, 2, 2, 0, 0, 0, 0, + : 1, 0, 0, 2, 1, 0, 0, 0, 0, + : 3, 0, 0, 0, 0, 0, 0, 0, 0, + : 0, 0, 2, -2, -1, 0, 0, 0, 0, + : 3, 0, 2, -2, 2, 0, 0, 0, 0, + : 0, 0, 4, -2, 2, 0, 0, 0, 0, + : -1, 0, 0, 4, 0, 0, 0, 0, 0, + : 0, 1, 2, 0, 1, 0, 0, 0, 0, + : 0, 0, 2, -2, 3, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=91,100 ) / + : -2, 0, 0, 4, 0, 0, 0, 0, 0, + : -1, -1, 0, 2, 1, 0, 0, 0, 0, + : -2, 0, 2, 0, -1, 0, 0, 0, 0, + : 0, 0, 2, 0, -1, 0, 0, 0, 0, + : 0, -1, 2, 0, 1, 0, 0, 0, 0, + : 0, 1, 0, 0, 2, 0, 0, 0, 0, + : 0, 0, 2, -1, 2, 0, 0, 0, 0, + : 2, 1, 0, -2, 0, 0, 0, 0, 0, + : 0, 0, 2, 4, 2, 0, 0, 0, 0, + : -1, -1, 0, 2, -1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=101,110 ) / + : -1, 1, 0, 2, 0, 0, 0, 0, 0, + : 1, -1, 0, 0, 1, 0, 0, 0, 0, + : 0, -1, 2, -2, 0, 0, 0, 0, 0, + : 0, 1, 0, 0, -2, 0, 0, 0, 0, + : 1, -1, 2, 2, 2, 0, 0, 0, 0, + : 1, 0, 0, 2, -1, 0, 0, 0, 0, + : -1, 1, 2, 2, 2, 0, 0, 0, 0, + : 3, 0, 2, 0, 1, 0, 0, 0, 0, + : 0, 1, 2, 2, 2, 0, 0, 0, 0, + : 1, 0, 2, -2, 0, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=111,120 ) / + : -1, 0, -2, 4, -1, 0, 0, 0, 0, + : -1, -1, 2, 2, 1, 0, 0, 0, 0, + : 0, -1, 2, 2, 1, 0, 0, 0, 0, + : 2, -1, 2, 0, 2, 0, 0, 0, 0, + : 0, 0, 0, 2, 2, 0, 0, 0, 0, + : 1, -1, 2, 0, 1, 0, 0, 0, 0, + : -1, 1, 2, 0, 2, 0, 0, 0, 0, + : 0, 1, 0, 2, 0, 0, 0, 0, 0, + : 0, 1, 2, -2, 0, 0, 0, 0, 0, + : 0, 3, 2, -2, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=121,130 ) / + : 0, 0, 0, 1, 1, 0, 0, 0, 0, + : -1, 0, 2, 2, 0, 0, 0, 0, 0, + : 2, 1, 2, 0, 2, 0, 0, 0, 0, + : 1, 1, 0, 0, 1, 0, 0, 0, 0, + : 2, 0, 0, 2, 0, 0, 0, 0, 0, + : 1, 1, 2, 0, 1, 0, 0, 0, 0, + : -1, 0, 0, 2, 2, 0, 0, 0, 0, + : 1, 0, -2, 2, 0, 0, 0, 0, 0, + : 0, -1, 0, 2, -1, 0, 0, 0, 0, + : -1, 0, 1, 0, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=131,140 ) / + : 0, 1, 0, 1, 0, 0, 0, 0, 0, + : 1, 0, -2, 2, -2, 0, 0, 0, 0, + : 0, 0, 0, 1, -1, 0, 0, 0, 0, + : 1, -1, 0, 0, -1, 0, 0, 0, 0, + : 0, 0, 0, 4, 0, 0, 0, 0, 0, + : 1, -1, 0, 2, 0, 0, 0, 0, 0, + : 1, 0, 2, 1, 2, 0, 0, 0, 0, + : 1, 0, 2, -1, 2, 0, 0, 0, 0, + : -1, 0, 0, 2, -2, 0, 0, 0, 0, + : 0, 0, 2, 1, 1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=141,150 ) / + : -1, 0, 2, 0, -1, 0, 0, 0, 0, + : -1, 0, 2, 4, 1, 0, 0, 0, 0, + : 0, 0, 2, 2, 0, 0, 0, 0, 0, + : 1, 1, 2, -2, 1, 0, 0, 0, 0, + : 0, 0, 1, 0, 1, 0, 0, 0, 0, + : -1, 0, 2, -1, 1, 0, 0, 0, 0, + : -2, 0, 2, 2, 1, 0, 0, 0, 0, + : 2, -1, 0, 0, 0, 0, 0, 0, 0, + : 4, 0, 2, 0, 2, 0, 0, 0, 0, + : 2, 1, 2, -2, 2, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=151,160 ) / + : 0, 1, 2, 1, 2, 0, 0, 0, 0, + : 1, 0, 4, -2, 2, 0, 0, 0, 0, + : 1, 1, 0, 0, -1, 0, 0, 0, 0, + : -2, 0, 2, 4, 1, 0, 0, 0, 0, + : 2, 0, 2, 0, 0, 0, 0, 0, 0, + : -1, 0, 1, 0, 0, 0, 0, 0, 0, + : 1, 0, 0, 1, 0, 0, 0, 0, 0, + : 0, 1, 0, 2, 1, 0, 0, 0, 0, + : -1, 0, 4, 0, 1, 0, 0, 0, 0, + : -1, 0, 0, 4, 1, 0, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=161,170 ) / + : 2, 0, 2, 2, 1, 0, 0, 0, 0, + : 2, 1, 0, 0, 0, 0, 0, 0, 0, + : 0, 0, 5, -5, 5, -3, 0, 0, 0, + : 0, 0, 0, 0, 0, 0, 0, 2, 0, + : 0, 0, 1, -1, 1, 0, 0, -1, 0, + : 0, 0, -1, 1, -1, 1, 0, 0, 0, + : 0, 0, -1, 1, 0, 0, 2, 0, 0, + : 0, 0, 3, -3, 3, 0, 0, -1, 0, + : 0, 0, -8, 8, -7, 5, 0, 0, 0, + : 0, 0, -1, 1, -1, 0, 2, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=171,180 ) / + : 0, 0, -2, 2, -2, 2, 0, 0, 0, + : 0, 0, -6, 6, -6, 4, 0, 0, 0, + : 0, 0, -2, 2, -2, 0, 8, -3, 0, + : 0, 0, 6, -6, 6, 0, -8, 3, 0, + : 0, 0, 4, -4, 4, -2, 0, 0, 0, + : 0, 0, -3, 3, -3, 2, 0, 0, 0, + : 0, 0, 4, -4, 3, 0, -8, 3, 0, + : 0, 0, -4, 4, -5, 0, 8, -3, 0, + : 0, 0, 0, 0, 0, 2, 0, 0, 0, + : 0, 0, -4, 4, -4, 3, 0, 0, 0 / + DATA ( ( NA(I,J), I=1,9 ), J=181,190 ) / + : 0, 1, -1, 1, -1, 0, 0, 1, 0, + : 0, 0, 0, 0, 0, 0, 0, 1, 0, + : 0, 0, 1, -1, 1, 1, 0, 0, 0, + : 0, 0, 2, -2, 2, 0, -2, 0, 0, + : 0, -1, -7, 7, -7, 5, 0, 0, 0, + : -2, 0, 2, 0, 2, 0, 0, -2, 0, + : -2, 0, 2, 0, 1, 0, 0, -3, 0, + : 0, 0, 2, -2, 2, 0, 0, -2, 0, + : 0, 0, 1, -1, 1, 0, 0, 1, 0, + : 0, 0, 0, 0, 0, 0, 0, 0, 2 / + DATA ( ( NA(I,J), I=1,9 ), J=191,NTERMS ) / + : 0, 0, 0, 0, 0, 0, 0, 0, 1, + : 2, 0, -2, 0, -2, 0, 0, 3, 0, + : 0, 0, 1, -1, 1, 0, 0, -2, 0, + : 0, 0, -7, 7, -7, 5, 0, 0, 0 / + +* Nutation series: longitude + DATA ( ( PSI(I,J), I=1,4 ), J=1,10 ) / + : 3341.5D0, 17206241.8D0, 3.1D0, 17409.5D0, + : -1716.8D0, -1317185.3D0, 1.4D0, -156.8D0, + : 285.7D0, -227667.0D0, 0.3D0, -23.5D0, + : -68.6D0, -207448.0D0, 0.0D0, -21.4D0, + : 950.3D0, 147607.9D0, -2.3D0, -355.0D0, + : -66.7D0, -51689.1D0, 0.2D0, 122.6D0, + : -108.6D0, 71117.6D0, 0.0D0, 7.0D0, + : 35.6D0, -38740.2D0, 0.1D0, -36.2D0, + : 85.4D0, -30127.6D0, 0.0D0, -3.1D0, + : 9.0D0, 21583.0D0, 0.1D0, -50.3D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=11,20 ) / + : 22.1D0, 12822.8D0, 0.0D0, 13.3D0, + : 3.4D0, 12350.8D0, 0.0D0, 1.3D0, + : -21.1D0, 15699.4D0, 0.0D0, 1.6D0, + : 4.2D0, 6313.8D0, 0.0D0, 6.2D0, + : -22.8D0, 5796.9D0, 0.0D0, 6.1D0, + : 15.7D0, -5961.1D0, 0.0D0, -0.6D0, + : 13.1D0, -5159.1D0, 0.0D0, -4.6D0, + : 1.8D0, 4592.7D0, 0.0D0, 4.5D0, + : -17.5D0, 6336.0D0, 0.0D0, 0.7D0, + : 16.3D0, -3851.1D0, 0.0D0, -0.4D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=21,30 ) / + : -2.8D0, 4771.7D0, 0.0D0, 0.5D0, + : 13.8D0, -3099.3D0, 0.0D0, -0.3D0, + : 0.2D0, 2860.3D0, 0.0D0, 0.3D0, + : 1.4D0, 2045.3D0, 0.0D0, 2.0D0, + : -8.6D0, 2922.6D0, 0.0D0, 0.3D0, + : -7.7D0, 2587.9D0, 0.0D0, 0.2D0, + : 8.8D0, -1408.1D0, 0.0D0, 3.7D0, + : 1.4D0, 1517.5D0, 0.0D0, 1.5D0, + : -1.9D0, -1579.7D0, 0.0D0, 7.7D0, + : 1.3D0, -2178.6D0, 0.0D0, -0.2D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=31,40 ) / + : -4.8D0, 1286.8D0, 0.0D0, 1.3D0, + : 6.3D0, 1267.2D0, 0.0D0, -4.0D0, + : -1.0D0, 1669.3D0, 0.0D0, -8.3D0, + : 2.4D0, -1020.0D0, 0.0D0, -0.9D0, + : 4.5D0, -766.9D0, 0.0D0, 0.0D0, + : -1.1D0, 756.5D0, 0.0D0, -1.7D0, + : -1.4D0, -1097.3D0, 0.0D0, -0.5D0, + : 2.6D0, -663.0D0, 0.0D0, -0.6D0, + : 0.8D0, -714.1D0, 0.0D0, 1.6D0, + : 0.4D0, -629.9D0, 0.0D0, -0.6D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=41,50 ) / + : 0.3D0, 580.4D0, 0.0D0, 0.6D0, + : -1.6D0, 577.3D0, 0.0D0, 0.5D0, + : -0.9D0, 644.4D0, 0.0D0, 0.0D0, + : 2.2D0, -534.0D0, 0.0D0, -0.5D0, + : -2.5D0, 493.3D0, 0.0D0, 0.5D0, + : -0.1D0, -477.3D0, 0.0D0, -2.4D0, + : -0.9D0, 735.0D0, 0.0D0, -1.7D0, + : 0.7D0, 406.2D0, 0.0D0, 0.4D0, + : -2.8D0, 656.9D0, 0.0D0, 0.0D0, + : 0.6D0, 358.0D0, 0.0D0, 2.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=51,60 ) / + : -0.7D0, 472.5D0, 0.0D0, -1.1D0, + : -0.1D0, -300.5D0, 0.0D0, 0.0D0, + : -1.2D0, 435.1D0, 0.0D0, -1.0D0, + : 1.8D0, -289.4D0, 0.0D0, 0.0D0, + : 0.6D0, -422.6D0, 0.0D0, 0.0D0, + : 0.8D0, -287.6D0, 0.0D0, 0.6D0, + : -38.6D0, -392.3D0, 0.0D0, 0.0D0, + : 0.7D0, -281.8D0, 0.0D0, 0.6D0, + : 0.6D0, -405.7D0, 0.0D0, 0.0D0, + : -1.2D0, 229.0D0, 0.0D0, 0.2D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=61,70 ) / + : 1.1D0, -264.3D0, 0.0D0, 0.5D0, + : -0.7D0, 247.9D0, 0.0D0, -0.5D0, + : -0.2D0, 218.0D0, 0.0D0, 0.2D0, + : 0.6D0, -339.0D0, 0.0D0, 0.8D0, + : -0.7D0, 198.7D0, 0.0D0, 0.2D0, + : -1.5D0, 334.0D0, 0.0D0, 0.0D0, + : 0.1D0, 334.0D0, 0.0D0, 0.0D0, + : -0.1D0, -198.1D0, 0.0D0, 0.0D0, + : -106.6D0, 0.0D0, 0.0D0, 0.0D0, + : -0.5D0, 165.8D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=71,80 ) / + : 0.0D0, 134.8D0, 0.0D0, 0.0D0, + : 0.9D0, -151.6D0, 0.0D0, 0.0D0, + : 0.0D0, -129.7D0, 0.0D0, 0.0D0, + : 0.8D0, -132.8D0, 0.0D0, -0.1D0, + : 0.5D0, -140.7D0, 0.0D0, 0.0D0, + : -0.1D0, 138.4D0, 0.0D0, 0.0D0, + : 0.0D0, 129.0D0, 0.0D0, -0.3D0, + : 0.5D0, -121.2D0, 0.0D0, 0.0D0, + : -0.3D0, 114.5D0, 0.0D0, 0.0D0, + : -0.1D0, 101.8D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=81,90 ) / + : -3.6D0, -101.9D0, 0.0D0, 0.0D0, + : 0.8D0, -109.4D0, 0.0D0, 0.0D0, + : 0.2D0, -97.0D0, 0.0D0, 0.0D0, + : -0.7D0, 157.3D0, 0.0D0, 0.0D0, + : 0.2D0, -83.3D0, 0.0D0, 0.0D0, + : -0.3D0, 93.3D0, 0.0D0, 0.0D0, + : -0.1D0, 92.1D0, 0.0D0, 0.0D0, + : -0.5D0, 133.6D0, 0.0D0, 0.0D0, + : -0.1D0, 81.5D0, 0.0D0, 0.0D0, + : 0.0D0, 123.9D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=91,100 ) / + : -0.3D0, 128.1D0, 0.0D0, 0.0D0, + : 0.1D0, 74.1D0, 0.0D0, -0.3D0, + : -0.2D0, -70.3D0, 0.0D0, 0.0D0, + : -0.4D0, 66.6D0, 0.0D0, 0.0D0, + : 0.1D0, -66.7D0, 0.0D0, 0.0D0, + : -0.7D0, 69.3D0, 0.0D0, -0.3D0, + : 0.0D0, -70.4D0, 0.0D0, 0.0D0, + : -0.1D0, 101.5D0, 0.0D0, 0.0D0, + : 0.5D0, -69.1D0, 0.0D0, 0.0D0, + : -0.2D0, 58.5D0, 0.0D0, 0.2D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=101,110 ) / + : 0.1D0, -94.9D0, 0.0D0, 0.2D0, + : 0.0D0, 52.9D0, 0.0D0, -0.2D0, + : 0.1D0, 86.7D0, 0.0D0, -0.2D0, + : -0.1D0, -59.2D0, 0.0D0, 0.2D0, + : 0.3D0, -58.8D0, 0.0D0, 0.1D0, + : -0.3D0, 49.0D0, 0.0D0, 0.0D0, + : -0.2D0, 56.9D0, 0.0D0, -0.1D0, + : 0.3D0, -50.2D0, 0.0D0, 0.0D0, + : -0.2D0, 53.4D0, 0.0D0, -0.1D0, + : 0.1D0, -76.5D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=111,120 ) / + : -0.2D0, 45.3D0, 0.0D0, 0.0D0, + : 0.1D0, -46.8D0, 0.0D0, 0.0D0, + : 0.2D0, -44.6D0, 0.0D0, 0.0D0, + : 0.2D0, -48.7D0, 0.0D0, 0.0D0, + : 0.1D0, -46.8D0, 0.0D0, 0.0D0, + : 0.1D0, -42.0D0, 0.0D0, 0.0D0, + : 0.0D0, 46.4D0, 0.0D0, -0.1D0, + : 0.2D0, -67.3D0, 0.0D0, 0.1D0, + : 0.0D0, -65.8D0, 0.0D0, 0.2D0, + : -0.1D0, -43.9D0, 0.0D0, 0.3D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=121,130 ) / + : 0.0D0, -38.9D0, 0.0D0, 0.0D0, + : -0.3D0, 63.9D0, 0.0D0, 0.0D0, + : -0.2D0, 41.2D0, 0.0D0, 0.0D0, + : 0.0D0, -36.1D0, 0.0D0, 0.2D0, + : -0.3D0, 58.5D0, 0.0D0, 0.0D0, + : -0.1D0, 36.1D0, 0.0D0, 0.0D0, + : 0.0D0, -39.7D0, 0.0D0, 0.0D0, + : 0.1D0, -57.7D0, 0.0D0, 0.0D0, + : -0.2D0, 33.4D0, 0.0D0, 0.0D0, + : 36.4D0, 0.0D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=131,140 ) / + : -0.1D0, 55.7D0, 0.0D0, -0.1D0, + : 0.1D0, -35.4D0, 0.0D0, 0.0D0, + : 0.1D0, -31.0D0, 0.0D0, 0.0D0, + : -0.1D0, 30.1D0, 0.0D0, 0.0D0, + : -0.3D0, 49.2D0, 0.0D0, 0.0D0, + : -0.2D0, 49.1D0, 0.0D0, 0.0D0, + : -0.1D0, 33.6D0, 0.0D0, 0.0D0, + : 0.1D0, -33.5D0, 0.0D0, 0.0D0, + : 0.1D0, -31.0D0, 0.0D0, 0.0D0, + : -0.1D0, 28.0D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=141,150 ) / + : 0.1D0, -25.2D0, 0.0D0, 0.0D0, + : 0.1D0, -26.2D0, 0.0D0, 0.0D0, + : -0.2D0, 41.5D0, 0.0D0, 0.0D0, + : 0.0D0, 24.5D0, 0.0D0, 0.1D0, + : -16.2D0, 0.0D0, 0.0D0, 0.0D0, + : 0.0D0, -22.3D0, 0.0D0, 0.0D0, + : 0.0D0, 23.1D0, 0.0D0, 0.0D0, + : -0.1D0, 37.5D0, 0.0D0, 0.0D0, + : 0.2D0, -25.7D0, 0.0D0, 0.0D0, + : 0.0D0, 25.2D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=151,160 ) / + : 0.1D0, -24.5D0, 0.0D0, 0.0D0, + : -0.1D0, 24.3D0, 0.0D0, 0.0D0, + : 0.1D0, -20.7D0, 0.0D0, 0.0D0, + : 0.1D0, -20.8D0, 0.0D0, 0.0D0, + : -0.2D0, 33.4D0, 0.0D0, 0.0D0, + : 32.9D0, 0.0D0, 0.0D0, 0.0D0, + : 0.1D0, -32.6D0, 0.0D0, 0.0D0, + : 0.0D0, 19.9D0, 0.0D0, 0.0D0, + : -0.1D0, 19.6D0, 0.0D0, 0.0D0, + : 0.0D0, -18.7D0, 0.0D0, 0.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=161,170 ) / + : 0.1D0, -19.0D0, 0.0D0, 0.0D0, + : 0.1D0, -28.6D0, 0.0D0, 0.0D0, + : 4.0D0, 178.8D0,-11.8D0, 0.3D0, + : 39.8D0, -107.3D0, -5.6D0, -1.0D0, + : 9.9D0, 164.0D0, -4.1D0, 0.1D0, + : -4.8D0, -135.3D0, -3.4D0, -0.1D0, + : 50.5D0, 75.0D0, 1.4D0, -1.2D0, + : -1.1D0, -53.5D0, 1.3D0, 0.0D0, + : -45.0D0, -2.4D0, -0.4D0, 6.6D0, + : -11.5D0, -61.0D0, -0.9D0, 0.4D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=171,180 ) / + : 4.4D0, -68.4D0, -3.4D0, 0.0D0, + : 7.7D0, -47.1D0, -4.7D0, -1.0D0, + : -42.9D0, -12.6D0, -1.2D0, 4.2D0, + : -42.8D0, 12.7D0, -1.2D0, -4.2D0, + : -7.6D0, -44.1D0, 2.1D0, -0.5D0, + : -64.1D0, 1.7D0, 0.2D0, 4.5D0, + : 36.4D0, -10.4D0, 1.0D0, 3.5D0, + : 35.6D0, 10.2D0, 1.0D0, -3.5D0, + : -1.7D0, 39.5D0, 2.0D0, 0.0D0, + : 50.9D0, -8.2D0, -0.8D0, -5.0D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=181,190 ) / + : 0.0D0, 52.3D0, 1.2D0, 0.0D0, + : -42.9D0, -17.8D0, 0.4D0, 0.0D0, + : 2.6D0, 34.3D0, 0.8D0, 0.0D0, + : -0.8D0, -48.6D0, 2.4D0, -0.1D0, + : -4.9D0, 30.5D0, 3.7D0, 0.7D0, + : 0.0D0, -43.6D0, 2.1D0, 0.0D0, + : 0.0D0, -25.4D0, 1.2D0, 0.0D0, + : 2.0D0, 40.9D0, -2.0D0, 0.0D0, + : -2.1D0, 26.1D0, 0.6D0, 0.0D0, + : 22.6D0, -3.2D0, -0.5D0, -0.5D0 / + DATA ( ( PSI(I,J), I=1,4 ), J=191,NTERMS ) / + : -7.6D0, 24.9D0, -0.4D0, -0.2D0, + : -6.2D0, 34.9D0, 1.7D0, 0.3D0, + : 2.0D0, 17.4D0, -0.4D0, 0.1D0, + : -3.9D0, 20.5D0, 2.4D0, 0.6D0 / + +* Nutation series: obliquity + DATA ( ( EPS(I,J), I=1,4 ), J=1,10 ) / + : 9205365.8D0, -1506.2D0, 885.7D0, -0.2D0, + : 573095.9D0, -570.2D0, -305.0D0, -0.3D0, + : 97845.5D0, 147.8D0, -48.8D0, -0.2D0, + : -89753.6D0, 28.0D0, 46.9D0, 0.0D0, + : 7406.7D0, -327.1D0, -18.2D0, 0.8D0, + : 22442.3D0, -22.3D0, -67.6D0, 0.0D0, + : -683.6D0, 46.8D0, 0.0D0, 0.0D0, + : 20070.7D0, 36.0D0, 1.6D0, 0.0D0, + : 12893.8D0, 39.5D0, -6.2D0, 0.0D0, + : -9593.2D0, 14.4D0, 30.2D0, -0.1D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=11,20 ) / + : -6899.5D0, 4.8D0, -0.6D0, 0.0D0, + : -5332.5D0, -0.1D0, 2.7D0, 0.0D0, + : -125.2D0, 10.5D0, 0.0D0, 0.0D0, + : -3323.4D0, -0.9D0, -0.3D0, 0.0D0, + : 3142.3D0, 8.9D0, 0.3D0, 0.0D0, + : 2552.5D0, 7.3D0, -1.2D0, 0.0D0, + : 2634.4D0, 8.8D0, 0.2D0, 0.0D0, + : -2424.4D0, 1.6D0, -0.4D0, 0.0D0, + : -123.3D0, 3.9D0, 0.0D0, 0.0D0, + : 1642.4D0, 7.3D0, -0.8D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=21,30 ) / + : 47.9D0, 3.2D0, 0.0D0, 0.0D0, + : 1321.2D0, 6.2D0, -0.6D0, 0.0D0, + : -1234.1D0, -0.3D0, 0.6D0, 0.0D0, + : -1076.5D0, -0.3D0, 0.0D0, 0.0D0, + : -61.6D0, 1.8D0, 0.0D0, 0.0D0, + : -55.4D0, 1.6D0, 0.0D0, 0.0D0, + : 856.9D0, -4.9D0, -2.1D0, 0.0D0, + : -800.7D0, -0.1D0, 0.0D0, 0.0D0, + : 685.1D0, -0.6D0, -3.8D0, 0.0D0, + : -16.9D0, -1.5D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=31,40 ) / + : 695.7D0, 1.8D0, 0.0D0, 0.0D0, + : 642.2D0, -2.6D0, -1.6D0, 0.0D0, + : 13.3D0, 1.1D0, -0.1D0, 0.0D0, + : 521.9D0, 1.6D0, 0.0D0, 0.0D0, + : 325.8D0, 2.0D0, -0.1D0, 0.0D0, + : -325.1D0, -0.5D0, 0.9D0, 0.0D0, + : 10.1D0, 0.3D0, 0.0D0, 0.0D0, + : 334.5D0, 1.6D0, 0.0D0, 0.0D0, + : 307.1D0, 0.4D0, -0.9D0, 0.0D0, + : 327.2D0, 0.5D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=41,50 ) / + : -304.6D0, -0.1D0, 0.0D0, 0.0D0, + : 304.0D0, 0.6D0, 0.0D0, 0.0D0, + : -276.8D0, -0.5D0, 0.1D0, 0.0D0, + : 268.9D0, 1.3D0, 0.0D0, 0.0D0, + : 271.8D0, 1.1D0, 0.0D0, 0.0D0, + : 271.5D0, -0.4D0, -0.8D0, 0.0D0, + : -5.2D0, 0.5D0, 0.0D0, 0.0D0, + : -220.5D0, 0.1D0, 0.0D0, 0.0D0, + : -20.1D0, 0.3D0, 0.0D0, 0.0D0, + : -191.0D0, 0.1D0, 0.5D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=51,60 ) / + : -4.1D0, 0.3D0, 0.0D0, 0.0D0, + : 130.6D0, -0.1D0, 0.0D0, 0.0D0, + : 3.0D0, 0.3D0, 0.0D0, 0.0D0, + : 122.9D0, 0.8D0, 0.0D0, 0.0D0, + : 3.7D0, -0.3D0, 0.0D0, 0.0D0, + : 123.1D0, 0.4D0, -0.3D0, 0.0D0, + : -52.7D0, 15.3D0, 0.0D0, 0.0D0, + : 120.7D0, 0.3D0, -0.3D0, 0.0D0, + : 4.0D0, -0.3D0, 0.0D0, 0.0D0, + : 126.5D0, 0.5D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=61,70 ) / + : 112.7D0, 0.5D0, -0.3D0, 0.0D0, + : -106.1D0, -0.3D0, 0.3D0, 0.0D0, + : -112.9D0, -0.2D0, 0.0D0, 0.0D0, + : 3.6D0, -0.2D0, 0.0D0, 0.0D0, + : 107.4D0, 0.3D0, 0.0D0, 0.0D0, + : -10.9D0, 0.2D0, 0.0D0, 0.0D0, + : -0.9D0, 0.0D0, 0.0D0, 0.0D0, + : 85.4D0, 0.0D0, 0.0D0, 0.0D0, + : 0.0D0, -88.8D0, 0.0D0, 0.0D0, + : -71.0D0, -0.2D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=71,80 ) / + : -70.3D0, 0.0D0, 0.0D0, 0.0D0, + : 64.5D0, 0.4D0, 0.0D0, 0.0D0, + : 69.8D0, 0.0D0, 0.0D0, 0.0D0, + : 66.1D0, 0.4D0, 0.0D0, 0.0D0, + : -61.0D0, -0.2D0, 0.0D0, 0.0D0, + : -59.5D0, -0.1D0, 0.0D0, 0.0D0, + : -55.6D0, 0.0D0, 0.2D0, 0.0D0, + : 51.7D0, 0.2D0, 0.0D0, 0.0D0, + : -49.0D0, -0.1D0, 0.0D0, 0.0D0, + : -52.7D0, -0.1D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=81,90 ) / + : -49.6D0, 1.4D0, 0.0D0, 0.0D0, + : 46.3D0, 0.4D0, 0.0D0, 0.0D0, + : 49.6D0, 0.1D0, 0.0D0, 0.0D0, + : -5.1D0, 0.1D0, 0.0D0, 0.0D0, + : -44.0D0, -0.1D0, 0.0D0, 0.0D0, + : -39.9D0, -0.1D0, 0.0D0, 0.0D0, + : -39.5D0, -0.1D0, 0.0D0, 0.0D0, + : -3.9D0, 0.1D0, 0.0D0, 0.0D0, + : -42.1D0, -0.1D0, 0.0D0, 0.0D0, + : -17.2D0, 0.1D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=91,100 ) / + : -2.3D0, 0.1D0, 0.0D0, 0.0D0, + : -39.2D0, 0.0D0, 0.0D0, 0.0D0, + : -38.4D0, 0.1D0, 0.0D0, 0.0D0, + : 36.8D0, 0.2D0, 0.0D0, 0.0D0, + : 34.6D0, 0.1D0, 0.0D0, 0.0D0, + : -32.7D0, 0.3D0, 0.0D0, 0.0D0, + : 30.4D0, 0.0D0, 0.0D0, 0.0D0, + : 0.4D0, 0.1D0, 0.0D0, 0.0D0, + : 29.3D0, 0.2D0, 0.0D0, 0.0D0, + : 31.6D0, 0.1D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=101,110 ) / + : 0.8D0, -0.1D0, 0.0D0, 0.0D0, + : -27.9D0, 0.0D0, 0.0D0, 0.0D0, + : 2.9D0, 0.0D0, 0.0D0, 0.0D0, + : -25.3D0, 0.0D0, 0.0D0, 0.0D0, + : 25.0D0, 0.1D0, 0.0D0, 0.0D0, + : 27.5D0, 0.1D0, 0.0D0, 0.0D0, + : -24.4D0, -0.1D0, 0.0D0, 0.0D0, + : 24.9D0, 0.2D0, 0.0D0, 0.0D0, + : -22.8D0, -0.1D0, 0.0D0, 0.0D0, + : 0.9D0, -0.1D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=111,120 ) / + : 24.4D0, 0.1D0, 0.0D0, 0.0D0, + : 23.9D0, 0.1D0, 0.0D0, 0.0D0, + : 22.5D0, 0.1D0, 0.0D0, 0.0D0, + : 20.8D0, 0.1D0, 0.0D0, 0.0D0, + : 20.1D0, 0.0D0, 0.0D0, 0.0D0, + : 21.5D0, 0.1D0, 0.0D0, 0.0D0, + : -20.0D0, 0.0D0, 0.0D0, 0.0D0, + : 1.4D0, 0.0D0, 0.0D0, 0.0D0, + : -0.2D0, -0.1D0, 0.0D0, 0.0D0, + : 19.0D0, 0.0D0, -0.1D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=121,130 ) / + : 20.5D0, 0.0D0, 0.0D0, 0.0D0, + : -2.0D0, 0.0D0, 0.0D0, 0.0D0, + : -17.6D0, -0.1D0, 0.0D0, 0.0D0, + : 19.0D0, 0.0D0, 0.0D0, 0.0D0, + : -2.4D0, 0.0D0, 0.0D0, 0.0D0, + : -18.4D0, -0.1D0, 0.0D0, 0.0D0, + : 17.1D0, 0.0D0, 0.0D0, 0.0D0, + : 0.4D0, 0.0D0, 0.0D0, 0.0D0, + : 18.4D0, 0.1D0, 0.0D0, 0.0D0, + : 0.0D0, 17.4D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=131,140 ) / + : -0.6D0, 0.0D0, 0.0D0, 0.0D0, + : -15.4D0, 0.0D0, 0.0D0, 0.0D0, + : -16.8D0, -0.1D0, 0.0D0, 0.0D0, + : 16.3D0, 0.0D0, 0.0D0, 0.0D0, + : -2.0D0, 0.0D0, 0.0D0, 0.0D0, + : -1.5D0, 0.0D0, 0.0D0, 0.0D0, + : -14.3D0, -0.1D0, 0.0D0, 0.0D0, + : 14.4D0, 0.0D0, 0.0D0, 0.0D0, + : -13.4D0, 0.0D0, 0.0D0, 0.0D0, + : -14.3D0, -0.1D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=141,150 ) / + : -13.7D0, 0.0D0, 0.0D0, 0.0D0, + : 13.1D0, 0.1D0, 0.0D0, 0.0D0, + : -1.7D0, 0.0D0, 0.0D0, 0.0D0, + : -12.8D0, 0.0D0, 0.0D0, 0.0D0, + : 0.0D0, -14.4D0, 0.0D0, 0.0D0, + : 12.4D0, 0.0D0, 0.0D0, 0.0D0, + : -12.0D0, 0.0D0, 0.0D0, 0.0D0, + : -0.8D0, 0.0D0, 0.0D0, 0.0D0, + : 10.9D0, 0.1D0, 0.0D0, 0.0D0, + : -10.8D0, 0.0D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=151,160 ) / + : 10.5D0, 0.0D0, 0.0D0, 0.0D0, + : -10.4D0, 0.0D0, 0.0D0, 0.0D0, + : -11.2D0, 0.0D0, 0.0D0, 0.0D0, + : 10.5D0, 0.1D0, 0.0D0, 0.0D0, + : -1.4D0, 0.0D0, 0.0D0, 0.0D0, + : 0.0D0, 0.1D0, 0.0D0, 0.0D0, + : 0.7D0, 0.0D0, 0.0D0, 0.0D0, + : -10.3D0, 0.0D0, 0.0D0, 0.0D0, + : -10.0D0, 0.0D0, 0.0D0, 0.0D0, + : 9.6D0, 0.0D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=161,170 ) / + : 9.4D0, 0.1D0, 0.0D0, 0.0D0, + : 0.6D0, 0.0D0, 0.0D0, 0.0D0, + : -87.7D0, 4.4D0, -0.4D0, -6.3D0, + : 46.3D0, 22.4D0, 0.5D0, -2.4D0, + : 15.6D0, -3.4D0, 0.1D0, 0.4D0, + : 5.2D0, 5.8D0, 0.2D0, -0.1D0, + : -30.1D0, 26.9D0, 0.7D0, 0.0D0, + : 23.2D0, -0.5D0, 0.0D0, 0.6D0, + : 1.0D0, 23.2D0, 3.4D0, 0.0D0, + : -12.2D0, -4.3D0, 0.0D0, 0.0D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=171,180 ) / + : -2.1D0, -3.7D0, -0.2D0, 0.1D0, + : -18.6D0, -3.8D0, -0.4D0, 1.8D0, + : 5.5D0, -18.7D0, -1.8D0, -0.5D0, + : -5.5D0, -18.7D0, 1.8D0, -0.5D0, + : 18.4D0, -3.6D0, 0.3D0, 0.9D0, + : -0.6D0, 1.3D0, 0.0D0, 0.0D0, + : -5.6D0, -19.5D0, 1.9D0, 0.0D0, + : 5.5D0, -19.1D0, -1.9D0, 0.0D0, + : -17.3D0, -0.8D0, 0.0D0, 0.9D0, + : -3.2D0, -8.3D0, -0.8D0, 0.3D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=181,190 ) / + : -0.1D0, 0.0D0, 0.0D0, 0.0D0, + : -5.4D0, 7.8D0, -0.3D0, 0.0D0, + : -14.8D0, 1.4D0, 0.0D0, 0.3D0, + : -3.8D0, 0.4D0, 0.0D0, -0.2D0, + : 12.6D0, 3.2D0, 0.5D0, -1.5D0, + : 0.1D0, 0.0D0, 0.0D0, 0.0D0, + : -13.6D0, 2.4D0, -0.1D0, 0.0D0, + : 0.9D0, 1.2D0, 0.0D0, 0.0D0, + : -11.9D0, -0.5D0, 0.0D0, 0.3D0, + : 0.4D0, 12.0D0, 0.3D0, -0.2D0 / + DATA ( ( EPS(I,J), I=1,4 ), J=191,NTERMS ) / + : 8.3D0, 6.1D0, -0.1D0, 0.1D0, + : 0.0D0, 0.0D0, 0.0D0, 0.0D0, + : 0.4D0, -10.8D0, 0.3D0, 0.0D0, + : 9.6D0, 2.2D0, 0.3D0, -1.2D0 / + + + +* Interval between fundamental epoch J2000.0 and given epoch (JC). + T = (DATE-DJM0)/DJC + +* Mean anomaly of the Moon. + EL = 134.96340251D0*DD2R+ + : MOD(T*(1717915923.2178D0+ + : T*( 31.8792D0+ + : T*( 0.051635D0+ + : T*( - 0.00024470D0)))),TURNAS)*DAS2R + +* Mean anomaly of the Sun. + ELP = 357.52910918D0*DD2R+ + : MOD(T*( 129596581.0481D0+ + : T*( - 0.5532D0+ + : T*( 0.000136D0+ + : T*( - 0.00001149D0)))),TURNAS)*DAS2R + +* Mean argument of the latitude of the Moon. + F = 93.27209062D0*DD2R+ + : MOD(T*(1739527262.8478D0+ + : T*( - 12.7512D0+ + : T*( - 0.001037D0+ + : T*( 0.00000417D0)))),TURNAS)*DAS2R + +* Mean elongation of the Moon from the Sun. + D = 297.85019547D0*DD2R+ + : MOD(T*(1602961601.2090D0+ + : T*( - 6.3706D0+ + : T*( 0.006539D0+ + : T*( - 0.00003169D0)))),TURNAS)*DAS2R + +* Mean longitude of the ascending node of the Moon. + OM = 125.04455501D0*DD2R+ + : MOD(T*( - 6962890.5431D0+ + : T*( 7.4722D0+ + : T*( 0.007702D0+ + : T*( - 0.00005939D0)))),TURNAS)*DAS2R + +* Mean longitude of Venus. + VE = 181.97980085D0*DD2R+MOD(210664136.433548D0*T,TURNAS)*DAS2R + +* Mean longitude of Mars. + MA = 355.43299958D0*DD2R+MOD( 68905077.493988D0*T,TURNAS)*DAS2R + +* Mean longitude of Jupiter. + JU = 34.35151874D0*DD2R+MOD( 10925660.377991D0*T,TURNAS)*DAS2R + +* Mean longitude of Saturn. + SA = 50.07744430D0*DD2R+MOD( 4399609.855732D0*T,TURNAS)*DAS2R + +* Geodesic nutation (Fukushima 1991) in microarcsec. + DP = -153.1D0*SIN(ELP)-1.9D0*SIN(2D0*ELP) + DE = 0D0 + +* Shirai & Fukushima (2001) nutation series. + DO J=NTERMS,1,-1 + THETA = DBLE(NA(1,J))*EL+ + : DBLE(NA(2,J))*ELP+ + : DBLE(NA(3,J))*F+ + : DBLE(NA(4,J))*D+ + : DBLE(NA(5,J))*OM+ + : DBLE(NA(6,J))*VE+ + : DBLE(NA(7,J))*MA+ + : DBLE(NA(8,J))*JU+ + : DBLE(NA(9,J))*SA + C = COS(THETA) + S = SIN(THETA) + DP = DP+(PSI(1,J)+PSI(3,J)*T)*C+(PSI(2,J)+PSI(4,J)*T)*S + DE = DE+(EPS(1,J)+EPS(3,J)*T)*C+(EPS(2,J)+EPS(4,J)*T)*S + END DO + +* Change of units, and addition of the precession correction. + DPSI = (DP*1D-6-0.042888D0-0.29856D0*T)*DAS2R + DEPS = (DE*1D-6-0.005171D0-0.02408D0*T)*DAS2R + +* Mean obliquity of date (Simon et al. 1994). + EPS0 = (84381.412D0+ + : (-46.80927D0+ + : (-0.000152D0+ + : (0.0019989D0+ + : (-0.00000051D0+ + : (-0.000000025D0)*T)*T)*T)*T)*T)*DAS2R + + END + SUBROUTINE sla_NUT (DATE, RMATN) +*+ +* - - - - +* N U T +* - - - - +* +* Form the matrix of nutation for a given date - Shirai & Fukushima +* 2001 theory (double precision) +* +* Reference: +* Shirai, T. & Fukushima, T., Astron.J. 121, 3270-3283 (2001). +* +* Given: +* DATE d TDB (loosely ET) as Modified Julian Date +* (=JD-2400000.5) +* Returned: +* RMATN d(3,3) nutation matrix +* +* Notes: +* +* 1 The matrix is in the sense v(true) = rmatn * v(mean) . +* where v(true) is the star vector relative to the true equator and +* equinox of date and v(mean) is the star vector relative to the +* mean equator and equinox of date. +* +* 2 The matrix represents forced nutation (but not free core +* nutation) plus corrections to the IAU~1976 precession model. +* +* 3 Earth attitude predictions made by combining the present nutation +* matrix with IAU~1976 precession are accurate to 1~mas (with +* respect to the ICRS) for a few decades around 2000. +* +* 4 The distinction between the required TDB and TT is always +* negligible. Moreover, for all but the most critical applications +* UTC is adequate. +* +* Called: sla_NUTC, sla_DEULER +* +* Last revision: 1 December 2005 +* +* Copyright P.T.Wallace. All rights reserved. +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE,RMATN(3,3) + + DOUBLE PRECISION DPSI,DEPS,EPS0 + + + +* Nutation components and mean obliquity + CALL sla_NUTC(DATE,DPSI,DEPS,EPS0) + +* Rotation matrix + CALL sla_DEULER('XZX',EPS0,-DPSI,-(EPS0+DEPS),RMATN) + + END + SUBROUTINE sla_PREBN (BEP0, BEP1, RMATP) +*+ +* - - - - - - +* P R E B N +* - - - - - - +* +* Generate the matrix of precession between two epochs, +* using the old, pre-IAU1976, Bessel-Newcomb model, using +* Kinoshita's formulation (double precision) +* +* Given: +* BEP0 dp beginning Besselian epoch +* BEP1 dp ending Besselian epoch +* +* Returned: +* RMATP dp(3,3) precession matrix +* +* The matrix is in the sense V(BEP1) = RMATP * V(BEP0) +* +* Reference: +* Kinoshita, H. (1975) 'Formulas for precession', SAO Special +* Report No. 364, Smithsonian Institution Astrophysical +* Observatory, Cambridge, Massachusetts. +* +* Called: sla_DEULER +* +* P.T.Wallace Starlink 23 August 1996 +* +* Copyright (C) 1996 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION BEP0,BEP1,RMATP(3,3) + +* Arc seconds to radians + DOUBLE PRECISION AS2R + PARAMETER (AS2R=0.484813681109535994D-5) + + DOUBLE PRECISION BIGT,T,TAS2R,W,ZETA,Z,THETA + + + +* Interval between basic epoch B1850.0 and beginning epoch in TC + BIGT = (BEP0-1850D0)/100D0 + +* Interval over which precession required, in tropical centuries + T = (BEP1-BEP0)/100D0 + +* Euler angles + TAS2R = T*AS2R + W = 2303.5548D0+(1.39720D0+0.000059D0*BIGT)*BIGT + + ZETA = (W+(0.30242D0-0.000269D0*BIGT+0.017996D0*T)*T)*TAS2R + Z = (W+(1.09478D0+0.000387D0*BIGT+0.018324D0*T)*T)*TAS2R + THETA = (2005.1125D0+(-0.85294D0-0.000365D0*BIGT)*BIGT+ + : (-0.42647D0-0.000365D0*BIGT-0.041802D0*T)*T)*TAS2R + +* Rotation matrix + CALL sla_DEULER('ZYZ',-ZETA,THETA,-Z,RMATP) + + END + SUBROUTINE sla_PRECES (SYSTEM, EP0, EP1, RA, DC) +*+ +* - - - - - - - +* P R E C E S +* - - - - - - - +* +* Precession - either FK4 (Bessel-Newcomb, pre IAU 1976) or +* FK5 (Fricke, post IAU 1976) as required. +* +* Given: +* SYSTEM char precession to be applied: 'FK4' or 'FK5' +* EP0,EP1 dp starting and ending epoch +* RA,DC dp RA,Dec, mean equator & equinox of epoch EP0 +* +* Returned: +* RA,DC dp RA,Dec, mean equator & equinox of epoch EP1 +* +* Called: sla_DRANRM, sla_PREBN, sla_PREC, sla_DCS2C, +* sla_DMXV, sla_DCC2S +* +* Notes: +* +* 1) Lowercase characters in SYSTEM are acceptable. +* +* 2) The epochs are Besselian if SYSTEM='FK4' and Julian if 'FK5'. +* For example, to precess coordinates in the old system from +* equinox 1900.0 to 1950.0 the call would be: +* CALL sla_PRECES ('FK4', 1900D0, 1950D0, RA, DC) +* +* 3) This routine will NOT correctly convert between the old and +* the new systems - for example conversion from B1950 to J2000. +* For these purposes see sla_FK425, sla_FK524, sla_FK45Z and +* sla_FK54Z. +* +* 4) If an invalid SYSTEM is supplied, values of -99D0,-99D0 will +* be returned for both RA and DC. +* +* P.T.Wallace Starlink 20 April 1990 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + CHARACTER SYSTEM*(*) + DOUBLE PRECISION EP0,EP1,RA,DC + + DOUBLE PRECISION PM(3,3),V1(3),V2(3) + CHARACTER SYSUC*3 + + DOUBLE PRECISION sla_DRANRM + + + + +* Convert to uppercase and validate SYSTEM + SYSUC=SYSTEM + IF (SYSUC(1:1).EQ.'f') SYSUC(1:1)='F' + IF (SYSUC(2:2).EQ.'k') SYSUC(2:2)='K' + IF (SYSUC.NE.'FK4'.AND.SYSUC.NE.'FK5') THEN + RA=-99D0 + DC=-99D0 + ELSE + +* Generate appropriate precession matrix + IF (SYSUC.EQ.'FK4') THEN + CALL sla_PREBN(EP0,EP1,PM) + ELSE + CALL sla_PREC(EP0,EP1,PM) + END IF + +* Convert RA,Dec to x,y,z + CALL sla_DCS2C(RA,DC,V1) + +* Precess + CALL sla_DMXV(PM,V1,V2) + +* Back to RA,Dec + CALL sla_DCC2S(V2,RA,DC) + RA=sla_DRANRM(RA) + + END IF + + END + SUBROUTINE sla_PREC (EP0, EP1, RMATP) +*+ +* - - - - - +* P R E C +* - - - - - +* +* Form the matrix of precession between two epochs (IAU 1976, FK5) +* (double precision) +* +* Given: +* EP0 dp beginning epoch +* EP1 dp ending epoch +* +* Returned: +* RMATP dp(3,3) precession matrix +* +* Notes: +* +* 1) The epochs are TDB (loosely ET) Julian epochs. +* +* 2) The matrix is in the sense V(EP1) = RMATP * V(EP0) +* +* 3) Though the matrix method itself is rigorous, the precession +* angles are expressed through canonical polynomials which are +* valid only for a limited time span. There are also known +* errors in the IAU precession rate. The absolute accuracy +* of the present formulation is better than 0.1 arcsec from +* 1960AD to 2040AD, better than 1 arcsec from 1640AD to 2360AD, +* and remains below 3 arcsec for the whole of the period +* 500BC to 3000AD. The errors exceed 10 arcsec outside the +* range 1200BC to 3900AD, exceed 100 arcsec outside 4200BC to +* 5600AD and exceed 1000 arcsec outside 6800BC to 8200AD. +* The SLALIB routine sla_PRECL implements a more elaborate +* model which is suitable for problems spanning several +* thousand years. +* +* References: +* Lieske,J.H., 1979. Astron.Astrophys.,73,282. +* equations (6) & (7), p283. +* Kaplan,G.H., 1981. USNO circular no. 163, pA2. +* +* Called: sla_DEULER +* +* P.T.Wallace Starlink 23 August 1996 +* +* Copyright (C) 1996 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION EP0,EP1,RMATP(3,3) + +* Arc seconds to radians + DOUBLE PRECISION AS2R + PARAMETER (AS2R=0.484813681109535994D-5) + + DOUBLE PRECISION T0,T,TAS2R,W,ZETA,Z,THETA + + + +* Interval between basic epoch J2000.0 and beginning epoch (JC) + T0 = (EP0-2000D0)/100D0 + +* Interval over which precession required (JC) + T = (EP1-EP0)/100D0 + +* Euler angles + TAS2R = T*AS2R + W = 2306.2181D0+(1.39656D0-0.000139D0*T0)*T0 + + ZETA = (W+((0.30188D0-0.000344D0*T0)+0.017998D0*T)*T)*TAS2R + Z = (W+((1.09468D0+0.000066D0*T0)+0.018203D0*T)*T)*TAS2R + THETA = ((2004.3109D0+(-0.85330D0-0.000217D0*T0)*T0) + : +((-0.42665D0-0.000217D0*T0)-0.041833D0*T)*T)*TAS2R + +* Rotation matrix + CALL sla_DEULER('ZYZ',-ZETA,THETA,-Z,RMATP) + + END + SUBROUTINE sla_PVOBS (P, H, STL, PV) +*+ +* - - - - - - +* P V O B S +* - - - - - - +* +* Position and velocity of an observing station (double precision) +* +* Given: +* P dp latitude (geodetic, radians) +* H dp height above reference spheroid (geodetic, metres) +* STL dp local apparent sidereal time (radians) +* +* Returned: +* PV dp(6) position/velocity 6-vector (AU, AU/s, true equator +* and equinox of date) +* +* Called: sla_GEOC +* +* IAU 1976 constants are used. +* +* P.T.Wallace Starlink 14 November 1994 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +*- + + IMPLICIT NONE + + DOUBLE PRECISION P,H,STL,PV(6) + + DOUBLE PRECISION R,Z,S,C,V + +* Mean sidereal rate (at J2000) in radians per (UT1) second + DOUBLE PRECISION SR + PARAMETER (SR=7.292115855306589D-5) + + + +* Geodetic to geocentric conversion + CALL sla_GEOC(P,H,R,Z) + +* Functions of ST + S=SIN(STL) + C=COS(STL) + +* Speed + V=SR*R + +* Position + PV(1)=R*C + PV(2)=R*S + PV(3)=Z + +* Velocity + PV(4)=-V*S + PV(5)=V*C + PV(6)=0D0 + + END diff --git a/lib/tm2.f90 b/lib/tm2.f90 deleted file mode 100644 index cf1107836..000000000 --- a/lib/tm2.f90 +++ /dev/null @@ -1,14 +0,0 @@ -subroutine tm2(day,xlat4,xlon4,xl4,b4) - - implicit real*8 (a-h,o-z) - parameter (RADS=0.0174532925199433d0) - - real*4 xlat4,xlon4,xl4,b4 - - glat=xlat4*RADS - glong=xlon4*RADS - call tmoonsub(day,glat,glong,el,rv,xl,b,pax) - xl4=xl - b4=b - -end subroutine tm2 diff --git a/lib/toxyz.f90 b/lib/toxyz.f90 deleted file mode 100644 index 646dd7a2a..000000000 --- a/lib/toxyz.f90 +++ /dev/null @@ -1,25 +0,0 @@ -subroutine toxyz(alpha,delta,r,vec) - - implicit real*8 (a-h,o-z) - real*8 vec(3) - - vec(1)=r*cos(delta)*cos(alpha) - vec(2)=r*cos(delta)*sin(alpha) - vec(3)=r*sin(delta) - - return -end subroutine toxyz - -subroutine fromxyz(vec,alpha,delta,r) - - implicit real*8 (a-h,o-z) - real*8 vec(3) - data twopi/6.283185307d0/ - - r=sqrt(vec(1)**2 + vec(2)**2 + vec(3)**2) - alpha=atan2(vec(2),vec(1)) - if(alpha.lt.0.d0) alpha=alpha+twopi - delta=asin(vec(3)/r) - - return -end subroutine fromxyz