Changes in / [09141ff:6ab64c9] in sasmodels

Files:

• doc/genapi.py (modified) (2 diffs)
• explore/angles.py (added)
• explore/asymint.py (added)
• explore/bccpy.py (added)
• explore/jitter.py (modified) (11 diffs)
• explore/sc.c (added)
• explore/sc.py (added)
• explore/symint.py (added)
• sasmodels/details.py (modified) (7 diffs)
• sasmodels/kernel_header.c (modified) (1 diff)
• sasmodels/kernel_iq.c (modified) (5 diffs)
• sasmodels/kernel_iq.cl (modified) (6 diffs)
• sasmodels/kernelpy.py (modified) (3 diffs)
• sasmodels/model_test.py (modified) (1 diff)
• sasmodels/models/_py_parallelepiped.py (added)
• sasmodels/models/barbell.c (modified) (5 diffs)
• sasmodels/models/bcc_paracrystal.c (modified) (2 diffs)
• sasmodels/models/capped_cylinder.c (modified) (5 diffs)
• sasmodels/models/core_shell_bicelle.c (modified) (3 diffs)
• sasmodels/models/core_shell_bicelle_elliptical.c (modified) (4 diffs)
• sasmodels/models/core_shell_cylinder.c (modified) (3 diffs)
• sasmodels/models/core_shell_ellipsoid.c (modified) (4 diffs)
• sasmodels/models/core_shell_ellipsoid.py (modified) (1 diff)
• sasmodels/models/core_shell_parallelepiped.c (modified) (7 diffs)
• sasmodels/models/cylinder.c (modified) (4 diffs)
• sasmodels/models/ellipsoid.c (modified) (3 diffs)
• sasmodels/models/elliptical_cylinder.c (modified) (4 diffs)
• sasmodels/models/fcc_paracrystal.c (modified) (2 diffs)
• sasmodels/models/hollow_cylinder.c (modified) (4 diffs)
• sasmodels/models/parallelepiped.c (modified) (2 diffs)
• sasmodels/models/sc_paracrystal.c (modified) (1 diff)
• sasmodels/models/sc_paracrystal.py (modified) (1 diff)
• sasmodels/models/stacked_disks.c (modified) (8 diffs)
• sasmodels/models/triaxial_ellipsoid.c (modified) (3 diffs)
• sasmodels/special.py (added)
• sasmodels/weights.py (modified) (1 diff)

doc/genapi.py

@@ -59 +59 @@
     #('alignment', 'GPU data alignment [unused]'),
     ('bumps_model', 'Bumps interface'),
+    ('compare', 'Compare models on different compute engines'),
     ('compare_many', 'Batch compare models on different compute engines'),
-    ('compare', 'Compare models on different compute engines'),
+    ('conversion_table', 'Model conversion table'),
     ('convert', 'Sasview to sasmodel converter'),
     ('core', 'Model access'),
@@ -82 +83 @@
     ('sasview_model', 'Sasview interface'),
     ('sesans', 'SESANS calculation routines'),
+    ('special', 'Special functions library'),
     ('weights', 'Distribution functions'),
 ]

explore/jitter.py

@@ -3 +3 @@
 dispersity and possible replacement algorithms.
 """
+import sys
+
 import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot
 import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider, CheckButtons
 from matplotlib import cm
-
 import numpy as np
 from numpy import pi, cos, sin, sqrt, exp, degrees, radians
 
-def draw_beam(ax):
+def draw_beam(ax, view=(0, 0)):
     #ax.plot([0,0],[0,0],[1,-1])
     #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8)
@@ -22 +23 @@
     x = r*np.outer(np.cos(u), np.ones_like(v))
     y = r*np.outer(np.sin(u), np.ones_like(v))
-    z = np.outer(np.ones_like(u), v)
+    z = 1.3*np.outer(np.ones_like(u), v)
+
+    theta, phi = view
+    shape = x.shape
+    points = np.matrix([x.flatten(), y.flatten(), z.flatten()])
+    points = Rz(phi)*Ry(theta)*points
+    x, y, z = [v.reshape(shape) for v in points]
 
     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5)
 
-def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi):
-    size=[0.1, 0.4, 1.0]
-    view=[theta, phi, psi]
-    shimmy=[0,0,0]
-    #draw_shape = draw_parallelepiped
-    draw_shape = draw_ellipsoid
+def draw_jitter(ax, view, jitter):
+    size = [0.1, 0.4, 1.0]
+    draw_shape = draw_parallelepiped
+    #draw_shape = draw_ellipsoid
 
     #np.random.seed(10)
@@ -64 +69 @@
         [ 1,  1,  1],
     ]
+    dtheta, dphi, dpsi = jitter
     if dtheta == 0:
         cloud = [v for v in cloud if v[0] == 0]
@@ -70 +76 @@
     if dpsi == 0:
         cloud = [v for v in cloud if v[2] == 0]
-    draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8)
+    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8)
     for point in cloud:
-        shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]]
-        draw_shape(ax, size, view, shimmy, alpha=0.8)
+        delta = [dtheta*point[0], dphi*point[1], dpsi*point[2]]
+        draw_shape(ax, size, view, delta, alpha=0.8)
     for v in 'xyz':
         a, b, c = size
@@ -80 +86 @@
         getattr(ax, v+'axis').label.set_text(v)
 
-def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1):
+def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1):
     a,b,c = size
-    theta, phi, psi = view
-    dtheta, dphi, dpsi = shimmy
-
     u = np.linspace(0, 2 * np.pi, steps)
     v = np.linspace(0, np.pi, steps)
@@ -90 +93 @@
     y = b*np.outer(np.sin(u), np.sin(v))
     z = c*np.outer(np.ones_like(u), np.cos(v))
-
-    shape = x.shape
-    points = np.matrix([x.flatten(),y.flatten(),z.flatten()])
-    points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points
-    points = Rz(phi)*Ry(theta)*Rz(psi)*points
-    x,y,z = [v.reshape(shape) for v in points]
+    x, y, z = transform_xyz(view, jitter, x, y, z)
 
     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
 
-def draw_parallelepiped(ax, size, view, shimmy, alpha=1):
+    draw_labels(ax, view, jitter, [
+         ('c+', [ 0, 0, c], [ 1, 0, 0]),
+         ('c-', [ 0, 0,-c], [ 0, 0,-1]),
+         ('a+', [ a, 0, 0], [ 0, 0, 1]),
+         ('a-', [-a, 0, 0], [ 0, 0,-1]),
+         ('b+', [ 0, b, 0], [-1, 0, 0]),
+         ('b-', [ 0,-b, 0], [-1, 0, 0]),
+    ])
+
+def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1):
     a,b,c = size
-    theta, phi, psi = view
-    dtheta, dphi, dpsi = shimmy
-
     x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
     y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
@@ -118 +122 @@
     ])
 
-    points = np.matrix([x,y,z])
-    points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points
-    points = Rz(phi)*Ry(theta)*Rz(psi)*points
-
-    x,y,z = [np.array(v).flatten() for v in points]
+    x, y, z = transform_xyz(view, jitter, x, y, z)
     ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
 
-def draw_sphere(ax, radius=10., steps=100):
-    u = np.linspace(0, 2 * np.pi, steps)
-    v = np.linspace(0, np.pi, steps)
-
-    x = radius * np.outer(np.cos(u), np.sin(v))
-    y = radius * np.outer(np.sin(u), np.sin(v))
-    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))
-    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
-
-def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'):
-    theta_center = radians(theta)
-    phi_center = radians(phi)
-    flow_center = radians(flow)
-    dtheta = radians(dtheta)
-    dphi = radians(dphi)
-
-    # 10 point 3-sigma gaussian weights
-    t = np.linspace(-3., 3., 11)
+    draw_labels(ax, view, jitter, [
+         ('c+', [ 0, 0, c], [ 1, 0, 0]),
+         ('c-', [ 0, 0,-c], [ 0, 0,-1]),
+         ('a+', [ a, 0, 0], [ 0, 0, 1]),
+         ('a-', [-a, 0, 0], [ 0, 0,-1]),
+         ('b+', [ 0, b, 0], [-1, 0, 0]),
+         ('b-', [ 0,-b, 0], [-1, 0, 0]),
+    ])
+
+def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gauss'):
+    theta, phi, psi = view
+    dtheta, dphi, dpsi = jitter
     if dist == 'gauss':
+        t = np.linspace(-3, 3, n)
         weights = exp(-0.5*t**2)
     elif dist == 'rect':
+        t = np.linspace(0, 1, n)
         weights = np.ones_like(t)
     else:
         raise ValueError("expected dist to be 'gauss' or 'rect'")
-    theta = dtheta*t
-    phi = dphi*t
-
-    x = radius * np.outer(cos(phi), cos(theta))
-    y = radius * np.outer(sin(phi), cos(theta))
-    z = radius * np.outer(np.ones_like(phi), sin(theta))
-    #w = np.outer(weights, weights*abs(cos(dtheta*t)))
-    w = np.outer(weights, weights*abs(cos(theta)))
-
-    x, y, z, w = [v.flatten() for v in (x,y,z,w)]
-    x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center)
-
-    ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.)
-
-def rotate(x, y, z, phi, theta, psi):
-    R = Rz(phi)*Ry(theta)*Rz(psi)
-    p = np.vstack([x,y,z])
-    return R*p
+
+    # mesh in theta, phi formed by rotating z
+    z = np.matrix([[0], [0], [radius]])
+    points = np.hstack([Rx(phi_i)*Ry(theta_i)*z
+                        for theta_i in dtheta*t
+                        for phi_i in dphi*t])
+    # rotate relative to beam
+    points = orient_relative_to_beam(view, points)
+
+    w = np.outer(weights, weights)
+
+    x, y, z = [np.array(v).flatten() for v in points]
+    ax.scatter(x, y, z, c=w.flatten(), marker='o', vmin=0., vmax=1.)
 
 def Rx(angle):
@@ -188 +179 @@
          [0., 0., 1.]]
     return np.matrix(R)
+
+def transform_xyz(view, jitter, x, y, z):
+    x, y, z = [np.asarray(v) for v in (x, y, z)]
+    shape = x.shape
+    points = np.matrix([x.flatten(),y.flatten(),z.flatten()])
+    points = apply_jitter(jitter, points)
+    points = orient_relative_to_beam(view, points)
+    x, y, z = [np.array(v).reshape(shape) for v in points]
+    return x, y, z
+
+def apply_jitter(jitter, points):
+    dtheta, dphi, dpsi = jitter
+    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
+    return points
+
+def orient_relative_to_beam(view, points):
+    theta, phi, psi = view
+    points = Rz(phi)*Ry(theta)*Rz(psi)*points
+    return points
+
+def draw_labels(ax, view, jitter, text):
+    labels, locations, orientations = zip(*text)
+    px, py, pz = zip(*locations)
+    dx, dy, dz = zip(*orientations)
+
+    px, py, pz = transform_xyz(view, jitter, px, py, pz)
+    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz)
+
+    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)):
+        zdir = np.asarray(zdir).flatten()
+        ax.text(p[0], p[1], p[2], label, zdir=zdir)
+
+def draw_sphere(ax, radius=10., steps=100):
+    u = np.linspace(0, 2 * np.pi, steps)
+    v = np.linspace(0, np.pi, steps)
+
+    x = radius * np.outer(np.cos(u), np.sin(v))
+    y = radius * np.outer(np.sin(u), np.sin(v))
+    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))
+    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
 
 def main():
@@ -206 +237 @@
 
     axcolor = 'lightgoldenrodyellow'
+
     axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor)
     axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor)
@@ -212 +244 @@
     sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi)
     spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi)
+
     axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor)
     axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor)
@@ -217 +250 @@
     sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta)
     sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi)
-    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi)
+    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dpsi)
 
     def update(val, axis=None):
-        theta, phi, psi = stheta.val, sphi.val, spsi.val
-        dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val
+        view = stheta.val, sphi.val, spsi.val
+        jitter = sdtheta.val, sdphi.val, sdpsi.val
         ax.cla()
-        draw_beam(ax)
-        draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi)
-        #if not axis.startswith('d'):
-        #    ax.view_init(elev=theta, azim=phi)
+        draw_beam(ax, (0, 0))
+        draw_jitter(ax, view, jitter)
+        #draw_jitter(ax, view, (0,0,0))
+        draw_mesh(ax, view, jitter)
         plt.gcf().canvas.draw()
 

sasmodels/details.py

@@ -15 +15 @@
 
 import numpy as np  # type: ignore
-from numpy import pi, cos, sin
+from numpy import pi, cos, sin, radians
 
 try:
@@ -29 +29 @@
 
 try:
-    from typing import List
+    from typing import List, Tuple, Sequence
 except ImportError:
     pass
 else:
     from .modelinfo import ModelInfo
+    from .modelinfo import ParameterTable
 
 
@@ -53 +54 @@
     coordinates, the total circumference decreases as latitude varies from
     pi r^2 at the equator to 0 at the pole, and the weight associated
-    with a range of phi values needs to be scaled by this circumference.
+    with a range of latitude values needs to be scaled by this circumference.
     This scale factor needs to be updated each time the theta value
     changes.  *theta_par* indicates which of the values in the parameter
@@ -231 +232 @@
     nvalues = kernel.info.parameters.nvalues
     scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs]
-    values, weights = zip(*pairs[2:npars+2]) if npars else ((),())
+    # skipping scale and background when building values and weights
+    values, weights = zip(*pairs[2:npars+2]) if npars else ((), ())
+    weights = correct_theta_weights(kernel.info.parameters, values, weights)
     length = np.array([len(w) for w in weights])
     offset = np.cumsum(np.hstack((0, length)))
@@ -244 +247 @@
     return call_details, data, is_magnetic
 
+def correct_theta_weights(parameters, values, weights):
+    # type: (ParameterTable, Sequence[np.ndarray], Sequence[np.ndarray]) -> Sequence[np.ndarray]
+    """
+    If there is a theta parameter, update the weights of that parameter so that
+    the cosine weighting required for polar integration is preserved.  Avoid
+    evaluation strictly at the pole, which would otherwise send the weight to
+    zero.
+
+    Note: values and weights do not include scale and background
+    """
+    # TODO: document code, explaining why scale and background are skipped
+    # given that we don't have scale and background in the list, we
+    # should be calling the variables something other than values and weights
+    # Apparently the parameters.theta_offset similarly skips scale and
+    # and background, so the indexing works out.
+    if parameters.theta_offset >= 0:
+        index = parameters.theta_offset
+        theta = values[index]
+        theta_weight = np.minimum(abs(cos(radians(theta))), 1e-6)
+        # copy the weights list so we can update it
+        weights = list(weights)
+        weights[index] = theta_weight*np.asarray(weights[index])
+        weights = tuple(weights)
+    return weights
+
 
 def convert_magnetism(parameters, values):
+    # type: (ParameterTable, Sequence[np.ndarray]) -> bool
     """
     Convert magnetism values from polar to rectangular coordinates.
@@ -255 +284 @@
     scale = mag[:,0]
     if np.any(scale):
-        theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180.
+        theta, phi = radians(mag[:, 1]), radians(mag[:, 2])
         cos_theta = cos(theta)
         mag[:, 0] = scale*cos_theta*cos(phi)  # mx
@@ -269 +298 @@
     """
     Create a mesh grid of dispersion parameters and weights.
+
+    pars is a list of pairs (values, weights), where the values are the
+    individual parameter values at which to evaluate the polydispersity
+    distribution and weights are the weights associated with each value.
+
+    Only the volume parameters should be included in this list.  Orientation
+    parameters do not affect the calculation of effective radius or volume
+    ratio.
 
     Returns [p1,p2,...],w where pj is a vector of values for parameter j

sasmodels/kernel_header.c

@@ -164 +164 @@
     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));  \
+    cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \
     sn = sqrt(1 - cn*cn); \
     } while (0)

sasmodels/kernel_iq.c

@@ -25 +25 @@
     int32_t num_weights;        // total length of the weights vector
     int32_t num_active;         // number of non-trivial pd loops
-    int32_t theta_par;          // id of spherical correction variable
+    int32_t theta_par;          // id of spherical correction variable (not used)
 } ProblemDetails;
 
@@ -173 +173 @@
 
 
-#if MAX_PD>0
-  const int theta_par = details->theta_par;
-  const int fast_theta = (theta_par == p0);
-  const int slow_theta = (theta_par >= 0 && !fast_theta);
-  double spherical_correction = 1.0;
-#else
-  // Note: if not polydisperse the weights cancel and we don't need the
-  // spherical correction.
-  const double spherical_correction = 1.0;
-#endif
-
   int step = pd_start;
 
@@ -220 +209 @@
 #endif
 #if MAX_PD>0
-  if (slow_theta) { // Theta is not in inner loop
-    spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6);
-  }
   while(i0 < n0) {
     local_values.vector[p0] = v0[i0];
     double weight0 = w0[i0] * weight1;
 //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0);
-    if (fast_theta) { // Theta is in inner loop
-      spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6);
-    }
 #else
     const double weight0 = 1.0;
@@ -244 +227 @@
       // Note: weight==0 must always be excluded
       if (weight0 > cutoff) {
-        // spherical correction is set at a minimum of 1e-6, otherwise there
-        // would be problems looking at models with theta=90.
-        const double weight = weight0 * spherical_correction;
-        pd_norm += weight * CALL_VOLUME(local_values.table);
+        pd_norm += weight0 * CALL_VOLUME(local_values.table);
 
         #ifdef USE_OPENMP
@@ -304 +284 @@
 #endif // !MAGNETIC
 //printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0);
-          result[q_index] += weight * scattering;
+          result[q_index] += weight0 * scattering;
         }
       }

sasmodels/kernel_iq.cl

@@ -25 +25 @@
     int32_t num_weights;        // total length of the weights vector
     int32_t num_active;         // number of non-trivial pd loops
-    int32_t theta_par;          // id of spherical correction variable
+    int32_t theta_par;          // id of spherical correction variable (not used)
 } ProblemDetails;
 
@@ -169 +169 @@
 
 
-#if MAX_PD>0
-  const int theta_par = details->theta_par;
-  const bool fast_theta = (theta_par == p0);
-  const bool slow_theta = (theta_par >= 0 && !fast_theta);
-  double spherical_correction = 1.0;
-#else
-  // Note: if not polydisperse the weights cancel and we don't need the
-  // spherical correction.
-  const double spherical_correction = 1.0;
-#endif
-
   int step = pd_start;
 
@@ -217 +206 @@
 #endif
 #if MAX_PD>0
-  if (slow_theta) { // Theta is not in inner loop
-    spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6);
-  }
   while(i0 < n0) {
     local_values.vector[p0] = v0[i0];
     double weight0 = w0[i0] * weight1;
 //if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0);
-    if (fast_theta) { // Theta is in inner loop
-      spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6);
-    }
 #else
     const double weight0 = 1.0;
@@ -232 +215 @@
 
 //if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n"); }
-//if (q_index == 0) printf("sphcor: %g\n", spherical_correction);
 
     #ifdef INVALID
@@ -241 +223 @@
       // Note: weight==0 must always be excluded
       if (weight0 > cutoff) {
-        // spherical correction is set at a minimum of 1e-6, otherwise there
-        // would be problems looking at models with theta=90.
-        const double weight = weight0 * spherical_correction;
-        pd_norm += weight * CALL_VOLUME(local_values.table);
+        pd_norm += weight0 * CALL_VOLUME(local_values.table);
 
 #if defined(MAGNETIC) && NUM_MAGNETIC > 0
@@ -296 +275 @@
         const double scattering = CALL_IQ(q, q_index, local_values.table);
 #endif // !MAGNETIC
-        this_result += weight * scattering;
+        this_result += weight0 * scattering;
       }
     }

sasmodels/kernelpy.py

@@ -201 +201 @@
 
     pd_norm = 0.0
-    spherical_correction = 1.0
     partial_weight = np.NaN
     weight = np.NaN
 
     p0_par = call_details.pd_par[0]
-    p0_is_theta = (p0_par == call_details.theta_par)
     p0_length = call_details.pd_length[0]
     p0_index = p0_length
@@ -223 +221 @@
             parameters[pd_par] = pd_value[pd_offset+pd_index]
             partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
-            if call_details.theta_par >= 0:
-                cor = sin(pi / 180 * parameters[call_details.theta_par])
-                spherical_correction = max(abs(cor), 1e-6)
             p0_index = loop_index%p0_length
 
         weight = partial_weight * pd_weight[p0_offset + p0_index]
         parameters[p0_par] = pd_value[p0_offset + p0_index]
-        if p0_is_theta:
-            cor = cos(pi/180 * parameters[p0_par])
-            spherical_correction = max(abs(cor), 1e-6)
         p0_index += 1
         if weight > cutoff:
@@ -244 +236 @@
 
             # update value and norm
-            weight *= spherical_correction
             total += weight * Iq
             pd_norm += weight * form_volume()

sasmodels/model_test.py

@@ -86 +86 @@
     suite = unittest.TestSuite()
 
-    if models[0] == 'all':
+    if models[0] in core.KINDS:
         skip = models[1:]
-        models = list_models()
+        models = list_models(models[0])
     else:
         skip = []

sasmodels/models/barbell.c

@@ -10 +10 @@
 //barbell kernel - same as dumbell
 static double
-_bell_kernel(double q, double h, double radius_bell,
-             double half_length, double sin_alpha, double cos_alpha)
+_bell_kernel(double qab, double qc, double h, double radius_bell,
+             double half_length)
 {
     // translate a point in [-1,1] to a point in [lower,upper]
@@ -26 +26 @@
     //    m = q R cos(alpha)
     //    b = q(L/2-h) cos(alpha)
-    const double m = q*radius_bell*cos_alpha; // cos argument slope
-    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept
-    const double qrst = q*radius_bell*sin_alpha; // Q*R*sin(theta)
+    const double m = radius_bell*qc; // cos argument slope
+    const double b = (half_length-h)*qc; // cos argument intercept
+    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;
         const double radical = 1.0 - t*t;
-        const double bj = sas_2J1x_x(qrst*sqrt(radical));
+        const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -44 +44 @@
 
 static double
-_fq(double q, double h,
-    double radius_bell, double radius, double half_length,
-    double sin_alpha, double cos_alpha)
+_fq(double qab, double qc, double h,
+    double radius_bell, double radius, double half_length)
 {
-    const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha);
-    const double bj = sas_2J1x_x(q*radius*sin_alpha);
-    const double si = sas_sinx_x(q*half_length*cos_alpha);
+    const double bell_fq = _bell_kernel(qab, qc, h, radius_bell, half_length);
+    const double bj = sas_2J1x_x(radius*qab);
+    const double si = sas_sinx_x(half_length*qc);
     const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si;
     const double Aq = bell_fq + cyl_fq;
@@ -84 +83 @@
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(alpha, sin_alpha, cos_alpha);
-        const double Aq = _fq(q, h, radius_bell, radius, half_length, 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;
     }
@@ -103 +102 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
     const double h = -sqrt(square(radius_bell) - square(radius));
-    const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha);
+    const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length);
 
     // Multiply by contrast^2 and convert to cm-1

sasmodels/models/bcc_paracrystal.c

@@ -1 @@
-double form_volume(double radius);
-double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double dnn,
-    double d_factor, double radius,double sld, double solvent_sld,
-    double theta, double phi, double psi);
+static double
+bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
+{
+#if 0  // Equations as written in Matsuoka
+    const double a1 = (+qa + qb + qc)/2.0;
+    const double a2 = (-qa - qb + qc)/2.0;
+    const double a3 = (-qa + qb - qc)/2.0;
+#else
+    const double a1 = (+qa + qb - qc)/2.0;
+    const double a2 = (+qa - qb + qc)/2.0;
+    const double a3 = (-qa + qb + qc)/2.0;
+#endif
 
-double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
-double _BCCeval(double Theta, double Phi, double temp1, double temp3);
-double _sphereform(double q, double radius, double sld, double solvent_sld);
+#if 1
+    // Numerator: (1 - exp(a)^2)^3
+    //         => (-(exp(2a) - 1))^3
+    //         => -expm1(2a)^3
+    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
+    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double exp_arg = exp(arg);
+    const double Zq = -cube(expm1(2.0*arg))
+        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
+#else
+    // Alternate form, which perhaps is more approachable
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double sinh_qd = sinh(arg);
+    const double cosh_qd = cosh(arg);
+    const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1))
+                    * sinh_qd/(cosh_qd - cos(dnn*a2))
+                    * sinh_qd/(cosh_qd - cos(dnn*a3));
+#endif
+
+    return Zq;
+}
 
 
-double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
-
-        const double Da = d_factor*dnn;
-        const double temp1 = q*q*Da*Da;
-        const double temp3 = q*dnn;
-
-        double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
-        return(retVal);
+// occupied volume fraction calculated from lattice symmetry and sphere radius
+static double
+bcc_volume_fraction(double radius, double dnn)
+{
+    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
 }
 
-double _BCCeval(double Theta, double Phi, double temp1, double temp3) {
-
-        double result;
-        double sin_theta,cos_theta,sin_phi,cos_phi;
-        SINCOS(Theta, sin_theta, cos_theta);
-        SINCOS(Phi, sin_phi, cos_phi);
-
-        const double temp6 =  sin_theta;
-        const double temp7 =  sin_theta*cos_phi + sin_theta*sin_phi + cos_theta;
-        const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta;
-        const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta;
-
-        const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9));
-        result = cube(1.0 - (temp10*temp10))*temp6
-            / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
-              * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
-              * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10));
-
-        return (result);
-}
-
-double form_volume(double radius){
+static double
+form_volume(double radius)
+{
     return sphere_volume(radius);
 }
 
 
-double Iq(double q, double dnn,
+static double Iq(double q, double dnn,
   double d_factor, double radius,
-  double sld, double solvent_sld){
+  double sld, double solvent_sld)
+{
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    const double phi_m = M_PI;
+    const double phi_b = M_PI;
+    // translate a point in [-1,1] to a point in [0, pi]
+    const double theta_m = M_PI_2;
+    const double theta_b = M_PI_2;
 
-        //Volume fraction calculated from lattice symmetry and sphere radius
-        const double s1 = dnn/sqrt(0.75);
-        const double latticescale = 2.0*sphere_volume(radius/s1);
-
-    const double va = 0.0;
-    const double vb = 2.0*M_PI;
-    const double vaj = 0.0;
-    const double vbj = M_PI;
-
-    double summ = 0.0;
-    double answer = 0.0;
-        for(int i=0; i<150; i++) {
-                //setup inner integral over the ellipsoidal cross-section
-                double summj=0.0;
-                const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi
-                for(int j=0;j<150;j++) {
-                        //20 gauss points for the inner integral
-                        double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta
-                        double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi);
-                        summj += yyy;
-                }
-                //now calculate the value of the inner integral
-                double answer = (vbj-vaj)/2.0*summj;
-
-                //now calculate outer integral
-                summ = summ+(Gauss150Wt[i] * answer);
-        }               //final scaling is done at the end of the function, after the NT_FP64 case
-
-        answer = (vb-va)/2.0*summ;
-        answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
-
-    return answer;
+    double outer_sum = 0.0;
+    for(int i=0; i<150; i++) {
+        double inner_sum = 0.0;
+        const double theta = Gauss150Z[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;
+            double sin_phi, cos_phi;
+            SINCOS(phi, sin_phi, cos_phi);
+            const double qa = qab*cos_phi;
+            const double qb = qab*sin_phi;
+            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * form;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m;
+    const double Zq = outer_sum/(4.0*M_PI);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
 }
 
 
-double Iqxy(double qx, double qy,
+static double Iqxy(double qx, double qy,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld,
@@ -92 +98 @@
     double q, zhat, yhat, xhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
-    const double a1 = +xhat - zhat + yhat;
-    const double a2 = +xhat + zhat - yhat;
-    const double a3 = -xhat + zhat + yhat;
-
-    const double qd = 0.5*q*dnn;
-    const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3);
-    const double tanh_qd = tanh(arg);
-    const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a2)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a3)/cosh_qd);
-
-    const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq;
-    //the occupied volume of the lattice
-    const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
-    return lattice_scale * Fq;
+    q = sqrt(qa*qa + qb*qb + qc*qc);
+    const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
 }

sasmodels/models/capped_cylinder.c

@@ -14 +14 @@
 //   radius_cap is the radius of the lens
 //   length is the cylinder length, or the separation between the lens halves
-//   alpha is the angle of the cylinder wrt q.
+//   theta is the angle of the cylinder wrt q.
 static double
-_cap_kernel(double q, double h, double radius_cap,
-                      double half_length, double sin_alpha, double cos_alpha)
+_cap_kernel(double qab, double qc, double h, double radius_cap,
+            double half_length)
 {
     // translate a point in [-1,1] to a point in [lower,upper]
@@ -26 +26 @@
 
     // cos term in integral is:
-    //    cos (q (R t - h + L/2) cos(alpha))
+    //    cos (q (R t - h + L/2) cos(theta))
     // so turn it into:
     //    cos (m t + b)
     // where:
-    //    m = q R cos(alpha)
-    //    b = q(L/2-h) cos(alpha)
-    const double m = q*radius_cap*cos_alpha; // cos argument slope
-    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept
-    const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta)
+    //    m = q R cos(theta)
+    //    b = q(L/2-h) cos(theta)
+    const double m = radius_cap*qc; // cos argument slope
+    const double b = (half_length-h)*qc; // cos argument intercept
+    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;
         const double radical = 1.0 - t*t;
-        const double bj = sas_2J1x_x(qrst*sqrt(radical));
+        const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -50 +50 @@
 
 static double
-_fq(double q, double h, double radius_cap, double radius, double half_length,
-    double sin_alpha, double cos_alpha)
+_fq(double qab, double qc, double h, double radius_cap, double radius, double half_length)
 {
-    const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha);
-    const double bj = sas_2J1x_x(q*radius*sin_alpha);
-    const double si = sas_sinx_x(q*half_length*cos_alpha);
+    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length);
+    const double bj = sas_2J1x_x(radius*qab);
+    const double si = sas_sinx_x(half_length*qc);
     const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si;
     const double Aq = cap_Fq + cyl_Fq;
@@ -101 +100 @@
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
-        const double alpha= Gauss76Z[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, h, radius_cap, radius, half_length, sin_alpha, cos_alpha);
-        // sin_alpha for spherical coord integration
-        total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
+        const double theta = Gauss76Z[i]*zm + zb;
+        double sin_theta, cos_theta; // slots to hold sincos function output
+        SINCOS(theta, sin_theta, cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        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;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -125 +125 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
     const double h = sqrt(radius_cap*radius_cap - radius*radius);
-    const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha);
+    const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length);
 
     // Multiply by contrast^2 and convert to cm-1

sasmodels/models/core_shell_bicelle.c

@@ -1 @@
-double form_volume(double radius, double thick_rim, double thick_face, double length);
-double Iq(double q,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld);
-
-
-double Iqxy(double qx, double qy,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi);
-
-
-double form_volume(double radius, double thick_rim, double thick_face, double length)
+static double
+form_volume(double radius, double thick_rim, double thick_face, double length)
 {
-    return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face);
+    return M_PI*square(radius+thick_rim)*(length+2.0*thick_face);
 }
 
 static double
-bicelle_kernel(double q,
-              double rad,
-              double radthick,
-              double facthick,
-              double halflength,
-              double rhoc,
-              double rhoh,
-              double rhor,
-              double rhosolv,
-              double sin_alpha,
-              double cos_alpha)
+bicelle_kernel(double qab,
+    double qc,
+    double radius,
+    double thick_radius,
+    double thick_face,
+    double halflength,
+    double sld_core,
+    double sld_face,
+    double sld_rim,
+    double sld_solvent)
 {
-    const double dr1 = rhoc-rhoh;
-    const double dr2 = rhor-rhosolv;
-    const double dr3 = rhoh-rhor;
-    const double vol1 = M_PI*square(rad)*2.0*(halflength);
-    const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick);
-    const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick);
+    const double dr1 = sld_core-sld_face;
+    const double dr2 = sld_rim-sld_solvent;
+    const double dr3 = sld_face-sld_rim;
+    const double vol1 = M_PI*square(radius)*2.0*(halflength);
+    const double vol2 = M_PI*square(radius+thick_radius)*2.0*(halflength+thick_face);
+    const double vol3 = M_PI*square(radius)*2.0*(halflength+thick_face);
 
-    const double be1 = sas_2J1x_x(q*(rad)*sin_alpha);
-    const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha);
-    const double si1 = sas_sinx_x(q*(halflength)*cos_alpha);
-    const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha);
+    const double be1 = sas_2J1x_x((radius)*qab);
+    const double be2 = sas_2J1x_x((radius+thick_radius)*qab);
+    const double si1 = sas_sinx_x((halflength)*qc);
+    const double si2 = sas_sinx_x((halflength+thick_face)*qc);
 
     const double t = vol1*dr1*si1*be1 +
@@ -58 +33 @@
                      vol3*dr3*si2*be1;
 
-    const double retval = t*t;
-
-    return retval;
-
+    return t;
 }
 
 static double
-bicelle_integration(double q,
-                   double rad,
-                   double radthick,
-                   double facthick,
-                   double length,
-                   double rhoc,
-                   double rhoh,
-                   double rhor,
-                   double rhosolv)
+Iq(double q,
+    double radius,
+    double thick_radius,
+    double thick_face,
+    double length,
+    double sld_core,
+    double sld_face,
+    double sld_rim,
+    double sld_solvent)
 {
     // set up the integration end points
@@ -79 +51 @@
     const double halflength = 0.5*length;
 
-    double summ = 0.0;
+    double total = 0.0;
     for(int i=0;i<N_POINTS_76;i++) {
-        double alpha = (Gauss76Z[i] + 1.0)*uplim;
-        double sin_alpha, cos_alpha; // slots to hold sincos function output
-        SINCOS(alpha, sin_alpha, cos_alpha);
-        double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick,
-                             halflength, rhoc, rhoh, rhor, rhosolv,
-                             sin_alpha, cos_alpha);
-        summ += yyy*sin_alpha;
+        double theta = (Gauss76Z[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;
     }
 
     // calculate value of integral to return
-    double answer = uplim*summ;
-    return answer;
+    double answer = total*uplim;
+    return 1.0e-4*answer;
 }
 
 static double
-bicelle_kernel_2d(double qx, double qy,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
+Iqxy(double qx, double qy,
+    double radius,
+    double thick_rim,
+    double thick_face,
+    double length,
+    double core_sld,
+    double face_sld,
+    double rim_sld,
+    double solvent_sld,
+    double theta,
+    double phi)
 {
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
-    double answer = bicelle_kernel(q, radius, thick_rim, thick_face,
+    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face,
                            0.5*length, core_sld, face_sld, rim_sld,
-                           solvent_sld, sin_alpha, cos_alpha);
-    return 1.0e-4*answer;
+                           solvent_sld);
+    return 1.0e-4*fq*fq;
 }
-
-double Iq(double q,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld)
-{
-    double intensity = bicelle_integration(q, radius, thick_rim, thick_face,
-                       length, core_sld, face_sld, rim_sld, solvent_sld);
-    return intensity*1.0e-4;
-}
-
-
-double Iqxy(double qx, double qy,
-          double radius,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
-{
-    double intensity = bicelle_kernel_2d(qx, qy,
-                      radius,
-                      thick_rim,
-                      thick_face,
-                      length,
-                      core_sld,
-                      face_sld,
-                      rim_sld,
-                      solvent_sld,
-                      theta,
-                      phi);
-
-    return intensity;
-}

sasmodels/models/core_shell_bicelle_elliptical.c

@@ -17 +17 @@
         double thick_face,
         double length,
-        double rhoc,
-        double rhoh,
-        double rhor,
-        double rhosolv)
+        double sld_core,
+        double sld_face,
+        double sld_rim,
+        double sld_solvent)
 {
-    double si1,si2,be1,be2;
      // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle
      // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle
-     //    const double uplim = M_PI_4;
     const double halfheight = 0.5*length;
-    //const double va = 0.0;
-    //const double vb = 1.0;
-    // inner integral limits
-    //const double vaj=0.0;
-    //const double vbj=M_PI;
-
     const double r_major = r_minor * x_core;
     const double r2A = 0.5*(square(r_major) + square(r_minor));
     const double r2B = 0.5*(square(r_major) - square(r_minor));
-    const double dr1 = (rhoc-rhoh)   *M_PI*r_minor*r_major*(2.0*halfheight);;
-    const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
-    const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
-    //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);
-    //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
-    //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
+    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);
+    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
+    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
+    const double dr1 = vol1*(sld_core-sld_face);
+    const double dr2 = vol2*(sld_rim-sld_solvent);
+    const double dr3 = vol3*(sld_face-sld_rim);
 
     //initialize integral
@@ -47 +39 @@
     for(int i=0;i<76;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 sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
-        double inner_sum=0;
-        double sinarg1 = q*halfheight*cos_alpha;
-        double sinarg2 = q*(halfheight+thick_face)*cos_alpha;
-        si1 = sas_sinx_x(sinarg1);
-        si2 = sas_sinx_x(sinarg2);
+        //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 sin_theta = sqrt(1.0 - cos_theta*cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        const double si1 = sas_sinx_x(halfheight*qc);
+        const double si2 = sas_sinx_x((halfheight+thick_face)*qc);
+        double inner_sum=0.0;
         for(int j=0;j<76;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 rr = sqrt(r2A - r2B*cos(beta));
-            double besarg1 = q*rr*sin_alpha;
-            double besarg2 = q*(rr+thick_rim)*sin_alpha;
-            be1 = sas_2J1x_x(besarg1);
-            be2 = sas_2J1x_x(besarg2);
-            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 +
-                                              dr2*si2*be2 +
-                                              dr3*si2*be1);
+            // 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 rr = sqrt(r2A - r2B*cos(phi));
+            const double be1 = sas_2J1x_x(rr*qab);
+            const double be2 = sas_2J1x_x((rr+thick_rim)*qab);
+            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1;
+
+            inner_sum += Gauss76Wt[j] * fq * fq;
         }
         //now calculate outer integral
@@ -83 +77 @@
           double thick_face,
           double length,
-          double rhoc,
-          double rhoh,
-          double rhor,
-          double rhosolv,
+          double sld_core,
+          double sld_face,
+          double sld_rim,
+          double sld_solvent,
           double theta,
           double phi,
           double psi)
 {
-       // THIS NEEDS TESTING
     double q, xhat, yhat, zhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
-    const double dr1 = rhoc-rhoh;
-    const double dr2 = rhor-rhosolv;
-    const double dr3 = rhoh-rhor;
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
+
+    const double dr1 = sld_core-sld_face;
+    const double dr2 = sld_rim-sld_solvent;
+    const double dr3 = sld_face-sld_rim;
     const double r_major = r_minor*x_core;
     const double halfheight = 0.5*length;
@@ -104 +101 @@
 
     // Compute effective radius in rotated coordinates
-    const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat));
-    const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat)
-                                   + square((r_minor+thick_rim)*yhat));
-    const double be1 = sas_2J1x_x( q*r_hat );
-    const double be2 = sas_2J1x_x( q*rshell_hat );
-    const double si1 = sas_sinx_x( q*halfheight*zhat );
-    const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat );
-    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1);
-    return 1.0e-4 * Aq;
+    const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb));
+    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa)
+                                   + square((r_minor+thick_rim)*qb));
+    const double be1 = sas_2J1x_x( qr_hat );
+    const double be2 = sas_2J1x_x( qrshell_hat );
+    const double si1 = sas_sinx_x( halfheight*qc );
+    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc );
+    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1;
+    return 1.0e-4 * fq*fq;
 }
-

sasmodels/models/core_shell_cylinder.c

@@ -1 @@
-double form_volume(double radius, double thickness, double length);
-double Iq(double q, double core_sld, double shell_sld, double solvent_sld,
-    double radius, double thickness, double length);
-double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld,
-    double radius, double thickness, double length, double theta, double phi);
-
 // vd = volume * delta_rho
-// besarg = q * R * sin(alpha)
-// siarg = q * L/2 * cos(alpha)
-double _cyl(double vd, double besarg, double siarg);
-double _cyl(double vd, double besarg, double siarg)
+// besarg = q * R * sin(theta)
+// siarg = q * L/2 * cos(theta)
+static double _cyl(double vd, double besarg, double siarg)
 {
     return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg);
 }
 
-double form_volume(double radius, double thickness, double length)
+static double
+form_volume(double radius, double thickness, double length)
 {
-    return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness);
+    return M_PI*square(radius+thickness)*(length+2.0*thickness);
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double core_sld,
     double shell_sld,
@@ -28 +23 @@
 {
     // precalculate constants
-    const double core_qr = q*radius;
-    const double core_qh = q*0.5*length;
+    const double core_r = radius;
+    const double core_h = 0.5*length;
     const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld);
-    const double shell_qr = q*(radius + thickness);
-    const double shell_qh = q*(0.5*length + thickness);
+    const double shell_r = (radius + thickness);
+    const double shell_h = (0.5*length + thickness);
     const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
     double total = 0.0;
-    // double lower=0, upper=M_PI_2;
     for (int i=0; i<76 ;i++) {
-        // translate a point in [-1,1] to a point in [lower,upper]
-        //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
-        double sn, cn;
-        const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2);
-        SINCOS(alpha, sn, cn);
-        const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn)
-            + _cyl(shell_vd, shell_qr*sn, shell_qh*cn);
-        total += Gauss76Wt[i] * fq * fq * sn;
+        // translate a point in [-1,1] to a point in [0, pi/2]
+        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        double sin_theta, cos_theta;
+        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4;
+        SINCOS(theta, sin_theta,  cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        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;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -64 +60 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
-    const double core_qr = q*radius;
-    const double core_qh = q*0.5*length;
+    const double core_r = radius;
+    const double core_h = 0.5*length;
     const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld);
-    const double shell_qr = q*(radius + thickness);
-    const double shell_qh = q*(0.5*length + thickness);
+    const double shell_r = (radius + thickness);
+    const double shell_h = (0.5*length + thickness);
     const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
 
-    const double fq = _cyl(core_vd, core_qr*sin_alpha, core_qh*cos_alpha)
-        + _cyl(shell_vd, shell_qr*sin_alpha, shell_qh*cos_alpha);
+    const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
+        + _cyl(shell_vd, shell_r*qab, shell_h*qc);
     return 1.0e-4 * fq * fq;
 }

sasmodels/models/core_shell_ellipsoid.c

@@ -1 @@
-double form_volume(double radius_equat_core,
-                   double polar_core,
-                   double equat_shell,
-                   double polar_shell);
-double Iq(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld);
 
+// Converted from Igor function gfn4, using the same pattern as ellipsoid
+// for evaluating the parts of the integral.
+//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
+//                       BY (53) & (58-59) IN CHEN AND
+//                       KOTLARCHYK REFERENCE
+//
+//       <OBLATE ELLIPSOID>
+static double
+_cs_ellipsoid_kernel(double qab, double qc,
+    double equat_core, double polar_core,
+    double equat_shell, double polar_shell,
+    double sld_core_shell, double sld_shell_solvent)
+{
+    const double qr_core = sqrt(square(equat_core*qab) + square(polar_core*qc));
+    const double si_core = sas_3j1x_x(qr_core);
+    const double volume_core = M_4PI_3*equat_core*equat_core*polar_core;
+    const double fq_core = si_core*volume_core*sld_core_shell;
 
-double Iqxy(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi);
+    const double qr_shell = sqrt(square(equat_shell*qab) + square(polar_shell*qc));
+    const double si_shell = sas_3j1x_x(qr_shell);
+    const double volume_shell = M_4PI_3*equat_shell*equat_shell*polar_shell;
+    const double fq_shell = si_shell*volume_shell*sld_shell_solvent;
 
+    return fq_core + fq_shell;
+}
 
-double form_volume(double radius_equat_core,
-                   double x_core,
-                   double thick_shell,
-                   double x_polar_shell)
+static double
+form_volume(double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell)
 {
     const double equat_shell = radius_equat_core + thick_shell;
@@ -37 +39 @@
 
 static double
-core_shell_ellipsoid_xt_kernel(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld)
+Iq(double q,
+    double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell,
+    double core_sld,
+    double shell_sld,
+    double solvent_sld)
 {
-    const double lolim = 0.0;
-    const double uplim = 1.0;
-
-
-    const double delpc = core_sld - shell_sld; //core - shell
-    const double delps = shell_sld - solvent_sld; //shell - solvent
-
+    const double sld_core_shell = core_sld - shell_sld;
+    const double sld_shell_solvent = shell_sld - solvent_sld;
 
     const double polar_core = radius_equat_core*x_core;
@@ -58 +55 @@
     const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell;
 
-    double summ = 0.0;   //initialize intergral
+    // translate from [-1, 1] => [0, 1]
+    const double m = 0.5;
+    const double b = 0.5;
+    double total = 0.0;     //initialize intergral
     for(int i=0;i<76;i++) {
-        double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim );
-        double yyy = gfn4(zi, radius_equat_core, polar_core, equat_shell,
-                          polar_shell, delpc, delps, q);
-        summ += Gauss76Wt[i] * yyy;
+        const double cos_theta = Gauss76Z[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,
+            radius_equat_core, polar_core,
+            equat_shell, polar_shell,
+            sld_core_shell, sld_shell_solvent);
+        total += Gauss76Wt[i] * fq * fq;
     }
-    summ *= 0.5*(uplim-lolim);
+    total *= m;
 
     // convert to [cm-1]
-    return 1.0e-4 * summ;
+    return 1.0e-4 * total;
 }
 
 static double
-core_shell_ellipsoid_xt_kernel_2d(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
+Iqxy(double qx, double qy,
+    double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell,
+    double core_sld,
+    double shell_sld,
+    double solvent_sld,
+    double theta,
+    double phi)
 {
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
-    const double sldcs = core_sld - shell_sld;
-    const double sldss = shell_sld- solvent_sld;
+    const double sld_core_shell = core_sld - shell_sld;
+    const double sld_shell_solvent = shell_sld - solvent_sld;
 
     const double polar_core = radius_equat_core*x_core;
@@ -93 +98 @@
     const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell;
 
-    // Call the IGOR library function to get the kernel:
-    // MUST use gfn4 not gf2 because of the def of params.
-    double answer = gfn4(cos_alpha,
-                  radius_equat_core,
-                  polar_core,
-                  equat_shell,
-                  polar_shell,
-                  sldcs,
-                  sldss,
-                  q);
+    double fq = _cs_ellipsoid_kernel(qab, qc,
+                  radius_equat_core, polar_core,
+                  equat_shell, polar_shell,
+                  sld_core_shell, sld_shell_solvent);
 
     //convert to [cm-1]
-    answer *= 1.0e-4;
-
-    return answer;
+    return 1.0e-4 * fq * fq;
 }
-
-double Iq(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld)
-{
-    double intensity = core_shell_ellipsoid_xt_kernel(q,
-           radius_equat_core,
-           x_core,
-           thick_shell,
-           x_polar_shell,
-           core_sld,
-           shell_sld,
-           solvent_sld);
-
-    return intensity;
-}
-
-
-double Iqxy(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
-{
-    double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy,
-                       radius_equat_core,
-                       x_core,
-                       thick_shell,
-                       x_polar_shell,
-                       core_sld,
-                       shell_sld,
-                       solvent_sld,
-                       theta,
-                       phi);
-
-    return intensity;
-}

sasmodels/models/core_shell_ellipsoid.py

@@ -141 +141 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c",
-          "core_shell_ellipsoid.c"]
+source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"]
 
 def ER(radius_equat_core, x_core, thick_shell, x_polar_shell):

sasmodels/models/core_shell_parallelepiped.c

@@ -1 @@
-double form_volume(double length_a, double length_b, double length_c, 
+double form_volume(double length_a, double length_b, double length_c,
                    double thick_rim_a, double thick_rim_b, double thick_rim_c);
 double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld,
@@ -9 +9 @@
             double thick_rim_c, double theta, double phi, double psi);
 
-double form_volume(double length_a, double length_b, double length_c, 
+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 + 
+    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;
@@ -34 +34 @@
     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
-    
+
     const double mu = 0.5 * q * length_b;
-    
+
     //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 + 
+    //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;
-    
+
     // Scale sides by B
     const double a_scaled = length_a / length_b;
@@ -101 +101 @@
             //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
             // This is the function called by csparallelepiped_analytical_2D_scaled,
-            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c        
-            
+            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c
+
             //  correct FF : sum of square of phase factors
             inner_total += Gauss76Wt[j] * form * form;
@@ -136 +136 @@
     double q, zhat, yhat, xhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
     // cspkernel in csparallelepiped recoded here
@@ -160 +163 @@
     double tc = length_a + 2.0*thick_rim_c;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    double siA = sas_sinx_x(0.5*q*length_a*xhat);
-    double siB = sas_sinx_x(0.5*q*length_b*yhat);
-    double siC = sas_sinx_x(0.5*q*length_c*zhat);
-    double siAt = sas_sinx_x(0.5*q*ta*xhat);
-    double siBt = sas_sinx_x(0.5*q*tb*yhat);
-    double siCt = sas_sinx_x(0.5*q*tc*zhat);
-    
+    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);
+
 
     // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
@@ -174 +177 @@
                + drB*siA*(siBt-siB)*siC*V2
                + drC*siA*siB*(siCt*siCt-siC)*V3);
-   
+
     return 1.0e-4 * f * f;
 }

sasmodels/models/cylinder.c

@@ -1 @@
-double form_volume(double radius, double length);
-double fq(double q, double sn, double cn,double radius, double length);
-double orient_avg_1D(double q, double radius, double length);
-double Iq(double q, double sld, double solvent_sld, double radius, double length);
-double Iqxy(double qx, double qy, double sld, double solvent_sld,
-    double radius, double length, double theta, double phi);
-
 #define INVALID(v) (v.radius<0 || v.length<0)
 
-double form_volume(double radius, double length)
+static double
+form_volume(double radius, double length)
 {
     return M_PI*radius*radius*length;
 }
 
-double fq(double q, double sn, double cn, double radius, double length)
+static double
+fq(double qab, double qc, double radius, double length)
 {
-    // precompute qr and qh to save time in the loop
-    const double qr = q*radius;
-    const double qh = q*0.5*length; 
-    return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn);
+    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length);
 }
 
-double orient_avg_1D(double q, double radius, double length)
+static double
+orient_avg_1D(double q, double radius, double length)
 {
     // translate a point in [-1,1] to a point in [0, pi/2]
     const double zm = M_PI_4;
-    const double zb = M_PI_4; 
+    const double zb = M_PI_4;
 
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
-        const double alpha = Gauss76Z[i]*zm + zb;
-        double sn, cn; // slots to hold sincos function output
-        // alpha(theta,phi) the projection of the cylinder on the detector plane
-        SINCOS(alpha, sn, cn);
-        total += Gauss76Wt[i] * square( fq(q, sn, cn, radius, length) ) * sn;
+        const double theta = Gauss76Z[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;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -39 +33 @@
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double solvent_sld,
@@ -49 +44 @@
 }
 
-
-double Iqxy(double qx, double qy,
+static double
+Iqxy(double qx, double qy,
     double sld,
     double solvent_sld,
@@ -60 +55 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
+
     const double s = (sld-solvent_sld) * form_volume(radius, length);
-    const double form = fq(q, sin_alpha, cos_alpha, radius, length);
+    const double form = fq(qab, qc, radius, length);
     return 1.0e-4 * square(s * form);
 }

sasmodels/models/ellipsoid.c

@@ -1 @@
-double form_volume(double radius_polar, double radius_equatorial);
-double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial);
-double Iqxy(double qx, double qy, double sld, double sld_solvent,
-    double radius_polar, double radius_equatorial, double theta, double phi);
-
-double form_volume(double radius_polar, double radius_equatorial)
+static double
+form_volume(double radius_polar, double radius_equatorial)
 {
     return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial;
 }
 
-double Iq(double q,
+static  double
+Iq(double q,
     double sld,
     double sld_solvent,
@@ -41 +38 @@
 }
 
-double Iqxy(double qx, double qy,
+static double
+Iqxy(double qx, double qy,
     double sld,
     double sld_solvent,
@@ -51 +49 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    const double r = sqrt(square(radius_equatorial*sin_alpha)
-                          + square(radius_polar*cos_alpha));
-    const double f = sas_3j1x_x(q*r);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
+
+    const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc));
+    const double f = sas_3j1x_x(qr);
     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
 
     return 1.0e-4 * square(f * s);
 }
-

sasmodels/models/elliptical_cylinder.c

@@ -1 @@
-double form_volume(double radius_minor, double r_ratio, double length);
-double Iq(double q, double radius_minor, double r_ratio, double length,
-          double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length,
-            double sld, double solvent_sld, double theta, double phi, double psi);
-
-
-double
+static double
 form_volume(double radius_minor, double r_ratio, double length)
 {
@@ -12 +5 @@
 }
 
-double
+static double
 Iq(double q, double radius_minor, double r_ratio, double length,
    double sld, double solvent_sld)
@@ -61 +54 @@
 
 
-double
+static double
 Iqxy(double qx, double qy,
      double radius_minor, double r_ratio, double length,
@@ -69 +62 @@
     double q, xhat, yhat, zhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
     // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2)
     // Given:    radius_major = r_ratio * radius_minor
-    const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat));
-    const double be = sas_2J1x_x(q*r);
-    const double si = sas_sinx_x(q*zhat*0.5*length);
+    const double qr = radius_minor*sqrt(square(r_ratio*qa) + square(qb));
+    const double be = sas_2J1x_x(qr);
+    const double si = sas_sinx_x(qc*0.5*length);
     const double Aq = be * si;
     const double delrho = sld - solvent_sld;

sasmodels/models/fcc_paracrystal.c

@@ -1 @@
-double form_volume(double radius);
-double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double dnn,
-    double d_factor, double radius,double sld, double solvent_sld,
-    double theta, double phi, double psi);
+static double
+fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
+{
+#if 0  // Equations as written in Matsuoka
+    const double a1 = ( qa + qb)/2.0;
+    const double a2 = (-qa + qc)/2.0;
+    const double a3 = (-qa + qb)/2.0;
+#else
+    const double a1 = ( qa + qb)/2.0;
+    const double a2 = ( qa + qc)/2.0;
+    const double a3 = ( qb + qc)/2.0;
+#endif
 
-double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
-double _FCCeval(double Theta, double Phi, double temp1, double temp3);
+    // Numerator: (1 - exp(a)^2)^3
+    //         => (-(exp(2a) - 1))^3
+    //         => -expm1(2a)^3
+    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
+    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double exp_arg = exp(arg);
+    const double Zq = -cube(expm1(2.0*arg))
+        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
+
+    return Zq;
+}
 
 
-double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
-
-        const double Da = d_factor*dnn;
-        const double temp1 = q*q*Da*Da;
-        const double temp3 = q*dnn;
-
-        double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
-        return(retVal);
+// occupied volume fraction calculated from lattice symmetry and sphere radius
+static double
+fcc_volume_fraction(double radius, double dnn)
+{
+    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn);
 }
 
-double _FCCeval(double Theta, double Phi, double temp1, double temp3) {
-
-        double result;
-        double sin_theta,cos_theta,sin_phi,cos_phi;
-        SINCOS(Theta, sin_theta, cos_theta);
-        SINCOS(Phi, sin_phi, cos_phi);
-
-        const double temp6 =  sin_theta;
-        const double temp7 =  sin_theta*sin_phi + cos_theta;
-        const double temp8 = -sin_theta*cos_phi + cos_theta;
-        const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi;
-
-        const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
-        result = cube(1.0-(temp10*temp10))*temp6
-            / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
-              * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
-              * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10));
-
-        return (result);
-}
-
-double form_volume(double radius){
+static double
+form_volume(double radius)
+{
     return sphere_volume(radius);
 }
 
 
-double Iq(double q, double dnn,
+static double Iq(double q, double dnn,
   double d_factor, double radius,
-  double sld, double solvent_sld){
+  double sld, double solvent_sld)
+{
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    const double phi_m = M_PI;
+    const double phi_b = M_PI;
+    // translate a point in [-1,1] to a point in [0, pi]
+    const double theta_m = M_PI_2;
+    const double theta_b = M_PI_2;
 
-        //Volume fraction calculated from lattice symmetry and sphere radius
-        const double s1 = dnn*sqrt(2.0);
-        const double latticescale = 4.0*sphere_volume(radius/s1);
+    double outer_sum = 0.0;
+    for(int i=0; i<150; i++) {
+        double inner_sum = 0.0;
+        const double theta = Gauss150Z[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;
+            double sin_phi, cos_phi;
+            SINCOS(phi, sin_phi, cos_phi);
+            const double qa = qab*cos_phi;
+            const double qb = qab*sin_phi;
+            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * form;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m;
+    const double Zq = outer_sum/(4.0*M_PI);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
 
-    const double va = 0.0;
-    const double vb = 2.0*M_PI;
-    const double vaj = 0.0;
-    const double vbj = M_PI;
-
-    double summ = 0.0;
-    double answer = 0.0;
-        for(int i=0; i<150; i++) {
-                //setup inner integral over the ellipsoidal cross-section
-                double summj=0.0;
-                const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi
-                for(int j=0;j<150;j++) {
-                        //20 gauss points for the inner integral
-                        double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta
-                        double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi);
-                        summj += yyy;
-                }
-                //now calculate the value of the inner integral
-                double answer = (vbj-vaj)/2.0*summj;
-
-                //now calculate outer integral
-                summ = summ+(Gauss150Wt[i] * answer);
-        }               //final scaling is done at the end of the function, after the NT_FP64 case
-
-        answer = (vb-va)/2.0*summ;
-        answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
-
-    return answer;
+    return fcc_volume_fraction(radius, dnn) * Pq * Zq;
+}
 
 
-}
-
-double Iqxy(double qx, double qy,
+static double Iqxy(double qx, double qy,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld,
@@ -92 +89 @@
     double q, zhat, yhat, xhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
-    const double a1 = yhat + xhat;
-    const double a2 = xhat + zhat;
-    const double a3 = yhat + zhat;
-    const double qd = 0.5*q*dnn;
-    const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3);
-    const double tanh_qd = tanh(arg);
-    const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a2)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*a3)/cosh_qd);
-
-    //if (isnan(Zq)) printf("q:(%g,%g) qd: %g a1: %g a2: %g a3: %g arg: %g\n", qx, qy, qd, a1, a2, a3, arg);
-
-    const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq;
-    //the occupied volume of the lattice
-    const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn);
-    return lattice_scale * Fq;
+    q = sqrt(qa*qa + qb*qb + qc*qc);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor);
+    return fcc_volume_fraction(radius, dnn) * Pq * Zq;
 }

sasmodels/models/hollow_cylinder.c

@@ -1 +1 @@
 double form_volume(double radius, double thickness, double length);
 double Iq(double q, double radius, double thickness, double length, double sld,
-        double solvent_sld);
+    double solvent_sld);
 double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld,
-        double solvent_sld, double theta, double phi);
+    double solvent_sld, double theta, double phi);
 
 //#define INVALID(v) (v.radius_core >= v.radius)
@@ -14 +14 @@
 }
 
-
 static double
-_hollow_cylinder_kernel(double q,
-    double radius, double thickness, double length, double sin_val, double cos_val)
+_fq(double qab, double qc,
+    double radius, double thickness, double length)
 {
-    const double qs = q*sin_val;
-    const double lam1 = sas_2J1x_x((radius+thickness)*qs);
-    const double lam2 = sas_2J1x_x(radius*qs);
+    const double lam1 = sas_2J1x_x((radius+thickness)*qab);
+    const double lam2 = sas_2J1x_x(radius*qab);
     const double gamma_sq = square(radius/(radius+thickness));
-    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qs)
-    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs)
-    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00
-    const double t2 = sas_sinx_x(0.5*q*length*cos_val);
+    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab)
+    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab)
+    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00
+    const double t2 = sas_sinx_x(0.5*length*qc);
     return psi*t2;
 }
@@ -43 +41 @@
 {
     const double lower = 0.0;
-    const double upper = 1.0;           //limits of numerical integral
+    const double upper = 1.0;        //limits of numerical integral
 
-    double summ = 0.0;                  //initialize intergral
+    double summ = 0.0;            //initialize intergral
     for (int i=0;i<76;i++) {
-        const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper );
-        const double sin_val = sqrt(1.0 - cos_val*cos_val);
-        const double inter = _hollow_cylinder_kernel(q, radius, thickness, length,
-                                                     sin_val, cos_val);
-        summ += Gauss76Wt[i] * inter * inter;
+        const double cos_theta = 0.5*( Gauss76Z[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;
     }
 
@@ -66 +64 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length,
-        sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
+
+    const double form = _fq(qab, qc, radius, thickness, length);
 
     const double vol = form_volume(radius, thickness, length);
-    return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol);
+    return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol);
 }
-

sasmodels/models/parallelepiped.c

@@ -20 +20 @@
 {
     const double mu = 0.5 * q * length_b;
-    
+
     // Scale sides by B
     const double a_scaled = length_a / length_b;
     const double c_scaled = length_c / length_b;
-        
+
     // outer integral (with gauss points), integration limits = 0, 1
     double outer_total = 0; //initialize integral
@@ -69 +69 @@
     double q, xhat, yhat, zhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
-    const double siA = sas_sinx_x(0.5*length_a*q*xhat);
-    const double siB = sas_sinx_x(0.5*length_b*q*yhat);
-    const double siC = sas_sinx_x(0.5*length_c*q*zhat);
+    const double siA = sas_sinx_x(0.5*length_a*qa);
+    const double siB = sas_sinx_x(0.5*length_b*qb);
+    const double siC = sas_sinx_x(0.5*length_c*qc);
     const double V = form_volume(length_a, length_b, length_c);
     const double drho = (sld - solvent_sld);

sasmodels/models/sc_paracrystal.c

@@ -1 @@
-double form_volume(double radius);
+static double
+sc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
+{
+    const double a1 = qa;
+    const double a2 = qb;
+    const double a3 = qc;
 
-double Iq(double q,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld);
+    // Numerator: (1 - exp(a)^2)^3
+    //         => (-(exp(2a) - 1))^3
+    //         => -expm1(2a)^3
+    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
+    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double exp_arg = exp(arg);
+    const double Zq = -cube(expm1(2.0*arg))
+        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
 
-double Iqxy(double qx, double qy,
-            double dnn,
-            double d_factor,
-            double radius,
-            double sphere_sld,
-            double solvent_sld,
-            double theta,
-            double phi,
-            double psi);
+    return Zq;
+}
 
-double form_volume(double radius)
+// occupied volume fraction calculated from lattice symmetry and sphere radius
+static double
+sc_volume_fraction(double radius, double dnn)
+{
+    return sphere_volume(radius/dnn);
+}
+
+static double
+form_volume(double radius)
 {
     return sphere_volume(radius);
 }
 
-static double
-sc_eval(double theta, double phi, double temp3, double temp4, double temp5)
+
+static double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld)
 {
-    double cnt, snt;
-    SINCOS(theta, cnt, snt);
+    // translate a point in [-1,1] to a point in [0, 2 pi]
+    const double phi_m = M_PI_4;
+    const double phi_b = M_PI_4;
+    // translate a point in [-1,1] to a point in [0, pi]
+    const double theta_m = M_PI_4;
+    const double theta_b = M_PI_4;
 
-    double cnp, snp;
-    SINCOS(phi, cnp, snp);
 
-        double temp6 = snt;
-        double temp7 = -1.0*temp3*snt*cnp;
-        double temp8 = temp3*snt*snp;
-        double temp9 = temp3*cnt;
-        double result = temp6/((1.0-temp4*cos((temp7))+temp5)*
-                               (1.0-temp4*cos((temp8))+temp5)*
-                               (1.0-temp4*cos((temp9))+temp5));
-        return (result);
+    double outer_sum = 0.0;
+    for(int i=0; i<150; i++) {
+        double inner_sum = 0.0;
+        const double theta = Gauss150Z[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;
+            double sin_phi, cos_phi;
+            SINCOS(phi, sin_phi, cos_phi);
+            const double qa = qab*cos_phi;
+            const double qb = qab*sin_phi;
+            const double form = sc_Zq(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * form;
+        }
+        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
+        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+    }
+    outer_sum *= theta_m;
+    const double Zq = outer_sum/M_PI_2;
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+
+    return sc_volume_fraction(radius, dnn) * Pq * Zq;
 }
 
-static double
-sc_integrand(double dnn, double d_factor, double qq, double xx, double yy)
-{
-    //Function to calculate integrand values for simple cubic structure
 
-        double da = d_factor*dnn;
-        double temp1 = qq*qq*da*da;
-        double temp2 = cube(-expm1(-temp1));
-        double temp3 = qq*dnn;
-        double temp4 = 2.0*exp(-0.5*temp1);
-        double temp5 = exp(-1.0*temp1);
-
-        double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2;
-
-        return(integrand);
-}
-
-double Iq(double q,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld)
-{
-        const double va = 0.0;
-        const double vb = M_PI_2; //orientation average, outer integral
-
-    double summ=0.0;
-    double answer=0.0;
-
-        for(int i=0;i<150;i++) {
-                //setup inner integral over the ellipsoidal cross-section
-                double summj=0.0;
-                double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;
-                for(int j=0;j<150;j++) {
-                        //150 gauss points for the inner integral
-                        double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0;
-                        double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij);
-                        summj += tmp;
-                }
-                //now calculate the value of the inner integral
-                answer = (vb-va)/2.0*summj;
-
-                //now calculate outer integral
-                double tmp = Gauss150Wt[i] * answer;
-                summ += tmp;
-        }               //final scaling is done at the end of the function, after the NT_FP64 case
-
-        answer = (vb-va)/2.0*summ;
-
-        //Volume fraction calculated from lattice symmetry and sphere radius
-        // NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3
-        const double latticeScale = sphere_volume(radius/dnn);
-
-        answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale;
-
-        return answer;
-}
-
-double Iqxy(double qx, double qy,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld,
-          double theta,
-          double phi,
-          double psi)
+static double Iqxy(double qx, double qy,
+    double dnn, double d_factor, double radius,
+    double sld, double solvent_sld,
+    double theta, double phi, double psi)
 {
     double q, zhat, yhat, xhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
-    const double qd = q*dnn;
-    const double arg = 0.5*square(qd*d_factor);
-    const double tanh_qd = tanh(arg);
-    const double cosh_qd = cosh(arg);
-    const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*yhat)/cosh_qd)
-                    * tanh_qd/(1. - cos(qd*xhat)/cosh_qd);
-
-    const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq;
-    //the occupied volume of the lattice
-    const double lattice_scale = sphere_volume(radius/dnn);
-    return lattice_scale * Fq;
+    q = sqrt(qa*qa + qb*qb + qc*qc);
+    const double Pq = sphere_form(q, radius, sld, solvent_sld);
+    const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor);
+    return sc_volume_fraction(radius, dnn) * Pq * Zq;
 }

sasmodels/models/sc_paracrystal.py

@@ -161 +161 @@
     [{'theta': 10.0, 'phi': 20, 'psi': 30.0}, (0.023, 0.045), 0.0177333171285],
     ]
-
-

sasmodels/models/stacked_disks.c

@@ -1 +1 @@
 static double stacked_disks_kernel(
-    double q,
+    double qab,
+    double qc,
     double halfheight,
     double thick_layer,
@@ -9 +10 @@
     double layer_sld,
     double solvent_sld,
-    double sin_alpha,
-    double cos_alpha,
     double d)
 
@@ -20 +19 @@
     // zi is the dummy variable for the integration (x in Feigin's notation)
 
-    const double besarg1 = q*radius*sin_alpha;
-    //const double besarg2 = q*radius*sin_alpha;
+    const double besarg1 = radius*qab;
+    //const double besarg2 = radius*qab;
 
-    const double sinarg1 = q*halfheight*cos_alpha;
-    const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha;
+    const double sinarg1 = halfheight*qc;
+    const double sinarg2 = (halfheight+thick_layer)*qc;
 
     const double be1 = sas_2J1x_x(besarg1);
@@ -43 +42 @@
 
     // loop for the structure factor S(q)
-    double qd_cos_alpha = q*d*cos_alpha;
+    double qd_cos_alpha = d*qc;
     //d*cos_alpha is the projection of d onto q (in other words the component
     //of d that is parallel to q.
@@ -84 +83 @@
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(zi, sin_alpha, cos_alpha);
-        double yyy = stacked_disks_kernel(q,
+        double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha,
                            halfheight,
                            thick_layer,
@@ -93 +92 @@
                            layer_sld,
                            solvent_sld,
-                           sin_alpha,
-                           cos_alpha,
                            d);
         summ += Gauss76Wt[i] * yyy * sin_alpha;
@@ -152 +149 @@
     double phi)
 {
-    int n_stacking = (int)(fp_n_stacking + 0.5);
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
+    int n_stacking = (int)(fp_n_stacking + 0.5);
     double d = 2.0 * thick_layer + thick_core;
     double halfheight = 0.5*thick_core;
-    double answer = stacked_disks_kernel(q,
+    double answer = stacked_disks_kernel(qab, qc,
                      halfheight,
                      thick_layer,
@@ -167 +166 @@
                      layer_sld,
                      solvent_sld,
-                     sin_alpha,
-                     cos_alpha,
                      d);
 

sasmodels/models/triaxial_ellipsoid.c

@@ -1 @@
-double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar);
-double Iq(double q, double sld, double sld_solvent,
-    double radius_equat_minor, double radius_equat_major, double radius_polar);
-double Iqxy(double qx, double qy, double sld, double sld_solvent,
-    double radius_equat_minor, double radius_equat_major, double radius_polar, double theta, double phi, double psi);
-
 //#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar)
 
-
-double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar)
+static double
+form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar)
 {
     return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar;
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double sld_solvent,
@@ -45 +40 @@
     // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi
     const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI);
-    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);
-    return 1.0e-4 * s * s * fqsq;
+    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar);
+    const double drho = (sld - sld_solvent);
+    return 1.0e-4 * square(vol*drho) * fqsq;
 }
 
-double Iqxy(double qx, double qy,
+static double
+Iqxy(double qx, double qy,
     double sld,
     double sld_solvent,
@@ -61 +58 @@
     double q, xhat, yhat, zhat;
     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    const double qa = q*xhat;
+    const double qb = q*yhat;
+    const double qc = q*zhat;
 
-    const double r = sqrt(square(radius_equat_minor*xhat)
-                          + square(radius_equat_major*yhat)
-                          + square(radius_polar*zhat));
-    const double fq = sas_3j1x_x(q*r);
-    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);
+    const double qr = sqrt(square(radius_equat_minor*qa)
+                           + square(radius_equat_major*qb)
+                           + square(radius_polar*qc));
+    const double fq = sas_3j1x_x(qr);
+    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar);
+    const double drho = (sld - sld_solvent);
 
-    return 1.0e-4 * square(s * fq);
+    return 1.0e-4 * square(vol * drho * fq);
 }
-

sasmodels/weights.py

@@ -55 +55 @@
         """
         sigma = self.width * center if relative else self.width
+        if not relative:
+            # For orientation, the jitter is relative to 0 not the angle
+            #center = 0
+            pass
         if sigma == 0 or self.npts < 2:
             if lb <= center <= ub: