1. 自由曲面描述方法对比：XY多项式与Zernike多项式

任务说明

本案例在 Python 中构建同一自由曲面（离轴抛物面），分别用 XY 多项式 与 Zernike 多项式 拟合其矢高，并对比拟合精度与近轴光线追迹结果。核心结论：XY 多项式在矩形孔径上对非对称面型拟合更灵活，Zernike 在圆形孔径上正交性更好但需足够项数才能逼近离轴面型。

环境与依赖

pip install numpy matplotlib

实现思路

完整代码

import numpy as np

# ---------- 1. 真实自由曲面：离轴抛物面 ----------
# z = (x^2 + y^2) / (4f)，孔径中心偏移 x0
f = 100.0
x0 = 15.0
def true_sag(x, y):
    return ((x - x0)**2 + y**2) / (4 * f)

# ---------- 2. 采样网格 ----------
N = 21
R = 10.0
# 矩形网格（XY 多项式适用）
xr = np.linspace(-R, R, N)
yr = np.linspace(-R, R, N)
Xr, Yr = np.meshgrid(xr, yr)
Zr = true_sag(Xr, Yr)

# 圆形网格（Zernike 适用）
r = np.linspace(0, R, N)
t = np.linspace(0, 2*np.pi, N)
Rr, Tt = np.meshgrid(r, t)
Xc = Rr * np.cos(Tt)
Yc = Rr * np.sin(Tt)
Zc = true_sag(Xc, Yc)

# ---------- 3. XY 多项式拟合 ----------
# 基函数：1, x, y, x^2, xy, y^2, x^3, x^2y, xy^2, y^3, x^4, ...
def xy_basis(x, y, n_terms):
    terms = []
    k = 0
    for deg in range(0, 10):
        for i in range(deg + 1):
            j = deg - i
            terms.append((x**i) * (y**j))
            k += 1
            if k >= n_terms:
                return np.array(terms)
    return np.array(terms)

n_xy = 15
A_xy = xy_basis(Xr.ravel(), Yr.ravel(), n_xy).T
coef_xy, *_ = np.linalg.lstsq(A_xy, Zr.ravel(), rcond=None)
Z_xy_fit = (A_xy @ coef_xy).reshape(Zr.shape)
res_xy = Zr - Z_xy_fit
rms_xy = np.sqrt(np.mean(res_xy**2))
pv_xy = np.max(np.abs(res_xy))

# ---------- 4. Zernike 多项式拟合 ----------
def zernike(n, m, rho, theta):
    """标准 Zernike 多项式（Noll 序简化实现）"""
    if m < 0:
        return zernike(n, -m, rho, theta) * np.sin(-m * theta) if False else \
               zernike_radial(n, -m, rho) * np.sin(-m * theta)
    return zernike_radial(n, m, rho) * np.cos(m * theta)

def zernike_radial(n, m, rho):
    m = abs(m)
    if (n - m) % 2 != 0:
        return np.zeros_like(rho)
    R = np.zeros_like(rho)
    for k in range((n - m) // 2 + 1):
        c = ((-1)**k * np.math.factorial(n - k) /
             (np.math.factorial(k) *
              np.math.factorial((n + m) // 2 - k) *
              np.math.factorial((n - m) // 2 - k)))
        R += c * rho**(n - 2*k)
    return R

# 归一化坐标
rho_c = Rr / R
theta_c = Tt

# 选取 Zernike 项（n 从 0 到 4）
zern_list = []
for n in range(0, 5):
    for m in range(-n, n + 1, 2):
        zern_list.append((n, m))

A_z = np.array([zernike(n, m, rho_c.ravel(), theta_c.ravel())
                for (n, m) in zern_list]).T
coef_z, *_ = np.linalg.lstsq(A_z, Zc.ravel(), rcond=None)
Z_z_fit = (A_z @ coef_z).reshape(Zc.shape)
res_z = Zc - Z_z_fit
rms_z = np.sqrt(np.mean(res_z**2))
pv_z = np.max(np.abs(res_z))

# ---------- 5. 近轴光线追迹对比 ----------
# 简化模型：光线沿 z 入射，反射后偏折角 ∝ 面型斜率
def slope_rms(Z, X, Y):
    dzdx = np.gradient(Z, axis=1) / np.gradient(X, axis=1)
    dzdy = np.gradient(Z, axis=0) / np.gradient(Y, axis=0)
    return np.sqrt(np.mean(dzdx**2 + dzdy**2))

slope_true = slope_rms(Zr, Xr, Yr)
slope_xy = slope_rms(Z_xy_fit, Xr, Yr)
slope_z = slope_rms(Z_z_fit, Xc, Yc)

# ---------- 6. 输出结果 ----------
print("=== 面型拟合精度 ===")
print(f"XY 多项式 ({n_xy}项):  RMS={rms_xy:.6f} mm,  PV={pv_xy:.6f} mm")
print(f"Zernike ({len(zern_list)}项): RMS={rms_z:.6f} mm,  PV={pv_z:.6f} mm")
print()
print("=== 近轴光线追迹（斜率 RMS）===")
print(f"真实面型:  {slope_true:.6f}")
print(f"XY 拟合:   {slope_xy:.6f}  (偏差 {abs(slope_xy-slope_true):.6f})")
print(f"Zernike:   {slope_z:.6f}  (偏差 {abs(slope_z-slope_true):.6f})")

运行与验证

python compare_sag.py

判定标准：

• XY 多项式 RMS 残差应显著小于 Zernike（因离轴抛物面在矩形孔径上更匹配 XY 基）。
• Zernike 需增加项数（n 到 6 以上）才能逼近同等精度，体现圆形正交基对非对称面型的局限。
• 斜率 RMS 偏差反映光线追迹差异：偏差越小，追迹结果越接近真实面型。

关键结论

• XY 多项式：适合矩形/非对称孔径，项数少即可高精度拟合离轴面型，但项间不正交，高阶易病态。
• Zernike：适合圆形孔径，项间正交、系数有像差物理意义，但描述离轴/非对称面型需更多项。
• 实际光学设计软件中，选择依据是孔径形状与面型对称性，而非单纯精度。