orb2cart.m

航天工程工具箱

function [R,V,P,dt]=orb2cart(OE,N,I,mu)
%ORB2CART  Orbital elements to Cartesian coordinates.
%   [R,V,P,DT] = ORB2CART(OE,N,I,MU) where N is the number of points along orbit
%      and I is the current position from these set of points (empty if all points).
%   [R,V,P] = ORB2CART(OE,M,MU) where M is the current mean anomaly.
%   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
%
%   R is the position vector (Nx3), V is the velocity vector (Nx3),
%   DT is the Delta t between two adjacent points of position in
%   the output, and P is the period time.
%   If I is included then R and V will be (1x3) vectors describing the
%   position vector at point I of the N points in the orbit.
%
%   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 CART2ORB, TLE2ORB, MA2TP, MSATTRACK.

% Copyright (c) 2002-12-07, B. Rasmus Anthin.
% Revision 2002-12-09, 2003-04-08, 2003-06-09.

%orb2cart(OE,M,mu)
%orb2cart(OE,N,I,mu)

d2r=pi/180;

if nargin<4, mu=I;end
M=N*d2r;
if nargin==4
   if N<1, error('N must be greater than 0.');end
   M=linspace(0,2*pi,N+1); M=M(1:N);
   if N>1
      dM=diff(M(1:2));
   end
   if I<1 | I>N, error('I must be between 1 and N.');end
   if ~isempty(I), M=M(I);end
end
a=OE(1);
e=OE(2);
i=OE(3)*d2r;
W=OE(4)*d2r;
w=OE(5)*d2r;

b=a*sqrt(1-e^2);
rp=a*(1-e);
ra=a*(1+e);
c=sqrt(a^2-b^2);
p=a*(1-e^2);

Ef=inline('M-E+e*sin(E)','E','e','M');
%eps=2/N;
warning off
E(1)=fzero(Ef,0,eps,0,e,M(1));
for k=2:length(M)
   E(k)=fzero(Ef,E(k-1),eps,0,e,M(k));
end
warning on
f=unwrap(2*atan(sqrt(ra/rp)*tan(E/2)));
r=p./(1+e*cos(f));
h=sqrt(mu*p);
P=2*pi*sqrt(a^3/mu);
n=2*pi/P;
if nargin==4 & N>1
   dt=dM/n;
elseif nargin==4
   dt=P;
end


x=r.*[cos(W)*cos(w+f)-sin(W)*sin(w+f)*cos(i)];
y=r.*[sin(W)*cos(w+f)+cos(W)*sin(w+f)*cos(i)];
z=r.*[sin(i)*sin(w+f)];
R=[x(:) y(:) z(:)];

vx=x*h*e./r/p.*sin(f)-h./r.*[cos(W)*sin(w+f)+sin(W)*cos(w+f)*cos(i)];
vy=y*h*e./r/p.*sin(f)-h./r.*[sin(W)*sin(w+f)-cos(W)*cos(w+f)*cos(i)];
vz=z*h*e./r/p.*sin(f)+h./r.*[sin(i)*cos(w+f)];
V=[vx(:) vy(:) vz(:)];