% This is a modified version of Nbody_test1 to investigate
% sensitivity to initial conditions in the 3-body problem.
% CCK 2009-Feb-24

function Nbody_test2

Nbodies = 3;

% 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 = 20; % 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 = 4000;
tolerance = [1e-4, 1e-4];
maximum_recursion_depth = 10;

% Initialize the state vector
q0 = randn(1,4*Nbodies);  % random initial state

% Modify the state vector by transforming into the center-of-mass frame
[K, U, E, Px, Py] = NbodyInvariants(q0, m); % calculate momentum & energy
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;
q0( (1:Nbodies) + (3*Nbodies) ) = q0( (1:Nbodies) + (3*Nbodies) ) - v_cm;

% Modify the state vector to put the system (probably) in a bound state.
[K, U, E, Px, Py] = NbodyInvariants(q0, m) % calculate momentum & energy
correction = sqrt(abs(0.4*U/K))
q0( (1:Nbodies) + (2*Nbodies) ) = q0( (1:Nbodies) + (2*Nbodies) ) * correction;
q0( (1:Nbodies) + (3*Nbodies) ) = q0( (1:Nbodies) + (3*Nbodies) ) * correction;
[K, U, E, Px, Py] = NbodyInvariants(q0, m)

% 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

% Perturb the initial conditions, but not by much!
perturbation = (1 + sqrt(tolerance(2))*randn(1,4*Nbodies));
q0p = q0 .* perturbation;

% Integrate the equations of motion again.
[qp,tp, greatest_depth_p] = adaptive_stepper(ODE, q0p, t0, tfinal, Nsteps,... 
   stepper, order, tolerance, maximum_recursion_depth);


% 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))
   hold on
   % plot(qp(1,body), qp(1,body+Nbodies), [style(body+Nbodies) '+'])
   plot(qp(:,body), qp(:,body+Nbodies), style(body+Nbodies))
end
title('Gravitational N-body simulation','fontsize',18)
legend('body 1','body 1 perturbed',...
       'body 2','body 2 perturbed',...
       'body 3','body 3 perturbed')
xlabel('x','fontsize',18)
ylabel('y','fontsize',18)
axis equal  % isotropy is desirable.
hold off
