% Simulation
% By:  Jay Borthen
% A multi-point-mass gravitational system

clear all
close all

% The Earth has a mass of 5.9736x10^24 kg
M1=550000000;  % kg
M2=125000000;  % kg
M3=315000000;  % kg
M4=75500000;   % kg
M5=1005000000;  % kg
G=6.67384*10^(-11);  %N*(m/kg)^2

P_M1=[23 20];  % Initial point in space of mass 1
%V_M1=[0.05 0.03];  % Initial velocity vector of mass 1
V_M1=[0 0];

P_M2=[10 15];  % Initial point in space of mass 2
%V_M2=[-0.02 0.1];  % Initial velocity vector of mass 2
V_M2=[0 0];

P_M3=[30 25];  % Initial point in space of mass 3
%V_M3=[-0.03 -0.2];  % Initial velocity vector of mass 3
V_M3=[0 0];

P_M4=[5 35];  % Initial point in space of mass 4
%V_M4=[0.1 -0.05];  % Initial velocity vector of mass 4
V_M4=[0 0];

P_M5=[18 40];  % Initial point in space of mass 5
%V_M5=[0.05 0.3];  % Initial velocity vector of mass 5
V_M5=[0 0];

for t=1:5000  % seconds

rVec12=P_M2-P_M1;
rVec13=P_M3-P_M1;
rVec14=P_M4-P_M1;
rVec15=P_M5-P_M1;

rVec21=P_M1-P_M2;
rVec23=P_M3-P_M2;
rVec24=P_M4-P_M2;
rVec25=P_M5-P_M2;

rVec31=P_M1-P_M3;
rVec32=P_M2-P_M3;
rVec34=P_M4-P_M3;
rVec35=P_M5-P_M3;

rVec41=P_M1-P_M4;
rVec42=P_M2-P_M4;
rVec43=P_M3-P_M4;
rVec45=P_M5-P_M4;

rVec51=P_M1-P_M5;
rVec52=P_M2-P_M5;
rVec53=P_M3-P_M5;
rVec54=P_M4-P_M5;



plot(P_M1(1),P_M1(2),'b.');
hold on
plot(P_M2(1),P_M2(2),'r.')
plot(P_M3(1),P_M3(2),'g.')
plot(P_M4(1),P_M4(2),'w.')
plot(P_M5(1),P_M5(2),'m.')
axis([-25 70 -10 70])
set(gca,'Color',[0,0,0]);
hold on





r12=sqrt((rVec12(1)^2)+(rVec12(2)^2));
F12=G.*(M1.*M2)/(r12.^2);

r13=sqrt((rVec13(1)^2)+(rVec13(2)^2));
F13=G.*(M1.*M3)/(r13.^2);

r14=sqrt((rVec14(1)^2)+(rVec14(2)^2));
F14=G.*(M1.*M4)/(r14.^2);

r15=sqrt((rVec15(1)^2)+(rVec15(2)^2));
F15=G.*(M1.*M5)/(r15.^2);




r21=sqrt((rVec21(1)^2)+(rVec21(2)^2));
F21=G.*(M2.*M1)/(r21.^2);

r23=sqrt((rVec23(1)^2)+(rVec23(2)^2));
F23=G.*(M2.*M3)/(r23.^2);

r24=sqrt((rVec24(1)^2)+(rVec24(2)^2));
F24=G.*(M2.*M4)/(r24.^2);

r25=sqrt((rVec25(1)^2)+(rVec25(2)^2));
F25=G.*(M2.*M5)/(r25.^2);





r31=sqrt((rVec31(1)^2)+(rVec31(2)^2));
F31=G.*(M3.*M1)/(r31.^2);

r32=sqrt((rVec32(1)^2)+(rVec32(2)^2));
F32=G.*(M3.*M2)/(r32.^2);

r34=sqrt((rVec34(1)^2)+(rVec34(2)^2));
F34=G.*(M3.*M4)/(r34.^2);

r35=sqrt((rVec35(1)^2)+(rVec35(2)^2));
F35=G.*(M3.*M5)/(r35.^2);





r41=sqrt((rVec41(1)^2)+(rVec41(2)^2));
F41=G.*(M4.*M1)/(r41.^2);

r42=sqrt((rVec42(1)^2)+(rVec42(2)^2));
F42=G.*(M4.*M2)/(r42.^2);

r43=sqrt((rVec43(1)^2)+(rVec43(2)^2));
F43=G.*(M4.*M3)/(r43.^2);

r45=sqrt((rVec45(1)^2)+(rVec45(2)^2));
F45=G.*(M4.*M5)/(r45.^2);



r51=sqrt((rVec51(1)^2)+(rVec51(2)^2));
F51=G.*(M5.*M1)/(r51.^2);

r52=sqrt((rVec52(1)^2)+(rVec52(2)^2));
F52=G.*(M5.*M2)/(r52.^2);

r53=sqrt((rVec53(1)^2)+(rVec53(2)^2));
F53=G.*(M5.*M3)/(r53.^2);

r54=sqrt((rVec54(1)^2)+(rVec54(2)^2));
F54=G.*(M5.*M4)/(r54.^2);



F1=F12.*rVec12+F13.*rVec13+F14.*rVec14+F15.*rVec15;
F2=F21.*rVec21+F23.*rVec23+F24.*rVec24+F25.*rVec25;
F3=F31.*rVec31+F32.*rVec32+F34.*rVec34+F35.*rVec35;
F4=F41.*rVec41+F42.*rVec42+F43.*rVec43+F45.*rVec45;
F5=F51.*rVec51+F52.*rVec52+F53.*rVec53+F54.*rVec54;


V_M1=V_M1+(F1/M1);
P_M1=P_M1+(V_M1)/2;

V_M2=V_M2+(F2/M2);
P_M2=P_M2+(V_M2)/2;

V_M3=V_M3+(F3/M3);
P_M3=P_M3+(V_M3)/2;

V_M4=V_M4+(F4/M4);
P_M4=P_M4+(V_M4)/2;

V_M5=V_M5+(F5/M5);
P_M5=P_M5+(V_M5)/2;


%grid on
pause(0.01)

end