Sympy计算圆孔PSF报错‘Ei未定义’，求新手技术指导

问题概述 #

• 编程新手尝试用Sympy通过惠更斯积分计算圆孔PSF，因菲涅尔数F>1无法使用夫琅禾费近似
• 运行代码时触发name 'Ei' is not defined错误
• 绘制图像未得到预期的衍射分布结果

错误根源分析 #

1. 积分模型错误：原代码未体现圆孔的孔径积分范围，仅对径向坐标r做无界积分，且被积函数不符合惠更斯-菲涅尔原理的物理形式
2. 特殊函数兼容性问题：Sympy符号积分结果生成指数积分Ei，但numpy无对应实现，lambdify转换时导致未定义错误
3. 物理过程缺失：惠更斯积分需对圆孔做面积分（极坐标下径向+角度积分），原代码缺少角度积分步骤，也未限定孔径的径向边界

修正方案与步骤 #

• 重构积分模型：采用极坐标对圆孔孔径（径向ρ从0到圆孔半径a，角度θ从0到2π）进行积分，利用贝塞尔函数简化角度积分
• 替换特殊函数实现：改用数值积分避免符号积分的复杂结果，规避特殊函数的兼容性问题
• 正确计算PSF：PSF为复振幅的模平方，修正共轭相乘的不当简化操作

修正后的可运行代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

# 定义物理参数（SI单位）
wavelength = 1e-3  # 波长，单位m
a = 5e-1 / 2       # 圆孔半径，单位m
L = 5              # 观察面到孔径的距离，单位m
k = 2 * np.pi / wavelength  # 波数

# 计算菲涅尔数
F = a**2 / (wavelength * L)
print(f"菲涅尔数F: {F:.2f}")

# 定义观察点径向坐标r对应的PSF计算函数
def calculate_psf(r):
    # 利用贝塞尔函数简化角度积分：∫₀²π exp(1j*kρr cosθ/L) dθ = 2π J₀(kρr/L)
    # 径向积分被积函数
    integrand = lambda rho: rho * np.exp(1j * k/(2*L) * (rho**2 + r**2)) * np.sinc((k * r * rho)/(np.pi * L)) * 2 * np.pi
    # 拆分实部虚部做数值积分
    real_part, _ = quad(lambda rho: np.real(integrand(rho)), 0, a)
    imag_part, _ = quad(lambda rho: np.imag(integrand(rho)), 0, a)
    # 计算复振幅的模平方得到PSF
    amplitude = real_part + 1j * imag_part
    return np.abs(amplitude)**2

# 生成观察点径向坐标序列（适配孔径尺寸调整范围）
r_vals = np.linspace(0, 0.1, 100)
# 批量计算PSF值
psf_vals = np.array([calculate_psf(r) for r in r_vals])

# 归一化后绘图
psf_vals /= psf_vals.max()
plt.plot(r_vals, psf_vals)
plt.xlabel("径向距离r (m)")
plt.ylabel("归一化PSF")
plt.title(f"圆孔点扩散函数 (F={F:.2f})")
plt.grid(True)
plt.show()

关键说明 #

• 数值积分优势：避免符号积分产生的特殊函数兼容性问题，scipy的quad数值积分更适合工程场景的衍射计算
• 物理准确性：严格遵循惠更斯-菲涅尔原理的极坐标积分过程，利用贝塞尔函数简化角度积分，符合圆孔衍射的物理规律
• 可视化优化：对PSF做归一化处理，便于清晰观察衍射斑的强度分布特征

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