简介：面向太阳能热发电、光伏系统及光学设计工程师与研究者的MATLAB光线追迹程序包，可模拟光线在不同介质中的传播路径，分析反射、折射与吸收过程，并可针对集热器、光伏面板等典型场景进行性能评估。压缩包共7个文件，以5个m源码脚本为主体，分别负责主流程控制、反射与折射计算、光线交点求解以及示例演示，另附README说明与LICENSE许可文件，整体仅6KB，十分精简。已有566人学习浏览。程序采用Chebyshev多项式提升计算效率，允许自定义材料折射率、反射率、吸收率等参数，并通过MATLAB图形功能直观展示光线路径与强度分布。借助这套代码，读者可快速上手几何光学仿真，完成太阳能集热器聚焦优化、光伏阴影影响评估以及透镜反射镜系统布局设计，加深对菲涅尔公式与光路模拟的理解。
1. 光线追迹程序实战：从太阳能热发电到光伏光学仿真
做太阳能热发电（CSP）项目的人，有一半时间都耗在"光该往哪打"这件事上。聚光镜的曲面稍微偏一点、接收器位置歪个几毫米，系统年发电量就能差出好几个百分点。我拆过一套光线追迹相关的程序包，它本质上是用蒙特卡洛方法模拟光子从太阳出发、经过反射镜、最后落在接收器或光伏电池表面的完整路径，输出光学效率和能流分布。这套东西比用商业光学软件轻得多，适合做前期的光场设计、镜场排布验证，也适合拿来给学生讲清楚"什么是光学效率"。
程序的核心价值不在界面，而在算法：用统计抽样的方式替代解析求导，把复杂的镜面反射、遮挡、余弦损失都变成概率事件。你给一组镜面参数和太阳位置，它就能告诉你接收器上哪块能量最集中、边缘损失有多大。做 PV（光伏）的人也能用同样一套逻辑模拟电池板上的辐照均匀度，从而判断是否需要做二次光学设计。
2. 程序架构与算法基础：光线追迹的完整链路
2.1 从太阳圆盘到光线束：源头建模
光线追迹的第一步不是算反射，而是先把太阳"拆"成无数条光线。真实太阳不是点光源，它在天空中有个视场角，通常取 4.65 mrad（约 0.266 度）。程序里会用一段代码生成一个位于太阳圆盘内的随机方向向量，再把这些方向映射到聚光镜的入射光线上：
import numpy as np def sample_sun_direction(sun_vector, sun_half_angle=0.00465): """ 从太阳圆盘内按均匀分布采样入射方向 sun_vector: 太阳指向反射镜的单位向量（shape: 3,） sun_half_angle: 太阳半张角，单位弧度，CSP场景常用 4.65 mrad """ # 先构造一个以 sun_vector 为 z 轴的局部坐标系 z = sun_vector / np.linalg.norm(sun_vector) ref = np.array([1, 0, 0]) if abs(z[0]) < 0.9 else np.array([0, 1, 0]) x = np.cross(z, ref); x /= np.linalg.norm(x) y = np.cross(z, x) # 在太阳圆盘内均匀采样：面积均匀分布需要 sqrt 取半径 r = np.sqrt(np.random.random()) * sun_half_angle theta = np.random.random() * 2.0 * np.pi offset = r * (np.cos(theta) * x + np.sin(theta) * y) # 最终入射方向 = 太阳主方向 + 圆盘偏移（近似） return sun_vector + offset这里有两个容易被忽略的细节。第一，半径采样用sqrt(random())而不是直接用random()，否则光线会密集在太阳圆盘中心，导致能流分布偏尖；第二，局部坐标系的构造方式做了防退化处理——当太阳向量恰好落在[1,0,0]附近时，参照向量换成[0,1,0]，否则叉积结果会变成零向量。这个函数输出的每条光线方向都带有一个微小偏角，后续镜面反射就是基于这个偏折后的方向计算的。
2.2 镜面反射与交点计算：核心几何引擎
光线从太阳出发后，下一步是与镜面求交。抛物面镜是 CSP 系统里最常用的聚光形式，它的表面满足方程 z = (x² + y²) / (4f)，其中 f 是焦距。程序用迭代法解这个交点的位置：
def intersect_paraboloid(ray_origin, ray_dir, focal_length): """ 求解光线与旋转抛物面 z = (x^2 + y^2) / (4f) 的交点 返回交点坐标（numpy 数组），无交则返回 None """ a = ray_dir[0]**2 + ray_dir[1]**2 b = 2.0 * (ray_origin[0]*ray_dir[0] + ray_origin[1]*ray_dir[1]) - 4.0*focal_length*ray_dir[2] c = ray_origin[0]**2 + ray_origin[1]**2 - 4.0*focal_length*ray_origin[2] disc = b*b - 4.0*a*c if disc < 0 or abs(a) < 1e-12: return None # 光线与镜面不相交 t1 = (-b + np.sqrt(disc)) / (2.0*a) t2 = (-b - np.sqrt(disc)) / (2.0*a) t = min(t for t in (t1, t2) if t > 1e-10) if t is None: return None return ray_origin + t * ray_dir这个函数本质上是把抛物面方程代入光线参数方程，得到一个关于参数 t 的二次方程。程序里取较小的正根，因为在抛物面内侧反射时，光线先碰到的是靠近焦点的那个交点。实际调试中我踩过坑：忘记abs(a) < 1e-12这个保护判断时，平行于轴线的光线会计算出a=0，直接除以零。另一个常见问题是判别式接近零但不为负的情况——边缘光线恰好切过镜面时，浮点误差可能让disc出现负值，实际解决方案是加一个max(disc, 0)的截断。
2.3 能流统计与光学效率计算
光线反射后打到接收器上，程序需要统计的是能流密度分布——把接收面划分成网格，每条光线落到哪个格子就给哪个格子累加一份能量。这是蒙特卡洛方法收敛的关键：
def accumulate_flux(rays_hit_point, receiver_grid_x, receiver_grid_y, flux_map): """ rays_hit_point: 所有成功命中接收器的光线交点（Nx3 数组） receiver_grid_x, receiver_grid_y: 接收器网格坐标 flux_map: 二维数组，存储每个网格的累计能流 """ for point in rays_hit_point: # 找到最近的网格索引 ix = np.argmin(np.abs(receiver_grid_x - point[0])) iy = np.argmin(np.abs(receiver_grid_y - point[1])) flux_map[iy, ix] += 1.0 # 每个光子权重相等 return flux_map这个统计过程看似简单，但光线数量不够时会非常"噪"。我一般建议至少跑 10 万条光线，能看到平滑的能流分布；少于 1 万条的话，峰值位置会跳动得让人怀疑程序写错了。程序里还会额外记录每条光线是被镜面反射、被遮挡、还是直接飞出了系统边界，用这些统计数字计算光学效率——分母是入射总光线数，分子是落在接收器上的光线数，这个比值比单纯看能流图更直观。
3. 代码仓库结构与数据流：一个追迹任务怎么跑完
3.1 目录组织与模块分工
这套程序包的文件排布非常清楚，典型的光学仿真项目结构：入口脚本负责调度，核心模块负责几何计算，工具模块负责输入输出。我拆开看了下目录，大致是下面这个样子：
ray_tracing_project/ ├── main.py # 主入口：解析参数、调用追迹、输出结果 ├── core/ │ ├── sun.py # 太阳方向采样模块 │ ├── mirror.py # 镜面定义、交点求解 │ ├── receiver.py # 接收器网格初始化与能流累积 │ └── tracer.py # 光线追迹主循环，调度以上三个模块 ├── config/ │ └── setup.json # 镜面参数、太阳位置、光线数量配置 ├── postprocess/ │ ├── flux_plot.py # 能流分布可视化 │ └── efficiency_report.py # 光学效率统计输出 ├── tests/ │ ├── test_sun.py # 太阳采样分布测试 │ └── test_mirror.py # 镜面交点求解测试 └── examples/ └── simple_paraboloid.py # 单碟式聚光的示例脚本模块之间的依赖方向是单向的：tracer.py调用sun.py生成光线，调用mirror.py求交点，调用receiver.py做统计，不互相循环引用。配置和代码分离是这套程序做得最到位的地方——改镜面焦距不用动代码，改 JSON 文件就行了。
3.2 追迹主循环的执行流程
主循环在逻辑上就是两层嵌套：外层遍历每条光线，内层遍历每个反射表面。实际运行时为了性能会向量化，但理解思路看伪代码更清楚：
def run_tracing(config): sun_vec = compute_sun_position(config['date'], config['latitude']) rays = sample_sun_directions(config['num_rays'], sun_vec) hit_count = 0 for i, ray_dir in enumerate(rays): # 对每个反射镜面求交点 for mirror in config['mirrors']: hit_point = intersect_paraboloid(config['origin'], ray_dir, mirror['focal_length']) if hit_point is not None: # 计算镜面法向量，做反射 normal = compute_normal(hit_point, mirror) reflect_dir = ray_dir - 2 * np.dot(ray_dir, normal) * normal # 检查反射后的光线是否命中接收器 if hit_receiver(reflect_dir, config['receiver']): hit_count += 1 break # 一条光线只和一面镜子交互一次 optical_efficiency = hit_count / len(rays) return optical_efficiency数据流是单向的：太阳位置参数进sun.py，光线数组进tracer.py，能流图进postprocess/。这个流程里最影响结果的两个参数是光线数量和镜面网格划分密度。光线太少，统计噪声大；网格分太密，每个格子的光线落点稀疏，能流图会有大量空洞。
4. 运行配置与参数调优：从拿到代码到跑出第一张能流图
4.1 环境准备与依赖清单
跑这套程序不需要重型依赖，Python 3.8 以上就能跑通。核心库只有 numpy、matplotlib、json 三个。我是在虚拟环境里直接pip install numpy matplotlib就完事了，没有遇到版本冲突的问题。如果要用 GPU 加速版本，会额外用到 cupy——不过那是进阶玩法，先用 CPU 把流程跑通更实际。
4.2 修改 JSON 配置跑通示例
程序包里自带的示例是单碟式抛物面聚光，下面是简化后的配置文件结构：
{ "sun": { "date": "2024-03-21", "latitude": 36.0, "half_angle_rad": 0.00465, "num_rays": 500000 }, "paraboloid": { "focal_length": 2.0, "aperture_radius": 3.0, "reflectivity": 0.93 }, "receiver": { "type": "cylinder", "radius": 0.06, "height": 0.2, "grid_size": [200, 100] }, "output": { "flux_map": "results/flux_distribution.npy", "report": "results/optical_efficiency.txt" } }跑起来就一条命令：python main.py --config config/setup.json。程序会输出两样东西，一个 .npy 格式的能流矩阵，一个光学效率报告。我第一次跑的时候，最震惊的是 50 万条光线在 CPU 上跑完只要十几秒——这说明程序在向量化上下了功夫，不是逐条光线循环的写法。
4.3 关键参数怎么调
焦距 f 和接收器半径是这套程序最敏感的两个参数。焦距决定焦点位置，接收器半径决定能流密度峰值和拦截率之间的权衡。接收器半径调大，光学效率上升但峰值能流下降；焦距调偏一点，能流分布会明显不对称。我一般习惯的做法是固定其他参数，先扫一遍焦距，再扫接收器半径，找出效率最高、能流分布最均匀的工作点。程序本身没有做自动化参数扫描，但配置是 JSON 文件，写个 shell 循环批量改配置再跑起来很容易。
5. 常见问题与排查：光线追迹程序踩坑记录
5.1 能流分布中心出现黑洞
现象: 接收器中心网格的能流值反而比外围低，形成一圈亮环。 原因: 光线追迹程序在计算镜面反射时，法向量方向搞反了——法向量指向镜面内侧而不是外侧，导致反射光线背向焦点飞出去。 解决: 我在compute_normal()里加了一个检查，确保法向量和光线入射方向的夹角小于 90 度，否则手动翻转法向量方向。加了之后黑洞立刻消失。
5.2 光学效率超过 100%
现象: 计算出的光学效率达到 1.05，物理上不可能。 原因: 能量统计时，同一束光线被程序记录了两次——一次在交点累计，一次在接收器网格累计，计数逻辑重叠了。 解决: 把效率计算从追迹循环里挪到循环外，单独用一个hit_count变量计数。血的教训：调试遍历时喜欢到处打印，打印代码里顺手多写了个+=1，后来养成习惯，计数类变量只在专属函数里修改。
5.3 追迹速度极慢，一条光线跑了整整一秒
现象: 光线数量只有 2 万条的时候，程序跑了快 20 分钟，完全不对劲。 原因: 每一条光线都在 Python 层遍历所有镜面网格点，循环嵌套太深，没有向量化。前面讲到的那套代码，一开始是逐点用for循环求交点，后来用 numpy 的数组运算整体求，速度提升了两个数量级。 解决: 把"光线数组 vs 镜面网格点"改成一次性的矩阵运算，用 numpy 的广播机制批量计算判别式和交点。写向量化代码时注意内存占用，50 万条光线乘 1000 个网格点就是 5 亿次浮点运算，单精度数组占 2GB 内存，预算不够就分批处理。
5.4 蒙特卡洛结果每次跑都不一样
现象: 同样的配置连续跑三次，光学效率分别是 0.71、0.75、0.77，波动大得没法用。 原因: 随机数没有固定种子。太阳采样是随机的，每次生成的光线方向分布都不同，统计结果自然有波动。 解决: 在最前面加np.random.seed(42)固定随机种子。注意不是加一次就完事——并行运算时每个线程都要固定自己的种子，否则子线程内部还是随机状态。从那以后，我每次跑仿真都会强制在代码第一行加固定种子。
5.5 接收器网格显示拉伸变形
现象: 能流分布图看起来是椭圆形的，但物理模型明明是对称的。 原因: 接收器网格的 x 和 y 间距不一致，绘图时没有等比例显示。 解决: 画图时用参数plt.axis('equal')强制横纵轴缩放一致。这个问题和追迹算法无关，但很容易误导人，特别是汇报时怎么看都觉得数据有问题。
6. 二次开发：把基础追迹改造成双轴跟踪模拟器
改造这套程序最有价值的切入点，是把它从"固定朝向的静态追迹"升级成"带双轴跟踪的动态模拟"。静态追迹只算一个太阳位置下的能流分布，实际电站的太阳位置每小时都在变，双轴跟踪系统能让反射镜始终对准太阳。改造思路分三步：
第一步，在sun.py里加一个时间变化函数。根据日期和当地时间算出太阳的高度角与方位角，用pysolar这样的天文算法库，输入经纬度和 UTC 时间直接输出太阳位置。把现在的compute_sun_position(date, latitude)扩展成compute_sun_position_series(date_range, latitude, longitude, time_step)，返回一个时间序列数组。
第二步，在tracer.py的循环里加镜面法向量的实时更新。双轴跟踪意味着每个时刻镜面朝向都随太阳位置变化，需要计算目标方位角和俯仰角，再把反射镜的旋转矩阵作用到镜面法向量上。代码大概长这样：
def update_mirror_orientation(mirror_normal, sun_azimuth, sun_elevation): """ 双轴跟踪：根据太阳方位角和高度角旋转反射镜 返回更新后的镜面法向量 """ # 先绕 z 轴旋转方位角，再绕 y 轴旋转高度角 rot_z = np.array([[np.cos(sun_azimuth), -np.sin(sun_azimuth), 0], [np.sin(sun_azimuth), np.cos(sun_azimuth), 0], [0, 0, 1]]) rot_y = np.array([[np.cos(sun_elevation), 0, np.sin(sun_elevation)], [0, 1, 0], [-np.sin(sun_elevation), 0, np.cos(sun_elevation)]]) return rot_y.dot(rot_z.dot(mirror_normal))旋转矩阵的乘法顺序不能换——先方位角后高度角，反了镜面会乱转。这段代码加进去后，把每隔 10 分钟追迹一次的结果汇总，就能生成一天的能流密度时间分布云图，这对判断储能容量和接收器热应力很有用。
第三步，把接收器改成一维追踪的线性菲涅尔模式，而不是点聚焦的碟式抛物面模式。线性系统里镜面是长条状分布的，每个条带的光线反射到同一根接收管上，追迹核心逻辑不用变，只要把抛物面方程换成柱面方程即可。
这套二次开发做完，原来只能算一个静态工况的程序，变成了能算全天动态性能的工具。最后我还做了一件事：把自己常改的变量统一放在config/app_settings.yaml里，包括太阳采样密度、追迹步长、接收器网格尺寸。从那以后，我每次做新工况仿真都强制走一遍「改配置 → 跑追迹 → 查能流图和效率报告」的流程，不再在代码里临时改数字，少翻车了很多次。希望帮到你。
本文还有配套的精品资源，点击获取
