tdehaeze/dcm-metrology
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Add ray tracing

Browse Source
This commit is contained in:
Thomas Dehaeze 2022-06-01 23:14:33 +02:00
parent 6214a64b21
commit e6ca5f96b3
1 changed files with 301 additions and 176 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

477
dcm-metrology.org
View File

@ -64,7 +64,7 @@ In order to increase the accuracy of the metrology system, two problems are to b
<<sec:metrology_concept>>
** 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
<<sec:relation_crystal_xray>>
** 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)
<<matlab-dir>>
#+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)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+end_src
#+begin_src matlab :tangle no :noweb yes
<<m-init-path>>
#+end_src
#+begin_src matlab :eval no :noweb yes
<<m-init-path-tangle>>
#+end_src
#+begin_src matlab :noweb yes
<<m-init-other>>
#+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
<<sec:frame_deformations>>
** 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:
<<sec:get_plane_reflection>>
#+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