% A test of adaptive rk23 (rk2.m, twostep.m, adaptive.m, adaptive_stepper.m)
% on the N-body problem with randomized initial conditions.
% CCK 2009-Feb-24

function Nbody_test1(Nbodies)

% Setup for the particular ODE problem
ODE = @Nbody;
m = exp( randn(1,Nbodies) );  % random masses, log-normally distributed
Nbody(0,0,m);
t0 = 0; % initial time
tfinal = 50; % final time

% Setup for adaptive solver,  using default tolerance and recursion settings:
stepper = @rk2;
order = 2;
% These next 3 parameters I keep changing to try and push to higher Nbodies...
Nsteps = 10000;
tolerance = [1e-3, 1e-3];
maximum_recursion_depth = 20;

% Initialize the state vector
q0 = randn(1,4*Nbodies);  % random initial state
[K, U, E, Px, Py] = NbodyInvariants(q0, m); % calculate momentum & energy

% Modify the state vector by transforming into the center-of-mass frame,
% and also halving the velocities to increase the likelihood of a bound state.
u_cm = Px / sum(m); % center of mass velocity (x component)
v_cm = Py / sum(m); % center of mass velocity (y component)
q0( (1:Nbodies) + (2*Nbodies) ) = (q0( (1:Nbodies) + (2*Nbodies) ) - u_cm)/2;
q0( (1:Nbodies) + (3*Nbodies) ) = (q0( (1:Nbodies) + (3*Nbodies) ) - v_cm)/2;


% Integrate the equations of motion.
[q,t, greatest_depth] = adaptive_stepper(ODE, q0, t0, tfinal, Nsteps,... 
   stepper, order, tolerance, maximum_recursion_depth);

% Work out time series of system energies and momenta.
for n=1:Nsteps+1  % includes initial conditions + subsequent time steps.
   [K(n), U(n), E(n), Px(n), Py(n)] = NbodyInvariants( q(n,:), m );
end

% Plots
figure(1)
hold off
style = ['b-';'k-';'r-';'g-';'m-';'c-';'y-']; % Array of line styles
for body = 1:Nbodies
   plot(q(1,body), q(1,body+Nbodies), [style(body) '+'])
   hold on
   plot(q(:,body), q(:,body+Nbodies), style(body))
end
title('Gravitational N-body simulation','fontsize',18)
% legend('10 quality-controlled RK23 steps','Analytic soln')
xlabel('x','fontsize',18)
ylabel('y','fontsize',18)
axis equal  % isotropy is desirable.
hold off

figure(2)
plot(t, U, 'b-', t, K, 'g-', t, E, 'r-')
title('Gravitational N-body simulation','fontsize',18)
xlabel('Time','fontsize',18)
ylabel('Energy','fontsize',18)
legend('Potential','Kinetic','Total')

figure(3)
plot(t, E, 'r-')
title('Gravitational N-body simulation','fontsize',18)
xlabel('Time','fontsize',18)
ylabel('Total Energy','fontsize',18)

figure(4)
plot(t, Px, 'b-', t, Py, 'g-')
title('Gravitational N-body simulation','fontsize',18)
xlabel('Time','fontsize',18)
ylabel('Momentum','fontsize',18)
legend('P_x','P_y')
