LG轨道角动量激光脉冲3D绘图异常：未呈现预期旋转效果

问题描述 #

尝试绘制带轨道角动量的LG激光脉冲场强图，场强方程依赖圆柱坐标系的r、theta、z三个变量，且脉冲场随z方向传播会产生旋转效果。使用这三个坐标创建3D meshgrid生成彩色图后，结果未呈现预期的旋转特性，需排查问题。

现有代码 #

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

# Defined Parameters
p = 1  # The order of the Laguerre polynomial
m = 1 #The azimuthal order of the pulse
theta0 = 0
r0 = 20e-6
E0 = 1.22e+12 #V/m
k0 = 536625
r = np.linspace(0, 20e-6, 100)
theta = np.linspace(0, np.pi, 100)
z = np.linspace(-3, 0, 100)
f = np.sqrt(np.math.factorial(p)/(np.math.factorial(m+p)))*((np.sqrt(2)*r/r0)**m)*(2*r**2/r0**2)*np.cos(m*(theta-theta0))
# Define the equation E = f(r, theta, z)
def equation(r, theta, z):
    return E0*f*np.exp(-r**2/r0**2)*np.cos(k0*z)

# Generate data for the plot

X, Y, Z = np.meshgrid(r, theta, z)
E = equation(X, Y, Z)

# Create a 3D density plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot the density using a colormap
density_plot = ax.scatter(Z, Y, X, c=E, cmap='plasma', s=1)

# Add a colorbar to the plot
cbar = plt.colorbar(density_plot)
cbar.set_label('E Value')

# Set axis labels
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

# Show the plot
plt.show()

效果对比 #

实际生成图 #

预期效果 #

问题分析与修复方案 #

核心问题点 #

1. 拉盖尔多项式实现错误：当前代码用(2*r**2/r0**2)代替正规拉盖尔多项式，不符合LG脉冲的幅值分布公式。
2. 方位角范围不足：theta仅取0到π，未覆盖完整的0到2π方位角范围，丢失了一半的旋转结构。
3. 相位项缺失耦合：场强方程仅包含cos(k0*z)，缺少m*theta与传播相位的耦合项，无法产生随z传播的旋转效果。
4. 坐标映射错误：直接用圆柱坐标的r、theta、z作为3D散点的坐标，未转换为直角坐标，无法直观展示空间旋转结构；且散点坐标顺序混乱，导致轴映射错误。
5. 幅值项未适配网格：f定义在一维数组上，未适配meshgrid生成的三维网格，无法对每个空间点计算正确的幅值。

修复后的代码 #

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

# 参数定义
p = 1  # 拉盖尔多项式阶数
m = 1  # 方位角阶数
theta0 = 0
r0 = 20e-6
E0 = 1.22e+12  # V/m
k0 = 536625

# 生成坐标网格：r范围0到r0，theta覆盖0-2π，z取一个波长的传播距离
r = np.linspace(0, r0, 100)
theta = np.linspace(0, 2*np.pi, 100)
z = np.linspace(0, 2*np.pi/k0, 100)

# 转换为三维网格
R, Theta, Z = np.meshgrid(r, theta, z)

# 圆柱坐标转直角坐标，用于空间结构展示
X = R * np.cos(Theta)
Y = R * np.sin(Theta)

# 计算LG脉冲场强：包含正确的拉盖尔多项式与相位耦合
def lg_field(R, Theta, Z):
    # 拉盖尔多项式项
    lag_poly = laguerre(p)(2 * R**2 / r0**2)
    # 幅值因子
    amp_factor = np.sqrt(2 * np.math.factorial(p) / (np.pi * np.math.factorial(p + np.abs(m)))) * \
                 (np.sqrt(2)*R/r0)**np.abs(m) * lag_poly * np.exp(-R**2 / r0**2)
    # 耦合相位项：实现随传播距离的旋转效果
    phase = np.cos(m * Theta - k0 * Z + theta0)
    return E0 * amp_factor * phase

# 计算全场强分布
E = lg_field(R, Theta, Z)

# 绘制3D散点图
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection='3d')

# 用直角坐标展示，颜色对应场强值
density_plot = ax.scatter(X, Y, Z, c=E, cmap='plasma', s=1)

# 配置颜色条与轴标签
cbar = plt.colorbar(density_plot, pad=0.1)
cbar.set_label('场强 E (V/m)')
ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_zlabel('传播距离 Z (m)')

plt.show()

修复说明 #

• 拉盖尔多项式：使用scipy.special.laguerre计算正规的拉盖尔多项式，匹配LG脉冲的幅值分布规律。
• 方位角范围：将theta扩展到0到2π，完整展示脉冲的环形旋转结构。
• 耦合相位：加入m*Theta - k0*Z相位项，使场强随z传播时，方位角方向的相位同步变化，呈现自然的旋转效果。
• 直角坐标转换：将圆柱坐标转为直角坐标(X,Y,Z)，直观展示脉冲在三维空间中的旋转形态。
• 网格适配：在函数内部对三维网格逐点计算幅值与相位，确保每个空间点的场强计算准确。

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