资源分类:Matlab 工具:MATLAB 7.9 (R2009b)


Orbital Elements from Position/Velocity Vectors

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

文件下载列表 (4.18KB)  
vec2orbElem.m  ._vec2orbElem.m  


用户反馈   关于我们
Copyright (©) ZHIHUISHI.COM 2013 All Rights Reserved.