圆柱对称磁场

Question

我想画一个正电荷在圆柱对称磁场中的运动。
假设一个圆柱体绕着z轴旋转，磁场是顺时针方向的。B场的震级为6T，距离z轴的距离R为3m。带电粒子沿z轴正向发射，能量为2 MeV。
我不确定如何正确模拟这个B场。我想在柱坐标系中创建B场，
从0到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
然后创建一个向量B=[Bx, By, Bz]，从中我可以用洛伦兹力计算时间跨度t的加速度。
但我想我会在这件事上兜圈子。有没有其他方法可以产生圆柱对称的磁场

Answer 1

与solve_ivp()一起：

只有两个函数，第一个函数B()负责磁场的几何结构（相对于z轴旋转不变，并且径向不变，在每个点上的大小始终相同），第二个函数f()负责物理，提供由磁场B()产生的速度和洛伦兹加速度，在每个点计算。后一个函数是运动微分方程的右侧，进入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()

也许这就是你想要的

'''
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()