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