function ar_freeform_prism_extended()
%% 默认参数（单位：长度单位与度）
D_default = 10; % 入射光束宽度（在 y-z 平面范围内）
R_default = -60; % 镜面B曲率半径（负值表示凹面镜）
Raynum_default = 5; % 光线数量（生成 Raynum×Raynum 个光线）
K_default = 0; % CONIC 系数
A_coeffs_default = [0, 0, 0]; % 非球面高阶系数
M_default = 10; % 镜面B尺寸（径向范围）
LensPosA_default = 10; % 平面A位置 (x 坐标)，固定为 10
LensApertureA_default = 30; % 平面A口径（半径）
MirrorRot_deg_default = -23; % 镜面B绕 Y 轴旋转角度，默认 -23°
LensRotA_deg_default = 10; % 平面A绕 Y 轴旋转角度，默认 10°
PlaneASphRad_default = 100; % 平面A球面半径（非0时采用球面；正值表示凸面，负值表示凹面）
% 平面C（出口面）的参数
LensPosC_default = 5; % 平面C位置 (x 坐标)
LensApertureC_default = 20; % 平面C口径（半径）
LensRotC_deg_default = 70; % 平面C绕 Y 轴旋转角度
PlaneCShiftZ_default = 30; % 平面C在Z轴的平移
PlaneCSphRad_default = 100; % 平面C球面半径（0 表示平面，否则为球面）
%% 创建主界面
f = figure('Name', 'AR Prism Extended: A-plane refraction → B-mirror reflection → A-plane internal reflection → C-plane exit', ...
'NumberTitle', 'off', 'Position', [100, 100, 1600, 700], 'Color', [1 1 1]);
%% 创建绘图坐标轴（3D视图）
ax = axes('Parent', f, 'Position', [0.35, 0.1, 0.6, 0.8]);
view(ax, 3);
%% 创建参数输入控件
% 基础参数
uicontrol(f, 'Style', 'text', 'Position', [20,680,140,20], 'String', 'Beam Width (D):', 'HorizontalAlignment', 'left');
D_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,680,80,25], 'String', num2str(D_default));
uicontrol(f, 'Style', 'text', 'Position', [20,650,140,20], 'String', 'Mirror Radius (R):', 'HorizontalAlignment', 'left');
R_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,650,80,25], 'String', num2str(R_default));
uicontrol(f, 'Style', 'text', 'Position', [20,620,140,20], 'String', 'Ray Number:', 'HorizontalAlignment', 'left');
Raynum_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,620,80,25], 'String', num2str(Raynum_default));
uicontrol(f, 'Style', 'text', 'Position', [20,590,140,20], 'String', 'Conic (K):', 'HorizontalAlignment', 'left');
K_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,590,80,25], 'String', num2str(K_default));
uicontrol(f, 'Style', 'text', 'Position', [20,560,140,20], 'String', 'A2,A3,A4:', 'HorizontalAlignment', 'left');
A_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,560,120,25], 'String', num2str(A_coeffs_default));
uicontrol(f, 'Style', 'text', 'Position', [20,530,140,20], 'String', 'Mirror Size (M):', 'HorizontalAlignment', 'left');
M_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,530,80,25], 'String', num2str(M_default));
%% 平面A参数
uicontrol(f, 'Style', 'text', 'Position', [20,500,140,20], 'String', 'Plane A Pos (x):', 'HorizontalAlignment', 'left');
LensPosA_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,500,80,25], 'String', num2str(LensPosA_default));
uicontrol(f, 'Style', 'text', 'Position', [20,470,140,20], 'String', 'Plane A Aperture:', 'HorizontalAlignment', 'left');
LensApertureA_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,470,80,25], 'String', num2str(LensApertureA_default));
uicontrol(f, 'Style', 'text', 'Position', [20,440,140,20], 'String', 'Plane A Rot (deg):', 'HorizontalAlignment', 'left');
LensRotA_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,440,80,25], 'String', num2str(LensRotA_deg_default));
uicontrol(f, 'Style', 'text', 'Position', [20,410,140,20], 'String', 'Plane A SphRad:', 'HorizontalAlignment', 'left');
PlaneASphRad_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,410,80,25], 'String', num2str(PlaneASphRad_default));
%% 平面C参数
uicontrol(f, 'Style', 'text', 'Position', [20,380,140,20], 'String', 'Plane C Pos (x):', 'HorizontalAlignment', 'left');
LensPosC_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,380,80,25], 'String', num2str(LensPosC_default));
uicontrol(f, 'Style', 'text', 'Position', [20,350,140,20], 'String', 'Plane C Aperture:', 'HorizontalAlignment', 'left');
LensApertureC_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,350,80,25], 'String', num2str(LensApertureC_default));
uicontrol(f, 'Style', 'text', 'Position', [20,320,140,20], 'String', 'Plane C Rot (deg):', 'HorizontalAlignment', 'left');
LensRotC_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,320,80,25], 'String', num2str(LensRotC_deg_default));
uicontrol(f, 'Style', 'text', 'Position', [20,290,140,20], 'String', 'Plane C Shift (Z):', 'HorizontalAlignment', 'left');
PlaneCShiftZ_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,290,80,25], 'String', num2str(PlaneCShiftZ_default));
uicontrol(f, 'Style', 'text', 'Position', [20,260,140,20], 'String', 'Plane C SphRad:', 'HorizontalAlignment', 'left');
PlaneCSphRad_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,260,80,25], 'String', num2str(PlaneCSphRad_default));
%% 其他参数
uicontrol(f, 'Style', 'text', 'Position', [20,230,140,20], 'String', 'Mirror Rot (deg):', 'HorizontalAlignment', 'left');
MirrorRot_edit = uicontrol(f, 'Style', 'edit', 'Position', [170,230,80,25], 'String', num2str(MirrorRot_deg_default));
%% 绘制按钮
uicontrol(f, 'Style', 'pushbutton', 'Position', [20,180,80,30], 'String', '绘制', ...
'Callback', @(~,~)draw_scene(D_edit, R_edit, Raynum_edit, K_edit, A_edit, ...
M_edit, LensPosA_edit, LensApertureA_edit, ...
MirrorRot_edit, LensRotA_edit, ...
LensPosC_edit, LensApertureC_edit, LensRotC_edit, PlaneCShiftZ_edit, PlaneASphRad_edit, PlaneCSphRad_edit, ax));
%% 初始绘图
draw_scene(D_edit, R_edit, Raynum_edit, K_edit, A_edit, ...
M_edit, LensPosA_edit, LensApertureA_edit, ...
MirrorRot_edit, LensRotA_edit, ...
LensPosC_edit, LensApertureC_edit, LensRotC_edit, PlaneCShiftZ_edit, PlaneASphRad_edit, PlaneCSphRad_edit, ax);
end
%% 光线追迹及系统绘制函数
function draw_scene(D_edit, R_edit, Raynum_edit, K_edit, A_edit, ...
M_edit, LensPosA_edit, LensApertureA_edit, ...
MirrorRot_edit, LensRotA_edit, ...
LensPosC_edit, LensApertureC_edit, LensRotC_edit, PlaneCShiftZ_edit, PlaneASphRad_edit, PlaneCSphRad_edit, ax)
%% 读取参数
D = str2double(get(D_edit, 'String'));
R = str2double(get(R_edit, 'String'));
Raynum = str2double(get(Raynum_edit, 'String'));
K = str2double(get(K_edit, 'String'));
A_coeffs = str2num(get(A_edit, 'String'));
mirrorSize = str2double(get(M_edit, 'String'));
LensPosA = str2double(get(LensPosA_edit, 'String'));
LensApertureA = str2double(get(LensApertureA_edit, 'String'));
MirrorRot_deg = str2double(get(MirrorRot_edit, 'String'));
LensRotA_deg = str2double(get(LensRotA_edit, 'String'));
PlaneASphRad = str2double(get(PlaneASphRad_edit, 'String'));
LensPosC = str2double(get(LensPosC_edit, 'String'));
LensApertureC = str2double(get(LensApertureC_edit, 'String'));
LensRotC_deg = str2double(get(LensRotC_edit, 'String'));
PlaneCShiftZ = str2double(get(PlaneCShiftZ_edit, 'String'));
PlaneCSphRad = str2double(get(PlaneCSphRad_edit, 'String'));
if isnan(D) || D<=0, D = 10; set(D_edit, 'String', '10'); end
if isnan(R) || R>=0, R = -60; set(R_edit, 'String', '-60'); end
if isnan(Raynum) || Raynum<1, Raynum = 5; set(Raynum_edit, 'String', '5'); end
if isnan(K), K = 0; set(K_edit, 'String', '0'); end
if isempty(A_coeffs), A_coeffs = [0, 0, 0]; set(A_edit, 'String', '0,0,0'); end
if isnan(mirrorSize) || mirrorSize<=0, mirrorSize = D/2; set(M_edit, 'String', num2str(D/2)); end
if isnan(LensPosA), LensPosA = 10; set(LensPosA_edit, 'String', '10'); end
if isnan(LensApertureA) || LensApertureA<=0, LensApertureA = 40; set(LensApertureA_edit, 'String', '40'); end
if isnan(MirrorRot_deg), MirrorRot_deg = -23; set(MirrorRot_edit, 'String', '-23'); end
if isnan(LensRotA_deg), LensRotA_deg = 10; set(LensRotA_edit, 'String', '10'); end
if isnan(PlaneASphRad), PlaneASphRad = 0; set(PlaneASphRad_edit, 'String', '0'); end
if isnan(LensPosC), LensPosC = 5; set(LensPosC_edit, 'String', '5'); end
if isnan(LensApertureC) || LensApertureC<=0, LensApertureC = 40; set(LensApertureC_edit, 'String', '40'); end
if isnan(LensRotC_deg), LensRotC_deg = -10; set(LensRotC_edit, 'String', '-10'); end
if isnan(PlaneCShiftZ), PlaneCShiftZ = 0; set(PlaneCShiftZ_edit, 'String', '0'); end
if isnan(PlaneCSphRad), PlaneCSphRad = 0; set(PlaneCSphRad_edit, 'String', '0'); end
cla(ax); hold(ax, 'on'); grid(ax, 'on');
xlabel(ax, 'X'); ylabel(ax, 'Y'); zlabel(ax, 'Z');
title(ax, 'AR Prism Extended: A-plane refraction → B-mirror reflection → A-plane internal reflection → C-plane exit');
axis(ax, 'equal'); view(ax, 3);
%% 1) 平面 A（入口面）：旋转矩阵 & 绘制（球面或平面）
LensRotA = deg2rad(LensRotA_deg);
R_y_lensA = rotationY(LensRotA);
P0_lensA = R_y_lensA * [LensPosA; 0; 0]; % 平面 A 中心
n_lensA = R_y_lensA * [1; 0; 0]; % 平面 A 法向量（默认指向外部）
if abs(PlaneASphRad) > 1e-6
[XA_A, YA_A, ZA_A, ~] = drawSphericalSurface(LensPosA, LensApertureA, R_y_lensA, PlaneASphRad);
else
[XA_A, YA_A, ZA_A] = drawPlane(ax, LensPosA, LensApertureA, R_y_lensA, 'g', 0.3);
end
surf(ax, XA_A, YA_A, ZA_A, 'FaceAlpha', 0.3, 'EdgeColor', 'none', 'FaceColor', 'g');
LensRotC = deg2rad(LensRotC_deg);
R_y_lensC = rotationY(LensRotC);
if abs(PlaneCSphRad) > 1e-6
[XC, YC, ZC] = drawSphericalSurface(LensPosC, LensApertureC, R_y_lensC, PlaneCSphRad);
else
[ThetaC, rC] = meshgrid(linspace(0, 2*pi, 50), linspace(0, LensApertureC, 50));
XC = LensPosC * ones(size(ThetaC));
YC = rC .* cos(ThetaC);
ZC = rC .* sin(ThetaC);
coordsC = R_y_lensC * [XC(:)'; YC(:)'; ZC(:)'];
XC = reshape(coordsC(1, :), size(XC));
YC = reshape(coordsC(2, :), size(YC));
ZC = reshape(coordsC(3, :), size(ZC));
end
ZC = ZC + PlaneCShiftZ; % 加上Z方向平移
surf(ax, XC, YC, ZC, 'FaceAlpha', 0.3, 'EdgeColor', 'none', 'FaceColor', 'm');
% 定义平面C中心及法向量（用于交点计算）
P0_lensC = R_y_lensC * [LensPosC; 0; 0];
P0_lensC(3) = P0_lensC(3) + PlaneCShiftZ;
n_lensC = R_y_lensC * [1; 0; 0];
%% 3) 镜面 B（自由曲面）：旋转矩阵 & 绘制（翻转后实现会聚效果）
MirrorRot = deg2rad(MirrorRot_deg);
R_y = rotationY(MirrorRot);
invR_y = R_y';
r_vals = linspace(0, mirrorSize, 50);
phi_vals = linspace(0, 2*pi, 100);
[RR, PP] = meshgrid(r_vals, phi_vals);
% 翻转后：sag取负，镜面采用 x = -evenAsphere(r)
baseX = -evenAsphere(RR, R, K, A_coeffs);
baseY = RR .* cos(PP);
baseZ = RR .* sin(PP);
mirror_coords = R_y * [baseX(:)'; baseY(:)'; baseZ(:)'];
X_m = reshape(mirror_coords(1, :), size(baseX));
Y_m = reshape(mirror_coords(2, :), size(baseY));
Z_m = reshape(mirror_coords(3, :), size(baseZ));
if ~isreal(X_m) || ~isreal(Y_m) || ~isreal(Z_m)
X_m = real(X_m); Y_m = real(Y_m); Z_m = real(Z_m);
end
surf(ax, X_m, Y_m, Z_m, 'FaceAlpha', 0.5, 'EdgeColor', 'none', 'FaceColor', [0.8 0.8 0.8]);
colormap(ax, 'gray');
%% 4) 入射光线生成
x_start = 50;
[Y0, Z0] = meshgrid(linspace(-D/2, D/2, Raynum), linspace(-D/2, D/2, Raynum));
Y0 = Y0(:); Z0 = Z0(:);
N_rays = length(Y0);
P_in = [x_start * ones(N_rays, 1), Y0, Z0];
I_dir = repmat([-1, 0, 0], N_rays, 1);
%% 5) 光线追迹（共五步）
for i = 1:N_rays
p0 = P_in(i, :)';
d0 = I_dir(i, :)';
% 第一步：平面 A 入射折射 (空气 → 棱镜)
if abs(PlaneASphRad) > 1e-6
[P_A1, validA1, N_A1] = intersectSagSurface(p0, d0, LensPosA, PlaneASphRad, LensApertureA, R_y_lensA);
else
[P_A1, validA1] = intersectPlane(p0, d0, P0_lensA, n_lensA, LensApertureA);
N_A1 = n_lensA;
end
if ~validA1, continue; end
plot3(ax, [p0(1) P_A1(1)], [p0(2) P_A1(2)], [p0(3) P_A1(3)], 'b-');
T_dir = refractRay(d0, N_A1, 1, 1.56);
if isempty(T_dir), continue; end
% 第二步：折射光线与镜面 B 交点 & 反射
[P_B, T_dir_reflect, validB] = reflectOnMirrorB(P_A1, T_dir, R_y, invR_y, R, K, A_coeffs);
if ~validB, continue; end
plot3(ax, [P_A1(1) P_B(1)], [P_A1(2) P_B(2)], [P_A1(3) P_B(3)], 'g-');
% 第三步：内部反射 —— 镜面B反射光线与平面A内侧交点
if abs(PlaneASphRad) > 1e-6
[P_A2, validA2, N_A2] = intersectSagSurface(P_B, T_dir_reflect, LensPosA, PlaneASphRad, LensApertureA, R_y_lensA);
if ~validA2, continue; end
N_A2 = -N_A2; % 内部反射取反法向量
else
[P_A2, validA2] = intersectPlane(P_B, T_dir_reflect, P0_lensA, -n_lensA, LensApertureA);
end
if ~validA2, continue; end
plot3(ax, [P_B(1) P_A2(1)], [P_B(2) P_A2(2)], [P_B(3) P_A2(3)], 'r-');
% 第四步：平面A内部反射：计算反射方向
if abs(PlaneASphRad) > 1e-6
T_dir_A2 = reflectRay(T_dir_reflect, -N_A2);
else
T_dir_A2 = reflectRay(T_dir_reflect, -n_lensA);
end
% 第五步：出口：光线从平面A内部反射后，与平面C交点并发生折射（出口）
% 注意：此处出口面采用反向法向量（n_exit = -n_lensC），确保交点在棱镜内部
[P_C, validC] = intersectPlane(P_A2, T_dir_A2, P0_lensC, -n_lensC, LensApertureC);
if ~validC, continue; end
plot3(ax, [P_A2(1) P_C(1)], [P_A2(2) P_C(2)], [P_A2(3) P_C(3)], 'c-');
% 出口折射：从棱镜内到空气（n1=1.56, n2=1），使用 n_exit = -n_lensC
T_dir_out = refractRay(T_dir_A2, n_lensC, 1.56, 1);
if isempty(T_dir_out), continue; end
% 为了显示效果，延长一定距离
P_final = P_C + 2*abs(R)*T_dir_out;
plot3(ax, [P_C(1) P_final(1)], [P_C(2) P_final(2)], [P_C(3) P_final(3)], 'c-');
end
% 为确保平面A显示在最上层，重新绘制平面A
if abs(PlaneASphRad) > 1e-6
[XA_A, YA_A, ZA_A, ~] = drawSphericalSurface(LensPosA, LensApertureA, R_y_lensA, PlaneASphRad);
else
[XA_A, YA_A, ZA_A] = drawPlane(ax, LensPosA, LensApertureA, R_y_lensA, 'g', 0.3);
end
surf(ax, XA_A, YA_A, ZA_A, 'FaceAlpha', 0.3, 'EdgeColor', 'none', 'FaceColor', 'g');
hold(ax, 'off');
end
%% --- 辅助函数 ---
function [P_int, valid] = intersectPlane(p0, d, P_plane, n_plane, aperture)
denom = dot(d, n_plane);
P_int = [];
valid = false;
if abs(denom) < 1e-9, return; end
t = dot(P_plane - p0, n_plane) / denom;
if t < 0, return; end
P_int = p0 + t*d;
if norm(P_int(2:3) - P_plane(2:3)) > aperture, return; end
valid = true;
end
function [P_int, valid, N_int] = intersectSagSurface(p0, d, LensPos, sphRad, aperture, R_mat)
p_local = R_mat \ p0;
d_local = R_mat \ d;
F = @(t) (p_local(1) + t*d_local(1)) - (LensPos + sphRad - sqrt(max(sphRad^2 - ((p_local(2)+t*d_local(2)).^2 + (p_local(3)+t*d_local(3)).^2), 0)));
t0 = (LensPos - p_local(1));
try
t_int = fzero(F, t0);
catch
valid = false; P_int = []; N_int = [];
return;
end
if isempty(t_int) || ~isfinite(t_int) || t_int < 0
valid = false; P_int = []; N_int = [];
return;
end
P_local_int = p_local + t_int*d_local;
r_local = norm(P_local_int(2:3));
if r_local > aperture
valid = false; P_int = []; N_int = [];
return;
end
valid = true;
P_int = R_mat * P_local_int;
if r_local < 1e-6
N_local = [1; 0; 0];
else
fprime = r_local / sqrt(max(sphRad^2 - r_local^2, 1e-6));
N_local = [1; -fprime*(P_local_int(2)/r_local); -fprime*(P_local_int(3)/r_local)];
N_local = N_local / norm(N_local);
end
N_int = R_mat * N_local;
end
function [X_rot, Y_rot, Z_rot, N_dummy] = drawSphericalSurface(planePos, aperture, R_mat, sphRad)
[Theta, r] = meshgrid(linspace(0, 2*pi, 50), linspace(0, aperture, 50));
if sphRad > 0
sag = sphRad - sqrt(max(sphRad^2 - r.^2, 0));
elseif sphRad < 0
sag = sphRad + sqrt(max(sphRad^2 - r.^2, 0));
else
sag = zeros(size(r));
end
X = planePos + sag;
Y = r .* cos(Theta);
Z = r .* sin(Theta);
coords = R_mat * [X(:)'; Y(:)'; Z(:)'];
X_rot = reshape(coords(1, :), size(X));
Y_rot = reshape(coords(2, :), size(Y));
Z_rot = reshape(coords(3, :), size(Z));
N_dummy = [];
end
function [X_rot, Y_rot, Z_rot] = drawPlane(ax, planePos, aperture, R_mat, faceColor, alphaVal)
[Theta, r] = meshgrid(linspace(0, 2*pi, 50), linspace(0, aperture, 50));
X = planePos * ones(size(Theta));
Y = r .* cos(Theta);
Z = r .* sin(Theta);
coords = R_mat * [X(:)'; Y(:)'; Z(:)'];
X_rot = reshape(coords(1, :), size(X));
Y_rot = reshape(coords(2, :), size(Y));
Z_rot = reshape(coords(3, :), size(Z));
end
function Rm = rotationY(angleRad)
Rm = [cos(angleRad), 0, sin(angleRad);
0, 1, 0;
-sin(angleRad), 0, cos(angleRad)];
end
function R = reflectRay(I, N)
I = I/norm(I);
N = N/norm(N);
R = I - 2*(dot(I, N))*N;
R = R/norm(R);
end
function T = refractRay(I, n, n1, n2)
I = I/norm(I);
n = n/norm(n);
eta = n1/n2;
cosi = -dot(n, I);
k = 1 - eta^2*(1 - cosi^2);
if k < 0
T = [];
else
T = eta*I + (eta*cosi - sqrt(k))*n;
T = T/norm(T);
end
end
function x = evenAsphere(r, R, K, A_coeffs)
if (K==0) && all(A_coeffs==0)
x = r.^2 ./ (R*(1+sqrt(max(1 - (r.^2)/(abs(R)^2), 0))));
else
R_eff = abs(R);
if (K==0)
L = sqrt(max(1 - (r.^2)/(R_eff^2), 0));
x = r.^2 ./ (R_eff*(1+L));
x = x * sign(R);
else
disc = 1 - (1+K)*(r.^2)/(R_eff^2);
disc = max(disc, 0);
L = sqrt(disc);
x_base = r.^2 ./ (R_eff*(1+L));
x = x_base * sign(R);
for i = 1:length(A_coeffs)
x = x + A_coeffs(i)*r.^(2*(i+1))*sign(R);
end
end
end
end
function dfdR = evenAsphereDerivative(r, R, K, A_coeffs)
R_eff = abs(R);
if (K==0) && all(A_coeffs==0)
L = sqrt(max(1 - (r.^2)/(R_eff^2), 0));
dfdR = (2*r.*(1+L) + r.^3./(R_eff^2.*L)) ./ (R_eff*(1+L).^2);
dfdR(r==0) = 0;
else
disc = 1 - (1+K)*(r.^2)/(R_eff^2);
disc = max(disc, 0);
L = sqrt(disc);
dfdR_base = (2*r.*(1+L) + r.^3./(R_eff^2.*L)) ./ (R_eff*(1+L).^2);
dfdR_extra = 0;
for i = 1:length(A_coeffs)
dfdR_extra = dfdR_extra + (2*(i+1))*A_coeffs(i)*r.^(2*i+1);
end
dfdR = dfdR_base + dfdR_extra;
dfdR(r==0) = 0;
end
dfdR = dfdR * sign(R);
end
function val = localIntersection(invR_y, Pstart, d, t, R, K, A_coeffs)
P_global = Pstart + t*d;
P_local = invR_y * P_global;
r_local = norm(P_local(2:3));
% 翻转后的镜面B交点条件：P_local(1) + evenAsphere(r)=0
val = P_local(1) + evenAsphere(r_local, R, K, A_coeffs);
end
function [P_m, R_dir, valid] = reflectOnMirrorB(p0, d, R_y, invR_y, R, K, A_coeffs)
F = @(t) real(localIntersection(invR_y, p0, d, t, R, K, A_coeffs));
p_local0 = invR_y * p0;
guess = (p_local0(1) + evenAsphere(norm(p_local0(2:3)), R, K, A_coeffs)) / max(abs(d(1)), 1e-6);
t_low = 0;
t_high = guess;
while sign(F(t_low)) == sign(F(t_high)) && t_high < 1e6
t_high = t_high * 2;
end
if sign(F(t_low)) == sign(F(t_high))
valid = false; P_m = []; R_dir = [];
return;
end
try
tB = fzero(F, [t_low, t_high]);
catch
valid = false; P_m = []; R_dir = [];
return;
end
if tB < 1e-6
valid = false; P_m = []; R_dir = [];
return;
end
P_m = p0 + tB*d;
P_localB = invR_y * P_m;
r_localB = norm(P_localB(2:3));
if R < 0
fp = -evenAsphereDerivative(r_localB, R, K, A_coeffs);
else
fp = evenAsphereDerivative(r_localB, R, K, A_coeffs);
end
if r_localB > 1e-9
N_local = [1; -fp*(P_localB(2)/r_localB); -fp*(P_localB(3)/r_localB)];
else
N_local = [1; 0; 0];
end
N_local = N_local / norm(N_local);
N_global = R_y * N_local;
if dot(N_global, d) > 0
N_global = -N_global;
end
R_dir = reflectRay(d, N_global);
valid = true;
end
这段MATLAB代码实现了一个用于模拟和绘制扩展增强现实（AR）棱镜系统中光线追迹的功能。该系统包括一个平面A（入口面）、一个镜面B（自由曲面）、另一个平面A（内部反射面）和一个平面C（出口面）。以下是对代码功能的详细描述：
-
默认参数设置：定义了系统的各种默认参数，如入射光束宽度、镜面曲率半径、光线数量、圆锥系数、非球面高阶系数等。
-
创建主界面：创建一个图形窗口，设置其名称、位置和大小，并在窗口中添加了一个3D视图的坐标轴。
-
创建参数输入控件：在图形窗口中添加了多个文本框和编辑框，用于用户输入和修改系统的各种参数。
-
绘制按钮：添加一个“绘制”按钮，当用户点击该按钮时，会调用
draw_scene函数进行光线追迹和图形绘制。 -
初始绘图：在程序启动时，调用
draw_scene函数进行一次初始的光线追迹和图形绘制。 -
光线追迹及系统绘制函数
draw_scene：- 读取参数：从用户输入的编辑框中读取各种参数，并进行有效性检查和默认值设置。
- 绘制平面A和平面C：根据输入的参数，绘制平面A和平面C，可以是平面或球面。
- 绘制镜面B：根据输入的参数，绘制镜面B，实现光线的反射效果。
- 生成入射光线：在指定位置生成一组入射光线。
- 光线追迹：对每条入射光线进行五步追迹：
- 第一步：光线在平面A的入射折射。
- 第二步：折射光线与镜面B的交点及反射。
- 第三步：反射光线与平面A内侧的交点。
- 第四步：平面A内部反射的反射方向计算。
- 第五步：光线从平面A内部反射后，与平面C的交点及折射（出口）。
- 重新绘制平面A：为确保平面A显示在最上层，重新绘制平面A。
-
辅助函数：
intersectPlane：计算光线与平面的交点。intersectSagSurface：计算光线与球面的交点。drawSphericalSurface：绘制球面。drawPlane：绘制平面。rotationY：计算绕Y轴的旋转矩阵。reflectRay：计算光线的反射方向。refractRay：计算光线的折射方向。evenAsphere：计算非球面的矢高。evenAsphereDerivative：计算非球面矢高的导数。localIntersection：计算局部交点。reflectOnMirrorB：计算光线在镜面B上的反射。
通过这些功能，代码能够实现对AR棱镜系统中光线追迹的模拟和可视化，帮助用户理解光线在系统中的传播路径和反射折射效果。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
