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 d0dc9a3eb51c366ad12aca86498d6593a798252a)
+++ 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
@@ -532,4 +529,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
@@ -796,13 +912,4 @@
 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.
-
 
 .. _Form_Factors:
@@ -1006,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
Index: explore/jitter.py
===================================================================
--- explore/jitter.py	(revision ff104792a8c7664b8c2b98af11a2db58f784b8d1)
+++ explore/jitter.py	(revision 8cfb486c14c2f481e9135c412b0ed9506517c141)
@@ -165,4 +165,5 @@
     # constants in kernel_iq.c
     'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance',
+    'azimuthal_equal_area',
 ]
 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian',
Index: explore/precision.py
===================================================================
--- explore/precision.py	(revision 2a602c705ae8385314cccd395073b6945524dbb0)
+++ explore/precision.py	(revision 2a7e20e9c4ba72da5c6088e743a64f3345ee554d)
@@ -345,5 +345,5 @@
 )
 add_function(
-    name="(1/2+(1-cos(x))/x^2-sin(x)/x)/x",
+    name="(1/2-sin(x)/x+(1-cos(x))/x^2)/x",
     mp_function=lambda x: (0.5 - mp.sin(x)/x + (1-mp.cos(x))/(x*x))/x,
     np_function=lambda x: (0.5 - np.sin(x)/x + (1-np.cos(x))/(x*x))/x,
@@ -609,5 +609,5 @@
     names = ", ".join(sorted(ALL_FUNCTIONS))
     print("""\
-usage: precision.py [-f/a/r] [-x<range>] name...
+usage: precision.py [-f/a/r] [-x<range>] "name" ...
 where
     -f indicates that the function value should be plotted,
@@ -620,5 +620,5 @@
       zoom indicates linear stepping in [1000, 1010]
       neg indicates linear stepping in [-100.1, 100.1]
-and name is "all [first]" or one of:
+and name is "all" or one of:
     """+names)
     sys.exit(1)
Index: explore/realspace.py
===================================================================
--- explore/realspace.py	(revision a1c32c237d9d08ea8ba721dbaf09ecf6b862c812)
+++ explore/realspace.py	(revision 8fb2a944dfdac4468e098855656ee605177840c0)
@@ -44,4 +44,28 @@
     """
     return Rx(phi)*Ry(theta)*Rz(psi)
+
+def apply_view(points, view):
+    """
+    Apply the view transform (theta, phi, psi) to a set of points.
+
+    Points are stored in a 3 x n numpy array.
+
+    View angles are in degrees.
+    """
+    theta, phi, psi = view
+    return np.asarray((Rz(phi)*Ry(theta)*Rz(psi))*np.matrix(points.T)).T
+
+
+def invert_view(qx, qy, view):
+    """
+    Return (qa, qb, qc) for the (theta, phi, psi) view angle at detector
+    pixel (qx, qy).
+
+    View angles are in degrees.
+    """
+    theta, phi, psi = view
+    q = np.vstack((qx.flatten(), qy.flatten(), 0*qx.flatten()))
+    return np.asarray((Rz(-psi)*Ry(-theta)*Rz(-phi))*np.matrix(q))
+
 
 class Shape:
@@ -219,4 +243,32 @@
         values = self.value.repeat(points.shape[0])
         return values, self._adjust(points)
+
+NUMBA = False
+if NUMBA:
+    from numba import njit
+    @njit("f8[:](f8[:],f8[:],f8[:],f8[:],f8[:],f8[:],f8[:])")
+    def _Iqxy(values, x, y, z, qa, qb, qc):
+        Iq = np.zeros_like(qa)
+        for j in range(len(Iq)):
+            total = 0. + 0j
+            for k in range(len(Iq)):
+                total += values[k]*np.exp(1j*(qa[j]*x[k] + qb[j]*y[k] + qc[j]*z[k]))
+            Iq[j] = abs(total)**2
+        return Iq
+
+def calc_Iqxy(qx, qy, rho, points, volume=1.0, view=(0, 0, 0)):
+    qx, qy = np.broadcast_arrays(qx, qy)
+    qa, qb, qc = invert_view(qx, qy, view)
+    rho, volume = np.broadcast_arrays(rho, volume)
+    values = rho*volume
+    x, y, z = points.T
+
+    # I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r)
+    if NUMBA:
+        Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten())
+    else:
+        Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2
+              for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)]
+    return np.asarray(Iq).reshape(qx.shape) / np.sum(volume)
 
 def _calc_Pr_nonuniform(r, rho, points):
@@ -239,23 +291,17 @@
             print("processing %d of %d"%(k, len(rho)-1))
     Pr = extended_Pr[1:-1]
-    return Pr / Pr.max()
-
-def calc_Pr(r, rho, points):
-    # P(r) with uniform steps in r is 3x faster; check if we are uniform
-    # before continuing
-    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01:
-        return _calc_Pr_nonuniform(r, rho, points)
-
+    return Pr
+
+def _calc_Pr_uniform(r, rho, points):
     # Make Pr a little be bigger than necessary so that only distances
     # min < d < max end up in Pr
-    n_max = len(r)
+    dr, n_max = r[0], len(r)
     extended_Pr = np.zeros(n_max+1, 'd')
     t0 = time.clock()
     t_next = t0 + 3
-    row_zero = np.zeros(len(rho), 'i')
     for k, rho_k in enumerate(rho[:-1]):
         distances = sqrt(np.sum((points[k] - points[k+1:])**2, axis=1))
         weights = rho_k * rho[k+1:]
-        index = np.minimum(np.asarray(distances/r[0], 'i'), n_max)
+        index = np.minimum(np.asarray(distances/dr, 'i'), n_max)
         # Note: indices may be duplicated, so "Pr[index] += w" will not work!!
         extended_Pr += np.bincount(index, weights, n_max+1)
@@ -269,5 +315,5 @@
     # we are only accumulating the upper triangular distances.
     #Pr = Pr * 2 / len(rho)**2
-    return Pr / Pr.max()
+    return Pr
 
     # Can get an additional 2x by going to C.  Cuda/OpenCL will allow even
@@ -291,4 +337,38 @@
 }
 """
+
+if NUMBA:
+    @njit("f8[:](f8[:], f8[:], f8[:,:])")
+    def _calc_Pr_uniform_jit(r, rho, points):
+        dr = r[0]
+        n_max = len(r)
+        Pr = np.zeros_like(r)
+        for j in range(len(rho) - 1):
+            x, y, z = points[j, 0], points[j, 1], points[j, 2]
+            for k in range(j+1, len(rho)):
+                distance = sqrt((x - points[k, 0])**2
+                                + (y - points[k, 1])**2
+                                + (z - points[k, 2])**2)
+                index = int(distance/dr)
+                if index < n_max:
+                    Pr[index] += rho[j] * rho[k]
+        # Make Pr independent of sampling density.  The factor of 2 comes because
+        # we are only accumulating the upper triangular distances.
+        #Pr = Pr * 2 / len(rho)**2
+        return Pr
+
+
+def calc_Pr(r, rho, points):
+    # P(r) with uniform steps in r is 3x faster; check if we are uniform
+    # before continuing
+    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01:
+        Pr = _calc_Pr_nonuniform(r, rho, points)
+    else:
+        if NUMBA:
+            Pr = _calc_Pr_uniform_jit(r, rho, points)
+        else:
+            Pr = _calc_Pr_uniform(r, rho, points)
+    return Pr / Pr.max()
+
 
 def j0(x):
@@ -333,4 +413,18 @@
         plt.legend()
 
+def plot_calc_2d(qx, qy, Iqxy, theory=None):
+    import matplotlib.pyplot as plt
+    qx, qy = bin_edges(qx), bin_edges(qy)
+    #qx, qy = np.meshgrid(qx, qy)
+    if theory is not None:
+        plt.subplot(121)
+    plt.pcolormesh(qx, qy, np.log10(Iqxy))
+    plt.xlabel('qx (1/A)')
+    plt.ylabel('qy (1/A)')
+    if theory is not None:
+        plt.subplot(122)
+        plt.pcolormesh(qx, qy, np.log10(theory))
+        plt.xlabel('qx (1/A)')
+
 def plot_points(rho, points):
     import mpl_toolkits.mplot3d
@@ -343,9 +437,56 @@
         pass
     n = len(points)
-    index = np.random.choice(n, size=1000) if n > 1000 else slice(None, None)
+    #print("len points", n)
+    index = np.random.choice(n, size=500) if n > 500 else slice(None, None)
     ax.scatter(points[index, 0], points[index, 1], points[index, 2], c=rho[index])
     #low, high = points.min(axis=0), points.max(axis=0)
     #ax.axis([low[0], high[0], low[1], high[1], low[2], high[2]])
     ax.autoscale(True)
+
+def check_shape(shape, fn=None):
+    rho_solvent = 0
+    q = np.logspace(-3, 0, 200)
+    r = shape.r_bins(q, r_step=0.01)
+    sampling_density = 6*5000 / shape.volume()
+    rho, points = shape.sample(sampling_density)
+    t0 = time.time()
+    Pr = calc_Pr(r, rho-rho_solvent, points)
+    print("calc Pr time", time.time() - t0)
+    Iq = calc_Iq(q, r, Pr)
+    theory = (q, fn(q)) if fn is not None else None
+
+    import pylab
+    #plot_points(rho, points); pylab.figure()
+    plot_calc(r, Pr, q, Iq, theory=theory)
+    pylab.show()
+
+def check_shape_2d(shape, fn=None, view=(0, 0, 0)):
+    rho_solvent = 0
+    nq, qmax = 100, 1.0
+    qx = np.linspace(0.0, qmax, nq)
+    qy = np.linspace(0.0, qmax, nq)
+    Qx, Qy = np.meshgrid(qx, qy)
+    sampling_density = 50000 / shape.volume()
+    #t0 = time.time()
+    rho, points = shape.sample(sampling_density)
+    #print("sample time", time.time() - t0)
+    t0 = time.time()
+    Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view)
+    print("calc time", time.time() - t0)
+    theory = fn(Qx, Qy) if fn is not None else None
+    Iqxy += 0.001 * Iqxy.max()
+    if theory is not None:
+        theory += 0.001 * theory.max()
+
+    import pylab
+    #plot_points(rho, points); pylab.figure()
+    plot_calc_2d(qx, qy, Iqxy, theory=theory)
+    pylab.show()
+
+def sas_sinx_x(x):
+    with np.errstate(all='ignore'):
+        retvalue = sin(x)/x
+    retvalue[x == 0.] = 1.
+    return retvalue
 
 def sas_2J1x_x(x):
@@ -373,4 +514,11 @@
     Iq = Iq/Iq[0]
     return Iq
+
+def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)):
+    qa, qb, qc = invert_view(qx, qy, view)
+    qab = np.sqrt(qa**2 + qb**2)
+    Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2)
+    Iq = Fq**2
+    return Iq.reshape(qx.shape)
 
 def sphere_Iq(q, radius):
@@ -415,18 +563,24 @@
     return Iq/Iq[0]
 
-def check_shape(shape, fn=None):
-    rho_solvent = 0
-    q = np.logspace(-3, 0, 200)
-    r = shape.r_bins(q, r_step=0.01)
-    sampling_density = 15000 / shape.volume()
-    rho, points = shape.sample(sampling_density)
-    Pr = calc_Pr(r, rho-rho_solvent, points)
-    Iq = calc_Iq(q, r, Pr)
-    theory = (q, fn(q)) if fn is not None else None
-
-    import pylab
-    #plot_points(rho, points); pylab.figure()
-    plot_calc(r, Pr, q, Iq, theory=theory)
-    pylab.show()
+def csbox_Iqxy(qx, qy, a, b, c, da, db, dc, slda, sldb, sldc, sld_core, view=(0,0,0)):
+    qa, qb, qc = invert_view(qx, qy, view)
+
+    sld_solvent = 0
+    overlapping = False
+    dr0 = sld_core - sld_solvent
+    drA, drB, drC = slda-sld_solvent, sldb-sld_solvent, sldc-sld_solvent
+    tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc
+    siA = a*sas_sinx_x(a*qa/2)
+    siB = b*sas_sinx_x(b*qb/2)
+    siC = c*sas_sinx_x(c*qc/2)
+    siAt = tA*sas_sinx_x(tA*qa/2)
+    siBt = tB*sas_sinx_x(tB*qb/2)
+    siCt = tC*sas_sinx_x(tC*qc/2)
+    Fq = (dr0*siA*siB*siC
+          + drA*(siAt-siA)*siB*siC
+          + drB*siA*(siBt-siB)*siC
+          + drC*siA*siB*(siCt-siC))
+    Iq = Fq**2
+    return Iq.reshape(qx.shape)
 
 def check_cylinder(radius=25, length=125, rho=2.):
@@ -434,4 +588,31 @@
     fn = lambda q: cylinder_Iq(q, radius, length)
     check_shape(shape, fn)
+
+def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)):
+    shape = EllipticalCylinder(radius, radius, length, rho)
+    fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view)
+    check_shape_2d(shape, fn, view=view)
+
+def check_cylinder_2d_lattice(radius=25, length=125, rho=2.,
+                              view=(0, 0, 0)):
+    nx, dx = 1, 2*radius
+    ny, dy = 30, 2*radius
+    nz, dz = 30, length
+    dx, dy, dz = 2*dx, 2*dy, 2*dz
+    def center(*args):
+        sigma = 0.333
+        space = 2
+        return [(space*n+np.random.randn()*sigma)*x for n, x in args]
+    shapes = [EllipticalCylinder(radius, radius, length, rho,
+                                 #center=(ix*dx, iy*dy, iz*dz)
+                                 orientation=np.random.randn(3)*0,
+                                 center=center((ix, dx), (iy, dy), (iz, dz))
+                                )
+              for ix in range(nx)
+              for iy in range(ny)
+              for iz in range(nz)]
+    shape = Composite(shapes)
+    fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view)
+    check_shape_2d(shape, fn, view=view)
 
 def check_sphere(radius=125, rho=2):
@@ -449,9 +630,18 @@
     side_c2 = copy(side_c).shift(0, 0, -c-dc)
     shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2))
-    fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core)
-    check_shape(shape, fn)
+    def fn(q):
+        return csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core)
+    #check_shape(shape, fn)
+
+    view = (20, 30, 40)
+    def fn_xy(qx, qy):
+        return csbox_Iqxy(qx, qy, a, b, c, da, db, dc,
+                          slda, sldb, sldc, sld_core, view=view)
+    check_shape_2d(shape, fn_xy, view=view)
 
 if __name__ == "__main__":
     check_cylinder(radius=10, length=40)
+    #check_cylinder_2d(radius=10, length=40, view=(90,30,0))
+    #check_cylinder_2d_lattice(radius=10, length=50, view=(90,30,0))
     #check_sphere()
     #check_csbox()
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/compare.py	(revision 2a7e20e9c4ba72da5c6088e743a64f3345ee554d)
@@ -92,4 +92,5 @@
     -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
     -neval=1 sets the number of evals for more accurate timing
+    -ngauss=0 overrides the number of points in the 1-D gaussian quadrature
 
     === precision options ===
Index: sasmodels/details.py
===================================================================
--- sasmodels/details.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/details.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -258,4 +258,7 @@
     # type: (...) -> Sequence[np.ndarray]
     """
+    **Deprecated** Theta weights will be computed in the kernel wrapper if
+    they are needed.
+
     If there is a theta parameter, update the weights of that parameter so that
     the cosine weighting required for polar integration is preserved.
@@ -272,5 +275,4 @@
     Returns updated weights vectors
     """
-    # TODO: explain in a comment why scale and background are missing
     # Apparently the parameters.theta_offset similarly skips scale and
     # and background, so the indexing works out, but they are still shipped
@@ -279,5 +281,4 @@
         index = parameters.theta_offset
         theta = dispersity[index]
-        # TODO: modify the dispersity vector to avoid the theta=-90,90,270,...
         theta_weight = abs(cos(radians(theta)))
         weights = tuple(theta_weight*w if k == index else w
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/generate.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -7,10 +7,17 @@
     particular dimensions averaged over all orientations.
 
-    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form
-    with particular dimensions for a single orientation.
-
-    *Imagnetic(qx, qy, result[], p1, p2, ...)* returns the scattering for the
-    polarized neutron spin states (up-up, up-down, down-up, down-down) for
-    a form with particular dimensions for a single orientation.
+    *Iqac(qab, qc, p1, p2, ...)* returns the scattering at qab, qc
+    for a rotationally symmetric form with particular dimensions.
+    qab, qc are determined from shape orientation and scattering angles.
+    This call is used if the shape has orientation parameters theta and phi.
+
+    *Iqabc(qa, qb, qc, p1, p2, ...)* returns the scattering at qa, qb, qc
+    for a form with particular dimensions.  qa, qb, qc are determined from
+    shape orientation and scattering angles. This call is used if the shape
+    has orientation parameters theta, phi and psi.
+
+    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy.  Use this
+    to create an arbitrary 2D theory function, needed for q-dependent
+    background functions and for models with non-uniform magnetism.
 
     *form_volume(p1, p2, ...)* returns the volume of the form with particular
@@ -31,13 +38,12 @@
 scale and background parameters for each model.
 
-*Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99
-functions written for OpenCL.  All functions need prototype declarations
-even if the are defined before they are used.  OpenCL does not support
-*#include* preprocessor directives, so instead the list of includes needs
-to be given as part of the metadata in the kernel module definition.
-The included files should be listed using a path relative to the kernel
-module, or if using "lib/file.c" if it is one of the standard includes
-provided with the sasmodels source.  The includes need to be listed in
-order so that functions are defined before they are used.
+C code should be stylized C-99 functions written for OpenCL. All functions
+need prototype declarations even if the are defined before they are used.
+Although OpenCL supports *#include* preprocessor directives, the list of
+includes should be given as part of the metadata in the kernel module
+definition. The included files should be listed using a path relative to the
+kernel module, or if using "lib/file.c" if it is one of the standard includes
+provided with the sasmodels source. The includes need to be listed in order
+so that functions are defined before they are used.
 
 Floating point values should be declared as *double*.  For single precision
@@ -107,9 +113,9 @@
     present, the volume ratio is 1.
 
-    *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the
-    C source code for the body of the volume, Iq, and Iqxy functions
+    *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing
+    the C source code for the body of the volume, Iq, and Iqac functions
     respectively.  These can also be defined in the last source file.
 
-    *Iq* and *Iqxy* also be instead be python functions defining the
+    *Iq*, *Iqac*, *Iqabc* also be instead be python functions defining the
     kernel.  If they are marked as *Iq.vectorized = True* then the
     kernel is passed the entire *q* vector at once, otherwise it is
@@ -168,4 +174,5 @@
 from zlib import crc32
 from inspect import currentframe, getframeinfo
+import logging
 
 import numpy as np  # type: ignore
@@ -181,4 +188,6 @@
     pass
 # pylint: enable=unused-import
+
+logger = logging.getLogger(__name__)
 
 # jitter projection to use in the kernel code.  See explore/jitter.py
@@ -626,6 +635,6 @@
 
 """
-def _gen_fn(name, pars, body, filename, line):
-    # type: (str, List[Parameter], str, str, int) -> str
+def _gen_fn(model_info, name, pars):
+    # type: (ModelInfo, str, List[Parameter]) -> str
     """
     Generate a function given pars and body.
@@ -639,7 +648,11 @@
     """
     par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void'
+    body = getattr(model_info, name)
+    filename = model_info.filename
+    # Note: if symbol is defined strangely in the module then default it to 1
+    lineno = model_info.lineno.get(name, 1)
     return _FN_TEMPLATE % {
         'name': name, 'pars': par_decl, 'body': body,
-        'filename': filename.replace('\\', '\\\\'), 'line': line,
+        'filename': filename.replace('\\', '\\\\'), 'line': lineno,
     }
 
@@ -656,34 +669,40 @@
 
 # type in IQXY pattern could be single, float, double, long double, ...
-_IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)",
+_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]",
                            flags=re.MULTILINE)
-def _have_Iqxy(sources):
+def find_xy_mode(source):
     # type: (List[str]) -> bool
     """
-    Return true if any file defines Iqxy.
+    Return the xy mode as qa, qac, qabc or qxy.
 
     Note this is not a C parser, and so can be easily confused by
     non-standard syntax.  Also, it will incorrectly identify the following
-    as having Iqxy::
+    as having 2D models::
 
         /*
-        double Iqxy(qx, qy, ...) { ... fill this in later ... }
+        double Iqac(qab, qc, ...) { ... fill this in later ... }
         */
 
-    If you want to comment out an Iqxy function, use // on the front of the
-    line instead.
-    """
-    for _path, code in sources:
-        if _IQXY_PATTERN.search(code):
-            return True
-    return False
-
-
-def _add_source(source, code, path):
+    If you want to comment out the function, use // on the front of the
+    line::
+
+        /*
+        // double Iqac(qab, qc, ...) { ... fill this in later ... }
+        */
+
+    """
+    for code in source:
+        m = _IQXY_PATTERN.search(code)
+        if m is not None:
+            return m.group('mode')
+    return 'qa'
+
+
+def _add_source(source, code, path, lineno=1):
     """
     Add a file to the list of source code chunks, tagged with path and line.
     """
     path = path.replace('\\', '\\\\')
-    source.append('#line 1 "%s"' % path)
+    source.append('#line %d "%s"' % (lineno, path))
     source.append(code)
 
@@ -716,9 +735,4 @@
     user_code = [(f, open(f).read()) for f in model_sources(model_info)]
 
-    # What kind of 2D model do we need?
-    xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str)
-               else 'qac' if not partable.is_asymmetric
-               else 'qabc')
-
     # Build initial sources
     source = []
@@ -727,19 +741,41 @@
         _add_source(source, code, path)
 
+    if model_info.c_code:
+        _add_source(source, model_info.c_code, model_info.filename,
+                    lineno=model_info.lineno.get('c_code', 1))
+
     # Make parameters for q, qx, qy so that we can use them in declarations
-    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')]
+    q, qx, qy, qab, qa, qb, qc \
+        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()]
     # Generate form_volume function, etc. from body only
     if isinstance(model_info.form_volume, str):
         pars = partable.form_volume_parameters
-        source.append(_gen_fn('form_volume', pars, model_info.form_volume,
-                              model_info.filename, model_info._form_volume_line))
+        source.append(_gen_fn(model_info, 'form_volume', pars))
     if isinstance(model_info.Iq, str):
         pars = [q] + partable.iq_parameters
-        source.append(_gen_fn('Iq', pars, model_info.Iq,
-                              model_info.filename, model_info._Iq_line))
+        source.append(_gen_fn(model_info, 'Iq', pars))
     if isinstance(model_info.Iqxy, str):
-        pars = [qx, qy] + partable.iqxy_parameters
-        source.append(_gen_fn('Iqxy', pars, model_info.Iqxy,
-                              model_info.filename, model_info._Iqxy_line))
+        pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters
+        source.append(_gen_fn(model_info, 'Iqxy', pars))
+    if isinstance(model_info.Iqac, str):
+        pars = [qab, qc] + partable.iq_parameters
+        source.append(_gen_fn(model_info, 'Iqac', pars))
+    if isinstance(model_info.Iqabc, str):
+        pars = [qa, qb, qc] + partable.iq_parameters
+        source.append(_gen_fn(model_info, 'Iqabc', pars))
+
+    # What kind of 2D model do we need?  Is it consistent with the parameters?
+    xy_mode = find_xy_mode(source)
+    if xy_mode == 'qabc' and not partable.is_asymmetric:
+        raise ValueError("asymmetric oriented models need to define Iqabc")
+    elif xy_mode == 'qac' and partable.is_asymmetric:
+        raise ValueError("symmetric oriented models need to define Iqac")
+    elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'):
+        raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode)
+    elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'):
+        if xy_mode == 'qxy':
+            logger.warn("oriented shapes should define Iqac or Iqabc")
+        else:
+            raise ValueError("Expected function Iqac or Iqabc for oriented shape")
 
     # Define the parameter table
@@ -767,14 +803,26 @@
     if xy_mode == 'qabc':
         pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
-        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars
+        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_ABC"
     elif xy_mode == 'qac':
         pars = ",".join(["_qa", "_qc"] + model_refs)
-        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars
+        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_AC"
-    else:  # xy_mode == 'qa'
+    elif xy_mode == 'qa':
         pars = ",".join(["_qa"] + model_refs)
         call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_A"
+    elif xy_mode == 'qxy':
+        orientation_refs = _call_pars("_v.", partable.orientation_parameters)
+        pars = ",".join(["_qx", "_qy"] + model_refs + orientation_refs)
+        call_iqxy = "#define CALL_IQ_XY(_qx,_qy,_v) Iqxy(%s)" % pars
+        clear_iqxy = "#undef CALL_IQ_XY"
+        if partable.orientation_parameters:
+            call_iqxy += "\n#define HAVE_THETA"
+            clear_iqxy += "\n#undef HAVE_THETA"
+        if partable.is_asymmetric:
+            call_iqxy += "\n#define HAVE_PSI"
+            clear_iqxy += "\n#undef HAVE_PSI"
+
 
     magpars = [k-2 for k, p in enumerate(partable.call_parameters)
Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision 8698a0ddea5a155fcc2e3f1085f06a618b46eebf)
+++ sasmodels/kernel_header.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -150,2 +150,69 @@
 inline double cube(double x) { return x*x*x; }
 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; }
+
+// CRUFT: support old style models with orientation received qx, qy and angles
+
+// To rotate from the canonical position to theta, phi, psi, first rotate by
+// psi about the major axis, oriented along z, which is a rotation in the
+// detector plane xy. Next rotate by theta about the y axis, aligning the major
+// axis in the xz plane. Finally, rotate by phi in the detector plane xy.
+// To compute the scattering, undo these rotations in reverse order:
+//     rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi
+// The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit
+// vector in the q direction.
+// To change between counterclockwise and clockwise rotation, change the
+// sign of phi and psi.
+
+#if 1
+//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017
+#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \
+    SINCOS(phi*M_PI_180, sn, cn); \
+    q = sqrt(qx*qx + qy*qy); \
+    cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \
+    sn = sqrt(1 - cn*cn); \
+    } while (0)
+#else
+// SasView 3.x definition of orientation
+#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \
+    SINCOS(theta*M_PI_180, sn, cn); \
+    q = sqrt(qx*qx + qy*qy);\
+    cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \
+    sn = sqrt(1 - cn*cn); \
+    } while (0)
+#endif
+
+#if 1
+#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \
+    q = sqrt(qx*qx + qy*qy); \
+    const double qxhat = qx/q; \
+    const double qyhat = qy/q; \
+    double sin_theta, cos_theta; \
+    double sin_phi, cos_phi; \
+    double sin_psi, cos_psi; \
+    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \
+    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \
+    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \
+    xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \
+         + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \
+    yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \
+         + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \
+    zhat = qxhat*(-sin_theta*cos_phi) \
+         + qyhat*(-sin_theta*sin_phi); \
+    } while (0)
+#else
+// SasView 3.x definition of orientation
+#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \
+    q = sqrt(qx*qx + qy*qy); \
+    const double qxhat = qx/q; \
+    const double qyhat = qy/q; \
+    double sin_theta, cos_theta; \
+    double sin_phi, cos_phi; \
+    double sin_psi, cos_psi; \
+    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \
+    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \
+    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \
+    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \
+    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \
+    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \
+    } while (0)
+#endif
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision ec8d4acf126bf6a5404bc968d4eea25959657b12)
+++ sasmodels/kernel_iq.c	(revision ec8d4acf126bf6a5404bc968d4eea25959657b12)
@@ -31,4 +31,5 @@
 //  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes
 //  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes
+//  CALL_IQ_XY(qx, qy, table) : call the Iqxy function for arbitrary models
 //  INVALID(table) : test if the current point is feesible to calculate.  This
 //      will be defined in the kernel definition file.
@@ -469,13 +470,50 @@
   #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc)
   #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
-#endif
-
-// Doing jitter projection code outside the previous if block so that we don't
-// need to repeat the identical logic in the IQ_AC and IQ_ABC branches.  This
-// will become more important if we implement more projections, or more
-// complicated projections.
-#if defined(CALL_IQ) || defined(CALL_IQ_A)
+#elif defined(CALL_IQ_XY)
+  // direct call to qx,qy calculator
+  double qx, qy;
+  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
+  #define BUILD_ROTATION() do {} while(0)
+  #define APPLY_ROTATION() do {} while(0)
+  #define CALL_KERNEL() CALL_IQ_XY(qx, qy, local_values.table)
+#endif
+
+// Define APPLY_PROJECTION depending on model symmetries. We do this outside
+// the previous if block so that we don't need to repeat the identical
+// logic in the IQ_AC and IQ_ABC branches.  This will become more important
+// if we implement more projections, or more complicated projections.
+#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation
   #define APPLY_PROJECTION() const double weight=weight0
-#else // !spherosymmetric projection
+#elif defined(CALL_IQ_XY) // pass orientation to the model
+  // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc
+  // Need to plug the values for the orientation angles back into parameter
+  // table in case they were overridden by the orientation offset.  This
+  // means that orientation dispersity will not work for these models, but
+  // it was broken anyway, so no matter.  Still want to provide Iqxy in case
+  // the user model wants full control of orientation/magnetism.
+  #if defined(HAVE_PSI)
+    const double theta = values[details->theta_par+2];
+    const double phi = values[details->theta_par+3];
+    const double psi = values[details->theta_par+4];
+    double weight;
+    #define APPLY_PROJECTION() do { \
+      local_values.table.theta = theta; \
+      local_values.table.phi = phi; \
+      local_values.table.psi = psi; \
+      weight=weight0; \
+    } while (0)
+  #elif defined(HAVE_THETA)
+    const double theta = values[details->theta_par+2];
+    const double phi = values[details->theta_par+3];
+    double weight;
+    #define APPLY_PROJECTION() do { \
+      local_values.table.theta = theta; \
+      local_values.table.phi = phi; \
+      weight=weight0; \
+    } while (0)
+  #else
+    #define APPLY_PROJECTION() const double weight=weight0
+  #endif
+#else // apply jitter and view before calling the model
   // Grab the "view" angles (theta, phi, psi) from the initial parameter table.
   const double theta = values[details->theta_par+2];
@@ -488,5 +526,5 @@
   // we go through the mesh.
   double dtheta, dphi, weight;
-  #if PROJECTION == 1
+  #if PROJECTION == 1 // equirectangular
     #define APPLY_PROJECTION() do { \
       dtheta = local_values.table.theta; \
@@ -494,5 +532,5 @@
       weight = fabs(cos(dtheta*M_PI_180)) * weight0; \
     } while (0)
-  #elif PROJECTION == 2
+  #elif PROJECTION == 2 // sinusoidal
     #define APPLY_PROJECTION() do { \
       dtheta = local_values.table.theta; \
@@ -504,5 +542,5 @@
     } while (0)
   #endif
-#endif // !spherosymmetric projection
+#endif // done defining APPLY_PROJECTION
 
 // ** define looping macros **
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/kernelpy.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -26,4 +26,6 @@
 # pylint: enable=unused-import
 
+logger = logging.getLogger(__name__)
+
 class PyModel(KernelModel):
     """
@@ -31,14 +33,13 @@
     """
     def __init__(self, model_info):
-        # Make sure Iq and Iqxy are available and vectorized
+        # Make sure Iq is available and vectorized
         _create_default_functions(model_info)
         self.info = model_info
         self.dtype = np.dtype('d')
-        logging.info("load python model " + self.info.name)
+        logger.info("load python model " + self.info.name)
 
     def make_kernel(self, q_vectors):
         q_input = PyInput(q_vectors, dtype=F64)
-        kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq
-        return PyKernel(kernel, self.info, q_input)
+        return PyKernel(self.info, q_input)
 
     def release(self):
@@ -89,5 +90,5 @@
     Callable SAS kernel.
 
-    *kernel* is the DllKernel object to call.
+    *kernel* is the kernel object to call.
 
     *model_info* is the module information
@@ -104,5 +105,5 @@
     Call :meth:`release` when done with the kernel instance.
     """
-    def __init__(self, kernel, model_info, q_input):
+    def __init__(self, model_info, q_input):
         # type: (callable, ModelInfo, List[np.ndarray]) -> None
         self.dtype = np.dtype('d')
@@ -110,5 +111,4 @@
         self.q_input = q_input
         self.res = np.empty(q_input.nq, q_input.dtype)
-        self.kernel = kernel
         self.dim = '2d' if q_input.is_2d else '1d'
 
@@ -159,6 +159,6 @@
         # Generate a closure which calls the form_volume if it exists.
         form_volume = model_info.form_volume
-        self._volume = ((lambda: form_volume(*volume_args)) if form_volume
-                        else (lambda: 1.0))
+        self._volume = ((lambda: form_volume(*volume_args)) if form_volume else
+                        (lambda: 1.0))
 
     def __call__(self, call_details, values, cutoff, magnetic):
@@ -261,6 +261,7 @@
     any functions that are not already marked as vectorized.
     """
+    # Note: must call create_vector_Iq before create_vector_Iqxy
     _create_vector_Iq(model_info)
-    _create_vector_Iqxy(model_info)  # call create_vector_Iq() first
+    _create_vector_Iqxy(model_info)
 
 
@@ -280,9 +281,10 @@
         model_info.Iq = vector_Iq
 
+
 def _create_vector_Iqxy(model_info):
     """
     Define Iqxy as a vector function if it exists, or default it from Iq().
     """
-    Iq, Iqxy = model_info.Iq, model_info.Iqxy
+    Iqxy = getattr(model_info, 'Iqxy', None)
     if callable(Iqxy):
         if not getattr(Iqxy, 'vectorized', False):
@@ -295,7 +297,8 @@
             vector_Iqxy.vectorized = True
             model_info.Iqxy = vector_Iqxy
-    elif callable(Iq):
+    else:
         #print("defaulting Iqxy")
         # Iq is vectorized because create_vector_Iq was already called.
+        Iq = model_info.Iq
         def default_Iqxy(qx, qy, *args):
             """
Index: sasmodels/modelinfo.py
===================================================================
--- sasmodels/modelinfo.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/modelinfo.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -37,5 +37,5 @@
 
 # assumptions about common parameters exist throughout the code, such as:
-# (1) kernel functions Iq, Iqxy, form_volume, ... don't see them
+# (1) kernel functions Iq, Iqxy, Iqac, Iqabc, form_volume, ... don't see them
 # (2) kernel drivers assume scale is par[0] and background is par[1]
 # (3) mixture models drop the background on components and replace the scale
@@ -256,7 +256,7 @@
 
     *type* indicates how the parameter will be used.  "volume" parameters
-    will be used in all functions.  "orientation" parameters will be used
-    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in
-    *Imagnetic* only.  If *type* is the empty string, the parameter will
+    will be used in all functions.  "orientation" parameters are not passed,
+    but will be used to convert from *qx*, *qy* to *qa*, *qb*, *qc* in calls to
+    *Iqxy* and *Imagnetic*.  If *type* is the empty string, the parameter will
     be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters
     can automatically be promoted to magnetic parameters, each of which
@@ -386,7 +386,4 @@
       with vector parameter p sent as p[].
 
-    * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...)
-      function, with vector parameter p sent as p[].
-
     * *form_volume_parameters* is the list of parameters to the form_volume(...)
       function, with vector parameter p sent as p[].
@@ -443,7 +440,6 @@
         self.iq_parameters = [p for p in self.kernel_parameters
                               if p.type not in ('orientation', 'magnetic')]
-        # note: orientation no longer sent to Iqxy, so its the same as
-        #self.iqxy_parameters = [p for p in self.kernel_parameters
-        #                        if p.type != 'magnetic']
+        self.orientation_parameters = [p for p in self.kernel_parameters
+                                       if p.type == 'orientation']
         self.form_volume_parameters = [p for p in self.kernel_parameters
                                        if p.type == 'volume']
@@ -490,9 +486,15 @@
                 if p.type != 'orientation':
                     raise TypeError("psi must be an orientation parameter")
+            elif p.type == 'orientation':
+                raise TypeError("only theta, phi and psi can be orientation parameters")
         if theta >= 0 and phi >= 0:
+            last_par = len(self.kernel_parameters) - 1
             if phi != theta+1:
                 raise TypeError("phi must follow theta")
             if psi >= 0 and psi != phi+1:
                 raise TypeError("psi must follow phi")
+            #if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par):
+            #    raise TypeError("orientation parameters must appear at the "
+            #                    "end of the parameter table")
         elif theta >= 0 or phi >= 0 or psi >= 0:
             raise TypeError("oriented shapes must have both theta and phi and maybe psi")
@@ -715,4 +717,7 @@
 
 
+#: Set of variables defined in the model that might contain C code
+C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code']
+
 def _find_source_lines(model_info, kernel_module):
     # type: (ModelInfo, ModuleType) -> None
@@ -720,37 +725,32 @@
     Identify the location of the C source inside the model definition file.
 
-    This code runs through the source of the kernel module looking for
-    lines that start with 'Iq', 'Iqxy' or 'form_volume'.  Clearly there are
-    all sorts of reasons why this might not work (e.g., code commented out
-    in a triple-quoted line block, code built using string concatenation,
-    or code defined in the branch of an 'if' block), but it should work
-    properly in the 95% case, and getting the incorrect line number will
-    be harmless.
-    """
-    # Check if we need line numbers at all
-    if callable(model_info.Iq):
-        return None
-
-    if (model_info.Iq is None
-            and model_info.Iqxy is None
-            and model_info.Imagnetic is None
-            and model_info.form_volume is None):
+    This code runs through the source of the kernel module looking for lines
+    that contain C code (because it is a c function definition).  Clearly
+    there are all sorts of reasons why this might not work (e.g., code
+    commented out in a triple-quoted line block, code built using string
+    concatenation, code defined in the branch of an 'if' block, code imported
+    from another file), but it should work properly in the 95% case, and for
+    the remainder, getting the incorrect line number will merely be
+    inconvenient.
+    """
+    # Only need line numbers if we are creating a C module and the C symbols
+    # are defined.
+    if (callable(model_info.Iq)
+            or not any(hasattr(model_info, s) for s in C_SYMBOLS)):
         return
 
-    # find the defintion lines for the different code blocks
+    # load the module source if we can
     try:
         source = inspect.getsource(kernel_module)
     except IOError:
         return
-    for k, v in enumerate(source.split('\n')):
-        if v.startswith('Imagnetic'):
-            model_info._Imagnetic_line = k+1
-        elif v.startswith('Iqxy'):
-            model_info._Iqxy_line = k+1
-        elif v.startswith('Iq'):
-            model_info._Iq_line = k+1
-        elif v.startswith('form_volume'):
-            model_info._form_volume_line = k+1
-
+
+    # look for symbol at the start of the line
+    for lineno, line in enumerate(source.split('\n')):
+        for name in C_SYMBOLS:
+            if line.startswith(name):
+                # Add 1 since some compilers complain about "#line 0"
+                model_info.lineno[name] = lineno + 1
+                break
 
 def make_model_info(kernel_module):
@@ -761,5 +761,5 @@
     Fill in default values for parts of the module that are not provided.
 
-    Note: vectorized Iq and Iqxy functions will be created for python
+    Note: vectorized Iq and Iqac/Iqabc functions will be created for python
     models when the model is first called, not when the model is loaded.
     """
@@ -790,4 +790,5 @@
     info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
     info.source = getattr(kernel_module, 'source', [])
+    info.c_code = getattr(kernel_module, 'c_code', None)
     # TODO: check the structure of the tests
     info.tests = getattr(kernel_module, 'tests', [])
@@ -797,4 +798,6 @@
     info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore
     info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore
+    info.Iqac = getattr(kernel_module, 'Iqac', None) # type: ignore
+    info.Iqabc = getattr(kernel_module, 'Iqabc', None) # type: ignore
     info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore
     info.profile = getattr(kernel_module, 'profile', None) # type: ignore
@@ -811,4 +814,8 @@
     info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore
 
+    if callable(info.Iq) and parameters.has_2d:
+        raise ValueError("oriented python models not supported")
+
+    info.lineno = {}
     _find_source_lines(info, kernel_module)
 
@@ -824,5 +831,5 @@
 
     The structure should be mostly static, other than the delayed definition
-    of *Iq* and *Iqxy* if they need to be defined.
+    of *Iq*, *Iqac* and *Iqabc* if they need to be defined.
     """
     #: Full path to the file defining the kernel, if any.
@@ -906,9 +913,9 @@
     structure_factor = None # type: bool
     #: List of C source files used to define the model.  The source files
-    #: should define the *Iq* function, and possibly *Iqxy*, though a default
-    #: *Iqxy = Iq(sqrt(qx**2+qy**2)* will be created if no *Iqxy* is provided.
-    #: Files containing the most basic functions must appear first in the list,
-    #: followed by the files that use those functions.  Form factors are
-    #: indicated by providing a :attr:`ER` function.
+    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the
+    #: model defines orientation parameters. Files containing the most basic
+    #: functions must appear first in the list, followed by the files that
+    #: use those functions.  Form factors are indicated by providing
+    #: an :attr:`ER` function.
     source = None           # type: List[str]
     #: The set of tests that must pass.  The format of the tests is described
@@ -935,4 +942,8 @@
     #: See :attr:`ER` for details on the parameters.
     VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]]
+    #: Arbitrary C code containing supporting functions, etc., to be inserted
+    #: after everything in source.  This can include Iq and Iqxy functions with
+    #: the full function signature, including all parameters.
+    c_code = None
     #: Returns the form volume for python-based models.  Form volume is needed
     #: for volume normalization in the polydispersity integral.  If no
@@ -955,6 +966,8 @@
     #: include the decimal point. See :mod:`generate` for more details.
     Iq = None               # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
-    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`.
-    Iqxy = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
+    #: Returns *I(qab, qc, a, b, ...)*.  The interface follows :attr:`Iq`.
+    Iqac = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
+    #: Returns *I(qa, qb, qc, a, b, ...)*.  The interface follows :attr:`Iq`.
+    Iqabc = None            # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
     #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`.
     Imagnetic = None        # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
@@ -972,12 +985,6 @@
     #: Returns a random parameter set for the model
     random = None           # type: Optional[Callable[[], Dict[str, float]]]
-
-    # line numbers within the python file for bits of C source, if defined
-    # NB: some compilers fail with a "#line 0" directive, so default to 1.
-    _Imagnetic_line = 1
-    _Iqxy_line = 1
-    _Iq_line = 1
-    _form_volume_line = 1
-
+    #: Line numbers for symbols defining C code
+    lineno = None           # type: Dict[str, int]
 
     def __init__(self):
Index: sasmodels/models/_spherepy.py
===================================================================
--- sasmodels/models/_spherepy.py	(revision ef07e950965796e2a3a8156fcade98ce5e403001)
+++ sasmodels/models/_spherepy.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -88,8 +88,4 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
-def Iqxy(qx, qy, sld, sld_solvent, radius):
-    return Iq(sqrt(qx ** 2 + qy ** 2), sld, sld_solvent, radius)
-Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values
-
 def sesans(z, sld, sld_solvent, radius):
     """
Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/barbell.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -90,5 +90,5 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double sld, double solvent_sld,
     double radius_bell, double radius, double length)
Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/bcc_paracrystal.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -107,5 +107,5 @@
 
 
-static double Iqxy(double qa, double qb, double qc,
+static double Iqabc(double qa, double qb, double qc,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld)
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/capped_cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -115,5 +115,5 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double sld, double solvent_sld, double radius,
     double radius_cap, double length)
Index: sasmodels/models/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/core_shell_bicelle.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -67,5 +67,5 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double radius,
     double thick_rim,
Index: sasmodels/models/core_shell_bicelle_elliptical.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/core_shell_bicelle_elliptical.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -71,5 +71,5 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double r_minor,
     double x_core,
Index: sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -79,5 +79,5 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
           double r_minor,
           double x_core,
@@ -114,3 +114,2 @@
     return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma));
 }
-
Index: sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py	(revision 110f69cfa1143196495c76df1d5cff08c6ce884e)
+++ sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -149,8 +149,8 @@
     ["sld_rim",        "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Cylinder rim scattering length density"],
     ["sld_solvent",    "1e-6/Ang^2", 6, [-inf, inf], "sld",         "Solvent scattering length density"],
+    ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"],
     ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"],
     ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"],
     ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about cylinder axis"],
-    ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"]
     ]
 
Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/core_shell_cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -48,5 +48,5 @@
 
 
-double Iqxy(double qab, double qc,
+double Iqac(double qab, double qc,
     double core_sld,
     double shell_sld,
Index: sasmodels/models/core_shell_ellipsoid.c
===================================================================
--- sasmodels/models/core_shell_ellipsoid.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/core_shell_ellipsoid.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -75,5 +75,5 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double radius_equat_core,
     double x_core,
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision 4493288f50c1e858772d5058d6ca027bfa72f50e)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision e0772316ffb087990c2e1b341c744594e3087097)
@@ -43,6 +43,6 @@
     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
-    // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
-    // Code rewritten (PAK)
+    // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker)
+    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK)
 
     const double half_q = 0.5*q;
@@ -102,5 +102,5 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double core_sld,
     double arim_sld,
@@ -121,6 +121,4 @@
     const double drC = crim_sld-solvent_sld;
 
-    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
-    // the scaling by B.
     const double tA = length_a + 2.0*thick_rim_a;
     const double tB = length_b + 2.0*thick_rim_b;
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision 10ee8389f2d5511546ae78ffddc6d16e36d2a6a8)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision 97be87700f140bc9317b791ad4691cc31fa30981)
@@ -136,6 +136,8 @@
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016
-* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017
-* **Currently Under review by:** Paul Butler
+* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
+* Cross-checked against hollow rectangular prism and rectangular prism for
+  equal thickness overlapping sides, and by Monte Carlo sampling of points
+  within the shape for non-uniform, non-overlapping sides.
 """
 
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -45,5 +45,5 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double sld,
     double solvent_sld,
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/ellipsoid.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -39,5 +39,5 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double sld,
     double sld_solvent,
Index: sasmodels/models/elliptical_cylinder.c
===================================================================
--- sasmodels/models/elliptical_cylinder.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/elliptical_cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -54,5 +54,5 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
      double radius_minor, double r_ratio, double length,
      double sld, double solvent_sld)
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/fcc_paracrystal.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -80,5 +80,5 @@
 
 
-static double Iqxy(double qa, double qb, double qc,
+static double Iqabc(double qa, double qb, double qc,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld)
Index: sasmodels/models/hollow_cylinder.c
===================================================================
--- sasmodels/models/hollow_cylinder.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/hollow_cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -52,5 +52,5 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double radius, double thickness, double length,
     double sld, double solvent_sld)
Index: sasmodels/models/hollow_rectangular_prism.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -85,5 +85,5 @@
 }
 
-double Iqxy(double qa, double qb, double qc,
+double Iqabc(double qa, double qb, double qc,
     double sld,
     double solvent_sld,
Index: sasmodels/models/line.py
===================================================================
--- sasmodels/models/line.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/line.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -57,4 +57,5 @@
 Iq.vectorized = True # Iq accepts an array of q values
 
+
 def Iqxy(qx, qy, *args):
     """
@@ -69,4 +70,12 @@
 
 Iqxy.vectorized = True  # Iqxy accepts an array of qx qy values
+
+# uncomment the following to test Iqxy in C models
+#del Iq, Iqxy
+#c_code = """
+#static double Iq(double q, double b, double m) { return m*q+b; }
+#static double Iqxy(double qx, double qy, double b, double m)
+#{ return (m*qx+b)*(m*qy+b); }
+#"""
 
 def random():
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/parallelepiped.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -53,5 +53,5 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double sld,
     double solvent_sld,
Index: sasmodels/models/rectangular_prism.c
===================================================================
--- sasmodels/models/rectangular_prism.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/rectangular_prism.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -71,5 +71,5 @@
 
 
-double Iqxy(double qa, double qb, double qc,
+double Iqabc(double qa, double qb, double qc,
     double sld,
     double solvent_sld,
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/sc_paracrystal.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -82,5 +82,5 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld)
Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/stacked_disks.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -142,5 +142,5 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double thick_core,
     double thick_layer,
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -46,5 +46,5 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double sld,
     double sld_solvent,
