clf
clear all
close all
%% global constants
global G b
G = +6.674e-20;                 % (kN*km^2/kg^2)
b = +4.185e-37;                 % (SpaceDensity/4) (kg/km^3)
%% Time
Time = 10760*24*60*60;          % (s) 167y for piriod of ne
t = 0;                          % (s) [1975-May-12 12:12:00.0000]
dt = (+3.156e+07)/256;         % (s)
%% initial positions & velocities
%% Sun (su)
r0  = [-3.869e+05;-1.516e+05;+9.562e+03];   % [x y z]*(km)
v0  = [+3.699e-03;-1.152e-02;-9.707e-05];   % [vx vy vz]*(km/s)
m   = +1.989e+30;                           % (kg)
rho = +1.408e+12;                           % (kg/km^3)
su  = [ r0 ; v0 ; m ; rho ];
sud = 2 * (3*m/(4*rho*pi))^(1/3);
%% Mercury (me)
r0  = [-5.296e+07;+1.435e+07;+6.023e+06];   % [x y z]*(km)
v0  = [-2.303e+01;-4.488e+01;-1.548e+00];   % [vx vy vz]*(km/s)
m   = +3.302e+23;                           % (kg)
rho = +5.427e+12;                           % (kg/km^3)
me = [ r0 ; v0 ; m ; rho ];
med = 2 * (3*m/(4*rho*pi))^(1/3);
%% Venus (ve)
r0  = [-1.032e+08;+3.085e+07;+6.368e+06];   % [x y z]*(km)
v0  = [-1.028e+01;-3.370e+01;+1.354e-01];   % [vx vy vz]*(km/s)
m   = +4.869e+24;                           % (kg)
rho = +5.204e+12;                           % (kg/km^3)
ve = [ r0 ; v0 ; m ; rho ];
ved = 2 * (3*m/(4*rho*pi))^(1/3);
%% Earth (ea)
r0  = [-9.446e+07;-1.184e+08;+2.493e+03];   % [x y z]*(km)
v0  = [+2.284e+01;-1.867e+01;+6.419e-05];   % [vx vy vz]*(km/s)
m   = +5.972e+24;                           % (kg)
rho = +5.515e+12;                           % (kg/km^3)
ea = [ r0 ; v0 ; m ; rho ];
ead = 2 * (3*m/(4*rho*pi))^(1/3);
%% Mars (ma)
r0  = [+1.481e+08;-1.453e+08;-6.688e+06];   % [x y z]*(km)
v0  = [+1.786e+01;+1.939e+01;-3.359e-02];   % [vx vy vz]*(km/s)
m   = +6.419e+23;                           % (kg)
rho = +3.933e+12;                           % (kg/km^3)
ma = [ r0 ; v0 ; m ; rho ];
mad = 2 * (3*m/(4*rho*pi))^(1/3);
%% Jupiter (ju)
r0  = [+7.368e+08;+7.566e+07;-1.682e+07];   % [x y z]*(km)
v0  = [-1.491e+00;+1.361e+01;-2.268e-02];   % [vx vy vz]*(km/s)
m   = +1.898e+27;                           % (kg)
rho = +1.326e+12;                           % (kg/km^3)
ju = [ r0 ; v0 ; m ; rho ];
jud = 2 * (3*m/(4*rho*pi))^(1/3);
%% Saturn (sa)
r0  = [-4.732e+08;+1.267e+09;-3.353e+06];   % [x y z]*(km)
v0  = [-9.571e+00;-3.405e+00;+4.405e-01];   % [vx vy vz]*(km/s)
m   = +5.683e+26;                           % (kg)
rho = +0.687e+12;                           % (kg/km^3)
sa = [ r0 ; v0 ; m ; rho ];
sad = 2 * (3*m/(4*rho*pi))^(1/3);
%% Uranus (ur)
r0  = [-2.368e+09;-1.426e+09;+2.544e+07];   % [x y z]*(km)
v0  = [+3.462e+00;-6.152e+00;-6.786e-02];   % [vx vy vz]*(km/s)
m   = +8.681e+25;                           % (kg)
rho = +1.318e+12;                           % (kg/km^3)
ur = [ r0 ; v0 ; m ; rho ];
urd = 2 * (3*m/(4*rho*pi))^(1/3);
%% Neptune (ne)
r0  = [-1.502e+09;-4.275e+09;+1.226e+08];   % [x y z]*(km)
v0  = [+5.094e+00;-1.771e+00;-8.092e-02];   % [vx vy vz]*(km/s)
m   = +1.024e+26;                           % (kg)
rho = +1.638e+12;                           % (kg/km^3)
ne = [ r0 ; v0 ; m ; rho ];
ned = 2 * (3*m/(4*rho*pi))^(1/3);
%%
Y = [ su ; me ; ve ; ea ; ma ; ju ; sa ; ur ; ne ];
%% marker sizes and colors
sums = 80;                      % Sun marker size
fact = 2;                       % Sun marker size (fact) times lager than Jupiter
mems = (sums/fact)*med/jud;
vems = (sums/fact)*ved/jud;
eams = (sums/fact)*(ead/jud);
mams = (sums/fact)*(mad/jud);
jums = (sums/fact);
sams = (sums/fact)*sad/jud;
urms = (sums/fact)*urd/jud;
nems = (sums/fact)*ned/jud;
%% Colors
suc = [1 1 0];
mec = [0.6 0.5 0.5];
vec = [1 1 0.7];
eac = [0.45 0.45 1];
mac = [0.9 0.5 0];
juc = [1 0.5 0.3];
sac = [0.7 0.7 0];
urc = [0.8 0.8 1];
nec = [0.6 0.6 1];
stc = [0.4 0.4 0.5];
%% Axis
axis1 = 3e8;
pov1 = [-axis1 -axis1 axis1]';
axis2 = 1.75e9;
pov2 = [-axis2 -axis2 axis2]';
axis3 = 4.75e9;
pov3 = [-axis3 -axis3 axis3]';
%% Backround Stars
nstars = 100;
BSx = ones(size(nstars,1));
BSy = ones(size(nstars,1));
BSz = ones(size(nstars,1));
f = 0;
for i=1:nstars
    BSx(i) = 2*(rand()-0.5);
    BSy(i) = 2*(rand()-0.5);
    BSz(i) = 2*(rand()-0.5);
    Wall = rand();
    test = 0;
    if Wall < 1/3
        BSx(i) = 1;
        test = test +1;
    end
    if Wall > 2/3
        BSy(i) = 1;
        test = test +1;
    end
    if test < 1
        BSz(i) = -1;
    end
end
%%
n = (Time-rem(Time,dt))/dt;
Yt = zeros(size(Y,1),n);

for i=1:n

    f1 = kepler_n_ode(t     , Y                );
    f2 = kepler_n_ode(t+dt/2, Y+f1 .*  (dt/2)  );
    f3 = kepler_n_ode(t+dt/2, Y+f2 .*  (dt/2)  );
    f4 = kepler_n_ode(t+dt  , Y+f3 .*   dt     );
    
    Y = Y + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);    
    t = t+dt;
    Yt(:,i) = Y;
    if rem(i-1,5) == 0
        subplot(3,1,1);plot3(BSx.*axis1,BSy.*axis1,BSz.*axis1,'.','color',stc);hold on
        subplot(3,1,1);plot3(Y(1,:),Y(2,:),Y(3,:),'.','markersize',sums,'markeredgecolor',suc);
        subplot(3,1,1);plot3(Y(9,:),Y(10,:),Y(11,:),'.','markersize',8*mems*(norm(pov1)/norm(Y(9:11,:)-pov1)),'markeredgecolor',mec);
        subplot(3,1,1);plot3(Y(17,:),Y(18,:),Y(19,:),'.','markersize',8*vems*(norm(pov1)/norm(Y(17:19,:)-pov1)),'markeredgecolor',vec);
        subplot(3,1,1);plot3(Y(25,:),Y(26,:),Y(27,:),'.','markersize',8*eams*(norm(pov1)/norm(Y(25:27,:)-pov1)),'markeredgecolor',eac);
        subplot(3,1,1);plot3(Y(33,:),Y(34,:),Y(35,:),'.','markersize',8*mams*(norm(pov1)/norm(Y(33:35,:)-pov1)),'markeredgecolor',mac);hold off;
        axis([-axis1 axis1 -axis1 axis1 -axis1 axis1]);
        subplot(3,1,2);plot3(BSx.*axis2,BSy.*axis2,BSz.*axis2,'.','color',stc);hold on;
        subplot(3,1,2);plot3(Y(1,:),Y(2,:),Y(3,:),'.','markersize',sums/1.5,'markeredgecolor',suc);
        subplot(3,1,2);plot3(Y(9,:),Y(10,:),Y(11,:),'.','markersize',mems*(norm(pov2)/norm(Y(9:11,:)-pov2)),'color',mec);
        subplot(3,1,2);plot3(Y(17,:),Y(18,:),Y(19,:),'.','markersize',vems*(norm(pov2)/norm(Y(17:19,:)-pov2)),'color',vec);
        subplot(3,1,2);plot3(Y(25,:),Y(26,:),Y(27,:),'.','markersize',eams*(norm(pov2)/norm(Y(25:27,:)-pov2)),'color',eac);
        subplot(3,1,2);plot3(Y(33,:),Y(34,:),Y(35,:),'.','markersize',mams*(norm(pov2)/norm(Y(33:35,:)-pov2)),'color',mac);
        subplot(3,1,2);plot3(Y(41,:),Y(42,:),Y(43,:),'.','markersize',jums*(norm(pov2)/norm(Y(41:43,:)-pov2)),'markeredgecolor',juc);
        subplot(3,1,2);plot3(Y(49,:),Y(50,:),Y(51,:),'.','markersize',sams*(norm(pov2)/norm(Y(49:51,:)-pov2)),'markeredgecolor',sac);
        subplot(3,1,2);plot3(Y(57,:),Y(58,:),Y(59,:),'.','markersize',urms*(norm(pov2)/norm(Y(57:59,:)-pov2)),'markeredgecolor',urc);
        subplot(3,1,2);plot3(Y(65,:),Y(66,:),Y(67,:),'.','markersize',nems*(norm(pov2)/norm(Y(65:67,:)-pov2)),'markeredgecolor',nec);hold off;
        axis([-axis2 axis2 -axis2 axis2 -axis2 axis2]);
        subplot(3,1,3);plot3(BSx.*axis3,BSy.*axis3,BSz.*axis3,'.','color',stc);hold on;
        subplot(3,1,3);plot3(Y(1,:),Y(2,:),Y(3,:),'.','markersize',sums/1.75,'markeredgecolor',suc);
        subplot(3,1,3);plot3(Y(41,:),Y(42,:),Y(43,:),'.','markersize',jums*(norm(pov3)/norm(Y(41:43,:)-pov3)),'markeredgecolor',juc);
        subplot(3,1,3);plot3(Y(49,:),Y(50,:),Y(51,:),'.','markersize',sams*(norm(pov3)/norm(Y(49:51,:)-pov3)),'markeredgecolor',sac);
        subplot(3,1,3);plot3(Y(57,:),Y(58,:),Y(59,:),'.','markersize',urms*(norm(pov3)/norm(Y(57:59,:)-pov3)),'markeredgecolor',urc);
        subplot(3,1,3);plot3(Y(65,:),Y(66,:),Y(67,:),'.','markersize',nems*(norm(pov3)/norm(Y(65:67,:)-pov3)),'markeredgecolor',nec);hold off;
        axis([-axis3 axis3 -axis3 axis3 -axis3 axis3]);
        whitebg(1,'k')
        set(gcf,'Color','k')
        pause(0.000001);
        
    end
end
