케플러 문제 (Kepler’s problem) - 4

우주비행체의 비행시간과 실제 비행각과의 함수 관계를 다루는 문제를 케플러 문제 (Kepler's problem)라고 한다. 케플러 문제는 비행시간(time of flight) 계산 문제와 예측(prediction) 문제로 나눌 수 있다.

비행시간 계산 문제는 시간 $t=t_0$ 에서 실제 비행각(true anomaly) ${\theta }_0$가 주어졌을 때 비행각이 $\mathrm{\Delta }\theta$ 만큼 변화하기까지 필요한 비행시간 $t-t_0$ 을 계산하는 문제다. 예측 문제는 비행시간 계산 문제의 역으로서 시간 $t=t_0$ 에서 실제 비행각 ${\theta }_0$ 과 비행시간 $t-t_0$ 이 주어졌을 때 실제 비행각 $\theta (t)$ 를 계산하는 문제다.

 

 

이전 게시글을 통해서 케플러 문제를 푸는데 필요한 케플러 방정식을 이미 유도한 바 있다. 여기서는 이 방정식을 종합하여 케플러 문제의 해를 구하는 알고리즘에 대해서 설명하고자 한다.

타원궤도에서 케플러 방정식을 다음과 같이 유도하였었다 (https://pasus.tistory.com/308).

 

$$\begin{array}{r}\text{(1)}& t-t_0& =\frac{h^3}{{\mu }^2}\frac{1}{(1-e^2{)}^{3/2}}(M_e-M_{e0})\\ & =\sqrt{\frac{a^3}{\mu }}(M_e-M_{e0})\end{array}$$

 

여기서 타원궤도의 평균 비행각(mean anomaly) $M_e$, 이심 비행각(eccentric anomaly) $E$, 실제 비행각 $\theta$ 의 관계식은 다음과 같다.

 

$$\begin{array}{r}\text{(2)}& & M_e=E-e\sin E\\ \text{(3)}& & \tan \frac{E}{2}=\sqrt{\frac{1-e}{1+e}}\tan \frac{\theta }{2}\end{array}$$

 

 

 

타원궤도를 돌고 있는 우주비행체의 위치벡터와 속도벡터 또는 고전 궤도요소가 시간 $t=t_0$ 에서 ${\overrightarrow{r}}_0,{\overrightarrow{v}}_0$ 또는 $(a,e,i,\mathrm{\Omega },\omega ,{\theta }_0)$ 로 주어졌을 때, 실제 비행각이 $\theta$ 에 도달하기 까지 걸리는 비행시간 $t-t_0$ 은 다음과 같은 순서로 계산할 수 있다.

[1] 먼저 식 (3)으로 주어진 ${\theta }_0$ 와 $\theta$ 에 대응하는 $E_0$ 와 $E$ 를 계산한다.
[2] 그리고 식 (2)로 $E_0$ 와 $E$ 에 대응하는 $M_{e0}$ 와 $M_e$ 를 계산한다.
[3] 마지막으로 식 (1)로 비행시간 $t-t_0$ 를 계산하면 된다.

타원궤도를 돌고 있는 우주비행체의 위치벡터와 속도벡터 또는 고전 궤도요소가 시간 $t=t_0$ 에서 ${\overrightarrow{r}}_0,{\overrightarrow{v}}_0$ 또는 $(a,e,i,\mathrm{\Omega },\omega ,{\theta }_0)$ 로 주어졌을 때, 비행시간이 $t-t_0$ 만큼 경과한 후의 실제 비행각 $\theta$ 는 다음과 같은 순서로 계산할 수 있다.

[1] 먼저 식 (3)과 (2)로 주어진 ${\theta }_0$ 에 대응하는 $E_0$ 와 $M_{e0}$ 를 계산한다.
[2] 식 (1)을 이용하여 $M_{e0}$ 와 비행시간 $t-t_0$ 에 대응하는 $M_e$ 를 계산한다.
[3] 식 (2)로 $M_e$ 에 대응하는 $E$ 를 계산한다.
[4] $E$ 를 이용하여 식 (3)으로 $\theta$ 를 계산한다.

케플러 방정식 (2)는 비선형 방정식이므로 해석적으로 해를 구하기가 어렵기 때문에 수치적인 반복계산법으로 계산해야 한다. 뉴톤-랩슨 방법(Newton-Raphson method)를 이용하면, 식 (2)를 다음과 같은 함수로 놓고,

 

$$\begin{array}{}\text{(4)}& f(E)=E-e\sin E-M_e\end{array}$$

 

$|E_{i+1}-E_i|<ϵ$ 의 조건이 수렴할 때까지 반복하여 업데이트 한다.

 

$$\begin{array}{}\text{(5)}& E_{i+1}=E_i-\frac{f(E_i)}{f^{\prime }(E_i)}\end{array}$$

 

여기서

 

$$\begin{array}{}\text{(6)}& f^{\prime }(E)=1-e\cos E\end{array}$$

 

이고 초기값 $E_0$ 는 다음 $M_e$ 와 $E$ 의 그래프를 참고하여 적절하게 선택한다.

 

 

 

 

포물선궤도에서의 케플러 방정식은 다음과 같았다 (https://pasus.tistory.com/307).

 

$$\begin{array}{}\text{(7)}& t-t_0=\frac{h^3}{{\mu }^2}(M_p-M_{p0})\end{array}$$

 

여기서 포물선궤도의 평균 비행각 $M_p$ 와 실제 비행각 $\theta$ 의 관계식은 다음과 같다.

 

$$\begin{array}{}\text{(8)}& M_p=\frac{1}{2}\tan \frac{\theta }{2}+\frac{1}{6}{\tan }^3\frac{\theta }{2}\end{array}$$

 

식 (8)에서 $\tan \frac{\theta }{2}$ 의 실수해는 다음과 같다.

 

$$\begin{array}{}\text{(9)}& \tan \frac{\theta }{2}={(3M_p+\sqrt{(3M_p{)}^2+1})}^{\frac{1}{3}}-{(3M_p+\sqrt{(3M_p{)}^2+1})}^{-\frac{1}{3}}\end{array}$$

 

 

 

포물선궤도에 있는 우주비행체의 위치벡터와 속도벡터 또는 고전 궤도요소가 시간 $t=t_0$ 에서 ${\overrightarrow{r}}_0,{\overrightarrow{v}}_0$ 또는 $(a,e,i,\mathrm{\Omega },\omega ,{\theta }_0)$ 로 주어졌을 때, 실제 비행각이 $\theta$ 에 도달하기 까지 걸리는 비행시간 $t-t_0$ 은 다음과 같은 순서로 계산할 수 있다.

[1] 식 (8)로 주어진 ${\theta }_0$ 와 $\theta$ 에 대응하는 $M_{p0}$ 와 $M_p$ 를 계산한다.
[2] 그리고 식 (7)로 비행시간 $t-t_0$ 를 계산하면 된다.

포물궤도에 있는 우주비행체의 위치벡터와 속도벡터 또는 고전 궤도요소가 시간 $t=t_0$ 에서 ${\overrightarrow{r}}_0,{\overrightarrow{v}}_0$ 또는 $(a,e,i,\mathrm{\Omega },\omega ,{\theta }_0)$ 로 주어졌을 때, 비행시간이 $t-t_0$ 만큼 경과한 후의 실제 비행각 $\theta$ 는 다음과 같은 순서로 계산할 수 있다.

[1] 먼저 식 (8) 로 주어진 ${\theta }_0$ 에 대응하는 $M_{p0}$ 를 계산한다.
[2] 식 (7)을 이용하여 $M_{p0}$ 와 비행시간 $t-t_0$ 에 대응하는 $M_p$ 를 계산한다.
[3] 식 (9)로 $M_p$ 에 대응하는 $\theta$ 를 계산한다.

쌍곡선궤도에서는 케플러 방정식이 다음과 같았다 (https://pasus.tistory.com/309).

 

$$\begin{array}{r}\text{(10)}& t-t_0& =\frac{h^3}{{\mu }^2}\frac{1}{(e^2-1{)}^{3/2}}(M_h-M_{h0})\\ & =\sqrt{\frac{-a^3}{\mu }}(M_h-M_{h0})\end{array}$$

 

여기서 쌍곡선궤도의 평균 비행각 $M_h$, 이심 비행각 $F$, 실제 비행각 $\theta$ 의 관계식은 다음과 같다.

 

$$\begin{array}{r}\text{(11)}& & M_h=e\sinh F-F\\ \text{(12)}& & \tanh \frac{F}{2}=\sqrt{\frac{e-1}{e+1}}\tan \frac{\theta }{2}\end{array}$$

 

 

 

쌍곡선궤도에 있는 우주비행체의 위치벡터와 속도벡터 또는 고전 궤도요소가 시간 $t=t_0$ 에서 ${\overrightarrow{r}}_0,{\overrightarrow{v}}_0$ 또는 $(a,e,i,\mathrm{\Omega },\omega ,{\theta }_0)$ 로 주어졌을 때, 실제 비행각이 $\theta$ 에 도달하기 까지 걸리는 비행시간 $t-t_0$ 는 다음과 같은 순서로 계산할 수 있다.

[1] 먼저 식 (12)로 주어진 ${\theta }_0$ 와 $\theta$ 에 대응하는 $F_0$ 와 $F$ 를 계산한다.
[2] 그리고 식 (11)로 $F_0$ 와 $F$ 에 대응하는 $M_{h0}$ 와 $M_h$ 를 계산한다.

[3] 마지막으로 식 (10)으로 비행시간 $t-t_0$ 를 계산하면 된다.

쌍곡선궤도에 있는 우주비행체의 위치벡터와 속도벡터 또는 고전 궤도요소가 시간 $t=t_0$ 에서 ${\overrightarrow{r}}_0,{\overrightarrow{v}}_0$ 또는 $(a,e,i,\mathrm{\Omega },\omega ,{\theta }_0)$ 로 주어졌을 때, 비행시간이 $t-t_0$ 만큼 경과한 후의 실제 비행각 $\theta$ 는 다음과 같은 순서로 계산할 수 있다.

[1] 먼저 식 (12)와 (11)로 주어진 ${\theta }_0$ 에 대응하는 $F_0$ 와 $M_{h0}$ 를 계산한다.
[2] 식 (10)을 이용하여 $M_{h0}$ 와 비행시간 $t-t_0$ 에 대응하는 $M_h$ 를 계산한다.
[3] 식 (11)로 $M_h$ 에 대응하는 $F$ 를 계산한다.
[4] $F$ 를 이용하여 식 (12)로 $\theta$ 를 계산한다.

쌍곡선궤도의 케플러 방정식 (11)은 비선형 방정식이므로 해석적으로 해를 구하기가 어렵기 때문에 수치적인 반복계산법으로 계산해야 한다. 마찬가지로 뉴톤-랩슨 방법을 이용하면, 식 (11)을 다음과 같은 함수로 놓고,

 

$$\begin{array}{}\text{(13)}& f(F)=e\sinh F-F-M_h\end{array}$$

 

$|F_{i+1}-F_i|<ϵ$ 의 조건이 수렴할 때까지 반복하여 업데이트 한다.

 

$$\begin{array}{}\text{(14)}& F_{i+1}=F_i-\frac{f(F_i)}{f^{\prime }(F_i)}\end{array}$$

 

여기서

 

$$\begin{array}{}\text{(15)}& f^{\prime }(F)=e\cosh F-1\end{array}$$

 

이고 초기값 $F_0$ 는 다음 $M_h$ 와 $F$ 의 그래프를 참고하여 적절하게 선택한다.