Cylindrically symmetric magnetic field

Question

I want to plot the motion of a positive charge in a cylindrically symmetric magnetic field.
I am assuming a cylinder around the z-axis, with the magnetic field going in clockwise direction. The B-field has magnitude of 6T and the distance R from the z-axis is 3m. The charged particle is launched in positive direction along the z-axis and has the energy 2 MeV.
I am uncertain of how to simulate this B-field correctly. I was thinking to create the B-field in cylindrical coordinates,
cylinder from 0 to 2pi:
theta=numpy.linspace(0, 2*numpy.pi, 360)
x=r*numpy.cos(theta)
y=r*numpy.sin(theta)

Bx=B0*(numpy.cos(numpy.arctan2(y,x)
By=B0*(-numpy.sin(numpy.arctan2(y,x)))
Bz=0
And then create a vector B=[Bx, By, Bz] from which I would calculate the acceleration using Lorentz force for a timespan t.
But I think I am going in circles with this. Is there another way to create a cylindrically symmetric magnetic field?

Answer 1

With solve_ivp():

There are only two functions, the first one, B() takes care of the geometry of the magnetic field (rotationally invariant relative to the z axis and radially invariant, with always the same magnitude at each point), and the second one f() takes care of the physics, providing the velocity and the Lorentz's acceleration generated by the magnetic field B(), calculated at each point. This latter function is the right hand side of the differential equations of motion, which goes into solve_ivp().

'''
Dynamics in cylindrical homogeneous magnetic field
'''

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def B(y, B_templet):
    r = np.array([y[0], y[1], 0]) # vector aligned with the x,y projection of the position vector y
    r = r / np.sqrt(r[0]**2 + r[1]**2) # unit vector aligned with the x,y projection of the position vector y
    r_perp = np.array([-r[1], r[0], 0]) # unit vector perpendicular to the x,y projection of position vector y
    B_field = B_templet[0] * r  +  B_templet[1] * r_perp
    B_field[2] = B_templet[2]
    return B_field

def f(t, y, parameters):   
    mass, charge, B_templet = parameters
    return np.concatenate( (y[3:6], (charge/mass)*np.cross(y[3:6], B(y[0:3], B_templet))) )


charge = 1
mass = 1
B_direction = np.array([0.3,1,0.1])
B_magnitude = 1
xv_start = np.array([3, 0, 0, 0, 0, 2])
time_step = 0.01
n_iter = 5000
t_span=[0,n_iter*time_step]

B_direction = B_magnitude * B_direction / np.sqrt(B_direction.dot(B_direction))

sol = solve_ivp(fun = lambda t, y : f(t, y, (mass, charge, B_direction)), t_span=t_span, y0=xv_start, max_step=time_step)


fig = plt.figure()
ax = fig.add_subplot(projection='3d')
r = 35
ax.set_xlim((-r, r))
ax.set_ylim((-r, r))
ax.set_zlim((-r, r))
ax.plot(sol.y[0,:], sol.y[1,:], sol.y[2,:], 'r-')
ax.plot(sol.y[0,0], sol.y[1,0], sol.y[2,0], 'bo')
ax.plot(sol.y[0,-1], sol.y[1,-1], sol.y[2,-1], 'go')
plt.show()

Maybe this is what you want?

'''
Dynamics in a cylindrical magnetic field
'''

import numpy as np
import matplotlib.pyplot as plt

def B(x, B_templet):
    r = np.array([x[0], x[1], 0]) # vector aligned with the x,y projection of the position vector x
    r = r / np.sqrt(r[0]**2 + r[1]**2) # unit vector aligned with the x,y projection of the position vector x
    r_perp = np.array([-r[1], r[0], 0]) # unit vector perpendicular to the x,y projection of the position vector x
    B_field = B_templet[0] * r  +  B_templet[1] * r_perp
    B_field[2] = B_templet[2]
    return B_field

def f(x, mass, charge, B_templet):    
    return np.concatenate( (x[3:6], (charge/mass)*np.cross(x[3:6], B(x[0:3], B_templet))) )

def time_step_f(x, mass, charge, B_templet, t_step):
    k1 = f(x, mass, charge, B_templet)
    k2 = f(x + t_step*k1/2, mass, charge, B_templet)
    k3 = f(x + t_step*k2/2, mass, charge, B_templet)
    k4 = f(x + t_step*k3, mass, charge, B_templet)
    return x + t_step * (k1 + 2*k2 + 2*k3 + k4) / 6

def flow_f(x_initial, mass, charge, B_templet, t_step, n_iter):
    traj = np.empty( (n_iter, x_initial.shape[0]),  dtype = float )
    traj[0, :] = x_initial
    for m in range(n_iter-1):
        traj[m+1,:] = time_step_f(traj[m,:], mass, charge, B_templet, t_step)
    return traj


charge = 1
mass = 1
B_direction = np.array([0.3,1,0.1])
B_magnitude = 3
xv_start = np.array([3, 0, 0, 0, 0, 2])
time_step = 0.01
n_iter = 5000

B_direction = B_magnitude * B_direction / np.sqrt(B_direction.dot(B_direction))


xv = flow_f(xv_start, mass, charge, B_direction, time_step, n_iter)


fig = plt.figure()
ax = fig.add_subplot(projection='3d')

r = 20

ax.set_xlim((-r, r))
ax.set_ylim((-r, r))
ax.set_zlim((-r, r))

ax.plot(xv[:,0], xv[:,1], xv[:,2], 'r-')
ax.plot(xv[0,0], xv[0,1], xv[0,2], 'bo')
ax.plot(xv[-1,0], xv[-1,1], xv[-1,2], 'go')
plt.show()