% 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

function [] = porbital

hold off;clf; clear; 

% Input Electric and Magnetic Fields


% Input parameters and initial conditions

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

% Book keeping

E = @EE;
B = @BB;

cc = q/m;

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

Xarray = X';
T = [t];


h = .001;
M = 10;
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),3,'o')
title('Particle Orbit');
grid;
hold on;

% this next part has a circle trace out the trajectory
% specifying a bigger nstep will make the "particle" traverse faster.

nstep = 10;
comet3(Xarray(1:nstep:size(Xarray,1),1), Xarray(1:nstep:size(Xarray,1),2), Xarray(1:nstep:size(Xarray,1),3))

% Define electric and magnetic fields

function Bout = BB(X,t)
rho = sqrt(X(1)^2 + X(2)^2);
Bout = [0,0,10*sin(t)];

function Eout = EE(X,t)
r = sqrt(dot(X,X));
Eout = [(-100*X(1)/r^2),(-100*X(2)/r^2),(-100*X(3)/r^2)];

% 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)));


%Another fun pair of fields to use is
%Bout = [-30*X(1)/(rho^2), 30*X(2)/(rho^2), 0];
%Eout = [0,0,0];