% N-body simulation of motion under gravity alone
% Author: Andri M. Gretarsson
% Last Modified: 9/28/09
%
% This is a routine for calculating the trajectories of the N masses 
% in an N-body system interacting under gravity alone.  


%-------------------------
% Set Initial Values etc.
%-------------------------

% Physical Constants
G=6.67300e-11;          % Universal Gravitational Constant (SI units)

% Number of Steps and Time Vector
tend=60*3600*24;        % Simulation end time (seconds):  days * 24 hrs/day * 3600 sec/hr
dt=0.001*3600*24;        % Time step (seconds):  days * 24 hrs/day * 3600 sec/hr
Nsteps=round(tend/dt);  % Number of simulation steps
t=1:dt:Nsteps*dt;       % Vector of times at which position etc will be calculated

% Masses of bodies in the system
M=[5.9736e24 7.3477e22];   % Currently using Earth and Moon masses (kg)

% Allocating memory for matrixes of pos, vel, accel. Columns hold values
% for individual masses. (Currently we have 2 masses and so we have 2
% columns). There is one row for each time step at which pos, vel, accel.
% are calculated
x=zeros(Nsteps,length(M));      y=zeros(Nsteps,length(M));      z=zeros(Nsteps,length(M)); 
vx=zeros(Nsteps,length(M));     vy=zeros(Nsteps,length(M));     vz=zeros(Nsteps,length(M));
ax=zeros(Nsteps-1,length(M));   ay=zeros(Nsteps-1,length(M));   az=zeros(Nsteps-1,length(M));


% Initial values occur at the first time step and so they go in the first
% row. I chose vy so that the net linear momentum is zero (that way the
% system  as a whole doesn't travel). 

% Initial Position (meters):
x(1,:)=[0 384399000]; 
y(1,:)=[0 0];                           
z(1,:)=[0 0];

% Initial Velocity (meters/second):
vx(1,:)=[0 0]*1e-5;             
vy(1,:)=[-M(2)/M(1)*1022 1022];        
vz(1,:)=[0 0];   



%-----------------------
% Iterate Calculation
%-----------------------

% On each pass through the outer for-loop, we calculate the new
% position of the bodies in the simulation and increment the time by dt.
% The number of steps is (Nsteps-1) since the first step corresponds to initial
% conditions. So, we start at s=2.
for s=2:Nsteps
    
    % The inner for-loop calculates the acceleration of the k'th 
    % mass during the (s-1)'th timestep
    for k=1:length(M)
        % xtemp, Mtemp etc. hold the x-pos, Mass, etc of all the bodies
        % EXCEPT the one for which we need to calculate the acceleration
        % components. 
        xtemp=x(s-1,[1:k-1,k+1:end]);   ytemp=y(s-1,[1:k-1,k+1:end]);   ztemp=z(s-1,[1:k-1,k+1:end]);
        Mtemp=M([1:k-1,k+1:end]);
        
        % The next three lines implement Newton's Law of Gravitation.
        % We find the accel. component for the k'th body by summing the
        % contributions to the accel from each of the other bodies. So,
        % we're summing over the "temp" variables that hold the position
        % and mass for all the other bodies.
        ax(s-1,k)=sum(-G*Mtemp.* (x(s-1,k)-xtemp) ./ ((x(s-1,k)-xtemp).^2+(y(s-1,k)-ytemp).^2+(z(s-1,k)-ztemp).^2).^(3/2) );
        ay(s-1,k)=sum(-G*Mtemp.* (y(s-1,k)-ytemp) ./ ((x(s-1,k)-xtemp).^2+(y(s-1,k)-ytemp).^2+(z(s-1,k)-ztemp).^2).^(3/2) );
        az(s-1,k)=sum(-G*Mtemp.* (z(s-1,k)-ztemp) ./ ((x(s-1,k)-xtemp).^2+(y(s-1,k)-ytemp).^2+(z(s-1,k)-ztemp).^2).^(3/2) );
    end
    
    % Integrate the (s-1)'th acceleration components above to find the s'th velocity components
    vx(s,:) =   ax(s-1,:)*dt    +   vx(s-1,:);   
    vy(s,:) =   ay(s-1,:)*dt    +   vy(s-1,:);
    vz(s,:) =   az(s-1,:)*dt    +   vz(s-1,:);
    
    % Integrate the (s-1)'th velocity components above to find the s'th
    % position components.
    % (The commented out portions of the calculation represent a trick for
    % much better accuracy for this particular problem. Try uncommenting
    % them if are interested, but complete the homework first.)
    x(s,:)  =   vx(s-1,:)*dt    +   x(s-1,:)   ;% +   dt^2*ax(s-1,:);   
    y(s,:)  =   vy(s-1,:)*dt    +   y(s-1,:)   ;% +   dt^2*ay(s-1,:);
    z(s,:)  =   vz(s-1,:)*dt    +   z(s-1,:)   ;% +   dt^2*az(s-1,:);
end

% Calculate magnitudes in case we want to look at them later
r=sqrt((x.^2+y.^2+z.^2));
v=sqrt(vx.^2+vy.^2+vz.^2);
a=sqrt(ax.^2+ay.^2+az.^2);




%-----------------------
% Plot the trajectories
%-----------------------

figure(1);
plot3(x/1000,y/1000,z/1000,'.','linewidth',2,'markersize',1);       
axis([-1 1 -1 1]*5e5);
axis square

figure(2);
plot(t/3600/24,r/1000,'.','linewidth',2,'markersize',1);
axis square
grid on;
shg    
figure(1);