% nbpRKN1210.m
%
% Numerically integrate a test particle in the circularly restricted
% three-body problem using the Runge-Kutta-Nystrom 12/10 integrator.
%
% Massive bodies M1 and M2 follow circular orbits around their center of 
% mass point, set to be the origin. The third test particle does affect the
% motion of the two massive bodies. 
%
% We assume G = 1 and R = 1 (distance between the primary bodies)
%
% This function is designed to be called by another function which sets 
% the initial conditions. 
%
% This function calls the following functions:
%   - rkn1210.m (written by Rody P.S. Oldenhuis, MATLAB File Exchange)
%   - LagrangePoints.m
%
% Passed parameters:
%   masses - 1x2 matrix containing m1 and m2
%   pos0   - 1x2N matrix cointaining pairs of (x,y) coordinates for
%            the initial conditions of the test particles. If one test
%            particle is used, pos0 will be a 1x2 matrix containing (x0, y0)
%            For 2 test particles, posq will contain (x10, y10, x20, y20)...
%   vel0   - 1x2N matrix containing pairs of (vx, vy) values for each of
%            the initial velocities of the test particles. 
%   times  - 1x2 array containing begin and end times [tBegin, tEnd]
%   flag   - 1x2 array containing plotting options (flags are optional):
%            flag(1) = true to animate the orbit as it integrates. Cool to
%                      watch, but it slows the calculate way down.
%            flag(2) = true to leave a trail behind the test particles when
%                      animation is turned on.
%
% Returned arrays:
%   t      - Tx1 array returning times used in integration, where T =
%            number of time steps
%   pos    - Tx2N array returning (x,y) coordinate pairs for each test
%            particle
%   vel    - Tx2N array returning (vx,vy) velocity component pairs for each 
%            test particle
%
%
% MATLAB-Monkey.com   10/6/2013


function returnVals = nbpRKN1210(m, pos0, vel0, times, flag)

if nargin<4 
    display('crtbpRKN1210 requires at least 4 parameters');
    exit;
end

    
%%%%%%%%%%                                                       %%%%%%%%%%
%%%%%%%%%%           Initialization and Parameters               %%%%%%%%%% 
%%%%%%%%%%                                                       %%%%%%%%%%


progressBar = false;     % if true, show progress bar

animateDelay = 0.0;      % delay (in seconds) between frames of animation.
                         % smaller the value, faster the animation will play
                        
leaveTrail = flag(2);    % if true and animation is on, past positions will
                         % not be erased.  If false and animation is on,
                         % no 'trail' of past positions will be left 
                        
trailLength = 100;      % if leaveTrail is true, this sets the length of
                         % trail


%%%%%%%%%%%%%%%%%%%%  Flags Controlling Output  %%%%%%%%%%%%%%%%%%%%%%%%%%%
 
if nargin >= 5
    animate = flag(1);       % if true, animate motion
                             % if false, just display final trajectory

    leaveTrail = flag(2);    % if true and animation is on, past positions will
                             % not be erased.  If false and animation is on,
                             % no 'trail' of past positions will be left 
end                 




%%%%%%%%%%%%%%%%%%%%%%%%%  System Parameters  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%


NP = length(pos0)/2;     % number of test particles.  Assumes 2D

G = 1;                   % Gravitational Constant

trails = ones(1,NP);     % if trails(i) == 1, trail is on
trails(1) = 0;           % turn off trail for particle 1
trails(2) = 0;           % turn off trail for particle 2


%%%%%%%%%%%%%%%%  Orbital properties of two massive bodies  %%%%%%%%%%%%%%%

plotLimits = 3.5;           % limit for plotting
tEnd = times(end);



%%%%%%%%%%%%%%%%%%%%%%%%%%%   Trail Properties  %%%%%%%%%%%%%%%%%%%%%%%%%%%

trail = ones(trailLength,2*NP);
for j = 1:2*NP
    trail(:,j) = pos0(j)*ones(trailLength,1);
end


%%%%%%%%%%                                                       %%%%%%%%%%
%%%%%%%%%%                      Integrate                        %%%%%%%%%% 
%%%%%%%%%%                                                       %%%%%%%%%%


% Use Runge-Kutta 45 integrator to solve the ODE

if animate
    options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12,'outputfcn',@OutputFcn2);
else
    % initialize waitbar
    if progressBar
        wait = waitbar(0, 'integrating...');
        options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12,'outputfcn',@OutputFcn1);
    else
        options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12);
    end
end

% Use Runge-Kutta-Nystrom integrator to solve the ODE
tic;
% [t,pos,vel] = rkn1210(@derivatives, times, pos0, vel0, options);
[t,pos,vel,TE,YE] = rkn1210(@derivatives, times, pos0, vel0, options);

DT = toc;


% returnVals = [t, pos, vel, YE];
returnVals = [t, pos, vel];




%%%%%%%%%%                                                       %%%%%%%%%%
%%%%%%%%%%                      Functions                        %%%%%%%%%% 
%%%%%%%%%%                                                       %%%%%%%%%%


    %%%%%%%%%%  Derivatives  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % Function defining derivatives dx/dt and dy/dt
    function dvdt = derivatives(tf,wf)
        NP = length(wf)/2;      % number of particles (assumes 2D)     
        
        dvdt = zeros(1,2*NP);
        
        for ii = 1:NP
            for jj = 1:NP
                if ii ~= jj
                    Xi = wf(2*ii-1);            % x component of ii^th particle
                    Yi = wf(2*ii);              % y component of ii^th particle

                    Xj = wf(2*jj-1);            % x component of jj^th particle
                    Yj = wf(2*jj);              % y component of jj^th particle

                    r1 = ((Xi-Xj)^2 + (Yi-Yj)^2)^(3/2);

                    dvdt(2*ii-1) = dvdt(2*ii-1) - m(jj) * (Xi-Xj) / r1;
                    dvdt(2*ii) =  dvdt(2*ii) - m(jj) * (Yi-Yj) / r1;
                end
            end
        end
    end



    %%%%%%%%%%  OutputFcn1:  Status Bar  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % the output function
    function stop = OutputFcn1(t,y,dy,flag)%#ok
        % don't stop
        stop = false;
        % only after sucessfull steps
        if ~isempty(flag), return
        else
            wait = waitbar(t/tEnd, wait);
        end
    end



    %%%%%%%%%%  OutputFcn2:  Animation  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % the output function for making the animation
    function stop = OutputFcn2(t,y,dy,flag)
        % don't stop
        stop = false;
        % only after sucessfull steps
        if ~isempty(flag), return
            
        else
            NP = length(y)/2;      % number of particles (assumes 2D)
            
            
            if leaveTrail
                for ii = 1:NP
                    if trails(ii)
                        trail(:,2*ii-1) = [y(2*ii-1); trail(1:end-1,2*ii-1)];
                        trail(:,2*ii) = [y(2*ii); trail(1:end-1,2*ii)];
                        
                        plot(trail(:,2*ii-1),trail(:,2*ii),'-');
                        hold all
                    end
                end
            end
            for ii = 1:NP
                plot(y(2*ii-1), y(2*ii), '.','MarkerSize',16)
                hold all
            end
        end
        
        axis equal
        axis(plotLimits*[-1 1 -1 1])
        hold off
        pause(animateDelay), drawnow % flush plotting commands
    end

end
