diff --git a/Doc/DevRef/DevRef.pdf b/Doc/DevRef/DevRef.pdf
index aceec5790929aaa1b7057791e92b70b3a2a78aa9..dc85f8c1e568d89cd791b39db9461ce4a7b68d72 100644
Binary files a/Doc/DevRef/DevRef.pdf and b/Doc/DevRef/DevRef.pdf differ
diff --git a/Doc/DevRef/Polarized.tex b/Doc/DevRef/Polarized.tex
index b7f469e2884b581530318653668b9411d8698461..5f5d83803d0b8afed30b346cc60d7a32a43b747f 100644
--- a/Doc/DevRef/Polarized.tex
+++ b/Doc/DevRef/Polarized.tex
@@ -143,6 +143,11 @@ By comparison with~\cref{EWaveZ2}, we see that the square of the operator~$\opka
 \begin{equation}\label{Dhp2}
   \opkappa^2 = \kappa_0^2-4\pi \hat \mv = \kappa_0^2-4\pi(\sldN + \sldM \eB \Pauli).
 \end{equation}
+
+%===============================================================================
+\subsection{Wavenumber operator}
+%===============================================================================
+
 Without derivation,\footnote
 {To verify, use standard properties of Pauli matrices.
 Square~\cref{Ehp1} to reproduce~\cref{Dhp2}.
@@ -214,7 +219,7 @@ with the $4\times4$ transfer matrix in place of~\cref{EMil}
 \begin{equation}\label{EMDSP}
   \WW{M}_l \coloneqq \WW{D}_{l-1}\, \WW{S}_{l}.
 \end{equation}
-The phase rotation matrix~\cref{DmatD} becomes a $4\times4$ block matrix
+The phase rotation matrix~\cref{DmatD} is replaced by the block matrix
 \begin{equation}
   \WW{D}_l
    \coloneqq
@@ -223,7 +228,8 @@ The phase rotation matrix~\cref{DmatD} becomes a $4\times4$ block matrix
      0&\hat{\delta_l}
    \end{pmatrix},
 \end{equation}
-and refraction matrix~\cref{DmatS} becomes
+to be discussed in the next section.
+The refraction matrix~\cref{DmatS} also is replaced by a block matrix,
 \begin{equation}\label{ESabP}
   \WW{S}_{l}
    \coloneqq
@@ -233,13 +239,17 @@ and refraction matrix~\cref{DmatS} becomes
        \hat s^-_l & \hat s^+_l
    \end{array}\right)
 \end{equation}
-with the coefficients
+with the coefficients\footnote
+{Currently (jun23), the matrix blocks $\hat s^+_l$ and $\hat s^-_l$,
+possibly modified by roughness factors (see below .... TODO),
+are computed through local function \T{refractionMatrixBlocks}
+in \SRC{Resample/Specular}{ComputeFluxMagnetic.cpp}.}
 \begin{equation}
   \hat s_l\pm \coloneqq 1 \pm \opkappa_{l-1}^{-1}\,\opkappa_l.
 \end{equation}
 
 %==================================================================================================%
-\subsection{Evaluation of the phase rotation matrix}\label{Sphase}
+\subsection{Phase rotation matrix}\label{Sphase}
 %==================================================================================================%
 
 The matrix $\eB\Pauli$ has the eigenvalues $\pm 1$ and the normalized eigenspinors
@@ -266,7 +276,8 @@ Then $\opkappa$ has the eigenvalue decomposition
   = \hat{Q} \, \begin{pmatrix}\evp_+ & 0 \\ 0 & \evp_-\end{pmatrix} \hat{Q}^\dagger,
 \end{equation}
 and the phase rotation matrix \cref{Dopdel} can be written\footnote
-{Currently (jun23) implemented in funtion \T{MatrixRTCoefficients::computeDeltaMatrix}.}
+{Currently (jun23) implemented in local function \T{PhaseRotationMatrix}
+in file \SRC{Resample/Flux}{MatrixFlux.cpp}.}
 \begin{equation}\label{EdP2}
   \hat\delta
   = \e^{i\opkappa d}
@@ -290,112 +301,6 @@ the critical factor $\e^{i\alpha d/2}$ may be drawn in front of $\hat{Q}$ in~\cr
 \section{From old document ``PolarizedImplementation''}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-%===============================================================================
-\subsection{Intensity analysis}\label{sec:intensity_analysis}
-%===============================================================================
-
-The wave function in the ambient material is given by
-\begin{equation}
-\U{\Phi_0} \left(z\right) = \underbrace{ \exp{i \UU{p_0}} \U{t_0} }_{\U{\Phi_i}(z)} +
- \underbrace{ \exp{-i \UU{p_0} } \UU{r_0} }_{\U{\Phi_r}(z)} \,.
-\end{equation}
-Here $\U{\Phi_i}$ is the incoming wave in a given polarization state $\U{t_0}$ and $\U{\Phi_r}$ is the reflected wave.
-The intensity measured on a detector is without polarization analysis given by
-\begin{equation} \label{eq:intensity_no_polarization_analysis}
-I_R = \left| \U{\Phi_r} \right|^2 = \Braket{\U{\Phi_r} | \U{\Phi_r}}
-\end{equation}
-This would allow for the direct computation of the reflected intensity if $\U{t_0}$ would describe the incoming polarization state.
-In case of an arbitrary but pure state of the incoming beam, the reflected wave can be described by a reflection matrix
-\begin{equation}
-\U{\Phi_r}(z) = \UU{\mathcal{R}} \, \U{\Phi_i} (z) \,,
-\end{equation}
-where $\UU{\mathcal{R}}$ is a $2 \times 2$ matrix:
-\begin{equation}
-\UU{\mathcal{R}}
-=
-\begin{pmatrix}
-r_{++} & r_{-+} \\
-r_{+-} & r_{--} \\
-\end{pmatrix}
-\end{equation}
-with its non-diagonal elements contributing to spin-flip reflections.
-If we consider the intensities right a the topmost interface of the sample at $z = z_0 = 0$, the phase factors drop out and we find the relation
-\begin{equation} \label{eq:R-matrix}
-\U{r_0} = \UU{\mathcal{R}} \U{t_0} \,.
-\end{equation}
-Therefore, as soon as $\UU{\mathcal{R}}$ is known, it is trivial to perform a calculation for any desired incoming polarization state.
-It is clear that \cref{eq:R-matrix} is a system of two equations with four unknown variables.
-Hence if two pairs of incoming and reflected waves $\U{t_0}, \U{r_0}$ and $\U{t_0'}, \U{r_0'}$ are known, the reflection matrix can be determined (if they are linearly independent).
-Writing these four equations as a matrix
-\begin{equation} \label{eq:R-matrix2}
-\left( \U{r_0}, \U{r_0}' \right) = \UU{\mathcal{R}} \left( \U{t_0}, \U{t_0}' \right) \,,
-\end{equation}
-one can see that the inversion of this equation becomes trivial if the two incoming waves are chosen such that they are in a $+$ and $-$ polarization state.
-If they are chosen differently, inverting equation \cref{eq:R-matrix2} corresponds to the rotation of the incoming polarization vectors such that they become pure $+$ and $-$ waves.
-
-
-If we perform polarization analysis, the analyzer will only pass a wave in the $\U{\Phi_f}$ polarization state.
-Hence the reflected wave needs to be projected onto this state to obtain the measured intensity
-\begin{equation} \label{eq:intensity_polarization_analysis}
-I_R
-=
-\left| \left\langle \U{\Phi_f} | \UU{\mathcal{R}} | \U{\Phi_i} \right\rangle \right|^2
-=
-\left\langle \U{\Phi_f} | \UU{\mathcal{R}} | \U{\Phi_i} \right\rangle
-\cdot
-\left\langle \U{\Phi_i} | \UU{\mathcal{R}}^{\dagger} | \U{\Phi_f} \right\rangle \,.
-\end{equation}
-Following \cite{Top15}, we introduce the density matrix $\UU{f_p}$ (polarizer) for an arbitrary mixed state of incoming beam and, correspondingly, $\UU{f_a}$ (analyzer) for an arbitrary mixed state passed through a polarization analyzer:
-\begin{align} \label{eq:density_operators}
-\UU{f_p} &= \frac{1}{2} \left(
-\UU{1} + \Pauli \cdot \boldsymbol{p}
-\right)
-&
-\UU{f_a} &= \frac{1}{2} \left(
-\UU{1} + \Pauli \cdot \boldsymbol{a}
-\right).
-\end{align}
-The beam polarization as well as analyzer direction and efficiency are described by the Bloch vectors $\boldsymbol{p}, \boldsymbol{a} \in \mathbb{R}^3$.
-$|\boldsymbol{p}| = 1$ corresponds to some pure state of beam polarization, while $|\boldsymbol{p}| < 1$ is for a state mixture (partial polarization).
-The same holds for non-perfect analyzers, where we call $|\boldsymbol{p}|$ the efficiency.
-\footnote{In PolarizedSpecular, Sec. 5.1, Dmitry claims this treatment of a non-perfect analyzer is not possible and suggests a different treatment.
-However i (rb) think that his argument is not correct.}
-
-In order to compute the reflection coefficient for a mixed-state beam, equation \cref{eq:intensity_polarization_analysis} needs to be rewritten in the density matrix formalism
-\begin{equation} \label{eq:intensity_trace}
-\Braket{ \U{\Phi_i} | \UU{\mathcal{R}}^{\dagger} | \Psi_f }
-\cdot
-\Braket{ \U{\Phi_f} | \UU{\mathcal{R}} | \U{\Phi_i}}
-=
-\text{Tr} \left(
-\Ket{ \U{\Phi_i} } \Bra{ \U{\Phi_i} } \UU{\mathcal{R}}^{\dagger} \Ket{ \U{\Phi_f} }  \Bra{ \U{\Phi_f}} \UU{\mathcal{R}}
-\right).
-\end{equation}
-Here Tr denotes trace operation.
-$\Ket{ \U{\Phi_f}} \Bra{ \U{\Phi_f}}$ and $ \Ket{ \U{\Phi_i}} \Bra{ \U{\Phi_i}}$ are respective outer products for the $\U{\Phi_f}$ and $\U{\Phi_i}$ pure states and coincide with the corresponding density matrices.
-To generalize expression \cref{eq:intensity_trace} to mixed states of the incoming beam and polarization analyzer, one has to replace the explicit outer products with the density matrices $\UU{f_p}$, $\UU{f_a}$ that describe the polarizer and analyzer as defined in \cref{eq:density_operators}.
-This will automatically take into account the averaging over all possible initial and final pure states of the system.
-Therefore, the final expression for $I_R$ reads
-\begin{equation} \label{eq:32}
-I_R
-=
-Tr \left(
-\UU{ f_p } \UU{\mathcal{R}}^{\dagger} \UU{f_a} \UU{\mathcal{R}}
-\right).
-\end{equation}
-This expression should work for both a perfect and imperfect polarizer and analyzer.
-It also seems to be consistent with Wildes \cite{Wil99,Wil06}, this paper was recommended as standard reference on this topic by Artur.
-
-It needs to be noted that the limit $\left| \vec a  \right| = 0$ does not correspond to no polarization analysis (i.e. a very common experiment without polarization analysis).
-Instead, if no polarization analysis is performed, writing \cref{eq:intensity_no_polarization_analysis} in the density matrix formalism yields
-\begin{equation}
-I_R = \left| \U{\Phi_r} \right|^2 = \Braket{ \U{\Phi_i} | \UU{\mathcal{R}}^\dagger \UU{\mathcal{R}} | \U{\Phi_i}}
-= \text{Tr}\left( \UU{f_p} \UU{\mathcal{R}}^\dagger \UU{\mathcal{R}} \right) \,,
-\end{equation}
-which corresponds to $\UU{f_a} = \UU{1}$.
-
-TODO: can we show the equivalence between the Wildes approach and our density matrix formulas?
-
 %===============================================================================
 \subsection{Numerically stable implementation}\label{sec:implementation}
 %===============================================================================
@@ -406,36 +311,33 @@ We combine the amplitude vectors for $+$ and $-$ polarization into a single matr
 \UU{T_j} \\ \UU{R_j}
 \end{pmatrix} &= \begin{pmatrix}
 \U{t_j}^+ & \U{t_j}^- \\ \U{r_j}^+ & \U{r_j}^-
-\end{pmatrix}\,,
-&
-\UU{T_j} &= \left( \U{t_j}^+, \U{t_j}^- \right) \,,
-&
-\UU{R_j} &= \left( \U{r_j}^+, \U{r_j}^-\right) \,,
+\end{pmatrix},
 \end{align}
-where the submatrices $\UU{T_j}$ and $\UU{R_j}$ are of dimension $2 \times 2$.
-The recursion equation \cref{eq:interface_transfer_matrix_re} can then be written simultaneously for both polarization states as
+where the submatrices
+\begin{equation}
+\UU{T_j} = \left( \U{t_j}^+, \U{t_j}^- \right) \quad\text{and}\quad
+\UU{R_j} = \left( \U{r_j}^+, \U{r_j}^-\right)
+\end{equation}
+are of dimension $2 \times 2$.
+The transfer matrix equation \cref{EcMcP}
+can then be written simultaneously for both polarization states as
 \begin{align}
 \begin{pmatrix}
 \UU{T_a} \\ \UU{R_a}
-\end{pmatrix} &= \frac{1}{2} \begin{pmatrix}
-\UU{\delta_a}^{-1} & 0 \\ 0 & \UU{\delta_a}
-\end{pmatrix}
-\begin{pmatrix}
-1 + \UU{P_{ab}} & 1 - \UU{P_{ab}}
-\\
-1 - \UU{P_{ab}} & 1 + \UU{P_{ab}}
-\end{pmatrix}
+\end{pmatrix} &= \WW{M}_{ab}
 \begin{pmatrix}
 \UU{T_{b}} \\ \UU{R_{b}}
 \end{pmatrix} \,.
 \end{align}
 Explicitly performing this multiplication yields two matrix recursion equations
 \begin{align}
-\UU{T_a}  &= \frac{1}{2} \UU{\delta_a}^{-1} \left( \left(1 + \UU{P_{ab}} \right) \UU{T_{b}} + \left(1 - \UU{P_{ab}} \right) \UU{R_{b}} \right)
+\UU{T_a}  &= \frac{1}{2} \hat\delta^{-1} \left( \hat{s}^+ \UU{T_{b}} + \hat{s}^- \UU{R_{b}} \right)
 \\
-\UU{R_a}  &= \frac{1}{2} \UU{\delta_a} \left( \left(1 - \UU{P_{ab}} \right) \UU{T_{b}} + \left(1 + \UU{P_{ab}} \right) \UU{R_{b}} \right)
+\UU{R_a}  &= \frac{1}{2} \hat\delta \left( \hat{s}^- \UU{T_{b}} + \hat{s}^+ \UU{R_{b}} \right)
 \end{align}
-After every step of the iteration \cref{eq:total_transfer_matrix}, we want to rotate and normalize the polarization, such that we obtain the bottom boundary condition $\UU{T_a} = \UU{1}$.
+After every step~\cref{EcMcP},
+we want to rotate and normalize the polarization,
+such that we obtain the bottom boundary condition $\UU{T_a} = \UU{1}$.
 For this purpose, we define the rotation matrix $\UU{S}$ as
 \begin{align}
 \UU{T'_a} &= \UU{T_a} \cdot \UU{S_a} = \UU{1} \,,
@@ -443,7 +345,8 @@ For this purpose, we define the rotation matrix $\UU{S}$ as
 and it must also be applied to rotate the reflected components as well
 \footnote{What happens here is a rotation of the wave function written as a superposition
  $\U{\Phi'_a}^\pm = a^\pm \U{\Phi_a}^+ + b^\pm \U{\Phi_a}^- $.
-This is just written as a single matrix equation and from this it is obvious that the same matrix must be applied to both $\UU{T_a}$ and $\UU{R_a}$.}
+This is just written as a single matrix equation
+and from this it is obvious that the same matrix must be applied to both $\UU{T_a}$ and $\UU{R_a}$.}
 \begin{align}
 \UU{R'_a} &= \UU{R_a} \cdot \UU{S_a} \,.
 \end{align}
@@ -452,9 +355,9 @@ Consequently, the recursion equations reduce to the primed version
 \begin{align}
 \UU{T'_a}  &= \UU{1}
 \\
-\UU{S_a}^{-1} &= \frac{1}{2} \UU{\delta_a}^{-1} \left( \left(1 + \UU{P_{ab}} \right) + \left(1 - \UU{P_{ab}} \right) \UU{R'_{b}} \right) \eqqcolon \UU{\delta_a}^{-1} \cdot \UU{\tilde{S}_a}^{-1}
+\UU{S_a}^{-1} &= \frac{1}{2} \hat\delta^{-1} \left( \hat{s}^+ + \hat{s}^- \UU{R'_{b}} \right) \eqqcolon \hat\delta^{-1} \cdot \UU{\tilde{S}_a}^{-1}
 \\
-\UU{R'_a}  &= \frac{1}{2} \UU{\delta_a} \left( \left(1 - \UU{P_{ab}} \right) + \left(1 + \UU{P_{ab}} \right) \UU{R'_{b}} \right) \cdot \UU{S_a} \,.
+\UU{R'_a}  &= \frac{1}{2} \hat\delta \left( \hat{s}^- + \hat{s}^+ \UU{R'_{b}} \right) \cdot \UU{S_a} \,.
 \end{align}
 \end{subequations}
 This process requires the inverse of $ \UU{S_a}^{-1}$, which is easy to obtain
@@ -463,11 +366,13 @@ This process requires the inverse of $ \UU{S_a}^{-1}$, which is easy to obtain
 \begin{pmatrix}
 \left( \UU{\tilde{S}_a}^{-1}\right)_{1, 1} & -\left(\UU{\tilde{S}_a}^{-1}\right)_{1, 0} \\
  -\left(\UU{\tilde{S}_a}^{-1}\right)_{0, 1} & \left(\UU{\tilde{S}_a}^{-1}\right)_{0, 0}
-\end{pmatrix} \UU{\delta_a}  \,,
+\end{pmatrix} \hat\delta  \,,
 \end{align}
-and since it does not contain the inverse $\delta^{-1}$-matrix anymore can be evaluated numerically stable.
+and since it does not contain the inverse $\delta^{-1}$-matrix anymore
+can be evaluated numerically stable.
 
-For the computation of reflectivity alone, we would be done at this point, however, for perturbation theory, we still need the amplitudes within the layer stack.
+For the computation of reflectivity alone, we would be done at this point,
+however, for perturbation theory, we still need the amplitudes within the layer stack.
 Therefore, all rotations of the polarization need to be forward-propagated according to
 \begin{align}
 \UU{T_a} &= \prod^{0}_{i = m-1} \UU{S_i}
@@ -477,7 +382,9 @@ Therefore, all rotations of the polarization need to be forward-propagated accor
 \text{for } & a \geq 1, \text{ where } \prod^{0}_{i = m-1} \UU{S_i} = \UU{S_{a-1}} \UU{S_{a-2}} \cdots \UU{S_0}
 \end{align}
 The proof of these relations can be done as follows.
-In order to obtain the correct amplitudes in each layer, after every step of the recursion \cref{eq:primed_recursion} the applied rotation needs to be propagated down through the bottom of the stack according to
+In order to obtain the correct amplitudes in each layer,
+after every step of the recursion \cref{eq:primed_recursion}
+the applied rotation needs to be propagated down through the bottom of the stack according to
 \begin{align}
 & \text{for } m = N - 1 \dots 0  \hspace{2em} \text{\small outer iteration from bottom to top}
 \\
@@ -494,13 +401,12 @@ In order to obtain the correct amplitudes in each layer, after every step of the
 
 
 \textbf{Remark}
-The defined $\UU{R'_m}$ is equal to the $\UU{X_m}$ in the Parratt formalism and equation \cref{eq:primed_recursion}c is almost identical to the usual recursion in $\UU{X_m}$ \cite{KeRT08,KeRT03}, apart from the treatment of the phase factor.
-The second recursion equation in \cite{KeRT08,KeRT03} that yields the amplitudes is replaced by storing the intermediate $\UU{S_m}$.
-
-  \T{Compute::SpecularMagnetic::topLayerR}
-  contains a lightweight implementation that does not store coefficients of intermediate layers and hence only computes the coefficient of the top layer and hence reflection.
-  This was introduced in order to speed up pure reflectometry computations.
-  This is similarly implemented for the scalar implementation.
+The defined $\UU{R'_m}$ is equal to the $\UU{X_m}$ in the Parratt formalism,
+and equation \cref{eq:primed_recursion}c is almost identical
+to the usual recursion in $\UU{X_m}$ \cite{KeRT08,KeRT03},
+apart from the treatment of the phase factor.
+The second recursion equation in \cite{KeRT08,KeRT03}
+that yields the amplitudes is replaced by storing the intermediate $\UU{S_m}$.
 
 %===============================================================================
 \subsection{Roughness}
@@ -798,17 +704,111 @@ Both of these conditions result in the $z$-component of the magnetic field (that
 \end{equation*}
 
 %===============================================================================
-\subsection{Further (potential) problems}
+\subsection{Intensity and density matrix}\label{sec:intensity_analysis}
 %===============================================================================
 
-\begin{itemize}
- \item $k_0^2 = 4 \pi \sldN$
- \item Many layers leading to infinity
-  \item Zero k-vector inside sample at critical angle? Probably not an issue, resolved due to checkforunderflow in kzcomputation.cpp
-  \item Branch cut from complex square-root?
-\end{itemize}
+The wave function in the ambient material is given by
+\begin{equation}
+\U{\Phi_0} \left(z\right) = \underbrace{ \exp{i \UU{p_0}} \U{t_0} }_{\U{\Phi_i}(z)} +
+ \underbrace{ \exp{-i \UU{p_0} } \UU{r_0} }_{\U{\Phi_r}(z)} \,.
+\end{equation}
+Here $\U{\Phi_i}$ is the incoming wave in a given polarization state $\U{t_0}$ and $\U{\Phi_r}$ is the reflected wave.
+The intensity measured on a detector is without polarization analysis given by
+\begin{equation} \label{eq:intensity_no_polarization_analysis}
+I_R = \left| \U{\Phi_r} \right|^2 = \Braket{\U{\Phi_r} | \U{\Phi_r}}
+\end{equation}
+This would allow for the direct computation of the reflected intensity if $\U{t_0}$ would describe the incoming polarization state.
+In case of an arbitrary but pure state of the incoming beam, the reflected wave can be described by a reflection matrix
+\begin{equation}
+\U{\Phi_r}(z) = \UU{\mathcal{R}} \, \U{\Phi_i} (z) \,,
+\end{equation}
+where $\UU{\mathcal{R}}$ is a $2 \times 2$ matrix:
+\begin{equation}
+\UU{\mathcal{R}}
+=
+\begin{pmatrix}
+r_{++} & r_{-+} \\
+r_{+-} & r_{--} \\
+\end{pmatrix}
+\end{equation}
+with its non-diagonal elements contributing to spin-flip reflections.
+If we consider the intensities right a the topmost interface of the sample at $z = z_0 = 0$, the phase factors drop out and we find the relation
+\begin{equation} \label{eq:R-matrix}
+\U{r_0} = \UU{\mathcal{R}} \U{t_0} \,.
+\end{equation}
+Therefore, as soon as $\UU{\mathcal{R}}$ is known, it is trivial to perform a calculation for any desired incoming polarization state.
+It is clear that \cref{eq:R-matrix} is a system of two equations with four unknown variables.
+Hence if two pairs of incoming and reflected waves $\U{t_0}, \U{r_0}$ and $\U{t_0'}, \U{r_0'}$ are known, the reflection matrix can be determined (if they are linearly independent).
+Writing these four equations as a matrix
+\begin{equation} \label{eq:R-matrix2}
+\left( \U{r_0}, \U{r_0}' \right) = \UU{\mathcal{R}} \left( \U{t_0}, \U{t_0}' \right) \,,
+\end{equation}
+one can see that the inversion of this equation becomes trivial if the two incoming waves are chosen such that they are in a $+$ and $-$ polarization state.
+If they are chosen differently, inverting equation \cref{eq:R-matrix2} corresponds to the rotation of the incoming polarization vectors such that they become pure $+$ and $-$ waves.
 
 
+If we perform polarization analysis, the analyzer will only pass a wave in the $\U{\Phi_f}$ polarization state.
+Hence the reflected wave needs to be projected onto this state to obtain the measured intensity
+\begin{equation} \label{eq:intensity_polarization_analysis}
+I_R
+=
+\left| \left\langle \U{\Phi_f} | \UU{\mathcal{R}} | \U{\Phi_i} \right\rangle \right|^2
+=
+\left\langle \U{\Phi_f} | \UU{\mathcal{R}} | \U{\Phi_i} \right\rangle
+\cdot
+\left\langle \U{\Phi_i} | \UU{\mathcal{R}}^{\dagger} | \U{\Phi_f} \right\rangle \,.
+\end{equation}
+Following \cite{Top15}, we introduce the density matrix $\UU{f_p}$ (polarizer) for an arbitrary mixed state of incoming beam and, correspondingly, $\UU{f_a}$ (analyzer) for an arbitrary mixed state passed through a polarization analyzer:
+\begin{align} \label{eq:density_operators}
+\UU{f_p} &= \frac{1}{2} \left(
+\UU{1} + \Pauli \cdot \boldsymbol{p}
+\right)
+&
+\UU{f_a} &= \frac{1}{2} \left(
+\UU{1} + \Pauli \cdot \boldsymbol{a}
+\right).
+\end{align}
+The beam polarization as well as analyzer direction and efficiency are described by the Bloch vectors $\boldsymbol{p}, \boldsymbol{a} \in \mathbb{R}^3$.
+$|\boldsymbol{p}| = 1$ corresponds to some pure state of beam polarization, while $|\boldsymbol{p}| < 1$ is for a state mixture (partial polarization).
+The same holds for non-perfect analyzers, where we call $|\boldsymbol{p}|$ the efficiency.
+\footnote{In PolarizedSpecular, Sec. 5.1, Dmitry claims this treatment of a non-perfect analyzer is not possible and suggests a different treatment.
+However i (rb) think that his argument is not correct.}
+
+In order to compute the reflection coefficient for a mixed-state beam, equation \cref{eq:intensity_polarization_analysis} needs to be rewritten in the density matrix formalism
+\begin{equation} \label{eq:intensity_trace}
+\Braket{ \U{\Phi_i} | \UU{\mathcal{R}}^{\dagger} | \Psi_f }
+\cdot
+\Braket{ \U{\Phi_f} | \UU{\mathcal{R}} | \U{\Phi_i}}
+=
+\text{Tr} \left(
+\Ket{ \U{\Phi_i} } \Bra{ \U{\Phi_i} } \UU{\mathcal{R}}^{\dagger} \Ket{ \U{\Phi_f} }  \Bra{ \U{\Phi_f}} \UU{\mathcal{R}}
+\right).
+\end{equation}
+Here Tr denotes trace operation.
+$\Ket{ \U{\Phi_f}} \Bra{ \U{\Phi_f}}$ and $ \Ket{ \U{\Phi_i}} \Bra{ \U{\Phi_i}}$ are respective outer products for the $\U{\Phi_f}$ and $\U{\Phi_i}$ pure states and coincide with the corresponding density matrices.
+To generalize expression \cref{eq:intensity_trace} to mixed states of the incoming beam and polarization analyzer, one has to replace the explicit outer products with the density matrices $\UU{f_p}$, $\UU{f_a}$ that describe the polarizer and analyzer as defined in \cref{eq:density_operators}.
+This will automatically take into account the averaging over all possible initial and final pure states of the system.
+Therefore, the final expression for $I_R$ reads
+\begin{equation} \label{eq:32}
+I_R
+=
+Tr \left(
+\UU{ f_p } \UU{\mathcal{R}}^{\dagger} \UU{f_a} \UU{\mathcal{R}}
+\right).
+\end{equation}
+This expression should work for both a perfect and imperfect polarizer and analyzer.
+It also seems to be consistent with Wildes \cite{Wil99,Wil06}, this paper was recommended as standard reference on this topic by Artur.
+
+It needs to be noted that the limit $\left| \vec a  \right| = 0$ does not correspond to no polarization analysis (i.e. a very common experiment without polarization analysis).
+Instead, if no polarization analysis is performed, writing \cref{eq:intensity_no_polarization_analysis} in the density matrix formalism yields
+\begin{equation}
+I_R = \left| \U{\Phi_r} \right|^2 = \Braket{ \U{\Phi_i} | \UU{\mathcal{R}}^\dagger \UU{\mathcal{R}} | \U{\Phi_i}}
+= \text{Tr}\left( \UU{f_p} \UU{\mathcal{R}}^\dagger \UU{\mathcal{R}} \right) \,,
+\end{equation}
+which corresponds to $\UU{f_a} = \UU{1}$.
+
+TODO: can we show the equivalence between the Wildes approach and our density matrix formulas?
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{From old document ``Reflectivity''}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Resample/Flux/MatrixFlux.cpp b/Resample/Flux/MatrixFlux.cpp
index ed9ee55df27bd50ccbd433aeaf780afac573ea55..0eb7db48c77d70d4e6bf8be18f5ae5e716154b76 100644
--- a/Resample/Flux/MatrixFlux.cpp
+++ b/Resample/Flux/MatrixFlux.cpp
@@ -26,6 +26,21 @@ complex_t GetImExponential(complex_t exponent)
     return std::exp(I * exponent);
 }
 
+//! See DevRef, chapter "Polarized", section "Phase rotation matrix"
+SpinMatrix PhaseRotationMatrix(const R3& B_direction, const Spinor& diagonal)
+{
+    const R3& e = B_direction; // direction of B field
+    const SpinMatrix M(diagonal.u, 0, 0, diagonal.v);
+
+    if (std::abs(e.mag() - 1.) < eps) {
+        const double factor1 = 2. * (1. + e.z());
+        SpinMatrix Q(1. + e.z(), I * e.y() - e.x(), e.x() + I * e.y(), e.z() + 1.);
+        return Q * M * Q.adjoint() / factor1;
+    }
+    ASSERT(e.mag() < eps); // remaining case: no field
+    return M;
+}
+
 } // namespace
 
 MatrixFlux::MatrixFlux(double kz_sign, const Spinor& eigenvalues, const R3& b, double magnetic_SLD)
@@ -40,74 +55,44 @@ MatrixFlux::MatrixFlux(double kz_sign, const Spinor& eigenvalues, const R3& b, d
 }
 
 
-SpinMatrix MatrixFlux::pMatrixHelper(double sign) const
+//! See DevRef, chapter "Polarized", section "Wavenumber operator".
+SpinMatrix MatrixFlux::computeKappa() const
 {
     const complex_t alpha = m_k_eigen.v + m_k_eigen.u;
     const complex_t beta = m_k_eigen.v - m_k_eigen.u;
-
-    SpinMatrix result(alpha + sign * beta * m_b.z(), sign * beta * (m_b.x() - I * m_b.y()),
-                      sign * beta * (m_b.x() + I * m_b.y()), alpha - sign * beta * m_b.z());
-
-    return m_kz_sign * result;
-}
-
-SpinMatrix MatrixFlux::computeKappa() const
-{
-    return pMatrixHelper(1.) / 2;
+    return SpinMatrix(alpha + beta * m_b.z(), beta * (m_b.x() - I * m_b.y()),
+                      beta * (m_b.x() + I * m_b.y()), alpha - beta * m_b.z())
+           / 2.;
 }
 
+//! See DevRef, chapter "Polarized", section "Wavenumber operator".
 SpinMatrix MatrixFlux::computeInverseKappa() const
 {
     const complex_t alpha = m_k_eigen.v + m_k_eigen.u;
     const complex_t beta = m_k_eigen.v - m_k_eigen.u;
-
     if (std::abs(alpha * alpha - beta * beta) == 0.)
-        return {};
+        throw std::runtime_error("Inverse wavenumber operator is singular for given B field."
+                                 " Add some absorption to the nuclear SLD.");
 
-    return 2. / (alpha * alpha - beta * beta) * pMatrixHelper(-1.);
+    return 2. / (alpha * alpha - beta * beta)
+           * SpinMatrix(alpha - beta * m_b.z(), beta * (-m_b.x() + I * m_b.y()),
+                        -beta * (m_b.x() + I * m_b.y()), alpha + beta * m_b.z());
 }
 
 SpinMatrix MatrixFlux::computeDeltaMatrix(double thickness) const
 {
-    const complex_t alpha = 0.5 * thickness * (m_k_eigen.v + m_k_eigen.u);
-
-    // Compute resulting phase matrix according to exp(i p_m d_m) = exp1 * Q * exp2 * Q.adjoint();
-    if (std::abs(m_b.mag() - 1.) < eps) {
-        const SpinMatrix exp2(GetImExponential(m_kz_sign * thickness * m_k_eigen.v), 0, 0,
-                              GetImExponential(m_kz_sign * thickness * m_k_eigen.u));
-        const double factor1 = 2. * (1. + m_b.z());
-        SpinMatrix Q((1. + m_b.z()), (I * m_b.y() - m_b.x()), (m_b.x() + I * m_b.y()),
-                     (m_b.z() + 1.));
-
-        return Q * exp2 * Q.adjoint() / factor1;
-    }
-    if (m_b.mag() < eps)
-        return SpinMatrix::One() * GetImExponential(m_kz_sign * alpha);
-
-    ASSERT(false); // invalid magnetic field vector
-}
-
-SpinMatrix MatrixFlux::TransformationMatrix(const Spinor& selection) const
-{
-    const SpinMatrix exp2(selection.u, 0, 0, selection.v);
-
-    if (std::abs(m_b.mag() - 1.) < eps) {
-        const double factor1 = 2. * (1. + m_b.z());
-        SpinMatrix Q(1. + m_b.z(), I * m_b.y() - m_b.x(), m_b.x() + I * m_b.y(), m_b.z() + 1.);
-        return Q * exp2 * Q.adjoint() / factor1;
-    }
-    ASSERT(m_b.mag() < eps); // otherwise the magnetic field vector is invalid
-    return exp2;
+    return PhaseRotationMatrix(m_b, {GetImExponential(m_kz_sign * thickness * m_k_eigen.v),
+                                     GetImExponential(m_kz_sign * thickness * m_k_eigen.u)});
 }
 
 SpinMatrix MatrixFlux::T1Matrix() const
 {
-    return TransformationMatrix({0., 1.});
+    return PhaseRotationMatrix(m_b, {0., 1.});
 }
 
 SpinMatrix MatrixFlux::T2Matrix() const
 {
-    return TransformationMatrix({1., 0.});
+    return PhaseRotationMatrix(m_b, {1., 0.});
 }
 
 Spinor MatrixFlux::T1plus() const
diff --git a/Resample/Flux/MatrixFlux.h b/Resample/Flux/MatrixFlux.h
index 421f3ab845bb6d869cbde8f870f23bf143000c13..be7617229d53564a9730a17c2c8d7920e0070926 100644
--- a/Resample/Flux/MatrixFlux.h
+++ b/Resample/Flux/MatrixFlux.h
@@ -67,11 +67,8 @@ private:
     double m_magnetic_SLD;
 
     // helper functions to compute DWBA compatible amplitudes used in the T1plus() etc. functions
-    SpinMatrix TransformationMatrix(const Spinor& selection) const;
     SpinMatrix T1Matrix() const;
     SpinMatrix T2Matrix() const;
-
-    SpinMatrix pMatrixHelper(double sign) const;
 };
 
 #endif // BORNAGAIN_RESAMPLE_FLUX_MATRIXFLUX_H
diff --git a/Resample/Specular/ComputeFluxMagnetic.cpp b/Resample/Specular/ComputeFluxMagnetic.cpp
index 5cb1bf6f743077e08cfe11e375c2d9d4cf0195f8..c7af127d873102b8e9cd2612294dcb138471c7c8 100644
--- a/Resample/Specular/ComputeFluxMagnetic.cpp
+++ b/Resample/Specular/ComputeFluxMagnetic.cpp
@@ -50,10 +50,9 @@ Spinor checkForUnderflow(const Spinor& eigenvs)
     return {k_eigen(eigenvs.u), k_eigen(eigenvs.v)};
 }
 
-std::pair<SpinMatrix, SpinMatrix> computeBackwardsSubmatrices(const MatrixFlux& coeff_i,
-                                                              const MatrixFlux& coeff_i1,
-                                                              double sigma,
-                                                              const RoughnessModel& r_model)
+std::pair<SpinMatrix, SpinMatrix> refractionMatrixBlocks(const MatrixFlux& coeff_i,
+                                                         const MatrixFlux& coeff_i1, double sigma,
+                                                         const RoughnessModel& r_model)
 {
     if (r_model == RoughnessModel::NEVOT_CROCE)
         return Compute::MagneticNevotCroceTransition::backwardsSubmatrices(coeff_i, coeff_i1,
@@ -78,7 +77,7 @@ void calculateUpwards(std::vector<MatrixFlux>& coeff, const SliceStack& slices,
             sigma = roughness->sigma();
 
         // compute the 2x2 submatrices in the write-up denoted as P+, P- and delta
-        const auto [mp, mm] = computeBackwardsSubmatrices(coeff[i - 1], coeff[i], sigma, r_model);
+        const auto [mp, mm] = refractionMatrixBlocks(coeff[i - 1], coeff[i], sigma, r_model);
 
         const SpinMatrix delta = coeff[i - 1].computeDeltaMatrix(slices[i - 1].thicknessOr0());
 
@@ -225,7 +224,7 @@ SpinMatrix Compute::SpecularMagnetic::topLayerR(const SliceStack& slices,
             sigma = roughness->sigma();
 
         // compute the 2x2 submatrices in the write-up denoted as P+, P- and delta
-        const auto [mp, mm] = ::computeBackwardsSubmatrices(c_i, c_i1, sigma, r_model);
+        const auto [mp, mm] = ::refractionMatrixBlocks(c_i, c_i1, sigma, r_model);
 
         const SpinMatrix delta = c_i.computeDeltaMatrix(slices[i - 1].thicknessOr0());
 
diff --git a/Resample/Specular/ComputeFluxScalar.cpp b/Resample/Specular/ComputeFluxScalar.cpp
index 6c436bed53de06491d3bfd66016145f7f919b4ba..8ea84621d5c03292380a112b52b8cc602c08f8c0 100644
--- a/Resample/Specular/ComputeFluxScalar.cpp
+++ b/Resample/Specular/ComputeFluxScalar.cpp
@@ -25,6 +25,7 @@ using std::numbers::pi;
 
 namespace {
 
+//! See DevRef, chapter "Scattering by rough interfaces", section "Interface with tanh profile".
 std::pair<complex_t, complex_t> transition(complex_t kzi, complex_t kzi1, double sigma,
                                            const RoughnessModel& r_model)
 {