-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathephMoon.m
108 lines (91 loc) · 3.56 KB
/
ephMoon.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
function [xP, vP] = ephMoon(MJD2000)
% ephMoon.m - Ephemerides (cartesian position and velocity) of the Moon.
%
% PROTOTYPE:
% [xP, vP] = ephMoon(MJD2000)
%
% DESCRIPTION:
% It gives the position and the velocity of the Moon at a given epoch in
% Geocentric Equatorial Reference Frame (IAU-76/FK5 J2000, mean equator,
% mean equinox frame) This frame {x,y,z} is characterised by:
% x-axis: on the equatorial plane, along the direction of the gamma
% point
% z-axis: direction of the north pole
% y-axis: on the equatorial plane, completes the reference frame
%
% INPUT:
% MJD2000[1] Epoch in Modified Julian Date 2000 (MJD2000 since 12:00
% noon 01/01/2000)
%
% OUTPUT:
% xP Position vector of the Moon in cartesian coordinates,
% expressed in the Geocentric Equatorial Reference Frame [km].
% vP Velocity vector of the Moon in cartesian coordinates,
% expressed in the Geocentric Equatorial Reference Frame
% [km/s]. The velocity is computed by numerical
% differentiation on a 1 second interval.
%
% CALLED FUNCTIONS:
% (none)
%
% REFERENCES:
% Algorithm taken from "Fundamentals of Astrodynamics and Applications"
% (3rd edition), D. A. Vallado, p.290 (algorithm 31).
%
% AUTHOR:
% Daniel Novak, 04/03/2008, MATLAB, ephMoon.m
%
% PREVIOUS VERSION:
% Daniel Novak, 04/03/2008, MATLAB, eph_moon.m
% - Header and function name in accordance with guidlines.
%
% CHANGELOG:
% 05/03/2008, REVISION, Matteo Ceriotti
% 06/05/2008, Nicolas Croisard: Use of COS and SIN instead of COSD and
% SIND for speed (15 times faster)
% 04/10/2010, Camilla Colombo: Header and function name in accordance
% with guidlines.
% 25/02/2013, Camilla Colombo: Header and function name in accordance
% with guidlines, fixed line 102.
%
% -------------------------------------------------------------------------
T_TDB = MJD2000/36525;
angles = [134.9 + 477198.85*T_TDB; ... % L_ecl 1 and p 1
259.2 - 413335.38*T_TDB; ... % L_ecl 2 and p 2
235.7 + 890534.23*T_TDB; ... % L_ecl 3 and p 3
269.9 + 954397.7*T_TDB; ... % L_ecl 4 and p 4
357.5 + 35999.05*T_TDB; ... % L_ecl 5
186.6 + 966404.05*T_TDB; ... % L_ecl 6
93.3 + 483202.03*T_TDB; ... % phi_ecl 1
228.2 + 960400.87*T_TDB; ... % phi_ecl 2
318.3 + 6003.18*T_TDB; ... % phi_ecl 3
217.6 - 407332.2*T_TDB; ... % phi_ecl 4
];
% Transform angles in radians
angles = angles * pi/180;
% Sinus of the first 10 angles
s = sin(angles(1:10));
% Cosinus of the first 4 angles
c = cos(angles(1:4));
% Compute L_ecl, phi_ecl, P and eps
L_ecl = (218.32 + 481267.883*T_TDB) ...
+ [+6.29, -1.27, +0.66, +0.21, -0.19, -0.11] * s(1:6);
phi_ecl = [+5.13, +0.28, -0.28, -0.17] * s(7:10);
P = 0.9508 ...
+ [+0.0518, +0.0095, +0.0078, +0.0028] * c;
eps = 23.439291 - 0.0130042*T_TDB - 1.64e-7*T_TDB^2 + 5.04e-7*T_TDB^3;
% Transform L_ecl, phi_ecl, P and eps in radians
L_ecl = L_ecl * pi/180;
phi_ecl = phi_ecl * pi/180;
P = P * pi/180;
eps = eps * pi/180;
r = 1/sin(P) * 6378.16; % Radius of the Earth copied from astro_constants for speed
% r = 1/sin(P) * 6371.01; % from Horizon (does not change much using one of
% the other...by Camilla
xP = r * [cos(L_ecl)*cos(phi_ecl), ...
cos(eps)*cos(phi_ecl)*sin(L_ecl) - sin(eps)*sin(phi_ecl), ...
sin(eps)*cos(phi_ecl)*sin(L_ecl) + cos(eps)*sin(phi_ecl)];
if nargout > 1
vP = (ephMoon(MJD2000+1/86400) - xP);
end
return