% Gravitational N-body equations of motion.
% CCK 2009-Feb-23
function dqdt = Nbody(q, t, mass)
% q is the state, t is time (t is not used, but is required for 
% compatibility with ODE solvers).
% The state vector is arranged as follows:
%     q = [ x1, x2,... xN,...
%           y1, y2,... yN,...
%           u1, u2,... uN,...
%           v1, v2,... vN ]
% The gravitational constant is taken to be G=1.
%
% mass is the N-element array of masses. It is typically specified
% once. When mass is specified, q and t are unused, and the function
% returns with a trivial (text confirmation message) result. Example:
%     octave:1> Nbody(0,0,[1,2,3])
%     mass =
%
%        1   2   3
%
%     ans = mass array loaded successfully.

persistent m; % persistent means m will be available next time Nbod is called.
if exist('mass','var')
   m = mass;
   dqdt = 'mass array loaded successfully.'; % Reassure the user.
   mass % Display the mass array that was loaded.
   return % return without evaluating any forces.
end
N = numel(m); % number of bodies


dqdt( 1:N ) = q( (1:N)+(2*N) );    % dx/dt = u
dqdt( (1:N)+N ) = q( (1:N)+(3*N) ); % dy/dt = v

for n = 1:N % Loop through all N bodies.
   % Find sum of forces on nth body:
   rx = q(n) - q(1:N);          % r vector x-components
   ry = q(n+N) - q( (1:N)+N );  % r vector y-components
   r3 = (rx.^2 + ry.^2).^(3/2); % |r|^3
   r3(n) = inf; % This will zero the self-force term.
   dqdt( n+(2*N) ) = -sum( m .* rx ./ r3 );  % dvx/dt
   dqdt( n+(3*N) ) = -sum( m .* ry ./ r3 );  % dvy/dt
end