궤도요소 (COE) 계산

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

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

 

 

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

각운동량 벡터 \(\vec{h}\):
각운동량 벡터 \(\vec{h}\) 는 정의에 의해서

 

\[ \vec{h}= \vec{r}(t_0 ) \times \vec{v}(t_0 ) \tag{1} \]

 

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

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

 

\[ \vec{n}= \hat{i}_3 \times \vec{h} \tag{2} \]

 

\(\vec{h}\) 가 상수벡터이므로 승교선 벡터도 상수벡터이다.

 

 

이심율 벡터 \(\vec{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)에 의하면 다음 식이 있었다.

 

\[ \vec{v} \times \vec{h}= \mu \frac{\vec{r}}{r}+\vec{B} \tag{3} \]

 

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

 

\[ \vec{e}= \frac{\vec{B}}{\mu} \tag{4} \]

 

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

 

\[ \vec{e}= \frac{1}{\mu} \left( \vec{v} \times \vec{h}- \mu \frac{\vec{r}}{r} \right) \tag{5} \]

 

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

 

\[ \vec{e}= \frac{1}{\mu} \left( \left(v^2-\frac{\mu}{r} \right) \vec{r}-( \vec{r} \cdot \vec{v} ) \vec{v} \right) \tag{6} \]

 

이 된다.

 

 

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

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

 

\[ \begin{align} & e=|\vec{e} | \tag{7} \\ \\ & p= \frac{h^2}{\mu}, \ \ \ \ \ a=\frac{p}{1-e^2 } \end{align} \]

 

여기서 \(h=|\vec{h}|\) 이다.

 

 

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

 

\[ \ i= \cos^{-1} \left( \frac{\hat{i}_3 \cdot \vec{h} }{h} \right) \tag{8} \]

 

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

 

\[ \ \Omega = \cos^{-1} \left( \frac{ \hat{i}_1 \cdot \vec{n} }{n} \right) \tag{9} \]

 

 

 

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

 

\[ \ \mbox{만약 } \ \hat{i}_2 \cdot \vec{n} \lt 0 \ \mbox{ 이면, } \ \Omega \gets -\Omega \]

 

 

 

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

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

 

\[ \ \omega = \cos^{-1} \left( \frac{ \vec{n} \cdot \vec{e} }{ne} \right) \tag{10} \]

 

 

 

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

 

\[ \ \mbox{만약 } \ \hat{i}_3 \cdot \vec{e} \lt 0 \ \mbox{ 이면, } \ \omega \gets -\omega \]

 

 

 

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

 

\[ \ \omega = \cos^{-1} \left( \frac{ \hat{i}_1 \cdot \vec{e} }{e} \right) \tag{11} \]

 

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

 

\[ \begin{align} & \mbox{만약 } \ \hat{i}_2 \cdot \vec{e} \lt 0 \ \mbox{ 이면, } \ \omega \gets -\omega \\ \\ & \mbox{만약 } \ i=180^0 \ \mbox{ 이면, } \ \omega \gets -\omega \end{align}\]

 

 

 

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

 

 

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

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

 

\[ \ \theta (t_0) = \cos^{-1} \left( \frac{ \vec{e} \cdot \vec{r} }{er} \right) \tag{12} \]

 

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

 

\[ \mbox{만약 } \ \vec{r} \cdot \vec{v} \lt 0 \ \mbox{ 이면, } \ \theta (t_0) \gets - \theta (t_0) \]

 

 

 

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

 

\[ \ \theta (t_0) = \cos^{-1} \left( \frac{ \vec{n} \cdot \vec{r} }{nr} \right) \tag{13} \]

 

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

 

\[ \mbox{만약 } \ \vec{n} \cdot \vec{v} \gt 0 \ \mbox{ 이면, } \ \theta (t_0) \gets - \theta (t_0) \]

 

 

 

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

 

\[ \ \theta (t_0) = \cos^{-1} \left( \frac{ \hat{i}_1 \cdot \vec{r}}{r} \right) \tag{14} \]

 

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

 

\[ \mbox{만약 } \ \hat{i}_1 \cdot \vec{v} \gt 0 \ \mbox{ 이면, } \ \theta (t_0) \gets - \theta (t_0) \]

 

 

 

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

 

\[ \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} \]

 

궤도요소 \(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} \]

 

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

 

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