在前一文章中，介绍了兰伯特方程的基本概念，并给出了无量纲飞行时间 T T T的具体的算法，本节给出上述逆过程，即由无量纲飞行时间 T T T求解自变量 x x x。
已知 T T T求解自变量 x x x采用函数求根方法，即求解下式的根：
T ( m , q , x ) ? T i = 0 T(m,q,x)-T_i=0 T(m,q,x)?Ti?=0
此处采用Halley迭代法(需要使用到 T T T对 x x x的2阶导数）。
自变量X的求解（XLAMB） 输入：
- m m m
- q q q
- 1 ? q 2 1-q^2 1?q2
- T i T_i Ti?,无量纲飞行时间
5.n n n,n=-1:非正常返回；n=0:无解；n=1:1个解；n=2:2个解
6.x x x,第1个解
7.x + x_+ x+?,第2个解（n=2时）
初值设定 x 0 x_0 x0?
为了使程序计算中迭代算法的高效快速，合适的初值选取是非常重要的。
首先给出本小节公式中一些常数值：
{ c 0 = 1.7 c 1 = 0.5 c 2 = 0.03 c 3 = 0.15 c 41 = 1 c 42 = 0.24 \begin{cases} c_0=1.7 \\ c_1=0.5 \\ c_2=0.03 \\ c_3=0.15 \\ c_{41}=1 \\ c_{42} = 0.24 \end{cases} ????????????????????c0?=1.7c1?=0.5c2?=0.03c3?=0.15c41?=1c42?=0.24?
转移角度 θ \theta θ
θ / 2 = a r c t a n 2 ( 1 ? q 2 , 2 q ) \theta/2 =arctan2(1-q^2,2q) θ/2=arctan2(1?q2,2q)
根据转移圈数 m m m的不同， x 0 x_0 x0?的初值也不同，下面分情况讨论。
m = 0 m=0 m=0时 此时对应转移角度不超过360°情形，由 T T T的图像可知， x x x必有解，且只有一个解。
因此，输出参数 n = 1 n=1 n=1
首先计算 x = 0 x=0 x=0时的无量纲时间 T 0 T_0 T0?,即： T 0 = T ( q , m = 0 , x = 0 ) T_0=T(q,m=0,x=0) T0?=T(q,m=0,x=0)
- T i ≤ T 0 T_i\le T_0 Ti?≤T0?
x 0 = ? T 0 ( T i ? T 0 ) 4 T i x_0=-\frac{T_0(T_i-T_0)}{4T_i} x0?=?4Ti?T0?(Ti??T0?)? - T i > T 0 T_i>T_0 Ti?>T0?
令
x 01 = ? T i ? T 0 T i ? T 0 + 4 x 02 = ? T i ? T 0 T i + T 0 / 2 W = x 01 + c 0 2 ? θ / π x_{01}=-\frac{T_i-T_0}{T_i-T_0+4} \\ x_{02}=-\sqrt{\frac{T_i-T_0}{T_i+T_0/2}} \\ W=x_{01}+c_0\sqrt{2-\theta/\pi} x01?=?Ti??T0?+4Ti??T0??x02?=?Ti?+T0?/2Ti??T0?? ?W=x01?+c0?2?θ/π ?
有
x 03 = { x 01 ifW ≥ 0 x 01 + ( ? W ) 1 / 16 ( x 02 ? x 01 ） ifW < 0 x_{03} = \begin{cases} x_{01} & \text{if $W\ge0$} \\ x_{01}+(-W)^{1/16}(x_{02}-x_{01}）& \text{if $W<0$} \end{cases} x03?={x01?x01?+(?W)1/16(x02??x01?）?if W≥0if W<0?
则有：
x 0 = λ x 03 x_0=\lambda x_{03} x0?=λx03?
其中：
λ = 1 + c 1 x 03 ( 1 + x 01 ) ? c 2 x 03 2 1 + x 01 \lambda = 1+c_1x_{03}(1+x_{01})-c_2x^2_{03}\sqrt{1+x_{01}} λ=1+c1?x03?(1+x01?)?c2?x032?1+x01? ?
最小时间 T m T_m Tm?及 x m x_m xm?的求解 首先求解出最小时间 T m T_m Tm?及其自变量 x m x_m xm?，也即寻找方程 T ′ ( q , m , x ) = 0 T'(q,m,x)=0 T′(q,m,x)=0 的根 ，其迭代过程采用Halley方法的迭代公式，因此在迭代求解过程中需要求解 T T T的三阶导数。
迭代时，仍需要一个合适的初值 x m x_m xm?，其由下式给出：
x m = { x m , π = 4 3 π ( 2 m + 1 ) ifθ = π x m , π ( θ π ) 1 / 8 ifθ < π x m , π ( 2 ? ( 2 ? θ π ) 1 / 8 ) ifθ > π x_m= \begin{cases} x_{m,\pi} =\frac{4}{3\pi(2m+1)} & \text{if $\theta =\pi$} \\ x_{m,\pi}(\frac{\theta}{\pi})^{1/8} & \text{if $\theta <\pi$} \\ x_{m,\pi}(2-(2-\frac{\theta}{\pi})^{1/8}) & \text{if $\theta >\pi$} \end{cases} xm?=??????xm,π?=3π(2m+1)4?xm,π?(πθ?)1/8xm,π?(2?(2?πθ?)1/8)?if θ=πif θ<πif θ>π?
有了初值后，采用halley迭代法求解方程 T ′ ( q , m , x ) = 0 T'(q,m,x)=0 T′(q,m,x)=0 的根。
迭代公式为：
T m , T m ′ , T m ′ ′ , T m ′ ′ ′ = T ( m , q , 1 ? q 2 , x m , 3 ) x m = x m ? T ′ T ′ ′ T ′ ′ T ′ ′ ? T ′ T ′ ′ ′ / 2 T_m,T'_m,T''_m,T'''_m=T(m,q,1-q^2,x_m,3) \\ x_m=x_m-\frac{T'T''}{T''T''-T'T'''/2} Tm?,Tm′?,Tm′′?,Tm′′′?=T(m,q,1?q2,xm?,3)xm?=xm??T′′T′′?T′T′′′/2T′T′′?
每一步迭代更新 x m x_m xm?前保留其值为 x m ′ = x m x'_m=x_m xm′?=xm?.
限定最大迭代次数为12次，迭代过程中，若 T ′ ′ = = 0 T''==0 T′′==0，则跳出循环；若 ∣ x m ′ / x m ? 1 ∣ < 3 × 1 0 ? 7 |x'_m/x_m-1|<3×10^{-7} ∣xm′?/xm??1∣<3×10?7也跳出循环。
若迭代次数超过12次，则令 n = ? 1 n=-1 n=?1，程序退出。
迭代完成后即可得到： T m , x m , T m ′ ′ T_m,x_m,T''_m Tm?,xm?,Tm′′?，注意，若最终求得的 T m ′ ′ = 0 T''_m=0 Tm′′?=0，则令 T m ′ ′ = 6 m π T''_m=6m\pi Tm′′?=6mπ。
得到 T m T_m Tm?后，则有：
- 若 T i < T m T_i
- 若 T i = T m T_i=T_m Ti?=Tm?，令 x = x m , n = 1 x=x_m,n=1 x=xm?,n=1，程序退出；
- 若 T i > T m T_i>T_m Ti?>Tm?，令 n = 2 n=2 n=2，有两解。由下面给出具体初值 x 0 + 、 x 0 ? x_{0+}、x_{0-} x0+?、x0??。
x 01 = T i ? T m T m ′ ′ 2 + T i ? T m ( 1 ? x m ) 2 x_{01}=\sqrt{\frac{T_i-T_m}{\frac{T''_m}{2}+\frac{T_i-T_m}{(1-x_m)^2}}} \\ x01?=2Tm′′??+(1?xm?)2Ti??Tm??Ti??Tm?? ?
则初值 x 0 + x_{0+} x0+?为：
x 0 + = x m + λ x 01 x_{0+}=x_m+\lambda x_{01} x0+?=xm?+λx01?
上式中，
{ W = 4 ( x m + x 01 ) 4 + T ? T m + ( 1 ? x m ? x 01 ) 2 λ ′ = ( 1 + m + c 41 ( θ / 2 π ? 0.5 ) ) / ( 1 + c 3 m ） λ = 1 ? x 01 λ ′ ( c 1 W + c 2 x 01 W ) \begin{cases} W=\frac{4(x_m+x_{01})}{4+T-T_m}+(1-x_m-x_{01})^2 \\ \lambda'=(1+m+c_{41}(\theta/2\pi-0.5))/(1+c_3m）\\ \lambda=1-x_{01}\lambda'(c_1W+c_2x_{01}\sqrt{W}) \end{cases} ??????W=4+T?Tm?4(xm?+x01?)?+(1?xm??x01?)2λ′=(1+m+c41?(θ/2π?0.5))/(1+c3?m）λ=1?x01?λ′(c1?W+c2?x01?W ?)?
x < x m x
- T i ≤ T 0 T_i\le T_0 Ti?≤T0?时
x 0 ? = x m ? T i ? T m T m ′ ′ / 2 ? ( T i ? T m ) ( T m ′ ′ 2 ( T 0 ? T m ) ? 1 x m 2 ) x_{0-}=x_m-\sqrt{\frac{T_i-T_m}{T''_m/2-(T_i-T_m)(\frac{T''_m}{2(T_0-T_m)}-\frac{1}{x_m^2})}} x0??=xm??Tm′′?/2?(Ti??Tm?)(2(T0??Tm?)Tm′′???xm2?1?)Ti??Tm?? ? - T i > T 0 T_i> T_0 Ti?>T0?时
令
x 01 = ? T i ? T 0 T i ? T 0 + 4 x 02 = ? T i ? T 0 T i + T 0 / 2 W = x 01 + c 0 2 ? θ / π x_{01}=-\frac{T_i-T_0}{T_i-T_0+4} \\ x_{02}=-\sqrt{\frac{T_i-T_0}{T_i+T_0/2}} \\ W=x_{01}+c_0\sqrt{2-\theta/\pi} x01?=?Ti??T0?+4Ti??T0??x02?=?Ti?+T0?/2Ti??T0?? ?W=x01?+c0?2?θ/π ?
有
x 03 = { x 01 ifW ≥ 0 x 01 + ( ? W ) 1 / 16 ( x 02 ? x 01 ） ifW < 0 x_{03} = \begin{cases} x_{01} & \text{if $W\ge0$} \\ x_{01}+(-W)^{1/16}(x_{02}-x_{01}）& \text{if $W<0$} \end{cases} x03?={x01?x01?+(?W)1/16(x02??x01?）?if W≥0if W<0?
则有：
x 0 ? = λ x 03 x_{0-}=\lambda x_{03} x0??=λx03?
其中：
{ λ = 1 + c 1 x 03 λ ′ ( 1 + x 01 ) ? c 2 x 03 2 λ ′ 1 + x 01 λ ′ = ( 1 + m + c 42 ( θ / 2 π ? 0.5 ) ) / ( 1 + c 3 m ） \begin{cases} \lambda = 1+c_1x_{03}\lambda'(1+x_{01})-c_2x^2_{03}\lambda'\sqrt{1+x_{01}} \\ \lambda'=(1+m+c_{42}(\theta/2\pi-0.5))/(1+c_3m） \end{cases} {λ=1+c1?x03?λ′(1+x01?)?c2?x032?λ′1+x01? ?λ′=(1+m+c42?(θ/2π?0.5))/(1+c3?m）?
有了初值后，采用Halley迭代法。迭代公式为：
T , T ′ , T ′ ′ = T ( m , q , 1 ? q 2 , x 0 , 2 ) x 0 = x 0 + ( T i ? T ) T ′ T ′ T ′ + ( T i ? T ) T ′ ′ / 2 T,T',T''=T(m,q,1-q^2,x_0,2) \\ x_0=x_0+\frac{(T_i-T)T'}{T'T'+(T_i-T)T''/2} T,T′,T′′=T(m,q,1?q2,x0?,2)x0?=x0?+T′T′+(Ti??T)T′′/2(Ti??T)T′?
迭代次数设定为3次即可。
迭代完成后，令 x = x 0 x=x_0 x=x0?，程序退出；
若 m > 0 m>0 m>0，即有两解的情况下，分别将初值 x 0 + 、 x 0 ? x_{0+}、x_{0-} x0+?、x0??带入上述迭代式，最终得到的解分别赋值给 x 、 x + x、x_+ x、x+?
C#源码
//Lambert 无量纲方程x的求解 (R.H.Gooding 方法)
//--------------------------------------------------------------------
//1已知tin,寻找下列方程的根x
//tin=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初值的选取是根据bilinear curve近似所得(分段分析)；求根迭代过程
//是根据Halley's method来iteration
//3算法引用的文献为：
//1Gooding,R.H.:: 1988a,'On the Solution of Lambert's Orbital
//Boundary-Value Problem',RAE Technical Report 88027
//2Gooding,R.H.:: 1989, 'A Procedure for the Solution of Lambert's
//Orbital Boundary
//--------------------------------------------------------------------///
/// Lambert 无量纲方程x的求解 (R.H.Gooding 方法)
///
/// 【兰伯特(Lambert)方程的求解算法2】飞行的圈数(0 for 0-2pi)
/// q=sqrt(r1*r2)/s*cos(0.5*theta)
/// (1-q*q): c/s
/// 无量纲飞行时间T
/// out:-1:异常;
0:无解(m>0,T0)
/// out:解1
/// out:解2
static void XLAMB(int m, double q, double qsqfm1, double tin, out int n, out double x, out double xpl)
{
double dt, d2t, d3t;
//系数
double tol = 3.0e-7;
double c0 = 1.7;
double c1 = 0.5;
double c2 = 0.03;
double c3 = 0.15;
double c41 = 1.0;
double c42 = 0.24;
//初始化
x = xpl = 0.0;
double t0 = 0.0;
//T(x=0)
double t = 0.0;
double tmin = 0.0;
//多圈时,T的最小值
double xm = 0.0;
//多圈时，T最小值的x
double tdiffm = 0;
//tin-tmin
double d2t2 = 0.0;
//T的二阶导数(x=xm)/2//theta/2pi
double thr2 = Math.Atan2(qsqfm1, 2.0 * q) / Math.PI;
#region1圈内转移,x的初值(x>-1;
可能为 椭圆，抛物线，双曲线)
if (m == 0)
{
//x仅有一个解
//Single-rev starter fromt (at x = 0) & bilinear (usually)
n = 1;
//求x=0时对应的T，以此判断x是否大于0，T随x单调减
TLAMB(m, q, qsqfm1, 0.0, 0, out t0, out dt, out d2t, out d3t);
double tdiff = tin - t0;
//当 x>0(tin <= t0) 时(bilinear curve拟合产生初始x0)
if (tdiff <= 0.0)
{
x = t0 * tdiff / (-4.0 * tin);
}
//当 x<0(tin > t0) 时 (bilinear curve(Need patch)拟合产生初始x0)
// (-4 is the value of dt, for x = 0)
else
{
x = -tdiff / (tdiff + 4.0);
double w = x + c0 * Math.Sqrt(2.0 * (1.0 - thr2));
if (w < 0.0)
x = x - Math.Sqrt(d8rt(-w)) * (x + Math.Sqrt(tdiff / (tdiff + 1.5 * t0)));
w = 4.0 / (4.0 + tdiff);
x = x * (1.0 + x * (c1 * w - c2 * x * Math.Sqrt(w)));
}
}
#endregion#region 多圈转移,x的初值(|x|<1,仅有椭圆情形)
else
{
//首先求出m圈转移中对应最小时间Tmin的Xm//xm初值的选取
xm = 1.0 / (1.5 * (m + 0.5) * Math.PI);
if (thr2 < 0.5) xm = d8rt(2.0 * thr2) * xm;
if (thr2 > 0.5) xm = (2.0 - d8rt(2.0 - 2.0 * thr2)) * xm;
#region 在12个循环内迭代找到tmin,及xm (Halley's method for iteration)
d2t = 0;
int i = 0;
while (i++ < 12)
{
//求解xm对应的tmin及其1-3阶导数
TLAMB(m, q, qsqfm1, xm, 3, out tmin, out dt, out d2t, out d3t);
if (d2t == 0) break;
//当q=1时,d2t=0,此时xm=0,停止迭代double xmold = xm;
xm = xm - dt * d2t / (d2t * d2t - dt * d3t / 2.0);
//Halley迭代,更新xm//若xm相对改变小于tol,则认为找到Xm,停止迭代
if (Math.Abs(xmold / xm - 1.0) <= tol) break;
}//找不到xm,令n=-1,然后返回!( 此种情况不应该发生)
if (i > 11)
{
n = -1;
return;
}
#endregiontdiffm = tin - tmin;
//当 tin < tmin 时，无解(N=0),程序退出
if (tdiffm < 0.0)
{
n = 0;
return;
}
//当 tin = tmin 时，仅有1解(N=1),程序退出
if (tdiffm == 0.0)
{
x = xm;
n = 1;
return;
}
//当 tin > tmin 时，有两解n = 2;
if (d2t == 0) d2t = 6.0 * m * Math.PI;
d2t2 = d2t / 2.0;
//T的二阶导数(x=xm时)/2#region 求出x>xm时的初值xpl
x = Math.Sqrt(tdiffm / (d2t / 2.0 + tdiffm / (1.0 - xm) / (1.0 - xm)));
double w = xm + x;
w = w * 4.0 / (4.0 + tdiffm) + (1.0 - w) * (1.0 - w);
x = x * (1.0 - (1.0 + m + c41 * (thr2 - 0.5)) / (1.0 + c3 * m) * x * (c1 * w + c2 * x * Math.Sqrt(w))) + xm;
xpl = x;
//若x>1,则x>xm时没有解
if (x >= 1.0)
n = 1;
#endregion#region 求出x0,x=0时的T
TLAMB(m, q, qsqfm1, 0, 0, out t0, out dt, out d2t, out d3t);
double tdiff0 = t0 - tmin;
double tdiff = tin - t0;
//tmin < tin t0 的情形
else
{
x = -tdiff / (tdiff + 4.0);
w = x + c0 * Math.Sqrt(2.0 * (1.0 - thr2));
if (w < 0.0) x = x - Math.Sqrt(d8rt(-w)) * (x + Math.Sqrt(tdiff / (tdiff + 1.5 * t0)));
w = 4.0 / (4.0 + tdiff);
x = x * (1.0 + (1.0 + m + c42 * (thr2 - 0.5)) / (1.0 + c3 * m) * x * (c1 * w - c2 * x * Math.Sqrt(w)));
}if (x <= -1.0)//若x<-1,则x 
