Finished discussion of S(a,b) and added description of probability tables in documentation.

docs/source/methods/index.rst

@@ -6,7 +6,7 @@ Theory and Methodology
 
 .. toctree::
     :numbered:
-    :maxdepth: 2
+    :maxdepth: 3
 
     introduction
     statistics

docs/source/methods/physics.rst

@@ -128,11 +128,17 @@ distribution function of the form :math:`p(x) = f_1(x) f_2(x)` with
 :math:`f_1(x)` bounded can be sampled by sampling :math:`x_s` from the
 distribution
 
-.. math:: \frac{f_2(x)}{\int f_2(x) dx}
+.. math::
+    :label: freegas7
+
+    \frac{f_2(x)}{\int f_2(x) dx}
 
 and accepting it with probability
 
-.. math:: \frac{f_1(x_s)}{\max f_1(x)}
+.. math::
+    :label: freegas8
+
+    \frac{f_1(x_s)}{\max f_1(x)}
 
 It is normally assumed that the velocity distribution of the target nucleus
 assumes a Maxwellian distribution in velocity.
@@ -181,6 +187,9 @@ scattering in hydrogenous solids such as polyethylene. As it occurs in ACE data,
 thermal inelastic scattering includes both coherent and incoherent effects and
 is dominant for most other materials including hydrogen in water.
 
+Calculating Integrated Cross Sections
+-------------------------------------
+
 The first aspect of using |sab| tables is calculating cross-sections to replace
 the data that would normally appear on the incident neutron data, which do not
 account for thermal binding effects. For incoherent elastic and inelastic
@@ -189,7 +198,7 @@ a specified energy grid. For coherent elastic data, the cross section can be
 expressed as
 
 .. math::
-    :label: bragg
+    :label: coherent-elastic-xs
 
     \sigma(E) = \frac{\sigma_c}{E} \sum_{E_i < E} f_i e^{-4WE_i}.
 
@@ -199,16 +208,170 @@ energies of the Bragg edges, and :math:`f_i` are related to crystallographic
 structure factors. Since the functional form of the cross-section is just 1/E
 and the proportionality constant changes only at Bragg edges, the
 proportionality constants are stored and then the cross-section can be
-calculated analytically based on equation :eq:`bragg`.
+calculated analytically based on equation :eq:`coherent-elastic-xs`.
+
+Outgoing Angle for Coherent Elastic Scattering
+----------------------------------------------
+
+The other aspect of using |sab| tables is determining the outgoing energy and
+angle of the neutron after scattering. For incoherent and coherent elastic
+scattering, the energy of the neutron does not actually change, but the angle
+does change. For coherent elastic scattering, the angle will depend on which
+Bragg edge scattered the neutron. The probability that edge :math:`i` will
+scatter then neutron is given by
+
+.. math::
+    :label: coherent-elastic-probability
+
+    \frac{f_i e^{-4WE_i}}{\sum_j f_j e^{-4WE_j}}.
+
+After a Bragg edge has been sampled, the cosine of the angle of scattering is
+given analytically by
+
+.. math::
+    :label: coherent-elastic-angle
+
+    \mu = 1 - \frac{E_i}{E}
+
+where :math:`E_i` is the energy of the Bragg edge that scattered the neutron. 
+
+Outgoing Angle for Incoherent Elastic Scattering
+------------------------------------------------
+
+For incoherent elastic scattering, the probability distribution for the cosine
+of the angle of scattering is represent as a series of equally-likely discrete
+cosines :math:`\mu_{i,j}` for each incoming energy :math:`E_i` on the thermal
+elastic energy grid. First the outgoing angle bin :math:`j` is sampled. Then, if
+the incoming energy of the neutron satisfies :math:`E_i < E < E_{i+1}` the final
+cosine is
+
+.. math::
+    :label: incoherent-elastic-angle
+
+    \mu = \mu_{i,j} + f (\mu_{i+1,j} - \mu_{i,j})
+
+where the interpolation factor is defined as
+
+.. math::
+    :label: sab-interpolation-factor
+
+    f = \frac{E - E_i}{E_{i+1} - E_i}.
+
+Outgoing Energy and Angle for Inelastic Scattering
+--------------------------------------------------
+
+On each |sab| table, there is a correlated angle-energy secondary distribution
+for neutron thermal inelastic scattering. While the documentation for the ACE
+format implies that there are a series of equiprobably outgoing energies, the
+outgoing energies may have non-uniform probability distribution. In particular,
+if the thermal data were processed with :math:`iwt = 0` in NJOY, then the first
+and last outgoing energies have a relative probability of 1, the second and
+second to last energies have a relative probability of 4, and all other energies
+have a relative probability of 10. The procedure to determine the outgoing
+energy and angle is as such. First, the inteprolation factor is determined from
+equation :eq:`sab-interpolation-factor`. Then, an outgoing energy bin is sampled
+either from a uniform distribution or from a skewed distribution as
+discussed. The outgoing energy is then interpolated between values corresponding
+to neighboring incoming energies:
+
+.. math::
+    :label: inelastic-energy
+
+    E = E_{i,j} + f (E_{i+1,j} - E_{i,j})
+
+where :math:`E_{i,j}` is the j-th outgoing energy corresponding to the i-th
+incoming energy. For each combination of incoming and outgoing energies, there
+is a series equiprobable outgoing cosines. An outgoing cosine bin is sampled
+uniformly and then the final cosine is interpolated on the incoming energy grid:
+
+.. math::
+    :label: inelastic-angle
+
+    \mu = \mu_{i,j,k} + f (\mu_{i+1,j,k} - \mu_{i,j,k})
+
+where :math:`\mu_{i,j,k}` is the k-th outgoing cosine corresponding to the j-th
+outgoing energy and the i-th incoming energy. 
 
 ----------------------------------------------
 Unresolved Resonance Region Probability Tables
 ----------------------------------------------
 
+In the unresolved resonance energy range, resonances may be so closely spaced
+that it is not possible for experimental measurements to resolve all
+resonances. To properly account for self-shielding in this energy range, OpenMC
+uses the probability table method [Levitt]_. For most thermal reactors, the use
+of probability tables will not significantly affect problem results. However,
+for some fast reactors and other problems with an appreciable flux spectrum in
+the unresolved resonance range, not using probability tables may lead to
+incorrect results.
+
+Probability tables in the ACE format are generated from the UNRESR module in
+NJOY following the method of Levitt. A similar method employed for the RACER and
+MC21_ Monte Carlo codes is described in [Sutton]_. For the discussion here, we
+will focus only on use of the probability table table as it appears in the ACE
+format.
+
+Each probability table for a nuclide contains the following information at a
+number of incoming energies within the unresolved resonance range:
+
+- Cumulative probabilities for cross section bands
+- Total cross section (or factor) in each band
+- Elastic scattering cross section (or factor) in each band
+- Fission cross section (or factor) in each band
+- :math:`(n,\gamma)` cross section (or factor) in each band
+- Neutron heating number (or factor) in each band
+
+It should be noted that unresolved resonance probability tables affect only
+integrated cross sections and no extra data need be given for secondary
+angle/energy distributions. Secondary distributions for elastic and inelastic
+scattering would be specified whether or not probability tables were present.
+
+The procedure for determining cross sections in the unresolved range using
+probability tables is as follows. First, the bounding incoming energies are
+determined, i.e. find :math:`i` such that :math:`E_i < E < E_{i+1}`. We then
+sample a cross section band :math:`j` using the cumulative probabilities for
+table :math:`i`. This allows us to then calculate the elastic, fission, and
+capture cross-sections from the probability tables interpolating between
+neighboring incoming energies. If interpolation is specified, then
+the cross sections are calculated as
+
+.. math::
+    :label: ptables-linlin
+
+    \sigma = \sigma_{i,j} + f (\sigma_{i+1,j} - \sigma{i,j})
+
+where :math:`f` is the interpolation factor defined in the same manner as
+:eq:`sab-interpolation-factor`. If logarithmic interpolation is specified, the
+cross sections are calculated as
+
+.. math::
+    :label: ptables-loglog
+
+    \sigma = \exp \left ( \log \sigma_{i,j} + f \log
+    \frac{\sigma_{i+1,j}}{\sigma_{i,j}} \right )
+
+where the interpolation factor is now defined as
+
+.. math::
+    :label: log-interpolation-factor
+
+    f = \frac{\log \frac{E}{E_i}}{\log \frac{E_{i+1}}{E_i}}
+
+A flag is also present in the probability table that specifies whether an
+inelastic cross section should be calculated. If so, this is done from a normal
+reaction cross section (either MT=51 or a special MT). Finally, if the
+cross-sections defined are above are specified to be factors and not true
+cross-sections, they are multiplied by the underlying smooth cross section in
+the unresolved range to get the actual cross sections. Lastly, the total cross
+section is calculated as the sum of the elastic, fission, capture, and inelastic
+cross sections.
+
 .. _NJOY: http://t2.lanl.gov/codes.shtml
 
 .. _PREPRO: http://www-nds.iaea.org/ndspub/endf/prepro/
 
+.. _MC21: http://www.osti.gov/bridge/servlets/purl/903083-HT5p1o/903083.pdf
+
 .. [SIGMA1] Dermett E. Cullen and Charles R. Weisbin, "Exact Doppler Broadening
    of Tabulated Cross Sections," *Nucl. Sci. Eng.*, **60**, pp. 199-229 (1976).
 
@@ -221,4 +384,13 @@ Unresolved Resonance Region Probability Tables
 .. [Squires] G. L. Squires, *Introduction to the Theory of Thermal Neutron
    Scattering*, Cambridge University Press (1978).
 
+.. [Levitt] Leo B. Levitt, "The Probability Table Method for Treating Unresolved
+   Neutron Resonances in Monte Carlo Calculations," *Nucl. Sci. Eng.*, **49**,
+   pp. 450-457 (1972).
+
+.. [Sutton] Thomas M. Sutton and Forrest B. Brown, "Implementation of
+   the Probability Table Method in a Continuous-Energy Monte Carlo Code System,"
+   *Proc. International Conf. on the Physics of Nucl. Sci. and Technology*,
+   October 5-8, Long Island, New York (1998).
+
 .. |sab| replace:: S(:math:`\alpha,\beta`)