📚 课程目录第4章基于光线追迹的初始结构生成

4. 基于光线追迹的初始结构生成

任务定义

本案例要求根据目标成像要求（物方点光源经反射后成像于指定像点），利用光线追迹方法求解自由曲面反射镜的离散面形点云，进而生成可用于后续优化的初始结构。输入为目标成像参数（物点坐标、像点坐标、入瞳口径），输出为自由曲面反射镜的面形点云数据文件。

设计原理与求解路径

对于单个反射面，给定物点 A 和理想像点 B，根据费马原理，反射面上任意点 P 需满足光程为常数的条件：|PA| + |PB| = 常数。该常数由轴上参考光线确定。对于离轴非球面或自由曲面，这一条件直接给出每个采样点的矢高解。本案例采用等光程条件逐点求解面形，而非逐光线迭代逼近，因此计算确定且无收敛问题。

具体实现：在入瞳平面（垂直于系统光轴）上按极坐标或矩形网格采样光线入射点坐标 (x, y)，每条光线对应反射点 P 的横向坐标与入射点相同（薄镜近似下忽略高度差），通过等光程方程解析求解 P 的 z 坐标。求解采用数值二分法，在给定 z 搜索区间内寻找满足光程差小于阈值的解。

核心代码实现


import numpy as np

def solve_sag(x, y, A, B, L0, z_min=-100.0, z_max=100.0, tol=1e-10):
    """
    对单个采样点 (x, y) 求解自由曲面矢高 z，满足等光程条件。
    
    参数:
        x, y : 入瞳面采样点坐标（即反射点横向坐标）
        A    : 物点坐标数组 [Ax, Ay, Az]
        B    : 像点坐标数组 [Bx, By, Bz]
        L0   : 参考光程（轴上光线光程）
        z_min, z_max : 矢高搜索区间
        tol  : 光程差收敛阈值
    
    返回:
        z    : 该采样点的矢高解
    """
    def optical_path(z):
        P = np.array([x, y, z])
        return np.linalg.norm(P - A) + np.linalg.norm(P - B) - L0
    
    # 二分法求解
    f_min = optical_path(z_min)
    f_max = optical_path(z_max)
    
    if f_min * f_max > 0:
        raise ValueError(f"采样点 ({x}, {y}) 在搜索区间内无解，请扩大 z 范围")
    
    z_lo, z_hi = z_min, z_max
    for _ in range(200):
        z_mid = 0.5 * (z_lo + z_hi)
        f_mid = optical_path(z_mid)
        if abs(f_mid) < tol:
            return z_mid
        if f_min * f_mid < 0:
            z_hi = z_mid
            f_max = f_mid
        else:
            z_lo = z_mid
            f_min = f_mid
    return 0.5 * (z_lo + z_hi)

def generate_freeform_surface(A, B, D, Nx=64, Ny=64, z_min=-100.0, z_max=100.0):
    """
    根据目标成像要求生成自由曲面反射镜面形点云。
    
    参数:
        A    : 物点坐标 [Ax, Ay, Az] (mm)
        B    : 像点坐标 [Bx, By, Bz] (mm)
        D    : 入瞳口径（方形孔径边长，mm）
        Nx, Ny : 采样网格数
        z_min, z_max : 矢高搜索范围 (mm)
    
    返回:
        points : (Nx*Ny, 3) 数组，每行为一个面形点 [x, y, z]
    """
    # 参考光线：取入瞳中心 (0, 0) 对应的反射点
    # 轴上光线从 A 经 (0,0,z_ref) 到 B，先求 z_ref
    def axial_op(z):
        P = np.array([0.0, 0.0, z])
        return np.linalg.norm(P - A) + np.linalg.norm(P - B)
    
    # 轴上参考点通过一维搜索确定（取光程最小处作为初始反射点）
    z_candidates = np.linspace(z_min, z_max, 2000)
    op_vals = [axial_op(z) for z in z_candidates]
    z_ref = z_candidates[int(np.argmin(op_vals))]
    L0 = axial_op(z_ref)  # 参考光程
    
    # 生成采样网格
    x_vals = np.linspace(-D/2, D/2, Nx)
    y_vals = np.linspace(-D/2, D/2, Ny)
    
    points = []
    for y in y_vals:
        for x in x_vals:
            z = solve_sag(x, y, A, B, L0, z_min, z_max)
            points.append([x, y, z])
    
    return np.array(points)

# ========== 案例参数与执行 ==========
if __name__ == "__main__":
    # 目标成像要求
    A = np.array([0.0, 0.0, -200.0])    # 物点：位于光轴上 z=-200mm
    B = np.array([0.0, 50.0, 300.0])    # 像点：离轴像点 (y=50mm, z=300mm)
    D = 60.0                             # 入瞳口径 60mm × 60mm
    
    # 生成面形点云（为演示取较小网格）
    surf_pts = generate_freeform_surface(A, B, D, Nx=32, Ny=32)
    
    # 保存为 CSV 文件（x, y, z 三列）
    output_file = "freeform_initial_surface.csv"
    np.savetxt(output_file, surf_pts, delimiter=",", 
               header="x_mm,y_mm,z_mm", comments="")
    
    # 输出面形统计信息
    z_min_actual = np.min(surf_pts[:, 2])
    z_max_actual = np.max(surf_pts[:, 2])
    print(f"面形点云已生成: {surf_pts.shape[0]} 个点")
    print(f"矢高范围: [{z_min_actual:.4f}, {z_max_actual:.4f}] mm")
    print(f"矢高峰谷值 (PV): {z_max_actual - z_min_actual:.4f} mm")
    print(f"数据已保存至: {output_file}")

落地与验证方法

运行方式：将上述代码保存为 gen_freeform.py，在装有 NumPy 的 Python 环境中执行 python gen_freeform.py。代码无外部依赖，仅需标准库与 NumPy。

验证步骤：

运行后检查终端输出，确认面形点数为 32×32 = 1024，矢高 PV 值在合理范围（通常为数十毫米量级，取决于物像距与离轴量）。
打开生成的 freeform_initial_surface.csv，确认包含 1024 行数据，每行三列（x, y, z），无 NaN 或异常跳变值。
独立验证：任取文件中一个点 (x₀, y₀, z₀)，手动计算 |PA| + |PB|，与参考光程 L₀ 比较，差值应小于 1×10⁻⁶ mm（即代码中 tol 量级）。

判定标准：若所有采样点的光程差均小于 1×10⁻⁶ mm，且面形点云在口径范围内连续无跳变，则判定初始结构生成正确。该点云可直接作为后续自由曲面优化迭代的初始面形。