% Simulation to study Child Langmuir current density limit
% For more see:
% http://www.particleincell.com/blog/2014/current-density

more off;			%turn off pagination in Octave

%physical constants
AMU = 1.6605e-27;	%atomic mass in kg
E = 1.6022e-19;		%elementary charge in Coulomb
EPS0 = 8.8542e-12;	%permittivity of free space F/m
ME = 9.1094e-31;	%electron mass

%set boundaries
phi_left = 0;
phi_right = -100;
x_left = 0;
x_right = 0.01;

%compute mesh parameters
nx = 100;							%number of mesh nodes
delta_x = (x_right-x_left)/(nx-1);	%cell spacing, there are nx-1 cells

%set matrix data
A = zeros(nx);
A(1,1) = 1;		%fixed potential at left
A(nx,nx) = 1;		%fixed potential at right
for i=2:nx-1
	%finite difference stencil for second derivative
	A(i,i-1) = 1;
	A(i,i) = -2;
	A(i,i+1) = 1;
end


%ion data
ion_mass = 131.3*AMU;
ion_charge = 1*E;

%time data
delta_t=0.5*delta_x/sqrt(2*ion_charge/ion_mass*(abs(phi_right-phi_left))) %time step
it_max = 3000;		%number of time steps


%injection current
j_CL = (4/9)*EPS0*sqrt(2*ion_charge/ion_mass)*(abs(phi_right-phi_left))^1.5/(x_right-x_left)^2;
I = 1.0*j_CL		%this is also current density due to unit length in y and z
np_inject = 10;	%number of particles to inject per time step
ion_spwt = I*delta_t/(ion_charge*np_inject);	%number of actual particles per macroparticle

%initialize
part_count = 0;		%number of particles
efield = zeros(nx,1);

%main loop
for it=1:it_max

	%*** 1) generate new particles ***    
	for p=1:np_inject		%inject particles
		%increment particle counter
		part_count=part_count+1;	
		
		part_x(part_count) = 0;	%initial position
		part_v(part_count) = 0;	%initial velocity
		part_new(part_count) = 1;	%mark new particles
		
		%compute electric field at particle position
		fi = (part_x(part_count)-x_left)/delta_x + 1;
		i = floor(fi);
		di = fi-i;
		part_ef = efield(i)*(1-di)+efield(i+1)*(di);	
		
		%rewind the velocity back by 0.5*delta_t, this part of the leapfrog scheme
		part_v(part_count) = part_v(part_count) - (ion_charge/ion_mass)*part_ef*0.5*delta_t; 
	end
	
	%*** 2) compute field quantities***
	u = zeros(nx,1);
	count = zeros(nx,1);		%this field holds number of macroparticles
	for p=1:part_count
		%get particle mesh index
		fi = (part_x(p)-x_left)/delta_x + 1;
		i = floor(fi);
		di = fi-i;
		count(i) = count(i)+(1-di);
		count(i+1) = count(i+1)+ (di);
		
		u(i) = u(i)+(1-di)*part_v(p);
		u(i+1) = u(i+1)+ (di)*part_v(p);
	end
	den = ion_spwt*count/delta_x;	%divide by "volume"
	den(1) = den(1)* 2.0;			%only half area at boundaries
	den(nx) = den(nx)*2.0;
	
	%average velocity, make sure we don't divide by zero
	u(count~=0) = u(count~=0)./count(count~=0);
	u(count==0) = 0;
	j = den.*u*ion_charge;	%current density
	
	%*** 3) set charge density
	rho = den*ion_charge;		
	
	%*** 4) solve potential ***
	b = -rho/EPS0*delta_x*delta_x;
	b(1) = phi_left;
	b(nx) = phi_right;
	phi = A\b;
	 
	%*** 5) update electric field ***
	efield = -gradient(phi,delta_x);
	
	%*** 6) move particles ***
	for p=1:part_count
		%obtain electric field at the particle position
		fi = (part_x(p)-x_left)/delta_x + 1;
		i = floor(fi);
		di = fi-i;
		part_ef = efield(i)*(1-di)+efield(i+1)*(di);	
		
		% set random delta_t for newly born particles
		if (part_new(p))
			dt = (1.0+1.0*rand())*delta_t;
			part_new(p) = 0;	%clear new particle flag			
		else
			dt = delta_t;
		end
		
		%update particle velocity from F=m(dv/delta_t)=qE
		part_v(p) = part_v(p)+(ion_charge/ion_mass)*part_ef*dt;
		
		%update particle position
		part_x(p) = part_x(p)+part_v(p)*dt;
		
		%remove particles that crossed boundaries
		if (part_x(p)<x_left || part_x(p)>=x_right)
			%remove particle by replacing the current array data with the last
			part_x(p) = part_x(part_count);
			part_v(p) = part_v(part_count);
			part_new(p) = part_new(part_count);
			part_count=part_count-1;		%decrement particle counter
			p=p-1;		%also decrement current index to rescan this position
		end		
	end
	
	%*** 7) output some diagnostics
	
	%plot results
	if (mod(it,25)==0 || it==it_max)
	disp(sprintf('it = %d\t j_ave = %4.3g',it,mean(j(j~=0))));
		x_pos = [1:nx]*delta_x+x_left;
 		figure(1);
 		plot(x_pos,phi);
 		title('Potential');
 		figure(2);
 		plot(x_pos,efield);
 		title('Electric Field');
 		figure(3);
 		plot(x_pos,j,'-*');
 		title('Current Density');
 		axis([x_left,x_right, 0, 1.2*I]);
 		figure(4);
		plot(x_pos,u,'-*');
		title('Velocity');
		axis([x_left,x_right, -2000, 12000]);
 		
		%for saving animation
		%filename=sprintf('v%05d.png',it);
		%print(filename);	
		drawnow;
	end

end



