function coe coe_from_sv This function computes the classical orbital

 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
% ˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜
function coe = coe_from_sv(R,V)
% ˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜
%
% This function computes the classical orbital elements (coe)
% from the state vector (R,V) using Algorithm 4.1.
%
% mu - gravitational parameter (kmˆ3/sˆ2)
% R - position vector in the geocentric equatorial frame
% (km)
% V - velocity vector in the geocentric equatorial frame
% (km)
% r, v - the magnitudes of R and V
% vr - radial velocity component (km/s)
% H - the angular momentum vector (kmˆ2/s)
% h - the magnitude of H (kmˆ2/s)
% incl - inclination of the orbit (rad)
% N - the node line vector (kmˆ2/s)
% n - the magnitude of N
% cp - cross product of N and R
% RA - right ascension of the ascending node (rad)
% E - eccentricity vector
% e - eccentricity (magnitude of E)
% eps - a small number below which the eccentricity is
% considered to be zero
% w - argument of perigee (rad)
% TA - true anomaly (rad)
% a - semimajor axis (km)
% pi - 3.1415926...
% coe - vector of orbital elements [h e RA incl w TA a]
%
% User M-functions required: None
% ------------------------------------------------------------
global mu;
eps = 1.e-10;
r = norm(R);
v = norm(V);
vr = dot(R,V)/r;
H = cross(R,V);
h = norm(H);
%...Equation 4.7:
incl = acos(H(3)/h);
%...Equation 4.8:
N = cross([0 0 1],H);
n = norm(N);
%...Equation 4.9:
if n ∼ = 0
RA = acos(N(1)/n);
if N(2) < 0
RA = 2*pi - RA;
end
else
RA = 0;
end
%...Equation 4.10:
E = 1/mu*((vˆ2 - mu/r)*R - r*vr*V);
e = norm(E);
%...Equation 4.12 (incorporating the case e = 0):
if n ∼ = 0
if e > eps
w = acos(dot(N,E)/n/e);
if E(3) < 0
w = 2*pi - w;
end
else
w = 0;
end
else
w = 0;
end
%...Equation 4.13a (incorporating the case e = 0):
if e > eps
TA = acos(dot(E,R)/e/r);
if vr < 0
TA = 2*pi - TA;
end
else
cp = cross(N,R);
if cp(3) >= 0
TA = acos(dot(N,R)/n/r);
else
TA = 2*pi - acos(dot(N,R)/n/r);
end
end
%...Equation 2.61 (a < 0 for a hyperbola):
a = hˆ2/mu/(1 - eˆ2);
coe = [h e RA incl w TA a];
% ˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜