兰伯特(Lambert)方程的求解算法1

本文针对兰伯特方程给出具体的算法，并不打算给出详细的过程。各位读者可参照此算法及相应的代码进行编程计算。

介绍

见下图，仅考虑中心天体C的万有引力，飞行器从P1P_1P1点飞行到P2P_2P2点，飞行时间为Δt\Delta tΔt。起点P1P_1P1的地心距为r1r_1r1,终点P2P_2P2的地心距为r2r_2r2,飞行路径所对应的地心角为θ\thetaθ。

兰伯特首先发现此两点轨道转移所需要的时间仅与始末的几何构型以及轨道半长径有关，而与其它轨道参数无关，其数学表述为：
μΔt=F(a,r1+r2,c)\sqrt\mu\Delta t =F(a,r_1+r_2,c) μΔt=F(a,r1+r2,c)

拉格朗日和高斯曾经都研究过这个问题，并给出相应的表达式，有兴趣的读者可直接参考相关教科书，此处不再详述。

下面介绍比较高效的求解算法，主要由R.H.Gooding在1990年的数篇文献中给出。想要详细了解，可参考他的两篇文献：

A procedure for the solution of Lambert’s orbital boundary-value problem
On the Solution of Lambert’s Orbital Boundary-Value Problem

本文基于上述两篇文献，直接给出算法。

初始化

由上图兰伯特转移轨道的几何构型，初始输入参数如下：

r1r_1r1，转移轨道起点P1P_1P1到引力中心C的距离 (m)；
r2r_2r2，转移轨道起点P2P_2P2到引力中心C的距离(m) ；
θ\thetaθ，转移轨道的地心夹角，不超过360°；
Δt\Delta tΔt，转移轨道的飞行时间(s)；
mmm，椭圆轨道时，转移轨道的飞行圈数，m=0m=0m=0时，表示转移轨道的地心角 θ\thetaθ不超过360°；m>0m>0m>0时，仅表示椭圆转移轨道，指转移轨道的地心角超过了360°。
μ\muμ，中心天体的引力常数；
由初始输入参数，可得部分参数：
c=(r1−r2)2+4r1r2sin2(θ/2)c=\sqrt{(r_1-r_2)^2+4r_1r_2sin^2(\theta/2)}c=(r1−r2)2+4r1r2sin2(θ/2)；
s=(r1+r2+c)/2s=(r_1+r_2+c)/2s=(r1+r2+c)/2 ；
q=r1r2cos(θ/2)/sq=\sqrt{r_1r_2}cos(\theta/2)/sq=r1r2cos(θ/2)/s；
1−q2=c/s1-q^2=c/s1−q2=c/s，为了避免计算中有效位数的损失，此参数单独计算；
T=8μ/s3ΔtT=\sqrt{8\mu/s^3}\Delta tT=8μ/s3Δt: 无量纲飞行时间；

转移时间T及其导数的求解（TLAMB）

输入：

mmm
qqq
1−q21-q^21−q2
xxx
nnn，n=0表示仅计算时间T，n=1表示计算一阶导数，n=2表示计算到二阶导数，n=3表示计算到三阶导数。

输出：T,T′,T′′,T′′′T,T',T'',T'''T,T′,T′′,T′′′
TTT的导数(T′,T′′,T′′′T',T'',T'''T′,T′′,T′′′)指TTT对xxx的1-3阶导数。

无量纲转移时间TTT可表达成如下形式：
T(m,q,x)(1)T(m,q,x)\tag1 T(m,q,x)(1)
其中，qqq与转移轨道几何构型相关，自变量xxx与轨道半长轴aaa相关，有：x2=1−s/(2a)x^2=1-s/(2a)x2=1−s/(2a)

根据参数m,q,xm,q,xm,q,x的初值，TTT的计算有两种方式：直接计算与级数计算。下面分别阐述。

直接计算方式

当下式成立时，，采用直接计算方式。
m>0或x<0或∣u∣>0.4(2)m>0 或 x<0 或|u|>0.4\tag2m>0或x<0或∣u∣>0.4(2)
其中
u=1−x2u=1-x^2u=1−x2
因此，针对xxx的具体范围为：x<0.6x<\sqrt{0.6}x<0.6或x>1.4x>\sqrt{1.4}x>1.4，即xxx的范围不在1附近（抛物线轨道）
定义
z=1−q2+q2x2z= \sqrt{1-q^2+q^2x^2}z=1−q2+q2x2
当qx≤0qx\le0qx≤0时，
{α=z−qxβ=qz−x\left\{ \begin{array}{l} \alpha=z-qx \\ \beta=qz-x \\ \end{array} \right. {α=z−qxβ=qz−x
当qx>0qx>0qx>0时，
{α=1−q2z+qxβ=(1−q2)(q2u−x2)qz+x\left\{ \begin{array}{l} \alpha=\frac{1-q^2}{z+qx} \\ \beta=\frac{(1-q^2)(q^2u-x^2)}{qz+x} \\ \end{array} \right. {α=z+qx1−q2β=qz+x(1−q2)(q2u−x2)
令
{y=∣u∣f=αy\left\{ \begin{array}{l} y=\sqrt{|u|} \\ f=\alpha y \end{array} \right. {y=∣u∣f=αy
且令
g={xz+qu,if qxu≥0x2−q2uxz−qu,if qxu<0g = \begin{cases} xz+qu, & \text{if $qxu\ge0$} \\ \frac{x^2-q^2u}{xz-qu} , & \text{if $qxu<0$} \end{cases} g={xz+qu,xz−qux2−q2u,if qxu≥0if qxu<0
则参数ddd为
d={mπ+arctan2(f,g)if x≤1tanh−1(f/g)=ln(f+g)if x>1d = \begin{cases} m\pi+arctan2(f,g) & \text{if $x\le1$} \\ tanh^{-1}(f/g)=ln(f+g) & \text{if $x>1$} \end{cases} d={mπ+arctan2(f,g)tanh−1(f/g)=ln(f+g)if x≤1if x>1
上式中，arctan2(f,g)arctan2(f,g)arctan2(f,g)为反正切函数tan−1(f/g)tan^{-1}(f/g)tan−1(f/g)的扩充，由fff,ggg 的符号决定了返回值的正确象限。
当x>1且f<0.4x>1且f<0.4x>1且f<0.4时，文献中使用级数方式计算ln(f+g)ln(f+g)ln(f+g)，经测试，误差不超过10−1510^{-15}10−15，因此此处仅使用ln(f+g)ln(f+g)ln(f+g)。

最后，直接给出 TTT及其导数的计算公式：
{T=2(d/y+β)/uT′=(3xT+4q3x/z−4)/uT′′=(3T+5xT′+4(q/z)3(1−q2))/uT′′′=(8T′+7xT′′−12x(q/z)5(1−q2))/u(3)\begin{cases} T=2(d/y+\beta)/u \\ T'=(3xT+4q^3x/z-4)/u \\ T''=(3T+5xT'+4(q/z)^3(1-q^2))/u \\ T'''=(8T'+7xT''-12x(q/z)^5(1-q^2))/u \end{cases}\tag3 ⎩⎪⎪⎪⎨⎪⎪⎪⎧T=2(d/y+β)/uT′=(3xT+4q3x/z−4)/uT′′=(3T+5xT′+4(q/z)3(1−q2))/uT′′′=(8T′+7xT′′−12x(q/z)5(1−q2))/u(3)
注意，上式中，参数uuu不可能为0；但参数zzz有可能为0，因此当z=0z=0z=0时，不计算T的导数(T′,T′′,T′′′T',T'',T'''T′,T′′,T′′′)，直接赋值为0。

级数计算方式

当直接计算方式的条件（式2）不成立时，则采用级数计算方式。此时x在1附近，接近抛物线。
此时无量纲时间 TTT的级数表达式为：
T=∑n=0∞Anbnun(4)T=\sum_{n=0}^\infty A_nb_nu^n\tag4T=n=0∑∞Anbnun(4)
其中
{An=(2n)!22n−2(n!)2(2n+3)bn=1−q2n+3\begin{cases} A_n =\frac{(2n)!}{2^{2n-2}(n!)^2(2n+3)} \\ b_n =1-q^{2n+3} \end{cases}{An=22n−2(n!)2(2n+3)(2n)!bn=1−q2n+3
在实际求解时，式(4)变为下式：
Tx2=∑n=0∞Cnun(5)Tx^2=\sum_{n=0}^\infty C_nu^n\tag5Tx2=n=0∑∞Cnun(5)
其中
Cn=−6n+14n2−1Anbn+An−1(bn−bn−1)C_n=-\frac{6n+1}{4n^2-1}A_nb_n+A_{n-1}(b_n-b_{n-1}) Cn=−4n2−16n+1Anbn+An−1(bn−bn−1)
上式中，AnA_nAn的递推公式为：
An=an2n+3an=2n−12nan−1A_n =\frac{a_n}{2n+3}\\ a_n =\frac{2n-1}{2n}a_{n-1} An=2n+3anan=2n2n−1an−1
bnb_nbn的递推公式为：
bn=bn−1+q2n+1(1−q2)b_n =b_{n-1}+q^{2n+1}(1-q^2)bn=bn−1+q2n+1(1−q2)
初值为：
C0=A0b0C_0 =A_0b_0C0=A0b0
其中，A0=4/3A_0 =4/3A0=4/3，b0b_0b0则由下式给出
b0={1−q3if q<0.5(q+1/(1+q))(1−q2)if q≥0.5b_0 = \begin{cases} 1-q^3 & \text{if $q<0.5$} \\ (q+1/(1+q))(1-q^2) & \text{if $q \ge0.5$} \end{cases} b0={1−q3(q+1/(1+q))(1−q2)if q<0.5if q≥0.5
最后，直接给出 TTT的导数计算公式：
{T′=−2xTu′T′′=−2Tu′+4x2Tu′′T′′′=12xTu′′−8x3Tu′′′(6)\begin{cases} T'=-2xT'_u \\ T''=-2T'_u+4x^2T''_u \\ T'''=12xT''_u-8x^3T'''_u \end{cases}\tag6 ⎩⎪⎨⎪⎧T′=−2xTu′T′′=−2Tu′+4x2Tu′′T′′′=12xTu′′−8x3Tu′′′(6)
也就是说，TTT对xxx的导数(T′,T′′,T′′′T',T'',T'''T′,T′′,T′′′)依赖TTT对uuu的导数(Tu′,Tu′′,Tu′′′T'_u,T''_u,T'''_uTu′,Tu′′,Tu′′′)

TTT对uuu的导数(Tu′,Tu′′,Tu′′′T'_u,T''_u,T'''_uTu′,Tu′′,Tu′′′)可由（4）式推导，得到：

Tu′=∑n=0∞Anbnnun−1Tu′′=∑n=0∞Anbnn(n−1)un−2Tu′′′=∑n=0∞Anbnn(n−1)(n−2)un−3(7)T'_u= \sum_{n=0}^\infty A_nb_nnu^{n-1} \\ T''_u= \sum_{n=0}^\infty A_nb_nn(n-1)u^{n-2} \\ T'''_u= \sum_{n=0}^\infty A_nb_nn(n-1)(n-2)u^{n-3} \tag 7 Tu′=n=0∑∞Anbnnun−1Tu′′=n=0∑∞Anbnn(n−1)un−2Tu′′′=n=0∑∞Anbnn(n−1)(n−2)un−3(7)
上式1-3阶导数的求和过程可与式（5）的求和过程同步。

下图给出了转移时间T（无量纲）的图像。横坐标为xxx,其范围为：x∈(−1,∞)x\in(-1,\infty)x∈(−1,∞)，图中只画到1.5。
x∈(−1,1)x\in(-1,1)x∈(−1,1)时，为椭圆轨道；
x=1x=1x=1时，为抛物线轨道；
x∈(1,∞)x\in(1,\infty)x∈(1,∞)时，为双曲轨道；
此外，图中仅画出了m=0,1m=0,1m=0,1的图像。

源码(C#)

//           Lambert 方程的无量纲时间 (R.H.Gooding 方法)
//--------------------------------------------------------------------
//   1   t=sqrt(8u/s^3)*Δt= 2*pi*m/(1-x*x)^1.5+
//           4/3*(F[3,1;2.5;0.5*(1-x)]-q^3*F[3,1;2.5;0.5*(1-y)])
//       其中:   y=sqrt(1-q*q+q^2*x^2)=sqrt(qsqfm1+q^2*x^2)
//   2   此算法根据不同情况，用直接计算法和级数法来计算
//       *与文献中的算法不同之处：f<0.4时(q接近1,x>1.2)，未采用级数循环计算，而直接采用
//          log(f+g)计算，经测试，误差小于1e-15。
//   3   算法引用的文献为：
//       1   Gooding,R.H.:: 1988a,'On the Solution of Lambert's Orbital
//           Boundary-Value Problem',RAE Technical Report 88027
//       2   Gooding,R.H.:: 1989, 'A Procedure for the Solution of Lambert's
//           Orbital Boundary
//   4   主要由子程序XLamb调用进行迭代寻根使用
//   5   注意输入参数n,当为-1时，只计算参数qz-x/qz+x/qx+z/// <summary>
/// Lambert 方程的无量纲时间 (R.H.Gooding 方法)
/// <para>主要由子程序XLamb调用进行迭代寻根使用</para>
/// </summary>
/// <param name="m">飞行的圈数(0 for 0-2pi)</param>
/// <param name="q">q=sqrt(r1*r2)/s*cos(0.5*theta)</param>
/// <param name="qsqfm1">(1-q*q): c/s</param>
/// <param name="x">自变量(x*x=1-am/a)</param>
/// <param name="n">0仅计算t; 2计算到d2t; 3计算到d3t  -1仅计算qz-x/qz+x/qx+z</param>
/// <param name="t">out:无量纲时间(T=sqrt(8u/s^3)*Δt)</param>
/// <param name="dt">out:T对x的1阶导数; n=-1时为qz-x</param>
/// <param name="d2t">out:T对x的2阶导数;n=-1时为qz+x</param>
/// <param name="d3t">out:T对x的3阶导数;n=-1时为qx+z</param>
public static void TLAMB(int m, double q, double qsqfm1, double x, int n, out double t, out double dt, out double d2t, out double d3t)
{dt = d2t = d3t = 0.0;           //  缺省值double sw = 0.4;                //  初值bool lm1 = (n == -1);           //是否计算qz-x/qz+x/qx+z//  导数计算的判断标志bool l1 = (n >= 1);bool l2 = (n >= 2);bool l3 = (n == 3);double qsq = q * q;             //q^2double xsq = x * x;             //x^2double u = (1 - x) * (1 + x);   //u=1-x^2#region 直接计算(非级数计算)，含n=-1时的情形if (lm1 || m > 0 || x < 0 || Math.Abs(u) > sw){double y = Math.Sqrt(Math.Abs(u));double z = Math.Sqrt(qsqfm1 + qsq * xsq);double qx = q * x;//  alpha,betadouble a = 0;       //z-qxdouble b = 0;       //qz-x  n=-1时使用dt返回double aa = 0;      //z+qx  n=-1时使用d2t返回double bb = 0;      //qz+x  n=-1时使用d3t返回if (qx <= 0){a = z - qx;b = q * z - x;}if (qx < 0 && lm1){aa = qsqfm1 / a;bb = qsqfm1 * (qsq * u - xsq) / b;}if (qx == 0.0 && lm1 || qx > 0){aa = z + qx;bb = q * z + x;}if (qx > 0.0){a = qsqfm1 / aa;b = qsqfm1 * (qsq * u - xsq) / bb;}//  如果n=-1，则返回，并通过dt/d2t/d3t返回:qz-x/qz+x/qx+zif (lm1){t = 0;dt = b;d2t = bb;d3t = aa;return;}double g;if (qx * u >= 0)g = x * z + q * u;elseg = (xsq - qsq * u) / (x * z - q * u);double f = a * y;//  椭圆情形if (x <= 1.0)t = m * Math.PI + Math.Atan2(f, g);//  双曲线情形else{if (f > sw)  //正常计算Log(f+g)t = Math.Log(f + g);#region 使用级数计算Log(f+g)else{double told;double fg1 = f / (g + 1.0);double term = 2.0 * fg1;double fg1sq = fg1 * fg1;t = term;double twoi1 = 1.0;do{twoi1 = twoi1 + 2.0;term = term * fg1sq;told = t;t = t + term / twoi1;} while (t != told);}#endregion}t = 2.0 * (t / y + b) / u;  //无量纲时间if (l1 && z != 0)           //1-3阶导数{double qz = q / z;double qz2 = qz * qz;double qz3 = qz2 * qz;dt = (3.0 * x * t - 4.0 * (a + qx * qsqfm1) / z) / u;if (l2)d2t = (3.0 * t + 5.0 * x * dt + 4.0 * qz3 * qsqfm1) / u;if (l3)d3t = (8.0 * dt + 7.0 * x * d2t - 12.0 * qz3 * qz2 * x * qsqfm1) / u;}}#endregion#region 级数计算( T*x^2=Sum(Cn*u^n) n=0,1,2...)else{//  求解T时级数中u的n次方double u0i = 1.0;//  求解T的1-3阶导数时级数中u的n次方double u1i = 0;double u2i = 0;double u3i = 0;if (l1) u1i = 1.0;if (l2) u2i = 1.0;if (l3) u3i = 1.0;double term = 4.0;          //a0double tq = q * qsqfm1;     //q*(1-q^2)double tqsum = 0.0;         //b0if (q < 0.5) tqsum = 1.0 - q * qsq;if (q >= 0.5) tqsum = (1.0 / (1.0 + q) + q) * qsqfm1;double ttmold = term / 3.0; //A0t = ttmold * tqsum;         //A0*b0: T*x^2的级数的0次方double told = t - 1.0;      //T*x^2的级数前一值，force t ~= told to get one pass throughint i = 0;double tterm, tqterm;       //An,An*bnwhile (i < n || t != told){i = i + 1;u0i = u0i * u;  //u^nif (l1 && i > 1) u1i = u1i * u;if (l2 && i > 2) u2i = u2i * u;if (l3 && i > 3) u3i = u3i * u;term = term * (i - 0.5) / i;    //antq = tq * qsq;                  //bn-b(n-1)tqsum = tqsum + tq;             //bn=b(n-1)+q^(2n)*q*(1-q^2)told = t;                       //保留T*x^2级数的当前值tterm = term / (2.0 * i + 3.0); //An=an/(2n+3)tqterm = tterm * tqsum;         //An*bn// T*x^2的级数: +当前的Cn*u^nt = t - u0i * ((1.5 * i + 0.25) * tqterm / (i * i - 0.25) - ttmold * tq);ttmold = tterm;                 //保留当前的Antqterm = tqterm * i;            //An*bn*n//  T对u的1-3阶级数:+当前的An*bn*n*..if (l1) dt = dt + tqterm * u1i;if (l2) d2t = d2t + tqterm * u2i * (i - 1.0);if (l3) d3t = d3t + tqterm * u3i * (i - 1.0) * (i - 2.0);}//  T对x的1-3阶导数if (l3) d3t = 8.0 * x * (1.5 * d2t - xsq * d3t);if (l2) d2t = 2.0 * (2.0 * xsq * d2t - dt);if (l1) dt = -2.0 * x * dt;//  T:最后求得的级数需要除以x^2t /= xsq;}#endregion
}