태양 주위의 가관측한 행성의 궤도

행성의 궤도를 그리는 스크립트를 작성해보았다.

지구의 회전 위상만 찾으면 특정 시간에서의 가관측한 행성을 알 수 있을 것이다.

 

 

전체 소스

행성의 율리우스력에 대한 상태를 계산

행성의 항성 주기에 맞추어서 위치를 계산

%% Get Julian day (JD)
clc;clear;
close all;
addpath('libs');
time = clock();
SD = 24*60*60; % Solar day
Day2Year = 365.256;
 
% time(1) = 2003;
% time(2) = 8;
% time(3) = 27;
% time(4) = 12;
% time(5) = 0;
% time(6) = 0;
 
fprintf("%4d.%02d.%02d %02d h %02d m %06.3f s\n",time);
UT = time(4) + time(5)/60 + time(6)/3600;
J0 = (367*time(1) - floor((7/4)*(time(1)+floor((time(2)+9)/12))) + floor(275*time(2)/9) + time(3) + 1721013.5);
JD = J0 + UT/24;
T0 = (JD - 2451545)/36525;
 
 
%% Planetary Properties
clc;
sun.mu = 132712E6;
% planet_list = 
Mercury.name = "Mercury";
Mercury.x0   = [0.38709893 0.20563069 7.00487 48.33167 77.45645 252.25084];
Mercury.dx   = [0.00000066 0.00002527 -23.51 -446.30 573.57 538101628.29];
Mercury.x    = Mercury.x0;
Mercury.mu   = 22030;
Mercury.osp  = 87.97; % Orbit Sidereal Period
Mercury      = trans(Mercury);
 
Venus.name = "Venus";
Venus.x0   = [0.72333199 0.00677323 3.39471 76.68069 131.53298 181.97973];
Venus.dx   = [0.00000092 -0.00004938 -2.86 -996.89 -108.80 210664136.06];
Venus.x    = Venus.x0;
Venus.mu   = 324900;
Venus.osp  = 224.7;
Venus      = trans(Venus);
 
Earth.name = "Earth";
Earth.x0 = [1.00000011 0.01671022 0.00005 -11.26064 102.94719 100.46435];
Earth.dx = [-0.00000005 -0.00003804 -46.94 -18228.25 1198.28 129597740.63];
Earth.x  = Earth.x0;
Earth.mu = 0.39860E6;
Earth.osp= 365.256;
Earth    = trans(Earth);
 
Mars.name = "Mars";
Mars.x0 = [1.52366231 0.09341233 1.85061 49.57854 336.04084 355.45332];
Mars.dx = [-0.00007221 0.00011902 -25.47 -1020.19 1560.78 68905103.78];
Mars.x  = Mars.x0;
Mars.mu = 0.042828E6;
Mars.osp= 1.881*Day2Year;
Mars = trans(Mars);
 
Jupiter.name = "Jupiter";
Jupiter.x0 = [ 5.20336301 0.04839266 1.30530 100.55615 14.75385 34.40438];
Jupiter.dx = [0.00060737 -0.00012880 -4.15 1217.17 839.93 10925078.35];
Jupiter.x  = Jupiter.x0;
Jupiter.mu = 126.687E6;
Jupiter.osp= 11.86*Day2Year;
Jupiter = trans(Jupiter);
 
Saturn.name = "Saturn";
Saturn.x0 = [ 9.53707032 0.05415060 2.48446 113.71504 92.43194 49.94432];
Saturn.dx = [-0.00301530 -0.00036762 6.11 -1591.05 -1948.89 4401052.95];
Saturn.x  = Saturn.x0;
Saturn.mu = 37.931E6;
Saturn.osp= 29.46*Day2Year;
Saturn = trans(Saturn);
 
RGBs = [[0.8500 0.3250 0.0980];
        [0.9290 0.6940 0.1250];
        [0.0000 0.4470 0.7410];
        [0.4940 0.1840 0.5560];
        [0.4660 0.6740 0.1880];
        [0.3010 0.7450 0.9330];
        [0.6350 0.0780 0.1840]];
 
%%
clc;
 
Planets = ["Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn"];
 
clear time;
time(1) = 2020;
time(2) = 1;
time(3) = 1;
time(4) = 0;
time(5) = 0;
time(6) = 0;
tictocs = datenum(time)+1800;
for i = 1 : 200
% tictocs = tictocs + 365.256/12/3;
tictocs = tictocs + 3;
time = datevec(tictocs);
% fprintf("%4d.%02d.%02d %02d h %02d m %06.3f s\n",time);
UT = time(4) + time(5)/60 + time(6)/3600;
J0 = (367*time(1) - floor((7/4)*(time(1)+floor((time(2)+9)/12))) + floor(275*time(2)/9) + time(3) + 1721013.5);
JD = J0 + UT/24;
T0 = (JD - 2451545)/36525;
 
% Get States of Planet @ specific time
for i = 1 : length(Planets)
    strings = sprintf("%s.x = %s.x0 + %s.dx * T0;",Planets{i},Planets{i},Planets{i});
    eval(strings);
    strings = sprintf("%s = planetary_properties(sun.mu, %s);",Planets{i},Planets{i});
    eval(strings);
end
 
figure(10);
set(gcf, "Position", [50 0 600 500]);
plot(0,0,'r*'); hold on; grid on;
for i = 1 : length(Planets)
    strings = sprintf("plot(%s.r(1)/1E6, %s.r(2)/1E6, 'o', 'Color', RGBs(%d,:));", Planets{i}, Planets{i}, i);
    eval(strings);
end
for i = 1 : length(Planets)
    strings = sprintf("plot(%s.circle(1,:)/1E6, %s.circle(2,:)/1E6, 'Color', RGBs(%d,:));",Planets{i}, Planets{i},i);
    eval(strings);
end
% fprintf("%4d.%02d.%02d %02d h %02d m %06.3f s\n",time);
axis([-1 1 -1 1]*1500);
strings = sprintf("%4d.%02d.%02d %02dh %02dm %06.3fs",time);
title(strings);
xlabel('X 10^6 [km]'); ylabel('Y 10^6 [km]');
hold off;
legend(["Sun", Planets]);
 
drawnow;
strings2 = sprintf("imgs\\%s.png",strings);
saveas(gcf, strings2);
end

 

행성의 위치와 관련 상태 계산 함수

function planet_out = planetary_properties(sun_mu, planet)
syms Eval;
 
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);
 
rp = (h*h/mu)/(1+e); vp = h/rp;
ra = (h*h/mu)/(1-e); va = h/ra;
%a  = (rp+ra)/2;
T  = (2*pi*power(a,3/2))/sqrt(mu);
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)];
try
    tmp = planet.circle(1,1);
catch
    degs        = linspace(0,360,60);
    const = (h.*h/mu)./(1+e*cosd(degs))';
    circle  = [cosd(degs)'.*const, sind(degs)'.*const, zeros(size(degs))'];
    planet_out.circle= QxX * circle';
end
 
% fprintf("%8s\t S-major axis %12.3f E6 km , h %12.3f E6 km2/s, e %9.8f \n",planet.name, a/1E6, h/1E6, e);
% fprintf("\t\t\t Perigee      %12.3f E6 km , %12.6f km/s\n", rp/1E6, vp);
% fprintf("\t\t\t Apogee       %12.3f E6 km , %12.6f km/s\n", ra/1E6, va);
% fprintf("\t\t\t Sidereal Orbital Period %10.5f years, %15.8f days - Calculated\n",T*Sec2Day*Day2Year, T*Sec2Day);
 
planet_out = planet;
planet_out.rp = rp;
planet_out.ra = ra;
planet_out.vp = vp;
planet_out.va = va;
planet_out.a  = a ;
planet_out.T  = T ;
 
planet_out.rx   = rx;
planet_out.vx   = vx;
planet_out.QxX  = QxX;
planet_out.r    = QxX * rx';
planet_out.v    = QxX * vx';
 
 
end

 

0~360도로 한정하는 함수

function y = warp_360(x)
mult = 360.0;
y  = mod(x, mult);
end

값 변환 함수

function out = trans(in)
AU = 1.49597871E8;
out = in;
out.x0(1) = out.x0(1) * AU;
out.dx(1) = out.dx(1) * AU;
out.x(1)  = out.x(1)  * AU;
 
out.dx(3) = out.dx(3) / 3600;
out.dx(4) = out.dx(4) / 3600;
out.dx(5) = out.dx(5) / 3600;
out.dx(6) = out.dx(6) / 3600;
 
end