램버트 문제 (Lambert’s problem)의 해

램버트 문제(Lambert's problem)는 이름에서 알 수 있듯이 18세기 수학자 Johann Heinrich Lambert에 의해 처음 제기된 문제다. 그는 이 문제를 해결하기 위해 램버트의 정리(https://pasus.tistory.com/315)를 고안하였다. 참고로 램버트 정리의 증명은 라그랑지(Lagrange)가 하였고 램버트 문제의 해를 처음으로 구한 사람은 가우스(Gauss)였다.

램버트 문제는 기본적으로 두 지점 사이를 비행하는 데 걸리는 시간이 주어졌을 때, 두 지점 사이를 연결하는 궤도를 찾는 것이다. 수학적으로 램버트 문제는 이체문제(two-body problem)에서 유도된 기본 궤도 미분 방정식에 대한 2점 경계값 문제(TPBVP, two-point boundary value problem)로 표현할 수 있다.

 

$$\begin{array}{r}\text{(1)}& & \frac{{}^i{d}^2\overrightarrow{r}}{d{t}^2}+\frac{\mu }{r^3}\overrightarrow{r}=0\\ & \overrightarrow{r}(t_1)={\overrightarrow{r}}_1,\overrightarrow{r}(t_2)={\overrightarrow{r}}_2\end{array}$$

 

여기서 $\mu$ 는 중력 파라미터, $\overrightarrow{r}$ 은 관성 좌표계의 원점에서 궤도 운동을 하는 물체까지의 위치 벡터다.

 

 

램버트 문제는 우주비행역학, 우주궤도역학 및 천체역학 분야에서 다루는 근본적인 문제에 관한 것이기 때문에 광범위하게 응용되고 있다. 응용 분야를 3가지 영역에서 정리하면 다음과 같다.

첫째, 위치벡터 ${\overrightarrow{r}}_1,{\overrightarrow{r}}_2$ 와 비행시간 $\mathrm{\Delta }t=t_2-t_1$ 을 궤도가 알려져 있지 않은 우주 물체의 관측값으로 해석하고, 그 물체의 궤도를 결정 (orbit determination)하는데 이용한다.

둘째, 위치벡터 ${\overrightarrow{r}}_1$ 을 현재의 궤도, ${\overrightarrow{r}}_2$ 를 최종 궤도상의 위치, 그리고 비행시간 $\mathrm{\Delta }t$ 를 설계변수로 해석하고, 우주비행체의 궤도를 변경 (orbit transfer)하는데 이용한다.

셋째, 위치벡터 ${\overrightarrow{r}}_1$ 을 인터셉터의 현재 위치, ${\overrightarrow{r}}_2$ 를 타깃의 위치, 그리고 $\mathrm{\Delta }t$ 를 비행시간으로 해석하고, 위성 요격 및 랑데뷰, 탄도 미사일 요격등에 이용한다.

실제로 램버트 문제는 보이저 탐사선을 비롯한 다양한 우주선의 궤도를 계산하는 데 사용되었으며 아폴로가 달 표면에 착륙하는데 유도 알고리즘으로도 사용되었다.

램버트 문제를 풀기 위한 알고리즘은 그동안 여러가지가 제시되어 왔지만, 대부분의 알고리즘은 반복적인 계산 방식을 사용한다. 반복적 계산 알고리즘에서는 특정 궤도에서 두 지점 사이를 비행하는 데 필요한 시간을 계산하고, 이를 주어진 비행시간과 비교한다. 비행시간이 맞지 않으면 궤도를 수정하여 비행시간을 다시 계산하고 이 과정을 비행시간이 일치할 때까지 반복하는 것이다. 다양한 알고리즘이 제시되어 왔음에도 불구하고 이심율이 매우 큰 타원궤도과 쌍곡선궤도, 비행각이 거의 180도에 달하는 경우, 그리고 출발 지점에서 수차례의 궤도 공전 (multiple revolution) 후에 목표지점으로 가는 경우 등에서 계산 상의 문제가 발생하기 때문에 아직 알고리즘이 완벽하다고 할 수는 없다.

램버트 문제를 풀기 위한 알고리즘은 기하학적 해석에 바탕을 둔 방법이 대부분이다. 물론 기하학적 방법외에도 2점 경계값 문제를 직접 푸는 방법도 있다 (https://pasus.tistory.com/297).

 

 

램버트 문제는 두 지점의 위치벡터 ${\overrightarrow{r}}_1$ 과 ${\overrightarrow{r}}_2$, 두 지점을 비행하는데 필요한 시간 $\mathrm{\Delta }t=t_2-t_1$, 그리고 비행방향이 주어졌을 때, 두 지점에서의 속도벡터 ${\overrightarrow{v}}_1$ 과 ${\overrightarrow{v}}_2$ 를 구하는 문제로 볼 수 있다. 왜냐하면 어느 지점에서 위치벡터와 속도벡터를 알면 그 궤도상의 모든 지점에서 위치벡터와 속도벡터를 시간의 함수로 계산할 수 있을 뿐만 아니라 궤도의 크기, 모양, 자세 등 궤도에 관한 모든 것을 알 수 있기 때문이다. 여기서 비행방향은 우주비행체가 ${\overrightarrow{r}}_1$ 에서 ${\overrightarrow{r}}_2$ 로 갈 때 작은각 경로(short way)로 가는지, 큰각 경로(long way)로 가는지를 의미한다.

 

 

만약 두 위치벡터가 평행하고 방향이 반대 ($\mathrm{\Delta }\theta ={180}^0$ ) 라면 궤도면을 유일하게 정의할 수 없기 때문에 속도벡터 ${\overrightarrow{v}}_1$ 과 ${\overrightarrow{v}}_2$ 를 유일하게 결정할 수 없는 문제가 있다. 또한 만약 두 벡터가 평행하고 방향이 같다면 ($\mathrm{\Delta }\theta =0^0$ ) 비행시간이 $0$ 이거나 혹은 이체문제가 아니기 때문에 이 두 위치벡터가 평행하게 주어진 경우는 제외하기로 한다.

 

 

네 개의 벡터 ${\overrightarrow{r}}_1,{\overrightarrow{r}}_2,{\overrightarrow{v}}_1,{\overrightarrow{v}}_2$ 사이의 관계는 라그랑지 계수(Lagrange coefficients) $f$ 와 $g$ 를 이용한 표현식에 나와있다 (https://pasus.tistory.com/311).

 

$$\begin{array}{r}\text{(2)}& & {\overrightarrow{r}}_2=f{\overrightarrow{r}}_1+g{\overrightarrow{v}}_1\\ & {\overrightarrow{v}}_2=\dot{f}{\overrightarrow{r}}_1+\dot{g}{\overrightarrow{v}}_1\end{array}$$

 

위 식에 의하면 ${\overrightarrow{r}}_1,{\overrightarrow{r}}_2$ 가 주어졌을 때 ${\overrightarrow{v}}_1,{\overrightarrow{v}}_2$ 는 다음과 같이 구할 수 있다.

 

$$\begin{array}{r}\text{(3)}& {\overrightarrow{v}}_1& =\frac{1}{g}({\overrightarrow{r}}_2-f{\overrightarrow{r}}_1)\\ \text{(4)}& {\overrightarrow{v}}_2& =\dot{f}{\overrightarrow{r}}_1+\frac{\dot{g}}{g}({\overrightarrow{r}}_2-f{\overrightarrow{r}}_1)& =\frac{\dot{g}}{g}{\overrightarrow{r}}_2-\frac{f\dot{g}-\dot{f}g}{g}{\overrightarrow{r}}_1\\ & =\frac{1}{g}(\dot{g}{\overrightarrow{r}}_2-{\overrightarrow{r}}_1)\end{array}$$

 

여기서 $f\dot{g}-\dot{f}g=1$ 임을 이용하였다. 따라서 라그랑지 계수 $f,g,\dot{f},\dot{g}$ 를 구할 수 있다면 램버트 문제를 풀 수 있게 된다.

램버트 문제를 풀기 위한 기하학적 기반 알고리즘의 거의 대부분은 주어진 조건과 라그랑지 계수의 관계에서 파생되었다고 보면 된다. 대표적으로 a-반복법(iteration) 계열, p-반복법 계열, 범용변수(universal variables)를 이용한 방법 등이 있다. 이 외에 램버트 정리 (https://pasus.tistory.com/315)를 이용한 방법 등이 있다. 램버트 정리는 두 지점 사이의 비행시간은 두 지점까지의 거리의 합, 코드길이 그리고 궤도의 장반경 $a$ 만의 함수라는 것인데, 여기서 주어진 비행시간으로 장반경 $a$ 를 계산할 수 있다면 램버트 문제를 푼 것이기 때문이다.

이 중에서 범용변수 방법은 a-반복법이나 p-반복법 보다도 우수하고, 또한 램버트 정리를 이용한 방법에서 처리해야 하는 여러가지 특수한 경우를 피할 수 있기 때문에 많이 사용되고 있는 알고리즘이다. 본 게시글에서는 범용변수를 이용한 램버트 문제의 해를 유도해 보기로 한다.

우선 라그랑지 계수를 실제 비행각(true anomaly)의 변화량 $\mathrm{\Delta }\theta$ 의 함수로 표현한 식을 다시 써보자 (https://pasus.tistory.com/311).

 

$$\begin{array}{r}\text{(5)}& & f=1-\frac{\mu r_2}{h^2}(1-\cos \mathrm{\Delta }\theta )\\ & g=\frac{r_1{r}_2}{h}\sin \mathrm{\Delta }\theta \\ \\ & \dot{f}=\frac{\mu }{h}\frac{1-\cos \mathrm{\Delta }\theta }{\sin \mathrm{\Delta }\theta }[\frac{\mu }{h^2}(1-\cos \mathrm{\Delta }\theta )-\frac{1}{r_1}-\frac{1}{r_2}]\\ \\ & \dot{g}=1-\frac{\mu r_1}{h^2}(1-\cos \mathrm{\Delta }\theta )\end{array}$$

 

또한 라그랑지 계수는 범용변수 또는 범용 비행각(universal anomaly) $\chi$ 를 이용하면 다음과 같이 쓸 수 있었다 (https://pasus.tistory.com/312).

 

$$\begin{array}{r}\text{(6)}& & f=1-\frac{{\chi }^2}{r_1}C(z)\\ & g=\mathrm{\Delta }t-\frac{{\chi }^3}{\sqrt{\mu }}S(z)\\ \\ & \dot{f}=\frac{\sqrt{\mu }}{r_1{r}_2}\chi [zS(z)-1]\\ \\ & \dot{g}=1-\frac{{\chi }^2}{r_2}C(z)\end{array}$$

 

여기서 $z={\chi }^{2}a$ 이고, $S(z)$ 와 $C(z)$ 는 Stumpff 함수이다.

식 (5)와 (6)에 의하면 라그랑지 계수를 구성하는 변수는 7개로서 ${\overrightarrow{r}}_1,{\overrightarrow{r}}_2,\mathrm{\Delta }\theta ,\mathrm{\Delta }t,h,z,\chi$ 가 있다는 것을 알 수 있다. 이 중에서 4개 변수 ${\overrightarrow{r}}_1,{\overrightarrow{r}}_2,\mathrm{\Delta }\theta ,\mathrm{\Delta }t$ 는 램버트 문제에서 주어지는 것이므로 3개의 미지수 $h,z,\chi$를 계산하면 된다. 식 (5)와 (6)에 의하면 3개의 식이 있으므로 이 식을 풀면 된다. 참고로 식 (5)와 (6)에는 식이 4개가 있지만 라그랑지 계수는 $f\dot{g}-\dot{f}g=1$ 를 만족해야 하므로 독립적인 식은 3개이다. 미지수가 3개이고 식이 3개이므로 풀 수는 있지만 삼각함수를 포함하고 있기 때문에 기본적으로 반복적인 풀이 방법을 이용해야 한다.

식 (6)은 범용변수를 이용하여 라그랑지 계수를 표현했지만, 이심 비행각의 변화량 $\mathrm{\Delta }E$ 와 $\mathrm{\Delta }F$ 로 표현하는 방법도 있다. 이 경우는 ${\overrightarrow{r}}_1,{\overrightarrow{r}}_2,\mathrm{\Delta }\theta ,\mathrm{\Delta }t,p,a,\mathrm{\Delta }E$ (또는 $\mathrm{\Delta }F$ ) 등 7개의 변수가 있고 이 중 $p,a,\mathrm{\Delta }E$ (또는 $\mathrm{\Delta }F$ ) 를 계산해야 한다.

반복적 풀이 방법의 일반적인 형식은 다음과 같다. 일단 3개의 미지수 중 하나를 선택하여 그 값을 추측한 다음에 나머지 2개의 미지수를 계산한다. 그런 후 선택한 미지수를 뉴톤 방법이나 이분법(bisection method) 등을 사용하여 업데이트 한다. 이 값으로 다시 나머지 미지수 2개를 계산한다. 이런 반복 작업을 미지수가 수렴할 때까지 계속한다.

여기서 반복 작업할 미지수로 $a$ 를 선택하면 a-반복법, $p$ 를 선택하면 p-반복법이라고 한다. 범용변수 방법에서는 미지수로 $z$ 를 선택한다.

 

 

식 (5)에서 실제 비행각의 변화량 $\mathrm{\Delta }\theta$ 는 두 위치벡터의 내적으로 계산할 수 있다.

 

$$\begin{array}{}\text{(7)}& \mathrm{\Delta }\theta ={\cos }^{-1}\left(\frac{{\overrightarrow{r}}_1\cdot {\overrightarrow{r}}_2}{r_1{r}_2}\right)\end{array}$$

 

여기서 arccos 은 $0^0\sim {180}^0$ 범위를 가지므로 큰각 경로(long way)인 경우는 다음 식으로 계산해야 한다.

 

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }\theta =2\pi -{\cos }^{-1}\left(\frac{{\overrightarrow{r}}_1\cdot {\overrightarrow{r}}_2}{r_1{r}_2}\right)\end{array}$$

 

$\mathrm{\Delta }t$ 와 $\mathrm{\Delta }\theta$ 의 관계식은 식 (5)와 (6)의 $g$ 식에서 얻을 수 있다.

 

$$\begin{array}{}\text{(9)}& g=\frac{r_1{r}_2}{h}\sin \mathrm{\Delta }\theta =\mathrm{\Delta }t-\frac{{\chi }^3}{\sqrt{\mu }}S(z)\end{array}$$

 

$f$ 식은 다음과 같으므로

 

$$\begin{array}{}\text{(10)}& f=1-\frac{\mu r_2}{h^2}(1-\cos \mathrm{\Delta }\theta )=1-\frac{{\chi }^2}{r_1}C(z)\end{array}$$

 

$h$ 를 다음과 같이 계산할 수 있다.

 

$$\begin{array}{}\text{(11)}& h=\sqrt{\frac{\mu r_1{r}_2(1-\cos \mathrm{\Delta }\theta )}{{\chi }^{2}C(z)}}\end{array}$$

 

위 식을 식 (9)에 대입하면,

 

$$\begin{array}{}\text{(12)}& \sqrt{\mu }\mathrm{\Delta }t={\chi }^{3}S(z)+\chi \sqrt{C(z)}(\sin \mathrm{\Delta }\theta \sqrt{\frac{r_1{r}_2}{1-\cos \mathrm{\Delta }\theta }})\end{array}$$

 

가 된다. 위 식의 괄호항은 상수인데. 이를 $A$ 로 놓자.

 

$$\begin{array}{}\text{(13)}& A=\sin \mathrm{\Delta }\theta \sqrt{\frac{r_1{r}_2}{1-\cos \mathrm{\Delta }\theta }}\end{array}$$

 

여기서 $A>0$ 이면 작은각 경로(short way), $A<0$ 이면 큰각 경로(long way)에 해당한다. 그러면 식 (12)는 다음과 같이 간단히 쓸 수 있다.

 

$$\begin{array}{}\text{(14)}& \sqrt{\mu }\mathrm{\Delta }t={\chi }^{3}S(z)+A\chi \sqrt{C(z)}\end{array}$$

 

$z$ 와 $\chi$ 의 관계식을 구하기 위해서 식 (5)와 (6)의 $\dot{f}$ 식을 이용한다.

 

$$\begin{array}{r}\text{(15)}& \dot{f}& =\frac{\mu }{h}\frac{1-\cos \mathrm{\Delta }\theta }{\sin \mathrm{\Delta }\theta }[\frac{\mu }{h^2}(1-\cos \mathrm{\Delta }\theta )-\frac{1}{r_1}-\frac{1}{r_2}]\\ & =\frac{\sqrt{\mu }}{r_1{r}_2}\chi [zS(z)-1]\end{array}$$

 

양변에 $r_1{r}_2$ 를 곱하고 식 (11)을 대입하면 다음과 같이 된다.

 

$$\begin{array}{r}\text{(16)}& & \frac{\mu }{\sqrt{\frac{\mu r_1{r}_2(1-\cos \mathrm{\Delta }\theta )}{{\chi }^{2}C(z))}}}\frac{1-\cos \mathrm{\Delta }\theta }{\sin \mathrm{\Delta }\theta }[\frac{\mu }{h^2}(1-\cos \mathrm{\Delta }\theta )-r_1-r_2]\\ & =\sqrt{\mu }\chi [zS(z)-1]\end{array}$$

 

위 식을 정리하면,

 

$$\begin{array}{}\text{(17)}& \frac{\sqrt{1-\cos \mathrm{\Delta }\theta }}{\sqrt{r_1{r}_2}\sin \mathrm{\Delta }\theta }\sqrt{C(z)}[{\chi }^{2}C(z)-r_1-r_2]=zS(z)-1\end{array}$$

 

이 되므로 결국 다음과 같이 $z$ 와 $\chi$ 의 관계식을 얻을 수 있다.

 

$$\begin{array}{}\text{(18)}& {\chi }^{2}C(z)=r_1+r_2+A\frac{zS(z)-1}{\sqrt{C(z)}}\end{array}$$

 

위 식의 오른쪽 항은 $z$ 만의 함수이므로

 

$$\begin{array}{}\text{(19)}& y(z)=r_1+r_2+A\frac{zS(z)-1}{\sqrt{C(z)}}\end{array}$$

 

로 놓으면, 식 (18)은

 

$$\begin{array}{}\text{(20)}& \chi =\sqrt{\frac{y(z)}{C(z)}}\end{array}$$

 

가 된다. 그러면 식 (14)는 다음과 같이 $z$ 만의 함수로 표현할 수 있다.

 

$$\begin{array}{}\text{(21)}& \sqrt{\mu }\mathrm{\Delta }t={\left[\frac{y(z)}{C(z)}\right]}^{\frac{3}{2}}S(z)+A\sqrt{y(z)}\end{array}$$

 

따라서 $\mathrm{\Delta }t$ 가 주어졌을 때, 식 (21)로 $z$ 를 계산할 수 있다. $z$ 를 구하면 식 (20)으로 $\chi$, 식 (13)으로 $h$ 를 계산할 수 있으므로 라그랑지 계수를 모두 계산할 수 있다.

뉴톤-랩슨 방법(Newton-Raphson method)을 이용하기 위해서 식 (21)을 다음과 같은 함수로 놓는다.

 

$$\begin{array}{}\text{(22)}& F(z)={\left[\frac{y(z)}{C(z)}\right]}^{\frac{3}{2}}S(z)+A\sqrt{y(z)}-\sqrt{\mu }\mathrm{\Delta }t\end{array}$$

 

$|z_{i+1}-z_i|<ϵ$ 의 조건이 수렴할 때까지 다음식을 반복하여 업데이트 하기 위해서는

 

$$\begin{array}{}\text{(23)}& z_{i+1}=z_i-\frac{F(z_i)}{F^{\prime }(z_i)}\end{array}$$

 

$F^{\prime }(z_i)$ 을 계산해야 한다.

 

$$\begin{array}{r}\text{(24)}& F^{\prime }(z)& =\frac{3}{2}{\left[\frac{y(z)}{C(z)}\right]}^{\frac{1}{2}}\left[\frac{y^{\prime }(z)C(z)-C^{\prime }(z)y(z)}{C^2(z)}\right]S(z)\\ & +{\left[\frac{y(z)}{C(z)}\right]}^{\frac{3}{2}}{S}^{\prime }(z)+\frac{1}{2}A\frac{y^{\prime }(z)}{\sqrt{y(z)}}\end{array}$$

 

여기서 $y^{\prime }(z),C^{\prime }(z),S^{\prime }(z)$ 를 계산해야 하는데, 먼저 $y^{\prime }(z)$ 를 계산해 보자.

 

$$\begin{array}{r}\text{(25)}& y^{\prime }(z)& =A\frac{\sqrt{C(z)}(S(z)+z{S}^{\prime }(z))-(zS(z)-1)\frac{1}{2\sqrt{C(z)}}{C}^{\prime }(z)}{C(z)}\\ & =\frac{A}{2{\left(\sqrt{C(z)}\right)}^3}[2(S(z)+z{S}^{\prime }(z))C(z)+(1-zS(z))C^{\prime }(z)]\end{array}$$

 

Stumpff 함수의 도함수 (https://pasus.tistory.com/314)를 이용하면,

 

$$\begin{array}{r}\text{(26)}& & C^{\prime }(z)=\frac{1}{2z}(1-zS(z)-2C(z))\\ & S^{\prime }(z)=\frac{1}{2z}(C(z)-3S(z))\end{array}$$

 

식 (25)는 다음과 같이 된다.

 

$$\begin{array}{rl}\text{(27)}& y^{\prime }(z)& =\frac{A}{2{\left(\sqrt{C(z)}\right)}^3}[2(S(z)+z(\frac{1}{2z}(C(z)-3S(z))))C(z)& +(1-zS(z))(\frac{1}{2z}(1-zS(z)-2C(z)))]\\ \\ & =\frac{A}{2{\left(\sqrt{C(z)}\right)}^3}[C^2(z)+\frac{1}{2z}-S(z)-\frac{C(z)}{z}+\frac{z{S}^2(z)}{2}]\end{array}$$

 

여기서 Stumpff 함수의 정의,

 

$$\begin{array}{r}\text{(28)}& & C(z)=\frac{1-\cos \sqrt{z}}{z}\\ & S(z)=\frac{\sqrt{z}-\sin \sqrt{z}}{{\left(\sqrt{z}\right)}^3}\end{array}$$

 

에 의하면, 식 (27)에서

 

$$\begin{array}{r}\text{(29)}& & z{S}^2(z)=\frac{{(\sqrt{z}-\sin \sqrt{z})}^2}{z^2}=\frac{1}{z^2}(z-2\sqrt{z}\sin \sqrt{z}+{\sin }^2\sqrt{z})\\ & C^2(z)=\frac{{(1-\cos \sqrt{z})}^2}{z^2}=\frac{1}{z^2}(1-2\cos \sqrt{z}+{\cos }^2\sqrt{z})\\ \\ & z{S}^2(z)+C^2(z)=2(\frac{1}{z^2}+\frac{1}{2z}-\frac{\sin \sqrt{z}}{(\sqrt{z}{)}^3}-\frac{\cos \sqrt{z}}{z^2})\\ \\ & =2(-\frac{1}{2z}+S(z)+\frac{C(z)}{z})\end{array}$$

 

이 성립하므로 식 (27)의 대괄호항은 다음과 같이 된다.

 

$$\begin{array}{rl}& C^2(z)+\frac{1}{2z}-S(z)-\frac{C(z)}{z}+\frac{z{S}^2(z)}{2}\\ \\ & =C^2(z)-\frac{C^2(z)}{2}\\ \\ & =\frac{C^2(z)}{2}\end{array}$$

 

따라서 y^' (z)은 다음과 같이 된다.

 

$$\begin{array}{r}\text{(30)}& y^{\prime }(z)& =\frac{A}{2{\left(\sqrt{C(z)}\right)}^3}\frac{C^2(z)}{2}\\ & =\frac{A}{4}\sqrt{C(z)}\end{array}$$

 

$z$ 의 정의에 의하면, $z<0$ 이면 쌍곡선, $z=0$ 이면 포물선, $z>0$ 이면 타원궤도다. $z=0$ 인 경우 식 (26)과 (28)에서 $z$ 가 분모에 있기 떼문에 계산상의 문제를 야기할 수 있는데 이 때는 다음 식으로 계산하면 된다. $z=0$ 일 때 Stumpff 함수는 다음과 같이 계산한다 (https://pasus.tistory.com/310).

 

$$\begin{array}{}\text{(31)}& C(0)=\frac{1}{2},S(0)=\frac{1}{6}\end{array}$$

 

Stumpff 함수의 정의를 사용하여 식 (26)을 전개하면 다음과 같다.

 

$$\begin{array}{rl}\text{(32)}& C^{\prime }(z)& =\frac{1}{2z}(1-zS(z)-2C(z))& =\frac{1}{2z}(1-z(\frac{1}{6}-\frac{z}{120}+\cdots )-2(\frac{1}{2}-\frac{z}{24}+\cdots ))\\ & =-\frac{3}{72}+\frac{z}{240}+\cdots \\ \\ S^{\prime }(z)& =\frac{1}{2z}(C(z)-3S(z))\\ & =\frac{1}{2z}((\frac{1}{2}-\frac{z}{24}+\cdots )-3(\frac{1}{6}-\frac{z}{120}+\cdots ))\\ & =-\frac{1}{80}+\cdots \end{array}$$

 

따라서 $z=0$ 일 때 Stumpff 함수의 도함수 값은 다음과 같다.

 

$$\begin{array}{}\text{(33)}& C^{\prime }(0)=-\frac{1}{24},S^{\prime }(0)=-\frac{1}{80}\end{array}$$

 

식 (28), (30), (31), (32), (33)을 식 (24)에 대입하면 뉴톤-랩슨 방법으로 $z$ 를 계산할 수 있다. 그런 후 식 (19)와 (20)을 식 (6)에 대입하면 라그랑지 계수를 모두 계산할 수 있다.

 

$$\begin{array}{r}\text{(34)}& & f=1-\frac{{\left[\sqrt{\frac{y(z)}{C(z)}}\right]}^2}{r_1}C(z)=1-\frac{y(z)}{r_1}\\ & g=\frac{1}{\sqrt{\mu }}[{\left[\frac{y(z)}{C(z)}\right]}^{\frac{3}{2}}S(z)+A\sqrt{y(z)}]-\frac{1}{\sqrt{\mu }}{\left[\frac{y(z)}{C(z)}\right]}^{\frac{3}{2}}S(z)\\ & =A\sqrt{\frac{y(z)}{\mu }}\\ \\ & \dot{f}=\frac{\sqrt{\mu }}{r_1{r}_2}\sqrt{\frac{y(z)}{C(z)}}[zS(z)-1]\\ \\ & \dot{g}=1-\frac{{\left[\sqrt{\frac{y(z)}{C(z)}}\right]}^2}{r_2}C(z)=1-\frac{y(z)}{r_2}\end{array}$$

 

라그랑지 계수를 계산한 후 식 (3)과 (4)에 대입하면 속도벡터 ${\overrightarrow{v}}_1,{\overrightarrow{v}}_2$ 를 구할 수 있다.

지금까지 논의한 내용을 정리하면 다음과 같다.

${\overrightarrow{r}}_1,{\overrightarrow{r}}_2,\mathrm{\Delta }t$ 와 비행방향(작은각 경로, 큰각 경로)이 주어졌을 때, 램버트 문제의 해는 다음과 같은 순서로 구할 수 있다.

[1] ${\overrightarrow{r}}_1,{\overrightarrow{r}}_2$ 에서 $r_1,r_2$ 를 계산한다.
[2] 비행방향에 따라 식 (7) 또는 (8)로 $\mathrm{\Delta }\theta$ 를 계산한다.
[3] 식 (13)으로 $A$ 를 계산한다.
[4] $z$ 를 계산한다.
     1. 식 (28), (31), (19), (22)로 $C(z),S(z),y(z),F(z)$ 를 계산한다.
     2. 식 (26), (33), (30), (24)로 $C^{\prime }(z),S^{\prime }(z),y^{\prime }(z),F^{\prime }(z)$ 를 계산한다.
     3. 식 (23)으로 $z$ 를 수렴할 때까지 업데이트 한다.
[5] 식 (34)로 $f,g,\dot{f},\dot{g}$ 을 계산한다.
[6] 식 (3)과 (4)로 ${\overrightarrow{v}}_1,{\overrightarrow{v}}_2$ 를 계산한다.

 

 

다음은 위 알고리즘을 매트랩 코드로 구현한 것이다.

 

lambert.m

 

function [v1_vec, v2_vec] = lambert(r1_vec, r2_vec, tof, dm)

% [v1_vec, v2_vec] = lambert(r1_vec, r2_vec, tof, dm)
%        given r1_vec, r2_vec, time-of-flight and direction
%        compute v1_vec and v2_vec
%       
% Input  r1_vec: initial position vector in km
%        r2_vec: final position vector in km
%        tof: time of flight in sec
%        dm: +1 short way (default), -1 long way
% Output v1_vec: initial velocity vector in km/sec
%        v2_vec: final velocity vector in km/sec
%
% all in inertial frame
%
% coded by st.watermelon

eps = 1e-8;
if nargin < 4
    dm = 1;
end

mu_e=398600; % gravitational parameter km^3/s^2

% step 1 : compute r1 and r2
r1 = norm(r1_vec);
r2 = norm(r2_vec);

% step 2 : compute d_theta
d_the = acos(dot(r1_vec,r2_vec)/(r1*r2));
if dm ~= 1
    d_the = 2*pi - d_the;
end

% step 3 : compute A
A = sin(d_the)*sqrt(r1*r2/(1-cos(d_the)));

% step 4 : compute z  
z = 0; % guess initial z

for ii = 1:100 
    % compute C(z) and S(z)
    [Cz, Sz] = stumpff(z);

    % compute y(z)
    yz = r1 + r2 + A*(z*Sz-1)/sqrt(Cz);

    % compute F(z)
    Fz = (yz/Cz)^(3/2)*Sz + A*sqrt(yz) - sqrt(mu_e)*tof;

    % compute derivatives
    yz_prime = A/4*sqrt(Cz);
    if (z<-eps) || (z>eps)
        Cz_prime = (1-z*Sz-2*Cz)/(2*z);
        Sz_prime = (Cz-3*Sz)/(2*z);
    else
        Cz_prime = -1/24;
        Sz_prime = -1/80;
    end

    Fz_prime = 3/2*sqrt(yz/Cz)*(yz_prime*Cz-Cz_prime*yz)/(Cz*Cz)*Sz ...
        + (yz/Cz)^(3/2)*Sz_prime + 1/2*A*yz_prime/sqrt(yz);

    % Newton-Raphson
    z_next = z - Fz/Fz_prime;
    if abs(z_next-z) < eps
        break;
    else
        z = z_next;
    end

end


% step 5 : compute Lagrange coefficients
%    compute C(z) and S(z)
[Cz, Sz] = stumpff(z);

%    compute y(z)
yz = r1 + r2 + A*(z*Sz-1)/sqrt(Cz);

%    compute f, g, g_dot
f = 1 - yz/r1;
g = A*sqrt(yz/mu_e);
g_dot = 1 - yz/r2;

% step 6 : compute v1_vec, v2_vec
v1_vec = (r2_vec-f*r1_vec)/g;
v2_vec = (g_dot*r2_vec-r1_vec)/g;

 

stumpff.m

 

function [Cz, Sz] = stumpff(z)

% [Cz, Sz] = stumpff(z)
%
% Stumpff functions
% coded by st.watermelon

eps = 1e-8;
if z > eps
    Cz = (1- cos(sqrt(z)))/z;
    Sz = (sqrt(z) - sin(sqrt(z)))/(sqrt(z))^3;

elseif z < -eps
    Cz = (cosh(sqrt(-z))-1)./(-z);
    Sz = (sinh(sqrt(-z))-sqrt(-z))/(sqrt(-z))^3;

else
    Cz = 1/2;
    Sz = 1/6;
    
end