궤도 계산 2D - 2 : 태양-지구 궤도

하루를 24시간이라 하자.

1 day = 24 hours = 1,440 mins = 86,400 mins

공전면에 대한 자전축 기울기 $\epsilon_{earth} = 23.44^{\circ}$

항성일 Sidereal day 86164.0905sec, 23.9344696 hrs

Orbital Parameters

근일점 Perihelion $r_p = 147.095 \cdot10^6 \text{km}$

원일점 Aphelion $r_a = 152.100 \cdot10^6\text{km}$

최고 궤도 속력 $v_p = 30.29 \text{km}$

최저 궤도 속력 $v_a = 29.29 \text{km}$

공전 궤도 이심률 Orbit eccentricity $e_{earth} = 0.0167282$

$$e=\frac{r_a - r_p}{r_a + r_p}$$

장반경 Semi-major axis $a = (r_p + r_a)/2 = 149.597.5 \cdot 10^6\text{km}$

단반경 Semi-minor axis $b = a \sqrt{1-e^2} = 149.576.6 \cdot 10^6\text{km}$

태양의 GM $GM_s = 132,712.0 \cdot 10^6 \text{km3/sec2}$

지구의 GM $GM_e = 0.39860 \cdot 10^6 \text{km3/sec2}$

$mu = GM_s + GM_e =  132.712.4 \cdot 10^6 \text{km3/sec2}$

비 각 운동량 Specific Angular Momentum $h= 4,455.5100 \cdot 10^6\text{km2/s}$

$$h=r_p \cdot v_p = r_a \cdot v_a \hspace{5mm} \text{or} \hspace{5mm} r = \frac{1}{\mu} \frac{h^2}{1+e \text{cos}\theta} \rightarrow h = \sqrt{r_p \cdot \mu \left( 1+ e \right) }$$

항성 공전 주기 Sidereal orbit period $T_e = 365.256 \text{days} = 8766.134 \text{hrs}$

$$T=\frac{2 \pi}{\sqrt{\mu}} a^{3/2} = \frac{2 \pi}{\mu^2}\left(\frac{h}{\sqrt{1-e^2}}\right)^3$$

항성 자전 주기(항성일) Sidereal rotation period, Revolution Period $t_e = 23.9345 \text{hrs}$

태양일 Solar day $1 \text{day}=\frac{1}{1/t_e - 1/T_e}=24.0000 \text{hrs}$

태양일은 회합 주기를 구하는 방법으로 구할 수 있다.

 

공전 궤도 경사 Orbit inclination $i_{moon} = 0.00^{\circ}$

 

Orbital Motion (1)

태양과 지구에 대해, 공전 궤도를 보면 다음과 같다.

항성 자전 주기 $t_e$와 태양일 $\text{SD}$에 대해서 애니메이션을 그려보면 다음과 같다.

지구의 항성 자전 주기 기준으로 볼 때, 지구는 항상 일정한 위상을 보인다.

항성 자전 주기 기준으로 태양에 대한 지구의 위치와 진행 방향, 그리고 지구의 방위각

태양일은 지구의 회합일이므로 지구의 위상은 항상 태양을 바라보는 형태로 나타나게 된다.

태양일 기준으로 태양에 대한 지구의 위치와 진행 방향, 그리고 지구의 방위각

애니메이션 코드는 다음과 같다.

%% 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;
Twe_body = Twe_sun;
Twe_sat  = Twe_earth;
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 = Twe_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");

 

끝.

 

 

[1] https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html