如何用Python生成贝塞尔光束3D图？代码问题求助

核心问题梳理 #

• 变量顺序颠倒：计算r、theta等变量时，x、y、z网格尚未生成，导致无效计算甚至报错。
• 取值范围无效：np.arange(0,50,100)的步长远大于区间长度，仅生成单个点，无法形成有效网格；且贝塞尔光束为对称分布，x/y需取正负区间。
• 光束公式错误：原代码中theta=k*r/n不符合贝塞尔光束的物理表达式，未区分径向与纵向波数分量。
• 绘图函数不匹配：plot_trisurf不适用于规则网格的3D强度可视化，无法呈现预期的光束结构。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv
from mpl_toolkits.mplot3d import Axes3D

# 定义物理参数
k = 2 * np.pi  # 波数
n = 1          # 贝塞尔光束阶数
alpha = np.pi / 12  # 圆锥角，控制光束扩散程度
k_rho = k * np.sin(alpha)  # 径向波数
k_z = k * np.cos(alpha)    # 纵向波数

# 定义取值范围（x/y对称分布，缩小步长提升分辨率）
x_range = np.arange(-25, 25, 0.5)
y_range = np.arange(-25, 25, 0.5)
z_range = np.arange(0, 10, 2)  # 选取多个z平面展示传播特性

# 生成x-y网格
x, y = np.meshgrid(x_range, y_range)
r = np.sqrt(x**2 + y**2)
phi = np.arctan2(y, x)

# 创建3D绘图对象
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# 遍历不同z平面绘制强度分布
for z in z_range:
    # 计算贝塞尔光束复振幅
    bessel_func = jv(n, k_rho * r)
    bessel_beam = bessel_func * np.exp(1j * (k_z * z - n * phi))
    intensity = np.abs(bessel_beam)**2
    
    # 将当前z平面的强度映射到3D空间
    z_plane = np.full_like(intensity, z)
    ax.plot_surface(x, y, z_plane, 
                    facecolors=plt.cm.viridis(intensity / np.max(intensity)),
                    rstride=2, cstride=2, alpha=0.7)

# 设置坐标轴与标题
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D Distribution of Bessel Beam Intensity')
plt.colorbar(plt.cm.ScalarMappable(cmap='viridis'), ax=ax, label='Intensity')
plt.show()

修正说明 #

1. 调整变量定义顺序，确保先初始化物理参数与网格，再进行光束计算。
2. 修正贝塞尔光束表达式，引入圆锥角区分径向/纵向波数，符合实际无衍射光束的传播规律。
3. 调整x/y为对称区间并缩小步长，保证光束的环形结构清晰可辨。
4. 采用多z平面的plot_surface绘图方式，直观展示贝塞尔光束在传播方向上的不变特性，匹配预期效果。

内容的提问来源于stack exchange，提问作者MV007