这里是一些轨道动力学中常见物理量的计算机算法的总结。

本文公式较多，在浏览器中将会花较长时间用于渲染公式。

偏近点角

已知：$\nu,e$

求：$EA$

如果$e>(1-10\times{-11})$，那么$EA=0$，否则

\begin{equation}\sin(EA)=\frac{\sqrt{1-e^2}\sin\nu}{1+e\cos\nu}\end{equation}

\begin{equation}\cos(EA)=\frac{e+\cos\nu}{1+e\cos\nu}\end{equation}

\begin{equation} \label{EA}

EA=atan2(\sin EA,\cos EA)\end{equation}

在资料里还查到了双曲线的Hyperbolic Anomaly(HA)，因为是双曲线轨道，先不写在这里了。

平近点角

已知：$\nu,e$

求：$MA$

对于椭圆轨道($e\leq 10\times{-11}$)，首先按照式(\ref{EA})算出偏近点角，然后

\begin{equation} MA=EA-e\sin EA\end{equation}

这个公式是和平均角速度的公式混合食用的：

\begin{equation}

n=\sqrt{\frac{\mu}{\pm a^3}}\end{equation}

正负号是因为双曲线的半长轴是负的。

偏近点角到真近点角

已知：$EA,e$

求：$\nu$

\begin{equation}\sin\nu=\frac{\sqrt{1-e^2}\sin EA}{1-e\cos EA} \end{equation}

\begin{equation}\cos\nu=\frac{\cos EA -e}{1-e\cos EA} \end{equation}

\begin{equation}\nu=atan2(\sin\nu,\cos\nu)\end{equation}

平近点角到偏近点角

已知：$MA,e$

求：$EA$

就是一个牛顿迭代法，这里直接给出实现的代码

double MA2EA(double MA, double ecc)

{

if (ecc < 1.0) {

// elliptic orbit case

double E;

if ((MA < .0 && MA > -pi) || (MA > pi))

E = MA - ecc;

else

E = MA + ecc;

double E_ = MA;

while (fabs(E - E_) > 1e-8) {

E_ = E;

E = E + (MA - E + ecc * sin(E)) / (1 - ecc * cos(E));

}

return E;

}

}

轨道外推的算法

根据上面各个角之间的转换关系，就可以实现轨道外推。写一个类，封装轨道六根数：

class KeplerianState {

public:

double SMA; // semimajor axis, a

double ECC; // eccentricity, e

double INC; // inclination, i

double AOP; // argument of periapsis, \omega

double RAAN; // right ascension of the ascending node, \Omega

double TA; // true anomaly, \phi

KeplerianState() {}

KeplerianState(double a, double e, double i, double omega, double Omega,

double phi, double mu = 3.986004415e14)

: SMA(a), ECC(e), INC(i), AOP(omega), RAAN(Omega), TA(phi),

gravityConst(mu) {}

~KeplerianState() {}

void toCartesian(vec3 *r, vec3 *vel);

void step(double t);

static KeplerianState fromR_V(const vec3 &r, const vec3 &v,

double mu = 3.98600445e14);

public:

double gravityConst;

};

轨道外推就是用平均角速度乘时间得到平近点角，然后再转化至真近点角：

void KeplerianState::step(double t)

{

double MA = TA2MA(TA, ECC);

double n = pow(gravityConst / SMA / SMA / SMA, .5);

MA += n * t; // t seconds later

TA = MA2TA(MA, ECC);

}

近地点、远地点速度

已知：$a,e,\mu$

求：$v_a,v_p$

如果$e > ( 1 - 10\times{−12} )$，那么$v_a=0$，否则

\begin{equation}v_a=\sqrt{\frac{\mu}{a}(\frac{1-e}{1+e})}

\end{equation}

\begin{equation}v_p=\sqrt{\frac{\mu}{a}(\frac{1+e}{1-e})}

\end{equation}

注意这些公式带入运算时都是要化成弧度的！！！