Index: .gitignore
===================================================================
--- .gitignore	(revision 2badeca29cc1a284f655f99079641fba3e5e5822)
+++ .gitignore	(revision 2badeca29cc1a284f655f99079641fba3e5e5822)
@@ -25,2 +25,3 @@
 /example/Fit_*/
 /example/batch_fit.csv
+/sasmodels/models/lib/gauss*.c
Index: doc/guide/plugin.rst
===================================================================
--- doc/guide/plugin.rst	(revision 0a9fcab08e03b08472a56fca8b7ff2a9d656a43b)
+++ doc/guide/plugin.rst	(revision c654160fb828fbdd130e49cbf46a616ce82875bb)
@@ -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,128 @@
 
     #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.  Shape orientation uses *a*, *b* and *c*
+axes, corresponding to the *x*, *y* and *z* axes in the laboratory coordinate
+system, with *z* along the beam and *x*-*y* in the detector plane, with *x*
+horizontal and *y* vertical.  The *psi* parameter rotates the shape
+about its *c* axis, the *theta* parameter then rotates the *c* axis toward
+the *x* axis of the detector, then *phi* rotates the shape in the detector
+plane.  (Prior to these rotations, orientation dispersity will be applied
+as roll-pitch-yaw, rotating *c*, then *b* then *a* in the shape coordinate
+system.)  A particular *qx*, *qy* point on the detector, then corresponds
+to *qa*, *qb*, *qc* with respect to the shape.
+
+The oriented C model is called as *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 +917,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 +1118,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/asymint.py
===================================================================
--- explore/asymint.py	(revision 1820208ff3a9fba87725a7037153e1cf974a72a9)
+++ explore/asymint.py	(revision a1c32c237d9d08ea8ba721dbaf09ecf6b862c812)
@@ -86,11 +86,44 @@
     a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
     def Fq(qa, qb, qc):
-        siA = env.sas_sinx_x(0.5*a*qa/2)
-        siB = env.sas_sinx_x(0.5*b*qb/2)
-        siC = env.sas_sinx_x(0.5*c*qc/2)
+        siA = env.sas_sinx_x(a*qa/2)
+        siB = env.sas_sinx_x(b*qb/2)
+        siC = env.sas_sinx_x(c*qc/2)
         return siA * siB * siC
     Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c)
     volume = a*b*c
     norm = CONTRAST**2*volume/10000
+    return norm, Fq
+
+def make_core_shell_parallelepiped(a, b, c, da, db, dc, slda, sldb, sldc, env=NPenv):
+    overlapping = False
+    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
+    da, db, dc = env.mpf(da), env.mpf(db), env.mpf(dc)
+    slda, sldb, sldc = env.mpf(slda), env.mpf(sldb), env.mpf(sldc)
+    dr0 = CONTRAST
+    drA, drB, drC = slda-SLD_SOLVENT, sldb-SLD_SOLVENT, sldc-SLD_SOLVENT
+    tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc
+    def Fq(qa, qb, qc):
+        siA = a*env.sas_sinx_x(a*qa/2)
+        siB = b*env.sas_sinx_x(b*qb/2)
+        siC = c*env.sas_sinx_x(c*qc/2)
+        siAt = tA*env.sas_sinx_x(tA*qa/2)
+        siBt = tB*env.sas_sinx_x(tB*qb/2)
+        siCt = tC*env.sas_sinx_x(tC*qc/2)
+        if overlapping:
+            return (dr0*siA*siB*siC
+                    + drA*(siAt-siA)*siB*siC
+                    + drB*siAt*(siBt-siB)*siC
+                    + drC*siAt*siBt*(siCt-siC))
+        else:
+            return (dr0*siA*siB*siC
+                    + drA*(siAt-siA)*siB*siC
+                    + drB*siA*(siBt-siB)*siC
+                    + drC*siA*siB*(siCt-siC))
+    Fq.__doc__ = "core-shell parallelepiped a=%g, b=%g c=%g"%(a, b, c)
+    if overlapping:
+        volume = a*b*c + 2*da*b*c + 2*tA*db*c + 2*tA*tB*dc
+    else:
+        volume = a*b*c + 2*da*b*c + 2*a*db*c + 2*a*b*dc
+    norm = 1/(volume*10000)
     return norm, Fq
 
@@ -184,4 +217,22 @@
     NORM, KERNEL = make_parallelepiped(A, B, C)
     NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv)
+elif shape == 'core_shell_parallelepiped':
+    #A, B, C = 4450, 14000, 47
+    #A, B, C = 445, 140, 47  # integer for the sake of mpf
+    A, B, C = 6800, 114, 1380
+    DA, DB, DC = 2300, 21, 58
+    SLDA, SLDB, SLDC = "5", "-0.3", "11.5"
+    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3
+    #SLD_SOLVENT,CONTRAST = 0, 4
+    if 1: # C shortest
+        B, C = C, B
+        DB, DC = DC, DB
+        SLDB, SLDC = SLDC, SLDB
+    elif 0: # C longest
+        A, C = C, A
+        DA, DC = DC, DA
+        SLDA, SLDC = SLDC, SLDA
+    NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
+    NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv)
 elif shape == 'paracrystal':
     LATTICE = 'bcc'
@@ -342,9 +393,12 @@
     print("gauss-150", *gauss_quad_2d(Q, n=150))
     print("gauss-500", *gauss_quad_2d(Q, n=500))
+    print("gauss-1025", *gauss_quad_2d(Q, n=1025))
+    print("gauss-2049", *gauss_quad_2d(Q, n=2049))
     #gridded_2d(Q, n=2**8+1)
     gridded_2d(Q, n=2**10+1)
-    #gridded_2d(Q, n=2**13+1)
+    #gridded_2d(Q, n=2**12+1)
     #gridded_2d(Q, n=2**15+1)
-    if shape != 'paracrystal':  # adaptive forms are too slow!
+    if shape not in ('paracrystal', 'core_shell_parallelepiped'):
+        # adaptive forms on models for which the calculations are fast enough
         print("dblquad", *scipy_dblquad_2d(Q))
         print("semi-romberg-100", *semi_romberg_2d(Q, n=100))
Index: explore/realspace.py
===================================================================
--- explore/realspace.py	(revision a1c32c237d9d08ea8ba721dbaf09ecf6b862c812)
+++ explore/realspace.py	(revision a1c32c237d9d08ea8ba721dbaf09ecf6b862c812)
@@ -0,0 +1,458 @@
+from __future__ import division, print_function
+
+import time
+from copy import copy
+
+import numpy as np
+from numpy import pi, radians, sin, cos, sqrt
+from numpy.random import poisson, uniform
+from numpy.polynomial.legendre import leggauss
+from scipy.integrate import simps
+from scipy.special import j1 as J1
+
+# Definition of rotation matrices comes from wikipedia:
+#    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
+def Rx(angle):
+    """Construct a matrix to rotate points about *x* by *angle* degrees."""
+    a = radians(angle)
+    R = [[1, 0, 0],
+         [0, +cos(a), -sin(a)],
+         [0, +sin(a), +cos(a)]]
+    return np.matrix(R)
+
+def Ry(angle):
+    """Construct a matrix to rotate points about *y* by *angle* degrees."""
+    a = radians(angle)
+    R = [[+cos(a), 0, +sin(a)],
+         [0, 1, 0],
+         [-sin(a), 0, +cos(a)]]
+    return np.matrix(R)
+
+def Rz(angle):
+    """Construct a matrix to rotate points about *z* by *angle* degrees."""
+    a = radians(angle)
+    R = [[+cos(a), -sin(a), 0],
+         [+sin(a), +cos(a), 0],
+         [0, 0, 1]]
+    return np.matrix(R)
+
+def rotation(theta, phi, psi):
+    """
+    Apply the jitter transform to a set of points.
+
+    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
+    """
+    return Rx(phi)*Ry(theta)*Rz(psi)
+
+class Shape:
+    rotation = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]])
+    center = np.array([0., 0., 0.])[:, None]
+    r_max = None
+
+    def volume(self):
+        # type: () -> float
+        raise NotImplementedError()
+
+    def sample(self, density):
+        # type: (float) -> np.ndarray[N], np.ndarray[N, 3]
+        raise NotImplementedError()
+
+    def rotate(self, theta, phi, psi):
+        self.rotation = rotation(theta, phi, psi) * self.rotation
+        return self
+
+    def shift(self, x, y, z):
+        self.center = self.center + np.array([x, y, z])[:, None]
+        return self
+
+    def _adjust(self, points):
+        points = np.asarray(self.rotation * np.matrix(points.T)) + self.center
+        return points.T
+
+    def r_bins(self, q, over_sampling=1, r_step=0.):
+        r_max = min(2 * pi / q[0], self.r_max)
+        if r_step == 0.:
+            r_step = 2 * pi / q[-1] / over_sampling
+        #r_step = 0.01
+        return np.arange(r_step, r_max, r_step)
+
+class Composite(Shape):
+    def __init__(self, shapes, center=(0, 0, 0), orientation=(0, 0, 0)):
+        self.shapes = shapes
+        self.rotate(*orientation)
+        self.shift(*center)
+
+        # Find the worst case distance between any two points amongst a set
+        # of shapes independent of orientation.  This could easily be a
+        # factor of two worse than necessary, e.g., a pair of thin rods
+        # end-to-end vs the same pair side-by-side.
+        distances = [((s1.r_max + s2.r_max)/2
+                      + sqrt(np.sum((s1.center - s2.center)**2)))
+                     for s1 in shapes
+                     for s2 in shapes]
+        self.r_max = max(distances + [s.r_max for s in shapes])
+
+    def volume(self):
+        return sum(shape.volume() for shape in self.shapes)
+
+    def sample(self, density):
+        values, points = zip(*(shape.sample(density) for shape in self.shapes))
+        return np.hstack(values), self._adjust(np.vstack(points))
+
+class Box(Shape):
+    def __init__(self, a, b, c,
+                 value, center=(0, 0, 0), orientation=(0, 0, 0)):
+        self.value = np.asarray(value)
+        self.rotate(*orientation)
+        self.shift(*center)
+        self.a, self.b, self.c = a, b, c
+        self._scale = np.array([a/2, b/2, c/2])[None, :]
+        self.r_max = sqrt(a**2 + b**2 + c**2)
+
+    def volume(self):
+        return self.a*self.b*self.c
+
+    def sample(self, density):
+        num_points = poisson(density*self.a*self.b*self.c)
+        points = self._scale*uniform(-1, 1, size=(num_points, 3))
+        values = self.value.repeat(points.shape[0])
+        return values, self._adjust(points)
+
+class EllipticalCylinder(Shape):
+    def __init__(self, ra, rb, length,
+                 value, center=(0, 0, 0), orientation=(0, 0, 0)):
+        self.value = np.asarray(value)
+        self.rotate(*orientation)
+        self.shift(*center)
+        self.ra, self.rb, self.length = ra, rb, length
+        self._scale = np.array([ra, rb, length/2])[None, :]
+        self.r_max = sqrt(4*max(ra, rb)**2 + length**2)
+
+    def volume(self):
+        return pi*self.ra*self.rb*self.length
+
+    def sample(self, density):
+        # density of the bounding box
+        num_points = poisson(density*4*self.ra*self.rb*self.length)
+        points = uniform(-1, 1, size=(num_points, 3))
+        radius = points[:, 0]**2 + points[:, 1]**2
+        points = self._scale*points[radius <= 1]
+        values = self.value.repeat(points.shape[0])
+        return values, self._adjust(points)
+
+class TriaxialEllipsoid(Shape):
+    def __init__(self, ra, rb, rc,
+                 value, center=(0, 0, 0), orientation=(0, 0, 0)):
+        self.value = np.asarray(value)
+        self.rotate(*orientation)
+        self.shift(*center)
+        self.ra, self.rb, self.rc = ra, rb, rc
+        self._scale = np.array([ra, rb, rc])[None, :]
+        self.r_max = 2*max(ra, rb, rc)
+
+    def volume(self):
+        return 4*pi/3 * self.ra * self.rb * self.rc
+
+    def sample(self, density):
+        # randomly sample from a box of side length 2*r, excluding anything
+        # not in the ellipsoid
+        num_points = poisson(density*8*self.ra*self.rb*self.rc)
+        points = uniform(-1, 1, size=(num_points, 3))
+        radius = np.sum(points**2, axis=1)
+        points = self._scale*points[radius <= 1]
+        values = self.value.repeat(points.shape[0])
+        return values, self._adjust(points)
+
+class Helix(Shape):
+    def __init__(self, helix_radius, helix_pitch, tube_radius, tube_length,
+                 value, center=(0, 0, 0), orientation=(0, 0, 0)):
+        self.value = np.asarray(value)
+        self.rotate(*orientation)
+        self.shift(*center)
+        self.helix_radius, self.helix_pitch = helix_radius, helix_pitch
+        self.tube_radius, self.tube_length = tube_radius, tube_length
+        helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2)
+        self.r_max = sqrt((2*helix_radius + 2*tube_radius)*2
+                          + (helix_length + 2*tube_radius)**2)
+
+    def volume(self):
+        # small tube radius approximation; for larger tubes need to account
+        # for the fact that the inner length is much shorter than the outer
+        # length
+        return pi*self.tube_radius**2*self.tube_length
+
+    def points(self, density):
+        num_points = poisson(density*4*self.tube_radius**2*self.tube_length)
+        points = uniform(-1, 1, size=(num_points, 3))
+        radius = points[:, 0]**2 + points[:, 1]**2
+        points = points[radius <= 1]
+
+        # Based on math stackexchange answer by Jyrki Lahtonen
+        #     https://math.stackexchange.com/a/461637
+        # with helix along z rather than x [so tuples in answer are (z, x, y)]
+        # and with random points in the cross section (p1, p2) rather than
+        # uniform points on the surface (cos u, sin u).
+        a, R = self.tube_radius, self.helix_radius
+        h = self.helix_pitch
+        scale = 1/sqrt(R**2 + h**2)
+        t = points[:, 3] * (self.tube_length * scale/2)
+        cos_t, sin_t = cos(t), sin(t)
+
+        # rx = R*cos_t
+        # ry = R*sin_t
+        # rz = h*t
+        # nx = -a * cos_t * points[:, 1]
+        # ny = -a * sin_t * points[:, 1]
+        # nz = 0
+        # bx = (a * h/scale) * sin_t * points[:, 2]
+        # by = (-a * h/scale) * cos_t * points[:, 2]
+        # bz = a*R/scale
+        # x = rx + nx + bx
+        # y = ry + ny + by
+        # z = rz + nz + bz
+        u, v = (R - a*points[:, 1]), (a * h/scale)*points[:, 2]
+        x = u * cos_t + v * sin_t
+        y = u * sin_t - v * cos_t
+        z = a*R/scale + h * t
+
+        points = np.hstack((x, y, z))
+        values = self.value.repeat(points.shape[0])
+        return values, self._adjust(points)
+
+def _calc_Pr_nonuniform(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)+1
+    extended_Pr = np.zeros(n_max+1, 'd')
+    # r refers to bin centers; find corresponding bin edges
+    bins = bin_edges(r)
+    t_next = time.clock() + 3
+    for k, rho_k in enumerate(rho[:-1]):
+        distance = sqrt(np.sum((points[k] - points[k+1:])**2, axis=1))
+        weights = rho_k * rho[k+1:]
+        index = np.searchsorted(bins, distance)
+        # Note: indices may be duplicated, so "Pr[index] += w" will not work!!
+        extended_Pr += np.bincount(index, weights, n_max+1)
+        t = time.clock()
+        if t > t_next:
+            t_next = t + 3
+            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)
+
+    # Make Pr a little be bigger than necessary so that only distances
+    # min < d < max end up in Pr
+    n_max = 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)
+        # Note: indices may be duplicated, so "Pr[index] += w" will not work!!
+        extended_Pr += np.bincount(index, weights, n_max+1)
+        t = time.clock()
+        if t > t_next:
+            t_next = t + 3
+            print("processing %d of %d"%(k, len(rho)-1))
+    #print("time py:", time.clock() - t0)
+    Pr = extended_Pr[:-1]
+    # 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 / Pr.max()
+
+    # Can get an additional 2x by going to C.  Cuda/OpenCL will allow even
+    # more speedup, though still bounded by the n^2 cose.
+    """
+void pdfcalc(int n, const double *pts, const double *rho,
+	     int nPr, double *Pr, double rstep)
+{
+  int i,j;
+
+  for (i=0; i<n-2; i++) {
+    for (j=i+1; j<=n-1; j++) {
+      const double dxx=pts[3*i]-pts[3*j];
+      const double dyy=pts[3*i+1]-pts[3*j+1];
+      const double dzz=pts[3*i+2]-pts[3*j+2];
+      const double d=sqrt(dxx*dxx+dyy*dyy+dzz*dzz);
+      const int k=rint(d/rstep);
+      if (k < nPr) Pr[k]+=rho[i]*rho[j];
+    }
+  }
+}
+"""
+
+def j0(x):
+    return np.sinc(x/np.pi)
+
+def calc_Iq(q, r, Pr):
+    Iq = np.array([simps(Pr * j0(qk*r), r) for qk in q])
+    #Iq = np.array([np.trapz(Pr * j0(qk*r), r) for qk in q])
+    Iq /= Iq[0]
+    return Iq
+
+# NOTE: copied from sasmodels/resolution.py
+def bin_edges(x):
+    """
+    Determine bin edges from bin centers, assuming that edges are centered
+    between the bins.
+
+    Note: this uses the arithmetic mean, which may not be appropriate for
+    log-scaled data.
+    """
+    if len(x) < 2 or (np.diff(x) < 0).any():
+        raise ValueError("Expected bins to be an increasing set")
+    edges = np.hstack([
+        x[0]  - 0.5*(x[1]  - x[0]),  # first point minus half first interval
+        0.5*(x[1:] + x[:-1]),        # mid points of all central intervals
+        x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval
+        ])
+    return edges
+
+def plot_calc(r, Pr, q, Iq, theory=None):
+    import matplotlib.pyplot as plt
+    plt.subplot(211)
+    plt.plot(r, Pr, '-', label="Pr")
+    plt.xlabel('r (A)')
+    plt.ylabel('Pr (1/A^2)')
+    plt.subplot(212)
+    plt.loglog(q, Iq, '-', label='from Pr')
+    plt.xlabel('q (1/A')
+    plt.ylabel('Iq')
+    if theory is not None:
+        plt.loglog(theory[0], theory[1], '-', label='analytic')
+        plt.legend()
+
+def plot_points(rho, points):
+    import mpl_toolkits.mplot3d
+    import matplotlib.pyplot as plt
+
+    ax = plt.axes(projection='3d')
+    try:
+        ax.axis('square')
+    except Exception:
+        pass
+    n = len(points)
+    index = np.random.choice(n, size=1000) if n > 1000 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 sas_2J1x_x(x):
+    with np.errstate(all='ignore'):
+        retvalue = 2*J1(x)/x
+    retvalue[x == 0] = 1.
+    return retvalue
+
+def sas_3j1x_x(x):
+    """return 3*j1(x)/x"""
+    with np.errstate(all='ignore'):
+        retvalue = 3*(sin(x) - x*cos(x))/x**3
+    retvalue[x == 0.] = 1.
+    return retvalue
+
+def cylinder_Iq(q, radius, length):
+    z, w = leggauss(76)
+    cos_alpha = (z+1)/2
+    sin_alpha = sqrt(1.0 - cos_alpha**2)
+    Iq = np.empty_like(q)
+    for k, qk in enumerate(q):
+        qab, qc = qk*sin_alpha, qk*cos_alpha
+        Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2)
+        Iq[k] = np.sum(w*Fq**2)
+    Iq = Iq/Iq[0]
+    return Iq
+
+def sphere_Iq(q, radius):
+    Iq = sas_3j1x_x(q*radius)**2
+    return Iq/Iq[0]
+
+def csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core):
+    z, w = leggauss(76)
+
+    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
+
+    outer_sum = np.zeros_like(q)
+    for cos_alpha, outer_w in zip((z+1)/2, w):
+        sin_alpha = sqrt(1.0-cos_alpha*cos_alpha)
+        qc = q*cos_alpha
+        siC = c*j0(c*qc/2)
+        siCt = tC*j0(tC*qc/2)
+        inner_sum = np.zeros_like(q)
+        for beta, inner_w in zip((z + 1)*pi/4, w):
+            qa, qb = q*sin_alpha*sin(beta), q*sin_alpha*cos(beta)
+            siA = a*j0(a*qa/2)
+            siB = b*j0(b*qb/2)
+            siAt = tA*j0(tA*qa/2)
+            siBt = tB*j0(tB*qb/2)
+            if overlapping:
+                Fq = (dr0*siA*siB*siC
+                      + drA*(siAt-siA)*siB*siC
+                      + drB*siAt*(siBt-siB)*siC
+                      + drC*siAt*siBt*(siCt-siC))
+            else:
+                Fq = (dr0*siA*siB*siC
+                      + drA*(siAt-siA)*siB*siC
+                      + drB*siA*(siBt-siB)*siC
+                      + drC*siA*siB*(siCt-siC))
+            inner_sum += inner_w * Fq**2
+        outer_sum += outer_w * inner_sum
+    Iq = outer_sum / 4  # = outer*um*zm*8.0/(4.0*M_PI)
+    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 check_cylinder(radius=25, length=125, rho=2.):
+    shape = EllipticalCylinder(radius, radius, length, rho)
+    fn = lambda q: cylinder_Iq(q, radius, length)
+    check_shape(shape, fn)
+
+def check_sphere(radius=125, rho=2):
+    shape = TriaxialEllipsoid(radius, radius, radius, rho)
+    fn = lambda q: sphere_Iq(q, radius)
+    check_shape(shape, fn)
+
+def check_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4):
+    core = Box(a, b, c, sld_core)
+    side_a = Box(da, b, c, slda, center=((a+da)/2, 0, 0))
+    side_b = Box(a, db, c, sldb, center=(0, (b+db)/2, 0))
+    side_c = Box(a, b, dc, sldc, center=(0, 0, (c+dc)/2))
+    side_a2 = copy(side_a).shift(-a-da, 0, 0)
+    side_b2 = copy(side_b).shift(0, -b-db, 0)
+    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)
+
+if __name__ == "__main__":
+    check_cylinder(radius=10, length=40)
+    #check_sphere()
+    #check_csbox()
+    #check_csbox(da=100, db=200, dc=300)
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/compare.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
@@ -42,5 +42,5 @@
 from .data import plot_theory, empty_data1D, empty_data2D, load_data
 from .direct_model import DirectModel, get_mesh
-from .generate import FLOAT_RE
+from .generate import FLOAT_RE, set_integration_size
 from .weights import plot_weights
 
@@ -693,5 +693,5 @@
         data = empty_data2D(q, resolution=res)
         data.accuracy = opts['accuracy']
-        set_beam_stop(data, 0.0004)
+        set_beam_stop(data, qmin)
         index = ~data.mask
     else:
@@ -706,5 +706,5 @@
     return data, index
 
-def make_engine(model_info, data, dtype, cutoff):
+def make_engine(model_info, data, dtype, cutoff, ngauss=0):
     # type: (ModelInfo, Data, str, float) -> Calculator
     """
@@ -714,4 +714,7 @@
     than OpenCL.
     """
+    if ngauss:
+        set_integration_size(model_info, ngauss)
+
     if dtype is None or not dtype.endswith('!'):
         return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff)
@@ -954,5 +957,5 @@
     'poly', 'mono', 'cutoff=',
     'magnetic', 'nonmagnetic',
-    'accuracy=',
+    'accuracy=', 'ngauss=',
     'neval=',  # for timing...
 
@@ -1089,4 +1092,5 @@
         'show_weights' : False,
         'sphere'    : 0,
+        'ngauss'    : '0',
     }
     for arg in flags:
@@ -1115,4 +1119,5 @@
         elif arg.startswith('-engine='):   opts['engine'] = arg[8:]
         elif arg.startswith('-neval='):    opts['count'] = arg[7:]
+        elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:]
         elif arg.startswith('-random='):
             opts['seed'] = int(arg[8:])
@@ -1169,4 +1174,5 @@
 
     comparison = any(PAR_SPLIT in v for v in values)
+
     if PAR_SPLIT in name:
         names = name.split(PAR_SPLIT, 2)
@@ -1181,4 +1187,10 @@
         return None
 
+    if PAR_SPLIT in opts['ngauss']:
+        opts['ngauss'] = [int(k) for k in opts['ngauss'].split(PAR_SPLIT, 2)]
+        comparison = True
+    else:
+        opts['ngauss'] = [int(opts['ngauss'])]*2
+
     if PAR_SPLIT in opts['engine']:
         opts['engine'] = opts['engine'].split(PAR_SPLIT, 2)
@@ -1199,7 +1211,9 @@
         opts['cutoff'] = [float(opts['cutoff'])]*2
 
-    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0])
+    base = make_engine(model_info[0], data, opts['engine'][0],
+                       opts['cutoff'][0], opts['ngauss'][0])
     if comparison:
-        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1])
+        comp = make_engine(model_info[1], data, opts['engine'][1],
+                           opts['cutoff'][1], opts['ngauss'][1])
     else:
         comp = None
@@ -1274,5 +1288,5 @@
         if model_info != model_info2:
             pars2 = randomize_pars(model_info2, pars2)
-            limit_dimensions(model_info, pars2, maxdim)
+            limit_dimensions(model_info2, pars2, maxdim)
             # Share values for parameters with the same name
             for k, v in pars.items():
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 db0340684d33368d63f9c39a84e1c2398bf1e953)
+++ 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
@@ -270,4 +279,22 @@
 """
 
+
+def set_integration_size(info, n):
+    # type: (ModelInfo, int) -> None
+    """
+    Update the model definition, replacing the gaussian integration with
+    a gaussian integration of a different size.
+
+    Note: this really ought to be a method in modelinfo, but that leads to
+    import loops.
+    """
+    if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)):
+        import os.path
+        from .gengauss import gengauss
+        path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n)
+        if not os.path.exists(path):
+            gengauss(n, path)
+        info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss')
+                        else lib for lib in info.source]
 
 def format_units(units):
@@ -608,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.
@@ -621,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,
     }
 
@@ -638,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)
 
@@ -698,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 = []
@@ -710,21 +742,40 @@
 
     if model_info.c_code:
-        source.append(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
@@ -752,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/gengauss.py
===================================================================
--- sasmodels/gengauss.py	(revision ff31782cc25cfbfea355acdd6d0204c959b41dcc)
+++ sasmodels/gengauss.py	(revision ff31782cc25cfbfea355acdd6d0204c959b41dcc)
@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+from __future__ import division, print_function
+
+import sys
+
+import numpy as np
+from numpy.polynomial.legendre import leggauss
+
+def gengauss(n, path):
+    z, w = leggauss(n)
+
+    # Make sure array size is a multiple of 4
+    if n%4:
+        array_size = n + (4 - n%4)
+        z, w = [np.hstack((v, [0.]*(4-n%4))) for v in (z, w)]
+    else:
+        array_size = n
+
+    with open(path, "w") as fid:
+        fid.write("""\
+// Generated by sasmodels.gengauss.gengauss(%d)
+
+#ifdef GAUSS_N
+# undef GAUSS_N
+# undef GAUSS_Z
+# undef GAUSS_W
+#endif
+#define GAUSS_N %d
+#define GAUSS_Z Gauss%dZ
+#define GAUSS_W Gauss%dWt
+
+"""%(n, n, n, n))
+
+        if array_size != n:
+            fid.write("// Note: using array size %d so that it is a multiple of 4\n\n"%array_size)
+
+        fid.write("constant double Gauss%dWt[%d]={\n"%(n, array_size))
+        fid.write(",\n".join("\t% .15e"%v for v in w))
+        fid.write("\n};\n")
+
+        fid.write("constant double Gauss%dZ[%d]={\n"%(n, array_size))
+        fid.write(",\n".join("\t% .15e"%v for v in z))
+        fid.write("\n};")
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 6aee3abf255884200f80256eb642d49c1e9c888e)
+++ sasmodels/kernel_iq.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -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,4 +470,11 @@
   #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc)
   #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
+#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
 
@@ -477,4 +485,34 @@
 #if defined(CALL_IQ) || defined(CALL_IQ_A)
   #define APPLY_PROJECTION() const double weight=weight0
+#elif defined(CALL_IQ_XY)
+  // 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 // !spherosymmetric projection
   // Grab the "view" angles (theta, phi, psi) from the initial parameter table.
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 db0340684d33368d63f9c39a84e1c2398bf1e953)
+++ sasmodels/modelinfo.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -42,5 +42,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
@@ -261,7 +261,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
@@ -391,7 +391,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[].
@@ -448,7 +445,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']
@@ -495,9 +491,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")
@@ -720,4 +722,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
@@ -725,37 +730,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):
@@ -766,5 +766,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.
     """
@@ -796,4 +796,5 @@
     info.c_code = getattr(kernel_module, 'c_code', None)
     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', [])
@@ -803,4 +804,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
@@ -817,4 +820,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)
     try:
@@ -835,5 +842,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.
@@ -917,9 +924,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
@@ -970,6 +977,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]]
@@ -987,12 +996,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 becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/barbell.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -23,10 +23,10 @@
     const double qab_r = radius_bell*qab; // Q*R*sin(theta)
     double total = 0.0;
-    for (int i = 0; i < 76; i++){
-        const double t = Gauss76Z[i]*zm + zb;
+    for (int i = 0; i < GAUSS_N; i++){
+        const double t = GAUSS_Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
         const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
-        total += Gauss76Wt[i] * Fq;
+        total += GAUSS_W[i] * Fq;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -73,10 +73,10 @@
     const double zb = M_PI_4;
     double total = 0.0;
-    for (int i = 0; i < 76; i++){
-        const double alpha= Gauss76Z[i]*zm + zb;
+    for (int i = 0; i < GAUSS_N; i++){
+        const double alpha= GAUSS_Z[i]*zm + zb;
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(alpha, sin_alpha, cos_alpha);
         const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length);
-        total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
+        total += GAUSS_W[i] * Aq * Aq * sin_alpha;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -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 ea60e084e3e40eeb3abfc5f728359135acaed020)
+++ sasmodels/models/bcc_paracrystal.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -81,13 +81,13 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<150; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
         double inner_sum = 0.0;
-        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        const double theta = GAUSS_Z[i]*theta_m + theta_b;
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
         const double qc = q*cos_theta;
         const double qab = q*sin_theta;
-        for(int j=0;j<150;j++) {
-            const double phi = Gauss150Z[j]*phi_m + phi_b;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double phi = GAUSS_Z[j]*phi_m + phi_b;
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -95,8 +95,8 @@
             const double qb = qab*sin_phi;
             const double form = bcc_Zq(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * form;
+            inner_sum += GAUSS_W[j] * form;
         }
         inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
-        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
     outer_sum *= theta_m;
@@ -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 becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/capped_cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -30,10 +30,10 @@
     const double qab_r = radius_cap*qab; // Q*R*sin(theta)
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
-        const double t = Gauss76Z[i]*zm + zb;
+    for (int i=0; i<GAUSS_N; i++) {
+        const double t = GAUSS_Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
         const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
-        total += Gauss76Wt[i] * Fq;
+        total += GAUSS_W[i] * Fq;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -95,6 +95,6 @@
     const double zb = M_PI_4;
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
-        const double theta = Gauss76Z[i]*zm + zb;
+    for (int i=0; i<GAUSS_N ;i++) {
+        const double theta = GAUSS_Z[i]*zm + zb;
         double sin_theta, cos_theta; // slots to hold sincos function output
         SINCOS(theta, sin_theta, cos_theta);
@@ -103,5 +103,5 @@
         const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length);
         // scale by sin_theta for spherical coord integration
-        total += Gauss76Wt[i] * Aq * Aq * sin_theta;
+        total += GAUSS_W[i] * Aq * Aq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -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 becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/core_shell_bicelle.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -52,11 +52,11 @@
 
     double total = 0.0;
-    for(int i=0;i<N_POINTS_76;i++) {
-        double theta = (Gauss76Z[i] + 1.0)*uplim;
+    for(int i=0;i<GAUSS_N;i++) {
+        double theta = (GAUSS_Z[i] + 1.0)*uplim;
         double sin_theta, cos_theta; // slots to hold sincos function output
         SINCOS(theta, sin_theta, cos_theta);
         double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
                                    halflength, sld_core, sld_face, sld_rim, sld_solvent);
-        total += Gauss76Wt[i]*fq*fq*sin_theta;
+        total += GAUSS_W[i]*fq*fq*sin_theta;
     }
 
@@ -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 82592da3721b7e2098cf6346cf224814189fadb3)
+++ sasmodels/models/core_shell_bicelle_elliptical.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -37,10 +37,10 @@
     //initialize integral
     double outer_sum = 0.0;
-    for(int i=0;i<76;i++) {
+    for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
         //const double va = 0.0;
         //const double vb = 1.0;
-        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
-        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0;
+        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0;
         const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
         const double qab = q*sin_theta;
@@ -49,11 +49,11 @@
         const double si2 = sas_sinx_x((halfheight+thick_face)*qc);
         double inner_sum=0.0;
-        for(int j=0;j<76;j++) {
+        for(int j=0;j<GAUSS_N;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
             // inner integral limits
             //const double vaj=0.0;
             //const double vbj=M_PI;
-            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
-            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2;
+            //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+            const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2;
             const double rr = sqrt(r2A - r2B*cos(phi));
             const double be1 = sas_2J1x_x(rr*qab);
@@ -61,8 +61,8 @@
             const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1;
 
-            inner_sum += Gauss76Wt[j] * fq * fq;
+            inner_sum += GAUSS_W[j] * fq * fq;
         }
         //now calculate outer integral
-        outer_sum += Gauss76Wt[i] * inner_sum;
+        outer_sum += GAUSS_W[i] * inner_sum;
     }
 
@@ -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 82592da3721b7e2098cf6346cf224814189fadb3)
+++ sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -7,5 +7,5 @@
         double length)
 {
-    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 
+    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length +
                  square(r_minor)*x_core*2.0*thick_face  );
 }
@@ -47,9 +47,9 @@
     //initialize integral
     double outer_sum = 0.0;
-    for(int i=0;i<76;i++) {
+    for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
         // since we generate these lots of times, why not store them somewhere?
-        //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
-        const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0;
+        //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0;
         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
         double inner_sum=0;
@@ -58,8 +58,8 @@
         si1 = sas_sinx_x(sinarg1);
         si2 = sas_sinx_x(sinarg2);
-        for(int j=0;j<76;j++) {
+        for(int j=0;j<GAUSS_N;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
-            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
-            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2;
+            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2;
             const double rr = sqrt(r2A - r2B*cos(beta));
             double besarg1 = q*rr*sin_alpha;
@@ -67,10 +67,10 @@
             be1 = sas_2J1x_x(besarg1);
             be2 = sas_2J1x_x(besarg2);
-            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 +
+            inner_sum += GAUSS_W[j] *square(dr1*si1*be1 +
                                               dr2*si1*be2 +
                                               dr3*si2*be1);
         }
         //now calculate outer integral
-        outer_sum += Gauss76Wt[i] * inner_sum;
+        outer_sum += GAUSS_W[i] * inner_sum;
     }
 
@@ -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 becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/core_shell_cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -30,9 +30,9 @@
     const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
+    for (int i=0; i<GAUSS_N ;i++) {
         // translate a point in [-1,1] to a point in [0, pi/2]
-        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0;
         double sin_theta, cos_theta;
-        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4;
+        const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4;
         SINCOS(theta, sin_theta,  cos_theta);
         const double qab = q*sin_theta;
@@ -40,5 +40,5 @@
         const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
             + _cyl(shell_vd, shell_r*qab, shell_h*qc);
-        total += Gauss76Wt[i] * fq * fq * sin_theta;
+        total += GAUSS_W[i] * fq * fq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -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 becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/core_shell_ellipsoid.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -59,6 +59,6 @@
     const double b = 0.5;
     double total = 0.0;     //initialize intergral
-    for(int i=0;i<76;i++) {
-        const double cos_theta = Gauss76Z[i]*m + b;
+    for(int i=0;i<GAUSS_N;i++) {
+        const double cos_theta = GAUSS_Z[i]*m + b;
         const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
         double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta,
@@ -66,5 +66,5 @@
             equat_shell, polar_shell,
             sld_core_shell, sld_shell_solvent);
-        total += Gauss76Wt[i] * fq * fq;
+        total += GAUSS_W[i] * fq * fq;
     }
     total *= m;
@@ -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 c69d6d6b7d1b4a8e1acc9c734f6ec3b26084cdd7)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -1,11 +1,28 @@
+// Set OVERLAPPING to 1 in order to fill in the edges of the box, with
+// c endcaps and b overlapping a.  With the proper choice of parameters,
+// (setting rim slds to sld, core sld to solvent, rim thickness to thickness
+// and subtracting 2*thickness from length, this should match the hollow
+// rectangular prism.)  Set it to 0 for the documented behaviour.
+#define OVERLAPPING 0
 static double
 form_volume(double length_a, double length_b, double length_c,
     double thick_rim_a, double thick_rim_b, double thick_rim_c)
 {
-    //return length_a * length_b * length_c;
-    return length_a * length_b * length_c +
-           2.0 * thick_rim_a * length_b * length_c +
-           2.0 * thick_rim_b * length_a * length_c +
-           2.0 * thick_rim_c * length_a * length_b;
+    return
+#if OVERLAPPING
+        // Hollow rectangular prism only includes the volume of the shell
+        // so uncomment the next line when comparing.  Solid rectangular
+        // prism, or parallelepiped want filled cores, so comment when
+        // comparing.
+        //-length_a * length_b * length_c +
+        (length_a + 2.0*thick_rim_a) *
+        (length_b + 2.0*thick_rim_b) *
+        (length_c + 2.0*thick_rim_c);
+#else
+        length_a * length_b * length_c +
+        2.0 * thick_rim_a * length_b * length_c +
+        2.0 * length_a * thick_rim_b * length_c +
+        2.0 * length_a * length_b * thick_rim_c;
+#endif
 }
 
@@ -24,90 +41,66 @@
     double thick_rim_c)
 {
-    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
+    // 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 is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
+    // Code rewritten (PAK)
 
-    const double mu = 0.5 * q * length_b;
+    const double half_q = 0.5*q;
 
-    //calculate volume before rescaling (in original code, but not used)
-    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
-    //double vol = length_a * length_b * length_c +
-    //       2.0 * thick_rim_a * length_b * length_c +
-    //       2.0 * thick_rim_b * length_a * length_c +
-    //       2.0 * thick_rim_c * length_a * length_b;
+    const double tA = length_a + 2.0*thick_rim_a;
+    const double tB = length_b + 2.0*thick_rim_b;
+    const double tC = length_c + 2.0*thick_rim_c;
 
-    // Scale sides by B
-    const double a_scaled = length_a / length_b;
-    const double c_scaled = length_c / length_b;
-
-    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b;
-    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b;
-    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present
-
-    double Vin = length_a * length_b * length_c;
-    //double Vot = (length_a * length_b * length_c +
-    //            2.0 * thick_rim_a * length_b * length_c +
-    //            2.0 * length_a * thick_rim_b * length_c +
-    //            2.0 * length_a * length_b * thick_rim_c);
-    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
-    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
-    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present
-
-    // Scale factors (note that drC is not used later)
-    const double drho0 = (core_sld-solvent_sld);
-    const double drhoA = (arim_sld-solvent_sld);
-    const double drhoB = (brim_sld-solvent_sld);
-    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot;
-
-
-    // Precompute scale factors for combining cross terms from the shape
-    const double scale23 = drhoA*V1;
-    const double scale14 = drhoB*V2;
-    const double scale24 = drhoC*V3;
-    const double scale11 = drho0*Vin;
-    const double scale12 = drho0*Vin - scale23 - scale14 - scale24;
+    // Scale factors
+    const double dr0 = (core_sld-solvent_sld);
+    const double drA = (arim_sld-solvent_sld);
+    const double drB = (brim_sld-solvent_sld);
+    const double drC = (crim_sld-solvent_sld);
 
     // outer integral (with gauss points), integration limits = 0, 1
-    double outer_total = 0; //initialize integral
+    double outer_sum = 0; //initialize integral
+    for( int i=0; i<GAUSS_N; i++) {
+        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 );
+        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha);
 
-    for( int i=0; i<76; i++) {
-        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
-        double mu_proj = mu * sqrt(1.0-sigma*sigma);
+        // inner integral (with gauss points), integration limits = 0, pi/2
+        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q);
+        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q);
+        double inner_sum = 0.0;
+        for(int j=0; j<GAUSS_N; j++) {
+            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 );
+            double sin_beta, cos_beta;
+            SINCOS(M_PI_2*beta, sin_beta, cos_beta);
+            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta);
+            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta);
+            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta);
+            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta);
 
-        // inner integral (with gauss points), integration limits = 0, 1
-        double inner_total = 0.0;
-        double inner_total_crim = 0.0;
-        for(int j=0; j<76; j++) {
-            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
-            double sin_uu, cos_uu;
-            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
-            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
-            const double si2 = sas_sinx_x(mu_proj * cos_uu );
-            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);
-            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb);
+#if OVERLAPPING
+            const double f = dr0*siA*siB*siC
+                + drA*(siAt-siA)*siB*siC
+                + drB*siAt*(siBt-siB)*siC
+                + drC*siAt*siBt*(siCt-siC);
+#else
+            const double f = dr0*siA*siB*siC
+                + drA*(siAt-siA)*siB*siC
+                + drB*siA*(siBt-siB)*siC
+                + drC*siA*siB*(siCt-siC);
+#endif
 
-            // Expression in libCylinder.c (neither drC nor Vot are used)
-            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4;
-            const double form_crim = scale11*si1*si2;
-
-            //  correct FF : sum of square of phase factors
-            inner_total += Gauss76Wt[j] * form * form;
-            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim;
+            inner_sum += GAUSS_W[j] * f * f;
         }
-        inner_total *= 0.5;
-        inner_total_crim *= 0.5;
+        inner_sum *= 0.5;
         // now sum up the outer integral
-        const double si = sas_sinx_x(mu * c_scaled * sigma);
-        const double si_crim = sas_sinx_x(mu * tc * sigma);
-        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim);
+        outer_sum += GAUSS_W[i] * inner_sum;
     }
-    outer_total *= 0.5;
+    outer_sum *= 0.5;
 
     //convert from [1e-12 A-1] to [cm-1]
-    return 1.0e-4 * outer_total;
+    return 1.0e-4 * outer_sum;
 }
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double core_sld,
     double arim_sld,
@@ -128,34 +121,27 @@
     const double drC = crim_sld-solvent_sld;
 
-    double Vin = length_a * length_b * length_c;
-    double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc)
-    double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
-    double V3 = 2.0 * length_a * length_b * thick_rim_c;
-    // As for the 1D case, Vot is not used
-    //double Vot = (length_a * length_b * length_c +
-    //              2.0 * thick_rim_a * length_b * length_c +
-    //              2.0 * length_a * thick_rim_b * length_c +
-    //              2.0 * length_a * length_b * thick_rim_c);
-
     // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
     // the scaling by B.
-    double ta = length_a + 2.0*thick_rim_a;
-    double tb = length_b + 2.0*thick_rim_b;
-    double tc = length_c + 2.0*thick_rim_c;
-    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    double siA = sas_sinx_x(0.5*length_a*qa);
-    double siB = sas_sinx_x(0.5*length_b*qb);
-    double siC = sas_sinx_x(0.5*length_c*qc);
-    double siAt = sas_sinx_x(0.5*ta*qa);
-    double siBt = sas_sinx_x(0.5*tb*qb);
-    double siCt = sas_sinx_x(0.5*tc*qc);
+    const double tA = length_a + 2.0*thick_rim_a;
+    const double tB = length_b + 2.0*thick_rim_b;
+    const double tC = length_c + 2.0*thick_rim_c;
+    const double siA = length_a*sas_sinx_x(0.5*length_a*qa);
+    const double siB = length_b*sas_sinx_x(0.5*length_b*qb);
+    const double siC = length_c*sas_sinx_x(0.5*length_c*qc);
+    const double siAt = tA*sas_sinx_x(0.5*tA*qa);
+    const double siBt = tB*sas_sinx_x(0.5*tB*qb);
+    const double siCt = tC*sas_sinx_x(0.5*tC*qc);
 
-
-    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
-    // in the 1D code, but should be checked!
-    double f = ( dr0*siA*siB*siC*Vin
-               + drA*(siAt-siA)*siB*siC*V1
-               + drB*siA*(siBt-siB)*siC*V2
-               + drC*siA*siB*(siCt-siC)*V3);
+#if OVERLAPPING
+    const double f = dr0*siA*siB*siC
+        + drA*(siAt-siA)*siB*siC
+        + drB*siAt*(siBt-siB)*siC
+        + drC*siAt*siBt*(siCt-siC);
+#else
+    const double f = dr0*siA*siB*siC
+        + drA*(siAt-siA)*siB*siC
+        + drB*siA*(siBt-siB)*siC
+        + drC*siA*siB*(siCt-siC);
+#endif
 
     return 1.0e-4 * f * f;
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision 10ee8389f2d5511546ae78ffddc6d16e36d2a6a8)
@@ -5,11 +5,5 @@
 Calculates the form factor for a rectangular solid with a core-shell structure.
 The thickness and the scattering length density of the shell or
-"rim" can be different on each (pair) of faces. However at this time the 1D
-calculation does **NOT** actually calculate a c face rim despite the presence
-of the parameter. Some other aspects of the 1D calculation may be wrong.
-
-.. note::
-   This model was originally ported from NIST IGOR macros. However, it is not
-   yet fully understood by the SasView developers and is currently under review.
+"rim" can be different on each (pair) of faces.
 
 The form factor is normalized by the particle volume $V$ such that
@@ -21,5 +15,4 @@
 where $\langle \ldots \rangle$ is an average over all possible orientations
 of the rectangular solid.
-
 
 The function calculated is the form factor of the rectangular solid below.
@@ -41,57 +34,77 @@
     V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
 
-**meaning that there are "gaps" at the corners of the solid.**  Again note that
-$t_C = 0$ currently.
+**meaning that there are "gaps" at the corners of the solid.**
 
 The intensity calculated follows the :ref:`parallelepiped` model, with the
 core-shell intensity being calculated as the square of the sum of the
-amplitudes of the core and shell, in the same manner as a core-shell model.
-
-.. math::
-
-    F_{a}(Q,\alpha,\beta)=
-    \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta)
-               }{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta}
-    - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta)
-           }{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right]
-    \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta)
-               }{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right]
-    \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta)
-               }{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right]
-
-.. note::
-
-    Why does t_B not appear in the above equation?
-    For the calculation of the form factor to be valid, the sides of the solid
-    MUST (perhaps not any more?) be chosen such that** $A < B < C$.
-    If this inequality is not satisfied, the model will not report an error,
-    but the calculation will not be correct and thus the result wrong.
+amplitudes of the core and the slabs on the edges.
+
+the scattering amplitude is computed for a particular orientation of the
+core-shell parallelepiped with respect to the scattering vector and then
+averaged over all possible orientations, where $\alpha$ is the angle between
+the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is
+the angle between projection of the particle in the $xy$ detector plane
+and the $y$ axis.
+
+.. math::
+
+    F(Q)
+    &= (\rho_\text{core}-\rho_\text{solvent})
+       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\
+    &+ (\rho_\text{A}-\rho_\text{solvent})
+        \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\
+    &+ (\rho_\text{B}-\rho_\text{solvent})
+        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\
+    &+ (\rho_\text{C}-\rho_\text{solvent})
+        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right]
+
+with
+
+.. math::
+
+    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L}
+
+and
+
+.. math::
+
+    Q_A &= \sin\alpha \sin\beta \\
+    Q_B &= \sin\alpha \cos\beta \\
+    Q_C &= \cos\alpha
+
+
+where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$
+are the scattering length of the parallelepiped core, and the rectangular
+slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$
+is the scattering length of the solvent.
 
 FITTING NOTES
+~~~~~~~~~~~~~
+
 If the scale is set equal to the particle volume fraction, $\phi$, the returned
-value is the scattered intensity per unit volume, $I(q) = \phi P(q)$.
-However, **no interparticle interference effects are included in this
-calculation.**
+value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However,
+**no interparticle interference effects are included in this calculation.**
 
 There are many parameters in this model. Hold as many fixed as possible with
 known values, or you will certainly end up at a solution that is unphysical.
 
-Constraints must be applied during fitting to ensure that the inequality
-$A < B < C$ is not violated. The calculation will not report an error,
-but the results will not be correct.
-
 The returned value is in units of |cm^-1|, on absolute scale.
 
 NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated
 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
-and length $(C+2t_C)$ values, after appropriately
-sorting the three dimensions to give an oblate or prolate particle, to give an
-effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied.
+and length $(C+2t_C)$ values, after appropriately sorting the three dimensions
+to give an oblate or prolate particle, to give an effective radius,
+for $S(Q)$ when $P(Q) * S(Q)$ is applied.
 
 For 2d data the orientation of the particle is required, described using
-angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
-of the calculation and angular dispersions see :ref:`orientation` .
+angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further
+details of the calculation and angular dispersions see :ref:`orientation`.
 The angle $\Psi$ is the rotational angle around the *long_c* axis. For example,
 $\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector.
+
+For 2d, constraints must be applied during fitting to ensure that the
+inequality $A < B < C$ is not violated, and hence the correct definition
+of angles is preserved. The calculation will not report an error,
+but the results may be not correct.
 
 .. figure:: img/parallelepiped_angle_definition.png
@@ -114,5 +127,5 @@
     Equations (1), (13-14). (in German)
 .. [#] D Singh (2009). *Small angle scattering studies of self assembly in
-   lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available
+   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available
    from Proquest <http://search.proquest.com/docview/304915826?accountid
    =26379>`_
@@ -175,12 +188,10 @@
         Return equivalent radius (ER)
     """
-
-    # surface average radius (rough approximation)
-    surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi)
-
-    height = length_c + 2.0*thick_rim_c
-
-    ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad))
-    return 0.5 * (ddd) ** (1. / 3.)
+    from .parallelepiped import ER as ER_p
+
+    a = length_a + 2*thick_rim_a
+    b = length_b + 2*thick_rim_b
+    c = length_c + 2*thick_rim_c
+    return ER_p(a, b, c)
 
 # VR defaults to 1.0
@@ -216,5 +227,6 @@
             psi_pd=10, psi_pd_n=1)
 
-# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
+# rkh 7/4/17 add random unit test for 2d, note make all params different,
+# 2d values not tested against other codes or models
 if 0:  # pak: model rewrite; need to update tests
     qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -21,11 +21,11 @@
 
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
-        const double theta = Gauss76Z[i]*zm + zb;
+    for (int i=0; i<GAUSS_N ;i++) {
+        const double theta = GAUSS_Z[i]*zm + zb;
         double sin_theta, cos_theta; // slots to hold sincos function output
         // theta (theta,phi) the projection of the cylinder on the detector plane
         SINCOS(theta , sin_theta, cos_theta);
         const double form = fq(q*sin_theta, q*cos_theta, radius, length);
-        total += Gauss76Wt[i] * form * form * sin_theta;
+        total += GAUSS_W[i] * form * form * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -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 becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/ellipsoid.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -22,13 +22,13 @@
 
     // translate a point in [-1,1] to a point in [0, 1]
-    // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2;
+    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2;
     const double zm = 0.5;
     const double zb = 0.5;
     double total = 0.0;
-    for (int i=0;i<76;i++) {
-        const double u = Gauss76Z[i]*zm + zb;
+    for (int i=0;i<GAUSS_N;i++) {
+        const double u = GAUSS_Z[i]*zm + zb;
         const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one);
         const double f = sas_3j1x_x(q*r);
-        total += Gauss76Wt[i] * f * f;
+        total += GAUSS_W[i] * f * f;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -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 82592da3721b7e2098cf6346cf224814189fadb3)
+++ sasmodels/models/elliptical_cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -22,16 +22,15 @@
     //initialize integral
     double outer_sum = 0.0;
-    for(int i=0;i<76;i++) {
+    for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
-        const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_val = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
         const double sin_val = sqrt(1.0 - cos_val*cos_val);
         //const double arg = radius_minor*sin_val;
         double inner_sum=0;
-        for(int j=0;j<76;j++) {
-            //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017
-            const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
             const double r = sin_val*sqrt(rA - rB*cos(theta));
             const double be = sas_2J1x_x(q*r);
-            inner_sum += Gauss76Wt[j] * be * be;
+            inner_sum += GAUSS_W[j] * be * be;
         }
         //now calculate the value of the inner integral
@@ -40,5 +39,5 @@
         //now calculate outer integral
         const double si = sas_sinx_x(q*0.5*length*cos_val);
-        outer_sum += Gauss76Wt[i] * inner_sum * si * si;
+        outer_sum += GAUSS_W[i] * inner_sum * si * si;
     }
     outer_sum *= 0.5*(vb-va);
@@ -55,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/elliptical_cylinder.py
===================================================================
--- sasmodels/models/elliptical_cylinder.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/elliptical_cylinder.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
@@ -121,6 +121,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c",
-          "elliptical_cylinder.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"]
 
 demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0,
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision f728001e783d4f2cad1f52bfa5374c457648c268)
+++ sasmodels/models/fcc_paracrystal.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -53,13 +53,13 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<150; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
         double inner_sum = 0.0;
-        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        const double theta = GAUSS_Z[i]*theta_m + theta_b;
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
         const double qc = q*cos_theta;
         const double qab = q*sin_theta;
-        for(int j=0;j<150;j++) {
-            const double phi = Gauss150Z[j]*phi_m + phi_b;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double phi = GAUSS_Z[j]*phi_m + phi_b;
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -67,8 +67,8 @@
             const double qb = qab*sin_phi;
             const double form = fcc_Zq(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * form;
+            inner_sum += GAUSS_W[j] * form;
         }
         inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
-        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
     outer_sum *= theta_m;
@@ -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/flexible_cylinder_elliptical.c
===================================================================
--- sasmodels/models/flexible_cylinder_elliptical.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/flexible_cylinder_elliptical.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -17,11 +17,11 @@
     double sum=0.0;
 
-    for(int i=0;i<N_POINTS_76;i++) {
-        const double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4;
+    for(int i=0;i<GAUSS_N;i++) {
+        const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4;
         double sn, cn;
         SINCOS(zi, sn, cn);
         const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn);
         const double yyy = sas_2J1x_x(arg);
-        sum += Gauss76Wt[i] * yyy * yyy;
+        sum += GAUSS_W[i] * yyy * yyy;
     }
     sum *= 0.5;
Index: sasmodels/models/hollow_cylinder.c
===================================================================
--- sasmodels/models/hollow_cylinder.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/hollow_cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -38,10 +38,10 @@
 
     double summ = 0.0;            //initialize intergral
-    for (int i=0;i<76;i++) {
-        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper );
+    for (int i=0;i<GAUSS_N;i++) {
+        const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper );
         const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
         const double form = _fq(q*sin_theta, q*cos_theta,
                                 radius, thickness, length);
-        summ += Gauss76Wt[i] * form * form;
+        summ += GAUSS_W[i] * form * form;
     }
 
@@ -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 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -1,4 +1,4 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness);
-double Iq(double q, double sld, double solvent_sld, double length_a, 
+double Iq(double q, double sld, double solvent_sld, double length_a,
           double b2a_ratio, double c2a_ratio, double thickness);
 
@@ -37,9 +37,9 @@
     const double v2a = 0.0;
     const double v2b = M_PI_2;  //phi integration limits
-    
+
     double outer_sum = 0.0;
-    for(int i=0; i<76; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
 
-        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
+        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
@@ -49,7 +49,7 @@
 
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
+        for(int j=0; j<GAUSS_N; j++) {
 
-            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
+            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -66,9 +66,9 @@
             const double AP2 = vol_core * termA2 * termB2 * termC2;
 
-            inner_sum += Gauss76Wt[j] * square(AP1-AP2);
+            inner_sum += GAUSS_W[j] * square(AP1-AP2);
         }
         inner_sum *= 0.5 * (v2b-v2a);
 
-        outer_sum += Gauss76Wt[i] * inner_sum * sin(theta);
+        outer_sum += GAUSS_W[i] * inner_sum * sin(theta);
     }
     outer_sum *= 0.5*(v1b-v1a);
@@ -84,2 +84,39 @@
     return 1.0e-4 * delrho * delrho * form;
 }
+
+double Iqabc(double qa, double qb, double qc,
+    double sld,
+    double solvent_sld,
+    double length_a,
+    double b2a_ratio,
+    double c2a_ratio,
+    double thickness)
+{
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double a_half = 0.5 * length_a;
+    const double b_half = 0.5 * length_b;
+    const double c_half = 0.5 * length_c;
+    const double vol_total = length_a * length_b * length_c;
+    const double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness);
+
+    // Amplitude AP from eqn. (13)
+
+    const double termA1 = sas_sinx_x(qa * a_half);
+    const double termA2 = sas_sinx_x(qa * (a_half-thickness));
+
+    const double termB1 = sas_sinx_x(qb * b_half);
+    const double termB2 = sas_sinx_x(qb * (b_half-thickness));
+
+    const double termC1 = sas_sinx_x(qc * c_half);
+    const double termC2 = sas_sinx_x(qc * (c_half-thickness));
+
+    const double AP1 = vol_total * termA1 * termB1 * termC1;
+    const double AP2 = vol_core * termA2 * termB2 * termC2;
+
+    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
+    const double delrho = sld - solvent_sld;
+
+    // Convert from [1e-12 A-1] to [cm-1]
+    return 1.0e-4 * square(delrho * (AP1-AP2));
+}
Index: sasmodels/models/hollow_rectangular_prism.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/hollow_rectangular_prism.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
@@ -5,5 +5,5 @@
 This model provides the form factor, $P(q)$, for a hollow rectangular
 parallelepiped with a wall of thickness $\Delta$.
-It computes only the 1D scattering, not the 2D.
+
 
 Definition
@@ -66,5 +66,25 @@
 (which is unitless).
 
-**The 2D scattering intensity is not computed by this model.**
+For 2d data the orientation of the particle is required, described using
+angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
+of the calculation and angular dispersions see :ref:`orientation` .
+The angle $\Psi$ is the rotational angle around the long *C* axis. For example,
+$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector.
+
+For 2d, constraints must be applied during fitting to ensure that the inequality
+$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error,
+but the results may be not correct.
+
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented hollow rectangular prism.
+    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then
+    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the prism.
+    The neutron or X-ray beam is along the $z$ axis.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for oriented hollow rectangular prisms against the
+    detector plane.
 
 
@@ -113,4 +133,10 @@
               ["thickness", "Ang", 1, [0, inf], "volume",
                "Thickness of parallelepiped"],
+              ["theta", "degrees", 0, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about c axis"],
              ]
 
Index: sasmodels/models/hollow_rectangular_prism_thin_walls.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision ab2aea81f7fb392e5f24f1c76dd989914d8f7022)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -1,4 +1,4 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio);
-double Iq(double q, double sld, double solvent_sld, double length_a, 
+double Iq(double q, double sld, double solvent_sld, double length_a,
           double b2a_ratio, double c2a_ratio);
 
@@ -29,8 +29,8 @@
     const double v2a = 0.0;
     const double v2b = M_PI_2;  //phi integration limits
-    
+
     double outer_sum = 0.0;
-    for(int i=0; i<76; i++) {
-        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
+    for(int i=0; i<GAUSS_N; i++) {
+        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
 
         double sin_theta, cos_theta;
@@ -44,6 +44,6 @@
 
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
-            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
+        for(int j=0; j<GAUSS_N; j++) {
+            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
 
             double sin_phi, cos_phi;
@@ -62,9 +62,9 @@
                 * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi );
 
-            inner_sum += Gauss76Wt[j] * square(AL+AT);
+            inner_sum += GAUSS_W[j] * square(AL+AT);
         }
 
         inner_sum *= 0.5 * (v2b-v2a);
-        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
 
Index: sasmodels/models/lib/gauss150.c
===================================================================
--- sasmodels/models/lib/gauss150.c	(revision 994d77f9b713691e41106d1f9af9661dfb2dc419)
+++ sasmodels/models/lib/gauss150.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -7,7 +7,17 @@
  *
  */
+ #ifdef GAUSS_N
+ # undef GAUSS_N
+ # undef GAUSS_Z
+ # undef GAUSS_W
+ #endif
+ #define GAUSS_N 150
+ #define GAUSS_Z Gauss150Z
+ #define GAUSS_W Gauss150Wt
+
+// Note: using array size 152 so that it is a multiple of 4
 
 // Gaussians
-constant double Gauss150Z[150]={
+constant double Gauss150Z[152]={
   	-0.9998723404457334,
   	-0.9993274305065947,
@@ -159,8 +169,10 @@
   	0.9983473449340834,
   	0.9993274305065947,
-  	0.9998723404457334
+  	0.9998723404457334,
+  	0.,
+  	0.
 };
 
-constant double Gauss150Wt[150]={
+constant double Gauss150Wt[152]={
   	0.0003276086705538,
   	0.0007624720924706,
@@ -312,4 +324,6 @@
   	0.0011976474864367,
   	0.0007624720924706,
-  	0.0003276086705538
+  	0.0003276086705538,
+  	0.,
+  	0.
 };
Index: sasmodels/models/lib/gauss20.c
===================================================================
--- sasmodels/models/lib/gauss20.c	(revision 994d77f9b713691e41106d1f9af9661dfb2dc419)
+++ sasmodels/models/lib/gauss20.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -7,4 +7,12 @@
  *
  */
+ #ifdef GAUSS_N
+ # undef GAUSS_N
+ # undef GAUSS_Z
+ # undef GAUSS_W
+ #endif
+ #define GAUSS_N 20
+ #define GAUSS_Z Gauss20Z
+ #define GAUSS_W Gauss20Wt
 
 // Gaussians
Index: sasmodels/models/lib/gauss76.c
===================================================================
--- sasmodels/models/lib/gauss76.c	(revision 66d119fbbe121774089e17e3086b4809c9cf2a80)
+++ sasmodels/models/lib/gauss76.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -7,8 +7,15 @@
  *
  */
-#define N_POINTS_76 76
+ #ifdef GAUSS_N
+ # undef GAUSS_N
+ # undef GAUSS_Z
+ # undef GAUSS_W
+ #endif
+ #define GAUSS_N 76
+ #define GAUSS_Z Gauss76Z
+ #define GAUSS_W Gauss76Wt
 
 // Gaussians
-constant double Gauss76Wt[N_POINTS_76]={
+constant double Gauss76Wt[76]={
 	.00126779163408536,		//0
 	.00294910295364247,
@@ -89,5 +96,5 @@
 };
 
-constant double Gauss76Z[N_POINTS_76]={
+constant double Gauss76Z[76]={
 	-.999505948362153,		//0
 	-.997397786355355,
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 9b7b23fffe93c3cd39b094e9e17d0bf62c39b2bc)
+++ sasmodels/models/parallelepiped.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -23,6 +23,6 @@
     double outer_total = 0; //initialize integral
 
-    for( int i=0; i<76; i++) {
-        const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
+    for( int i=0; i<GAUSS_N; i++) {
+        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
         const double mu_proj = mu * sqrt(1.0-sigma*sigma);
 
@@ -30,16 +30,16 @@
         // corresponding to angles from 0 to pi/2.
         double inner_total = 0.0;
-        for(int j=0; j<76; j++) {
-            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
+        for(int j=0; j<GAUSS_N; j++) {
+            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
             double sin_uu, cos_uu;
             SINCOS(M_PI_2*uu, sin_uu, cos_uu);
             const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
             const double si2 = sas_sinx_x(mu_proj * cos_uu);
-            inner_total += Gauss76Wt[j] * square(si1 * si2);
+            inner_total += GAUSS_W[j] * square(si1 * si2);
         }
         inner_total *= 0.5;
 
         const double si = sas_sinx_x(mu * c_scaled * sigma);
-        outer_total += Gauss76Wt[i] * inner_total * si * si;
+        outer_total += GAUSS_W[i] * inner_total * si * si;
     }
     outer_total *= 0.5;
@@ -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/pringle.c
===================================================================
--- sasmodels/models/pringle.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
+++ sasmodels/models/pringle.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -29,11 +29,11 @@
     double sumC = 0.0;		// initialize integral
     double r;
-    for (int i=0; i < 76; i++) {
-        r = Gauss76Z[i]*zm + zb;
+    for (int i=0; i < GAUSS_N; i++) {
+        r = GAUSS_Z[i]*zm + zb;
 
         const double qrs = r*q_sin_psi;
         const double qrrc = r*r*q_cos_psi;
 
-        double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);
+        double y = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);
         double S, C;
         SINCOS(alpha*qrrc, S, C);
@@ -86,6 +86,6 @@
 
     double sum = 0.0;
-    for (int i = 0; i < 76; i++) {
-        double psi = Gauss76Z[i]*zm + zb;
+    for (int i = 0; i < GAUSS_N; i++) {
+        double psi = GAUSS_Z[i]*zm + zb;
         double sin_psi, cos_psi;
         SINCOS(psi, sin_psi, cos_psi);
@@ -93,5 +93,5 @@
         double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0));
         double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term;
-        sum += Gauss76Wt[i] * pringle_kernel;
+        sum += GAUSS_W[i] * pringle_kernel;
     }
 
Index: sasmodels/models/rectangular_prism.c
===================================================================
--- sasmodels/models/rectangular_prism.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
+++ sasmodels/models/rectangular_prism.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -1,4 +1,4 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio);
-double Iq(double q, double sld, double solvent_sld, double length_a, 
+double Iq(double q, double sld, double solvent_sld, double length_a,
           double b2a_ratio, double c2a_ratio);
 
@@ -26,8 +26,8 @@
     const double v2a = 0.0;
     const double v2b = M_PI_2;  //phi integration limits
-    
+
     double outer_sum = 0.0;
-    for(int i=0; i<76; i++) {
-        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
+    for(int i=0; i<GAUSS_N; i++) {
+        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
@@ -36,6 +36,6 @@
 
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
-            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
+        for(int j=0; j<GAUSS_N; j++) {
+            double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -45,16 +45,16 @@
             const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi);
             const double AP = termA * termB * termC;
-            inner_sum += Gauss76Wt[j] * AP * AP;
+            inner_sum += GAUSS_W[j] * AP * AP;
         }
         inner_sum = 0.5 * (v2b-v2a) * inner_sum;
-        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
 
     double answer = 0.5*(v1b-v1a)*outer_sum;
 
-    // Normalize by Pi (Eqn. 16). 
-    // The term (ABC)^2 does not appear because it was introduced before on 
+    // Normalize by Pi (Eqn. 16).
+    // The term (ABC)^2 does not appear because it was introduced before on
     // the definitions of termA, termB, termC.
-    // The factor 2 appears because the theta integral has been defined between 
+    // The factor 2 appears because the theta integral has been defined between
     // 0 and pi/2, instead of 0 to pi.
     answer /= M_PI_2; //Form factor P(q)
@@ -64,7 +64,37 @@
     answer *= square((sld-solvent_sld)*volume);
 
-    // Convert from [1e-12 A-1] to [cm-1] 
+    // Convert from [1e-12 A-1] to [cm-1]
     answer *= 1.0e-4;
 
     return answer;
 }
+
+
+double Iqabc(double qa, double qb, double qc,
+    double sld,
+    double solvent_sld,
+    double length_a,
+    double b2a_ratio,
+    double c2a_ratio)
+{
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double a_half = 0.5 * length_a;
+    const double b_half = 0.5 * length_b;
+    const double c_half = 0.5 * length_c;
+    const double volume = length_a * length_b * length_c;
+
+    // Amplitude AP from eqn. (13)
+
+    const double termA = sas_sinx_x(qa * a_half);
+    const double termB = sas_sinx_x(qb * b_half);
+    const double termC = sas_sinx_x(qc * c_half);
+
+    const double AP = termA * termB * termC;
+
+    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
+    const double delrho = sld - solvent_sld;
+
+    // Convert from [1e-12 A-1] to [cm-1]
+    return 1.0e-4 * square(volume * delrho * AP);
+}
Index: sasmodels/models/rectangular_prism.py
===================================================================
--- sasmodels/models/rectangular_prism.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/rectangular_prism.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
@@ -12,6 +12,4 @@
 the prism (e.g. setting $b/a = 1$ and $c/a = 1$ and applying polydispersity
 to *a* will generate a distribution of cubes of different sizes).
-Note also that, contrary to :ref:`parallelepiped`, it does not compute
-the 2D scattering.
 
 
@@ -26,6 +24,5 @@
 that reference), with $\theta$ corresponding to $\alpha$ in that paper,
 and not to the usual convention used for example in the
-:ref:`parallelepiped` model. As the present model does not compute
-the 2D scattering, this has no further consequences.
+:ref:`parallelepiped` model.
 
 In this model the scattering from a massive parallelepiped with an
@@ -65,5 +62,26 @@
 units) *scale* represents the volume fraction (which is unitless).
 
-**The 2D scattering intensity is not computed by this model.**
+For 2d data the orientation of the particle is required, described using
+angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
+of the calculation and angular dispersions see :ref:`orientation` .
+The angle $\Psi$ is the rotational angle around the long *C* axis. For example,
+$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector.
+
+For 2d, constraints must be applied during fitting to ensure that the inequality
+$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error,
+but the results may be not correct.
+
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented core-shell parallelepipeds.
+    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then
+    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder.
+    The neutron or X-ray beam is along the $z$ axis.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for oriented rectangular prisms against the
+    detector plane.
+
 
 
@@ -108,4 +126,10 @@
               ["c2a_ratio", "", 1, [0, inf], "volume",
                "Ratio sides c/a"],
+              ["theta", "degrees", 0, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about c axis"],
              ]
 
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision f728001e783d4f2cad1f52bfa5374c457648c268)
+++ sasmodels/models/sc_paracrystal.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -54,13 +54,13 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<150; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
         double inner_sum = 0.0;
-        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        const double theta = GAUSS_Z[i]*theta_m + theta_b;
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
         const double qc = q*cos_theta;
         const double qab = q*sin_theta;
-        for(int j=0;j<150;j++) {
-            const double phi = Gauss150Z[j]*phi_m + phi_b;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double phi = GAUSS_Z[j]*phi_m + phi_b;
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -68,8 +68,8 @@
             const double qb = qab*sin_phi;
             const double form = sc_Zq(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * form;
+            inner_sum += GAUSS_W[j] * form;
         }
         inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
-        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
     outer_sum *= theta_m;
@@ -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 becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/stacked_disks.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -81,6 +81,6 @@
     double halfheight = 0.5*thick_core;
 
-    for(int i=0; i<N_POINTS_76; i++) {
-        double zi = (Gauss76Z[i] + 1.0)*M_PI_4;
+    for(int i=0; i<GAUSS_N; i++) {
+        double zi = (GAUSS_Z[i] + 1.0)*M_PI_4;
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(zi, sin_alpha, cos_alpha);
@@ -95,5 +95,5 @@
                            solvent_sld,
                            d);
-        summ += Gauss76Wt[i] * yyy * sin_alpha;
+        summ += GAUSS_W[i] * yyy * sin_alpha;
     }
 
@@ -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 becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
@@ -21,7 +21,7 @@
     const double zb = M_PI_4;
     double outer = 0.0;
-    for (int i=0;i<76;i++) {
-        //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2;
-        const double phi = Gauss76Z[i]*zm + zb;
+    for (int i=0;i<GAUSS_N;i++) {
+        //const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2;
+        const double phi = GAUSS_Z[i]*zm + zb;
         const double pa_sinsq_phi = pa*square(sin(phi));
 
@@ -29,12 +29,12 @@
         const double um = 0.5;
         const double ub = 0.5;
-        for (int j=0;j<76;j++) {
+        for (int j=0;j<GAUSS_N;j++) {
             // translate a point in [-1,1] to a point in [0, 1]
-            const double usq = square(Gauss76Z[j]*um + ub);
+            const double usq = square(GAUSS_Z[j]*um + ub);
             const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq);
             const double fq = sas_3j1x_x(q*r);
-            inner += Gauss76Wt[j] * fq * fq;
+            inner += GAUSS_W[j] * fq * fq;
         }
-        outer += Gauss76Wt[i] * inner;  // correcting for dx later
+        outer += GAUSS_W[i] * inner;  // correcting for dx later
     }
     // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi
@@ -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,
