This chapter deals with charge particle motion.

Particles experience three different types of interactions

- collisions (direct or indirect): due to direct interaction between particles.
- grativity: when particles have mass
- electromagnetic forces: when particles have charge

In this chapter, we ignore collisions. They will be treated in a different chapter.

Gravitation is always present in any experiment in a laboratory or in astrophysical settings. However, it can be ignored when charged particles are surrounded by large electric fields or magnetic fields. However, one should always estimate the impact of any force to assess which are negligible and which is not. Ultimately, developing a model where all forces are taken into account will lead to unwanted complexity.

Overall, the best approach is *back of the envelope* calculation before taking on big computational problems and apply the full force of numerical analysis to generate large 3-dimensional datasets.

When the change in gravitation is small compared to the particle trajectory, the force applied on a particle of mass $m$ is simply its weight. $$\vec{F}=m\vec{g}\tag{2.1}$$ We now need to define the trajectory as a function of time using Newton's second law $$\frac{d\vec{p}}{dt}=\sum_i{\vec{F_i}}$$ $\vec{p}$ is the particle momentum given by $\vec{p}=m\vec{v}$ using the particle mass and it velocity $\vec{v}$. Clearly, we have only one force here! So, we shall simplify it right away to $$ \frac{d\vec{p}}{dt}=m\vec{g}\tag{2.2} $$ It is possible to integrate Eq. $(2.2)$ using a very simple integration method for Ordinary Differential Equations.

Newton's second law, i.e. Eq. $(2.1)$, can be discretized. Discretization means in this context that $dt$, representing an infinitely small amount of time, is turned into $\Delta t$, a finite amount of time. We can rewrite Newton's law as
$$\frac{\Delta p}{\Delta t}=\vec{F}\tag{2.3}.$$
Clearly, we get from this $p_{t+\Delta t}-p_t=\vec{F}\Delta t$. If we know the momentum and the force at $t$, then we can compute the new momentum at $\Delta t$. To extract the trajectory from this equation, we have to realize that $\vec{v}=\frac{d\vec{r}}{dt}$ where $\vec{r}(t)=(x(t),y(t),z(t))$, the coordinate of the particle at the time $t$. We can rewrite Eq. $(2.3)$ as $$m\vec{v}_{t+\Delta t}-m\vec{v}_t=m\vec{g}\Delta t$$ or $$\vec{v}_{t+\Delta t}-\vec v_t=\vec g\Delta t$$
It seems that **this trajectory does not dependent on the particle mass!** This is a remarkable result. An electron inside a gravitational field falls as fast as a proton, which is 2000 times more massive. We need to expand this equation to our coordinate system $(x,y,z)$

We are going to solve equation $(2.1)$ using the following conditions. First, $$\vec{g}=(0,g_y,0)$$ with $g_y=-9.81m/s^2$. Our initial conditions are $$v_x=0\\v_y=0\\v_z=0$$ and $$x=0\\y=0\\z=0$$ Clearly, nothing happens in the z direction. So, we will ignore it for now.

Let's start by importing our tools into the python code, math, numpy and matplotlib.

In [1]:

```
import math
from math import sqrt as sqrt
import numpy as np
import matplotlib.pyplot as plt
#to display inside the notebook!
%matplotlib inline
```

Then we define our main (i.e. global) variables, all in mks

In [2]:

```
#initial time
t_initial=0
#final time
t_final=10
#time step
dt=0.5
#number of time steps declared as integer
n = int((t_final-t_initial)/dt)
#electron mass
m=9.10938356e-31
#gravity at sea level
gx=0
gy=-9.81
#initial velocity
vx_initial=0
vy_initial=0
#we use the momemtum to be physically sound, even if using velocity is more computationally effecient
px_initial=m*vx_initial
py_initial=m*vy_initial
#initial position
x_initial=0
y_initial=0
```

**Avoid using np.empty(...) as this function allocates memory but does not set it to a given value (like 0s).**

In [3]:

```
x=np.zeros(1)
y=np.zeros(1)
```

In [4]:

```
px_old=px_initial
py_old=py_initial
x_old=x_initial
y_old=y_initial
```

`for`

loop.
`for i in range(1,n,1):`

This loop goes from `1`

to `n`

with a step of `1`

. Note that there is a column `:`

at the end that indicates to python that the code included in the loop is below. Indentations are used to tell python which piece of code is in the loop. So **indentations are important in python**. To put a piece of code outside of the loop, you need to indent back to the left.

In [5]:

```
for i in range(1,n,1):
fx=m*gx #this is the force along x
fy=m*gy #this is the force along y
px_new=fx*dt+px_old #new momentum
py_new=fy*dt+py_old
vx=px_new/m #new velocity
vy=py_new/m
x_new=vx*dt+x_old #new position
y_new=vy*dt+y_old
x=np.append(x,[x_new]) #the positions are appended at the end of the x and y arrays.
y=np.append(y,[y_new])
x_old=x_new # we now replace the old data with the new ones
y_old=y_new
px_old=px_new
py_old=py_new
# and the loop starts again.
```

`i`

is the step counter going from `1`

to `n`

, with s step of `1`

. Note the indent to indicate python that this piece of code is inside the `for`

loop. Note that we note only compute the components of the force $\vec{F}$ `fx`

and `fy`

but both variables are also declared at the same time. This goes for the other variables also. It is time to plot the results now

In [6]:

```
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
# Now, we draw our points with a gradient of colors.
ax.scatter(x, y, c=range(n), linewidths=0,marker='o', s=15, cmap=plt.cm.jet,)
ax.axis('equal')
ax.set_axis_off()
```

`x`

and `y`

.

In [7]:

```
k = 100
x_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, x)
y_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, y)
```

And see what we get

In [8]:

```
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
# Now, we draw our points with a gradient of colors.
ax.scatter(x_interpolated, y_interpolated, c=range(n*k), linewidths=0,
marker='o', s=15, cmap=plt.cm.jet,)
ax.axis('equal')
ax.set_axis_off()
```

We can make a simple plotting function by rewriting the code above slightly.

In [9]:

```
def plot_2D(x,y,n_interpolated,plot_size,plot_axes):
n = len(x)
k = n_interpolated
x_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, x)
y_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, y)
fig, ax = plt.subplots(1, 1, figsize=(plot_size, plot_size))
# Now, we draw our points with a gradient of colors.
ax.scatter(x_interpolated, y_interpolated, c=range(n*k), linewidths=0,
marker='o', s=2*plot_size, cmap=plt.cm.jet,) #s is the size of the plotted symbol 'o'
ax.axis('equal')
if (plot_axes!="axes"):#so we can switch the axis on or off
ax.set_axis_off()
else:
plt.xlabel("x (m)")
plt.ylabel("y (m)")
```

And we can print easily now.

In [10]:

```
plot_2D(x_interpolated,y_interpolated,10,4,"no")
plot_2D(x_interpolated,y_interpolated,10,4,"axes")
```

Let's take a break from physics and see if we can use a function to integrate Eq. $(2.3)$. We can reuse the `for`

loop above to generate an ordinary differential equation integrator using $x_0$, $y_0$, $v_{x,0}$, $v_{y,0}$, $t_{initial}$, $t_{final}$ and the number of steps $n_{steps}$ as inputs. We define it as

```
def OED_integration(x_initial,y_initial,vx_initial,vy_initial,t_initial,t_final,n_steps):
```

The other variables, like `m`

, `t`

, `dt`

, ... are kept as is and do not have to be part of the input of OED_integration.

In [11]:

```
def ODE_integration(x_initial,y_initial,vx_initial,vy_initial,t_initial,t_final,n_steps):
x=np.full(1,x_initial) # we do this to initialize the x variable as a numpy array
y=np.full(1,y_initial)
px_old=vx_initial*m
py_old=vy_initial*m
dt=(t_final-t_initial)/n_steps
for i in range(1,n_steps,1):
fx=m*gx #this is the force along x
fy=m*gy #this is the force along y
px_new=fx*dt+px_old #new momentum
py_new=fy*dt+py_old
vx=px_new/m #new velocity
vy=py_new/m
x_new=vx*dt+x[i-1] #new position
y_new=vy*dt+y[i-1] # in python arrays start at index 0
x=np.append(x,x_new) #the positions are appended at the end of the x and y arrays.
y=np.append(y,y_new)
"""The ones below are not needed. We use the ones from the previous step inside the array
x_old=x_new
y_old=y_new"""
px_old=px_new
py_old=py_new
return x,y
# and the loop starts again.
```

`ODE_integration(x_initial,y_initial,vx_initial,vy_initial,t_initial,t_final,n_steps)`

In [12]:

```
x,y=ODE_integration(0,0,10,10,0,3,80)
plot_2D(x,y,10,8,"axes")
```

It is now time to move on to electrostatic forces. A particle with charge $q$ placed inside an electric field $\vec{E}$ is submitted to a force $$\vec F=q\vec E\tag{2.4}$$ Does this formula remind you of another one seen before?

Indeed, the particle feels a force that looks like gravity (i.e. Eq. $2.1$). However, there is a subtle difference when we use Newton's second law $$ m\frac{d\vec v}{dt}=q\vec E$$ This time around the mass does not go away and we expect the trajectory to depend on the mass. Here is the new system of equations we need to solve for. $$v_{x,t+\Delta t}=\frac{q}{m}E_x\Delta t+v_{x,t}\\ v_{y,t+\Delta t}=\frac{q}{m}E_y\Delta t+v_{y,t}\\ v_{z,t+\Delta t}=\frac{q}{m}E_z\Delta t+v_{z,t}$$ And since $\vec{v}=\frac{d\vec{r}}{dt}$ we get the same system as before for the coordinates, i.e. $$x_{t+\Delta t}=v_{x,t+\Delta t}\Delta t+x_t\\ y_{t+\Delta t}=v_{y,t+\Delta t}\Delta t+y_t\\ z_{t+\Delta t}=v_{z,t+\Delta t}\Delta t+z_t$$

We are going to solve equation $(2.1)$ using the following conditions. First, $$\vec E=(0,E_y,0)$$ with $E_y=-1 V/m$. Our initial conditions are $$v_x=0\\v_y=0\\v_z=0$$ and $$x=0\\y=0\\z=0$$ Clearly, nothing happens in the z direction. So, we will ignore it for now.

In [13]:

```
q=-1.60217662e-19
Ex=0
Ey=1
vx_final=0
vy_final=0
clight=299792458
```

And then we make our new function, named `OED_integration_E`

.

In [14]:

```
def ODE_integration_E(x_initial,y_initial,vx_initial,vy_initial,t_initial,t_final,n_steps):
x=np.full(1,x_initial) # we do this to initialize the x variable as a numpy array
y=np.full(1,y_initial)
px_old=vx_initial*m
py_old=vy_initial*m
dt=(t_final-t_initial)/n_steps
for i in range(1,n_steps,1):
fx=q*Ex #this is the force along x
fy=q*Ey #this is the force along y
px_new=fx*dt+px_old #new momentum
py_new=fy*dt+py_old
vx=px_new/m #new velocity
vy=py_new/m
x_new=vx*dt+x[i-1] #new position
y_new=vy*dt+y[i-1] # in python arrays start at index 0
x=np.append(x,x_new) #the positions are appended at the end of the x and y arrays using numpy.
y=np.append(y,y_new)
px_old=px_new
py_old=py_new
print(vx/clight,vy/clight)
print(x[n_steps-1],"m",y[n_steps-1],"m")
return x,y
# and the loop starts again.
```

Can you guess what changed??

Besides the name (which we changed to avoid overriding the previous function), we changed only two lines

```
fx=q*Ex #this is the force along x
fy=q*Ey #this is the force along y
```

We could do this because the previous code conversed physical information and separated force and momentum. While being less efficient computationally, we were able to reuse the code, mostly unchanged.

So let's run it!

In [15]:

```
x,y=ODE_integration_E(0,0,0,0,0,3,2000)
plot_2D(x,y,1,8,"axes")
```

0.0 -1759.1576096745716 0.0 m -791073275720.6312 m

In [16]:

```
print(q*Ey/(2.*m)*3.**2,"m")
```

-791469010225.7589 m

In [17]:

```
x,y=ODE_integration_E(0,0,0,0,0,3e-6,20)
plot_2D(x,y,1,8,"axes")
```

0.0 -0.0016720357470644369 0.0 m -0.7518955597144708 m

In [18]:

```
x,y=ODE_integration_E(0,0,2e5,2e5,0,5e-6,60)
plot_2D(x,y,1,8,"axes")
```

0.0006671281903963041 -0.002217377922960473 0.9833333333333347 m -1.1785496112648044 m

**Figure 1: The geometry used to compute the electric field generated by a charge q at a remote location (the cross).**

In [19]:

```
def ODE_integration_E_with_charge(charge_location,charge_value,initial_position,initial_velocity,time_frame,n_steps):
loc=np.zeros((3,n_steps)) #index 0 is x, index 1 is y and index 2 is z
loc[:,0]=initial_position
""" Let's start to use vectors and combine notations
the column : is equivalent to writing
loc[0,0]=initial_position[0]
loc[1,0]=initial_position[1]
loc[2,0]=initial_position[2]
"""
v_old=initial_velocity
dt=(time_frame[1]-time_frame[0])/n_steps
E=[0,0,0]
f=np.zeros(3)
v_new=np.zeros(3)
for i in range(1,n_steps,1):
r=sqrt((loc[0,i-1]-charge_location[0])**2+(loc[1,i-1]-charge_location[1])**2+(loc[2,i-1]-charge_location[2])**2)
if (r<1e-7):
r=1e-7 # we won't divide by 0
r_unit=(loc[:,i-1]-charge_location)/r
E=1./(4*3.1416*8.85e-12)*charge_value/r**2*r_unit
f=q*E
v_new=f/m*dt+v_old #new velocity
loc[:,i]=v_new*dt+loc[:,i-1] #new position
v_old=v_new
return loc
```

Now let's define a new plotting function that indicates where the motionless charge is located

In [20]:

```
def plot_2D_with_charge(rc,x,y,n_interpolated,plot_size,plot_axes):
n = len(x)
k = n_interpolated
x_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, x)
y_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, y)
fig, ax = plt.subplots(1, 1, figsize=(plot_size, plot_size))
# Now, we draw our points with a gradient of colors.
ax.scatter(x_interpolated, y_interpolated, c=range(n*k), linewidths=0,
marker='o', s=2*plot_size, cmap=plt.cm.jet,) #s is the size of the plotted symbol 'o'
ax.scatter(rc[0], rc[1], c='black', linewidths=0,
marker='x', s=3*plot_size,)
ax.axis('equal')
if (plot_axes!="axes"):#so we can switch the axis on or off
ax.set_axis_off()
else:
plt.xlabel("x (m)")
plt.ylabel("y (m)")
```

In [21]:

```
rc=[4e-3,0,0] #charge location
qc=1.6e-19 #charge value
loc=ODE_integration_E_with_charge(rc,qc,(-5e-3,10e-4,0),(500,0,0),(0,50e-6),1000)
plot_2D_with_charge(rc,loc[0,:],loc[1,:],1,7,"axes")
```

After looking at electrostatic forces, we now study magnetic forces. A particle of charge $q$ and traveling at a velocity $\vec v$ inside a magnetic field $\vec{B}$ is given by $$\vec F=q\vec v\times \vec B$$ In this case, the equation of motion is given by $$m\frac{d\vec v}{dt}=q\vec v\times \vec B\tag{2.7}$$

We now have a more complex behavior since the system $$\frac{dv_x}{dt}=\frac{q}{m}(v_yB_z-v_zB_y)\\ \frac{dv_y}{dt}=\frac{q}{m}(v_zB_x-v_xB_z)\\ \frac{dv_z}{dt}=\frac{q}{m}(v_xB_y-v_yB_x)$$ formed a set of coupled ODEs.

In [22]:

```
r_L=1e3*9e-31/(1.6e-19*1)
print(r_L,"m")
```

5.6250000000000016e-09 m

Now let's integrate our ODE numerically

In [23]:

```
m=9e-31
q=-1.6e-19
def ODE_integration_B(initial_position,initial_velocity,time_frame,n_steps):
loc=np.zeros((3,n_steps)) #3 for x,y,z
loc[:,0]=initial_position
v_old=initial_velocity
dt=(time_frame[1]-time_frame[0])/n_steps
B=[0,0,1]
f=np.zeros(3)
v_new=np.zeros(3)
for i in range(1,n_steps,1):
f=q*np.cross(v_old,B)
v_new=f/m*dt+v_old #new momentum
loc[:,i]=v_new*dt+loc[:,i-1] #new position
v_old=v_new
return loc
```

Let's try our code now!

In [24]:

```
r_initial=np.array([0,0,0])
v_initial=np.array([1000,0,0])
time_frame=[0,1e-9]
r=ODE_integration_B(r_initial,v_initial,time_frame,2000)
plot_2D(r[0,:],r[1,:],1,3,"axes")
```

The autoscale default from matplotlib does not always work and this is a good example. The distances, like the orbit of an electron inside a field of 1 T, are too small. We need to write an even better plotting function than `plot_2D`

. This function is called

```
def plot_2D_autoscale(x,y,n_interpolated=1,plot_size=3,plot_axes="none"):
```

Note that this function has three optional arguments. `n_interpolated`

corresponds to the number of interpolation points between time steps, as in `plot_2D`

. It's default value is set to `1`

, i.e. no interpolation. `plot_size`

sets the plot size to `3`

by default. Finally, `plot_axes`

is set to `"none"`

In [25]:

```
def plot_2D_autoscale(x,y,n_interpolated=1,plot_size=3,plot_axes="none"):
n = len(x)
k = n_interpolated
#interpolation
x_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, x)
y_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, y)
#generate figure and axes
fig, ax = plt.subplots(1, 1, figsize=(plot_size, plot_size))
ax.scatter(x_interpolated, y_interpolated, c=range(n*k), linewidths=0,
marker='o', s=2*plot_size, cmap=plt.cm.jet,) #s is the size of the plotted symbol 'o'
ax.axis('equal')
#compute autoscale parameters
xc=(x.max()+x.min())/2.
x_low=xc-(x.max()-x.min())/2.*1.1
x_high=xc+(x.max()-x.min())/2.*1.1
yc=(y.max()+y.min())/2.
y_low=yc-(y.max()-y.min())/2.*1.1
y_high=yc+(y.max()-y.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"):#so we can switch the axis on or off
ax.set_axis_off()
else:
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
```

In [26]:

```
plot_2D_autoscale(r[0,:],r[1,:],plot_axes="axes")
```

Clearly, this is not a circle. Using the same strategy to solve this ODE did not work well here, as we are trying to approximate a second order ODE using a first order approach... Before, the underlying structure of motion was close to a straight line (hence first order). Now we have circular motion, which is second order, and cannot be approximated by a straight line. Here again we need to develop a model that preserves the underlying nature of motion. Doing a second order integration scheme is possible with python.

However, rather than using this method we are going to use something quite exotic called exponential integrators....

`Lorentz_factor(v)`

that computes the Lorentz factor as a function of a velocity $\vec v=(v_x,v_y,v_z)$

In [27]:

```
def Lorentz_factor(v):
v2=norm(v)**2
c=299792458
return 1./sqrt(1-v2/c**2)
```

In [28]:

```
def get_E(r,t):
return np.array([0,0,0])
def get_B(r,t):
return np.array([0,0,1])
```

Let's add another function that computes the electromagnetic field tensor $F_{\mu\nu}$

In [29]:

```
def EM_field_tensor(E,B):
c=299792458
F_mu_nu=np.array([0,E[0]/c,E[1]/c,E[2]/c,E[0]/c,0,-B[2],B[1],E[1]/c,B[2],0,-B[0],E[2]/c,-B[1],B[0],0])
F_mu_nu.shape=(4,4)
return F_mu_nu
```

Let's now add two new functions and a library to measure the computational time

In [30]:

```
from scipy.linalg import expm
from numpy.linalg import norm
import time
```

In [31]:

```
m=9e-31
q=-1.6e-19
c=299792458
```

In [32]:

```
def ODE_integration_exponential_integrator(initial_position,initial_velocity,time_frame,n_steps,n_integrate):
start_time = time.time()
loc=np.zeros((3,n_steps)) #3 for x,y,z
loc[:,0]=initial_position
v=initial_velocity
gamma=Lorentz_factor(v)
p_old=np.array([m*c,m*v[0],m*v[1],m*v[2]])*gamma
p_new=np.zeros(4)
t=0
dt=(time_frame[1]-time_frame[0])/n_steps
ds=dt/n_integrate
for i in range(1,n_steps,1): #solution of the momentum equation
B=get_B(loc[:,i-1],t)
E=get_E(loc[:,i-1],t)
F_mu_nu=EM_field_tensor(E,B)
p_new=np.matmul(expm(q/(m*gamma)*F_mu_nu*dt),p_old)
gamma=p_new[0]/(m*c)
P=np.zeros((4,4))
s=0.
dPds=np.eye(4)*ds/2
x=np.zeros(4)
for j in range(0,n_integrate,1): #trapezoidal integration to get the trajectory
s+=ds
if (i>0):
P+=dPds
dPds=expm(q/(m*gamma)*F_mu_nu*s)*ds/2
if (i==0):
P+=dPds
P+=dPds
x=np.matmul(P,p_old)/(m*gamma)
loc[:,i]=x[1:4]+loc[:,i-1] #new position
p_old=p_new
t+=dt
#print("--- %s seconds ---" % (time.time() - start_time))
return loc
```

In [33]:

```
def ODE_integration_first_order_exponential_integrator(initial_position,initial_velocity,time_frame,n_steps):
start_time = time.time()
loc=np.zeros((3,n_steps)) #3 for x,y,z
loc[:,0]=initial_position
v_old=initial_velocity
p_old=np.zeros(4)
p_old[0]=m*c
p_old[1:4]=m*v_old
p_old*=Lorentz_factor(initial_velocity)
p_new=np.zeros(4)
dt=(time_frame[1]-time_frame[0])/n_steps
f=np.zeros(3)
v_new=np.zeros(3)
gamma=Lorentz_factor(v_old)
t=0
for i in range(1,n_steps,1):
B=get_B(loc[:,i-1],t)
E=get_E(loc[:,i-1],t)
F_mu_nu=EM_field_tensor(E,B)
p_new=np.matmul(expm(q/(m*gamma)*F_mu_nu*dt),p_old)
gamma=p_new[0]/(m*c)
v_new=p_new[1:4]/(m*gamma)
loc[:,i]=v_new*dt+loc[:,i-1] #new position
p_old=p_new
t+=dt
#print("--- %s seconds ---" % (time.time() - start_time))
return loc
```

In [34]:

```
r_initial=np.array([0,0,0])
v_initial=np.array([1000,0,0])
t_gyration=m/(abs(q)*norm(get_B(r_initial,0)))*(2*math.pi)
time_frame=[0,t_gyration]
print('Gyration time =',t_gyration,'s')
print('Larmor radius =',norm(v_initial)*m/(abs(q)*norm(get_B(r_initial,time_frame[0]))),'m')
r=ODE_integration_exponential_integrator(r_initial,v_initial,time_frame,150,6)
print(r[0,:].min(),r[0,:].max())
plot_2D_autoscale(r[0,:],r[1,:],n_interpolated=1,plot_size=5,plot_axes="axes")
```

In [35]:

```
r=ODE_integration_first_order_exponential_integrator(r_initial,v_initial,time_frame,150)
print('Larmor radius =',norm(v_initial)*m/(abs(q)*norm(get_B(r_initial,time_frame[0]))),'m')
print(r[0,:].min(),r[0,:].max())
plot_2D_autoscale(r[0,:],r[1,:],n_interpolated=1,plot_size=5,plot_axes="axes")
```

Larmor radius = 5.6250000000000016e-09 m -5.743220979103771e-09 5.507601530084524e-09

Let's find out.

In [36]:

```
time_frame=[0,109*t_gyration]
r=ODE_integration_first_order_exponential_integrator(r_initial,v_initial,time_frame,150)
print('Larmor radius =',norm(v_initial)*m/(abs(q)*norm(get_B(r_initial,time_frame[0]))),'m')
print(r[0,:].min(),r[0,:].max())
plot_2D_autoscale(r[0,:],r[1,:],n_interpolated=1,plot_size=5,plot_axes="axes")
```

Larmor radius = 5.6250000000000016e-09 m -2.980472608078362e-08 4.122206137687232e-09

How about the trapezoid integration?

In [37]:

```
r=ODE_integration_exponential_integrator(r_initial,v_initial,time_frame,150,6)
print('Larmor radius =',norm(v_initial)*m/(abs(q)*norm(get_B(r_initial,time_frame[0]))),'m')
print(r[0,:].min(),r[0,:].max())
plot_2D_autoscale(r[0,:],r[1,:],n_interpolated=1,plot_size=5,plot_axes="axes")
```

Larmor radius = 5.6250000000000016e-09 m -5.349733178899195e-09 5.349733178685792e-09

`ODE_integration_first_order_exponential_integrator`

had for one turn! You can check that both functions are actually giving good results, when the time step is small enough. However, you will discover that `ODE_integration_exponential_integrator`

does much better for larger time steps, especially if you are increasing `n_integrate`

. So, again, beware of time stepping! In fact, for inhomegeneous fields, this is `n_steps`

that has to be increased.

Now we combine the effect of both electric fields and magnetic fields. If a magnetic field is given by $\vec B=(B_x,B_y,B_z)$, then the particle will gyrate in the $(x,y)$ plane. But what happens if there is also and electric field along the z-axis? Note that by **parallel electric fields** we mean parallel to $\vec B$. In this case, we need to solve the following system of differentiual equations.
$$\frac{d^2v_x}{dt^2}=\frac{q}{m}\frac{dv_y}{dt}B_z\\
\frac{d^2v_y}{dt^2}=-\frac{q}{m}\frac{dv_x}{dt}B_z\\
\frac{dv_z}{dt}=\frac{q}{m}E_z$$
As before, we have a circular trajectory in the $(x,y)$ plane and a "free" fall along the $z$ axis.

In [38]:

```
def plot_3D_autoscale(loc,n_interpolated=1,plot_size=3,plot_axes="none"):
from mpl_toolkits.mplot3d import Axes3D
n = len(loc[0])
k = n_interpolated
#interpolation
x_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, loc[0,:])
y_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, loc[1,:])
z_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, loc[2,:])
#generate figure and axes
fig = plt.figure(figsize=(plot_size,plot_size))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x_interpolated, y_interpolated, z_interpolated, c=range(n*k), linewidths=0,
marker='o', s=2*plot_size, cmap=plt.cm.jet,) #s is the size of the plotted symbol 'o'
#compute autoscale parameters
xc=(loc[0,:].max()+loc[0,:].min())/2.
x_low=xc-(loc[0,:].max()-loc[0,:].min())/2.*1.1-1e-12
x_high=xc+(loc[0,:].max()-loc[0,:].min())/2.*1.1+1e-12
yc=(loc[1,:].max()+loc[1,:].min())/2.
y_low=yc-(loc[1,:].max()-loc[1,:].min())/2.*1.1-1e-12
y_high=yc+(loc[1,:].max()-loc[1,:].min())/2.*1.1+1e-12
zc=(loc[2,:].max()+loc[2,:].min())/2.
z_low=zc-(loc[2,:].max()-loc[2,:].min())/2.*1.1-1e-12
z_high=zc+(loc[2,:].max()-loc[2,:].min())/2.*1.1+1e-12
#set autoscale parameters
ax.set_xlim(min(x_low,y_low,z_low),max(x_high,y_high,z_high))
ax.set_ylim(min(x_low,y_low,z_low),max(x_high,y_high,z_high))
ax.set_zlim(min(x_low,y_low,z_low),max(x_high,y_high,z_high))
if (plot_axes!="axes"):#so we can switch the axis on or off
ax.set_axis_off()
else:
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
ax.set_zlabel("z (m)")
```

In [39]:

```
def get_E(r,t):
return np.array([0,300,10])
def get_B(r,t):
return np.array([0,0,-1])
r_initial=np.array([0,0,0])
v_initial=np.array([1000,0,0])
t_gyration=m/(abs(q)*norm(get_B(r_initial,time_frame[0])))*(2*math.pi)
time_frame=[0,3*t_gyration]
r=ODE_integration_exponential_integrator(r_initial,v_initial,time_frame,150,6)
plot_3D_autoscale(r,n_interpolated=10,plot_size=8,plot_axes="axes")
```