Orbital Elements from Position/Velocity Vectors

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

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