mirror of
https://github.com/eddyem/BTA_lib.git
synced 2025-12-06 02:35:20 +03:00
380 lines
15 KiB
Fortran
380 lines
15 KiB
Fortran
SUBROUTINE sla_MOON (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 sla_CALYD or
|
|
* sla_CLYD.
|
|
*
|
|
* 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 sla_DMOON).
|
|
*
|
|
* Reference:
|
|
* Meeus, l'Astronomie, June 1984, p348.
|
|
*
|
|
* Called: sla_CS2C6
|
|
*
|
|
* 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., 59 Temple Place, Suite 330,
|
|
* Boston, MA 02111-1307 USA
|
|
*
|
|
*-
|
|
|
|
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 sla_CS2C6(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
|