Working ray tracing with dz/dry/drx

matlab/src/getBeamPath.m

@@ -0,0 +1,72 @@
+function [results] = getBeamPath(theta, args)
+% getBeamPath -
+%
+% Syntax: [in_data] = getBeamPath(drx, dry, dz, theta, )
+%
+% Inputs:
+%    - drx, dry, dz, theta,  -
+%
+% Outputs:
+%    - in_data -
+
+arguments
+    theta
+    args.drx    (1,1) double {mustBeNumeric} = 0 % [rad]
+    args.dry    (1,1) double {mustBeNumeric} = 0 % [rad]
+    args.dz     (1,1) double {mustBeNumeric} = 0 % [m]
+end
+
+% Rotation matrix for drx
+udrx = [cos(theta), 0, -sin(theta)];
+Rdrx = cos(args.drx)*eye(3)+sin(args.drx)*[0, -udrx(3), udrx(2); udrx(3), 0, -udrx(1); -udrx(2), udrx(1), 0] + (1-cos(args.drx))*(udrx'*udrx);
+
+% Rotation matrix for dry
+Rdry = [ cos(args.dry), 0, sin(args.dry); ...
+         0,             1, 0; ...
+        -sin(args.dry), 0, cos(args.dry)];
+
+%% Input Beam
+p1 = [-0.1, 0, 0]; % [m]
+s1 = [ 1, 0, 0];
+
+%% Primary Mirror
+pp = [0, 0, 0]; % [m]
+np = [cos(pi/2-theta), 0, sin(pi/2-theta)];
+
+%% Reflected beam
+[p2, s2] = get_plane_reflection(p1, s1, pp, np);
+
+%% Secondary Mirror
+ps = pp ...
+     + 0.07*[cos(theta), 0, -sin(theta)] ... % x offset (does not matter a lot)
+     - np*10e-3./(2*cos(theta)) ... % Theoretical offset
+     + np*args.dz; % Add error in distance
+
+ns = [Rdry*Rdrx*[cos(pi/2-theta), 0, sin(pi/2-theta)]']'; % Normal
+
+%% Output Beam
+[p3, s3] = get_plane_reflection(p2, s2, ps, ns);
+
+%% Detector
+pd = [1, 0, 0]; % [m]
+nd = [-1, 0, 0];
+
+%% Get beam position on the decector
+p4 = get_plane_reflection(p3, s3, pd, nd);
+
+results = struct();
+% Beam position and orientations
+results.p1 = p1;
+results.s1 = s1;
+results.p2 = p2;
+results.s2 = s2;
+results.p3 = p3;
+results.s3 = s3;
+results.p4 = p4;
+% Mirrors/detector positions and orientations
+results.pp = pp;
+results.np = np;
+results.ps = ps;
+results.ns = ns;
+results.pd = pd;
+results.nd = nd;

matlab/src/get_plane_reflection.m

@@ -0,0 +1,16 @@
+function [p_reflect, s_reflect] = get_plane_reflection(p_in, s_in, p_plane, n_plane)
+% get_plane_reflection -
+%
+% Syntax: [p_reflect, s_reflect] = get_plane_reflection(p_in, s_in, p_plane, n_plane)
+%
+% Inputs:
+%    - p_in, s_in, p_plane, n_plane -
+%
+% Outputs:
+%    - p_reflect, s_reflect -
+
+p_reflect = p_in + s_in*dot(p_plane - p_in, n_plane)/dot(s_in, n_plane);
+s_reflect = s_in-2*n_plane*dot(s_in, n_plane);
+s_reflect = s_reflect./norm(s_reflect);
+
+end

matlab/src/plotBeamPath.m

@@ -0,0 +1,14 @@
+function [] = plotBeamPath(results)
+% plotBeamPath -
+%
+% Syntax: [in_data] = plotBeamPath(drx, dry, dz, theta, )
+%
+% Inputs:
+%    - drx, dry, dz, theta,  -
+%
+% Outputs:
+%    - in_data -
+
+plot3([results.p1(1), results.p2(1)],[results.p1(2), results.p2(2)], [results.p1(3), results.p2(3)])
+plot3([results.p2(1), results.p3(1)],[results.p2(2), results.p3(2)], [results.p2(3), results.p3(3)])
+plot3([results.p3(1), results.p4(1)],[results.p3(2), results.p4(2)], [results.p3(3), results.p4(3)])