View License

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

» Watch video

Highlights from
Orbital Elements from Position/Velocity Vectors

Join the 15-year community celebration.

Play games and win prizes!

» Learn more

Be the first to rate this file! 5 Downloads (last 30 days) File Size: 4.2 KB File ID: #31333 Version: 1.0

Orbital Elements from Position/Velocity Vectors



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

| Watch this File

File Information

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

 rs: 3n x 1 stacked initial position vectors:
          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
 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 ( 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.

 %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)));

MATLAB release MATLAB 7.9 (R2009b)
Other requirements Tested on MATLAB 7.8 - 7.11
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (1)
28 Oct 2015 Carlo Convenevole

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

Comment only

Contact us