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

타원궤도, 포물선궤도, 쌍곡선궤도의 케플러 방정식을 이용하여 케플러 문제 (Kepler's problem)를 풀어보았는데(https://pasus.tistory.com/313). 이번에는 범용(universal) 케플러 방정식을 이용하여 케플러 문제를 풀어보도록 하겠다.

범용 케플러 방정식은 다음과 같았다. (https://pasus.tistory.com/310).

 

\[ \begin{align} \sqrt{\mu} (t-t_0 )=S(z) \chi^3+ \frac{ \vec{r}_0 \cdot \vec{v}_0}{\sqrt{\mu}} \chi^2 C(z)+r_0 \chi \left( 1-zS(z) \right) \tag{1} \end{align} \]

 

여기서 \(z= \frac{\chi^2}{a}\) 이다.

범용변수(universal variable) \(\chi\) 와 타원궤도의 이심 비행각 \(E\), 포물선궤도의 실제 비행각 \(\theta\), 그리고 쌍곡선궤도의 이심 비행각 \(F\) 와의 관계식은 다음과 같다.

 

\[ \begin{align} \chi &= \sqrt{a} (E-E_0) \tag{2} \\ \\ \chi &= \frac{h}{\sqrt{\mu}} \left( \tan \frac{\theta}{2}- \tan \frac{\theta_0}{2} \right) \\ \\ \chi &= \sqrt{-a} (F-F_0) \end{align} \]

 

식 (2)에서 보듯이 범용변수 \(\chi\) 는 각 궤도의 이심 비행각과 연결되어 있기 때문에 범용 '비행각' (universal anomaly)라고도 한다. 참고로 \(\chi\) 의 단위는 \(\sqrt{km}\) 이다.

 

 

우주비행체의 위치벡터와 속도벡터가 시간 \(t=t_0\) 에서 \(\vec{r}_0, \vec{v}_0\) 로 주어졌을 때, 비행시간이 \(t-t_0\) 만큼 경과한 후의 범용 비행각 \(\chi\) 는 다음과 같은 순서로 계산할 수 있다.

[1] \(\vec{r}_0\) 와 \(\vec{v}_0\) 로부터 \(r_0\) 와 장반경 \(a\) 를 계산한다. 여기서 포물선궤도의 경우 \(a=\infty\) 이기 때문에 장반경 \(a\) 의 역수인 \(\alpha=\frac{1}{a}\) 를 사용하는 것이 편리하다.
[2] 식 (1)을 이용하여 비행시간 \(t-t_0\) 에 대응하는 \(\chi\) 를 계산한다.

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

 

\[ \begin{align} f(\chi) &=S(z) \chi^3+\frac{\vec{r}_0 \cdot \vec{v}_0}{\sqrt{\mu}} \chi^2 C(z)\tag{3} \\ \\ & \ \ \ \ \ \ \ +r_0 \chi (1-zS(z)) - \sqrt{\mu} (t-t_0 ) \end{align} \]

 

\(z= \frac{\chi^2}{a}= \alpha \chi^2\) 의 관계식을 이용하면 위 식을 다음과 같이 쓸 수 있다.

 

\[ \begin{align} f(\chi) &= S(z) \chi^3+ \frac{\vec{r}_0 \cdot \vec{v}_0}{ \sqrt{\mu}} \chi^2 C(z) \tag{4} \\ & \ \ \ \ \ +r_0 \chi (1-\alpha \chi^2 S(z))-\sqrt{\mu} (t-t_0 ) \\ \\ &= \frac{\vec{r}_0 \cdot \vec{v}_0}{ \sqrt{\mu} } \chi^2 C(z)+(1-\alpha r_0 ) \chi^3 S(z)+r_0 \chi- \sqrt{\mu} (t-t_0 ) \end{align} \]

 

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

 

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

 

\(f^\prime (\chi_i )\) 을 계산해야 한다. \(f^\prime (\chi_i )\) 는 다음과 같이 계산할 수 있다.

 

\[ \begin{align} \frac{f(\chi)}{d\chi} &= 2 \frac{\vec{r}_0 \cdot \vec{v}_0}{\sqrt{\mu}} \chi C(z)+\frac{\vec{r}_0 \cdot \vec{v}_0}{\sqrt{\mu}} \chi^2 \frac{dC(z)}{dz} \frac{dz}{d\chi} \tag{6} \\ \\ & \ \ +3(1-\alpha r_0 ) \chi^2 S(z)+(1-\alpha r_0 ) \chi^3 \frac{dS(z)}{dz} \frac{dz}{d\chi}+r_0 \end{align} \]

 

여기서

 

\[ \begin{align} \frac{dz}{d \chi} &= 2 \alpha \chi \tag{7} \\ \\ \frac{dC(z)}{dz} &= - \frac{1}{z^2} -\frac{ -z \left( \frac{1}{2 \sqrt{z}} \right) \sin \sqrt{z}- \cos \sqrt{z} }{z^2} \tag{8} \\ \\ & = \frac{1}{z^2} \left( -1+ \frac{\sqrt{z}}{2} \sin \sqrt{z} + \cos \sqrt{z} \right) \\ \\ &= \frac{1}{z^2} \left( -1+ \frac{\sqrt{z}}{2} \left( \sqrt{z}-(\sqrt{z})^3 S(z) \right)+(1-zC(z)) \right) \\ \\ & = \frac{1}{2z} (1-zS(z)-2C(z)) \\ \\ \frac{dS(z)}{dz} &= -\frac{1}{z^2} - \frac{ (\sqrt{z})^3 \cos \sqrt{z} \left( \frac{1}{2 \sqrt{z}} \right)-\left( \frac{3}{2} \sqrt{z} \right) \sin \sqrt{z} }{z^3} \tag{9} \\ \\ &= \frac{1}{z^3} \left( -z-\frac{z}{2} \cos \sqrt{z}+ \frac{3}{2} \sqrt{z} \sin \sqrt{z} \right) \\ \\ &= \frac{1}{z^3} \left( -z- \frac{z}{2} (1-zC(z)) + \frac{3}{2} \sqrt{z} \left( \sqrt{z}-(\sqrt{z})^3 S(z) \right) \right) \\ \\ &= \frac{1}{2z} (C(z)-3S(z)) \end{align} \]

 

이므로 식 (6)은 \(z= \alpha \chi^2\) 임을 이용하면 다음과 같이 쓸 수 있다.

 

\[ \begin{align} \frac{f(\chi)}{d\chi} &= 2 \frac{\vec{r}_0 \cdot \vec{v}_0}{ \sqrt{\mu}} \chi C(z) \tag{10} \\ \\ & \ \ \ \ \ +\frac{\vec{r}_0 \cdot \vec{v}_0}{ \sqrt{\mu}} \chi^2 \left( \frac{1}{2 \alpha \chi^2 } (1- \alpha \chi^2 S(z)-2C(z)) \right) 2 \alpha \chi \\ \\ & \ \ \ \ \ +3(1-\alpha r_0 ) \chi^2 S(z) \\ \\ & \ \ \ \ \ +(1-\alpha r_0 ) \chi^3 \left( \frac{1}{2 \alpha \chi^2 } (C(z)-3S(z)) \right) 2 \alpha \chi +r_0 \\ \\ &= \frac{\vec{r}_0 \cdot \vec{v}_0}{ \sqrt{\mu}} \chi (1-\alpha \chi^2 S(z))+(1- \alpha r_0 ) \chi^2 C(z)+r_0 \end{align} \]

 

식 (5)에서 \(\chi\) 의 초기값을 \(\chi_0= \sqrt{\mu} \vert \alpha \vert (t-t_0)\) 로 하는 것이 일반적이다.

 

 

만약 범용 비행각으로부터 실제 비행각 \(\theta\) 를 계산하고자 한다면, \(\alpha\) 또는 \(a\) 로 궤도의 모양을 판단하여 식 (2)로 궤도의 이심 비행각의 변화량을 계산해야 한다. 그런 후 다음 관계식을 이용하여 \(\theta\) 를 계산하면 된다.

 

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

 

\(\alpha\) 는 다음과 같이 에너지 방정식으로부터 유도할 수 있다.

 

\[ \begin{align} \alpha = \frac{1}{a}= \frac{2}{r_0} -\frac{v_0^2}{\mu} \tag{12} \end{align} \]

 

\(\alpha \gt 0\) 이면 타원궤도, \(\alpha=0\) 이면 포물선궤도, \(\alpha \lt 0\) 이면 쌍곡선궤도이다.

초기 비행각 \(\theta_0\) 는 \(\vec{r}_0\) 와 \(\vec{v}_0\) 를 이용하여 계산할 수 있다.

다음은 \(\chi\) 의 계산 과정을 매트랩 코드로 구현한 것이다.

 

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

 

universal_chi.m

 

function chi = universal_chi(r0_vec, v0_vec, tof)

% chi = universal_chi(r0_vec, v0_vec, tof)
%        given r0_vec, r0_vec and time-of-flight
%        compute universal alomaly chi
%       
% Input  r0_vec: initial position vector in km
%        v0_vec: initial velocity vector in km/sec
%        tof: time of flight in sec
% Output chi: universal anomaly in sqrt(km)
%
% all in inertial frame
%
% coded by st.watermelon

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


% step 1: compute r0 and a ---------------
r0 = sqrt(r0_vec'*r0_vec);
v0 = sqrt(v0_vec'*v0_vec);

alpha = 2/r0 - v0^2/mu_e; % alpha = 1/a 

% step 2 ----------------------------------

eps = 1e-8;

% compute chi
chi = sqrt(mu_e)*abs(alpha)*tof;
for ii = 1:100

    % Stumpff
    z = alpha * chi^2;
    [Cz, Sz] = stumpff(z);

    % Newton-Raphson
    f_chi = (r0_vec'*v0_vec)/sqrt(mu_e)*chi^2*Cz ...
             + (1-alpha*r0)*chi^3*Sz + r0*chi - sqrt(mu_e)*tof;
    dfdchi = (r0_vec'*v0_vec)/sqrt(mu_e)*chi*(1-alpha*chi^2*Sz) ...
             + (1-alpha*r0)*chi^2*Cz + r0;

    chi_next = chi - f_chi/dfdchi;

    if abs(chi_next-chi) < eps
        break;
    else
        chi = chi_next;
    end
   
end

 

 

 

만약 우주비행체의 초기 위치벡터 \(\vec{r}_0\) 와 속도벡터 \(\vec{v}_0\) 가 주어졌을 때, 비행시간이 \(t-t_0\) 만큼 변화한 후의 위치벡터와 속도벡터 \(\vec{r}, \vec{v}\) 를 구하고자 한다면 범용변수와 라그랑지 계수(Lagrange coefficients)와의 관계식을 이용하면 된다.

라그랑지 계수는 다음과 같았다 (https://pasus.tistory.com/312).

 

\[ \begin{align} f &= 1- \frac{\chi^2}{r_0} C(z) \tag{13} \\ \\ g &= (t-t_0 )- \frac{\chi^3}{\sqrt{\mu}} S(z) \\ \\ \dot{f} &= \frac{ \sqrt{\mu}}{rr_0 } \chi (zS(z)-1) \tag{14} \\ \\ \dot{g} &=1- \frac{\chi^2}{r} C(z) \end{align} \]

 

위치벡터 및 속도벡터와 라그랑지 계수와의 관계식은 다음과 같았다 (https://pasus.tistory.com/311).

 

\[ \begin{align} \vec{r} &= f \vec{r}_0+g \vec{v}_0 \tag{15} \\ \\ \vec{v} &= \dot{f} \vec{r}_0+ \dot{g} \vec{v}_0 \end{align} \]

 

계산 순서는 다음과 같다.

[1] \(\vec{r}_0\) 와 \(\vec{v}_0\) 로부터 \(r_0\) 와 장반경 \(a\) 를 계산한다.
[2] \(t-t_0\) 로부터 범용변수 \(\chi\) 를 계산한다.
[3] 식 (13) 으로 라그랑지 계수 \(f, g\) 를 계산한다.
[4] 식 (15) 로 \(\vec{r}\) 를 계산한다. 그리고 \(\vec{r}\) 의 크기인 \(r\) 도 계산한다.
[5] 식 (14) 로 라그랑지 계수 \(\dot{f}, \dot{g}\) 를 계산한다.
[6] 식 (15) 로 \(\vec{v}\) 를 계산한다.

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

 

kepler_chi.m

 

function [r_vec, v_vec] = kepler_chi(r0_vec, v0_vec, tof)

% [r_vec, v_vec] = kepler_chi(r0_vec, v0_vec, tof)
%        given r0_vec, r0_vec and time-of-flight
%        compute r_vec and v_vec
%       
% Input  r0_vec: initial position vector in km
%        v0_vec: initial velocity vector in km/sec
%        tof: time of flight in sec
% Output r_vec: final position vector in km
%        v_vec: final velocity vector in km/sec
%
% all in inertial frame
%
% coded by st.watermelon

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

% step 1 : compute r0 and v0
r0 = sqrt(r0_vec'*r0_vec);
v0 = sqrt(v0_vec'*v0_vec);

alpha = 2/r0 - v0^2/mu_e; % alpha = 1/a 

% step 2 : compute chi
chi = universal_chi(r0_vec, v0_vec, tof);

% compute C(z), S(z) and z
z = alpha * chi^2;
[Cz, Sz] = stumpff(z);

% step 3 : Lagrange coefficients f and g
f = 1 - chi^2/r0 * Cz;
g = tof - chi^3/sqrt(mu_e)*Sz;

% step 4 : compute r_vec
r_vec = f*r0_vec + g*v0_vec;
r = sqrt(r_vec'*r_vec);

% step 5 : Lagrange coefficients f_dot and g_dot
f_dot = sqrt(mu_e)/(r*r0) * chi * (z*Sz-1);
g_dot = 1 - chi^2/r*Cz;

v_vec = f_dot*r0_vec + g_dot*v0_vec;