%modelling orbits of 4 gas giants only without planet 9

clf
clear all
close all
global mS G mJ mSat mN mU
mS = (2e+30);
G  = 6.7e-11; %gravitational constant
mJ=1.8986e+27;
mSat=5.6836e+26;
mN=1.0243e+26;
mU=8.681e+25;

f=0; %to be able to make a movie

RJ=778330000000; %jupiter-Sun
RSat=1.436e+12; %Saturn-Sun
RN=4498252900000; %Neptune-Sun
RU=3.2e12; %Uranus-Sun

theta7=0.32; %Jupiter
theta8=0.93; %saturn
theta9=0.72; %Neptune
theta10=1.02; %Uranus

rJx=RJ*cosd(theta7);
rJy=RJ*sind(theta7);
rSatx=RSat*cosd(theta8);
rSaty=RSat*sind(theta8);
rNx=RN*cosd(theta9);
rNy=RN*sind(theta9);
rUx=RU*cosd(theta10);
rUy=RU*sind(theta10);

cmX=(0  + mJ*rJx + mSat*rSatx + mN*rNx + mU*rUx)/(mS+mJ+mSat+mN+mU);
cmY=(0  + mJ*rJy + mSat*rSaty + mN*rNy + mU*rUy)/(mS+mJ+mSat+mN+mU);
rCM=[cmX cmY]; %center of mass

rS  =[0 0]-rCM; %position Sun relative to the center of mass
rJ  =[rJx rJy]-rCM;
rSat =[rSatx rSaty]-rCM;
rN  =[rNx rNy]-rCM;
rU  =[rUx rUy]-rCM;

%calculating angular velocities
%mu
muJ=(mJ*mS)/(mJ+mS);
muSat=(mSat*mS)/(mSat+mS);
muN=(mN*mS)/(mN+mS);
muU=(mU*mS)/(mU+mS);

%angular velocities
wJ=sqrt((G*mJ*mS)/(muJ*norm(rJ)^3));
wSat=sqrt((G*mSat*mS)/(muSat*norm(rSat)^3));
wN=sqrt((G*mN*mS)/(muN*norm(rN)^3));
wU=sqrt((G*mU*mS)/(muU*norm(rU)^3));

vS  =[0;0];
vJ  =[0;(wJ*norm(rJ))];
vSat=[0;(wSat*norm(rSat))];
vN  =[0;(wN*norm(rN))];
vU  =[0;(wU*norm(rU))];

Y=[rS(1);rS(2);vS(1);vS(2);rJ(1);rJ(2);vJ(1);vJ(2);rSat(1);rSat(2);vSat(1);vSat(2);rN(1);rN(2);vN(1);vN(2);rU(1);rU(2);vU(1);vU(2)];
f =  0;

t =  0.0;
dt= 2*pi/wJ / 10000;
n = 110000; % number of steps
tt = zeros(n+1,1);
Yt = zeros(n+1,size(Y,1));
tt(1)   = t;
Yt(1,:) = Y';
image= imread('space.jpg');

for i=1:n
    f1 = kepler_ode_giant(t     , Y );
    f2 = kepler_ode_giant(t+dt/2, Y+f1 .*  (dt/2));
    f3 = kepler_ode_giant(t+dt/2, Y+f2 .*  (dt/2));
    f4 = kepler_ode_giant(t+dt  , Y+f3 .*   dt);
    
    
    Y = Y + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    
    t = t+dt;
    tt(i+1)   = t;  % this records the time at each step
    
    if mod(i,1500)==0
        Yt(i+1,:) = Y';
        f=f+1
        imagesc([-5e+12 5e+12], [-5e+12 5e+12],image);
        hold on;
        plot(0,0,'X')
        hold on
        plot(Yt(i+1,1),Yt(i+1,2),'y*','markersize',15) %plots Sun position
        hold on
        plot(Yt(i+1,5),Yt(i+1,6),'g.', 'markersize',30) %plots Jupiter
        hold on
        plot(Yt(i+1,9),Yt(i+1,10),'m.','markersize',25) %plots Saturn
        hold on
        plot(Yt(i+1,13),Yt(i+1,14),'b.','markersize',15) %plots Neptune
        hold on
        plot(Yt(i+1,17),Yt(i+1,18),'c.','markersize',15) %plots Uranus
        axis([-5e+12 5e+12 -5e+12 5e+12])
        title('Orbits of the four gas giants around the Sun','Fontsize',9)
        legend('Centre of mass','Sun','Jupiter','Saturn','Neptune','Uranus','Location','northeastoutside','Orientation','vertical')
        legend('boxoff')
        set(gca,'xtick',[])
        set(gca,'xticklabel',[])
        set(gca,'ytick',[])
        set(gca,'yticklabel',[])
        drawnow
        saveas(gcf, ['Project_giant_reduced',num2str(f,'%03d'),'.png'])
        
    end
end

