cart2orb.m

航天工程工具箱

function oe=cart2orb(rr,vv,mu)
%CART2ORB  Cartesian coordinates to orbital elements.
%   OE = ORB2CART(R,V,MU) where R is the position vector (1x3)
%   and V is the velocity vector (1x3) at a given time point (epoch).
%   MU is G*M for the central object,
%   and OE is the orbital elements which has the form:
%
%   OE = [a, e, i, W, w, M]
%
%   However, only the first five elements will be used by this function.
%
%   Where:
%      a [m]   : semi-major axis
%      e []    : eccentricity
%      i [deg] : inclination
%      W [deg] : longitude of the ascending node
%      w [deg] : argument of periapsis
%      M [deg] : mean anomaly at epoch
%
%   M is the mean anomaly at epoch, and is here the same as the mean anomaly
%   for the position R and velocity V given above. If R is the location of the
%   periapsis, then the mean anomaly will be M = 0.
%
%   The xy-plane is the reference plane and the zero point of longitude
%   (or direction of vernal equinox) is in the direction of the x-axis.
%
%   See also ORB2CART, TLE2ORB, MA2TP, MSATTRACK.

% Copyright (c) 2003-04-18, B. Rasmus Anthin.

r2d=180/pi;
rr=rr(1,:);
vv=vv(1,:);
hh=cross(rr,vv);
nn=cross([0 0 1],hh);
r=norm(rr);
v=norm(vv);
h=norm(hh);
n=norm(nn);
ee=((v^2-mu/r)*rr-(rr*vv')*vv)/mu;
e=norm(ee);

W=acos(nn(1)/n);
if nn(2)<0, W=2*pi-W;end
i=atan(n/hh(3));
a=mu*r/(2*mu-r*v^2);
E=atan2(r*v/(e*sqrt(mu*a)),(a-r)/(a*e));
M=E-e*sin(E);
w=acos(nn*ee'/n/e);
if ee(3)<0, w=2*pi-w;end

oe=[a e i*r2d W*r2d w*r2d M*r2d];