% SPM: Driver for Single Particle Motion in given E and B fields
% EXAMPLE OF SWITCH-CASE STRUCTURE and FUNCTION CALLS
% Numerical Method choices are made by setting variable (meth):
%    Possible selections are:
%    1 Euler's method
%    2 Leapfrog method
% Input x0, v0, B0, E0,
% Compute x and v as a function of time
'======================================================================='

%Specify method for integration
meth=1;

%Setting TimeSpan and Initial Conditions
tspan= [0. 10.*2.*pi]';
nn=1000;
%Initial Conditions
x0=zeros(3,1);
v0=zeros(3,1);
x0 = [0. 1. 0. ]';  % Initial position (Column vector)
v0 = [1. 0. 0. ]';  % Initial velocity (Column vector)

%Parameters
q=1.;  %Charge
m=1.;  %Mass

%Electromagnetic Field Input
B0=[0. 0. 1.]';   % Initial Magnetic Field (Column vector)
E0=[0. 0.1 0.]';   % Initial Electric Field (Column vector)
Bmag=norm(B0);    % Magnitude of B0 (2-norm)

% Compute timestep
dt=(tspan(end)-tspan(1))/nn;
%Calculate row vector of times
t=(tspan(1):dt:tspan(end))';
%Compute total number of times
nt=size(t,1)
%Initial position and velocity vectors to be solved
x=zeros(3,nt);
v=zeros(3,nt);
%Set initial conditions into x and v
x(1:3,1)=x0(1:3);
v(1:3,1)=v0(1:3);

% Initialize time derivatives of x and v
%dxdt=zeros(3,1);
%dvdt=zeros(3,1);

switch(meth)
 case 1
  %Euler's method (first-order)
  method='Euler Method'
  % Loop over timesteps
  for i= 2:nt
    %Calculate time rate of change of x and v
    [dxdt,dvdt]=lorentz(q,m,x(1:3,i-1),v(1:3,i-1),E0,B0);
    %Update x and v using dxdt and dvdt
    x(1:3,i)= x(1:3,i-1) + dt*dxdt(1:3);
    v(1:3,i)= v(1:3,i-1) + dt*dvdt(1:3);
  end
case 2
  %Leapfrog method (second-order)
  method='Leapfrog Method'

  %Update with Leapfrog method, calling lorentz function as needed

end

%Set Figure Size
h1=figure();
%Plot trajectory
plot(x(1,:),x(2,:),'LineWidth',2.0,'Color','b')
hold on
% Plot Axes
xL = xlim;
yL = ylim;
line([0 0], yL,'LineWidth',1.0,'Color','k');  %y-axis
line(xL, [0 0],'LineWidth',1.0,'Color','k');  %x-axis
%Plot Labels
xlabel('Position, $x$','Interpreter','latex','FontName','TimesNewRoman','FontSize',20,'FontWeight','bold')
ylabel('Position, $y$','Interpreter','latex','FontName','TimesNewRoman','FontSize',20,'FontWeight','bold')
set(gca,'FontSize',16,'FontName','TimesNewRoman','FontWeight','bold','LineWidth',2)
title('Trajectory');


