Lambert's problem은 초기 위치 $r_1$ 에서 목표 위치 $r_2$로 전이하기 위한 초기속도를 구하는 문제이다. 전이 궤도의 방정식 혹은 궤도 요소를 구하면, 그 위치에서의 속도를 구할 수 있다. 그렇기에 Lambert's problem의 핵심은 전이 궤도의 궤도 요소를 구하는 것이 된다.
개요
Lambert's problem의 인풋은 아래와 같다.
- 초기 위치 $\vec{r_1}$
- 목표 위치 $\vec{r_2}$
- 목표 전이 시간 $TOF$
이 정보로부터 궤도의 6요소를 유도해 보자. 궤도의 6요소는 $a, e, i, \omega, \Omega, \nu, $이다.
- Inclination, raan
두 위치벡터를 외적하면 각 운동량 벡터의 방향을 알 수 있다. 이로부터 궤도의 $i, \Omega$를 찾을 수 있다.
- 장반경
장반경은 잠시 보류하자. 장반경을 구하는 것이 Lambert's problem의 핵심이다. 아래에서 자세하게 다룰 것이다. 여기서는 장반경을 구했다고 가정하자.
- 이심률, True anomaly
궤도 방정식 $ r_1 = \frac{a(1-e^2)}{1+ecos(\nu)} \\ r_2 = \frac{a(1-e^2)}{1+ecos(\nu+\Delta \nu)}$ 에서 $r_1, r_2, a, \Delta \nu$를 알고있으므로 식 2개, 미지수 2개에서 $e, \nu$를 찾을 수 있다.
- Argument of perigee
$\vec{r_1}$과 이때의 각도 $\nu$를 알고있으므로 $\omega$를 구할 수 있다.
이런 방법으로 주어진 조건(인풋)으로부터 궤도의 6요소를 찾아낼 수 있다. 자세한 구현 방법은 지금부터 다룰 것이다.
Recap Kepler's Equation (elliptic)
$$ r = a(1 - ecosE) $$
$$ n(t-T_0) = E - esinE + 2M\pi $$
$$ rcos\nu = a(cosE - e) $$
$$ rsin\nu = a\sqrt{1-e^2}sinE $$
Function of semi-axis
위 케플러 방정식으로부터 궤도 장반경을 구해보자.
$$n = \sqrt{\frac{\mu}{a^3}}$$
$$\sqrt\mu(t_1-T_0) = a^{3/2}(E_1-esinE_1+2M_1\pi)$$
$$\sqrt\mu(t_2-T_0) = a^{3/2}(E_2-esinE_2+2M_2\pi)$$
$$\downarrow$$
$$\sqrt\mu(t_2-t_1) = a^{3/2}(E_2-E_1-esinE_2+esinE_1+2M\pi) \cdots (1)$$
$$M = M_2 - M_1 : \textit{Complete rovolutions made during the transfer}$$
수식 (1)은 전이 시간 $t_2 - t_1$과 궤도 장반경 $a$의 관계를 나타내는 식이다. $M$은 전이할 대 궤도를 몇바퀴 회전하였는지 나타내는 순 회전수로, 입력값이다.
수식 안에 $E_1, E_2$과 같은 문자들이 있기에 아직 풀 수 는 없다. 풀기 위해서는 2개의 추가적인 관계식을 찾아야 한다.
- Define two new quantities
$$\psi = \frac{E_2 - E_1}{2}, \qquad cos\phi = ecos\frac{E_2+E_1}{2}$$
여기서 기하학적인 부분을 생각해보자.
$$E_2-E_1$$은 원으로 확장한(케플러 원) 두 위치의 사잇각이다. 그러므로 $E_2-E_1 \in [0, 2\pi], \psi \in [0,\pi]$일 것이다.
$\phi$는 ambiguity를 피하기 위해 $\phi \in [0, \pi]$인 변수라고 정하자. 범위는 매우 중요한 조건이 되므로 기억하자.
그러면 수식 (1)은 아래와 같이 정리할 수 있다.
$$\sqrt\mu(t_2-t_1) = a^{3/2}(\phi - e(sinE_2 - sinE_1)+2M\pi) $$
$$ \small{(\textit{Using,} \quad sin(\alpha+\beta)-sin(\alpha-\beta) = 2cos\alpha sin\beta)} $$
$$ = a^{3/2}(\phi - 2ecos\frac{E_2+E_1}{2}sin\frac{E_2-E_1}{2}+2M\pi))$$
$$\therefore \sqrt\mu(t_2-t_1) =2a^{3/2}(\phi-cos\phi sin\psi +M\pi) \cdots (2)$$
- Geometry relation 1
이어서 기하학적인 관계식을 유도하자
$$r_1 = a(1-ecosE_1), \quad r_2=a(1-ecosE_2)$$
$$r_1 + r_2 = a(2 - e(cosE_2+cosE_1))$$
$$\small{(\textit{Using,} \quad cos(\alpha+\beta)+cos(\alpha-\beta) = 2cos\alpha cos\beta)}$$
$$ = 2a(1 - cos\phi cos\psi) \cdots (3) $$
- Geometry relation 2
빗변의 길이를 c라고 하자.
$$c^2 = |\vec{r_2} - \vec{r_1}|^2 = (r_2cosf_2 - r_1cosf_1)^2 + (r_2sinf_2 - r_1sinf_1)^2 $$
$$ = a^2(cosE_2-cosE_1)^2 + a^2(1-e^2)(sinE_2-sinE_1)^2 = a^2((cosE_2-cosE_1)^2+(1-e^2)(sinE_2-sinE_1)^2))$$
$$ = a^2(4sin^2\frac{E_2+E_1}{2}sin^2\frac{E_2-E_1}{2}+(1-e^2)4cos^2\frac{E_2+E_1}{2}sin^2\frac{E_2-E_1}{2} ) $$
$$ = 4a^2sin^2\psi (sin^2\frac{E_2+E_1}{2}+(1-e^2)cos^2\frac{E_2+E_1}{2}) $$
$$ = 4a^2sin^2\psi (1-cos^2\phi) = 4a^2sin^2\psi sin^2\phi$$
$$\therefore c = 2asin\psi sin\phi \cdots (4) $$
또한
$$c^2 = r_1^2 + r_2^2 - 2r_1r_2cos\Delta\nu \cdots (5) $$
이다.
식 (2),(3),(4)를 보자. $t_2-t_1$는 $a, \psi, \phi$에 대한 함수이고, $ \psi, \phi $는 $r_1+r_1, c$에 대한 함수 이므로 결국 $t_2-t_1$ 는 $a, r_1+r_2, c$에 대한 함수이다. 그런데, 이미 $r_1,r_2,c$는 알고 있다. 즉 $t_2-t_1$는 오로지 $a$에 대한 함수가 된다. 따라서 이제 전이 궤도의 장반경 $a$를 구할 수 있다.
Solution. Find semi-axis
궤도 장반경을 구하기 위한 모든 관계식을 찾았으니, 이제 구해보자. 쉽게 구하기 위해 새로운 변수를 도입한다.
- Define three new quantities ($\alpha, \beta, s$)
$$\alpha = \phi + \psi, \qquad \beta=\phi-\psi$$
$$\alpha \in [0,\pi], \beta \in [-pi, pi]$$
$$s = \frac{r_1+r_2+c}{2} $$
$$\downarrow$$
$$\begin{cases} (2) : \sqrt\mu(t_2-t_1) = a^{3/2}((\alpha-sin\alpha) - (\beta - sin\beta) + 2M\pi) \\
(3)+(4) : \frac{s}{2a} = sin^2\frac{\alpha}{2} \\
(3) -(4) : \frac{s-c}{2a} = sin^2\frac{\beta}{2}
\end{cases}$$
이 세가지 식을 풀면 a를 구할 수 있다.
그런데 이때 아래의 두 식에서 $\alpha, \beta$에 대해 quadrand ambiguity가 존재한다. 이에 따라 $\alpha$와 $\beta$가 두개 씩 존재한다. 해가 4개가 존재하는 것이다.
앞 글에서 Lambert's problem의 solution이 4개가 존재한다는 것을 기하학적으로 살펴 보았다! 즉, 기하학적인 조건에 따라 $\alpha, \beta$를 설정해 주어야 한다.
우선 $\beta$의 부호를 정해보자.
$$s(s-c) = \frac{(r_1+r_2+c)(r_1+r_2-c)}{4} = \frac{(r_1+r_2)^2-c^2}{4} $$
$$ = \frac{2r_1r_2+2r_1r_2cos\theta}{4} = r_1r_2\frac{1+cos\theta}{2}=r_1r_2cos^2\frac{\theta}{2}$$
또한
$$ s(s-c) = 4a^2sin^2\frac{\alpha}{2}sin^2\frac{\beta}{2} $$
그러므로
$$\sqrt{r_1r_2}cos\frac{\theta}{2}=\pm\sqrt{s(s-c)}=\pm2asin\frac{\alpha}{2}sin\frac{\beta}{2}$$
이때 $a, sin\frac{\alpha}{2}$는 양수이므로 $\beta$는 $\theta$와 어떤 관계가 있다는 것을 알 수 있다. 여기서 수식 (2)를 보자. $\beta$가 음수이면 같은 $a, \alpha$에 대해 TOF가 크다. 이를 통상적으로, $\theta \in [\pi, 2\pi]$일때 $\beta \in [-\pi,0]$로 정한다. 같은 의미로, $\theta \in [0, \pi]$ 일때 $\beta \in [0, \pi]$로 한다.
즉 $\beta$는 예각 전파, 둔각 전파를 결정하는 요소이다.
그렇다면 $\alpha$는 무엇을 의미할까. 전파 시간을 의미한다. 두 궤도 중 TOF가 작은 궤도와 큰 궤도를 나누는 것이 된다.
아래 그림을 통해 확인해 보자.
Solution Curve
이 그래프는 두 위치 $[8000,0,0], [7500,1000,500]$에 대해 푼 Lambert problem의 솔루션 이다. $x$축은 궤도 장반경을 최소 장반경으로 나눈 값(무차원값) 이고 y축은 전파 시간(s)이다.
그래프에서 알 수 있듯이, 전이 궤도의 최소 장반경이 존재한다. 해당 장반경보다 작은 장반경으로 전파할 수 없다는 의미이다. 이때의 시간을 $T_{amin}$이라고 하자. $T_{amin}$에 대한 자세한 내용은 아래를 참고하길 바란다.
이제 $\alpha$에 대한 이야기로 돌아오자. $\alpha \in [0,\pi]$인 경우 전파 시간은 $T_{amin}$보다 작다. 반대로 $\alpha \in [\pi, 2\pi]$인 경우 전파 시간은 $T_{amin}$보다 크다.
총 4개의 전이 궤도 중, $\beta$의 범위에 따라 둔각 혹은 예각 전파 궤도를 선택할 수 있다. 그리고 남은 두 궤도 중 한 궤도는 $T_{amin}$보다 크고 나머지 궤도는 $T_{amin}$보다 작다는 의미이다. 이제 설계자가 원하는 궤도에 따라 $\alpha, \beta$의 범위를 정하면 된다. 그렇게 선택된 하나의 궤도는 Lambert's problem의 해가 된다.
사실 원하는 것은 속도이기 때문에, 해를 구했다고 하기에는 아직 이르다. 지금까지 전이 궤도의 6요소를 구했다면 이제 남은 것은 이 요소들을 통해 속도를 계산하는 것이다.
$T_{amin}$
$T_{amin}$을 구해보자.
$$a_{min} = \frac{s}{2}$$
$$\downarrow$$
$$sin^2(\frac{\alpha}{2})=\frac{s}{2a_{min}}=1, \alpha = \pi$$
$$\frac{s-c}{2a} = sin^2(\frac{\beta}{2})$$
$$\downarrow$$
$$\therefore T_{amin} = \sqrt{\frac{a^3}{\mu}}(\pi-(\beta-sin\beta)))$$
Summary
$$\textit{Input : }\vec{r}_1, \vec{r}_2, TOF(=t_2-t_1), M$$
$$\downarrow$$
- Step 1
$\theta = cos^{-1}\frac{\vec{r}_1 \cdot \vec{r}_2}{ |\vec{r}_1| |\vec{r}_2|}$ (예각, 둔각 구분 필요)
$ c = | \vec{r}_1 - \vec{r_2} |$
$ s = \frac{|\vec{r}_1|+|\vec{r}_2|+c}{2}$
- Step 2
$\sqrt\mu(t_2-t_1) = a^{3/2}((\alpha-sin\alpha) - (\beta - sin\beta) + 2M\pi)$
$\frac{s}{2a} = sin^2\frac{\alpha}{2}, \textit{ (if TOF } > T_{amin}, \alpha = 2pi - \alpha \textit{)} $
$\frac{s-c}{2a} = sin^2\frac{\beta}{2}, \textit{ (if } \theta > \pi, \beta = - \beta \textit{)} $
- Step 3 (Find a)
Step 2의 함수들은 입력은 a이고 출력이 TOF인 Object function이다. Root finding method를 통해 a를 찾을 수 있다.
2023.10.22 - [수치해석] - [수치해석] 1. 근사해 구하기. bisect, newton raphson
여기서 주의할 점이 있다. 입력된 a의 범위를 {a_min} 이상으로 설정해 주어야 한다. 그렇지 않으면 수학적 오류가 발생한다.
여기까지가 이 글에서 다룬 부분이며 Lambert's problem의 핵심이다.
이제 남은 것은 하나의 과정뿐이다.
- Step 4
Find $e, i, \omega, \Omega, \nu$
Find $\vec{v_1}, \vec{v_2}$
Step 4는 지금까지 다룬 내용과는 또 다른 내용이다. 지금까지는 전이 시간에 따른 궤도 장반경을 구하는 것이 목적이었다면, 이제는 궤도의 장반경과 두 벡터를 통해 궤도의 6요소를 구하는 것이 목적이다. 해당 내용은 다음 글에서 자세하게 다룰 것이다.
