궤도요소 (COE) 계산

고전 궤도요소 (COE, classical orbital elements)의 6개 파라미터는 우주비행체의 위치벡터 및 속도벡터와 함수관계에 있다. 따라서 임의의 시간 $t=t_0$ 에서 주어진 위치벡터와 속도벡터를 궤도요소로 변환하면 궤도의 크기, 모양, 자세 등을 알 수 있다. 시간 $t=t_0$ 에서 주어진 위치벡터와 속도벡터 $\overrightarrow{r}(t_0),\overrightarrow{v}(t_0)$ 에서 궤도요소를 구하는 과정은 두 단계로 나누어진다.

우선 위치벡터와 속도벡터로부터 각운동량 벡터 $\overrightarrow{h}$, 승교선 벡터(ascending node vector) $\overrightarrow{n}$, 이심율 벡터(eccentricity vector) $\overrightarrow{e}$ 를 구하는 단계와 이들 벡터로부터 궤도요소를 구하는 단계이다.

 

 

다음은 위치벡터와 속도벡터로부터 각운동량 벡터, 승교선 벡터, 이심율 벡터를 구하는 방법에 대한 설명이다.

각운동량 벡터 $\overrightarrow{h}$:
각운동량 벡터 $\overrightarrow{h}$ 는 정의에 의해서

 

$$\begin{array}{}\text{(1)}& \overrightarrow{h}=\overrightarrow{r}(t_0)\times \overrightarrow{v}(t_0)\end{array}$$

 

로 구하면 된다. 각운동량 벡터는 전 궤도상에서 상수로 유지되므로 모든 $t\ge t_0$ 에서 같다. 즉, $\overrightarrow{h}=\overrightarrow{r}(t)\times \overrightarrow{v}(t),t\ge t_0$ 이다.

승교선 벡터 $\overrightarrow{n}$:
승교선 벡터는 적도면과 궤도면의 교차선과 평행하며 크기는 각운동량과 같고 방향은 지구중심에서 우주비행체가 지구의 남반구에서 북반구로 올라가면서 만나는 점을 향하는 벡터이다. 적도면과 직각인 벡터는 ${\hat{i}}_3$ 이고 궤도면과 직각인 벡터는 $\overrightarrow{h}$ 이므로 승교선 벡터는 두 벡터에 모두 직각이 되도록 다음과 같이 정의한다.

 

$$\begin{array}{}\text{(2)}& \overrightarrow{n}={\hat{i}}_3\times \overrightarrow{h}\end{array}$$

 

$\overrightarrow{h}$ 가 상수벡터이므로 승교선 벡터도 상수벡터이다.

 

 

이심율 벡터 $\overrightarrow{e}$:
이심율 벡터는 크기가 이심율 $e$ 과 같고 방향은 근지점을 향하는, 즉 궤도중심좌표계(perifocal frame)의 ${\hat{p}}_1$ 축과 같은 방향을 가지는 벡터로 정의한다. 궤도중심좌표계 $\{p\}$ 는 우주비행체의 위치를 궤도상에 표시하기 위한 기준좌표계로서 지구의 중심에 원점이 위치한다. ${\hat{p}}_1$ 축은 근지점을 향하며, ${\hat{p}}_3$ 는 궤도면의 직각 방향 또는 각운동량 벡터 방향을 향하고 ${\hat{p}}_2$ 는 오른손 법칙에 의해서 정해진다. 원궤도의 경우 근지점이 없으므로 승교선 벡터 방향을 ${\hat{p}}_1$ 축으로 정의하며, 적도궤도이면서 원궤도인 경우는 ECI좌표계의 ${\hat{i}}_1$ 축을 ${\hat{p}}_1$ 축으로 정의한다. 궤도중심좌표계는 ECI좌표계에 고정되어 있으므로 관성좌표계의 하나로 볼 수 있다.

 

 

궤적방정식 유도 과정 (https://pasus.tistory.com/102)에 의하면 다음 식이 있었다.

 

$$\begin{array}{}\text{(3)}& \overrightarrow{v}\times \overrightarrow{h}=\mu \frac{\overrightarrow{r}}{r}+\overrightarrow{B}\end{array}$$

 

여기서 $\overrightarrow{B}$ 는 적분상수벡터이다. 실제비행각 $\theta$ 는 위치벡터 $\overrightarrow{r}$ 과 적분상수벡터 $\overrightarrow{B}$ 의 사잇각이므로 (https://pasus.tistory.com/171) $\overrightarrow{B}$ 는 근지점을 향하는, 또는 궤도중심좌표계의 ${\hat{p}}_1$ 과 같은 방향을 가지는 벡터이다. 이심율은 $e=\frac{B}{\mu }$ 이므로 이심율 벡터는

 

$$\begin{array}{}\text{(4)}& \overrightarrow{e}=\frac{\overrightarrow{B}}{\mu }\end{array}$$

 

로 놓을 수 있다. $\overrightarrow{B}$ 가 상수벡터이므로 이심율 벡터도 상수벡터임을 알 수 있다. 식 (3)에서 $\overrightarrow{B}$ 를 구하고 식 (4)에 대입하면

 

$$\begin{array}{}\text{(5)}& \overrightarrow{e}=\frac{1}{\mu }(\overrightarrow{v}\times \overrightarrow{h}-\mu \frac{\overrightarrow{r}}{r})\end{array}$$

 

가 된다. $\overrightarrow{h}=\overrightarrow{r}\times \overrightarrow{v}$ 이고, $\overrightarrow{a}\times (\overrightarrow{b}\times \overrightarrow{c})=(\overrightarrow{a}\cdot \overrightarrow{c})\overrightarrow{b}-(\overrightarrow{a}\cdot \overrightarrow{b})\overrightarrow{c}$ 의 관계를 이용하여 오른쪽 항을 다시 쓰면,

 

$$\begin{array}{}\text{(6)}& \overrightarrow{e}=\frac{1}{\mu }((v^2-\frac{\mu }{r})\overrightarrow{r}-(\overrightarrow{r}\cdot \overrightarrow{v})\overrightarrow{v})\end{array}$$

 

이 된다.

 

 

다음은 각운동량 벡터 $\overrightarrow{h}$, 승교선 벡터 $\overrightarrow{n}$, 이심율 벡터 $\overrightarrow{e}$ 로부터 궤도요소를 구하는 단계이다.

통반경 $p$ 와 이심율 $e$:
이심율은 이심율 벡터의 크기이고, 통반경은 정의에 의하여 구하면 된다. 장반경은 포물선궤도를 제외하고는 $p=a(1-e^2)$ 의 관계가 성립하므로 다음과 같이 구할 수 있다.

 

$$\begin{array}{r}\text{(7)}& & e=|\overrightarrow{e}|\\ & p=\frac{h^2}{\mu },a=\frac{p}{1-e^2}\end{array}$$

 

여기서 $h=|\overrightarrow{h}|$ 이다.

 

 

경사각 $i$:
경사각은 ${\hat{i}}_3$ 와 $\overrightarrow{h}$ 의 사잇각이므로 ${\hat{i}}_3$ 와 $\overrightarrow{h}$ 의 내적을 이용하여 구할 수 있다. 경사각의 범위는 $0$도에서 $180$도 이다.

 

$$\begin{array}{}\text{(8)}& i={\cos }^{-1}\left(\frac{{\hat{i}}_3\cdot \overrightarrow{h}}{h}\right)\end{array}$$

 

RAAN $\mathrm{\Omega }$:
RAAN은 ${\hat{i}}_1$ 과 $\overrightarrow{n}$ 의 사잇각이므로 ${\hat{i}}_1$ 과 $\overrightarrow{n}$ 의 내적을 이용하여 계산할 수 있다.

 

$$\begin{array}{}\text{(9)}& \mathrm{\Omega }={\cos }^{-1}\left(\frac{{\hat{i}}_1\cdot \overrightarrow{n}}{n}\right)\end{array}$$

 

 

 

RAAN은 $-180$도에서 $180$도의 범위를 가지나 식 (9)의 계산으로는 $0$도에서 $180$도의 범위밖에 계산하지 못하므로 부호를 구별할 수 있는 정보가 더 필요하다. 다음 그림의 A영역에서는 RAAN 이 $0$도에서 $180$도 사이지만 B영역에서는 RAAN 이 $-180$도에서 $0$도 사이가 된다. 만약 $\overrightarrow{n}$ 이 B영역에 있다면 식 (9)로 계산한 RAAN 을 음수로 바꾸어 줘야 한다. 두 영역의 구별은 ${\hat{i}}_2$ 와 $\overrightarrow{n}$ 의 사잇각이 예각인지 둔각인지를 판별함으로써 가능하므로, RAAN 부호에 관한 다음 교정식을 도입한다.

 

$$\text{만약 }{\hat{i}}_2\cdot \overrightarrow{n}<0\text{ 이면, }\mathrm{\Omega }\leftarrow -\mathrm{\Omega }$$

 

 

 

적도궤도인 경우 $\overrightarrow{n}$ 이 정의되지 않으므로 RAAN 역시 정의되지 않는다. 이 경우 RAAN을 $0$도로 간주한다.

근점편각 $\omega$:
근점편각은 $\overrightarrow{n}$ 과 $\overrightarrow{e}$ 의 사잇각이므로 $\overrightarrow{n}$ 과 $\overrightarrow{e}$ 의 내적을 이용하여 구할 수 있다.

 

$$\begin{array}{}\text{(10)}& \omega ={\cos }^{-1}\left(\frac{\overrightarrow{n}\cdot \overrightarrow{e}}{ne}\right)\end{array}$$

 

 

 

근점편각은 궤도의 운동방향으로 측정해야 한다. RAAN 과 마찬가지로 식 (10)으로 계산한 근점편각도 0도에서 180도의 범위밖에 계산하지 못하므로 부호를 구별할 수 있는 정보가 더 필요하다. 아래 그림의 A영역에서는 근점편각이 $0$도에서 $180$사이지만 B영역에서는 근점편각이 $-180$도에서 $0$도 사이가 된다. 만약 $\overrightarrow{e}$ 가 B영역에 있다면 식 (10)으로 계산한 근점편각을 음수로 바꾸어 줘야 한다. 두 영역의 구별은 ${\hat{i}}_3$ 와 $\overrightarrow{e}$ 의 사잇각이 예각인지 둔각인지를 판별함으로써 가능하므로, 근점편각 부호에 관한 다음 교정식을 도입한다.

 

$$\text{만약 }{\hat{i}}_3\cdot \overrightarrow{e}<0\text{ 이면, }\omega \leftarrow -\omega$$

 

 

 

적도궤도인 경우 $\overrightarrow{n}$ 이 정의되지 않으므로 근점편각은 ${\hat{i}}_1$ 과 $\overrightarrow{e}$ 의 사잇각으로 정의하며 ${\hat{i}}_1$ 과 $\overrightarrow{e}$ 의 내적을 이용하여 구할 수 있다.

 

$$\begin{array}{}\text{(11)}& \omega ={\cos }^{-1}\left(\frac{{\hat{i}}_1\cdot \overrightarrow{e}}{e}\right)\end{array}$$

 

다음 그림의 A영역에서는 근점편각이 $0$도에서 $180$도 사이지만 B영역에서는 근점편각이 $-180$도에서 $0$도 사이가 되므로, 만약 $\overrightarrow{e}$ 가 B영역에 있다면 식 (11)로 계산한 근점편각을 음수로 바꾸어 줘야 한다. 두 영역의 구별은 ${\hat{i}}_2$ 와 $\overrightarrow{e}$ 의 사잇각이 예각인지 둔각인지를 판별함으로써 가능하다. 한편, 경사각이 $180$도이면 궤도의 회전방향이 반대이므로 근점편각의 부호를 반대로 해줘야 한다. 따라서 근점편각 부호에 관한 다음 교정식을 도입하기로 한다.

 

$$\begin{array}{rl}& \text{만약 }{\hat{i}}_2\cdot \overrightarrow{e}<0\text{ 이면, }\omega \leftarrow -\omega \\ \\ & \text{만약 }i={180}^0\text{ 이면, }\omega \leftarrow -\omega \end{array}$$

 

 

 

원궤도인 경우 $\overrightarrow{e}$ 가 정의되지 않으므로 근점편각 역시 정의되지 않는다. 이 경우 근점편각을 $0$도로 간주한다.

 

 

실제비행각 $\theta$:
실제비행각은 위치벡터 $\overrightarrow{r}$ 의 함수이므로 다른 궤도요소가 모두 상수이었던 것과는 달리 시간의 함수이다. 따라서 이를 명확히 하기 위하여 $\theta (t)$ 라고 표기하기도 한다. 이 경우는 위치벡터가 $t=t_0$ 에서 주어졌으므로 실제비행각은 $\theta (t_0)$ 라고 표기한다.

실제비행각은 $\overrightarrow{e}$ 와 $\overrightarrow{r}$ 의 사잇각이므로 $\overrightarrow{e}$ 와 $\overrightarrow{r}$ 의 내적을 이용하여 구할 수 있다.

 

$$\begin{array}{}\text{(12)}& \theta (t_0)={\cos }^{-1}\left(\frac{\overrightarrow{e}\cdot \overrightarrow{r}}{er}\right)\end{array}$$

 

실제비행각은 궤도의 운동방향으로 측정해야 한다. 다음 그림의 A영역에서는 실제비행각이 $0$도에서 $180$도 사이지만 B영역에서는 실제비행각이 $-180$도에서 $0$도 사이가 되므로, 만약 $\overrightarrow{r}$ 이 B영역에 있다면 식 (12)로 계산한 실제비행각을 음수로 바꾸어 줘야 한다. 두 영역의 구별은 $\overrightarrow{r}$ 과 $\overrightarrow{v}$ 의 사잇각이 예각인지 둔각인지를 판별함으로써 가능하므로, 실제비행각 부호에 관한 다음 교정식을 도입한다.

 

$$\text{만약 }\overrightarrow{r}\cdot \overrightarrow{v}<0\text{ 이면, }\theta (t_0)\leftarrow -\theta (t_0)$$

 

 

 

원궤도인 경우 $\overrightarrow{e}$ 가 정의되지 않으므로 실제비행각은 $\overrightarrow{n}$ 과 $\overrightarrow{r}$ 의 사잇각으로 정의하며 $\overrightarrow{n}$ 과 $\overrightarrow{r}$ 의 내적을 이용하여 구할 수 있다.

 

$$\begin{array}{}\text{(13)}& \theta (t_0)={\cos }^{-1}\left(\frac{\overrightarrow{n}\cdot \overrightarrow{r}}{nr}\right)\end{array}$$

 

다음 그림의 A영역에서는 실제비행각이 $0$도에서 $180$도 사이지만 B영역에서는 실제비행각이 $-180$도에서 $0$도 사이가 되므로, 만약 $\overrightarrow{r}$ 이 B영역에 있다면 식 (13)으로 계산한 실제비행각을 음수로 바꾸어 줘야 한다. 두 영역의 구별은 $\overrightarrow{n}$ 과 $\overrightarrow{v}$ 의 사잇각이 예각인지 둔각인지를 판별함으로써 가능하므로, 실제비행각 부호에 관한 다음 교정식을 도입한다.

 

$$\text{만약 }\overrightarrow{n}\cdot \overrightarrow{v}>0\text{ 이면, }\theta (t_0)\leftarrow -\theta (t_0)$$

 

 

 

원궤도이면서 적도궤도인 경우 $\overrightarrow{e}$ 뿐만 아니라 $\overrightarrow{n}$ 도 정의되지 않으므로 실제비행각은 ${\hat{i}}_1$ 과 $\overrightarrow{r}$ 의 사잇각으로 정의하며 ${\hat{i}}_1$ 과 $\overrightarrow{r}$ 의 내적을 이용하여 구할 수 있다.

 

$$\begin{array}{}\text{(14)}& \theta (t_0)={\cos }^{-1}\left(\frac{{\hat{i}}_1\cdot \overrightarrow{r}}{r}\right)\end{array}$$

 

다음 그림의 A영역에서는 실제 비행각이 $0$도에서 $180$도 사이지만 B영역에서는 실제비행각이 $-180$도에서 $0$도 사이가 되므로, 만약 $\overrightarrow{r}$ 이 B영역에 있다면 식 (14)로 계산한 실제비행각을 음수로 바꾸어 줘야 한다. 두 영역의 구별은 ${\hat{i}}_1$ 과 $\overrightarrow{v}$ 의 사잇각이 예각인지 둔각인지를 판별함으로써 가능하므로, 실제 비행각 부호에 관한 다음 교정식을 도입한다.

 

$$\text{만약 }{\hat{i}}_1\cdot \overrightarrow{v}>0\text{ 이면, }\theta (t_0)\leftarrow -\theta (t_0)$$

 

 

 

예를 들어서, 어떤 싯점에서 우주비행체의 위치벡터와 속도벡터가 다음과 같이 주어졌다면,

 

$$\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}$$

 

궤도요소 $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}$$

 

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

 

function [a,es,ii,W,w,theta] = rv2coe(r,v)

% [a,e,i,W,w,theta] = rv2coe(r,v)
%        Converting vectors R and V in inertial frame to orbit elements,
%        where V=dR/dt in inertial frame.
% Input  (R, V) in km, km/sec
% Output (a, e, i, W, w, theta) in (km, sec) and degrees  (p=a(1-e^2))
%        equatorial orbit: W=0, w is angle between (i_1) and perigee
%        circular orbit:   w=0, theta is angle between ascending node and R
%        equatorial and circular orbit: W=w=0, theta is angle between (i_1) and R
%
% coded by st.watermelon

Req = 6378.137; % km Earth radius
e_er = 0.0818191908426; % Earth eccentricity


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

% angular momentum, ascending node and eccentricity vectors
h = cross(r,v);
n = cross([0 0 1]',h);
ev = ((v'*v-mu/norm(r))*r-(r'*v)*v)/mu;

% 
p = h'*h/mu;   % semi-latus rectum
es = norm(ev); % eccentricity
ns = norm(n);
ii = acos([0 0 1]*h/(norm(h)))*180/pi;  % inclination
a = p/(1-es^2);

eps = 1e-8;
 

if abs(ii)<=eps | abs(ii-180)<=eps % equatorial orbit

    W=0;
    
    if abs(es)>eps % equatorial, but NOT circular orbit
    
        w = acos([1 0 0]*ev/es)*180/pi;
        if ([0 1 0]*ev<0) w = -w; end
        if (abs((ii-180))<=eps) w = -w; end
        theta = acos(ev'*r/(es*norm(r)))*180/pi;
        if (r'*v<0) theta = -theta; end
        
    else % euqatorial and circular orbit
        w = 0;
        theta = acos([1 0 0]*r/(nomr(r)))*180/pi;
        if ([1 0 0]*v>0) theta = -theta; end
    end
    
else  % NOT equatorial orbit
    
    W = acos([1 0 0]*n/ns)*180/pi;
    if ([0 1 0]*n<0) W = -W; end
    
    if abs(es)>eps % NOT circular orbit
    
        w = acos(n'*ev/(ns*es))*180/pi;
        if ([0 0 1]*ev<0) w = -w; end
        theta=acos(ev'*r/(es*norm(r)))*180/pi;
        if (r'*v<0) theta = -theta; end
        
    else % circular orbit
        w = 0;
        theta = acos(n'*r/(ns*norm(r)))*180/pi;
        if (n'*v>0) theta = -theta; end
    end
    
end