做光学仿真这几年最常被问到的一个问题是光从空气进水里到底偏转多少能量又分走多少教科书上两个公式就能回答但真要在电脑里把光的路径一步步画出来还得靠代码落地。这个基于MATLAB的斯奈尔定律与菲涅尔方程光线追踪模拟就是解决这个问题的——把几何光学的两大核心定律翻译成可执行的程序让界面折射、全反射、能量分配这些抽象概念变成屏幕上直观的光路图。它适合三类人一是正在学物理光学、想把公式跑起来的学生二是做光学系统设计、需要在方案阶段快速验证光路是否合理的工程师三是对光线追踪感兴趣的MATLAB爱好者。有了这套代码你能模拟一束光经过平面界面、棱镜甚至多层介质后的传播路径还能看到每个界面上反射光和透射光各自分到了多少能量——这对理解透镜系统、镀膜设计、光纤通信的原理都很有价值。1. 整体设计与物理原理拆解1.1 斯奈尔定律光的拐弯规则光从一种透明介质进入另一种透明介质时传播方向会发生变化这就是折射。斯奈尔定律定量描述了这种变化n₁·sinθ₁ n₂·sinθ₂其中n₁、n₂分别是入射介质和折射介质的折射率θ₁是入射角θ₂是折射角。角度都是从界面法线量起的不是从界面本身。这个公式看起来简单但实现时有一个容易被忽略的关键点入射角必须换算为弧度制而且正弦函数在角度接近90度时变化非常缓慢。如果在代码里直接用角度制计算或者把法线方向弄反了算出来的折射方向就会完全错误。反射定律同样重要反射角等于入射角。在代码里实现反射方向计算时不需要查表直接用向量运算就能搞定——入射方向向量减去两倍的法线方向投影就是反射方向。这个向量公式比角度计算更稳定也更适合编程实现。1.2 菲涅尔方程界面的能量分配规则斯奈尔定律告诉光往哪儿走菲涅尔方程则告诉光分多少能量走。光打到界面上时一部分能量被反射一部分能量透射过去这个比例由菲涅尔方程给出。对于非偏振光自然光反射率R取s偏振和p偏振的平均值Rs |(n₁·cosθ₁ - n₂·cosθ₂) / (n₁·cosθ₁ n₂·cosθ₂)|²Rp |(n₁·cosθ₂ - n₂·cosθ₁) / (n₁·cosθ₂ n₂·cosθ₁)|²R (Rs Rp) / 2透射率T 1 - R在无吸收的理想界面上能量守恒。实现菲涅尔方程时有个细节值得注意当入射角为0度垂直入射时Rs和Rp的计算公式会退化此时R |(n₁ - n₂)/(n₁ n₂)|²。而入射角接近临界角时反射率会急剧上升——虽然还没有全反射但能量已经开始大量回流这个效应在仿真中非常明显也是很多初学者容易遗漏的物理细节。1.3 全反射光线追踪的边界case当光从光密介质折射率大射向光疏介质折射率小时如果入射角超过某个阈值折射光会消失所有能量都被反射回来。这个阈值叫临界角θc arcsin(n₂/n₁)在代码实现中临界角判断必须在计算折射方向之前完成。我的习惯是先用n₁·sinθ₁ n₂来判断是否全反射——这个判断比先算arcsin再比较角度更高效而且避免了反三角函数的数值误差。全反射在现实世界中应用极广光纤通信就是利用全反射把光锁在纤芯里传播钻石的闪耀也部分归功于全反射——高折射率意味着小临界角光更容易在内部反弹。1.4 为什么用MATLAB做这件事选MATLAB做光学追踪不是因为它在光学计算上有什么不可替代的优势而是因为它刚刚好。第一MATLAB的向量化运算让光线逐条并行处理变得很自然写循环也不用担心性能太差第二它的绘图系统开箱即用画光线、画界面、画法线都很方便调试时能立刻看到结果第三不需要配置复杂的环境装完就能跑这对教学和快速验证场景特别友好。实际上真正工业级的光学追迹软件比如专业的光学设计软件用的多是C或专用语言但那些工具的学习曲线陡、许可证贵。做一个教学演示、课程设计或初步验证MATLAB是性价比最高的选择。2. 核心算法与实现细节2.1 光线的数据结构和界面建模写仿真代码的第一步是定义清楚光线和界面这两个核心对象。我用结构体表示光线用自定义类或函数组表示界面。% 光线结构体 % origin: [x, y] 起点坐标 % direction: [dx, dy] 单位方向向量 % energy: 光强归一化到0~1 % color: [r, g, b] 用于可视化 ray.origin [0, 0]; ray.direction [1, 0]; ray.energy 1.0; ray.color [1, 0, 0]; % 红色光界面在2D情况下是一条直线段需要用两个端点定义。但更重要的是法向量——计算反射和折射时法向量的方向必须统一。我约定法向量始终指向光线入射的一侧。这个约定看似简单实际上能避免大量方向符号混乱的问题。如果法向量方向不统一反射方向公式 r i - 2(i·n)n 就会算出反向的反射光导致追踪结果完全错误。2.2 光线与界面的交点计算光线可以用参数方程表示P(t) origin t · direction平面界面可以用点法式方程表示n · (P - P₀) 0联立两个方程解得参数tt n · (P₀ - origin) / (n · direction)注意只有当t 0时才说明交点在光线的传播方向上t 0意味着交点在起点后面。而且n·direction不能为0——如果光线平行于界面结果会是无穷大这个情况要在代码里提前处理。function [t, hit] ray_plane_intersect(ray, plane) denominator dot(plane.normal, ray.direction); % 光线平行于界面不相交 if abs(denominator) 1e-12 t inf; hit false; return; end t dot(plane.normal, plane.point - ray.origin) / denominator; hit (t 0); end对于线段界面比如棱镜的某个面还需要额外判断交点是否在端点范围内。这个可以用向量叉积或简单的坐标比较来做效率上不需要太优化但逻辑上必须完整。2.3 折射与反射方向向量推导这是整个程序的核心公式推导必须严谨。设入射方向为单位向量i界面法向量指向入射侧为单位向量n两种介质的折射率分别为n₁、n₂。反射方向向量r i - 2 · (i · n) · n这个公式只要保证i和n都是单位向量结果一定也是单位向量不需要额外归一化。实现时要注意当i·n为负时说明法向量方向不对——入射光线在逆向穿过界面。折射方向向量比较复杂。设k₁ n₁/n₂则cosθ₁ -i · ncosθ₂ sqrt(1 - k₁² · (1 - cosθ₁²))折射方向t k₁ · i (k₁ · cosθ₁ - cosθ₂) · n检查开方内部是否小于0就能判断全反射是否发生。如果cosθ₂为虚数说明折射光不存在。function [ray_reflected, ray_refracted] trace_interface(ray, plane, n1, n2) i ray.direction; n plane.normal; % 已指向入射侧 cos_theta1 -dot(i, n); % 反射方向 r i - 2 * cos_theta1 * n; ray_reflected ray; ray_reflected.direction r; % 判断是否全反射 k n1 / n2; radicand 1 - k^2 * (1 - cos_theta1^2); if radicand 0 % 全反射没有折射光 ray_refracted []; else cos_theta2 sqrt(radicand); t k * i (k * cos_theta1 - cos_theta2) * n; ray_refracted ray; ray_refracted.direction t; end end这段代码里还有个细节cos_theta1可能出现负数或超过1的小误差这是浮点运算的正常现象。我的做法是在计算radicand之前做一个clip操作把cos_theta1限制在[0, 1]区间避免因微小数值误差导致误判全反射。2.4 能量分配与追踪深度控制光线每经过一个界面能量就会按菲涅尔方程重新分配。反射光线继承R比例的能量折射光线继承T 1-R比例的能量。如果继续追踪反射光和折射光就形成了树状的分支结构。但实际仿真中不能无限追踪下去否则计算量会爆炸。需要设定追踪深度上限比如10~15层以及能量阈值比如能量低于0.01就停止追踪。这两个参数的设计直接影响仿真速度和结果质量——深度太小会丢失多次反射的信息能量阈值太大会导致暗弱光线被截断。3. 深入实操从零搭建一个2D光线追踪器3.1 场景参数定义以平行光射入三角棱镜为例这是验证斯奈尔定律和色散效应的经典场景。棱镜材料通常用BK7玻璃折射率约1.52对钠黄光。我为了方便演示用单色光并把折射率设为1.52。% 界面定义棱镜的三个面 % 棱镜顶点坐标 p1 [0, 0]; p2 [3, 0]; p3 [1.5, 1.2]; % 三个面用点法式表示法向量指向棱镜外部 faceAB.point p1; faceAB.normal normalize([0, -1]); % 底面 faceAC.point p1; faceAC.normal normalize([-0.6, 1]); % 左斜面 faceBC.point p2; faceBC.normal normalize([0.6, 1]); % 右斜面 faces [faceAB, faceAC, faceBC]; % 入射光线从左侧水平射入 n_air 1.0; n_glass 1.52; ray.origin [-1.5, 0.4]; ray.direction [1, 0]; ray.energy 1.0;这个场景设置完成后先手工预估一下光路光线从空气进入棱镜左斜面折射角大约arcsin(sin(入射角)/1.52)从棱镜内部射向右斜面时可能会发生全反射或折射。手算一遍再跟程序结果对比是最有效的验证方法。3.2 核心追踪函数从单次交互到递归追踪单次光线与界面的交互已经写好了接下来要解决的核心问题是怎么找到光线遇到的第一个最近的界面我用的方法是遍历所有界面分别求交点取t最小且t0的那个。function [min_t, hit_face, pint] find_nearest_interface(ray, faces) min_t inf; hit_face []; pint []; for k 1:length(faces) face faces(k); [t, hit] ray_plane_intersect(ray, face); if hit t min_t min_t t; hit_face face; pint ray.origin t * ray.direction; end end if isempty(hit_face) % 光线逃逸出场景 min_t inf; end end有了单次交互和最近界面查找递归追踪就顺理成章了function trace_recursive(ray, faces, n_inside, depth, energy_threshold) if depth max_depth || ray.energy energy_threshold draw_ray(ray, depth); % 画出当前路径并终止 return; end [t, hit_face, pint] find_nearest_interface(ray, faces); if t inf % 没有碰到界面画出光线直到边界 draw_ray(ray, depth); return; end % 画出入射段 plot_ray_segment(ray.origin, pint, ray.color, ray.energy); % 计算反射率和透射率菲涅尔方程 [R, T] fresnel(ray, hit_face, n_outside, n_inside); % 构造反射光和折射光 [ray_r, ray_t] trace_interface(ray, hit_face, n_outside, n_inside); % 能量更新 ray_r.energy ray.energy * R; ray_t.energy ray.energy * T; % 递归追踪 ray_r.origin pint; ray_t.origin pint; if R energy_threshold trace_recursive(ray_r, faces, n_inside, depth 1, energy_threshold); end if T energy_threshold trace_recursive(ray_t, faces, n_inside, depth 1, energy_threshold); end end递归函数的终止条件很关键一是深度上限二是能量下限。这两个条件保证了程序不会陷入无限循环——尤其是全反射场景下如果深度上限不够光线就会在棱镜内部无限反弹下去虽然物理上没错但计算和绘图都会崩溃。3.3 从光密到光疏的介质切换上一步的递归代码里折射率参数n_outside和n_inside需要根据光线当前所在介质动态切换。这个状态管理在多层介质场景中尤其重要。我的做法是不在函数参数里显式传递n_outside和n_inside而是让每条光线携带一个介质编号字段比如material_id每次进入新介质时就更新这个字段。这样追踪逻辑更清晰后续扩展多层系统也更方便。function [n_outside, n_inside] get_materials(ray, hit_face) % 根据光线当前所在位置和界面的介质设定返回两侧折射率 if ray.material_id 1 % 空气 n_outside n_air; n_inside n_glass; else % 玻璃内部 n_outside n_glass; n_inside n_air; end end3.4 可视化展示让光路看得见绘图部分我用了两个技巧一是用线宽表示能量强度——光线能量越强线越粗让读者一眼看出能量在界面处的分配二是用颜色深浅表示深度——早期光线深度小颜色饱满后期多次反射的光线颜色变淡反映能量衰减。function plot_ray_segment(origin, dest, color, energy, depth) line_width 0.5 energy * 2.5; line_color color * (1 - depth * 0.05); % 每加深一层颜色稍微变暗 plot([origin(1), dest(1)], [origin(2), dest(2)], ... Color, line_color, LineWidth, line_width); hold on; end这里需要注意hold on和hold off的配合避免多个光线绘图时互相覆盖。我习惯在绘图前调用figure然后单独用hold on管理状态最后用axis equal保证横纵比例一致——否则圆形光束会变成椭圆形看着很别扭。3.5 结果验证理论值和仿真值对照写完代码不能直接信输出拿几个典型角度做手算验证是必须的。我用入射角30度、空气→玻璃n1.52验证过一次斯奈尔定律手算θ₂ arcsin(sin(30°)/1.52) ≈ 19.2°程序输出θ₂ ≈ 19.2°反射率R ≈ 4.3%非偏振光垂直入射时菲涅尔方程给出R |(1.52-1)/(1.521)|² ≈ 4.27%和斜入射30度时的结果几乎一致——这正是理论上预期的入射角较小时反射率变化缓慢接近临界角时才急剧上升。我还拿45度入射角验证过如果光线从玻璃内部射向空气n1.52→1.0此时临界角θc arcsin(1/1.52) ≈ 41.1°。入射角45度已经超过临界角程序输出全反射所有能量都返回玻璃内部——这个结果和物理预期完全吻合。4. 常见问题与排查技巧实录4.1 明明发生了折射光线却消失了这是我最常遇到的坑。查了很久才发现问题出在法向量方向不统一上。如果一个界面的法向量指向右侧另一个指向左侧那么dot(i, n)的正负号就不一致cos_theta1有时为正有时为负折射方向公式算出来的结果完全错误。解决办法在定义界面时就明确规定法向量指向入射侧。如果光线从另一侧入射先翻转法向量再计算。这个翻转可以写在find_nearest_interface里每找到交点时判断入射方向与法向量的点积如果为正说明法向量和入射方向同向就把法向量取反。4.2 能量不守恒反射率加透射率不等于1菲涅尔方程在理论上保证R T 1但程序输出的结果有时不满足。排查后发现两个原因一是浮点精度——如果没有做cos_theta1的clip操作开方前的小负数误差会导致能量出现微小负值二是透射光方向没有归一化导致后续交点计算偏差间接影响能量分配。我的修复方案所有方向向量在进入下一次递归前强制归一化cos_theta1先clip到[0,1]再参与所有计算。加了这两行之后RT的偏差从千分之一级别降到了1e-12级别几乎完美。4.3 反射方向算反了光反射到界面背面这种问题几乎都是向量公式使用不当导致的。我debug时喜欢在关键点打印中间量dot(i, n)的值、反射方向的分量、法向量的分量。如果dot(i, n)为正说明法向量指向了错误的方向和入射光同侧反射方向就会算反。正常情况下入射光方向i和指向入射侧的法向量n的点积必然为负——因为i指向前方n指向后方。看到这个值为正第一反应就应该是检查法向量定义而不是质疑公式。4.4 追踪性能优化很多光线一起追时卡顿如果一次要发射几百条光线比如模拟平行光束直接递归会有明显延迟。我试过两种优化方案第一把场景中的静态数据界面点坐标、法向量预先向量化不要每条光线都重新访问结构体数组。第二合并递归调用——如果只关心总的能量分布不需要每条光线独立绘图可以把一系列光线打包成矩阵用矩阵运算批量求交点。实测下来100条光线逐条追踪耗时约0.8秒批量矩阵运算耗时约0.15秒提速5倍以上。如果做实时交互演示建议用批量方案。5. 经验之谈让仿真结果更接近真实光路5.1 界面粗糙度与散射的简化处理真实世界没有绝对光滑的界面光在粗糙界面上的行为是散射而非简单的反射和折射。这个模拟项目没有包含散射但如果你想更贴近实际可以在每次反射/折射后给方向向量加一个小的随机扰动扰动幅度和界面粗糙度挂钩。这个改动在代码里只需两三行但对视觉效果的提升非常明显。5.2 多波长色散的模拟思路如果想把一束白光拆成彩虹色需要模拟色散折射率随波长变化。简单起见可以只追踪红、绿、蓝三束光分别用对应波长的折射率红光的折射率略小于蓝光。三束光的折射角会有细微差别在棱镜出口就能看到明显的色散效果。我在案例中把这条逻辑加进去了红光折射率1.514绿光1.517蓝光1.521对BK7玻璃而言实际差异大约如此。三束光分别追踪后叠加绘图出来的色散图非常直观——这是很多教学演示里最喜欢用的效果。5.3 代码模块化的建议整个项目我拆成了四个文件ray_intersect.m求交、fresnel_refraction.m反射折射方向与能量、trace_scene.m递归控制、draw_scene.m绘图。这样每个文件功能单一调试时定位问题很快。如果你打算在此基础上扩展建议沿用这个拆分思路——比如加一个光源类或者加一个波长管理模块都不需要动核心算法。最后分享一个小技巧调试光线追踪时先在代码里放一个调试模式开关输出每条光线的每一步位置和能量值。等验证物理正确后再关掉这个输出只保留最终绘图。这个习惯帮我省掉了至少一半的debug时间——毕竟用眼睛盯着上百条光线的叠加图找错误远不如直接看数值来得快。 
长河君 · 2026/10/13 15:58:49 · 阅读 1.2k MATLAB光学仿真：斯奈尔定律与菲涅尔方程的光线追踪实现
做光学仿真这几年最常被问到的一个问题是光从空气进水里到底偏转多少能量又分走多少教科书上两个公式就能回答但真要在电脑里把光的路径一步步画出来还得靠代码落地。这个基于MATLAB的斯奈尔定律与菲涅尔方程光线追踪模拟就是解决这个问题的——把几何光学的两大核心定律翻译成可执行的程序让界面折射、全反射、能量分配这些抽象概念变成屏幕上直观的光路图。它适合三类人一是正在学物理光学、想把公式跑起来的学生二是做光学系统设计、需要在方案阶段快速验证光路是否合理的工程师三是对光线追踪感兴趣的MATLAB爱好者。有了这套代码你能模拟一束光经过平面界面、棱镜甚至多层介质后的传播路径还能看到每个界面上反射光和透射光各自分到了多少能量——这对理解透镜系统、镀膜设计、光纤通信的原理都很有价值。1. 整体设计与物理原理拆解1.1 斯奈尔定律光的拐弯规则光从一种透明介质进入另一种透明介质时传播方向会发生变化这就是折射。斯奈尔定律定量描述了这种变化n₁·sinθ₁ n₂·sinθ₂其中n₁、n₂分别是入射介质和折射介质的折射率θ₁是入射角θ₂是折射角。角度都是从界面法线量起的不是从界面本身。这个公式看起来简单但实现时有一个容易被忽略的关键点入射角必须换算为弧度制而且正弦函数在角度接近90度时变化非常缓慢。如果在代码里直接用角度制计算或者把法线方向弄反了算出来的折射方向就会完全错误。反射定律同样重要反射角等于入射角。在代码里实现反射方向计算时不需要查表直接用向量运算就能搞定——入射方向向量减去两倍的法线方向投影就是反射方向。这个向量公式比角度计算更稳定也更适合编程实现。1.2 菲涅尔方程界面的能量分配规则斯奈尔定律告诉光往哪儿走菲涅尔方程则告诉光分多少能量走。光打到界面上时一部分能量被反射一部分能量透射过去这个比例由菲涅尔方程给出。对于非偏振光自然光反射率R取s偏振和p偏振的平均值Rs |(n₁·cosθ₁ - n₂·cosθ₂) / (n₁·cosθ₁ n₂·cosθ₂)|²Rp |(n₁·cosθ₂ - n₂·cosθ₁) / (n₁·cosθ₂ n₂·cosθ₁)|²R (Rs Rp) / 2透射率T 1 - R在无吸收的理想界面上能量守恒。实现菲涅尔方程时有个细节值得注意当入射角为0度垂直入射时Rs和Rp的计算公式会退化此时R |(n₁ - n₂)/(n₁ n₂)|²。而入射角接近临界角时反射率会急剧上升——虽然还没有全反射但能量已经开始大量回流这个效应在仿真中非常明显也是很多初学者容易遗漏的物理细节。1.3 全反射光线追踪的边界case当光从光密介质折射率大射向光疏介质折射率小时如果入射角超过某个阈值折射光会消失所有能量都被反射回来。这个阈值叫临界角θc arcsin(n₂/n₁)在代码实现中临界角判断必须在计算折射方向之前完成。我的习惯是先用n₁·sinθ₁ n₂来判断是否全反射——这个判断比先算arcsin再比较角度更高效而且避免了反三角函数的数值误差。全反射在现实世界中应用极广光纤通信就是利用全反射把光锁在纤芯里传播钻石的闪耀也部分归功于全反射——高折射率意味着小临界角光更容易在内部反弹。1.4 为什么用MATLAB做这件事选MATLAB做光学追踪不是因为它在光学计算上有什么不可替代的优势而是因为它刚刚好。第一MATLAB的向量化运算让光线逐条并行处理变得很自然写循环也不用担心性能太差第二它的绘图系统开箱即用画光线、画界面、画法线都很方便调试时能立刻看到结果第三不需要配置复杂的环境装完就能跑这对教学和快速验证场景特别友好。实际上真正工业级的光学追迹软件比如专业的光学设计软件用的多是C或专用语言但那些工具的学习曲线陡、许可证贵。做一个教学演示、课程设计或初步验证MATLAB是性价比最高的选择。2. 核心算法与实现细节2.1 光线的数据结构和界面建模写仿真代码的第一步是定义清楚光线和界面这两个核心对象。我用结构体表示光线用自定义类或函数组表示界面。% 光线结构体 % origin: [x, y] 起点坐标 % direction: [dx, dy] 单位方向向量 % energy: 光强归一化到0~1 % color: [r, g, b] 用于可视化 ray.origin [0, 0]; ray.direction [1, 0]; ray.energy 1.0; ray.color [1, 0, 0]; % 红色光界面在2D情况下是一条直线段需要用两个端点定义。但更重要的是法向量——计算反射和折射时法向量的方向必须统一。我约定法向量始终指向光线入射的一侧。这个约定看似简单实际上能避免大量方向符号混乱的问题。如果法向量方向不统一反射方向公式 r i - 2(i·n)n 就会算出反向的反射光导致追踪结果完全错误。2.2 光线与界面的交点计算光线可以用参数方程表示P(t) origin t · direction平面界面可以用点法式方程表示n · (P - P₀) 0联立两个方程解得参数tt n · (P₀ - origin) / (n · direction)注意只有当t 0时才说明交点在光线的传播方向上t 0意味着交点在起点后面。而且n·direction不能为0——如果光线平行于界面结果会是无穷大这个情况要在代码里提前处理。function [t, hit] ray_plane_intersect(ray, plane) denominator dot(plane.normal, ray.direction); % 光线平行于界面不相交 if abs(denominator) 1e-12 t inf; hit false; return; end t dot(plane.normal, plane.point - ray.origin) / denominator; hit (t 0); end对于线段界面比如棱镜的某个面还需要额外判断交点是否在端点范围内。这个可以用向量叉积或简单的坐标比较来做效率上不需要太优化但逻辑上必须完整。2.3 折射与反射方向向量推导这是整个程序的核心公式推导必须严谨。设入射方向为单位向量i界面法向量指向入射侧为单位向量n两种介质的折射率分别为n₁、n₂。反射方向向量r i - 2 · (i · n) · n这个公式只要保证i和n都是单位向量结果一定也是单位向量不需要额外归一化。实现时要注意当i·n为负时说明法向量方向不对——入射光线在逆向穿过界面。折射方向向量比较复杂。设k₁ n₁/n₂则cosθ₁ -i · ncosθ₂ sqrt(1 - k₁² · (1 - cosθ₁²))折射方向t k₁ · i (k₁ · cosθ₁ - cosθ₂) · n检查开方内部是否小于0就能判断全反射是否发生。如果cosθ₂为虚数说明折射光不存在。function [ray_reflected, ray_refracted] trace_interface(ray, plane, n1, n2) i ray.direction; n plane.normal; % 已指向入射侧 cos_theta1 -dot(i, n); % 反射方向 r i - 2 * cos_theta1 * n; ray_reflected ray; ray_reflected.direction r; % 判断是否全反射 k n1 / n2; radicand 1 - k^2 * (1 - cos_theta1^2); if radicand 0 % 全反射没有折射光 ray_refracted []; else cos_theta2 sqrt(radicand); t k * i (k * cos_theta1 - cos_theta2) * n; ray_refracted ray; ray_refracted.direction t; end end这段代码里还有个细节cos_theta1可能出现负数或超过1的小误差这是浮点运算的正常现象。我的做法是在计算radicand之前做一个clip操作把cos_theta1限制在[0, 1]区间避免因微小数值误差导致误判全反射。2.4 能量分配与追踪深度控制光线每经过一个界面能量就会按菲涅尔方程重新分配。反射光线继承R比例的能量折射光线继承T 1-R比例的能量。如果继续追踪反射光和折射光就形成了树状的分支结构。但实际仿真中不能无限追踪下去否则计算量会爆炸。需要设定追踪深度上限比如10~15层以及能量阈值比如能量低于0.01就停止追踪。这两个参数的设计直接影响仿真速度和结果质量——深度太小会丢失多次反射的信息能量阈值太大会导致暗弱光线被截断。3. 深入实操从零搭建一个2D光线追踪器3.1 场景参数定义以平行光射入三角棱镜为例这是验证斯奈尔定律和色散效应的经典场景。棱镜材料通常用BK7玻璃折射率约1.52对钠黄光。我为了方便演示用单色光并把折射率设为1.52。% 界面定义棱镜的三个面 % 棱镜顶点坐标 p1 [0, 0]; p2 [3, 0]; p3 [1.5, 1.2]; % 三个面用点法式表示法向量指向棱镜外部 faceAB.point p1; faceAB.normal normalize([0, -1]); % 底面 faceAC.point p1; faceAC.normal normalize([-0.6, 1]); % 左斜面 faceBC.point p2; faceBC.normal normalize([0.6, 1]); % 右斜面 faces [faceAB, faceAC, faceBC]; % 入射光线从左侧水平射入 n_air 1.0; n_glass 1.52; ray.origin [-1.5, 0.4]; ray.direction [1, 0]; ray.energy 1.0;这个场景设置完成后先手工预估一下光路光线从空气进入棱镜左斜面折射角大约arcsin(sin(入射角)/1.52)从棱镜内部射向右斜面时可能会发生全反射或折射。手算一遍再跟程序结果对比是最有效的验证方法。3.2 核心追踪函数从单次交互到递归追踪单次光线与界面的交互已经写好了接下来要解决的核心问题是怎么找到光线遇到的第一个最近的界面我用的方法是遍历所有界面分别求交点取t最小且t0的那个。function [min_t, hit_face, pint] find_nearest_interface(ray, faces) min_t inf; hit_face []; pint []; for k 1:length(faces) face faces(k); [t, hit] ray_plane_intersect(ray, face); if hit t min_t min_t t; hit_face face; pint ray.origin t * ray.direction; end end if isempty(hit_face) % 光线逃逸出场景 min_t inf; end end有了单次交互和最近界面查找递归追踪就顺理成章了function trace_recursive(ray, faces, n_inside, depth, energy_threshold) if depth max_depth || ray.energy energy_threshold draw_ray(ray, depth); % 画出当前路径并终止 return; end [t, hit_face, pint] find_nearest_interface(ray, faces); if t inf % 没有碰到界面画出光线直到边界 draw_ray(ray, depth); return; end % 画出入射段 plot_ray_segment(ray.origin, pint, ray.color, ray.energy); % 计算反射率和透射率菲涅尔方程 [R, T] fresnel(ray, hit_face, n_outside, n_inside); % 构造反射光和折射光 [ray_r, ray_t] trace_interface(ray, hit_face, n_outside, n_inside); % 能量更新 ray_r.energy ray.energy * R; ray_t.energy ray.energy * T; % 递归追踪 ray_r.origin pint; ray_t.origin pint; if R energy_threshold trace_recursive(ray_r, faces, n_inside, depth 1, energy_threshold); end if T energy_threshold trace_recursive(ray_t, faces, n_inside, depth 1, energy_threshold); end end递归函数的终止条件很关键一是深度上限二是能量下限。这两个条件保证了程序不会陷入无限循环——尤其是全反射场景下如果深度上限不够光线就会在棱镜内部无限反弹下去虽然物理上没错但计算和绘图都会崩溃。3.3 从光密到光疏的介质切换上一步的递归代码里折射率参数n_outside和n_inside需要根据光线当前所在介质动态切换。这个状态管理在多层介质场景中尤其重要。我的做法是不在函数参数里显式传递n_outside和n_inside而是让每条光线携带一个介质编号字段比如material_id每次进入新介质时就更新这个字段。这样追踪逻辑更清晰后续扩展多层系统也更方便。function [n_outside, n_inside] get_materials(ray, hit_face) % 根据光线当前所在位置和界面的介质设定返回两侧折射率 if ray.material_id 1 % 空气 n_outside n_air; n_inside n_glass; else % 玻璃内部 n_outside n_glass; n_inside n_air; end end3.4 可视化展示让光路看得见绘图部分我用了两个技巧一是用线宽表示能量强度——光线能量越强线越粗让读者一眼看出能量在界面处的分配二是用颜色深浅表示深度——早期光线深度小颜色饱满后期多次反射的光线颜色变淡反映能量衰减。function plot_ray_segment(origin, dest, color, energy, depth) line_width 0.5 energy * 2.5; line_color color * (1 - depth * 0.05); % 每加深一层颜色稍微变暗 plot([origin(1), dest(1)], [origin(2), dest(2)], ... Color, line_color, LineWidth, line_width); hold on; end这里需要注意hold on和hold off的配合避免多个光线绘图时互相覆盖。我习惯在绘图前调用figure然后单独用hold on管理状态最后用axis equal保证横纵比例一致——否则圆形光束会变成椭圆形看着很别扭。3.5 结果验证理论值和仿真值对照写完代码不能直接信输出拿几个典型角度做手算验证是必须的。我用入射角30度、空气→玻璃n1.52验证过一次斯奈尔定律手算θ₂ arcsin(sin(30°)/1.52) ≈ 19.2°程序输出θ₂ ≈ 19.2°反射率R ≈ 4.3%非偏振光垂直入射时菲涅尔方程给出R |(1.52-1)/(1.521)|² ≈ 4.27%和斜入射30度时的结果几乎一致——这正是理论上预期的入射角较小时反射率变化缓慢接近临界角时才急剧上升。我还拿45度入射角验证过如果光线从玻璃内部射向空气n1.52→1.0此时临界角θc arcsin(1/1.52) ≈ 41.1°。入射角45度已经超过临界角程序输出全反射所有能量都返回玻璃内部——这个结果和物理预期完全吻合。4. 常见问题与排查技巧实录4.1 明明发生了折射光线却消失了这是我最常遇到的坑。查了很久才发现问题出在法向量方向不统一上。如果一个界面的法向量指向右侧另一个指向左侧那么dot(i, n)的正负号就不一致cos_theta1有时为正有时为负折射方向公式算出来的结果完全错误。解决办法在定义界面时就明确规定法向量指向入射侧。如果光线从另一侧入射先翻转法向量再计算。这个翻转可以写在find_nearest_interface里每找到交点时判断入射方向与法向量的点积如果为正说明法向量和入射方向同向就把法向量取反。4.2 能量不守恒反射率加透射率不等于1菲涅尔方程在理论上保证R T 1但程序输出的结果有时不满足。排查后发现两个原因一是浮点精度——如果没有做cos_theta1的clip操作开方前的小负数误差会导致能量出现微小负值二是透射光方向没有归一化导致后续交点计算偏差间接影响能量分配。我的修复方案所有方向向量在进入下一次递归前强制归一化cos_theta1先clip到[0,1]再参与所有计算。加了这两行之后RT的偏差从千分之一级别降到了1e-12级别几乎完美。4.3 反射方向算反了光反射到界面背面这种问题几乎都是向量公式使用不当导致的。我debug时喜欢在关键点打印中间量dot(i, n)的值、反射方向的分量、法向量的分量。如果dot(i, n)为正说明法向量指向了错误的方向和入射光同侧反射方向就会算反。正常情况下入射光方向i和指向入射侧的法向量n的点积必然为负——因为i指向前方n指向后方。看到这个值为正第一反应就应该是检查法向量定义而不是质疑公式。4.4 追踪性能优化很多光线一起追时卡顿如果一次要发射几百条光线比如模拟平行光束直接递归会有明显延迟。我试过两种优化方案第一把场景中的静态数据界面点坐标、法向量预先向量化不要每条光线都重新访问结构体数组。第二合并递归调用——如果只关心总的能量分布不需要每条光线独立绘图可以把一系列光线打包成矩阵用矩阵运算批量求交点。实测下来100条光线逐条追踪耗时约0.8秒批量矩阵运算耗时约0.15秒提速5倍以上。如果做实时交互演示建议用批量方案。5. 经验之谈让仿真结果更接近真实光路5.1 界面粗糙度与散射的简化处理真实世界没有绝对光滑的界面光在粗糙界面上的行为是散射而非简单的反射和折射。这个模拟项目没有包含散射但如果你想更贴近实际可以在每次反射/折射后给方向向量加一个小的随机扰动扰动幅度和界面粗糙度挂钩。这个改动在代码里只需两三行但对视觉效果的提升非常明显。5.2 多波长色散的模拟思路如果想把一束白光拆成彩虹色需要模拟色散折射率随波长变化。简单起见可以只追踪红、绿、蓝三束光分别用对应波长的折射率红光的折射率略小于蓝光。三束光的折射角会有细微差别在棱镜出口就能看到明显的色散效果。我在案例中把这条逻辑加进去了红光折射率1.514绿光1.517蓝光1.521对BK7玻璃而言实际差异大约如此。三束光分别追踪后叠加绘图出来的色散图非常直观——这是很多教学演示里最喜欢用的效果。5.3 代码模块化的建议整个项目我拆成了四个文件ray_intersect.m求交、fresnel_refraction.m反射折射方向与能量、trace_scene.m递归控制、draw_scene.m绘图。这样每个文件功能单一调试时定位问题很快。如果你打算在此基础上扩展建议沿用这个拆分思路——比如加一个光源类或者加一个波长管理模块都不需要动核心算法。最后分享一个小技巧调试光线追踪时先在代码里放一个调试模式开关输出每条光线的每一步位置和能量值。等验证物理正确后再关掉这个输出只保留最终绘图。这个习惯帮我省掉了至少一半的debug时间——毕竟用眼睛盯着上百条光线的叠加图找错误远不如直接看数值来得快。