라그랑지 계수 (Lagrange coefficients) - 1

어느 우주비행체의 초기 위치벡터 \(\vec{r}_0\) 와 속도벡터 \(\vec{v}_0\) 가 주어졌을 때, 실제 비행각(true anomaly)이 \(\Delta \theta\) 만큼 변화한 후, 변화된 위치벡터와 속도벡터 \(\vec{r}, \vec{v}\) 를 초기 위치벡터 및 속도벡터, 그리고 \(\Delta \theta\) 의 함수로 표현하고자 한다.

 

 

우주비행체는 궤도면(orbital plane) 상에서만 운동하므로(https://pasus.tistory.com/96) 위치벡터와 속도벡터 \(\vec{r}, \vec{v}\) 는 항상 궤도면 상에 존재한다. 따라서 임의의 시간에서의 위치벡터와 속도벡터는 초기 위치벡터와 속도벡터 \(\vec{r}_0, \vec{v}_0\) 의 선형 조합으로 표현할 수 있다. 즉,

 

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

 

여기서 \(f\) 와 \(g\) 는 시간의 함수로서 라그랑지 계수(Lagrange coefficients) 또는 '\(f\)와 \(g\) 함수' 라고 한다. 따라서 \(f\) 와 \(g\) 를 \(\Delta \theta\) 의 함수로 나타낼 수 있다면 목적을 달성할 수 있다.

먼저 궤도중심좌표계(perifocal frame)에서 위치벡터와 속도벡터를 구해 보자 (https://pasus.tistory.com/288). 참고로 궤도중심좌표계는 관성좌표계이므로 3차원 공간상에 고정되어 있다.

 

\[ \begin{align} & \vec{r}=x \ \hat{p}_1+y \ \hat{p}_2 \tag{2} \\ \\ & \vec{v}= \dot{x} \ \hat{p}_1+ \dot{y} \ \hat{p}_2 \end{align} \]

 

여기서

 

\[ \begin{align} & x=r \cos \theta, \ \ \ \ \ y=r \sin \theta \tag{3} \\ \\ & \dot{x}= -\frac{\mu}{h} \sin \theta, \ \ \ \ \ \dot{y}= \frac{\mu}{h} (e+ \cos \theta ) \end{align} \]

 

이다. 식 (2)에 의해서 초기 위치벡터와 속도벡터도 다음과 같이 표현할 수 있다.

 

\[ \begin{align} & \vec{r}_0=x_0 \ \hat{p}_1+y_0 \ \hat{p}_2 \tag{4} \\ \\ & \vec{v}_0= \dot{x}_0 \ \hat{p}_1+ \dot{y}_0 \ \hat{p}_2 \end{align} \]

 

식 (4)로부터 궤도중심좌표계의 \(\hat{p}_1\) 축과 \(\hat{p}_2\) 축을 \(\vec{r}_0\) 와 \(\vec{v}_0\) 의 함수로 나타낼 수 있다면 식 (2)에 의해서 식 (1)과 같은 관계식을 얻을 수 있을 것이다.

식 (4)에서 \(\hat{p}_2\) 를 계산하면,

 

\[ \hat{p}_2= \frac{1}{y_0} \vec{r}_0- \frac{x_0}{y_0} \hat{p}_1 \tag{5} \]

 

인데 이를 \(\vec{v}_0\) 식에 대입하면 다음과 같이 된다.

 

\[ \begin{align} \vec{v}_0 &= \dot{x}_0 \ \hat{p}_1+ \dot{y}_0 \left( \frac{1}{y_0} \vec{r}_0- \frac{x_0}{y_0} \hat{p}_1 \right) \tag{6} \\ \\ &= \left( \frac{\dot{x}_0 y_0- \dot{y}_0 x_0}{y_0} \right) \hat{p}_1+ \frac{ \dot{y}_0}{y_0} \vec{r}_0 \end{align} \]

 

한편 각운동량 벡터는 정의에 의해서 다음과 같이 계산할 수있다.

 

\[ \begin{align} \vec{h}_0 &= \vec{r}_0 \times \vec{v}_0=(x_0 \dot{y}_0-y_0 \dot{x}_0 ) \hat{p}_3 \tag{7} \\ \\ &= h \ \hat{p}_3 \end{align} \]

 

여기서

 

\[ h=x_0 \dot{y}_0-y_0 \dot{x}_0 \]

 

이다. \(h\) 를 식(6)에 대입하면 \(\vec{v}_0\) 는 다음과 같이 된다.

 

\[ \vec{v}_0=- \frac{h}{y_0} \hat{p}_1+ \frac{ \dot{y}_0}{y_0} \vec{r}_0 \tag{8} \]

 

식 (8)에서 \(\hat{p}_1\) 을 계산하고,

 

\[ \hat{p}_1= \frac{ \dot{y}_0}{h} \vec{r}_0- \frac{y_0}{h} \vec{v}_0 \tag{9} \]

 

이를 식 (5)에 대입하면,

 

\[ \begin{align} \hat{p}_2 &= \frac{1}{y_0} \vec{r}_0- \frac{x_0}{y_0} \left(\frac{ \dot{y}_0}{h} \vec{r}_0- \frac{y_0}{h} \vec{v}_0 \right) \tag{10} \\ \\ &= \frac{h-x_0 \dot{y}_0}{y_0 h} \vec{r}_0+ \frac{x_0}{h} \vec{v}_0 \\ \\ &=- \frac{\dot{x}_0}{h} \vec{r}_0+ \frac{x_0}{h} \vec{v}_0 \end{align} \]

 

이 된다. 이제 \(\hat{p}_1\) 과 \(\hat{p}_2\) 를 \(\vec{r}_0\) 와 \(\vec{v}_0\) 의 함수로 구했으므로 식 (9)와 (10)을 식 (2)에 대입한다.

 

\[ \begin{align} \vec{r} &=x \left( \frac{ \dot{y}_0}{h} \vec{r}_0 - \frac{y_0}{h} \vec{v}_0 \right) + y \left( -\frac{\dot{x}_0}{h} \vec{r}_0+ \frac{x_0}{h} \vec{v}_0 \right) \tag{11} \\ \\ &= \frac{x \dot{y}_0- y \dot{x}_0}{h} \vec{r}_0+ \frac{-xy_0+yx_0}{h} \vec{v}_0 \\ \\ \vec{v} &= \dot{x} \left( \frac{ \dot{y}_0}{h} \vec{r}_0 - \frac{y_0}{h} \vec{v}_0 \right) + \dot{y} \left( -\frac{ \dot{x}_0}{h} \vec{r}_0+ \frac{x_0}{h} \vec{v}_0 \right) \\ \\ &= \frac{\dot{x} \dot{y}_0- \dot{y} \dot{x}_0}{h} \vec{r}_0+ \frac{-\dot{x}y_0+ \dot{y}x_0}{h} \vec{v}_0 \end{align} \]

 

식 (11)과 (1)을 비교하면 라그랑지 계수를 다음과 같이 구할 수 있다.

 

\[ \begin{align} f &= \frac{x \dot{y}_0- y \dot{x}_0}{h}, \ \ \ \ \ g=\frac{-xy_0+yx_0}{h} \tag{12} \\ \\ \dot{f} &= \frac{\dot{x} \dot{y}_0- \dot{y} \dot{x}_0}{h}, \ \ \ \ \ \dot{g}= \frac{-\dot{x}y_0+ \dot{y}x_0}{h} \end{align} \]

 

다시 식 (1)을 이용하며 각운동량 벡터를 계산해 보면,

 

\[ \begin{align} \vec{h} &= \vec{r} \times \vec{v}= (f \vec{r}_0+g \vec{v}_0 ) \times (\dot{f} \vec{r}_0+ \dot{g} \vec{v}_0 ) \tag{13} \\ \\ &=(f \dot{g}-\dot{f}g)( \vec{r}_0 \times \vec{v}_0 ) \\ \\ &=(f \dot{g}- \dot{f}g) \vec{h}_0 \end{align} \]

 

이 되는데, 각운동량 벡터는 상수벡터이므로 \(\vec{h}= \vec{h}_0\) 이기 때문에

 

\[ f \dot{g}- \dot{f}g = 1 \tag{14} \]

 

이 됨을 알 수 있다. 즉 라그랑지 계수 4개는 서로 독립적이지 않고 4개 중 3개만 알면 나머지는 식 (14)를 이용하여 계산할 수 있다.

한편 식 (3)에 의해서 초기 위치와 속도의 성분은 다음과 같이 계산할 수 있으므로,

 

\[ \begin{align} & x_0=r_0 \cos \theta_0, \ \ \ \ \ y_0=r_0 \sin \theta_0 \tag{15} \\ \\ & \dot{x}_0= -\frac{\mu}{h} \sin \theta_0, \ \ \ \ \ \dot{y}_0= \frac{\mu}{h} (e+ \cos \theta_0 ) \end{align} \]

 

이를 식 (3)과 함께 식 (12)에 대입하면 라그랑지 계수를 계산할 수 있다.

 

 

먼저 \(f\) 를 계산해 보자.

 

\[ \begin{align} f &= \frac{x\dot{y}_0-y \dot{x}_0}{h} \tag{16} \\ \\ &= \frac{1}{h} \left[ r \cos \theta \frac{\mu}{h} (e+ \cos \theta_0 )+r \sin \theta \frac{\mu}{h} \sin \theta_0 \right] \\ \\ &= \frac{\mu r}{h^2} \left[ e \cos \theta + \cos \theta \cos \theta_0+ \sin \theta \sin \theta_0 \right] \\ \\ &= \frac{\mu r}{h^2} \left( e \cos \theta + \cos⁡(\theta - \theta_0 ) \right) \\ \\ &= \frac{\mu r}{h^2} (e \cos \theta + \cos \Delta \theta ) \end{align} \]

 

여기서 \(\Delta \theta = \theta- \theta_0\) 이다. 한편

 

\[ r= \frac{h^2}{\mu} \frac{1}{1+e \cos \theta} \tag{17} \]

 

에서 \(e \cos \theta = \frac{h^2}{\mu r}-1\) 이기 때문에 식 (16)은 다음과 같이 정리된다.

 

\[ f=1- \frac{\mu r}{h^2} (1- \cos \Delta \theta ) \tag{18} \]

 

같은 방법으로 \(g\) 를 계산해 보면 다음과 같다.

 

\[ \begin{align} g &= \frac{-xy_0+yx_0}{h} \tag{19} \\ \\ &= \frac{1}{h} \left[ (-r \cos \theta ) r_0 \sin \theta_0 +r \sin \theta r_0 \cos \theta_0 \right] \\ \\ &= \frac{rr_0}{h} (- \cos \theta \sin \theta_0 + \sin \theta \cos \theta_0 ) \\ \\ &= \frac{rr_0}{h} \sin⁡ (\theta - \theta_0 ) \\ \\ &= \frac{rr_0}{h} \sin \Delta \theta \end{align} \]

 

\(\dot{g}\) 은 다음과 같이 계산된다.

 

\[ \begin{align} \dot{g} &= \frac{-\dot{x}y_0+ \dot{y}x_0}{h} \tag{20} \\ \\ &= \frac{1}{h} \left[ \frac{\mu}{h} \sin \theta r_0 \sin \theta_0 + \frac{\mu}{h} (e+ \cos \theta ) r_0 \cos⁡ \theta_0 \right] \\ \\ &= \frac{ \mu r_0}{h^2} (\sin \theta \sin \theta_0 + \cos \theta \cos \theta_0 +e \cos \theta_0 ) \\ \\ &= \frac{\mu r_0 }{h^2} \left( \cos (\theta - \theta_0 )+ \left( \frac{h^2}{μr_0}-1 \right) \right) \\ \\ &= 1- \frac{ \mu r_0}{h^2} (1- \cos \Delta \theta ) \end{align} \]

 

마지막으로 \(\dot{f}\) 은 식 (14)를 이용하여 계산할 수 있다.

 

\[ \begin{align} \dot{f} &= \frac{1}{g} (f \dot{g}-1) \tag{21} \\ \\ &= \frac{h}{rr_0 \sin \Delta \theta } \left[ \left( 1- \frac{\mu r}{h^2} (1- \cos \Delta \theta ) \right) \left( 1- \frac{\mu r_0 }{h^2} (1- \cos \Delta \theta ) \right) -1 \right] \\ \\ &= \frac{h}{rr_0 \sin \Delta \theta } \left[ - \frac{\mu r_0 }{h^2} (1- \cos \Delta \theta )- \frac{\mu r}{h^2} (1- \cos \Delta \theta )+ \frac{\mu^2 rr_0 }{h^4} (1- \cos \Delta \theta )^2 \right] \\ \\ &= \frac{\mu}{h} \frac{(1- \cos \Delta \theta )}{\sin \Delta \theta} \left[\frac{\mu}{h^2} (1-\cos \Delta \theta )- \frac{1}{r_0} -\frac{1}{r} \right] \end{align} \]

 

식 (18)~(21)에 \(r\) 이 있으므로 이것도 \(\vec{r}_0, \vec{v}_0\) 와 \(\Delta \theta \) 의 함수로 표현할 필요가 있다. 식 (17)로부터 다음식을 얻을 수 있다.

 

\[ \begin{align} r &= \frac{h^2}{\mu} \frac{1}{1+e \cos⁡ (\theta_0+ \Delta \theta) } \tag{22} \\ \\ &= \frac{h^2}{\mu} \frac{1}{1+e \cos \theta_0 \cos \Delta \theta -e \sin \theta_0 \sin \Delta \theta} \end{align} \]

 

여기서 \(e \cos \theta_0 =\frac{h^2}{\mu r_0 }-1\) 이고, \(\dot{r}= \frac{\mu}{h} e \sin \theta\) 이므로(https://pasus.tistory.com/288), \(r\dot{r}= \vec{r} \cdot \vec{v}\) 의 관계식을 이용하면,

 

\[ e \sin \theta_0 = \frac{h}{\mu} \dot{r}_0= \frac{h}{\mu} \frac{\vec{r}_0 \cdot \vec{v}_0}{r_0} \tag{23} \]

 

이 된다. 이를 식 (22)에 대입하면 \(r\) 은 다음과 같이 된다.

 

\[ r = \frac{h^2}{\mu} \frac{1}{ 1+\left( \frac{h^2}{\mu r_0 }-1 \right) \cos \Delta \theta -\left( \frac{h}{\mu} \frac{\vec{r}_0 \cdot \vec{v}_0}{r_0} \right) \sin \Delta \theta} \tag{24} \]

 

정리하면, 어느 우주비행체의 초기 위치벡터 \(\vec{r}_0\) 와 속도벡터 \(\vec{v}_0\) 가 주어졌을 때, 실제 비행각이 \(\Delta \theta\) 만큼 변화한 후의 위치벡터와 속도벡터 \(\vec{r}, \vec{v}\) 를 구하는 방법은 다음과 같다.

1. \(\vec{r}_0\) 와 \(\vec{v}_0\) 로부터 \(r_0\) 와 각운동량 \(h\) 를 계산한다.
2. 식 (24)로 \(r\) 을 계산한다.
3. 식 (18)~(21)로 라그랑지 계수를 계산한다.
4. 식 (1)로 \(\vec{r}, \vec{v}\) 를 계산한다.

 

 

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

 

function [r_vec, v_vec] = kepler_the(r0_vec, v0_vec, d_the)

% [r_vec, v_vec] = kepler_the(r0_vec, v0_vec, d_the)
%        given r0_vec, r0_vec and true anomaly increments
%        compute r_vec and v_vec
%       
% Input  r0_vec: initial position vector in km
%        v0_vec: initial velocity vector in km/sec
%        d_the: true anomaly change (degrees)
% 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
r0 = sqrt(r0_vec'*r0_vec);
h_vec = cross(r0_vec, v0_vec);
h = sqrt(h_vec'*h_vec);

% step 2
d_the = d_the * pi/180;
r_tmp = 1 + ( h^2/(mu_e*r0) - 1 )*cos(d_the) ...
         - (h/mu_e)*r0_vec'*v0_vec/r0*sin(d_the);
r = h^2/mu_e/r_tmp;

% step 3 : Lagrange coefficients
f = 1 - mu_e*r/h^2 * (1-cos(d_the));
g = r*r0/h * sin(d_the);
f_dot = mu_e/h*(1-cos(d_the))/sin(d_the) * (mu_e/h^2*(1-cos(d_the))-1/r0-1/r);
g_dot = 1 - mu_e*r0/h^2*(1-cos(d_the));

% step 4
r_vec = f*r0_vec + g*v0_vec;
v_vec = f_dot*r0_vec + g_dot*v0_vec;