# Plasma Model

### Plasma Model (Under Construction 2/7/2012)

This is currently a work in progress.  The goal is to monitor the position of charged particles in a box (meaning the particles must not escape) as a function of time.  A plasma is ionized, overall electrically neutral gas. It is composed of ionized atoms and free electrons.  The problem I'm encountering now is that since the interaction forces are weak, the particles don't interact too much.  I need to add more particles.  This however makes computation time insanely long on my computer.  I will continue working on this a bit every day and after testing a simple 1-D case, will post the final results.  Below is the code as it stands now (generic 3-D model). There are no comments in it right now (to be added later).  The computation is long because the computer has to sum forces of each of the N particles and then advance the time to do it again. There are 3*(N-1) interactions on each particle.  The 3 is due to 3 dimensions and there are N-1 interactions per dimension.  I have to compute the Coulomb force on each particle due to all of the other particles (not counting itself).
clc
clear
format 'long'
X1 = [];
dt = 10^-15;
t = 100*dt;
N = 500;
T0 = 20;
mp = 1.66*10^-27;
me = 9.11*10^-31;
Lx = 10^-9;
Ly = 10^-9;
Lz = 10^-9;
e = -1.6*10^-19;
p = 1.6*10^-19;
k = 8.9876*10^9;
kb = 1.38*10^-23;
xp = [];
yp = [];
zp = [];
vpx = [];
vpy = [];
vpz = [];
xe = [];
ye = [];
ze = [];
vex = [];
vey = [];
vez = [];
for i = 1:N
xp = [xp,Lx*rand(1)];
yp = [yp,Ly*rand(1)];
zp = [zp,Lz*rand(1)];
vpx = [vpx,(-1)^i*sqrt(kb*T0/mp)];
vpy = [vpy,(-1)^i*sqrt(kb*T0/mp)];
vpz = [vpz,(-1)^i*sqrt(kb*T0/mp)];
xe = [xe,Lx*rand(1)];
ye = [ye,Ly*rand(1)];
ze = [ze,Lz*rand(1)];
vex = [vex,(-1)^i*sqrt(kb*T0/me)];
vey = [vey,(-1)^i*sqrt(kb*T0/me)];
vez = [vez,(-1)^i*sqrt(kb*T0/me)];
end
for tau = 1:t/dt
Fex = [];
NFex = [];
Fey = [];
NFey = [];
Fez = [];
NFez = [];
Fpx = [];
NFpx = [];
Fpy = [];
NFpy = [];
Fpz = [];
NFpz = [];
X1 = [X1,xe(1)];
for j = 1:N
for k = 1:N
if j ~= k
Ree = sqrt((xe(j)-xe(k))^2+(ye(j)-ye(k))^2+(ze(j)-ze(k))^2);
Rep = sqrt((xe(j)-xp(k))^2+(ye(j)-yp(k))^2+(ze(j)-zp(k))^2);
Fex = [Fex,k*e*e*(xe(j)-xe(k))/Ree^3 + k*e*p*(xe(j)-xp(k))/Rep^3];
Fey = [Fey,k*e*e*(ye(j)-ye(k))/Ree^3 + k*e*p*(ye(j)-yp(k))/Rep^3];
Fez = [Fez,k*e*e*(ze(j)-ze(k))/Ree^3 + k*e*p*(ze(j)-zp(k))/Rep^3];
else
Rep = sqrt((xe(j)-xp(k))^2+(ye(j)-yp(k))^2+(ze(j)-zp(k))^2);
Fex = [Fex,k*e*p*(xe(j)-xp(k))/Rep^3];
Fey = [Fey,k*e*p*(ye(j)-yp(k))/Rep^3];
Fez = [Fez,k*e*p*(ze(j)-zp(k))/Rep^3];
end
if k == N
Fex = sum(Fex);
Fey = sum(Fey);
Fez = sum(Fez);
end
end
NFex = [NFex;Fex];
Fex = [];
NFey = [NFey;Fey];
Fey = [];
NFez = [NFez;Fez];
Fez = [];
end
for j = 1:N
for k = 1:N
if j ~= k
Rpp = sqrt((xp(j)-xp(k))^2+(yp(j)-yp(k))^2+(zp(j)-zp(k))^2);
Rpe = sqrt((xp(j)-xe(k))^2+(yp(j)-ye(k))^2+(zp(j)-ze(k))^2);
Fpx = [Fpx,k*p*p*(xp(j)-xp(k))/Rpp^3 + k*p*e*(xp(j)-xe(k))/Rpe^3];
Fpy = [Fpy,k*p*p*(yp(j)-yp(k))/Rpp^3 + k*p*e*(yp(j)-ye(k))/Rpe^3];
Fpz = [Fpz,k*p*p*(zp(j)-zp(k))/Rpp^3 + k*p*e*(zp(j)-ze(k))/Rpe^3];
else
Rpe = sqrt((xp(j)-xe(k))^2+(yp(j)-ye(k))^2+(zp(j)-ze(k))^2);
Fpx = [Fpx,k*p*e*(xp(j)-xe(k))/Rpe^3];
Fpy = [Fpy,k*p*e*(yp(j)-ye(k))/Rpe^3];
Fpz = [Fpz,k*p*e*(zp(j)-ze(k))/Rpe^3];
end
if k == N
Fpx = sum(Fpx);
Fpy = sum(Fpy);
Fpz = sum(Fpz);
end
end
NFpx = [NFpx;Fpx];
Fpx = [];
NFpy = [NFpy;Fpy];
Fpy = [];
NFpz = [NFpz;Fpz];
Fpz = [];
end
for i = 1:N
xe(i) = (NFex(i)/me)*dt^2 + vex(i)*dt + xe(i);
vex(i) = (NFex(i)/me)*dt + vex(i);
xp(i) = (NFpx(i)/mp)*dt^2 + vpx(i)*dt + xp(i);
vpx(i) = (NFpx(i)/mp)*dt + vpx(i);
ye(i) = (NFey(i)/me)*dt^2 + vey(i)*dt + ye(i);
vey(i) = (NFey(i)/me)*dt + vey(i);
yp(i) = (NFpy(i)/mp)*dt^2 + vpy(i)*dt + yp(i);
vpy(i) = (NFpy(i)/mp)*dt + vpy(i);
ze(i) = (NFez(i)/me)*dt^2 + vez(i)*dt + ze(i);
vez(i) = (NFez(i)/me)*dt + vez(i);
zp(i) = (NFpz(i)/mp)*dt^2 + vpz(i)*dt + zp(i);
vpz(i) = (NFpz(i)/mp)*dt + vpz(i);
if xe(i) <= 0 || xe(i) >= Lx
vex(i) = -vex(i);
end
if ye(i) <= 0 || ye(i) >= Ly
vey(i) = -vey(i);
end
if ze(i) <= 0 || ze(i) >= Lz
vez(i) = -vez(i);
end
if xp(i) <= 0 || xp(i) >= Lx
vpx(i) = -vpx(i);
end
if yp(i) <= 0 || yp(i) >= Ly
vpy(i) = -vpy(i);
end
if zp(i) <= 0  || zp(i) >= Lz
vpz(i) = -vpz(i);
end
end
disp(tau)
end

plot(1:tau,X1)
disp(max(vex))
rho = N*(me + mp)/(Lx*Ly*Lz)