% Created by Ethan Deyle
% January 31, 2008
% Calculates motion of a particle of charge q
% in an electric field E and magnetic field B
% created for Plasma Physics Seminar, Swarthmore College
% solves for Bellan problem 3.2

function [] = orbital2

hold off;clf; clear; 

% Input Electric and Magnetic Fields


% Input parameters and initial conditions

q = 5;
m = 1;
X0 = [0,30,0];
vX0 = [0,0,0];

% Book keeping

E = @EE;
B = @BB;

cc = q/m;

X = X0;
vX = vX0;
t = 0;

Xarray = X';
T = [t];


h = .0003;
M = 1.0;
figure(1)
hold on

counter = 0;
while counter <= M
fX1 = vX;
g1 = FF(X,vX, E, B, t, cc);
fX2 = vX + h*(g1)/2;
g2 = FF((X + (1/2)*fX1), fX2, E, B, t, cc);
fX3 = vX + h*(g2)/2;
g3 = FF((X + h*(1/2)*fX2), fX3, E, B, t, cc);
fX4 = vX + h*g3;
g4 = FF((X+h*fX3), fX4, E, B, t, cc);
	X = X + 1/6*(fX1 + 2*fX2 + 2*fX3 + fX4);
    vX = vX + h*(g1+2*g2+2*g3+g4)/6; 
	t = t+h;
	T = [T,t];
	Xarray = [Xarray,X'];
counter = counter + h;
end

Xarray = Xarray';

scatter3(Xarray(:,1), Xarray(:,2), Xarray(:,3),10,'o')
title('Particle Orbit');
grid;

% Define electric and magnetic fields
function Bout = BB(X,t)
rho = sqrt(X(1)^2 + X(2)^2);
%Bout = [0,0,10*sin(t)];
Bout = [0,0,X(2)/15];

function Eout = EE(X,t)
r = sqrt(dot(X,X));
Eout = [1,0,0];

% define the function v'(t) = FF
function Fout = FF(X, vX, E, B, t, cc)
Fout = cc*(E(X,t) + cross(vX,B(X,t)));