第6章：如何计算折射方向
一束光线现在到达了表面。
这句话听起来简单，只是因为我们早已完成了那些隐秘的工作。
在第4章中，我们为光线赋予了位置、方向、波长、强度以及光程。在第5章中，我们让这束光线与平面或球面发生碰撞。我们计算了碰撞点，检查了通光孔径，又计算了局部表面法线。
现在，这束光线就位于表面处。
它有一个入射方向：
d_in
表面有一个局部法线：
N
材料模型给出了两侧的折射率：
n1 = 表面前的折射率
n2 = 表面后的折射率
下一个问题是每个光学学子都期待的：
光线经过折射后，会朝哪个方向前进？
本章就将回答这个问题。
我们会从标量形式的斯涅耳定律开始，但不会止步于此。光学软件通常不会用量角器画角度来追迹光线，它用的是方向向量和表面法线。因此，有用的形式是矢量折射。
这正是许多错误发生的地方。
公式本身并不算长。困难之处在于：
法线方向
单位向量
哪个折射率是 n1，哪个是 n2
全内反射
数值容差
我们会慢慢来。这一章是本书中至关重要的章节之一。
在这一章之后，光线最终会做出镜头图纸上一直展示的那些动作：
撞击一个表面
根据局部法线及折射率发生偏折
继续沿着新的方向前进
只是现在，偏折将通过计算得出。
从标量斯涅耳定律开始
我们熟悉的斯涅耳定律形式是：
n1 sin θ1 = n2 sin θ2
其中：
n1 = 入射介质的折射率
n2 = 透射介质的折射率
θ1 = 从表面法线量起的入射角
θ2 = 从表面法线量起的折射角
这条定律指出，当光线穿过折射率不同的介质时，其方向会发生改变。
如果光从空气进入玻璃：
n1 ≈ 1.0
n2 ≈ 1.5
那么光线会向法线方向偏折。
如果光从玻璃进入空气：
n1 ≈ 1.5
n2 ≈ 1.0
那么光线会偏离法线。
当光线以 θ1=0 垂直入射时，其几何方向保持不变。虽然光在介质中的速度和波长可能会改变，但它的几何方向保持不变。
我们先写一个标量角度函数，仅仅是为了观察定律的数值表现。
import math
def snell_angle(theta1_deg, n1, n2):
"""返回用标量斯涅耳定律算出的折射角（度）。
theta1_deg 是从表面法线量起的角度。
如果发生全内反射，则返回 None。
"""
theta1 = math.radians(theta1_deg)
sin_theta2 = (n1 / n2) * math.sin(theta1)
if abs(sin_theta2) > 1.0:
return None
theta2 = math.asin(sin_theta2)
return math.degrees(theta2)
再试一下空气到玻璃的情况：
for theta in [0, 10, 20, 30, 40]:
theta2 = snell_angle(theta, n1=1.0, n2=1.5)
print(theta, "→", theta2)
折射角更小，光线向法线方向偏折。
再试试玻璃到空气：
for theta in [0, 10, 20, 30, 40, 50]:
theta2 = snell_angle(theta, n1=1.5, n2=1.0)
print(theta, "→", theta2)
在某个入射角，函数会返回 None，这就是全内反射。
标量形式有助于建立直觉，但对于真正的光线追迹来说，它还不够。
光线追迹程序通常不直接存储“与法线的夹角”，它存储的是：
入射方向向量
表面法线向量
所以我们需要矢量版本。
为什么标量角度不够用
想象一束光线偏离光轴射到一个球面上。
局部法线不再是 [0, 0, 1]，它从球心指向碰撞点。它的 x 分量和 y 分量可能不再为零。
在标量斯涅耳定律中，你可以这样写：
入射角是 θ1
折射角是 θ2
但在三维代码中，这并不能告诉你出射向量是什么。除非同时知道入射面，否则有无数个方向能以同样的角度与法线相交。
入射面由以下两者定义：
入射方向
表面法线
折射光线必须位于该平面内。用矢量公式可以自然而然地处理这一点。
所以计算任务就是：
给定：
d_in, normal, n1, n2
计算：
d_out
其中两个方向向量都是单位长度。
这正是我们想要的函数：
d_out = refract_direction(d_in, normal, n1, n2)
但在编写这个函数之前，我们需要处理一个看似琐碎却极其重要的细节。
法线可能有两个方向
表面法线可以指向两个完全相反的方向。
对于一个平面而言，下面这两个法线都是有效的：
[0, 0, 1]
[0, 0, -1]
对于一个球面，其几何法线通常指向球心外侧。那是一个完全合理的几何法线。但对于折射来说，公式要求法线相对于入射光线的方向保持一致。
这正是许多光线追迹程序最先出错的地方。
让我们定义一条约定：
在应用矢量折射公式之前，需要先调整法线的方向，使其与入射光线的方向相反。
换句话说，调整之后应满足：
dot(d_in, normal) <= 0
如果入射光线大致沿着 +z 方向，调整后的法线应大致指向入射侧的后方。
以下是辅助函数：
import numpy as np
def normalize(v):
"""返回向量 v 的单位长度副本。"""
v = np.asarray(v, dtype=float)
norm = np.linalg.norm(v)
if norm == 0:
raise ValueError("无法对零向量进行归一化。")
return v / norm
def orient_normal_against_ray(direction, normal):
"""调整法线方向，使其与入射方向相反。"""
d = normalize(direction)
n = normalize(normal)
if np.dot(d, n) > 0:
n = -n
return n
测试一下：
d_in = np.array([0.0, 0.0, 1.0])
print(orient_normal_against_ray(d_in, [0, 0, 1]))
print(orient_normal_against_ray(d_in, [0, 0, -1]))
两种情形的输出结果都应为：
[ 0. 0. -1.]
这正是我们想要的。折射函数不应关心表面给出的是 [0,0,1] 还是 [0,0,-1]，它应该在使用法线之前先调整其方向。
这个小小的函数消除了大量由符号错误引起的 bug。
矢量折射公式
设：
i = 入射单位方向
n = 调整后与入射方向相反的单位法线
η = n1 / n2
入射角的余弦是：
cosθ1 = - i · n
透射方向可以写作：
t = η i + (η cosθ1 - cosθ2) n
其中：
cosθ2 = sqrt(1 - η²(1 - cos²θ1))
平方根内的项至关重要：
k = 1 - η²(1 - cos²θ1)
如果 k < 0，则不存在实数的折射光线，即发生全内反射。
以上就是完整的计算逻辑。
归一化入射方向
归一化并调整法线方向
计算 η = n1 / n2
计算 cosθ1
计算 k
如果 k < 0：全内反射
否则计算出射方向
归一化出射方向
现在用 Python 来写。
def refract_direction(direction, normal, n1, n2, eps=1e-12):
"""用矢量斯涅耳定律计算折射光线方向。
参数
----------
direction:
入射光线方向。不需要预先归一化。
normal:
碰撞点的表面法线。它可以指向任意方向；
本函数会将其调整为与入射光线相反的方向。
n1:
入射介质的折射率。
n2:
透射介质的折射率。
返回值
-------
单位出射方向向量；如果发生全内反射，则返回 None。
"""
i = normalize(direction)
n = orient_normal_against_ray(i, normal)
eta = n1 / n2
cos_theta1 = -float(np.dot(i, n))
# 数值安全：确保余弦值在有效范围内。
cos_theta1 = max(-1.0, min(1.0, cos_theta1))
sin2_theta2 = eta * eta * (1.0 - cos_theta1 * cos_theta1)
if sin2_theta2 > 1.0 + eps:
return None
# 钳制微小数值溢出。
sin2_theta2 = min(1.0, max(0.0, sin2_theta2))
cos_theta2 = math.sqrt(1.0 - sin2_theta2)
t = eta * i + (eta * cos_theta1 - cos_theta2) * n
return normalize(t)
这个函数就是本章的核心。
它不是一个黑盒，它是斯涅耳定律的直接计算形式。
入射光线方向和表面法线决定了入射面。折射率比值控制着偏折程度。平方根测试则检测全内反射。
现在我们必须仔细测试它。
一个测量角度的辅助函数
为了与标量斯涅耳定律进行比较，我们需要一种方法来测量方向与法线之间的夹角。
法线有两个方向，即 normal 和 -normal，因此我们使用点积的绝对值。
def angle_to_normal_line(direction, normal):
"""返回方向与法线之间较小的夹角。"""
d = normalize(direction)
n = normalize(normal)
c = abs(float(np.dot(d, n)))
c = max(-1.0, min(1.0, c))
return math.degrees(math.acos(c))
现在我们可以将矢量折射与标量斯涅耳定律进行比较测试。
测试1：垂直入射
沿法线传播的光线不应改变方向。
d_in = np.array([0.0, 0.0, 1.0])
normal = np.array([0.0, 0.0, -1.0])
d_out = refract_direction(d_in, normal, n1=1.0, n2=1.5)
print("入射方向:", d_in)
print("出射方向:", d_out)
print("角度:", angle_to_normal_line(d_out, normal))
出射方向仍应为：
[0, 0, 1]
与法线的夹角应为 0。
这是一个基本的完整性检查。如果这里失败了，不要继续，先修正法线方向或公式。
测试2：空气到玻璃
现在使用与法线成30度角的光线。设z轴为法线方向，并让光线在x-z平面内传播。
def direction_from_angle_to_z(theta_deg):
"""在 x-z 平面内、大致沿 +z 方向传播的方向。"""
theta = math.radians(theta_deg)
return normalize([math.sin(theta), 0.0, math.cos(theta)])
d_in = direction_from_angle_to_z(30.0)
normal = np.array([0.0, 0.0, -1.0])
d_out = refract_direction(d_in, normal, n1=1.0, n2=1.5)
print("入射角:", angle_to_normal_line(d_in, normal))
print("出射角:", angle_to_normal_line(d_out, normal))
print("标量斯涅耳:", snell_angle(30.0, n1=1.0, n2=1.5))
print("出射方向:", d_out)
矢量的出射角度应与标量斯涅耳定律相符。
从空气到玻璃会向法线偏折，因此出射角应小于30度。
这个测试证实了在简单平面情况下矢量公式与标量定律是一致的。
测试3：玻璃到空气
现在交换折射率。
d_in = direction_from_angle_to_z(30.0)
normal = np.array([0.0, 0.0, -1.0])
d_out = refract_direction(d_in, normal, n1=1.5, n2=1.0)
print("入射角:", angle_to_normal_line(d_in, normal))
print("出射角:", angle_to_normal_line(d_out, normal))
print("标量斯涅耳:", snell_angle(30.0, n1=1.5, n2=1.0))
print("出射方向:", d_out)
从玻璃到空气会偏离法线，因此出射角应大于30度。
同样，矢量结果应与标量结果相符。
测试4：全内反射
当光试图以过于陡峭的角度从较高折射率介质进入较低折射率介质时，会发生全内反射。
临界角满足：
sin θc = n2 / n1
其中 n1 > n2。
对于玻璃到空气：
θc = arcsin(1.0 / 1.5) ≈ 41.8 度
因此，50度的入射角不应产生透射光线。
d_in = direction_from_angle_to_z(50.0)
normal = np.array([0.0, 0.0, -1.0])
d_out = refract_direction(d_in, normal, n1=1.5, n2=1.0)
print(d_out)
结果应为：
None
这并不代表光线在物理上消失了，而是意味着不存在折射的透射光线。此时界面会对光线进行内部反射。在折射透镜追迹路径中，根据系统模型的不同，我们可以终止该光线，也可以处理反射。
对于本书中基础的折射光线追迹器，我们将会终止该光线并记录原因。
全内反射并非错误
很容易将 None 当作代码的失败，但事实并非如此。
全内反射是一种真实的物理状态。代码通过检测到不存在折射方向来完成其职责。
光线追迹器接下来做什么，取决于光学模型。
对于纯折射系统，可以将该光线标记为无效：
终止光线：全内反射
对于支持反射的系统，则可以改为计算反射方向：
d_reflect = d_in - 2(d_in · n)n
在此我们不会构建反射式光线追迹器，但值得了解的是，全内反射是一个分支事件，而非数值意外。
在我们基础透镜光线追迹器中，首先停止追迹光线是最合适的行为。
用于对比的反射函数
尽管本章的重点是折射，但反射的公式足够简短，值得展示一下。它也有助于理解法线方向的问题。
def reflect_direction(direction, normal):
"""将方向关于表面法线进行反射。"""
i = normalize(direction)
n = normalize(normal)
return normalize(i - 2.0 * np.dot(i, n) * n)
如果发生全内反射，更完整的追迹器可以使用此函数。目前，当折射无法进行时，我们仅会终止光线的折射追迹。
本书的主线是透镜成像，而非反射镜系统。
将波长重新纳入计算
到目前为止，我们使用的都是固定值，例如：
n1 = 1.0
n2 = 1.5
但在第3章中，我们已经知道真实材料的折射率是随波长变化的：
n = n(λ)
光线携带了波长信息：
ray.wavelength_um
因此，在每个表面上，追迹代码都应该向材料库查询：
n1 = materials.n(material_before, ray.wavelength_um)
n2 = materials.n(material_after, ray.wavelength_um)
这就是色散折射进入光线追迹器的方式。
一束蓝光和一束红光可以有相同的入射方向，并撞击到相同的表面点。但由于玻璃在不同波长下的折射率不同，它们的出射方向会略有差异。
让我们用一个平面界面来演示一下。
假设我们有第3章中的 N-BK7 材料模型：
# N_BK7.n(wavelength_um) 返回折射率。
# 在本示例中，假设 N_BK7 已经定义好了。
现在对比蓝光、绿光和红光以30度角进入玻璃的情况。
F_LINE_UM = 0.4861327
D_LINE_UM = 0.5875618
C_LINE_UM = 0.6562725
normal = np.array([0.0, 0.0, -1.0])
d_in = direction_from_angle_to_z(30.0)
for wl in [F_LINE_UM, D_LINE_UM, C_LINE_UM]:
n_glass = N_BK7.n(wl)
d_out = refract_direction(d_in, normal, n1=1.0, n2=n_glass)
angle = angle_to_normal_line(d_out, normal)
print(f"{wl:.7f} µm n={n_glass:.7f} 出射角={angle:.6f} deg")
虽然差异很小，但却真实存在。
在正常色散材料中，较短波长通常对应更高的折射率，因此蓝光会比红光向法线偏折得更厉害一些。
这就是通过一行代码体现出的色散行为：
n_glass = N_BK7.n(wl)
这就是第3章不可或缺的原因。如果没有与波长相关的材料折射率，光线追迹器就无法展现色差。
曲面上的折射
现在，让我们使用第5章所给出的内容：一个球面上的碰撞点和局部法线。
想象一束平行光线射入一个 N-BK7 球面透镜的前表面。
我们已经有了：
光线
表面
intersect_surface(ray, surface)
碰撞记录中包含：
hit.point
hit.normal
现在我们计算出射方向。
surface = IntersectableSurface(
radius=50.0,
z_vertex=0.0,
semi_diameter=10.0,
comment="前表面球面"
)
ray = Ray(
position=[5.0, 0.0, -20.0],
direction=[0.0, 0.0, 1.0],
wavelength_um=0.5875618
)
hit = intersect_surface(ray, surface)
n_air = 1.0
n_glass = N_BK7.n(ray.wavelength_um)
d_out = refract_direction(
direction=ray.direction,
normal=hit.normal,
n1=n_air,
n2=n_glass
)
print("碰撞点:", hit.point)
print("表面法线:", hit.normal)
print("入射方向:", ray.direction)
print("出射方向:", d_out)
对于离轴光线，表面法线具有 x 分量。出射方向也将获得 x 分量。这就是表面在聚焦光线。
这正是透镜发挥作用的局部原因。
镜头并不知道像面在哪里。在每个表面上，光线只看到：
局部法线
入射方向
折射率 n1
折射率 n2
整体的聚焦行为是由许多局部折射和传播汇聚而成的。
这是一个强有力的理念，也正是光线追迹软件所计算的内容。
在表面上更新光线
折射函数返回一个方向，但我们的光线对象必须相应地更新。
一次表面交互包括以下步骤：
1. 找到表面碰撞点
2. 将光线通过当前介质传播到碰撞点
3. 检查孔径
4. 计算折射率
5. 计算折射方向
6. 更新光线方向
7. 继续进入下一个介质
第3、4、5章几乎提供了所有所需的部件，现在我们将它们组合起来。
首先，定义一个小的结果对象：
from dataclasses import dataclass
@dataclass
class RefractionResult:
success: bool
direction: np.ndarray | None
reason: str = ""
如果你的 Python 版本不支持 np.ndarray | None，可以从 typing 导入 Optional[np.ndarray]。
现在包装折射调用：
def refract_at_surface(ray, hit, n1, n2):
"""返回光线在表面碰撞处的折射结果。"""
d_out = refract_direction(
direction=ray.direction,
normal=hit.normal,
n1=n1,
n2=n2
)
if d_out is None:
return RefractionResult(
success=False,
direction=None,
reason="total internal reflection"
)
return RefractionResult(
success=True,
direction=d_out,
reason=""
)
然后，一次表面交互就可以更新光线：
def interact_refractive_surface(ray, surface, n1, n2):
"""将光线传播到一个表面并进行折射。
此函数假设 n1 是表面前的介质折射率，
n2 是表面后的介质折射率。
"""
if not ray.alive:
return None
hit = intersect_surface(ray, surface)
if hit is None:
ray.stop(f"missed surface: {surface.comment}")
return None
# 在入射介质中传播至碰撞点。
ray.propagate(hit.t, refractive_index=n1)
if not hit.inside_aperture:
ray.stop(f"blocked by aperture: {surface.comment}")
return hit
result = refract_at_surface(ray, hit, n1=n1, n2=n2)
if not result.success:
ray.stop(result.reason)
return hit
old_direction = ray.direction.copy()
ray.direction = result.direction
ray.history.append({
"event": "refract",
"surface": surface.comment,
"point": ray.position.copy(),
"normal": hit.normal.copy(),
"n1": n1,
"n2": n2,
"direction_before": old_direction,
"direction_after": ray.direction.copy(),
})
return hit
这个函数是一个里程碑。
现在光线可以：
找到一个表面
移动到该表面
检查孔径
偏折进入下一介质
以新方向继续前进
这就是真正的光线追迹。
虽然仍然简单，但确实是真实的。
穿过单个球面界面的追迹
让我们追迹几束平行光线，观察它们如何穿过一个球面空气-玻璃界面的。
surface = IntersectableSurface(
radius=50.0,
z_vertex=0.0,
semi_diameter=10.0,
comment="air-to-glass spherical surface"
)
rays = [
Ray(position=[x, 0.0, -20.0], direction=[0.0, 0.0, 1.0], wavelength_um=0.5875618)
for x in np.linspace(-8.0, 8.0, 5)
]
n_air = 1.0
n_glass = N_BK7.n(0.5875618)
for ray in rays:
interact_refractive_surface(ray, surface, n1=n_air, n2=n_glass)
print(
"位置:", ray.position,
"方向:", ray.direction,
"存活:", ray.alive
)
轴上光线应保持在 z 方向上。离轴光线应向内侧偏折。相对两侧的光线应呈对称偏折。
这种对称性是一项有用的健全性检查。如果两侧的离轴光线都向横向的同一方向偏折，那一定是哪里出了问题。
最可能的原因包括：
法线符号错误
球面中心符号错误
n1/n2 顺序错误
交点分支选取错误
由于我们保持了各函数间的独立性，调试路线已经十分清晰。
可视化偏折
让我们绘制围绕球面的入射与出射线段。
import matplotlib.pyplot as plt
def line_points(start, direction, length):
start = np.asarray(start, dtype=float)
direction = normalize(direction)
end = start + length * direction
return start, end
surface = IntersectableSurface(
radius=50.0,
z_vertex=0.0,
semi_diameter=10.0,
comment="前表面球面"
)
rays = [
Ray(position=[x, 0.0, -20.0], direction=[0.0, 0.0, 1.0], wavelength_um=0.5875618)
for x in np.linspace(-8.0, 8.0, 7)
]
plt.figure()
# 在 x-z 平面中绘制表面轮廓。
xs = np.linspace(-10.0, 10.0, 300)
zs = [spherical_sag(x, 0.0, surface) for x in xs]
plt.plot(zs, xs, label="spherical surface")
for ray in rays:
start = ray.position.copy()
incoming_dir = ray.direction.copy()
hit = interact_refractive_surface(ray, surface, n1=1.0, n2=N_BK7.n(ray.wavelength_um))
if hit is None:
continue
hit_point = ray.position.copy()
outgoing_dir = ray.direction.copy()
# 入射线段。
plt.plot([start[2], hit_point[2]], [start[0], hit_point[0]], linewidth=0.8)
# 出射线段的预览。
end = hit_point + 25.0 * outgoing_dir
plt.plot([hit_point[2], end[2]], [hit_point[0], end[0]], linewidth=0.8)
plt.xlabel("z [mm]")
plt.ylabel("x [mm]")
plt.title("球面空气-玻璃界面的折射")
plt.gca().set_aspect("equal", adjustable="box")
plt.legend()
plt.show()
这幅图是真正计算得来的。
它并非一张装饰性的草图，而是对计算出的交点和计算出的折射方向的可视化呈现。
这正是我们在全书中希望推动的那种转变。
n1 与 n2 的顺序至关重要
在透镜的前表面：
空气 → 玻璃
因此：
n1 = n_air
n2 = n_glass
在透镜的后表面：
玻璃 → 空气
因此：
n1 = n_glass
n2 = n_air
如果你把它们弄反了，光线就会向错误的方向偏折。
这就是第3章中材料过渡逻辑的重要性所在。追迹代码不应随意猜测。
像下面这样的函数会很有帮助：
def material_before_after(lens, surface_index, starting_material="air"):
"""返回顺序表面之前和之后的材料名称。"""
if surface_index == 0:
before = starting_material
else:
before = lens.surfaces[surface_index - 1].material_after
after = lens.surfaces[surface_index].material_after
return before, after
def index_transition(lens, materials, surface_index, wavelength_um):
before_name, after_name = material_before_after(lens, surface_index)
n_before = materials.n(before_name, wavelength_um)
n_after = materials.n(after_name, wavelength_um)
return n_before, n_after
然后在某个表面上：
n1, n2 = index_transition(lens, materials, surface_index, ray.wavelength_um)
这一行代码将正确的折射率顺序编码了进去。
没有这套规范，多表面追迹就会变得很脆弱。
将交点求取与材料过渡结合
在顺序系统中，一个完整的表面步骤需要同时用到几何数据和材料数据。
表面几何告诉我们：
光线在何处碰撞
法线是什么
光线是否在孔径内
而透镜规格和材料库告诉我们：
n_before
n_after
折射公式则告诉我们：
新的方向
让我们勾勒一下组合循环。
def trace_refractive_surfaces(ray, intersectable_surfaces, lens, materials, start_index=1):
"""将一束光线追迹通过多个折射表面。
这是一个简化的骨架。它假设 intersectable_surfaces 与
经过对象行处理后的 lens.surfaces 对应相同的顺序表面。
"""
for local_i, surface in enumerate(intersectable_surfaces):
if not ray.alive:
break
# 跳过对象面（如果存在）。
if surface.comment.lower().startswith("object"):
continue
# 像面是一个评估面，而非折射边界。
if surface.comment.lower().startswith("image"):
hit = intersect_surface(ray, surface)
if hit is None:
ray.stop("missed image plane")
break
ray.propagate(hit.t, refractive_index=1.0)
break
prescription_index = start_index + local_i
n1, n2 = index_transition(
lens,
materials,
prescription_index,
ray.wavelength_um
)
interact_refractive_surface(ray, surface, n1=n1, n2=n2)
return ray
这仍然是简化版。它假设了一个直接的顺序系统。它不处理坐标转折、反射镜、反射镀膜、浸没介质、倾斜表面或更复杂的孔径定义。
这样就可以了。
在同一个地方集中表达了计算链条：
光线波长
→ 材料折射率
→ 表面碰撞点
→ 法线
→ 矢量斯涅耳折射
→ 更新后的光线方向
这就是真实光线追迹的核心。
双面单透镜追迹
现在我们来做一个基本的双凸单透镜。
我们将追迹光线通过：
前表面：空气 → N-BK7
后表面：N-BK7 → 空气
像面： 仅用于评估
使用前面章节中的表面：
front = IntersectableSurface(
radius=50.0,
z_vertex=0.0,
semi_diameter=10.0,
comment="前表面"
)
back = IntersectableSurface(
radius=-50.0,
z_vertex=5.0,
semi_diameter=10.0,
comment="后表面"
)
image = IntersectableSurface(
radius=float("inf"),
z_vertex=50.0,
semi_diameter=None,
surface_type="plane",
comment="像面"
)
追迹一束光线：
ray = Ray(
position=[5.0, 0.0, -20.0],
direction=[0.0, 0.0, 1.0],
wavelength_um=0.5875618
)
n_air = 1.0
n_glass = N_BK7.n(ray.wavelength_um)
interact_refractive_surface(ray, front, n1=n_air, n2=n_glass)
interact_refractive_surface(ray, back, n1=n_glass, n2=n_air)
# 在空气中传播到像面。
hit_image = intersect_surface(ray, image)
ray.propagate(hit_image.t, refractive_index=n_air)
print("像点:", ray.position)
print("最终方向:", ray.direction)
print("OPL:", ray.opl)
print("存活:", ray.alive)
print("终止原因:", ray.stop_reason)
这是首次完成一个简单透镜的完整折射追迹。
它还不是一个成熟的光学设计程序。它没有从适当的物面视场发射光线，也没有计算最佳焦点，也不处理所有边缘情况。但它执行了核心操作：
空气传播
前表面交点求取
空气到玻璃的折射
玻璃传播
后表面交点求取
玻璃到空气的折射
空气传播至像面
这正是许多光学设计工具底层的骨架。
从一束光线到一个小光扇
让我们追迹一个平行光扇穿过单透镜，并标出它们的落点。
def trace_simple_singlet_ray(ray, front, back, image, glass_material):
n_air = 1.0
n_glass = glass_material.n(ray.wavelength_um)
interact_refractive_surface(ray, front, n1=n_air, n2=n_glass)
if ray.alive:
interact_refractive_surface(ray, back, n1=n_glass, n2=n_air)
if ray.alive:
hit_image = intersect_surface(ray, image)
if hit_image is None:
ray.stop("missed image plane")
else:
ray.propagate(hit_image.t, refractive_index=n_air)
return ray
现在追迹若干条光线。
rays = [
Ray(position=[x, 0.0, -20.0], direction=[0.0, 0.0, 1.0], wavelength_um=0.5875618)
for x in np.linspace(-8.0, 8.0, 9)
]
traced = [
trace_simple_singlet_ray(ray, front, back, image, N_BK7)
for ray in rays
]
for ray in traced:
print(ray.position, ray.alive, ray.stop_reason)
如果像面没有恰好位于最佳焦点处，光线可能不会汇聚到一个点上。这很正常。事实上，这正是图像分析的起点。
一个真正的镜头设计程序会帮助寻找焦点，计算近轴量，并在选定的视场和波长下生成点列图。我们还没有到达那一步。但现在光线确实已经穿过了折射表面。
追迹机制已经运转起来了。
绘制追迹出的光线
由于我们的光线存储了历史记录，我们可以画出它们的路径。
def path_points_from_ray(ray):
points = []
for event in ray.history:
if event["event"] == "propagate":
if not points:
points.append(event["from"])
points.append(event["to"])
if len(points) == 0:
return np.empty((0, 3))
return np.array(points)
plt.figure()
# 绘制大致的表面轮廓。
for s in [front, back]:
xs = np.linspace(-s.semi_diameter, s.semi_diameter, 300)
zs = [spherical_sag(x, 0.0, s) for x in xs]
plt.plot(zs, xs, color="black", linewidth=1.0)
# 绘制像面。
plt.axvline(image.z_vertex, linestyle="--", linewidth=1.0)
# 绘制光线路径。
for ray in traced:
path = path_points_from_ray(ray)
if len(path) > 0:
plt.plot(path[:, 2], path[:, 0], linewidth=0.8)
plt.xlabel("z [mm]")
plt.ylabel("x [mm]")
plt.title("穿过折射单透镜的小光扇")
plt.gca().set_aspect("equal", adjustable="box")
plt.show()
这幅图依然简单，但它已经跨越了一条重要的界限。
光线在表面处偏折，是因为代码计算了以下内容：
表面碰撞点
局部法线
n1 和 n2
矢量斯涅耳方向
这已经不仅仅是几何学了，它已经是折射光线追迹了。
为何不同高度的光线偏折不同
一束平行的轴上光线会以不同的高度打到透镜上。
在球面的中心处，法线与光轴对齐。离轴越远，法线就越倾斜。
斯涅耳定律使用的是入射光线与局部法线之间的夹角。因此，不同的光线高度会对应不同的入射角。
这意味着球面并不会将所有平行光线偏折完全相同的量。这就是球差存在的原因之一。我们稍后会研究它，但目前这一原因已可见端倪。
在每一个光线高度：
相同的入射方向
不同的表面法线
不同的折射方向
这就是边缘光线和近轴光线可能不会汇聚在同一焦点上的计算层面原因。
透镜并不是以某个整体动作去“聚焦光线”的，它是在整个表面上局部地改变光线的方向。而像面正是这些局部作用的集体结果。
这是一种更好的思维模型。
为何不同波长的光线偏折不同
现在追迹同一光线在三个不同波长下的情形。
for wl in [F_LINE_UM, D_LINE_UM, C_LINE_UM]:
ray = Ray(
position=[5.0, 0.0, -20.0],
direction=[0.0, 0.0, 1.0],
wavelength_um=wl
)
trace_simple_singlet_ray(ray, front, back, image, N_BK7)
print(
f"λ={wl:.7f} µm "
f"像面 x={ray.position[0]:.6f} "
f"z={ray.position[2]:.3f} "
f"存活={ray.alive}"
)
像面 x 位置可能略有不同，最终方向也可能不同。
这一微小的差异就是实际追迹中色差产生的根源。
再次强调，这并不神秘。它源于：
N_BK7.n(wavelength)
这改变了 n1/n2，进而改变了折射方向。
整个链条都是明确的。
Optiland 在高层次上做了什么
此时，将我们的教学代码与真实库的工作流程进行比较会很有帮助。
在我们的代码中，我们手动编写了：
求表面交点
计算法线
查找折射率
应用矢量斯涅耳定律
更新光线
在成熟的光学设计库中，这些步骤都被封装在追迹方法之中。用户可能只需几行代码来定义镜头并绘制追迹出的光线。
一个类似 Optiland 风格的单透镜工作流程可能如下所示：
from optiland import optic
lens = optic.Optic()
lens.surfaces.add(index=0, radius=float("inf"), thickness=float("inf"))
lens.surfaces.add(index=1, radius=50.0, thickness=5.0, material="N-BK7", is_stop=True)
lens.surfaces.add(index=2, radius=-50.0, thickness=45.0)
lens.surfaces.add(index=3)
lens.set_aperture(aperture_type="EPD", value=10.0)
lens.fields.set_type("angle")
lens.fields.add(y=0.0)
lens.wavelengths.add(value=0.5876, is_primary=True)
lens.draw(num_rays=7)
根据所安装的 Optiland 版本，具体的光线访问方法可能会有所不同。这里进行比较的目的并不是记忆某个特定的 API 调用，而是为了认识到，当库绘制折射光线时，它在底层必须完成哪些工作。
对于每束光线以及每个折射表面，都会发生某种形式的以下过程：
找到碰撞点
找到法线
计算折射率
应用斯涅耳定律
继续追迹
库可能对其进行了向量化处理，可能使用了后端的数组系统，可能支持更多种类的表面类型，也可能包含了鲁棒的光线瞄准和瞳孔管理功能，但它无法跳过物理定律。
软件的命令可以隐藏这些步骤。
但它无法废除这些步骤。
如何将我们的结果与库的结果对比
一次良好的对比要保持适度。
不要一开始就去对比整个 MTF 曲线，那中间隔着太多环节。先对比一个表面或一个简单的单透镜。
一个实用的对比计划是：
1. 在教学代码和 Optiland 中构建同一个简单透镜。
2. 使用相同的波长。
3. 如果 API 暴露了光线层级的追迹，使用相同的入射光线高度和方向。
4. 比较经过第一表面之后的光线位置和方向。
5. 然后比较经过第二表面之后的情况。
6. 只有在这之后，才去比较像面的截距。
如果库不能方便地暴露同层级的光线数据，可以先比较布局图或最终光线截距。但概念上，最好的对比方式是逐表面进行。
当结果出现差异时，不要立刻认为库有误或自己的代码有误。首先检查约定：
表面半径符号
表面顶点位置
材料分配
波长单位
法线方向
孔径定义
物面光线设置
像面位置
大多数早期的差异都来自约定不匹配，而非物理错误。
这就是我们保持代码显式的原因。它让我们有东西可以检查。
一个小的诊断打印输出
当折射后的光线看起来不对劲时，可以在表面处打印出相关数值。
def debug_refraction(ray, surface, n1, n2):
hit = intersect_surface(ray, surface)
if hit is None:
print("无碰撞点。")
return
d_out = refract_direction(ray.direction, hit.normal, n1, n2)
print("表面:", surface.comment)
print("碰撞点:", hit.point)
print("入射方向:", normalize(ray.direction))
print("几何法线:", normalize(hit.normal))
print("调整后法线:", orient_normal_against_ray(ray.direction, hit.normal))
print("n1, n2:", n1, n2)
print("入射角:", angle_to_normal_line(ray.direction, hit.normal))
if d_out is None:
print("全内反射")
else:
print("出射方向:", d_out)
print("出射角:", angle_to_normal_line(d_out, hit.normal))
使用它：
ray = Ray(position=[5.0, 0.0, -20.0], direction=[0.0, 0.0, 1.0])
debug_refraction(ray, front, n1=1.0, n2=N_BK7.n(0.5875618))
这种诊断输出虽然不那么光鲜，但正是理解光线追迹器的途径。
如果光线向错误的方向偏折，诊断输出通常能揭示原因。
可能是法线方向反了。
可能是 n1 和 n2 被交换了。
可能是碰撞点与你想象的不一样。
可能是光线打到了球面的错误分支上。
数据会讲述出真相。
矢量折射中的常见错误
让我们清楚地列举一下这些常见错误。
错误1：未对向量进行归一化
公式要求使用单位方向向量和单位法线。
如果 direction 或 normal 没有归一化，余弦计算就会出错。务必归一化它们。
错误2：在未调整方向的情况下直接使用法线
几何法线可能指向任意一侧。在应用公式之前，需调整其方向，使其与入射方向相反。
应满足如下条件：
dot(direction, adjusted_normal) <= 0
错误3：交换 n1 和 n2
在空气到玻璃的界面，应使用：
n1 = air
n2 = glass
在玻璃到空气的界面，应使用：
n1 = glass
n2 = air
若交换二者，光线会向错误方向偏折。
错误4：将全内反射视为崩溃
当平方根项变为负数时，物理上意味着不存在折射光线。需明确处理这种情况。
错误5：从表面而不是法线量测角度
斯涅耳定律使用的是与法线之间的夹角，而不是与表面切平面的夹角。
这是最容易出现的概念性错误之一。
错误6：忘记波长
对于真实玻璃，n2 取决于波长。如果所有波长都使用相同的折射率，模型中的色差效应就会消失。
错误7：在检查孔径之前进行折射
如果光线落在通光孔径之外，它应被阻挡。不要将其视为通过玻璃而进行折射。
错误8：在错误的时刻增加光程
光线在介质中传播时会累积光程。在理想零厚度的界面处进行折射时，改变的是方向，而非几何路程长度。
在我们的简化模型中：
传播到表面 → 更新光程
在表面折射 → 更新方向
错误9：忘记数值容差
在临界角附近，极小的浮点误差就可能决定计算出的 sin²θ2 是略低于还是略高于 1。需使用一个较小的容差。
矢量公式的物理意义
人们很容易把矢量公式当作一段可以复制粘贴的代码，让我们把它与物理联系起来。
相对于表面法线，入射方向可以分解为两个部分：
切向分量
法向分量
在理想折射界面上，切向分量受到斯涅耳定律的约束，而法向分量会发生变化，以使出射方向为单位长度并位于透射一侧。
公式：
t = η i + (η cosθ1 - cosθ2) n
恰恰完成了这一点。它根据折射率比值缩放切向行为，然后选择合适的法向分量。
这就是出射光线保持在入射面内的原因。
这也正是当切向要求将迫使出射方向具有不可能的长度时，全内反射就会出现的原因。
因此，这段代码并非魔法，它是以矢量形式写出的斯涅耳定律的几何体现。
一个完整的单面教学示例
以下是整章核心计算的一个紧凑版本。
# 一束光线在空气中射向一个球面 N-BK7 表面。
ray = Ray(
position=[5.0, 0.0, -20.0],
direction=[0.0, 0.0, 1.0],
wavelength_um=0.5875618
)
surface = IntersectableSurface(
radius=50.0,
z_vertex=0.0,
semi_diameter=10.0,
comment="前表面球面"
)
# 1. 找到碰撞点。
hit = intersect_surface(ray, surface)
if hit is None:
ray.stop("missed surface")
elif not hit.inside_aperture:
ray.stop("blocked by aperture")
else:
# 2. 在空气中移动到碰撞点。
ray.propagate(hit.t, refractive_index=1.0)
# 3. 计算折射率。
n1 = 1.0
n2 = N_BK7.n(ray.wavelength_um)
# 4. 折射。
d_out = refract_direction(ray.direction, hit.normal, n1, n2)
if d_out is None:
ray.stop("total internal reflection")
else:
ray.direction = d_out
print("位置:", ray.position)
print("方向:", ray.direction)
print("OPL:", ray.opl)
print("存活:", ray.alive)
print("原因:", ray.stop_reason)
这就是折射事件的全部，完全暴露在外。
没有按钮，没有隐藏的光线追迹，没有绘图中难以解释的偏折。
只有计算本身。
本章如何改变了整本书
本章之前，我们拥有的是零散的部件：
透镜规格
材料模型
光线对象
表面交点求取
表面法线
本章之后，这些部件协同作用起来。
一束光线现在能够使用物理折射穿过表面了。
这意味着我们可以开始产出几何光学的首批真正成果：
穿过透镜的光线路径
像面截距
点列图
边缘光线与主光线的行为
实际光线偏差
色散焦点漂移
这些还不是最终的像质评价指标，但它们是基础。
下游的一切都依赖这一步骤。
点列图是由许多折射光线构成的。
光扇图也是通过比较不同瞳孔坐标下的折射光线而构建的。
波前计算需要光线路径和光程长度。
PSF 和 MTF 最终依赖于波前信息。
因此，这一束光线在这一表面上的偏折并不是一个小事件。它是链条中的一环，这个链条一直通回到第 1 章提出那个问题：
MTF 曲线到底是怎么来的？
答案的一部分是：
它来自于许多由矢量斯涅耳定律计算出的局部折射。
下一步是什么
一束折射光线固然有用，但光学分析通常不会就此止步。
一个透镜不会仅凭一条光线就形成图像。我们需要一束光线：
遍布瞳孔的光线
来自不同视场点的光线
不同波长的光线
当这些光线到达像面时，它们的截距就成了数据。标出这些截距，你就得到了点列图。比较不同瞳孔坐标，你就能开始看出像差的结构。追迹不同波长，颜色误差就会显现出来。
这就是第 7 章的内容。
我们将有意识地生成光线束，而不是让软件包替我们隐式地选择它们。我们将对瞳孔进行采样，赋予视场方向，将光线追迹穿过单透镜，并标出它们的落点。
到那时，一幅熟悉的光学分析图形就不再只是一张图片了。
它将变成一系列计算出的光线事件的集合。
生成的验证图
