Changeset ac60a39 – SasView

sasmodels/product.py

-                      rce99754
+                      rac60a39
     def __init__(self, model_info, P, S):
         # type: (ModelInfo, KernelModel, KernelModel) -> None
+        #: Combined info plock for the product model
         self.info = model_info
+        #: Form factor modelling individual particles.
         self.P = P
+        #: Structure factor modelling interaction between particles.
         self.S = S
+        self.dtype = P.dtype
+        #: Model precision. This is not really relevant, since it is the
+        #: individual P and S models that control the effective dtype,
+        #: converting the q-vectors to the correct type when the kernels
+        #: for each are created. Ideally this should be set to the more
+        #: precise type to avoid loss of precision, but precision in q is
+        #: not critical (single is good enough for our purposes), so it just
+        #: uses the precision of the form factor.
+        self.dtype = P.dtype  # type: np.dtype
     def make_kernel(self, q_vectors):

LICENSE.txt

-                      rf68e2a5
+                      r5f8b72b
 All rights reserved.
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+. Redistributions of source code must retain the above copyright notice,
+this list of conditions and the following disclaimer.
+. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+. Redistributions in binary form must reproduce the above copyright notice,
+this list of conditions and the following disclaimer in the documentation
+and/or other materials provided with the distribution.
+. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+. Neither the name of the copyright holder nor the names of its contributors
+may be used to endorse or promote products derived from this software without
+specific prior written permission.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.

doc/_extensions/dollarmath.py

r103ea45	rdd4d95d
12	12	import re
13	13
14		_dollar = re.compile(r"(?:^\|(?<=\s\|[~~(]))[$]([^\n]*?)(?<![\\])[$](?:$\|(?=\s\|[.,;)\\~~]))")
	14	_dollar = re.compile(r"(?:^\|(?<=\s\|[-(]))[$]([^\n]*?)(?<![\\])[$](?:$\|(?=\s\|[-.,;:?\\)]))")
15	15	_notdollar = re.compile(r"\\[$]")
16	16

doc/developer/calculator.rst

-                      r870a2f4
+                      r2c108a3
 This document describes the layer between the form factor kernels and the
 model calculator which implements the polydispersity and magnetic SLD
+model calculator which implements the dispersity and magnetic SLD
 calculations.  There are three separate implementations of this layer,
 :mod:`kernelcl` for OpenCL, which operates on a single Q value at a time,
 …
 Each implementation provides three different calls *Iq*, *Iqxy* and *Imagnetic*
+for 1-D, 2-D and 2-D magnetic kernels respectively.
+The C code is defined in *kernel_iq.c* and *kernel_iq.cl* for DLL and OpenCL
+respectively.  The kernel call looks as follows::
+for 1-D, 2-D and 2-D magnetic kernels respectively. The C code is defined
+in *kernel_iq.c*, with the minor differences between OpenCL and DLL handled
+by #ifdef statements.
+The kernel call looks as follows::
   kernel void KERNEL_NAME(
       int nq,                  // Number of q values in the q vector
       int pd_start,            // Starting position in the polydispersity loop
       int pd_stop,             // Ending position in the polydispersity loop
       ProblemDetails *details, // Polydispersity info
+      int pd_start,            // Starting position in the dispersity loop
+      int pd_stop,             // Ending position in the dispersity loop
+      ProblemDetails *details, // dispersity info
       double *values,          // Value and weights vector
       double *q,               // q or (qx,qy) vector
       double *result,          // returned I(q), with result[nq] = pd_weight
       double cutoff)           // polydispersity weight cutoff
+      double cutoff)           // dispersity weight cutoff
 The details for OpenCL and the python loop are slightly different, but these
 …
 *nq* indicates the number of q values that will be calculated.
 The *pd_start* and *pd_stop* parameters set the range of the polydispersity
 loop to compute for the current kernel call.   Give a polydispersity
+The *pd_start* and *pd_stop* parameters set the range of the dispersity
+loop to compute for the current kernel call.   Give a dispersity
 calculation with 30 weights for length and 30 weights for radius for example,
 there are a total of 900 calls to the form factor required to compute the
 …
 the length index to 3 and the radius index to 10 for a position of 3*30+10=100,
 and could then proceed to position 200.  This allows us to interrupt the
 calculation in the middle of a long polydispersity loop without having to
+calculation in the middle of a long dispersity loop without having to
 do special tricks with the C code.  More importantly, it stops the OpenCL
 kernel in a reasonable time; because the GPU is used by the operating
 …
 The *ProblemDetails* structure is a direct map of the
 :class:`details.CallDetails` buffer.  This indicates which parameters are
 polydisperse, and where in the values vector the values and weights can be
 found.  For each polydisperse parameter there is a parameter id, the length
 of the polydispersity loop for that parameter, the offset of the parameter
+:class:`details.CallDetails` buffer.  This indicates which parameters have
+dispersity, and where in the values vector the values and weights can be
+found.  For each parameter with dispersity there is a parameter id, the length
+of the dispersity loop for that parameter, the offset of the parameter
 values in the pd value and pd weight vectors and the 'stride' from one index
 to the next, which is used to translate between the position in the
 polydispersity loop and the particular parameter indices.  The *num_eval*
 field is the total size of the polydispersity loop.  *num_weights* is the
+dispersity loop and the particular parameter indices.  The *num_eval*
+field is the total size of the dispersity loop.  *num_weights* is the
 number of elements in the pd value and pd weight vectors.  *num_active* is
 the number of non-trivial pd loops (polydisperse parameters should be ordered
 by decreasing pd vector length, with a length of 1 meaning no polydispersity).
+the number of non-trivial pd loops (parameters with dispersity should be ordered
+by decreasing pd vector length, with a length of 1 meaning no dispersity).
 Oriented objects in 2-D need a cos(theta) spherical correction on the angular
 variation in order to preserve the 'surface area' of the weight distribution.
 …
 *(Mx, My, Mz)*.  Sample magnetization is translated from *(M, theta, phi)*
 to *(Mx, My, Mz)* before the kernel is called.   After the fixed values comes
 the pd value vector, with the polydispersity values for each parameter
+the pd value vector, with the dispersity values for each parameter
 stacked one after the other.  The order isn't important since the location
 for each parameter is stored in the *pd_offset* field of the *ProblemDetails*
 …
 values, the pd weight vector is stored, with the same configuration as the
 pd value vector.  Note that the pd vectors can contain values that are not
 in the polydispersity loop; this is used by :class:`mixture.MixtureKernel`
+in the dispersity loop; this is used by :class:`mixture.MixtureKernel`
 to make it easier to call the various component kernels.
 …
 The *results* vector contains one slot for each of the *nq* values, plus
 one extra slot at the end for the current polydisperse normalization.  This
+is required when the polydispersity loop is broken across several kernel
 calls.
+one extra slot at the end for the weight normalization accumulated across
+all points in the dispersity mesh.  This is required when the dispersity
+loop is broken across several kernel calls.
 *cutoff* is a importance cutoff so that points which contribute negligibly
 …
 - USE_OPENCL is defined if running in opencl
 - MAX_PD is the maximum depth of the polydispersity loop [model specific]
+- MAX_PD is the maximum depth of the dispersity loop [model specific]
 - NUM_PARS is the number of parameter values in the kernel.  This may be
   more than the number of parameters if some of the parameters are vector
   values.
 - NUM_VALUES is the number of fixed values, which defines the offset in the
   value list to the polydisperse value and weight vectors.
+  value list to the dispersity value and weight vectors.
 - NUM_MAGNETIC is the number of magnetic SLDs
 - MAGNETIC_PARS is a comma separated list of the magnetic SLDs, indicating
   their locations in the values vector.
-- MAGNETIC_PAR0 to MAGNETIC_PAR2 are the first three magnetic parameter ids
-  so we can hard code the setting of magnetic values if there are only a
-  few of them.
 - KERNEL_NAME is the name of the function being declared
 - PARAMETER_TABLE is the declaration of the parameters to the kernel:
 …
     Cylinder2D::
         #define CALL_IQ(q, i, var) Iqxy(q[2*i], q[2*i+1], \
+        #define CALL_IQ(q, i, var) Iqxy(qa, qc, \
         var.length, \
         var.radius, \
         var.sld, \
+        var.sld_solvent, \
+        var.theta, \
+        var.phi)
+        var.sld_solvent)
 - CALL_VOLUME(var) is similar, but for calling the form volume::
 …
         #define INVALID(var) constrained(var.p1, var.p2, var.p3)
 Our design supports a limited number of polydispersity loops, wherein
 we need to cycle through the values of the polydispersity, calculate
+Our design supports a limited number of dispersity loops, wherein
+we need to cycle through the values of the dispersity, calculate
 the I(q, p) for each combination of parameters, and perform a normalized
 weighted sum across all the weights.  Parameters may be passed to the
 underlying calculation engine as scalars or vectors, but the polydispersity
+underlying calculation engine as scalars or vectors, but the dispersity
 calculator treats the parameter set as one long vector.
 Let's assume we have 8 parameters in the model, with two polydisperse.  Since
 this is a 1-D model the orientation parameters won't be used::
+Let's assume we have 8 parameters in the model, two of which allow dispersity.
+Since this is a 1-D model the orientation parameters won't be used::
 : scale        {scl = constant}
 …
 : radius       {r = vector of 10pts}
 : length       {l = vector of 30pts}
 : sld          {s = constant/(radius**2*length)}
+: sld          {s1 = constant/(radius**2*length)}
 : sld_solvent  {s2 = constant}
 : theta        {not used}
 …
 This generates the following call to the kernel.  Note that parameters 4 and
 are treated as polydisperse even though they are not --- this is because
+are treated as having dispersity even though they don't --- this is because
 it is harmless to do so and simplifies the looping code::
 …
     NUM_MAGNETIC = 2      // two parameters might be magnetic
     MAGNETIC_PARS = 4, 5  // they are sld and sld_solvent
-    MAGNETIC_PAR0 = 4     // sld index
-    MAGNETIC_PAR1 = 5     // solvent index
     details {
 …
         pd_offset = {10, 0, 31, 32}   // *length* starts at index 10 in weights
         pd_stride = {1, 30, 300, 300} // cumulative product of pd length
         num_eval = 300   // 300 values in the polydispersity loop
+        num_eval = 300   // 300 values in the dispersity loop
         num_weights = 42 // 42 values in the pd vector
         num_active = 2   // only the first two pd are active
 …
     values = { scl, bkg,                                  // universal
                r, l, s, s2, theta, phi,                   // kernel pars
+               r, l, s1, s2, theta, phi,                  // kernel pars
                in spin, out spin, spin angle,             // applied magnetism
                mx s, my s, mz s, mx s2, my s2, mz s2,     // magnetic slds
+               mx s1, my s1, mz s1, mx s2, my s2, mz s2,  // magnetic slds
                r0, .., r9, l0, .., l29, s, s2,            // pd values
                r0, .., r9, l0, .., l29, s, s2}            // pd weights
 …
     result = {r1, ..., r130, pd_norm, x }
 The polydisperse parameters are stored in as an array of parameter
 indices, one for each polydisperse parameter, stored in pd_par[n].
+Non-polydisperse parameters do not appear in this array. Each polydisperse
+The dispersity parameters are stored in as an array of parameter
+indices, one for each parameter, stored in pd_par[n]. Parameters which do
+not support dispersity do not appear in this array. Each dispersity
 parameter has a weight vector whose length is stored in pd_length[n].
 The weights are stored in a contiguous vector of weights for all
 …
 in pd_offset[n].  The values corresponding to the weights are stored
 together in a separate weights[] vector, with offset stored in
 par_offset[pd_par[n]]. Polydisperse parameters should be stored in
+par_offset[pd_par[n]]. Dispersity parameters should be stored in
 decreasing order of length for highest efficiency.
 We limit the number of polydisperse dimensions to MAX_PD (currently 4),
 though some models may have fewer if they have fewer polydisperse
+We limit the number of dispersity dimensions to MAX_PD (currently 4),
+though some models may have fewer if they have fewer dispersity
 parameters.  The main reason for the limit is to reduce code size.
+Each additional polydisperse parameter requires a separate polydispersity
+loop.  If more than 4 levels of polydispersity are needed, then kernel_iq.c
+and kernel_iq.cl will need to be extended.
+Each additional dispersity parameter requires a separate dispersity
+loop.  If more than 4 levels of dispersity are needed, then we need to
+switch to a monte carlo importance sampling algorithm with better
+performance for high-dimensional integrals.
 Constraints between parameters are not supported.  Instead users will
 …
 theta since the polar coordinates normalization is tied to this parameter.
 If there is no polydispersity we pretend that it is polydisperisty with one
+parameter, pd_start=0 and pd_stop=1.  We may or may not short circuit the
 calculation in this case, depending on how much time it saves.
+If there is no dispersity we pretend that we have a disperisty mesh over
+a single parameter with a single point in the distribution, giving
+pd_start=0 and pd_stop=1.
 The problem details structure could be allocated and sent in as an integer
 array using the read-only flag.  This would allow us to copy it once per fit
 along with the weights vector, since features such as the number of
+polydisperity elements per pd parameter or the coordinated won't change
+between function evaluations.  A new parameter vector must be sent for
+each I(q) evaluation.  This is not currently implemented, and would require
 some resturcturing of the :class:`sasview_model.SasviewModel` interface.
 The results array will be initialized to zero for polydispersity loop
+disperity points per pd parameter won't change between function evaluations.
+A new parameter vector must be sent for each I(q) evaluation.  This is
+not currently implemented, and would require some resturcturing of
+the :class:`sasview_model.SasviewModel` interface.
+The results array will be initialized to zero for dispersity loop
 entry zero, and preserved between calls to [start, stop] so that the
 results accumulate by the time the loop has completed.  Background and
 …
 This will require accumulated error for each I(q) value to be preserved
 between kernel calls to implement this fully.  The kernel_iq.c code, which
 loops over q for each parameter set in the polydispersity loop, will need
 also need the accumalation vector.
+between kernel calls to implement this fully.  The *kernel_iq.c* code, which
+loops over q for each parameter set in the dispersity loop, will also need
+the accumulation vector.

doc/developer/overview.rst

-                      r870a2f4
+                      r3d40839
 :ref:`Calculator_Interface`
+.. _orientation_developer:
+Orientation and Numerical Integration
+-------------------------------------
+For 2d data from oriented anisotropic particles, the mean particle
+orientation is defined by angles $\theta$, $\phi$ and $\Psi$, which are not
+in general the same as similarly named angles in many form factors. The
+wikipedia page on Euler angles (https://en.wikipedia.org/wiki/Euler_angles)
+lists the different conventions available. To quote: "Different authors may
+use different sets of rotation axes to define Euler angles, or different
+names for the same angles. Therefore, any discussion employing Euler angles
+should always be preceded by their definition."
+We are using the $z$-$y$-$z$ convention with extrinsic rotations
+$\Psi$-$\theta$-$\phi$ for the particle orientation and $x$-$y$-$z$
+convention with extrinsic rotations $\Psi$-$\theta$-$\phi$ for jitter, with
+jitter applied before particle orientation.
+For numerical integration within form factors etc. sasmodels is mostly using
+Gaussian quadrature with 20, 76 or 150 points depending on the model. It also
+makes use of symmetries such as calculating only over one quadrant rather
+than the whole sphere. There is often a U-substitution replacing $\theta$
+with $cos(\theta)$ which changes the limits of integration from 0 to $\pi/2$
+to 0 to 1 and also conveniently absorbs the $sin(\theta)$ scale factor in the
+integration. This can cause confusion if checking equations to include in a
+paper or thesis! Most models use the same core kernel code expressed in terms
+of the rotated view ($q_a$, $q_b$, $q_c$) for both the 1D and the 2D models,
+but there are also historical quirks such as the parallelepiped model, which
+has a useless transformation representing $j_0(a q_a)$ as $j_0(b q_a a/b)$.
+Useful testing routines include:
+:mod:`asymint` a direct implementation of the surface integral for certain
+models to get a more trusted value for the 1D integral using a
+reimplementation of the 2D kernel in python and mpmath (which computes math
+functions to arbitrary precision). It uses $\theta$ ranging from 0 to $\pi$
+and $\phi$ ranging from 0 to $2\pi$. It perhaps would benefit from including
+the U-substitution for $\theta$.
+:mod:`check1d` uses sasmodels 1D integration and compares that with a
+rectangle distribution in $\theta$ and $\phi$, with $\theta$ limits set to
+$\pm 90/\sqrt(3)$ and $\phi$ limits set to $\pm 180/\sqrt(3)$ [The rectangle
+weight function uses the fact that the distribution width column is labelled
+sigma to decide that the 1-$\sigma$ width of a rectangular distribution needs to
+be multiplied by $\sqrt(3)$ to get the corresponding gaussian equivalent, or
+similar reasoning.] This should rotate the sample through the entire
+$\theta$-$\phi$ surface according to the pattern that you see in jitter.py when
+you modify it to use 'rectangle' rather than 'gaussian' for its distribution
+without changing the viewing angle. In order to match the 1-D pattern for
+an arbitrary viewing angle on triaxial shapes, we need to integrate
+over $\theta$, $\phi$ and $\Psi$.
+When computing the dispersity integral, weights are scaled by
+$|\cos(\delta \theta)|$ to account for the points in $\phi$ getting closer
+together as $\delta \theta$ increases.
+[This will probably change so that instead of adjusting the weights, we will
+adjust $\delta\theta$-$\delta\phi$ mesh so that the point density in
+$\delta\phi$ is lower at larger $\delta\theta$. The flag USE_SCALED_PHI in
+*kernel_iq.c* selects an alternative algorithm.]
+The integrated dispersion is computed at a set of $(qx, qy)$ points $(q
+\cos(\alpha), q \sin(\alpha))$ at some angle $\alpha$ (currently angle=0) for
+each $q$ used in the 1-D integration. The individual $q$ points should be
+equivalent to asymint rect-n when the viewing angle is set to
+$(\theta,\phi,\Psi) = (90, 0, 0)$. Such tests can help to validate that 2d
+models are consistent with 1d models.
+:mod:`sascomp -sphere=n` uses the same rectangular distribution as check1d to
+compute the pattern of the $q_x$-$q_y$ grid.
+The :mod:`sascomp` utility can be used for 2d as well as 1d calculations to
+compare results for two sets of parameters or processor types, for example
+these two oriented cylinders here should be equivalent.
+:mod:`\./sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10`
 Testing

doc/genapi.py

-                      r2e66ef5
+                      r706f466
     #('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'),
 …
     ('sasview_model', 'Sasview interface'),
     ('sesans', 'SESANS calculation routines'),
+    ('special', 'Special functions library'),
     ('weights', 'Distribution functions'),
+]

doc/guide/index.rst

rc0d7ab3	rda5536f
13	13	resolution.rst
14	14	magnetism/magnetism.rst
	15	orientation/orientation.rst
15	16	sesans/sans_to_sesans.rst
16	17	sesans/sesans_fitting.rst

doc/guide/magnetism/magnetism.rst

-                      r1f058ea
+                      r4f5afc9
 Models which define a scattering length density parameter can be evaluated
  as magnetic models. In general, the scattering length density (SLD =
  $\beta$) in each region where the SLD is uniform, is a combination of the
  nuclear and magnetic SLDs and, for polarised neutrons, also depends on the
  spin states of the neutrons.
+as magnetic models. In general, the scattering length density (SLD =
+$\beta$) in each region where the SLD is uniform, is a combination of the
+nuclear and magnetic SLDs and, for polarised neutrons, also depends on the
+spin states of the neutrons.
 For magnetic scattering, only the magnetization component $\mathbf{M_\perp}$
 perpendicular to the scattering vector $\mathbf{Q}$ contributes to the magnetic
 scattering length.
 .. figure::
 …
 is the Pauli spin.
 Assuming that incident neutrons are polarized parallel (+) and anti-parallel (-)
+to the $x'$ axis, the possible spin states after the sample are then
+Assuming that incident neutrons are polarized parallel $(+)$ and anti-parallel
+$(-)$ to the $x'$ axis, the possible spin states after the sample are then:
+No spin-flips (+ +) and (- -)
+* Non spin-flip $(+ +)$ and $(- -)$
+Spin-flips    (+ -) and (- +)
+* Spin-flip $(+ -)$ and $(- +)$
+Each measurement is an incoherent mixture of these spin states based on the
+fraction of $+$ neutrons before ($u_i$) and after ($u_f$) the sample,
+with weighting:
+.. math::
+    -- &= ((1-u_i)(1-u_f))^{1/4} \\
+    -+ &= ((1-u_i)(u_f))^{1/4} \\
+    +- &= ((u_i)(1-u_f))^{1/4} \\
+    ++ &= ((u_i)(u_f))^{1/4}
+Ideally the experiment would measure the pure spin states independently and
+perform a simultaneous analysis of the four states, tying all the model
+parameters together except $u_i$ and $u_f$.
 .. figure::
 …
 $\phi$ and $\theta_{up}$, respectively, then, depending on the spin state of the
 neutrons, the scattering length densities, including the nuclear scattering
 length density ($\beta{_N}$) are
+length density $(\beta{_N})$ are
 .. math::
     \beta_{\pm\pm} =  \beta_N \mp D_M M_{\perp x'}
     \text{ when there are no spin-flips}
+    \text{ for non spin-flip states}
 and
 …
 .. math::
     \beta_{\pm\mp} =  -D_M (M_{\perp y'} \pm iM_{\perp z'})
     \text{ when there are}
+    \text{ for spin-flip states}
 where
 .. math::
     M_{\perp x'} = M_{0q_x}\cos(\theta_{up})+M_{0q_y}\sin(\theta_{up}) \\
     M_{\perp y'} = M_{0q_y}\cos(\theta_{up})-M_{0q_x}\sin(\theta_{up}) \\
     M_{\perp z'} = M_{0z} \\
     M_{0q_x} = (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\
     M_{0q_y} = (M_{0y}\sin\phi - M_{0x}\cos\phi)\sin\phi
+    M_{\perp x'} &= M_{0q_x}\cos(\theta_{up})+M_{0q_y}\sin(\theta_{up}) \\
+    M_{\perp y'} &= M_{0q_y}\cos(\theta_{up})-M_{0q_x}\sin(\theta_{up}) \\
+    M_{\perp z'} &= M_{0z} \\
+    M_{0q_x} &= (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\
+    M_{0q_y} &= (M_{0y}\sin\phi - M_{0x}\cos\phi)\sin\phi
 Here, $M_{0x}$, $M_{0x}$, $M_{0z}$ are the x, y and z components
 …
 .. math::
     M_{0x} = M_0\cos\theta_M\cos\phi_M \\
     M_{0y} = M_0\sin\theta_M \\
     M_{0z} = -M_0\cos\theta_M\sin\phi_M
+    M_{0x} &= M_0\cos\theta_M\cos\phi_M \\
+    M_{0y} &= M_0\sin\theta_M \\
+    M_{0z} &= -M_0\cos\theta_M\sin\phi_M
 and the magnetization angles $\theta_M$ and $\phi_M$ are defined in
 …
 ===========   ================================================================
  M0_sld        = $D_M M_0$
  Up_theta      = $\theta_\mathrm{up}$
  M_theta       = $\theta_M$
  M_phi         = $\phi_M$
  Up_frac_i     = (spin up)/(spin up + spin down) neutrons *before* the sample
  Up_frac_f     = (spin up)/(spin up + spin down) neutrons *after* the sample
+ M0:sld       $D_M M_0$
+ mtheta:sld   $\theta_M$
+ mphi:sld     $\phi_M$
+ up:angle     $\theta_\mathrm{up}$
+ up:frac_i    $u_i$ = (spin up)/(spin up + spin down) *before* the sample
+ up:frac_f    $u_f$ = (spin up)/(spin up + spin down) *after* the sample
 ===========   ================================================================
 .. note::
     The values of the 'Up_frac_i' and 'Up_frac_f' must be in the range 0 to 1.
+    The values of the 'up:frac_i' and 'up:frac_f' must be in the range 0 to 1.
 *Document History*
 | 2015-05-02 Steve King
 | 2017-05-08 Paul Kienzle
+| 2017-11-15 Paul Kienzle

doc/guide/pd/polydispersity.rst

-                      r1f058ea
+                      reda8b30
 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
+.. _polydispersityhelp:
 Polydispersity Distributions
 ----------------------------
 With some models in sasmodels we can calculate the average form factor for a
+With some models in sasmodels we can calculate the average intensity for a
 population of particles that exhibit size and/or orientational
 polydispersity. The resultant form factor is normalized by the average
+polydispersity. The resultant intensity is normalized by the average
 particle volume such that

doc/guide/resolution.rst

-                      r1f058ea
+                      r0db85af
 $x'_0 = x_0 \cos(\theta) + y_0 \sin(\theta)$ and
 $y'_0 = -x_0 \sin(\theta) + y_0 \cos(\theta)$.
 Note that the rotation angle is zero for a $x$\ -\ $y$ symmetric
+Note that the rotation angle is zero for a $x$-$y$ symmetric
 elliptical Gaussian distribution. The $A$ is a normalization factor.
 …
 Since the weighting factor on each of the bins is known, it is convenient to
 transform $x'$\ -\ $y'$ back to $x$\ -\ $y$ coordinates (by rotating it
+transform $x'$-$y'$ back to $x$-$y$ coordinates (by rotating it
 by $-\theta$ around the $z$\ -axis).
 …
     y'_0 &= 0
 while for a $x$\ -\ $y$ symmetric smear
+while for a $x$-$y$ symmetric smear
 .. math::

explore/jitter.py

Property mode changed from 100644 to 100755

-                      r85190c2
+                      rff10479
+#!/usr/bin/env python
 """
 Application to explore the difference between sasview 3.x orientation
 dispersity and possible replacement algorithms.
 """
+from __future__ import division, print_function
+import sys, os
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
+import argparse
 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)):
+    """
+    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1)
+    """
     #ax.plot([0,0],[0,0],[1,-1])
     #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8)
 …
     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, dist='gaussian', size=(0.1, 0.4, 1.0)):
+    """
+    Represent jitter as a set of shapes at different orientations.
+    """
+    # set max diagonal to 0.95
+    scale = 0.95/sqrt(sum(v**2 for v in size))
+    size = tuple(scale*v for v in size)
+    draw_shape = draw_parallelepiped
+    #draw_shape = draw_ellipsoid
     #np.random.seed(10)
 …
         [ 1,  1,  1],
+    ]
+    dtheta, dphi, dpsi = jitter
     if dtheta == 0:
         cloud = [v for v in cloud if v[0] == 0]
 …
     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)
+    scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist]
     for point in cloud:
         shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]]
         draw_shape(ax, size, view, shimmy, alpha=0.8)
+        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]]
+        draw_shape(ax, size, view, delta, alpha=0.8)
     for v in 'xyz':
         a, b, c = size
 …
         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):
+    """Draw an ellipsoid."""
     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)
 …
     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):
+    """Draw a parallelepiped."""
     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])
 …
     ])
+    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)
+    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_sphere(ax, radius=10., steps=100):
+    """Draw a sphere"""
     u = np.linspace(0, 2 * np.pi, steps)
     v = np.linspace(0, np.pi, steps)
 …
     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)
+    if dist == 'gauss':
+PROJECTIONS = [
+    # in order of PROJECTION number; do not change without updating the
+    # constants in kernel_iq.c
+    'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance',
+]
+def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian',
+              projection='equirectangular'):
+    """
+    Draw the dispersion mesh showing the theta-phi orientations at which
+    the model will be evaluated.
+    jitter projections
+    <https://en.wikipedia.org/wiki/List_of_map_projections>
+    equirectangular (standard latitude-longitude mesh)
+        <https://en.wikipedia.org/wiki/Equirectangular_projection>
+        Allows free movement in phi (around the equator), but theta is
+        limited to +/- 90, and points are cos-weighted. Jitter in phi is
+        uniform in weight along a line of latitude.  With small theta and
+        phi ranging over +/- 180 this forms a wobbling disk.  With small
+        phi and theta ranging over +/- 90 this forms a wedge like a slice
+        of an orange.
+    azimuthal_equidistance (Postel)
+        <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>
+        Preserves distance from center, and so is an excellent map for
+        representing a bivariate gaussian on the surface.  Theta and phi
+        operate identically, cutting wegdes from the antipode of the viewing
+        angle.  This unfortunately does not allow free movement in either
+        theta or phi since the orthogonal wobble decreases to 0 as the body
+        rotates through 180 degrees.
+    sinusoidal (Sanson-Flamsteed, Mercator equal-area)
+        <https://en.wikipedia.org/wiki/Sinusoidal_projection>
+        Preserves arc length with latitude, giving bad behaviour at
+        theta near +/- 90.  Theta and phi operate somewhat differently,
+        so a system with a-b-c dtheta-dphi-dpsi will not give the same
+        value as one with b-a-c dphi-dtheta-dpsi, as would be the case
+        for azimuthal equidistance.  Free movement using theta or phi
+        uniform over +/- 180 will work, but not as well as equirectangular
+        phi, with theta being slightly worse.  Computationally it is much
+        cheaper for wide theta-phi meshes since it excludes points which
+        lie outside the sinusoid near theta +/- 90 rather than packing
+        them close together as in equirectangle.  Note that the poles
+        will be slightly overweighted for theta > 90 with the circle
+        from theta at 90+dt winding backwards around the pole, overlapping
+        the circle from theta at 90-dt.
+    Guyou (hemisphere-in-a-square) **not weighted**
+        <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>
+        With tiling, allows rotation in phi or theta through +/- 180, with
+        uniform spacing.  Both theta and phi allow free rotation, with wobble
+        in the orthogonal direction reasonably well behaved (though not as
+        good as equirectangular phi). The forward/reverse transformations
+        relies on elliptic integrals that are somewhat expensive, so the
+        behaviour has to be very good to justify the cost and complexity.
+        The weighting function for each point has not yet been computed.
+        Note: run the module *guyou.py* directly and it will show the forward
+        and reverse mappings.
+    azimuthal_equal_area  **incomplete**
+        <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>
+        Preserves the relative density of the surface patches.  Not that
+        useful and not completely implemented
+    Gauss-Kreuger **not implemented**
+        <https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Ellipsoidal_transverse_Mercator>
+        Should allow free movement in theta, but phi is distorted.
+    """
+    theta, phi, psi = view
+    dtheta, dphi, dpsi = jitter
+    t = np.linspace(-1, 1, n)
+    weights = np.ones_like(t)
+    if dist == 'gaussian':
+        t *= 3
         weights = exp(-0.5*t**2)
+    elif dist == 'rect':
+        weights = np.ones_like(t)
+    elif dist == 'rectangle':
+        # Note: uses sasmodels ridiculous definition of rectangle width
+        t *= sqrt(3)
+    elif dist == 'uniform':
+        pass
     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)
+        raise ValueError("expected dist to be gaussian, rectangle or uniform")
+    if projection == 'equirectangular':  #define PROJECTION 1
+        def rotate(theta_i, phi_j):
+            return Rx(phi_j)*Ry(theta_i)
+        def weight(theta_i, phi_j, wi, wj):
+            return wi*wj*abs(cos(radians(theta_i)))
+    elif projection == 'sinusoidal':  #define PROJECTION 2
+        def rotate(theta_i, phi_j):
+            latitude = theta_i
+            scale = cos(radians(latitude))
+            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0
+            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
+            return Rx(longitude)*Ry(latitude)
+        def weight(theta_i, phi_j, wi, wj):
+            latitude = theta_i
+            scale = cos(radians(latitude))
+            w = 1 if abs(phi_j) < abs(scale)*180 else 0
+            return w*wi*wj
+    elif projection == 'guyou':  #define PROJECTION 3  (eventually?)
+        def rotate(theta_i, phi_j):
+            from guyou import guyou_invert
+            #latitude, longitude = guyou_invert([theta_i], [phi_j])
+            longitude, latitude = guyou_invert([phi_j], [theta_i])
+            return Rx(longitude[0])*Ry(latitude[0])
+        def weight(theta_i, phi_j, wi, wj):
+            return wi*wj
+    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry
+        def rotate(theta_i, phi_j):
+            latitude = sqrt(theta_i**2 + phi_j**2)
+            longitude = degrees(np.arctan2(phi_j, theta_i))
+            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
+            return Rz(longitude)*Ry(latitude)
+        def weight(theta_i, phi_j, wi, wj):
+            # Weighting for each point comes from the integral:
+            #     \int\int I(q, lat, log) sin(lat) dlat dlog
+            # We are doing a conformal mapping from disk to sphere, so we need
+            # a change of variables g(theta, phi) -> (lat, long):
+            #     lat, long = sqrt(theta^2 + phi^2), arctan(phi/theta)
+            # giving:
+            #     dtheta dphi = det(J) dlat dlong
+            # where J is the jacobian from the partials of g. Using
+            #     R = sqrt(theta^2 + phi^2),
+            # then
+            #     J = [[x/R, Y/R], -y/R^2, x/R^2]]
+            # and
+            #     det(J) = 1/R
+            # with the final integral being:
+            #    \int\int I(q, theta, phi) sin(R)/R dtheta dphi
+            #
+            # This does approximately the right thing, decreasing the weight
+            # of each point as you go farther out on the disk, but it hasn't
+            # yet been checked against the 1D integral results. Prior
+            # to declaring this "good enough" and checking that integrals
+            # work in practice, we will examine alternative mappings.
+            #
+            # The issue is that the mapping does not support the case of free
+            # rotation about a single axis correctly, with a small deviation
+            # in the orthogonal axis independent of the first axis.  Like the
+            # usual polar coordiates integration, the integrated sections
+            # form wedges, though at least in this case the wedge cuts through
+            # the entire sphere, and treats theta and phi identically.
+            latitude = sqrt(theta_i**2 + phi_j**2)
+            w = sin(radians(latitude))/latitude if latitude != 0 else 1
+            return w*wi*wj if latitude < 180 else 0
+    elif projection == 'azimuthal_equal_area':
+        def rotate(theta_i, phi_j):
+            R = min(1, sqrt(theta_i**2 + phi_j**2)/180)
+            latitude = 180-degrees(2*np.arccos(R))
+            longitude = degrees(np.arctan2(phi_j, theta_i))
+            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
+            return Rz(longitude)*Ry(latitude)
+        def weight(theta_i, phi_j, wi, wj):
+            latitude = sqrt(theta_i**2 + phi_j**2)
+            w = sin(radians(latitude))/latitude if latitude != 0 else 1
+            return w*wi*wj if latitude < 180 else 0
+    else:
+        raise ValueError("unknown projection %r"%projection)
+    # mesh in theta, phi formed by rotating z
+    z = np.matrix([[0], [0], [radius]])
+    points = np.hstack([rotate(theta_i, phi_j)*z
+                        for theta_i in dtheta*t
+                        for phi_j in dphi*t])
+    # select just the active points (i.e., those with phi < 180
+    w = np.array([weight(theta_i, phi_j, wi, wj)
+                  for wi, theta_i in zip(weights, dtheta*t)
+                  for wj, phi_j in zip(weights, dphi*t)])
+    #print(max(w), min(w), min(w[w>0]))
+    points = points[:, w>0]
+    w = w[w>0]
+    w /= max(w)
+    if 0: # Kent distribution
+        points = np.hstack([Rx(phi_j)*Ry(theta_i)*z for theta_i in 30*t for phi_j in 60*t])
+        xp, yp, zp = [np.array(v).flatten() for v in points]
+        kappa = max(1e6, radians(dtheta)/(2*pi))
+        beta = 1/max(1e-6, radians(dphi)/(2*pi))/kappa
+        w = exp(kappa*zp) #+ beta*(xp**2 + yp**2)
+        print(kappa, dtheta, radians(dtheta), min(w), max(w), sum(w))
+        #w /= abs(cos(radians(
+        #w /= sum(w)
+    # rotate relative to beam
+    points = orient_relative_to_beam(view, points)
+    x, y, z = [np.array(v).flatten() for v in points]
+    #plt.figure(2); plt.clf(); plt.hist(z, bins=np.linspace(-1, 1, 51))
     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
+def draw_labels(ax, view, jitter, text):
+    """
+    Draw text at a particular location.
+    """
+    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)
+    # TODO: zdir for labels is broken, and labels aren't appearing.
+    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)
+# Definition of rotation matrices comes from wikipedia:
+#    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
 def Rx(angle):
+    """Construct a matrix to rotate points about *x* by *angle* degrees."""
     a = radians(angle)
     R = [[1., 0., 0.],
          [0.,  cos(a), sin(a)],
          [0., -sin(a), cos(a)]]
+    R = [[1, 0, 0],
+         [0, +cos(a), -sin(a)],
+         [0, +sin(a), +cos(a)]]
     return np.matrix(R)
 def Ry(angle):
+    """Construct a matrix to rotate points about *y* by *angle* degrees."""
     a = radians(angle)
     R = [[cos(a), 0., -sin(a)],
          [0., 1., 0.],
          [sin(a), 0.,  cos(a)]]
+    R = [[+cos(a), 0, +sin(a)],
+         [0, 1, 0],
+         [-sin(a), 0, +cos(a)]]
     return np.matrix(R)
 def Rz(angle):
+    """Construct a matrix to rotate points about *z* by *angle* degrees."""
     a = radians(angle)
     R = [[cos(a), -sin(a), 0.],
          [sin(a),  cos(a), 0.],
          [0., 0., 1.]]
+    R = [[+cos(a), -sin(a), 0],
+         [+sin(a), +cos(a), 0],
+         [0, 0, 1]]
     return np.matrix(R)
+def main():
+def transform_xyz(view, jitter, x, y, z):
+    """
+    Send a set of (x,y,z) points through the jitter and view transforms.
+    """
+    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):
+    """
+    Apply the jitter transform to a set of points.
+    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
+    """
+    dtheta, dphi, dpsi = jitter
+    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
+    return points
+def orient_relative_to_beam(view, points):
+    """
+    Apply the view transform to a set of points.
+    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
+    """
+    theta, phi, psi = view
+    points = Rz(phi)*Ry(theta)*Rz(psi)*points
+    return points
+# translate between number of dimension of dispersity and the number of
+# points along each dimension.
+PD_N_TABLE = {
+    (0, 0, 0): (0, 0, 0),     # 0
+    (1, 0, 0): (100, 0, 0),   # 100
+    (0, 1, 0): (0, 100, 0),
+    (0, 0, 1): (0, 0, 100),
+    (1, 1, 0): (30, 30, 0),   # 900
+    (1, 0, 1): (30, 0, 30),
+    (0, 1, 1): (0, 30, 30),
+    (1, 1, 1): (15, 15, 15),  # 3375
+}
+def clipped_range(data, portion=1.0, mode='central'):
+    """
+    Determine range from data.
+    If *portion* is 1, use full range, otherwise use the center of the range
+    or the top of the range, depending on whether *mode* is 'central' or 'top'.
+    """
+    if portion == 1.0:
+        return data.min(), data.max()
+    elif mode == 'central':
+        data = np.sort(data.flatten())
+        offset = int(portion*len(data)/2 + 0.5)
+        return data[offset], data[-offset]
+    elif mode == 'top':
+        data = np.sort(data.flatten())
+        offset = int(portion*len(data) + 0.5)
+        return data[offset], data[-1]
+def draw_scattering(calculator, ax, view, jitter, dist='gaussian'):
+    """
+    Plot the scattering for the particular view.
+    *calculator* is returned from :func:`build_model`.  *ax* are the 3D axes
+    on which the data will be plotted.  *view* and *jitter* are the current
+    orientation and orientation dispersity.  *dist* is one of the sasmodels
+    weight distributions.
+    """
+    if dist == 'uniform':  # uniform is not yet in this branch
+        dist, scale = 'rectangle', 1/sqrt(3)
+    else:
+        scale = 1
+    # add the orientation parameters to the model parameters
+    theta, phi, psi = view
+    theta_pd, phi_pd, psi_pd = [scale*v for v in jitter]
+    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd>0, phi_pd>0, psi_pd>0)]
+    ## increase pd_n for testing jitter integration rather than simple viz
+    #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)]
+    pars = dict(
+        theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n,
+        phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n,
+        psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n,
+    )
+    pars.update(calculator.pars)
+    # compute the pattern
+    qx, qy = calculator._data.x_bins, calculator._data.y_bins
+    Iqxy = calculator(**pars).reshape(len(qx), len(qy))
+    # scale it and draw it
+    Iqxy = np.log(Iqxy)
+    if calculator.limits:
+        # use limits from orientation (0,0,0)
+        vmin, vmax = calculator.limits
+    else:
+        vmax = Iqxy.max()
+        vmin = vmax*10**-7
+        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top')
+    #print("range",(vmin,vmax))
+    #qx, qy = np.meshgrid(qx, qy)
+    if 0:
+        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i')
+        level[level<0] = 0
+        colors = plt.get_cmap()(level)
+        ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors)
+    elif 1:
+        ax.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1,
+                    levels=np.linspace(vmin, vmax, 24))
+    else:
+        ax.pcolormesh(qx, qy, Iqxy)
+def build_model(model_name, n=150, qmax=0.5, **pars):
+    """
+    Build a calculator for the given shape.
+    *model_name* is any sasmodels model.  *n* and *qmax* define an n x n mesh
+    on which to evaluate the model.  The remaining parameters are stored in
+    the returned calculator as *calculator.pars*.  They are used by
+    :func:`draw_scattering` to set the non-orientation parameters in the
+    calculation.
+    Returns a *calculator* function which takes a dictionary or parameters and
+    produces Iqxy.  The Iqxy value needs to be reshaped to an n x n matrix
+    for plotting.  See the :class:`sasmodels.direct_model.DirectModel` class
+    for details.
+    """
+    from sasmodels.core import load_model_info, build_model
+    from sasmodels.data import empty_data2D
+    from sasmodels.direct_model import DirectModel
+    model_info = load_model_info(model_name)
+    model = build_model(model_info) #, dtype='double!')
+    q = np.linspace(-qmax, qmax, n)
+    data = empty_data2D(q, q)
+    calculator = DirectModel(data, model)
+    # stuff the values for non-orientation parameters into the calculator
+    calculator.pars = pars.copy()
+    calculator.pars.setdefault('backgound', 1e-3)
+    # fix the data limits so that we can see if the pattern fades
+    # under rotation or angular dispersion
+    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars)
+    Iqxy = np.log(Iqxy)
+    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top')
+    calculator.limits = vmin, vmax+1
+    return calculator
+def select_calculator(model_name, n=150, size=(10,40,100)):
+    """
+    Create a model calculator for the given shape.
+    *model_name* is one of sphere, cylinder, ellipsoid, triaxial_ellipsoid,
+    parallelepiped or bcc_paracrystal. *n* is the number of points to use
+    in the q range.  *qmax* is chosen based on model parameters for the
+    given model to show something intersting.
+    Returns *calculator* and tuple *size* (a,b,c) giving minor and major
+    equitorial axes and polar axis respectively.  See :func:`build_model`
+    for details on the returned calculator.
+    """
+    a, b, c = size
+    if model_name == 'sphere':
+        calculator = build_model('sphere', n=n, radius=c)
+        a = b = c
+    elif model_name == 'bcc_paracrystal':
+        calculator = build_model('bcc_paracrystal', n=n, dnn=c,
+                                  d_factor=0.06, radius=40)
+        a = b = c
+    elif model_name == 'cylinder':
+        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c)
+        a = b
+    elif model_name == 'ellipsoid':
+        calculator = build_model('ellipsoid', n=n, qmax=1.0,
+                                 radius_polar=c, radius_equatorial=b)
+        a = b
+    elif model_name == 'triaxial_ellipsoid':
+        calculator = build_model('triaxial_ellipsoid', n=n, qmax=0.5,
+                                 radius_equat_minor=a,
+                                 radius_equat_major=b,
+                                 radius_polar=c)
+    elif model_name == 'parallelepiped':
+        calculator = build_model('parallelepiped', n=n, a=a, b=b, c=c)
+    else:
+        raise ValueError("unknown model %s"%model_name)
+    return calculator, (a, b, c)
+SHAPES = [
+    'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal',
+    'cylinder', 'ellipsoid',
+    'sphere',
+ ]
+DISTRIBUTIONS = [
+    'gaussian', 'rectangle', 'uniform',
+]
+DIST_LIMITS = {
+    'gaussian': 30,
+    'rectangle': 90/sqrt(3),
+    'uniform': 90,
+}
+def run(model_name='parallelepiped', size=(10, 40, 100),
+        dist='gaussian', mesh=30,
+        projection='equirectangular'):
+    """
+    Show an interactive orientation and jitter demo.
+    *model_name* is one of the models available in :func:`select_model`.
+    """
+    # projection number according to 1-order position in list, but
+    # only 1 and 2 are implemented so far.
+    from sasmodels import generate
+    generate.PROJECTION = PROJECTIONS.index(projection) + 1
+    if generate.PROJECTION > 2:
+        print("*** PROJECTION %s not implemented in scattering function ***"%projection)
+        generate.PROJECTION = 2
+    # set up calculator
+    calculator, size = select_calculator(model_name, n=150, size=size)
+    ## uncomment to set an independent the colour range for every view
+    ## If left commented, the colour range is fixed for all views
+    calculator.limits = None
+    ## initial view
+    #theta, dtheta = 70., 10.
+    #phi, dphi = -45., 3.
+    #psi, dpsi = -45., 3.
+    theta, phi, psi = 0, 0, 0
+    dtheta, dphi, dpsi = 0, 0, 0
+    ## create the plot window
     #plt.hold(True)
+    plt.subplots(num=None, figsize=(5.5, 5.5))
     plt.set_cmap('gist_earth')
     plt.clf()
+    plt.gcf().canvas.set_window_title(projection)
     #gs = gridspec.GridSpec(2,1,height_ratios=[4,1])
     #ax = plt.subplot(gs[0], projection='3d')
     ax = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d')
+    theta, dtheta = 70., 10.
+    phi, dphi = -45., 3.
+    psi, dpsi = -45., 3.
+    theta, phi, psi = 0, 0, 0
+    dtheta, dphi, dpsi = 0, 0, 0
+    #dist = 'rect'
+    dist = 'gauss'
+    try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection
+        ax.axis('square')
+    except Exception:
+        pass
     axcolor = 'lightgoldenrodyellow'
+    ## add control widgets to plot
     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)
 …
     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)
     axdpsi= plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor)
+    sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta)
+    sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi)
+    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi)
+    # Note: using ridiculous definition of rectangle distribution, whose width
+    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep
+    # the maximum width to 90.
+    dlimit = DIST_LIMITS[dist]
+    sdtheta = Slider(axdtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta)
+    sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi)
+    sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi)
+    ## callback to draw the new view
     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
+        # set small jitter as 0 if multiple pd dims
+        dims = sum(v > 0 for v in jitter)
+        limit = [0, 0, 0.5, 5][dims]
+        jitter = [0 if v < limit else v for v in jitter]
         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, dist=dist, size=size)
+        #draw_jitter(ax, view, (0,0,0))
+        draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection)
+        draw_scattering(calculator, ax, view, jitter, dist=dist)
         plt.gcf().canvas.draw()
+    ## bind control widgets to view updater
     stheta.on_changed(lambda v: update(v,'theta'))
     sphi.on_changed(lambda v: update(v, 'phi'))
 …
     sdpsi.on_changed(lambda v: update(v, 'dpsi'))
+    ## initialize view
     update(None, 'phi')
+    ## go interactive
     plt.show()
+def main():
+    parser = argparse.ArgumentParser(
+        description="Display jitter",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+        )
+    parser.add_argument('-p', '--projection', choices=PROJECTIONS, default=PROJECTIONS[0], help='coordinate projection')
+    parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths')
+    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='jitter distribution')
+    parser.add_argument('-m', '--mesh', type=int, default=30, help='#points in theta-phi mesh')
+    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape')
+    opts = parser.parse_args()
+    size = tuple(int(v) for v in opts.size.split(','))
+    run(opts.shape, size=size,
+        mesh=opts.mesh, dist=opts.distribution,
+        projection=opts.projection)
 if __name__ == "__main__":

explore/precision.py

-                      r237c9cf
+                      ra1c5758
 #!/usr/bin/env python
 r"""
+Show numerical precision of $2 J_1(x)/x$.
+Show numerical precision of various expressions.
+Evaluates the same function(s) in single and double precision and compares
+the results to 500 digit mpmath evaluation of the same function.
+Note: a quick way to generation C and python code for taylor series
+expansions from sympy:
+    import sympy as sp
+    x = sp.var("x")
+    f = sp.sin(x)/x
+    t = sp.series(f, n=12).removeO()  # taylor series with no O(x^n) term
+    p = sp.horner(t)   # Horner representation
+    p = p.replace(x**2, sp.var("xsq")  # simplify if alternate terms are zero
+    p = p.n(15)  # evaluate coefficients to 15 digits (optional)
+    c_code = sp.ccode(p, assign_to=sp.var("p"))  # convert to c code
+    py_code = c[:-1]  # strip semicolon to convert c to python
+    # mpmath has pade() rational function approximation, which might work
+    # better than the taylor series for some functions:
+    P, Q = mp.pade(sp.Poly(t.n(15),x).coeffs(), L, M)
+    P = sum(a*x**n for n,a in enumerate(reversed(P)))
+    Q = sum(a*x**n for n,a in enumerate(reversed(Q)))
+    c_code = sp.ccode(sp.horner(P)/sp.horner(Q), assign_to=sp.var("p"))
+    # There are richardson and shanks series accelerators in both sympy
+    # and mpmath that may be helpful.
 """
 from __future__ import division, print_function
 …
     np_function=scipy.special.erfc,
     ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]),
+    limits=(-5., 5.),
+)
+add_function(
+    name="expm1(x)",
+    mp_function=mp.expm1,
+    np_function=np.expm1,
+    ocl_function=make_ocl("return expm1(q);", "sas_expm1"),
     limits=(-5., 5.),
+)
 …
+)
+replacement_expm1 = """\
+      double x = (double)q;  // go back to float for single precision kernels
+      // Adapted from the cephes math library.
+      // Copyright 1984 - 1992 by Stephen L. Moshier
+      if (x != x || x == 0.0) {
+         return x; // NaN and +/- 0
+      } else if (x < -0.5 || x > 0.5) {
+         return exp(x) - 1.0;
+      } else {
+         const double xsq = x*x;
+         const double p = (((
+            +1.2617719307481059087798E-4)*xsq
+            +3.0299440770744196129956E-2)*xsq
+            +9.9999999999999999991025E-1);
+         const double q = ((((
+            +3.0019850513866445504159E-6)*xsq
+            +2.5244834034968410419224E-3)*xsq
+            +2.2726554820815502876593E-1)*xsq
+            +2.0000000000000000000897E0);
+         double r = x * p;
+         r =  r / (q - r);
+         return r+r;
+       }
+"""
+add_function(
+    name="sas_expm1(x)",
+    mp_function=mp.expm1,
+    np_function=np.expm1,
+    ocl_function=make_ocl(replacement_expm1, "sas_expm1"),
+)
 # Alternate versions of 3 j1(x)/x, for posterity
 def taylor_3j1x_x(x):

sasmodels/compare.py

-                      rced5bd2
+                      r3bfd924
 from . import exception
 from .data import plot_theory, empty_data1D, empty_data2D, load_data
 from .direct_model import DirectModel
+from .direct_model import DirectModel, get_mesh
 from .convert import revert_name, revert_pars, constrain_new_to_old
 from .generate import FLOAT_RE
+from .weights import plot_weights
 try:
 …
     -pars/-nopars* prints the parameter set or not
     -default/-demo* use demo vs default parameters
+    -sphere[=150] set up spherical integration over theta/phi using n points
     === calculation options ===
     -mono*/-poly force monodisperse or allow polydisperse demo parameters
+    -mono*/-poly force monodisperse or allow polydisperse random parameters
     -cutoff=1e-5* cutoff value for including a point in polydispersity
     -magnetic/-nonmagnetic* suppress magnetism
 …
     === precision options ===
     -calc=default uses the default calcution precision
+    -engine=default uses the default calcution precision
     -single/-double/-half/-fast sets an OpenCL calculation engine
     -single!/-double!/-quad! sets an OpenMP calculation engine
 …
     -abs/-rel* plot relative or absolute error
     -title="note" adds note to the plot title, after the model name
+    -weights shows weights plots for the polydisperse parameters
     === output options ===
 …
 vary from 64-bit to 128-bit, with 80-bit floats being common (1e-19 precision).
 On unix and mac you may need single quotes around the DLL computation
 engines, such as -calc='single!,double!' since !, is treated as a history
+engines, such as -engine='single!,double!' since !, is treated as a history
 expansion request in the shell.
 …
     # compare single and double precision calculation for a barbell
     sascomp barbell -calc=single,double
+    sascomp barbell -engine=single,double
     # generate 10 random lorentz models, with seed=27
 …
     # model timing test requires multiple evals to perform the estimate
     sascomp pringle -calc=single,double -timing=100,100 -noplot
+    sascomp pringle -engine=single,double -timing=100,100 -noplot
 """
 …
         # orientation in [-180,180], orientation pd in [0,45]
         if p.endswith('_pd'):
             return 0., 45.
+            return 0., 180.
         else:
             return -180., 180.
 …
+def limit_dimensions(model_info, pars, maxdim):
+    # type: (ModelInfo, ParameterSet, float) -> None
+    """
+    Limit parameters of units of Ang to maxdim.
+    """
+    for p in model_info.parameters.call_parameters:
+        value = pars[p.name]
+        if p.units == 'Ang' and value > maxdim:
+            pars[p.name] = maxdim*10**np.random.uniform(-3,0)
 def constrain_pars(model_info, pars):
     # type: (ModelInfo, ParameterSet) -> None
 …
     Format the parameter list for printing.
     """
+    is2d = True
     lines = []
     parameters = model_info.parameters
 …
 def compare(opts, limits=None):
+def compare(opts, limits=None, maxdim=np.inf):
     # type: (Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
     """
 …
     as the values are adjusted, making it easier to see the effects of the
     parameters.
+    """
+    limits = np.Inf, -np.Inf
+    *maxdim* is the maximum value for any parameter with units of Angstrom.
+    """
     for k in range(opts['sets']):
+        opts['pars'] = parse_pars(opts)
+        if k > 1:
+            # print a separate seed for each dataset for better reproducibility
+            new_seed = np.random.randint(1000000)
+            print("Set %d uses -random=%i"%(k+1,new_seed))
+            np.random.seed(new_seed)
+        opts['pars'] = parse_pars(opts, maxdim=maxdim)
         if opts['pars'] is None:
             return
 …
         if opts['plot']:
             limits = plot_models(opts, result, limits=limits, setnum=k)
+        if opts['show_weights']:
+            base, _ = opts['engines']
+            base_pars, _ = opts['pars']
+            model_info = base._kernel.info
+            dim = base._kernel.dim
+            plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim))
     if opts['plot']:
         import matplotlib.pyplot as plt
         plt.show()
+    return limits
 def run_models(opts, verbose=False):
 …
 def plot_models(opts, result, limits=(np.Inf, -np.Inf), setnum=0):
+def plot_models(opts, result, limits=None, setnum=0):
     # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
+    import matplotlib.pyplot as plt
     base_value, comp_value = result['base_value'], result['comp_value']
     base_time, comp_time = result['base_time'], result['comp_time']
 …
     # Plot if requested
     view = opts['view']
     import matplotlib.pyplot as plt
     vmin, vmax = limits
     if have_base:
         vmin = min(vmin, base_value.min())
         vmax = max(vmax, base_value.max())
     if have_comp:
         vmin = min(vmin, comp_value.min())
         vmax = max(vmax, comp_value.max())
     limits = vmin, vmax
+    if limits is None:
+        vmin, vmax = np.inf, -np.inf
+        if have_base:
+            vmin = min(vmin, base_value.min())
+            vmax = max(vmax, base_value.max())
+        if have_comp:
+            vmin = min(vmin, comp_value.min())
+            vmax = max(vmax, comp_value.max())
+        limits = vmin, vmax
     if have_base:
 …
             err[err > cutoff] = cutoff
         #err,errstr = base/comp,"ratio"
+        plot_theory(data, None, resid=err, view=view, use_data=use_data)
+        plt.yscale(errview)
+        plot_theory(data, None, resid=err, view=errview, use_data=use_data)
+        plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear')
+        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best')
         plt.title("max %s = %.3g"%(errstr, abs(err).max()))
         #cbar_title = errstr if errview=="linear" else "log "+errstr
 …
 OPTIONS = [
     # Plotting
     'plot', 'noplot',
+    'plot', 'noplot', 'weights',
     'linear', 'log', 'q4',
     'rel', 'abs',
 …
     'demo', 'default',  # TODO: remove demo/default
     'nopars', 'pars',
+    'sphere', 'sphere=', # integrate over a sphere in 2d with n points
     # Calculation options
 …
     # Precision options
     'calc=',
+    'engine=',
     'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!',
     'sasview',  # TODO: remove sasview 3.x support
 …
         'sets'      : 0,
         'engine'    : 'default',
+        'evals'     : '1',
+        'count'     : '1',
+        'show_weights' : False,
+        'sphere'    : 0,
+    }
     for arg in flags:
 …
         elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
         elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:]
-        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:])
         elif arg.startswith('-title='):    opts['title'] = arg[7:]
         elif arg.startswith('-data='):     opts['datafile'] = arg[6:]
+        elif arg.startswith('-calc='):     opts['engine'] = arg[6:]
+        elif arg.startswith('-neval='):    opts['evals'] = arg[7:]
+        elif arg == '-random':  opts['seed'] = np.random.randint(1000000)
+        elif arg.startswith('-engine='):   opts['engine'] = arg[8:]
+        elif arg.startswith('-neval='):    opts['count'] = arg[7:]
+        elif arg.startswith('-random='):
+            opts['seed'] = int(arg[8:])
+            opts['sets'] = 0
+        elif arg == '-random':
+            opts['seed'] = np.random.randint(1000000)
+            opts['sets'] = 0
+        elif arg.startswith('-sphere'):
+            opts['sphere'] = int(arg[8:]) if len(arg) > 7 else 150
+            opts['is2d'] = True
         elif arg == '-preset':  opts['seed'] = -1
         elif arg == '-mono':    opts['mono'] = True
 …
         elif arg == '-demo':    opts['use_demo'] = True
         elif arg == '-default': opts['use_demo'] = False
+        elif arg == '-weights': opts['show_weights'] = True
         elif arg == '-html':    opts['html'] = True
         elif arg == '-help':    opts['html'] = True
 …
     if PAR_SPLIT in opts['engine']:
         engine_types = opts['engine'].split(PAR_SPLIT, 2)
+        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2)
         comparison = True
     else:
         engine_types = [opts['engine']]*2
     if PAR_SPLIT in opts['evals']:
         evals = [int(k) for k in opts['evals'].split(PAR_SPLIT, 2)]
+        opts['engine'] = [opts['engine']]*2
+    if PAR_SPLIT in opts['count']:
+        opts['count'] = [int(k) for k in opts['count'].split(PAR_SPLIT, 2)]
         comparison = True
     else:
         evals = [int(opts['evals'])]*2
+        opts['count'] = [int(opts['count'])]*2
     if PAR_SPLIT in opts['cutoff']:
         cutoff = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)]
+        opts['cutoff'] = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)]
         comparison = True
     else:
         cutoff = [float(opts['cutoff'])]*2
     base = make_engine(model_info[0], data, engine_types[0], cutoff[0])
+        opts['cutoff'] = [float(opts['cutoff'])]*2
+    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0])
     if comparison:
         comp = make_engine(model_info[1], data, engine_types[1], cutoff[1])
+        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1])
     else:
         comp = None
 …
         'data'      : data,
         'name'      : names,
+        'def'       : model_info,
+        'count'     : evals,
+        'info'      : model_info,
         'engines'   : [base, comp],
         'values'    : values,
 …
     # pylint: enable=bad-whitespace
+    # Set the integration parameters to the half sphere
+    if opts['sphere'] > 0:
+        set_spherical_integration_parameters(opts, opts['sphere'])
     return opts
+def parse_pars(opts):
+    model_info, model_info2 = opts['def']
+def set_spherical_integration_parameters(opts, steps):
+    """
+    Set integration parameters for spherical integration over the entire
+    surface in theta-phi coordinates.
+    """
+    # Set the integration parameters to the half sphere
+    opts['values'].extend([
+        #'theta=90',
+        'theta_pd=%g'%(90/np.sqrt(3)),
+        'theta_pd_n=%d'%steps,
+        'theta_pd_type=rectangle',
+        #'phi=0',
+        'phi_pd=%g'%(180/np.sqrt(3)),
+        'phi_pd_n=%d'%(2*steps),
+        'phi_pd_type=rectangle',
+        #'background=0',
+    ])
+    if 'psi' in opts['info'][0].parameters:
+        opts['values'].extend([
+            #'psi=0',
+            'psi_pd=%g'%(180/np.sqrt(3)),
+            'psi_pd_n=%d'%(2*steps),
+            'psi_pd_type=rectangle',
+        ])
+        pass
+def parse_pars(opts, maxdim=np.inf):
+    model_info, model_info2 = opts['info']
     # Get demo parameters from model definition, or use default parameters
 …
     if opts['seed'] > -1:
         pars = randomize_pars(model_info, pars)
+        limit_dimensions(model_info, pars, maxdim)
         if model_info != model_info2:
             pars2 = randomize_pars(model_info2, pars2)
+            limit_dimensions(model_info, pars2, maxdim)
             # Share values for parameters with the same name
             for k, v in pars.items():
 …
     from . import rst2html
     info = opts['def'][0]
+    info = opts['info'][0]
     html = make_html(info)
     path = os.path.dirname(info.filename)
 …
         opts['pars'] = list(opts['pars'])
         p1, p2 = opts['pars']
         m1, m2 = opts['def']
+        m1, m2 = opts['info']
         self.fix_p2 = m1 != m2 or p1 != p2
         model_info = m1
 …
         self.starting_values = dict((k, v.value) for k, v in pars.items())
         self.pd_types = pd_types
         self.limits = np.Inf, -np.Inf
+        self.limits = None
     def revert_values(self):

sasmodels/core.py

-                      ra85a569
+                      r9e771a3
     Possible types include 'half', 'single', 'double' and 'quad'.  If the
     type is 'fast', then this is equivalent to dtype 'single' but using
+    fast native functions rather than those with the precision level guaranteed
+    by the OpenCL standard.
+    fast native functions rather than those with the precision level
+    guaranteed by the OpenCL standard.  'default' will choose the appropriate
+    default for the model and platform.
     Platform preference can be specfied ("ocl" vs "dll"), with the default

sasmodels/data.py

-                      rba7302a
+                      ra1c5758
     if hasattr(data, 'isSesans') and data.isSesans:
         _plot_result_sesans(data, None, None, use_data=True, limits=limits)
     elif hasattr(data, 'qx_data'):
+    elif hasattr(data, 'qx_data') and not getattr(data, 'radial', False):
         _plot_result2D(data, None, None, view, use_data=True, limits=limits)
     else:
 …
     if hasattr(data, 'isSesans') and data.isSesans:
         _plot_result_sesans(data, theory, resid, use_data=True, limits=limits)
     elif hasattr(data, 'qx_data'):
+    elif hasattr(data, 'qx_data') and not getattr(data, 'radial', False):
         _plot_result2D(data, theory, resid, view, use_data, limits=limits)
     else:
 …
     import matplotlib.pyplot as plt  # type: ignore
     from numpy.ma import masked_array, masked  # type: ignore
+    if getattr(data, 'radial', False):
+        radial_data.x = radial_data.q_data
+        radial_data.y = radial_data.data
     use_data = use_data and data.y is not None
 …
     ymin, ymax = min(data.qy_data), max(data.qy_data)
     if view == 'log':
+        vmin, vmax = np.log10(vmin), np.log10(vmax)
+        vmin_scaled, vmax_scaled= np.log10(vmin), np.log10(vmax)
+    else:
+        vmin_scaled, vmax_scaled = vmin, vmax
     plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)),
                interpolation='nearest', aspect=1, origin='lower',
+               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
+               extent=[xmin, xmax, ymin, ymax],
+               vmin=vmin_scaled, vmax=vmax_scaled)
     plt.xlabel("$q_x$/A$^{-1}$")
     plt.ylabel("$q_y$/A$^{-1}$")

sasmodels/details.py

-                      rccd5f01
+                      rce99754
 import numpy as np  # type: ignore
 from numpy import pi, cos, sin
+from numpy import pi, cos, sin, radians
 try:
 …
 try:
     from typing import List
+    from typing import List, Tuple, Sequence
 except ImportError:
     pass
 else:
     from .modelinfo import ModelInfo
+    from .modelinfo import ParameterTable
 …
     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
 …
     #print("off:",offset)
     # Check that we arn't using too many polydispersity loops
+    # Check that we aren't using too many polydispersity loops
     num_active = np.sum(length > 1)
     max_pd = model_info.parameters.max_pd
 …
 ZEROS = tuple([0.]*31)
 def make_kernel_args(kernel, pairs):
+def make_kernel_args(kernel, mesh):
     # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool]
     """
     Converts (value, weight) pairs into parameters for the kernel call.
+    Converts (value, dispersity, weight) for each parameter into kernel pars.
     Returns a CallDetails object indicating the polydispersity, a data object
 …
     npars = kernel.info.parameters.npars
     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 ((),())
+    scalars = [value for value, dispersity, weight in mesh]
+    # skipping scale and background when building values and weights
+    values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ())
+    #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights)
     length = np.array([len(w) for w in weights])
     offset = np.cumsum(np.hstack((0, length)))
     call_details = make_details(kernel.info, length, offset[:-1], offset[-1])
     # Pad value array to a 32 value boundaryd
     data_len = nvalues + 2*sum(len(v) for v in values)
+    data_len = nvalues + 2*sum(len(v) for v in dispersity)
     extra = (32 - data_len%32)%32
     data = np.hstack((scalars,) + values + weights + ZEROS[:extra])
+    data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra])
     data = data.astype(kernel.dtype)
     is_magnetic = convert_magnetism(kernel.info.parameters, data)
 …
     return call_details, data, is_magnetic
+def correct_theta_weights(parameters, dispersity, 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.  This is probably not a problem in practice (if dispersity
+    is +/- 90, then you probably should be using a 1-D model of the circular
+    average).
+    Note: scale and background parameters are not include in the tuples for
+    dispersity and weights, so index is parameters.theta_offset, not
+    parameters.theta_offset+2
+    Returns updated weights vectors
+    """
+    # TODO: explain in a comment why scale and background are missing
+    # Apparently the parameters.theta_offset similarly skips scale and
+    # and background, so the indexing works out, but they are still shipped
+    # to the kernel, so we need to add two there.
+    if parameters.theta_offset >= 0:
+        index = parameters.theta_offset
+        theta = dispersity[index]
+        # TODO: modify the dispersity vector to avoid the theta=-90,90,270,...
+        theta_weight = abs(cos(radians(theta)))
+        weights = tuple(theta_weight*w if k == index else w
+                        for k, w in enumerate(weights))
+    return weights
 def convert_magnetism(parameters, values):
+    # type: (ParameterTable, Sequence[np.ndarray]) -> bool
     """
     Convert magnetism values from polar to rectangular coordinates.
 …
     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
 …
 def dispersion_mesh(model_info, pars):
+def dispersion_mesh(model_info, mesh):
     # type: (ModelInfo) -> Tuple[List[np.ndarray], List[np.ndarray]]
     """
     Create a mesh grid of dispersion parameters and weights.
+    *mesh* is a list of (value, dispersity, weights), where the values
+    are the individual parameter values, and (dispersity, weights) is
+    the distribution of parameter values.
+    Only the volume parameters should be included in this list.  Orientation
+    parameters do not affect the calculation of effective radius or volume
+    ratio.  This is convenient since it avoids the distinction between
+    value and dispersity that is present in orientation parameters but not
+    shape parameters.
     Returns [p1,p2,...],w where pj is a vector of values for parameter j
 …
     parameter set in the vector.
     """
     value, weight = zip(*pars)
+    _, dispersity, weight = zip(*mesh)
     #weight = [w if len(w)>0 else [1.] for w in weight]
     weight = np.vstack([v.flatten() for v in meshgrid(*weight)])
     weight = np.prod(weight, axis=0)
     value = [v.flatten() for v in meshgrid(*value)]
+    dispersity = [v.flatten() for v in meshgrid(*dispersity)]
     lengths = [par.length for par in model_info.parameters.kernel_parameters
                if par.type == 'volume']
 …
         offset = 0
         for n in lengths:
             pars.append(np.vstack(value[offset:offset+n])
                         if n > 1 else value[offset])
+            pars.append(np.vstack(dispersity[offset:offset+n])
+                        if n > 1 else dispersity[offset])
             offset += n
         value = pars
     return value, weight
+        dispersity = pars
+    return dispersity, weight

sasmodels/direct_model.py

-                      rd1ff3a5
+                      r9e771a3
     *mono* is True if polydispersity should be set to none on all parameters.
     """
+    parameters = calculator.info.parameters
+    if mono:
+        active = lambda name: False
+    elif calculator.dim == '1d':
+        active = lambda name: name in parameters.pd_1d
+    elif calculator.dim == '2d':
+        active = lambda name: name in parameters.pd_2d
+    else:
+        active = lambda name: True
+    #print("pars",[p.id for p in parameters.call_parameters])
+    vw_pairs = [(get_weights(p, pars) if active(p.name)
+                 else ([pars.get(p.name, p.default)], [1.0]))
+                for p in parameters.call_parameters]
+    call_details, values, is_magnetic = make_kernel_args(calculator, vw_pairs)
+    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono)
+    #print("pars", list(zip(*mesh))[0])
+    call_details, values, is_magnetic = make_kernel_args(calculator, mesh)
     #print("values:", values)
     return calculator(call_details, values, cutoff, is_magnetic)
 def call_ER(model_info, pars):
 …
     return x, y, model_info.profile_axes
+def get_weights(parameter, values):
+    # type: (Parameter, Dict[str, float]) -> Tuple[np.ndarray, np.ndarray]
+def get_mesh(model_info, values, dim='1d', mono=False):
+    # type: (ModelInfo, Dict[str, float], str, bool) -> List[Tuple[float, np.ndarray, np.ndarry]]
+    """
+    Retrieve the dispersity mesh described by the parameter set.
+    Returns a list of *(value, dispersity, weights)* with one tuple for each
+    parameter in the model call parameters.  Inactive parameters return the
+    default value with a weight of 1.0.
+    """
+    parameters = model_info.parameters
+    if mono:
+        active = lambda name: False
+    elif dim == '1d':
+        active = lambda name: name in parameters.pd_1d
+    elif dim == '2d':
+        active = lambda name: name in parameters.pd_2d
+    else:
+        active = lambda name: True
+    #print("pars",[p.id for p in parameters.call_parameters])
+    mesh = [_get_par_weights(p, values, active(p.name))
+            for p in parameters.call_parameters]
+    return mesh
+def _get_par_weights(parameter, values, active=True):
+    # type: (Parameter, Dict[str, float]) -> Tuple[float, np.ndarray, np.ndarray]
     """
     Generate the distribution for parameter *name* given the parameter values
 …
     """
     value = float(values.get(parameter.name, parameter.default))
-    relative = parameter.relative_pd
-    limits = parameter.limits
-    disperser = values.get(parameter.name+'_pd_type', 'gaussian')
     npts = values.get(parameter.name+'_pd_n', 0)
     width = values.get(parameter.name+'_pd', 0.0)
+    nsigma = values.get(parameter.name+'_pd_nsigma', 3.0)
+    if npts == 0 or width == 0:
+        return [value], [1.0]
+    value, weight = weights.get_weights(
+        disperser, npts, width, nsigma, value, limits, relative)
+    return value, weight / np.sum(weight)
+def _vol_pars(model_info, pars):
+    relative = parameter.relative_pd
+    if npts == 0 or width == 0.0 or not active:
+        # Note: orientation parameters have the viewing angle as the parameter
+        # value and the jitter in the distribution, so be sure to set the
+        # empty pd for orientation parameters to 0.
+        pd = [value if relative or not parameter.polydisperse else 0.0], [1.0]
+    else:
+        limits = parameter.limits
+        disperser = values.get(parameter.name+'_pd_type', 'gaussian')
+        nsigma = values.get(parameter.name+'_pd_nsigma', 3.0)
+        pd = weights.get_weights(disperser, npts, width, nsigma,
+                                value, limits, relative)
+    return value, pd[0], pd[1]
+def _vol_pars(model_info, values):
     # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray]
     vol_pars = [get_weights(p, pars)
+    vol_pars = [_get_par_weights(p, values)
                 for p in model_info.parameters.call_parameters
                 if p.type == 'volume']
     #import pylab; pylab.plot(vol_pars[0][0],vol_pars[0][1]); pylab.show()
     value, weight = dispersion_mesh(model_info, vol_pars)
     return value, weight
+    dispersity, weight = dispersion_mesh(model_info, vol_pars)
+    return dispersity, weight

sasmodels/generate.py

-                      r6db17bd
+                      rff10479
 import string
 from zlib import crc32
+from inspect import currentframe, getframeinfo
 import numpy as np  # type: ignore
 …
 except ImportError:
     pass
+# jitter projection to use in the kernel code.  See explore/jitter.py
+# for details.  To change it from a program, set generate.PROJECTION.
+PROJECTION = 1
 def get_data_path(external_dir, target_file):
 …
     # Load templates and user code
     kernel_header = load_template('kernel_header.c')
+    dll_code = load_template('kernel_iq.c')
+    ocl_code = load_template('kernel_iq.cl')
+    #ocl_code = load_template('kernel_iq_local.cl')
+    kernel_code = load_template('kernel_iq.c')
     user_code = [(f, open(f).read()) for f in model_sources(model_info)]
+    # What kind of 2D model do we need?
+    xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str)
+               else 'qac' if not partable.is_asymmetric
+               else 'qabc')
     # Build initial sources
 …
     # Define the parameter table
+    # TODO: plug in current line number
+    source.append('#line 542 "sasmodels/generate.py"')
+    lineno = getframeinfo(currentframe()).lineno + 2
+    source.append('#line %d "sasmodels/generate.py"'%lineno)
+    #source.append('introduce breakage in generate to test lineno reporting')
     source.append("#define PARAMETER_TABLE \\")
     source.append("\\\n".join(p.as_definition()
 …
     source.append(call_volume)
+    refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters)
+    call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs))
+    if _have_Iqxy(user_code) or isinstance(model_info.Iqxy, str):
+        # Call 2D model
+        refs = ["_q[2*_i]", "_q[2*_i+1]"] + _call_pars("_v.", partable.iqxy_parameters)
+        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs))
+    else:
+        # Call 1D model with sqrt(qx^2 + qy^2)
+        #warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))")
+        # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters)
+        pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:]
+        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt))
+    model_refs = _call_pars("_v.", partable.iq_parameters)
+    pars = ",".join(["_q"] + model_refs)
+    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
+    if xy_mode == 'qabc':
+        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
+        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars
+        clear_iqxy = "#undef CALL_IQ_ABC"
+    elif xy_mode == 'qac':
+        pars = ",".join(["_qa", "_qc"] + model_refs)
+        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars
+        clear_iqxy = "#undef CALL_IQ_AC"
+    else:  # xy_mode == 'qa'
+        pars = ",".join(["_qa"] + model_refs)
+        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
+        clear_iqxy = "#undef CALL_IQ_A"
     magpars = [k-2 for k, p in enumerate(partable.call_parameters)
 …
     source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic)
     source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
+    for k, v in enumerate(magpars[:3]):
+        source.append("#define MAGNETIC_PAR%d %d"%(k+1, v))
+    source.append("#define PROJECTION %d"%PROJECTION)
     # TODO: allow mixed python/opencl kernels?
     ocl = kernels(ocl_code, call_iq, call_iqxy, model_info.name)
     dll = kernels(dll_code, call_iq, call_iqxy, model_info.name)
+    ocl = kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
+    dll = kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
     result = {
         'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]),
 …
 def kernels(kernel, call_iq, call_iqxy, name):
+def kernels(kernel, call_iq, call_iqxy, clear_iqxy, name):
     # type: ([str,str], str, str, str) -> List[str]
     code = kernel[0]
 …
         '#line 1 "%s Iqxy"' % path,
         code,
         "#undef CALL_IQ",
+        clear_iqxy,
         "#undef KERNEL_NAME",
+    ]
 …
         '#line 1 "%s Imagnetic"' % path,
         code,
+        clear_iqxy,
         "#undef MAGNETIC",
-        "#undef CALL_IQ",
         "#undef KERNEL_NAME",
+    ]

sasmodels/kernel_header.c

-                      rbb4b509
+                      r8698a0d
 #if defined(NEED_EXPM1)
+   // TODO: precision is a half digit lower than numpy on mac in [1e-7, 0.5]
+   // Run "explore/precision.py sas_expm1" to see this (may have to fiddle
+   // the xrange for log to see the complete range).
    static SAS_DOUBLE expm1(SAS_DOUBLE x_in) {
       double x = (double)x_in;  // go back to float for single precision kernels
 …
 inline double cube(double x) { return x*x*x; }
 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; }
-// To rotate from the canonical position to theta, phi, psi, first rotate by
-// psi about the major axis, oriented along z, which is a rotation in the
-// detector plane xy. Next rotate by theta about the y axis, aligning the major
-// axis in the xz plane. Finally, rotate by phi in the detector plane xy.
-// To compute the scattering, undo these rotations in reverse order:
-//     rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi
-// The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit
-// vector in the q direction.
-// To change between counterclockwise and clockwise rotation, change the
-// sign of phi and psi.
-#if 1
-//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017
-#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \
-    SINCOS(phi*M_PI_180, sn, cn); \
-    q = sqrt(qx*qx + qy*qy); \
-    cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \
-    sn = sqrt(1 - cn*cn); \
-    } while (0)
-#else
-// SasView 3.x definition of orientation
-#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \
-    SINCOS(theta*M_PI_180, sn, cn); \
-    q = sqrt(qx*qx + qy*qy);\
-    cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \
-    sn = sqrt(1 - cn*cn); \
-    } while (0)
-#endif
-#if 1
-#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \
-    q = sqrt(qx*qx + qy*qy); \
-    const double qxhat = qx/q; \
-    const double qyhat = qy/q; \
-    double sin_theta, cos_theta; \
-    double sin_phi, cos_phi; \
-    double sin_psi, cos_psi; \
-    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \
-    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \
-    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \
-    xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \
-         + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \
-    yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \
-         + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \
-    zhat = qxhat*(-sin_theta*cos_phi) \
-         + qyhat*(-sin_theta*sin_phi); \
-    } while (0)
-#else
-// SasView 3.x definition of orientation
-#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \
-    q = sqrt(qx*qx + qy*qy); \
-    const double qxhat = qx/q; \
-    const double qyhat = qy/q; \
-    double sin_theta, cos_theta; \
-    double sin_phi, cos_phi; \
-    double sin_psi, cos_psi; \
-    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \
-    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \
-    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \
-    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \
-    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \
-    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \
-    } while (0)
-#endif

sasmodels/kernel_iq.c

-                      rbde38b5
+                      r6aee3ab
 /*
     ##########################################################
 …
 */
+// NOTE: the following macros are defined in generate.py:
+//
+//  MAX_PD : the maximum number of dispersity loops allowed for this model,
+//      which will be at most modelinfo.MAX_PD.
+//  NUM_PARS : the number of parameters in the parameter table
+//  NUM_VALUES : the number of values to skip at the start of the
+//      values array before you get to the dispersity values.
+//  PARAMETER_TABLE : list of parameter declarations used to create the
+//      ParameterTable type.
+//  KERNEL_NAME : model_Iq, model_Iqxy or model_Imagnetic.  This code is
+//      included three times, once for each kernel type.
+//  MAGNETIC : defined when the magnetic kernel is being instantiated
+//  NUM_MAGNETIC : the number of magnetic parameters
+//  MAGNETIC_PARS : a comma-separated list of indices to the sld
+//      parameters in the parameter table.
+//  CALL_VOLUME(table) : call the form volume function
+//  CALL_IQ(q, table) : call the Iq function for 1D calcs.
+//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data.
+//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes
+//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes
+//  INVALID(table) : test if the current point is feesible to calculate.  This
+//      will be defined in the kernel definition file.
+//  PROJECTION : equirectangular=1, sinusoidal=2
+//      see explore/jitter.py for definitions.
 #ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
 #define _PAR_BLOCK_
 …
 typedef struct {
 #if MAX_PD > 0
     int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable
     int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector
+    int32_t pd_par[MAX_PD];     // id of the nth dispersity variable
+    int32_t pd_length[MAX_PD];  // length of the nth dispersity weight vector
     int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
     int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
 …
     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 first orientation variable
 } ProblemDetails;
 …
 #endif // _PAR_BLOCK_
+#if defined(MAGNETIC) && NUM_MAGNETIC>0
+#if defined(MAGNETIC) && NUM_MAGNETIC > 0
+// ===== Helper functions for magnetism =====
 // Return value restricted between low and high
 …
 //     uu * (sld - m_sigma_x);
 //     dd * (sld + m_sigma_x);
+//     ud * (m_sigma_y + 1j*m_sigma_z);
+//     du * (m_sigma_y - 1j*m_sigma_z);
+static void set_spins(double in_spin, double out_spin, double spins[4])
+//     ud * (m_sigma_y - 1j*m_sigma_z);
+//     du * (m_sigma_y + 1j*m_sigma_z);
+// weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag
+static void set_spin_weights(double in_spin, double out_spin, double spins[4])
+{
   in_spin = clip(in_spin, 0.0, 1.0);
   out_spin = clip(out_spin, 0.0, 1.0);
   spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd
   spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du
   spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud
+  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real
+  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real
   spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu
+}
+static double mag_sld(double qx, double qy, double p,
+                       double mx, double my, double sld)
+{
+  spins[4] = spins[1]; // du imag
+  spins[5] = spins[2]; // ud imag
+}
+// Compute the magnetic sld
+static double mag_sld(
+  const int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag
+  const double qx, const double qy,
+  const double px, const double py,
+  const double sld,
+  const double mx, const double my, const double mz
+)
+{
+  if (xs < 4) {
     const double perp = qy*mx - qx*my;
+    return sld + perp*p;
+}
+#endif // MAGNETIC
+    switch (xs) {
+      case 0: // uu => sld - D M_perpx
+          return sld - px*perp;
+      case 1: // ud real => -D M_perpy
+          return py*perp;
+      case 2: // du real => -D M_perpy
+          return py*perp;
+      case 3: // dd real => sld + D M_perpx
+          return sld + px*perp;
+    }
+  } else {
+    if (xs== 4) {
+      return -mz;  // ud imag => -D M_perpz
+    } else { // index == 5
+      return mz;   // du imag => D M_perpz
+    }
+  }
+}
+#endif
+// ===== Helper functions for orientation and jitter =====
+// To change the definition of the angles, run explore/angles.py, which
+// uses sympy to generate the equations.
+#if !defined(_QAC_SECTION) && defined(CALL_IQ_AC)
+#define _QAC_SECTION
+typedef struct {
+    double R31, R32;
+} QACRotation;
+// Fill in the rotation matrix R from the view angles (theta, phi) and the
+// jitter angles (dtheta, dphi).  This matrix can be applied to all of the
+// (qx, qy) points in the image to produce R*[qx,qy]' = [qa,qc]'
+static void
+qac_rotation(
+    QACRotation *rotation,
+    double theta, double phi,
+    double dtheta, double dphi)
+{
+    double sin_theta, cos_theta;
+    double sin_phi, cos_phi;
+    // reverse view matrix
+    SINCOS(theta*M_PI_180, sin_theta, cos_theta);
+    SINCOS(phi*M_PI_180, sin_phi, cos_phi);
+    const double V11 = cos_phi*cos_theta;
+    const double V12 = sin_phi*cos_theta;
+    const double V21 = -sin_phi;
+    const double V22 = cos_phi;
+    const double V31 = sin_theta*cos_phi;
+    const double V32 = sin_phi*sin_theta;
+    // reverse jitter matrix
+    SINCOS(dtheta*M_PI_180, sin_theta, cos_theta);
+    SINCOS(dphi*M_PI_180, sin_phi, cos_phi);
+    const double J31 = sin_theta;
+    const double J32 = -sin_phi*cos_theta;
+    const double J33 = cos_phi*cos_theta;
+    // reverse matrix
+    rotation->R31 = J31*V11 + J32*V21 + J33*V31;
+    rotation->R32 = J31*V12 + J32*V22 + J33*V32;
+}
+// Apply the rotation matrix returned from qac_rotation to the point (qx,qy),
+// returning R*[qx,qy]' = [qa,qc]'
+static double
+qac_apply(
+    QACRotation rotation,
+    double qx, double qy,
+    double *qa_out, double *qc_out)
+{
+    const double dqc = rotation.R31*qx + rotation.R32*qy;
+    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2
+    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy);
+    *qa_out = dqa;
+    *qc_out = dqc;
+}
+#endif // _QAC_SECTION
+#if !defined(_QABC_SECTION) && defined(CALL_IQ_ABC)
+#define _QABC_SECTION
+typedef struct {
+    double R11, R12;
+    double R21, R22;
+    double R31, R32;
+} QABCRotation;
+// Fill in the rotation matrix R from the view angles (theta, phi, psi) and the
+// jitter angles (dtheta, dphi, dpsi).  This matrix can be applied to all of the
+// (qx, qy) points in the image to produce R*[qx,qy]' = [qa,qb,qc]'
+static void
+qabc_rotation(
+    QABCRotation *rotation,
+    double theta, double phi, double psi,
+    double dtheta, double dphi, double dpsi)
+{
+    double sin_theta, cos_theta;
+    double sin_phi, cos_phi;
+    double sin_psi, cos_psi;
+    // reverse view matrix
+    SINCOS(theta*M_PI_180, sin_theta, cos_theta);
+    SINCOS(phi*M_PI_180, sin_phi, cos_phi);
+    SINCOS(psi*M_PI_180, sin_psi, cos_psi);
+    const double V11 = -sin_phi*sin_psi + cos_phi*cos_psi*cos_theta;
+    const double V12 = sin_phi*cos_psi*cos_theta + sin_psi*cos_phi;
+    const double V21 = -sin_phi*cos_psi - sin_psi*cos_phi*cos_theta;
+    const double V22 = -sin_phi*sin_psi*cos_theta + cos_phi*cos_psi;
+    const double V31 = sin_theta*cos_phi;
+    const double V32 = sin_phi*sin_theta;
+    // reverse jitter matrix
+    SINCOS(dtheta*M_PI_180, sin_theta, cos_theta);
+    SINCOS(dphi*M_PI_180, sin_phi, cos_phi);
+    SINCOS(dpsi*M_PI_180, sin_psi, cos_psi);
+    const double J11 = cos_psi*cos_theta;
+    const double J12 = sin_phi*sin_theta*cos_psi + sin_psi*cos_phi;
+    const double J13 = sin_phi*sin_psi - sin_theta*cos_phi*cos_psi;
+    const double J21 = -sin_psi*cos_theta;
+    const double J22 = -sin_phi*sin_psi*sin_theta + cos_phi*cos_psi;
+    const double J23 = sin_phi*cos_psi + sin_psi*sin_theta*cos_phi;
+    const double J31 = sin_theta;
+    const double J32 = -sin_phi*cos_theta;
+    const double J33 = cos_phi*cos_theta;
+    // reverse matrix
+    rotation->R11 = J11*V11 + J12*V21 + J13*V31;
+    rotation->R12 = J11*V12 + J12*V22 + J13*V32;
+    rotation->R21 = J21*V11 + J22*V21 + J23*V31;
+    rotation->R22 = J21*V12 + J22*V22 + J23*V32;
+    rotation->R31 = J31*V11 + J32*V21 + J33*V31;
+    rotation->R32 = J31*V12 + J32*V22 + J33*V32;
+}
+// Apply the rotation matrix returned from qabc_rotation to the point (qx,qy),
+// returning R*[qx,qy]' = [qa,qb,qc]'
+static double
+qabc_apply(
+    QABCRotation rotation,
+    double qx, double qy,
+    double *qa_out, double *qb_out, double *qc_out)
+{
+    *qa_out = rotation.R11*qx + rotation.R12*qy;
+    *qb_out = rotation.R21*qx + rotation.R22*qy;
+    *qc_out = rotation.R31*qx + rotation.R32*qy;
+}
+#endif // _QABC_SECTION
+// ==================== KERNEL CODE ========================
 kernel
 void KERNEL_NAME(
     int32_t nq,                 // number of q values
     const int32_t pd_start,     // where we are in the polydispersity loop
     const int32_t pd_stop,      // where we are stopping in the polydispersity loop
+    const int32_t pd_start,     // where we are in the dispersity loop
+    const int32_t pd_stop,      // where we are stopping in the dispersity loop
     global const ProblemDetails *details,
     global const double *values,
     global const double *q, // nq q values, with padding to boundary
     global double *result,  // nq+1 return values, again with padding
     const double cutoff     // cutoff in the polydispersity weight product
+    const double cutoff     // cutoff in the dispersity weight product
+    )
+{
+  // Storage for the current parameter values.  These will be updated as we
+  // walk the polydispersity cube.
+#ifdef USE_OPENCL
+  // who we are and what element we are working with
+  const int q_index = get_global_id(0);
+  if (q_index >= nq) return;
+#else
+  // Define q_index here so that debugging statements can be written to work
+  // for both OpenCL and DLL using:
+  //    if (q_index == 0) {printf(...);}
+  int q_index = 0;
+#endif
+  // ** Fill in the local values table **
+  // Storage for the current parameter values.
+  // These will be updated as we walk the dispersity mesh.
   ParameterBlock local_values;
-#if defined(MAGNETIC) && NUM_MAGNETIC>0
-  // Location of the sld parameters in the parameter vector.
-  // These parameters are updated with the effective sld due to magnetism.
-  #if NUM_MAGNETIC > 3
-  const int32_t slds[] = { MAGNETIC_PARS };
-  #endif
-  // TODO: could precompute these outside of the kernel.
-  // Interpret polarization cross section.
-  //     up_frac_i = values[NUM_PARS+2];
-  //     up_frac_f = values[NUM_PARS+3];
-  //     up_angle = values[NUM_PARS+4];
-  double spins[4];
-  double cos_mspin, sin_mspin;
-  set_spins(values[NUM_PARS+2], values[NUM_PARS+3], spins);
-  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
-#endif // MAGNETIC
-  // Fill in the initial variables
   //   values[0] is scale
   //   values[1] is background
 …
   for (int i=0; i < NUM_PARS; i++) {
     local_values.vector[i] = values[2+i];
 //printf("p%d = %g\n",i, local_values.vector[i]);
+    //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]);
+  }
+//printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
+//printf("start:%d  stop:%d\n", pd_start, pd_stop);
+  double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+  if (pd_start == 0) {
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
+  //if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
+  //if (q_index==0) printf("start:%d  stop:%d\n", pd_start, pd_stop);
+  // ** Precompute magnatism values **
+#if defined(MAGNETIC) && NUM_MAGNETIC>0
+  // Location of the sld parameters in the parameter vector.
+  // These parameters are updated with the effective sld due to magnetism.
+  const int32_t slds[] = { MAGNETIC_PARS };
+  // Interpret polarization cross section.
+  //     up_frac_i = values[NUM_PARS+2];
+  //     up_frac_f = values[NUM_PARS+3];
+  //     up_angle = values[NUM_PARS+4];
+  // TODO: could precompute more magnetism parameters before calling the kernel.
+  double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill
+  double cos_mspin, sin_mspin;
+  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins);
+  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
+#endif // MAGNETIC
+  // ** Fill in the initial results **
+  // If pd_start is zero that means that we are starting a new calculation,
+  // and must initialize the result to zero.  Otherwise, we are restarting
+  // the calculation from somewhere in the middle of the dispersity mesh,
+  // and we update the value rather than reset it. Similarly for the
+  // normalization factor, which is stored as the final value in the
+  // results vector (one past the number of q values).
+  //
+  // The code differs slightly between opencl and dll since opencl is only
+  // seeing one q value (stored in the variable "this_result") while the dll
+  // version must loop over all q.
+  #ifdef USE_OPENCL
+    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+    double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
+  #else // !USE_OPENCL
+    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+    if (pd_start == 0) {
+      #ifdef USE_OPENMP
+      #pragma omp parallel for
+      #endif
+      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
+    }
+    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
+#endif // !USE_OPENCL
+// ====== macros to set up the parts of the loop =======
+/*
+Based on the level of the loop, uses C preprocessor magic to construct
+level-specific looping variables, including these from loop level 3:
+  int n3 : length of loop for mesh level 3
+  int i3 : current position in the loop for level 3, which is calculated
+       from a combination of pd_start, pd_stride[3] and pd_length[3].
+  int p3 : is the index into the parameter table for mesh level 3
+  double v3[] : pointer into dispersity array to values for loop 3
+  double w3[] : pointer into dispersity array to weights for loop 3
+  double weight3 : the product of weights from levels 3 and up, computed
+       as weight5*weight4*w3[i3].  Note that we need an outermost
+       value weight5 set to 1.0 for this to work properly.
+After expansion, the loop struction will look like the following:
+  // --- PD_INIT(4) ---
+  const int n4 = pd_length[4];
+  const int p4 = pd_par[4];
+  global const double *v4 = pd_value + pd_offset[4];
+  global const double *w4 = pd_weight + pd_offset[4];
+  int i4 = (pd_start/pd_stride[4])%n4;  // position in level 4 at pd_start
+  // --- PD_INIT(3) ---
+  const int n3 = pd_length[3];
+  ...
+  int i3 = (pd_start/pd_stride[3])%n3;  // position in level 3 at pd_start
+  PD_INIT(2)
+  PD_INIT(1)
+  PD_INIT(0)
+  // --- PD_OUTERMOST_WEIGHT(5) ---
+  const double weight5 = 1.0;
+  // --- PD_OPEN(4,5) ---
+  while (i4 < n4) {
+    parameter[p4] = v4[i4];  // set the value for pd parameter 4 at this mesh point
+    const double weight4 = w4[i4] * weight5;
+    // from PD_OPEN(3,4)
+    while (i3 < n3) {
+      parameter[p3] = v3[i3];  // set the value for pd parameter 3 at this mesh point
+      const double weight3 = w3[i3] * weight4;
+      PD_OPEN(3,2)
+      PD_OPEN(2,1)
+      PD_OPEN(0,1)
+      // ... main loop body ...
+      APPLY_PROJECTION    // convert jitter values to spherical coords
+      BUILD_ROTATION      // construct the rotation matrix qxy => qabc
+      for each q
+          FETCH_Q         // set qx,qy from the q input vector
+          APPLY_ROTATION  // convert qx,qy to qa,qb,qc
+          CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...)
+      ++step;  // increment counter representing position in dispersity mesh
+      PD_CLOSE(0)
+      PD_CLOSE(1)
+      PD_CLOSE(2)
+      // --- PD_CLOSE(3) ---
+      if (step >= pd_stop) break;
+      ++i3;
+    }
+    i3 = 0; // reset loop counter for next round through the loop
+    // --- PD_CLOSE(4) ---
+    if (step >= pd_stop) break;
+    ++i4;
+  }
+//printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
+  i4 = 0; // reset loop counter even though no more rounds through the loop
+*/
+// ** prepare inner loops **
+// Depending on the shape type (radial, axial, triaxial), the variables
+// and calling parameters in the loop body will be slightly different.
+// Macros capture the differences in one spot so the rest of the code
+// is easier to read. The code below both declares variables for the
+// inner loop and defines the macros that use them.
+#if defined(CALL_IQ)
+  // unoriented 1D
+  double qk;
+  #define FETCH_Q() do { qk = q[q_index]; } while (0)
+  #define BUILD_ROTATION() do {} while(0)
+  #define APPLY_ROTATION() do {} while(0)
+  #define CALL_KERNEL() CALL_IQ(qk, local_values.table)
+#elif defined(CALL_IQ_A)
+  // unoriented 2D
+  double qx, qy;
+  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
+  #define BUILD_ROTATION() do {} while(0)
+  #define APPLY_ROTATION() do {} while(0)
+  #define CALL_KERNEL() CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table)
+#elif defined(CALL_IQ_AC)
+  // oriented symmetric 2D
+  double qx, qy;
+  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
+  double qa, qc;
+  QACRotation rotation;
+  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
+  #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi);
+  #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc)
+  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table)
+#elif defined(CALL_IQ_ABC)
+  // oriented asymmetric 2D
+  double qx, qy;
+  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
+  double qa, qb, qc;
+  QABCRotation rotation;
+  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
+  // psi and dpsi are only for IQ_ABC, so they are processed here.
+  const double psi = values[details->theta_par+4];
+  local_values.table.psi = 0.;
+  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi)
+  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc)
+  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
+#endif
+// Doing jitter projection code outside the previous if block so that we don't
+// need to repeat the identical logic in the IQ_AC and IQ_ABC branches.  This
+// will become more important if we implement more projections, or more
+// complicated projections.
+#if defined(CALL_IQ) || defined(CALL_IQ_A)
+  #define APPLY_PROJECTION() const double weight=weight0
+#else // !spherosymmetric projection
+  // Grab the "view" angles (theta, phi, psi) from the initial parameter table.
+  const double theta = values[details->theta_par+2];
+  const double phi = values[details->theta_par+3];
+  // Make sure jitter angle defaults to zero if there is no jitter distribution
+  local_values.table.theta = 0.;
+  local_values.table.phi = 0.;
+  // The "jitter" angles (dtheta, dphi, dpsi) are stored with the
+  // dispersity values and copied to the local parameter table as
+  // we go through the mesh.
+  double dtheta, dphi, weight;
+  #if PROJECTION == 1
+    #define APPLY_PROJECTION() do { \
+      dtheta = local_values.table.theta; \
+      dphi = local_values.table.phi; \
+      weight = fabs(cos(dtheta*M_PI_180)) * weight0; \
+    } while (0)
+  #elif PROJECTION == 2
+    #define APPLY_PROJECTION() do { \
+      dtheta = local_values.table.theta; \
+      dphi = local_values.table.phi; \
+      weight = weight0; \
+      if (dtheta != 90.0) dphi /= cos(dtheta*M_PI_180); \
+      else if (dphi != 0.0) weight = 0.; \
+      if (fabs(dphi) >= 180.) weight = 0.; \
+    } while (0)
+  #endif
+#endif // !spherosymmetric projection
+// ** define looping macros **
+// Define looping variables
+#define PD_INIT(_LOOP) \
+  const int n##_LOOP = details->pd_length[_LOOP]; \
+  const int p##_LOOP = details->pd_par[_LOOP]; \
+  global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
+  global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
+  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP;
+// Jump into the middle of the dispersity loop
+#define PD_OPEN(_LOOP,_OUTER) \
+  while (i##_LOOP < n##_LOOP) { \
+    local_values.vector[p##_LOOP] = v##_LOOP[i##_LOOP]; \
+    const double weight##_LOOP = w##_LOOP[i##_LOOP] * weight##_OUTER;
+// create the variable "weight#=1.0" where # is the outermost level+1 (=MAX_PD).
+#define _PD_OUTERMOST_WEIGHT(_n) const double weight##_n = 1.0;
+#define PD_OUTERMOST_WEIGHT(_n) _PD_OUTERMOST_WEIGHT(_n)
+// Close out the loop
+#define PD_CLOSE(_LOOP) \
+    if (step >= pd_stop) break; \
+    ++i##_LOOP; \
+  } \
+  i##_LOOP = 0;
+// ====== construct the loops =======
+// Pointers to the start of the dispersity and weight vectors, if needed.
 #if MAX_PD>0
   global const double *pd_value = values + NUM_VALUES;
 …
 #endif
+  // Jump into the middle of the polydispersity loop
+// The variable "step" is the current position in the dispersity loop.
+// It will be incremented each time a new point in the mesh is accumulated,
+// and used to test whether we have reached pd_stop.
+int step = pd_start;
+// *** define loops for each of 0, 1, 2, ..., modelinfo.MAX_PD-1 ***
+// define looping variables
 #if MAX_PD>4
+  int n4=details->pd_length[4];
+  int i4=(pd_start/details->pd_stride[4])%n4;
+  const int p4=details->pd_par[4];
+  global const double *v4 = pd_value + details->pd_offset[4];
+  global const double *w4 = pd_weight + details->pd_offset[4];
+  PD_INIT(4)
 #endif
 #if MAX_PD>3
+  int n3=details->pd_length[3];
+  int i3=(pd_start/details->pd_stride[3])%n3;
+  const int p3=details->pd_par[3];
+  global const double *v3 = pd_value + details->pd_offset[3];
+  global const double *w3 = pd_weight + details->pd_offset[3];
+//printf("offset %d: %d %d\n", 3, details->pd_offset[3], NUM_VALUES);
+  PD_INIT(3)
 #endif
 #if MAX_PD>2
+  int n2=details->pd_length[2];
+  int i2=(pd_start/details->pd_stride[2])%n2;
+  const int p2=details->pd_par[2];
+  global const double *v2 = pd_value + details->pd_offset[2];
+  global const double *w2 = pd_weight + details->pd_offset[2];
+  PD_INIT(2)
 #endif
 #if MAX_PD>1
+  int n1=details->pd_length[1];
+  int i1=(pd_start/details->pd_stride[1])%n1;
+  const int p1=details->pd_par[1];
+  global const double *v1 = pd_value + details->pd_offset[1];
+  global const double *w1 = pd_weight + details->pd_offset[1];
+  PD_INIT(1)
 #endif
 #if MAX_PD>0
+  int n0=details->pd_length[0];
+  int i0=(pd_start/details->pd_stride[0])%n0;
+  const int p0=details->pd_par[0];
+  global const double *v0 = pd_value + details->pd_offset[0];
+  global const double *w0 = pd_weight + details->pd_offset[0];
+//printf("w0:%p, values:%p, diff:%ld, %d\n",w0,values,(w0-values), NUM_VALUES);
+#endif
+  PD_INIT(0)
+#endif
+// open nested loops
+PD_OUTERMOST_WEIGHT(MAX_PD)
+#if MAX_PD>4
+  PD_OPEN(4,5)
+#endif
+#if MAX_PD>3
+  PD_OPEN(3,4)
+#endif
+#if MAX_PD>2
+  PD_OPEN(2,3)
+#endif
+#if MAX_PD>1
+  PD_OPEN(1,2)
+#endif
 #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;
+#if MAX_PD>4
+  const double weight5 = 1.0;
+  while (i4 < n4) {
+    local_values.vector[p4] = v4[i4];
+    double weight4 = w4[i4] * weight5;
+//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, local_values.vector[p4], weight4);
+#elif MAX_PD>3
+    const double weight4 = 1.0;
+#endif
+#if MAX_PD>3
+  while (i3 < n3) {
+    local_values.vector[p3] = v3[i3];
+    double weight3 = w3[i3] * weight4;
+//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, local_values.vector[p3], weight3);
+#elif MAX_PD>2
+    const double weight3 = 1.0;
+#endif
+#if MAX_PD>2
+  while (i2 < n2) {
+    local_values.vector[p2] = v2[i2];
+    double weight2 = w2[i2] * weight3;
+//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, local_values.vector[p2], weight2);
+#elif MAX_PD>1
+    const double weight2 = 1.0;
+#endif
+#if MAX_PD>1
+  while (i1 < n1) {
+    local_values.vector[p1] = v1[i1];
+    double weight1 = w1[i1] * weight2;
+//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, local_values.vector[p1], weight1);
+#elif MAX_PD>0
+    const double weight1 = 1.0;
+#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;
+#endif
+//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");
+//printf("sphcor: %g\n", spherical_correction);
+    #ifdef INVALID
+    if (!INVALID(local_values.table))
+    #endif
+    {
+      // Accumulate I(q)
+      // 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);
+        #ifdef USE_OPENMP
+        #pragma omp parallel for
+        #endif
+        for (int q_index=0; q_index<nq; q_index++) {
+#if defined(MAGNETIC) && NUM_MAGNETIC > 0
+          const double qx = q[2*q_index];
+          const double qy = q[2*q_index+1];
+  PD_OPEN(0,1)
+#endif
+//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");}
+  // ====== loop body =======
+  #ifdef INVALID
+  if (!INVALID(local_values.table))
+  #endif
+  {
+     APPLY_PROJECTION();
+    // Accumulate I(q)
+    // Note: weight==0 must always be excluded
+    if (weight > cutoff) {
+      pd_norm += weight * CALL_VOLUME(local_values.table);
+      BUILD_ROTATION();
+#ifndef USE_OPENCL
+      // DLL needs to explicitly loop over the q values.
+      #ifdef USE_OPENMP
+      #pragma omp parallel for
+      #endif
+      for (q_index=0; q_index<nq; q_index++)
+#endif // !USE_OPENCL
+      {
+        FETCH_Q();
+        APPLY_ROTATION();
+        // ======= COMPUTE SCATTERING ==========
+        #if defined(MAGNETIC) && NUM_MAGNETIC > 0
+          // Compute the scattering from the magnetic cross sections.
+          double scattering = 0.0;
           const double qsq = qx*qx + qy*qy;
-          // Constant across orientation, polydispersity for given qx, qy
-          double scattering = 0.0;
-          // TODO: what is the magnetic scattering at q=0
           if (qsq > 1.e-16) {
+            double p[4];  // dd, du, ud, uu
+            p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq;
+            p[3] = -p[0];
+            p[1] = p[2] = (qy*sin_mspin - qx*cos_mspin)/qsq;
+            for (int index=0; index<4; index++) {
+              const double xs = spins[index];
+              if (xs > 1.e-8) {
+                const int spin_flip = (index==1) || (index==2);
+                const double pk = p[index];
+                for (int axis=0; axis<=spin_flip; axis++) {
+                  #define M1 NUM_PARS+5
+                  #define M2 NUM_PARS+8
+                  #define M3 NUM_PARS+13
+                  #define SLD(_M_offset, _sld_offset) \
+                      local_values.vector[_sld_offset] = xs * (axis \
+                      ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \
+                      : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \
+                                (spin_flip ? 0.0 : values[_sld_offset+2])))
+                  #if NUM_MAGNETIC==1
+                      SLD(M1, MAGNETIC_PAR1);
+                  #elif NUM_MAGNETIC==2
+                      SLD(M1, MAGNETIC_PAR1);
+                      SLD(M2, MAGNETIC_PAR2);
+                  #elif NUM_MAGNETIC==3
+                      SLD(M1, MAGNETIC_PAR1);
+                      SLD(M2, MAGNETIC_PAR2);
+                      SLD(M3, MAGNETIC_PAR3);
+                  #else
+                  for (int sk=0; sk<NUM_MAGNETIC; sk++) {
+                      SLD(M1+3*sk, slds[sk]);
+                  }
+                  #endif
+                  scattering += CALL_IQ(q, q_index, local_values.table);
+            // TODO: what is the magnetic scattering at q=0
+            const double px = (qy*cos_mspin + qx*sin_mspin)/qsq;
+            const double py = (qy*sin_mspin - qx*cos_mspin)/qsq;
+            // loop over uu, ud real, du real, dd, ud imag, du imag
+            for (int xs=0; xs<6; xs++) {
+              const double xs_weight = spins[xs];
+              if (xs_weight > 1.e-8) {
+                // Since the cross section weight is significant, set the slds
+                // to the effective slds for this cross section, call the
+                // kernel, and add according to weight.
+                for (int sk=0; sk<NUM_MAGNETIC; sk++) {
+                  const int32_t mag_index = NUM_PARS+5 + 3*sk;
+                  const int32_t sld_index = slds[sk];
+                  const double mx = values[mag_index];
+                  const double my = values[mag_index+1];
+                  const double mz = values[mag_index+2];
+                  local_values.vector[sld_index] =
+                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz);
+                }
+                scattering += xs_weight * CALL_KERNEL();
+              }
+            }
+          }
+#else  // !MAGNETIC
+          const double scattering = CALL_IQ(q, q_index, local_values.table);
+#endif // !MAGNETIC
+//printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0);
+        #else  // !MAGNETIC
+          const double scattering = CALL_KERNEL();
+        #endif // !MAGNETIC
+//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0);
+        #ifdef USE_OPENCL
+          this_result += weight * scattering;
+        #else // !USE_OPENCL
           result[q_index] += weight * scattering;
+        }
+        #endif // !USE_OPENCL
+      }
+    }
+    ++step;
+  }
+// close nested loops
+++step;
 #if MAX_PD>0
+    if (step >= pd_stop) break;
+    ++i0;
+  }
+  i0 = 0;
+  PD_CLOSE(0)
 #endif
 #if MAX_PD>1
+    if (step >= pd_stop) break;
+    ++i1;
+  }
+  i1 = 0;
+  PD_CLOSE(1)
 #endif
 #if MAX_PD>2
+    if (step >= pd_stop) break;
+    ++i2;
+  }
+  i2 = 0;
+  PD_CLOSE(2)
 #endif
 #if MAX_PD>3
+    if (step >= pd_stop) break;
+    ++i3;
+  }
+  i3 = 0;
+  PD_CLOSE(3)
 #endif
 #if MAX_PD>4
+    if (step >= pd_stop) break;
+    ++i4;
+  }
+  i4 = 0;
+#endif
+  PD_CLOSE(4)
+#endif
+// Remember the current result and the updated norm.
+#ifdef USE_OPENCL
+  result[q_index] = this_result;
+  if (q_index == 0) result[nq] = pd_norm;
+//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm);
+#else // !USE_OPENCL
+  result[nq] = pd_norm;
 //printf("res: %g/%g\n", result[0], pd_norm);
+  // Remember the updated norm.
+  result[nq] = pd_norm;
+}
+#endif // !USE_OPENCL
+// ** clear the macros in preparation for the next kernel **
+#undef PD_INIT
+#undef PD_OPEN
+#undef PD_CLOSE
+#undef FETCH_Q
+#undef APPLY_PROJECTION
+#undef BUILD_ROTATION
+#undef APPLY_ROTATION
+#undef CALL_KERNEL
+}

sasmodels/kernelpy.py

-                      r3b32f8d
+                      r8698a0d
         partable = model_info.parameters
+        kernel_parameters = (partable.iqxy_parameters if q_input.is_2d
+                             else partable.iq_parameters)
+        #kernel_parameters = (partable.iqxy_parameters if q_input.is_2d
+        #                     else partable.iq_parameters)
+        kernel_parameters = partable.iq_parameters
         volume_parameters = partable.form_volume_parameters
 …
     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
 …
             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:
 …
             # update value and norm
-            weight *= spherical_correction
             total += weight * Iq
             pd_norm += weight * form_volume()

sasmodels/model_test.py

-                      r65314f7
+                      r20fe0cd
     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 = []
 …
+                ]
             tests = smoke_tests
+            #tests = []
             if self.info.tests is not None:
                 tests += self.info.tests

sasmodels/modelinfo.py

-                      ra85a569
+                      r6aee3ab
     TestCondition = Tuple[ParameterSetUser, TestInput, TestValue]
+MAX_PD = 4 #: Maximum number of simultaneously polydisperse parameters
+# If MAX_PD changes, need to change the loop macros in kernel_iq.c
+MAX_PD = 5 #: Maximum number of simultaneously polydisperse parameters
 # assumptions about common parameters exist throughout the code, such as:
 …
       with vector parameter p sent as p[].
     * *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...)
+    * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...)
       function, with vector parameter p sent as p[].
 …
         # type: (List[Parameter]) -> None
         self.kernel_parameters = parameters
+        self._check_angles()
         self._set_vector_lengths()
 …
         self.iq_parameters = [p for p in self.kernel_parameters
                               if p.type not in ('orientation', 'magnetic')]
+        self.iqxy_parameters = [p for p in self.kernel_parameters
+                                if p.type != 'magnetic']
+        # note: orientation no longer sent to Iqxy, so its the same as
+        #self.iqxy_parameters = [p for p in self.kernel_parameters
+        #                        if p.type != 'magnetic']
         self.form_volume_parameters = [p for p in self.kernel_parameters
                                        if p.type == 'volume']
 …
         self.has_2d = any(p.type in ('orientation', 'magnetic')
                           for p in self.kernel_parameters)
+        self.is_asymmetric = any(p.name == 'psi' for p in self.kernel_parameters)
         self.magnetism_index = [k for k, p in enumerate(self.call_parameters)
                                 if p.id.startswith('M0:')]
 …
                          if p.polydisperse and p.type not in ('orientation', 'magnetic'))
         self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse)
+    def _check_angles(self):
+        theta = phi = psi = -1
+        for k, p in enumerate(self.kernel_parameters):
+            if p.name == 'theta':
+                theta = k
+                if p.type != 'orientation':
+                    raise TypeError("theta must be an orientation parameter")
+            elif p.name == 'phi':
+                phi = k
+                if p.type != 'orientation':
+                    raise TypeError("phi must be an orientation parameter")
+            elif p.name == 'psi':
+                psi = k
+                if p.type != 'orientation':
+                    raise TypeError("psi must be an orientation parameter")
+        if theta >= 0 and phi >= 0:
+            if phi != theta+1:
+                raise TypeError("phi must follow theta")
+            if psi >= 0 and psi != phi+1:
+                raise TypeError("psi must follow phi")
+        elif theta >= 0 or phi >= 0 or psi >= 0:
+            raise TypeError("oriented shapes must have both theta and phi and maybe psi")
     def __getitem__(self, key):
 …
             raise KeyError("unknown parameter %r"%key)
         return par
+    def __contains__(self, key):
+        for par in self.call_parameters:
+            if par.name == key:
+                return True
+        else:
+            return False
     def _set_vector_lengths(self):

sasmodels/models/barbell.c

-                      r592343f
+                      rbecded3
-double form_volume(double radius_bell, double radius, double length);
-double Iq(double q, double sld, double solvent_sld,
-        double radius_bell, double radius, double length);
-double Iqxy(double qx, double qy, double sld, double solvent_sld,
-        double radius_bell, double radius, double length,
-        double theta, double phi);
 #define INVALID(v) (v.radius_bell < v.radius)
 //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]
 …
     //    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;
 …
 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;
 …
+}
 double form_volume(double radius_bell,
         double radius,
         double length)
+static double
+form_volume(double radius_bell,
+    double radius,
+    double length)
+{
     // bell radius should never be less than radius when this is called
 …
+}
+double Iq(double q, double sld, double solvent_sld,
+          double radius_bell, double radius, double length)
+static double
+Iq(double q, double sld, double solvent_sld,
+    double radius_bell, double radius, double length)
+{
     const double h = -sqrt(radius_bell*radius_bell - radius*radius);
 …
         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;
+    }
 …
+double Iqxy(double qx, double qy,
         double sld, double solvent_sld,
         double radius_bell, double radius, double length,
         double theta, double phi)
+static double
+Iqxy(double qab, double qc,
+    double sld, double solvent_sld,
+    double radius_bell, double radius, double length)
+{
-    double q, sin_alpha, cos_alpha;
-    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, 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

-                      r50beefe
+                      rea60e08
+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)
+{
+    // Equations from Matsuoka 26-27-28, multiplied by |q|
+    const double a1 = (-qa + qb + qc)/2.0;
+    const double a2 = (+qa - qb + qc)/2.0;
+    const double a3 = (+qa + qb - qc)/2.0;
+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
+    // Matsuoka 29-30-31
+    //     Z_k numerator: 1 - exp(a)^2
+    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a)
+    // Rewriting numerator
+    //         => -(exp(2a) - 1)
+    //         => -expm1(2a)
+    // Rewriting denominator
+    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => (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));
+#elif 0
+    // ** Alternate form, which perhaps is more approachable
+    //     Z_k numerator   => -[(exp(2a) - 1) / 2.exp(a)] 2.exp(a)
+    //                     => -[sinh(a)] exp(a)
+    //     Z_k denominator => [(exp(2a) + 1) / 2.exp(a) - cos(d a_k)] 2.exp(a)
+    //                     => [cosh(a) - cos(d a_k)] 2.exp(a)
+    //     => Z_k = -sinh(a) / [cosh(a) - cos(d a_k)]
+    //            = sinh(-a) / [cosh(-a) - cos(d a_k)]
+    //
+    // One more step leads to the form in sasview 3.x for 2d models
+    //            = tanh(-a) / [1 - cos(d a_k)/cosh(-a)]
+    //
+    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));
+#else
+    const double arg = 0.5*square(dnn*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.0 - cos(dnn*a1)/cosh_qd)
+                    * tanh_qd/(1.0 - cos(dnn*a2)/cosh_qd)
+                    * tanh_qd/(1.0 - cos(dnn*a3)/cosh_qd);
+#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,
+  double d_factor, double radius,
+  double sld, double solvent_sld){
+static double Iq(double q, double dnn,
+    double d_factor, double radius,
+    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 qa, double qb, double qc,
     double dnn, double d_factor, double radius,
+    double sld, double solvent_sld,
+    double theta, double phi, double psi)
+    double sld, double solvent_sld)
+{
+    double q, zhat, yhat, xhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, 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;
+    const double 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/bcc_paracrystal.py

-                      r8f04da4
+                      reda8b30
     \end{array}
+**NB**: The calculation of $Z(q)$ is a double numerical integral that must
+be carried out with a high density of points to properly capture the sharp
+peaks of the paracrystalline scattering. So be warned that the calculation
+is SLOW. Go get some coffee. Fitting of any experimental data must be
+resolution smeared for any meaningful fit. This makes a triple integral.
+Very, very slow. Go get lunch!
+.. note::
+  The calculation of $Z(q)$ is a double numerical integral that
+  must be carried out with a high density of points to properly capture
+  the sharp peaks of the paracrystalline scattering.
+  So be warned that the calculation is slow. Fitting of any experimental data
+  must be resolution smeared for any meaningful fit. This makes a triple integral
+  which may be very slow.
 This example dataset is produced using 200 data points,
 *qmin* = 0.001 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above default values.
 …
 The 2D (Anisotropic model) is based on the reference below where $I(q)$ is
 approximated for 1d scattering. Thus the scattering pattern for 2D may not
+be accurate.
+be accurate, particularly at low $q$. For general details of the calculation and angular
+dispersions for oriented particles see :ref:`orientation` .
+Note that we are not responsible for any incorrectness of the 2D model computation.
 .. figure:: img/parallelepiped_angle_definition.png

sasmodels/models/capped_cylinder.c

-                      r592343f
+                      rbecded3
-double form_volume(double radius, double radius_cap, double length);
-double Iq(double q, double sld, double solvent_sld,
-    double radius, double radius_cap, double length);
-double Iqxy(double qx, double qy, double sld, double solvent_sld,
-    double radius, double radius_cap, double length, double theta, double phi);
 #define INVALID(v) (v.radius_cap < v.radius)
 …
 //   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]
 …
     // 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;
 …
 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;
 …
+}
+double form_volume(double radius, double radius_cap, double length)
+static double
+form_volume(double radius, double radius_cap, double length)
+{
     // cap radius should never be less than radius when this is called
 …
+}
+double Iq(double q, double sld, double solvent_sld,
+          double radius, double radius_cap, double length)
+static double
+Iq(double q, double sld, double solvent_sld,
+    double radius, double radius_cap, double length)
+{
     const double h = sqrt(radius_cap*radius_cap - radius*radius);
 …
     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]
 …
+double Iqxy(double qx, double qy,
+static double
+Iqxy(double qab, double qc,
     double sld, double solvent_sld, double radius,
+    double radius_cap, double length,
+    double theta, double phi)
+    double radius_cap, double length)
+{
-    double q, sin_alpha, cos_alpha;
-    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, 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

-                      rb260926
+                      rbecded3
+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 +
 …
                      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
 …
     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 qab, double qc,
+    double radius,
+    double thick_rim,
+    double thick_face,
+    double length,
+    double core_sld,
+    double face_sld,
+    double rim_sld,
+    double solvent_sld)
+{
+    double q, sin_alpha, cos_alpha;
+    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    double answer = bicelle_kernel(q, radius, thick_rim, thick_face,
+    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face,
 .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

-                      rdedcf34
+                      r82592da
 static double
 form_volume(double r_minor,
         double x_core,
         double thick_rim,
         double thick_face,
         double length)
+    double x_core,
+    double thick_rim,
+    double thick_face,
+    double length)
+{
     return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face);
 …
 static double
 Iq(double q,
         double r_minor,
         double x_core,
         double thick_rim,
         double thick_face,
         double length,
         double rhoc,
         double rhoh,
         double rhor,
         double rhosolv)
+    double r_minor,
+    double x_core,
+    double thick_rim,
+    double thick_face,
+    double length,
+    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
 …
     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
 …
 static double
+Iqxy(double qx, double qy,
+          double r_minor,
+          double x_core,
+          double thick_rim,
+          double thick_face,
+          double length,
+          double rhoc,
+          double rhoh,
+          double rhor,
+          double rhosolv,
+          double theta,
+          double phi,
+          double psi)
+Iqxy(double qa, double qb, double qc,
+    double r_minor,
+    double x_core,
+    double thick_rim,
+    double thick_face,
+    double length,
+    double sld_core,
+    double sld_face,
+    double sld_rim,
+    double sld_solvent)
+{
+       // 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 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;
 …
     // 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*qb) + square(r_minor*qa));
+    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qb)
+                                   + square((r_minor+thick_rim)*qa));
+    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

-                      r592343f
+                      rbecded3
-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,
 …
+{
     // 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]
 …
 double Iqxy(double qx, double qy,
+double Iqxy(double qab, double qc,
     double core_sld,
     double shell_sld,
 …
     double radius,
     double thickness,
+    double length,
+    double theta,
+    double phi)
+    double length)
+{
+    double q, sin_alpha, cos_alpha;
+    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, 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

-                      r0a3d9b2
+                      rbecded3
-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;
 …
 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;
 …
     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 qab, double qc,
+    double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell,
+    double core_sld,
+    double shell_sld,
+    double solvent_sld)
+{
+    double q, sin_alpha, cos_alpha;
+    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, 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;
 …
     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

-                      r30b60d2
+                      r8db25bf
 # 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

-                      rc69d6d6
+                      r904cd9c
+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,
+          double solvent_sld, double length_a, double length_b, double length_c,
+          double thick_rim_a, double thick_rim_b, double thick_rim_c);
+double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld,
+            double crim_sld, double solvent_sld, double length_a, double length_b,
+            double length_c, double thick_rim_a, double thick_rim_b,
+            double thick_rim_c, double theta, double phi, double psi);
+double form_volume(double length_a, double length_b, double length_c,
+                   double thick_rim_a, double thick_rim_b, double thick_rim_c)
+static double
+form_volume(double length_a, double length_b, double length_c,
+    double thick_rim_a, double thick_rim_b, double thick_rim_c)
+{
     //return length_a * length_b * length_c;
 …
+}
+double Iq(double q,
+static double
+Iq(double q,
     double core_sld,
     double arim_sld,
 …
     double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
     double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present
-    double Vot = Vin + V1 + V2 + V3;
     // Scale factors (note that drC is not used later)
 …
             const double form_crim = scale11*si1*si2;
             //  correct FF : sum of square of phase factors
             inner_total += Gauss76Wt[j] * form * form;
 …
+}
+double Iqxy(double qx, double qy,
+static double
+Iqxy(double qa, double qb, double qc,
     double core_sld,
     double arim_sld,
 …
     double thick_rim_a,
     double thick_rim_b,
+    double thick_rim_c,
+    double theta,
+    double phi,
+    double psi)
+    double thick_rim_c)
+{
-    double q, zhat, yhat, xhat;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
     // cspkernel in csparallelepiped recoded here
     const double dr0 = core_sld-solvent_sld;
 …
     double tc = length_c + 2.0*thick_rim_c;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
     double siA = sas_sinx_x(0.5*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);

sasmodels/models/core_shell_parallelepiped.py

-                      r8f04da4
+                      r904cd9c
 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied.
+To provide easy access to the orientation of the parallelepiped, we define the
+axis of the cylinder using three angles $\theta$, $\phi$ and $\Psi$.
+(see :ref:`cylinder orientation <cylinder-angle-definition>`).
+The angle $\Psi$ is the rotational angle around the *long_c* axis against the
+$q$ plane. For example, $\Psi = 0$ when the *short_b* axis is parallel to the
+*x*-axis of the detector.
+For 2d data the orientation of the particle is required, described using
+angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
+of the calculation and angular dispersions see :ref:`orientation` .
+The angle $\Psi$ is the rotational angle around the *long_c* axis. For example,
+$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector.
 .. figure:: img/parallelepiped_angle_definition.png
     Definition of the angles for oriented core-shell parallelepipeds.
+    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then
+    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder.
+    The neutron or X-ray beam is along the $z$ axis.
 .. figure:: img/parallelepiped_angle_projection.png

sasmodels/models/cylinder.c

-                      r592343f
+                      rbecded3
-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]
 …
+}
+double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double solvent_sld,
 …
+}
 double Iqxy(double qx, double qy,
+static double
+Iqxy(double qab, double qc,
     double sld,
     double solvent_sld,
     double radius,
+    double length,
+    double theta,
+    double phi)
+    double length)
+{
-    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 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/cylinder.py

-                      r31df0c9
+                      reda8b30
 when $P(q) \cdot S(q)$ is applied.
 For oriented cylinders, we define the direction of the
+For 2d scattering from oriented cylinders, we define the direction of the
 axis of the cylinder using two angles $\theta$ (note this is not the
 same as the scattering angle used in q) and $\phi$. Those angles
 are defined in :numref:`cylinder-angle-definition` .
+are defined in :numref:`cylinder-angle-definition` , for further details see :ref:`orientation` .
 .. _cylinder-angle-definition:
 …
 .. figure:: img/cylinder_angle_definition.png
+    Definition of the $\theta$ and $\phi$ orientation angles for a cylinder relative
+    to the beam line coordinates, plus an indication of their orientation distributions
+    which are described as rotations about each of the perpendicular axes $\delta_1$ and $\delta_2$
+    Angles $\theta$ and $\phi$ orient the cylinder relative
+    to the beam line coordinates, where the beam is along the $z$ axis. Rotation $\theta$, initially
+    in the $xz$ plane, is carried out first, then rotation $\phi$ about the $z$ axis. Orientation distributions
+    are described as rotations about two perpendicular axes $\delta_1$ and $\delta_2$
     in the frame of the cylinder itself, which when $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes.
 …
 The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.
-On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
-appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, which when $\theta = \phi = 0$ are parallel
-to the $Y$ and $X$ axes of the instrument respectively. Some experimentation may be required to understand the 2d patterns fully.
-(Earlier implementations had numerical integration issues in some circumstances when orientation distributions passed through 90 degrees, such
-situations, with very broad distributions, should still be approached with care.)
 Validation

sasmodels/models/ellipsoid.c

-                      r3b571ae
+                      rbecded3
+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,
 …
+}
+double Iqxy(double qx, double qy,
+static double
+Iqxy(double qab, double qc,
     double sld,
     double sld_solvent,
     double radius_polar,
+    double radius_equatorial,
+    double theta,
+    double phi)
+    double radius_equatorial)
+{
+    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 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/ellipsoid.py

-                      r92708d8
+                      reda8b30
     r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2}
+To provide easy access to the orientation of the ellipsoid, we define
+the rotation axis of the ellipsoid using two angles $\theta$ and $\phi$.
+These angles are defined in the
+For 2d data from oriented ellipsoids the direction of the rotation axis of
+the ellipsoid is defined using two angles $\theta$ and $\phi$ as for the
 :ref:`cylinder orientation figure <cylinder-angle-definition>`.
 For the ellipsoid, $\theta$ is the angle between the rotational axis
 and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$
+in the $xy$ plane.
+in the $xy$ plane, for further details of the calculation and angular
+dispersions see :ref:`orientation` .
 NB: The 2nd virial coefficient of the solid ellipsoid is calculated based

sasmodels/models/elliptical_cylinder.c

-                      r61104c8
+                      r82592da
+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)
+{
 …
+}
 double
+static double
 Iq(double q, double radius_minor, double r_ratio, double length,
    double sld, double solvent_sld)
 …
         //const double arg = radius_minor*sin_val;
         double inner_sum=0;
         for(int j=0;j<20;j++) {
             //20 gauss points for the inner integral
             const double theta = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+        for(int j=0;j<76;j++) {
+            //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017
+            const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
             const double r = sin_val*sqrt(rA - rB*cos(theta));
             const double be = sas_2J1x_x(q*r);
             inner_sum += Gauss20Wt[j] * be * be;
+            inner_sum += Gauss76Wt[j] * be * be;
+        }
         //now calculate the value of the inner integral
 …
 double
 Iqxy(double qx, double qy,
+static double
+Iqxy(double qa, double qb, double qc,
      double radius_minor, double r_ratio, double length,
+     double sld, double solvent_sld,
+     double theta, double phi, double psi)
+     double sld, double solvent_sld)
+{
-    double q, xhat, yhat, zhat;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, 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*qb) + square(qa));
+    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/elliptical_cylinder.py

-                      rd9ec8f9
+                      reda8b30
 # pylint: disable=line-too-long
 r"""
-Definition for 2D (orientated system)
--------------------------------------
-The angles $\theta$ and $\phi$ define the orientation of the axis of the
-cylinder. The angle $\Psi$ is defined as the orientation of the major
-axis of the ellipse with respect to the vector $Q$. A gaussian polydispersity
-can be added to any of the orientation angles, and also for the minor
-radius and the ratio of the ellipse radii.
 .. figure:: img/elliptical_cylinder_geometry.png
 …
+Definition for 1D (no preferred orientation)
+--------------------------------------------
+The form factor is averaged over all possible orientation before normalized
+For 1D scattering, with no preferred orientation, the form factor is averaged over all possible orientations and normalized
 by the particle volume
 …
     P(q) = \text{scale}  <F^2> / V
+To provide easy access to the orientation of the elliptical cylinder, we
+define the axis of the cylinder using two angles $\theta$, $\phi$ and $\Psi$
+(see :ref:`cylinder orientation <cylinder-angle-definition>`). The angle
+$\Psi$ is the rotational angle around its own long_c axis.
+For 2d data the orientation of the particle is required, described using a different set
+of angles as in the diagrams below, for further details of the calculation and angular
+dispersions  see :ref:`orientation` .
-All angle parameters are valid and given only for 2D calculation; ie, an
-oriented system.
 .. figure:: img/elliptical_cylinder_angle_definition.png
+    Definition of angles for oriented elliptical cylinder, where axis_ratio is drawn >1,
+    and angle $\Psi$ is now a rotation around the axis of the cylinder.
+    Note that the angles here are not the same as in the equations for the scattering function.
+    Rotation $\theta$, initially in the $xz$ plane, is carried out first, then
+    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder.
+    The neutron or X-ray beam is along the $z$ axis.
 .. figure:: img/elliptical_cylinder_angle_projection.png
 …
 The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.
+On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
+appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the
+cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.)
+The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to
+understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation
+distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)
 NB: The 2nd virial coefficient of the cylinder is calculated based on the

sasmodels/models/fcc_paracrystal.c

-                      r50beefe
+                      rf728001
+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)
+{
+    // Equations from Matsuoka 17-18-19, multiplied by |q|
+    const double a1 = ( qa + qb)/2.0;
+    const double a2 = ( qa + qc)/2.0;
+    const double a3 = ( qb + qc)/2.0;
+double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
+double _FCCeval(double Theta, double Phi, double temp1, double temp3);
+    // Matsuoka 23-24-25
+    //     Z_k numerator: 1 - exp(a)^2
+    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a)
+    // Rewriting numerator
+    //         => -(exp(2a) - 1)
+    //         => -expm1(2a)
+    // Rewriting denominator
+    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => (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;
+}
+static double Iqxy(double qa, double qb, double qc,
+    double dnn, double d_factor, double radius,
+    double sld, double solvent_sld)
+{
+    const double 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;
+}
-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 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;
+}

sasmodels/models/fcc_paracrystal.py

-                      r8f04da4
+                      reda8b30
     \end{array}
+**NB**: The calculation of $Z(q)$ is a double numerical integral that
+must be carried out with a high density of points to properly capture
+the sharp peaks of the paracrystalline scattering. So be warned that the
+calculation is SLOW. Go get some coffee. Fitting of any experimental data
+must be resolution smeared for any meaningful fit. This makes a triple
+integral. Very, very slow. Go get lunch!
+.. note::
+  The calculation of $Z(q)$ is a double numerical integral that
+  must be carried out with a high density of points to properly capture
+  the sharp peaks of the paracrystalline scattering.
+  So be warned that the calculation is slow. Fitting of any experimental data
+  must be resolution smeared for any meaningful fit. This makes a triple integral
+  which may be very slow.
 The 2D (Anisotropic model) is based on the reference below where $I(q)$ is
 approximated for 1d scattering. Thus the scattering pattern for 2D may not
+be accurate. Note that we are not responsible for any incorrectness of the
+be accurate particularly at low $q$. For general details of the calculation
+and angular dispersions for oriented particles see :ref:`orientation` .
+Note that we are not responsible for any incorrectness of the
 D model computation.

sasmodels/models/hollow_cylinder.c

-                      r592343f
+                      rbecded3
-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 Iqxy(double qx, double qy, double radius, double thickness, double length, double sld,
-        double solvent_sld, double theta, double phi);
 //#define INVALID(v) (v.radius_core >= v.radius)
 …
+}
 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;
+}
 double
+static double
 form_volume(double radius, double thickness, double length)
+{
 …
 double
+static double
 Iq(double q, double radius, double th

Changeset ac60a39 in sasmodels

Legend: