ColumX 3D 相对论光线追踪实战指南

电子在电磁场中的运动是电子光学的基础。ColumX 的 ray3d 模块提供了一个完整的三维相对论 Lorentz 力光线追踪器，使用 SciPy 自适应步长 RK45 积分器，绝对容差可设至 10^-14。本文将从基础原理出发，通过 6 个实战案例，手把手教你使用这个模块。

一、物理原理：相对论 Lorentz 力方程

ColumX 的 ray3d 求解的是完整的相对论 Lorentz 力方程：

展开为速度导数形式：

其中 γ = 1/√(1 - v²/c²) 是 Lorentz 因子。这个公式包含了所有相对论修正项，在 200-300 kV 的 TEM 工作电压下，电子速度已达光速的 70-78%，相对论效应不可忽略。

积分器使用 SciPy 的 solve_ivp，默认 Dormand-Prince RK45 自适应步长方法，容差设置为 rtol=1e-10, atol=1e-14，确保能量守恒精度优于 10^-8。

二、快速上手：5 分钟追踪第一条光线

安装 ColumX 后，只需 10 行代码就能追踪一条电子：

from columx.ray3d import trace_ray, uniform_B, zero_E
import numpy as np

# 初始状态：[x, y, z, vx, vy, vz] (SI 单位)
state0 = [0, 0, 0,       # 位置：原点
          1e6, 0, 2e7]   # 速度：x 方向 1e6 m/s，z 方向 2e7 m/s

# 均匀磁场 Bz = 0.5 T
B_func = uniform_B(0.5)

# 追踪 2 纳秒
result = trace_ray(zero_E, B_func, state0, t_span=(0, 2e-9))

# 提取轨迹
x_mm = result.y[0] * 1e3  # 转换为 mm
z_mm = result.y[2] * 1e3
print(f"终点 z = {z_mm[-1]:.2f} mm, 共 {len(result.t)} 步")

输出类似：终点 z = 42.38 mm, 共 187 步。自适应步长积分器自动调整步长，在曲率大的地方加密，在直线段稀疏，既保证精度又高效。

三、案例 1：均匀磁场中的螺旋运动

电子在均匀磁场 B_z 中会做螺旋运动（Larmor 回旋），这是磁透镜聚焦的物理基础。让我们追踪一个 200 kV 的电子，初始速度有轴向分量和径向分量：

from columx.ray3d import trace_ray, uniform_B, zero_E
from columx.constants import q, m, c_light
import numpy as np

V = 200e3  # 200 kV
gamma = 1 + q * V / (m * c_light**2)  # 相对论因子
v_total = c_light * np.sqrt(1 - 1/gamma**2)
vz = v_total * 0.95   # 95% 轴向
vr = v_total * np.sqrt(1 - 0.95**2)  # 径向分量

state0 = [0, 0, -0.01, vr, 0, vz]  # 从 z=-10mm 出发

result = trace_ray(zero_E, uniform_B(0.5), state0, t_span=(0, 2e-9))

四、案例 2：平行电子束在不同磁场中的行为

TEM 中的电子束通常近似为平行束。ColumX 提供了 generate_parallel_beam 工具函数，自动生成同心圆环分布的初始条件（"乌克兰坟墓"分布）：

from columx.ray3d import generate_parallel_beam, trace_ray, uniform_B, zero_E

# 12 条光线，束半径 1mm，200 kV
states0 = generate_parallel_beam(n_rays=12, beam_radius=1e-3, V=200e3)

# 在不同 Bz 下追踪
for Bz in [0.0, 0.1, 0.5, 1.0, 2.0]:
    B_func = uniform_B(Bz)
    for state in states0:
        result = trace_ray(zero_E, B_func, state, t_span=(0, 1e-9))
        # 绘图 ...

五、案例 3：点源发射与锥形束发散

场发射电子源可以近似为点源，电子在锥角内均匀发射。generate_point_beam 生成锥形束初始条件：

from columx.ray3d import generate_point_beam, trace_ray, uniform_B, zero_E
import numpy as np

# 20 条光线，200 kV，半角 15 度
states0 = generate_point_beam(n_rays=20, V=200e3, half_angle=np.radians(15))

# 无磁场（纯漂移）
for state in states0:
    result = trace_ray(zero_E, uniform_B(0.0), state, t_span=(0, 0.5e-9))

六、案例 4：能量守恒验证

数值积分的精度可以通过能量守恒来验证。在纯磁场中（E=0），Lorentz 力不做功，电子动能应严格守恒。让我们验证不同电压下的能量守恒精度：

from columx.ray3d import trace_ray, uniform_B, zero_E
from columx.constants import q, m, c_light
import numpy as np

V = 200e3
gamma0 = 1 + q * V / (m * c_light**2)
v = c_light * np.sqrt(1 - 1/gamma0**2)

# 初始状态：三个方向都有速度分量
state0 = [0, 0, 0, v*0.3, v*0.2, v*0.9]

result = trace_ray(zero_E, uniform_B(0.8), state0,
                   t_span=(0, 3e-9), rtol=1e-10, atol=1e-14)

# 计算每个时间步的动能
vx, vy, vz = result.y[3], result.y[4], result.y[5]
speed = np.sqrt(vx**2 + vy**2 + vz**2)
gamma_t = 1.0 / np.sqrt(1 - speed**2 / c_light**2)
ke_eV = (gamma_t - 1) * m * c_light**2 / q

# 相对误差
error = np.abs(ke_eV - ke_eV[0]) / ke_eV[0]
print(f"最大相对误差: {error.max():.2e}")

七、案例 5：束流截面演化分析

观察束流在不同 z 平面的截面分布，是评估束流品质（均匀性、旋转对称性）的重要手段：

from columx.ray3d import generate_parallel_beam, trace_ray, uniform_B, zero_E
import numpy as np

# 30 条光线，束半径 1.5mm，2 度角发散
states0 = generate_parallel_beam(30, beam_radius=1.5e-3,
                                 V=200e3, angular_spread=np.radians(2))

all_results = []
for s in states0:
    res = trace_ray(zero_E, uniform_B(0.5), s, t_span=(0, 1.5e-9))
    all_results.append(res)

# 在 z=0, 2, 5, 10 mm 处取截面
for z_target in [0, 2, 5, 10]:  # mm
    xs, ys = [], []
    for res in all_results:
        z_arr = res.y[2] * 1e3
        i = np.argmin(np.abs(z_arr - z_target))
        xs.append(res.y[0, i] * 1e3)
        ys.append(res.y[1, i] * 1e3)
    # 绘制 scatter(xs, ys) ...

八、案例 6：不同加速电压的轨迹对比

加速电压直接影响电子的相对论因子 γ，进而改变螺旋轨迹的半径和频率。对比不同电压下的轨迹，可以直观理解 TEM 电压选择对光学性能的影响：

voltages = [20e3, 50e3, 100e3, 200e3, 300e3]
for V in voltages:
    gamma = 1 + q * V / (m * c_light**2)
    v = c_light * np.sqrt(1 - 1/gamma**2)
    state0 = [1e-3, 0, 0, 0, v*0.15, v*0.98]
    result = trace_ray(zero_E, uniform_B(0.5), state0, t_span=(0, 2e-9))
    # 绘制轨迹 ...

九、自定义场函数：接入 FEM 磁场

ray3d 的强大之处在于它接受任意场函数。你可以将 ColumX 的 FEM 求解结果作为输入，追踪电子在真实磁透镜场中的轨迹：

from columx.fem2d import AxisymmetricLens, bh_mild_steel
from columx.ray3d import trace_ray, zero_E
import numpy as np

# 1. 用 FEM 求解磁透镜场
solver, coil_info = AxisymmetricLens.magnetic_lens_full(
    NI=3000, r_pole=0.005, gap=0.004,
    z_positions=(-0.04, 0.04), r_max=0.05,
    iron_mu_r=2000.0, iron_bh_curve=bh_mild_steel,
    iron_r_outer=0.035, pole_length=0.012, yoke_width=0.006,
    nr=25, nz=50,
)
A = solver.solve_magnetic_newton_raphson(
    A_func=coil_info['J_func'], max_iter=30, tol=1e-4,
)

# 2. 构建插值场函数（从 FEM 结果）
from scipy.interpolate import RegularGridInterpolator

# 在网格上计算 B_r 和 B_z（使用 get_field_2d 批量计算）
r_grid = np.linspace(0.001, 0.01, 50)  # 避开 r=0 奇点
z_grid = np.linspace(-0.03, 0.03, 100)
R, Z = np.meshgrid(r_grid, z_grid, indexing='ij')
Br_flat, Bz_flat = solver.get_field_2d(
    R.ravel(), Z.ravel(), field_type='magnetic'
)
Br = Br_flat.reshape(R.shape)
Bz = Bz_flat.reshape(R.shape)

Br_interp = RegularGridInterpolator((r_grid, z_grid), Br,
                                     bounds_error=False, fill_value=0)
Bz_interp = RegularGridInterpolator((r_grid, z_grid), Bz,
                                     bounds_error=False, fill_value=0)

# 3. 定义 ray3d 格式的场函数
def B_func_lens(t, r):
    """B(t, [x,y,z]) -> [Bx, By, Bz]"""
    x, y, z = r
    rho = np.sqrt(x**2 + y**2)
    if rho < 1e-10:
        return np.array([0.0, 0.0, Bz_interp([[0.001, z]])[0]])
    Br_val = Br_interp([[rho, z]])[0]
    Bz_val = Bz_interp([[rho, z]])[0]
    # 将柱坐标分量转换为笛卡尔
    cos_phi = x / rho
    sin_phi = y / rho
    return np.array([Br_val * cos_phi, Br_val * sin_phi, Bz_val])

# 4. 追踪电子通过磁透镜
V = 200e3
gamma = 1 + 1.602e-19 * V / (9.109e-31 * (3e8)**2)
v = 3e8 * np.sqrt(1 - 1/gamma**2)
state0 = [0.5e-3, 0, -0.03, 0, 0, v]  # 偏离轴 0.5mm

result = trace_ray(zero_E, B_func_lens, state0, t_span=(0, 5e-10))
z_mm = result.y[2] * 1e3
x_mm = result.y[0] * 1e3
print(f"电子从 z={z_mm[0]:.1f}mm 到 z={z_mm[-1]:.1f}mm")
print(f"最大偏转 x = {np.max(np.abs(x_mm)):.3f} mm")

十、API 速查表

场函数接口规范

场函数签名为 f(t, r) -> array[3]，其中 t 是时间（秒），r = [x, y, z] 是位置（米），返回值是 [Ex/Bx, Ey/By, Ez/Bz]（V/m 或 Tesla）。

积分器参数

安装

pip install columx==2.3.6

ColumX v2.3.6 · Shuimu Scientific · columx.com · PyPI