clf
clear all
close all

global mS mX G mJ mSat mN mU
mS = (2e+30);
mX=(6e+25); %hypothetical mass of Planet X
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
%Planet X distance from Sun
RX=2.244e+13;
%RX=5.236e+13;
RJ=778330000000; %jupiter-Sun
RSat=1.436e+12; %Saturn-Sun
RN=4498252900000; %Neptune-Sun
RU=3.2e12; %Uranus-Sun

theta6=30; %angle of orbit of Planet X
theta7=0.32; %Jupiter
theta8=0.93; %saturn
theta9=0.72; %Neptune
theta10=1.02; %Uranus

rSx=0;
rSy=0;
rXx=RX*cosd(theta6);
rXy=RX*sind(theta6);
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 + mX*rXx + mJ*rJx + mSat*rSatx + mN*rNx + mU*rUx)/(mS+mX+mJ+mSat+mN+mU);
cmY=(0 + mX*rXy + mJ*rJy + mSat*rSaty + mN*rNy + mU*rUy)/(mS+mX+mJ+mSat+mN+mU);
rCM=[cmX cmY]; %center of mass

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

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

%angular velocities
wX=sqrt((G*mX*mS)/(muX*norm(rX)^3));
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));

% initial velocities
vS  =[0;0];
vX  =[0;(wX*norm(rX))]; 
vJ  =[0;(wJ*norm(rJ))];
vSat=[0;(wSat*norm(rSat))];
vN  =[0;(wN*norm(rN))];
vU  =[0;(wU*norm(rU))];


Y=[rX(1);rX(2);vX(1);vX(2);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 / 100;
n = 17000*2; % 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_planet9(t     , Y );
    f2 = kepler_ode_planet9(t+dt/2, Y+f1 .*  (dt/2));
    f3 = kepler_ode_planet9(t+dt/2, Y+f2 .*  (dt/2));
    f4 = kepler_ode_planet9(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,1000)==0
        Yt(i+1,:) = Y';
        %f=f+1
         imagesc([-1e+14 1e+14],[-1e+14 1e+14],image);
         hold on;
        set(gca,'ydir','normal');
        plot(0,0,'X')
        hold on
        plot(Yt(i+1,1),Yt(i+1,2),'.r','markersize',30) %plots planet X trajectory
        hold on
        plot(Yt(i+1,5),Yt(i+1,6),'y*','markersize',5) %plots Sun position
        hold on
%         plot(Yt(i+1,9),Yt(i+1,10),'g.','markersize',15) %plots Jupiter
%         hold on
%         plot(Yt(i+1,13),Yt(i+1,14),'m.','markersize',15) %plots Saturn
%         hold on
        plot(Yt(i+1,17),Yt(i+1,18),'b.','markersize',15) %plots Neptune
        hold on
%         plot(Yt(i+1,21),Yt(i+1,22),'c.','markersize',15) %plots Uranus
        legend('Centre of mass','Planet 9','Sun','Neptune','Location','northeastoutside','Orientation','vertical')
        legend('boxoff')
%         axis([-1.5e+14 1.5e+14 -10e+13 30e+13])
        title('Orbits of Neptune and Planet 9 around the Sun','Fontsize',9)
         axis([-1e+14 1e+14 -1e+14 1e+14]);
        set(gca,'xtick',[])
        set(gca,'xticklabel',[])
        set(gca,'ytick',[])
        set(gca,'yticklabel',[])
        drawnow
        %saveas(gcf, ['Project_Neptune_close_reduced',num2str(f,'%03d'),'.png'])
    end
end


