星历表-matlab

用于测量八大行星在任一时刻在测站坐标系下所处的位置和分速度。

function coe = coe_from_sv（RV）
% computes the classical orbital elements （coe） from the state vector （RV）
% ------------------------------------------------------------
%
% This function computes the classical orbital elements （coe）
% from the state vector （RV） 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
% ------------------------------------------------------------
%   This .m file was from Appendix D of the book:
%        （Howard D. Curtis）
%   You could get appendix D from: http://books.elsevier.com/companions
% ------------------------------------------------------------
%   Last Edit by:           Li yunfei   2008/07/31
% ------------------------------------------------------------

global mu;
eps = 1.e-10;

r = norm（R）;
v = norm（V）;

vr = dot（RV）/r;

H = cross（RV）;
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（NE）/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（ER）/e/r）;
    if vr < 0
        TA = 2*pi - TA;
    end
else
    cp = cross（NR）;
    if cp（3） >= 0
        TA = acos（dot（NR）/n/r）;
    else
        TA = 2*pi - acos（dot（NR）/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];

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2757  2018-06-18 13:57  星历表\coe_from_sv.m

     文件       1388  2017-04-26 20:00  星历表\fDot_and_gDot.m

     文件       1193  2017-04-26 00:05  星历表\f_and_g.m

     文件       3124  2017-04-26 20:24  星历表\interplanetary.m

     文件        972  2008-07-31 09:15  星历表\J0.m

     文件        121  2018-06-18 11:25  星历表\kepler_E.m

     文件       4069  2018-06-18 17:13  星历表\lambert.m

     文件        430  2018-06-18 17:10  星历表\main.m

     文件       1770  2008-07-31 09:10  星历表\month_planet_names.m

     文件       8442  2018-06-18 14:13  星历表\planet_elements_and_sv.m

     文件       1954  2018-06-18 14:57  星历表\skmain.m

     文件        931  2008-07-31 09:06  星历表\stumpC.m

     文件        893  2008-07-31 09:05  星历表\stumpS.m

     文件       2527  2018-06-18 14:13  星历表\sv_from_coe.m

     目录          0  2018-07-04 10:41  星历表

----------- ---------  ---------- -----  ----

                30571                    15