램버트 문제 (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{align} & \frac{ ^i d^2 \vec{r}}{dt^2 }+ \frac{\mu}{r^3} \vec{r}=0 \tag{1} \\ \\ & \vec{r}(t_1 )= \vec{r}_1, \ \ \ \ \ \vec{r}(t_2 )= \vec{r}_2 \end{align} \]

 

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

 

 

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

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

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

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

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

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

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

 

 

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

 

 

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

 

 

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

 

\[ \begin{align} & \vec{r}_2=f \vec{r}_1+g \vec{v}_1 \tag{2} \\ \\ & \vec{v}_2= \dot{f} \vec{r}_1+ \dot{g} \vec{v}_1 \end{align} \]

 

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

 

\[ \begin{align} \vec{v}_1 &= \frac{1}{g} (\vec{r}_2-f \vec{r}_1) \tag{3} \\ \\ \vec{v}_2 &= \dot{f} \vec{r}_1+ \frac{ \dot{g}}{g} ( \vec{r}_2-f \vec{r}_1 ) \tag{4} \\ &= \frac{ \dot{g}}{g} \vec{r}_2- \frac{f \dot{g}- \dot{f}g}{g} \vec{r}_1 \\ & = \frac{1}{g} ( \dot{g} \vec{r}_2- \vec{r}_1) \end{align} \]

 

여기서 \(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)의 변화량 \(\Delta \theta\) 의 함수로 표현한 식을 다시 써보자 (https://pasus.tistory.com/311).

 

\[ \begin{align} & f= 1-\frac{\mu r_2}{h^2} (1- \cos \Delta \theta ) \tag{5} \\ \\ & g= \frac{r_1 r_2}{h} \sin \Delta \theta \\ \\ & \dot{f} = \frac{\mu}{h} \frac{1- \cos \Delta \theta}{ \sin⁡ \Delta \theta} \left[ \frac{\mu}{h^2} (1-\cos \Delta \theta )- \frac{1}{r_1} - \frac{1}{r_2} \right] \\ \\ & \dot{g} =1- \frac{\mu r_1}{h^2} (1- \cos \Delta \theta ) \end{align} \]

 

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

 

\[ \begin{align} & f=1- \frac{\chi^2}{r_1} C(z) \tag{6} \\ \\ & g= \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{align} \]

 

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

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

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

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

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

 

 

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

 

\[ \begin{align} \Delta \theta = \cos^{-1} \left( \frac{ \vec{r}_1 \cdot \vec{r}_2}{ r_1 r_2 } \right) \tag{7} \end{align} \]

 

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

 

\[ \begin{align} \Delta \theta = 2\pi - \cos^{-1} \left( \frac{ \vec{r}_1 \cdot \vec{r}_2}{ r_1 r_2 } \right) \tag{8} \end{align} \]

 

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

 

\[ \begin{align} g= \frac{r_1 r_2}{h} \sin \Delta \theta = \Delta t- \frac{\chi^3}{ \sqrt{\mu} } S(z) \tag{9} \end{align} \]

 

\(f\) 식은 다음과 같으므로

 

\[ \begin{align} f=1-\frac{\mu r_2}{h^2} (1- \cos \Delta \theta )=1- \frac{\chi^2}{r_1} C(z) \tag{10} \end{align} \]

 

\(h\) 를 다음과 같이 계산할 수 있다.

 

\[ \begin{align} h= \sqrt{ \frac{ \mu r_1 r_2 (1- \cos \Delta \theta )}{ \chi^2 C(z) } } \tag{11} \end{align} \]

 

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

 

\[ \begin{align} \sqrt{\mu} \Delta t= \chi^3 S(z)+ \chi \sqrt{C(z)} \left( \sin \Delta \theta \sqrt{ \frac{r_1 r_2}{1-\cos \Delta \theta} } \ \right) \tag{12} \end{align} \]

 

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

 

\[ \begin{align} A= \sin \Delta \theta \sqrt{ \frac{r_1 r_2}{1-\cos \Delta \theta} } \tag{13} \end{align} \]

 

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

 

\[ \begin{align} \sqrt{\mu} \Delta t= \chi^3 S(z)+ A \chi \sqrt{C(z)} \tag{14} \end{align} \]

 

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

 

\[ \begin{align} \dot{f} &= \frac{\mu}{h} \frac{1- \cos \Delta \theta}{ \sin \Delta \theta} \left[ \frac{\mu}{h^2} (1-\cos \Delta \theta )- \frac{1}{r_1} - \frac{1}{r_2} \right] \tag{15} \\ \\ &= \frac{\sqrt{\mu}}{r_1 r_2 } \chi [zS(z)-1] \end{align} \]

 

양변에 \(r_1 r_2\) 를 곱하고 식 (11)을 대입하면 다음과 같이 된다.

 

\[ \begin{align} & \frac{\mu}{ \sqrt{ \frac{ \mu r_1 r_2 (1- \cos \Delta \theta )}{\chi^2 C(z)) } } } \frac{1-\cos \Delta \theta}{ \sin \Delta \theta} \left[ \frac{\mu}{h^2} (1-\cos \Delta \theta )-r_1-r_2 \right] \tag{16} \\ \\ & \ \ \ \ \ = \sqrt{\mu} \chi [zS(z)-1] \end{align} \]

 

위 식을 정리하면,

 

\[ \begin{align} \frac{ \sqrt{1- \cos \Delta \theta }}{ \sqrt{r_1 r_2 } \sin \Delta \theta } \sqrt{C(z)} \ \left[ \chi^2 C(z)-r_1-r_2 \right]=zS(z)-1 \tag{17} \end{align} \]

 

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

 

\[ \begin{align} \chi^2 C(z)=r_1+r_2+A \frac{zS(z)-1}{\sqrt{C(z)} } \tag{18} \end{align} \]

 

위 식의 오른쪽 항은 \(z\) 만의 함수이므로

 

\[ \begin{align} y(z)=r_1+r_2+A \frac{zS(z)-1}{\sqrt{C(z)} } \tag{19} \end{align} \]

 

로 놓으면, 식 (18)은

 

\[ \begin{align} \chi =\sqrt{ \frac{y(z)}{C(z)} } \tag{20} \end{align} \]

 

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

 

\[ \begin{align} \sqrt{\mu} \Delta t= \left[ \frac{y(z)}{C(z)} \right]^{\frac{3}{2}} S(z)+A \sqrt{y(z)} \tag{21} \end{align} \]

 

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

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

 

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

 

\(\vert z_{i+1}-z_i \vert \lt \epsilon\) 의 조건이 수렴할 때까지 다음식을 반복하여 업데이트 하기 위해서는

 

\[ \begin{align} z_{i+1}=z_i- \frac{F(z_i )}{F^\prime (z_i ) } \tag{23} \end{align} \]

 

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

 

\[ \begin{align} 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) \tag{24} \\ \\ & \ \ \ + \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{align} \]

 

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

 

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

 

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

 

\[ \begin{align} & C^\prime(z)= \frac{1}{2z} (1-zS(z)-2C(z)) \tag{26} \\ \\ & S^\prime (z)= \frac{1}{2z} (C(z)-3S(z)) \end{align} \]

 

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

 

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

 

여기서 Stumpff 함수의 정의,

 

\[ \begin{align} & C(z)= \frac{1-\cos \sqrt{z}}{z} \tag{28} \\ \\ & S(z)= \frac{\sqrt{z}- \sin \sqrt{z}}{ \left( \sqrt{z} \right)^3} \end{align} \]

 

에 의하면, 식 (27)에서

 

\[ \begin{align} & zS^2 (z)= \frac{ \left(\sqrt{z}- \sin \sqrt{z} \right)^2}{z^2} = \frac{1}{z^2} (z-2 \sqrt{z} \sin \sqrt{z}+ \sin^2 \sqrt{z}) \tag{29} \\ \\ & C^2 (z)= \frac{ \left( 1- \cos \sqrt{z} \right)^2}{z^2} = \frac{1}{z^2} (1-2 \cos \sqrt{z}+ \cos^2 \sqrt{z}) \\ \\ & zS^2 (z)+C^2 (z)=2 \left( \frac{1}{z^2} + \frac{1}{2z}- \frac{ \sin \sqrt{z} }{(\sqrt{z})^3} - \frac{\cos \sqrt{z}}{z^2} \right) \\ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = 2 \left( -\frac{1}{2z}+ S(z)+ \frac{C(z)}{z} \right) \end{align} \]

 

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

 

\[ \begin{align} & C^2 (z)+ \frac{1}{2z}-S(z)-\frac{C(z)}{z}+ \frac{zS^2 (z)}{2} \\ \\ & \ \ \ \ \ =C^2 (z)-\frac{C^2 (z)}{2} \\ \\ & \ \ \ \ \ =\frac{C^2 (z)}{2} \end{align} \]

 

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

 

\[ \begin{align} y^\prime (z) & = \frac{A}{2 \left( \sqrt{C(z) } \right)^3 } \frac{C^2 (z)}{2} \tag{30} \\ \\ & = \frac{A}{4} \sqrt{C(z)} \end{align} \]

 

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

 

\[ \begin{align} C(0)= \frac{1}{2}, \ \ \ \ \ S(0)=\frac{1}{6} \tag{31} \end{align} \]

 

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

 

\[ \begin{align} C^\prime (z) &= \frac{1}{2z} (1-zS(z)-2C(z)) \tag{32} \\ &= \frac{1}{2z} \left( 1-z \left( \frac{1}{6}- \frac{z}{120} + \cdots \right)- 2 \left( \frac{1}{2}- \frac{z}{24} + \cdots \right) \right) \\ & =- \frac{3}{72}+ \frac{z}{240} + \cdots \\ \\ S^\prime (z) & = \frac{1}{2z} (C(z)-3S(z)) \\ &= \frac{1}{2z} \left( \left( \frac{1}{2}- \frac{z}{24} + \cdots \right)-3 \left( \frac{1}{6}- \frac{z}{120} + \cdots \right) \right) \\ &= -\frac{1}{80} + \cdots \end{align} \]

 

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

 

\[ \begin{align} C^\prime (0)=- \frac{1}{24}, \ \ \ \ \ S^\prime (0) =- \frac{1}{80} \tag{33} \end{align} \]

 

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

 

\[ \begin{align} & f=1- \frac{ \left[ \sqrt{ \frac{y(z)}{C(z)} } \right]^2}{r_1} C(z)=1- \frac{y(z)}{r_1} \tag{34} \\ \\ & g= \frac{1}{\sqrt{\mu}} \left[ \left[ \frac{y(z)}{C(z)} \right]^{\frac{3}{2}} S(z)+A \sqrt{y(z)} \right] - \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{align} \]

 

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

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

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

[1] \( \vec{r}_1, \ \vec{r}_2\) 에서 \(r_1, \ r_2\) 를 계산한다.
[2] 비행방향에 따라 식 (7) 또는 (8)로 \(\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)로 \(\vec{v}_1, \ \vec{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