# How to determine the minimum grid length | Numerical Plasma physics

I am trying to understand the attached python 3X code.

The following dispersion relation is given

$$\epsilon (k, \omega) = 1 - \frac 1 2 \left[ \frac{\omega_p^2}{(\omega-kv_0)^2} + \frac{\omega_p^2}{(\omega+kv_0)^2}\right]$$

What I am trying to understand is how to find the minimum grid length $$L_{min}$$ (as a function of $$\frac{v_0}{w}$$). $$L_{min}$$ indicates the needed minimum grid length to support such unstable modes.

To do so the code below should be checked in order to understand the plasma behavior either side of $$L_{min}$$ threshold, ie $$L_{min} < L$$ and $$L < L_{min}$$. To do so we should: 1) Adjust the number of simulation particles to grid points to improve the statistics, 2) Fix the number of particles per cell (i.e. npart/ngrid) and keep it $$>> 1$$, in order to reduce numerical noise and the runtime needed to observe the instability can be estimated from the maximum growth rate.

The theory is OK but I am struggling to actually find the minimum grid length $$L_{min}$$ once given the code. Might you please shed some light and guide me through? Thank you.

#! /usr/bin/python
#
#  Python script for computing and plotting single charged particle
#  trajectories in prescribed electric and magnetic fields.
#  Roughly equivalent to boris.m matlab program

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.widgets import Slider, Button, RadioButtons
from mpl_toolkits.mplot3d import Axes3D
import os
import os.path
import sys
from sys import exit
from time import sleep

# ===================================
#
# Function to integrate particle trajectory
# in given E, B fields
#
# ===================================

def integrate(E0, B0, vz0):
global dt, v0, x0, xp, yp, zp, qom, larmor, nsteps
wc=qom*B0 # cyclotron frequency
larmor=vperp/wc
print ("Cyclotron frequency =",wc)
print ("Perpendicular velocity v_p=",vperp)

norm = 1.  # choose whether to normalise plot axes dimensions to Larmor radius
trun=5*2*np.pi/wc  # total runtime
dt=.1/wc  # timestep - adjust to current B-field

nsteps=int(trun/dt)  # timesteps
E=np.array([0.,E0,0.])  # initial E-field
B=np.array([0.,0.,B0])  # initial B-field
u=np.array([0.,0.,0.])  # intermediate velocity
h=np.array([0.,0.,0.])  # normalized B-field
xp[0]=x0[0]
yp[0]=x0[1]
zp[0]=x0[2]
v0[2]=vz0 # z-component

v=v0+.5*dt*qom*(E+np.cross(v0,B)) # shift initial velocity back 1/2 step
x=x0

for itime in range(1,nsteps):
x=x+dt*v
xp[itime]=x[0] /norm
yp[itime]=x[1] /norm
zp[itime]=x[2] /norm
tp[itime]=itime*dt
#
# Boris mover: solves dv/dt = q/m*(E + vxB) to 2nd order accuracy in dt
#
qomdt2 = dt*qom/2
h = qomdt2*B
s=2*h/(1+np.dot(h,h))
u = v + qomdt2*E
up=u+np.cross(u+np.cross(u,h),s)
v=up+qomdt2*E

#     vxp[itime] = v[0]

# ===================================

# Make 2D plots of particle orbit
#
# ===================================

def plot_track2D():
global xp,yp,nsteps,ax1

fig = plt.figure(figsize=(8,8)) # initialize plot
xmin=np.min(xp)
xmax=np.max(xp)
ymin=np.min(yp)
ymax=np.max(yp)
fig.add_subplot(221) # 1st subplot in 2x2 arrangement
plt.cla()
plt.grid(True, which='both')
plt.xlim( (xmin, xmax) )
plt.ylim( (ymin, ymax) )
plt.xlabel('$$x$$')
plt.ylabel('$$y$$')
plt.plot(xp[0:nsteps],yp[0:nsteps],c='b')

plt.draw()
plt.savefig('./particle_orbit.png') # Save plot to file

# ===================================
#
# Make 3D plot of particle orbit
#
# ===================================

def plot_track3D():
global xp,yp,zp,nsteps,ax1
xmin=np.min(xp)
xmax=np.max(xp)
ymin=np.min(yp)
ymax=np.max(yp)
zmin=np.min(zp)
zmax=np.max(zp)
ax1.cla()

plt.ion()
plt.grid(True, which='both')
ax1.set_xlim( (xmin, xmax) )
ax1.set_ylim( (ymin, ymax) )
ax1.set_zlim( (zmin, zmax) )
ax1.set_xlabel('$$x$$ [m]')
ax1.set_ylabel('$$y$$ [m]')
ax1.set_zlabel('$$z$$ [m]')
#ax1.set_aspect(1.)
ax1.scatter(xp,yp,zp,c=tp,marker='o') # tracks coloured by elapsed time since start
plt.draw()

# =============================================
#
#  Main program
#
# =============================================

print ("Charged particle orbit solver")
plotboxsize   = 8.
animated = True

x0=np.array([0.,0.,0.])     # initial coords
vz0=0.
v0=np.array([-1e2,0.,vz0]) # initial velocity
vperp = np.sqrt(v0[0]**2+v0[2]**2)
E0=0.
B0=.1

e=1.602176e-19 # electron charge
m=9.109e-31 # electron mass
qom=e/m  # charge/mass ratio

wc=qom*B0 # cyclotron frequency
larmor=vperp/wc
print (wc,vperp,larmor)

trun=5*2*np.pi/wc  # total runtime
dt=.1/wc  # timestep - adjust to current B-field

nsteps=int(trun/dt)  # timesteps

#wc=qom*np.linalg.norm(B) # cyclotron frequency

#nsteps=2
tp = np.zeros(nsteps)  # variables to store particle tracks
xp = np.zeros(nsteps)
yp = np.zeros(nsteps)
zp = np.zeros(nsteps)
vxp = np.zeros(nsteps)
vyp = np.zeros(nsteps)
vzp = np.zeros(nsteps)

# Compute orbit
integrate(E0, B0, vz0)

# 2D orbit plotter
plot_track2D()

exit(0) # Quit script before 3D plot - comment out to continue!

# Start 3D interactive mode with sliders for B, E and v0

plt.ion() # Turn on interactive plot display
fig = plt.figure(figsize=(8,8))
# Get instance of Axis3D

# Get current rotation angle
print (ax1.azim)

# Set initial view to x-y plane
ax1.view_init(elev=90,azim=0)
ax1.set_xlabel('$$x$$[microns]')
ax1.set_ylabel('$$y$$[microns]')
ax1.set_zlabel('$$z$$[microns]')
plot_track3D()

#filename = 'a0_45/parts_p0000.%0*d'%(6, ts)
#plot_from_file(filename):
axcolor = 'lightgoldenrodyellow'
axe0 = fig.add_axes([0.1, 0.95, 0.3, 0.03])#, facecolor=axcolor) # box position, color & size
axb0  = fig.add_axes([0.5, 0.95, 0.3, 0.03])#, facecolor=axcolor)
axv0  = fig.add_axes([0.1, 0.9, 0.3, 0.03])#, facecolor=axcolor)

sefield = Slider(axe0, 'Ey [V/m]', -5.0,5.0, valinit=E0)
sbfield = Slider(axb0, 'Bz [T]', -1.0, 1.0, valinit=B0)
svz = Slider(axv0, 'vz [m/s]', 0.0, 1.0, valinit=0.)

def update(val):
E0 = sefield.val
B0 = sbfield.val
vz0 = svz.val

integrate(E0,B0,vz0)
plot_track3D()
plt.draw()

sefield.on_changed(update)
sbfield.on_changed(update)
svz.on_changed(update)

resetax = fig.add_axes([0.8, 0.025, 0.1, 0.04])
button = Button(resetax, 'Reset', color=axcolor, hovercolor='0.975')
def reset(event):
global ax1
sefield.reset()
sbfield.reset()
svz.reset()
ax1.cla()
ax1.set_xlabel('$$x$$[microns]')
ax1.set_ylabel('$$y$$[microns]')
ax1.set_xlim( (0., 10.) )
#    ax1.set_ylim( (-sigma, sigma) )
ax1.grid(True, which='both')
plt.draw()
button.on_clicked(reset)

#plt.show()
plt.show(block=False)
• Unless you bump into someone who is an expert on plasma physics simulations it is very hard to help you. It is not clear what numerical method you use to simulate what particular equation / model. Furthermore, providing a bunch of code without any comments or explanations in the text is also not very helpful. Try to link your questions to concepts the community is more familiar with. Jan 2 at 18:16
• @DanDoe you are absolutely right, I am thinking how to edit my question in the most efficient way to make it more to the point. Jan 2 at 21:05