diff --git a/Doc/DevRef/DevRef.pdf b/Doc/DevRef/DevRef.pdf
index 6e3be484bb4d95f557b758a473bedf16d349c7ad..eaf81aa82f6409347db63dee74c5ecf7612623ed 100644
Binary files a/Doc/DevRef/DevRef.pdf and b/Doc/DevRef/DevRef.pdf differ
diff --git a/Doc/DevRef/DevRef.tex b/Doc/DevRef/DevRef.tex
index ae2fe73404b0d778e5edbff720734458383973d6..783f24c2f4be0ad5cf5502ff5225d64114a4c047 100644
--- a/Doc/DevRef/DevRef.tex
+++ b/Doc/DevRef/DevRef.tex
@@ -26,6 +26,9 @@
\documentclass[a4paper,11pt,fleqn,draft]{report}\usepackage[final]{graphicx}
%\documentclass[a4paper,11pt,fleqn,draft]{report}\usepackage[draft]{graphicx}
+\usepackage{titlesec}
+\titleformat*{\section}{\bfseries\boldmath}
+
\def\shorttitle{BornAgain Developers Reference}
\def\longtitle{\shorttitle\\ Physics and Numerics}
diff --git a/Doc/DevRef/Macros.tex b/Doc/DevRef/Macros.tex
index 8080cf6c113b0eb0f37acc6cbd7021722c4fefbf..bdde67e5ff6f93766bdbc59be27a0184d6909bd8 100644
--- a/Doc/DevRef/Macros.tex
+++ b/Doc/DevRef/Macros.tex
@@ -23,6 +23,9 @@
% Tensors
\def\TENS#1{\ensuremath{\mathsf{#1}}}
\def\ONE{\TENS{1}}
+\def\U#1{\underline{#1}}
+\def\UU#1{\underline{\underline{#1}}}
+\def\WW#1{\mathbb{#1}}
% Math operators
\newcommand{\ggll}{\mathrel{\substack{\ll\\[-.05em]\gg}}}
@@ -62,7 +65,8 @@
\def\eps{\epsilon}
\def\etal{\E{et al.}}
\def\e{\mathrm e}
-\def\eB{\v{e}_{\v{B}}}
+\def\eB{\breve{\v{B}}}
+\def\eb{\breve{B}} % component of \eB
\def\evp{c} % eigenvalue of \hat p
\def\hx{\v{\hat x}}
\def\hy{\v{\hat y}}
diff --git a/Doc/DevRef/Polarized.tex b/Doc/DevRef/Polarized.tex
index 638cb933a08dbd1162393feaf8a339657b515304..12e471033b42835882af731d99ecd905cb3aa18e 100644
--- a/Doc/DevRef/Polarized.tex
+++ b/Doc/DevRef/Polarized.tex
@@ -24,12 +24,8 @@ to polarized neutrons.
We therefore need to study spinor wave equations,
in contrast to the scalar theory of the previous chapters.
-\def\U#1{\underline{#1}}
-\def\UU#1{\underline{\underline{#1}}}
-\def\WW#1{\mathbb{#1}}
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{Polarized neutron propagation}\label{Snpol0}
+\section{Polarized neutrons in 2+1 dimensions}\label{Snpol0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%===============================================================================
@@ -195,6 +191,51 @@ It does not exist if $\sldN$ is real and $\sldM= \kappa_0^2/(4\pi)-\sldN$.
If $\sldM$ is even larger, then $\opkappa$ becomes pure imaginary,
causing evanescent waves, to be discussed later~(\cref{Seva}).
+%==================================================================================================%
+\subsection{Eigenvalue decomposition}\label{Skeigen}
+%==================================================================================================%
+
+To evaluate functions of the operator~$\opkappa$, we will need its eigenvalue decomposition.
+We start with the matrix $\eB\Pauli$,
+which has the eigenvalues $\pm 1$ and the normalized eigenspinors
+\begin{equation}\label{Ev1v2}
+ V_+ = \frac{1}{\sqrt{2(1+\eb_z)}} \begin{pmatrix}
+ 1 +\eb_z \\ \eb_x+i\eb_y \end{pmatrix},\quad
+ V_- = \frac{1}{\sqrt{2(1+\eb_z)}} \begin{pmatrix}
+ \eb_x-i\eb_y\\ -1 -\eb_z \end{pmatrix}.
+\end{equation}
+For readability, we have omitted the subscript $\v{B}$ from the components of~$\eB$.
+and the same eigenvectors as $\eB\Pauli$.
+We introduce the eigenvector matrix
+\begin{equation}\label{D0QofB}
+ \hat{Q}(\v{B}) \coloneqq \left(V_+, V_-\right)
+ = \frac{1}{\sqrt{2(1+\eb_z)}}
+ \begin{pmatrix} 1 +\eb_z & \eb_x-i\eb_y \\ \eb_x+i\eb_y & -1 -\eb_z \end{pmatrix}.
+\end{equation}
+The normalization factor becomes singular for $\eb_z=-1$.
+In this case, the matrix $\eB\Pauli$ is just $\hat\sigma_z$
+and has eigenvectors $V_+=(0,1)^\dagger$ and $V_-=(1,0)^\dagger$.
+Furthermore, we need to take care of the case $\v{B}=0$.
+Altogether, we let
+\begin{equation}\label{DQofB}
+ \hat{Q}(\v{B}) \coloneqq \left\{\begin{array}{ll}
+ \hat{1} &\text{ if } B=0, \\
+ \hat\sigma_x &\text{ if } B_z=-B, \\
+ (\eB+\breve{\v{z}})\Pauli/\sqrt{2(1+\eb_z)} &\text{ else.}
+ \end{array}\right.
+\end{equation}
+The matrix $\opkappa$ has the eigenvalues $\evp_\pm$,
+and the same eigenvectors as~$\eB\Pauli$.
+Accordingly, it has the eigenvalue decomposition
+\begin{equation}
+ \opkappa
+ = \hat{Q} \, \begin{pmatrix}\evp_+ & 0 \\ 0 & \evp_-\end{pmatrix} \hat{Q}^\dagger,
+\end{equation}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Refraction and reflection at interfaces}\label{Spolif}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
%===============================================================================
\subsection{Transfer matrix}
%===============================================================================
@@ -252,30 +293,7 @@ in \SRC{Resample/Specular}{ComputeFluxMagnetic.cpp}.}
\subsection{Phase rotation matrix}\label{Sphase}
%==================================================================================================%
-The matrix $\eB\Pauli$ has the eigenvalues $\pm 1$ and the normalized eigenspinors
-\begin{equation}\label{Ev1v2}
- V_+ = \frac{1}{\sqrt{2(1+e_z)}} \begin{pmatrix}
- 1 +e_z \\ e_x+ie_y \end{pmatrix},\quad
- V_- = \frac{1}{\sqrt{2(1+e_z)}} \begin{pmatrix}
- e_x-ie_y\\ -1 -e_z \end{pmatrix}.
-\end{equation}
-For readability, we have omitted the subscript $\v{B}$ from the components of~$\eB$.
-The normalization factor becomes singular for $e_z=-1$.
-Therefore for $e_z<1+10^{-8}$, we simply take the known
-eigenvalues of $\hat\sigma_z$, namely $V_+=(0,1)^\dagger$ and $V_-=(1,0)^\dagger$.
-
-The matrix $\opkappa$ has the eigenvalues $\evp_\pm$,
-and the same eigenvectors as $\eB\Pauli$.
-We introduce the eigenvector matrix
-\begin{equation}\label{DQofB}
- \hat{Q} \coloneqq \left(V_+, V_-\right).
-\end{equation}
-Then $\opkappa$ has the eigenvalue decomposition
-\begin{equation}
- \opkappa
- = \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
+The phase rotation matrix \cref{Dopdel} can be written\footnote
{Currently (jun23) implemented in local function \T{PhaseRotationMatrix}
in file \SRC{Resample/Flux}{MatrixFlux.cpp}.}
\begin{equation}\label{EdP2}
diff --git a/Doc/DevRef/Title.tex b/Doc/DevRef/Title.tex
index 946b6dc3f1158eb07f3cac7c1b81cf04c4f027ee..706046194b3351e4c5187c824fbf3c1dd989a073 100644
--- a/Doc/DevRef/Title.tex
+++ b/Doc/DevRef/Title.tex
@@ -47,7 +47,14 @@ see the reference paper Pospelov et al 2020 \cite{PoHB20}
and the web docs at \url{https://www.bornagainproject.org}.
This reference provides some of the theory behind the code.
-It is principally targetted at fellow developers.
+\\[2ex]
+
+Notation: Bold symbols ($\r$, $\v{B}$) are real or complex 3-vectors.
+A hat ($\hat{v}$, $\hat\sigma_z$) denotes an operator in spin space,
+represented by a complex $2\times2$ matrix.
+Spinors are represented by uppercase letters ($\Psi$, $\Phi$, $T$, $R$).
+Black-board bold ($\WW{M}$) denotes $4\times4$ matrices.
+
~\vfill
\noindent
\begin{tabular}{@{}p{.25\textwidth}@{}p{.75\textwidth}@{}}