Index: doc/developer/overview.rst
===================================================================
--- doc/developer/overview.rst	(revision 3d4083969c45f693ad16ca85a10880859ed708e7)
+++ doc/developer/overview.rst	(revision 2a7e20e9c4ba72da5c6088e743a64f3345ee554d)
@@ -185,4 +185,12 @@
 jitter applied before particle orientation.
 
+When computing the orientation dispersity integral, the weights for
+the individual points depends on the map projection used to translate jitter
+angles into latitude/longitude.  The choice of projection is set by
+*sasmodels.generate.PROJECTION*, with the default *PROJECTION=1* for
+equirectangular and *PROJECTION=2* for sinusoidal.  The more complicated
+Guyou and Postel projections are not implemented. See explore.jitter.draw_mesh
+for details.
+
 For numerical integration within form factors etc. sasmodels is mostly using
 Gaussian quadrature with 20, 76 or 150 points depending on the model. It also
@@ -199,12 +207,41 @@
 Useful testing routines include:
 
-:mod:`asymint` a direct implementation of the surface integral for certain
-models to get a more trusted value for the 1D integral using a
-reimplementation of the 2D kernel in python and mpmath (which computes math
-functions to arbitrary precision). It uses $\theta$ ranging from 0 to $\pi$
-and $\phi$ ranging from 0 to $2\pi$. It perhaps would benefit from including
-the U-substitution for $\theta$.
-
-:mod:`check1d` uses sasmodels 1D integration and compares that with a
+The *sascomp* utility is used to view and compare models with different
+parameters and calculation engines. The usual case is to simply plot a
+model that you are developing::
+
+    python sascomp path/to/model.py
+
+Once the obvious problems are addressed, check the numerical precision
+across a variety of randomly generated inputs::
+
+    python sascomp -engine=single,double path/to/model.py -sets=10
+
+You can compare different parameter values for the same or different models.
+For example when looking along the long axis of a cylinder ($\theta=0$),
+dispersity in $\theta$ should be equivalent to dispersity in $\phi$
+when $\phi=90$::
+
+    python sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle \\
+    phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10
+
+It turns out that they are not because the equirectangular map projection
+weights the points by $\cos(\theta)$ so $\Delta\theta$ is not identical
+to $\Delta\phi$.  Setting *PROJECTION=2* in :mod:`sasmodels.generate` helps
+somewhat.  Postel would help even more in this case, though leading
+to distortions elsewhere.  See :mod:`sasmodels.compare` for many more details.
+
+*sascomp -ngauss=n* allows you to set the number of quadrature points used
+for the 1D integration for any model.  For example, a carbon nanotube with
+length 10 $\mu$\ m and radius 1 nm is not computed correctly at high $q$::
+
+    python sascomp cylinder length=100000 radius=10 -ngauss=76,500 -double -highq
+
+Note: ticket 702 gives some forms for long rods and thin disks that may
+be used in place of the integration, depending on $q$, radius and length; if
+the cylinder model is updated with these corrections then above call will show
+no difference.
+
+*explore/check1d.py* uses sasmodels 1D integration and compares that with a
 rectangle distribution in $\theta$ and $\phi$, with $\theta$ limits set to
 $\pm 90/\sqrt(3)$ and $\phi$ limits set to $\pm 180/\sqrt(3)$ [The rectangle
@@ -214,33 +251,57 @@
 similar reasoning.] This should rotate the sample through the entire
 $\theta$-$\phi$ surface according to the pattern that you see in jitter.py when
-you modify it to use 'rectangle' rather than 'gaussian' for its distribution
-without changing the viewing angle. In order to match the 1-D pattern for
-an arbitrary viewing angle on triaxial shapes, we need to integrate
+you use 'rectangle' rather than 'gaussian' for its distribution without
+changing the viewing angle. In order to match the 1-D pattern for an arbitrary
+viewing angle on triaxial shapes, we need to integrate
 over $\theta$, $\phi$ and $\Psi$.
 
-When computing the dispersity integral, weights are scaled by
-$|\cos(\delta \theta)|$ to account for the points in $\phi$ getting closer
-together as $\delta \theta$ increases.
-[This will probably change so that instead of adjusting the weights, we will
-adjust $\delta\theta$-$\delta\phi$ mesh so that the point density in
-$\delta\phi$ is lower at larger $\delta\theta$. The flag USE_SCALED_PHI in
-*kernel_iq.c* selects an alternative algorithm.]
-
-The integrated dispersion is computed at a set of $(qx, qy)$ points $(q
-\cos(\alpha), q \sin(\alpha))$ at some angle $\alpha$ (currently angle=0) for
-each $q$ used in the 1-D integration. The individual $q$ points should be
-equivalent to asymint rect-n when the viewing angle is set to
-$(\theta,\phi,\Psi) = (90, 0, 0)$. Such tests can help to validate that 2d
-models are consistent with 1d models.
-
-:mod:`sascomp -sphere=n` uses the same rectangular distribution as check1d to
-compute the pattern of the $q_x$-$q_y$ grid.
-
-The :mod:`sascomp` utility can be used for 2d as well as 1d calculations to
-compare results for two sets of parameters or processor types, for example
-these two oriented cylinders here should be equivalent.
-
-:mod:`\./sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10`
-
+*sascomp -sphere=n* uses the same rectangular distribution as check1d to
+compute the pattern of the $q_x$-$q_y$ grid.  This ought to produce a
+spherically symmetric pattern.
+
+*explore/precision.py* investigates the accuracy of individual functions
+on the different execution platforms.  This includes the basic special
+functions as well as more complex expressions used in scattering.  In many
+cases the OpenCL function in sasmodels will use a polynomial approximation
+over part of the range to improve accuracy, as shown in::
+
+    python explore/precision.py 3j1/x
+
+*explore/realspace.py* allows you to set up a Monte Carlo simulation of your
+model by sampling random points within and computing the 1D and 2D scattering
+patterns.  This was used to check the core shell parallelepiped example.  This
+is similar to the general sas calculator in sasview, though it uses different
+code.
+
+*explore/jitter.py* is for exploring different options for handling
+orientation and orientation dispersity.  It uses *explore/guyou.py* to
+generate the Guyou projection.
+
+*explore/asymint.py* is a direct implementation of the 1D integration for
+a number of different models.  It calculates the integral for a particular
+$q$ using several different integration schemes, including mpmath with 100
+bits of precision (33 digits), so you can use it to check the target values
+for the 1D model tests.  This is not a general purpose tool; you will need to
+edit the file to change the model parameters. It does not currently
+apply the $u=cos(\theta)$ substitution that many models are using
+internally so the 76-point gaussian quadrature may not match the value
+produced by the eqivalent model from sasmodels.
+
+*explore/symint.py* is for rotationally symmetric models (just cylinder for
+now), so it can compute an entire curve rather than a single $q$ point.  It
+includes code to compare the long cylinder approximation to cylinder.
+
+*explore/rpa.py* is for checking different implementations of the RPA model
+against calculations for specific blends.  This is a work in (suspended)
+progress.
+
+There are a few modules left over in *explore* that can probably be removed.
+*angular_pd.py* was an early version of *jitter.py*.  *sc.py* and *sc.c*
+was used to explore different calculations for paracrystal models; it has
+been absorbed into *asymint.py*. *transform_angles.py* translates orientation
+parameters from the SasView 3.x definition to sasmodels.
+
+*explore/angles.py* generates code for the view and jitter transformations.
+Keep this around since it may be needed if we add new projections.
 
 Testing
Index: doc/guide/orientation/orientation.rst
===================================================================
--- doc/guide/orientation/orientation.rst	(revision 82592da3721b7e2098cf6346cf224814189fadb3)
+++ doc/guide/orientation/orientation.rst	(revision 5fb0634b45bba451ca824eecf40346b1ed5305a7)
@@ -4,16 +4,21 @@
 ==================
 
-With two dimensional small angle diffraction data SasView will calculate
+With two dimensional small angle diffraction data sasmodels will calculate
 scattering from oriented particles, applicable for example to shear flow
 or orientation in a magnetic field.
 
 In general we first need to define the reference orientation
-of the particles with respect to the incoming neutron or X-ray beam. This
-is done using three angles: $\theta$ and $\phi$ define the orientation of
-the axis of the particle, angle $\Psi$ is defined as the orientation of
-the major axis of the particle cross section with respect to its starting
-position along the beam direction. The figures below are for an elliptical
-cross section cylinder, but may be applied analogously to other shapes of
-particle.
+of the particle's $a$-$b$-$c$ axes with respect to the incoming
+neutron or X-ray beam. This is done using three angles: $\theta$ and $\phi$
+define the orientation of the $c$-axis of the particle, and angle $\Psi$ is
+defined as the orientation of the major axis of the particle cross section
+with respect to its starting position along the beam direction (or
+equivalently, as rotation about the $c$ axis). There is an unavoidable
+ambiguity when $c$ is aligned with $z$ in that $\phi$ and $\Psi$ both
+serve to rotate the particle about $c$, but this symmetry is destroyed
+when $\theta$ is not a multiple of 180.
+
+The figures below are for an elliptical cross section cylinder, but may
+be applied analogously to other shapes of particle.
 
 .. note::
@@ -29,8 +34,8 @@
 
     Definition of angles for oriented elliptical cylinder, where axis_ratio
-    b/a is shown >1, Note that rotation $\theta$, initially in the $x$-$z$
+    b/a is shown >1. Note that rotation $\theta$, initially in the $x$-$z$
     plane, is carried out first, then rotation $\phi$ about the $z$-axis,
     finally rotation $\Psi$ is around the axis of the cylinder. The neutron
-    or X-ray beam is along the $z$ axis.
+    or X-ray beam is along the $-z$ axis.
 
 .. figure::
@@ -40,57 +45,154 @@
     with $\Psi$ = 0.
 
-Having established the mean direction of the particle we can then apply
-angular orientation distributions. This is done by a numerical integration
-over a range of angles in a similar way to particle size dispersity.
-In the current version of sasview the orientational dispersity is defined
-with respect to the axes of the particle.
+Having established the mean direction of the particle (the view) we can then
+apply angular orientation distributions (jitter). This is done by a numerical
+integration over a range of angles in a similar way to particle size
+dispersity. The orientation dispersity is defined with respect to the
+$a$-$b$-$c$ axes of the particle, with roll angle $\Psi$ about the $c$-axis,
+yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the
+$a$-axis.
+
+More formally, starting with axes $a$-$b$-$c$ of the particle aligned
+with axes $x$-$y$-$z$ of the laboratory frame, the orientation dispersity
+is applied first, using the
+`Tait-Bryan <https://en.wikipedia.org/wiki/Euler_angles#Conventions_2>`_
+$x$-$y'$-$z''$ convention with angles $\Delta\phi$-$\Delta\theta$-$\Delta\Psi$.
+The reference orientation then follows, using the
+`Euler angles <https://en.wikipedia.org/wiki/Euler_angles#Conventions>`_
+$z$-$y'$-$z''$ with angles $\phi$-$\theta$-$\Psi$.  This is implemented
+using rotation matrices as
+
+.. math::
+
+    R = R_z(\phi)\, R_y(\theta)\, R_z(\Psi)\,
+        R_x(\Delta\phi)\, R_y(\Delta\theta)\, R_z(\Delta\Psi)
+
+To transform detector $(q_x, q_y)$ values into $(q_a, q_b, q_c)$ for the
+shape in its canonical orientation, use
+
+.. math::
+
+    [q_a, q_b, q_c]^T = R^{-1} \, [q_x, q_y, 0]^T
+
+
+The inverse rotation is easily calculated by rotating the opposite directions
+in the reverse order, so
+
+.. math::
+
+    R^{-1} = R_z(-\Delta\Psi)\, R_y(-\Delta\theta)\, R_x(-\Delta\phi)\,
+             R_z(-\Psi)\, R_y(-\theta)\, R_z(-\phi)
+
 
 The $\theta$ and $\phi$ orientation parameters for the cylinder only appear
-when fitting 2d data. On introducing "Orientational Distribution" in
-the angles, "distribution of theta" and "distribution of phi" parameters will
+when fitting 2d data. On introducing "Orientational Distribution" in the
+angles, "distribution of theta" and "distribution of phi" parameters will
 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$
-of the cylinder, the $b$ and $a$ axes of the cylinder cross section. (When
-$\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the
-instrument.) The third orientation distribution, in $\Psi$, is about the $c$
-axis of the particle. Some experimentation may be required to understand the
-2d patterns fully. A number of different shapes of distribution are
-available, as described for polydispersity, see :ref:`polydispersityhelp` .
+of the cylinder, which correspond to the $b$ and $a$ axes of the cylinder
+cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and
+$X$ axes of the instrument.) The third orientation distribution, in $\Psi$,
+is about the $c$ axis of the particle. Some experimentation may be required
+to understand the 2d patterns fully. A number of different shapes of
+distribution are available, as described for size dispersity, see
+:ref:`polydispersityhelp`.
 
-Earlier versions of SasView had numerical integration issues in some
-circumstances when distributions passed through 90 degrees. The distributions
-in particle coordinates are more robust, but should still be approached with
-care for large ranges of angle.
+Given that the angular dispersion distribution is defined in cartesian space,
+over a cube defined by
+
+.. math::
+
+    [-\Delta \theta, \Delta \theta] \times
+    [-\Delta \phi, \Delta \phi] \times
+    [-\Delta \Psi, \Delta \Psi]
+
+but the orientation is defined over a sphere, we are left with a
+`map projection <https://en.wikipedia.org/wiki/List_of_map_projections>`_
+problem, with different tradeoffs depending on how values in $\Delta\theta$
+and $\Delta\phi$ are translated into latitude/longitude on the sphere.
+
+Sasmodels is using the
+`equirectangular projection <https://en.wikipedia.org/wiki/Equirectangular_projection>`_.
+In this projection, square patches in angular dispersity become wedge-shaped
+patches on the sphere. To correct for the changing point density, there is a
+scale factor of $\sin(\Delta\theta)$ that applies to each point in the
+integral. This is not enough, though. Consider a shape which is tumbling
+freely around the $b$ axis, with $\Delta\theta$ uniform in $[-180, 180]$. At
+$\pm 90$, all points in $\Delta\phi$ map to the pole, so the jitter will have
+a distinct angular preference. If the spin axis is along the beam (which
+will be the case for $\theta=90$ and $\Psi=90$) the scattering pattern
+should be circularly symmetric, but it will go to zero at $q_x = 0$ due to the
+$\sin(\Delta\theta)$ correction. This problem does not appear for a shape
+that is tumbling freely around the $a$ axis, with $\Delta\phi$ uniform in
+$[-180, 180]$, so swap the $a$ and $b$ axes so $\Delta\theta < \Delta\phi$
+and adjust $\Psi$ by 90. This works with the current sasmodels shapes due to
+symmetry.
+
+Alternative projections were considered.
+The `sinusoidal projection <https://en.wikipedia.org/wiki/Sinusoidal_projection>`_
+works by scaling $\Delta\phi$ as $\Delta\theta$ increases, and dropping those
+points outside $[-180, 180]$. The distortions are a little less for middle
+ranges of $\Delta\theta$, but they are still severe for large $\Delta\theta$
+and the model is much harder to explain.
+The `azimuthal equidistance projection <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>`_
+also improves on the equirectangular projection by extending the range of
+reasonable values for the $\Delta\theta$ range, with $\Delta\phi$ forming a
+wedge that cuts to the opposite side of the sphere rather than cutting to the
+pole. This projection has the nice property that distance from the center are
+preserved, and that $\Delta\theta$ and $\Delta\phi$ act the same.
+The `azimuthal equal area projection <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_
+is like the azimuthal equidistance projection, but it preserves area instead
+of distance. It also has the same behaviour for $\Delta\theta$ and $\Delta\phi$.
+The `Guyou projection <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>`_
+has an excellent balance with reasonable distortion in both $\Delta\theta$
+and $\Delta\phi$, as well as preserving small patches. However, it requires
+considerably more computational overhead, and we have not yet derived the
+formula for the distortion correction, measuring the degree of stretch at
+the point $(\Delta\theta, \Delta\phi)$ on the map.
 
 .. note::
-    Note that the form factors for oriented particles are also performing
-    numerical integrations over one or more variables, so care should be taken,
-    especially with very large particles or more extreme aspect ratios. In such 
-    cases results may not be accurate, particularly at very high Q, unless the model
-    has been specifically coded to use limiting forms of the scattering equations.
-    
-    For best numerical results keep the $\theta$ distribution narrower than the $\phi$ 
-    distribution. Thus for asymmetric particles, such as elliptical_cylinder, you may 
-    need to reorder the sizes of the three axes to acheive the desired result. 
-    This is due to the issues of mapping a rectangular distribution onto the 
-    surface of a sphere.
+    Note that the form factors for oriented particles are performing
+    numerical integrations over one or more variables, so care should be
+    taken, especially with very large particles or more extreme aspect
+    ratios. In such cases results may not be accurate, particularly at very
+    high Q, unless the model has been specifically coded to use limiting
+    forms of the scattering equations.
 
-Users can experiment with the values of *Npts* and *Nsigs*, the number of steps 
-used in the integration and the range spanned in number of standard deviations. 
-The standard deviation is entered in units of degrees. For a "rectangular" 
-distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations. 
-The new "uniform" distribution avoids this by letting you directly specify the 
+    For best numerical results keep the $\theta$ distribution narrower than
+    the $\phi$ distribution. Thus for asymmetric particles, such as
+    elliptical_cylinder, you may need to reorder the sizes of the three axes
+    to acheive the desired result. This is due to the issues of mapping a
+    rectanglar distribution onto the surface of a sphere.
+
+Users can experiment with the values of *Npts* and *Nsigs*, the number of steps
+used in the integration and the range spanned in number of standard deviations.
+The standard deviation is entered in units of degrees. For a "rectangular"
+distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations.
+The new "uniform" distribution avoids this by letting you directly specify the
 half width.
 
-The angular distributions will be truncated outside of the range -180 to +180 
-degrees, so beware of using saying a broad Gaussian distribution with large value
-of *Nsigs*, as the array of *Npts* may be truncated to many fewer points than would 
-give a good integration,as well as becoming rather meaningless. (At some point 
-in the future the actual polydispersity arrays may be made available to the user 
-for inspection.)
+The angular distributions may be truncated outside of the range -180 to +180
+degrees, so beware of using saying a broad Gaussian distribution with large
+value of *Nsigs*, as the array of *Npts* may be truncated to many fewer
+points than would give a good integration,as well as becoming rather
+meaningless. (At some point in the future the actual dispersion arrays may be
+made available to the user for inspection.)
 
 Some more detailed technical notes are provided in the developer section of
 this manual :ref:`orientation_developer` .
 
+This definition of orientation is new to SasView 4.2.  In earlier versions,
+the orientation distribution appeared as a distribution of view angles.
+This led to strange effects when $c$ was aligned with $z$, where changes
+to the $\phi$ angle served only to rotate the shape about $c$, rather than
+having a consistent interpretation as the pitch of the shape relative to
+the flow field defining the reference orientation.  Prior to SasView 4.1,
+the reference orientation was defined using a Tait-Bryan convention, making
+it difficult to control.  Now, rotation in $\theta$ modifies the spacings
+in the refraction pattern, and rotation in $\phi$ rotates it in the detector
+plane.
+
+
 *Document History*
 
 | 2017-11-06 Richard Heenan
+| 2017-12-20 Paul Kienzle
Index: doc/guide/plugin.rst
===================================================================
--- doc/guide/plugin.rst	(revision 3048ec65c991d5665c396efe7de93e067fce4707)
+++ doc/guide/plugin.rst	(revision 7e6bc45eddf06819ea47d06f521d0c27c582996e)
@@ -292,6 +292,7 @@
 **Note: The order of the parameters in the definition will be the order of the
 parameters in the user interface and the order of the parameters in Iq(),
-Iqxy() and form_volume(). And** *scale* **and** *background* **parameters are
-implicit to all models, so they do not need to be included in the parameter table.**
+Iqac(), Iqabc() and form_volume(). And** *scale* **and** *background*
+**parameters are implicit to all models, so they do not need to be included
+in the parameter table.**
 
 - **"name"** is the name of the parameter shown on the FitPage.
@@ -362,9 +363,11 @@
     scattered intensity.
 
-  - "volume" parameters are passed to Iq(), Iqxy(), and form_volume(), and
-    have polydispersity loops generated automatically.
-
-  - "orientation" parameters are only passed to Iqxy(), and have angular
-    dispersion.
+  - "volume" parameters are passed to Iq(), Iqac(), Iqabc() and form_volume(),
+    and have polydispersity loops generated automatically.
+
+  - "orientation" parameters are not passed, but instead are combined with
+    orientation dispersity to translate *qx* and *qy* to *qa*, *qb* and *qc*.
+    These parameters should appear at the end of the table with the specific
+    names *theta*, *phi* and for asymmetric shapes *psi*, in that order.
 
 Some models will have integer parameters, such as number of pearls in the
@@ -419,7 +422,5 @@
 That is, the individual models do not need to include polydispersity
 calculations, but instead rely on numerical integration to compute the
-appropriately smeared pattern.   Angular dispersion values over polar angle
-$\theta$ requires an additional $\cos \theta$ weighting due to decreased
-arc length for the equatorial angle $\phi$ with increasing latitude.
+appropriately smeared pattern.
 
 Python Models
@@ -468,8 +469,4 @@
 barbell).  If I(q; pars) is NaN for any $q$, then those parameters will be
 ignored, and not included in the calculation of the weighted polydispersity.
-
-Similar to *Iq*, you can define *Iqxy(qx, qy, par1, par2, ...)* where the
-parameter list includes any orientation parameters.  If *Iqxy* is not defined,
-then it will default to *Iqxy = Iq(sqrt(qx**2+qy**2), par1, par2, ...)*.
 
 Models should define *form_volume(par1, par2, ...)* where the parameter
@@ -497,16 +494,16 @@
     }
 
-*Iqxy* is similar to *Iq*, except it uses parameters *qx, qy* instead of *q*,
-and it includes orientation parameters.
-
 *form_volume* defines the volume of the shape. As in python models, it
 includes only the volume parameters.
 
-*Iqxy* will default to *Iq(sqrt(qx**2 + qy**2), par1, ...)* and
-*form_volume* will default to 1.0.
-
 **source=['fn.c', ...]** includes the listed C source files in the
-program before *Iq* and *Iqxy* are defined. This allows you to extend the
-library of C functions available to your model.
+program before *Iq* and *form_volume* are defined. This allows you to
+extend the library of C functions available to your model.
+
+*c_code* includes arbitrary C code into your kernel, which can be
+handy for defining helper functions for *Iq* and *form_volume*. Note that
+you can put the full function definition for *Iq* and *form_volume*
+(include function declaration) into *c_code* as well, or put them into an
+external C file and add that file to the list of sources.
 
 Models are defined using double precision declarations for the
@@ -533,4 +530,123 @@
     #define INVALID(v) (v.bell_radius < v.radius)
 
+The INVALID define can go into *Iq*, or *c_code*, or an external C file
+listed in *source*.
+
+Oriented Shapes
+...............
+
+If the scattering is dependent on the orientation of the shape, then you
+will need to include *orientation* parameters *theta*, *phi* and *psi*
+at the end of the parameter table.  As described in the section
+:ref:`orientation`, the individual $(q_x, q_y)$ points on the detector will
+be rotated into $(q_a, q_b, q_c)$ points relative to the sample in its
+canonical orientation with $a$-$b$-$c$ aligned with $x$-$y$-$z$ in the
+laboratory frame and beam travelling along $-z$.
+
+The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where
+*par1*, etc. are the parameters to the model.  If the shape is rotationally
+symmetric about *c* then *psi* is not needed, and the model is called
+as *Iqac(qab, qc, par1, par2, ...)*.  In either case, the orientation
+parameters are not included in the function call.
+
+For 1D oriented shapes, an integral over all angles is usually needed for
+the *Iq* function. Given symmetry and the substitution $u = \cos(\alpha)$,
+$du = -\sin(\alpha)\,d\alpha$ this becomes
+
+.. math::
+
+    I(q) &= \frac{1}{4\pi} \int_{-\pi/2}^{pi/2} \int_{-pi}^{pi}
+            F(q_a, q_b, q_c)^2 \sin(\alpha)\,d\beta\,d\alpha \\
+        &= \frac{8}{4\pi} \int_{0}^{pi/2} \int_{0}^{\pi/2}
+            F^2 \sin(\alpha)\,d\beta\,d\alpha \\
+        &= \frac{8}{4\pi} \int_1^0 \int_{0}^{\pi/2} - F^2 \,d\beta\,du \\
+        &= \frac{8}{4\pi} \int_0^1 \int_{0}^{\pi/2} F^2 \,d\beta\,du
+
+for
+
+.. math::
+
+    q_a &= q \sin(\alpha)\sin(\beta) = q \sqrt{1-u^2} \sin(\beta) \\
+    q_b &= q \sin(\alpha)\cos(\beta) = q \sqrt{1-u^2} \cos(\beta) \\
+    q_c &= q \cos(\alpha) = q u
+
+Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the
+numerical integration is then::
+
+    double outer_sum = 0.0;
+    for (int i = 0; i < GAUSS_N; i++) {
+        const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5;
+        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
+        const double qc = cos_alpha * q;
+        double inner_sum = 0.0;
+        for (int j = 0; j < GAUSS_N; j++) {
+            const double beta = M_PI_4 * GAUSS_Z[j] + M_PI_4;
+            double sin_beta, cos_beta;
+            SINCOS(beta, sin_beta, cos_beta);
+            const double qa = sin_alpha * sin_beta * q;
+            const double qb = sin_alpha * cos_beta * q;
+            const double form = Fq(qa, qb, qc, ...);
+            inner_sum += GAUSS_W[j] * form * form;
+        }
+        outer_sum += GAUSS_W[i] * inner_sum;
+    }
+    outer_sum *= 0.25; // = 8/(4 pi) * outer_sum * (pi/2) / 4
+
+The *z* values for the Gauss-Legendre integration extends from -1 to 1, so
+the double sum of *w[i]w[j]* explains the factor of 4.  Correcting for the
+average *dz[i]dz[j]* gives $(1-0) \cdot (\pi/2-0) = \pi/2$.  The $8/(4 \pi)$
+factor comes from the integral over the quadrant.  With less symmetry (eg.,
+in the bcc and fcc paracrystal models), then an integral over the entire
+sphere may be necessary.
+
+For simpler models which are rotationally symmetric a single integral
+suffices:
+
+.. math::
+
+    I(q) &= \frac{1}{\pi}\int_{-\pi/2}^{\pi/2}
+            F(q_{ab}, q_c)^2 \sin(\alpha)\,d\alpha/\pi \\
+        &= \frac{2}{\pi} \int_0^1 F^2\,du
+
+for
+
+.. math::
+
+    q_{ab} &= q \sin(\alpha) = q \sqrt{1 - u^2} \\
+    q_c &= q \cos(\alpha) = q u
+
+
+with integration loop::
+
+    double sum = 0.0;
+    for (int i = 0; i < GAUSS_N; i++) {
+        const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5;
+        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
+        const double qab = sin_alpha * q;
+        const double qc = cos_alpha * q;
+        const double form = Fq(qab, qc, ...);
+        sum += GAUSS_W[j] * form * form;
+    }
+    sum *= 0.5; // = 2/pi * sum * (pi/2) / 2
+
+Magnetism
+.........
+
+Magnetism is supported automatically for all shapes by modifying the
+effective SLD of particle according to the Halpern-Johnson vector
+describing the interaction between neutron spin and magnetic field.  All
+parameters marked as type *sld* in the parameter table are treated as
+possibly magnetic particles with magnitude *M0* and direction
+*mtheta* and *mphi*.  Polarization parameters are also provided
+automatically for magnetic models to set the spin state of the measurement.
+
+For more complicated systems where magnetism is not uniform throughout
+the individual particles, you will need to write your own models.
+You should not mark the nuclear sld as type *sld*, but instead leave
+them unmarked and provide your own magnetism and polarization parameters.
+For 2D measurements you will need $(q_x, q_y)$ values for the measurement
+to compute the proper magnetism and orientation, which you can implement
+using *Iqxy(qx, qy, par1, par2, ...)*.
+
 Special Functions
 .................
@@ -543,8 +659,8 @@
     M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E:
         $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$
-    exp, log, pow(x,y), expm1, sqrt:
-        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$.
-        The function expm1(x) is accurate across all $x$, including $x$
-        very close to zero.
+    exp, log, pow(x,y), expm1, log1p, sqrt, cbrt:
+        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$,
+        $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x)
+        are accurate across all $x$, including $x$ very close to zero.
     sin, cos, tan, asin, acos, atan:
         Trigonometry functions and inverses, operating on radians.
@@ -557,10 +673,11 @@
         atan(y/x) would return a value in quadrant I. Similarly for
         quadrants II and IV when $x$ and $y$ have opposite sign.
-    fmin(x,y), fmax(x,y), trunc, rint:
+    fabs(x), fmin(x,y), fmax(x,y), trunc, rint:
         Floating point functions.  rint(x) returns the nearest integer.
     NAN:
         NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that
         you cannot use :code:`x == NAN` to test for NaN values since that
-        will always return false.  NAN does not equal NAN!
+        will always return false.  NAN does not equal NAN!  The alternative,
+        :code:`x != x` may fail if the compiler optimizes the test away.
     INFINITY:
         $\infty, 1/0$.  Use isinf(x) to test for infinity, or isfinite(x)
@@ -734,4 +851,7 @@
         Similar arrays are available in :code:`gauss20.c` for 20-point
         quadrature and in :code:`gauss150.c` for 150-point quadrature.
+        The macros :code:`GAUSS_N`, :code:`GAUSS_Z` and :code:`GAUSS_W` are
+        defined so that you can change the order of the integration by
+        selecting an different source without touching the C code.
 
         :code:`source = ["lib/gauss76.c", ...]`
@@ -791,13 +911,4 @@
 can be a model that runs 1000x faster on a good card.  Even your laptop may
 show a 50x improvement or more over the equivalent pure python model.
-
-External C Models
-.................
-
-External C models are very much like embedded C models, except that
-*Iq*, *Iqxy* and *form_volume* are defined in an external source file
-loaded using the *source=[...]* statement. You need to supply the function
-declarations for each of these that you need instead of building them
-automatically from the parameter table.
 
 
@@ -1002,5 +1113,5 @@
 variable name *Rg* for example because $R_g$ is the right name for the model
 parameter then ignore the lint errors.  Also, ignore *missing-docstring*
-for standard model functions *Iq*, *Iqxy*, etc.
+for standard model functions *Iq*, *Iqac*, etc.
 
 We will have delinting sessions at the SasView Code Camps, where we can
