简介该压缩包提供一份基于MATLAB编写的厚透镜光路模拟程序面向光学工程、应用物理及图像处理方向的初学者与科研人员用于在不采用薄透镜近似的条件下追踪真实光线传播路径。整个资源仅含1个thicklens.m源代码文件压缩包大小811B结构精简便于快速阅读与二次开发。已有639人学习下载适合作为光学仿真入门练习或教学演示素材。程序定义了透镜几何参数与光线初始条件依据斯涅尔定律逐面计算折射并利用MATLAB绘图函数绘制光线轨迹与聚焦点同时涉及光学坐标建立、菲涅尔公式和矩阵光学等关键概念可帮助读者对比厚透镜与薄透镜模型在像差、焦点位置上的差异。通过修改代码中的参数或追加表面定义还能扩展为双胶合、非球面等复杂系统具有很好的可定制性与教学价值。1. 厚透镜光路模拟为什么薄透镜近似在这里不适用在几何光学里处理透镜成像最顺手的方式是薄透镜公式 (1/f (n-1)(1/R_1 - 1/R_2))配合高斯成像公式就能完成大部分设计和计算。但一旦透镜厚度相对焦距不可忽略这套近似就开始偏离真实光路主平面位置会移动焦点不再与透镜几何中心重合边缘光线和近轴光线的会聚点也不一致。thicklens.m 这套 MATLAB 程序做的就是去掉“薄透镜”和“近轴”两个假设用真实光线追迹逐点计算光线经过两个球面的折射路径得到厚透镜的焦点、主平面和实际光束会聚行为。适合的场景包括光学设计课程中对透镜参数的定量验证、需要精确理解厚透镜主点位置的工程计算以及作为自定义光学系统追迹代码的起点。下面直接从追迹原理、代码结构和结果验证三个层次拆开讲。2. 光线追迹原理斯涅尔定律与球面折射的向量形式2.1 薄透镜近似到底丢掉了什么薄透镜近似本质上是把透镜的两个球面合并到同一个平面上处理认为光线在通过该平面的瞬间发生偏折忽略了光线在玻璃内部的实际平移。对于焦距远大于厚度的透镜这种简化带来的误差在工程公差范围内但对于曲率半径较小、厚度较大的透镜比如平凸或双凸短焦透镜忽略厚度会导致主平面位置计算错误进而影响系统组装时的间隔设计。另一个被约掉的因素是球差的高阶项。近轴近似假设所有光线与光轴的夹角足够小sin 和 tan 直接替换为角度本身。实际光线以较大孔径角入射时边缘光线比近轴光线会聚得更早形成焦点弥散。thicklens.m 中的“无小角度近似”正是针对这一点每条光线都按实际入射角计算折射方向而不是做线性化处理。2.2 球面折射的向量计算光线追迹的核心是求光线与球面的交点然后在该交点处应用斯涅尔定律。设光线起点为 (P_0)方向为单位向量 (D)参数化表示为 (P(t) P_0 tD)。球心在 (C)、半径为 (R) 的球面满足 (|P - C|^2 R^2)代入后得到一个关于 (t) 的二次方程[ |D|^2 t^2 2D \cdot (P_0 - C) t |P_0 - C|^2 - R^2 0 ]求解得到交点后法向量为 (N (P - C)/R)。设入射介质折射率为 (n_1)出射介质折射率为 (n_2)根据斯涅尔定律的向量形式[ n_1 (D \times N) n_2 (D \times N) ]其中 (D) 是折射后的方向向量。实际计算中更常用的是先求折射向量在法线方向上的分量再求出切向分量。令 (cos\theta_1 -D \cdot N)则k n1 / n2 cosθ2 sqrt(1 - k^2 * (1 - cosθ1^2)) D k * D (k * cosθ1 - cosθ2) * N注意当根号内为负值时说明发生了全反射需要特殊标记。2.3 厚透镜的表面顺序与坐标系约定厚透镜包含两个球面第一个面光入射面和第二个面光出射面它们的球心分别在光轴的不同侧。thicklens.m 中我沿用的是从左到右的传播约定光轴为 x 轴垂直方向为 y 轴透镜中心厚度为 (t)两个球面的曲率半径分别为 (R_1) 和 (R_2)符号按照几何光学惯例——球心在出射方向一侧时曲率半径为正反之为负。追迹过程中每一条光线依次执行四个步骤从物方出发直线传播到第一个球面在第一个球面处折射在玻璃内部直线传播到第二个球面在第二个球面处再次折射后进入空气。这套顺序和光线追迹专业软件中的“surface by surface”逻辑一致只不过这里直接写在 MATLAB 脚本里方便逐行检查和修改。3. thicklens.m 代码拆解从光线生成到光路可视化3.1 全局参数区与输入约定打开 thicklens.m 后首先看到的是参数区。这里定义了透镜的几何尺寸、折射率和入射光束的参数。以下代码是典型的参数定义方式% 透镜几何参数 R1 50; % 第一面曲率半径单位 mm正值表示球心在右侧 R2 -50; % 第二面曲率半径负值表示球心在左侧 thickness 10; % 透镜中心厚度单位 mm n_glass 1.5168; % 玻璃折射率BK7 在 587.6nm 下的典型值 n_air 1.0; % 空气折射率 % 光线参数 num_rays 21; % 光线数量奇数条保证包含光轴上的中心光线 source_distance -100; % 光源到透镜第一面的距离单位 mm负号表示在左侧 beam_half_height 20; % 光束半高度决定边缘光线的入射高度参数说明R1和R2的符号需要特别注意程序中约定正值表示球心位于该面右侧负值表示球心在左侧。对于双凸透镜R1 0、R2 0是标准配置。num_rays设为奇数是为了让中间那根光线严格沿光轴传播便于后续验证焦点位置是否落在轴上。3.2 光线与球面求交的函数实现核心求交逻辑建议封装成独立函数方便复用和调试。以下是光线与球面求交的代码function [t_intersect, is_hit] ray_sphere_intersection(P0, D, C, R) % 光线与球面求交返回沿光线方向的参数 t % P0: 光线起点坐标1x2 行向量 [x, y] % D: 光线方向向量1x2 行向量需要归一化 % C: 球心坐标1x2 行向量 % R: 球面半径标量带符号正负决定球心方向 D D / norm(D); oc P0 - C; a dot(D, D); % 应为 1因为 D 已归一化 b 2 * dot(oc, D); c dot(oc, oc) - R^2; discriminant b^2 - 4 * a * c; if discriminant 0 t_intersect NaN; is_hit false; return; end t1 (-b - sqrt(discriminant)) / (2 * a); t2 (-b sqrt(discriminant)) / (2 * a); % 取正数且较小的那个交点表示光线前进方向上的第一个交点 candidates [t1, t2]; candidates candidates(candidates 0); if isempty(candidates) t_intersect NaN; is_hit false; else t_intersect min(candidates); is_hit true; end end逻辑说明光线参数方程 (P_0 tD) 代入球面方程后得到一个一元二次方程判别式小于零说明光线与球面无交点。在有两个正根时取较小的 t 表示进入球面的那个交点一正一负时取正根说明光线起点在球体内部往外传播遇到球面。这个函数是整个程序中最基础也最关键的部分后续所有折射计算都依赖它返回的交点坐标。3.3 折射方向的实现与边界条件处理求解折射方向时需要注意光线方向和法向量方向的夹角关系。MATLAB 实现如下function D_refracted refract_ray(D_incident, N_normal, n1, n2) % 斯涅尔定律向量形式 % D_incident: 入射光方向单位向量指向传播方向 % N_normal: 球面法向量单位向量指向球心即指向透镜内部 % n1: 入射介质折射率 % n2: 出射介质折射率 cos_theta1 -dot(D_incident, N_normal); % 防止浮点误差导致 cos 超出 [-1,1] cos_theta1 max(-1, min(1, cos_theta1)); sin_theta1_sq 1 - cos_theta1^2; ratio n1 / n2; sin_theta2_sq ratio^2 * sin_theta1_sq; % 检查全反射条件 if sin_theta2_sq 1 D_refracted [NaN, NaN]; % 全反射标记 return; end cos_theta2 sqrt(1 - sin_theta2_sq); % 折射向量合成 D_refracted ratio * D_incident (ratio * cos_theta1 - cos_theta2) * N_normal; D_refracted D_refracted / norm(D_refracted); end关键点cos_theta1的计算公式中负号来自光线方向与法向量的方向约定。入射光线指向球面法向量指向球心两者夹角如果小于 90 度cos_theta1为正值对应入射角小于 90 度。全反射判断不能省略特别是在处理大角度入射或从玻璃射向空气的光线时否则会出现复数折射方向导致绘图异常。3.4 主循环与光路画图主循环的基本流程对每条入射光线做四次追迹每到一个交点就往结果数组中记录一段线段。以下是核心追迹和绘图代码figure(Position, [100, 100, 1200, 500]); hold on; % 绘制透镜轮廓两个球面的弧线 theta1 linspace(-60, 60, 100) * pi/180; x1 R1 * sin(theta1); y1 thickness/2 - R1 * cos(theta1) R1; % 左表面 % 这里仅示意实际应计算球面在透镜厚度范围内的有效弧度 % 对每条光线进行追迹 for i 1:num_rays % 生成入射光线从物方出发方向平行于光轴 y0 linspace(-beam_half_height, beam_half_height, num_rays); P_start [source_distance, y0(i)]; D_inc [1, 0]; % 平行于光轴 % 追迹第一个面 C1 [0, 0]; % 第一面球心位置需要根据透镜位置修正 [t1, hit1] ray_sphere_intersection(P_start, D_inc, C1, R1); P1 P_start t1 * D_inc; % 计算第一面法向量并折射 N1 (P1 - C1) / R1; D_mid refract_ray(D_inc, N1, n_air, n_glass); % 追迹第二个面 C2 [thickness, 0]; % 第二面球心位置 [t2, hit2] ray_sphere_intersection(P1, D_mid, C2, R2); P2 P1 t2 * D_mid; % 第二次折射 N2 (P2 - C2) / R2; D_out refract_ray(D_mid, N2, n_glass, n_air); % 将出射光线延长到绘图区域边界 t_extend 200; P3 P2 t_extend * D_out; % 绘制光线路径 plot([P_start(1), P1(1)], [P_start(2), P1(2)], b-, LineWidth, 1.2); plot([P1(1), P2(1)], [P1(2), P2(2)], r-, LineWidth, 1.2); plot([P2(1), P3(1)], [P2(2), P3(2)], g-, LineWidth, 1.2); end xlabel(x (mm)); ylabel(y (mm)); axis equal; grid on;关于球心位置的计算第一面球心的 x 坐标需要根据球面顶点位置修正。设第一面顶点在 (x0)曲率半径 (R_150)则球心在 (x50) 处。但第二面的顶点在 (xthickness) 处其球心在 (x thickness R_2)。调试时最常见的错误就是把两个球心都放在光轴的同一位置上导致光线追迹结果看起来像是穿过了一个平板而不是透镜。3.5 透镜轮廓绘制的小技巧绘制透镜轮廓时不能直接画完整圆需要截取透镜厚度范围内的那段弧。对第一面来说y 坐标范围受限于透镜的口径而 x 坐标范围受限于厚度。一个简洁的做法是用参数方程按角度采样再筛选落在透镜范围内的点% 绘制第一面轮廓 angle1 linspace(-asin(beam_half_height / R1), asin(beam_half_height / R1), 100); x_surf1 R1 * sin(angle1); y_surf1 R1 - R1 * cos(angle1); % 球面顶点在 (0,0) valid1 x_surf1 0 x_surf1 thickness; plot(x_surf1(valid1), y_surf1(valid1), k-, LineWidth, 2);这段代码的关键在于y_surf1的表达式。球面顶点在原点球心在 ((R_1, 0))球面上一点的坐标满足 ((x - R_1)^2 y^2 R_1^2)取 y 为正值那一支就得到 (y \sqrt{R_1^2 - (x - R_1)^2})。按参数方程展开后即上面代码的形式。第二面同理以 ((thickness R_2, 0)) 为球心、以厚度位置为顶点做镜像处理。4. 光路模拟的数值稳定性与参数敏感性分析4.1 光线从玻璃射向空气时的全反射边界这是厚透镜模拟中最容易踩坑的地方。光线在第二个球面从玻璃进入空气折射率从 1.5168 降到 1.0对应的临界角约为 41.8 度。对于曲率半径较小的透镜边缘光线在第二面的入射角很容易超过临界角此时程序里必须明确处理入射角范围物理现象程序行为(\theta_1 \theta_{crit})正常折射光线出射计算折射方向并延长绘制(\theta_1 \theta_{crit})掠射出射角 90 度折射光线沿界面切向(\theta_1 \theta_{crit})全反射返回 NaN停止该光线追迹全反射在真实透镜中并不是错误而是大孔径光学系统设计时需要考虑的重要约束。比如在数值孔径很高的显微物镜中边缘光线在透镜内部的全反射会导致光能损失和杂散光。thicklens.m 中对这类光线直接跳过并做标记是合理的默认行为。如果希望可视化全反射路径可以把refract_ray中的返回值改为反射方向向量% 全反射时的反射方向 if sin_theta2_sq 1 D_reflect D_incident - 2 * dot(D_incident, N_normal) * N_normal; D_refracted D_reflect; return; end4.2 平行光入射的焦点位置计算用平行于光轴的光束入射时所有出射光线与光轴的交点就是焦点位置。但厚透镜存在球差边缘光线和近轴光线的焦点位置并不重合。可以从程序输出中提取每个交点的 x 坐标绘制成“入射高度-焦点位置”曲线% 提取每条光线的焦点位置 focus_x zeros(num_rays, 1); for i 1:num_rays % 光线追迹后得到出射方向 D_out 和出射点 P2 % 求与光轴 y0 的交点 if abs(D_out(2)) 1e-8 focus_x(i) NaN; % 平行于光轴无交点 else t_axis -P2(2) / D_out(2); focus_x(i) P2(1) t_axis * D_out(1); end end y_input linspace(-beam_half_height, beam_half_height, num_rays); plot(focus_x, y_input, o-); xlabel(焦点位置 x (mm)); ylabel(入射光线高度 (mm));这个曲线的形状直接反映了球差的大小。如果所有点落在同一条竖直线上说明系统无球差实际的厚透镜通常呈现边缘光线焦点偏左球差为负或偏右球差为正的趋势。这里得到的结果可以用球差系数做定量分析纵向球差定义为近轴焦点与边缘光线焦点之间的差值。4.3 曲率半径符号设置错误导致的怪异结果程序调试中频率最高的错误是三组参数之间的符号逻辑不一致。以双凸透镜为例第一面曲率半径为正球心在右第二面曲率半径为负球心在左。如果对两个面都填正曲率半径追迹出的光线会在第二个面处向错误方向折射表现为主光线不再沿光轴传播。排查时可以先用一根沿光轴的光线做单光线测试理想情况下它应该直接穿过透镜不改变方向。若发生偏折优先检查第二面法向量的指向。法向量方向的选择同样影响结果。ray_sphere_intersection返回的 t 是正根对应的法向量取 ( (P - C) / R )此时 R 的符号会影响法向量方向。对于第一面R1 为正入射光线从左侧进入交点的法向量指向球心右侧即光从空气进入玻璃时法向量与入射光近似同向保证 (cos\theta_1 \approx 1)。如果这个惯例被打破折射方向会出现 180 度的偏差。4.4 大孔径角入射时的追迹精度验证当入射光线与光轴夹角很大时参数的数值误差会被放大。一个简单的验证方法是使用能量守恒约束折射前后光线的方向向量模长必须保持为 1。程序每次调用refract_ray后都对结果做了归一化但如果输入方向向量未归一化这个约束会被破坏导致交点计算出现漂移。我建议在ray_sphere_intersection开头强制归一化而不是依赖外部传入的单位向量。另一个验证手段是逆向追迹。将从焦点发出的光线反向追迹回物方如果光路可逆最终出射光线应与初始入射光线平行。将两条光线的夹角作为误差指标可以在不引入额外光学理论的情况下验证追迹代码的正确性。5. 焦点、主平面与焦距的提取方法5.1 用焦距公式验证追迹结果厚透镜的焦距由透镜制造者公式精确描述[ \frac{1}{f} (n-1)\left[\frac{1}{R_1} - \frac{1}{R_2} \frac{(n-1)t}{n R_1 R_2}\right] ]将参数 (R_150) mm、(R_2-50) mm、(t10) mm、(n1.5168) 代入时可以同时用 MATLAB 计算理论值和追迹得到的数值解n 1.5168; R1 50; R2 -50; t 10; % 理论焦距 f_theory 1 / ((n-1) * (1/R1 - 1/R2 (n-1)*t/(n*R1*R2))); % 追迹焦距取近轴光线的焦点位置 % 程序输出 f_trace 后与理论值比较 fprintf(理论焦距: %.3f mm\n, f_theory); fprintf(追迹焦距: %.3f mm\n, f_trace); fprintf(偏差: %.3f mm\n, abs(f_theory - f_trace));注意透镜制造者公式算出的焦距是“焦距”即主平面到焦点的距离不是透镜后表面顶点到焦点的距离。两者之差正是主平面位移。程序中通过追迹平行光得到的焦点位置是“后焦距”back focal length与理论焦距之间存在一个系统偏差这个偏差等于主平面偏离透镜顶点的距离。5.2 主平面位置的数值提取主平面定义为延长入射光线和反向延长出射光线两者延长线的交点所在的垂直于光轴的平面。对平行入射光来说这条延长线交点的高度与入射高度一致因此可以更简单地提取对每一条光线计算入射光线延长线和出射光线反向延长线的交点全部交点的 x 坐标平均就是主平面位置。% 计算主平面位置入射和出射方向延长线交点 % 入射光线方程: y y0 (平行光入射) % 出射光线反向延长: P2 s * (-D_out) % 交点在 x x2 s * (-D_out_x)其中 s (y0 - y2) / (-D_out_y) H_plane_x zeros(num_rays, 1); for i 1:num_rays P2_i P2_array(i, :); D_out_i D_out_array(i, :); y0_i y0_array(i); s (y0_i - P2_i(2)) / (-D_out_i(2)); H_plane_x(i) P2_i(1) s * (-D_out_i(1)); end H_plane mean(H_plane_x);这里的思路是入射平行光线的延长线是水平直线出射光线的反向延长线与这根水平线的交点就在主平面上。把所有光线求得的 x 坐标取平均可以减少数值噪声。当球差较大时各条光线算出的主平面位置可能不重合这正是非近轴效应的表现。5.3 球差量化与光斑半径估计球差的直观表现是焦点在轴向上的弥散。将光线追迹得到的焦点位置列表排序最大值与最小值之差就是纵向球差longitudinal_aberration max(focus_x) - min(focus_x);如果关心的是像面上的弥散斑尺寸需要先把出射光线延长到一个给定的像面位置 (x_{image})计算光线在该位置的 y 坐标然后取最大差值作为几何光斑直径y_at_image P2(:, 2) (x_image - P2(:, 1)) .* (D_out(:, 2) ./ D_out(:, 1)); spot_diameter max(y_at_image) - min(y_at_image);当 (x_{image}) 取在近轴焦点处得到的是“最小弥散圆”直径取在边缘光线焦点处则反映边缘光线的会聚点位置。对于工程设计来说最小弥散圆位置通常更接近最佳像面。6. 从平行光到任意物点物距、像距与光路可逆性验证6.1 修改光源参数实现点光源追迹thicklens.m 默认处理的是平行光入射但真实光学系统中更多遇到的是有限物距的点光源。修改方式很简单将入射光线的起点改为物点坐标方向从物点指向透镜上不同的采样点。以下是点光源入射的实现改法% 物点位置 object_point [-80, 5]; % 物点x-80mmy5mm离轴物点 % 对每条光线从物点出发到透镜第一面上的采样点 for i 1:num_rays % 在第一面上采样 y_target linspace(-beam_half_height, beam_half_height, num_rays); x_target 0; % 第一面顶点位置 % 入射方向 D_inc [x_target - object_point(1), y_target(i) - object_point(2)]; D_inc D_inc / norm(D_inc); P_start object_point; % 后续追迹流程与平行光一致 % ... end物点离轴时需要特别留意边缘光线在第二面的入射角。离轴物点的光线到达透镜时的入射角比轴上物点更大更容易触发全反射条件。如果追迹过程中某一根光线返回 NaN不要直接跳过先检查是不是物点本身已经超出透镜的视场角。6.2 像点聚集性与像差分析点光源经过透镜后在像方形成的不是理想像点而是一个弥散斑。收集所有出射光线与指定像面的交点坐标绘制散点图可以直观看到几何像差分布figure; scatter(y_at_image, zeros(size(y_at_image)), 20, filled); xlabel(像面上的 y 坐标 (mm)); title(点光源成像弥散斑分布);如果追迹光线数量足够多散点图的分布范围反映了彗差和球差的综合效果。对于离轴物点弥散斑的重心会偏离理想像点位置这个偏移量就是畸变的几何表现。将物点高度逐步增加并重复追迹可以画出畸变曲线。6.3 光路可逆性检查作为正确性判据光线追迹代码最容易出错的环节是法向量方向或多介质折射率配对。一个高效的调试方法是做全程可逆性验证从像方焦点反向追迹光线看是否回到物点位置。具体做法是把物方和像方参数互换然后比较两个方向的追迹结果正向追迹物点到像方反向追迹把像方点作为“物点”用相反的传播方向追迹两条路径在透镜内部的交点坐标差应小于 (10^{-6}) mm这个验证不需要任何额外的光学理论纯粹是几何一致性检查。如果差值过大说明折射向量计算中存在方向约定问题优先检查refract_ray中的cos_theta1符号。6.4 点列图半径与衍射极限的对照对于成像系统几何追迹得到的弥散斑半径需要与衍射极限做对比才能判断系统是衍射受限还是像差受限。艾里斑半径公式为 (r_{airy} 1.22 \cdot \lambda \cdot f / D)其中 (D) 为光束口径。用 MATLAB 快速估算lambda 0.5876e-3; % 波长单位 mm587.6nm 绿光 f f_trace; % 追迹得到的焦距 D_aperture 2 * beam_half_height; r_airy 1.22 * lambda * f / D_aperture; fprintf(衍射极限半径: %.4f mm\n, r_airy); fprintf(几何弥散斑半直径: %.4f mm\n, spot_diameter / 2);当几何弥散斑远大于艾里斑半径时系统处于像差主导区增加口径只会加剧弥散当两者接近时系统已接近衍射极限继续优化透镜参数的空间不大。这个对比为修改R1、R2和thickness提供了量化目标。本文还有配套的精品资源点击获取