Changes in / [d5014e4:6ceca44] in sasmodels

Files:

• .gitignore (modified) (1 diff)
• .travis.yml (modified) (2 diffs)
• appveyor.yml (modified) (1 diff)
• doc/developer/overview.rst (modified) (3 diffs)
• explore/precision.py (modified) (3 diffs)
• explore/realspace.py (modified) (11 diffs)
• sasmodels/compare.py (modified) (1 diff)
• sasmodels/kernel_iq.c (modified) (9 diffs)
• sasmodels/models/core_shell_parallelepiped.c (modified) (2 diffs)
• sasmodels/models/core_shell_parallelepiped.py (modified) (1 diff)
• sasmodels/models/lib/gauss150.c (modified) (4 diffs)
• sasmodels/models/lib/gauss20.c (modified) (1 diff)
• sasmodels/models/lib/gauss76.c (modified) (1 diff)

.gitignore

@@ -24 +24 @@
 /example/Fit_*/
 /example/batch_fit.csv
+# Note: changes to gauss20, gauss76 and gauss150 are still tracked since they
+# are part of the repo and required by some models.
 /sasmodels/models/lib/gauss*.c

.travis.yml

@@ -6 +6 @@
     env:
     - PY=2.7
+    - DEPLOY=True
   - os: linux
     env:
@@ -51 +52 @@
   on:
     branch: master
+    condition: $DEPLOY = True
 notifications:
   slack:

appveyor.yml

@@ -67 +67 @@
     - cmd: conda install --yes --quiet obvious-ci
     - cmd: conda install --yes --quiet numpy toolchain scipy cython cffi
-    - cmd: conda install --yes --channel conda-forge pyopencl
+    #- cmd: conda install --yes --channel conda-forge pyopencl
     - cmd: pip install bumps unittest-xml-reporting tinycc
 

doc/developer/overview.rst

@@ -185 +185 @@
 jitter applied before particle orientation.
 
+When computing the orientation dispersity integral, the weights for
+the individual points depends on the map projection used to translate jitter
+angles into latitude/longitude.  The choice of projection is set by
+*sasmodels.generate.PROJECTION*, with the default *PROJECTION=1* for
+equirectangular and *PROJECTION=2* for sinusoidal.  The more complicated
+Guyou and Postel projections are not implemented. See explore.jitter.draw_mesh
+for details.
+
 For numerical integration within form factors etc. sasmodels is mostly using
 Gaussian quadrature with 20, 76 or 150 points depending on the model. It also
@@ -199 +207 @@
 Useful testing routines include:
 
-:mod:`asymint` a direct implementation of the surface integral for certain
-models to get a more trusted value for the 1D integral using a
-reimplementation of the 2D kernel in python and mpmath (which computes math
-functions to arbitrary precision). It uses $\theta$ ranging from 0 to $\pi$
-and $\phi$ ranging from 0 to $2\pi$. It perhaps would benefit from including
-the U-substitution for $\theta$.
-
-:mod:`check1d` uses sasmodels 1D integration and compares that with a
+The *sascomp* utility is used to view and compare models with different
+parameters and calculation engines. The usual case is to simply plot a
+model that you are developing::
+
+    python sascomp path/to/model.py
+
+Once the obvious problems are addressed, check the numerical precision
+across a variety of randomly generated inputs::
+
+    python sascomp -engine=single,double path/to/model.py -sets=10
+
+You can compare different parameter values for the same or different models.
+For example when looking along the long axis of a cylinder ($\theta=0$),
+dispersity in $\theta$ should be equivalent to dispersity in $\phi$
+when $\phi=90$::
+
+    python sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle \\
+    phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10
+
+It turns out that they are not because the equirectangular map projection
+weights the points by $\cos(\theta)$ so $\Delta\theta$ is not identical
+to $\Delta\phi$.  Setting *PROJECTION=2* in :mod:`sasmodels.generate` helps
+somewhat.  Postel would help even more in this case, though leading
+to distortions elsewhere.  See :mod:`sasmodels.compare` for many more details.
+
+*sascomp -ngauss=n* allows you to set the number of quadrature points used
+for the 1D integration for any model.  For example, a carbon nanotube with
+length 10 $\mu$\ m and radius 1 nm is not computed correctly at high $q$::
+
+    python sascomp cylinder length=100000 radius=10 -ngauss=76,500 -double -highq
+
+Note: ticket 702 gives some forms for long rods and thin disks that may
+be used in place of the integration, depending on $q$, radius and length; if
+the cylinder model is updated with these corrections then above call will show
+no difference.
+
+*explore/check1d.py* uses sasmodels 1D integration and compares that with a
 rectangle distribution in $\theta$ and $\phi$, with $\theta$ limits set to
 $\pm 90/\sqrt(3)$ and $\phi$ limits set to $\pm 180/\sqrt(3)$ [The rectangle
@@ -214 +251 @@
 similar reasoning.] This should rotate the sample through the entire
 $\theta$-$\phi$ surface according to the pattern that you see in jitter.py when
-you modify it to use 'rectangle' rather than 'gaussian' for its distribution
-without changing the viewing angle. In order to match the 1-D pattern for
-an arbitrary viewing angle on triaxial shapes, we need to integrate
+you use 'rectangle' rather than 'gaussian' for its distribution without
+changing the viewing angle. In order to match the 1-D pattern for an arbitrary
+viewing angle on triaxial shapes, we need to integrate
 over $\theta$, $\phi$ and $\Psi$.
 
-When computing the dispersity integral, weights are scaled by
-$|\cos(\delta \theta)|$ to account for the points in $\phi$ getting closer
-together as $\delta \theta$ increases.
-[This will probably change so that instead of adjusting the weights, we will
-adjust $\delta\theta$-$\delta\phi$ mesh so that the point density in
-$\delta\phi$ is lower at larger $\delta\theta$. The flag USE_SCALED_PHI in
-*kernel_iq.c* selects an alternative algorithm.]
-
-The integrated dispersion is computed at a set of $(qx, qy)$ points $(q
-\cos(\alpha), q \sin(\alpha))$ at some angle $\alpha$ (currently angle=0) for
-each $q$ used in the 1-D integration. The individual $q$ points should be
-equivalent to asymint rect-n when the viewing angle is set to
-$(\theta,\phi,\Psi) = (90, 0, 0)$. Such tests can help to validate that 2d
-models are consistent with 1d models.
-
-:mod:`sascomp -sphere=n` uses the same rectangular distribution as check1d to
-compute the pattern of the $q_x$-$q_y$ grid.
-
-The :mod:`sascomp` utility can be used for 2d as well as 1d calculations to
-compare results for two sets of parameters or processor types, for example
-these two oriented cylinders here should be equivalent.
-
-:mod:`\./sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10`
-
+*sascomp -sphere=n* uses the same rectangular distribution as check1d to
+compute the pattern of the $q_x$-$q_y$ grid.  This ought to produce a
+spherically symmetric pattern.
+
+*explore/precision.py* investigates the accuracy of individual functions
+on the different execution platforms.  This includes the basic special
+functions as well as more complex expressions used in scattering.  In many
+cases the OpenCL function in sasmodels will use a polynomial approximation
+over part of the range to improve accuracy, as shown in::
+
+    python explore/precision.py 3j1/x
+
+*explore/realspace.py* allows you to set up a Monte Carlo simulation of your
+model by sampling random points within and computing the 1D and 2D scattering
+patterns.  This was used to check the core shell parallelepiped example.  This
+is similar to the general sas calculator in sasview, though it uses different
+code.
+
+*explore/jitter.py* is for exploring different options for handling
+orientation and orientation dispersity.  It uses *explore/guyou.py* to
+generate the Guyou projection.
+
+*explore/asymint.py* is a direct implementation of the 1D integration for
+a number of different models.  It calculates the integral for a particular
+$q$ using several different integration schemes, including mpmath with 100
+bits of precision (33 digits), so you can use it to check the target values
+for the 1D model tests.  This is not a general purpose tool; you will need to
+edit the file to change the model parameters. It does not currently
+apply the $u=cos(\theta)$ substitution that many models are using
+internally so the 76-point gaussian quadrature may not match the value
+produced by the eqivalent model from sasmodels.
+
+*explore/symint.py* is for rotationally symmetric models (just cylinder for
+now), so it can compute an entire curve rather than a single $q$ point.  It
+includes code to compare the long cylinder approximation to cylinder.
+
+*explore/rpa.py* is for checking different implementations of the RPA model
+against calculations for specific blends.  This is a work in (suspended)
+progress.
+
+There are a few modules left over in *explore* that can probably be removed.
+*angular_pd.py* was an early version of *jitter.py*.  *sc.py* and *sc.c*
+was used to explore different calculations for paracrystal models; it has
+been absorbed into *asymint.py*. *transform_angles.py* translates orientation
+parameters from the SasView 3.x definition to sasmodels.
+
+*explore/angles.py* generates code for the view and jitter transformations.
+Keep this around since it may be needed if we add new projections.
 
 Testing

explore/precision.py

@@ -345 +345 @@
 )
 add_function(
-    name="(1/2+(1-cos(x))/x^2-sin(x)/x)/x",
+    name="(1/2-sin(x)/x+(1-cos(x))/x^2)/x",
     mp_function=lambda x: (0.5 - mp.sin(x)/x + (1-mp.cos(x))/(x*x))/x,
     np_function=lambda x: (0.5 - np.sin(x)/x + (1-np.cos(x))/(x*x))/x,
@@ -609 +609 @@
     names = ", ".join(sorted(ALL_FUNCTIONS))
     print("""\
-usage: precision.py [-f/a/r] [-x<range>] name...
+usage: precision.py [-f/a/r] [-x<range>] "name" ...
 where
     -f indicates that the function value should be plotted,
@@ -620 +620 @@
       zoom indicates linear stepping in [1000, 1010]
       neg indicates linear stepping in [-100.1, 100.1]
-and name is "all [first]" or one of:
+and name is "all" or one of:
     """+names)
     sys.exit(1)

explore/realspace.py

@@ -44 +44 @@
     """
     return Rx(phi)*Ry(theta)*Rz(psi)
+
+def apply_view(points, view):
+    """
+    Apply the view transform (theta, phi, psi) to a set of points.
+
+    Points are stored in a 3 x n numpy array.
+
+    View angles are in degrees.
+    """
+    theta, phi, psi = view
+    return np.asarray((Rz(phi)*Ry(theta)*Rz(psi))*np.matrix(points.T)).T
+
+
+def invert_view(qx, qy, view):
+    """
+    Return (qa, qb, qc) for the (theta, phi, psi) view angle at detector
+    pixel (qx, qy).
+
+    View angles are in degrees.
+    """
+    theta, phi, psi = view
+    q = np.vstack((qx.flatten(), qy.flatten(), 0*qx.flatten()))
+    return np.asarray((Rz(-psi)*Ry(-theta)*Rz(-phi))*np.matrix(q))
+
 
 class Shape:
@@ -219 +243 @@
         values = self.value.repeat(points.shape[0])
         return values, self._adjust(points)
+
+NUMBA = False
+if NUMBA:
+    from numba import njit
+    @njit("f8[:](f8[:],f8[:],f8[:],f8[:],f8[:],f8[:],f8[:])")
+    def _Iqxy(values, x, y, z, qa, qb, qc):
+        Iq = np.zeros_like(qa)
+        for j in range(len(Iq)):
+            total = 0. + 0j
+            for k in range(len(Iq)):
+                total += values[k]*np.exp(1j*(qa[j]*x[k] + qb[j]*y[k] + qc[j]*z[k]))
+            Iq[j] = abs(total)**2
+        return Iq
+
+def calc_Iqxy(qx, qy, rho, points, volume=1.0, view=(0, 0, 0)):
+    qx, qy = np.broadcast_arrays(qx, qy)
+    qa, qb, qc = invert_view(qx, qy, view)
+    rho, volume = np.broadcast_arrays(rho, volume)
+    values = rho*volume
+    x, y, z = points.T
+
+    # I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r)
+    if NUMBA:
+        Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten())
+    else:
+        Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2
+              for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)]
+    return np.asarray(Iq).reshape(qx.shape) / np.sum(volume)
 
 def _calc_Pr_nonuniform(r, rho, points):
@@ -239 +291 @@
             print("processing %d of %d"%(k, len(rho)-1))
     Pr = extended_Pr[1:-1]
-    return Pr / Pr.max()
-
-def calc_Pr(r, rho, points):
-    # P(r) with uniform steps in r is 3x faster; check if we are uniform
-    # before continuing
-    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01:
-        return _calc_Pr_nonuniform(r, rho, points)
-
+    return Pr
+
+def _calc_Pr_uniform(r, rho, points):
     # Make Pr a little be bigger than necessary so that only distances
     # min < d < max end up in Pr
-    n_max = len(r)
+    dr, n_max = r[0], len(r)
     extended_Pr = np.zeros(n_max+1, 'd')
     t0 = time.clock()
     t_next = t0 + 3
-    row_zero = np.zeros(len(rho), 'i')
     for k, rho_k in enumerate(rho[:-1]):
         distances = sqrt(np.sum((points[k] - points[k+1:])**2, axis=1))
         weights = rho_k * rho[k+1:]
-        index = np.minimum(np.asarray(distances/r[0], 'i'), n_max)
+        index = np.minimum(np.asarray(distances/dr, 'i'), n_max)
         # Note: indices may be duplicated, so "Pr[index] += w" will not work!!
         extended_Pr += np.bincount(index, weights, n_max+1)
@@ -269 +315 @@
     # we are only accumulating the upper triangular distances.
     #Pr = Pr * 2 / len(rho)**2
-    return Pr / Pr.max()
+    return Pr
 
     # Can get an additional 2x by going to C.  Cuda/OpenCL will allow even
@@ -291 +337 @@
 }
 """
+
+if NUMBA:
+    @njit("f8[:](f8[:], f8[:], f8[:,:])")
+    def _calc_Pr_uniform_jit(r, rho, points):
+        dr = r[0]
+        n_max = len(r)
+        Pr = np.zeros_like(r)
+        for j in range(len(rho) - 1):
+            x, y, z = points[j, 0], points[j, 1], points[j, 2]
+            for k in range(j+1, len(rho)):
+                distance = sqrt((x - points[k, 0])**2
+                                + (y - points[k, 1])**2
+                                + (z - points[k, 2])**2)
+                index = int(distance/dr)
+                if index < n_max:
+                    Pr[index] += rho[j] * rho[k]
+        # Make Pr independent of sampling density.  The factor of 2 comes because
+        # we are only accumulating the upper triangular distances.
+        #Pr = Pr * 2 / len(rho)**2
+        return Pr
+
+
+def calc_Pr(r, rho, points):
+    # P(r) with uniform steps in r is 3x faster; check if we are uniform
+    # before continuing
+    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01:
+        Pr = _calc_Pr_nonuniform(r, rho, points)
+    else:
+        if NUMBA:
+            Pr = _calc_Pr_uniform_jit(r, rho, points)
+        else:
+            Pr = _calc_Pr_uniform(r, rho, points)
+    return Pr / Pr.max()
+
 
 def j0(x):
@@ -333 +413 @@
         plt.legend()
 
+def plot_calc_2d(qx, qy, Iqxy, theory=None):
+    import matplotlib.pyplot as plt
+    qx, qy = bin_edges(qx), bin_edges(qy)
+    #qx, qy = np.meshgrid(qx, qy)
+    if theory is not None:
+        plt.subplot(121)
+    plt.pcolormesh(qx, qy, np.log10(Iqxy))
+    plt.xlabel('qx (1/A)')
+    plt.ylabel('qy (1/A)')
+    if theory is not None:
+        plt.subplot(122)
+        plt.pcolormesh(qx, qy, np.log10(theory))
+        plt.xlabel('qx (1/A)')
+
 def plot_points(rho, points):
     import mpl_toolkits.mplot3d
@@ -343 +437 @@
         pass
     n = len(points)
-    index = np.random.choice(n, size=1000) if n > 1000 else slice(None, None)
+    #print("len points", n)
+    index = np.random.choice(n, size=500) if n > 500 else slice(None, None)
     ax.scatter(points[index, 0], points[index, 1], points[index, 2], c=rho[index])
     #low, high = points.min(axis=0), points.max(axis=0)
     #ax.axis([low[0], high[0], low[1], high[1], low[2], high[2]])
     ax.autoscale(True)
+
+def check_shape(shape, fn=None):
+    rho_solvent = 0
+    q = np.logspace(-3, 0, 200)
+    r = shape.r_bins(q, r_step=0.01)
+    sampling_density = 6*5000 / shape.volume()
+    rho, points = shape.sample(sampling_density)
+    t0 = time.time()
+    Pr = calc_Pr(r, rho-rho_solvent, points)
+    print("calc Pr time", time.time() - t0)
+    Iq = calc_Iq(q, r, Pr)
+    theory = (q, fn(q)) if fn is not None else None
+
+    import pylab
+    #plot_points(rho, points); pylab.figure()
+    plot_calc(r, Pr, q, Iq, theory=theory)
+    pylab.show()
+
+def check_shape_2d(shape, fn=None, view=(0, 0, 0)):
+    rho_solvent = 0
+    nq, qmax = 100, 1.0
+    qx = np.linspace(0.0, qmax, nq)
+    qy = np.linspace(0.0, qmax, nq)
+    Qx, Qy = np.meshgrid(qx, qy)
+    sampling_density = 50000 / shape.volume()
+    #t0 = time.time()
+    rho, points = shape.sample(sampling_density)
+    #print("sample time", time.time() - t0)
+    t0 = time.time()
+    Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view)
+    print("calc time", time.time() - t0)
+    theory = fn(Qx, Qy) if fn is not None else None
+    Iqxy += 0.001 * Iqxy.max()
+    if theory is not None:
+        theory += 0.001 * theory.max()
+
+    import pylab
+    #plot_points(rho, points); pylab.figure()
+    plot_calc_2d(qx, qy, Iqxy, theory=theory)
+    pylab.show()
+
+def sas_sinx_x(x):
+    with np.errstate(all='ignore'):
+        retvalue = sin(x)/x
+    retvalue[x == 0.] = 1.
+    return retvalue
 
 def sas_2J1x_x(x):
@@ -373 +514 @@
     Iq = Iq/Iq[0]
     return Iq
+
+def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)):
+    qa, qb, qc = invert_view(qx, qy, view)
+    qab = np.sqrt(qa**2 + qb**2)
+    Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2)
+    Iq = Fq**2
+    return Iq.reshape(qx.shape)
 
 def sphere_Iq(q, radius):
@@ -415 +563 @@
     return Iq/Iq[0]
 
-def check_shape(shape, fn=None):
-    rho_solvent = 0
-    q = np.logspace(-3, 0, 200)
-    r = shape.r_bins(q, r_step=0.01)
-    sampling_density = 15000 / shape.volume()
-    rho, points = shape.sample(sampling_density)
-    Pr = calc_Pr(r, rho-rho_solvent, points)
-    Iq = calc_Iq(q, r, Pr)
-    theory = (q, fn(q)) if fn is not None else None
-
-    import pylab
-    #plot_points(rho, points); pylab.figure()
-    plot_calc(r, Pr, q, Iq, theory=theory)
-    pylab.show()
+def csbox_Iqxy(qx, qy, a, b, c, da, db, dc, slda, sldb, sldc, sld_core, view=(0,0,0)):
+    qa, qb, qc = invert_view(qx, qy, view)
+
+    sld_solvent = 0
+    overlapping = False
+    dr0 = sld_core - sld_solvent
+    drA, drB, drC = slda-sld_solvent, sldb-sld_solvent, sldc-sld_solvent
+    tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc
+    siA = a*sas_sinx_x(a*qa/2)
+    siB = b*sas_sinx_x(b*qb/2)
+    siC = c*sas_sinx_x(c*qc/2)
+    siAt = tA*sas_sinx_x(tA*qa/2)
+    siBt = tB*sas_sinx_x(tB*qb/2)
+    siCt = tC*sas_sinx_x(tC*qc/2)
+    Fq = (dr0*siA*siB*siC
+          + drA*(siAt-siA)*siB*siC
+          + drB*siA*(siBt-siB)*siC
+          + drC*siA*siB*(siCt-siC))
+    Iq = Fq**2
+    return Iq.reshape(qx.shape)
 
 def check_cylinder(radius=25, length=125, rho=2.):
@@ -434 +588 @@
     fn = lambda q: cylinder_Iq(q, radius, length)
     check_shape(shape, fn)
+
+def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)):
+    shape = EllipticalCylinder(radius, radius, length, rho)
+    fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view)
+    check_shape_2d(shape, fn, view=view)
+
+def check_cylinder_2d_lattice(radius=25, length=125, rho=2.,
+                              view=(0, 0, 0)):
+    nx, dx = 1, 2*radius
+    ny, dy = 30, 2*radius
+    nz, dz = 30, length
+    dx, dy, dz = 2*dx, 2*dy, 2*dz
+    def center(*args):
+        sigma = 0.333
+        space = 2
+        return [(space*n+np.random.randn()*sigma)*x for n, x in args]
+    shapes = [EllipticalCylinder(radius, radius, length, rho,
+                                 #center=(ix*dx, iy*dy, iz*dz)
+                                 orientation=np.random.randn(3)*0,
+                                 center=center((ix, dx), (iy, dy), (iz, dz))
+                                )
+              for ix in range(nx)
+              for iy in range(ny)
+              for iz in range(nz)]
+    shape = Composite(shapes)
+    fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view)
+    check_shape_2d(shape, fn, view=view)
 
 def check_sphere(radius=125, rho=2):
@@ -449 +630 @@
     side_c2 = copy(side_c).shift(0, 0, -c-dc)
     shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2))
-    fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core)
-    check_shape(shape, fn)
+    def fn(q):
+        return csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core)
+    #check_shape(shape, fn)
+
+    view = (20, 30, 40)
+    def fn_xy(qx, qy):
+        return csbox_Iqxy(qx, qy, a, b, c, da, db, dc,
+                          slda, sldb, sldc, sld_core, view=view)
+    check_shape_2d(shape, fn_xy, view=view)
 
 if __name__ == "__main__":
     check_cylinder(radius=10, length=40)
+    #check_cylinder_2d(radius=10, length=40, view=(90,30,0))
+    #check_cylinder_2d_lattice(radius=10, length=50, view=(90,30,0))
     #check_sphere()
     #check_csbox()

sasmodels/compare.py

@@ -92 +92 @@
     -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
     -neval=1 sets the number of evals for more accurate timing
+    -ngauss=0 overrides the number of points in the 1-D gaussian quadrature
 
     === precision options ===

sasmodels/kernel_iq.c

@@ -174 +174 @@
 static double
 qac_apply(
-    QACRotation rotation,
+    QACRotation *rotation,
     double qx, double qy,
     double *qa_out, double *qc_out)
 {
-    const double dqc = rotation.R31*qx + rotation.R32*qy;
+    const double dqc = rotation->R31*qx + rotation->R32*qy;
     // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2
     const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy);
@@ -247 +247 @@
 static double
 qabc_apply(
-    QABCRotation rotation,
+    QABCRotation *rotation,
     double qx, double qy,
     double *qa_out, double *qb_out, double *qc_out)
 {
-    *qa_out = rotation.R11*qx + rotation.R12*qy;
-    *qb_out = rotation.R21*qx + rotation.R22*qy;
-    *qc_out = rotation.R31*qx + rotation.R32*qy;
+    *qa_out = rotation->R11*qx + rotation->R12*qy;
+    *qb_out = rotation->R21*qx + rotation->R22*qy;
+    *qc_out = rotation->R31*qx + rotation->R32*qy;
 }
 
@@ -454 +454 @@
   // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
   #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi);
-  #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc)
+  #define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc)
   #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table)
 
@@ -468 +468 @@
   local_values.table.psi = 0.;
   #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi)
-  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc)
+  #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)
@@ -479 +479 @@
 #endif
 
-// Doing jitter projection code outside the previous if block so that we don't
-// need to repeat the identical logic in the IQ_AC and IQ_ABC branches.  This
-// will become more important if we implement more projections, or more
-// complicated projections.
-#if defined(CALL_IQ) || defined(CALL_IQ_A)
+// Define APPLY_PROJECTION depending on model symmetries. We do this outside
+// the previous if block so that we don't need to repeat the identical
+// logic in the IQ_AC and IQ_ABC branches.  This will become more important
+// if we implement more projections, or more complicated projections.
+#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation
   #define APPLY_PROJECTION() const double weight=weight0
-#elif defined(CALL_IQ_XY)
+#elif defined(CALL_IQ_XY) // pass orientation to the model
   // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc
   // Need to plug the values for the orientation angles back into parameter
@@ -515 +515 @@
     #define APPLY_PROJECTION() const double weight=weight0
   #endif
-#else // !spherosymmetric projection
+#else // apply jitter and view before calling the model
   // Grab the "view" angles (theta, phi, psi) from the initial parameter table.
   const double theta = values[details->theta_par+2];
@@ -526 +526 @@
   // we go through the mesh.
   double dtheta, dphi, weight;
-  #if PROJECTION == 1
+  #if PROJECTION == 1 // equirectangular
     #define APPLY_PROJECTION() do { \
       dtheta = local_values.table.theta; \
@@ -532 +532 @@
       weight = fabs(cos(dtheta*M_PI_180)) * weight0; \
     } while (0)
-  #elif PROJECTION == 2
+  #elif PROJECTION == 2 // sinusoidal
     #define APPLY_PROJECTION() do { \
       dtheta = local_values.table.theta; \
@@ -542 +542 @@
     } while (0)
   #endif
-#endif // !spherosymmetric projection
+#endif // done defining APPLY_PROJECTION
 
 // ** define looping macros **

sasmodels/models/core_shell_parallelepiped.c

@@ -43 +43 @@
     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
-    // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
-    // Code rewritten (PAK)
+    // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker)
+    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK)
 
     const double half_q = 0.5*q;
@@ -121 +121 @@
     const double drC = crim_sld-solvent_sld;
 
-    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
-    // the scaling by B.
     const double tA = length_a + 2.0*thick_rim_a;
     const double tB = length_b + 2.0*thick_rim_b;

sasmodels/models/core_shell_parallelepiped.py

@@ -136 +136 @@
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016
-* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017
-* **Currently Under review by:** Paul Butler
+* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
+* Cross-checked against hollow rectangular prism and rectangular prism for
+  equal thickness overlapping sides, and by Monte Carlo sampling of points
+  within the shape for non-uniform, non-overlapping sides.
 """
 

sasmodels/models/lib/gauss150.c

@@ -1 @@
-/*
- *  GaussWeights.c
- *  SANSAnalysis
- *
- *  Created by Andrew Jackson on 4/23/07.
- *  Copyright 2007 __MyCompanyName__. All rights reserved.
- *
- */
+// Created by Andrew Jackson on 4/23/07
+
  #ifdef GAUSS_N
  # undef GAUSS_N
@@ -16 +10 @@
  #define GAUSS_W Gauss150Wt
 
-// Note: using array size 152 so that it is a multiple of 4
 
-// Gaussians
+// Note: using array size 152 rather than 150 so that it is a multiple of 4.
+// Some OpenCL devices prefer that vectors start and end on nice boundaries.
 constant double Gauss150Z[152]={
         -0.9998723404457334,
@@ -170 +164 @@
         0.9993274305065947,
         0.9998723404457334,
-        0.,
-        0.
+        0., // zero padding is ignored
+        0.  // zero padding is ignored
 };
 
@@ -325 +319 @@
         0.0007624720924706,
         0.0003276086705538,
-        0.,
-        0.
+        0., // zero padding is ignored
+        0.  // zero padding is ignored
 };

sasmodels/models/lib/gauss20.c

@@ -1 @@
-/*
- *  GaussWeights.c
- *  SANSAnalysis
- *
- *  Created by Andrew Jackson on 4/23/07.
- *  Copyright 2007 __MyCompanyName__. All rights reserved.
- *
- */
+// Created by Andrew Jackson on 4/23/07
+
  #ifdef GAUSS_N
  # undef GAUSS_N

sasmodels/models/lib/gauss76.c

@@ -1 @@
-/*
- *  GaussWeights.c
- *  SANSAnalysis
- *
- *  Created by Andrew Jackson on 4/23/07.
- *  Copyright 2007 __MyCompanyName__. All rights reserved.
- *
- */
+// Created by Andrew Jackson on 4/23/07
+
  #ifdef GAUSS_N
  # undef GAUSS_N