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

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

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

 

 

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

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

 

 

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

 

$$\begin{array}{}\text{(1)}& \overrightarrow{r}=r{\hat{e}}_r\end{array}$$

 

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

 

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

 

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

 

$$\begin{array}{r}\text{(3)}& \overrightarrow{v}=\frac{{}^{i}d\overrightarrow{r}}{dt}& =\frac{{}^{p}d\overrightarrow{r}}{dt}\\ & =\frac{{}^{q}d\overrightarrow{r}}{dt}+{}^p{\overrightarrow{\omega }}^q\times \overrightarrow{r}\end{array}$$

 

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

 

$$\begin{array}{}\text{(4)}& \overrightarrow{v}=\dot{r}{\hat{e}}_r+r\dot{\theta }{\hat{e}}_{\theta }\end{array}$$

 

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

 

$$\begin{array}{}\text{(5)}& r\dot{\theta }=\frac{h}{r}\end{array}$$

 

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

 

$$\begin{array}{}\text{(6)}& r\dot{\theta }=\sqrt{\frac{\mu }{p}}(1+e\cos \theta )\end{array}$$

 

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

 

$$\begin{array}{}\text{(7)}& \dot{r}=\frac{pe\dot{\theta }\sin \theta }{(1+e\cos \theta {)}^2}=\frac{r\dot{\theta }e\sin \theta }{1+e\cos \theta }\end{array}$$

 

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

 

$$\begin{array}{}\text{(8)}& \dot{r}=\sqrt{\frac{\mu }{p}}e\sin \theta \end{array}$$

 

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

 

$$\begin{array}{}\text{(9)}& \overrightarrow{v}=\sqrt{\frac{\mu }{p}}e\sin \theta {\hat{e}}_r+\sqrt{\frac{\mu }{p}}(1+e\cos \theta ){\hat{e}}_{\theta }\end{array}$$

 

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

 

$$\begin{array}{}\text{(10)}& C_q^p=C(z,\theta )=\left[\begin{array}{ccc}\cos \theta & -\sin \theta & 0\\ \sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right]\end{array}$$

 

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

 

$$\begin{array}{r}\text{(11)}& \overrightarrow{r}& =r\cos \theta {\hat{p}}_1+r\sin \theta {\hat{p}}_2\\ \overrightarrow{v}& =-\sqrt{\frac{\mu }{p}}\sin \theta {\hat{p}}_1+\sqrt{\frac{\mu }{p}}(e+\cos \theta ){\hat{p}}_2\end{array}$$

 

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

 

 

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

 

 

따라서

 

$$\begin{array}{r}\text{(12)}& C_p^i& =C(z,\mathrm{\Omega })C(x,i)C(z,\omega )\\ & =\left[\begin{array}{ccc}c\mathrm{\Omega }c\omega -s\mathrm{\Omega }cis\omega & -c\mathrm{\Omega }s\omega -s\mathrm{\Omega }cic\omega & s\mathrm{\Omega }si\\ s\mathrm{\Omega }c\omega +c\mathrm{\Omega }cis\omega & -s\mathrm{\Omega }s\omega +c\mathrm{\Omega }cic\omega & -c\mathrm{\Omega }si\\ sis\omega & sic\omega & ci\end{array}\right]\end{array}$$

 

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

 

$$\begin{array}{r}\text{(13)}& {\mathbf{r}}^i& =C_p^i{\mathbf{r}}^p\\ {\mathbf{v}}^i& =C_p^i{\mathbf{v}}^p\end{array}$$

 

여기서 벡터 ${\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{array}{r}\text{(14)}& \overrightarrow{r}& =r_1{\hat{i}}_1+r_2{\hat{i}}_2+r_3{\hat{i}}_3\\ \overrightarrow{v}& =v_1{\hat{i}}_1+v_2{\hat{i}}_2+v_3{\hat{i}}_3\end{array}$$

 

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

 

$$\begin{array}{rl}& a=8788.1km,e=0.1712,i={153.25}^0\\ \\ & \mathrm{\Omega }={255.30}^0,\omega ={20.07}^0,\theta ={28.45}^0\end{array}$$

 

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

 

$$\begin{array}{rl}\overrightarrow{r}& =-6044.2{\hat{i}}_1-3491.6{\hat{i}}_2+2500.2{\hat{i}}_3(km)\\ \\ \overrightarrow{v}& =-3.4587{\hat{i}}_1+6.6171{\hat{i}}_2+2.5326{\hat{i}}_3(km/sec)\end{array}$$

 

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

 

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