Changes in / [6ceca44:d5014e4] 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:
@@ -52 +51 @@
   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
@@ -207 +199 @@
 Useful testing routines include:
 
-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
+: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
 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
@@ -251 +214 @@
 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 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 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
 over $\theta$, $\phi$ and $\Psi$.
 
-*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.
+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`
+
 
 Testing

explore/precision.py

@@ -345 +345 @@
 )
 add_function(
-    name="(1/2-sin(x)/x+(1-cos(x))/x^2)/x",
+    name="(1/2+(1-cos(x))/x^2-sin(x)/x)/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" or one of:
+and name is "all [first]" 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:
@@ -243 +219 @@
         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):
@@ -291 +239 @@
             print("processing %d of %d"%(k, len(rho)-1))
     Pr = extended_Pr[1:-1]
-    return Pr
-
-def _calc_Pr_uniform(r, rho, points):
+    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
-    dr, n_max = r[0], len(r)
+    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/dr, 'i'), n_max)
+        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)
@@ -315 +269 @@
     # we are only accumulating the upper triangular distances.
     #Pr = Pr * 2 / len(rho)**2
-    return Pr
+    return Pr / Pr.max()
 
     # Can get an additional 2x by going to C.  Cuda/OpenCL will allow even
@@ -337 +291 @@
 }
 """
-
-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):
@@ -413 +333 @@
         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
@@ -437 +343 @@
         pass
     n = len(points)
-    #print("len points", n)
-    index = np.random.choice(n, size=500) if n > 500 else slice(None, None)
+    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 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):
@@ -514 +373 @@
     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):
@@ -563 +415 @@
     return Iq/Iq[0]
 
-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_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.):
@@ -588 +434 @@
     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):
@@ -630 +449 @@
     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))
-    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)
+    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_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
 
-// 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
+// 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() const double weight=weight0
-#elif defined(CALL_IQ_XY) // pass orientation to the model
+#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
@@ -515 +515 @@
     #define APPLY_PROJECTION() const double weight=weight0
   #endif
-#else // apply jitter and view before calling the model
+#else // !spherosymmetric projection
   // 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 // equirectangular
+  #if PROJECTION == 1
     #define APPLY_PROJECTION() do { \
       dtheta = local_values.table.theta; \
@@ -532 +532 @@
       weight = fabs(cos(dtheta*M_PI_180)) * weight0; \
     } while (0)
-  #elif PROJECTION == 2 // sinusoidal
+  #elif PROJECTION == 2
     #define APPLY_PROJECTION() do { \
       dtheta = local_values.table.theta; \
@@ -542 +542 @@
     } while (0)
   #endif
-#endif // done defining APPLY_PROJECTION
+#endif // !spherosymmetric 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 Singh's thesis now (Dirk Honecker)
-    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK)
+    // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
+    // Code rewritten (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:** 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.
+* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017
+* **Currently Under review by:** Paul Butler
 """
 

sasmodels/models/lib/gauss150.c

@@ -1 @@
-// Created by Andrew Jackson on 4/23/07
-
+/*
+ *  GaussWeights.c
+ *  SANSAnalysis
+ *
+ *  Created by Andrew Jackson on 4/23/07.
+ *  Copyright 2007 __MyCompanyName__. All rights reserved.
+ *
+ */
  #ifdef GAUSS_N
  # undef GAUSS_N
@@ -10 +16 @@
  #define GAUSS_W Gauss150Wt
 
+// Note: using array size 152 so that it is a multiple of 4
 
-// 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.
+// Gaussians
 constant double Gauss150Z[152]={
         -0.9998723404457334,
@@ -164 +170 @@
         0.9993274305065947,
         0.9998723404457334,
-        0., // zero padding is ignored
-        0.  // zero padding is ignored
+        0.,
+        0.
 };
 
@@ -319 +325 @@
         0.0007624720924706,
         0.0003276086705538,
-        0., // zero padding is ignored
-        0.  // zero padding is ignored
+        0.,
+        0.
 };

sasmodels/models/lib/gauss20.c

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

sasmodels/models/lib/gauss76.c

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