Orbital Mechanics, Ch.8.10 행성의 천체력 계산

 상태 벡터 $R$과 $V$는 태양 중심 황도 좌표계 (Heliocentric ecliptic frame)에 상대적으로 정의된다. 이는 지구중심 적도 좌표계 (Geo-centric equatorial frame)과 매우 유사하다. 인력 중심(The center of attraction)으로써 지구를 태양으로 대체하고, 황도 평면이 지구 적도 평면을 대체하면 된다. 춘분선은 관성 X 축으로 동일하게 정의된다.

 실제적인 항성간 임무를 설계하기 위해서는, 어떤 주어진 시간에서 행성의 상태 벡터를 결정할수 있어야만 한다. 그래서 테이블 [1]과 같은 궤도 성분이 제공된다.

 

 

주어지는 일자와 시간에서 행성의 상태 벡터를 결정하는 알고리즘은 다음과 같다.

 

1) 율리우스일 JD 계산하기

$$\text{JD} = J_0 + \frac{\text{UT}}{24}$$

$$J_0 = 367y - \text{INT} \left[\frac{7}{4}\left[ y+ \text{INT} \left( \frac{m+9}{12}\right) \right] \right] + \text{INT} \left( \frac{275m}{9} \right) + d + 1,721,013.5$$

 

2) J2000 시대로부터의 율리우스 일 $T_0$ 계산하기

$$T_0 = \frac{J_0 - 2,451,545}{36,525}$$

 

3) 성분 변화량 계산하기

$$x = x_0 + \dot{x} T_0$$

여기서 다음과 같은 어떤 시간에서의 행성의 궤도 성분 (Planetary orbital element)을 구하게 된다.

• $a$ 장축단, semi-major axis
• $e$ 이심률, eccentricity
• $i$ 궤도 기울기, the inclination to the ecliptic plane
• $\Omega$ 승교점의 적경, the right ascension of the ascending node
• $\tilde{\omega}$ 근일점의 경도, the longitude of perihelion
• $L$ 평균 경도, the mean longitude

 

4) 장축 반지름 semi-major axis $a$와 이심률 $e$으로부터 비 각 운동량 specific angular momentum $h$ 구하기

$$h = \sqrt{\mu a \left( 1-e ^2 \right) }$$

 

5) 근일점의 매개변수 argument of perihelion $\omega$와 평균 이각 $M$ 구하기

$$\omega = \tilde{\omega} - \Omega$$

$$M = L - \tilde{\omega}$$

 

6) 편심 이각 $E$와 전근점 이각 $\theta$ 구하기

평균 근점 이각 $M_e$와 편심 이각 $E$의 관계는 다음과 같다. 이를 케플러 방정식이라 부른다.

$$M_e = E-e\text{sin} E$$

편심 이각 $E$와 진 근점 이각 $\theta$간의 관계는 다음과 같다.

$$\frac{E}{2} = \sqrt{\frac{1-e}{1+e}} \text{tan}\frac{\theta}{2}$$

$$\theta = 2 \text{tan}^{-1}\left( \sqrt{\frac{1+e}{1-e}} \text{tan} \frac{E}{2} \right)$$

 

7) 다른 상태량을 통해서 태양 중심의 위치 벡터와 속도 벡터를 얻는다.

초점으로부터의 거리 $r$
$$r= \frac{h^2}{\mu} \frac{1}{1+e \text{cos}\theta }$$
시선 속도 Radial velocity $v_r$
$$v_r = \frac{\mu}{h} e \text{sin}\theta$$
횡속도 Transverse velocity $v_{\perp}$
$$v_r = \frac{\mu}{h} \left(1+ e \text{cos}\theta \right)$$

 

8) 태양좌표계 혹은 지구 좌표계에 고정하기

$$r^x = \frac{h^2}{\mu} \frac{1}{1+e \text{cos}\theta} \left[ \matrix{\text{cos}\theta \\ \text{sin}\theta \\ 0 } \right]$$

$$v^x = \frac{\mu}{h} \left[ \matrix{-\text{sin}\theta \\ e+\text{cos}\theta \\ 0} \right]$$

 

위에서 계산한 위치 벡터 $r^x$와 속도 벡터 $v^x$는 궤도 평면 상에서의 위치이다.

기준이 되는 고정 좌표계로 변환은 다음과 같다.

$$Q_x^X = \left[ \matrix{ \text{cos}\Omega \text{cos}\omega - \text{sin}\Omega \text{sin}\omega \text{cos}i & -\text{cos}\Omega \text{sin}\omega -\text{sin}\Omega \text{cos}i \text{cos}\omega  & \text{sin}\Omega \text{sin}i \\ \text{sin}\Omega \text{cos}\omega +\text{cos}\Omega \text{sin}\omega \text{cos}i & -\text{sin}\Omega \text{sin}\omega +\text{cos}\Omega \text{cos}i \text{cos}\omega  & -\text{cos}\Omega \text{sin}i \\ \text{sin}\omega \text{sin}i & \text{sin}i \text{cos}\omega  & \text{cos}i } \right]$$

 

태양중심, 혹은 지구 중심 좌표계로의 변환은 다음과 같다.

$$r^X = Q_x^X r^x$$

$$v^X = Q_x^X v^x$$

 

전체 코드는 다음과 같다,

a           =          planet.x(1);
e           =          planet.x(2);
i           = warp_360(planet.x(3));
Omega       = warp_360(planet.x(4));
tilde_omega = warp_360(planet.x(5));
L           = warp_360(planet.x(6));
mu          = planet.mu + sun_mu;
h           = sqrt(mu*a*(1-e.*e));
 
Sec2Day = 1/(24*60*60);
Day2Year= 1/(365.256);
D2R     = pi/180;
R2D     = 180/pi;
 
omega   = warp_360(tilde_omega - Omega);
M       = warp_360(L - tilde_omega);
E       = warp_360(R2D*double(vpasolve((Eval - e * sin(Eval)) - M * D2R, Eval)));
theta   = warp_360(2 * atan(sqrt((1-e)/(1+e))*tand(E/2)) * R2D);
 
rx  = (h.*h/mu)/(1+e*cosd(theta))*[cosd(theta), sind(theta), 0];
vx  = (h/mu) * [-sind(theta), e + cosd(theta), 0];
QxX = [ cosd(Omega)*cosd(omega)-sind(Omega)*sind(omega)*cosd(i), -cosd(Omega)*sind(omega)-sind(Omega)*cosd(i)*cosd(omega),  sind(Omega)*sind(i);
        sind(Omega)*cosd(omega)+cosd(Omega)*sind(omega)*cosd(i), -sind(Omega)*sind(omega)+cosd(Omega)*cosd(i)*cosd(omega), -cosd(Omega)*sind(i);
                                            sind(omega)*sind(i),                                      sind(i)*cosd(omega),              cosd(i)];

 

예제 1) 2003년 8월 27일의 12hr UT의 지구와 화성의 거리를 계산해보자.

 

예제 2) 오늘 밤 자정, 경기 평택 어딘가에서의 달의 모습을 나타내어 보자.

오늘 밤 자정은 2023년 03월 13일 0hr KST

경기 평택 어딘가는 N 37.00000, E 127.000000 이라 하자.

 

[1]  "Planetary Orbital Elements - Planetary Mean Orbits (J2000)" http://www.met.rdg.ac.uk/~ross/Astronomy/Planets.html