4. 光束模式仿真基础：Python环境搭建与高斯光束传播

说实话，做光束模式识别，光有理论不行。你得亲手把光束“跑”出来看看。我见过不少同学，公式背得滚瓜烂熟，一打开IDE就懵了。嗯，这一章我们就来解决这个问题。

我会带你搭好Python环境，然后从零开始模拟一束高斯光束的传播。你想想看，当你在屏幕上看到那个光斑随着传播距离慢慢散开的时候，很多抽象的概念一下子就通了。

4.1 环境搭建：别小看这一步

我个人习惯用Anaconda来管理Python环境。为什么？因为它把NumPy、SciPy、Matplotlib这些科学计算库打包好了，省得你一个一个装，版本还容易打架。

如果你还没装，可以去Anaconda官网下载个Miniconda，轻量够用。装完之后，打开终端（Windows下是Anaconda Prompt），敲下面这行：

conda create -n beam_sim python=3.9
conda activate beam_sim
pip install numpy scipy matplotlib

等进度条跑完，环境就搭好了。我在项目中遇到过好几次因为库版本不对导致计算结果异常的情况——比如SciPy的傅里叶变换函数接口变了，代码跑出来全是NaN。所以建议你锁定版本：

pip install numpy==1.23.5 scipy==1.9.3 matplotlib==3.6.2

4.2 核心库速览：三剑客

我们这次主要用三个库。说白了，它们的分工很明确：

你不需要把每个函数都背下来。记住一个原则：NumPy造数据，SciPy算数据，Matplotlib看数据。这就够了。

4.3 高斯光束的数学表达

先回顾一下公式。基模高斯光束在z=0处的光场分布是：

E(x, y, 0) = exp(-(x² + y²) / w₀²)

其中w₀是束腰半径。传播到距离z处，光场变成：

E(x, y, z) = (w₀ / w(z)) · exp(-(x² + y²) / w(z)²) · exp(-i·k·z - i·k·(x²+y²)/(2R(z)) + i·ζ(z))

看着复杂？别怕。我们一步一步拆开写代码。我在项目中习惯先把参数定义清楚：

import numpy as np

# 物理参数
wavelength = 532e-9       # 波长 532nm（绿光）
w0 = 1e-3                 # 束腰半径 1mm
z = 0.5                   # 传播距离 0.5m
k = 2 * np.pi / wavelength # 波数

# 计算光束参数
zR = np.pi * w0**2 / wavelength  # 瑞利距离
wz = w0 * np.sqrt(1 + (z/zR)**2) # z处的光斑半径
Rz = z * (1 + (zR/z)**2)         # 波前曲率半径
zeta = np.arctan(z / zR)         # 古伊相位

4.4 构建二维网格与光场

要模拟二维光场，先得建一个网格。我建议网格范围取到束腰的3~5倍，不然边界效应会干扰结果。

# 网格参数
N = 512                   # 网格点数
L = 5 * w0                # 网格半宽
x = np.linspace(-L, L, N)
y = np.linspace(-L, L, N)
X, Y = np.meshgrid(x, y)
R2 = X**2 + Y**2

# z=0处的光场（复振幅）
E0 = np.exp(-R2 / w0**2)

# z处的光场
E = (w0/wz) * np.exp(-R2/wz**2) * \
    np.exp(-1j * k * z - 1j * k * R2/(2*Rz) + 1j * zeta)

这里用np.meshgrid生成二维坐标网格，R2是径向距离的平方。注意E是复数——振幅和相位都包含在里面了。

4.5 可视化：让光束“现形”

光场算出来了，得画出来看看。我一般同时看三个东西：光强分布、一维截面、相位分布。

import matplotlib.pyplot as plt

# 光强
I = np.abs(E)**2

plt.figure(figsize=(12, 4))

# 子图1：二维光强
plt.subplot(131)
plt.imshow(I, extent=[-L, L, -L, L], cmap='hot')
plt.colorbar(label='Intensity')
plt.title('Intensity at z=0.5m')

# 子图2：一维截面（取y=0）
plt.subplot(132)
plt.plot(x, I[N//2, :])
plt.xlabel('x (m)')
plt.ylabel('Intensity')
plt.title('Cross-section (y=0)')
plt.grid(True)

# 子图3：相位
plt.subplot(133)
phase = np.angle(E)
plt.imshow(phase, extent=[-L, L, -L, L], cmap='hsv', vmin=-np.pi, vmax=np.pi)
plt.colorbar(label='Phase (rad)')
plt.title('Phase at z=0.5m')

plt.tight_layout()
plt.show()

运行这段代码，你会看到：光强是高斯形状的圆斑，截面是一条钟形曲线，相位则呈现弯曲的波前——这就是高斯光束的“指纹”。

4.6 传播模拟：角谱法

上面我们用了解析公式。但如果光束不是理想高斯形，或者遇到复杂孔径，解析解就不存在了。这时候要用角谱传播法。

核心思想很简单：把光场分解成平面波，每个平面波传播一段距离，再合起来。数学上就是两次傅里叶变换。

from scipy.fft import fft2, ifft2, fftshift, ifftshift

def angular_spectrum_propagation(E0, wavelength, dx, z):
    """
    角谱法传播
    E0: 初始光场
    dx: 网格间距
    z:  传播距离
    """
    N = E0.shape[0]
    # 频率坐标
    fx = np.fft.fftfreq(N, d=dx)
    fy = np.fft.fftfreq(N, d=dx)
    FX, FY = np.meshgrid(fx, fy)
    F2 = FX**2 + FY**2
    
    # 角谱传递函数
    k = 2 * np.pi / wavelength
    H = np.exp(1j * k * z * np.sqrt(1 - wavelength**2 * F2))
    # 注意：当1 - wavelength**2 * F2 < 0时，对应倏逝波，这里直接置零
    H[np.imag(H) > 0] = 0
    
    # 传播
    A0 = fftshift(fft2(E0))
    A = A0 * H
    E = ifft2(ifftshift(A))
    
    return E

调用也很简单：

dx = 2*L / N
E_prop = angular_spectrum_propagation(E0, wavelength, dx, z=0.5)

对比一下解析解和角谱法的结果——你会发现它们几乎完全重合。这就是数值仿真的魅力：用通用方法验证特解。

4.7 本章知识结构

下面这张图把本章的核心逻辑串起来了。你可以把它当作一个“思维导图”来用：

4.8 避坑指南

最后，分享几个我踩过的坑：

• 网格间距太小或太大： 太小导致计算量爆炸，太大导致采样不足。我一般取dx = w0/10作为起点。
• 忘记归一化： 角谱法传播后，总能量可能会漂移。记得检查np.sum(np.abs(E0)**2)和np.sum(np.abs(E)**2)是否一致。
• 相位缠绕： np.angle返回的是主值(-π, π]，如果相位变化超过2π，画出来会有跳变。可以用np.unwrap解缠绕。
• 边界反射： 如果网格边界处光场没有衰减到0，FFT会产生伪影。解决办法是加一个孔径光阑或者扩大网格范围。

嗯，环境搭好了，代码能跑了，光束也画出来了。接下来你就可以拿着这套工具，去分析各种光束模式了。记住，仿真不是目的，理解物理才是。