File Exchange

image thumbnail

Orbital Elements from Position/Velocity Vectors

version 1.0 (4.2 KB) by

Convert positions and velocity state vectors to osculating Keplerian orbital elements.

15 Downloads

Updated

View License

vec2orbElem(rs,vs,mus) converts positions (rs) and velocities (vs)
of bodies with gravitational parameters (mus) to Keplerian orbital elements.

INPUTS:
 rs: 3n x 1 stacked initial position vectors:
      [r1(1);r1(2);r1(3);r2(1);r2(2)r2(3);...;rn(1);rn(2);rn(3)]
          or 3 x n matrix of position vectors.
 vs: 3n x 1 stacked initial velocity vectors or 3 x n matrix
 mus: gravitational parameters (G*m_i) where G is the
          gravitational constant and m_i is the mass of the ith body.
          if all vectors represent the same body, mus may be a
          scalar.
 OUTPUTS:
 a: semi-major axes
 e: eccentricities
 E: eccentric anomalies
 I: inclinations
 omega: arguments of pericenter
 Omega: longitudes of ascending nodes
 P: orbital periods
 tau: time of periapsis crossing
 A, B: orientation matrices (see Vinti, 1998)

 All units must be complementary, i.e., if positions are in AUs, and
 time is in days, dx0 must be in AU/day, mus must be in
 AU^3/day^2 (these are the units in solarSystemData.mat).

 The data in solarSystemData.mat was downloaded from JPL's System Web
 Interface (http://ssd.jpl.nasa.gov/?horizons). It includes
 positions for the planets, the sun and pluto (because I went to
 grade school before 2006). Positions for planets with moons are for
 the barycenters.

 Example:
 %solar system oribtal elements
 ssdat = load('solarSystemData.mat');
 rs = ssdat.p0(1:end-3) - repmat(ssdat.p0(end-2:end),9,1);
 vs = ssdat.v0(1:end-3) - repmat(ssdat.v0(end-2:end),9,1);
 mus = ssdat.mus(1:9) + ssdat.mus(10);
 [a,e,E,I,omega,Omega,P,tau,A,B] = vec2orbElem(rs,vs,mus);
 %convert back:
 r = A*diag(cos(E) - e) + B*diag(sin(E));
 rdot = (-A*diag(sin(E))+B*diag(cos(E)))*...
         diag(sqrt(mus(:).'.*a.^-3)./(1 - e.*cos(E)));

Comments and Ratings (1)

How did you solve the problems with equatorial orbit and circular inclinated orbit?

MATLAB Release
MATLAB 7.9 (R2009b)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video