梯度折射率介质内的光线追迹算法

借助Gemini 2.5 pro，综合了ChatGPT和Gemini的意见。为Gemini 2.5 pro撰写

目的

• 梯度折射率介质内光线追迹的数学模型

• 以Gradient3面型为例，建立数学模型。这是一个在眼视光学中很常用的面型，可用于多种模型眼中。该介质的折射率定义为：

  $$n(x,y,z) = n_0 + n_{r2}(x^2+y^2) + n_{r4}(x^2+y^2)^2 + n_{r6}(x^2+y^2)^3 + n_{z1}z + n_{z2}z^2 + n_{z3}z^3$$

• 使用4阶Runge-Kutta 算法(RK4)，对Gradient3面型内的光线追迹给出数值算法。

教科书内容

《梯度折射率光学》ISBN 7-03-002306-4/O.433

在设计光学仪器时，按照光学理论，准确地得到通过光学系统的光线径迹的过程，人们通常称为光线追迹。我们知道，由任何一点所发出的光线有无数条，我们只能选取其中具有代表性的几条光线进行追迹，以便根据这些光线在离断面上交点的分布来“标明”这个光学系统的完善程度，因此准确地定义应该称为有限光线追迹。

在均匀介质（光学玻璃、塑料、镜等）中光线是直线传播的，而在梯度折射率介质中，光线是曲线传播的，而且只有在少数情况下可以求得光线方程的解析解，这就懂得光线追迹显得复杂而重要。梯度折射率介质的光线追迹必须在电子计算机上进行。目前为了追迹梯度折射率介质中的光线，已经研究了一些方法，但是这些追迹方法和程序编制的正确与否，常常需要判别。这两个判别公式是光线传播理论一章导出的：

和

前者对于任意情况都是适用的，后者对检验旋转对称系统中的空间光线追迹是正确的。同时，一些呈特殊函数分布介质的解析解不仅描述了那种介质本身的特性，而且为各种数值计算方法提供了检验和判断的依据。当一个计算程序通过时，首先必须用来追迹一些具有解析解的介质内的光线，以判断其计算的可靠性。光线传播理论一章已为设计者提供了考核光线追迹程序的多种典型分布函数，光线追迹解析式及其光学特性。

光线追迹这一章，首先讨论梯度折射率介质中求解光线方程的一些基本方法，也就是求光线方程数值解的一些方法。这不仅是光线理论用于波导中追迹计算所必须的，而且是光学系统中利用梯度折射率材料时，进行光线追迹的基础。其次将讨论含有梯度折射率元件的光学系统的光线追迹问题，目的在于解决光学设计中的像差计算问题。

Runge-Kutta 法

高于一阶的任意高阶微分方程，从理论上说，都能化为一阶常微分方程组。光线微分方程是一个二阶微分方程，我们可以把它化为一对一阶微分方程组，解含有两个一阶微分方程组的方程组的一种典型方法就是 Runge-Kutta 法。

我们令光学光线矢

式中 t 是一个引入的新参量，由式中可以看出

利用该参量可将光线方程 [式(4.16)]

写为

这个方程有三个分量，它们可以用三元一维数组同时解得。令

则该方程又可写为

这个方程在初始条件为$R=R_0(x_0, y_0, z_0)$，$T=T_0$时可解得，由已知点 $(R_0, T_0)$ 出发，可逐次得到$(R_1, T_1), (R_2, T_2), \dots, (R_n, T_n)$ ，也就是利用 Runge-Kutta 算法追迹通过梯度折射率介质的光线，其算法可以写作

这里矩阵

式中 $\Delta t$ 是外推距离，是 $t$ 的增量，我们在计算中用$\Delta t$ ，而不用 t。由式(4.33)可以看出，外推距离是几何步长 $\Delta s$ 与折射率的商：$\Delta t =\frac{\Delta s}{n}$，外推距离 $\Delta t$ 小，则精度好。此法与其他方法比较，一些学者认为在达到同样精度的情况下，这种外推法花费的运算量要少得多。

【所有的AI都认为4.37和4.38有问题】
【建议修订为：】

把二阶方程 $\dfrac{d^2R}{dt^2}=D(R)$ 化为一阶系统

$$\begin{cases} \dfrac{dR}{dt}=T,\\[4pt] \dfrac{dT}{dt}=D(R). \end{cases}$$

标准 RK4为：

$$\begin{aligned} K_{1,R}&=T_n, &\quad K_{1,T}&=D(R_n),\\[4pt] K_{2,R}&=T_n+\dfrac{\Delta t}{2}K_{1,T}, &\quad K_{2,T}&=D\!\Big(R_n+\dfrac{\Delta t}{2}K_{1,R}\Big),\\[4pt] K_{3,R}&=T_n+\dfrac{\Delta t}{2}K_{2,T}, &\quad K_{3,T}&=D\!\Big(R_n+\dfrac{\Delta t}{2}K_{2,R}\Big),\\[4pt] K_{4,R}&=T_n+\Delta t\,K_{3,T}, &\quad K_{4,T}&=D\!\big(R_n+\Delta t\,K_{3,R}\big). \end{aligned}$$

修正后的 4.37：

$$\boxed{\;R_{n+1}=R_n+\dfrac{\Delta t}{6}\big(K_{1,R}+2K_{2,R}+2K_{3,R}+K_{4,R}\big)\; }$$

修正后的 4.38：

$$\boxed{\;T_{n+1}=T_n+\dfrac{\Delta t}{6}\big(K_{1,T}+2K_{2,T}+2K_{3,T}+K_{4,T}\big)\; }$$

如需用书中原来用的 A,B,C 符号表示，可定义

$$A=\Delta t\,K_{1,T},\quad B=\Delta t\,K_{2,T},\quad C=\Delta t\,K_{3,T},\quad D=\Delta t\,K_{4,T},$$

则等价写为

$$T_{n+1}=T_n+\dfrac{1}{6}\,(A+2B+2C+D),$$

使用 Runge-Kutta 法，也可以把光线微分方程作如下改造。

仍令

则光线方程

可写为

或

因为

上式独立变量仅有四个，且利用哈密顿算符子

将改写后的光线方程化为四个一阶微分方程

这样一来就使计算量大为减少了。

近轴情况下，令

则由近轴光线微分方程

得到

为了计算介质内光程差，需对光程

微分，利用关系式(4.42)和(4.43)得到

利用四阶 Runge-Kutta 法计算式(4.44)和式(4.48)是很方便的，共有五个变量，使用五元一维矩阵即可同时求解。

令

【公式4.49与公式4.44有矛盾，第3，4项缺了/H】
【建议修订为：】

$$r = \begin{pmatrix} x \\ y \\ T_x \\ T_y \\ [OP] \end{pmatrix}, \quad F = -\begin{pmatrix} T_x/H \\ T_y/H \\ n\frac{\partial n / \partial x}{H} \\ n\frac{\partial n / \partial y}{H} \\ n^2/H \end{pmatrix} \quad \tag{4.49}$$

则数值计算公式为

式中

$\Delta z$ 为 Z 轴方向增量，即外推步长。

Gradient3介质

“Gradient3”介质其折射率 n 和梯度 ∇n 的计算公式如下：

1. 折射率 n(x, y, z):

首先计算径向距离的平方： $r^2 = x^2 + y^2$

$$n(x,y,z) = n_0 + n_{r2}(x^2+y^2) + n_{r4}(x^2+y^2)^2 + n_{r6}(x^2+y^2)^3 + n_{z1}z + n_{z2}z^2 + n_{z3}z^3$$

2. 折射率梯度 ∇n = (∂n/∂x, ∂n/∂y, ∂n/∂z):

在Gradient3介质内光线追迹的完整数学模型与数值算法

第1步：建立通用数学模型（以z为自变量）

1. 基本光线方程:

光线在梯度折射率介质中的传播路径 $\vec{r}(s) = (x(s), y(s), z(s))$ 遵循以下微分方程，其中

2. 引入光学方向矢量:

定义光学方向矢量 $\vec{T} = n \frac{d\vec{r}}{ds}$ ，其分量为$(T_x, T_y, T_z)$ 。该矢量的模长平方为$|\vec{T}|^2 = T_x^2 + T_y^2 + T_z^2 = n^2$ 。基本方程可化为一阶方程组：

3. 更换自变量为 z:

在光学系统中，我们通常希望追踪光线随光轴坐标z的变化。利用链式法则$\frac{d}{dz} = \frac{ds}{dz} \frac{d}{ds}$ ，以及关系式$\frac{dz}{ds} = \frac{T_z}{n}$ ，可得$\frac{ds}{dz} = \frac{n}{T_z}$ 。

为了与教科书统一并处理光线传播方向，定义哈密顿量 H：

4. 最终模型:

将自变量替换后，我们得到一个由五个常微分方程组成的方程组，它描述了光线的状态（位置x, y，方向Tx, Ty，以及光程OP）如何随z演化：

第2步：代入Gradient3面型公式

将Gradient3的折射率及其偏导数公式代入上述模型的函数 $\vec{F}$ 中。

1. 折射率 n(x, y, z):

   $$n(x,y,z) = n_0 + n_{r2}(x^2+y^2) + n_{r4}(x^2+y^2)^2 + n_{r6}(x^2+y^2)^3 + n_{z1}z + n_{z2}z^2 + n_{z3}z^3$$

2. 折射率偏导数:

   $$\frac{\partial n}{\partial x} = x \left[ 2n_{r2} + 4n_{r4}(x^2+y^2) + 6n_{r6}(x^2+y^2)^2 \right]$$

   $$\frac{\partial n}{\partial y} = y \left[ 2n_{r2} + 4n_{r4}(x^2+y^2) + 6n_{r6}(x^2+y^2)^2 \right]$$

在每次计算导数函数 $\vec{F}$ 时，都需要根据当前的状态 $(x, y, z)$ 计算出相应的$n$,$\frac{\partial n}{\partial x}$, $\frac{\partial n}{\partial y}$ 和 H。

第3步：使用4阶龙格库塔(RK4)的数值算法

我们使用RK4方法求解上述微分方程组。从初始状态 $\vec{v}_0$ 在 $z_0$ 处开始，以步长 $\Delta z$ 迭代计算，得到下一位置 $z_{i+1} = z_i + \Delta z$ 处的状态 $\vec{v}\_{i+1}$。

迭代公式:

其中：

第4步：算法伪代码（综合优化版）

以下伪代码结合了AI_2的清晰结构和AI_1的高效计算逻辑。

// -------------------------------------------------------------------------
// 函数: F_derivative(z, v, params)
// 功能: 计算在给定点(z, v)的状态导数向量 F = dv/dz
// 输入:
//   z: 当前的z坐标 (浮点数)
//   v: 当前的状态向量 [x, y, Tx, Ty, OP] (5维数组)
//   params: Gradient3介质的系数 (结构体或对象) 
// 输出:
//   导数向量 F (5维数组)，如果发生全反射则返回错误标记（如 null）
// -------------------------------------------------------------------------
FUNCTION F_derivative(z, v, params):
    // 1. 解包状态向量
    x = v[0], y = v[1], Tx = v[2], Ty = v[3]

    // 2. 高效计算 n 和其偏导数
    r_sq = x*x + y*y
    r_sq_2 = r_sq * r_sq
    r_sq_3 = r_sq_2 * r_sq
    
    n = params.n0 + params.nr2*r_sq + params.nr4*r_sq_2 + params.nr6*r_sq_3 + 
        params.nz1*z + params.nz2*z*z + params.nz3*z*z*z

    g_r = 2*params.nr2 + 4*params.nr4*r_sq + 6*params.nr6*r_sq_2
    dn_dx = x * g_r
    dn_dy = y * g_r

    // 3. 计算哈密顿量 H 并检查合法性
    n_sq = n*n
    T_sq_transverse = Tx*Tx + Ty*Ty
    
    radicand = n_sq - T_sq_transverse // 根号内的值
    IF radicand <= 0 THEN
        // 发生全反射或光线与z轴垂直，无法继续以z为步长追迹
        RETURN null // 返回错误标记
    END IF
    
    H = -sqrt(radicand)

    // 4. 计算并返回导数向量 F
    inv_H = 1.0 / H // 预先计算倒数，避免多次除法
    F = [
        -Tx * inv_H,
        -Ty * inv_H,
        -n * dn_dx * inv_H,
        -n * dn_dy * inv_H,
        -n_sq * inv_H
    ]
    RETURN F
END FUNCTION


// -------------------------------------------------------------------------
// 主程序：RayTrace_Gradient3_RK4
// -------------------------------------------------------------------------
PROCEDURE RayTrace_Gradient3_RK4:
    // 1. 初始化
    // 定义介质参数
    gradient_params = {n0, nr2, nr4, nr6, nz1, nz2, nz3} 
    
    // 定义光线初始状态
    z_start = 0.0
    z_end = 10.0
    step_size = 0.1
    
    initial_pos = [x0, y0]
    initial_dir = [dir_x0, dir_y0, dir_z0] // 单位方向矢量
    
    // 根据初始条件计算初始状态向量
    n_start = calculate_n(z_start, initial_pos, gradient_params) // 需要辅助函数
    Tx0 = n_start * initial_dir[0]
    Ty0 = n_start * initial_dir[1]
    
    z = z_start
    v = [initial_pos[0], initial_pos[1], Tx0, Ty0, 0.0] // [x, y, Tx, Ty, OP]

    // 存储轨迹
    ray_path = []
    ADD {z: z, state: v} TO ray_path

    // 2. RK4 主循环
    WHILE z < z_end:
        // 确保最后一步正好落在 z_end 上
        current_step = min(step_size, z_end - z)
        
        // 计算 K1
        K1 = F_derivative(z, v, gradient_params)
        IF K1 IS null THEN BREAK // 错误处理

        // 计算 K2
        v_temp2 = v + (current_step/2.0) * K1
        K2 = F_derivative(z + current_step/2.0, v_temp2, gradient_params)
        IF K2 IS null THEN BREAK // 错误处理

        // 计算 K3
        v_temp3 = v + (current_step/2.0) * K2
        K3 = F_derivative(z + current_step/2.0, v_temp3, gradient_params)
        IF K3 IS null THEN BREAK // 错误处理
        
        // 计算 K4
        v_temp4 = v + current_step * K3
        K4 = F_derivative(z + current_step, v_temp4, gradient_params)
        IF K4 IS null THEN BREAK // 错误处理

        // 3. 更新状态向量和z坐标
        v = v + (current_step / 6.0) * (K1 + 2*K2 + 2*K3 + K4)
        z = z + current_step

        ADD {z: z, state: v} TO ray_path
    
    END WHILE

    // 4. 输出结果
    IF z >= z_end THEN
      PRINT "光线追迹成功完成。"
      PRINT "最终状态在 z = ", z, ":"
      PRINT "  位置 (x, y) = (", v[0], ", ", v[1], ")"
      PRINT "  方向 (Tx, Ty) = (", v[2], ", ", v[3], ")"
      PRINT "  光程 OP = ", v[4]
    ELSE
      PRINT "光线追迹因错误（如全反射）提前终止于 z = ", z
    END IF

END PROCEDURE