离轴反射系统真实光线追迹

在学习离轴反射式初始结构生成时，发现matlab和CodeV的交互速度很慢，这严重影响了算法运行速度。并且对于变焦系统来说，使用API获取到的镜面边界位置不准，所以希望自己编写一套光线追迹程序以解决这一问题。

以下内容均为以主镜作为光阑推导

变量定义

定义四面反射镜分别为$\text{M}_j$，$j$表示反射镜序号。反射镜中心点定义为$\text{O}_j$，其坐标为$(x_{O_j},y_{O_j})$。反射镜曲率半径为$R_j$，球面球心定义为$Q_j$，坐标$(x_{Q_j},y_{Q_j})$。

规定$O_1$为坐标原点，坐标为$(0,0)$。向右为x轴正方向，向上为y轴正方向。

求解主镜入射点及出射光线方向

首先计算球心位置$Q_j$。主镜的球心位置可以用$(R_1 \cdot \sin{\alpha_1},R_1 \cdot \cos{\alpha_1})+\vec{O_1}$来表示。其余反射镜需要使用到线性代数中旋转矩阵的知识。规定笛卡尔系中，逆时针旋转为正方向，则有：

$$\begin{align} \vec R(\theta) \vec r= \begin{bmatrix} \cos{\theta}&-\sin{\theta}\\ \sin{\theta}&\cos{\theta} \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} \end{align}$$

得到旋转后的向量。将$\overrightarrow {E_1Q_1} = -(\overrightarrow{OE_1}-\overrightarrow{OQ_1})$进行旋转后，乘以$d_j$，再加上$O_j$即可。

然后求解入射点，这一过程实际可以理解为求直线和圆的交点。规定入射光线的直线方程为$L_j: Ax+By+C=0$，圆的方程为$Q_j:(x-x_{Q_j})^2+(y-y_{Q_j})^2=R_j^2$

求解方法根据(99+ 封私信 / 80 条消息) 如何优雅的求直线与圆的交点？ - 知乎给出的公式计算得出。

1. 首先计算光线与镜片是否存在交点，如果存在则继续，如果不存在则退出计算，使用CodeV返回的光线追迹值。

   $\dfrac{|Ax_{Q_j}+By_{Q_j}+C|}{\sqrt{A^2+B^2}}\ge R_j$则认为不存在交点。

2. 计算圆心到直线的距离$\delta=\dfrac{Ax_{Q_j}+By_{Q_j}+C}{\sqrt{A^2+B^2}}$。

   计算直线的一个单位法向量$\vec {e_n} = \left ( \dfrac{A}{\sqrt{A^2+B^2}},\dfrac{B}{\sqrt{A^2+B^2}}\right )$

   如果O在$\vec {e_n}$指向的一侧，$\delta>0$，反之小于0，如果等于0则判定异常，退出调用codev。

   最后交点坐标为：

   $$\begin{align} \overrightarrow{OA} = \overrightarrow{OO_1}-\delta \cdot \vec e_n+\sqrt{R_j^2-|\delta|^2} \cdot \vec a \\ \overrightarrow{OB} = \overrightarrow{OO_1}-\delta \cdot \vec e_n-\sqrt{R_j^2-|\delta|^2} \cdot \vec a \end{align}$$

   $$\overrightarrow{OO_1} = (x_{Q_j},y_{Q_j}),\delta =\dfrac{Ax_{Q_j}+By_{Q_j}+C}{\sqrt{A^2+B^2}},\vec {e_n} = \left ( \dfrac{A}{\sqrt{A^2+B^2}},\dfrac{B}{\sqrt{A^2+B^2}}\right )$$

   $$\vec {a} = \left ( \dfrac{B}{\sqrt{A^2+B^2}},-\dfrac{A}{\sqrt{A^2+B^2}}\right )$$

得到的两个交点中，距离$O_j$的欧式距离最近的即为入射点。

最后求解出射光线。

这里将主镜单独拿出来讨论主要是为了定义第一条入射光线的方向。直线的方向矢量为$\vec v = (-B,A)$，规定$-B$必须为正数（即入射光线方向必须向右）。

入射光线的方向单位向量可以表示为$\vec v= \left ( \dfrac{-B}{\sqrt{A^2+B^2}},\dfrac{A}{\sqrt{A^2+B^2}}\right )$

定义$\overrightarrow {E_1Q_1} = -(\overrightarrow{OE_1}-\overrightarrow{OQ_1})$，规定$\overrightarrow {E_1Q_1}$的单位向量为$\vec e$

对于入射光线存在两种情况，第一种是$\vec v$和$\vec{e}$夹角的余弦值小于0（如上图），另一种是大于0（如下图）。

定义$\vec f$为出射光线单位矢量。

使用向量形式的反射定律$\vec f = \vec v-2(\vec v \cdot \vec e)\vec e$。

求解其余反射镜入射点及出射光线方向

入射点求解与主镜一样，反射光线求解将入射光线方向矢量替换为出射光线方向矢量。

Matlab代码编写

整体思路：

先求出每个反射镜的中心点$O_j$，然后定义光阑前的参考面位置$O_{ref}$，计算像面中心位置$O_{im}$。

然后计算主镜球心坐标，定义入射光线。

接着计算每个面上的入射点和反射光线方向矢量，重复计算上述过程。

最后输出每个镜片上对应光线的位置得到光线边缘。

function result = trace_single_config_offaxis(R, alpha, d, epd, fieldAnglesDeg, r)
%TRACE_SINGLE_CONFIG_OFFAXIS Single-configuration off-axis mirror ray tracing (YOZ plane).
% Inputs:
%   R              1x4 signed curvature radii [R1 R2 R3 R4]
%   alpha          1x4 mirror tilt increments in deg [a1 a2 a3 a4]
%   d              1x3 distances [d1 d2 d3]
%   epd            entrance pupil diameter
%   fieldAnglesDeg 1x3 field angles in deg
%
% Notes:
% - One configuration only.
% - O-point chain follows the user's specified equations.
% - If any ray has no intersection with a mirror, this function throws an
%   error as requested (instead of using CodeV returned ray data).

    validateattributes(R, {'numeric'}, {'vector', 'numel', 4, 'finite', 'real'});
    validateattributes(alpha, {'numeric'}, {'vector', 'numel', 4, 'finite', 'real'});
    validateattributes(d, {'numeric'}, {'vector', 'finite', 'real'});
    validateattributes(epd, {'numeric'}, {'scalar', 'positive', 'finite', 'real'});
    validateattributes(fieldAnglesDeg, {'numeric'}, {'vector', 'numel', 3, 'finite', 'real'});
    if nargin < 6
        r = [];
    end

    R = R(:).';
    alpha = alpha(:).';
    d = d(:).';
    fieldAnglesDeg = fieldAnglesDeg(:).';
    absR = abs(R);

    if ~(numel(d) == 3 || numel(d) == 4)
        error('trace_single_config_offaxis:InvalidD', ...
            'd must have 3 elements [d1 d2 d3] or 4 elements [d1 d2 d3 dImage].');
    end

    result = struct();
    result.valid = false;
    result.errorId = '';
    result.errorMsg = '';

    % Mirror vertex points from user-provided off-axis geometry chain.
    O1 = [0, 400];
    a12 = alpha(1) + alpha(2);
    a123 = a12 + alpha(3);
    O2 = O1 + d(1) * [sind(2 * alpha(1)), cosd(2 * alpha(1))];
    O3 = O2 + d(2) * [sind(2 * a12), cosd(2 * a12)];
    O4 = O3 + d(3) * [sind(2 * a123), cosd(2 * a123)];
    
    % 光阑前参考面.
    dRef = 400;

    Oref = O1 + [-dRef * sind(0), -dRef * cosd(0)];
    if numel(d) == 4
        dImage = d(4);
    else
        % d4 is replaced by BFL computed from radii and d1~d3.
        if ~isempty(r) && numel(r) >= 7
            dImage = calc_bfl_from_r(R, r(5:7));
        else
            dImage = calc_bfl_from_r(R, d(1:3));
        end
        if ~isfinite(dImage)
            result.errorId = 'trace_single_config_offaxis:InvalidBFL';
            result.errorMsg = 'Computed BFL is not finite.';
            return;
        end
    end
    alphaSum = sum(alpha);
    Oim = O4 + dImage * [sind(2 * alphaSum), cosd(2 * alphaSum)];
    O = [O1; O2; O3; O4];

    % Circle centers from signed radius and normal direction.
    % Signed radius controls which side the center lies on.
    Q = zeros(4, 2);
    if R(1) < 0
        Q(1, :) = [R(1) * sind(alpha(1)), R(1) * cosd(alpha(1))] + O1;
    else
        Q(1, :) = [-R(1) * sind(alpha(1)), -R(1) * cosd(alpha(1))] + O1;
    end

    for k = 2:4
        n_vec = (O(k - 1, :) - O(k, :))';
        n_vec = n_vec / norm(n_vec);

        R_vec = [cosd(alpha(k)), -sind(alpha(k)); sind(alpha(k)), cosd(alpha(k))];
        C_base = (R_vec * [n_vec(2); n_vec(1)]) * R(k);

        if k == 2 || k == 4
            signFactor = 1;
            if R(k) <= 0
                signFactor = -1;
            end
        else
            signFactor = -1;
            if R(k) <= 0
                signFactor = 1;
            end
        end

        C_vec = signFactor * C_base;
        Q(k, :) = O(k, :) + [C_vec(2), C_vec(1)];
    end

    % Use 3 representative pupil samples for each field: lower, chief, upper.
    pupilY = [-epd / 2, 0, epd / 2];
    zStart = O1(2);

    rayIdx = 0;
    totalRays = numel(fieldAnglesDeg) * numel(pupilY);
    rays(1, totalRays) = struct('fieldDeg', 0, 'pupilY', 0, 'points', zeros(6, 2), 'dirs', zeros(4, 2), 'imagePoint', [NaN, NaN], 'refPoint', [NaN, NaN], 'surfacePoints', zeros(6, 2));

    for fi = 1:numel(fieldAnglesDeg)
        fld = fieldAnglesDeg(fi);
        % 求入射光方向矢量
        d0 = normalize_vec([sind(fld), cosd(fld)]);

        for pi = 1:numel(pupilY)
            rayIdx = rayIdx + 1;

            % 入射光上的一个点
            p = [pupilY(pi), zStart];
            v = d0;
            refHit = caculate_plane_hit(p, v, O1, Oref);

            pts = zeros(6, 2);
            dirs = zeros(4, 2);
            pts(1, :) = p;

            for k = 1:4
                [hit, ok] = ray_circle_hit(p, v, Q(k, :), absR(k), O(k, :));
                if ~ok
                    result.errorId = 'trace_single_config_offaxis:NoIntersection';
                    result.errorMsg = sprintf(['No ray-mirror intersection at mirror %d (field=%.6f deg, pupilY=%.6f). ', ...
                        'Stop here and use CodeV returned ray tracing value.'], k, fld, pupilY(pi));
                    return;
                end
                
                E1Q1_vec = -([hit(2),hit(1)] - [Q(k,2),Q(k,1)]);

                v = caculate_reflect_ray(E1Q1_vec, v);
                pts(k + 1, :) = hit;
                dirs(k, :) = v;
                p = hit;
            end

            imageHit = caculate_plane_hit(p, v, O4, Oim);
            pts(6, :) = imageHit;

            % Six surface points aligned with ObtainData semantics:
            % S1=reference plane, S2~S5=M1~M4, S6=image plane.
            surfacePts = [refHit; pts(2:6, :)];

            rays(rayIdx).fieldDeg = fld;
            rays(rayIdx).pupilY = pupilY(pi);
            rays(rayIdx).points = pts;
            rays(rayIdx).dirs = dirs;
            rays(rayIdx).imagePoint = imageHit;
            rays(rayIdx).refPoint = refHit;
            rays(rayIdx).surfacePoints = surfacePts;
        end
    end

    result.valid = true;
    result.O = O;
    result.centers = Q;
    result.radii = absR;
    result.referencePlaneRef = [O1; Oref];
    result.imagePlaneRef = [O4; Oim];
    result.rays = rays;
    result.P_S = build_point_data(rays, fieldAnglesDeg, epd);

    % plot_trace_result(result);
end

function reflect_ray = caculate_reflect_ray(E1Q1_vec, v)
    nv = norm(v);
    ne = norm(E1Q1_vec);
    
    if ~isfinite(nv) || ~isfinite(ne) || nv < 1e-15 || ne < 1e-15
        reflect_ray = [NaN, NaN];
        return;
    end
    
    vin_xy = [v(2), v(1)] / nv;     % [x, y]
    n_xy = E1Q1_vec / ne;           % 单位法向 [x, y]
    
    % 镜面反射通式：v_out = v_in - 2*(v_in·n)*n
    vout_xy = vin_xy - 2 * dot(vin_xy, n_xy) * n_xy;
    
    no = norm(vout_xy);
    if ~isfinite(no) || no < 1e-15
        reflect_ray = [NaN, NaN];
        return;
    end
    
    vout_xy = vout_xy / no;
    reflect_ray = [vout_xy(2), vout_xy(1)];   % 转回 [y, x]
end

function plot_trace_result(result)
    fig = figure('Name', 'Single Config Off-Axis Trace', 'Color', 'w');
    ax = axes('Parent', fig);
    hold(ax, 'on');
    axis(ax, 'equal');
    grid(ax, 'on');
    xlabel(ax, 'Z (mm)');
    ylabel(ax, 'Y (mm)');
    title(ax, 'Off-axis single-configuration ray tracing (z-right, y-up)');

    O = result.O;
    Q = result.centers;
    radii = result.radii;
    referencePlaneRef = result.referencePlaneRef;
    imagePlaneRef = result.imagePlaneRef;

    % Draw mirror parent circles and key points.
    t = linspace(0, 2 * pi, 360);
    for k = 1:4
        yc = Q(k, 1) + radii(k) * cos(t);
        zc = Q(k, 2) + radii(k) * sin(t);
        plot(ax, zc, yc, '-', 'Color', [0.75, 0.75, 0.75], 'LineWidth', 0.8);
        plot(ax, O(k, 2), O(k, 1), 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 5);
        plot(ax, Q(k, 2), Q(k, 1), 'o', 'Color', [0.2, 0.5, 0.9], 'MarkerSize', 4);
        text(O(k, 2), O(k, 1), sprintf('  O%d', k), 'FontSize', 9, 'Color', 'k');
    end

    % Draw ray paths with color by field.
    % Draw segment-by-segment so reflected paths are always visible.
    fields = unique([result.rays.fieldDeg]);
    cmap = lines(numel(fields));
    for i = 1:numel(result.rays)
        fld = result.rays(i).fieldDeg;
        ci = find(fields == fld, 1);
        ptsAll = result.rays(i).points;

        refPt = result.rays(i).refPoint;
        pFirst = ptsAll(1, :);
        if all(isfinite(refPt)) && all(isfinite(pFirst))
            plot(ax, [refPt(2), pFirst(2)], [refPt(1), pFirst(1)], '--', 'Color', cmap(ci, :), 'LineWidth', 1.2);
            plot(ax, refPt(2), refPt(1), 'd', 'Color', cmap(ci, :), 'MarkerFaceColor', cmap(ci, :), 'MarkerSize', 4);
        end

        % Incoming + four reflected segments (p0->M1->M2->M3->M4).
        for s = 1:5
            pA = ptsAll(s, :);
            pB = ptsAll(s + 1, :);
            if all(isfinite(pA)) && all(isfinite(pB))
                plot(ax, [pA(2), pB(2)], [pA(1), pB(1)], '-', 'Color', cmap(ci, :), 'LineWidth', 1.8);
            end
        end

        % Mark mirror hit points (M1~M4).
        mirrorHits = ptsAll(2:5, :);
        validMirrorHits = all(isfinite(mirrorHits), 2);
        mirrorHits = mirrorHits(validMirrorHits, :);
        if ~isempty(mirrorHits)
            plot(ax, mirrorHits(:, 2), mirrorHits(:, 1), 'o', ...
                'Color', cmap(ci, :), 'MarkerFaceColor', cmap(ci, :), 'MarkerSize', 4);
        end

        imgPt = result.rays(i).imagePoint;
        if all(isfinite(imgPt))
            % Final segment from M4 to image plane.
            pLast = ptsAll(5, :);
            if all(isfinite(pLast))
                plot(ax, [pLast(2), imgPt(2)], [pLast(1), imgPt(1)], '--', 'Color', cmap(ci, :), 'LineWidth', 1.5);
            end
            plot(ax, imgPt(2), imgPt(1), 'x', 'Color', cmap(ci, :), 'LineWidth', 1.6, 'MarkerSize', 8);
        end
    end

    % Draw reference plane through Oref and normal to (O1->Oref) direction.
    O1ref = referencePlaneRef(1, :);
    Oref = referencePlaneRef(2, :);
    dirRef = [Oref(2) - O1ref(2), -(Oref(1) - O1ref(1))];
    nRef = norm(dirRef);
    if nRef > 1e-12
        dirRef = dirRef / nRef;
        allPts = [O; Q];
        for i = 1:numel(result.rays)
            rp = result.rays(i).points;
            allPts = [allPts; rp(all(isfinite(rp), 2), :)]; %#ok<AGROW>
            ip = result.rays(i).imagePoint;
            if all(isfinite(ip))
                allPts = [allPts; ip]; %#ok<AGROW>
            end
            rf = result.rays(i).refPoint;
            if all(isfinite(rf))
                allPts = [allPts; rf]; %#ok<AGROW>
            end
        end
        yRange = max(allPts(:, 1)) - min(allPts(:, 1));
        zRange = max(allPts(:, 2)) - min(allPts(:, 2));
        spanRef = max([yRange, zRange, 1]) * 0.8;
        pr1 = Oref - spanRef * dirRef;
        pr2 = Oref + spanRef * dirRef;
        plot(ax, [pr1(2), pr2(2)], [pr1(1), pr2(1)], '--', 'Color', [0.15, 0.45, 0.85], 'LineWidth', 1.8);
        plot(ax, Oref(2), Oref(1), 's', 'Color', [0.15, 0.45, 0.85], 'MarkerFaceColor', [0.15, 0.45, 0.85], 'MarkerSize', 5);
        text(Oref(2), Oref(1), '  Reference Plane', 'Color', [0.15, 0.45, 0.85], 'FontSize', 9, 'FontWeight', 'bold');
    end

    % Draw image plane through Oim and normal to (O4->Oim) direction.
    O4ref = imagePlaneRef(1, :);
    Oim = imagePlaneRef(2, :);
    dirImage = [Oim(2) - O4ref(2), -(Oim(1) - O4ref(1))];
    nDir = norm(dirImage);
    if nDir > 1e-12
        dirImage = dirImage / nDir;
        % Use plotted data extent instead of radius scale so the image plane is always visible.
        allPts = [O; Q];
        for i = 1:numel(result.rays)
            rp = result.rays(i).points;
            allPts = [allPts; rp(all(isfinite(rp), 2), :)]; %#ok<AGROW>
            ip = result.rays(i).imagePoint;
            if all(isfinite(ip))
                allPts = [allPts; ip]; %#ok<AGROW>
            end
        end
        yRange = max(allPts(:, 1)) - min(allPts(:, 1));
        zRange = max(allPts(:, 2)) - min(allPts(:, 2));
        span = max([yRange, zRange, 1]) * 0.8;
        p1 = Oim - span * dirImage;
        p2 = Oim + span * dirImage;
        plot(ax, [p1(2), p2(2)], [p1(1), p2(1)], '--', 'Color', [0.85, 0.1, 0.1], 'LineWidth', 2.0);
        plot(ax, Oim(2), Oim(1), 's', 'Color', [0.85, 0.1, 0.1], 'MarkerFaceColor', [0.85, 0.1, 0.1], 'MarkerSize', 5);
        text(Oim(2), Oim(1), '  Oim / Image Plane', 'Color', [0.85, 0.1, 0.1], 'FontSize', 10, 'FontWeight', 'bold');
    end

    legendLabels = arrayfun(@(v) sprintf('Field %.3f deg', v), fields, 'UniformOutput', false);
    dummy = gobjects(numel(fields), 1);
    for i = 1:numel(fields)
        dummy(i) = plot(ax, NaN, NaN, '-', 'Color', cmap(i, :), 'LineWidth', 1.5);
    end
    refPlaneDummy = plot(ax, NaN, NaN, '--', 'Color', [0.15, 0.45, 0.85], 'LineWidth', 1.8);
    imagePlaneDummy = plot(ax, NaN, NaN, '--', 'Color', [0.8, 0.2, 0.2], 'LineWidth', 1.4);
    imageSpotDummy = plot(ax, NaN, NaN, 'x', 'Color', [0.8, 0.2, 0.2], 'LineWidth', 1.4, 'MarkerSize', 8);
    legend(ax, [dummy; refPlaneDummy; imagePlaneDummy; imageSpotDummy], [legendLabels, {'Reference Plane', 'Image Plane', 'Image Points'}], 'Location', 'best');
end

function hit_plane = caculate_plane_hit(hit, v, P1, P2)
    % Ray line: p = hit + t * v (y-z coordinates).
    % Plane line in YOZ: through P2, normal to (P1 -> P2),
    % direction dImg = [dz, -dy].
    dImg = [P2(2) - P1(2), -(P2(1) - P1(1))];
    nv = norm(v);
    nd = norm(dImg);

    if ~isfinite(nv) || ~isfinite(nd) || nv < 1e-12 || nd < 1e-12
        hit_plane = [NaN, NaN];
        return;
    end

    rv = v / nv;
    rd = dImg / nd;

    % Solve 2x2 system hit + t*rv = P2 + s*rd explicitly.
    b = P2 - hit;
    detA = rv(2) * rd(1) - rv(1) * rd(2);
    if abs(detA) < 1e-12
        hit_plane = [NaN, NaN];
        return;
    end

    t = (b(2) * rd(1) - b(1) * rd(2)) / detA;
    hit_plane = hit + t * rv;
end

function P_S = build_point_data(rays, fieldAnglesDeg, epd)
    % Build macro-compatible point data vector (48x1) in this order for each S:
    % [Y R2 F1, Y R2 F3, Y R3 F1, Y R3 F3, Z R2 F1, Z R2 F3, Z R3 F1, Z R3 F3].
    P_S = NaN(48, 1);

    if numel(fieldAnglesDeg) < 3
        return;
    end

    f1 = fieldAnglesDeg(1);
    f3 = fieldAnglesDeg(3);
    r2 = epd / 2;   % +Y ray (upper edge)
    r3 = -epd / 2;  % -Y ray (lower edge)
    tol = max(1e-9, abs(epd) * 1e-12);

    rayR2F1 = pick_ray(rays, f1, r2, tol);
    rayR2F3 = pick_ray(rays, f3, r2, tol);
    rayR3F1 = pick_ray(rays, f1, r3, tol);
    rayR3F3 = pick_ray(rays, f3, r3, tol);

    idx = 1;
    for s = 1:6
        P_S(idx) = get_coord(rayR2F1, s, 1); idx = idx + 1; % Y R2 F1
        P_S(idx) = get_coord(rayR2F3, s, 1); idx = idx + 1; % Y R2 F3
        P_S(idx) = get_coord(rayR3F1, s, 1); idx = idx + 1; % Y R3 F1
        P_S(idx) = get_coord(rayR3F3, s, 1); idx = idx + 1; % Y R3 F3

        P_S(idx) = get_coord(rayR2F1, s, 2); idx = idx + 1; % Z R2 F1
        P_S(idx) = get_coord(rayR2F3, s, 2); idx = idx + 1; % Z R2 F3
        P_S(idx) = get_coord(rayR3F1, s, 2); idx = idx + 1; % Z R3 F1
        P_S(idx) = get_coord(rayR3F3, s, 2); idx = idx + 1; % Z R3 F3
    end
end

function raySel = pick_ray(rays, fldTarget, pupilTarget, tol)
    raySel = [];
    for i = 1:numel(rays)
        if abs(rays(i).fieldDeg - fldTarget) <= tol && abs(rays(i).pupilY - pupilTarget) <= tol
            raySel = rays(i);
            return;
        end
    end
end

function val = get_coord(raySel, s, coordIdx)
    if isempty(raySel)
        val = NaN;
        return;
    end
    sp = raySel.surfacePoints;
    if size(sp, 1) < s || size(sp, 2) < coordIdx
        val = NaN;
        return;
    end
    val = sp(s, coordIdx);
end

function [hit, ok] = ray_circle_hit(p0, v, Q, radius, O)
    % Solve |p0 + t*v - c|^2 = radius^2 for t > 0.
    [a,b,cc] = lineFromPointSlope(v(1) / v(2), p0(2), p0(1));
    denom = sqrt(a * a + b * b);
    delta = (a * Q(2) + b * Q(1) + cc) / denom;
    abs_delta = abs(delta);

    if abs_delta - abs(radius) >= 0
        hit = [NaN, NaN];
        ok = false;
        return;
    end

    if delta == 0
        hit = [NaN, NaN];
        ok = false;
        return;
    end

    a_vec = [b / denom, -a / denom];
    oo1_vec = [Q(2), Q(1)];
    en_vec = [a / denom, b / denom];

    halfChord = sqrt(radius ^ 2 - abs_delta ^ 2);
    OA_vec = oo1_vec - delta * en_vec + halfChord * a_vec;
    OB_vec = oo1_vec - delta * en_vec - halfChord * a_vec;

    if (norm(OA_vec - [O(2),O(1)]) < norm(OB_vec - [O(2),O(1)]))
        hit = [OA_vec(2),OA_vec(1)];
        ok = all(isfinite(hit));
    else
        hit = [OB_vec(2),OB_vec(1)];
        ok = all(isfinite(hit));
    end
end


function [a, b, c] = lineFromPointSlope(k, x0, y0)
% lineFromPointSlope 返回直线标准方程 ax + by + c = 0 的系数
%   给定直线的斜率 k 和直线上一点 (x0, y0)，计算标准形式的系数。
%   处理斜率无穷大（垂直线）的情况。
%   输入：
%       k  - 斜率（实数，可为 Inf 或 -Inf）
%       x0, y0 - 点的坐标
%   输出：
%       a, b, c - 直线方程 ax + by + c = 0 的系数

    if isinf(k)
        % 垂直线：x = x0  →  1*x + 0*y - x0 = 0
        a = 1;
        b = 0;
        c = -x0;
    else
        % 斜线：y - y0 = k(x - x0)  →  kx - y + (y0 - k*x0) = 0
        a = k;
        b = -1;
        c = y0 - k * x0;
    end
end

function u = normalize_vec(v)
    n = norm(v);
    if ~isfinite(n) || n < 1e-15
        u = [NaN, NaN];
    else
        u = v / n;
    end
    if(u(2)<0)
        u=-u;
    end
end

function BFL = calc_bfl_from_r(R, d123)
    R_BFL = d123(:).';
    if numel(R_BFL) ~= 3
        BFL = NaN;
        return;
    end

    den = 4 * R(2) * (R_BFL(1) - R_BFL(2)) * (R(4) + 2 * R_BFL(3)) ...
        - 2 * R(2) * R(3) * (R(4) + 2 * (R_BFL(1) - R_BFL(2) + R_BFL(3))) ...
        - 4 * R_BFL(1) * (2 * R_BFL(2) * (R(4) + 2 * R_BFL(3)) + R(3) * (R(4) - 2 * R_BFL(2) + 2 * R_BFL(3))) ...
        + 2 * R(1) * (R(2) * (R(3) - R(4) - 2 * R_BFL(3)) + 2 * R_BFL(2) * (R(4) + 2 * R_BFL(3)) + R(3) * (R(4) - 2 * R_BFL(2) + 2 * R_BFL(3)));

    num = R(4) * (R(1) * (-2 * R(3) * R_BFL(2) + R(2) * (R(3) - 2 * R_BFL(3)) + 2 * R(3) * R_BFL(3) + 4 * R_BFL(2) * R_BFL(3)) ...
        - 2 * (4 * R_BFL(1) * R_BFL(2) * R_BFL(3) + 2 * R(2) * (-R_BFL(1) + R_BFL(2)) * R_BFL(3) ...
        + 2 * R(3) * R_BFL(1) * (-R_BFL(2) + R_BFL(3)) + R(2) * R(3) * (R_BFL(1) - R_BFL(2) + R_BFL(3))));

    if abs(den) < 1e-12
        BFL = NaN;
    else
        BFL = num / den;
    end
end