python - 圆柱对称磁场
问题描述
我想绘制圆柱对称磁场中正电荷的运动。
我假设一个围绕 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 的加速度。
但我想我正在绕圈子。是否有另一种方法来创建圆柱对称磁场?
解决方案
与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()
推荐阅读
- swift - 如何在 Swift 中向 NSString 添加方法以设置 Self 的值?
- json - 如何添加从 forEach 获得的值?(反应)
- angular - 调试器在我的 Ionic(angular) 项目 (VSCODE) 中不起作用
- django - 如何插入到与Django ORM中的外键连接的多个postgres表中
- apache-kafka - 数据不一致的 kafka 代理 - NotLeaderForPartitionError
- c# - 为什么从 .net core 3.1 迁移到 .net 5 时 JSON 返回值发生了变化
- sass - Angular 11=>12 材质更新 SCSS
- javascript - 在加载到 Firebase 数据库之前读取本地文件有问题吗?
- python - TensorFlow 重复成功消息和 NUMA 节点读取警告
- mongodb - 如何打印 mongoDB 查询说明和在 Scala 中执行查询所需的时间?