 % ALL UNITS ARE METRIC  (kg, m, s)
 % total number of big objects in our solar system
 Nsphere  = 10;

 % Gravitational constant in kg^-2*m^3*s^-2
 global Grav;
 Grav     = 6.67259*(10^(-11));
 
 % mass of planets in kg
 global Mass;
 Mass     = zeros(Nsphere,1);
 Mass(1)  = 0.33*(10^24);       % Mercury
 Mass(2)  = 4.87*(10^24);       % Venus  
 Mass(3)  = 5.97*(10^24);       % Earth  
 Mass(4)  = 0.642*(10^24);      % Mars   
 Mass(5)  = 1899*(10^24);       % Jupiter
 Mass(6)  = 568*(10^24);        % Saturn 
 Mass(7)  = 86.8*(10^24);       % Uranus 
 Mass(8)  = 102*(10^24);        % Neptune
 Mass(9)  = 0.0125*(10^24);     % Pluto  
 Mass(10) = 1047.3486*Mass(5);  % The sun 

 % initial position in m
 X        = zeros(Nsphere,1);
 Y        = zeros(Nsphere,1);
 X(1)     =   57.9*(10^9);      % Mercury
 X(2)     =  108.2*(10^9); 	% Venus  
 X(3)     =  149.6*(10^9);	% Earth  
 X(4)     =  227.9*(10^9); 	% Mars   
 X(5)     =  778.6*(10^9);	% Jupiter
 X(6)     = 1433.5*(10^9);	% Saturn 
 X(7)     = 2872.5*(10^9);	% Uranus 
 X(8)     = 4495.1*(10^9);	% Neptune
 X(9)     = 5870.0*(10^9);	% Pluto  
 X(10)    = 0;			% The sun

 % initial velocity in m*s^-1 (using mean orbital velocity)
 VX       = zeros(Nsphere,1);
 VY       = zeros(Nsphere,1);
 VY(1)    = 47.89*(10^3); 
 VY(2)    = 35.03*(10^3); 
 VY(3)    = 29.79*(10^3); 
 VY(4)    = 24.13*(10^3); 
 VY(5)    = 13.06*(10^3); 
 VY(6)    =  9.64*(10^3); 
 VY(7)    =  6.80*(10^3); 
 VY(8)    =  5.40*(10^3); 
 VY(9)    =  4.70*(10^3);
 VY(10)   =  0;

 % final time is 3 years
 TimeLimits = [0 3*(365*24*60*60)];

 % initial position and velocity
 CoordVel = zeros(Nsphere, 4);
 CoordVel(:,1) = X;
 CoordVel(:,2) = Y;
 CoordVel(:,3) = VX;
 CoordVel(:,4) = VY;

 [T,CoordVel] = ode23('forcing', TimeLimits, CoordVel(:));

 % find the number of time steps taken
 [Tsteps,dummy] = size(T);
 
 % plot first 4 planet motions
 CoordVel =  reshape(CoordVel, Tsteps, Nsphere, 4);
 X  = CoordVel(:,1:4,1);
 Y  = CoordVel(:,1:4,2);
 VX = CoordVel(:,1:4,3);
 VY = CoordVel(:,1:4,4);

 title(sprintf('day=%lf', TimeLimits(2)));
 Tcol = zeros(Tsteps, 4);
 for i=1:4
   Tcol(:,i) = T(:);
 end
 scatter(X(:),Y(:),'o',Tcol(:));  
 







