due 10/03/2018 at 11:59pm
binary_Coulomb_collisions
that can handle any arbitrary number of particles.plot_2D
function so that you can print the results.Let's setup the main parameters and objects for this homework.
import math
from math import sqrt as sqrt
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
class charged_particle:
m=9.10938356e-31 #charge
q=-1.60217662e-19 #mass
x=np.zeros(3) #location
v=np.zeros(3) #speed
We use here an array of object to pass all the particles to the function. Much faster and more general than the function binary_Coulomb_collision
in the text of Ch03. Note that we added a cut-off radius r_cutoff
to limit large acceleration between time steps. This allows to increase the time steps between iterations but it also increase the error on the trajectory.
def Coulomb_collisions(charged_particles,time_frame,n_steps,r_cutoff):
N=len(charged_particles)
loc=np.zeros((3,n_steps,N))
dt=(time_frame[1]-time_frame[0])/n_steps
for i in range(n_steps):
for j in range(N):
E=np.zeros(3)
for k in range(N):
R=charged_particles[j].x-charged_particles[k].x
r=np.linalg.norm(R)
if (r>r_cutoff):
E+=charged_particles[k].q/r**2*R/r
charged_particles[j].v+=1/(4*np.pi*8.85e-12)*E*charged_particles[j].q/charged_particles[j].m*dt
for j in range(N):
charged_particles[j].x+=charged_particles[j].v*dt # we compute all the velocities then we compute the locations
loc[:,i,j]=charged_particles[j].x
return loc
That's plot_2D
from the text.
def plot_2D(r,n_interpolated=1,plot_size=3,plot_axes='axes'):
n = len(r[0,:,0])
k = n_interpolated
fig, ax = plt.subplots(1, 1, figsize=(plot_size, plot_size))
# Now, we draw our points with a gradient of colors.
for i in range(len(r[0,0,:])):
x_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, r[0,:,i])
y_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, r[1,:,i])
ax.scatter(x_interpolated, y_interpolated, c=range(n*k), linewidths=0,
marker='o', s=2*plot_size, cmap=plt.cm.jet,)
#compute autoscale parameters
xc=(r[0,:,:].max()+r[0,:,:].min())/2.
x_low=xc-(r[0,:,:].max()-r[0,:,:].min())/2.*1.1
x_high=xc+(r[0,:,:].max()-r[0,:,:].min())/2.*1.1
yc=(r[1,:,:].max()+r[1,:,:].min())/2.
y_low=yc-(r[1,:,:].max()-r[1,:,:].min())/2.*1.1
y_high=yc+(r[1,:,:].max()-r[1,:,:].min())/2.*1.1
#set autoscale parameters
ax.set_xlim(x_low,x_high)
ax.set_ylim(y_low,y_high)
if (plot_axes!="axes"):
ax.set_axis_off()
else:
plt.xlabel("x (m)")
plt.ylabel("y (m)")
We now set the initial conditions
N=25
cp=np.empty(N,dtype=object)
for i in range (N):
cp[i]=charged_particle()
cp[i].x=np.array([(np.random.random_sample()-0.5)*30e-3,(np.random.random_sample()-0.5)*30e-3,0])
cp[i].v=np.array([100*(np.random.random_sample()-0.5),100*(np.random.random_sample()-0.5),0])
time_frame=[0,4e-5]
n_steps=200
And then we run the code and plot the results.
r=Coulomb_collisions(cp,time_frame,n_steps,1e-10)
plot_2D(r,n_interpolated=1,plot_size=8)
It worked!
This is a difficult problem to adjust to get the right oscillations. The particles have to be far apart so the electrostatic forces are weak enough to allow for large time steps. We need an acceptable cut-off radius when computing the electrostatic force so the particles do not acquire large velocities between two consecutive time steps. Note that this is really to allow for fast computations rather than enabling oscillations to occur. Here we use $n_0=10^{13}m^{-3}$ which gives an oscillation period of 40ns, a time scale that has worked for us in previous problems. We also try to keep the size of the lattice at a scale that we have also used successfully. We used a length of $100 \mu m$. We take a radius cut-off on the order of $1\mu m$, a thousand smaller than the scale length of our lattice.
N=3
length=1e-4
n0=N**2/length**3
omega=sqrt(n0*1.6e-19**2/(8.85e-12*9e-31))
f=omega/(2*3.1416)
T=1/f
print('Density',n0,'m^-3 Frequency',f,'Hz Period',T,'s')
Density 8999999999999.998 m^-3 Frequency 27068704.021020308 Hz Period 3.694303204259229e-08 s
Now, starting small, we slowly increase the velocity so that oscillations start to appear. It is important to realize that these oscillations correspond to motion between two consecutive ions, not motion around one single ion. Check the difference between the velocity where oscillations between two ions occur ($v_x=6km/s$) and revolutions around a single ion ($v_x=4km/s$), or complete escape ($v_x=8 km/s$). It is important to note that complete escape would not happen if our lattice was repeating to infinity in all directions.
With these parameters we get particles that start at a given position $x$ and get back top that position in 60ns, which is close to the plasma frequency. We do not expect to find the exact period due to the limitation in lattice size.
delta_x=length/10
delta_y=length/10
lattice=np.linspace(-length,length,N)
lattice=np.linspace(-length,length,N)
cp=np.empty(0,dtype=object)
l=0
for i in range (N):
for j in range (N):
cp=np.append(cp,charged_particle())
cp[l].x=np.array([lattice[i],lattice[j],0])
cp[l].m=1.67e-27
cp[l].q*=-1
cp[l].v=np.array([0.,0,0])
l+=1
for i in range (N):
for j in range (N):
cp=np.append(cp,charged_particle())
cp[l].x=np.array([lattice[i]+delta_x,lattice[j]+delta_y,0])
cp[l].v=np.array([6000.,0,0])
l+=1
time_frame=[0,6e-8]
n_steps=200
r=Coulomb_collisions(cp,time_frame,n_steps,3e-6)
plot_2D(r,n_interpolated=3,plot_size=8)