% SPM: Driver for Single Particle Motion in given E and B fields
% Variables to set to specify the problem are:
%    tspan is time interval 
%    nn is number of steps to take (unless using adaptive RK45
%         method)
%    y0 are initial conditions [x0 y0 z0 vx0 vy0 vz0]
% Derivatives of [x0 y0 z0 vx0 vy0 vz0] are computed using the
%    function lorentz
%    Within lorentz charge-to-mass ratio qoverm is specified
%    lorentz calls functions electric and magnetic to specify
%        fields as a function of [t y(6)]
% Numerical Method choices are made by setting variable meth:
%    Possible selections are:
%    1 Euler's method
%    2 Leapfrog method
%    3 Adaptive Runge-Kutta 45 method
%    4 Adaptive Runge-Kutta 45 method with RelTol specified
%

'======================================================================='
%Setting TimeSpan and Initial Conditions
tspan= [0. 2.*pi*10.]';
nn=1000;
y0 = [0. 0. 0. 1. 0. 0. ]';  %Specify column vector of initial conditions 

dt=(tspan(end)-tspan(1))/nn;
t=tspan(1):dt:tspan(end);%Calculate row vector of times
t=t';  %Convert to column vector of times

%Setting function handle to Lorentz Force Law function
lorhandle= @lorentz ;

%Specify method for integration
meth=1;

switch(meth)
 case 1
  %Calling Euler's method (first-order)
  method='Euler Method'
  [T,Y]=euler1(lorhandle,t,y0);
 case 2
  %Calling Leapfrog method (second-order)
  method='Leapfrog Method'
  [T,Y]=leapfrog2(lorhandle,t,y0);
 case 3
  %Calling ode45 for Runge-Kutta
  method='RK45 Method'
  [T,Y]=ode45(lorhandle,tspan,y0);
 case 4
  %Calling ode45 for Runge-Kutta with option for setting error tolerance
  options = odeset('RelTol',1.0e-5);
  method='RK45 Method with specified RelTol'
  [T,Y]=ode45(lorhandle,tspan,y0,options);
end

nsteps=size(T,1)

terr=( (Y(1,1)-Y(end,1))^2. + (Y(1,2)-Y(end,2))^2.)^0.5

plot(Y(:,1),Y(:,2))