SMS多表面同步透镜设计

一、设计原理

1、Snell定律的矢量形式

折射定律又称Snell定律，主要包括两个方面：一是入射光线、法线和折射光线共面，二是入射角和折射角满足以下关系：

$$n_{1}sin{\theta }_1=n_{2}sin{\theta }_2$$

式中，$n_1$是入射光所在介质的折射率，$n_2$是折射光所在介质的折射率，${\theta }_1$入射角，${\theta }_2$是折射角。在LED二次光学设计中，斯涅尔（Snell）定律是最为简单也是最为重要的定律自由曲面模型的计算主要是通过采用折射光线、入射光线和法线矢量之间的等量关系建立方程，并联立角度与自由曲面间的关系，求解出一系列离散点，最终确定曲面的模型。

所有光线都遵循Snell定律，接下来讨论一下Snell定律的矢量形式,如下图所示：

原始的折射定律可以写成如下的矢量等式：

$$n^{\prime }=(Out\times N)=n\ast (In\times N)$$

式中$Out$和$In$分別是折射率为$n^{\prime }$和$n$时入射光线和折射光线的单位矢量，对上式进行化简，则有

$$Out\times N=In\times N$$

对上式进行变换可得

$$(Out-In)\ast N=0$$

在上式中，因为$(Ou{t}^{\prime }-I{n}^{\prime })$与$N$不可能为零，所以则有这两个矢量是共线或是平行的，于是上式可以被写成如下形式：

$$Out-In=K\ast N$$

式中$K$为待定常数。

给上式的等式两边同时点乘$N$，则有

$$K=Out\ast N-In\ast N=n^{\prime }cos({\theta }_1)-ncos({\theta }_2)$$

对上式分析可知：当$n^{\prime }>n$时，$K>0$，说明矢量$Ou{t}^{\prime }$和$I{n}^{\prime }$是同向的，反之，则是反向。通常情况下，在己知两种介质的折射率及入射角需要我们求解折射角时，此时上式可以化为：

$$K=\sqrt{n^{\prime 2}-n^2+ncos{{\theta }_1}^2}-ncos{\theta }_1$$

推导获得矢量形式：

$$n^{\prime }In-nOut=(n^{\prime }In-nOut)\ast N^2$$

将上式展开：

$$n^{\prime }In-nOut=\sqrt{(n^{\prime }|In|\ast |N|cos{\theta }_2-n|Out|\ast |N|cos{\theta }_1{)}^2}\ast N$$

通过进一步化简，得到最终形式：

$$\sqrt{n^{\prime 2}+n^2-2n^{\prime }n(Out\ast In)}\ast N=nOu{t}^{\prime }-n^{\prime }In$$

上式即为Snell定律的矢量形式。

2、多表面同步设计法

同步多表面设计法又称 Minano-Benitez设计方法、SMS(simultaneously multiple surface)设计法，来源于其允许对多个光学表面同时进行设计的性质。

一般的非成像光学设计方法都是以点光源为基础，不考虑光源的尺寸，这可以大大降低设计的难度，但在实际应用中多采用扩展光源，从而导致实际效果与设计预期效果存在偏差。但是在 SMS 设计方法中，光学表面按顺序的反射或折射全部来自于扩展光源端点发出的边缘光线。在光源尺寸确定的情况，保证光展恒定的前提下，由起点顺序地按照法线方向应用反射定律或者折射定律，在计算机程序的辅助下可以得到待定光学表面上的参数，由样条曲线的拟合最终确定曲面形状。

根据透镜的上顶点和光线经过该顶点后的路径，利用折射定律可以确定光线在透镜下表面的折射点及其法线方向。

接着，根据等光程原则，从接收面的另一侧作为光线的起点，以透镜下表面的折射点为参照，求解出透镜上表面的另一个关键点。重复这一操作，可以确定透镜的多个关键节点。最终，通过使用NURBS曲线拟合，确定透镜的表面形状。

根据透镜下表面的第一个点及其法线方向，用圆弧来拟合透镜下表面。然后，使用等光程方法及关键节点的计算方法，确定透镜上表面的形状。

3、等光程原理

在确定光线从透镜上顶点穿过透镜并经过折射到达接收面边缘后，可以计算出整个过程的光程。SMS的一个重要原则是等光程，后续的光线计算，如光线在透镜中的光程和光线在另一表面的折射点等，都需要依据等光程原则进行。

$$光程=\{\begin{array}{l}A{P}_1+n{P}_1{Q}_0+Q_{0}D\\ A{P}_2+n{P}_2{Q}_1+Q_{1}D\\ B{P}_1+n{P}_1{Q}_1+Q_{1}C\\ B{P}_0+n{P}_0{Q}_0+Q_{0}C\end{array}$$

二、设计目标

1. 通过SMS多表面同步设计方法设计一个能将扩展光源半径为$a=300mm$，在目标面距离为$H=3000mm$处，打出一个目标面光斑半径为$r=500mm$的透镜。

三、设计步骤

1、设计参数

$$\begin{array}{c}H=3000mm\\ r=500mm\\ a=300mm\end{array}$$

2、编写matlab程序，计算光学母线

main.m

%% 初始化
clc
clear

%% 初始参数设定
H = 3000; % 接收面到光源距离
r = 500; % 光斑半径
a = 300; % 光源半径
d = 2000; % 初始高度
n = 5; % 分段数量
l1 = 250; % B穿过PO后折线光线长度

%% 计算透镜上下表面坐标
[Q, P] = SMSchar(d, l1, H, r, a, n);

%% 绘制图像
plot(Q(:, 1), Q(:, 2), '-r');
axis equal;
grid on;
hold on;
plot(P(:, 1), P(:, 2), '-b');

%% 处理坐标数据并保存
z_r = zeros(length(Q), 1);
Q = [Q z_r];
z_r = zeros(length(P), 1);
P = [P z_r];
save('SMS上表面数据.txt', 'P', '-ascii');
save('SMS下表面数据.txt', 'Q', '-ascii');

SMSchar.m

function [Q, P] = SMSchar(d, l1, H, r, a, n)
    % 计算透镜上下表面的坐标
    % d 是 p0 的高度
    % l1 是 B 穿过 PO 后折线光线长度
    % H 是接收面到光源距离
    % a 是光源半宽
    % n 是分段数量
    % 圆心设在接收面中点
    
    % 参数定义
    p = zeros(n, 2);
    p(1, :) = [0 d];
    A = [-a H];
    B = [a H];
    C = [-r 0];
    D = [r 0];
    r1 = p(1, :) - B;
    n1 = 1; % 介质 1 的折射率
    n2 = 1.5; % 介质 2 的折射率
    
    F = zeros(n, 2);
    lb = zeros(n, 2);
    rb = zeros(n, 2);
    N = zeros(n, 2);
    
    rb(1, :) = r1 / norm(r1); % rb 表示从 B 点到 P 点的入射光线单位向量
    N(1, :) = [0 1]; % 法线单位向量，奇数为 P 点法线，偶数为 Q 点法线
    lb(1, :) = ref(rb(1, :), N(1, :));
    
    Q = zeros(n, 2);
    Q(1, :) = p(1, :) + l1 * lb(1, :); % Q1 点坐标
    K = n1 * (norm(B - p(1, :)) + norm(C - Q(1, :))) + n2 * norm(p(1, :) - Q(1, :)); % 求光程
    
    % 计算循环
    m = 1;
    for i = 1:n
        F(m, :) = C - Q(i, :); % QC 向量
        F(m, :) = F(m, :) / norm(F(m, :));
        k = -n1 * F(m, :) + n2 * lb(m, :);
        m = m + 1;
    
        rb(m, :) = (Q(i, :) - D) / norm(Q(i, :) - D); % DQ 光线
        N(m, :) = k;
        N(m, :) = N(m, :) / norm(N(m, :));
        lb(m, :) = ref(rb(m, :), N(m, :));
        L = (H - Q(i, 2)) / abs(lb(m, 2)); % lb 已归一化其 y 坐标即为角度 sin 值或者 -sin 值
        b = Q(i, :);
        a = Q(i, :) + L * lb(m, :);
        iter_count = 0;
        
        while true % 二分法求解 P 坐标
            iter_count = iter_count + 1;
            if iter_count > 1000
                break;
            end
            p(i + 1, :) = (a + b) / 2;
            Kn = n1 * (norm(A - p(i + 1, :)) + norm(D - Q(i, :))) + n2 * norm(p(i + 1, :) - Q(i, :));
            if Kn - K < -1e-6
                b = (a + b) / 2;
            elseif Kn - K > 1e-6
                a = (a + b) / 2;
            else
                break;
            end
        end
    
        rb(m + 1, :) = (p(i + 1, :) - B) / norm(p(i + 1, :) - B);
        F(m, :) = (A - p(i + 1, :)) / norm(A - p(i + 1));
        k = -n1 * F(m, :) + n2 * lb(m, :);
        N(m + 1, :) = k;
        N(m + 1, :) = N(m + 1, :) / norm(N(m + 1, :));
        lb(m + 1, :) = ref(rb(m + 1, :), N(m + 1, :));
        L = p(i + 1, 2) / abs(lb(m + 1, 2));
        b = p(i + 1, :);
        a = p(i + 1, :) + L * lb(m + 1, :);
        iter_count = 0;
        
        while true % 二分法求解 Q 坐标
            iter_count = iter_count + 1;
            if iter_count > 1000
                break;
            end
            Q(i + 1, :) = (a + b) / 2;
            Kn = n1 * (norm(B - p(i + 1, :)) + norm(C - Q(i + 1, :))) + n2 * norm(p(i + 1, :) - Q(i + 1, :));
            if Kn - K < -1e-6
                b = (a + b) / 2;
            elseif Kn - K > 1e-6
                a = (a + b) / 2;
            else
                break;
            end
        end
        m = m + 1;
    
        if Q(i + 1, 1) - Q(i, 1) > 0
            o = 1:n;
            ind = (o > i);
            Q(ind, :) = [];
            p(ind, :) = [];
            break;
        end
    end
    
    clear P
    num = 90;
    Qx = [];
    Qy = [];
    P = zeros(n, 2);
    rq = [];
    Oqy = [];
    x = [];
    y = [];
    dx = abs(2 * Q(1, 1)) / num;
    
    % 只保存 Q1 到 0 位置上的作为光学母线
    i = 1;
    x = 0:-dx:Q(i, 1);
    rq(1) = abs(Q(i, 1)) / abs(N(2, 1));
    Oqy(i) = Q(i, 2) + abs(N(2, 2)) * rq(i);
    y = -(rq(i)^2 - x.^2).^(1/2) + Oqy(i);
    Qx = [Qx x];
    Qy = [Qy y];
    
    h1 = length(x);
    
    % 为了求出 P1 到 P2 间所有点需要对 Q1 做对称 Q1' 来求出初始点
    i = 1;
    x = -Q(1, 1) - dx:-dx:Q(1, 1); % 为了防止重复，右端点不取
    rq(1) = -Q(1, 1) / -abs(N(2, 1));
    Oqy(1) = Q(1, 2) + N(2, 2) * rq(1);
    y = -(rq(1)^2 - x.^2).^(1/2) + Oqy(1);
    
    h = length(x);
    clear lb Q
    lb = zeros(n, 2);
    Q = zeros(n, 2);
    Q = [Qx', Qy'];
    
    % 计算从 Q1' 到 Q1 透射到 P1 到 P2 部分
    for g = 1:h
        nQg = [x(g) y(g)] - [0 Oqy(i)];
        nQg = nQg / norm(nQg);
        rbg = [x(g) y(g)] - D; % 入射光线为
        rbg = rbg / norm(rbg);
    
        lb(g, :) = ref(rbg, nQg);
    
        L = (H - y(g)) / abs(lb(g, 2));
        b = [x(g) y(g)];
        a = [x(g) y(g)] + L * lb(g, :);
        iter_count = 0;
        
        while true % 二分法求解 P 坐标
            iter_count = iter_count + 1;
            if iter_count > 1000
                disp(Kn - K);
                break;
            end
            Pi = (a + b) / 2;
            Kn = n1 * (norm(A - Pi) + norm(D - [x(g) y(g)])) + n2 * norm(Pi - [x(g) y(g)]);
            if Kn - K < -1e-6
                b = (a + b) / 2;
            elseif Kn - K > 1e-6
                a = (a + b) / 2;
            else
                break;
            end
        end
        
        % Px 和 Py 用来存储最终的 P 点坐标
        P(g, :) = Pi;
    end

    clear F N rbig
    F = zeros(n, 2);
    N = zeros(n, 2);
    rbig = zeros(n, 2);
    
    m = 1;
    
    % 重新使用 F 储存 PA 矢量，计算法线
    % 清空 N 储存法线
    for v = 0:n
        for i = 1:h
            rbig(i, :) = P(v * h + i, :) - B;
            rbig(i, :) = rbig(i, :) / norm(rbig(i, :));
            F(i, :) = (A - P(v * h + i, :)) / norm(A - P(v * h + i, :));
            N(i, :) = -n1 * F(i, :) + n2 * lb(i, :);
            N(i, :) = N(i, :) / norm(N(i, :));
            
            % lb((v + 1) * h + i, :) 表示第 v+1 次的第 i 个折射光线计算
            lb(i, :) = ref(rbig(i, :), N(i, :));
            
            L = P(v * h + i, 2) / abs(lb(i, 2));
            b = P(v * h + i, :);
            a = P(v * h + i, :) + L * lb(i, :);
            iter_count = 0;
            
            while true % 二分法求解 Q 坐标
                iter_count = iter_count + 1;
                if iter_count > 1000
                    break;
                end
                Q(h1 + v * h + i, :) = (a + b) / 2;
                Kn = n1 * (norm(B - P(v * h + i, :)) + norm(C - Q(h1 + v * h + i, :))) + n2 * norm(P(v * h + i, :) - Q(h1 + v * h + i, :));
                if Kn - K < -1e-6
                    b = (a + b) / 2;
                elseif Kn - K > 1e-6
                    a = (a + b) / 2;
                else
                    break;
                end
            end
    
            m = m + 1;
    
            F(i, :) = C - Q(h1 + v * h + i, :); % QC 向量
            F(i, :) = F(i, :) / norm(F(i, :));
            N(i, :) = -n1 * F(i, :) + n2 * lb(i, :);
            N(i, :) = N(i, :) / norm(N(i, :));
            rbig(i, :) = (Q(h1 + v * h + i, :) - D) / norm(Q(h1 + v * h + i, :) - D); % DQ 光线
            lb(i, :) = ref(rbig(i, :), N(i, :));
    
            L = (H - Q(h1 + v * h + i, 2)) / abs(lb(i, 2));
            b = Q(h1 + v * h + i, :);
            a = Q(h1 + v * h + i, :) + L * lb(i, :);
            iter_count = 0;
    
            while true % 二分法求解 P 坐标
                iter_count = iter_count + 1;
                if iter_count > 1000
                    break;
                end
                P((v + 1) * h + i, :) = (a + b) / 2;
                Kn = n1 * (norm(A - P((v + 1) * h + i, :)) + norm(D - Q(h1 + v * h + i, :))) + n2 * norm(P((v + 1) * h + i, :) - Q(h1 + v * h + i, :));
                if Kn - K < -1e-6
                    b = (a + b) / 2;
                elseif Kn - K > 1e-6
                    a = (a + b) / 2;
                else
                    break;
                end
            end
    
            m = m + 1;
        end
    end

ref.m

function r3 = ref(incidence, N)
     % 计算折射光线
     % incidence 是入射光线单位矢量，指向界面的方向
     % N 是法线单位矢量
     % 默认从 n1=1 的介质入射到 n2=1.5 的介质
     
     n1 = 1;
     n2 = 1.5;
     r1 = -incidence;
     n = N;
     sintheta1 = r1(1) * n(2) - r1(2) * n(1);
     sintheta2 = n1 * sintheta1 / n2;
     theta2 = -asin(sintheta2);
     r3 = -n;
     R = [cos(theta2), -sin(theta2); sin(theta2), cos(theta2)];
     r3 = R * r3';
     r3 = r3';
end

计算得到光学母线：

3、将计算好的光学母线数据导入SolidWorks，建立透镜模型

4、在SolidWorks中保存零件为.sat(R20)格式，并导入TracePro中

从辐照度分析图可以看出，该透镜的光效为$91.95\%$左右，整体均匀度大致在$50\%$以上，在中心位置处均匀度在$95\%$左右。

四、误差分析

目标光斑在边缘处呈现的光强较弱，或许是因为边缘处的光线根据等光程原理，在透镜中经过的光程更长，被透镜吸收的光强相比中心位置处更多。

五、总结

整体而言，本次设计成功实现了对扩展光源的SMS同步多表面透镜的设计，该透镜对于拓展光源打在目标面的光斑大小、光效表现良好，均满足设计目标。通过对折射定律、SMS同步多表面生成以及等光程原理的深入分析，对拓展光源的光学设计有了一定的理论基础。建立模型后，通过matlab计算光学母线，solidworks建立透镜模型以及tracepro进行光学仿真，展示了设计从理论到实践的转化过程。

六、参考

1. 张航, 严金华. 非成像光学设计[M]. 北京: 科学出版社, 2016.
2. 胡甜甜. 自由曲面菲涅尔面TIR透镜的设计及其应用[D].苏州大学,2020.DOI:10.27351/d.cnki.gszhu.2019.001237.