Orbital Elements from Position/Velocity Vectors 1.0

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

Requirements:	No special requirements
Platforms:	Matlab
Keyword:	Grade,  Includes,  Interface,  Moons,  Planets,  Pluto,  Positions,  School
Users rating:	0/10