From e6ca5f96b3e31c3e54875a3d09ab946068372a5d Mon Sep 17 00:00:00 2001 From: Thomas Dehaeze Date: Wed, 1 Jun 2022 23:14:33 +0200 Subject: [PATCH] Add ray tracing --- dcm-metrology.org | 477 +++++++++++++++++++++++++++++----------------- 1 file changed, 301 insertions(+), 176 deletions(-) diff --git a/dcm-metrology.org b/dcm-metrology.org index 9d82b4a..ebb337e 100644 --- a/dcm-metrology.org +++ b/dcm-metrology.org @@ -64,7 +64,7 @@ In order to increase the accuracy of the metrology system, two problems are to b <> ** Introduction :ignore: -The goal of the metrology system is to measure the distance and default of parallelism orientation between the first and second crystals +The goal of the metrology system is to measure the distance and default of parallelism between the first and second crystals. Only 3 degrees of freedom are of interest: - $d_z$ @@ -86,198 +86,111 @@ In order to measure the relative pose of the two crystals, instead of performing Three interferometers are used to measured the 3dof of interest for each crystals. Three additional interferometers are used to measured the relative motion of the metrology frame. + +In total, there are 15 interferometers represented in Figure [[fig:metrology_schematic]]. +The measurements are summarized in Table [[tab:interferometer_list]]. + #+name: tab:metrology_notations #+caption: Notations for the metrology frame #+attr_latex: :environment tabularx :width 0.4\linewidth :align cX #+attr_latex: :center t :booktabs t -| Notation | Meaning | -|----------+--------------------------| -| =d= | "Downstream": Positive X | -| =u= | "Upstream": Negative X | -| =h= | "Hall": Positive Y | -| =r= | "Ring": Negative Y | -| =f= | "Frame" | -| =1= | "First Crystals" | -| =2= | "Second Crystals" | +| *Notation* | *Meaning* | +|------------+--------------------------| +| =d= | "Downstream": Positive X | +| =u= | "Upstream": Negative X | +| =h= | "Hall": Positive Y | +| =r= | "Ring": Negative Y | +| =f= | "Frame" | +| =1= | "First Crystals" | +| =2= | "Second Crystals" | + +#+name: tab:interferometer_list +#+caption: List of Interferometer measurements +#+attr_latex: :environment tabularx :width 0.7\linewidth :align ccl +#+attr_latex: :center t :booktabs t +| *Number* | *Measurement* | *Description* | +|----------+---------------+-------------------------------------| +| 1 | $z_{1r,u}$ | First "Ring" Crystal, "upstream" | +| 2 | $z_{1r,c}$ | First "Ring" Crystal, "center" | +| 3 | $z_{1r,d}$ | First "Ring" Crystal, "downstream" | +| 4 | $z_{1h,u}$ | First "Hall" Crystal, "upstream" | +| 5 | $z_{1h,c}$ | First "Hall" Crystal, "center" | +| 6 | $z_{1h,d}$ | First "Hall" Crystal, "downstream" | +| 7 | $z_{2h,u}$ | Second "Hall" Crystal, "upstream" | +| 8 | $z_{2h,c}$ | Second "Hall" Crystal, "center" | +| 9 | $z_{2h,d}$ | Second "Hall" Crystal, "downstream" | +| 10 | $z_{2r,u}$ | Second "Ring" Crystal, "upstream" | +| 11 | $z_{2r,c}$ | Second "Ring" Crystal, "center" | +| 12 | $z_{2r,d}$ | Second "Ring" Crystal, "downstream" | +| 13 | $z_{mf,u}$ | Metrology Frame, "upstream" | +| 14 | $z_{mf,dr}$ | Metrology Frame, "downstream-ring" | +| 15 | $z_{mf,dh}$ | Metrology Frame, "downstream-hall" | #+name: fig:metrology_schematic #+caption: Schematic of the Metrology System [[file:figs/metrology_schematic.png]] -** Crystal's motion computation +** Computation of the relative pose between first and second crystals -From the raw interferometric measurements, the pose between the first and second crystals can be computed. +To understand how the relative pose between the crystals is computed from the interferometer signals, have a look at [[https://gitlab.esrf.fr/dehaeze/dcm-kinematics][this repository]] (=https://gitlab.esrf.fr/dehaeze/dcm-kinematics=). -First, Jacobian matrices can be used to convert raw interferometer measurements to axial displacement and orientation of the crystals and metrology frame. +Basically, Jacobian matrices are derived from the geometry and are used to convert the 15 interferometer signals to the *relative pose* of the primary and secondary crystals $[d_{h,z},\ r_{h,y},\ r_{h,x}]$ or $[d_{r,z},\ r_{r,y},\ r_{r,x}]$. -For the 311 crystals: +#+begin_note +The sign conventions for the relative crystal pose are: +- An increase of $d_{h,z}$ means the two crystals are further apart +- An increase of $r_{h,}$ means that the second crystals experiences a rotation around $y$ with respect to the primary crystal +- An increase of $r_{h,x}$ means that the second crystals experiences a rotation around $x$ with respect to the primary crystal +#+end_note -#+name: tab:311_crystals_notations_raw_interf -#+caption: Table caption -#+attr_latex: :environment tabularx :width 0.5\linewidth :align lX -#+attr_latex: :center t :booktabs t -| Notation | Description | -|----------+-----------------------------------| -| =um= | Metrology Frame - Upstream | -| =dhm= | Metrology Frame - Downstream Hall | -| =drm= | Metrology Frame - Downstream Ring | -|----------+-----------------------------------| -| =ur1= | First Crystal - Upstream Ring | -| =h1= | First Crystal - Hall | -| =dr1= | First Crystal - Downstream Ring | -|----------+-----------------------------------| -| =ur2= | First Crystal - Upstream Ring | -| =h2= | First Crystal - Hall | -| =dr2= | First Crystal - Downstream Ring | +The relative pose can be expressed as a function of the interferometers using the Jacobian matrices for the "hall" crystals: +\begin{equation} +\boxed{ +\begin{bmatrix} d_{h,z} \\ r_{h,y} \\ r_{h,x} \end{bmatrix} += +\bm{J}_{2h,s}^{-1} +\begin{bmatrix} z_{2h,u} \\ z_{2h,c} \\ z_{2h,d} \end{bmatrix} +- +\bm{J}_{1h,s}^{-1} +\begin{bmatrix} z_{1h,u} \\ z_{1h,c} \\ z_{1h,d} \end{bmatrix} +- +\bm{J}_{mf,s}^{-1} +\begin{bmatrix} z_{mf,u} \\ z_{mf,dh} \\ z_{mf,dr} \end{bmatrix} +} +\end{equation} -#+name: tab:311_crystals_notations_converted_interf -#+caption: Table caption -#+attr_latex: :environment tabularx :width 0.5\linewidth :align lX -#+attr_latex: :center t :booktabs t -| Notation | Description | -|----------+--------------------------------| -| =dzm= | Positive: increase of distance | -| =rym= | | -| =rxm= | | -|----------+--------------------------------| -| =dz1= | Positive: decrease of distance | -| =ry1= | | -| =rx1= | | -|----------+--------------------------------| -| =dz2= | Positive: increase of distance | -| =ry2= | | -| =rx2= | | +As well as for the "ring" crystals: +\begin{equation} +\boxed{ +\begin{bmatrix} d_{r,z} \\ r_{r,y} \\ r_{r,x} \end{bmatrix} += +\bm{J}_{2r,s}^{-1} +\begin{bmatrix} z_{2r,u} \\ z_{2r,c} \\ z_{2r,d} \end{bmatrix} +- +\bm{J}_{1r,s}^{-1} +\begin{bmatrix} z_{1r,u} \\ z_{1r,c} \\ z_{1r,d} \end{bmatrix} +- +\bm{J}_{mf,s}^{-1} +\begin{bmatrix} z_{mf,u} \\ z_{mf,dr} \\ z_{mf,dr} \end{bmatrix} +} +\end{equation} -#+begin_src latex :file schematic_sensor_jacobian_forward_kinematics_m.pdf -\begin{tikzpicture} - % Blocs - \node[block] (Js_inv) {$\bm{J}_{s,m}^{-1}$}; - - % Connections and labels - \draw[->] ($(Js_inv.west)+(-1.5,0)$) node[above right]{$\begin{bmatrix} u_{m} \\ dh_{m} \\ dr_{m} \end{bmatrix}$} -- (Js_inv.west); - \draw[->] (Js_inv.east) -- ++(1.5, 0) node[above left]{$\begin{bmatrix} d_{zm} \\ r_{ym} \\ r_{xm} \end{bmatrix}$}; -\end{tikzpicture} -#+end_src - -#+name: fig:schematic_sensor_jacobian_forward_kinematics_m -#+caption: Forward Kinematics for the Metrology frame -#+RESULTS: -[[file:figs/schematic_sensor_jacobian_forward_kinematics_m.png]] - -#+begin_src latex :file schematic_sensor_jacobian_forward_kinematics_1.pdf -\begin{tikzpicture} - % Blocs - \node[block] (Js_inv) {$\bm{J}_{s,1}^{-1}$}; - - % Connections and labels - \draw[->] ($(Js_inv.west)+(-1.5,0)$) node[above right]{$\begin{bmatrix} u_{r1} \\ h_1 \\ d_{r1} \end{bmatrix}$} -- (Js_inv.west); - \draw[->] (Js_inv.east) -- ++(1.5, 0) node[above left]{$\begin{bmatrix} d_{z1} \\ r_{y1} \\ r_{x1} \end{bmatrix}$}; -\end{tikzpicture} -#+end_src - -#+name: fig:schematic_sensor_jacobian_forward_kinematics_1 -#+caption: Forward Kinematics for the 1st crystal -#+RESULTS: -[[file:figs/schematic_sensor_jacobian_forward_kinematics_1.png]] - -#+begin_src latex :file schematic_sensor_jacobian_forward_kinematics_2.pdf -\begin{tikzpicture} - % Blocs - \node[block] (Js_inv) {$\bm{J}_{s,2}^{-1}$}; - - % Connections and labels - \draw[->] ($(Js_inv.west)+(-1.5,0)$) node[above right]{$\begin{bmatrix} u_{r2} \\ h_2 \\ d_{r2} \end{bmatrix}$} -- (Js_inv.west); - \draw[->] (Js_inv.east) -- ++(1.5, 0) node[above left]{$\begin{bmatrix} d_{z2} \\ r_{y2} \\ r_{x2} \end{bmatrix}$}; -\end{tikzpicture} -#+end_src - -#+name: fig:schematic_sensor_jacobian_forward_kinematics_2 -#+caption: Forward Kinematics for the 2nd crystal -#+RESULTS: -[[file:figs/schematic_sensor_jacobian_forward_kinematics_2.png]] - -Then, the displacement and orientations can be combined as follows: -\begin{align} - d_{z} &= + d_{z1} - d_{z2} + d_{zm} \\ - d_{r_y} &= - r_{y1} + r_{y2} - r_{ym} \\ - d_{r_x} &= - r_{x1} + r_{x2} - r_{xm} -\end{align} - -Therefore: -- $d_z$ represents the distance between the two crystals -- $d_{r_y}$ represents the rotation of the second crystal w.r.t. the first crystal around $y$ axis -- $d_{r_x}$ represents the rotation of the second crystal w.r.t. the first crystal around $x$ axis - -If $d_{r_y}$ is positive, the second crystal has a positive rotation around $y$ w.r.t. the first crystal. -Therefore, the second crystal should be actuated such that it is making a negative rotation around $y$ w.r.t. metrology frame. - -The Jacobian matrices are defined as follow: -#+begin_src matlab -%% Sensor Jacobian matrix for the metrology frame -J_m = [1, 0.102, 0 - 1, -0.088, 0.1275 - 1, -0.088, -0.1275]; - -%% Sensor Jacobian matrix for 1st "111" crystal -J_s_111_1 = [-1, -0.036, -0.015 - -1, 0, 0.015 - -1, 0.036, -0.015]; - -%% Sensor Jacobian matrix for 2nd "111" crystal -J_s_111_2 = [1, 0.07, 0.015 - 1, 0, -0.015 - 1, -0.07, 0.015]; -#+end_src - -Therefore, the matrix that gives the relative pose of the crystal from the 9 interferometers is: -#+begin_src matlab -%% Compute the transformation matrix -G_111_t = [-inv(J_s_111_1), inv(J_s_111_2), -inv(J_m)]; - -% Sign convention for the axial motion -G_111_t(1,:) = -G_111_t(1,:); -#+end_src - -#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*) -data2orgtable(G_111_t, ... - {'=dz= [nm]', '=rx= [nrad]', '=ry= [nrad]'}, ... - {'=ur1= [nm]', '=h1= [nm]', '=dr1= [nm]', '=ur2= [nm]', '=h2= [nm]', '=dr1= [nm]', '=um= [nm]', '=dhm= [nm]', '=drm= [nm]'}, ... - ' %.3f '); -#+end_src - -#+name: tab:transformation_matrix -#+caption: Transformation Matrix -#+attr_latex: :environment tabularx :width \linewidth :align cccccccccc -#+attr_latex: :center t :booktabs t :font \scriptsize -#+RESULTS: -| | =ur1= [nm] | =h1= [nm] | =dr1= [nm] | =ur2= [nm] | =h2= [nm] | =dr1= [nm] | =um= [nm] | =dhm= [nm] | =drm= [nm] | -|-------------+------------+-----------+------------+------------+-----------+------------+-----------+------------+------------| -| =dz= [nm] | -0.25 | -0.5 | -0.25 | -0.25 | -0.5 | -0.25 | 0.463 | 0.268 | 0.268 | -| =rx= [nrad] | 13.889 | 0.0 | -13.889 | 7.143 | 0.0 | -7.143 | -5.263 | 2.632 | 2.632 | -| =ry= [nrad] | 16.667 | -33.333 | 16.667 | 16.667 | -33.333 | 16.667 | 0.0 | -3.922 | 3.922 | - -From table [[fig:schematic_sensor_jacobian_forward_kinematics_2]], we can determine the effect of each interferometer on the estimated relative pose between the crystals. -For instance, an error on =dr1= will have much greater impact on =ry= than an error on =drm=. +Values of the matrices can be found in the document describing the kinematics of the DCM (see =https://gitlab.esrf.fr/dehaeze/dcm-kinematics=). * Relation Between Crystal position and X-ray measured displacement -** Setup +<> +** Introduction :ignore: -#+name: fig:calibration_setup -#+caption: Schematic of the setup -[[file:figs/calibration_setup.png]] +In this section, the impact of an error in the relative pose between the first and second crystals on the output X-ray beam is studied. -Detector: +This is very important in order to: +- link a measurement of the x-ray beam position to a default in the crystal position +- understand which pose default will have large impact on the output beam position/orientation +- calibrate the deformations of the metrology frame using and external metrology measuring the x-ray beam position/orientation -https://www.baslerweb.com/en/products/cameras/area-scan-cameras/ace/aca1920-40gc/ +In order to simplify the problem, the first crystal is supposed to be fixed (i.e. ideally positioned), and only the motion of the second crystal is studied. -Pixel size depends on the magnification used (1x, 6x, 12x). - -Pixel size of camera is 5.86 um x 5.86 um. -With typical magnification of 6x, pixel size is ~1.44um x 1.44um - -Frame rate is: 42 fps - -** Matlab Init :noexport:ignore: +** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) <> #+end_src @@ -290,8 +203,7 @@ Frame rate is: 42 fps bragg = pi/180*linspace(5,75,100); #+end_src -** Relation between second crystal motion and beam motion -*** Axial motion of second crystal +** Axial motion of second crystal Let's consider the relation between the $[y, z]$ motion of the beam and the motion of the second crystal $[z^\prime, R_{y^\prime}, R_{x^\prime}]$. #+name: fig:relation_dz_output_beam @@ -319,7 +231,7 @@ exportFig('figs/relation_vert_motion_crystal_beam.pdf', 'width', 'wide', 'height #+RESULTS: [[file:figs/relation_vert_motion_crystal_beam.png]] -*** Ry motion of second crystal +** Ry motion of second crystal \begin{equation} d_z = D_{\text{vlm}} d_{R_y^\prime} @@ -327,7 +239,7 @@ d_z = D_{\text{vlm}} d_{R_y^\prime} with $D_{\text{vlm}} \approx 10\,m$. -*** Rx motion of second crystal +** Rx motion of second crystal \begin{equation} d_y = 2 D_{\text{vlm}} \sin \theta \cdot d_{R_x^\prime} @@ -335,6 +247,173 @@ d_y = 2 D_{\text{vlm}} \sin \theta \cdot d_{R_x^\prime} +* Ray Tracing +** Matlab Init :noexport:ignore: +#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) +<> +#+end_src + +#+begin_src matlab :exports none :results silent :noweb yes +<> +#+end_src + +#+begin_src matlab :tangle no :noweb yes +<> +#+end_src + +#+begin_src matlab :eval no :noweb yes +<> +#+end_src + +#+begin_src matlab :noweb yes +<> +#+end_src + +** Definition of frame + +| | Position | Orientation | +|------------------+------------+-------------| +| input beam | =x1,y1,z1= | =s1= | +| primary mirror | =xp,yp,zp= | =np= | +| reflected beam | =x2,y2,z2= | =s2= | +| secondary mirror | =xs,ys,zz= | =ns= | +| output beam | =x3,y3,z3= | =s3= | + +#+begin_src matlab +theta = 5*pi/180; % [rad] +#+end_src + +#+begin_src matlab +%% Secondary crystal defaults +drx = 5*pi/180; % [rad] +dry = 0*pi/180; % [rad] +dz = 0; % [m] + +% Rotation matrix for drx +udrx = [cos(theta), 0, -sin(theta)]; +Rdrx = cos(drx)*eye(3)+sin(drx)*[0, -udrx(3), udrx(2); udrx(3), 0, -udrx(1); -udrx(2), udrx(1), 0] + (1-cos(drx))*(udrx'*udrx); + +% Rotation matrix for dry +Rdry = [ cos(dry), 0, sin(dry); ... + 0, 1, 0; ... + -sin(dry), 0, cos(dry)]; +#+end_src + +#+begin_src matlab +%% 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*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); +#+end_src + +#+begin_src matlab +%% Primary crystal plane +z = np; +y = [0,1,0]; +x = cross(y,z); +xtal1_rectangle = [pp + 0.02*y + 0.07*x; + pp - 0.02*y + 0.07*x; + pp - 0.02*y - 0.07*x; + pp + 0.02*y - 0.07*x]; + +%% Secondary crystal plane +z = ns; +y = [0,cos(drx),sin(drx)]; +x = cross(y,z); +xtal2_rectangle = [ps + 0.02*y + 0.07*x; + ps - 0.02*y + 0.07*x; + ps - 0.02*y - 0.07*x; + ps + 0.02*y - 0.07*x]; +#+end_src + +#+begin_src matlab +figure; +tiledlayout(2, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile(); +hold on; +plot3([p1(1), p2(1)],[p1(2), p2(2)], [p1(3), p2(3)]) +plot3([p2(1), p3(1)],[p2(2), p3(2)], [p2(3), p3(3)]) +plot3([p3(1), p3(1)+0.3*s3(1)],[p3(2), p3(2)+0.3*s3(2)], [p3(3), p3(3)+0.3*s3(3)]) +patch(xtal1_rectangle(:,1), xtal1_rectangle(:,2), xtal1_rectangle(:,3), 'k') +patch(xtal2_rectangle(:,1), xtal2_rectangle(:,2), xtal2_rectangle(:,3), 'k') +hold off; +view(0,0) +axis equal +xlim([-0.1, 0.15]) +zlim([-0.02, 0.01]) +grid off; +xlabel('X') +ylabel('Y') +zlabel('Z') + +ax2 = nexttile(); +hold on; +plot3([p1(1), p2(1)],[p1(2), p2(2)], [p1(3), p2(3)]) +plot3([p2(1), p3(1)],[p2(2), p3(2)], [p2(3), p3(3)]) +plot3([p3(1), p3(1)+0.3*s3(1)],[p3(2), p3(2)+0.3*s3(2)], [p3(3), p3(3)+0.3*s3(3)]) +patch(xtal1_rectangle(:,1), xtal1_rectangle(:,2), xtal1_rectangle(:,3), 'k') +patch(xtal2_rectangle(:,1), xtal2_rectangle(:,2), xtal2_rectangle(:,3), 'k') +hold off; +view(0,90) +axis equal +xlim([-0.1, 0.15]) +zlim([-0.02, 0.01]) +grid off; +xlabel('X') +ylabel('Y') +zlabel('Z') +#+end_src + +#+begin_src matlab +figure; +hold on; +plot3([p1(1), p2(1)],[p1(2), p2(2)], [p1(3), p2(3)]) +plot3([p2(1), p3(1)],[p2(2), p3(2)], [p2(3), p3(3)]) +plot3([p3(1), p3(1)+0.3*s3(1)],[p3(2), p3(2)+0.3*s3(2)], [p3(3), p3(3)+0.3*s3(3)]) +surf(xtal1_x, xtal1_y, xtal1_z) +patch(xtal2_rectangle(:,1), xtal2_rectangle(:,2), xtal2_rectangle(:,3), 'k') +hold off; +view(0,0) +axis equal +xlim([-0.1, 0.15]) +zlim([-0.02, 0.01]) +grid off; +#+end_src + +#+begin_src matlab +figure; +hold on; +plot3([p1(1), p2(1)],[p1(2), p2(2)], [p1(3), p2(3)]) +plot3([p2(1), p3(1)],[p2(2), p3(2)], [p2(3), p3(3)]) +plot3([p3(1), p3(1)+0.3*s3(1)],[p3(2), p3(2)+0.3*s3(2)], [p3(3), p3(3)+0.3*s3(3)]) +surf(xtal1_x, xtal1_y, xtal1_z) +patch(xtal2_rectangle(:,1), xtal2_rectangle(:,2), xtal2_rectangle(:,3), 'k') +hold off; +view(90,90) +axis equal +xlim([-0.1, 0.15]) +zlim([-0.02, 0.01]) +grid off; +#+end_src + * Deformations of the Metrology Frame <> ** Introduction :ignore: @@ -356,6 +435,22 @@ For each Bragg angle, the Fast Jacks are actuated to that the beam is at the cen Then, then position of the crystals as measured by the interferometers is recorded. This position is the wanted position for a given Bragg angle. +#+name: fig:calibration_setup +#+caption: Schematic of the setup +[[file:figs/calibration_setup.png]] + +Detector: + +https://www.baslerweb.com/en/products/cameras/area-scan-cameras/ace/aca1920-40gc/ + +Pixel size depends on the magnification used (1x, 6x, 12x). + +Pixel size of camera is 5.86 um x 5.86 um. +With typical magnification of 6x, pixel size is ~1.44um x 1.44um + +Frame rate is: 42 fps + + ** Simulations The deformations of the metrology frame and therefore the expected interferometric measurements can be computed as a function of the Bragg angle. @@ -1050,6 +1145,36 @@ colors = colororder; freqs = logspace(1, 3, 1000); #+END_SRC +** Intersection between Ray and Plane +:PROPERTIES: +:header-args:matlab+: :tangle matlab/src/get_plane_reflection.m +:header-args:matlab+: :comments none :mkdirp yes :eval no +:END: +<> + +#+begin_src matlab +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 - +#+end_src + +#+begin_src matlab +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_src + +#+begin_src matlab +end +#+end_src + * Bibliography :PROPERTIES: :UNNUMBERED: t