4. 自由曲面光线追迹基础

任务与实现边界

本案例实现一个最小可运行的自由曲面单面折射光线追迹器：给定入射光线（起点、方向）、自由曲面（用高度函数 z=f(x,y) 表示）和折射率，计算光线与曲面的交点、曲面法线，并按斯涅尔定律的矢量形式求折射方向。曲面用旋转对称的偶次非球面作为自由曲面的最小实例，便于解析求法线并验证。

核心原理

自由曲面法线由隐函数 F(x,y,z)=z-f(x,y)=0 的梯度给出：N=(-∂f/∂x, -∂f/∂y, 1)，归一化后使用。斯涅尔定律矢量形式：

t = (n₁/n₂)·i + [ (n₁/n₂)·cosθᵢ − cosθₜ ]·N，其中 cosθᵢ = −i·N（i 为入射单位方向，N 指向入射侧），cosθₜ = √(1−(n₁/n₂)²(1−cos²θᵢ))。当根号内为负时发生全反射。

完整代码

import numpy as np

# 自由曲面：旋转对称偶次非球面 z = f(x,y)
# 参数：曲率 c、圆锥常数 k、高阶系数 A4,A6
def surface_f(x, y, c=0.02, k=-0.5, A4=1e-5, A6=1e-7):
    r2 = x*x + y*y
    denom = 1.0 - (1.0 + k) * c * c * r2
    z = c * r2 / (1.0 + np.sqrt(np.maximum(denom, 1e-12)))
    z += A4 * r2**2 + A6 * r2**3
    return z

def surface_grad(x, y, c=0.02, k=-0.5, A4=1e-5, A6=1e-7, h=1e-6):
    # 数值梯度，通用且对任意自由曲面有效
    fx = (surface_f(x+h, y, c, k, A4, A6) - surface_f(x-h, y, c, k, A4, A6)) / (2*h)
    fy = (surface_f(x, y+h, c, k, A4, A6) - surface_f(x, y-h, c, k, A4, A6)) / (2*h)
    return fx, fy

def intersect(P0, d, c=0.02, k=-0.5, A4=1e-5, A6=1e-7, t0=0.0, tol=1e-10, itmax=60):
    # 求光线 P0 + t*d 与 z=f(x,y) 的交点，牛顿迭代
    t = t0
    for _ in range(itmax):
        P = P0 + t * d
        F = P[2] - surface_f(P[0], P[1], c, k, A4, A6)
        fx, fy = surface_grad(P[0], P[1], c, k, A4, A6)
        dF = -fx * d[0] - fy * d[1] + d[2]
        if abs(dF) < 1e-14:
            break
        dt = -F / dF
        t += dt
        if abs(dt) < tol:
            break
    return P0 + t * d, t

def refract(i, N, n1, n2):
    # 斯涅尔定律矢量形式；i 为入射单位方向，N 为指向入射侧的单位法线
    i = i / np.linalg.norm(i)
    N = N / np.linalg.norm(N)
    cos_i = -np.dot(i, N)
    if cos_i < 0:
        N = -N
        cos_i = -np.dot(i, N)
    eta = n1 / n2
    sin2_t = eta * eta * (1.0 - cos_i * cos_i)
    if sin2_t > 1.0:
        return None  # 全反射
    cos_t = np.sqrt(1.0 - sin2_t)
    t = eta * i + (eta * cos_i - cos_t) * N
    return t / np.linalg.norm(t)

def trace(P0, d, n1=1.0, n2=1.5, **kw):
    P0 = np.asarray(P0, float)
    d = np.asarray(d, float)
    d = d / np.linalg.norm(d)
    P, t = intersect(P0, d, **kw)
    fx, fy = surface_grad(P[0], P[1], **kw)
    N = np.array([-fx, -fy, 1.0])
    N = N / np.linalg.norm(N)
    if np.dot(d, N) > 0:
        N = -N  # 保证 N 指向入射侧
    t_dir = refract(d, N, n1, n2)
    return P, N, t_dir

if __name__ == "__main__":
    P0 = np.array([0.0, 0.0, -10.0])
    d = np.array([0.0, 0.0, 1.0])
    P, N, t_dir = trace(P0, d, n1=1.0, n2=1.5)
    print("交点 P =", np.round(P, 6))
    print("法线 N =", np.round(N, 6))
    print("折射方向 t =", np.round(t_dir, 6))
    print("|t| =", round(float(np.linalg.norm(t_dir)), 10))
    # 验证斯涅尔定律：n1*sin(θi) = n2*sin(θt)
    i = d / np.linalg.norm(d)
    cos_i = abs(np.dot(i, N))
    sin_i = np.sqrt(1 - cos_i**2)
    cos_t = abs(np.dot(t_dir, N))
    sin_t = np.sqrt(1 - cos_t**2)
    print("n1*sin_i =", round(1.0*sin_i, 8), " n2*sin_t =", round(1.5*sin_t, 8))

运行与验证

环境：Python 3 + NumPy（pip install numpy）。直接运行脚本：python trace.py。

通过条件：

• 输出 |t| 等于 1.0（方向向量已归一化）；
• n1*sin_i 与 n2*sin_t 数值相等（斯涅尔定律成立）；
• 对轴上入射光线（P0=(0,0,-10)，d=(0,0,1)），交点应为曲面顶点 (0,0,0)，法线为 (0,0,1)，折射方向仍为 (0,0,1)。

可进一步验证：把入射方向改为斜入射（如 d=(0.1,0,1)），观察折射方向偏离法线的角度随 n2 增大而增大；把 n1 设为 1.5、n2 设为 1.0 并增大入射角，可触发全反射（函数返回 None）。