Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/generate.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -258,5 +258,5 @@
 %(parameters)s
 
-The returned value is scaled to units of |cm^-1|.
+%(returns)s
 
 %(docs)s
@@ -609,9 +609,13 @@
     Return the documentation for the model.
     """
+    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale."
+    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$."
     info = make_info(kernel_module)
+    is_Sq = ("structure-factor" in info['category'])
     subst = dict(id=info['id'].replace('_', '-'),
                  name=info['name'],
                  title=info['title'],
                  parameters=make_partable(info['parameters']),
+                 returns=Sq_units if is_Sq else Iq_units,
                  docs=kernel_module.__doc__)
     return DOC_HEADER % subst
Index: sasmodels/models/HayterMSAsq.py
===================================================================
--- sasmodels/models/HayterMSAsq.py	(revision cd3dba0c60af7479ce1bae901d065d3d5786c8ea)
+++ sasmodels/models/HayterMSAsq.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -1,38 +1,33 @@
 # Note: model title and parameter table are inserted automatically
 r"""
-This calculates the structure factor (the Fourier transform of the pair correlation function *g(r)*) for a system of
-charged, spheroidal objects in a dielectric medium. When combined with an appropriate form factor (such as sphere,
-core+shell, ellipsoid, etc), this allows for inclusion of the interparticle interference effects due to screened coulomb
-repulsion between charged particles.
+This calculates the structure factor (the Fourier transform of the pair
+correlation function $g(r)$) for a system of charged, spheroidal objects
+in a dielectric medium. When combined with an appropriate form factor
+(such as sphere, core+shell, ellipsoid, etc), this allows for inclusion
+of the interparticle interference effects due to screened coulomb repulsion
+between charged particles.
 
-**This routine only works for charged particles**. If the charge is set to zero the routine will self-destruct!
-For non-charged particles use a hard sphere potential.
+**This routine only works for charged particles**. If the charge is set to
+zero the routine will self-destruct! For non-charged particles use a hard
+sphere potential.
 
-The salt concentration is used to compute the ionic strength of the solution which in turn is used to compute the Debye
-screening length. At present there is no provision for entering the ionic strength directly nor for use of any
-multivalent salts. The counterions are also assumed to be monovalent.
+The salt concentration is used to compute the ionic strength of the solution
+which in turn is used to compute the Debye screening length. At present
+there is no provision for entering the ionic strength directly nor for use
+of any multivalent salts. The counterions are also assumed to be monovalent.
 
-For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as
+For 2D data, the scattering intensity is calculated in the same way as 1D,
+where the $q$ vector is defined as
 
 .. math::
 
-    Q = \sqrt{Q_x^2 + Q_y^2}
+    q = \sqrt{q_x^2 + q_y^2}
 
-==============  ========  =============
-Parameter name  Units     Default value
-==============  ========  =============
-effect_radius   |Ang|     20.8
-charge          *e*       19
-volfraction     None      0.2
-temperature     K         318
-salt conc       M         0
-dielectconst    None      71.1
-==============  ========  =============
+.. figure:: img/HayterMSAsq_227.jpg
 
-.. image:: img/HayterMSAsq_227.jpg
+    1D plot using the default values (in linear scale).
 
-*Figure. 1D plot using the default values (in linear scale).*
-
-REFERENCE
+References
+----------
 
 J B Hayter and J Penfold, *Molecular Physics*, 42 (1981) 109-118
Index: sasmodels/models/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision 5ef0633b151126a2a8e376666e05ccbcc0341316)
+++ sasmodels/models/barbell.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -1,59 +1,47 @@
 #barbell model
-# cylinder model
 # Note: model title and parameter table are inserted automatically
 r"""
-
-Calculates the scattering from a barbell-shaped cylinder (This model simply
-becomes the DumBellModel when the length of the cylinder, *L*, is set to zero).
-That is, a sphereocylinder with spherical end caps that have a radius larger
-than that of the cylinder and the center of the end cap radius lies outside
-of the cylinder. All dimensions of the BarBell are considered to be
-monodisperse. See the diagram for the details of the geometry and restrictions
-on parameter values.
+Calculates the scattering from a barbell-shaped cylinder.  Like
+:ref:`capped-cylinder`, this is a sphereocylinder with spherical end
+caps that have a radius larger than that of the cylinder, but with the center
+of the end cap radius lying outside of the cylinder. See the diagram for
+the details of the geometry and restrictions on parameter values.
 
 Definition
 ----------
 
-The returned value is scaled to units of |cm^-1|\ |sr^-1|, absolute scale.
+.. figure:: img/barbell_geometry.jpg
 
-The barbell geometry is defined as
+    Barbell geometry, where $r$ is *radius*, $R$ is *bell_radius* and
+    $L$ is *length*. Since the end cap radius $R \geq r$ and by definition
+    for this geometry $h < 0$, $h$ is then defined by $r$ and $R$ as
+    $h = - \sqrt{R^2 - r^2}$
 
-.. image:: img/barbell_geometry.jpg
-
-where *r* is the radius of the cylinder. All other parameters are as defined
-in the diagram.
-
-Since the end cap radius
-*R* >= *r* and by definition for this geometry *h* < 0, *h* is then
-defined by *r* and *R* as
-
-*h* = -1 \* sqrt(*R*\ :sup:`2` - *r*\ :sup:`2`)
-
-The scattered intensity *I(q)* is calculated as
+The scattered intensity $I(q)$ is calculated as
 
 .. math::
 
-    I(Q) = \frac{(\Delta \rho)^2}{V} \left< A^2(Q)\right>
+    I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q)\right>
 
-where the amplitude *A(q)* is given as
+where the amplitude $A(q)$ is given as
 
 .. math::
 
-    A(Q) =&\ \pi r^2L
-        {\sin\left(\tfrac12 QL\cos\theta\right)
-            \over \tfrac12 QL\cos\theta}
-        {2 J_1(Qr\sin\theta) \over Qr\sin\theta} \\
+    A(q) =&\ \pi r^2L
+        \frac{\sin\left(\tfrac12 qL\cos\theta\right)}
+             {\tfrac12 qL\cos\theta}
+        \frac{2 J_1(qr\sin\theta)}{qr\sin\theta} \\
         &\ + 4 \pi R^3 \int_{-h/R}^1 dt
-        \cos\left[ Q\cos\theta
+        \cos\left[ q\cos\theta
             \left(Rt + h + {\tfrac12} L\right)\right]
         \times (1-t^2)
-        {J_1\left[QR\sin\theta \left(1-t^2\right)^{1/2}\right]
-             \over QR\sin\theta \left(1-t^2\right)^{1/2}}
+        \frac{J_1\left[qR\sin\theta \left(1-t^2\right)^{1/2}\right]}
+             {qR\sin\theta \left(1-t^2\right)^{1/2}}
 
-The < > brackets denote an average of the structure over all orientations.
-<*A* :sup:`2`\ *(q)*> is then the form factor, *P(q)*. The scale factor is
-equivalent to the volume fraction of cylinders, each of volume, *V*. Contrast
-is the difference of scattering length densities of the cylinder and the
-surrounding solvent.
+The $\left<\ldots\right>$ brackets denote an average of the structure over
+all orientations. $\left<A^2(q)\right>$ is then the form factor, $P(q)$.
+The scale factor is equivalent to the volume fraction of cylinders, each of
+volume, $V$. Contrast $\Delta\rho$ is the difference of scattering length
+densities of the cylinder and the surrounding solvent.
 
 The volume of the barbell is
@@ -64,5 +52,5 @@
 
 
-and its radius-of-gyration is
+and its radius of gyration is
 
 .. math::
@@ -76,38 +64,34 @@
         \left( 4R^3 6R^2h - 2h^3 + 3r^2L \right)^{-1}
 
-**The requirement that** *R* >= *r* **is not enforced in the model!** It is
-up to you to restrict this during analysis.
+.. note::
+    The requirement that $R \geq r$ is not enforced in the model! It is
+    up to you to restrict this during analysis.
 
-This example dataset is produced by running the Macro PlotBarbell(),
-using 200 data points, *qmin* = 0.001 |Ang^-1|, *qmax* = 0.7 |Ang^-1|,
-*sld* = 4e-6 |Ang^-2| and the default model values.
+.. figure:: img/barbell_1d.jpg
 
-.. image:: img/barbell_1d.jpg
+    1D plot using the default values (w/256 data point).
 
-*Figure. 1D plot using the default values (w/256 data point).*
+For 2D data, the scattering intensity is calculated similar to the 2D
+cylinder model.
 
-For 2D data: The 2D scattering intensity is calculated similar to the 2D
-cylinder model. For example, for |theta| = 45 deg and |phi| = 0 deg with
-default values for other parameters
+.. figure:: img/barbell_2d.jpg
 
-.. image:: img/barbell_2d.jpg
+    2D plot (w/(256X265) data points) for $\theta = 45^\circ$ and
+    $\phi = 0^\circ$ with default values for the remaining parameters.
 
-*Figure. 2D plot (w/(256X265) data points).*
+.. figure:: img/orientation.jpg
 
-.. image:: img/orientation.jpg
+    Definition of the angles for oriented 2D barbells.
 
-Figure. Definition of the angles for oriented 2D barbells.
+.. figure:: img/orientation2.jpg
 
-.. image:: img/orientation2.jpg
+    Examples of the angles for oriented pp against the detector plane.
 
-*Figure. Examples of the angles for oriented pp against the detector plane.*
-
-REFERENCE
----------
+References
+----------
 
 H Kaya, *J. Appl. Cryst.*, 37 (2004) 37 223-230
 
 H Kaya and N R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda and errata)
-
 """
 from numpy import inf
@@ -124,5 +108,5 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "", "Barbell scattering length density"],
+parameters = [["sld", "4e-6/Ang^2", 4, [-inf, inf], "", "Barbell scattering length density"],
               ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", "Solvent scattering length density"],
               ["bell_radius", "Ang", 40, [0, inf], "volume", "Spherical bell radius"],
Index: sasmodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/models/bcc.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -13,16 +13,16 @@
 The scattering intensity $I(q)$ is calculated as
 
-.. math:
+.. math::
 
-    I(q) = \frac{\text{scale}}{V_P} V_\text{lattice} P(q) Z(q)
+    I(q) = \frac{\text{scale}}{V_p} V_\text{lattice} P(q) Z(q)
 
 
-where *scale* is the volume fraction of spheres, *Vp* is the volume of the
-primary particle, *V(lattice)* is a volume correction for the crystal
+where *scale* is the volume fraction of spheres, $V_p$ is the volume of the
+primary particle, $V_\text{lattice}$ is a volume correction for the crystal
 structure, $P(q)$ is the form factor of the sphere (normalized), and $Z(q)$
 is the paracrystalline structure factor for a body-centered cubic structure.
 
 Equation (1) of the 1990 reference is used to calculate $Z(q)$, using
-equations (29)-(31) from the 1987 paper for *Z1*\ , *Z2*\ , and *Z3*\ .
+equations (29)-(31) from the 1987 paper for $Z1$, $Z2$, and $Z3$.
 
 The lattice correction (the occupied volume of the lattice) for a
@@ -30,5 +30,5 @@
 separation $D$ is
 
-.. math:
+.. math::
 
     V_\text{lattice} = \frac{16\pi}{3} \frac{R^3}{\left(D\sqrt{2}\right)^3}
@@ -38,5 +38,5 @@
 in the calculation of $Z(q)$
 
-.. math:
+.. math::
 
     \Delta a = g D
@@ -51,5 +51,5 @@
 For a crystal, diffraction peaks appear at reduced q-values given by
 
-.. math:
+.. math::
 
     \frac{qD}{2\pi} = \sqrt{h^2 + k^2 + l^2}
@@ -60,10 +60,10 @@
 correspond to (just the first 5)
 
-.. math:
+.. math::
 
-    \begin{eqnarray}
-    &q/q_o&&\quad 1&& \ \sqrt{2} && \ \sqrt{3} && \ \sqrt{4} && \ \sqrt{5} \\
-    &\text{Indices}&& (110) && (200) && (211) && (220) && (310)
-    \end{eqnarray}
+    \begin{array}{lccccc}
+    q/q_o          &   1   & \sqrt{2} & \sqrt{3} & \sqrt{4} & \sqrt{5} \\
+    \text{Indices} & (110) &    (200) & (211)    & (220)    & (310)    \\
+    \end{array}
 
 **NB**: The calculation of $Z(q)$ is a double numerical integral that must
@@ -86,5 +86,5 @@
 model computation.
 
-.. figure:: img/crystal_orientation.gif
+.. figure:: img/crystal_orientation.png
 
     Orientation of the crystal with respect to the scattering plane.
@@ -94,6 +94,6 @@
     2D plot using the default values (w/200X200 pixels).*
 
-Reference
----------
+References
+----------
 
 Hideki Matsuoka et. al. *Physical Review B*, 36 (1987) 1754-1765
Index: sasmodels/models/broad_peak.py
===================================================================
--- sasmodels/models/broad_peak.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/models/broad_peak.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -10,21 +10,19 @@
 spherical morphologies, or for bicontinuous structures).
 
-The returned value is scaled to units of |cm^-1|, absolute scale.
-
 Definition
 ----------
 
-The scattering intensity *I(q)* is calculated as
+The scattering intensity $I(q)$ is calculated as
 
-.. math:
+.. math::
 
-    I(q) = \frac{A}{Q^n} + \frac{C}{1 + (Q\xi}^m} + B
+    I(q) = \frac{A}{q^n} + \frac{C}{1 + (q\xi)^m} + B
 
-Here the peak position is related to the d-spacing as *Q0* = 2|pi| / *d0*.
+Here the peak position is related to the d-spacing as $q_o = 2\pi / d_o$.
 
-For 2D data: The 2D scattering intensity is calculated in the same way as 1D,
-where the *q* vector is defined as
+For 2D data the scattering intensity is calculated in the same way as 1D,
+where the $q$ vector is defined as
 
-.. math:
+.. math::
 
     q = \sqrt{q_x^2 + q_y^2}
@@ -35,6 +33,6 @@
     1D plot using the default values (w/200 data point).
 
-REFERENCE
----------
+References
+----------
 
 None.
@@ -75,9 +73,9 @@
              / (1.0 + (abs(q - peak_pos) * lorentz_length) ** lorentz_exp))
     return inten
-Iq.vectorized = True  # Iq accepts an array of Q values
+Iq.vectorized = True  # Iq accepts an array of q values
 
 def Iqxy(qx, qy, *args):
     return Iq(sqrt(qx ** 2 + qy ** 2), *args)
-Iqxy.vectorized = True # Iqxy accepts an array of Qx, Qy values
+Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
 
 
Index: sasmodels/models/capped_cylinder.py
===================================================================
--- sasmodels/models/capped_cylinder.py	(revision 485aee225e23a0dad14b9acfc70c0c7072e650b9)
+++ sasmodels/models/capped_cylinder.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -1,53 +1,46 @@
 r"""
 Calculates the scattering from a cylinder with spherical section end-caps.
-This model simply becomes the a convex lens when the length of the cylinder
-$L=0$, that is, a sphereocylinder with end caps that have a radius larger
-than that of the cylinder and the center of the end cap radius lies within
-the cylinder. See the diagram for the details of the geometry and
-restrictions on parameter values.
+Like :ref:`barbell`, this is a sphereocylinder with end caps that have a
+radius larger than that of the cylinder, but with the center of the end cap
+radius lying within the cylinder. This model simply becomes the a convex
+lens when the length of the cylinder $L=0$. See the diagram for the details
+of the geometry and restrictions on parameter values.
 
 Definitions
 -----------
 
-The returned value is scaled to units of |cm^-1|\ |sr^-1|, absolute scale.
+.. figure:: img/capped_cylinder_geometry.jpg
 
-The capped cylinder geometry is defined as
+    Capped cylinder geometry, where $r$ is *radius*, $R$ is *bell_radius* and
+    $L$ is *length*. Since the end cap radius $R \geq r$ and by definition
+    for this geometry $h < 0$, $h$ is then defined by $r$ and $R$ as
+    $h = - \sqrt{R^2 - r^2}$
 
-.. image:: img/capped_cylinder_geometry.jpg
-
-where $r$ is the radius of the cylinder. All other parameters are as defined
-in the diagram. Since the end cap radius $R \ge r$ and by definition for this
-geometry $h < 0$, $h$ is then defined by $r$ and $R$ as
+The scattered intensity $I(q)$ is calculated as
 
 .. math::
 
-    h = - \sqrt{R^2 - r^2}
+    I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q)\right>
 
-The scattered intensity $I(Q)$ is calculated as
+where the amplitude $A(q)$ is given as
 
 .. math::
 
-    I(Q) = \frac{(\Delta \rho)^2}{V} \left< A^2(Q)\right>
-
-where the amplitude $A(Q)$ is given as
-
-.. math::
-
-    A(Q) =&\ \pi r^2L
-        {\sin\left(\tfrac12 QL\cos\theta\right)
-            \over \tfrac12 QL\cos\theta}
-        {2 J_1(Qr\sin\theta) \over Qr\sin\theta} \\
+    A(q) =&\ \pi r^2L
+        \frac{\sin\left(\tfrac12 qL\cos\theta\right)}
+            {\tfrac12 qL\cos\theta}
+        \frac{2 J_1(qr\sin\theta)}{qr\sin\theta} \\
         &\ + 4 \pi R^3 \int_{-h/R}^1 dt
-        \cos\left[ Q\cos\theta
+        \cos\left[ q\cos\theta
             \left(Rt + h + {\tfrac12} L\right)\right]
         \times (1-t^2)
-        {J_1\left[QR\sin\theta \left(1-t^2\right)^{1/2}\right]
-             \over QR\sin\theta \left(1-t^2\right)^{1/2}}
+        \frac{J_1\left[qR\sin\theta \left(1-t^2\right)^{1/2}\right]}
+             {qR\sin\theta \left(1-t^2\right)^{1/2}}
 
-The $\left< \ldots \right>$ brackets denote an average of the structure over
-all orientations. $\left< A^2(Q)\right>$ is then the form factor, $P(Q)$.
+The $\left<\ldots\right>$ brackets denote an average of the structure over
+all orientations. $\left< A^2(q)\right>$ is then the form factor, $P(q)$.
 The scale factor is equivalent to the volume fraction of cylinders, each of
-volume, $V$. Contrast is the difference of scattering length densities of
-the cylinder and the surrounding solvent.
+volume, $V$. Contrast $\Delta\rho$ is the difference of scattering length
+densities of the cylinder and the surrounding solvent.
 
 The volume of the capped cylinder is (with $h$ as a positive value here)
@@ -58,5 +51,5 @@
 
 
-and its radius-of-gyration is
+and its radius of gyration is
 
 .. math::
@@ -73,5 +66,5 @@
 .. note::
 
-    The requirement that $R \ge r$ is not enforced in the model!
+    The requirement that $R \geq r$ is not enforced in the model!
     It is up to you to restrict this during analysis.
 
@@ -104,5 +97,6 @@
     Examples of the angles for oriented pp against the detector plane.
 
-REFERENCE
+References
+----------
 
 H Kaya, *J. Appl. Cryst.*, 37 (2004) 223-230
Index: sasmodels/models/core_shell_cylinder.py
===================================================================
--- sasmodels/models/core_shell_cylinder.py	(revision 485aee225e23a0dad14b9acfc70c0c7072e650b9)
+++ sasmodels/models/core_shell_cylinder.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -12,5 +12,5 @@
 .. math::
 
-    P(Q,\alpha) = {\text{scale} \over V_s} F^2(Q) + \text{background}
+    P(q,\alpha) = \frac{\text{scale}}{V_s} F^2(q) + \text{background}
 
 where
@@ -18,14 +18,14 @@
 .. math::
 
-    F(Q) = &\ (\rho_c - \rho_s) V_c
-           {\sin \left( Q \tfrac12 L\cos\alpha \right)
-               \over Q \tfrac12 L\cos\alpha }
-           {2 J_1 \left( QR\sin\alpha \right)
-               \over QR\sin\alpha } \\
+    F(q) = &\ (\rho_c - \rho_s) V_c
+           \frac{\sin \left( q \tfrac12 L\cos\alpha \right)}
+                {q \tfrac12 L\cos\alpha}
+           \frac{2 J_1 \left( qR\sin\alpha \right)}
+                {qR\sin\alpha} \\
          &\ + (\rho_s - \rho_\text{solv}) V_s
-           {\sin \left( Q \left(\tfrac12 L+T\right) \cos\alpha \right)
-               \over Q \left(\tfrac12 L +T \right) \cos\alpha }
-           { 2 J_1 \left( Q(R+T)\sin\alpha \right)
-               \over Q(R+T)\sin\alpha }
+           \frac{\sin \left( q \left(\tfrac12 L+T\right) \cos\alpha \right)}
+                {q \left(\tfrac12 L +T \right) \cos\alpha}
+           \frac{ 2 J_1 \left( q(R+T)\sin\alpha \right)}
+                {q(R+T)\sin\alpha}
 
 and
@@ -58,9 +58,7 @@
 NB: The 2nd virial coefficient of the cylinder is calculated based on
 the radius and 2 length values, and used as the effective radius for
-$S(Q)$ when $P(Q) \cdot S(Q)$ is applied.
+$S(q)$ when $P(q) \cdot S(q)$ is applied.
 
-The $\theta$ and $\phi$ parameters are not used for the 1D output. Our
-implementation of the scattering kernel and the 1D scattering intensity
-use the c-library from NIST.
+The $\theta$ and $\phi$ parameters are not used for the 1D output.
 
 Validation
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/models/cylinder.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -12,5 +12,5 @@
 .. math::
 
-    P(Q,\alpha) = {\text{scale} \over V} F^2(Q) + \text{background}
+    P(q,\alpha) = \frac{\text{scale}}{V} F^2(q) + \text{background}
 
 where
@@ -18,8 +18,8 @@
 .. math::
 
-    F(Q) = 2 (\Delta \rho) V
-           {\sin \left(Q\tfrac12 L\cos\alpha \right)
-               \over Q\tfrac12 L \cos \alpha}
-           {J_1 \left(Q R \sin \alpha\right) \over Q R \sin \alpha}
+    F(q) = 2 (\Delta \rho) V
+           \frac{\sin \left(q\tfrac12 L\cos\alpha \right)}
+                {q\tfrac12 L \cos \alpha}
+           \frac{J_1 \left(q R \sin \alpha\right)}{q R \sin \alpha}
 
 and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V$
@@ -44,6 +44,6 @@
 
 NB: The 2nd virial coefficient of the cylinder is calculated based on the
-radius and length values, and used as the effective radius for $S(Q)$
-when $P(Q) \cdot S(Q)$ is applied.
+radius and length values, and used as the effective radius for $S(q)$
+when $P(q) \cdot S(q)$ is applied.
 
 The output of the 1D scattering intensity function for randomly oriented
@@ -52,10 +52,8 @@
 .. math::
 
-    P(Q) = {\text{scale} \over V}
-        \int_0^{\pi/2} F^2(Q,\alpha) \sin \alpha\ d\alpha + \text{background}
+    P(q) = \frac{\text{scale}}{V}
+        \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background}
 
-The *theta* and *phi* parameters are not used for the 1D output. Our
-implementation of the scattering kernel and the 1D scattering intensity
-use the c-library from NIST.
+The $\theta$ and $\phi$ parameters are not used for the 1D output.
 
 Validation
@@ -82,10 +80,10 @@
 .. math::
 
-    P(Q) = \int_0^{\pi/2} d\phi
-        \int_0^\pi p(\theta, \phi) P_0(Q,\alpha) \sin \theta\ d\theta
+    P(q) = \int_0^{\pi/2} d\phi
+        \int_0^\pi p(\theta, \phi) P_0(q,\alpha) \sin \theta\ d\theta
 
 
 where $p(\theta,\phi)$ is the probability distribution for the orientation
-and $P_0(Q,\alpha)$ is the scattering intensity for the fully oriented
+and $P_0(q,\alpha)$ is the scattering intensity for the fully oriented
 system. Since we have no other software to compare the implementation of
 the intensity for fully oriented cylinders, we can compare the result of
@@ -129,5 +127,5 @@
 
 #             [ "name", "units", default, [lower, upper], "type", "description"],
-parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "",
+parameters = [["sld", "4e-6/Ang^2", 4, [-inf, inf], "",
                "Cylinder scattering length density"],
               ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
Index: sasmodels/models/dab.py
===================================================================
--- sasmodels/models/dab.py	(revision f224873e62472b942530cb10cfcab1f4e13ebb87)
+++ sasmodels/models/dab.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -1,8 +1,3 @@
 r"""
-
-Scattering model class for the DAB (Debye Anderson Brumberger) Model
-
-Definition
-----------
 
 Calculates the scattering from a randomly distributed, two-phase system based on
@@ -11,31 +6,25 @@
 measure of the average spacing between regions of phase 1 and phase 2. **The
 model also assumes smooth interfaces between the phases** and hence exhibits
-Porod behavior (I ~ *q*\ :sup:`-4`) at large *q* (*QL* >> 1).
+Porod behavior $(I \sim q^{-4})$ at large $q$, $(qL \gg 1)$.
 
 The DAB model is ostensibly a development of the earlier Debye-Bueche model.
 
-*Definition*
+Definition
+----------
 
-.. math:: I(q) = \text{scale}\cdot\frac{L^3}{(1 + (q \cdot L)^2)^2} + \text{background}
+.. math::
+
+    I(q) = \text{scale}\cdot\frac{L^3}{(1 + (q\cdot L)^2)^2} + \text{background}
 
 where scale is
 
-.. math:: \text{scale} = 8 \pi \phi (1-\phi)(\Delta \rho)^2
+.. math:: \text{scale} = 8 \pi \phi (1-\phi) \Delta\rho^2
 
-and the parameter *L* is the correlation length.
+and the parameter $L$ is the correlation length.
 
-For 2D data: The 2D scattering intensity is calculated in the same way as 1D,
-where the *q* vector is defined as
+For 2D data the scattering intensity is calculated in the same way as 1D,
+where the $q$ vector is defined as
 
 .. math:: q = \sqrt{q_x^2 + q_y^2}
-
-==============  ========  =============
-Parameter name  Units     Default value
-==============  ========  =============
-scale           None      1.0
-corr length L   |Ang|     50.0
-background      |cm^-1|   0.0
-==============  ========  =============
-
 
 .. figure:: img/dab_1d.jpg
@@ -44,6 +33,6 @@
 
 
-Reference
----------
+References
+----------
 
 P Debye, H R Anderson, H Brumberger, *Scattering by an Inhomogeneous Solid. II.
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision 3e428ec139734cd9698fc2d854c41eb3498c4c47)
+++ sasmodels/models/ellipsoid.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -12,5 +12,5 @@
 .. math::
 
-    P(Q,\alpha) = {\text{scale} \over V} F^2(Q) + \text{background}
+    P(q,\alpha) = \frac{\text{scale}}{V} F^2(q) + \text{background}
 
 where
@@ -18,7 +18,7 @@
 .. math::
 
-    F(Q) = {3 (\Delta rho)) V (\sin[Qr(R_p,R_e,\alpha)]
-                - \cos[Qr(R_p,R_e,\alpha)])
-            \over [Qr(R_p,R_e,\alpha)]^3 }
+    F(q) = \frac{3 \Delta \rho V (\sin[qr(R_p,R_e,\alpha)]
+                - \cos[qr(R_p,R_e,\alpha)])}
+                {[qr(R_p,R_e,\alpha)]^3}
 
 and
@@ -45,9 +45,9 @@
 NB: The 2nd virial coefficient of the solid ellipsoid is calculated based
 on the $R_p$ and $R_e$ values, and used as the effective radius for
-$S(Q)$ when $P(Q) \cdot S(Q)$ is applied.
+$S(q)$ when $P(q) \cdot S(q)$ is applied.
 
 .. _ellipsoid-1d:
 
-.. figure:: img/ellipsoid_1d.JPG
+.. figure:: img/ellipsoid_1d.jpg
 
     The output of the 1D scattering intensity function for randomly oriented
@@ -55,11 +55,9 @@
 
 
-The $\theta$ and $\phi$ parameters are not used for the 1D output. Our
-implementation of the scattering kernel and the 1D scattering intensity
-use the c-library from NIST.
+The $\theta$ and $\phi$ parameters are not used for the 1D output.
 
 .. _ellipsoid-geometry:
 
-.. figure:: img/ellipsoid_geometry.JPG
+.. figure:: img/ellipsoid_geometry.jpg
 
     The angles for oriented ellipsoid.
@@ -100,15 +98,16 @@
     *contrast* = 3e-6 |Ang^-2|, and *background* = 0.0 |cm^-1|.
 
-The discrepancy above *q* = 0.3 |cm^-1| is due to the way the form factors
+The discrepancy above $q$ = 0.3 |cm^-1| is due to the way the form factors
 are calculated in the c-library provided by NIST. A numerical integration
-has to be performed to obtain $P(Q)$ for randomly oriented particles.
+has to be performed to obtain $P(q)$ for randomly oriented particles.
 The NIST software performs that integration with a 76-point Gaussian
-quadrature rule, which will become imprecise at high $Q$ where the amplitude
-varies quickly as a function of $Q$. The SasView result shown has been
+quadrature rule, which will become imprecise at high $q$ where the amplitude
+varies quickly as a function of $q$. The SasView result shown has been
 obtained by summing over 501 equidistant points. Our result was found
-to be stable over the range of $Q$ shown for a number of points higher
+to be stable over the range of $q$ shown for a number of points higher
 than 500.
 
-REFERENCE
+References
+----------
 
 L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum,
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/models/fcc.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -16,8 +16,8 @@
 .. math::
 
-    I(q) = \frac{\text{scale}}{V_p} P(q) Z(q)
+    I(q) = \frac{\text{scale}}{V_p} V_\text{lattice} P(q) Z(q)
 
 where *scale* is the volume fraction of spheres, $V_p$ is the volume of
-the primary particle, *V(lattice)* is a volume correction for the crystal
+the primary particle, $V_\text{lattice}$ is a volume correction for the crystal
 structure, $P(q)$ is the form factor of the sphere (normalized), and $Z(q)$
 is the paracrystalline structure factor for a face-centered cubic structure.
@@ -27,5 +27,5 @@
 
 The lattice correction (the occupied volume of the lattice) for a
-face-centered cubic structure of particles of radius *R* and nearest
+face-centered cubic structure of particles of radius $R$ and nearest
 neighbor separation $D$ is
 
@@ -78,10 +78,10 @@
     1D plot in the linear scale using the default values (w/200 data point).
 
-The 2D (Anisotropic model) is based on the reference below where *I(q)* is
+The 2D (Anisotropic model) is based on the reference below where $I(q)$ is
 approximated for 1d scattering. Thus the scattering pattern for 2D may not
 be accurate. Note that we are not responsible for any incorrectness of the
 2D model computation.
 
-.. figure:: img/crystal_orientation.gif
+.. figure:: img/crystal_orientation.png
 
     Orientation of the crystal with respect to the scattering plane.
@@ -91,6 +91,6 @@
     2D plot using the default values (w/200X200 pixels).
 
-Reference
----------
+References
+----------
 
 Hideki Matsuoka et. al. *Physical Review B*, 36 (1987) 1754-1765
Index: sasmodels/models/gaussian_peak.py
===================================================================
--- sasmodels/models/gaussian_peak.py	(revision 0e9048fc51455645fd5530e83ead6393f51bbf74)
+++ sasmodels/models/gaussian_peak.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -26,6 +26,6 @@
     1D plot using the default values (w/500 data points).
 
-Reference
----------
+References
+----------
 
 None.
Index: sasmodels/models/guinier.py
===================================================================
--- sasmodels/models/guinier.py	(revision 66ebdd668828275e87ada92bedafaa22888a51db)
+++ sasmodels/models/guinier.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -1,5 +1,3 @@
 r"""
-Guinier (Model)
-
 Definition
 ----------
@@ -7,15 +5,19 @@
 This model fits the Guinier function
 
-.. math:: Q_1=\frac{1}{R_g}\sqrt{\frac{(m-s)(3-s)}{2}}
+.. math:: q_1=\frac{1}{R_g}\sqrt{\frac{(m-s)(3-s)}{2}}
 
-to the data directly without any need for linearisation (*cf*. Ln *I(q)* vs *q*\ :sup:`2`).
+to the data directly without any need for linearisation
+(*cf*. $\ln I(q)$ vs $q^2$\ ).
 
-For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as
+For 2D data the scattering intensity is calculated in the same way as 1D,
+where the $q$ vector is defined as
 
 .. math:: q=\sqrt{q_x^2 + q_y^2}
 
-REFERENCE
+References
+----------
 
-A Guinier and G Fournet, *Small-Angle Scattering of X-Rays*, John Wiley & Sons, New York (1955)
+A Guinier and G Fournet, *Small-Angle Scattering of X-Rays*,
+John Wiley & Sons, New York (1955)
 """
 
Index: sasmodels/models/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision 3e428ec139734cd9698fc2d854c41eb3498c4c47)
+++ sasmodels/models/hardsphere.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -7,5 +7,5 @@
 potential is
 
-.. math:
+.. math::
 
     U(r) = \begin{cases}
@@ -14,5 +14,5 @@
     \end{cases}
 
-where *r* is the distance from the center of the sphere of a radius *R*.
+where $r$ is the distance from the center of the sphere of a radius $R$.
 
 For a 2D plot, the wave transfer is defined as
@@ -23,9 +23,10 @@
 
 
-.. image:: img/HardSphere_1d.jpg
+.. figure:: img/HardSphere_1d.jpg
 
-*Figure. 1D plot using the default values (in linear scale).*
+    1D plot using the default values (in linear scale).
 
-REFERENCE
+References
+----------
 
 J K Percus, J Yevick, *J. Phys. Rev.*, 110, (1958) 1
Index: sasmodels/models/hollow_cylinder.py
===================================================================
--- sasmodels/models/hollow_cylinder.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/models/hollow_cylinder.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -6,7 +6,7 @@
 .. math::
 
-    P(q) = \text{scale} \langle F^2 \rangle/V_\text{shell} + \text{background}
+    P(q) = \text{scale} \left<F^2\right>/V_\text{shell} + \text{background}
 
-where the averaging $\langle \rangle$ is applied only for the 1D calculation.
+where the averaging $\left<\ldots\right>$ is applied only for the 1D calculation.
 
 The inside and outside of the hollow cylinder are assumed have the same SLD.
@@ -19,17 +19,18 @@
 .. math::
 
-    \begin{eqnarray}
-    P(q)           &=& (\text{scale})V_\text{shell}\Delta\rho^2
+    %\begin{align*} % isn't working with pdflatex
+    \begin{array}{rl}
+    P(q)           &= (\text{scale})V_\text{shell}\Delta\rho^2
             \int_0^{1}\Psi^2
             \left[q_z, R_\text{shell}(1-x^2)^{1/2},
                        R_\text{core}(1-x^2)^{1/2}\right]
             \left[\frac{\sin(qHx)}{qHx}\right]^2 dx \\
-    \Psi[q,y,z]    &=& \frac{1}{1-\gamma^2}
+    \Psi[q,y,z]    &= \frac{1}{1-\gamma^2}
             \left[ \Lambda(qy) - \gamma^2\Lambda(qz) \right] \\
-    \Lambda(a)     &=& 2 J_1(a) / a \\
-    \gamma         &=& R_\text{core} / R_\text{shell} \\
-    V_\text{shell} &=& \pi \left(R_\text{shell}^2 - R_\text{core}^2 \right)L \\
-    J_1(x)         &=& \frac{(\sin(x)-x\cdot \cos(x))}{x^2} \\
-    \end{eqnarray}
+    \Lambda(a)     &= 2 J_1(a) / a \\
+    \gamma         &= R_\text{core} / R_\text{shell} \\
+    V_\text{shell} &= \pi \left(R_\text{shell}^2 - R_\text{core}^2 \right)L \\
+    J_1(x)         &= (\sin(x)-x\cdot \cos(x)) / x^2 \\
+    \end{array}
 
 where *scale* is a scale factor and $J_1$ is the 1st order
@@ -42,5 +43,5 @@
 **NB**: The 2nd virial coefficient of the cylinder is calculated
 based on the radius and 2 length values, and used as the effective radius
-for $S(Q)$ when $P(Q) * S(Q)$ is applied.
+for $S(q)$ when $P(q) \cdot S(q)$ is applied.
 
 In the parameters, the contrast represents SLD :sub:`shell` - SLD :sub:`solvent`
@@ -59,6 +60,6 @@
     Examples of the angles for oriented pp against the detector plane.
 
-Reference
----------
+References
+----------
 
 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and
@@ -88,6 +89,6 @@
               ["sld", "1/Ang^2", 6.3, [-inf, inf], "", "Cylinder sld"],
               ["solvent_sld", "1/Ang^2", 1, [-inf, inf], "", "Solvent sld"],
-              ["theta", "[deg]", 90, [-360, 360], "orientation", "Theta angle"],
-              ["phi", "[deg]", 0, [-360, 360], "orientation", "Phi angle"],
+              ["theta", "degrees", 90, [-360, 360], "orientation", "Theta angle"],
+              ["phi", "degrees", 0, [-360, 360], "orientation", "Phi angle"],
               ]
 
Index: sasmodels/models/lamellar.py
===================================================================
--- sasmodels/models/lamellar.py	(revision 38d8774e2ac4fb70034f42ecb95289ad2fba4613)
+++ sasmodels/models/lamellar.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -5,9 +5,9 @@
 ----------
 
-The scattering intensity *I(Q)* is
+The scattering intensity $I(q)$ is
 
 .. math::
 
-    I(Q) = 2\pi{P(Q) \over \delta Q^2}
+    I(q) = \frac{2\pi P(q)}{\delta q^2}
 
 
@@ -16,19 +16,16 @@
 .. math::
 
-    P(Q) = {2\Delta\rho^2 \over Q^2}(1-cos(Q\delta))
+    P(q) = \frac{2\Delta\rho^2}{q^2}(1-cos(q\delta))
 
 
-where |delta| = bilayer thickness.
+where $\delta$ is the bilayer thickness.
 
-The 2D scattering intensity is calculated in the same way as 1D, where the $Q$ vector is defined as
+The 2D scattering intensity is calculated in the same way as 1D, where
+the $q$ vector is defined as
 
 .. math::
 
-    Q = \sqrt{Q_x^2 + Q_y^2}
+    q = \sqrt{q_x^2 + q_y^2}
 
-
-
-Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research
-(Kline, 2006).
 
 .. figure:: img/lamellar_1d.jpg
@@ -37,6 +34,6 @@
 
 
-Reference
----------
+References
+----------
 
 F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502
Index: sasmodels/models/lamellarCaille.py
===================================================================
--- sasmodels/models/lamellarCaille.py	(revision 3e428ec139734cd9698fc2d854c41eb3498c4c47)
+++ sasmodels/models/lamellarCaille.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -3,9 +3,12 @@
 This model provides the scattering intensity, $I(q) = P(q) S(q)$, for a
 lamellar phase where a random distribution in solution are assumed.
-Here a Caille $S(Q)$ is used for the lamellar stacks.
+Here a Caille $S(q)$ is used for the lamellar stacks.
+
+Definition
+----------
 
 The scattering intensity $I(q)$ is
 
-.. math:
+.. math::
 
     I(q) = 2\pi \frac{P(q)S(q)}{\delta q^2}
@@ -13,5 +16,5 @@
 The form factor is
 
-.. math:
+.. math::
 
     P(q) = \frac{2\Delta\rho^2}{q^2}\left(1-\cos q\delta \right)
@@ -19,5 +22,5 @@
 and the structure factor is
 
-.. math:
+.. math::
 
     S(q) = 1 + 2 \sum_1^{N-1}\left(1-\frac{n}{N}\right)
@@ -26,11 +29,16 @@
 where
 
-.. math:
+.. math::
 
-    \begin{eqnarray}
-    \alpha(n) &=& \frac{\eta_{cp}}{4\pi^2} \left(\ln(\pi n)+\gamma_E\right)  \\
-    \gamma_E &=& 0.5772156649 && \text{Euler's constant} \\
-    \eta_{cp} &=& \frac{q_o^2k_B T}{8\pi\sqrt{K\overline{B}}} && \text{Caille constant}
-    \end{eqnarray}
+    %\begin{align*} % isn't working with pdflatex
+    \begin{array}{rll}
+    \alpha(n) &= \frac{\eta_{cp}}{4\pi^2} \left(\ln(\pi n)+\gamma_E\right)
+              & \\
+    \gamma_E  &= 0.5772156649
+              & \text{Euler's constant} \\
+    \eta_{cp} &= \frac{q_o^2k_B T}{8\pi\sqrt{K\overline{B}}}
+              & \text{Caille constant} \\
+              &
+    \end{array}
 
 Here $d$ = (repeat) spacing, $\delta$ = bilayer thickness,
@@ -54,15 +62,10 @@
     q = \sqrt{q_x^2 + q_y^2}
 
-The returned value is in units of |cm^-1|, on absolute scale.
+.. figure:: img/lamellarCaille_1d.jpg
 
-.. image:: img/lamellarCaille_1d.jpg
+    1D plot using the default values (w/6000 data point).
 
-*Figure. 1D plot using the default values (w/6000 data point).*
-
-Our model uses the form factor calculations as implemented in a c library
-provided by the NIST Center for Neutron Research (Kline, 2006).
-
-REFERENCE
----------
+References
+----------
 
 F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502
Index: sasmodels/models/lamellarCailleHG.py
===================================================================
--- sasmodels/models/lamellarCailleHG.py	(revision 3e428ec139734cd9698fc2d854c41eb3498c4c47)
+++ sasmodels/models/lamellarCailleHG.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -2,5 +2,5 @@
 r"""
 This model provides the scattering intensity, $I(q) = P(q)S(q)$, for a lamellar
-phase where a random distribution in solution are assumed. Here a Caille $S(Q)$
+phase where a random distribution in solution are assumed. Here a Caille $S(q)$
 is used for the lamellar stacks.
 
@@ -31,9 +31,14 @@
 .. math::
 
-    \begin{eqnarray}
-    \alpha(n) &=& \frac{\eta_{cp}}{4\pi^2} \left(\ln(\pi n)+\gamma_E\right)  \\
-    \gamma_E &=& 0.5772156649&&\text{Euler's constant} \\
-    \eta_{cp} &=& \frac{q_o^2k_B T}{8\pi\sqrt{K\overline{B}}} && \text{Caille constant}
-    \end{eqnarray}
+    %\begin{align*} % isn't working with pdflatex
+    \begin{array}{rll}
+    \alpha(n) &= \frac{\eta_{cp}}{4\pi^2} \left(\ln(\pi n)+\gamma_E\right)
+              &  \\
+    \gamma_E  &= 0.5772156649
+              & \text{Euler's constant} \\
+    \eta_{cp} &= \frac{q_o^2k_B T}{8\pi\sqrt{K\overline{B}}}
+              & \text{Caille constant} \\
+              &
+    \end{array}
 
 
@@ -59,14 +64,10 @@
     q = \sqrt{q_x^2 + q_y^2}
 
-The returned value is in units of |cm^-1|, on absolute scale.
+.. figure:: img/lamellarCailleHG_1d.jpg
 
-.. image:: img/lamellarCailleHG_1d.jpg
+    1D plot using the default values (w/6000 data point).
 
-*Figure. 1D plot using the default values (w/6000 data point).*
-
-Our model uses the form factor calculations implemented in a C library provided
-by the NIST Center for Neutron Research (Kline, 2006).
-
-REFERENCE
+References
+----------
 
 F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502
Index: sasmodels/models/lamellarFFHG.py
===================================================================
--- sasmodels/models/lamellarFFHG.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/models/lamellarFFHG.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -1,5 +1,5 @@
 # Note: model title and parameter table are inserted automatically
 r"""
-This model provides the scattering intensity, *I(q)*, for a lyotropic lamellar
+This model provides the scattering intensity, $I(q)$, for a lyotropic lamellar
 phase where a random distribution in solution are assumed. The SLD of the head
 region is taken to be different from the SLD of the tail region.
@@ -25,10 +25,10 @@
         \right\rbrace^2
 
-where $\delta_T$ is *tail_length*, \delta_H is *head_length*,
+where $\delta_T$ is *tail_length*, $\delta_H$ is *head_length*,
 $\Delta\rho_H$ is the head contrast (*head_sld* $-$ *solvent_sld*),
 and $\Delta\rho_T$ is tail contrast (*sld* $-$ *solvent_sld*).
 
 The 2D scattering intensity is calculated in the same way as 1D, where
-the *q* vector is defined as
+the $q$ vector is defined as
 
 .. math::
Index: sasmodels/models/lamellarPC.py
===================================================================
--- sasmodels/models/lamellarPC.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/models/lamellarPC.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -23,8 +23,8 @@
 
 - *spacing* is the average distance between adjacent layers
-  $\langle D \rangle$, and
+  $\left<D\right>$, and
 
 - *spacing_polydisp* is the relative standard deviation of the Gaussian
-  layer distance distribution $\sigma_D / \langle D \rangle$.
+  layer distance distribution $\sigma_D / \left<D\right>$.
 
 
@@ -49,5 +49,5 @@
 
 
-    Z_N(q) = \frac{1 - w^2}{1 + w^2 - 2w \cos(q \langle D \rangle)}
+    Z_N(q) = \frac{1 - w^2}{1 + w^2 - 2w \cos(q \left<D\right>)}
         + x_N S_N + (1 - x_N) S_{N+1}
 
@@ -56,5 +56,5 @@
 .. math::
 
-    S_N(q) = \frac{a_N}{N}[1 + w^2 - 2 w \cos(q \langle D \rangle)]^2
+    S_N(q) = \frac{a_N}{N}[1 + w^2 - 2 w \cos(q \left<D\right>)]^2
 
 and
@@ -62,8 +62,10 @@
 .. math::
 
-    a_N = 4w^2 - 2(w^3 + w) \cos(q \langle D \rangle)
-        - 4w^{N+2}\cos(Nq \langle D \rangle)
-        + 2 w^{N+3}\cos[(N-1)q \langle D \rangle]
-        + 2w^{N+1}\cos[(N+1)q \langle D \rangle]
+    \begin{array}{rcl}
+    a_N &=& 4w^2 - 2(w^3 + w) \cos(q \left<D\right>)
+        - 4w^{N+2}\cos(Nq \left<D\right>) \\
+        &&{} + 2 w^{N+3}\cos[(N-1)q \left<D\right>]
+        + 2w^{N+1}\cos[(N+1)q \left<D\right>]
+    \end{array}
 
 for the layer spacing distribution $w = \exp(-\sigma_D^2 q^2/2)$.
Index: sasmodels/models/lorentz.py
===================================================================
--- sasmodels/models/lorentz.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/models/lorentz.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -7,9 +7,9 @@
 The Ornstein-Zernicke model is defined by
 
-.. math:: I(q)=\frac{\text{scale}}{1+(qL)^2}+\text{bkg}
+.. math:: I(q)=\frac{\text{scale}}{1+(qL)^2}+\text{background}
 
-The parameter *L* is the screening length.
+The parameter $L$ is the screening length *cor_length*.
 
-For 2D data: The 2D scattering intensity is calculated in the same way as 1D, 
+For 2D data the scattering intensity is calculated in the same way as 1D,
 where the $q$ vector is defined as
 
@@ -20,6 +20,6 @@
     1D plot using the default values (w/200 data point).
 
-Reference
----------
+References
+----------
 
 L.S. Qrnstein and F. Zernike, *Proc. Acad. Sci. Amsterdam* 17, 793 (1914), and
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/models/parallelepiped.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -7,5 +7,5 @@
 ----------
 
-This model provides the form factor, *P(q)*, for a rectangular parallelepiped
+This model provides the form factor, $P(q)$, for a rectangular parallelepiped
 (below) where the form factor is normalized by the volume of the
 parallelepiped. If you need to apply polydispersity, see also
@@ -16,8 +16,8 @@
 .. math::
 
-    P(Q) = {\text{scale} \over V} F^2(Q) + \text{background}
+    P(q) = \frac{\text{scale}}{V} F^2(q) + \text{background}
 
-where the volume *V* = *A B C* and the averaging < > is applied over all
-orientations for 1D.
+where the volume $V = A B C$ and the averaging $\left<\ldots\right>$ is
+applied over all orientations for 1D.
 
 .. figure:: img/parallelepiped.jpg
@@ -25,12 +25,13 @@
    Parallelepiped with the corresponding definition of sides.
 
-The edge of the solid must satisfy the condition that** *A* < *B* < *C*.
-Then, assuming *a* = *A* / *B* < 1, *b* = *B* / *B* = 1, and
-*c* = *C* / *B* > 1, the form factor is
+The edge of the solid must satisfy the condition that $A < B < C$.
+Then, assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the
+form factor is
 
 .. math::
 
-    P(q) = \frac{\textstyle{scale}}{V}\int_0^1 \phi(\mu \sqrt{1-\sigma^2},a)
-    [S(\mu c \sigma/2)]^2 d\sigma
+    P(q) = \frac{\text{scale}}{V}\int_0^1
+        \phi\left(\mu \sqrt{1-\sigma^2},a\right)
+        \left[S(\mu c \sigma/2)\right]^2 d\sigma
 
 with
@@ -38,6 +39,8 @@
 .. math::
 
-    \phi(\mu,a) = \int_0^1 \{S[\frac{\mu}{2}\cos(\frac{\pi}{2}u)]
-    S[\frac{\mu a}{2}\sin(\frac{\pi}{2}u)]\}^2 du
+    \phi(\mu,a) = \int_0^1
+        \left\{S\left[\frac{\mu}{2}\cos(\frac{\pi}{2}u)\right]
+               S\left[\frac{\mu a}{2}\sin(\frac{\pi}{2}u)\right]
+               \right\}^2 du
 
     S(x) = \frac{\sin x}{x}
@@ -52,17 +55,17 @@
 
 The scattering intensity per unit volume is returned in units of |cm^-1|;
-ie, *I(q)* = |phi| *P(q)*\ .
+i.e., $I(q) = \phi P(q)$.
 
 NB: The 2nd virial coefficient of the parallelpiped is calculated based on
-the averaged effective radius (= sqrt(*short_a* \* *short_b* / |pi|)) and
-length(= *long_c*) values, and used as the effective radius for
-*S(Q)* when *P(Q)* \* *S(Q)* is applied.
+the averaged effective radius $(=\sqrt{A B / \pi})$ and
+length $(= C)$ values, and used as the effective radius for
+$S(q)$ when $P(q) \cdot S(q)$ is applied.
 
 To provide easy access to the orientation of the parallelepiped, we define
-three angles |theta|, |phi| and |bigpsi|. The definition of |theta| and |phi|
-is the same as for the cylinder model (see also figures below).
-The angle |bigpsi| is the rotational angle around the *long_c* axis against
-the *q* plane. For example, |bigpsi| = 0 when the *short_b* axis is parallel
-to the *x*-axis of the detector.
+three angles $\theta$, $\phi$ and $\Psi$. The definition of $\theta$ and
+$\phi$ is the same as for the cylinder model (see also figures below).
+The angle $\Psi$ is the rotational angle around the $C$ axis against
+the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel
+to the $x$-axis of the detector.
 
 
@@ -85,10 +88,10 @@
 angles. The Figure below shows the comparison where the solid dot refers to
 averaged 2D while the line represents the result of the 1D calculation (for
-the averaging, 76, 180, 76 points are taken for the angles of |theta|, |phi|,
-and |psi| respectively).
+the averaging, 76, 180, 76 points are taken for the angles of $\theta$,
+$\phi$, and $\Psi$ respectively).
 
 .. _parallelepiped-compare:
 
-.. figure:: img/parallelepiped_compare.gif
+.. figure:: img/parallelepiped_compare.png
 
    Comparison between 1D and averaged 2D.
Index: sasmodels/models/power_law.py
===================================================================
--- sasmodels/models/power_law.py	(revision 525f3a9ece3be16c08905b2a1525c19bd9db33de)
+++ sasmodels/models/power_law.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -10,5 +10,5 @@
 .. math::
 
-    I(q) = \text{scale} \cdot q^{-power} + \text{background}
+    I(q) = \text{scale} \cdot q^{-\text{power}} + \text{background}
 
 Note the minus sign in front of the exponent. The exponent *power* 
@@ -19,9 +19,10 @@
 combining this model with other models.
 
-.. image:: img/power_law_1d.jpg
+.. figure:: img/power_law_1d.jpg
 
-*Figure. 1D plot using the default values (w/200 data point).*
+    1D plot using the default values (w/200 data point).
 
-REFERENCE
+References
+----------
 
 None.
@@ -47,9 +48,9 @@
     inten = (q**-power)
     return inten
-Iq.vectorized = True  # Iq accepts an array of Q values
+Iq.vectorized = True  # Iq accepts an array of q values
 
 def Iqxy(qx, qy, *args):
     return Iq(sqrt(qx ** 2 + qy ** 2), *args)
-Iqxy.vectorized = True # Iqxy accepts an array of Qx, Qy values
+Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
 
 demo = dict(scale=1.0,
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 0fa687db9fc44ec8f0be13c81ef29ce8071d8a81)
+++ sasmodels/models/sphere.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -11,10 +11,10 @@
 .. math::
 
-    I(Q) = \frac{\text{scale}}{V} \cdot \left[ \
-        3V(\Delta\rho) \cdot \frac{\sin(QR) - QR\cos(QR))}{(QR)^3} \
+    I(q) = \frac{\text{scale}}{V} \cdot \left[
+        3V(\Delta\rho) \cdot \frac{\sin(qr) - qr\cos(qr))}{(qr)^3}
         \right]^2 + \text{background}
 
 where *scale* is a volume fraction, $V$ is the volume of the scatterer,
-$R$ is the radius of the sphere, *background* is the background level and
+$r$ is the radius of the sphere, *background* is the background level and
 *sld* and *solvent_sld* are the scattering length densities (SLDs) of the
 scatterer and the solvent respectively.
@@ -27,7 +27,4 @@
 The 2D scattering intensity is the same as above, regardless of the
 orientation of $\vec q$.
-
-Our model uses the form factor calculations as defined in the IGOR
-package provided by the NIST Center for Neutron Research (Kline, 2006).
 
 Validation
@@ -49,6 +46,6 @@
 
 
-Reference
----------
+References
+----------
 
 A Guinier and G. Fournet, *Small-Angle Scattering of X-Rays*,
@@ -63,7 +60,7 @@
 title = "Spheres with uniform scattering length density"
 description = """\
-P(q)=(scale/V)*[3V(sld-solvent_sld)*(sin(qR)-qRcos(qR))
-                /(qR)^3]^2 + background
-    R: radius of sphere
+P(q)=(scale/V)*[3V(sld-solvent_sld)*(sin(qr)-qr cos(qr))
+                /(qr)^3]^2 + background
+    r: radius of sphere
     V: The volume of the scatter
     sld: the SLD of the sphere
@@ -93,5 +90,5 @@
     double sn, cn;
     SINCOS(qr, sn, cn);
-    // Use taylor series for low Q to avoid cancellation error.  Tested against
+    // Use taylor series for low q to avoid cancellation error.  Tested against
     // the following expression in quad precision:
     //     3.0*(sn-qr*cn)/(qr*qr*qr);
Index: sasmodels/models/spherepy.py
===================================================================
--- sasmodels/models/spherepy.py	(revision d138d4377bdeeb2f6748c00827e90f3421165989)
+++ sasmodels/models/spherepy.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -11,10 +11,10 @@
 .. math::
 
-    I(Q) = \frac{\text{scale}}{V} \cdot \left[ \
-        3V(\Delta\rho) \cdot \frac{\sin(QR) - QR\cos(QR))}{(QR)^3} \
+    I(q) = \frac{\text{scale}}{V} \cdot \left[
+        3V(\Delta\rho) \cdot \frac{\sin(qr) - qr\cos(qr))}{(qr)^3}
         \right]^2 + \text{background}
 
 where *scale* is a volume fraction, $V$ is the volume of the scatterer,
-$R$ is the radius of the sphere, *background* is the background level and
+$r$ is the radius of the sphere, *background* is the background level and
 *sld* and *solvent_sld* are the scattering length densities (SLDs) of the
 scatterer and the solvent respectively.
@@ -27,7 +27,4 @@
 The 2D scattering intensity is the same as above, regardless of the
 orientation of $\vec q$.
-
-Our model uses the form factor calculations as defined in the IGOR
-package provided by the NIST Center for Neutron Research (Kline, 2006).
 
 Validation
@@ -49,6 +46,6 @@
 
 
-Reference
----------
+References
+----------
 
 A Guinier and G. Fournet, *Small-Angle Scattering of X-Rays*,
@@ -64,7 +61,7 @@
 title = "Spheres with uniform scattering length density"
 description = """\
-P(q)=(scale/V)*[3V(sld-solvent_sld)*(sin(qR)-qRcos(qR))
-                /(qR)^3]^2 + background
-    R: radius of sphere
+P(q)=(scale/V)*[3V(sld-solvent_sld)*(sin(qr)-qr cos(qr))
+                /(qr)^3]^2 + background
+    r: radius of sphere
     V: The volume of the scatter
     sld: the SLD of the sphere
@@ -101,9 +98,9 @@
     fq = bes * (sld - solvent_sld) * form_volume(radius)
     return 1.0e-4 * fq ** 2
-Iq.vectorized = True  # Iq accepts an array of Q values
+Iq.vectorized = True  # Iq accepts an array of q values
 
 def Iqxy(qx, qy, sld, solvent_sld, radius):
     return Iq(sqrt(qx ** 2 + qy ** 2), sld, solvent_sld, radius)
-Iqxy.vectorized = True  # Iqxy accepts arrays of Qx, Qy values
+Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values
 
 def sesans(z, sld, solvent_sld, radius):
Index: sasmodels/models/stickyhardsphere.py
===================================================================
--- sasmodels/models/stickyhardsphere.py	(revision 3e428ec139734cd9698fc2d854c41eb3498c4c47)
+++ sasmodels/models/stickyhardsphere.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -4,22 +4,33 @@
 with a narrow attractive well. A perturbative solution of the Percus-Yevick
 closure is used. The strength of the attractive well is described in terms
-of "stickiness" as defined below. The returned value is a dimensionless
-structure factor, *S(q)*.
+of "stickiness" as defined below.
 
-The perturb (perturbation parameter), |epsilon|, should be held between 0.01
+The perturb (perturbation parameter), $\epsilon$, should be held between 0.01
 and 0.1. It is best to hold the perturbation parameter fixed and let
 the "stickiness" vary to adjust the interaction strength. The stickiness,
-|tau|, is defined in the equation below and is a function of both the
-perturbation parameter and the interaction strength. |tau| and |epsilon|
-are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ), the
-width of the square well, |bigdelta| (same units as *R*), and the depth of
-the well, *Uo*, in units of kT. From the definition, it is clear that
-smaller |tau| means stronger attraction.
+$\tau$, is defined in the equation below and is a function of both the
+perturbation parameter and the interaction strength. $\tau$ and $\epsilon$
+are defined in terms of the hard sphere diameter $(\sigma = 2 R)$, the
+width of the square well, $\Delta$ (same units as $R$\ ), and the depth of
+the well, $U_o$, in units of $kT$. From the definition, it is clear that
+smaller $\tau$ means stronger attraction.
 
-.. image:: img/stickyhardsphere_228.PNG
+.. math::
+
+    %\begin{align*} % isn't working with pdflatex
+    \begin{array}{rl}
+    \tau     &= \frac{1}{12\epsilon} \exp(u_o / kT) \\
+    \epsilon &= \Delta / (\sigma + \Delta) \\
+    \end{array}
 
 where the interaction potential is
 
-.. image:: img/stickyhardsphere_229.PNG
+.. math::
+
+    U(r) = \begin{cases}
+        \infty & r < \sigma \\
+        -U_o   & \sigma \leq r \leq \sigma + \Delta \\
+        0      & r > \sigma + \Delta
+        \end{cases}
 
 The Percus-Yevick (PY) closure was used for this calculation, and is an
@@ -28,5 +39,5 @@
 good agreement.
 
-The true particle volume fraction, |phi|, is not equal to *h*, which appears
+The true particle volume fraction, $\phi$, is not equal to $h$, which appears
 in most of the reference. The two are related in equation (24) of the
 reference. The reference also describes the relationship between this
@@ -34,7 +45,7 @@
 sphere) model by Baxter.
 
-NB: The calculation can go haywire for certain combinations of the input
+**NB**: The calculation can go haywire for certain combinations of the input
 parameters, producing unphysical solutions - in this case errors are
-reported to the command window and the *S(q)* is set to -1 (so it will
+reported to the command window and the $S(q)$ is set to -1 (so it will
 disappear on a log-log plot). Use tight bounds to keep the parameters to
 values that you know are physical (test them) and keep nudging them until
@@ -42,27 +53,19 @@
 
 In sasview the effective radius will be calculated from the parameters
-used in the form factor P(Q) that this S(Q) is combined with.
+used in the form factor $P(q)$ that this $S(q)$ is combined with.
 
-For 2D data: The 2D scattering intensity is calculated in the same way
-as 1D, where the *q* vector is defined as
+For 2D data the scattering intensity is calculated in the same way
+as 1D, where the $q$ vector is defined as
 
 .. math::
 
-    Q = \sqrt{Q_x^2 + Q_y^2}
+    q = \sqrt{q_x^2 + q_y^2}
 
-==============  ========  =============
-Parameter name  Units     Default value
-==============  ========  =============
-effect_radius   |Ang|     50
-perturb         None      0.05
-volfraction     None      0.1
-stickiness      K         0.2
-==============  ========  =============
+.. figure:: img/stickyhardsphere_1d.jpg
 
-.. image:: img/stickyhardsphere_230.jpg
+    1D plot using the default values (in linear scale).
 
-*Figure. 1D plot using the default values (in linear scale).*
-
-REFERENCE
+References
+----------
 
 S V G Menon, C Manohar, and K S Rao, *J. Chem. Phys.*, 95(12) (1991) 9186-9190
Index: sasmodels/models/teubner_strey.py
===================================================================
--- sasmodels/models/teubner_strey.py	(revision 0e9048fc51455645fd5530e83ead6393f51bbf74)
+++ sasmodels/models/teubner_strey.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -40,6 +40,6 @@
     1D plot using the default values (w/200 data point).
 
-Reference
----------
+References
+----------
 
 M Teubner, R Strey, *J. Chem. Phys.*, 87 (1987) 3195
@@ -77,9 +77,9 @@
 def Iq(q, a2, c1, c2):
     return 1. / np.polyval([c2, c1, a2], q**2)
-Iq.vectorized = True  # Iq accepts an array of Q values
+Iq.vectorized = True  # Iq accepts an array of q values
 
 def Iqxy(qx, qy, a2, c1, c2):
     return Iq(sqrt(qx**2 + qy**2), a2, c1, c2)
-Iqxy.vectorized = True  # Iqxy accepts arrays of Qx, Qy values
+Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values
 
 # ER defaults to 0.0
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision 3e428ec139734cd9698fc2d854c41eb3498c4c47)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
@@ -2,19 +2,17 @@
 # Note: model title and parameter table are inserted automatically
 r"""
-All three axes are of different lengths with $R_a \le R_b <= R_c$
+All three axes are of different lengths with $R_a \leq R_b \leq R_c$
 **Users should maintain this inequality for all calculations**.
 
 .. math::
 
-    P(Q) = \text{scale} V \left< F^2(Q) \right> + \text{background}
+    P(q) = \text{scale} V \left< F^2(q) \right> + \text{background}
 
 where the volume $V = 4/3 \pi R_a R_b R_c$, and the averaging
-$\left< \cdots \right>$ is applied over all orientations for 1D.
+$\left<\ldots\right>$ is applied over all orientations for 1D.
 
 .. figure:: img/triaxial_ellipsoid_geometry.jpg
 
     Ellipsoid schematic.
-
-The returned value is in units of |cm^-1|, on absolute scale.
 
 Definition
@@ -25,6 +23,6 @@
 .. math::
 
-    P(Q) = \frac{\text{scale}}{V}\int_0^1\int_0^1
-        \Phi^2(QR_a^2\cos^2( \pi x/2) + QR_b^2\sin^2(\pi y/2)(1-y^2) + c^2y^2)
+    P(q) = \frac{\text{scale}}{V}\int_0^1\int_0^1
+        \Phi^2(qR_a^2\cos^2( \pi x/2) + qR_b^2\sin^2(\pi y/2)(1-y^2) + R_c^2y^2)
         dx dy
 
@@ -40,5 +38,5 @@
 :num:`figure #triaxial-ellipsoid-angles`.
 The angle $\psi$ is the rotational angle around its own $c$ axis
-against the $Q$ plane. For example, $\psi = 0$ when the
+against the $q$ plane. For example, $\psi = 0$ when the
 $a$ axis is parallel to the $x$ axis of the detector.
 
@@ -52,11 +50,11 @@
 
 The contrast is defined as SLD(ellipsoid) - SLD(solvent).  In the
-parameters, *a* is the minor equatorial radius, *b* is the major
-equatorial radius, and c is the polar radius of the ellipsoid.
+parameters, $R_a$ is the minor equatorial radius, $R_b$ is the major
+equatorial radius, and $R_c$ is the polar radius of the ellipsoid.
 
 NB: The 2nd virial coefficient of the triaxial solid ellipsoid is
 calculated based on the polar radius $R_p = R_c$ and equatorial
 radius $R_e = \sqrt{R_a R_b}$, and used as the effective radius for
-$S(Q)$ when $P(Q) \cdot S(Q)$ is applied.
+$S(q)$ when $P(q) \cdot S(q)$ is applied.
 
 .. figure:: img/triaxial_ellipsoid_1d.jpg
@@ -81,8 +79,6 @@
     Comparison between 1D and averaged 2D.
 
-Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research
-(Kline, 2006)
-
-REFERENCE
+References
+----------
 
 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum,
