function xyzfinal = nbody(n, steps)
% NBODY : Simulate n-body problem
%
% xyzfinal = nbody(n, steps)
%
% We simulate Newton's equations of motion in steps of dt seconds,
% drawing a picture after every time step.
%
% n = 0 or omitted simulates the solar system (sun and 9 planets)
% n = 3 simulates a stable reel of 3 bodies
% other 1 <= n <= 10 simulates sun and n-1 planets
% otherwise generate |n| bodies at random to simulate

% John R. Gilbert, January 2010, last updated 30 January 2017

if nargin < 2
    steps = 200;           % number of time steps to simulate
end;
if nargin < 1 || n == 0
    n = 10;
end;
dim = 1e13;        % size of the area to plot
if n == 3
    [mass,xyz,vxyz] = gen_3body;    % generate the reel of 3
elseif 0 < n && n <= 10
    [mass,xyz,vxyz] = gen_solar(n); % generate the solar system
    dim = 1.5*norm(max(abs(xyz)));
else
    n = abs(n);
    [mass,xyz,vxyz] = gen_nbody(n); % generate random initial data
end;

% mass(i)     = mass of i-th body (kg)
% xyz(i,1:3)  = initial position of i-th body (m in x,y,z)
% vxyz(i,1:3) = initial velocity of i-th body (m/s in x,y,z
                          
G = 6.674e-11;             % gravitational constant, m^3 kg^-1 s^-2
week = 60 * 60 * 24 * 7;   % one week, in seconds
dt = 1 * week;             % size of a time step
pix = 1;                   % 1 to draw pictures, 0 not to

if pix > 0, plot_pix(xyz,0); end;

for s = 1:steps
    axyz = zeros(n,3);   % acceleration of each body at this time step
    for i = 1:n
        for j = 1:n
            dxyz = xyz(i,:) - xyz(j,:);  % (x,y,z) vector from body j to i
            r = norm(dxyz);              % distance between bodies i and j
            if r == 0, continue, end;    % no force on a body from itself
            if ~isfinite(r), continue, end;   % ignore bodies at infinite dist
            f = G * mass(i) * mass(j) / r^2;  % force as a scalar
            fxyz = f * dxyz/r;                % force as a vector
            axyz(i,:) = axyz(i,:) - fxyz / mass(i);  % update to acceleration
            if ~all(isfinite(axyz(i,:)))
                axyz(i,:) = [0 0 0];      % don't update from infs and nans
            end;
        end;
        
    end;
    vxyz = vxyz + dt * axyz;    % update to velocity
    xyz = xyz + dt * vxyz;      % update to position
    if pix > 0, plot_pix(xyz,s,dim); end;
end;
xyzfinal = xyz;
end



function plot_pix(v,t,dim)
view = [1 0 0 ; 0 1 0]; % projection matrix for 2D pictures
if nargin < 3, dim = 1e13; end;
plot(v(:,1),v(:,2),'r.','MarkerSize',15);
hold on;
plot(v(1,1),v(1,2),'y.','MarkerSize',15);
axis([-dim dim -dim dim]);
axis square
title(sprintf('step %d',t));
hold off
shg;
pause(.02);
end
