Matlab光路仿真：PQ向量光线追迹源码与工程实践

简介：这是一份基于PQ分解法的MATLAB潮流计算源码，主要面向电力系统分析初学者、电气工程专业学生以及需要快速掌握潮流计算原理的工程师。源码通过清晰的代码结构展示如何建立电力网络拓扑模型，区分PQ节点（负荷节点）与PV节点（发电机节点），并依据基尔霍夫定律构造功率平衡方程组，再采用牛顿-拉弗森迭代法求解非线性方程，最终得到各节点电压幅值、相角及功率分布。压缩包中仅包含一个PQ.m文件，文件大小仅2KB，代码轻量、结构紧凑，便于逐行阅读和调试。该资源目前已获得160人次的下载学习，学习热度稳定。通过实际运行此程序，您可以深入理解电力系统稳态分析中节点分类、雅可比矩阵、迭代收敛等关键知识点，同时提高使用MATLAB进行矩阵运算和算法实现的编程能力，也为后续学习牛顿法、高斯-塞德尔法等其他潮流算法提供了可直接修改与扩展的基础。

1. 为什么光路仿真要把一条光线拆成 PQ 两个向量

多数人第一次在代码库拿到“PQ,matlab光路程序源码”这类资源时，以为PQ是某个算法的缩写；其实它就是光线状态的两种描述：P是光线与参考面的交点坐标，Q是光线的传播方向。用Matlab做光路仿真，核心工作就是更新每一条光线的(P,Q)，让它在空气中直线传播、在透镜表面按Snell定律折射、在反射镜上转向，最后落在接收面上。这个表示看似简单，却比用角度追迹健壮得多，也方便与商用光学软件对齐。

光学设计、激光系统、课程仿真的从业者都会遇到符号越调越乱的情况：角度在界面上要看象限、判断正负；换成PQ之后，问题变成向量运算，全反射也能自然暴露出来。下文从PQ的几何定义开始，给出一个可运行的Matlab追迹源码框架，用三片式透镜组跑出点列图，再把追迹结果接到优化工具箱和外部软件上，解决源码“能跑但不敢用”的问题。

2. 从PQ坐标到传输矩阵：光路源码背后的几何模型

2.1 PQ光线状态的定义与坐标约定

一条光线在某个参考坐标系里可以用两个三维向量表示：位置向量P=(x,y,z)和光学方向余弦向量Q=(L,M,N)。其中Q不是简单的单位方向矢量，而是介质折射率与方向余弦的乘积，即Q=n·(cosα, cosβ, cosγ)，α、β、γ是光线与三个坐标轴的夹角。这样定义有一个很重要的性质：在均匀介质里，|Q|恒等于所在介质的折射率n，光线经过折射面时，Q的切向分量连续，法向分量按折射率跳变。

在傍轴光学里，如果光轴沿z方向，通常只保留四维子空间(px,py,qx,qy)，并把向量长度归一化近似为1。很多源码包在开头的坐标转换函数里做的事情，就是把用户输入的物距、视场角换算成这四维量。这里建议保留一个约定：P和Q永远使用同一个右手直角坐标系，z朝系统出射方向，而不是让每个表面都建局部系。这样做能让后续的光线交点计算省掉大量坐标变换，代价是到了离轴反射系统时，需要额外维护一个从全局到局部的旋转矩阵。

2.2 常见光学元件的近轴矩阵与PQ的边界条件

当系统满足傍轴条件时，光线在一个元件前后的状态可以用2×2矩阵连接起来：(x', qx')=M·(x, qx)。下面是几类常见元件的ABCD矩阵，使用角度量与方向余弦混合的傍轴形式：

这张表对应的是近轴一阶量。源码里如果直接拿这些矩阵去计算大视场或者短焦距系统，点列图会明显偏离真实结果。成熟的PQ光路程序会把“近轴矩阵模式”和“精确向量追迹模式”分开：前者用于初始设计和快速评估，后者用于像差分析。实现时不需要两个引擎，只要在每个表面保留解析求交，就已经是精确模式；矩阵模式只是对入射光线做快速变换，不真正计算交点。

2.3 为什么源码里用方向余弦而不是角度来存Q

角度表示在界面处必须做sin、cos和象限判断，而在向量形式下，折射定律可以写成紧凑的迭代式。设界面法线N指向入射介质一侧，入射光学方向余弦向量为Q1，折射方向为：

Q2 = Q1 − (n2·cosI2 − n1·cosI1)·N

其中cosI1 = −(Q1·N)/n1，cosI2由折射定律的平方根形式给出。这个公式在代码里实现只要几步：计算cosI1，判断是否小于0表示光线从背面入射，再取平方根得到cosI2，最后合成Q2。当表达式根号内出现负数时就是全反射，返回值可以直接标记为ray.valid=false，而不用像角度法那样先把角转到“看得懂”的范围。常见做法是把无效光线留在数组里只做标记，后面评价函数统一过滤，避免在追迹循环里动态缩减数组导致索引错位。

3. 把PQ光路源码整理成可维护的Matlab工程

3.1 用struct组织光线与光学表面

拿到一个光路程序源码包之后，先别急着读主脚本，读数据结构更快。我一般会用以下定义来组织光线和光学表面：

% ray: 3xN 矩阵存储，1条光线占1列 rays.P = zeros(3, N); % 位置向量 (x,y,z)，单位mm rays.Q = zeros(3, N); % 光学方向余弦 (L,M,N)，无量纲 rays.valid = true(1, N);% 标记全反射或被口径截断的光线 % 光学表面：1个元素代表1个表面顶点 surf.z = 0; % 表面顶点z坐标，mm surf.R = 20; % 曲率半径，mm，R=inf代表平面 surf.n1 = 1.0; % 入射侧折射率 surf.n2 = 1.5168; % 出射侧折射率，典型BK7在587nm surf.radius = 10; % 半口径，mm surf.type = "refract"; % 或 "reflect"

用3×N而不是1×N的结构体，是为了让后续交点计算能直接利用Matlab的向量化点乘。N不必提前预知，追迹过程也尽量不要改变数组长度，把失效光线用valid标记下来是源码可维护性的关键。注意surf里存的是顶点的全局z坐标，而不是相邻面的间距；相邻间隔在装配时被换算成z，可以避免“厚度加错符号”这类低级错误。

3.2 传播、折射与反射三个基础函数

自由空间传播是把每条光线沿自身方向移动已知距离d。由于Q带有折射率长度，先归一化再乘d：

function rays = advance(rays, d) nrm = vecnorm(rays.Q, 2, 1); % 每个方向余弦矢量的长度，即折射率 dir = rays.Q ./ nrm; % 单位方向向量 rays.P = rays.P + dir .* d; % 3xN 矩阵整体平移 end

球面求交是光路程序里最容易写错的一步。设球心在C=(0,0,z_s+R)，半径R，光线从参考点P出发沿单位方向u前进，则二次方程|P+u·t−C|²=R²决定交点：

function [t, P1] = sphereIntersect(P0, u, C, R) d = P0 - C; b = 2 * dot(u, d); c = dot(d, d) - R^2; disc = b^2 - 4 * c; t = Inf; P1 = P0; if disc >= 0 s = sqrt(disc); t1 = (-b - s) / 2; % 正值解，靠近光源侧 t2 = (-b + s) / 2; t = min([t1 t2]); if t <= 0 t = max([t1 t2]); % 起点在球内时取远交点 end P1 = P0 + u * t; end end

这里的t是沿u方向的传播距离，不是z方向增量。计算时先取较小正根，如果两个根都小于等于0，说明光线起点在球内，这种情况少见，但处理带厚度透镜时会出现。折射更新采用2.3里的向量式Snell定律：

function Q2 = snell(Q1, N, n1, n2) % N是单位法线，指向入射介质一侧 cosI1 = -dot(Q1, N) / n1; sinI2_sq = (n1 / n2)^2 * (1 - cosI1^2); if sinI2_sq > 1 Q2 = []; % 全反射，由调用方标记valid return; end cosI2 = sqrt(1 - sinI2_sq); Q2 = Q1 - (n2 * cosI2 - n1 * cosI1) .* N; end

注意n1/n2传的是绝对值；如果光线从玻璃射向空气，调用方要自己把n1、n2对调，并且把N反向，方向余弦向量Q2的模长会自动变成n2。反射情况更简单：Q_ref = Q1 − 2·dot(Q1,N)·N，同样要求N为单位法线。三个函数都操作3×N矩阵，但snell里的cosI1是1×N向量，dot需要逐列计算；上面代码为了可读性用了Matlab的向量点乘，实际批量运算时把Q1和N分别转成3行N列，用Q1.*N再按行求和更高效。

3.3 用光学元件列表驱动整条光路

把3.1的surface结构体放进数组，就得到一条光路。表面之间的厚度不是预先传播的距离，而是通过顶点z坐标对曲面定位；球面求交天然给出传播距离。装配循环可以这样写：

function rays = traceThrough(sys, rays) for k = 1:length(sys) e = sys(k); u = rays.Q ./ vecnorm(rays.Q, 2, 1); % 单位方向 if isinf(e.R) d = (e.z - rays.P(3,:)) ./ u(3,:); rays.P = rays.P + u .* d; % 平面直接走到顶点平面 else C = [0; 0; e.z + e.R]; [~, rays.P] = sphereIntersect(rays.P, u, C, e.R); end N = surfaceNormal(e, rays.P); % 外法线 N = N .* sign(-dot(u, N)); % 翻转到指向入射介质 if e.type == "refract" Q2 = snell(rays.Q, N, e.n1, e.n2); else Q2 = reflectRay(rays.Q, N); end if isempty(Q2) rays.valid = false; rays.Q(:, ~rays.valid) = NaN; continue; end rays.Q = Q2; r = sqrt(rays.P(1,:).^2 + rays.P(2,:).^2); rays.valid = rays.valid & (r <= e.radius); end end

这个循环里最关键的一行是N = N .* sign(-dot(u, N))。球面外法线由球心指向表面，如果光线从左侧射入凸球面，外法线与光线同向，必须翻转才能让N指向入射介质；凹面则不一定。逐光线翻转比在surfaceNormal里写死方向可靠，也适配反射镜。追迹过程中不要删除失效列，置NaN和valid=false即可，否则后面评价函数的下标会乱。

3.4 追迹结果的自检手段

追迹完成后要先验证再画图。三条内存检查通常够用：第一，每条光线的|Q|应该等于光线当前所在介质的折射率，如果中途某次snell写错n1/n2的次序，这里立刻暴露；第二，P的更新必须来自二次方程的正根，若出现负根，表现是光线在表面处倒退，检查反射镜是否误写成了折射分支；第三，口径截断后的光线数要等于valid的true数量。这些检查写成assert放在traceThrough末尾，会比任何注释都可靠。尤其从网上下载的matlab光路源码，先加三段assert再跑，能省下很多查符号的时间。

4. 用PQ追迹跑通一套三片式透镜组：从参数到点列图

4.1 输入参数与系统装配

用一个三片式透镜组当测试对象：两侧是BK7，中间是SF2。表面参数如下表，波长587nm，单位mm：

注意S2在空气侧曲率半径是负号，表示球心在光轴左侧；S3的负号同理。装配时把上表逐行转成3.1的surface结构体，并由调用者检查顶点z的单调性。不要使用面间距作为输入，那样一旦有人在两个面之间插入坐标变换，整个系统就会错位。

4.2 平行光入射与追迹循环

入射光束在z=0处生成，覆盖半口径5mm：

[xg, yg] = meshgrid(linspace(-5, 5, 21)); N = numel(xg); rays.P = [xg(:)'; yg(:)'; zeros(1, N)]; rays.Q = [zeros(2, N); ones(1, N)]; % 沿z正向，|Q|=1 rays.valid = true(1, N); rays = traceThrough(sys, rays);

21×21=441条光线，用来评估几何光斑已经足够密。追迹逻辑直接复用traceThrough；如果用的是3.3的平面-球面混合版，像面是Inf平面，光线会被advance函数带到像面位置。跑完后检查sum(rays.valid)，出现大量false说明有光线在半口径边缘被截断，或者全反射判断太激进，要回到第4.4节找原因。

4.3 点列图与RMS光斑半径

成像质量的快速评估不依赖官方光学工具箱，只要最后一步光线在像面上的分布。绘制点列图和计算RMS半径：

figure; hold on; axis equal; plot(rays.P(1, rays.valid), rays.P(2, rays.valid), '.'); xlabel('x / mm'); ylabel('y / mm'); title('Spot Diagram'); grid on; r2 = rays.P(1,:).^2 + rays.P(2,:).^2; r2 = r2(rays.valid); rmsR = sqrt(mean(r2)); maxR = sqrt(max(r2)); fprintf('RMS radius=%.4f mm, max radius=%.4f mm\n', rmsR, maxR);

点列图整体呈圆形且RMS半径远小于艾里斑直径时，可以认为几何像差不是限制因素，后面的优化更多是调整像面位置。如果只关心子午面，把入射网格改成y=0的单行采样即可，但计算RMS时建议保留二维网格，避免漏掉彗差的方向信息。输出评价数字时，除了RMS半径，P-V（最远点与主光线的距离）也要一起看，两个指标可能给出相反的排序。想做成三维动画或放到3D大屏案例里展示，把P的三行直接plot3出来逐帧旋转视角即可；441条光线没问题，几万条就要先抽稀再画。

4.4 追迹源码最常见的四类坑

追迹结果出现NaN，最常见原因是透镜半口径小于光线需要经过的入射高度，导致某条光线的传播距离为负并越过球面顶点；此时不要急着删光线，把valid和NaN同步置位。第二类坑是曲率半径符号与坐标轴约定冲突：本文令球心在出射方向一侧为正；如果你拿到的源码包习惯相反，整个系统的符号要在装配脚本开头统一翻转，不要在追迹函数里打补丁。第三类是折射率引起的方向余弦长度误差，尤其从玻璃进入空气时，若忘记把|Q2|归一化到空气的1，后面的advance会放大位置误差。第四类是口径截断只压了口径参数，没有考虑光阑引起的边缘遮挡，导致点列图边缘出现一圈不该有的点；把光阑写成一个口径极小的独立表面即可。

5. 把追迹源码扩展成波前评估与自动优化

5.1 从光程累计重建波前

追迹除了给出交点和方向，还能顺带累积光程。在每段传播的返回值里增加一个光程字段，传播时累加d * n_medium，到像面时每条光线相对于主光线的OPD就组成了波前图。重建时把光线网格展开成二维矩阵：

W = reshape(OPD - OPD(mid), size(xg)); figure; surf(xg, yg, W*1e6, 'EdgeColor', 'none'); xlabel('x/mm'); ylabel('y/mm'); zlabel('OPD/um');

这个波前用于判断是球差还是像散：沿径向对称的圆环条纹对应球差，沿45度方向的花样对应像散。当点列图RMS被衍射极限掩盖时，波前图比点列图更早暴露问题。

5.2 用优化工具箱调整最后一块间隔

把像面的轴向位置当作变量，用matlab优化工具箱里的fminunc自动收敛到最小RMS。目标函数需要固定入射光束网格，否则变量变化时采样点位移会造成数值噪声：

function loss = evalImagePos(zImg) sys(end).z = zImg; % 像面顶点z rays = traceThrough(sys, raysIn); P = rays.P(:, rays.valid); loss = mean(sum(P.^2, 1)); % 平方光斑半径均值 end x0 = 80; zOpt = fminunc(@evalImagePos, x0, ... optimoptions('fminunc', 'Display', 'iter', 'FiniteDifferenceStepSize', 1e-4));

FiniteDifferenceStepSize建议显式给出；默认值在追迹这类内部带有坐标判断的目标函数上，可能跳过边界造成loss不光滑。常见做法是先粗扫一维z从70到100，记录loss曲线，再以最小值附近为初值进入fminunc，比直接跑优化快得多。

5.3 PQ数据的导出与外部软件联调

需要把追迹结果交给外部软件时，PQ格式是天然的交换格式。导出每一束光线的状态为文本行(x,y,z,L,M,N,valid)，用writematrix一次写入即可。如果目标是类似matlab调用tracepro那样的联合仿真，常见做法是让Matlab生成入射光线文件，外部软件改变接收面或光源属性后返回照度图，再用matlab图像处理里的regionprops统计照度图质心，形成闭环。这个流程里要注意坐标单位：Matlab里用mm，外部软件默认也要设置成mm；PQ的方向余弦依赖介质折射率，传给某些软件前需要乘上所在介质的折射率，否则材质折射率不为1时整体偏差。保存时用单精度可以把文本体积砍掉一半，代价是返回数据二次读取时有效位数下降，这个取舍在大批量离线分析时很划算。

本文还有配套的精品资源，点击获取