Python中四变量绘图：LG光束电场三维可视化求助

首先需要修正数据维度问题：你生成的r、theta、z是一维数组，直接计算E会出现维度不匹配，必须先将其转换为三维网格坐标，才能逐点计算每个位置的电场值。

一、数据网格化与电场计算示例 #

import numpy as np

# 生成一维变量数组
r = np.linspace(0, 10, 100)
theta = np.linspace(0, np.pi, 100)
z = np.linspace(0, 80, 100)

# 转换为三维网格（indexing='ij'保证维度与输入顺序对应）
R, Z, Theta = np.meshgrid(r, z, theta, indexing='ij')

# 替换为你自己的LG光束电场计算函数f(r,z,theta)
def calculate_lg_beam(r, z, theta, l=1, p=0, w0=1, wavelength=1):
    k = 2 * np.pi / wavelength
    zR = np.pi * w0**2 / wavelength  # 瑞利长度
    w_z = w0 * np.sqrt(1 + (z/zR)**2)
    R_z = z * (1 + (zR/z)**2) if z != 0 else np.inf
    psi_z = np.arctan(z/zR)
    
    # LG光束电场模值计算
    E = (np.sqrt(2 * np.math.factorial(p)) / (np.sqrt(np.pi) * np.math.factorial(p + l) * w0)) * \
        (np.sqrt(2) * r / w_z)**l * np.exp(-r**2 / w_z**2) * \
        np.abs(np.polynomial.hermite.hermval(2 * r**2 / w_z**2, [0]*p + [1])) * \
        np.exp(-1j * k * r**2 / (2 * R_z)) * np.exp(-1j * (2*p + l + 1) * psi_z) * np.exp(1j * l * theta)
    return np.abs(E)  # 取模值用于可视化

# 计算所有网格点的电场值
E = calculate_lg_beam(R, Z, Theta)

二、可视化方案选择 #

1. 三维散点图（Matplotlib） #

由于100×100×100的网格数据量达1e6，直接绘制所有点会严重卡顿，建议采样后再绘制：

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

# 采样间隔，可根据性能调整
sample_rate = 5
R_sample = R[::sample_rate, ::sample_rate, ::sample_rate].flatten()
Z_sample = Z[::sample_rate, ::sample_rate, ::sample_rate].flatten()
Theta_sample = Theta[::sample_rate, ::sample_rate, ::sample_rate].flatten()
E_sample = E[::sample_rate, ::sample_rate, ::sample_rate].flatten()

# 创建3D散点图
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
scatter = ax.scatter(R_sample, Z_sample, Theta_sample, c=E_sample, cmap='viridis', alpha=0.6)

# 设置标签与标题
plt.colorbar(scatter, label='电场模值 |E|')
ax.set_xlabel('r (横向距离)')
ax.set_ylabel('z (传播距离)')
ax.set_zlabel('theta (旋转参数)')
ax.set_title('LG光束电场三维散点图')
plt.show()

2. 切片曲面图（Matplotlib） #

直接展示三维全量数据难度大，可固定一个参数，绘制另外两个参数的电场分布：

# 固定theta=π/2，绘制r-z平面的电场分布
theta_fixed_idx = np.argmin(np.abs(theta - np.pi/2))
E_slice = E[:, :, theta_fixed_idx]

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
R_slice, Z_slice = np.meshgrid(r, z)
surf = ax.plot_surface(R_slice, Z_slice, E_slice, cmap='viridis', edgecolor='none')

plt.colorbar(surf, label='电场模值 |E|')
ax.set_xlabel('r')
ax.set_ylabel('z')
ax.set_zlabel('|E|')
ax.set_title(f'LG光束电场分布 (theta={np.pi/2:.2f})')
plt.show()

# 同理可固定z或r，绘制对应二维平面的电场分布

3. 体可视化与等值面（Mayavi） #

Mayavi对三维体数据的渲染更专业，适合展示电场的空间分布：

from mayavi import mlab

# 体绘制：展示电场的密度分布
mlab.figure(size=(800, 600))
vol = mlab.pipeline.volume(mlab.pipeline.scalar_field(R, Z, Theta, E), 
                           vmin=np.min(E), vmax=np.max(E))
mlab.xlabel('r')
mlab.ylabel('z')
mlab.zlabel('theta')
mlab.title('LG光束电场体可视化')
mlab.show()

# 等值面绘制：展示|E|等于特定值的空间曲面
mlab.figure(size=(800, 600))
mlab.contour3d(R, Z, Theta, E, contours=[np.max(E)*0.5], colormap='viridis')
mlab.xlabel('r')
mlab.ylabel('z')
mlab.zlabel('theta')
mlab.title('LG光束电场等值面')
mlab.show()

4. 交互式三维图（Plotly） #

Plotly生成的交互式图支持旋转、缩放、hover查看数值，适合精细分析：

import plotly.graph_objects as go

# 采样数据
sample_rate = 5
R_sample = R[::sample_rate, ::sample_rate, ::sample_rate].flatten()
Z_sample = Z[::sample_rate, ::sample_rate, ::sample_rate].flatten()
Theta_sample = Theta[::sample_rate, ::sample_rate, ::sample_rate].flatten()
E_sample = E[::sample_rate, ::sample_rate, ::sample_rate].flatten()

# 创建交互式散点图
fig = go.Figure(data=[go.Scatter3d(
    x=R_sample, y=Z_sample, z=Theta_sample,
    mode='markers',
    marker=dict(
        size=3,
        color=E_sample,
        colorscale='Viridis',
        opacity=0.7
    )
)])

# 更新布局
fig.update_layout(
    scene=dict(
        xaxis_title='r (横向距离)',
        yaxis_title='z (传播距离)',
        zaxis_title='theta (旋转参数)'
    ),
    title='LG光束电场交互式三维散点图'
)

fig.show()

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