# 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]; % ˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜˜ ```