%run this one
clear
clc
clf
global Mj mio meuropa mganymede G Yt
Mj=1.8986e27;
mio=8.93e22;
meuropa=4.8e22;
mganymede=1.48e23;
G=6.67e-11;
%%
%Semi-Major axes in meters
io_axis=4.2e8;
europa_axis=6.709e8;
ganymede_axis=1.0692e9;
%%
%Initial positions
r1_io=[-io_axis 0 0];
r1_eu=[europa_axis 0 0];
r1_gy=[ganymede_axis 0 0];
r1_j=[0 0 0];
rcm=((mio*r1_io)+(meuropa*r1_eu)+(mganymede*r1_gy))/(Mj+mio+meuropa+mganymede);
r_io=r1_io-rcm;
r_eu=r1_eu-rcm;
r_gy=r1_gy-rcm;
r_j=r1_j-rcm;

omega_io=[0 0 sqrt(G*(Mj+mio+meuropa+mganymede)/norm(r_io)^3)];
omega_eu=[0 0 sqrt(G*(Mj+mio+meuropa+mganymede)/norm(r_eu)^3)];
omega_gy=[0 0 sqrt(G*(Mj+mio+meuropa+mganymede)/norm(r_gy)^3)];
omega_j=[0 0 sqrt(G*(Mj+mio+meuropa+mganymede)/norm(r_j)^3)];
w_io=omega_io(3);

%%
%Initial velocities
v_io=cross(omega_io,r_io);
v_eu=cross(omega_eu,r_eu);
v_gy=cross(omega_gy,r_gy);
%v_j=cross(omega_gy,r_j)+[0 0.4273 0]; %.4273 keeps the system from flying away
v_j=[0 -0.400908147163966 0]; %This value keeps Jupiter's orbit from flying away
%%
%Setting up Y vector
y=[r_io(1),r_io(2),v_io(1),v_io(2),r_eu(1),r_eu(2),v_eu(1),v_eu(2),...
    r_gy(1),r_gy(2),v_gy(1),v_gy(2),r_j(1),r_j(2),v_j(1),v_j(2)]';
t=0;

P=2*pi/w_io;
N=1000;
dt=P/N;
n=8*N; %# of Io orbits
%%
%Loop
for i=1:n

    f1 = moon_ode(t     , y                );
    f2 = moon_ode(t+dt/2, y+f1 .*  (dt/2)  );
    f3 = moon_ode(t+dt/2, y+f2 .*  (dt/2)  );
    f4 = moon_ode(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
    Yt(i+1,:) = y'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    
end
ecc=eccentricity(Yt);

for i=1:n
    %%
    %Plots
    q=n/100;
    if mod(i,q)==0
        j=ceil(i/1000);
    subplot(1,3,1)
    
    plot(Yt(i,1),Yt(i,2),'og',Yt(i,5),Yt(i,6),'ob',Yt(i,9),...
        Yt(i,10),'or',Yt(i,13),Yt(i,14),'ok',0,0,'X',Yt(:,1),Yt(:,2),'-g',...
        Yt(:,5),Yt(:,6),'-b',Yt(:,9),Yt(:,10),'-r',Yt(:,13),Yt(:,14),'-k');
    xlabel('Galilean Moon Orbits','fontsize',10);
    legend('Io','Europa','Ganymede','Jupiter','Location','North')
    axis equal
    
    subplot(1,3,2)
    
    plot(Yt(i,13),Yt(i,14),'ok',Yt(:,13),Yt(:,14),'-k',0,0,'X')
    xlabel('Jupiter''s Orbit','fontsize',10);
    legend('Jupiter','Location','North');
    axis equal
    
    subplot(1,3,3)
    
    plot(1:j,ecc(1,1:j),'-go',1:2:j,ecc(2,1:2:j),'-bo',1:4:j,ecc(3,1:4:j),'-ro')
    xlabel('Eccentricities of Moon Orbits','fontsize',10);
    legend('Io','Europa','Ganymede','Location','North')
    axis([0 8 0 .003])
    grid on
    
    k=i/(q);
    saveas(gcf,['zzzz',sprintf('%03d',k),'.png'])
    end
end