diff --git a/docs/tex/sammyrefs.bib b/docs/tex/sammyrefs.bib
index 8552a67091a40da45069d094910c1704272cfac4..491bd21593301f813169e6ba179a9e351130e8d9 100644
--- a/docs/tex/sammyrefs.bib
+++ b/docs/tex/sammyrefs.bib
@@ -1228,6 +1228,11 @@ abstract = "The response of the RPI 16-segment NaI(Tl) multiplicity detector sys
publisher={APS}
}
+@misc{frohner_lecture_notes_02,
+ author={Fr\"{o}hner, F.H.},
+ date = {2002},
+ howpublished = {{lecture notes on R-matrix theory, presented at the SAMMY workshop of February 18-22, 2002, at the NEA Headquarters in Paris, France (2002)}}
+}
@Comment{jabref-meta: databaseType:bibtex;}
diff --git a/docs/tex/scattering-theory.tex b/docs/tex/scattering-theory.tex
index 9af8110105abac34bb04f7ab1463219efc5d0198..6285c1548bf25c7275d940ca2740992d1b56e56a 100644
--- a/docs/tex/scattering-theory.tex
+++ b/docs/tex/scattering-theory.tex
@@ -228,7 +228,7 @@ The R-matrix is not the only possibility for parameterization of the scattering
\noindent
To see the relationship of the A-matrix to the R-matrix, we begin by multiplying both
-sides of Eq. (II A1.4) by $A$ and summing over $\lambda$:
+sides of Eq. (II A1.4)\ref{eq:a-matrix} by $A$ and summing over $\lambda$:
\begin{equation}
\begin{aligned}
@@ -241,7 +241,7 @@ sides of Eq. (II A1.4) by $A$ and summing over $\lambda$:
\noindent
Dividing by $(E_\mu - E)$, multiplying on the left by $\gamma_{\mu c'}$ and on the right by $\gamma_{\nu c''}$, and summing over $\mu$ puts this equation into the form
-\begin{equation}
+\begin{equation}\label{a-matrix-transform-1}
\begin{aligned}
\sum_\mu \gamma_{\mu c'}\left(E_{\mu}-E\right)^{-1}\delta_{\mu\nu}\gamma_{\nu c''} = & \sum_\mu \gamma_{\mu c'}\left(E_\mu-E\right)^{-1}\left(E_\mu-E\right)A_{\mu\nu}\gamma_{\nu c''} \\
& - \sum_\mu \gamma_{\mu c'} \left(E_{\mu}-E\right)^{-1} \sum_c \gamma_{\mu c} L_c \sum_\lambda \gamma_{\lambda c} A_{\lambda \nu}\gamma_{\nu c''} ,
@@ -251,7 +251,7 @@ Dividing by $(E_\mu - E)$, multiplying on the left by $\gamma_{\mu c'}$ and on t
\noindent
which can be reduced to
-\begin{equation}
+\begin{equation}\label{a-matrix-transform-2}
\begin{aligned}
\gamma_{\nu c'}\left(E_{\nu}-E\right)^{-1}\gamma_{\nu c''} = & \sum_{\mu} \gamma_{\mu c'} A_{\mu\nu}\gamma_{\nu c''} \\
& - \sum_c \left[ \sum_\mu \gamma_{\mu c'}\left(E_{\mu}-E\right)^{-1}\gamma_{\mu c} \right] L_c \sum_\lambda \gamma_{\lambda c} A_{\lambda\nu} \gamma_{\nu c''}.
@@ -281,11 +281,212 @@ in which we can replace the quantities in square brackets by the R-matrix, givin
\noindent
Solving for the summation, this equation can be rewritten as
-\begin{equation}
+\begin{equation}\label{eq:r-matrix-to-a-matrix}
\left[(I-RL)^{-1}R\right]_{cc''} = \sum_{\lambda\nu} \gamma_{\lambda c} A_{\lambda\nu} \gamma_{\nu c''}.
\end{equation}
+\noindent
+To relate this to the scattering matrix, we note that Eq. (II A.6)\ref{eq:W-matrix} can be rewritten using Eq. (II A.7)\ref{eq:L-matrix} into the form
+
+\begin{equation}\label{eq:w-matrix-transform}
+ \begin{aligned}
+ W & = P^{1/2}\left(I-RL\right)^{-1}\left(I-RL^*\right)P^{-1/2} \\
+ & = P^{1/2}\left(I-RL\right)^{-1}\left(I-RL+2iRP\right)P^{-1/2} \\
+ & = P^{1/2}\left[\left(I-RL\right)^{-1}\left(I-RL\right) + 2i\left(I-RL\right)^{-1}RP\right]P^{-1/2} \\
+ & = P^{1/2}P^{-1/2} + 2iP^{1/2}\left(I-RL\right)^{-1}RP P^{-1/2} \\
+ & = I + 2iP^{1/2}\left(I-RL\right)^{-1}RP^{1/2} .
+ \end{aligned}
+\end{equation}
+
+\noindent
+Comparing Eq. (II A1.10)\ref{eq:r-matrix-to-a-matrix} to Eq. (II A1.11)\ref{eq:w-matrix-transform} gives, in matrix form,
+
+\begin{equation}
+ W = I + 2iP^{1/2}\gamma A\gamma P^{1/2} .
+\end{equation}
+
+\noindent
+These equations are exact; no approximations have been made.
+
+One common approximation should be discussed here: the ``eliminated channel'' approximation, for which one particular type of channel is treated in aggregate and assumed to not interfere from level to level. This is most easily understood in the A-matrix definition, Eq. (II A1.4)\ref{eq:a-matrix}; assuming no level-level interference for the gamma channels (for example), this equation can be approximated as
+
+\begin{equation}\label{eq:a-matrix-rm-approx}
+ A_{\mu\lambda}^{-1} \approx (E_\lambda-E)\delta_{\mu\lambda} - \left[\sum_{\substack{\gamma=gamma \\ channels}} \gamma_{\mu\gamma}L_{\gamma}\gamma_{\lambda\gamma}\right]\delta_{\mu\lambda} - \sum_{\substack{c=particle \\channels}} \gamma_{\mu c}L_c\gamma_{\lambda c} .
+\end{equation}
+
+\noindent
+The quantity in square brackets corresponds to those channels for which the level-level
+interference is to be neglected; that is, only the interactions within one level are important. For
+gamma channels, $L=S+iP$ reduces to $L=i$, so Eq. (II A1.13)\ref{eq:a-matrix-rm-approx} becomes
+
+\begin{equation}\label{eq:a-matrix-rm-reduced}
+ A_{\mu\lambda}^{-1} \approx \left(E_\lambda-E-i\bar{\Gamma}_{\lambda\gamma}/2\right)\delta_{\mu\lambda} - \sum_{\substack{c=particle \\ channels}} \gamma_{\mu c}L_c\gamma_{\lambda c} .
+\end{equation}
+
+\noindent
+The bar over $\bar{\Gamma}_{\lambda\gamma}$ is used to indicate the special treatment for this channel.
+
+In this form, our expression for $A$ is analogous to the exact expression in Eq. (II A1.4)\ref{a-matrix} with two modifications: the additional imaginary term is added to the energy difference, and the sum over the channels includes only the ``particle channels'' (non-eliminated channels). It is therefore possible to immediately write the R-matrix formula for the eliminated-channel approximation as
+
+\begin{equation}
+ R_{cc'} = \sum_{\lambda}\frac{\gamma_{\lambda c}\gamma_{\lambda c'}}{E_\lambda-E-i\bar{\Gamma}_{\lambda\gamma}/2} \delta_{JJ'} ,
+\end{equation}
+
+where the channel indices c and c ' refer only to particle channels, not to the gamma channels.
+This formula for the R-matrix is the Reich-Moore approximation and is the form which is used in
+the SAMMY code. See Section II.B.1 for more about this formulation of R-matrix theory.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Derivation of Scattering Theory Equations}\label{subsec:derivation-of-scat-theory}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Many authors have given derivations of the equations for the scattering matrix in terms of the R-matrix. Sources for the derivation shown here are unpublished lecture notes of Fr\"{o}hner [FF02]\cite{frohner_lecture_notes_02}, presented at the SAMMY workshop in Paris in 2002, and Foderaro [AF71]\cite{foderaro_1971}. This derivation is valid for only the simple case of spinless projectiles and target nuclei, assuming only elastic scattering and absorption. For the general case, the reader is referred to Lane and Thomas [AL58]\cite{lane_thomas_1958}.
+
+\noindent
+\textbf{Schr\"{o}dinger equation}
+
+The Schr\"{o}dinger equation with a complex potential is
+
+\begin{equation}
+ \left( \frac{-\hbar^2}{2m}\nabla^2 + V + iW \right)\psi = E\psi,
+\end{equation}
+
+\noindent
+in which one can consider that $V$ causes scattering and $W$ causes absorption. The wave function
+can be expanded in the usual fashion,
+
+\begin{equation}
+ \psi\left(r,\text{cos}(\theta)\right) = \sum_{l=0}^\infty \frac{u_l(r)}{r} P_l(\text{cos}(\theta)) ,
+\end{equation}
+
+\noindent
+for which the radial portion obeys the equation
+
+\begin{equation}
+ \frac{d^2u_l}{dr^2} + \left[ k^2-\frac{2m}{\hbar^2}(V+iW)-\frac{l(l+1)}{r^2} \right]u_l = 0,
+\end{equation}
+
+\noindent
+subject to the conditions that $|\psi|^2$ is everywhere finite and that
+
+\begin{equation}
+ u_l(r=0) = 0.
+\end{equation}
+
+In the external region, $r>a$, the nuclear forces are zero $(V=W=0)$, so the solution has
+the form
+
+\begin{equation}
+ u_l(r) = I_l(r) - U_l O_l(r) .
+\end{equation}
+
+$I_l$ represents an incoming free wave, and $O_l$ represents an outgoing free wave. $U_l$ is the ``collision function'' or ``S function'' that describes the effects of the nuclear interaction, giving both the attenuation and the phase shift of the outgoing wave:
+
+\begin{equation}
+ \begin{aligned}
+ |U_l|^2 = 1 & \text{ for } W = 0, \\
+ |U_l|^2 < 1 & \text{ for } W \neq 0.
+ \end{aligned}
+\end{equation}
+
+\noindent
+Our goal is to determine an appropriate analytic form for $U_l$.
+
+\noindent
+\textbf{Orthogonal eigenvectors in interior region}
+
+For the interior region $r<a$ , we define eigenfunctions $w_{\lambda l}(r)$ and eigenvalues $E_\lambda$,
+
+\begin{equation}
+ E_\lambda = \frac{\hbar^2k_\lambda^2}{2m},
+\end{equation}
+
+\noindent
+for the wave equation without absorption $(W = 0)$,
+
+\begin{equation}\label{eq:wave-func-without-abs}
+ \frac{d^2w_{\lambda l}}{dr^2} + \left[k_\lambda^2-\frac{2m}{\hbar^2}V-\frac{l(l+1)}{r^2}\right]w_{\lambda l} = 0,
+\end{equation}
+
+\noindent
+for which the boundary conditions are
+
+\begin{equation}\label{eq:schrod-bound-conditions}
+ w_{\lambda l}(r=0) = 0 \;\;\;\;\; \text{and} \;\;\;\;\; \frac{a}{w_{\lambda l}(a)}\frac{dw_{\lambda l}}{dr}\Bigr|_{r=a} = B_l.
+\end{equation}
+
+\noindent
+Note that $w_{\lambda l}(r)$ is real if the boundary parameter $B_l$ is chosen to be real. The eigenfunctions are orthogonal, since
+
+\begin{equation}\label{eq:eigen-func-orthog}
+ \begin{aligned}
+ \int_0^a\left( \frac{d^2w_{\lambda l}}{dr^2}w_{\mu l} - w_{\lambda l}\frac{d^2w_{\mu l}}{dr^2} \right) dr & = \int_0^a \frac{d}{dr} \left( \frac{dw_{\lambda l}}{dr}w_{\mu l} - w_{\lambda l}\frac{dw_{\mu l}}{dr} \right) dr \\
+ & = \left[\frac{dw_{\lambda l}}{dr}w_{\mu l} - w_{\lambda l}\frac{dw_{\mu l}}{dr}\right]_0^a \\
+ & = \frac{dw_{\lambda l}}{dr}\Bigr|_{r=a} w_{\mu l}(a) - w_{\lambda l}(a)\frac{dw_{\mu l}}{dr}\Bigr|_{r=a} - [0] \\
+ & = \frac{B_l}{a}\left[ w_{\lambda l}(a)w_{\mu l}(a)-w_{\lambda l}(a)w_{\mu l}(a) \right] = 0,
+ \end{aligned}
+\end{equation}
+
+\noindent
+in which both equations of (II A2.9)\ref{eq:schrod-bound-conditions} have been invoked. The integral in Eq. (II A2.10)\ref{eq:eigen-func-orthog} can also be evaluated using Eq. (II A2.8)\ref{eq:wave-func-without-abs}, giving
+
+\begin{equation}\label{eq:reform-orthog}
+ \begin{aligned}
+ \int_0^a\left( \frac{d^2w_{\lambda l}}{dr^2}w_{\mu l} - w_{\lambda l}\frac{d^2w_{\mu l}}{dr^2} \right) dr & = \int_0^a\left( \left[-k_\lambda^2-\frac{2mV}{\hbar^2}\right]w_{\lambda l}w_{\mu l} - w_{\lambda l}\left[-k_\lambda^2-\frac{2mV}{\hbar^2}\right]w_{\lambda l} \right)dr \\
+ & = \int_0^a\left( -k_\lambda^2w_{\lambda l}w_{\mu l} + k_\mu^2w_{\lambda l}w_{\mu l} \right)dr \\
+ & = -(k_\lambda^2-k_\mu^2)\int_0^a w_{\lambda l}w_{\mu l}dr .
+ \end{aligned}
+\end{equation}
+
+\noindent
+Equating Eq. (II A2.10)\ref{eq:eigen-func-orthog} to Eq. (II A2.11)\ref{eq:reform-orthog} gives
+
+\begin{equation}
+ (k_\lambda^2-k_\mu^2)\int_0^a w_{\lambda l}w_{\mu l}dr = 0.
+\end{equation}
+
+\noindent
+For $\lambda \neq \mu$, assuming no degenerate states, it therefore follows that
+
+\begin{equation}
+ \int_0^a w_{\lambda l}w_{\mu l}dr = 0 \;\;\;\;\; \text{if} \;\;\;\;\; \lambda\neq\mu.
+\end{equation}
+
+The orthogonality of the eigenvectors is therefore established. We assume that these wave
+functions are normalized such that
+
+\begin{equation}\label{eq:normalized-wave-func}
+ \int_0^a w_{\lambda l}w_{\mu l}dr = \delta_{\lambda\mu}.
+\end{equation}
+
+\noindent
+\textbf{Matching at the surface}
+
+The internal wave function for the true potential (including the imaginary part $iW$ ) can be
+expanded in terms of the eigenfunctions as
+
+\begin{equation}\label{eq:internal-wave-func}
+ u_l(r) = \sum_{\lambda} c_{\lambda l}w_{\lambda l}(r) \;\;\;\;\; \text{for} \;\; r\leq a,
+\end{equation}
+
+\noindent
+with
+
+\begin{equation}
+ c_{\lambda l} = \int_0^a u_l w_{\lambda l} dr.
+\end{equation}
+
+\noindent
+This equation for $c_{\lambda l}$ is derived by multiplying Eq. (II A2.15)\ref{eq:internal-wave-func} by $u_{\lambda l}(r)$ , integrating, and applying Eq. (II A2.14)\ref{eq:normalized-wave-func}.
+
+\noindent
+Consider now the integral
+
+\begin{equation}
+ \int_0^a\left( \frac{d^2w_{\lambda l}}{dr^2}w_{\mu l} - w_{\lambda l}\frac{d^2w_{\mu l}}{dr^2} \right) dr,
+\end{equation}
+which can be expanded by use of Eqs. (II A2.3) and (II A2.8) to give