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

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

 

 

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

 

$$\begin{array}{r}\text{(1)}& & \overrightarrow{r}=f{\overrightarrow{r}}_0+g{\overrightarrow{v}}_0\\ & \overrightarrow{v}=\dot{f}{\overrightarrow{r}}_0+\dot{g}{\overrightarrow{v}}_0\end{array}$$

 

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

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

 

$$\begin{array}{r}\text{(2)}& & \overrightarrow{r}=x{\hat{p}}_1+y{\hat{p}}_2\\ & \overrightarrow{v}=\dot{x}{\hat{p}}_1+\dot{y}{\hat{p}}_2\end{array}$$

 

여기서

 

$$\begin{array}{r}\text{(3)}& & x=r\cos \theta ,y=r\sin \theta \\ & \dot{x}=-\frac{\mu }{h}\sin \theta ,\dot{y}=\frac{\mu }{h}(e+\cos \theta )\end{array}$$

 

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

 

$$\begin{array}{r}\text{(4)}& & {\overrightarrow{r}}_0=x_0{\hat{p}}_1+y_0{\hat{p}}_2\\ & {\overrightarrow{v}}_0={\dot{x}}_0{\hat{p}}_1+{\dot{y}}_0{\hat{p}}_2\end{array}$$

 

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

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

 

$$\begin{array}{}\text{(5)}& {\hat{p}}_2=\frac{1}{y_0}{\overrightarrow{r}}_0-\frac{x_0}{y_0}{\hat{p}}_1\end{array}$$

 

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

 

$$\begin{array}{r}\text{(6)}& {\overrightarrow{v}}_0& ={\dot{x}}_0{\hat{p}}_1+{\dot{y}}_0(\frac{1}{y_0}{\overrightarrow{r}}_0-\frac{x_0}{y_0}{\hat{p}}_1)\\ & =\left(\frac{{\dot{x}}_0{y}_0-{\dot{y}}_0{x}_0}{y_0}\right){\hat{p}}_1+\frac{{\dot{y}}_0}{y_0}{\overrightarrow{r}}_0\end{array}$$

 

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

 

$$\begin{array}{r}\text{(7)}& {\overrightarrow{h}}_0& ={\overrightarrow{r}}_0\times {\overrightarrow{v}}_0=(x_0{\dot{y}}_0-y_0{\dot{x}}_0){\hat{p}}_3\\ & =h{\hat{p}}_3\end{array}$$

 

여기서

 

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

 

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

 

$$\begin{array}{}\text{(8)}& {\overrightarrow{v}}_0=-\frac{h}{y_0}{\hat{p}}_1+\frac{{\dot{y}}_0}{y_0}{\overrightarrow{r}}_0\end{array}$$

 

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

 

$$\begin{array}{}\text{(9)}& {\hat{p}}_1=\frac{{\dot{y}}_0}{h}{\overrightarrow{r}}_0-\frac{y_0}{h}{\overrightarrow{v}}_0\end{array}$$

 

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

 

$$\begin{array}{r}\text{(10)}& {\hat{p}}_2& =\frac{1}{y_0}{\overrightarrow{r}}_0-\frac{x_0}{y_0}(\frac{{\dot{y}}_0}{h}{\overrightarrow{r}}_0-\frac{y_0}{h}{\overrightarrow{v}}_0)\\ & =\frac{h-x_0{\dot{y}}_0}{y_{0}h}{\overrightarrow{r}}_0+\frac{x_0}{h}{\overrightarrow{v}}_0\\ \\ & =-\frac{{\dot{x}}_0}{h}{\overrightarrow{r}}_0+\frac{x_0}{h}{\overrightarrow{v}}_0\end{array}$$

 

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

 

$$\begin{array}{r}\text{(11)}& \overrightarrow{r}& =x(\frac{{\dot{y}}_0}{h}{\overrightarrow{r}}_0-\frac{y_0}{h}{\overrightarrow{v}}_0)+y(-\frac{{\dot{x}}_0}{h}{\overrightarrow{r}}_0+\frac{x_0}{h}{\overrightarrow{v}}_0)\\ & =\frac{x{\dot{y}}_0-y{\dot{x}}_0}{h}{\overrightarrow{r}}_0+\frac{-x{y}_0+y{x}_0}{h}{\overrightarrow{v}}_0\\ \\ \overrightarrow{v}& =\dot{x}(\frac{{\dot{y}}_0}{h}{\overrightarrow{r}}_0-\frac{y_0}{h}{\overrightarrow{v}}_0)+\dot{y}(-\frac{{\dot{x}}_0}{h}{\overrightarrow{r}}_0+\frac{x_0}{h}{\overrightarrow{v}}_0)\\ \\ & =\frac{\dot{x}{\dot{y}}_0-\dot{y}{\dot{x}}_0}{h}{\overrightarrow{r}}_0+\frac{-\dot{x}{y}_0+\dot{y}{x}_0}{h}{\overrightarrow{v}}_0\end{array}$$

 

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

 

$$\begin{array}{r}\text{(12)}& f& =\frac{x{\dot{y}}_0-y{\dot{x}}_0}{h},g=\frac{-x{y}_0+y{x}_0}{h}\\ \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{array}$$

 

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

 

$$\begin{array}{r}\text{(13)}& \overrightarrow{h}& =\overrightarrow{r}\times \overrightarrow{v}=(f{\overrightarrow{r}}_0+g{\overrightarrow{v}}_0)\times (\dot{f}{\overrightarrow{r}}_0+\dot{g}{\overrightarrow{v}}_0)\\ & =(f\dot{g}-\dot{f}g)({\overrightarrow{r}}_0\times {\overrightarrow{v}}_0)\\ \\ & =(f\dot{g}-\dot{f}g){\overrightarrow{h}}_0\end{array}$$

 

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

 

$$\begin{array}{}\text{(14)}& f\dot{g}-\dot{f}g=1\end{array}$$

 

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

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

 

$$\begin{array}{r}\text{(15)}& & x_0=r_0\cos {\theta }_0,y_0=r_0\sin {\theta }_0\\ & {\dot{x}}_0=-\frac{\mu }{h}\sin {\theta }_0,{\dot{y}}_0=\frac{\mu }{h}(e+\cos {\theta }_0)\end{array}$$

 

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

 

 

먼저 $f$ 를 계산해 보자.

 

$$\begin{array}{r}\text{(16)}& f& =\frac{x{\dot{y}}_0-y{\dot{x}}_0}{h}\\ & =\frac{1}{h}[r\cos \theta \frac{\mu }{h}(e+\cos {\theta }_0)+r\sin \theta \frac{\mu }{h}\sin {\theta }_0]\\ \\ & =\frac{\mu r}{h^2}[e\cos \theta +\cos \theta \cos {\theta }_0+\sin \theta \sin {\theta }_0]\\ \\ & =\frac{\mu r}{h^2}(e\cos \theta +\cos (\theta -{\theta }_0))\\ \\ & =\frac{\mu r}{h^2}(e\cos \theta +\cos \mathrm{\Delta }\theta )\end{array}$$

 

여기서 $\mathrm{\Delta }\theta =\theta -{\theta }_0$ 이다. 한편

 

$$\begin{array}{}\text{(17)}& r=\frac{h^2}{\mu }\frac{1}{1+e\cos \theta }\end{array}$$

 

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

 

$$\begin{array}{}\text{(18)}& f=1-\frac{\mu r}{h^2}(1-\cos \mathrm{\Delta }\theta )\end{array}$$

 

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

 

$$\begin{array}{r}\text{(19)}& g& =\frac{-x{y}_0+y{x}_0}{h}\\ & =\frac{1}{h}[(-r\cos \theta )r_0\sin {\theta }_0+r\sin \theta r_0\cos {\theta }_0]\\ \\ & =\frac{r{r}_0}{h}(-\cos \theta \sin {\theta }_0+\sin \theta \cos {\theta }_0)\\ \\ & =\frac{r{r}_0}{h}\sin (\theta -{\theta }_0)\\ \\ & =\frac{r{r}_0}{h}\sin \mathrm{\Delta }\theta \end{array}$$

 

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

 

$$\begin{array}{r}\text{(20)}& \dot{g}& =\frac{-\dot{x}{y}_0+\dot{y}{x}_0}{h}\\ & =\frac{1}{h}[\frac{\mu }{h}\sin \theta r_0\sin {\theta }_0+\frac{\mu }{h}(e+\cos \theta )r_0\cos {\theta }_0]\\ \\ & =\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}(\cos (\theta -{\theta }_0)+(\frac{h^2}{\mu r_0}-1))\\ \\ & =1-\frac{\mu r_0}{h^2}(1-\cos \mathrm{\Delta }\theta )\end{array}$$

 

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

 

$$\begin{array}{r}\text{(21)}& \dot{f}& =\frac{1}{g}(f\dot{g}-1)\\ & =\frac{h}{r{r}_0\sin \mathrm{\Delta }\theta }[(1-\frac{\mu r}{h^2}(1-\cos \mathrm{\Delta }\theta ))(1-\frac{\mu r_0}{h^2}(1-\cos \mathrm{\Delta }\theta ))-1]\\ \\ & =\frac{h}{r{r}_0\sin \mathrm{\Delta }\theta }[-\frac{\mu r_0}{h^2}(1-\cos \mathrm{\Delta }\theta )-\frac{\mu r}{h^2}(1-\cos \mathrm{\Delta }\theta )+\frac{{\mu }^{2}r{r}_0}{h^4}(1-\cos \mathrm{\Delta }\theta {)}^2]\\ \\ & =\frac{\mu }{h}\frac{(1-\cos \mathrm{\Delta }\theta )}{\sin \mathrm{\Delta }\theta }[\frac{\mu }{h^2}(1-\cos \mathrm{\Delta }\theta )-\frac{1}{r_0}-\frac{1}{r}]\end{array}$$

 

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

 

$$\begin{array}{r}\text{(22)}& r& =\frac{h^2}{\mu }\frac{1}{1+e\cos ({\theta }_0+\mathrm{\Delta }\theta )}\\ & =\frac{h^2}{\mu }\frac{1}{1+e\cos {\theta }_0\cos \mathrm{\Delta }\theta -e\sin {\theta }_0\sin \mathrm{\Delta }\theta }\end{array}$$

 

여기서 $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}=\overrightarrow{r}\cdot \overrightarrow{v}$ 의 관계식을 이용하면,

 

$$\begin{array}{}\text{(23)}& e\sin {\theta }_0=\frac{h}{\mu }{\dot{r}}_0=\frac{h}{\mu }\frac{{\overrightarrow{r}}_0\cdot {\overrightarrow{v}}_0}{r_0}\end{array}$$

 

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

 

$$\begin{array}{}\text{(24)}& r=\frac{h^2}{\mu }\frac{1}{1+(\frac{h^2}{\mu r_0}-1)\cos \mathrm{\Delta }\theta -\left(\frac{h}{\mu }\frac{{\overrightarrow{r}}_0\cdot {\overrightarrow{v}}_0}{r_0}\right)\sin \mathrm{\Delta }\theta }\end{array}$$

 

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

1. ${\overrightarrow{r}}_0$ 와 ${\overrightarrow{v}}_0$ 로부터 $r_0$ 와 각운동량 $h$ 를 계산한다.
2. 식 (24)로 $r$ 을 계산한다.
3. 식 (18)~(21)로 라그랑지 계수를 계산한다.
4. 식 (1)로 $\overrightarrow{r},\overrightarrow{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;