%EA III
%
%---------------------------------------------------------------------------
%This M-file,calculates the motion of an electron in an EM-field
% 	Force vector F=q{E+vXB}
% 	q...charge of particle
%	E...electric field vector
%	B...magnetic field vector
%	v...velocity of particle
%	X...vector cross-product
%---------------------------------------------------------------------------

clear all;
close all;			%% forget previous MATLAB runs									


disp('Motion of a charged particle in an EM-field');


%%Parameters
m	=	9.11*10^(-31);		%mass of electron
q	=	-1.6*10^(-19);		%charge of electron

%%Initial conditions:
k	=	1;					% time index
%% initial position; 
%% since we want to plot the paths of the particles we will
%% store the x, y and z positions in vectors rx, ry and rz respectively.
rx(k)	=	0; ry(k)	=	0; rz(k)	=	0;
%% Input initial velocity; we are not going to keep track of velocity for
%% all times, so we will use only one vector v to store vx, vy and vz
%% at current time.
v	=	[2.2e7	0 0];

%% Time parameters
start_time = 0;
dt = 0.5e-11;
end_time = 8e-9;
time(1) = start_time;	% let us keep a list of time values
w=2e9; 					%the omega that appears in parts d and e of Pset1/3

%% integrate for motion

for t = start_time:dt:end_time

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Describe the EM-field
	% Reason this is inside the loop is because this could be time varying

		E=[0 15e3 0]; B=[0 0 0];
							
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	

	rx(k+1) = rx(k) + v(1)*dt;	% calculate next x-position
	ry(k+1) = ry(k) + v(2)*dt;
	rz(k+1) = rz(k) + v(3)*dt;
	
	a = (q/m)*(E+cross(v,B));	% calculate acceleration and then next velocity
	v = v + a*dt;				% note we are replacing current v with next v
	
	time(k+1) = time(k) + dt;	% store time values for plotting
	k = k+1;
	
end



%%some fancy plotting code to better view motion
%% Determine if motion is 1D, 2D or 3D

motion_plane=0;	% default is 3D plot


if ( (max(rz)==min(rz))	& (min(ry)==max(ry)) )
	motion_plane = 1;	% motion is 1D along x-direction only
	disp('motion is 1D along x-direction only')
end
if ( (max(rz)==min(rz))	& (min(rx)==max(rx)) )
	motion_plane = 2;	% motion is 1D along y-direction only
		disp('motion is 1D along y-direction only')
end
if ( (max(rx)==min(rx))	& (min(ry)==max(ry)) )
	motion_plane = 3;	% motion is 1D along z-direction only
		disp('motion is 1D along z-direction only')
end
if ( (max(rz)==min(rz))	& (min(ry)~=max(ry)) & (min(rx)~=max(rx)) )
	motion_plane = 4;	% motion is 2D in xy-plane
		disp('motion is 2D in xy-plane only')
end
if ( (max(rx)==min(rx))	& (min(ry)~=max(ry)) & (min(rz)~=max(rz)) )
	motion_plane = 5;	% motion is 2D in yz-plane
			disp('motion is 2D in yz-plane only')
end
if (max(ry)==min(ry)	& min(rx)~=max(rx) & min(rz)~=max(rz) )
	motion_plane = 6;	% motion is 2D in xz-plane
			disp('motion is 2D in zx-plane only')
end
if ( max(rz)~=min(rz)	& min(ry)~=max(ry) & min(rx)~=max(rx) )
	motion_plane = 0;	% motion is 3D in xyz space
			disp('motion is 3D')
end

switch motion_plane

	case 1,
		{plot(time,rx)					% so plot rx vs time
		title('Plot of x as a function of time');
		xlabel('t (s) ');
		ylabel('x location ');
		}
	case 2,
		{plot(time,ry)					% so plot yx vs time
		title('Plot of y as a function of time');
		xlabel('t (s) ');
		ylabel('y location ');
		}
	case 3,
		{plot(time,rz)					% so plot rz vs time
		title('Plot of z as a function of time');
		xlabel('t (s) ');
		ylabel('z location ');
		}
	case 4
		comet(rx,ry)					% animated motion
		pause;
		plot(rx,ry);					% fixed plot
		title('Motion of an electron in an electromagnetic field');
		xlabel('x'); ylabel('y ');
		grid
	case 5
		comet(ry,rz)					% animated motion
		pause;
		plot(ry,rz);					% fixed plot
		title('Motion of an electron in an electromagnetic field');
		zlabel('z'); ylabel('y ');
		grid
	case 6
		comet(rx,rz)					% animated motion
		pause;
		plot(rx,rz);					% fixed plot
		title('Motion of an electron in an electromagnetic field');
		xlabel('x'); zlabel('z ');
		grid
	case 0
		comet3(rx,ry,rz);				% animated motion
		pause;
		plot3(rx,ry,rz);				%fiexd plot
		title('Motion of an electron in an electromagnetic field');
		xlabel('x'); ylabel('y'); zlabel('z');
		grid
end