aboutsummaryrefslogtreecommitdiff
path: root/math/slalib/moon.f
diff options
context:
space:
mode:
Diffstat (limited to 'math/slalib/moon.f')
-rw-r--r--math/slalib/moon.f380
1 files changed, 380 insertions, 0 deletions
diff --git a/math/slalib/moon.f b/math/slalib/moon.f
new file mode 100644
index 00000000..c77395ee
--- /dev/null
+++ b/math/slalib/moon.f
@@ -0,0 +1,380 @@
+ SUBROUTINE slMOON (IY, ID, FD, PV)
+*+
+* - - - - -
+* M O O N
+* - - - - -
+*
+* Approximate geocentric position and velocity of the Moon
+* (single precision).
+*
+* Given:
+* IY i year
+* ID i day in year (1 = Jan 1st)
+* FD r fraction of day
+*
+* Returned:
+* PV r(6) Moon position & velocity vector
+*
+* Notes:
+*
+* 1 The date and time is TDB (loosely ET) in a Julian calendar
+* which has been aligned to the ordinary Gregorian
+* calendar for the interval 1900 March 1 to 2100 February 28.
+* The year and day can be obtained by calling slCAYD or
+* slCLYD.
+*
+* 2 The Moon 6-vector is Moon centre relative to Earth centre,
+* mean equator and equinox of date. Position part, PV(1-3),
+* is in AU; velocity part, PV(4-6), is in AU/sec.
+*
+* 3 The position is accurate to better than 0.5 arcminute
+* in direction and 1000 km in distance. The velocity
+* is accurate to better than 0.5"/hour in direction and
+* 4 m/s in distance. (RMS figures with respect to JPL DE200
+* for the interval 1960-2025 are 14 arcsec and 0.2 arcsec/hour in
+* longitude, 9 arcsec and 0.2 arcsec/hour in latitude, 350 km and
+* 2 m/s in distance.) Note that the distance accuracy is
+* comparatively poor because this routine is principally intended
+* for computing topocentric direction.
+*
+* 4 This routine is only a partial implementation of the original
+* Meeus algorithm (reference below), which offers 4 times the
+* accuracy in direction and 30 times the accuracy in distance
+* when fully implemented (as it is in slDMON).
+*
+* Reference:
+* Meeus, l'Astronomie, June 1984, p348.
+*
+* Called: slS2C6
+*
+* P.T.Wallace Starlink 8 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
+*
+* Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
+*-
+
+ IMPLICIT NONE
+
+ INTEGER IY,ID
+ REAL FD,PV(6)
+
+ INTEGER ITP(4,4),ITL(4,39),ITB(4,29),I,IY4,N
+ REAL D2R,RATCON,ERADAU
+ REAL ELP0,ELP1,ELP1I,ELP1F
+ REAL EM0,EM1,EM1F
+ REAL EMP0,EMP1,EMP1I,EMP1F
+ REAL D0,D1,D1I,D1F
+ REAL F0,F1,F1I,F1F
+ REAL TL(39)
+ REAL TB(29)
+ REAL TP(4)
+ REAL YI,YF,T,ELP,EM,EMP,D,F,EL,ELD,COEFF,CEM,CEMP
+ REAL CD,CF,THETA,THETAD,B,BD,P,PD,SP,R,RD
+ REAL V(6),EPS,SINEPS,COSEPS
+
+* Degrees to radians
+ PARAMETER (D2R=1.745329252E-2)
+
+* Rate conversion factor: D2R**2/(86400*365.25)
+ PARAMETER (RATCON=9.652743551E-12)
+
+* Earth radius in AU: 6378.137/149597870
+ PARAMETER (ERADAU=4.2635212653763E-5)
+
+*
+* Coefficients for fundamental arguments
+*
+* Fixed term (deg), term in T (deg & whole revs + fraction per year)
+*
+* Moon's mean longitude
+ DATA ELP0,ELP1,ELP1I,ELP1F /
+ : 270.434164, 4812.678831, 4680., 132.678831 /
+*
+* Sun's mean anomaly
+ DATA EM0,EM1,EM1F /
+ : 358.475833, 359.990498, 359.990498 /
+*
+* Moon's mean anomaly
+ DATA EMP0,EMP1,EMP1I,EMP1F /
+ : 296.104608, 4771.988491, 4680., 91.988491 /
+*
+* Moon's mean elongation
+ DATA D0,D1,D1I,D1F /
+ : 350.737486, 4452.671142, 4320., 132.671142 /
+*
+* Mean distance of the Moon from its ascending node
+ DATA F0,F1,F1I,F1F /
+ : 11.250889, 4832.020251, 4680., 152.020251 /
+
+*
+* Coefficients for Moon position
+*
+* T(N) = coefficient of term (deg)
+* IT(N,1-4) = coefficients of M, M', D, F in argument
+*
+* Longitude
+* M M' D F
+ DATA TL( 1)/ +6.288750 /,
+ : (ITL(I, 1),I=1,4)/ 0, +1, 0, 0 /
+ DATA TL( 2)/ +1.274018 /,
+ : (ITL(I, 2),I=1,4)/ 0, -1, +2, 0 /
+ DATA TL( 3)/ +0.658309 /,
+ : (ITL(I, 3),I=1,4)/ 0, 0, +2, 0 /
+ DATA TL( 4)/ +0.213616 /,
+ : (ITL(I, 4),I=1,4)/ 0, +2, 0, 0 /
+ DATA TL( 5)/ -0.185596 /,
+ : (ITL(I, 5),I=1,4)/ +1, 0, 0, 0 /
+ DATA TL( 6)/ -0.114336 /,
+ : (ITL(I, 6),I=1,4)/ 0, 0, 0, +2 /
+ DATA TL( 7)/ +0.058793 /,
+ : (ITL(I, 7),I=1,4)/ 0, -2, +2, 0 /
+ DATA TL( 8)/ +0.057212 /,
+ : (ITL(I, 8),I=1,4)/ -1, -1, +2, 0 /
+ DATA TL( 9)/ +0.053320 /,
+ : (ITL(I, 9),I=1,4)/ 0, +1, +2, 0 /
+ DATA TL(10)/ +0.045874 /,
+ : (ITL(I,10),I=1,4)/ -1, 0, +2, 0 /
+ DATA TL(11)/ +0.041024 /,
+ : (ITL(I,11),I=1,4)/ -1, +1, 0, 0 /
+ DATA TL(12)/ -0.034718 /,
+ : (ITL(I,12),I=1,4)/ 0, 0, +1, 0 /
+ DATA TL(13)/ -0.030465 /,
+ : (ITL(I,13),I=1,4)/ +1, +1, 0, 0 /
+ DATA TL(14)/ +0.015326 /,
+ : (ITL(I,14),I=1,4)/ 0, 0, +2, -2 /
+ DATA TL(15)/ -0.012528 /,
+ : (ITL(I,15),I=1,4)/ 0, +1, 0, +2 /
+ DATA TL(16)/ -0.010980 /,
+ : (ITL(I,16),I=1,4)/ 0, -1, 0, +2 /
+ DATA TL(17)/ +0.010674 /,
+ : (ITL(I,17),I=1,4)/ 0, -1, +4, 0 /
+ DATA TL(18)/ +0.010034 /,
+ : (ITL(I,18),I=1,4)/ 0, +3, 0, 0 /
+ DATA TL(19)/ +0.008548 /,
+ : (ITL(I,19),I=1,4)/ 0, -2, +4, 0 /
+ DATA TL(20)/ -0.007910 /,
+ : (ITL(I,20),I=1,4)/ +1, -1, +2, 0 /
+ DATA TL(21)/ -0.006783 /,
+ : (ITL(I,21),I=1,4)/ +1, 0, +2, 0 /
+ DATA TL(22)/ +0.005162 /,
+ : (ITL(I,22),I=1,4)/ 0, +1, -1, 0 /
+ DATA TL(23)/ +0.005000 /,
+ : (ITL(I,23),I=1,4)/ +1, 0, +1, 0 /
+ DATA TL(24)/ +0.004049 /,
+ : (ITL(I,24),I=1,4)/ -1, +1, +2, 0 /
+ DATA TL(25)/ +0.003996 /,
+ : (ITL(I,25),I=1,4)/ 0, +2, +2, 0 /
+ DATA TL(26)/ +0.003862 /,
+ : (ITL(I,26),I=1,4)/ 0, 0, +4, 0 /
+ DATA TL(27)/ +0.003665 /,
+ : (ITL(I,27),I=1,4)/ 0, -3, +2, 0 /
+ DATA TL(28)/ +0.002695 /,
+ : (ITL(I,28),I=1,4)/ -1, +2, 0, 0 /
+ DATA TL(29)/ +0.002602 /,
+ : (ITL(I,29),I=1,4)/ 0, +1, -2, -2 /
+ DATA TL(30)/ +0.002396 /,
+ : (ITL(I,30),I=1,4)/ -1, -2, +2, 0 /
+ DATA TL(31)/ -0.002349 /,
+ : (ITL(I,31),I=1,4)/ 0, +1, +1, 0 /
+ DATA TL(32)/ +0.002249 /,
+ : (ITL(I,32),I=1,4)/ -2, 0, +2, 0 /
+ DATA TL(33)/ -0.002125 /,
+ : (ITL(I,33),I=1,4)/ +1, +2, 0, 0 /
+ DATA TL(34)/ -0.002079 /,
+ : (ITL(I,34),I=1,4)/ +2, 0, 0, 0 /
+ DATA TL(35)/ +0.002059 /,
+ : (ITL(I,35),I=1,4)/ -2, -1, +2, 0 /
+ DATA TL(36)/ -0.001773 /,
+ : (ITL(I,36),I=1,4)/ 0, +1, +2, -2 /
+ DATA TL(37)/ -0.001595 /,
+ : (ITL(I,37),I=1,4)/ 0, 0, +2, +2 /
+ DATA TL(38)/ +0.001220 /,
+ : (ITL(I,38),I=1,4)/ -1, -1, +4, 0 /
+ DATA TL(39)/ -0.001110 /,
+ : (ITL(I,39),I=1,4)/ 0, +2, 0, +2 /
+*
+* Latitude
+* M M' D F
+ DATA TB( 1)/ +5.128189 /,
+ : (ITB(I, 1),I=1,4)/ 0, 0, 0, +1 /
+ DATA TB( 2)/ +0.280606 /,
+ : (ITB(I, 2),I=1,4)/ 0, +1, 0, +1 /
+ DATA TB( 3)/ +0.277693 /,
+ : (ITB(I, 3),I=1,4)/ 0, +1, 0, -1 /
+ DATA TB( 4)/ +0.173238 /,
+ : (ITB(I, 4),I=1,4)/ 0, 0, +2, -1 /
+ DATA TB( 5)/ +0.055413 /,
+ : (ITB(I, 5),I=1,4)/ 0, -1, +2, +1 /
+ DATA TB( 6)/ +0.046272 /,
+ : (ITB(I, 6),I=1,4)/ 0, -1, +2, -1 /
+ DATA TB( 7)/ +0.032573 /,
+ : (ITB(I, 7),I=1,4)/ 0, 0, +2, +1 /
+ DATA TB( 8)/ +0.017198 /,
+ : (ITB(I, 8),I=1,4)/ 0, +2, 0, +1 /
+ DATA TB( 9)/ +0.009267 /,
+ : (ITB(I, 9),I=1,4)/ 0, +1, +2, -1 /
+ DATA TB(10)/ +0.008823 /,
+ : (ITB(I,10),I=1,4)/ 0, +2, 0, -1 /
+ DATA TB(11)/ +0.008247 /,
+ : (ITB(I,11),I=1,4)/ -1, 0, +2, -1 /
+ DATA TB(12)/ +0.004323 /,
+ : (ITB(I,12),I=1,4)/ 0, -2, +2, -1 /
+ DATA TB(13)/ +0.004200 /,
+ : (ITB(I,13),I=1,4)/ 0, +1, +2, +1 /
+ DATA TB(14)/ +0.003372 /,
+ : (ITB(I,14),I=1,4)/ -1, 0, -2, +1 /
+ DATA TB(15)/ +0.002472 /,
+ : (ITB(I,15),I=1,4)/ -1, -1, +2, +1 /
+ DATA TB(16)/ +0.002222 /,
+ : (ITB(I,16),I=1,4)/ -1, 0, +2, +1 /
+ DATA TB(17)/ +0.002072 /,
+ : (ITB(I,17),I=1,4)/ -1, -1, +2, -1 /
+ DATA TB(18)/ +0.001877 /,
+ : (ITB(I,18),I=1,4)/ -1, +1, 0, +1 /
+ DATA TB(19)/ +0.001828 /,
+ : (ITB(I,19),I=1,4)/ 0, -1, +4, -1 /
+ DATA TB(20)/ -0.001803 /,
+ : (ITB(I,20),I=1,4)/ +1, 0, 0, +1 /
+ DATA TB(21)/ -0.001750 /,
+ : (ITB(I,21),I=1,4)/ 0, 0, 0, +3 /
+ DATA TB(22)/ +0.001570 /,
+ : (ITB(I,22),I=1,4)/ -1, +1, 0, -1 /
+ DATA TB(23)/ -0.001487 /,
+ : (ITB(I,23),I=1,4)/ 0, 0, +1, +1 /
+ DATA TB(24)/ -0.001481 /,
+ : (ITB(I,24),I=1,4)/ +1, +1, 0, +1 /
+ DATA TB(25)/ +0.001417 /,
+ : (ITB(I,25),I=1,4)/ -1, -1, 0, +1 /
+ DATA TB(26)/ +0.001350 /,
+ : (ITB(I,26),I=1,4)/ -1, 0, 0, +1 /
+ DATA TB(27)/ +0.001330 /,
+ : (ITB(I,27),I=1,4)/ 0, 0, -1, +1 /
+ DATA TB(28)/ +0.001106 /,
+ : (ITB(I,28),I=1,4)/ 0, +3, 0, +1 /
+ DATA TB(29)/ +0.001020 /,
+ : (ITB(I,29),I=1,4)/ 0, 0, +4, -1 /
+*
+* Parallax
+* M M' D F
+ DATA TP( 1)/ +0.051818 /,
+ : (ITP(I, 1),I=1,4)/ 0, +1, 0, 0 /
+ DATA TP( 2)/ +0.009531 /,
+ : (ITP(I, 2),I=1,4)/ 0, -1, +2, 0 /
+ DATA TP( 3)/ +0.007843 /,
+ : (ITP(I, 3),I=1,4)/ 0, 0, +2, 0 /
+ DATA TP( 4)/ +0.002824 /,
+ : (ITP(I, 4),I=1,4)/ 0, +2, 0, 0 /
+
+
+
+* Whole years & fraction of year, and years since J1900.0
+ YI=FLOAT(IY-1900)
+ IY4=MOD(MOD(IY,4)+4,4)
+ YF=(FLOAT(4*(ID-1/(IY4+1))-IY4-2)+4.0*FD)/1461.0
+ T=YI+YF
+
+* Moon's mean longitude
+ ELP=D2R*MOD(ELP0+ELP1I*YF+ELP1F*T,360.0)
+
+* Sun's mean anomaly
+ EM=D2R*MOD(EM0+EM1F*T,360.0)
+
+* Moon's mean anomaly
+ EMP=D2R*MOD(EMP0+EMP1I*YF+EMP1F*T,360.0)
+
+* Moon's mean elongation
+ D=D2R*MOD(D0+D1I*YF+D1F*T,360.0)
+
+* Mean distance of the moon from its ascending node
+ F=D2R*MOD(F0+F1I*YF+F1F*T,360.0)
+
+* Longitude
+ EL=0.0
+ ELD=0.0
+ DO N=39,1,-1
+ COEFF=TL(N)
+ CEM=FLOAT(ITL(1,N))
+ CEMP=FLOAT(ITL(2,N))
+ CD=FLOAT(ITL(3,N))
+ CF=FLOAT(ITL(4,N))
+ THETA=CEM*EM+CEMP*EMP+CD*D+CF*F
+ THETAD=CEM*EM1+CEMP*EMP1+CD*D1+CF*F1
+ EL=EL+COEFF*SIN(THETA)
+ ELD=ELD+COEFF*COS(THETA)*THETAD
+ END DO
+ EL=EL*D2R+ELP
+ ELD=RATCON*(ELD+ELP1/D2R)
+
+* Latitude
+ B=0.0
+ BD=0.0
+ DO N=29,1,-1
+ COEFF=TB(N)
+ CEM=FLOAT(ITB(1,N))
+ CEMP=FLOAT(ITB(2,N))
+ CD=FLOAT(ITB(3,N))
+ CF=FLOAT(ITB(4,N))
+ THETA=CEM*EM+CEMP*EMP+CD*D+CF*F
+ THETAD=CEM*EM1+CEMP*EMP1+CD*D1+CF*F1
+ B=B+COEFF*SIN(THETA)
+ BD=BD+COEFF*COS(THETA)*THETAD
+ END DO
+ B=B*D2R
+ BD=RATCON*BD
+
+* Parallax
+ P=0.0
+ PD=0.0
+ DO N=4,1,-1
+ COEFF=TP(N)
+ CEM=FLOAT(ITP(1,N))
+ CEMP=FLOAT(ITP(2,N))
+ CD=FLOAT(ITP(3,N))
+ CF=FLOAT(ITP(4,N))
+ THETA=CEM*EM+CEMP*EMP+CD*D+CF*F
+ THETAD=CEM*EM1+CEMP*EMP1+CD*D1+CF*F1
+ P=P+COEFF*COS(THETA)
+ PD=PD-COEFF*SIN(THETA)*THETAD
+ END DO
+ P=(P+0.950724)*D2R
+ PD=RATCON*PD
+
+* Transform parallax to distance (AU, AU/sec)
+ SP=SIN(P)
+ R=ERADAU/SP
+ RD=-R*PD/SP
+
+* Longitude, latitude to x,y,z (AU)
+ CALL slS2C6(EL,B,R,ELD,BD,RD,V)
+
+* Mean obliquity
+ EPS=D2R*(23.45229-0.00013*T)
+ SINEPS=SIN(EPS)
+ COSEPS=COS(EPS)
+
+* Rotate Moon position and velocity into equatorial system
+ PV(1)=V(1)
+ PV(2)=V(2)*COSEPS-V(3)*SINEPS
+ PV(3)=V(2)*SINEPS+V(3)*COSEPS
+ PV(4)=V(4)
+ PV(5)=V(5)*COSEPS-V(6)*SINEPS
+ PV(6)=V(5)*SINEPS+V(6)*COSEPS
+
+ END
generated by cgit (git 2.34.1) at 2025-09-20 10:16:01 -0400