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

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

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

 

 

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

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

 

\[ \begin{align} t-t_0 &= \frac{h^3}{\mu^2} \frac{1}{(1-e^2 )^{3/2}} (M_e-M_{e0} ) \tag{1} \\ \\ &= \sqrt{ \frac{a^3}{\mu}} (M_e-M_{e0} ) \end{align} \]

 

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

 

\[ \begin{align} & M_e=E-e \sin E \tag{2} \\ \\ & \tan \frac{E}{2} = \sqrt{ \frac{1-e}{1+e}} \tan \frac{\theta}{2} \tag{3} \end{align} \]

 

 

 

타원궤도를 돌고 있는 우주비행체의 위치벡터와 속도벡터 또는 고전 궤도요소가 시간 \(t=t_0\) 에서 \(\vec{r}_0, \vec{v}_0\) 또는 \( (a, e, i, \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\) 에서 \(\vec{r}_0, \vec{v}_0\) 또는 \( (a, e, i, \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{align} f(E)=E-e \sin E-M_e \tag{4} \end{align} \]

 

\(\vert E_{i+1}-E_i \vert \lt \epsilon\) 의 조건이 수렴할 때까지 반복하여 업데이트 한다.

 

\[ \begin{align} E_{i+1}=E_i- \frac{f(E_i )}{f^\prime (E_i ) } \tag{5} \end{align} \]

 

여기서

 

\[ \begin{align} f^\prime (E)=1-e \cos E \tag{6} \end{align} \]

 

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

 

 

 

 

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

 

\[ \begin{align} t-t_0 = \frac{h^3}{\mu^2} (M_p-M_{p0} ) \tag{7} \end{align} \]

 

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

 

\[ \begin{align} M_p= \frac{1}{2} \tan \frac{\theta}{2} + \frac{1}{6} \tan^3 \frac{\theta}{2} \tag{8} \end{align} \]

 

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

 

\[ \begin{align} \tan \frac{\theta}{2}= \left( 3M_p+ \sqrt{ (3M_p )^2+1} \right)^{\frac{1}{3}}- \left( 3M_p+ \sqrt{ (3M_p )^2+1} \right)^{-\frac{1}{3}} \tag{9} \end{align} \]

 

 

 

포물선궤도에 있는 우주비행체의 위치벡터와 속도벡터 또는 고전 궤도요소가 시간 \(t=t_0\) 에서 \(\vec{r}_0, \vec{v}_0\) 또는 \( (a, e, i, \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\) 에서 \(\vec{r}_0, \vec{v}_0\) 또는 \( (a, e, i, \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{align} t-t_0 &= \frac{h^3}{\mu^2} \frac{1}{(e^2-1)^{3/2} } (M_h-M_{h0})\tag{10} \\ \\ &= \sqrt{ \frac{-a^3}{\mu}} (M_h-M_{h0} ) \end{align} \]

 

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

 

\[ \begin{align} & M_h=e \sinh F-F \tag{11} \\ \\ & \tanh⁡ \frac{F}{2}= \sqrt{ \frac{e-1}{e+1}} \tan \frac{\theta}{2} \tag{12} \end{align} \]

 

 

 

쌍곡선궤도에 있는 우주비행체의 위치벡터와 속도벡터 또는 고전 궤도요소가 시간 \(t=t_0\) 에서 \(\vec{r}_0, \vec{v}_0\) 또는 \( (a, e, i, \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\) 에서 \(\vec{r}_0, \vec{v}_0\) 또는 \((a, e, i, \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{align} f(F)=e \sinh F-F-M_h \tag{13} \end{align} \]

 

\(\vert F_{i+1}-F_i \vert \lt \epsilon\) 의 조건이 수렴할 때까지 반복하여 업데이트 한다.

 

\[ \begin{align} F_{i+1}=F_i- \frac{f(F_i )}{f^\prime (F_i ) } \tag{14} \end{align} \]

 

여기서

 

\[ \begin{align} f^\prime (F)=e \cosh F-1 \tag{15} \end{align} \]

 

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