궤도 계산 2D - 1 : 지구-달

식 정리

타원 궤도를 가정하면..

강체의 단위 질량 당 상대 각운동량(relative angular momentum of the body per unit mass) $h=\sqrt{\mu a \left (1-e^2 \right)}$ (m2/s)

극좌표계의 반경 $r=\frac{h^2}{\mu} \frac{1}{1+e \text{cos}\theta}$ (m)

근지점 반경 Perigee radius $r_p = \frac{h^2}{\mu}\frac{1}{1+e}$ (m)

원지점 반경 Apogee radius $r_a = \frac{h^2}{\mu}\frac{1}{1-e}$ (m)

장반경 semi-major axis $a=\frac{h^2}{\mu} \frac{1}{1-e^2} = \left(r_p + r_a \right)/2$ (m)

단반경 semi-minor axis $b=a\sqrt{1-e^2}$ (m)

이심률 eccentricity $e=\left(r_a - r_p \right) / \left( r_a - r_p \right) $

공전 주기 $T=\frac{2\pi}{\sqrt{\mu}} a^{3/2}$

공전 주기 $T$ 와 시간 $t$ 에 대한 평균 이각 $M_e$ (degree)

$$M_e = \frac{2\pi}{T} t $$

평균 이각 $M_e$으로부터 편심 이각 $E$ 

$$E=\text{arg}_E \left( M_e =  E-e \text{sin}(E) \right)$$

편심 이각 $E$으로부터 전근점 이각 $\theta$

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

달의 궤도 특성

Orbit of the Moon [2]

지구의 공전 궤도(황도)에 대한 지구의 자전축 기울기는 $\epsilon_{earth} = 23.4^{ \circ}$

지구의 공전 궤도(황도)에 대한 달의 공전 궤도 경사 $i_{moon}=5.145^{\circ}$

달의 공전 궤도에 대한 달의 자전축 기울기는 $\epsilon_{moon}=6.687^{\circ}$

지구의 공전 궤도(황도)에 대한 달의 자전축 기울기는 $1.543^{\circ}=\epsilon_{moon} - i_{moon}$

Perigee $r_p=363,229 \text{km}$

Apogee $r_a=405,400 \text{km}$

비각운동량 Angular momentum $h = \sqrt{r_p \mu + \left( 1+e \text{cos}0\right) }=0.39424 \text{E6 m2/s} $

Sidereal Rotation Period of the Earth	Sidereal Orbital Period of Moon about the Earth	Sidereal Rotation Period of Moon
23.9345 hrs	654.833 hrs (계산값)	655.72 hrs
0.997271 days	27.2847 days (계산값)	27.321667 days

 

애니메이션으로 만든 지구-달 궤도의 시간에 대한 위상과 위치

아래 애니메이션은 1시간 간격으로 시뮬레이션을 했을 때, 지구의 위상과 달의 위상과 위치를 나타낸 것이다.

달은 지구에 대한 동기 궤도를 가진다고 한다. 이는 달의 공전 속력과 자전 속력이 동일하기에 지구의 입장에서 볼 때, 항상 같은 면만 보이는 현상이다. 이는 위의 표로 보아도 알가 쉽다.

 

코드로 구현한 애니메이션

코드는 다음과 같다.

%% Basic Parameters
clc;clear;
 
SD = 24*60*60; % Solar day
Twe_sun   = 609.12  * 60*60;
Twe_earth = 23.9345 * 60*60; % Sidereal rotation period 
Twe_moon  = 655.720 * 60*60; % Sidereal rotation period
 
% Earth - Moon
fprintf("Earth - Moon Orbit\n");
relation = "Earth-Moon";
ra = 0.4055E6;
rp = 0.3633E6;
Twe_body = Twe_earth;
Twe_sat  = Twe_moon;
i = 5.14;
GM_m = 0.004900E6;
GM_e = 0.398600E6;
mu = GM_m + GM_e;
 
% Sun - Earth
% fprintf("Sun - Earth Orbit\n");
% relation = "Sun-Earth";
% rp = 147.095E6;
% ra = 152.100E6;
% GM_s = 132712E6;
% GM_e = 0.39860E6;
% te = 23.9345 * 60*60;
% i = 0.000;
% mu = GM_s + GM_e;
 
e = (ra - rp)/(ra + rp);
h = sqrt(rp * mu * (1 + e * cos(0)));
T = (2 * pi / power(mu,2)) * power(h/sqrt(1-e*e),3);
 
fprintf("GM        \t:%15.6f 10^6 km2/sec \n",mu/1E6);
fprintf("Semi-major\t:%15.6f 10^6 km \n",ra/1E6);
fprintf("Semi-minor\t:%15.6f 10^6 km \n",rp/1E6);
fprintf("Periods\n");
fprintf("- Sidereal Rotation ------ Sidereal Orbit ------ Sidereal Rotation ---\n"); 
fprintf("%15.2f sec - %15.2f sec - %15.2f sec \n", Twe_body, T, Twe_sat);
fprintf("%15.6f hrs - %15.6f hrs - %15.6f hrs \n", Twe_body/3600, T/3600, Twe_sat/3600);
fprintf("%15.6f day - %15.6f day - %15.6f day \n", Twe_body/(24*60*60), T/(24*60*60), Twe_sat/(24*60*60));
 
%% Calculate states
syms Eval;
 
we_sat  = 2*pi / Twe_sat;
we_body = 2*pi / Twe_body;
% dt = we_sat;
% dt = SD;
dt = 60*60;
scale_bar = rp * 0.01;
for i = 0 : dt : floor(T)
    t = i; % Seconds
    
    % Rotation Angle from Revolution Period
    Ae_body = we_body * t;
    Ae_sat  = we_sat  * t;
    
    % Anomaly as a function of time
    Me = warp_2pi(2 * pi * t / T);
    E = double(vpasolve((Eval - e * sin(Eval)) - Me, Eval));
    theta = 2 * atan(sqrt((1+e)/(1-e))*tan(E/2));
%     fprintf("T %8.2f hrs : E %6.1f D / t %6.1f D \n",t/(3600), E*57.296,theta*57.296);
    fprintf("T %8.2f days : E %6.1f D / t %6.1f D \n",t/(24*3600), E*57.296,theta*57.296);
    
    % Properties from true anomaly
    r  = (h*h/mu)/(1+e*cos(theta));
    r_vec = [ra * (cos(theta) - e),   rp * sin(theta)];
    ur_vec = r_vec/norm(r_vec);
    vr = (mu/h)*e*sin(theta);     % Radial component of velocity
    vt = (mu/h)*(1+e*cos(theta)); % Transverse component of velocity
    v  = norm([vr vt]);
    vr_vec = vr * ur_vec;
    vt_vec = vt * [-ur_vec(2) ur_vec(1)];
    v_vec = vr_vec + vt_vec;
    
    % Show me
    figure(1);
    angles = [0:10:360];
    % Position of Planet
    plot(ra * cos(theta),   rp * sin(theta),'xr');hold on; grid on;
    % Line from origin to planet
    plot([ra*e, ra * cos(theta)], [0 rp * sin(theta)],':r');
    % Velocity vector of Planet
    plot(ra * cos(theta)+[0, v_vec(1)] *scale_bar,...
         rp * sin(theta)+[0, v_vec(2)] *scale_bar,'b');
    plot(ra * cos(theta)+[0, vr_vec(1)]*scale_bar,...
         rp * sin(theta)+[0, vr_vec(2)]*scale_bar,'m');    % Vr
    plot(ra * cos(theta)+[0, vt_vec(1)]*scale_bar,...
         rp * sin(theta)+[0, vt_vec(2)]*scale_bar,'m');    % Vt
    % Focus of eclipse
    plot([ ra*e],[0],'or'); % 1st Focus
    plot([-ra*e],[0],'ok'); % 2nd Focus
    % Phase of Smaller planet
    plot(ra * cos(theta)+[0, cos(pi+Ae_sat)]*scale_bar*30.0,...
         rp * sin(theta)+[0, sin(pi+Ae_sat)]*scale_bar*30.0,'g');
    % Phase of Bigger planet
    plot(ra*e+[0, cos(pi+Ae_body)]*scale_bar*30.0,...
         0   +[0, sin(pi+Ae_body)]*scale_bar*30.0,'g');    % 0, 0
    
    plot(ra,0,'kx');
    plot(0,rp,'kx');
    plot(ra * cosd(angles), rp * sind(angles),'--r');   % Main Path
    plot(ra * cosd(angles), ra * sind(angles),':k');    % Mean Circle
     
    str = sprintf("Rp %6.4f E6 km", rp/1E6); text(ra*0.6, rp*1.3, str);
    str = sprintf("Ra %6.4f E6 km", ra/1E6); text(ra*0.6, rp*1.2, str);
%     str = sprintf("T  %6.3f hrs", T/(60*60)); text(ra*0.6, rp*1.1, str);
    str = sprintf("T  %6.3f days", T/(24*60*60)); text(ra*0.6, rp*1.1, str);
%     title_string = sprintf("%s : Time %8.2f hrs : E %7.2f D / t %7.2f D",relation, t/(60*60), E*57.296, theta*57.296);
    title_string = sprintf("Time %8.2f days : E %7.2f D / t %7.2f D",t/(24*60*60), E*57.296, theta*57.296);
    title(title_string);
%     legend({"Position","Course","Heading"},"Location","Best");
    hold off;
    % break;
    drawnow;
end
fprintf("End\n");

asdfasdf

 

 

 

 

[1] Moon Fact Sheet, https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html

[2] Orbit of the Moon, https://en.wikipedia.org/wiki/Orbit_of_the_Moon