궤도요소 (COE)로 부터 위치 및 속도벡터 계산

우주비행체의 위치벡터 및 속도벡터를 궤도요소(COE, classical orbital elements)로 변환할 수 있었다 (https://pasus.tistory.com/287). 이번에는 이와 반대로 궤도요소를 위치벡터와 속도벡터로 변환하는 방법에 대해서 알아보기로 하자.

시간 \(t=t_0\) 에서 주어진 궤도요소 \((a, \ e, \ i, \ \Omega, \ \omega, \ \theta (t_0 ))\) 에서 위치벡터 \(\vec{r}\) 과 속도벡터 \(\vec{v}\) 를 구하는 과정은 두 단계로 나누어진다. 궤도중심좌표계(perifocal frame)에서 위치벡터와 속도벡터를 구하는 단계와 좌표변환을 통하여 ECI좌표계로 이들 벡터를 변환하는 단계이다.

 

 

먼저 궤도중심좌표계에서 위치벡터와 속도벡터를 구하는 단계이다.

다음 그림과 같이 \((\hat{e}_r, \ \hat{e}_\theta , \ \hat{e}_z )\) 를 기저벡터로 하는 새로운 극좌표계 \(\{q\}\) 를 도입하기로 한다. \(\hat{e}_r\) 축은 벡터 \(\vec{r}\) 을 향하며 \(\hat{e}_z\) 는 궤도중심좌표계 \(\{p\}\) 의 \(\hat{p}_3\) 와 동일하다.

 

 

극좌표계 \(\{q\}\) 에서 우주비행체의 위치벡터를 표현하면 다음과 같다.

 

\[ \vec{r}=r \ \hat{e}_r \tag{1} \]

 

여기서 거리 \(r\) 은 다음과 같이 궤도요소의 함수로 표현된다.

 

\[ r= \frac{p}{1+e \cos \theta } \tag{2} \]

 

속도벡터는 식 (1)을 ECI좌표계에서 미분하면 된다. 궤도중심좌표계 \(\{p\}\) 는 ECI좌표계 \(\{i\}\) 대해서 고정이므로 속도벡터는 기본운동학방정식(BKE)을 이용하면 다음과 같이 계산된다 (https://pasus.tistory.com/121).

 

\[ \begin{align} \vec{v} = \frac{ ^i d \vec{r}}{dt} & = \frac{ ^p d \vec{r} }{dt} \tag{3} \\ \\ &= \frac{^q d \vec{r} }{dt} + {^p {\vec{\omega}}^q } \times \vec{r} \end{align} \]

 

여기서 \( \ ^p {\vec{\omega}^q} = \dot{\theta} \ \hat{p}_3= \dot{\theta} \ \hat{e}_z\) 이므로 위 식은 다음과 같이 된다.

 

\[ \vec{v}= \dot{r} \ \hat{e}_r+ r \dot{\theta} \ \hat{e}_\theta \tag{4} \]

 

한편 각운동량 벡터는 \(\vec{h}= \vec{r} \times \vec{v}\) 이므로 식 (1)과 (4)를 이용하면 \(\vec{h}=r^2 \dot{\theta} \ \hat{p}_3\) 이기 때문에 각운동량의 크기는 \(h=r^2 \dot{\theta}\) 가 된다. 따라서

 

\[ r \dot{\theta} = \frac{h}{r} \tag{5} \]

 

가 된다. \(h= \sqrt{\mu p}\) 의 관계식과 식 (2)를 위 식에 대입하면

 

\[ r \dot{\theta} = \sqrt{ \frac{\mu}{p} } (1+ e \cos \theta ) \tag{6} \]

 

가 된다. 식 (2)를 미분하면

 

\[ \dot{r} = \frac{ pe \dot{\theta} \sin \theta }{(1+ e \cos \theta )^2} = \frac{r \dot{\theta} e \sin \theta}{1+e \cos \theta } \tag{7} \]

 

가 되므로 식 (6)을 위 식에 대입하면

 

\[ \dot{r}= \sqrt{ \frac{ \mu}{p} } e \sin \theta \tag{8} \]

 

가 된다. 식 (6)과 (8)을 식 (4)에 대입하면 속도벡터는

 

\[ \vec{v} = \sqrt{ \frac{\mu}{p} } e \sin \theta \ \hat{e}_r+ \sqrt{ \frac{\mu}{p} } (1+ e \cos \theta ) \ \hat{e}_\theta \tag{9} \]

 

가 된다. 한편 좌표계 \(\{p\}\) 에서 \(\{q\}\) 로의 방향코사인행렬 (DCM)은

 

\[ C_q^p = C( z, \theta) = \begin{bmatrix} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{bmatrix} \tag{10} \]

 

이므로 식 (1)과 (9)를 궤도중심좌표계로 표현하면 다음과 같다 (https://pasus.tistory.com/83).

 

\[ \begin{align} \vec{r} &= r \cos \theta \ \hat{p}_1+ r \sin \theta \ \hat{p}_2 \tag{11} \\ \\ \vec{v} &= - \sqrt{ \frac{\mu}{p} } \sin \theta \ \hat{p}_1+ \sqrt{ \frac{\mu}{p} } (e+ \cos \theta ) \ \hat{p}_2 \end{align} \]

 

식 (11)은 우주비행체의 위치벡터와 속도벡터를 궤도중심좌표계로 표현한 식이다.

 

 

이제 식 (11)을 ECI좌표계로 변환하면 된다. ECI좌표계 \(\{i\}\) 에서 궤도중심좌표계 \(\{p\}\) 로의 DCM은 다음 그림과 같이 연쇄적인 좌표축 회전으로 표현할 수 있다. 우선 ECI좌표계의 \(\hat{i}_3\) 축으로 \(\Omega\) 만큼 회전한 후, 돌아간 좌표계의 \(\hat{i}'_1\) 축으로 \(i\) 만큼 회전하고, 다시 돌아간 좌표계의 \(\hat{i}''_3\) 축으로 \(\omega\) 만큼 회전하면 궤도중심좌표계로의 변환이 완성된다.

 

 

따라서

 

\[ \begin{align} C_p^i &= C(z, \Omega) \ C(x,i) \ C(z, \omega) \tag{12} \\ \\ &= \begin{bmatrix} c \Omega c \omega - s \Omega c ⁡i s \omega & -c \Omega s \omega - s \Omega c i c \omega & s \Omega s i \\ s \Omega c \omega + c \Omega c⁡ i s \omega & -s \Omega s \omega + c \Omega c i c \omega & -c \Omega s i \\ s i s \omega & s i c \omega & c i \end{bmatrix} \end{align} \]

 

가 된다. 여기서 \(c= \cos , \ s= \sin \) 을 의미한다. 식 (11)과 (12)를 이용하여 위치벡터와 속도벡터를 ECI좌표계로 표현하면 다음과 같다.

 

\[ \begin{align} \mathbf{r}^i &= C_p^i \ \mathbf{r}^p \tag{13} \\ \\ \mathbf{v}^i &= C_p^i \ \mathbf{v}^p \end{align} \]

 

여기서 벡터 \(\mathbf{r}^p\) 와 \(\mathbf{v}^p\) 의 성분은 식 (11)의 축 성분이며 \(\mathbf{r}^i=[r_1 \ \ r_2 \ \ r_3 ]^T\) 와 \(\mathbf{v}^i=[v_1 \ \ v_2 \ \ v_3 ]^T\) 는 식 (14)의 축 성분을 성분으로 하는 벡터이다.

 

\[ \begin{align} \vec{r} &= r_1 \ \hat{i}_1+r_2 \ \hat{i}_2+r_3 \ \hat{i}_3 \tag{14} \\ \\ \vec{v} &=v_1 \ \hat{i}_1+v_2 \ \hat{i}_2+v_3 \ \hat{i}_3 \end{align} \]

 

예를 들어서, 우주비행체의 궤도요소 \((a,\ e, \ i, \ \Omega, \ \omega , \ \theta )\) 가 다음과 같이 주어졌다면,

 

\[ \begin{align} & a=8788.1 \ km, \ \ \ \ \ e=0.1712, \ \ \ \ \ i=153.25^0 \\ \\ & \Omega=255.30^0, \ \ \ \ \ \omega=20.07^0, \ \ \ \ \ \theta=28.45^0 \end{align} \]

 

위치벡터와 속도벡터는 각각 다음과 같이 계산된다.

 

\[ \begin{align} \vec{r} &= -6044.2 \ \hat{i}_1-3491.6 \ \hat{i}_2+2500.2 \ \hat{i}_3 \ (km) \\ \\ \vec{v} &= -3.4587 \ \hat{i}_1+6.6171 \ \hat{i}_2+2.5326 \ \hat{i}_3 \ (km / sec) \end{align} \]

 

다음은 궤도요소를 위치벡터와 속도벡터로 변환해 주는 매트랩 코드다.

 

function [r,v] = coe2rv(a,e,i,W,w,theta)

% [r,v]=coe2rv(a,e,i,W,w,theta)
% Converting orbital elements to vectors R and V in inertial frame,
%        where V=dR/dt in inertial frame.
% Input  (a, e, i, W, w, theta) in (km, sec) and degrees (p=a(1-e^2))
%        put undefined values as 0.
% Output (R, V) in km, km/s
%
% coded by st.watermelon

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

p = a*(1-e^2);
theta = theta*pi/180;
sthe = sin(theta); 
cthe = cos(theta);
mup = sqrt(mu/p);

s_r = p/(1+e*cthe);
r_o = [s_r*cthe s_r*sthe 0]'; % in orbital frame
v_o = [-(mup*sthe) (mup*(e+cthe)) 0]';

C = eci2perifocal(i,W,w);
r = C*r_o;
v = C*v_o;

end

% DCM ECI to Perifocal
function C = eci2perifocal(i,W,w)

i = i*pi/180;
W = W*pi/180;
w = w*pi/180;

ci = cos(i); si = sin(i);
cW = cos(W); sW = sin(W);
cw = cos(w); sw = sin(w);

C1 = [cW -sW 0; sW cW 0; 0 0 1];
C2 = [1 0 0; 0 ci -si; 0 si ci];
C3 = [cw -sw 0; sw cw 0; 0 0 1];

C = C1*C2*C3;
end