diff --git a/Doc/UserManual/Assemblies.tex b/Doc/UserManual/Assemblies.tex
index dbc9276c37b0f79ba9a32d6f32fe2f85ae3db0bd..2e30a03d298cb09b64e4512432ebf06439eaf0a0 100644
--- a/Doc/UserManual/Assemblies.tex
+++ b/Doc/UserManual/Assemblies.tex
@@ -35,13 +35,15 @@ of a mesoscopic size (nanometer to micrometer).
 In the following, all such inhomogeneities will be described as
 \E{particles} that are embedded in a material layer.
 
+\iffalse % minimal version
 \Work
 {Documentation on inter-particle correlations is under preparation.
 \BornAgain\ currently offers the same choice of models as \IsGISAXS.
 Therefore, for the moment we refer to the \IsGISAXS\ manual \cite{Laz08}.}
+\fi
 
 %\iffalse
-Our derivation of the DWBA (\cref{Sdwba}) has shown
+Our derivation of the DWBA (\cref{SGisas}) has shown
 that the distinction between refraction and reflection on the one hand
 and scattering on the other hand is not a fundamental one,
 but to a certain degree arbitrary:
@@ -150,7 +152,7 @@ The first one is the shape transform of the entire layer,
   = \chi_\text{m} \int_{z_{\il+1}}^{z_{\il}}\!\d z \int\!\d^2r_\plll\, \e^{i\v{q}\,\r}
   = (2\pi)^2 \chi_\text{m} d_\il \e^{i q_\perp \overline{z}}  \delta(\q_\plll) \sinc\left( q_\perp d_\il /2 \right),
 \end{equation}
-with $\overline{z} \coloneqq \frac{z_\il + z_{\il+1}}{2}$ and $d_\il \coloneqq z_\il - z_{\il+1}$.
+with $\overline{z} \coloneqq (z_\il + z_{\il+1})/2$ and $d_\il \coloneqq z_\il - z_{\il+1}$.
 Thanks to the delta function,
 this term only contributes to the direct beam and specular peak.
 \Warn{Since $\chi_\text{m}=0$ in BornAgain,
@@ -262,7 +264,7 @@ we encode the particle coordinates in a distribution function
 \index{Statistics!particle distribution}%
 \begin{equation}\label{EPDdeltas}
   \PD_{l}(\r,\v\tau,\v\tau')
-  = \sum_{ij} \delta(\r-(\R_{j\plll}-\R_{i\plll}))
+  \coloneqq \sum_{ij} \delta(\r-(\R_{j\plll}-\R_{i\plll}))
   \delta(\v\tau-\v{T}_i)\delta(\v\tau'-\v{T}_j).
 \end{equation}
 Since in the following we are concerned with only one layer,
@@ -279,7 +281,7 @@ The differential cross section~\cref{Eass3} becomes
 \end{equation}
 where the integration over $\r$ is restricted to the horizontal plane.
 
-The normalization of the distribution function $\PD$ is such that it's total
+The normalization of the distribution function $\PD$ is such that its total
 integral equals the square of the number of particles that contribute to the scattering:
 \begin{equation}
   \int\!\d^2r\, \int\! \d^{d_T}\tau\, \d^{d_T}\tau' \PD(\r,\v\tau,\v\tau')
@@ -534,32 +536,38 @@ particles is the coherent scattering intensity of~\cref{EICoherent}.
 
 The conditional interference function is defined as
 \begin{equation}
-  \SD(\q |\v\tau,\v\tau') \coloneqq 1 + \rho_S \int\! d^2r\, e^{i\q_\plll\r}\! \GD(\r|\v\tau,\v\tau').
+  \SD(\q |\v\tau,\v\tau')
+  \coloneqq 1 + \rho_S \int\!\d^2r\, \e^{i\q_\plll\r}\! \GD(\r|\v\tau,\v\tau').
 \end{equation}
 With the chosen normalization, $\rho_S \GD(\r|\v\tau,\v\tau') d^2r$ gives the probability of
-finding a particle of type $\v\tau'$ in the infinitesimal area $d^2r$ at the relative position $\r$ from
+finding a particle of type $\v\tau'$ in the infinitesimal area~$d^2r$
+at the relative position~$\r$ from
 a given (different) particle of type $\v\tau$.
 
-In the following subsections, the supported interference functions in \BornAgain\ will be discussed.
-With one exception, they all assume the R-T decoupling approximation, making the interference function
-independent of the particle types ($\v\tau$ and $\v\tau'$).
+In the following subsections,
+the supported interference functions in \BornAgain\ will be discussed.
+With one exception,
+they all assume the R-T decoupling approximation, making the interference function
+independent of the particle types ($\v\tau$ and~$\v\tau'$).
 
 %===============================================================================
 \subsection{One-dimensional lattice} \label{sec:sect:1dlattice}
 %===============================================================================
+
 For a perfect one-dimensional lattice along the x-axis with period $a$, the position
 correlation function is given by:
 \begin{equation}
   \rho_S\GD(\r) = \sum_{n\neq 0} \delta(x-na)\delta(y).
 \end{equation}
-The corresponding interference function then becomes
+Using standard relations for the Dirac comb,
+one obtains the interference function
 \begin{equation}
-  \SD(\q) = \frac{2\pi}{a}\sum_k \delta(q_x - \frac{2\pi k}{a}),
+ \SD(\q) = \frac{2\pi}{a}\sum_k \delta\left(q_x - \frac{2\pi k}{a}\right),
 \end{equation}
-where $2\pi /a$ is a basis vector for the reciprocal lattice.
+which has the reciprocal lattice period $2\pi/a$.
 
 For computational reasons in \BornAgain, the delta functions appearing in the interference function
-are replaced by distributions of a finite width $H(q_x-2\pi k/a)$. This ammounts to convoluting the
+are replaced by distributions of a finite width $H(q_x-2\pi k/a)$. This amounts to convoluting the
 previously given interference function with $H(q_x)$ or, equivalently, multiplying
 the position correlation function by the inverse Fourier image of $H(q_x)$, called the
 \index{Decay function}
@@ -567,7 +575,7 @@ the position correlation function by the inverse Fourier image of $H(q_x)$, call
 
 The interference function can then be written as:
 \begin{equation}
-  \SD(q_x) = \frac{1}{a}\sum_k H(q_x - \frac{2\pi k}{a}).
+  \SD(q_x) = \frac{1}{a}\sum_k H\left(q_x - \frac{2\pi k}{a}\right).
 \end{equation}
 
 \BornAgain\ currently supports the following types of one-dimensional decay functions
@@ -577,21 +585,24 @@ The interference function can then be written as:
     \hline
     \text{Name} & h(x) & H(q_x) \\
     \hline
-    \text{Cauchy} & e^{-|x|/\lambda} & \frac{2\lambda}{1+q_x^2\lambda^2} \\
+    \text{Cauchy} & \e^{-|x|/\lambda} & \frac{2\lambda}{1+q_x^2\lambda^2} \\
     \hline
-    \text{Gauss} & e^{-x^2/2\lambda^2} & \sqrt{2\pi}\lambda e^{-q_x^2\lambda^2/2} \\
+    \text{Gauss} & \e^{-x^2/2\lambda^2} & \sqrt{2\pi}\lambda \e^{-q_x^2\lambda^2/2} \\
     \hline
     \text{Triangle} & 1-|x|/\lambda \quad \text{for} \quad |x|<\lambda & \lambda \sinc^2(q_x\lambda/2) \\
     \hline
   \end{array}
 \end{equation}
 
-In addition, a pseudo-Voigt decay function is available, which is a convex combination of the Cauchy and Gauss decay functions.
+In addition, a pseudo-Voigt decay function is available,
+which is a convex combination of the Cauchy and Gauss decay functions.
 \Emph{
-The decay functions are all normalized such that $\int\!dq_x H(q_x)=2\pi$, which is equivalent to $h(0)=1$.
+The decay functions are all normalized such that $\int\!\d q_x H(q_x)=2\pi$,
+which is equivalent to $h(0)=1$.
 }
 \Emph{
-If the one-dimensional lattice is rotated with respect to the x-axis by an angle $\xi$, the corresponding interference function is calculated
+If the one-dimensional lattice is rotated with respect to the x-axis by an angle $\xi$,
+the corresponding interference function is calculated
 by using the correct projection of the in-plane reciprocal vector instead of $q_x$:
 \begin{equation}
   q_\xi = q_x \cos\xi + q_y \sin\xi.
@@ -601,6 +612,7 @@ by using the correct projection of the in-plane reciprocal vector instead of $q_
 %===============================================================================
 \subsection{Two-dimensional lattice} \label{sec:sect:2dlattice}
 %===============================================================================
+
 For a perfect two-dimensional lattice with lattice basis $(\v a,\v b)$, the position
 correlation function is given by:
 \begin{equation}
@@ -618,9 +630,10 @@ decay functions that are defined in the radial variable
 \begin{equation}
   \phi \coloneqq \sqrt{X^2/\lambda_X^2 + Y^2/\lambda_Y^2},
 \end{equation}
-where $(X,Y)$ are the coordinates in an orthonormal coordinate system, where the $X$-axis is rotated
+where $(X,Y)$ are the coordinates in an orthonormal coordinate system,
+where the $X$-axis is rotated
 by an angle $\gamma$ with respect to the first lattice vector $\v a$.
-This ammounts to convoluting the
+This amounts to convoluting the
 previously given interference function with $H(\q)$ or, equivalently, multiplying
 the position correlation function by the inverse Fourier image of $H(\q)$, called the
 \index{Decay function}
@@ -638,28 +651,33 @@ The interference function can then be written as:
     \hline
     \text{Name} & h(\phi) & H(Q) \\
     \hline
-    \text{Cauchy} & e^{-\phi} & \frac{2\pi\lambda_X\lambda_Y}{(1+Q^2)^{3/2}} \\
+    \text{Cauchy} & \e^{-\phi} & \frac{2\pi\lambda_X\lambda_Y}{(1+Q^2)^{3/2}} \\
     \hline
-    \text{Gauss} & e^{-\phi^2/2} & 2\pi\lambda_X\lambda_Y e^{-Q^2/2} \\
+    \text{Gauss} & \e^{-\phi^2/2} & 2\pi\lambda_X\lambda_Y \e^{-Q^2/2} \\
     \hline
   \end{array}
 \end{equation}
 with $Q \coloneqq \sqrt{q_X^2\lambda_X^2 + q_Y^2\lambda_Y^2}$.
 
-In addition, a pseudo-Voigt decay function is available, which is a convex combination of the Cauchy and Gauss decay functions.
+In addition, a pseudo-Voigt decay function is available,
+which is a convex combination of the Cauchy and Gauss decay functions.
 \Emph{
-The decay functions are all normalized such that $\int\!d\q_\plll H(\q)= 4\pi^2$, which is equivalent to $h(0)=1$.
+The decay functions are all normalized such that $\int\!\d^2q_\plll H(\q)= 4\pi^2$,
+which is equivalent to $h(0)=1$.
 }
 
 %===============================================================================
 \subsection{The one-dimensional paracrystal} \label{sec:sect:1dparacrystal}
 %===============================================================================
-A paracrystal, originally developed by Hosemann\cite{Hos51}, models the cumulative disorder of
-the interparticle distances. Although the paracrystal in one dimension is not directly implemented in \BornAgain, it forms the
-basis for the paracrystal models in \BornAgain\ and will thus be discussed first.
+
+A paracrystal, originally developed by Hosemann \cite{Hos51},
+models the cumulative disorder of the interparticle distances.
+Although the paracrystal in one dimension is not directly implemented in \BornAgain,
+it forms the basis for the paracrystal models in \BornAgain\ and will thus be discussed first.
 
 In one dimension, the paracrystal is parameterized by the
-position distribution of the nearest neighbour $p(x)$, centered at a peak distance $D$. The probablility
+position distribution of the nearest neighbour $p(x)$, centered at a peak distance $D$.
+The probablility
 distribution of the position of the nearest neighbour to the right of the particle at the origin then becomes
 \begin{equation}
   p_1(x) = p(x-D).
@@ -685,7 +703,7 @@ Since the Fourier transformed distribution $P(q_x)$ will always be $1$ at $q_x=0
 This singularity is caused by the forward scattering contribution of an infinitely long lattice and can be removed in two ways in \BornAgain.
 The simplest way is to multiply $P(q_x)$ by a fixed factor that is very close, but not equal to, one:
 \begin{equation}
-  P'(q_x) \coloneqq P(q_x) e^{-D/\Lambda},
+  P'(q_x) \coloneqq P(q_x) \e^{-D/\Lambda},
 \end{equation}
 where $\Lambda$ has dimension of length and is called the damping length.
 
@@ -769,7 +787,7 @@ The following distributions are currently supported:
     \hline
     \text{Gate} & \frac{1}{\pi\omega_X\omega_Y}  \quad \text{for} \quad \phi < 1 & \frac{2J_1(Q)}{Q} \\
     \hline
-    \text{Cone} & \frac{3}{\pi\omega_X\omega_Y}\left(1 - \phi \right)  \quad \text{for} \quad \phi < 1 & 6\left[ \frac{J_1(Q)}{Q} - \frac{1}{Q^3}\int_0^Q\!du u^2 J_0(u) \right] \\
+    \text{Cone} & \frac{3}{\pi\omega_X\omega_Y}\left(1 - \phi \right)  \quad \text{for} \quad \phi < 1 & 6\left[ \frac{J_1(Q)}{Q} - \frac{1}{Q^3}\int_0^Q\!\d u u^2 J_0(u) \right] \\
     \hline
   \end{array}
 \end{equation}
@@ -799,7 +817,7 @@ we present particle assembly models that are supported in \BornAgain.
 %\index{Particles!intersecting layer boundaries}}
 
 The nanoparticles are characterized by their form factors
-(\idest the Fourier transform of the shape function - see \cref{app:ff} for a list of form factors implemented in \BornAgain) and the composing material.
+(\cref{SFF}) and the composing material.
 
 
 \Note{\indent
@@ -842,13 +860,16 @@ where
   I_d(\q) &
   \equiv {\bra\left\rvert F_\alpha (\q) \right\rvert^2\ket}_{\alpha}
        - \left\lvert {\bra F_\alpha (\q)\ket}_{\alpha} \right\rvert^2, \\
-  S_{\alpha\beta} (\q) &\equiv 1 + \rho_V \int_V d^3\r\mathcal{G}_{\alpha\beta}(\r) \exp \left[ i\q\cdot \r \right].
+  S_{\alpha\beta} (\q) &\equiv 1 + \rho_V \int_V \d^3\r\mathcal{G}_{\alpha\beta}(\r)
+                       \exp \left[ i\q\cdot \r \right].
 \end{align*}
-$S_{\alpha\beta} (\q)$ is called the \emph{interference function} and $\bra\dotso\ket_\alpha$ is the expectation value over the classes $\lbrace \alpha\rbrace$.
+$S_{\alpha\beta} (\q)$ is called the \emph{interference function},
+and $\bra\dotso\ket_\alpha$ is the expectation value over the classes $\lbrace \alpha\rbrace$.
 
 
 TO MERGE IN:
-More precisely, the probability density for finding a particle $\alpha$ at position $\v{R}_\alpha$ and another one of type $\beta$ at $\v{R}_\beta$ is given by:
+More precisely, the probability density for finding a particle~$\alpha$
+at position $\v{R}_\alpha$ and another one of type $\beta$ at $\v{R}_\beta$ is given by:
 \begin{equation*}
   \mathcal{P}(\alpha, \v{R}_\alpha ; \beta , \v{R}_\beta ) \equiv \rho_V^2 p_\alpha p_\beta \ppcf{\alpha}{\beta}{R}.
 \end{equation*}
@@ -858,7 +879,10 @@ More precisely, the probability density for finding a particle $\alpha$ at posit
 \section{Horizontal particle distributions}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-To proceed further, when the morphology and topology are not exactly known, some hypotheses need to be made since the correlation between the kinds of scatterers and their relative positions included in the pair correlation functions are difficult to estimate. Several options are available:
+To proceed further, when the morphology and topology are not exactly known,
+some hypotheses need to be made since the correlation between the kinds of scatterers
+and their relative positions included in the pair correlation functions are difficult to estimate.
+Several options are available:
 
 
 
@@ -878,10 +902,21 @@ The particles are positioned at regular intervals generating a layout characteri
 This lattice can be two or one-dimensional depending on the characteristics of the particles. For example when they are infinitely long, the implementation can be simplified and reduced to a "pseudo" 1D system.
 
 \index{Paracrystal}
-A paracrystal, whose notion was developed by Hosemann\cite{Hos51}, allows fluctuations of the lengths and orientations of lattice vectors. Paracrystals can be defined as distorted crystals in which the crystalline order has not disappeared and for which the behavior of the interference functions  at small angles is coherent.
+A paracrystal, whose notion was developed by Hosemann \cite{Hos51},
+allows fluctuations of the lengths and orientations of lattice vectors.
+Paracrystals can be defined as distorted crystals
+in which the crystalline order has not disappeared
+and for which the behavior of the interference functions at small angles is coherent.
 It is a transition between the regular lattice and the disordered state.\\
 
-For example, in one dimension, a paracrystal is generated using the following method. First we place a particle at the origin. The second particle is put at a distance $x$ with a density probability $p(x)$ that is peaked at a mean value $D$: $\int_{-\infty} ^{\infty}p(x)dx=1$ and $\int_{-\infty}^{\infty}xp(x)dx=D$. The third one is added at a distance $y$ from the second site using the same rule with a density probability $p_2(y)= \int_{-\infty}^{\infty}p(x)p(y-x)dx=p\otimes p(y)$.\\ With such a method, the pair correlation function $g(x)$ is built step by step. Its expression and the one of its Fourier transform, which is the interference function are
+For example, in one dimension, a paracrystal is generated using the following method. First we place a particle at the origin. The second particle is put at a distance $x$ with a density probability $p(x)$ that is peaked at a mean value $D$: $\int_{-\infty} ^{\infty}p(x)\d x=1$
+and $\int_{-\infty}^{\infty}xp(x)\d x=D$.
+The third one is added at a distance $y$ from the second site
+using the same rule with a density probability
+$p_2(y)= \int_{-\infty}^{\infty}p(x)p(y-x)dx=p\otimes p(y)$.
+
+With such a method, the pair correlation function $g(x)$ is built step by step.
+Its expression and the one of its Fourier transform, which is the interference function are
 \begin{equation*}
 g(x)=\delta(x)+ p(x)+ p(x)\otimes p(x)+\ldots + p(-x)+\ldots \: \mathrm{and}\:\, S(q)=\Re\left(\dfrac{1+P(q)}{1-P(q)}\right),
 \end{equation*}
@@ -908,8 +943,16 @@ where $\Lambda$ is a damping length used in order to introduce some finite-size
 \label{fig:1dparas_q}
 \end{figure}
 
-In two dimensions, the paracrystal is constructed on a pseudo-regular lattice with base vectors $\v{a}$ and $\v{b}$ using the following conditions for the densities of probabilities:\\ $\int p_{\v{a}}(\r)d^2r=\int p_{\v{b}}(\r)d^2r=1$, $\int \r p_{\v{a}}(\r)d^2r=\v{a}$, $\int \r p_{\v{b}}(\r)d^2r=\v{b}$.\\
-In the ideal case the deformations along the two axes are decoupled and each unit cell should retain a parallelogram shape. The interference function is given by\\ $S(q_{\plll})=\prod_{k=a,b}\Re\left(\dfrac{1+P_k(q_{\plll})}{1-P_k(q_{\plll})} \right)$ with $P_k$ the Fourier transform of $p_k$, $k=a, b$.
+In two dimensions, the paracrystal is constructed on a pseudo-regular lattice with base vectors $\v{a}$ and $\v{b}$ using the following conditions for the densities of probabilities:\\
+$\int p_{\v{a}}(\r)\d^2r=\int p_{\v{b}}(\r)\d^2r=1$,
+$\int \r p_{\v{a}}(\r)\d^2r=\v{a}$,
+$\int \r p_{\v{b}}(\r)\d^2r=\v{b}$.
+
+In the ideal case the deformations along the two axes are decoupled
+and each unit cell should retain a parallelogram shape.
+The interference function is given by
+$S(q_{\plll})=\prod_{k=a,b}\Re\left(\dfrac{1+P_k(q_{\plll})}{1-P_k(q_{\plll})} \right)$
+with $P_k$ the Fourier transform of $p_k$, $k=a, b$.
 
 \paragraph{Probability distributions} \mbox{}\\
 The scattering by an ordered lattice gives rise to a series of Bragg peaks situated at the nodes of the reciprocal lattice. Any divergence from the ideal crystalline case modifies the output spectrum by, for example, widening or attenuating the Bragg peaks. The influence of these "defects" can be accounted for
@@ -993,18 +1036,20 @@ I(\q ) = {\bra\left| F_\alpha(\q ) \right| ^2\ket}_{\alpha}
 \end{equation}
 with
 \begin{align*}
-  \tilde{p}_\kappa(\q ) &= \int d\alpha\; p(\alpha) e^{i\kappa q_x \Delta R(\alpha)}  \\
-  \Omega_\kappa(\q ) &= \tilde{p}_{2\kappa}(\q ) \phi(\q) e^{i q_x D_0}  \\
-  \widetilde{\curlf_\kappa}(\q ) &= \int d\alpha\; p(\alpha)F_\alpha (\q ) e^{i\kappa q_x \Delta R(\alpha)},
+  \tilde{p}_\kappa(\q ) &= \int \d\alpha\; p(\alpha) \e^{i\kappa q_x \Delta R(\alpha)}  \\
+  \Omega_\kappa(\q ) &= \tilde{p}_{2\kappa}(\q ) \phi(\q) \e^{i q_x D_0}  \\
+  \widetilde{\curlf_\kappa}(\q ) &=
+       \int \d\alpha\; p(\alpha)F_\alpha (\q ) \e^{i\kappa q_x \Delta R(\alpha)},
 \end{align*}
 and the Fourier transform of $P_1(x|\alpha_0\alpha_1)$ is
 \begin{equation*}
-  \curlp (\q ) = \phi (\q )e^{i q_x D_0} e^{i \kappa q_x \left[ \Delta R(\alpha_0) + \Delta R(\alpha_1) \right] }.
+  \curlp (\q ) = \phi (\q )\e^{i q_x D_0}
+                    \e^{i \kappa q_x \left[ \Delta R(\alpha_0) + \Delta R(\alpha_1) \right] }.
 \end{equation*}
 
 Using the result from the one-dimensional analysis, one can apply this formula ad hoc for distributions of particles in a plane, where the coordinate $x$ will now be replaced with $(x,y)$, while the $s$ coordinate encodes a position in the remaining orthogonal direction. One must be aware however that this constitutes a further approximation, since this type of correlation does not have a general solution in more than one dimension.
 
-The intensity in \cref{Esscainf} will contain a Dirac delta function contribution, caused by taking an infinite sum of terms that are perfectly correlated at $\q = 0$. One can leverage this behaviour by multiplying the nearest neighbor distribution by a constant factor $e^{-D/\Lambda}$, which removes the division by zero in \cref{Esscainf}.
+The intensity in \cref{Esscainf} will contain a Dirac delta function contribution, caused by taking an infinite sum of terms that are perfectly correlated at $\q = 0$. One can leverage this behaviour by multiplying the nearest neighbor distribution by a constant factor $\e^{-D/\Lambda}$, which removes the division by zero in \cref{Esscainf}.
 Another way of dealing with this infinity at $\q =0$ consists of taking only a finite number of terms, in which case the geometric series still has an analytic solution, but becomes a bit more cumbersome:
 \begin{equation*}
 \begin{split}
@@ -1119,10 +1164,10 @@ The system considered in this section consists of particles encapsulated in a la
 
 \begin{align}
   F_{\rm{DWBA}}(q_{\plll}, k_{i,z}, k_{f,z})
-  &= T_\text{i} T_\text{f} F_{\rm{BA}}(q_{\plll}, k_{i,z}-k_{f,z})e^{i(k_{i,z}-k_{f,z})d}\nonumber \\
-  &+ R_\text{i} T_\text{f} F_{\rm{BA}}(q_{\plll}, -k_{i,z}-k_{f,z})e^{i(-k_{i,z}-k_{f,z})d} \nonumber \\
-  &+ R_\text{f} T_\text{i} F_{\rm{BA}}(q_{\plll}, k_{i,z}+k_{f,z}) e^{i(k_{i,z}+k_{f,z})d}\nonumber \\
-  &+ R_\text{f} R_\text{i}F_{\rm{BA}}(q_{\plll},-k_{i,z}+k_{f,z})e^{i(-k_{i,z}+k_{f,z})d}, \label{Edwbaburied}
+  &= T_\text{i} T_\text{f} F_{\rm{BA}}(q_{\plll}, k_{i,z}-k_{f,z})\e^{i(k_{i,z}-k_{f,z})d}\nonumber \\
+  &+ R_\text{i} T_\text{f} F_{\rm{BA}}(q_{\plll}, -k_{i,z}-k_{f,z})\e^{i(-k_{i,z}-k_{f,z})d} \nonumber \\
+  &+ R_\text{f} T_\text{i} F_{\rm{BA}}(q_{\plll}, k_{i,z}+k_{f,z}) \e^{i(k_{i,z}+k_{f,z})d}\nonumber \\
+  &+ R_\text{f} R_\text{i}F_{\rm{BA}}(q_{\plll},-k_{i,z}+k_{f,z})\e^{i(-k_{i,z}+k_{f,z})d}, \label{Edwbaburied}
 \end{align}
 
 \begin{equation*}
@@ -1454,7 +1499,7 @@ where ($L_1$, $L_2$, $\alpha$, $\xi$) are shown in \cref{fig:2dlattice} with
 \begin{itemize}
 \item[]$L_1$, $L_2$ the lengths of the lattice cell,
 \item[]$\alpha$ the angle between the lattice basis vectors $\v{a}, \v{b}$ in direct space,
-\item[] $\xi$ is the angle defining the lattice orientation (set to $0$ by default); it is taken as the angle between the $\v{a}$ vector of the lattice basis and the $\v{x}$ axis of the reference Cartesian frame (as shown in \cref{fig:multil3d}).
+\item[] $\xi$ is the angle defining the lattice orientation (set to $0$ by default); it is taken as the angle between the $\v{a}$ vector of the lattice basis and the $\v{x}$ axis of the reference Cartesian frame (as shown in \cref{fig:2dlattice}).
 \end{itemize}
 
 \begin{figure}[tb]
@@ -1591,5 +1636,5 @@ Function  & Parameters & Comments\\
 \caption{List of interference functions implemented in \BornAgain. pdf : probability distribution function, $\v{a}, \v{b}$ are the lattice base vectors, and $\v{x}$ is the axis vector perpendicular to the detector plane.}
 \end{table}
 
-\index{Particle assemblies)}%
+\index{Particle assemblies|)}%
 %\fi
diff --git a/Doc/UserManual/BornAgainManual.tex b/Doc/UserManual/BornAgainManual.tex
index 14b73ad9aa9eaa398941add8c8eaa8bcb4312986..99412df92ebba5d7bf813f51d56f850e42295d38 100644
--- a/Doc/UserManual/BornAgainManual.tex
+++ b/Doc/UserManual/BornAgainManual.tex
@@ -44,21 +44,17 @@ Walter Van Herck, Joachim Wuttke}
 \tableofcontents\cleardoublepage
 
 \include{Introduction}
-
 \include{SAS}
 \include{GisasFoundations}
-\include{PolarizedScattering}
 \include{Multilayers}
-\include{polMultilayers}
 \include{Assemblies}
-%void \include{Domains}
-%void \include{Roughness}
-%void \include{Experiment}
+%%\include{Roughness}
+%%\include{PolarizedScattering}
 
 \include{ThreeAPIs}
 \include{Usage}
-%\include{Fitting}
-%\include{IntensityData}
+%%\include{Fitting}
+%%\include{IntensityData}
 
 \include{FormFactors}
 
diff --git a/Doc/UserManual/Domains.tex b/Doc/UserManual/Domains.tex
deleted file mode 100644
index d4d783953af825bdb33b131651775d6a0b4e7c6e..0000000000000000000000000000000000000000
--- a/Doc/UserManual/Domains.tex
+++ /dev/null
@@ -1,19 +0,0 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%
-%%   BornAgain User Manual
-%%
-%%   homepage:   http://www.bornagainproject.org
-%%
-%%   copyright:  Forschungszentrum Jülich GmbH 2015
-%%
-%%   license:    Creative Commons CC-BY-SA
-%%   
-%%   authors:    Scientific Computing Group at MLZ Garching
-%%               C. Durniak, M. Ganeva, G. Pospelov, W. Van Herck, J. Wuttke
-%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-
-\chapter{Scattering by domain structures}  \label{sec:Domains}
-
-\MissingSection
diff --git a/Doc/UserManual/Experiment.tex b/Doc/UserManual/Experiment.tex
deleted file mode 100644
index 5df921e5784aec8aa3284c6ba825164e499826e1..0000000000000000000000000000000000000000
--- a/Doc/UserManual/Experiment.tex
+++ /dev/null
@@ -1,36 +0,0 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%
-%%   BornAgain User Manual
-%%
-%%   homepage:   http://www.bornagainproject.org
-%%
-%%   copyright:  Forschungszentrum Jülich GmbH 2015
-%%
-%%   license:    Creative Commons CC-BY-SA
-%%   
-%%   authors:    Scientific Computing Group at MLZ Garching
-%%               C. Durniak, M. Ganeva, G. Pospelov, W. Van Herck, J. Wuttke
-%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\chapter{Simulating a GISAS experiment}  \label{sec:Exp}
-
-To fit experimental data, it is not sufficient to simulate
-the scattering of plane waves by the sample.
-To account for actual detector images,
-we must account for the finite instrument resolution,
-and compute detector images in appropriate coordinates.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{Instrumental resolution}\label{Sresolution}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\MissingSection
-
-%In this section we describe the conversion between the scattering vector~$\q$
-%and the experimental coordinates.
-
-%The experimental geometry has already been introduced in \cref{fig:expgeom}.
-
-\MissingSection
-
diff --git a/Doc/UserManual/Introduction.tex b/Doc/UserManual/Introduction.tex
index b15a9c427fac123a5c1312d9894bd604d6f1c397..677069713a9c0cfe5c33172a773b7d21bd167450 100644
--- a/Doc/UserManual/Introduction.tex
+++ b/Doc/UserManual/Introduction.tex
@@ -209,10 +209,6 @@ with the software or the documentation.}}
   an \textbf{implementation note} that explains
   how the theory exposed in this manual is actually used in \BornAgain.}}
 
-\medskip
-\demobox{\Tuto{\indent This is a link to a usage example in the online tutorial.
-  To be merged into this manual \ldots}}
-
 \medskip
 \demobox{\Link{\indent This is a link to the online docs.}}
 
diff --git a/Doc/UserManual/Macros.tex b/Doc/UserManual/Macros.tex
index 4e520068feb148ee2dc95d1a1979caf1285634c6..6e26b3a869cafdc9cbf13ab86405fd27b07bc48a 100644
--- a/Doc/UserManual/Macros.tex
+++ b/Doc/UserManual/Macros.tex
@@ -81,4 +81,4 @@
 %	HYPHENATION
 %-------------------------------------------------------------------------------
 
-\hyphenation{ MacOS Schrö-ding-er nano-par-ti-cle nano-par-ti-cles }
+\hyphenation{ MacOS Schrö-ding-er nano-par-ti-cle nano-par-ti-cles para-crys-tal wave-num-ber}
diff --git a/Doc/UserManual/Multilayers.tex b/Doc/UserManual/Multilayers.tex
index 2c87a81686ebc94901f7f619e882f562d524f030..2f6d9c56304e66bd9976f7a4a7a074b9c49e6488 100644
--- a/Doc/UserManual/Multilayers.tex
+++ b/Doc/UserManual/Multilayers.tex
@@ -441,3 +441,5 @@ The following cases are treated seperately:
     representing a slight absorption. However, this should be inconsequential because the index of refraction
     of non-vacuum layer always contains an absorptive component.
 \end{itemize}
+
+\index{Multilayer|)}%
diff --git a/Doc/UserManual/PolarizedScattering.tex b/Doc/UserManual/PolarizedScattering.tex
index b67f20c8ef2dfd282f71baa3d114f708f03450cc..fe40ce8f81dea948777e4cfdac46267a48e55bcb 100644
--- a/Doc/UserManual/PolarizedScattering.tex
+++ b/Doc/UserManual/PolarizedScattering.tex
@@ -7,7 +7,7 @@
 %%   copyright:  Forschungszentrum Jülich GmbH 2015
 %%
 %%   license:    Creative Commons CC-BY-SA
-%%   
+%%
 %%   authors:    Scientific Computing Group at MLZ Garching
 %%               C. Durniak, M. Ganeva, G. Pospelov, W. Van Herck, J. Wuttke
 %%
@@ -40,3 +40,154 @@ in contrast to the scalar theory of the previous chapter.
 
 \MissingSection
 
+
+%===============================================================================
+\subsection{Wave equation and propagation within one layer}
+%===============================================================================
+
+%\cite{Deak_ppt, PhysRevB.76.224420, Deak2001113, PhysRevB.53.6158}.
+%\cite{RevModPhys.23.287}
+
+To allow for polarization-dependent interactions,
+we replace the squared index of refraction $n^2$
+by $1+\uu\chi$, where $\uu\chi$ is a $2\times 2$ susceptibility matrix.
+The wave equation \cref{Escalar_wave} for layer~$\il$ becomes
+\begin{equation}\label{Ewaveqp}
+(\Delta +K^2 +K^2 \uu\chi_\il) \u\psi(\r)= 0,
+\end{equation}
+where $\u\psi(\r)$ is a two-component spinor wavefunction,
+with components $\psi_\UP(\r)$ and~$\psi_\DN(\r)$.
+At interfaces between layers,
+both spinor components of $\u\psi(\r)$ and $\Nabla\u\psi(\r)$
+must evolve continuously.
+
+The reasons for the factorization \cref{Ewave3} still apply,
+and so we can write
+\begin{equation}\label{Ewave3p}
+\u\psi(\r) = \u\psi(z) \e^{i \k_\parallel\r_\parallel}.
+\end{equation}
+As before, $\k_\parallel$ is constant across layers.
+The wave equation~\cref{Ewaveqp} reduces to
+\begin{equation}\label{Ewavezp}
+\left(\partial_z^2 + K^2 + K^2\uu\chi_\il - k_\parallel^2 \right) \u\psi(z) = 0.
+\end{equation}
+We abbreviate
+\begin{equation}
+  \uu H_\il \coloneqq  K^2(1+\uu\chi_\il)-k_\parallel^2
+\end{equation}
+so that the wave equation becomes simply
+\begin{equation}\label{Ewaveqp2}
+  \left(\partial_z^2 + \uu H_\il\right) \u\psi(z) = 0.
+\end{equation}
+The solution is
+\begin{equation}\label{Epsizp}
+  \u\psi_\il(z)
+  = \sum_{k=1}^2 \u x_{\il k}\left(\alpha_{\il k}\e^{i p_{\il k}(z-z_k)}
+                            + \beta_{\il k}\e^{-i p_{\il k}(z-z_k)}\right),
+\end{equation}
+where the $\u x_{\il k}$ are eigenvectors of $\uu H_\il$
+with eigenvalues $p_{\il k}^2$:
+\begin{equation}
+  \left( -p_{\il k}^2 + \uu H_\il \right) \u x_{\il k} = 0
+   \;\text{ for }\;\il=1,2.
+\end{equation}
+In a reproducible algorithm,
+the eigenvectors $\u x_{\il k}$ must be chosen according to some arbitrary
+normalization rule,
+for instance
+\begin{equation}
+  |\u x_{\il k}|=1,\quad x_{i\il\UP} \text{ real and nonnegative}.
+\end{equation}
+Similarly,
+a rule is needed how to handle the case of one degenerate eigenvalue,
+which includes in particular the case of scalar interactions.
+
+
+%-------------------------------------------------------------------------------
+\subsection{Wave propagation across layers}
+%-------------------------------------------------------------------------------
+
+Generalizing \cref{Evecc},
+we introduce the coefficient vector
+\begin{equation}
+  c_\il \coloneqq  {(\alpha_{\il1}, \alpha_{\il2}, \beta_{\il1}, \beta_{\il2})}^\text{T}.
+\end{equation}
+To match solutions for neighboring layers,
+continuity is requested for both spinorial components
+of $\u\psi$ and $\Nabla\u\psi$.
+As before \cref{EFcFDc}, we have at the bottom of layer~$\il$
+\begin{equation}\label{EFcFDcp}
+  F_\il c_\il = F_{\il+1} D_{\il+1} c_{\il+1},
+\end{equation}
+where the matrices are now
+\begin{equation}
+  F_\il \coloneqq  \left(\begin{array}{cccc}
+    x_{i1\UP}      &x_{i2\UP}     &x_{i1\UP}       &x_{i2\UP}       \\
+    x_{i1\DN}      &x_{i2\DN}     &x_{i1\DN}       &x_{i2\DN}       \\
+    x_{i1\UP}p_{\il1}&x_{i2\UP}p_{\il2}&-x_{i1\UP}p_{\il1}&-x_{i2\UP}p_{\il2}\\
+    x_{i1\DN}p_{\il1}&x_{i2\DN}p_{\il2}&-x_{i1\DN}p_{\il1}&-x_{i2\DN}p_{\il2}
+  \end{array}\right)
+\end{equation}
+and
+\begin{equation}
+  D_\il \coloneqq  \text{diag}(\delta_{\il1}, \delta_{\il2}, \delta_{\il1}^*, \delta_{\il2}^*)
+\end{equation}
+with the phase factor
+\begin{equation}
+   \delta_{\il k} \coloneqq  \e^{ip_{\il k}d_k}.
+\end{equation}
+Note that matrix $F_\il$ has the block form
+\begin{equation}
+  F_\il
+  =\left(\begin{array}{ll}\uu x_\il&\hphantom{-}\uu x_\il\\[1ex]
+    \uu x_\il\; \uu P_\il&-\uu x_\il\; \uu P_\il\end{array}\right)
+    = \uu x_\il \cdot
+    \left(\begin{array}{cc}\uu 1&\uu 1\\[1ex]
+    \uu P_\il&-\uu P_\il\end{array}\right),
+\end{equation}
+with
+\begin{equation}
+  \uu x_\il \coloneqq
+  \left(\u x_{\il1}, \u x_{\il2}\right),
+  \quad
+  \uu P_\il \coloneqq
+  \text{diag}\left(p_{\il1},p_{\il2}\right).
+\end{equation}
+This facilitates the computation of the inverse
+\begin{equation}
+  F_\il^{-1}
+    = \frac{1}{2}
+    \left(\begin{array}{cc}\uu 1&\hphantom{-}\uu P_\il^{-1}\\[1.2ex]
+      \uu 1 &-\uu P_\il^{-1}\end{array}\right)
+      \cdot\uu x_\il^{-1},
+\end{equation}
+which is needed for the transfer matrix $M_\il$,
+defined as in \cref{Edef_M}.
+With the new meaning of $c_\il$ and $M_\il$,
+the recursion \cref{EcMc} and the explicit solution~\cref{Eci}
+hold as derived above.
+To resolve~\cref{Eci} for the reflected amplitudes $\alpha_{0\il}$
+as function of the incident amplitudes $\beta_{0\il}$,
+we choose the notations
+\begin{equation}
+  \u\alpha_\il
+  \coloneqq \left(\begin{array}{c}\alpha_{\il1}\\\alpha_{\il2}\end{array}\right),\quad
+  \u\beta_\il
+  \coloneqq \left(\begin{array}{c}\beta_{\il1}\\\beta_{\il2}\end{array}\right),\quad
+  M\coloneqq M_1 ... M_N % TODO restore \cdots
+  \eqqcolon \left(\begin{array}{cc}\uu m_{11}&\uu m_{12}\\
+                           \uu m_{21}&\uu m_{22}\end{array}\right),
+\end{equation}
+where the $\uu m_{\il k}$ are $2\times2$ matrices.
+Eq.~\cref{Eci} then takes the form
+\begin{equation}
+  \left(\begin{array}{c}\u\alpha_{0}\\\u\beta_{0}\end{array}\right)
+  =
+  \left(\begin{array}{cc}\uu m_{11}&\uu m_{12}\\
+    \uu m_{21}&\uu m_{22}\end{array}\right)
+  \left(\begin{array}{c}\u{0}\\\u\beta_{N}\end{array}\right),
+\end{equation}
+which immediately yields
+\begin{equation}
+  \u\alpha_0 = \uu m_{12}\,\uu m_{22}^{-1}\,\u\beta_0.
+\end{equation}
diff --git a/Doc/UserManual/SAS.tex b/Doc/UserManual/SAS.tex
index ee70a48f6e255aff81e433b87d16ae723ae175fb..e3cb67496a0186090b22e04a19a69084856438f7 100644
--- a/Doc/UserManual/SAS.tex
+++ b/Doc/UserManual/SAS.tex
@@ -36,7 +36,7 @@ needed for grazing-incidence small-angle scattering (\Cref{SBornApprox}).
 
 Finally, we discuss how the physical scattering law
 (the double differential cross section as function of the scattering vector~$\q$)
-relates to the experimental detector image (Sec.\ref{SdetImg}).
+relates to the experimental detector image (\Cref{SdetImg}).
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -45,6 +45,8 @@ relates to the experimental detector image (Sec.\ref{SdetImg}).
 \index{Wave propagation!neutrons|(}%
 \index{Neutrons!wave propagation|(}%
 
+\def\Vmac{\tilde{V}}
+
 \index{Schrodinger@Schrödinger equation!microscopic}%
 The scalar wavefunction $\psi(\r,t)$
 \nomenclature[2t020]{$t$}{Time}%
@@ -137,58 +139,75 @@ The macroscopic field equation
 has still the form of a stationary Schrödinger equation,
 \index{Schrodinger@Schrödinger equation!macroscopic}%
 \begin{equation}\label{EmacrSchrodi}
-  \left\{-\frac{\hbar^2}{2m}\Nabla^2+v(\r)-\hbar\omega\right\} \psi(\r) = 0,
+  \left\{-\frac{\hbar^2}{2m}\Nabla^2+\Vmac(\r)-\hbar\omega\right\} \psi(\r) = 0,
 \end{equation}
 \nomenclature[1ψ030 2r040 2t020]{$\psi(\r,t)$}{Coherent wavefunction}%
-\nomenclature[2v020 2r040]{$v(\r)$}{Macrosopic optical potential}%
+\nomenclature[2v131 2r040]{$\Vmac(\r)$}{Macroscopic optical potential}%
 where $\psi$ now stands for the \E{coherent wavefunction}
 \index{Coherent wavefunction}%
 \index{Wave propagation!coherent}%
 obtained by superposition of
 incident and forward scattered states,
-and $v(\r)$ is the \E{macroscopic optical potential}.
+and $\Vmac(\r)$ is the \E{macroscopic optical potential}.
 \index{Optical potential!macroscopic}%
 This potential is weak, and slowly varying compared to atomic length scales.
-It can be rewritten in a number of ways,
-especially in terms of a
+In the following it shall be expressed through the
 \E{bound scattering length density}
 \index{Scattering length density}%
-$\rho_s(\r)$ \cite[eq.\ 2.8.37]{Sea89},
-\nomenclature[1ρ 034 2s000 2r040]{$\rho_s(\r)$}{Scattering length density}%
+\cite[eq.\ 2.8.37]{Sea89},
+\nomenclature[2v020 2r040]{$v(\r)$}{Scattering length density}%
 \begin{equation}
-  v(\r)=\frac{2\pi \hbar^2}{m}\rho_s(\r),
-\end{equation}
-or of a \E{refractive index}~$n(\r)$
-\nomenclature[2n020 2r040]{$n(\r)$}{Refractive index}%
-\index{Refractive index} % !vs scattering length density}%
+  v(\r)\coloneqq\frac{m}{2\pi \hbar^2}\Vmac(\r),
+\end{equation}
+which shall be decomposed into a constant average and a fluctuating part,
+\begin{equation}\label{Edecompose_v}
+  v(\r) = \overline{v} + \delta v(\r).
+\end{equation}
+\nomenclature[1d030 2v230 2r040]{$\delta v(\r)$}{Fluctuating part of the scattering length density}%
+\nomenclature[2v021]{$\overline{v}$}{Average scattering length density}%
+The macroscopic Schrödinger equation~\cref{EmacrSchrodi} becomes
+\begin{equation}\label{ESchrodi3}
+  \left\{ \Nabla^2 + \frac{2m\omega}{\hbar} - 4\pi\overline{v} \right\}\psi(\r)
+  = 4\pi\delta v(\r)\psi(\r).
+\end{equation}
+In a homogeneous medium with $\delta v=0$,
+\cref{ESchrodi3} reduces to the Helmholtz wave equation.
+This motivates us to rewrite~\cref{ESchrodi3} as
+\Emph{
+\begin{equation}\label{ESchrodiK}
+  \left\{ \Nabla^2 + K^2 \right\}\psi(\r)
+  = 4\pi\delta v(\r)\psi(\r),
+\end{equation}}
+where the wavenumber~$K$ is given by the dispersion relation
+\nomenclature[2k120]{$K$}{Wavenumber}%
+\index{Dispersion!neutron in homogeneous medium}%
+\index{Wavenumber!neutron}%
+\begin{equation}
+  K^2 = \frac{2m\omega}{\hbar} - 4\pi\overline{v}
+      = K_0^2\left(1-\frac{4\pi}{K_0^2}\overline{v}\right)
+      = K_0^2 n^2.
+\end{equation}
+To arrive at the final compact expression $K=K_0 n$
+we introduced the vacuum wavenumber~$K_0$
+\nomenclature[2k122]{$K_0$}{Wavenumber in vacuum}%
+and the \E{refractive index}
+\nomenclature[2n020]{$n$}{Refractive index}%
+\index{Refractive index}% !vs scattering length density}%
 \index{Index of refraction|see {Refractive index}}%
-defined by
 \begin{equation}\label{EnRefrIndx}
-  n(\r)^2\coloneqq 1-\frac{4\pi}{K^2}\rho_s(\r) = 1 -\frac{2m}{\hbar^2 K^2}v(\r).
-\end{equation}
-In the latter expression,
-we introduced the \E{vacuum wavenumber}~$K$,
-\nomenclature[2k120]{$K$}{Vacuum wavenumber, corresponding to the frequency~$\omega$}%
-which is connected with the frequency~$\omega$ through the
-\E{dispersion relation}
-\begin{equation}
-  \frac{\hbar^2 K^2}{2m} = \hbar\omega.
-\end{equation}
-Since we only consider stationary solutions~\cref{Estationarywave},
-$\omega$ will not appear any further in our derivations.
-Instead, we use~$K$ as the given parameter that characterizes the
-incoming radiation.
-In terms of $K$ and $n$,
-the macroscopic Schrödinger equation \cref{EmacrSchrodi}
-can be rewritten as
-\Emph{
-\begin{equation}\label{EnSchrodi}
-  \left\{\Nabla^2+K^2n(\r)^2\right\}\psi(\r) = 0.
-\end{equation}\vspace*{-10pt}
-}
-This equation is the starting point for the analysis of all
-small-angle scattering experiments,
-whether under grazing incidence (GISAS) or not (regular SAS).
+  n\coloneqq \sqrt{1-\frac{4\pi}{K_0^2}\overline{v}}.
+\end{equation}
+
+The Schrödinger equation \cref{ESchrodiK}
+is the starting point for the analysis of small-angle scattering (SAS) experiments
+in Born approximation (\cref{SBornApprox}).
+Later (\cref{Swave21}), for the special case of grazing incidence,
+the distorted wave Born approximation (DWBA)
+\index{Distorted-wave Born approximation}%
+will start with a slight modification of the decomposition~\cref{Edecompose_v}:
+vertical variations of the average scattering length density
+will be reassigned from $\delta v(\r)$ to $n(z)$
+in order to account for refraction and reflection in layered samples.
 \index{SAS|see {Small-angle scattering}}%
 \index{Small-angle scattering}%
 
@@ -206,7 +225,7 @@ whether under grazing incidence (GISAS) or not (regular SAS).
 \index{Born approximation|(}%
 
 To describe an elastic scattering experiment,
-we need to solve the Schrödinger equation~\cref{EnSchrodi}
+we need to solve the Schrödinger equation~\cref{ESchrodiK}
 under the asymptotic boundary condition
 \begin{equation}\label{Escabouco}
   \psi(\r)
@@ -222,49 +241,30 @@ the outgoing scattered wave
 that carries information in form of the angular distribution
 $f(\vartheta,\varphi)$.
 
-For thermal or cold neutrons,
-as for X-rays, the refractive index~$n$ is almost always
-very close to~1.
+For thermal or cold neutrons, as for X-rays,
+the refractive index~$n$ is almost always very close to~1,
+fluctuations of the scattering length density are at most of the order of $1-n$,
+and therefore the right-hand side of~\cref{ESchrodiK}
+is but a weak perturbation.
 This suggests a solution of the Schrödinger equation
-by means of a perturbation expansion in powers of $n^2-1$.
-This expansion is named after Max Born
-who introduced it in quantum mechanics.\footnote
-{It goes back to Lord Rayleigh
+by means of a perturbation expansion in powers of $\delta v$,
+known as the \E{Born expansion} or \E{Born series}.\footnote{
+Named after Max Born
+who introduced it in quantum mechanics.
+The idea goes back to Lord Rayleigh
 who devised it for sound,
 and later also applied it to electromagnetic waves,
 which resulted in his famous explanation of the blue sky.}
 
-To carry out this idea, we rewrite the Schrödinger equation
-once more so that it takes the form of a Helmholtz equation
-\index{Helmholtz equation}%
-with a perturbation term on the right side:
-\begin{equation}\label{ESchrodiHelmholtz}
-  \left(\Nabla^2+K^2\right)\psi(\r)
-  = 4\pi\chi(\r)\psi(\r)
-\end{equation}
-with
-\begin{equation}\label{EChiDef}
-  \chi(\r) \coloneqq  \frac{K^2}{4\pi}\left(1-n^2(\r)\right).
-\end{equation}
-\nomenclature[1χ030 2r040]{$\chi(\r)$}{Perturbative potential,
-for neutrons equal to the scattering-length density~$\rho_s$}%
-This definition just compensates \cref{EnRefrIndx} so that $\chi=\rho_s$.
-In the following, we prefer the notation~$\chi$
-and the appellation \E{perturbative potential}
-\index{Potential|see {Perturbation}}%
-\index{Perturbation}%
-over the scattering length density~$\rho_s$
-to prepare for the generalization to the electromagnetic case.
-
-Equation~\cref{ESchrodiHelmholtz} looks
+Equation~\cref{ESchrodiK} looks
 like an inhomogeneous differential equation ---
-provided we neglect for a moment that the unknown function~$\psi$
-reappears on the right side.
-The homogeneous equation
+except that the right-hand side contains the unknown function~$\psi$.
+In the Born expansion this problem is solved by iteration.
+We start with the homogeneous equation
 \begin{equation}\label{EHelmholtzHomog}
-  \left(\Nabla^2+K^2\right)\psi(\r) = 0
+  \left(\Nabla^2+K^2\right)\psi(\r) = 0,
 \end{equation}
-is solved by plane waves and superpositions thereof.
+which is solved by plane waves and superpositions thereof.
 It applies in particular to the incident wave~$\psi_\ti$.
 
 For an isolated inhomogeneity,
@@ -284,40 +284,40 @@ This explains the appearance of the factor $4\pi$.}
   G(\r,\r') = \frac{\e^{iK|\r-\r'|}}{4\pi |\r-\r'|},
 \end{equation}
 which is an outgoing spherical wave centered at $\r'$.
-Convoluting this function with the given inhomogeneity $4\pi\chi\psi$,
+Convoluting this function with the given inhomogeneity $4\pi\delta v\psi$,
 we obtain what is known as the Lippmann-Schwinger equation,
 \index{Lippmann-Schwinger equation}%the formal solution
 \begin{equation}\label{EPsiFormal}
   \psi(\r)
   = \psi_\ti(\r)
-  + \int\!\d^3r'\, G(\r,\r') 4\pi\chi(\r')\psi(\r').
+  + \int\!\d^3r'\, G(\r,\r') 4\pi\delta v(\r')\psi(\r').
 \end{equation}
 This integral equation for $\psi(\r)$ improves
-upon the original stationary Schrödinger equation \cref{ESchrodiHelmholtz}
+upon the original stationary Schrödinger equation \cref{ESchrodiK}
 in that it ensures the boundary condition~\cref{Escabouco}.
 It can be resolved into an infinite series
 by iteratively substituting the full right-hand side of~\cref{EPsiFormal}
 into the integrand.
-Successive terms in this series contain rising powers of $\chi$.
-Since $\chi$ is assumed to be small, the series is likely to converge.
+Successive terms in this series contain rising powers of $\delta v$.
+Since $\delta v$ is assumed to be small, the series is likely to converge.
 In \E{first-order Born approximation},
-only the linear order in $\chi$ is retained,
+only the linear order in $\delta v$ is retained,
 \begin{equation}\label{EBorn}
   \psi(\r)
   \doteq \psi_\ti(\r)
-  + 4\pi \int\!\d^3r'\, G(\r,\r') \chi(\r') \psi_\ti(\r').
+  + 4\pi \int\!\d^3r'\, G(\r,\r') \delta v(\r') \psi_\ti(\r').
 \end{equation}
 This is practically always adequate for
 material investigations with X-rays or neutrons,
 where the aim is to
-deduce $\chi(\r')$ from the scattered intensity ${|\psi(\r)|}^2$.
+deduce $\delta v(\r')$ from the scattered intensity ${|\psi(\r)|}^2$.
 Since detectors are always placed at positions $\r$
 that are not illuminated by the incident beam,
 we are only interested in the scattered wave field
 \begin{equation}\label{EBornS}
   \psi_\text{s}(\r)
   \coloneqq
-  4\pi \int\!\d^3r'\, G(\r,\r') \chi(\r') \psi_\ti(\r').
+  4\pi \int\!\d^3r'\, G(\r,\r') \delta v(\r') \psi_\ti(\r').
 \end{equation}
 \nomenclature[1ψ034 2s000 0 2r040]{$\psi_\text{s}(\r)$}{Scattered wavefunction}%
 \nomenclature[2s000 0]{s}{Subscript ``scattered''}%
@@ -349,8 +349,8 @@ This allows us to expand
   \equiv r - \frac{\k_\tf \r'}{K},
 \end{equation}
 \nomenclature[2f000]{f}{Subscript ``final'',
-for outgoing waves scattered into the direction of the detector}%
-\nomenclature[2k040]{$\k$}{wavevector}
+   for outgoing waves scattered into the direction of the detector}%
+\nomenclature[2k040]{$\k$}{wavevector}%
 where we have introduced the outgoing wavevector
 \begin{equation}
   \k_\tf\coloneqq K\frac{\r}{r}.
@@ -362,14 +362,14 @@ and obtain in leading order the far-field Green function
   G_\text{far}(\r,\r')
   = \frac{\e^{iKr}}{4\pi r}\psi^*_\tf(\r')
 \end{equation}
-\nomenclature[2g134 2far]{$G_\text{far}(\r,\r')$}
-  {Far-field approximation to the Green function $G(\r,\r')$}
+\nomenclature[2g134 2far]{$G_\text{far}(\r,\r')$}{Far-field
+   approximation to the Green function $G(\r,\r')$}%
 where
 \begin{equation}
   \psi_\tf(\r) \coloneqq  \e^{i\k_\tf \r}
 \end{equation}
-\nomenclature[1ψ034 2f000 2r040]{$\psi_\tf(\r)$}
-  {Plane wave propagating from the sample towards the detector}%
+\nomenclature[1ψ034 2f000 2r040]{$\psi_\tf(\r)$}{Plane
+  wave propagating from the sample towards the detector}%
 is a plane wave propagating towards the detector,
 and $\psi^*$ designates the complex conjugate of $\psi$.
 With respect to $\r$, $G_\text{far}$ is an outgoing spherical wave.
@@ -379,21 +379,21 @@ becomes in the far-field approximation
 \begin{equation}\label{EsandwichC}
   \psi_\text{s,far}(\r)
   = \frac{\e^{iKr}}{r}
-    \bra \psi_\tf|\chi|\psi_\ti\ket,
+    \bra \psi_\tf|\delta v|\psi_\ti\ket,
 \end{equation}
-\nomenclature[1ψ034 2s000 2far]{$\psi_\text{s,far}(\r)$}
-  {Far-field approximation to the scattered wavefunction $\psi_\text{s}(\r)$}%
+\nomenclature[1ψ034 2s000 2far]{$\psi_\text{s,far}(\r)$}{Far-field
+   approximation to the scattered wavefunction $\psi_\text{s}(\r)$}%
 where we used Dirac notation for the transition matrix element
 \index{Transition matrix}%
 \begin{equation}\label{Etrama}
-  \bra \psi_\tf|\chi|\psi_\ti\ket
-  \coloneqq  \int\!\d^3r\, \psi^*_\tf(\r)\chi(\r)\psi_\ti(\r).
+  \bra \psi_\tf|\delta v|\psi_\ti\ket
+  \coloneqq  \int\!\d^3r\, \psi^*_\tf(\r)\delta v(\r)\psi_\ti(\r).
 \end{equation}
-\nomenclature[0$\langle$0]{{$\bra\ldots\vert\ldots\vert\ldots\ket$}}
-  {Matrix element, defined as a volume integral}%
+\nomenclature[0$\langle$0]{{$\bra\ldots\vert\ldots\vert\ldots\ket$}}{Matrix
+  element, defined as a volume integral}%
 In order to reconcile conflicting sign conventions,
 we will in the following rather use its complex conjugate
-$\bra \psi_\ti|\chi|\psi_\tf\ket = \bra \psi_\tf|\chi|\psi_\ti\ket^*$.
+$\bra \psi_\ti|\delta v|\psi_\tf\ket = \bra \psi_\tf|\delta v|\psi_\ti\ket^*$.
 Under the standard assumption
 that the incident radiation is a plane wave
 \begin{equation}\label{EPsi0Plane}
@@ -402,13 +402,13 @@ that the incident radiation is a plane wave
 with $k_\ti=K$,
 the matrix element takes the form
 \begin{equation}\label{Echiq}
-  \bra \psi_\ti|\chi|\psi_\tf\ket
-  = \int\!\d^3r\, {\rm e}^{-i\k_\ti\r}\chi(\r){\rm e}^{i\k_\tf\r}
-  = \int\!\d^3r\, {\rm e}^{i\q\r}\chi(\r)
-  \eqqcolon \chi(\q),
+  \bra \psi_\ti|\delta v|\psi_\tf\ket
+  = \int\!\d^3r\, {\rm e}^{-i\k_\ti\r}\delta v(\r){\rm e}^{i\k_\tf\r}
+  = \int\!\d^3r\, {\rm e}^{i\q\r}\delta v(\r)
+  \eqqcolon v(\q),
 \end{equation}
-\nomenclature[1χ030 2q040]{$\chi(\v{q})$}
-  {Fourier transform of the perturbation potential $\chi(\r)$}%
+\nomenclature[1χ030 2q040]{$v(\v{q})$}{Fourier
+   transform of the perturbation potential $\delta v(\r)$}%
 where we have introduced the \E{scattering vector}\footnote
 {With this choice of sign,
 \index{Sign convention!scattering vector}%
@@ -429,7 +429,7 @@ which can only be achieved by the convention~\cref{Eq}.}
   \q\coloneqq \k_\tf-\k_\ti
 \end{equation}
 \nomenclature[2q040]{$\q$}{Scattering vector}%
-and the notation $\chi(\q)$ for
+and the notation $v(\q)$ for
 the Fourier transform of the perturbative potential,
 \index{Optical potential!Fourier transform}%
 which is what small-angle neutron scattering basically measures.
@@ -450,7 +450,7 @@ We define it in arbitrary relative units as
 \begin{equation}\label{EdefJ}
   \v{J}(\r) \coloneqq  \psi^*\frac{\Nabla}{2i}\psi - \psi\frac{\Nabla}{2i}\psi^*.
 \end{equation}
-\nomenclature[2j150 2r040]{$\v{J}(\r)$}{Probability flux}
+\nomenclature[2j150 2r040]{$\v{J}(\r)$}{Probability flux}%
 \index{Flux!incident and scattered}%
 The ratio of the scattered flux hitting an infinitesimal detector area
 $r^2\d\Omega$ to the incident flux is expressed as a
@@ -475,7 +475,7 @@ With \cref{EsandwichC}, the scattered flux at the detector is
 \begin{equation}\label{EJr}
   \v{J}(\r)
   = \v{\hat r}\frac{K}{r^2}
-    {\left|\bra\psi_\ti|\chi|\psi_\tf\ket\right|}^2.
+    {\left|\bra\psi_\ti|\delta v|\psi_\tf\ket\right|}^2.
 \end{equation}
 From \cref{Exsectiondef} we obtain
 the generic differential cross section
@@ -483,7 +483,7 @@ of elastic scattering in first order Born approximation,
 \Emph{
 \begin{equation}\label{Exsection}
   \frac{\d\sigma}{\d\Omega}
-  =  {\left|\bra\psi_\ti|\chi|\psi_\tf\ket\right|}^2.
+  =  {\left|\bra\psi_\ti|\delta v|\psi_\tf\ket\right|}^2.
 \end{equation}\vspace*{-5pt}
 }
 As we shall see below,
@@ -496,7 +496,7 @@ the differential cross section is just the squared modulus
 of the Fourier transform of the perturbative potential,
 \begin{equation}\label{Ecross1}
   \frac{\d\sigma}{\d\Omega}
-  = {\left| \chi(\q) \right|}^2.
+  = {\left| v(\q) \right|}^2.
 \end{equation}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -670,7 +670,7 @@ that $q_x$ goes with the square of the two other components of the scattering ve
 \end{equation}
 Therefore, under typical small angle conditions $y,z\to L$
 the dependence of the scattering signal on $q_x$ is unimportant:
-one basically measures $\chi(\q)\simeq \chi(0,q_y,q_z)$.
+one basically measures $v(\q)\simeq v(0,q_y,q_z)$.
 The exception, for sample structures with long correlations in $x$~direction,
 is illustrated in~\cref{Fdetbox}.
 
diff --git a/Doc/UserManual/Setup.tex b/Doc/UserManual/Setup.tex
index 31f5b1ee9a6a26685fd0ca6d02a535ffcd065c7f..4f30a65b17fc6f92b3e77be6435f83201d8a8d3c 100644
--- a/Doc/UserManual/Setup.tex
+++ b/Doc/UserManual/Setup.tex
@@ -198,7 +198,7 @@
 \makenomenclature
 \renewcommand{\nomname}{List of Symbols}
   % see nomencl.txt for how to force the ordering of symbols
-
+\def\nompageref#1{,~\hyperpage{#1}\nomentryend\endgroup}
 %-------------------------------------------------------------------------------
 %  Improve LaTeX basics
 %-------------------------------------------------------------------------------
@@ -288,8 +288,6 @@
 \defineBox{boxWarn}{boxxWarn}{magenta!40}{magenta}
   {\marginSymbolLarge{fig/icons/Achtung.png}{WARN}}
 \defineBox{boxNote}{boxxNote}{yellow!33}{yellow}{{}}
-\defineBox{boxTuto}{boxxTuto}{blue!25}{red}
-  {\marginSymbolLarge{fig/icons/Arbeiten.png}{TODO}}
 \defineBox{boxEmph}{boxxEmph}{green!20}{green}{{}}
 \defineBox{boxLink}{boxxLink}{blue!25}{blue}
   {\marginSymbolLarge{fig/icons/Weblink.png}{LINK}}
@@ -297,8 +295,8 @@
 \def\Warn#1{\begin{boxWarn}#1\end{boxWarn}}
 \def\Work#1{\begin{boxWork}#1\end{boxWork}}
 \def\Note#1{\begin{boxNote}#1\end{boxNote}}
-\def\Tuto#1{\begin{boxTuto}#1\end{boxTuto}}
 \def\Link#1{\begin{boxLink}#1\end{boxLink}}
+\let\Tuto=\Link
 \def\Emph#1{\begin{boxEmph}#1\end{boxEmph}}
 \def\Emphc#1{\begin{boxEmph}#1\vskip -5pt\end{boxEmph}}
 
diff --git a/Doc/UserManual/Usage.tex b/Doc/UserManual/Usage.tex
index d83da4839cae7a1f462c36f794f737e8e9c1fc26..b0808f032849b2a21de9ac8826dc311103bbe39e 100644
--- a/Doc/UserManual/Usage.tex
+++ b/Doc/UserManual/Usage.tex
@@ -130,7 +130,7 @@ We implement two different shapes of particles: cylinders and
 prisms (\idest  elongated particles with a constant equilateral triangular cross section).
 
 All particles implemented in \BornAgain\ are defined by their
-form factors (see \cref{app:ff}), their sizes and the material
+form factors (see \cref{SFF}), their sizes and the material
 they are made of. Here, for the
 cylindrical particle, we input its radius and height.  For the prism,
 the possible inputs are the length of one side of its equilateral triangular
@@ -237,9 +237,8 @@ The first stage is to create the \Code{Simulation()} object (line~\ref{run2}). T
 parameters (line~\ref{runbeam}). %, which are associated with the
 %sample previously defined (line~\ref{runsample}). Finally we run
 %the simulation (line~\ref{runsimul}).
-Those functions are part of the Simulation
-class.  The different incident and exit angles are
-shown in \cref{fig:multil3d}.
+Those functions are part of the Simulation class.
+The different incident and exit angles are shown in \cref{fig:multil3d}.
 
 The detector parameters are set using ranges of angles via
 the function:
diff --git a/Doc/UserManual/polMultilayers.tex b/Doc/UserManual/polMultilayers.tex
deleted file mode 100644
index cb94cc327c920150bdbd2648726dbd88d2c719e0..0000000000000000000000000000000000000000
--- a/Doc/UserManual/polMultilayers.tex
+++ /dev/null
@@ -1,164 +0,0 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% \newpage % TEMPORARY
-\section{Reflection with polarization-dependent interactions [TO REVISE]}\label{s:pol}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-\MissingSection
-\index{Multilayer|)}%
-
-\Work{\indent To be reinserted once
- the chapter on polarized waves propagation is written.}
-%\iffalse
-
-%\cite{Deak_ppt, PhysRevB.76.224420, Deak2001113, PhysRevB.53.6158}.
-%\cite{RevModPhys.23.287}
-
-%-------------------------------------------------------------------------------
-\subsection{Wave equation and propagation within one layer}
-%-------------------------------------------------------------------------------
-
-To allow for polarization-dependent interactions,
-we replace the squared index of refraction $n^2$
-by $1+\uu\chi$, where $\uu\chi$ is a $2\times 2$ susceptibility matrix.
-The wave equation \cref{Escalar_wave} for layer~$\il$ becomes
-\begin{equation}\label{Ewaveqp}
-(\Delta +K^2 +K^2 \uu\chi_\il) \u\psi(\r)= 0,
-\end{equation}
-where $\u\psi(\r)$ is a two-component spinor wavefunction,
-with components $\psi_\UP(\r)$ and~$\psi_\DN(\r)$.
-At interfaces between layers,
-both spinor components of $\u\psi(\r)$ and $\Nabla\u\psi(\r)$
-must evolve continuously.
-
-The reasons for the factorization \cref{Ewave3} still apply,
-and so we can write
-\begin{equation}\label{Ewave3p}
-\u\psi(\r) = \u\psi(z) \e^{i \k_\parallel\r_\parallel}.
-\end{equation}
-As before, $\k_\parallel$ is constant across layers.
-The wave equation~\cref{Ewaveqp} reduces to
-\begin{equation}\label{Ewavezp}
-\left(\partial_z^2 + K^2 + K^2\uu\chi_\il - k_\parallel^2 \right) \u\psi(z) = 0.
-\end{equation}
-We abbreviate
-\begin{equation}
-  \uu H_\il \coloneqq  K^2(1+\uu\chi_\il)-k_\parallel^2
-\end{equation}
-so that the wave equation becomes simply
-\begin{equation}\label{Ewaveqp2}
-  \left(\partial_z^2 + \uu H_\il\right) \u\psi(z) = 0.
-\end{equation}
-The solution is
-\begin{equation}\label{Epsizp}
-  \u\psi_\il(z)
-  = \sum_{k=1}^2 \u x_{\il k}\left(\alpha_{\il k}\e^{i p_{\il k}(z-z_k)}
-                            + \beta_{\il k}\e^{-i p_{\il k}(z-z_k)}\right),
-\end{equation}
-where the $\u x_{\il k}$ are eigenvectors of $\uu H_\il$
-with eigenvalues $p_{\il k}^2$:
-\begin{equation}
-  \left( -p_{\il k}^2 + \uu H_\il \right) \u x_{\il k} = 0
-   \;\text{ for }\;\il=1,2.
-\end{equation}
-In a reproducible algorithm,
-the eigenvectors $\u x_{\il k}$ must be chosen according to some arbitrary
-normalization rule,
-for instance
-\begin{equation}
-  |\u x_{\il k}|=1,\quad x_{i\il\UP} \text{ real and nonnegative}.
-\end{equation}
-Similarly,
-a rule is needed how to handle the case of one degenerate eigenvalue,
-which includes in particular the case of scalar interactions.
-
-
-%-------------------------------------------------------------------------------
-\subsection{Wave propagation across layers}
-%-------------------------------------------------------------------------------
-
-Generalizing \cref{Evecc},
-we introduce the coefficient vector
-\begin{equation}
-  c_\il \coloneqq  {(\alpha_{\il1}, \alpha_{\il2}, \beta_{\il1}, \beta_{\il2})}^\text{T}.
-\end{equation}
-To match solutions for neighboring layers,
-continuity is requested for both spinorial components
-of $\u\psi$ and $\Nabla\u\psi$.
-As before \cref{EFcFDc}, we have at the bottom of layer~$\il$
-\begin{equation}\label{EFcFDcp}
-  F_\il c_\il = F_{\il+1} D_{\il+1} c_{\il+1},
-\end{equation}
-where the matrices are now
-\begin{equation}
-  F_\il \coloneqq  \left(\begin{array}{cccc}
-    x_{i1\UP}      &x_{i2\UP}     &x_{i1\UP}       &x_{i2\UP}       \\
-    x_{i1\DN}      &x_{i2\DN}     &x_{i1\DN}       &x_{i2\DN}       \\
-    x_{i1\UP}p_{\il1}&x_{i2\UP}p_{\il2}&-x_{i1\UP}p_{\il1}&-x_{i2\UP}p_{\il2}\\
-    x_{i1\DN}p_{\il1}&x_{i2\DN}p_{\il2}&-x_{i1\DN}p_{\il1}&-x_{i2\DN}p_{\il2}
-  \end{array}\right)
-\end{equation}
-and
-\begin{equation}
-  D_\il \coloneqq  \text{diag}(\delta_{\il1}, \delta_{\il2}, \delta_{\il1}^*, \delta_{\il2}^*)
-\end{equation}
-with the phase factor
-\begin{equation}
-   \delta_{\il k} \coloneqq  \e^{ip_{\il k}d_k}.
-\end{equation}
-Note that matrix $F_\il$ has the block form
-\begin{equation}
-  F_\il
-  =\left(\begin{array}{ll}\uu x_\il&\hphantom{-}\uu x_\il\\[1ex]
-    \uu x_\il\; \uu P_\il&-\uu x_\il\; \uu P_\il\end{array}\right)
-    = \uu x_\il \cdot
-    \left(\begin{array}{cc}\uu 1&\uu 1\\[1ex]
-    \uu P_\il&-\uu P_\il\end{array}\right),
-\end{equation}
-with
-\begin{equation}
-  \uu x_\il \coloneqq
-  \left(\u x_{\il1}, \u x_{\il2}\right),
-  \quad
-  \uu P_\il \coloneqq
-  \text{diag}\left(p_{\il1},p_{\il2}\right).
-\end{equation}
-This facilitates the computation of the inverse
-\begin{equation}
-  F_\il^{-1}
-    = \frac{1}{2}
-    \left(\begin{array}{cc}\uu 1&\hphantom{-}\uu P_\il^{-1}\\[1.2ex]
-      \uu 1 &-\uu P_\il^{-1}\end{array}\right)
-      \cdot\uu x_\il^{-1},
-\end{equation}
-which is needed for the transfer matrix $M_\il$,
-defined as in \cref{Edef_M}.
-With the new meaning of $c_\il$ and $M_\il$,
-the recursion \cref{EcMc} and the explicit solution~\cref{Eci}
-hold as derived above.
-To resolve~\cref{Eci} for the reflected amplitudes $\alpha_{0\il}$
-as function of the incident amplitudes $\beta_{0\il}$,
-we choose the notations
-\begin{equation}
-  \u\alpha_\il
-  \coloneqq \left(\begin{array}{c}\alpha_{\il1}\\\alpha_{\il2}\end{array}\right),\quad
-  \u\beta_\il
-  \coloneqq \left(\begin{array}{c}\beta_{\il1}\\\beta_{\il2}\end{array}\right),\quad
-  M\coloneqq M_1 ... M_N % TODO restore \cdots
-  \eqqcolon \left(\begin{array}{cc}\uu m_{11}&\uu m_{12}\\
-                           \uu m_{21}&\uu m_{22}\end{array}\right),
-\end{equation}
-where the $\uu m_{\il k}$ are $2\times2$ matrices.
-Eq.~\cref{Eci} then takes the form
-\begin{equation}
-  \left(\begin{array}{c}\u\alpha_{0}\\\u\beta_{0}\end{array}\right)
-  =
-  \left(\begin{array}{cc}\uu m_{11}&\uu m_{12}\\
-    \uu m_{21}&\uu m_{22}\end{array}\right)
-  \left(\begin{array}{c}\u{0}\\\u\beta_{N}\end{array}\right),
-\end{equation}
-which immediately yields
-\begin{equation}
-  \u\alpha_0 = \uu m_{12}\,\uu m_{22}^{-1}\,\u\beta_0.
-\end{equation}
-
-%\fi