xyz.m

matlab编写的轨道坐标系变换程序 地心直角坐标系与轨道参数互换

function cart=orb2cart(OE)
%ORB2CART  Orbital elements to Cartesian coordinates.
%  Mfile for converting orbital elements ==> r, v vectors
%  input is vector z (6 x 1) of classical orbital elements
%       a is in DU, and all angles are in deg
%
%   OE is the orbital elements which has the form:
%   OE = [a, e, i, W, w, M]
%
%   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
%
%  output is vector x_ECI (r,v vectors in ECI)
%      r is in DU, and v is in DU/TU

function x_ECI = xyz(z);

r2d = 180/pi;
d2r = 1/r2d;
mu=1.0

sma = z(1);
ecc = z(2);
incld = z(3);
raand = z(4);
argpd = z(5);
thetad = z(6);

incl = incld*d2r;
raan = raand*d2r;
argp = argpd*d2r;
theta = thetad*d2r;

%  parameter

p = sma*( 1 - ecc^2);
cth = cos(theta);
sth = sin(theta);

%  radius

r = p/(1 + ecc*cth);

%  compute r vector in perifocal (PQW) frame

rp = r*cth;
rq = r*sth;
rw = 0;
r_PQW = [ rp ; rq ; rw ];

%  compute v vector in PQW frame

vp = -sqrt(1/p)*sth;
vq = sqrt(1/p)*(ecc + cth);
vw = 0;
v_PQW = [ vp ; vq ; vw ];

%  compute the r & v vectors in the ECI frame

cosa = cos(raan);
sina = sin(raan);
cosi = cos(incl);
sini = sin(incl);
cosw = cos(argp);
sinw = sin(argp);

%  rotation matrix R

r11 = cosa*cosw - sina*sinw*cosi;
r12 = -cosa*sinw - sina*cosw*cosi;
r13 = sina*sini;
r21 = sina*cosw + cosa*sinw*cosi;
r22 = -sina*sinw + cosa*cosw*cosi;
r23 = -cosa*sini;
r31 = sinw*sini;
r32 = cosw*sini;
r33 = cosi;

R = [ r11 r12 r13 ; r21 r22 r23 ; r31 r32 r33 ];

r_ECI = R*r_PQW;

v_ECI = R*v_PQW;

r_mag = norm(r_ECI);
v_mag = norm(v_ECI);

E = v_mag^2/2 - mu/r_mag;

Esma = -mu/(2*sma);

rp = p/(1 + ecc);
ra = p/(1 - ecc);

%  store state vector for output

x_ECI = [ r_ECI' v_ECI' ]';