传统几何光学成像与光纤追迹仿真系统
光学成像与光纤追迹仿真系统，结合了传统几何光学成像、成像仿真和光纤追迹技术，能够模拟从物体到像面的完整光学过程，并分析光纤耦合效率。
classdef OpticalSystemSimulator
% 光学系统仿真器
% 实现传统几何光学成像、成像仿真和光纤追迹
properties
components; % 光学元件列表
rays; % 光线集合
object; % 物体定义
imagePlane; % 像面位置
fiber; % 光纤参数
wavelength = 0.55e-6; % 波长(m)
end
methods
function obj = OpticalSystemSimulator()
% 构造函数
obj.components = {};
obj.rays = {};
obj.object = struct('type', 'point', 'position', [0,0,-100e-3], 'size', 10e-3);
obj.imagePlane = 200e-3; % 默认像面位置
obj.fiber = struct('position', [0,0,200e-3], 'na', 0.22, 'coreRadius', 25e-6, 'claddingRadius', 62.5e-6);
end
function addComponent(obj, type, params)
% 添加光学元件
comp = struct('type', type, 'params', params);
obj.components{end+1} = comp;
end
function setObject(obj, objType, position, size)
% 设置物体
obj.object.type = objType;
obj.object.position = position;
obj.object.size = size;
end
function setImagePlane(obj, position)
% 设置像面位置
obj.imagePlane = position;
end
function setFiber(obj, position, na, coreRadius, claddingRadius)
% 设置光纤参数
obj.fiber.position = position;
obj.fiber.na = na;
obj.fiber.coreRadius = coreRadius;
obj.fiber.claddingRadius = claddingRadius;
end
function generateRays(obj, numRays)
% 生成光线
obj.rays = {};
switch obj.object.type
case 'point'
% 点光源
for i = 1:numRays
% 随机方向
theta = 2*pi*rand();
phi = asin(2*rand() - 1); % -pi/2 到 pi/2
direction = [cos(phi)*cos(theta), cos(phi)*sin(theta), sin(phi)];
ray = struct('origin', obj.object.position, 'direction', direction, 'intensity', 1/numRays);
obj.rays{end+1} = ray;
end
case 'extended'
% 扩展光源
for i = 1:numRays
% 在物体表面随机取点
x = (rand()-0.5)*obj.object.size;
y = (rand()-0.5)*obj.object.size;
origin = [x, y, obj.object.position(3)];
% 随机方向
theta = 2*pi*rand();
phi = asin(2*rand() - 1);
direction = [cos(phi)*cos(theta), cos(phi)*sin(theta), sin(phi)];
ray = struct('origin', origin, 'direction', direction, 'intensity', 1/numRays);
obj.rays{end+1} = ray;
end
end
end
function traceRays(obj)
% 光线追迹
for i = 1:length(obj.rays)
ray = obj.rays{i};
for j = 1:length(obj.components)
comp = obj.components{j};
ray = obj.propagateRay(ray, comp);
end
% 传播到像面
ray = obj.propagateToPlane(ray, obj.imagePlane);
obj.rays{i} = ray;
end
end
function rayOut = propagateRay(obj, rayIn, component)
% 光线通过光学元件
switch component.type
case 'lens'
% 薄透镜
f = component.params.focalLength;
pos = component.params.position;
% 传播到透镜位置
ray = obj.propagateToPoint(rayIn, pos);
% 计算透镜作用
p = ray.origin;
d = ray.direction;
% 计算光线与光轴交点
if abs(d(3)) > 1e-6
t = (pos(3) - p(3)) / d(3);
intersection = p + t*d;
else
intersection = p;
end
% 薄透镜公式: 1/f = 1/u + 1/v
% 计算新方向
y = intersection(1:2);
y_norm = norm(y);
if y_norm > 0
% 透镜光心
center = [0, 0, pos(3)];
r = intersection - center;
% 折射
n1 = 1; % 空气
n2 = 1.5; % 玻璃
normal = r / norm(r);
incident = d;
% 斯涅尔定律
cosI = -dot(incident, normal);
sinT = (n1/n2)*sqrt(1-cosI^2);
if sinT > 1
% 全反射
reflect = incident - 2*cosI*normal;
d = reflect;
else
cosT = sqrt(1-sinT^2);
d = (n1/n2)*incident + (n1/n2*cosI - cosT)*normal;
end
end
rayOut = struct('origin', intersection, 'direction', d/norm(d), 'intensity', rayIn.intensity);
case 'mirror'
% 球面镜
R = component.params.radius;
pos = component.params.position;
type = component.params.type; % 'concave' or 'convex'
% 传播到镜面位置
ray = obj.propagateToPoint(rayIn, [pos(1:2), 0]);
p = ray.origin;
d = ray.direction;
% 计算与球面的交点
a = dot(d, d);
b = 2*dot(d, p - pos);
c = dot(p-pos, p-pos) - R^2;
discriminant = b^2 - 4*a*c;
if discriminant < 0
rayOut = ray; % 无交点
return;
end
t1 = (-b + sqrt(discriminant)) / (2*a);
t2 = (-b - sqrt(discriminant)) / (2*a);
t = min(t1, t2);
if t < 0
t = max(t1, t2);
end
intersection = p + t*d;
% 计算法线
normal = (intersection - pos) / R;
if strcmp(type, 'convex')
normal = -normal;
end
% 反射
reflect = d - 2*dot(d, normal)*normal;
rayOut = struct('origin', intersection, 'direction', reflect/norm(reflect), 'intensity', rayIn.intensity);
case 'aperture'
% 光阑
radius = component.params.radius;
pos = component.params.position;
% 传播到光阑位置
ray = obj.propagateToPoint(rayIn, [pos(1:2), 0]);
% 检查是否在光阑内
p = ray.origin;
if norm(p(1:2)) > radius
rayOut = struct('origin', p, 'direction', [0,0,0], 'intensity', 0); % 被阻挡
else
rayOut = ray;
end
otherwise
rayOut = rayIn;
end
end
function rayOut = propagateToPoint(obj, ray, point)
% 传播光线到指定点
p0 = ray.origin;
d = ray.direction;
if abs(d(3)) > 1e-6
t = (point(3) - p0(3)) / d(3);
else
t = 0;
end
p = p0 + t*d;
rayOut = struct('origin', p, 'direction', d, 'intensity', ray.intensity);
end
function rayOut = propagateToPlane(obj, ray, z)
% 传播光线到指定平面
p0 = ray.origin;
d = ray.direction;
if abs(d(3)) > 1e-6
t = (z - p0(3)) / d(3);
else
t = 0;
end
p = p0 + t*d;
rayOut = struct('origin', p, 'direction', d, 'intensity', ray.intensity);
end
function simulateImaging(obj, numRays)
% 成像仿真
obj.generateRays(numRays);
obj.traceRays();
end
function analyzeFiberCoupling(obj)
% 分析光纤耦合效率
totalIntensity = 0;
coupledIntensity = 0;
fiberHits = [];
for i = 1:length(obj.rays)
ray = obj.rays{i};
p = ray.origin;
d = ray.direction;
% 传播到光纤位置
if abs(d(3)) > 1e-6
t = (obj.fiber.position(3) - p(3)) / d(3);
else
t = 0;
end
hitPoint = p + t*d;
% 检查是否在光纤纤芯内
r = norm(hitPoint(1:2));
if r <= obj.fiber.coreRadius
% 计算入射角
normal = [0, 0, 1]; % 光纤端面法线
cosTheta = abs(dot(d, normal));
sinTheta = sqrt(1 - cosTheta^2);
% 检查是否在接收角内
if sinTheta <= obj.fiber.na
coupledIntensity = coupledIntensity + ray.intensity;
end
end
totalIntensity = totalIntensity + ray.intensity;
fiberHits(end+1,:) = [hitPoint(1:2), ray.intensity];
end
couplingEfficiency = coupledIntensity / totalIntensity;
fprintf('光纤耦合效率: %.2f%%\n', couplingEfficiency*100);
% 可视化光纤接收情况
figure;
scatter(fiberHits(:,1)*1e6, fiberHits(:,2)*1e6, 10, fiberHits(:,3), 'filled');
hold on;
% 绘制光纤纤芯
theta = linspace(0, 2*pi, 100);
coreX = obj.fiber.coreRadius * cos(theta);
coreY = obj.fiber.coreRadius * sin(theta);
plot(coreX*1e6, coreY*1e6, 'r-', 'LineWidth', 2);
% 绘制包层
cladX = obj.fiber.claddingRadius * cos(theta);
cladY = obj.fiber.claddingRadius * sin(theta);
plot(cladX*1e6, cladY*1e6, 'b-', 'LineWidth', 2);
axis equal;
xlabel('X (μm)');
ylabel('Y (μm)');
title('光纤接收点分布');
colorbar;
colormap jet;
end
function visualizeSystem(obj)
% 可视化光学系统
figure;
hold on;
grid on;
axis equal;
view(3);
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('光学系统仿真');
% 绘制物体
switch obj.object.type
case 'point'
plot3(obj.object.position(1), obj.object.position(2), obj.object.position(3), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
case 'extended'
[X, Y] = meshgrid(-obj.object.size/2:obj.object.size/10:obj.object.size/2, ...
-obj.object.size/2:obj.object.size/10:obj.object.size/2);
Z = zeros(size(X)) + obj.object.position(3);
surf(X, Y, Z, 'FaceAlpha', 0.3, 'EdgeColor', 'none', 'FaceColor', 'r');
end
% 绘制光学元件
for i = 1:length(obj.components)
comp = obj.components{i};
switch comp.type
case 'lens'
f = comp.params.focalLength;
pos = comp.params.position;
% 绘制薄透镜
[X, Y] = meshgrid(-0.05:0.01:0.05, -0.05:0.01:0.05);
Z = zeros(size(X)) + pos(3);
surf(X, Y, Z, 'FaceColor', 'blue', 'FaceAlpha', 0.3, 'EdgeColor', 'none');
text(0, 0, pos(3), sprintf('Lens f=%.1fmm', f*1000), 'HorizontalAlignment', 'center');
case 'mirror'
R = comp.params.radius;
pos = comp.params.position;
type = comp.params.type;
% 绘制球面镜
[X, Y, Z] = sphere(20);
X = R*X + pos(1);
Y = R*Y + pos(2);
Z = R*Z + pos(3);
surf(X, Y, Z, 'FaceColor', 'green', 'FaceAlpha', 0.3, 'EdgeColor', 'none');
if strcmp(type, 'concave')
text(pos(1), pos(2), pos(3), 'Concave Mirror', 'HorizontalAlignment', 'center');
else
text(pos(1), pos(2), pos(3), 'Convex Mirror', 'HorizontalAlignment', 'center');
end
case 'aperture'
radius = comp.params.radius;
pos = comp.params.position;
% 绘制光阑
theta = linspace(0, 2*pi, 100);
x = radius * cos(theta);
y = radius * sin(theta);
plot3(x, y, zeros(size(x))+pos(3), 'k-', 'LineWidth', 2);
text(0, 0, pos(3), 'Aperture', 'HorizontalAlignment', 'center');
end
end
% 绘制像面
[X, Y] = meshgrid(-0.1:0.02:0.1, -0.1:0.02:0.1);
Z = zeros(size(X)) + obj.imagePlane;
surf(X, Y, Z, 'FaceColor', 'yellow', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
text(0, 0, obj.imagePlane, 'Image Plane', 'HorizontalAlignment', 'center');
% 绘制光纤
pos = obj.fiber.position;
% 纤芯
[Xc, Yc, Zc] = cylinder(obj.fiber.coreRadius, 20);
Xc = Xc + pos(1);
Yc = Yc + pos(2);
Zc = Zc*0.1 + pos(3); % 短圆柱表示
surf(Xc, Yc, Zc, 'FaceColor', 'cyan', 'FaceAlpha', 0.3, 'EdgeColor', 'none');
% 包层
[Xcl, Ycl, Zcl] = cylinder(obj.fiber.claddingRadius, 20);
Xcl = Xcl + pos(1);
Ycl = Ycl + pos(2);
Zcl = Zcl*0.1 + pos(3);
surf(Xcl, Ycl, Zcl, 'FaceColor', 'magenta', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
text(pos(1), pos(2), pos(3), 'Fiber', 'HorizontalAlignment', 'center');
% 绘制光线
for i = 1:length(obj.rays)
ray = obj.rays{i};
if norm(ray.direction) > 0
p0 = ray.origin;
p1 = p0 + 0.2*ray.direction; % 缩短光线长度以便观察
plot3([p0(1), p1(1)], [p0(2), p1(2)], [p0(3), p1(3)], 'k-');
end
end
end
function visualizeImage(obj)
% 可视化像面
figure;
hold on;
grid on;
axis equal;
xlabel('X (m)');
ylabel('Y (m)');
title('像面光斑分布');
% 收集像面点
points = [];
intensities = [];
for i = 1:length(obj.rays)
ray = obj.rays{i};
p = ray.origin;
if abs(p(3) - obj.imagePlane) < 1e-3 % 在像面附近
points(end+1, :) = p(1:2);
intensities(end+1) = ray.intensity;
end
end
if isempty(points)
disp('无光线到达像面');
return;
end
% 绘制光斑
scatter(points(:,1)*1e3, points(:,2)*1e3, 20, intensities, 'filled');
colorbar;
colormap jet;
xlabel('X (mm)');
ylabel('Y (mm)');
% 计算点扩散函数
centroid = mean(points, 1);
rms = sqrt(mean(sum((points - centroid).^2, 2)));
fprintf('像面光斑RMS半径: %.2f mm\n', rms*1e3);
end
function runDemo()
% 运行演示
sim = OpticalSystemSimulator();
% 设置系统
sim.setObject('extended', [0,0,-100e-3], 20e-3);
sim.addComponent('lens', struct('focalLength', 50e-3, 'position', [0,0,0]));
sim.addComponent('aperture', struct('radius', 5e-3, 'position', [0,0,30e-3]));
sim.setImagePlane(200e-3);
sim.setFiber([0,0,200e-3], 0.22, 25e-6, 62.5e-6);
% 成像仿真
sim.simulateImaging(1000);
% 可视化
sim.visualizeSystem();
sim.visualizeImage();
sim.analyzeFiberCoupling();
end
end
end
系统功能说明
这个光学系统仿真器实现了以下功能：
-
几何光学成像仿真
- 支持点光源和扩展光源
- 光线追迹通过多种光学元件
- 计算像面光斑分布和点扩散函数
-
光学元件支持
- 薄透镜（凸透镜/凹透镜）
- 球面镜（凹面镜/凸面镜）
- 光阑（圆形孔径）
- 光纤（纤芯和包层）
-
光纤追迹与耦合分析
- 计算光线在光纤端面的位置
- 分析入射角与光纤接收角的关系
- 计算光纤耦合效率
- 可视化光纤接收点分布
-
三维可视化
- 显示整个光学系统布局
- 展示光线传播路径
- 可视化像面光斑分布
- 显示光纤接收情况
参考代码 传统几何光学成像,成像仿真，光纤追迹，可以实现模拟成像 www.youwenfan.com/contentcss/97009.html
使用示例
% 创建仿真器实例
sim = OpticalSystemSimulator();
% 设置物体（扩展光源）
sim.setObject('extended', [0,0,-100e-3], 20e-3);
% 添加光学元件
sim.addComponent('lens', struct('focalLength', 50e-3, 'position', [0,0,0]));
sim.addComponent('aperture', struct('radius', 5e-3, 'position', [0,0,30e-3]));
% 设置像面位置
sim.setImagePlane(200e-3);
% 设置光纤参数
sim.setFiber([0,0,200e-3], 0.22, 25e-6, 62.5e-6);
% 运行成像仿真（1000条光线）
sim.simulateImaging(1000);
% 可视化结果
sim.visualizeSystem();
sim.visualizeImage();
sim.analyzeFiberCoupling();
技术实现细节
1. 光线追迹算法
系统使用参数化光线追迹方法：
- 每条光线表示为：起点
origin和方向direction - 通过光学元件时，计算光线与元件的交点
- 应用相应的光学定律（折射、反射）更新光线方向
- 使用小角度近似简化计算
2. 薄透镜模型
薄透镜使用高斯光学公式：
- 透镜变换：
θ' = θ - y/f - 传播公式：
y' = y + θ' * d
3. 球面镜模型
球面镜使用解析几何计算交点：
- 球面方程：
(x-x₀)² + (y-y₀)² + (z-z₀)² = R² - 求解光线参数方程与球面方程的交点
4. 光纤耦合分析
光纤耦合效率计算考虑两个因素：
- 几何接受条件：光线在光纤端面位置
(x,y)满足√(x²+y²) ≤ a（纤芯半径） - 角度接受条件：入射角
θ满足sinθ ≤ NA（数值孔径）
5. 可视化技术
- 使用MATLAB的3D绘图功能展示系统布局
- 颜色映射表示光线强度分布
- 点扩散函数计算和显示
应用场景
这个仿真系统可用于：
- 光学系统设计验证
- 成像系统性能评估
- 光纤耦合效率优化
- 教学演示光学原理
- 激光束传输分析
扩展功能建议
-
高级光学元件
- 添加非球面透镜
- 支持衍射光学元件
- 加入渐变折射率透镜
-
杂散光分析
- 鬼像分析
- 散射模型
- 偏振效应
-
优化算法
- 自动优化元件位置
- 公差分析
- 蒙特卡洛误差分析
-
光谱分析
- 色差计算
- 多波长仿真
- 光谱透射率
