Changeset 167d0f1 in sasmodels


Ignore:
Timestamp:
Nov 28, 2017 9:54:33 AM (5 years ago)
Author:
GitHub <noreply@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
110f69c
Parents:
fa70e04 (diff), 8c7d5d5 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Paul Kienzle <pkienzle@…> (11/28/17 09:54:33)
git-committer:
GitHub <noreply@…> (11/28/17 09:54:33)
Message:

Merge pull request #52 from SasView?/ticket-776-orientation

Ticket 776 orientation

Files:
16 added
1 deleted
53 edited

Legend:

Unmodified
Added
Removed
  • LICENSE.txt

    rf68e2a5 r5f8b72b  
    22All rights reserved. 
    33 
    4 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 
     4Redistribution and use in source and binary forms, with or without 
     5modification, are permitted provided that the following conditions are met: 
    56 
    6 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 
     71. Redistributions of source code must retain the above copyright notice, 
     8this list of conditions and the following disclaimer. 
    79 
    8 2. 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. 
     102. Redistributions in binary form must reproduce the above copyright notice, 
     11this list of conditions and the following disclaimer in the documentation 
     12and/or other materials provided with the distribution. 
    913 
    10 3. 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. 
     143. Neither the name of the copyright holder nor the names of its contributors 
     15may be used to endorse or promote products derived from this software without 
     16specific prior written permission. 
    1117 
    12 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. 
     18THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
     19AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
     20IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
     21ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 
     22LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
     23CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
     24SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
     25INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
     26CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
     27ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
     28POSSIBILITY OF SUCH DAMAGE. 
  • doc/_extensions/dollarmath.py

    r103ea45 rdd4d95d  
    1212import re 
    1313 
    14 _dollar = re.compile(r"(?:^|(?<=\s|[(]))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;)\\]))") 
     14_dollar = re.compile(r"(?:^|(?<=\s|[-(]))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[-.,;:?\\)]))") 
    1515_notdollar = re.compile(r"\\[$]") 
    1616 
  • doc/developer/calculator.rst

    r870a2f4 r2c108a3  
    77 
    88This document describes the layer between the form factor kernels and the 
    9 model calculator which implements the polydispersity and magnetic SLD 
     9model calculator which implements the dispersity and magnetic SLD 
    1010calculations.  There are three separate implementations of this layer, 
    1111:mod:`kernelcl` for OpenCL, which operates on a single Q value at a time, 
     
    1414 
    1515Each implementation provides three different calls *Iq*, *Iqxy* and *Imagnetic* 
    16 for 1-D, 2-D and 2-D magnetic kernels respectively. 
    17  
    18 The C code is defined in *kernel_iq.c* and *kernel_iq.cl* for DLL and OpenCL 
    19 respectively.  The kernel call looks as follows:: 
     16for 1-D, 2-D and 2-D magnetic kernels respectively. The C code is defined 
     17in *kernel_iq.c*, with the minor differences between OpenCL and DLL handled 
     18by #ifdef statements. 
     19 
     20The kernel call looks as follows:: 
    2021 
    2122  kernel void KERNEL_NAME( 
    2223      int nq,                  // Number of q values in the q vector 
    23       int pd_start,            // Starting position in the polydispersity loop 
    24       int pd_stop,             // Ending position in the polydispersity loop 
    25       ProblemDetails *details, // Polydispersity info 
     24      int pd_start,            // Starting position in the dispersity loop 
     25      int pd_stop,             // Ending position in the dispersity loop 
     26      ProblemDetails *details, // dispersity info 
    2627      double *values,          // Value and weights vector 
    2728      double *q,               // q or (qx,qy) vector 
    2829      double *result,          // returned I(q), with result[nq] = pd_weight 
    29       double cutoff)           // polydispersity weight cutoff 
     30      double cutoff)           // dispersity weight cutoff 
    3031 
    3132The details for OpenCL and the python loop are slightly different, but these 
     
    3435*nq* indicates the number of q values that will be calculated. 
    3536 
    36 The *pd_start* and *pd_stop* parameters set the range of the polydispersity 
    37 loop to compute for the current kernel call.   Give a polydispersity 
     37The *pd_start* and *pd_stop* parameters set the range of the dispersity 
     38loop to compute for the current kernel call.   Give a dispersity 
    3839calculation with 30 weights for length and 30 weights for radius for example, 
    3940there are a total of 900 calls to the form factor required to compute the 
     
    4243the length index to 3 and the radius index to 10 for a position of 3*30+10=100, 
    4344and could then proceed to position 200.  This allows us to interrupt the 
    44 calculation in the middle of a long polydispersity loop without having to 
     45calculation in the middle of a long dispersity loop without having to 
    4546do special tricks with the C code.  More importantly, it stops the OpenCL 
    4647kernel in a reasonable time; because the GPU is used by the operating 
     
    4950 
    5051The *ProblemDetails* structure is a direct map of the 
    51 :class:`details.CallDetails` buffer.  This indicates which parameters are 
    52 polydisperse, and where in the values vector the values and weights can be 
    53 found.  For each polydisperse parameter there is a parameter id, the length 
    54 of the polydispersity loop for that parameter, the offset of the parameter 
     52:class:`details.CallDetails` buffer.  This indicates which parameters have 
     53dispersity, and where in the values vector the values and weights can be 
     54found.  For each parameter with dispersity there is a parameter id, the length 
     55of the dispersity loop for that parameter, the offset of the parameter 
    5556values in the pd value and pd weight vectors and the 'stride' from one index 
    5657to the next, which is used to translate between the position in the 
    57 polydispersity loop and the particular parameter indices.  The *num_eval* 
    58 field is the total size of the polydispersity loop.  *num_weights* is the 
     58dispersity loop and the particular parameter indices.  The *num_eval* 
     59field is the total size of the dispersity loop.  *num_weights* is the 
    5960number of elements in the pd value and pd weight vectors.  *num_active* is 
    60 the number of non-trivial pd loops (polydisperse parameters should be ordered 
    61 by decreasing pd vector length, with a length of 1 meaning no polydispersity). 
     61the number of non-trivial pd loops (parameters with dispersity should be ordered 
     62by decreasing pd vector length, with a length of 1 meaning no dispersity). 
    6263Oriented objects in 2-D need a cos(theta) spherical correction on the angular 
    6364variation in order to preserve the 'surface area' of the weight distribution. 
     
    7273*(Mx, My, Mz)*.  Sample magnetization is translated from *(M, theta, phi)* 
    7374to *(Mx, My, Mz)* before the kernel is called.   After the fixed values comes 
    74 the pd value vector, with the polydispersity values for each parameter 
     75the pd value vector, with the dispersity values for each parameter 
    7576stacked one after the other.  The order isn't important since the location 
    7677for each parameter is stored in the *pd_offset* field of the *ProblemDetails* 
     
    7879values, the pd weight vector is stored, with the same configuration as the 
    7980pd value vector.  Note that the pd vectors can contain values that are not 
    80 in the polydispersity loop; this is used by :class:`mixture.MixtureKernel` 
     81in the dispersity loop; this is used by :class:`mixture.MixtureKernel` 
    8182to make it easier to call the various component kernels. 
    8283 
     
    8788 
    8889The *results* vector contains one slot for each of the *nq* values, plus 
    89 one extra slot at the end for the current polydisperse normalization.  This 
    90 is required when the polydispersity loop is broken across several kernel 
    91 calls. 
     90one extra slot at the end for the weight normalization accumulated across 
     91all points in the dispersity mesh.  This is required when the dispersity 
     92loop is broken across several kernel calls. 
    9293 
    9394*cutoff* is a importance cutoff so that points which contribute negligibly 
     
    9798 
    9899- USE_OPENCL is defined if running in opencl 
    99 - MAX_PD is the maximum depth of the polydispersity loop [model specific] 
     100- MAX_PD is the maximum depth of the dispersity loop [model specific] 
    100101- NUM_PARS is the number of parameter values in the kernel.  This may be 
    101102  more than the number of parameters if some of the parameters are vector 
    102103  values. 
    103104- NUM_VALUES is the number of fixed values, which defines the offset in the 
    104   value list to the polydisperse value and weight vectors. 
     105  value list to the dispersity value and weight vectors. 
    105106- NUM_MAGNETIC is the number of magnetic SLDs 
    106107- MAGNETIC_PARS is a comma separated list of the magnetic SLDs, indicating 
    107108  their locations in the values vector. 
    108 - MAGNETIC_PAR0 to MAGNETIC_PAR2 are the first three magnetic parameter ids 
    109   so we can hard code the setting of magnetic values if there are only a 
    110   few of them. 
    111109- KERNEL_NAME is the name of the function being declared 
    112110- PARAMETER_TABLE is the declaration of the parameters to the kernel: 
     
    152150    Cylinder2D:: 
    153151 
    154         #define CALL_IQ(q, i, var) Iqxy(q[2*i], q[2*i+1], \ 
     152        #define CALL_IQ(q, i, var) Iqxy(qa, qc, \ 
    155153        var.length, \ 
    156154        var.radius, \ 
    157155        var.sld, \ 
    158         var.sld_solvent, \ 
    159         var.theta, \ 
    160         var.phi) 
     156        var.sld_solvent) 
    161157 
    162158- CALL_VOLUME(var) is similar, but for calling the form volume:: 
     
    182178        #define INVALID(var) constrained(var.p1, var.p2, var.p3) 
    183179 
    184 Our design supports a limited number of polydispersity loops, wherein 
    185 we need to cycle through the values of the polydispersity, calculate 
     180Our design supports a limited number of dispersity loops, wherein 
     181we need to cycle through the values of the dispersity, calculate 
    186182the I(q, p) for each combination of parameters, and perform a normalized 
    187183weighted sum across all the weights.  Parameters may be passed to the 
    188 underlying calculation engine as scalars or vectors, but the polydispersity 
     184underlying calculation engine as scalars or vectors, but the dispersity 
    189185calculator treats the parameter set as one long vector. 
    190186 
    191 Let's assume we have 8 parameters in the model, with two polydisperse.  Since 
    192 this is a 1-D model the orientation parameters won't be used:: 
     187Let's assume we have 8 parameters in the model, two of which allow dispersity. 
     188Since this is a 1-D model the orientation parameters won't be used:: 
    193189 
    194190    0: scale        {scl = constant} 
     
    196192    2: radius       {r = vector of 10pts} 
    197193    3: length       {l = vector of 30pts} 
    198     4: sld          {s = constant/(radius**2*length)} 
     194    4: sld          {s1 = constant/(radius**2*length)} 
    199195    5: sld_solvent  {s2 = constant} 
    200196    6: theta        {not used} 
     
    202198 
    203199This generates the following call to the kernel.  Note that parameters 4 and 
    204 5 are treated as polydisperse even though they are not --- this is because 
     2005 are treated as having dispersity even though they don't --- this is because 
    205201it is harmless to do so and simplifies the looping code:: 
    206202 
     
    210206    NUM_MAGNETIC = 2      // two parameters might be magnetic 
    211207    MAGNETIC_PARS = 4, 5  // they are sld and sld_solvent 
    212     MAGNETIC_PAR0 = 4     // sld index 
    213     MAGNETIC_PAR1 = 5     // solvent index 
    214208 
    215209    details { 
     
    218212        pd_offset = {10, 0, 31, 32}   // *length* starts at index 10 in weights 
    219213        pd_stride = {1, 30, 300, 300} // cumulative product of pd length 
    220         num_eval = 300   // 300 values in the polydispersity loop 
     214        num_eval = 300   // 300 values in the dispersity loop 
    221215        num_weights = 42 // 42 values in the pd vector 
    222216        num_active = 2   // only the first two pd are active 
     
    225219 
    226220    values = { scl, bkg,                                  // universal 
    227                r, l, s, s2, theta, phi,                   // kernel pars 
     221               r, l, s1, s2, theta, phi,                  // kernel pars 
    228222               in spin, out spin, spin angle,             // applied magnetism 
    229                mx s, my s, mz s, mx s2, my s2, mz s2,     // magnetic slds 
     223               mx s1, my s1, mz s1, mx s2, my s2, mz s2,  // magnetic slds 
    230224               r0, .., r9, l0, .., l29, s, s2,            // pd values 
    231225               r0, .., r9, l0, .., l29, s, s2}            // pd weights 
     
    235229    result = {r1, ..., r130, pd_norm, x } 
    236230 
    237 The polydisperse parameters are stored in as an array of parameter 
    238 indices, one for each polydisperse parameter, stored in pd_par[n]. 
    239 Non-polydisperse parameters do not appear in this array. Each polydisperse 
     231The dispersity parameters are stored in as an array of parameter 
     232indices, one for each parameter, stored in pd_par[n]. Parameters which do 
     233not support dispersity do not appear in this array. Each dispersity 
    240234parameter has a weight vector whose length is stored in pd_length[n]. 
    241235The weights are stored in a contiguous vector of weights for all 
     
    243237in pd_offset[n].  The values corresponding to the weights are stored 
    244238together in a separate weights[] vector, with offset stored in 
    245 par_offset[pd_par[n]]. Polydisperse parameters should be stored in 
     239par_offset[pd_par[n]]. Dispersity parameters should be stored in 
    246240decreasing order of length for highest efficiency. 
    247241 
    248 We limit the number of polydisperse dimensions to MAX_PD (currently 4), 
    249 though some models may have fewer if they have fewer polydisperse 
     242We limit the number of dispersity dimensions to MAX_PD (currently 4), 
     243though some models may have fewer if they have fewer dispersity 
    250244parameters.  The main reason for the limit is to reduce code size. 
    251 Each additional polydisperse parameter requires a separate polydispersity 
    252 loop.  If more than 4 levels of polydispersity are needed, then kernel_iq.c 
    253 and kernel_iq.cl will need to be extended. 
     245Each additional dispersity parameter requires a separate dispersity 
     246loop.  If more than 4 levels of dispersity are needed, then we need to 
     247switch to a monte carlo importance sampling algorithm with better 
     248performance for high-dimensional integrals. 
    254249 
    255250Constraints between parameters are not supported.  Instead users will 
     
    262257theta since the polar coordinates normalization is tied to this parameter. 
    263258 
    264 If there is no polydispersity we pretend that it is polydisperisty with one 
    265 parameter, pd_start=0 and pd_stop=1.  We may or may not short circuit the 
    266 calculation in this case, depending on how much time it saves. 
     259If there is no dispersity we pretend that we have a disperisty mesh over 
     260a single parameter with a single point in the distribution, giving 
     261pd_start=0 and pd_stop=1. 
    267262 
    268263The problem details structure could be allocated and sent in as an integer 
    269264array using the read-only flag.  This would allow us to copy it once per fit 
    270265along with the weights vector, since features such as the number of 
    271 polydisperity elements per pd parameter or the coordinated won't change 
    272 between function evaluations.  A new parameter vector must be sent for 
    273 each I(q) evaluation.  This is not currently implemented, and would require 
    274 some resturcturing of the :class:`sasview_model.SasviewModel` interface. 
    275  
    276 The results array will be initialized to zero for polydispersity loop 
     266disperity points per pd parameter won't change between function evaluations. 
     267A new parameter vector must be sent for each I(q) evaluation.  This is 
     268not currently implemented, and would require some resturcturing of 
     269the :class:`sasview_model.SasviewModel` interface. 
     270 
     271The results array will be initialized to zero for dispersity loop 
    277272entry zero, and preserved between calls to [start, stop] so that the 
    278273results accumulate by the time the loop has completed.  Background and 
     
    295290 
    296291This will require accumulated error for each I(q) value to be preserved 
    297 between kernel calls to implement this fully.  The kernel_iq.c code, which 
    298 loops over q for each parameter set in the polydispersity loop, will need 
    299 also need the accumalation vector. 
     292between kernel calls to implement this fully.  The *kernel_iq.c* code, which 
     293loops over q for each parameter set in the dispersity loop, will also need 
     294the accumulation vector. 
  • doc/developer/overview.rst

    r870a2f4 r3d40839  
    166166:ref:`Calculator_Interface` 
    167167 
     168.. _orientation_developer: 
     169 
     170Orientation and Numerical Integration 
     171------------------------------------- 
     172 
     173For 2d data from oriented anisotropic particles, the mean particle 
     174orientation is defined by angles $\theta$, $\phi$ and $\Psi$, which are not 
     175in general the same as similarly named angles in many form factors. The 
     176wikipedia page on Euler angles (https://en.wikipedia.org/wiki/Euler_angles) 
     177lists the different conventions available. To quote: "Different authors may 
     178use different sets of rotation axes to define Euler angles, or different 
     179names for the same angles. Therefore, any discussion employing Euler angles 
     180should always be preceded by their definition." 
     181 
     182We are using the $z$-$y$-$z$ convention with extrinsic rotations 
     183$\Psi$-$\theta$-$\phi$ for the particle orientation and $x$-$y$-$z$ 
     184convention with extrinsic rotations $\Psi$-$\theta$-$\phi$ for jitter, with 
     185jitter applied before particle orientation. 
     186 
     187For numerical integration within form factors etc. sasmodels is mostly using 
     188Gaussian quadrature with 20, 76 or 150 points depending on the model. It also 
     189makes use of symmetries such as calculating only over one quadrant rather 
     190than the whole sphere. There is often a U-substitution replacing $\theta$ 
     191with $cos(\theta)$ which changes the limits of integration from 0 to $\pi/2$ 
     192to 0 to 1 and also conveniently absorbs the $sin(\theta)$ scale factor in the 
     193integration. This can cause confusion if checking equations to include in a 
     194paper or thesis! Most models use the same core kernel code expressed in terms 
     195of the rotated view ($q_a$, $q_b$, $q_c$) for both the 1D and the 2D models, 
     196but there are also historical quirks such as the parallelepiped model, which 
     197has a useless transformation representing $j_0(a q_a)$ as $j_0(b q_a a/b)$. 
     198 
     199Useful testing routines include: 
     200 
     201:mod:`asymint` a direct implementation of the surface integral for certain 
     202models to get a more trusted value for the 1D integral using a 
     203reimplementation of the 2D kernel in python and mpmath (which computes math 
     204functions to arbitrary precision). It uses $\theta$ ranging from 0 to $\pi$ 
     205and $\phi$ ranging from 0 to $2\pi$. It perhaps would benefit from including 
     206the U-substitution for $\theta$. 
     207 
     208:mod:`check1d` uses sasmodels 1D integration and compares that with a 
     209rectangle distribution in $\theta$ and $\phi$, with $\theta$ limits set to 
     210$\pm 90/\sqrt(3)$ and $\phi$ limits set to $\pm 180/\sqrt(3)$ [The rectangle 
     211weight function uses the fact that the distribution width column is labelled 
     212sigma to decide that the 1-$\sigma$ width of a rectangular distribution needs to 
     213be multiplied by $\sqrt(3)$ to get the corresponding gaussian equivalent, or 
     214similar reasoning.] This should rotate the sample through the entire 
     215$\theta$-$\phi$ surface according to the pattern that you see in jitter.py when 
     216you modify it to use 'rectangle' rather than 'gaussian' for its distribution 
     217without changing the viewing angle. In order to match the 1-D pattern for 
     218an arbitrary viewing angle on triaxial shapes, we need to integrate 
     219over $\theta$, $\phi$ and $\Psi$. 
     220 
     221When computing the dispersity integral, weights are scaled by 
     222$|\cos(\delta \theta)|$ to account for the points in $\phi$ getting closer 
     223together as $\delta \theta$ increases. 
     224[This will probably change so that instead of adjusting the weights, we will 
     225adjust $\delta\theta$-$\delta\phi$ mesh so that the point density in 
     226$\delta\phi$ is lower at larger $\delta\theta$. The flag USE_SCALED_PHI in 
     227*kernel_iq.c* selects an alternative algorithm.] 
     228 
     229The integrated dispersion is computed at a set of $(qx, qy)$ points $(q 
     230\cos(\alpha), q \sin(\alpha))$ at some angle $\alpha$ (currently angle=0) for 
     231each $q$ used in the 1-D integration. The individual $q$ points should be 
     232equivalent to asymint rect-n when the viewing angle is set to 
     233$(\theta,\phi,\Psi) = (90, 0, 0)$. Such tests can help to validate that 2d 
     234models are consistent with 1d models. 
     235 
     236:mod:`sascomp -sphere=n` uses the same rectangular distribution as check1d to 
     237compute the pattern of the $q_x$-$q_y$ grid. 
     238 
     239The :mod:`sascomp` utility can be used for 2d as well as 1d calculations to 
     240compare results for two sets of parameters or processor types, for example 
     241these two oriented cylinders here should be equivalent. 
     242 
     243: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` 
     244 
    168245 
    169246Testing 
  • doc/genapi.py

    r2e66ef5 r706f466  
    5959    #('alignment', 'GPU data alignment [unused]'), 
    6060    ('bumps_model', 'Bumps interface'), 
     61    ('compare', 'Compare models on different compute engines'), 
    6162    ('compare_many', 'Batch compare models on different compute engines'), 
    62     ('compare', 'Compare models on different compute engines'), 
     63    ('conversion_table', 'Model conversion table'), 
    6364    ('convert', 'Sasview to sasmodel converter'), 
    6465    ('core', 'Model access'), 
     
    8283    ('sasview_model', 'Sasview interface'), 
    8384    ('sesans', 'SESANS calculation routines'), 
     85    ('special', 'Special functions library'), 
    8486    ('weights', 'Distribution functions'), 
    8587] 
  • doc/guide/index.rst

    rc0d7ab3 rda5536f  
    1313   resolution.rst 
    1414   magnetism/magnetism.rst 
     15   orientation/orientation.rst 
    1516   sesans/sans_to_sesans.rst 
    1617   sesans/sesans_fitting.rst 
  • doc/guide/magnetism/magnetism.rst

    r1f058ea r4f5afc9  
    55 
    66Models which define a scattering length density parameter can be evaluated 
    7  as magnetic models. In general, the scattering length density (SLD = 
    8  $\beta$) in each region where the SLD is uniform, is a combination of the 
    9  nuclear and magnetic SLDs and, for polarised neutrons, also depends on the 
    10  spin states of the neutrons. 
     7as magnetic models. In general, the scattering length density (SLD = 
     8$\beta$) in each region where the SLD is uniform, is a combination of the 
     9nuclear and magnetic SLDs and, for polarised neutrons, also depends on the 
     10spin states of the neutrons. 
    1111 
    1212For magnetic scattering, only the magnetization component $\mathbf{M_\perp}$ 
    1313perpendicular to the scattering vector $\mathbf{Q}$ contributes to the magnetic 
    1414scattering length. 
    15  
    1615 
    1716.. figure:: 
     
    2827is the Pauli spin. 
    2928 
    30 Assuming that incident neutrons are polarized parallel (+) and anti-parallel (-) 
    31 to the $x'$ axis, the possible spin states after the sample are then 
     29Assuming that incident neutrons are polarized parallel $(+)$ and anti-parallel 
     30$(-)$ to the $x'$ axis, the possible spin states after the sample are then: 
    3231 
    33 No spin-flips (+ +) and (- -) 
     32* Non spin-flip $(+ +)$ and $(- -)$ 
    3433 
    35 Spin-flips    (+ -) and (- +) 
     34* Spin-flip $(+ -)$ and $(- +)$ 
     35 
     36Each measurement is an incoherent mixture of these spin states based on the 
     37fraction of $+$ neutrons before ($u_i$) and after ($u_f$) the sample, 
     38with weighting: 
     39 
     40.. math:: 
     41    -- &= ((1-u_i)(1-u_f))^{1/4} \\ 
     42    -+ &= ((1-u_i)(u_f))^{1/4} \\ 
     43    +- &= ((u_i)(1-u_f))^{1/4} \\ 
     44    ++ &= ((u_i)(u_f))^{1/4} 
     45 
     46Ideally the experiment would measure the pure spin states independently and 
     47perform a simultaneous analysis of the four states, tying all the model 
     48parameters together except $u_i$ and $u_f$. 
    3649 
    3750.. figure:: 
     
    4154$\phi$ and $\theta_{up}$, respectively, then, depending on the spin state of the 
    4255neutrons, the scattering length densities, including the nuclear scattering 
    43 length density ($\beta{_N}$) are 
     56length density $(\beta{_N})$ are 
    4457 
    4558.. math:: 
    4659    \beta_{\pm\pm} =  \beta_N \mp D_M M_{\perp x'} 
    47     \text{ when there are no spin-flips} 
     60    \text{ for non spin-flip states} 
    4861 
    4962and 
     
    5164.. math:: 
    5265    \beta_{\pm\mp} =  -D_M (M_{\perp y'} \pm iM_{\perp z'}) 
    53     \text{ when there are} 
     66    \text{ for spin-flip states} 
    5467 
    5568where 
    5669 
    5770.. math:: 
    58     M_{\perp x'} = M_{0q_x}\cos(\theta_{up})+M_{0q_y}\sin(\theta_{up}) \\ 
    59     M_{\perp y'} = M_{0q_y}\cos(\theta_{up})-M_{0q_x}\sin(\theta_{up}) \\ 
    60     M_{\perp z'} = M_{0z} \\ 
    61     M_{0q_x} = (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\ 
    62     M_{0q_y} = (M_{0y}\sin\phi - M_{0x}\cos\phi)\sin\phi 
     71    M_{\perp x'} &= M_{0q_x}\cos(\theta_{up})+M_{0q_y}\sin(\theta_{up}) \\ 
     72    M_{\perp y'} &= M_{0q_y}\cos(\theta_{up})-M_{0q_x}\sin(\theta_{up}) \\ 
     73    M_{\perp z'} &= M_{0z} \\ 
     74    M_{0q_x} &= (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\ 
     75    M_{0q_y} &= (M_{0y}\sin\phi - M_{0x}\cos\phi)\sin\phi 
    6376 
    6477Here, $M_{0x}$, $M_{0x}$, $M_{0z}$ are the x, y and z components 
     
    6679 
    6780.. math:: 
    68     M_{0x} = M_0\cos\theta_M\cos\phi_M \\ 
    69     M_{0y} = M_0\sin\theta_M \\ 
    70     M_{0z} = -M_0\cos\theta_M\sin\phi_M 
     81    M_{0x} &= M_0\cos\theta_M\cos\phi_M \\ 
     82    M_{0y} &= M_0\sin\theta_M \\ 
     83    M_{0z} &= -M_0\cos\theta_M\sin\phi_M 
    7184 
    7285and the magnetization angles $\theta_M$ and $\phi_M$ are defined in 
     
    7689 
    7790===========   ================================================================ 
    78  M0_sld        = $D_M M_0$ 
    79  Up_theta      = $\theta_\mathrm{up}$ 
    80  M_theta       = $\theta_M$ 
    81  M_phi         = $\phi_M$ 
    82  Up_frac_i     = (spin up)/(spin up + spin down) neutrons *before* the sample 
    83  Up_frac_f     = (spin up)/(spin up + spin down) neutrons *after* the sample 
     91 M0:sld      $D_M M_0$ 
     92 mtheta:sld   $\theta_M$ 
     93 mphi:sld     $\phi_M$ 
     94 up:angle     $\theta_\mathrm{up}$ 
     95 up:frac_i    $u_i$ = (spin up)/(spin up + spin down) *before* the sample 
     96 up:frac_f    $u_f$ = (spin up)/(spin up + spin down) *after* the sample 
    8497===========   ================================================================ 
    8598 
    8699.. note:: 
    87     The values of the 'Up_frac_i' and 'Up_frac_f' must be in the range 0 to 1. 
     100    The values of the 'up:frac_i' and 'up:frac_f' must be in the range 0 to 1. 
    88101 
    89102*Document History* 
    90103 
    91104| 2015-05-02 Steve King 
    92 | 2017-05-08 Paul Kienzle 
     105| 2017-11-15 Paul Kienzle 
  • doc/guide/pd/polydispersity.rst

    r1f058ea reda8b30  
    66.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
    77 
     8.. _polydispersityhelp: 
     9 
    810Polydispersity Distributions 
    911---------------------------- 
    1012 
    11 With some models in sasmodels we can calculate the average form factor for a 
     13With some models in sasmodels we can calculate the average intensity for a 
    1214population of particles that exhibit size and/or orientational 
    13 polydispersity. The resultant form factor is normalized by the average 
     15polydispersity. The resultant intensity is normalized by the average 
    1416particle volume such that 
    1517 
  • doc/guide/resolution.rst

    r1f058ea r0db85af  
    209209$x'_0 = x_0 \cos(\theta) + y_0 \sin(\theta)$ and 
    210210$y'_0 = -x_0 \sin(\theta) + y_0 \cos(\theta)$. 
    211 Note that the rotation angle is zero for a $x$\ -\ $y$ symmetric 
     211Note that the rotation angle is zero for a $x$-$y$ symmetric 
    212212elliptical Gaussian distribution. The $A$ is a normalization factor. 
    213213 
     
    233233 
    234234Since the weighting factor on each of the bins is known, it is convenient to 
    235 transform $x'$\ -\ $y'$ back to $x$\ -\ $y$ coordinates (by rotating it 
     235transform $x'$-$y'$ back to $x$-$y$ coordinates (by rotating it 
    236236by $-\theta$ around the $z$\ -axis). 
    237237 
     
    254254    y'_0 &= 0 
    255255 
    256 while for a $x$\ -\ $y$ symmetric smear 
     256while for a $x$-$y$ symmetric smear 
    257257 
    258258.. math:: 
  • explore/jitter.py

    • Property mode changed from 100644 to 100755
    r85190c2 rff10479  
     1#!/usr/bin/env python 
    12""" 
    23Application to explore the difference between sasview 3.x orientation 
    34dispersity and possible replacement algorithms. 
    45""" 
     6from __future__ import division, print_function 
     7 
     8import sys, os 
     9sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
     10 
     11import argparse 
     12 
    513import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
    614import matplotlib.pyplot as plt 
    715from matplotlib.widgets import Slider, CheckButtons 
    816from matplotlib import cm 
    9  
    1017import numpy as np 
    1118from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1219 
    13 def draw_beam(ax): 
     20def draw_beam(ax, view=(0, 0)): 
     21    """ 
     22    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 
     23    """ 
    1424    #ax.plot([0,0],[0,0],[1,-1]) 
    1525    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
     
    2232    x = r*np.outer(np.cos(u), np.ones_like(v)) 
    2333    y = r*np.outer(np.sin(u), np.ones_like(v)) 
    24     z = np.outer(np.ones_like(u), v) 
     34    z = 1.3*np.outer(np.ones_like(u), v) 
     35 
     36    theta, phi = view 
     37    shape = x.shape 
     38    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
     39    points = Rz(phi)*Ry(theta)*points 
     40    x, y, z = [v.reshape(shape) for v in points] 
    2541 
    2642    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    2743 
    28 def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi): 
    29     size=[0.1, 0.4, 1.0] 
    30     view=[theta, phi, psi] 
    31     shimmy=[0,0,0] 
    32     #draw_shape = draw_parallelepiped 
    33     draw_shape = draw_ellipsoid 
     44def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0)): 
     45    """ 
     46    Represent jitter as a set of shapes at different orientations. 
     47    """ 
     48    # set max diagonal to 0.95 
     49    scale = 0.95/sqrt(sum(v**2 for v in size)) 
     50    size = tuple(scale*v for v in size) 
     51    draw_shape = draw_parallelepiped 
     52    #draw_shape = draw_ellipsoid 
    3453 
    3554    #np.random.seed(10) 
     
    6483        [ 1,  1,  1], 
    6584    ] 
     85    dtheta, dphi, dpsi = jitter 
    6686    if dtheta == 0: 
    6787        cloud = [v for v in cloud if v[0] == 0] 
     
    7090    if dpsi == 0: 
    7191        cloud = [v for v in cloud if v[2] == 0] 
    72     draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8) 
     92    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
     93    scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] 
    7394    for point in cloud: 
    74         shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]] 
    75         draw_shape(ax, size, view, shimmy, alpha=0.8) 
     95        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
     96        draw_shape(ax, size, view, delta, alpha=0.8) 
    7697    for v in 'xyz': 
    7798        a, b, c = size 
     
    80101        getattr(ax, v+'axis').label.set_text(v) 
    81102 
    82 def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1): 
     103def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
     104    """Draw an ellipsoid.""" 
    83105    a,b,c = size 
    84     theta, phi, psi = view 
    85     dtheta, dphi, dpsi = shimmy 
    86  
    87106    u = np.linspace(0, 2 * np.pi, steps) 
    88107    v = np.linspace(0, np.pi, steps) 
     
    90109    y = b*np.outer(np.sin(u), np.sin(v)) 
    91110    z = c*np.outer(np.ones_like(u), np.cos(v)) 
    92  
    93     shape = x.shape 
    94     points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
    95     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    96     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    97     x,y,z = [v.reshape(shape) for v in points] 
     111    x, y, z = transform_xyz(view, jitter, x, y, z) 
    98112 
    99113    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    100114 
    101 def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 
     115    draw_labels(ax, view, jitter, [ 
     116         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     117         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     118         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     119         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     120         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     121         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     122    ]) 
     123 
     124def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
     125    """Draw a parallelepiped.""" 
    102126    a,b,c = size 
    103     theta, phi, psi = view 
    104     dtheta, dphi, dpsi = shimmy 
    105  
    106127    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    107128    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
     
    118139    ]) 
    119140 
    120     points = np.matrix([x,y,z]) 
    121     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    122     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    123  
    124     x,y,z = [np.array(v).flatten() for v in points] 
     141    x, y, z = transform_xyz(view, jitter, x, y, z) 
    125142    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    126143 
     144    draw_labels(ax, view, jitter, [ 
     145         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     146         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     147         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     148         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     149         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     150         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     151    ]) 
     152 
    127153def draw_sphere(ax, radius=10., steps=100): 
     154    """Draw a sphere""" 
    128155    u = np.linspace(0, 2 * np.pi, steps) 
    129156    v = np.linspace(0, np.pi, steps) 
     
    134161    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    135162 
    136 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    137     theta_center = radians(theta) 
    138     phi_center = radians(phi) 
    139     flow_center = radians(flow) 
    140     dtheta = radians(dtheta) 
    141     dphi = radians(dphi) 
    142  
    143     # 10 point 3-sigma gaussian weights 
    144     t = np.linspace(-3., 3., 11) 
    145     if dist == 'gauss': 
     163PROJECTIONS = [ 
     164    # in order of PROJECTION number; do not change without updating the 
     165    # constants in kernel_iq.c 
     166    'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance', 
     167] 
     168def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', 
     169              projection='equirectangular'): 
     170    """ 
     171    Draw the dispersion mesh showing the theta-phi orientations at which 
     172    the model will be evaluated. 
     173 
     174    jitter projections 
     175    <https://en.wikipedia.org/wiki/List_of_map_projections> 
     176 
     177    equirectangular (standard latitude-longitude mesh) 
     178        <https://en.wikipedia.org/wiki/Equirectangular_projection> 
     179        Allows free movement in phi (around the equator), but theta is 
     180        limited to +/- 90, and points are cos-weighted. Jitter in phi is 
     181        uniform in weight along a line of latitude.  With small theta and 
     182        phi ranging over +/- 180 this forms a wobbling disk.  With small 
     183        phi and theta ranging over +/- 90 this forms a wedge like a slice 
     184        of an orange. 
     185    azimuthal_equidistance (Postel) 
     186        <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection> 
     187        Preserves distance from center, and so is an excellent map for 
     188        representing a bivariate gaussian on the surface.  Theta and phi 
     189        operate identically, cutting wegdes from the antipode of the viewing 
     190        angle.  This unfortunately does not allow free movement in either 
     191        theta or phi since the orthogonal wobble decreases to 0 as the body 
     192        rotates through 180 degrees. 
     193    sinusoidal (Sanson-Flamsteed, Mercator equal-area) 
     194        <https://en.wikipedia.org/wiki/Sinusoidal_projection> 
     195        Preserves arc length with latitude, giving bad behaviour at 
     196        theta near +/- 90.  Theta and phi operate somewhat differently, 
     197        so a system with a-b-c dtheta-dphi-dpsi will not give the same 
     198        value as one with b-a-c dphi-dtheta-dpsi, as would be the case 
     199        for azimuthal equidistance.  Free movement using theta or phi 
     200        uniform over +/- 180 will work, but not as well as equirectangular 
     201        phi, with theta being slightly worse.  Computationally it is much 
     202        cheaper for wide theta-phi meshes since it excludes points which 
     203        lie outside the sinusoid near theta +/- 90 rather than packing 
     204        them close together as in equirectangle.  Note that the poles 
     205        will be slightly overweighted for theta > 90 with the circle 
     206        from theta at 90+dt winding backwards around the pole, overlapping 
     207        the circle from theta at 90-dt. 
     208    Guyou (hemisphere-in-a-square) **not weighted** 
     209        <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection> 
     210        With tiling, allows rotation in phi or theta through +/- 180, with 
     211        uniform spacing.  Both theta and phi allow free rotation, with wobble 
     212        in the orthogonal direction reasonably well behaved (though not as 
     213        good as equirectangular phi). The forward/reverse transformations 
     214        relies on elliptic integrals that are somewhat expensive, so the 
     215        behaviour has to be very good to justify the cost and complexity. 
     216        The weighting function for each point has not yet been computed. 
     217        Note: run the module *guyou.py* directly and it will show the forward 
     218        and reverse mappings. 
     219    azimuthal_equal_area  **incomplete** 
     220        <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection> 
     221        Preserves the relative density of the surface patches.  Not that 
     222        useful and not completely implemented 
     223    Gauss-Kreuger **not implemented** 
     224        <https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Ellipsoidal_transverse_Mercator> 
     225        Should allow free movement in theta, but phi is distorted. 
     226    """ 
     227    theta, phi, psi = view 
     228    dtheta, dphi, dpsi = jitter 
     229 
     230    t = np.linspace(-1, 1, n) 
     231    weights = np.ones_like(t) 
     232    if dist == 'gaussian': 
     233        t *= 3 
    146234        weights = exp(-0.5*t**2) 
    147     elif dist == 'rect': 
    148         weights = np.ones_like(t) 
     235    elif dist == 'rectangle': 
     236        # Note: uses sasmodels ridiculous definition of rectangle width 
     237        t *= sqrt(3) 
     238    elif dist == 'uniform': 
     239        pass 
    149240    else: 
    150         raise ValueError("expected dist to be 'gauss' or 'rect'") 
    151     theta = dtheta*t 
    152     phi = dphi*t 
    153  
    154     x = radius * np.outer(cos(phi), cos(theta)) 
    155     y = radius * np.outer(sin(phi), cos(theta)) 
    156     z = radius * np.outer(np.ones_like(phi), sin(theta)) 
    157     #w = np.outer(weights, weights*abs(cos(dtheta*t))) 
    158     w = np.outer(weights, weights*abs(cos(theta))) 
    159  
    160     x, y, z, w = [v.flatten() for v in (x,y,z,w)] 
    161     x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center) 
    162  
     241        raise ValueError("expected dist to be gaussian, rectangle or uniform") 
     242 
     243    if projection == 'equirectangular':  #define PROJECTION 1 
     244        def rotate(theta_i, phi_j): 
     245            return Rx(phi_j)*Ry(theta_i) 
     246        def weight(theta_i, phi_j, wi, wj): 
     247            return wi*wj*abs(cos(radians(theta_i))) 
     248    elif projection == 'sinusoidal':  #define PROJECTION 2 
     249        def rotate(theta_i, phi_j): 
     250            latitude = theta_i 
     251            scale = cos(radians(latitude)) 
     252            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0 
     253            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
     254            return Rx(longitude)*Ry(latitude) 
     255        def weight(theta_i, phi_j, wi, wj): 
     256            latitude = theta_i 
     257            scale = cos(radians(latitude)) 
     258            w = 1 if abs(phi_j) < abs(scale)*180 else 0 
     259            return w*wi*wj 
     260    elif projection == 'guyou':  #define PROJECTION 3  (eventually?) 
     261        def rotate(theta_i, phi_j): 
     262            from guyou import guyou_invert 
     263            #latitude, longitude = guyou_invert([theta_i], [phi_j]) 
     264            longitude, latitude = guyou_invert([phi_j], [theta_i]) 
     265            return Rx(longitude[0])*Ry(latitude[0]) 
     266        def weight(theta_i, phi_j, wi, wj): 
     267            return wi*wj 
     268    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry 
     269        def rotate(theta_i, phi_j): 
     270            latitude = sqrt(theta_i**2 + phi_j**2) 
     271            longitude = degrees(np.arctan2(phi_j, theta_i)) 
     272            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
     273            return Rz(longitude)*Ry(latitude) 
     274        def weight(theta_i, phi_j, wi, wj): 
     275            # Weighting for each point comes from the integral: 
     276            #     \int\int I(q, lat, log) sin(lat) dlat dlog 
     277            # We are doing a conformal mapping from disk to sphere, so we need 
     278            # a change of variables g(theta, phi) -> (lat, long): 
     279            #     lat, long = sqrt(theta^2 + phi^2), arctan(phi/theta) 
     280            # giving: 
     281            #     dtheta dphi = det(J) dlat dlong 
     282            # where J is the jacobian from the partials of g. Using 
     283            #     R = sqrt(theta^2 + phi^2), 
     284            # then 
     285            #     J = [[x/R, Y/R], -y/R^2, x/R^2]] 
     286            # and 
     287            #     det(J) = 1/R 
     288            # with the final integral being: 
     289            #    \int\int I(q, theta, phi) sin(R)/R dtheta dphi 
     290            # 
     291            # This does approximately the right thing, decreasing the weight 
     292            # of each point as you go farther out on the disk, but it hasn't 
     293            # yet been checked against the 1D integral results. Prior 
     294            # to declaring this "good enough" and checking that integrals 
     295            # work in practice, we will examine alternative mappings. 
     296            # 
     297            # The issue is that the mapping does not support the case of free 
     298            # rotation about a single axis correctly, with a small deviation 
     299            # in the orthogonal axis independent of the first axis.  Like the 
     300            # usual polar coordiates integration, the integrated sections 
     301            # form wedges, though at least in this case the wedge cuts through 
     302            # the entire sphere, and treats theta and phi identically. 
     303            latitude = sqrt(theta_i**2 + phi_j**2) 
     304            w = sin(radians(latitude))/latitude if latitude != 0 else 1 
     305            return w*wi*wj if latitude < 180 else 0 
     306    elif projection == 'azimuthal_equal_area': 
     307        def rotate(theta_i, phi_j): 
     308            R = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
     309            latitude = 180-degrees(2*np.arccos(R)) 
     310            longitude = degrees(np.arctan2(phi_j, theta_i)) 
     311            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
     312            return Rz(longitude)*Ry(latitude) 
     313        def weight(theta_i, phi_j, wi, wj): 
     314            latitude = sqrt(theta_i**2 + phi_j**2) 
     315            w = sin(radians(latitude))/latitude if latitude != 0 else 1 
     316            return w*wi*wj if latitude < 180 else 0 
     317    else: 
     318        raise ValueError("unknown projection %r"%projection) 
     319 
     320    # mesh in theta, phi formed by rotating z 
     321    z = np.matrix([[0], [0], [radius]]) 
     322    points = np.hstack([rotate(theta_i, phi_j)*z 
     323                        for theta_i in dtheta*t 
     324                        for phi_j in dphi*t]) 
     325    # select just the active points (i.e., those with phi < 180 
     326    w = np.array([weight(theta_i, phi_j, wi, wj) 
     327                  for wi, theta_i in zip(weights, dtheta*t) 
     328                  for wj, phi_j in zip(weights, dphi*t)]) 
     329    #print(max(w), min(w), min(w[w>0])) 
     330    points = points[:, w>0] 
     331    w = w[w>0] 
     332    w /= max(w) 
     333 
     334    if 0: # Kent distribution 
     335        points = np.hstack([Rx(phi_j)*Ry(theta_i)*z for theta_i in 30*t for phi_j in 60*t]) 
     336        xp, yp, zp = [np.array(v).flatten() for v in points] 
     337        kappa = max(1e6, radians(dtheta)/(2*pi)) 
     338        beta = 1/max(1e-6, radians(dphi)/(2*pi))/kappa 
     339        w = exp(kappa*zp) #+ beta*(xp**2 + yp**2) 
     340        print(kappa, dtheta, radians(dtheta), min(w), max(w), sum(w)) 
     341        #w /= abs(cos(radians( 
     342        #w /= sum(w) 
     343 
     344    # rotate relative to beam 
     345    points = orient_relative_to_beam(view, points) 
     346 
     347    x, y, z = [np.array(v).flatten() for v in points] 
     348    #plt.figure(2); plt.clf(); plt.hist(z, bins=np.linspace(-1, 1, 51)) 
    163349    ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
    164350 
    165 def rotate(x, y, z, phi, theta, psi): 
    166     R = Rz(phi)*Ry(theta)*Rz(psi) 
    167     p = np.vstack([x,y,z]) 
    168     return R*p 
    169  
     351def draw_labels(ax, view, jitter, text): 
     352    """ 
     353    Draw text at a particular location. 
     354    """ 
     355    labels, locations, orientations = zip(*text) 
     356    px, py, pz = zip(*locations) 
     357    dx, dy, dz = zip(*orientations) 
     358 
     359    px, py, pz = transform_xyz(view, jitter, px, py, pz) 
     360    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz) 
     361 
     362    # TODO: zdir for labels is broken, and labels aren't appearing. 
     363    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 
     364        zdir = np.asarray(zdir).flatten() 
     365        ax.text(p[0], p[1], p[2], label, zdir=zdir) 
     366 
     367# Definition of rotation matrices comes from wikipedia: 
     368#    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations 
    170369def Rx(angle): 
     370    """Construct a matrix to rotate points about *x* by *angle* degrees.""" 
    171371    a = radians(angle) 
    172     R = [[1., 0., 0.], 
    173          [0.,  cos(a), sin(a)], 
    174          [0., -sin(a), cos(a)]] 
     372    R = [[1, 0, 0], 
     373         [0, +cos(a), -sin(a)], 
     374         [0, +sin(a), +cos(a)]] 
    175375    return np.matrix(R) 
    176376 
    177377def Ry(angle): 
     378    """Construct a matrix to rotate points about *y* by *angle* degrees.""" 
    178379    a = radians(angle) 
    179     R = [[cos(a), 0., -sin(a)], 
    180          [0., 1., 0.], 
    181          [sin(a), 0.,  cos(a)]] 
     380    R = [[+cos(a), 0, +sin(a)], 
     381         [0, 1, 0], 
     382         [-sin(a), 0, +cos(a)]] 
    182383    return np.matrix(R) 
    183384 
    184385def Rz(angle): 
     386    """Construct a matrix to rotate points about *z* by *angle* degrees.""" 
    185387    a = radians(angle) 
    186     R = [[cos(a), -sin(a), 0.], 
    187          [sin(a),  cos(a), 0.], 
    188          [0., 0., 1.]] 
     388    R = [[+cos(a), -sin(a), 0], 
     389         [+sin(a), +cos(a), 0], 
     390         [0, 0, 1]] 
    189391    return np.matrix(R) 
    190392 
    191 def main(): 
     393def transform_xyz(view, jitter, x, y, z): 
     394    """ 
     395    Send a set of (x,y,z) points through the jitter and view transforms. 
     396    """ 
     397    x, y, z = [np.asarray(v) for v in (x, y, z)] 
     398    shape = x.shape 
     399    points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     400    points = apply_jitter(jitter, points) 
     401    points = orient_relative_to_beam(view, points) 
     402    x, y, z = [np.array(v).reshape(shape) for v in points] 
     403    return x, y, z 
     404 
     405def apply_jitter(jitter, points): 
     406    """ 
     407    Apply the jitter transform to a set of points. 
     408 
     409    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     410    """ 
     411    dtheta, dphi, dpsi = jitter 
     412    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     413    return points 
     414 
     415def orient_relative_to_beam(view, points): 
     416    """ 
     417    Apply the view transform to a set of points. 
     418 
     419    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     420    """ 
     421    theta, phi, psi = view 
     422    points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     423    return points 
     424 
     425# translate between number of dimension of dispersity and the number of 
     426# points along each dimension. 
     427PD_N_TABLE = { 
     428    (0, 0, 0): (0, 0, 0),     # 0 
     429    (1, 0, 0): (100, 0, 0),   # 100 
     430    (0, 1, 0): (0, 100, 0), 
     431    (0, 0, 1): (0, 0, 100), 
     432    (1, 1, 0): (30, 30, 0),   # 900 
     433    (1, 0, 1): (30, 0, 30), 
     434    (0, 1, 1): (0, 30, 30), 
     435    (1, 1, 1): (15, 15, 15),  # 3375 
     436} 
     437 
     438def clipped_range(data, portion=1.0, mode='central'): 
     439    """ 
     440    Determine range from data. 
     441 
     442    If *portion* is 1, use full range, otherwise use the center of the range 
     443    or the top of the range, depending on whether *mode* is 'central' or 'top'. 
     444    """ 
     445    if portion == 1.0: 
     446        return data.min(), data.max() 
     447    elif mode == 'central': 
     448        data = np.sort(data.flatten()) 
     449        offset = int(portion*len(data)/2 + 0.5) 
     450        return data[offset], data[-offset] 
     451    elif mode == 'top': 
     452        data = np.sort(data.flatten()) 
     453        offset = int(portion*len(data) + 0.5) 
     454        return data[offset], data[-1] 
     455 
     456def draw_scattering(calculator, ax, view, jitter, dist='gaussian'): 
     457    """ 
     458    Plot the scattering for the particular view. 
     459 
     460    *calculator* is returned from :func:`build_model`.  *ax* are the 3D axes 
     461    on which the data will be plotted.  *view* and *jitter* are the current 
     462    orientation and orientation dispersity.  *dist* is one of the sasmodels 
     463    weight distributions. 
     464    """ 
     465    if dist == 'uniform':  # uniform is not yet in this branch 
     466        dist, scale = 'rectangle', 1/sqrt(3) 
     467    else: 
     468        scale = 1 
     469 
     470    # add the orientation parameters to the model parameters 
     471    theta, phi, psi = view 
     472    theta_pd, phi_pd, psi_pd = [scale*v for v in jitter] 
     473    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd>0, phi_pd>0, psi_pd>0)] 
     474    ## increase pd_n for testing jitter integration rather than simple viz 
     475    #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)] 
     476 
     477    pars = dict( 
     478        theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n, 
     479        phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n, 
     480        psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n, 
     481    ) 
     482    pars.update(calculator.pars) 
     483 
     484    # compute the pattern 
     485    qx, qy = calculator._data.x_bins, calculator._data.y_bins 
     486    Iqxy = calculator(**pars).reshape(len(qx), len(qy)) 
     487 
     488    # scale it and draw it 
     489    Iqxy = np.log(Iqxy) 
     490    if calculator.limits: 
     491        # use limits from orientation (0,0,0) 
     492        vmin, vmax = calculator.limits 
     493    else: 
     494        vmax = Iqxy.max() 
     495        vmin = vmax*10**-7 
     496        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top') 
     497    #print("range",(vmin,vmax)) 
     498    #qx, qy = np.meshgrid(qx, qy) 
     499    if 0: 
     500        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
     501        level[level<0] = 0 
     502        colors = plt.get_cmap()(level) 
     503        ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
     504    elif 1: 
     505        ax.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
     506                    levels=np.linspace(vmin, vmax, 24)) 
     507    else: 
     508        ax.pcolormesh(qx, qy, Iqxy) 
     509 
     510def build_model(model_name, n=150, qmax=0.5, **pars): 
     511    """ 
     512    Build a calculator for the given shape. 
     513 
     514    *model_name* is any sasmodels model.  *n* and *qmax* define an n x n mesh 
     515    on which to evaluate the model.  The remaining parameters are stored in 
     516    the returned calculator as *calculator.pars*.  They are used by 
     517    :func:`draw_scattering` to set the non-orientation parameters in the 
     518    calculation. 
     519 
     520    Returns a *calculator* function which takes a dictionary or parameters and 
     521    produces Iqxy.  The Iqxy value needs to be reshaped to an n x n matrix 
     522    for plotting.  See the :class:`sasmodels.direct_model.DirectModel` class 
     523    for details. 
     524    """ 
     525    from sasmodels.core import load_model_info, build_model 
     526    from sasmodels.data import empty_data2D 
     527    from sasmodels.direct_model import DirectModel 
     528 
     529    model_info = load_model_info(model_name) 
     530    model = build_model(model_info) #, dtype='double!') 
     531    q = np.linspace(-qmax, qmax, n) 
     532    data = empty_data2D(q, q) 
     533    calculator = DirectModel(data, model) 
     534 
     535    # stuff the values for non-orientation parameters into the calculator 
     536    calculator.pars = pars.copy() 
     537    calculator.pars.setdefault('backgound', 1e-3) 
     538 
     539    # fix the data limits so that we can see if the pattern fades 
     540    # under rotation or angular dispersion 
     541    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars) 
     542    Iqxy = np.log(Iqxy) 
     543    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top') 
     544    calculator.limits = vmin, vmax+1 
     545 
     546    return calculator 
     547 
     548def select_calculator(model_name, n=150, size=(10,40,100)): 
     549    """ 
     550    Create a model calculator for the given shape. 
     551 
     552    *model_name* is one of sphere, cylinder, ellipsoid, triaxial_ellipsoid, 
     553    parallelepiped or bcc_paracrystal. *n* is the number of points to use 
     554    in the q range.  *qmax* is chosen based on model parameters for the 
     555    given model to show something intersting. 
     556 
     557    Returns *calculator* and tuple *size* (a,b,c) giving minor and major 
     558    equitorial axes and polar axis respectively.  See :func:`build_model` 
     559    for details on the returned calculator. 
     560    """ 
     561    a, b, c = size 
     562    if model_name == 'sphere': 
     563        calculator = build_model('sphere', n=n, radius=c) 
     564        a = b = c 
     565    elif model_name == 'bcc_paracrystal': 
     566        calculator = build_model('bcc_paracrystal', n=n, dnn=c, 
     567                                  d_factor=0.06, radius=40) 
     568        a = b = c 
     569    elif model_name == 'cylinder': 
     570        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c) 
     571        a = b 
     572    elif model_name == 'ellipsoid': 
     573        calculator = build_model('ellipsoid', n=n, qmax=1.0, 
     574                                 radius_polar=c, radius_equatorial=b) 
     575        a = b 
     576    elif model_name == 'triaxial_ellipsoid': 
     577        calculator = build_model('triaxial_ellipsoid', n=n, qmax=0.5, 
     578                                 radius_equat_minor=a, 
     579                                 radius_equat_major=b, 
     580                                 radius_polar=c) 
     581    elif model_name == 'parallelepiped': 
     582        calculator = build_model('parallelepiped', n=n, a=a, b=b, c=c) 
     583    else: 
     584        raise ValueError("unknown model %s"%model_name) 
     585 
     586    return calculator, (a, b, c) 
     587 
     588SHAPES = [ 
     589    'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal', 
     590    'cylinder', 'ellipsoid', 
     591    'sphere', 
     592 ] 
     593 
     594DISTRIBUTIONS = [ 
     595    'gaussian', 'rectangle', 'uniform', 
     596] 
     597DIST_LIMITS = { 
     598    'gaussian': 30, 
     599    'rectangle': 90/sqrt(3), 
     600    'uniform': 90, 
     601} 
     602 
     603def run(model_name='parallelepiped', size=(10, 40, 100), 
     604        dist='gaussian', mesh=30, 
     605        projection='equirectangular'): 
     606    """ 
     607    Show an interactive orientation and jitter demo. 
     608 
     609    *model_name* is one of the models available in :func:`select_model`. 
     610    """ 
     611    # projection number according to 1-order position in list, but 
     612    # only 1 and 2 are implemented so far. 
     613    from sasmodels import generate 
     614    generate.PROJECTION = PROJECTIONS.index(projection) + 1 
     615    if generate.PROJECTION > 2: 
     616        print("*** PROJECTION %s not implemented in scattering function ***"%projection) 
     617        generate.PROJECTION = 2 
     618 
     619    # set up calculator 
     620    calculator, size = select_calculator(model_name, n=150, size=size) 
     621 
     622    ## uncomment to set an independent the colour range for every view 
     623    ## If left commented, the colour range is fixed for all views 
     624    calculator.limits = None 
     625 
     626    ## initial view 
     627    #theta, dtheta = 70., 10. 
     628    #phi, dphi = -45., 3. 
     629    #psi, dpsi = -45., 3. 
     630    theta, phi, psi = 0, 0, 0 
     631    dtheta, dphi, dpsi = 0, 0, 0 
     632 
     633    ## create the plot window 
    192634    #plt.hold(True) 
     635    plt.subplots(num=None, figsize=(5.5, 5.5)) 
    193636    plt.set_cmap('gist_earth') 
    194637    plt.clf() 
     638    plt.gcf().canvas.set_window_title(projection) 
    195639    #gs = gridspec.GridSpec(2,1,height_ratios=[4,1]) 
    196640    #ax = plt.subplot(gs[0], projection='3d') 
    197641    ax = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d') 
    198  
    199     theta, dtheta = 70., 10. 
    200     phi, dphi = -45., 3. 
    201     psi, dpsi = -45., 3. 
    202     theta, phi, psi = 0, 0, 0 
    203     dtheta, dphi, dpsi = 0, 0, 0 
    204     #dist = 'rect' 
    205     dist = 'gauss' 
     642    try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection 
     643        ax.axis('square') 
     644    except Exception: 
     645        pass 
    206646 
    207647    axcolor = 'lightgoldenrodyellow' 
     648 
     649    ## add control widgets to plot 
    208650    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    209651    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
     
    212654    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    213655    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
     656 
    214657    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    215658    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    216659    axdpsi= plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
    217     sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta) 
    218     sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi) 
    219     sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi) 
    220  
     660    # Note: using ridiculous definition of rectangle distribution, whose width 
     661    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
     662    # the maximum width to 90. 
     663    dlimit = DIST_LIMITS[dist] 
     664    sdtheta = Slider(axdtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
     665    sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
     666    sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
     667 
     668    ## callback to draw the new view 
    221669    def update(val, axis=None): 
    222         theta, phi, psi = stheta.val, sphi.val, spsi.val 
    223         dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val 
     670        view = stheta.val, sphi.val, spsi.val 
     671        jitter = sdtheta.val, sdphi.val, sdpsi.val 
     672        # set small jitter as 0 if multiple pd dims 
     673        dims = sum(v > 0 for v in jitter) 
     674        limit = [0, 0, 0.5, 5][dims] 
     675        jitter = [0 if v < limit else v for v in jitter] 
    224676        ax.cla() 
    225         draw_beam(ax) 
    226         draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi) 
    227         #if not axis.startswith('d'): 
    228         #    ax.view_init(elev=theta, azim=phi) 
     677        draw_beam(ax, (0, 0)) 
     678        draw_jitter(ax, view, jitter, dist=dist, size=size) 
     679        #draw_jitter(ax, view, (0,0,0)) 
     680        draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection) 
     681        draw_scattering(calculator, ax, view, jitter, dist=dist) 
    229682        plt.gcf().canvas.draw() 
    230683 
     684    ## bind control widgets to view updater 
    231685    stheta.on_changed(lambda v: update(v,'theta')) 
    232686    sphi.on_changed(lambda v: update(v, 'phi')) 
     
    236690    sdpsi.on_changed(lambda v: update(v, 'dpsi')) 
    237691 
     692    ## initialize view 
    238693    update(None, 'phi') 
    239694 
     695    ## go interactive 
    240696    plt.show() 
     697 
     698def main(): 
     699    parser = argparse.ArgumentParser( 
     700        description="Display jitter", 
     701        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
     702        ) 
     703    parser.add_argument('-p', '--projection', choices=PROJECTIONS, default=PROJECTIONS[0], help='coordinate projection') 
     704    parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths') 
     705    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='jitter distribution') 
     706    parser.add_argument('-m', '--mesh', type=int, default=30, help='#points in theta-phi mesh') 
     707    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape') 
     708    opts = parser.parse_args() 
     709    size = tuple(int(v) for v in opts.size.split(',')) 
     710    run(opts.shape, size=size, 
     711        mesh=opts.mesh, dist=opts.distribution, 
     712        projection=opts.projection) 
    241713 
    242714if __name__ == "__main__": 
  • explore/precision.py

    r237c9cf ra1c5758  
    11#!/usr/bin/env python 
    22r""" 
    3 Show numerical precision of $2 J_1(x)/x$. 
     3Show numerical precision of various expressions. 
     4 
     5Evaluates the same function(s) in single and double precision and compares 
     6the results to 500 digit mpmath evaluation of the same function. 
     7 
     8Note: a quick way to generation C and python code for taylor series 
     9expansions from sympy: 
     10 
     11    import sympy as sp 
     12    x = sp.var("x") 
     13    f = sp.sin(x)/x 
     14    t = sp.series(f, n=12).removeO()  # taylor series with no O(x^n) term 
     15    p = sp.horner(t)   # Horner representation 
     16    p = p.replace(x**2, sp.var("xsq")  # simplify if alternate terms are zero 
     17    p = p.n(15)  # evaluate coefficients to 15 digits (optional) 
     18    c_code = sp.ccode(p, assign_to=sp.var("p"))  # convert to c code 
     19    py_code = c[:-1]  # strip semicolon to convert c to python 
     20 
     21    # mpmath has pade() rational function approximation, which might work 
     22    # better than the taylor series for some functions: 
     23    P, Q = mp.pade(sp.Poly(t.n(15),x).coeffs(), L, M) 
     24    P = sum(a*x**n for n,a in enumerate(reversed(P))) 
     25    Q = sum(a*x**n for n,a in enumerate(reversed(Q))) 
     26    c_code = sp.ccode(sp.horner(P)/sp.horner(Q), assign_to=sp.var("p")) 
     27 
     28    # There are richardson and shanks series accelerators in both sympy 
     29    # and mpmath that may be helpful. 
    430""" 
    531from __future__ import division, print_function 
     
    284310    np_function=scipy.special.erfc, 
    285311    ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]), 
     312    limits=(-5., 5.), 
     313) 
     314add_function( 
     315    name="expm1(x)", 
     316    mp_function=mp.expm1, 
     317    np_function=np.expm1, 
     318    ocl_function=make_ocl("return expm1(q);", "sas_expm1"), 
    286319    limits=(-5., 5.), 
    287320) 
     
    448481) 
    449482 
     483replacement_expm1 = """\ 
     484      double x = (double)q;  // go back to float for single precision kernels 
     485      // Adapted from the cephes math library. 
     486      // Copyright 1984 - 1992 by Stephen L. Moshier 
     487      if (x != x || x == 0.0) { 
     488         return x; // NaN and +/- 0 
     489      } else if (x < -0.5 || x > 0.5) { 
     490         return exp(x) - 1.0; 
     491      } else { 
     492         const double xsq = x*x; 
     493         const double p = ((( 
     494            +1.2617719307481059087798E-4)*xsq 
     495            +3.0299440770744196129956E-2)*xsq 
     496            +9.9999999999999999991025E-1); 
     497         const double q = (((( 
     498            +3.0019850513866445504159E-6)*xsq 
     499            +2.5244834034968410419224E-3)*xsq 
     500            +2.2726554820815502876593E-1)*xsq 
     501            +2.0000000000000000000897E0); 
     502         double r = x * p; 
     503         r =  r / (q - r); 
     504         return r+r; 
     505       } 
     506""" 
     507add_function( 
     508    name="sas_expm1(x)", 
     509    mp_function=mp.expm1, 
     510    np_function=np.expm1, 
     511    ocl_function=make_ocl(replacement_expm1, "sas_expm1"), 
     512) 
     513 
    450514# Alternate versions of 3 j1(x)/x, for posterity 
    451515def taylor_3j1x_x(x): 
  • sasmodels/compare.py

    rced5bd2 r3bfd924  
    4242from . import exception 
    4343from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    44 from .direct_model import DirectModel 
     44from .direct_model import DirectModel, get_mesh 
    4545from .convert import revert_name, revert_pars, constrain_new_to_old 
    4646from .generate import FLOAT_RE 
     47from .weights import plot_weights 
    4748 
    4849try: 
     
    8384    -pars/-nopars* prints the parameter set or not 
    8485    -default/-demo* use demo vs default parameters 
     86    -sphere[=150] set up spherical integration over theta/phi using n points 
    8587 
    8688    === calculation options === 
    87     -mono*/-poly force monodisperse or allow polydisperse demo parameters 
     89    -mono*/-poly force monodisperse or allow polydisperse random parameters 
    8890    -cutoff=1e-5* cutoff value for including a point in polydispersity 
    8991    -magnetic/-nonmagnetic* suppress magnetism 
     
    9294 
    9395    === precision options === 
    94     -calc=default uses the default calcution precision 
     96    -engine=default uses the default calcution precision 
    9597    -single/-double/-half/-fast sets an OpenCL calculation engine 
    9698    -single!/-double!/-quad! sets an OpenMP calculation engine 
     
    103105    -abs/-rel* plot relative or absolute error 
    104106    -title="note" adds note to the plot title, after the model name 
     107    -weights shows weights plots for the polydisperse parameters 
    105108 
    106109    === output options === 
     
    111114vary from 64-bit to 128-bit, with 80-bit floats being common (1e-19 precision). 
    112115On unix and mac you may need single quotes around the DLL computation 
    113 engines, such as -calc='single!,double!' since !, is treated as a history 
     116engines, such as -engine='single!,double!' since !, is treated as a history 
    114117expansion request in the shell. 
    115118 
     
    123126 
    124127    # compare single and double precision calculation for a barbell 
    125     sascomp barbell -calc=single,double 
     128    sascomp barbell -engine=single,double 
    126129 
    127130    # generate 10 random lorentz models, with seed=27 
     
    132135 
    133136    # model timing test requires multiple evals to perform the estimate 
    134     sascomp pringle -calc=single,double -timing=100,100 -noplot 
     137    sascomp pringle -engine=single,double -timing=100,100 -noplot 
    135138""" 
    136139 
     
    278281        # orientation in [-180,180], orientation pd in [0,45] 
    279282        if p.endswith('_pd'): 
    280             return 0., 45. 
     283            return 0., 180. 
    281284        else: 
    282285            return -180., 180. 
     
    433436 
    434437 
     438def limit_dimensions(model_info, pars, maxdim): 
     439    # type: (ModelInfo, ParameterSet, float) -> None 
     440    """ 
     441    Limit parameters of units of Ang to maxdim. 
     442    """ 
     443    for p in model_info.parameters.call_parameters: 
     444        value = pars[p.name] 
     445        if p.units == 'Ang' and value > maxdim: 
     446            pars[p.name] = maxdim*10**np.random.uniform(-3,0) 
     447 
    435448def constrain_pars(model_info, pars): 
    436449    # type: (ModelInfo, ParameterSet) -> None 
     
    495508    Format the parameter list for printing. 
    496509    """ 
     510    is2d = True 
    497511    lines = [] 
    498512    parameters = model_info.parameters 
     
    815829 
    816830 
    817 def compare(opts, limits=None): 
     831def compare(opts, limits=None, maxdim=np.inf): 
    818832    # type: (Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float] 
    819833    """ 
     
    826840    as the values are adjusted, making it easier to see the effects of the 
    827841    parameters. 
    828     """ 
    829     limits = np.Inf, -np.Inf 
     842 
     843    *maxdim* is the maximum value for any parameter with units of Angstrom. 
     844    """ 
    830845    for k in range(opts['sets']): 
    831         opts['pars'] = parse_pars(opts) 
     846        if k > 1: 
     847            # print a separate seed for each dataset for better reproducibility 
     848            new_seed = np.random.randint(1000000) 
     849            print("Set %d uses -random=%i"%(k+1,new_seed)) 
     850            np.random.seed(new_seed) 
     851        opts['pars'] = parse_pars(opts, maxdim=maxdim) 
    832852        if opts['pars'] is None: 
    833853            return 
     
    835855        if opts['plot']: 
    836856            limits = plot_models(opts, result, limits=limits, setnum=k) 
     857        if opts['show_weights']: 
     858            base, _ = opts['engines'] 
     859            base_pars, _ = opts['pars'] 
     860            model_info = base._kernel.info 
     861            dim = base._kernel.dim 
     862            plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim)) 
    837863    if opts['plot']: 
    838864        import matplotlib.pyplot as plt 
    839865        plt.show() 
     866    return limits 
    840867 
    841868def run_models(opts, verbose=False): 
     
    912939 
    913940 
    914 def plot_models(opts, result, limits=(np.Inf, -np.Inf), setnum=0): 
     941def plot_models(opts, result, limits=None, setnum=0): 
    915942    # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float] 
     943    import matplotlib.pyplot as plt 
     944 
    916945    base_value, comp_value = result['base_value'], result['comp_value'] 
    917946    base_time, comp_time = result['base_time'], result['comp_time'] 
     
    925954    # Plot if requested 
    926955    view = opts['view'] 
    927     import matplotlib.pyplot as plt 
    928     vmin, vmax = limits 
    929     if have_base: 
    930         vmin = min(vmin, base_value.min()) 
    931         vmax = max(vmax, base_value.max()) 
    932     if have_comp: 
    933         vmin = min(vmin, comp_value.min()) 
    934         vmax = max(vmax, comp_value.max()) 
    935     limits = vmin, vmax 
     956    if limits is None: 
     957        vmin, vmax = np.inf, -np.inf 
     958        if have_base: 
     959            vmin = min(vmin, base_value.min()) 
     960            vmax = max(vmax, base_value.max()) 
     961        if have_comp: 
     962            vmin = min(vmin, comp_value.min()) 
     963            vmax = max(vmax, comp_value.max()) 
     964        limits = vmin, vmax 
    936965 
    937966    if have_base: 
     
    962991            err[err > cutoff] = cutoff 
    963992        #err,errstr = base/comp,"ratio" 
    964         plot_theory(data, None, resid=err, view=view, use_data=use_data) 
    965         plt.yscale(errview) 
     993        plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
     994        plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear') 
     995        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best') 
    966996        plt.title("max %s = %.3g"%(errstr, abs(err).max())) 
    967997        #cbar_title = errstr if errview=="linear" else "log "+errstr 
     
    9951025OPTIONS = [ 
    9961026    # Plotting 
    997     'plot', 'noplot', 
     1027    'plot', 'noplot', 'weights', 
    9981028    'linear', 'log', 'q4', 
    9991029    'rel', 'abs', 
     
    10101040    'demo', 'default',  # TODO: remove demo/default 
    10111041    'nopars', 'pars', 
     1042    'sphere', 'sphere=', # integrate over a sphere in 2d with n points 
    10121043 
    10131044    # Calculation options 
     
    10181049 
    10191050    # Precision options 
    1020     'calc=', 
     1051    'engine=', 
    10211052    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!', 
    10221053    'sasview',  # TODO: remove sasview 3.x support 
     
    11451176        'sets'      : 0, 
    11461177        'engine'    : 'default', 
    1147         'evals'     : '1', 
     1178        'count'     : '1', 
     1179        'show_weights' : False, 
     1180        'sphere'    : 0, 
    11481181    } 
    11491182    for arg in flags: 
     
    11681201        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:] 
    11691202        elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:] 
    1170         elif arg.startswith('-random='):   opts['seed'] = int(arg[8:]) 
    11711203        elif arg.startswith('-title='):    opts['title'] = arg[7:] 
    11721204        elif arg.startswith('-data='):     opts['datafile'] = arg[6:] 
    1173         elif arg.startswith('-calc='):     opts['engine'] = arg[6:] 
    1174         elif arg.startswith('-neval='):    opts['evals'] = arg[7:] 
    1175         elif arg == '-random':  opts['seed'] = np.random.randint(1000000) 
     1205        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
     1206        elif arg.startswith('-neval='):    opts['count'] = arg[7:] 
     1207        elif arg.startswith('-random='): 
     1208            opts['seed'] = int(arg[8:]) 
     1209            opts['sets'] = 0 
     1210        elif arg == '-random': 
     1211            opts['seed'] = np.random.randint(1000000) 
     1212            opts['sets'] = 0 
     1213        elif arg.startswith('-sphere'): 
     1214            opts['sphere'] = int(arg[8:]) if len(arg) > 7 else 150 
     1215            opts['is2d'] = True 
    11761216        elif arg == '-preset':  opts['seed'] = -1 
    11771217        elif arg == '-mono':    opts['mono'] = True 
     
    11961236        elif arg == '-demo':    opts['use_demo'] = True 
    11971237        elif arg == '-default': opts['use_demo'] = False 
     1238        elif arg == '-weights': opts['show_weights'] = True 
    11981239        elif arg == '-html':    opts['html'] = True 
    11991240        elif arg == '-help':    opts['html'] = True 
     
    12321273 
    12331274    if PAR_SPLIT in opts['engine']: 
    1234         engine_types = opts['engine'].split(PAR_SPLIT, 2) 
     1275        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
    12351276        comparison = True 
    12361277    else: 
    1237         engine_types = [opts['engine']]*2 
    1238  
    1239     if PAR_SPLIT in opts['evals']: 
    1240         evals = [int(k) for k in opts['evals'].split(PAR_SPLIT, 2)] 
     1278        opts['engine'] = [opts['engine']]*2 
     1279 
     1280    if PAR_SPLIT in opts['count']: 
     1281        opts['count'] = [int(k) for k in opts['count'].split(PAR_SPLIT, 2)] 
    12411282        comparison = True 
    12421283    else: 
    1243         evals = [int(opts['evals'])]*2 
     1284        opts['count'] = [int(opts['count'])]*2 
    12441285 
    12451286    if PAR_SPLIT in opts['cutoff']: 
    1246         cutoff = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)] 
     1287        opts['cutoff'] = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)] 
    12471288        comparison = True 
    12481289    else: 
    1249         cutoff = [float(opts['cutoff'])]*2 
    1250  
    1251     base = make_engine(model_info[0], data, engine_types[0], cutoff[0]) 
     1290        opts['cutoff'] = [float(opts['cutoff'])]*2 
     1291 
     1292    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
    12521293    if comparison: 
    1253         comp = make_engine(model_info[1], data, engine_types[1], cutoff[1]) 
     1294        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
    12541295    else: 
    12551296        comp = None 
     
    12601301        'data'      : data, 
    12611302        'name'      : names, 
    1262         'def'       : model_info, 
    1263         'count'     : evals, 
     1303        'info'      : model_info, 
    12641304        'engines'   : [base, comp], 
    12651305        'values'    : values, 
     
    12671307    # pylint: enable=bad-whitespace 
    12681308 
     1309    # Set the integration parameters to the half sphere 
     1310    if opts['sphere'] > 0: 
     1311        set_spherical_integration_parameters(opts, opts['sphere']) 
     1312 
    12691313    return opts 
    12701314 
    1271 def parse_pars(opts): 
    1272     model_info, model_info2 = opts['def'] 
     1315def set_spherical_integration_parameters(opts, steps): 
     1316    """ 
     1317    Set integration parameters for spherical integration over the entire 
     1318    surface in theta-phi coordinates. 
     1319    """ 
     1320    # Set the integration parameters to the half sphere 
     1321    opts['values'].extend([ 
     1322        #'theta=90', 
     1323        'theta_pd=%g'%(90/np.sqrt(3)), 
     1324        'theta_pd_n=%d'%steps, 
     1325        'theta_pd_type=rectangle', 
     1326        #'phi=0', 
     1327        'phi_pd=%g'%(180/np.sqrt(3)), 
     1328        'phi_pd_n=%d'%(2*steps), 
     1329        'phi_pd_type=rectangle', 
     1330        #'background=0', 
     1331    ]) 
     1332    if 'psi' in opts['info'][0].parameters: 
     1333        opts['values'].extend([ 
     1334            #'psi=0', 
     1335            'psi_pd=%g'%(180/np.sqrt(3)), 
     1336            'psi_pd_n=%d'%(2*steps), 
     1337            'psi_pd_type=rectangle', 
     1338        ]) 
     1339        pass 
     1340 
     1341def parse_pars(opts, maxdim=np.inf): 
     1342    model_info, model_info2 = opts['info'] 
    12731343 
    12741344    # Get demo parameters from model definition, or use default parameters 
     
    12811351    if opts['seed'] > -1: 
    12821352        pars = randomize_pars(model_info, pars) 
     1353        limit_dimensions(model_info, pars, maxdim) 
    12831354        if model_info != model_info2: 
    12841355            pars2 = randomize_pars(model_info2, pars2) 
     1356            limit_dimensions(model_info, pars2, maxdim) 
    12851357            # Share values for parameters with the same name 
    12861358            for k, v in pars.items(): 
     
    13591431    from . import rst2html 
    13601432 
    1361     info = opts['def'][0] 
     1433    info = opts['info'][0] 
    13621434    html = make_html(info) 
    13631435    path = os.path.dirname(info.filename) 
     
    14101482        opts['pars'] = list(opts['pars']) 
    14111483        p1, p2 = opts['pars'] 
    1412         m1, m2 = opts['def'] 
     1484        m1, m2 = opts['info'] 
    14131485        self.fix_p2 = m1 != m2 or p1 != p2 
    14141486        model_info = m1 
     
    14291501        self.starting_values = dict((k, v.value) for k, v in pars.items()) 
    14301502        self.pd_types = pd_types 
    1431         self.limits = np.Inf, -np.Inf 
     1503        self.limits = None 
    14321504 
    14331505    def revert_values(self): 
  • sasmodels/core.py

    ra85a569 r9e771a3  
    272272    Possible types include 'half', 'single', 'double' and 'quad'.  If the 
    273273    type is 'fast', then this is equivalent to dtype 'single' but using 
    274     fast native functions rather than those with the precision level guaranteed 
    275     by the OpenCL standard. 
     274    fast native functions rather than those with the precision level 
     275    guaranteed by the OpenCL standard.  'default' will choose the appropriate 
     276    default for the model and platform. 
    276277 
    277278    Platform preference can be specfied ("ocl" vs "dll"), with the default 
  • sasmodels/data.py

    rba7302a ra1c5758  
    363363    if hasattr(data, 'isSesans') and data.isSesans: 
    364364        _plot_result_sesans(data, None, None, use_data=True, limits=limits) 
    365     elif hasattr(data, 'qx_data'): 
     365    elif hasattr(data, 'qx_data') and not getattr(data, 'radial', False): 
    366366        _plot_result2D(data, None, None, view, use_data=True, limits=limits) 
    367367    else: 
     
    391391    if hasattr(data, 'isSesans') and data.isSesans: 
    392392        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 
    393     elif hasattr(data, 'qx_data'): 
     393    elif hasattr(data, 'qx_data') and not getattr(data, 'radial', False): 
    394394        _plot_result2D(data, theory, resid, view, use_data, limits=limits) 
    395395    else: 
     
    425425    import matplotlib.pyplot as plt  # type: ignore 
    426426    from numpy.ma import masked_array, masked  # type: ignore 
     427 
     428    if getattr(data, 'radial', False): 
     429        radial_data.x = radial_data.q_data 
     430        radial_data.y = radial_data.data 
    427431 
    428432    use_data = use_data and data.y is not None 
     
    651655    ymin, ymax = min(data.qy_data), max(data.qy_data) 
    652656    if view == 'log': 
    653         vmin, vmax = np.log10(vmin), np.log10(vmax) 
     657        vmin_scaled, vmax_scaled= np.log10(vmin), np.log10(vmax) 
     658    else: 
     659        vmin_scaled, vmax_scaled = vmin, vmax 
    654660    plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)), 
    655661               interpolation='nearest', aspect=1, origin='lower', 
    656                extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
     662               extent=[xmin, xmax, ymin, ymax], 
     663               vmin=vmin_scaled, vmax=vmax_scaled) 
    657664    plt.xlabel("$q_x$/A$^{-1}$") 
    658665    plt.ylabel("$q_y$/A$^{-1}$") 
  • sasmodels/details.py

    rccd5f01 rce99754  
    1515 
    1616import numpy as np  # type: ignore 
    17 from numpy import pi, cos, sin 
     17from numpy import pi, cos, sin, radians 
    1818 
    1919try: 
     
    2929 
    3030try: 
    31     from typing import List 
     31    from typing import List, Tuple, Sequence 
    3232except ImportError: 
    3333    pass 
    3434else: 
    3535    from .modelinfo import ModelInfo 
     36    from .modelinfo import ParameterTable 
    3637 
    3738 
     
    5354    coordinates, the total circumference decreases as latitude varies from 
    5455    pi r^2 at the equator to 0 at the pole, and the weight associated 
    55     with a range of phi values needs to be scaled by this circumference. 
     56    with a range of latitude values needs to be scaled by this circumference. 
    5657    This scale factor needs to be updated each time the theta value 
    5758    changes.  *theta_par* indicates which of the values in the parameter 
     
    192193    #print("off:",offset) 
    193194 
    194     # Check that we arn't using too many polydispersity loops 
     195    # Check that we aren't using too many polydispersity loops 
    195196    num_active = np.sum(length > 1) 
    196197    max_pd = model_info.parameters.max_pd 
     
    218219 
    219220ZEROS = tuple([0.]*31) 
    220 def make_kernel_args(kernel, pairs): 
     221def make_kernel_args(kernel, mesh): 
    221222    # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool] 
    222223    """ 
    223     Converts (value, weight) pairs into parameters for the kernel call. 
     224    Converts (value, dispersity, weight) for each parameter into kernel pars. 
    224225 
    225226    Returns a CallDetails object indicating the polydispersity, a data object 
     
    230231    npars = kernel.info.parameters.npars 
    231232    nvalues = kernel.info.parameters.nvalues 
    232     scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs] 
    233     values, weights = zip(*pairs[2:npars+2]) if npars else ((),()) 
     233    scalars = [value for value, dispersity, weight in mesh] 
     234    # skipping scale and background when building values and weights 
     235    values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 
     236    #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights) 
    234237    length = np.array([len(w) for w in weights]) 
    235238    offset = np.cumsum(np.hstack((0, length))) 
    236239    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
    237240    # Pad value array to a 32 value boundaryd 
    238     data_len = nvalues + 2*sum(len(v) for v in values) 
     241    data_len = nvalues + 2*sum(len(v) for v in dispersity) 
    239242    extra = (32 - data_len%32)%32 
    240     data = np.hstack((scalars,) + values + weights + ZEROS[:extra]) 
     243    data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra]) 
    241244    data = data.astype(kernel.dtype) 
    242245    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
     
    244247    return call_details, data, is_magnetic 
    245248 
     249def correct_theta_weights(parameters, dispersity, weights): 
     250    # type: (ParameterTable, Sequence[np.ndarray], Sequence[np.ndarray]) -> Sequence[np.ndarray] 
     251    """ 
     252    If there is a theta parameter, update the weights of that parameter so that 
     253    the cosine weighting required for polar integration is preserved. 
     254 
     255    Avoid evaluation strictly at the pole, which would otherwise send the 
     256    weight to zero.  This is probably not a problem in practice (if dispersity 
     257    is +/- 90, then you probably should be using a 1-D model of the circular 
     258    average). 
     259 
     260    Note: scale and background parameters are not include in the tuples for 
     261    dispersity and weights, so index is parameters.theta_offset, not 
     262    parameters.theta_offset+2 
     263 
     264    Returns updated weights vectors 
     265    """ 
     266    # TODO: explain in a comment why scale and background are missing 
     267    # Apparently the parameters.theta_offset similarly skips scale and 
     268    # and background, so the indexing works out, but they are still shipped 
     269    # to the kernel, so we need to add two there. 
     270    if parameters.theta_offset >= 0: 
     271        index = parameters.theta_offset 
     272        theta = dispersity[index] 
     273        # TODO: modify the dispersity vector to avoid the theta=-90,90,270,... 
     274        theta_weight = abs(cos(radians(theta))) 
     275        weights = tuple(theta_weight*w if k == index else w 
     276                        for k, w in enumerate(weights)) 
     277    return weights 
     278 
    246279 
    247280def convert_magnetism(parameters, values): 
     281    # type: (ParameterTable, Sequence[np.ndarray]) -> bool 
    248282    """ 
    249283    Convert magnetism values from polar to rectangular coordinates. 
     
    255289    scale = mag[:,0] 
    256290    if np.any(scale): 
    257         theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. 
     291        theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
    258292        cos_theta = cos(theta) 
    259293        mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
     
    265299 
    266300 
    267 def dispersion_mesh(model_info, pars): 
     301def dispersion_mesh(model_info, mesh): 
    268302    # type: (ModelInfo) -> Tuple[List[np.ndarray], List[np.ndarray]] 
    269303    """ 
    270304    Create a mesh grid of dispersion parameters and weights. 
     305 
     306    *mesh* is a list of (value, dispersity, weights), where the values 
     307    are the individual parameter values, and (dispersity, weights) is 
     308    the distribution of parameter values. 
     309 
     310    Only the volume parameters should be included in this list.  Orientation 
     311    parameters do not affect the calculation of effective radius or volume 
     312    ratio.  This is convenient since it avoids the distinction between 
     313    value and dispersity that is present in orientation parameters but not 
     314    shape parameters. 
    271315 
    272316    Returns [p1,p2,...],w where pj is a vector of values for parameter j 
     
    274318    parameter set in the vector. 
    275319    """ 
    276     value, weight = zip(*pars) 
     320    _, dispersity, weight = zip(*mesh) 
    277321    #weight = [w if len(w)>0 else [1.] for w in weight] 
    278322    weight = np.vstack([v.flatten() for v in meshgrid(*weight)]) 
    279323    weight = np.prod(weight, axis=0) 
    280     value = [v.flatten() for v in meshgrid(*value)] 
     324    dispersity = [v.flatten() for v in meshgrid(*dispersity)] 
    281325    lengths = [par.length for par in model_info.parameters.kernel_parameters 
    282326               if par.type == 'volume'] 
     
    285329        offset = 0 
    286330        for n in lengths: 
    287             pars.append(np.vstack(value[offset:offset+n]) 
    288                         if n > 1 else value[offset]) 
     331            pars.append(np.vstack(dispersity[offset:offset+n]) 
     332                        if n > 1 else dispersity[offset]) 
    289333            offset += n 
    290         value = pars 
    291     return value, weight 
     334        dispersity = pars 
     335    return dispersity, weight 
  • sasmodels/direct_model.py

    rd1ff3a5 r9e771a3  
    5555    *mono* is True if polydispersity should be set to none on all parameters. 
    5656    """ 
    57     parameters = calculator.info.parameters 
    58     if mono: 
    59         active = lambda name: False 
    60     elif calculator.dim == '1d': 
    61         active = lambda name: name in parameters.pd_1d 
    62     elif calculator.dim == '2d': 
    63         active = lambda name: name in parameters.pd_2d 
    64     else: 
    65         active = lambda name: True 
    66  
    67     #print("pars",[p.id for p in parameters.call_parameters]) 
    68     vw_pairs = [(get_weights(p, pars) if active(p.name) 
    69                  else ([pars.get(p.name, p.default)], [1.0])) 
    70                 for p in parameters.call_parameters] 
    71  
    72     call_details, values, is_magnetic = make_kernel_args(calculator, vw_pairs) 
     57    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
     58    #print("pars", list(zip(*mesh))[0]) 
     59    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
    7360    #print("values:", values) 
    7461    return calculator(call_details, values, cutoff, is_magnetic) 
    75  
    7662 
    7763def call_ER(model_info, pars): 
     
    129115    return x, y, model_info.profile_axes 
    130116 
    131  
    132 def get_weights(parameter, values): 
    133     # type: (Parameter, Dict[str, float]) -> Tuple[np.ndarray, np.ndarray] 
     117def get_mesh(model_info, values, dim='1d', mono=False): 
     118    # type: (ModelInfo, Dict[str, float], str, bool) -> List[Tuple[float, np.ndarray, np.ndarry]] 
     119    """ 
     120    Retrieve the dispersity mesh described by the parameter set. 
     121 
     122    Returns a list of *(value, dispersity, weights)* with one tuple for each 
     123    parameter in the model call parameters.  Inactive parameters return the 
     124    default value with a weight of 1.0. 
     125    """ 
     126    parameters = model_info.parameters 
     127    if mono: 
     128        active = lambda name: False 
     129    elif dim == '1d': 
     130        active = lambda name: name in parameters.pd_1d 
     131    elif dim == '2d': 
     132        active = lambda name: name in parameters.pd_2d 
     133    else: 
     134        active = lambda name: True 
     135 
     136    #print("pars",[p.id for p in parameters.call_parameters]) 
     137    mesh = [_get_par_weights(p, values, active(p.name)) 
     138            for p in parameters.call_parameters] 
     139    return mesh 
     140 
     141 
     142def _get_par_weights(parameter, values, active=True): 
     143    # type: (Parameter, Dict[str, float]) -> Tuple[float, np.ndarray, np.ndarray] 
    134144    """ 
    135145    Generate the distribution for parameter *name* given the parameter values 
     
    140150    """ 
    141151    value = float(values.get(parameter.name, parameter.default)) 
    142     relative = parameter.relative_pd 
    143     limits = parameter.limits 
    144     disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
    145152    npts = values.get(parameter.name+'_pd_n', 0) 
    146153    width = values.get(parameter.name+'_pd', 0.0) 
    147     nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
    148     if npts == 0 or width == 0: 
    149         return [value], [1.0] 
    150     value, weight = weights.get_weights( 
    151         disperser, npts, width, nsigma, value, limits, relative) 
    152     return value, weight / np.sum(weight) 
    153  
    154  
    155 def _vol_pars(model_info, pars): 
     154    relative = parameter.relative_pd 
     155    if npts == 0 or width == 0.0 or not active: 
     156        # Note: orientation parameters have the viewing angle as the parameter 
     157        # value and the jitter in the distribution, so be sure to set the 
     158        # empty pd for orientation parameters to 0. 
     159        pd = [value if relative or not parameter.polydisperse else 0.0], [1.0] 
     160    else: 
     161        limits = parameter.limits 
     162        disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
     163        nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
     164        pd = weights.get_weights(disperser, npts, width, nsigma, 
     165                                value, limits, relative) 
     166    return value, pd[0], pd[1] 
     167 
     168 
     169def _vol_pars(model_info, values): 
    156170    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray] 
    157     vol_pars = [get_weights(p, pars) 
     171    vol_pars = [_get_par_weights(p, values) 
    158172                for p in model_info.parameters.call_parameters 
    159173                if p.type == 'volume'] 
    160174    #import pylab; pylab.plot(vol_pars[0][0],vol_pars[0][1]); pylab.show() 
    161     value, weight = dispersion_mesh(model_info, vol_pars) 
    162     return value, weight 
     175    dispersity, weight = dispersion_mesh(model_info, vol_pars) 
     176    return dispersity, weight 
    163177 
    164178 
  • sasmodels/generate.py

    r6db17bd rff10479  
    167167import string 
    168168from zlib import crc32 
     169from inspect import currentframe, getframeinfo 
    169170 
    170171import numpy as np  # type: ignore 
     
    178179except ImportError: 
    179180    pass 
     181 
     182# jitter projection to use in the kernel code.  See explore/jitter.py 
     183# for details.  To change it from a program, set generate.PROJECTION. 
     184PROJECTION = 1 
    180185 
    181186def get_data_path(external_dir, target_file): 
     
    689694    # Load templates and user code 
    690695    kernel_header = load_template('kernel_header.c') 
    691     dll_code = load_template('kernel_iq.c') 
    692     ocl_code = load_template('kernel_iq.cl') 
    693     #ocl_code = load_template('kernel_iq_local.cl') 
     696    kernel_code = load_template('kernel_iq.c') 
    694697    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
     698 
     699    # What kind of 2D model do we need? 
     700    xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str) 
     701               else 'qac' if not partable.is_asymmetric 
     702               else 'qabc') 
    695703 
    696704    # Build initial sources 
     
    717725 
    718726    # Define the parameter table 
    719     # TODO: plug in current line number 
    720     source.append('#line 542 "sasmodels/generate.py"') 
     727    lineno = getframeinfo(currentframe()).lineno + 2 
     728    source.append('#line %d "sasmodels/generate.py"'%lineno) 
     729    #source.append('introduce breakage in generate to test lineno reporting') 
    721730    source.append("#define PARAMETER_TABLE \\") 
    722731    source.append("\\\n".join(p.as_definition() 
     
    734743    source.append(call_volume) 
    735744 
    736     refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 
    737     call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 
    738     if _have_Iqxy(user_code) or isinstance(model_info.Iqxy, str): 
    739         # Call 2D model 
    740         refs = ["_q[2*_i]", "_q[2*_i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 
    741         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 
    742     else: 
    743         # Call 1D model with sqrt(qx^2 + qy^2) 
    744         #warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
    745         # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
    746         pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:] 
    747         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
     745    model_refs = _call_pars("_v.", partable.iq_parameters) 
     746    pars = ",".join(["_q"] + model_refs) 
     747    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     748    if xy_mode == 'qabc': 
     749        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
     750        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
     751        clear_iqxy = "#undef CALL_IQ_ABC" 
     752    elif xy_mode == 'qac': 
     753        pars = ",".join(["_qa", "_qc"] + model_refs) 
     754        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 
     755        clear_iqxy = "#undef CALL_IQ_AC" 
     756    else:  # xy_mode == 'qa' 
     757        pars = ",".join(["_qa"] + model_refs) 
     758        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
     759        clear_iqxy = "#undef CALL_IQ_A" 
    748760 
    749761    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
     
    756768    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic) 
    757769    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    758     for k, v in enumerate(magpars[:3]): 
    759         source.append("#define MAGNETIC_PAR%d %d"%(k+1, v)) 
     770    source.append("#define PROJECTION %d"%PROJECTION) 
    760771 
    761772    # TODO: allow mixed python/opencl kernels? 
    762773 
    763     ocl = kernels(ocl_code, call_iq, call_iqxy, model_info.name) 
    764     dll = kernels(dll_code, call_iq, call_iqxy, model_info.name) 
     774    ocl = kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
     775    dll = kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    765776    result = { 
    766777        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
     
    771782 
    772783 
    773 def kernels(kernel, call_iq, call_iqxy, name): 
     784def kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
    774785    # type: ([str,str], str, str, str) -> List[str] 
    775786    code = kernel[0] 
     
    791802        '#line 1 "%s Iqxy"' % path, 
    792803        code, 
    793         "#undef CALL_IQ", 
     804        clear_iqxy, 
    794805        "#undef KERNEL_NAME", 
    795806    ] 
     
    802813        '#line 1 "%s Imagnetic"' % path, 
    803814        code, 
     815        clear_iqxy, 
    804816        "#undef MAGNETIC", 
    805         "#undef CALL_IQ", 
    806817        "#undef KERNEL_NAME", 
    807818    ] 
  • sasmodels/kernel_header.c

    rbb4b509 r8698a0d  
    8787 
    8888#if defined(NEED_EXPM1) 
     89   // TODO: precision is a half digit lower than numpy on mac in [1e-7, 0.5] 
     90   // Run "explore/precision.py sas_expm1" to see this (may have to fiddle 
     91   // the xrange for log to see the complete range). 
    8992   static SAS_DOUBLE expm1(SAS_DOUBLE x_in) { 
    9093      double x = (double)x_in;  // go back to float for single precision kernels 
     
    147150inline double cube(double x) { return x*x*x; } 
    148151inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    149  
    150 // To rotate from the canonical position to theta, phi, psi, first rotate by 
    151 // psi about the major axis, oriented along z, which is a rotation in the 
    152 // detector plane xy. Next rotate by theta about the y axis, aligning the major 
    153 // axis in the xz plane. Finally, rotate by phi in the detector plane xy. 
    154 // To compute the scattering, undo these rotations in reverse order: 
    155 //     rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi 
    156 // The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit 
    157 // vector in the q direction. 
    158 // To change between counterclockwise and clockwise rotation, change the 
    159 // sign of phi and psi. 
    160  
    161 #if 1 
    162 //think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
    163 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
    164     SINCOS(phi*M_PI_180, sn, cn); \ 
    165     q = sqrt(qx*qx + qy*qy); \ 
    166     cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
    167     sn = sqrt(1 - cn*cn); \ 
    168     } while (0) 
    169 #else 
    170 // SasView 3.x definition of orientation 
    171 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
    172     SINCOS(theta*M_PI_180, sn, cn); \ 
    173     q = sqrt(qx*qx + qy*qy);\ 
    174     cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \ 
    175     sn = sqrt(1 - cn*cn); \ 
    176     } while (0) 
    177 #endif 
    178  
    179 #if 1 
    180 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \ 
    181     q = sqrt(qx*qx + qy*qy); \ 
    182     const double qxhat = qx/q; \ 
    183     const double qyhat = qy/q; \ 
    184     double sin_theta, cos_theta; \ 
    185     double sin_phi, cos_phi; \ 
    186     double sin_psi, cos_psi; \ 
    187     SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
    188     SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
    189     SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
    190     xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \ 
    191          + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \ 
    192     yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \ 
    193          + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \ 
    194     zhat = qxhat*(-sin_theta*cos_phi) \ 
    195          + qyhat*(-sin_theta*sin_phi); \ 
    196     } while (0) 
    197 #else 
    198 // SasView 3.x definition of orientation 
    199 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
    200     q = sqrt(qx*qx + qy*qy); \ 
    201     const double qxhat = qx/q; \ 
    202     const double qyhat = qy/q; \ 
    203     double sin_theta, cos_theta; \ 
    204     double sin_phi, cos_phi; \ 
    205     double sin_psi, cos_psi; \ 
    206     SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
    207     SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
    208     SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
    209     cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ 
    210     cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ 
    211     cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
    212     } while (0) 
    213 #endif 
  • sasmodels/kernel_iq.c

    rbde38b5 r6aee3ab  
    1  
    21/* 
    32    ########################################################## 
     
    1211*/ 
    1312 
     13// NOTE: the following macros are defined in generate.py: 
     14// 
     15//  MAX_PD : the maximum number of dispersity loops allowed for this model, 
     16//      which will be at most modelinfo.MAX_PD. 
     17//  NUM_PARS : the number of parameters in the parameter table 
     18//  NUM_VALUES : the number of values to skip at the start of the 
     19//      values array before you get to the dispersity values. 
     20//  PARAMETER_TABLE : list of parameter declarations used to create the 
     21//      ParameterTable type. 
     22//  KERNEL_NAME : model_Iq, model_Iqxy or model_Imagnetic.  This code is 
     23//      included three times, once for each kernel type. 
     24//  MAGNETIC : defined when the magnetic kernel is being instantiated 
     25//  NUM_MAGNETIC : the number of magnetic parameters 
     26//  MAGNETIC_PARS : a comma-separated list of indices to the sld 
     27//      parameters in the parameter table. 
     28//  CALL_VOLUME(table) : call the form volume function 
     29//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
     30//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
     31//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 
     32//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes 
     33//  INVALID(table) : test if the current point is feesible to calculate.  This 
     34//      will be defined in the kernel definition file. 
     35//  PROJECTION : equirectangular=1, sinusoidal=2 
     36//      see explore/jitter.py for definitions. 
     37 
    1438#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
    1539#define _PAR_BLOCK_ 
     
    1741typedef struct { 
    1842#if MAX_PD > 0 
    19     int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable 
    20     int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector 
     43    int32_t pd_par[MAX_PD];     // id of the nth dispersity variable 
     44    int32_t pd_length[MAX_PD];  // length of the nth dispersity weight vector 
    2145    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector 
    2246    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level 
     
    2549    int32_t num_weights;        // total length of the weights vector 
    2650    int32_t num_active;         // number of non-trivial pd loops 
    27     int32_t theta_par;          // id of spherical correction variable 
     51    int32_t theta_par;          // id of first orientation variable 
    2852} ProblemDetails; 
    2953 
     
    3862#endif // _PAR_BLOCK_ 
    3963 
    40  
    41 #if defined(MAGNETIC) && NUM_MAGNETIC>0 
     64#if defined(MAGNETIC) && NUM_MAGNETIC > 0 
     65// ===== Helper functions for magnetism ===== 
    4266 
    4367// Return value restricted between low and high 
     
    5175//     uu * (sld - m_sigma_x); 
    5276//     dd * (sld + m_sigma_x); 
    53 //     ud * (m_sigma_y + 1j*m_sigma_z); 
    54 //     du * (m_sigma_y - 1j*m_sigma_z); 
    55 static void set_spins(double in_spin, double out_spin, double spins[4]) 
     77//     ud * (m_sigma_y - 1j*m_sigma_z); 
     78//     du * (m_sigma_y + 1j*m_sigma_z); 
     79// weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag 
     80static void set_spin_weights(double in_spin, double out_spin, double spins[4]) 
    5681{ 
    5782  in_spin = clip(in_spin, 0.0, 1.0); 
    5883  out_spin = clip(out_spin, 0.0, 1.0); 
    5984  spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd 
    60   spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du 
    61   spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud 
     85  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real 
     86  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real 
    6287  spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu 
    63 } 
    64  
    65 static double mag_sld(double qx, double qy, double p, 
    66                        double mx, double my, double sld) 
    67 { 
     88  spins[4] = spins[1]; // du imag 
     89  spins[5] = spins[2]; // ud imag 
     90} 
     91 
     92// Compute the magnetic sld 
     93static double mag_sld( 
     94  const int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 
     95  const double qx, const double qy, 
     96  const double px, const double py, 
     97  const double sld, 
     98  const double mx, const double my, const double mz 
     99) 
     100{ 
     101  if (xs < 4) { 
    68102    const double perp = qy*mx - qx*my; 
    69     return sld + perp*p; 
    70 } 
    71  
    72 #endif // MAGNETIC 
     103    switch (xs) { 
     104      case 0: // uu => sld - D M_perpx 
     105          return sld - px*perp; 
     106      case 1: // ud real => -D M_perpy 
     107          return py*perp; 
     108      case 2: // du real => -D M_perpy 
     109          return py*perp; 
     110      case 3: // dd real => sld + D M_perpx 
     111          return sld + px*perp; 
     112    } 
     113  } else { 
     114    if (xs== 4) { 
     115      return -mz;  // ud imag => -D M_perpz 
     116    } else { // index == 5 
     117      return mz;   // du imag => D M_perpz 
     118    } 
     119  } 
     120} 
     121 
     122 
     123#endif 
     124 
     125// ===== Helper functions for orientation and jitter ===== 
     126 
     127// To change the definition of the angles, run explore/angles.py, which 
     128// uses sympy to generate the equations. 
     129 
     130#if !defined(_QAC_SECTION) && defined(CALL_IQ_AC) 
     131#define _QAC_SECTION 
     132 
     133typedef struct { 
     134    double R31, R32; 
     135} QACRotation; 
     136 
     137// Fill in the rotation matrix R from the view angles (theta, phi) and the 
     138// jitter angles (dtheta, dphi).  This matrix can be applied to all of the 
     139// (qx, qy) points in the image to produce R*[qx,qy]' = [qa,qc]' 
     140static void 
     141qac_rotation( 
     142    QACRotation *rotation, 
     143    double theta, double phi, 
     144    double dtheta, double dphi) 
     145{ 
     146    double sin_theta, cos_theta; 
     147    double sin_phi, cos_phi; 
     148 
     149    // reverse view matrix 
     150    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     151    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     152    const double V11 = cos_phi*cos_theta; 
     153    const double V12 = sin_phi*cos_theta; 
     154    const double V21 = -sin_phi; 
     155    const double V22 = cos_phi; 
     156    const double V31 = sin_theta*cos_phi; 
     157    const double V32 = sin_phi*sin_theta; 
     158 
     159    // reverse jitter matrix 
     160    SINCOS(dtheta*M_PI_180, sin_theta, cos_theta); 
     161    SINCOS(dphi*M_PI_180, sin_phi, cos_phi); 
     162    const double J31 = sin_theta; 
     163    const double J32 = -sin_phi*cos_theta; 
     164    const double J33 = cos_phi*cos_theta; 
     165 
     166    // reverse matrix 
     167    rotation->R31 = J31*V11 + J32*V21 + J33*V31; 
     168    rotation->R32 = J31*V12 + J32*V22 + J33*V32; 
     169} 
     170 
     171// Apply the rotation matrix returned from qac_rotation to the point (qx,qy), 
     172// returning R*[qx,qy]' = [qa,qc]' 
     173static double 
     174qac_apply( 
     175    QACRotation rotation, 
     176    double qx, double qy, 
     177    double *qa_out, double *qc_out) 
     178{ 
     179    const double dqc = rotation.R31*qx + rotation.R32*qy; 
     180    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
     181    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
     182 
     183    *qa_out = dqa; 
     184    *qc_out = dqc; 
     185} 
     186#endif // _QAC_SECTION 
     187 
     188#if !defined(_QABC_SECTION) && defined(CALL_IQ_ABC) 
     189#define _QABC_SECTION 
     190 
     191typedef struct { 
     192    double R11, R12; 
     193    double R21, R22; 
     194    double R31, R32; 
     195} QABCRotation; 
     196 
     197// Fill in the rotation matrix R from the view angles (theta, phi, psi) and the 
     198// jitter angles (dtheta, dphi, dpsi).  This matrix can be applied to all of the 
     199// (qx, qy) points in the image to produce R*[qx,qy]' = [qa,qb,qc]' 
     200static void 
     201qabc_rotation( 
     202    QABCRotation *rotation, 
     203    double theta, double phi, double psi, 
     204    double dtheta, double dphi, double dpsi) 
     205{ 
     206    double sin_theta, cos_theta; 
     207    double sin_phi, cos_phi; 
     208    double sin_psi, cos_psi; 
     209 
     210    // reverse view matrix 
     211    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     212    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     213    SINCOS(psi*M_PI_180, sin_psi, cos_psi); 
     214    const double V11 = -sin_phi*sin_psi + cos_phi*cos_psi*cos_theta; 
     215    const double V12 = sin_phi*cos_psi*cos_theta + sin_psi*cos_phi; 
     216    const double V21 = -sin_phi*cos_psi - sin_psi*cos_phi*cos_theta; 
     217    const double V22 = -sin_phi*sin_psi*cos_theta + cos_phi*cos_psi; 
     218    const double V31 = sin_theta*cos_phi; 
     219    const double V32 = sin_phi*sin_theta; 
     220 
     221    // reverse jitter matrix 
     222    SINCOS(dtheta*M_PI_180, sin_theta, cos_theta); 
     223    SINCOS(dphi*M_PI_180, sin_phi, cos_phi); 
     224    SINCOS(dpsi*M_PI_180, sin_psi, cos_psi); 
     225    const double J11 = cos_psi*cos_theta; 
     226    const double J12 = sin_phi*sin_theta*cos_psi + sin_psi*cos_phi; 
     227    const double J13 = sin_phi*sin_psi - sin_theta*cos_phi*cos_psi; 
     228    const double J21 = -sin_psi*cos_theta; 
     229    const double J22 = -sin_phi*sin_psi*sin_theta + cos_phi*cos_psi; 
     230    const double J23 = sin_phi*cos_psi + sin_psi*sin_theta*cos_phi; 
     231    const double J31 = sin_theta; 
     232    const double J32 = -sin_phi*cos_theta; 
     233    const double J33 = cos_phi*cos_theta; 
     234 
     235    // reverse matrix 
     236    rotation->R11 = J11*V11 + J12*V21 + J13*V31; 
     237    rotation->R12 = J11*V12 + J12*V22 + J13*V32; 
     238    rotation->R21 = J21*V11 + J22*V21 + J23*V31; 
     239    rotation->R22 = J21*V12 + J22*V22 + J23*V32; 
     240    rotation->R31 = J31*V11 + J32*V21 + J33*V31; 
     241    rotation->R32 = J31*V12 + J32*V22 + J33*V32; 
     242} 
     243 
     244// Apply the rotation matrix returned from qabc_rotation to the point (qx,qy), 
     245// returning R*[qx,qy]' = [qa,qb,qc]' 
     246static double 
     247qabc_apply( 
     248    QABCRotation rotation, 
     249    double qx, double qy, 
     250    double *qa_out, double *qb_out, double *qc_out) 
     251{ 
     252    *qa_out = rotation.R11*qx + rotation.R12*qy; 
     253    *qb_out = rotation.R21*qx + rotation.R22*qy; 
     254    *qc_out = rotation.R31*qx + rotation.R32*qy; 
     255} 
     256 
     257#endif // _QABC_SECTION 
     258 
     259 
     260// ==================== KERNEL CODE ======================== 
    73261 
    74262kernel 
    75263void KERNEL_NAME( 
    76264    int32_t nq,                 // number of q values 
    77     const int32_t pd_start,     // where we are in the polydispersity loop 
    78     const int32_t pd_stop,      // where we are stopping in the polydispersity loop 
     265    const int32_t pd_start,     // where we are in the dispersity loop 
     266    const int32_t pd_stop,      // where we are stopping in the dispersity loop 
    79267    global const ProblemDetails *details, 
    80268    global const double *values, 
    81269    global const double *q, // nq q values, with padding to boundary 
    82270    global double *result,  // nq+1 return values, again with padding 
    83     const double cutoff     // cutoff in the polydispersity weight product 
     271    const double cutoff     // cutoff in the dispersity weight product 
    84272    ) 
    85273{ 
    86   // Storage for the current parameter values.  These will be updated as we 
    87   // walk the polydispersity cube. 
     274#ifdef USE_OPENCL 
     275  // who we are and what element we are working with 
     276  const int q_index = get_global_id(0); 
     277  if (q_index >= nq) return; 
     278#else 
     279  // Define q_index here so that debugging statements can be written to work 
     280  // for both OpenCL and DLL using: 
     281  //    if (q_index == 0) {printf(...);} 
     282  int q_index = 0; 
     283#endif 
     284 
     285  // ** Fill in the local values table ** 
     286  // Storage for the current parameter values. 
     287  // These will be updated as we walk the dispersity mesh. 
    88288  ParameterBlock local_values; 
    89  
    90 #if defined(MAGNETIC) && NUM_MAGNETIC>0 
    91   // Location of the sld parameters in the parameter vector. 
    92   // These parameters are updated with the effective sld due to magnetism. 
    93   #if NUM_MAGNETIC > 3 
    94   const int32_t slds[] = { MAGNETIC_PARS }; 
    95   #endif 
    96  
    97   // TODO: could precompute these outside of the kernel. 
    98   // Interpret polarization cross section. 
    99   //     up_frac_i = values[NUM_PARS+2]; 
    100   //     up_frac_f = values[NUM_PARS+3]; 
    101   //     up_angle = values[NUM_PARS+4]; 
    102   double spins[4]; 
    103   double cos_mspin, sin_mspin; 
    104   set_spins(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
    105   SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
    106 #endif // MAGNETIC 
    107  
    108   // Fill in the initial variables 
    109289  //   values[0] is scale 
    110290  //   values[1] is background 
     
    114294  for (int i=0; i < NUM_PARS; i++) { 
    115295    local_values.vector[i] = values[2+i]; 
    116 //printf("p%d = %g\n",i, local_values.vector[i]); 
     296    //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]); 
    117297  } 
    118 //printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD); 
    119 //printf("start:%d  stop:%d\n", pd_start, pd_stop); 
    120  
    121   double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    122   if (pd_start == 0) { 
    123     #ifdef USE_OPENMP 
    124     #pragma omp parallel for 
    125     #endif 
    126     for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     298  //if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD); 
     299  //if (q_index==0) printf("start:%d  stop:%d\n", pd_start, pd_stop); 
     300 
     301  // ** Precompute magnatism values ** 
     302#if defined(MAGNETIC) && NUM_MAGNETIC>0 
     303  // Location of the sld parameters in the parameter vector. 
     304  // These parameters are updated with the effective sld due to magnetism. 
     305  const int32_t slds[] = { MAGNETIC_PARS }; 
     306 
     307  // Interpret polarization cross section. 
     308  //     up_frac_i = values[NUM_PARS+2]; 
     309  //     up_frac_f = values[NUM_PARS+3]; 
     310  //     up_angle = values[NUM_PARS+4]; 
     311  // TODO: could precompute more magnetism parameters before calling the kernel. 
     312  double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
     313  double cos_mspin, sin_mspin; 
     314  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
     315  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
     316#endif // MAGNETIC 
     317 
     318  // ** Fill in the initial results ** 
     319  // If pd_start is zero that means that we are starting a new calculation, 
     320  // and must initialize the result to zero.  Otherwise, we are restarting 
     321  // the calculation from somewhere in the middle of the dispersity mesh, 
     322  // and we update the value rather than reset it. Similarly for the 
     323  // normalization factor, which is stored as the final value in the 
     324  // results vector (one past the number of q values). 
     325  // 
     326  // The code differs slightly between opencl and dll since opencl is only 
     327  // seeing one q value (stored in the variable "this_result") while the dll 
     328  // version must loop over all q. 
     329  #ifdef USE_OPENCL 
     330    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     331    double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
     332  #else // !USE_OPENCL 
     333    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     334    if (pd_start == 0) { 
     335      #ifdef USE_OPENMP 
     336      #pragma omp parallel for 
     337      #endif 
     338      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     339    } 
     340    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
     341#endif // !USE_OPENCL 
     342 
     343 
     344// ====== macros to set up the parts of the loop ======= 
     345/* 
     346Based on the level of the loop, uses C preprocessor magic to construct 
     347level-specific looping variables, including these from loop level 3: 
     348 
     349  int n3 : length of loop for mesh level 3 
     350  int i3 : current position in the loop for level 3, which is calculated 
     351       from a combination of pd_start, pd_stride[3] and pd_length[3]. 
     352  int p3 : is the index into the parameter table for mesh level 3 
     353  double v3[] : pointer into dispersity array to values for loop 3 
     354  double w3[] : pointer into dispersity array to weights for loop 3 
     355  double weight3 : the product of weights from levels 3 and up, computed 
     356       as weight5*weight4*w3[i3].  Note that we need an outermost 
     357       value weight5 set to 1.0 for this to work properly. 
     358 
     359After expansion, the loop struction will look like the following: 
     360 
     361  // --- PD_INIT(4) --- 
     362  const int n4 = pd_length[4]; 
     363  const int p4 = pd_par[4]; 
     364  global const double *v4 = pd_value + pd_offset[4]; 
     365  global const double *w4 = pd_weight + pd_offset[4]; 
     366  int i4 = (pd_start/pd_stride[4])%n4;  // position in level 4 at pd_start 
     367 
     368  // --- PD_INIT(3) --- 
     369  const int n3 = pd_length[3]; 
     370  ... 
     371  int i3 = (pd_start/pd_stride[3])%n3;  // position in level 3 at pd_start 
     372 
     373  PD_INIT(2) 
     374  PD_INIT(1) 
     375  PD_INIT(0) 
     376 
     377  // --- PD_OUTERMOST_WEIGHT(5) --- 
     378  const double weight5 = 1.0; 
     379 
     380  // --- PD_OPEN(4,5) --- 
     381  while (i4 < n4) { 
     382    parameter[p4] = v4[i4];  // set the value for pd parameter 4 at this mesh point 
     383    const double weight4 = w4[i4] * weight5; 
     384 
     385    // from PD_OPEN(3,4) 
     386    while (i3 < n3) { 
     387      parameter[p3] = v3[i3];  // set the value for pd parameter 3 at this mesh point 
     388      const double weight3 = w3[i3] * weight4; 
     389 
     390      PD_OPEN(3,2) 
     391      PD_OPEN(2,1) 
     392      PD_OPEN(0,1) 
     393 
     394      // ... main loop body ... 
     395      APPLY_PROJECTION    // convert jitter values to spherical coords 
     396      BUILD_ROTATION      // construct the rotation matrix qxy => qabc 
     397      for each q 
     398          FETCH_Q         // set qx,qy from the q input vector 
     399          APPLY_ROTATION  // convert qx,qy to qa,qb,qc 
     400          CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...) 
     401 
     402      ++step;  // increment counter representing position in dispersity mesh 
     403 
     404      PD_CLOSE(0) 
     405      PD_CLOSE(1) 
     406      PD_CLOSE(2) 
     407 
     408      // --- PD_CLOSE(3) --- 
     409      if (step >= pd_stop) break; 
     410      ++i3; 
     411    } 
     412    i3 = 0; // reset loop counter for next round through the loop 
     413 
     414    // --- PD_CLOSE(4) --- 
     415    if (step >= pd_stop) break; 
     416    ++i4; 
    127417  } 
    128 //printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    129  
     418  i4 = 0; // reset loop counter even though no more rounds through the loop 
     419 
     420*/ 
     421 
     422 
     423// ** prepare inner loops ** 
     424 
     425// Depending on the shape type (radial, axial, triaxial), the variables 
     426// and calling parameters in the loop body will be slightly different. 
     427// Macros capture the differences in one spot so the rest of the code 
     428// is easier to read. The code below both declares variables for the 
     429// inner loop and defines the macros that use them. 
     430 
     431#if defined(CALL_IQ) 
     432  // unoriented 1D 
     433  double qk; 
     434  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     435  #define BUILD_ROTATION() do {} while(0) 
     436  #define APPLY_ROTATION() do {} while(0) 
     437  #define CALL_KERNEL() CALL_IQ(qk, local_values.table) 
     438 
     439#elif defined(CALL_IQ_A) 
     440  // unoriented 2D 
     441  double qx, qy; 
     442  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     443  #define BUILD_ROTATION() do {} while(0) 
     444  #define APPLY_ROTATION() do {} while(0) 
     445  #define CALL_KERNEL() CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table) 
     446 
     447#elif defined(CALL_IQ_AC) 
     448  // oriented symmetric 2D 
     449  double qx, qy; 
     450  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     451  double qa, qc; 
     452  QACRotation rotation; 
     453  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code. 
     454  #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi); 
     455  #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc) 
     456  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table) 
     457 
     458#elif defined(CALL_IQ_ABC) 
     459  // oriented asymmetric 2D 
     460  double qx, qy; 
     461  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     462  double qa, qb, qc; 
     463  QABCRotation rotation; 
     464  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code. 
     465  // psi and dpsi are only for IQ_ABC, so they are processed here. 
     466  const double psi = values[details->theta_par+4]; 
     467  local_values.table.psi = 0.; 
     468  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi) 
     469  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
     470  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     471#endif 
     472 
     473// Doing jitter projection code outside the previous if block so that we don't 
     474// need to repeat the identical logic in the IQ_AC and IQ_ABC branches.  This 
     475// will become more important if we implement more projections, or more 
     476// complicated projections. 
     477#if defined(CALL_IQ) || defined(CALL_IQ_A) 
     478  #define APPLY_PROJECTION() const double weight=weight0 
     479#else // !spherosymmetric projection 
     480  // Grab the "view" angles (theta, phi, psi) from the initial parameter table. 
     481  const double theta = values[details->theta_par+2]; 
     482  const double phi = values[details->theta_par+3]; 
     483  // Make sure jitter angle defaults to zero if there is no jitter distribution 
     484  local_values.table.theta = 0.; 
     485  local_values.table.phi = 0.; 
     486  // The "jitter" angles (dtheta, dphi, dpsi) are stored with the 
     487  // dispersity values and copied to the local parameter table as 
     488  // we go through the mesh. 
     489  double dtheta, dphi, weight; 
     490  #if PROJECTION == 1 
     491    #define APPLY_PROJECTION() do { \ 
     492      dtheta = local_values.table.theta; \ 
     493      dphi = local_values.table.phi; \ 
     494      weight = fabs(cos(dtheta*M_PI_180)) * weight0; \ 
     495    } while (0) 
     496  #elif PROJECTION == 2 
     497    #define APPLY_PROJECTION() do { \ 
     498      dtheta = local_values.table.theta; \ 
     499      dphi = local_values.table.phi; \ 
     500      weight = weight0; \ 
     501      if (dtheta != 90.0) dphi /= cos(dtheta*M_PI_180); \ 
     502      else if (dphi != 0.0) weight = 0.; \ 
     503      if (fabs(dphi) >= 180.) weight = 0.; \ 
     504    } while (0) 
     505  #endif 
     506#endif // !spherosymmetric projection 
     507 
     508// ** define looping macros ** 
     509 
     510// Define looping variables 
     511#define PD_INIT(_LOOP) \ 
     512  const int n##_LOOP = details->pd_length[_LOOP]; \ 
     513  const int p##_LOOP = details->pd_par[_LOOP]; \ 
     514  global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \ 
     515  global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \ 
     516  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP; 
     517 
     518// Jump into the middle of the dispersity loop 
     519#define PD_OPEN(_LOOP,_OUTER) \ 
     520  while (i##_LOOP < n##_LOOP) { \ 
     521    local_values.vector[p##_LOOP] = v##_LOOP[i##_LOOP]; \ 
     522    const double weight##_LOOP = w##_LOOP[i##_LOOP] * weight##_OUTER; 
     523 
     524// create the variable "weight#=1.0" where # is the outermost level+1 (=MAX_PD). 
     525#define _PD_OUTERMOST_WEIGHT(_n) const double weight##_n = 1.0; 
     526#define PD_OUTERMOST_WEIGHT(_n) _PD_OUTERMOST_WEIGHT(_n) 
     527 
     528// Close out the loop 
     529#define PD_CLOSE(_LOOP) \ 
     530    if (step >= pd_stop) break; \ 
     531    ++i##_LOOP; \ 
     532  } \ 
     533  i##_LOOP = 0; 
     534 
     535// ====== construct the loops ======= 
     536 
     537// Pointers to the start of the dispersity and weight vectors, if needed. 
    130538#if MAX_PD>0 
    131539  global const double *pd_value = values + NUM_VALUES; 
     
    133541#endif 
    134542 
    135   // Jump into the middle of the polydispersity loop 
     543// The variable "step" is the current position in the dispersity loop. 
     544// It will be incremented each time a new point in the mesh is accumulated, 
     545// and used to test whether we have reached pd_stop. 
     546int step = pd_start; 
     547 
     548// *** define loops for each of 0, 1, 2, ..., modelinfo.MAX_PD-1 *** 
     549 
     550// define looping variables 
    136551#if MAX_PD>4 
    137   int n4=details->pd_length[4]; 
    138   int i4=(pd_start/details->pd_stride[4])%n4; 
    139   const int p4=details->pd_par[4]; 
    140   global const double *v4 = pd_value + details->pd_offset[4]; 
    141   global const double *w4 = pd_weight + details->pd_offset[4]; 
     552  PD_INIT(4) 
    142553#endif 
    143554#if MAX_PD>3 
    144   int n3=details->pd_length[3]; 
    145   int i3=(pd_start/details->pd_stride[3])%n3; 
    146   const int p3=details->pd_par[3]; 
    147   global const double *v3 = pd_value + details->pd_offset[3]; 
    148   global const double *w3 = pd_weight + details->pd_offset[3]; 
    149 //printf("offset %d: %d %d\n", 3, details->pd_offset[3], NUM_VALUES); 
     555  PD_INIT(3) 
    150556#endif 
    151557#if MAX_PD>2 
    152   int n2=details->pd_length[2]; 
    153   int i2=(pd_start/details->pd_stride[2])%n2; 
    154   const int p2=details->pd_par[2]; 
    155   global const double *v2 = pd_value + details->pd_offset[2]; 
    156   global const double *w2 = pd_weight + details->pd_offset[2]; 
     558  PD_INIT(2) 
    157559#endif 
    158560#if MAX_PD>1 
    159   int n1=details->pd_length[1]; 
    160   int i1=(pd_start/details->pd_stride[1])%n1; 
    161   const int p1=details->pd_par[1]; 
    162   global const double *v1 = pd_value + details->pd_offset[1]; 
    163   global const double *w1 = pd_weight + details->pd_offset[1]; 
     561  PD_INIT(1) 
    164562#endif 
    165563#if MAX_PD>0 
    166   int n0=details->pd_length[0]; 
    167   int i0=(pd_start/details->pd_stride[0])%n0; 
    168   const int p0=details->pd_par[0]; 
    169   global const double *v0 = pd_value + details->pd_offset[0]; 
    170   global const double *w0 = pd_weight + details->pd_offset[0]; 
    171 //printf("w0:%p, values:%p, diff:%ld, %d\n",w0,values,(w0-values), NUM_VALUES); 
    172 #endif 
    173  
    174  
     564  PD_INIT(0) 
     565#endif 
     566 
     567// open nested loops 
     568PD_OUTERMOST_WEIGHT(MAX_PD) 
     569#if MAX_PD>4 
     570  PD_OPEN(4,5) 
     571#endif 
     572#if MAX_PD>3 
     573  PD_OPEN(3,4) 
     574#endif 
     575#if MAX_PD>2 
     576  PD_OPEN(2,3) 
     577#endif 
     578#if MAX_PD>1 
     579  PD_OPEN(1,2) 
     580#endif 
    175581#if MAX_PD>0 
    176   const int theta_par = details->theta_par; 
    177   const int fast_theta = (theta_par == p0); 
    178   const int slow_theta = (theta_par >= 0 && !fast_theta); 
    179   double spherical_correction = 1.0; 
    180 #else 
    181   // Note: if not polydisperse the weights cancel and we don't need the 
    182   // spherical correction. 
    183   const double spherical_correction = 1.0; 
    184 #endif 
    185  
    186   int step = pd_start; 
    187  
    188 #if MAX_PD>4 
    189   const double weight5 = 1.0; 
    190   while (i4 < n4) { 
    191     local_values.vector[p4] = v4[i4]; 
    192     double weight4 = w4[i4] * weight5; 
    193 //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); 
    194 #elif MAX_PD>3 
    195     const double weight4 = 1.0; 
    196 #endif 
    197 #if MAX_PD>3 
    198   while (i3 < n3) { 
    199     local_values.vector[p3] = v3[i3]; 
    200     double weight3 = w3[i3] * weight4; 
    201 //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); 
    202 #elif MAX_PD>2 
    203     const double weight3 = 1.0; 
    204 #endif 
    205 #if MAX_PD>2 
    206   while (i2 < n2) { 
    207     local_values.vector[p2] = v2[i2]; 
    208     double weight2 = w2[i2] * weight3; 
    209 //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); 
    210 #elif MAX_PD>1 
    211     const double weight2 = 1.0; 
    212 #endif 
    213 #if MAX_PD>1 
    214   while (i1 < n1) { 
    215     local_values.vector[p1] = v1[i1]; 
    216     double weight1 = w1[i1] * weight2; 
    217 //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); 
    218 #elif MAX_PD>0 
    219     const double weight1 = 1.0; 
    220 #endif 
    221 #if MAX_PD>0 
    222   if (slow_theta) { // Theta is not in inner loop 
    223     spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6); 
    224   } 
    225   while(i0 < n0) { 
    226     local_values.vector[p0] = v0[i0]; 
    227     double weight0 = w0[i0] * weight1; 
    228 //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); 
    229     if (fast_theta) { // Theta is in inner loop 
    230       spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6); 
    231     } 
    232 #else 
    233     const double weight0 = 1.0; 
    234 #endif 
    235  
    236 //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"); 
    237 //printf("sphcor: %g\n", spherical_correction); 
    238  
    239     #ifdef INVALID 
    240     if (!INVALID(local_values.table)) 
    241     #endif 
    242     { 
    243       // Accumulate I(q) 
    244       // Note: weight==0 must always be excluded 
    245       if (weight0 > cutoff) { 
    246         // spherical correction is set at a minimum of 1e-6, otherwise there 
    247         // would be problems looking at models with theta=90. 
    248         const double weight = weight0 * spherical_correction; 
    249         pd_norm += weight * CALL_VOLUME(local_values.table); 
    250  
    251         #ifdef USE_OPENMP 
    252         #pragma omp parallel for 
    253         #endif 
    254         for (int q_index=0; q_index<nq; q_index++) { 
    255 #if defined(MAGNETIC) && NUM_MAGNETIC > 0 
    256           const double qx = q[2*q_index]; 
    257           const double qy = q[2*q_index+1]; 
     582  PD_OPEN(0,1) 
     583#endif 
     584 
     585//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");} 
     586 
     587  // ====== loop body ======= 
     588  #ifdef INVALID 
     589  if (!INVALID(local_values.table)) 
     590  #endif 
     591  { 
     592     APPLY_PROJECTION(); 
     593 
     594    // Accumulate I(q) 
     595    // Note: weight==0 must always be excluded 
     596    if (weight > cutoff) { 
     597      pd_norm += weight * CALL_VOLUME(local_values.table); 
     598      BUILD_ROTATION(); 
     599 
     600#ifndef USE_OPENCL 
     601      // DLL needs to explicitly loop over the q values. 
     602      #ifdef USE_OPENMP 
     603      #pragma omp parallel for 
     604      #endif 
     605      for (q_index=0; q_index<nq; q_index++) 
     606#endif // !USE_OPENCL 
     607      { 
     608 
     609        FETCH_Q(); 
     610        APPLY_ROTATION(); 
     611 
     612        // ======= COMPUTE SCATTERING ========== 
     613        #if defined(MAGNETIC) && NUM_MAGNETIC > 0 
     614          // Compute the scattering from the magnetic cross sections. 
     615          double scattering = 0.0; 
    258616          const double qsq = qx*qx + qy*qy; 
    259  
    260           // Constant across orientation, polydispersity for given qx, qy 
    261           double scattering = 0.0; 
    262           // TODO: what is the magnetic scattering at q=0 
    263617          if (qsq > 1.e-16) { 
    264             double p[4];  // dd, du, ud, uu 
    265             p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq; 
    266             p[3] = -p[0]; 
    267             p[1] = p[2] = (qy*sin_mspin - qx*cos_mspin)/qsq; 
    268  
    269             for (int index=0; index<4; index++) { 
    270               const double xs = spins[index]; 
    271               if (xs > 1.e-8) { 
    272                 const int spin_flip = (index==1) || (index==2); 
    273                 const double pk = p[index]; 
    274                 for (int axis=0; axis<=spin_flip; axis++) { 
    275                   #define M1 NUM_PARS+5 
    276                   #define M2 NUM_PARS+8 
    277                   #define M3 NUM_PARS+13 
    278                   #define SLD(_M_offset, _sld_offset) \ 
    279                       local_values.vector[_sld_offset] = xs * (axis \ 
    280                       ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \ 
    281                       : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \ 
    282                                 (spin_flip ? 0.0 : values[_sld_offset+2]))) 
    283                   #if NUM_MAGNETIC==1 
    284                       SLD(M1, MAGNETIC_PAR1); 
    285                   #elif NUM_MAGNETIC==2 
    286                       SLD(M1, MAGNETIC_PAR1); 
    287                       SLD(M2, MAGNETIC_PAR2); 
    288                   #elif NUM_MAGNETIC==3 
    289                       SLD(M1, MAGNETIC_PAR1); 
    290                       SLD(M2, MAGNETIC_PAR2); 
    291                       SLD(M3, MAGNETIC_PAR3); 
    292                   #else 
    293                   for (int sk=0; sk<NUM_MAGNETIC; sk++) { 
    294                       SLD(M1+3*sk, slds[sk]); 
    295                   } 
    296                   #endif 
    297                   scattering += CALL_IQ(q, q_index, local_values.table); 
     618            // TODO: what is the magnetic scattering at q=0 
     619            const double px = (qy*cos_mspin + qx*sin_mspin)/qsq; 
     620            const double py = (qy*sin_mspin - qx*cos_mspin)/qsq; 
     621 
     622            // loop over uu, ud real, du real, dd, ud imag, du imag 
     623            for (int xs=0; xs<6; xs++) { 
     624              const double xs_weight = spins[xs]; 
     625              if (xs_weight > 1.e-8) { 
     626                // Since the cross section weight is significant, set the slds 
     627                // to the effective slds for this cross section, call the 
     628                // kernel, and add according to weight. 
     629                for (int sk=0; sk<NUM_MAGNETIC; sk++) { 
     630                  const int32_t mag_index = NUM_PARS+5 + 3*sk; 
     631                  const int32_t sld_index = slds[sk]; 
     632                  const double mx = values[mag_index]; 
     633                  const double my = values[mag_index+1]; 
     634                  const double mz = values[mag_index+2]; 
     635                  local_values.vector[sld_index] = 
     636                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz); 
    298637                } 
     638                scattering += xs_weight * CALL_KERNEL(); 
    299639              } 
    300640            } 
    301641          } 
    302 #else  // !MAGNETIC 
    303           const double scattering = CALL_IQ(q, q_index, local_values.table); 
    304 #endif // !MAGNETIC 
    305 //printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0); 
     642        #else  // !MAGNETIC 
     643          const double scattering = CALL_KERNEL(); 
     644        #endif // !MAGNETIC 
     645//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 
     646 
     647        #ifdef USE_OPENCL 
     648          this_result += weight * scattering; 
     649        #else // !USE_OPENCL 
    306650          result[q_index] += weight * scattering; 
    307         } 
     651        #endif // !USE_OPENCL 
    308652      } 
    309653    } 
    310     ++step; 
     654  } 
     655 
     656// close nested loops 
     657++step; 
    311658#if MAX_PD>0 
    312     if (step >= pd_stop) break; 
    313     ++i0; 
    314   } 
    315   i0 = 0; 
     659  PD_CLOSE(0) 
    316660#endif 
    317661#if MAX_PD>1 
    318     if (step >= pd_stop) break; 
    319     ++i1; 
    320   } 
    321   i1 = 0; 
     662  PD_CLOSE(1) 
    322663#endif 
    323664#if MAX_PD>2 
    324     if (step >= pd_stop) break; 
    325     ++i2; 
    326   } 
    327   i2 = 0; 
     665  PD_CLOSE(2) 
    328666#endif 
    329667#if MAX_PD>3 
    330     if (step >= pd_stop) break; 
    331     ++i3; 
    332   } 
    333   i3 = 0; 
     668  PD_CLOSE(3) 
    334669#endif 
    335670#if MAX_PD>4 
    336     if (step >= pd_stop) break; 
    337     ++i4; 
    338   } 
    339   i4 = 0; 
    340 #endif 
    341  
     671  PD_CLOSE(4) 
     672#endif 
     673 
     674// Remember the current result and the updated norm. 
     675#ifdef USE_OPENCL 
     676  result[q_index] = this_result; 
     677  if (q_index == 0) result[nq] = pd_norm; 
     678//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
     679#else // !USE_OPENCL 
     680  result[nq] = pd_norm; 
    342681//printf("res: %g/%g\n", result[0], pd_norm); 
    343   // Remember the updated norm. 
    344   result[nq] = pd_norm; 
    345 } 
     682#endif // !USE_OPENCL 
     683 
     684// ** clear the macros in preparation for the next kernel ** 
     685#undef PD_INIT 
     686#undef PD_OPEN 
     687#undef PD_CLOSE 
     688#undef FETCH_Q 
     689#undef APPLY_PROJECTION 
     690#undef BUILD_ROTATION 
     691#undef APPLY_ROTATION 
     692#undef CALL_KERNEL 
     693} 
  • sasmodels/kernelpy.py

    r3b32f8d r8698a0d  
    113113 
    114114        partable = model_info.parameters 
    115         kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 
    116                              else partable.iq_parameters) 
     115        #kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 
     116        #                     else partable.iq_parameters) 
     117        kernel_parameters = partable.iq_parameters 
    117118        volume_parameters = partable.form_volume_parameters 
    118119 
     
    201202 
    202203    pd_norm = 0.0 
    203     spherical_correction = 1.0 
    204204    partial_weight = np.NaN 
    205205    weight = np.NaN 
    206206 
    207207    p0_par = call_details.pd_par[0] 
    208     p0_is_theta = (p0_par == call_details.theta_par) 
    209208    p0_length = call_details.pd_length[0] 
    210209    p0_index = p0_length 
     
    223222            parameters[pd_par] = pd_value[pd_offset+pd_index] 
    224223            partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
    225             if call_details.theta_par >= 0: 
    226                 cor = sin(pi / 180 * parameters[call_details.theta_par]) 
    227                 spherical_correction = max(abs(cor), 1e-6) 
    228224            p0_index = loop_index%p0_length 
    229225 
    230226        weight = partial_weight * pd_weight[p0_offset + p0_index] 
    231227        parameters[p0_par] = pd_value[p0_offset + p0_index] 
    232         if p0_is_theta: 
    233             cor = cos(pi/180 * parameters[p0_par]) 
    234             spherical_correction = max(abs(cor), 1e-6) 
    235228        p0_index += 1 
    236229        if weight > cutoff: 
     
    244237 
    245238            # update value and norm 
    246             weight *= spherical_correction 
    247239            total += weight * Iq 
    248240            pd_norm += weight * form_volume() 
  • sasmodels/model_test.py

    r65314f7 r20fe0cd  
    8686    suite = unittest.TestSuite() 
    8787 
    88     if models[0] == 'all': 
     88    if models[0] in core.KINDS: 
    8989        skip = models[1:] 
    90         models = list_models() 
     90        models = list_models(models[0]) 
    9191    else: 
    9292        skip = [] 
     
    202202                ] 
    203203            tests = smoke_tests 
     204            #tests = [] 
    204205            if self.info.tests is not None: 
    205206                tests += self.info.tests 
  • sasmodels/modelinfo.py

    ra85a569 r6aee3ab  
    3030    TestCondition = Tuple[ParameterSetUser, TestInput, TestValue] 
    3131 
    32 MAX_PD = 4 #: Maximum number of simultaneously polydisperse parameters 
     32# If MAX_PD changes, need to change the loop macros in kernel_iq.c 
     33MAX_PD = 5 #: Maximum number of simultaneously polydisperse parameters 
    3334 
    3435# assumptions about common parameters exist throughout the code, such as: 
     
    382383      with vector parameter p sent as p[]. 
    383384 
    384     * *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 
     385    * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 
    385386      function, with vector parameter p sent as p[]. 
    386387 
     
    420421        # type: (List[Parameter]) -> None 
    421422        self.kernel_parameters = parameters 
     423        self._check_angles() 
    422424        self._set_vector_lengths() 
    423425 
     
    438440        self.iq_parameters = [p for p in self.kernel_parameters 
    439441                              if p.type not in ('orientation', 'magnetic')] 
    440         self.iqxy_parameters = [p for p in self.kernel_parameters 
    441                                 if p.type != 'magnetic'] 
     442        # note: orientation no longer sent to Iqxy, so its the same as 
     443        #self.iqxy_parameters = [p for p in self.kernel_parameters 
     444        #                        if p.type != 'magnetic'] 
    442445        self.form_volume_parameters = [p for p in self.kernel_parameters 
    443446                                       if p.type == 'volume'] 
     
    461464        self.has_2d = any(p.type in ('orientation', 'magnetic') 
    462465                          for p in self.kernel_parameters) 
     466        self.is_asymmetric = any(p.name == 'psi' for p in self.kernel_parameters) 
    463467        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    464468                                if p.id.startswith('M0:')] 
     
    467471                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    468472        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
     473 
     474    def _check_angles(self): 
     475        theta = phi = psi = -1 
     476        for k, p in enumerate(self.kernel_parameters): 
     477            if p.name == 'theta': 
     478                theta = k 
     479                if p.type != 'orientation': 
     480                    raise TypeError("theta must be an orientation parameter") 
     481            elif p.name == 'phi': 
     482                phi = k 
     483                if p.type != 'orientation': 
     484                    raise TypeError("phi must be an orientation parameter") 
     485            elif p.name == 'psi': 
     486                psi = k 
     487                if p.type != 'orientation': 
     488                    raise TypeError("psi must be an orientation parameter") 
     489        if theta >= 0 and phi >= 0: 
     490            if phi != theta+1: 
     491                raise TypeError("phi must follow theta") 
     492            if psi >= 0 and psi != phi+1: 
     493                raise TypeError("psi must follow phi") 
     494        elif theta >= 0 or phi >= 0 or psi >= 0: 
     495            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
    469496 
    470497    def __getitem__(self, key): 
     
    476503            raise KeyError("unknown parameter %r"%key) 
    477504        return par 
     505 
     506    def __contains__(self, key): 
     507        for par in self.call_parameters: 
     508            if par.name == key: 
     509                return True 
     510        else: 
     511            return False 
    478512 
    479513    def _set_vector_lengths(self): 
  • sasmodels/models/barbell.c

    r592343f rbecded3  
    1 double form_volume(double radius_bell, double radius, double length); 
    2 double Iq(double q, double sld, double solvent_sld, 
    3         double radius_bell, double radius, double length); 
    4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    5         double radius_bell, double radius, double length, 
    6         double theta, double phi); 
    7  
    81#define INVALID(v) (v.radius_bell < v.radius) 
    92 
    103//barbell kernel - same as dumbell 
    114static double 
    12 _bell_kernel(double q, double h, double radius_bell, 
    13              double half_length, double sin_alpha, double cos_alpha) 
     5_bell_kernel(double qab, double qc, double h, double radius_bell, 
     6             double half_length) 
    147{ 
    158    // translate a point in [-1,1] to a point in [lower,upper] 
     
    2619    //    m = q R cos(alpha) 
    2720    //    b = q(L/2-h) cos(alpha) 
    28     const double m = q*radius_bell*cos_alpha; // cos argument slope 
    29     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    30     const double qrst = q*radius_bell*sin_alpha; // Q*R*sin(theta) 
     21    const double m = radius_bell*qc; // cos argument slope 
     22    const double b = (half_length-h)*qc; // cos argument intercept 
     23    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    3124    double total = 0.0; 
    3225    for (int i = 0; i < 76; i++){ 
    3326        const double t = Gauss76Z[i]*zm + zb; 
    3427        const double radical = 1.0 - t*t; 
    35         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     28        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    3629        const double Fq = cos(m*t + b) * radical * bj; 
    3730        total += Gauss76Wt[i] * Fq; 
     
    4437 
    4538static double 
    46 _fq(double q, double h, 
    47     double radius_bell, double radius, double half_length, 
    48     double sin_alpha, double cos_alpha) 
     39_fq(double qab, double qc, double h, 
     40    double radius_bell, double radius, double half_length) 
    4941{ 
    50     const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 
    51     const double bj = sas_2J1x_x(q*radius*sin_alpha); 
    52     const double si = sas_sinx_x(q*half_length*cos_alpha); 
     42    const double bell_fq = _bell_kernel(qab, qc, h, radius_bell, half_length); 
     43    const double bj = sas_2J1x_x(radius*qab); 
     44    const double si = sas_sinx_x(half_length*qc); 
    5345    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5446    const double Aq = bell_fq + cyl_fq; 
     
    5648} 
    5749 
    58  
    59 double form_volume(double radius_bell, 
    60         double radius, 
    61         double length) 
     50static double 
     51form_volume(double radius_bell, 
     52    double radius, 
     53    double length) 
    6254{ 
    6355    // bell radius should never be less than radius when this is called 
     
    7062} 
    7163 
    72 double Iq(double q, double sld, double solvent_sld, 
    73           double radius_bell, double radius, double length) 
     64static double 
     65Iq(double q, double sld, double solvent_sld, 
     66    double radius_bell, double radius, double length) 
    7467{ 
    7568    const double h = -sqrt(radius_bell*radius_bell - radius*radius); 
     
    8477        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8578        SINCOS(alpha, sin_alpha, cos_alpha); 
    86         const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 
     79        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    8780        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8881    } 
     
    9689 
    9790 
    98 double Iqxy(double qx, double qy, 
    99         double sld, double solvent_sld, 
    100         double radius_bell, double radius, double length, 
    101         double theta, double phi) 
     91static double 
     92Iqxy(double qab, double qc, 
     93    double sld, double solvent_sld, 
     94    double radius_bell, double radius, double length) 
    10295{ 
    103     double q, sin_alpha, cos_alpha; 
    104     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    105  
    10696    const double h = -sqrt(square(radius_bell) - square(radius)); 
    107     const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 
     97    const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length); 
    10898 
    10999    // Multiply by contrast^2 and convert to cm-1 
  • sasmodels/models/bcc_paracrystal.c

    r50beefe rea60e08  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    // Equations from Matsuoka 26-27-28, multiplied by |q| 
     5    const double a1 = (-qa + qb + qc)/2.0; 
     6    const double a2 = (+qa - qb + qc)/2.0; 
     7    const double a3 = (+qa + qb - qc)/2.0; 
    68 
    7 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _BCCeval(double Theta, double Phi, double temp1, double temp3); 
    9 double _sphereform(double q, double radius, double sld, double solvent_sld); 
     9#if 1 
     10    // Matsuoka 29-30-31 
     11    //     Z_k numerator: 1 - exp(a)^2 
     12    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a) 
     13    // Rewriting numerator 
     14    //         => -(exp(2a) - 1) 
     15    //         => -expm1(2a) 
     16    // Rewriting denominator 
     17    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     18    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
     19    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     20    const double exp_arg = exp(arg); 
     21    const double Zq = -cube(expm1(2.0*arg)) 
     22        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     24          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     25 
     26#elif 0 
     27    // ** Alternate form, which perhaps is more approachable 
     28    //     Z_k numerator   => -[(exp(2a) - 1) / 2.exp(a)] 2.exp(a) 
     29    //                     => -[sinh(a)] exp(a) 
     30    //     Z_k denominator => [(exp(2a) + 1) / 2.exp(a) - cos(d a_k)] 2.exp(a) 
     31    //                     => [cosh(a) - cos(d a_k)] 2.exp(a) 
     32    //     => Z_k = -sinh(a) / [cosh(a) - cos(d a_k)] 
     33    //            = sinh(-a) / [cosh(-a) - cos(d a_k)] 
     34    // 
     35    // One more step leads to the form in sasview 3.x for 2d models 
     36    //            = tanh(-a) / [1 - cos(d a_k)/cosh(-a)] 
     37    // 
     38    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     39    const double sinh_qd = sinh(arg); 
     40    const double cosh_qd = cosh(arg); 
     41    const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1)) 
     42                    * sinh_qd/(cosh_qd - cos(dnn*a2)) 
     43                    * sinh_qd/(cosh_qd - cos(dnn*a3)); 
     44#else 
     45    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     46    const double tanh_qd = tanh(arg); 
     47    const double cosh_qd = cosh(arg); 
     48    const double Zq = tanh_qd/(1.0 - cos(dnn*a1)/cosh_qd) 
     49                    * tanh_qd/(1.0 - cos(dnn*a2)/cosh_qd) 
     50                    * tanh_qd/(1.0 - cos(dnn*a3)/cosh_qd); 
     51#endif 
     52 
     53    return Zq; 
     54} 
    1055 
    1156 
    12 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    13  
    14         const double Da = d_factor*dnn; 
    15         const double temp1 = q*q*Da*Da; 
    16         const double temp3 = q*dnn; 
    17  
    18         double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    19         return(retVal); 
     57// occupied volume fraction calculated from lattice symmetry and sphere radius 
     58static double 
     59bcc_volume_fraction(double radius, double dnn) 
     60{ 
     61    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
    2062} 
    2163 
    22 double _BCCeval(double Theta, double Phi, double temp1, double temp3) { 
    23  
    24         double result; 
    25         double sin_theta,cos_theta,sin_phi,cos_phi; 
    26         SINCOS(Theta, sin_theta, cos_theta); 
    27         SINCOS(Phi, sin_phi, cos_phi); 
    28  
    29         const double temp6 =  sin_theta; 
    30         const double temp7 =  sin_theta*cos_phi + sin_theta*sin_phi + cos_theta; 
    31         const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 
    32         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 
    33  
    34         const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 
    35         result = cube(1.0 - (temp10*temp10))*temp6 
    36             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    38               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    39  
    40         return (result); 
    41 } 
    42  
    43 double form_volume(double radius){ 
     64static double 
     65form_volume(double radius) 
     66{ 
    4467    return sphere_volume(radius); 
    4568} 
    4669 
    4770 
    48 double Iq(double q, double dnn, 
    49   double d_factor, double radius, 
    50   double sld, double solvent_sld){ 
     71static double Iq(double q, double dnn, 
     72    double d_factor, double radius, 
     73    double sld, double solvent_sld) 
     74{ 
     75    // translate a point in [-1,1] to a point in [0, 2 pi] 
     76    const double phi_m = M_PI; 
     77    const double phi_b = M_PI; 
     78    // translate a point in [-1,1] to a point in [0, pi] 
     79    const double theta_m = M_PI_2; 
     80    const double theta_b = M_PI_2; 
    5181 
    52         //Volume fraction calculated from lattice symmetry and sphere radius 
    53         const double s1 = dnn/sqrt(0.75); 
    54         const double latticescale = 2.0*sphere_volume(radius/s1); 
    55  
    56     const double va = 0.0; 
    57     const double vb = 2.0*M_PI; 
    58     const double vaj = 0.0; 
    59     const double vbj = M_PI; 
    60  
    61     double summ = 0.0; 
    62     double answer = 0.0; 
    63         for(int i=0; i<150; i++) { 
    64                 //setup inner integral over the ellipsoidal cross-section 
    65                 double summj=0.0; 
    66                 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi 
    67                 for(int j=0;j<150;j++) { 
    68                         //20 gauss points for the inner integral 
    69                         double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta 
    70                         double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi); 
    71                         summj += yyy; 
    72                 } 
    73                 //now calculate the value of the inner integral 
    74                 double answer = (vbj-vaj)/2.0*summj; 
    75  
    76                 //now calculate outer integral 
    77                 summ = summ+(Gauss150Wt[i] * answer); 
    78         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    79  
    80         answer = (vb-va)/2.0*summ; 
    81         answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 
    82  
    83     return answer; 
     82    double outer_sum = 0.0; 
     83    for(int i=0; i<150; i++) { 
     84        double inner_sum = 0.0; 
     85        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     86        double sin_theta, cos_theta; 
     87        SINCOS(theta, sin_theta, cos_theta); 
     88        const double qc = q*cos_theta; 
     89        const double qab = q*sin_theta; 
     90        for(int j=0;j<150;j++) { 
     91            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     92            double sin_phi, cos_phi; 
     93            SINCOS(phi, sin_phi, cos_phi); 
     94            const double qa = qab*cos_phi; 
     95            const double qb = qab*sin_phi; 
     96            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     97            inner_sum += Gauss150Wt[j] * form; 
     98        } 
     99        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     100        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     101    } 
     102    outer_sum *= theta_m; 
     103    const double Zq = outer_sum/(4.0*M_PI); 
     104    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     105    return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
    84106} 
    85107 
    86108 
    87 double Iqxy(double qx, double qy, 
     109static double Iqxy(double qa, double qb, double qc, 
    88110    double dnn, double d_factor, double radius, 
    89     double sld, double solvent_sld, 
    90     double theta, double phi, double psi) 
     111    double sld, double solvent_sld) 
    91112{ 
    92     double q, zhat, yhat, xhat; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    94  
    95     const double a1 = +xhat - zhat + yhat; 
    96     const double a2 = +xhat + zhat - yhat; 
    97     const double a3 = -xhat + zhat + yhat; 
    98  
    99     const double qd = 0.5*q*dnn; 
    100     const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
    101     const double tanh_qd = tanh(arg); 
    102     const double cosh_qd = cosh(arg); 
    103     const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 
    104                     * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 
    105                     * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 
    106  
    107     const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
    108     //the occupied volume of the lattice 
    109     const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
    110     return lattice_scale * Fq; 
     113    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
     114    const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     115    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     116    return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
    111117} 
  • sasmodels/models/bcc_paracrystal.py

    r8f04da4 r1f159bd  
    6565    \end{array} 
    6666 
    67 **NB**: The calculation of $Z(q)$ is a double numerical integral that must 
    68 be carried out with a high density of points to properly capture the sharp 
    69 peaks of the paracrystalline scattering. So be warned that the calculation 
    70 is SLOW. Go get some coffee. Fitting of any experimental data must be 
    71 resolution smeared for any meaningful fit. This makes a triple integral. 
    72 Very, very slow. Go get lunch! 
     67.. note:: 
     68 
     69  The calculation of $Z(q)$ is a double numerical integral that 
     70  must be carried out with a high density of points to properly capture 
     71  the sharp peaks of the paracrystalline scattering. 
     72  So be warned that the calculation is slow. Fitting of any experimental data 
     73  must be resolution smeared for any meaningful fit. This makes a triple integral 
     74  which may be very slow. 
    7375 
    7476This example dataset is produced using 200 data points, 
     
    7779The 2D (Anisotropic model) is based on the reference below where $I(q)$ is 
    7880approximated for 1d scattering. Thus the scattering pattern for 2D may not 
    79 be accurate. 
     81be accurate, particularly at low $q$. For general details of the calculation and angular 
     82dispersions for oriented particles see :ref:`orientation` . 
     83Note that we are not responsible for any incorrectness of the 2D model computation. 
    8084 
    8185.. figure:: img/parallelepiped_angle_definition.png 
     
    154158# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    155159# add 2d test later 
     160# TODO: fix the 2d tests 
    156161q = 4.*pi/220. 
    157162tests = [ 
    158163    [{}, [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]], 
    159     [{'theta': 20.0, 'phi': 30, 'psi': 40.0}, (-0.017, 0.035), 2082.20264399], 
    160     [{'theta': 20.0, 'phi': 30, 'psi': 40.0}, (-0.081, 0.011), 0.436323144781], 
     164    #[{'theta': 20.0, 'phi': 30, 'psi': 40.0}, (-0.017, 0.035), 2082.20264399], 
     165    #[{'theta': 20.0, 'phi': 30, 'psi': 40.0}, (-0.081, 0.011), 0.436323144781], 
    161166    ] 
  • sasmodels/models/capped_cylinder.c

    r592343f rbecded3  
    1 double form_volume(double radius, double radius_cap, double length); 
    2 double Iq(double q, double sld, double solvent_sld, 
    3     double radius, double radius_cap, double length); 
    4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    5     double radius, double radius_cap, double length, double theta, double phi); 
    6  
    71#define INVALID(v) (v.radius_cap < v.radius) 
    82 
     
    148//   radius_cap is the radius of the lens 
    159//   length is the cylinder length, or the separation between the lens halves 
    16 //   alpha is the angle of the cylinder wrt q. 
     10//   theta is the angle of the cylinder wrt q. 
    1711static double 
    18 _cap_kernel(double q, double h, double radius_cap, 
    19                       double half_length, double sin_alpha, double cos_alpha) 
     12_cap_kernel(double qab, double qc, double h, double radius_cap, 
     13    double half_length) 
    2014{ 
    2115    // translate a point in [-1,1] to a point in [lower,upper] 
     
    2620 
    2721    // cos term in integral is: 
    28     //    cos (q (R t - h + L/2) cos(alpha)) 
     22    //    cos (q (R t - h + L/2) cos(theta)) 
    2923    // so turn it into: 
    3024    //    cos (m t + b) 
    3125    // where: 
    32     //    m = q R cos(alpha) 
    33     //    b = q(L/2-h) cos(alpha) 
    34     const double m = q*radius_cap*cos_alpha; // cos argument slope 
    35     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    36     const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta) 
     26    //    m = q R cos(theta) 
     27    //    b = q(L/2-h) cos(theta) 
     28    const double m = radius_cap*qc; // cos argument slope 
     29    const double b = (half_length-h)*qc; // cos argument intercept 
     30    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3731    double total = 0.0; 
    3832    for (int i=0; i<76 ;i++) { 
    3933        const double t = Gauss76Z[i]*zm + zb; 
    4034        const double radical = 1.0 - t*t; 
    41         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     35        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    4236        const double Fq = cos(m*t + b) * radical * bj; 
    4337        total += Gauss76Wt[i] * Fq; 
     
    5044 
    5145static double 
    52 _fq(double q, double h, double radius_cap, double radius, double half_length, 
    53     double sin_alpha, double cos_alpha) 
     46_fq(double qab, double qc, double h, double radius_cap, double radius, double half_length) 
    5447{ 
    55     const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 
    56     const double bj = sas_2J1x_x(q*radius*sin_alpha); 
    57     const double si = sas_sinx_x(q*half_length*cos_alpha); 
     48    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length); 
     49    const double bj = sas_2J1x_x(radius*qab); 
     50    const double si = sas_sinx_x(half_length*qc); 
    5851    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5952    const double Aq = cap_Fq + cyl_Fq; 
     
    6154} 
    6255 
    63 double form_volume(double radius, double radius_cap, double length) 
     56static double 
     57form_volume(double radius, double radius_cap, double length) 
    6458{ 
    6559    // cap radius should never be less than radius when this is called 
     
    9084} 
    9185 
    92 double Iq(double q, double sld, double solvent_sld, 
    93           double radius, double radius_cap, double length) 
     86static double 
     87Iq(double q, double sld, double solvent_sld, 
     88    double radius, double radius_cap, double length) 
    9489{ 
    9590    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
     
    10196    double total = 0.0; 
    10297    for (int i=0; i<76 ;i++) { 
    103         const double alpha= Gauss76Z[i]*zm + zb; 
    104         double sin_alpha, cos_alpha; // slots to hold sincos function output 
    105         SINCOS(alpha, sin_alpha, cos_alpha); 
    106  
    107         const double Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha); 
    108         // sin_alpha for spherical coord integration 
    109         total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
     98        const double theta = Gauss76Z[i]*zm + zb; 
     99        double sin_theta, cos_theta; // slots to hold sincos function output 
     100        SINCOS(theta, sin_theta, cos_theta); 
     101        const double qab = q*sin_theta; 
     102        const double qc = q*cos_theta; 
     103        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
     104        // scale by sin_theta for spherical coord integration 
     105        total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
    110106    } 
    111107    // translate dx in [-1,1] to dx in [lower,upper] 
     
    118114 
    119115 
    120 double Iqxy(double qx, double qy, 
     116static double 
     117Iqxy(double qab, double qc, 
    121118    double sld, double solvent_sld, double radius, 
    122     double radius_cap, double length, 
    123     double theta, double phi) 
     119    double radius_cap, double length) 
    124120{ 
    125     double q, sin_alpha, cos_alpha; 
    126     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    127  
    128121    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    129     const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 
     122    const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length); 
    130123 
    131124    // Multiply by contrast^2 and convert to cm-1 
  • sasmodels/models/core_shell_bicelle.c

    rb260926 rbecded3  
    1 double form_volume(double radius, double thick_rim, double thick_face, double length); 
    2 double Iq(double q, 
    3           double radius, 
    4           double thick_rim, 
    5           double thick_face, 
    6           double length, 
    7           double core_sld, 
    8           double face_sld, 
    9           double rim_sld, 
    10           double solvent_sld); 
    11  
    12  
    13 double Iqxy(double qx, double qy, 
    14           double radius, 
    15           double thick_rim, 
    16           double thick_face, 
    17           double length, 
    18           double core_sld, 
    19           double face_sld, 
    20           double rim_sld, 
    21           double solvent_sld, 
    22           double theta, 
    23           double phi); 
    24  
    25  
    26 double form_volume(double radius, double thick_rim, double thick_face, double length) 
     1static double 
     2form_volume(double radius, double thick_rim, double thick_face, double length) 
    273{ 
    28     return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face); 
     4    return M_PI*square(radius+thick_rim)*(length+2.0*thick_face); 
    295} 
    306 
    317static double 
    32 bicelle_kernel(double q, 
    33               double rad, 
    34               double radthick, 
    35               double facthick, 
    36               double halflength, 
    37               double rhoc, 
    38               double rhoh, 
    39               double rhor, 
    40               double rhosolv, 
    41               double sin_alpha, 
    42               double cos_alpha) 
     8bicelle_kernel(double qab, 
     9    double qc, 
     10    double radius, 
     11    double thick_radius, 
     12    double thick_face, 
     13    double halflength, 
     14    double sld_core, 
     15    double sld_face, 
     16    double sld_rim, 
     17    double sld_solvent) 
    4318{ 
    44     const double dr1 = rhoc-rhoh; 
    45     const double dr2 = rhor-rhosolv; 
    46     const double dr3 = rhoh-rhor; 
    47     const double vol1 = M_PI*square(rad)*2.0*(halflength); 
    48     const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick); 
    49     const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick); 
     19    const double dr1 = sld_core-sld_face; 
     20    const double dr2 = sld_rim-sld_solvent; 
     21    const double dr3 = sld_face-sld_rim; 
     22    const double vol1 = M_PI*square(radius)*2.0*(halflength); 
     23    const double vol2 = M_PI*square(radius+thick_radius)*2.0*(halflength+thick_face); 
     24    const double vol3 = M_PI*square(radius)*2.0*(halflength+thick_face); 
    5025 
    51     const double be1 = sas_2J1x_x(q*(rad)*sin_alpha); 
    52     const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha); 
    53     const double si1 = sas_sinx_x(q*(halflength)*cos_alpha); 
    54     const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha); 
     26    const double be1 = sas_2J1x_x((radius)*qab); 
     27    const double be2 = sas_2J1x_x((radius+thick_radius)*qab); 
     28    const double si1 = sas_sinx_x((halflength)*qc); 
     29    const double si2 = sas_sinx_x((halflength+thick_face)*qc); 
    5530 
    5631    const double t = vol1*dr1*si1*be1 + 
     
    5833                     vol3*dr3*si2*be1; 
    5934 
    60     const double retval = t*t; 
    61  
    62     return retval; 
    63  
     35    return t; 
    6436} 
    6537 
    6638static double 
    67 bicelle_integration(double q, 
    68                    double rad, 
    69                    double radthick, 
    70                    double facthick, 
    71                    double length, 
    72                    double rhoc, 
    73                    double rhoh, 
    74                    double rhor, 
    75                    double rhosolv) 
     39Iq(double q, 
     40    double radius, 
     41    double thick_radius, 
     42    double thick_face, 
     43    double length, 
     44    double sld_core, 
     45    double sld_face, 
     46    double sld_rim, 
     47    double sld_solvent) 
    7648{ 
    7749    // set up the integration end points 
     
    7951    const double halflength = 0.5*length; 
    8052 
    81     double summ = 0.0; 
     53    double total = 0.0; 
    8254    for(int i=0;i<N_POINTS_76;i++) { 
    83         double alpha = (Gauss76Z[i] + 1.0)*uplim; 
    84         double sin_alpha, cos_alpha; // slots to hold sincos function output 
    85         SINCOS(alpha, sin_alpha, cos_alpha); 
    86         double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 
    87                              halflength, rhoc, rhoh, rhor, rhosolv, 
    88                              sin_alpha, cos_alpha); 
    89         summ += yyy*sin_alpha; 
     55        double theta = (Gauss76Z[i] + 1.0)*uplim; 
     56        double sin_theta, cos_theta; // slots to hold sincos function output 
     57        SINCOS(theta, sin_theta, cos_theta); 
     58        double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
     59                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
     60        total += Gauss76Wt[i]*fq*fq*sin_theta; 
    9061    } 
    9162 
    9263    // calculate value of integral to return 
    93     double answer = uplim*summ; 
    94     return answer; 
     64    double answer = total*uplim; 
     65    return 1.0e-4*answer; 
    9566} 
    9667 
    9768static double 
    98 bicelle_kernel_2d(double qx, double qy, 
    99           double radius, 
    100           double thick_rim, 
    101           double thick_face, 
    102           double length, 
    103           double core_sld, 
    104           double face_sld, 
    105           double rim_sld, 
    106           double solvent_sld, 
    107           double theta, 
    108           double phi) 
     69Iqxy(double qab, double qc, 
     70    double radius, 
     71    double thick_rim, 
     72    double thick_face, 
     73    double length, 
     74    double core_sld, 
     75    double face_sld, 
     76    double rim_sld, 
     77    double solvent_sld) 
    10978{ 
    110     double q, sin_alpha, cos_alpha; 
    111     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    112  
    113     double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
     79    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face, 
    11480                           0.5*length, core_sld, face_sld, rim_sld, 
    115                            solvent_sld, sin_alpha, cos_alpha); 
    116     return 1.0e-4*answer; 
     81                           solvent_sld); 
     82    return 1.0e-4*fq*fq; 
    11783} 
    118  
    119 double Iq(double q, 
    120           double radius, 
    121           double thick_rim, 
    122           double thick_face, 
    123           double length, 
    124           double core_sld, 
    125           double face_sld, 
    126           double rim_sld, 
    127           double solvent_sld) 
    128 { 
    129     double intensity = bicelle_integration(q, radius, thick_rim, thick_face, 
    130                        length, core_sld, face_sld, rim_sld, solvent_sld); 
    131     return intensity*1.0e-4; 
    132 } 
    133  
    134  
    135 double Iqxy(double qx, double qy, 
    136           double radius, 
    137           double thick_rim, 
    138           double thick_face, 
    139           double length, 
    140           double core_sld, 
    141           double face_sld, 
    142           double rim_sld, 
    143           double solvent_sld, 
    144           double theta, 
    145           double phi) 
    146 { 
    147     double intensity = bicelle_kernel_2d(qx, qy, 
    148                       radius, 
    149                       thick_rim, 
    150                       thick_face, 
    151                       length, 
    152                       core_sld, 
    153                       face_sld, 
    154                       rim_sld, 
    155                       solvent_sld, 
    156                       theta, 
    157                       phi); 
    158  
    159     return intensity; 
    160 } 
  • sasmodels/models/core_shell_bicelle_elliptical.c

    rdedcf34 r82592da  
    22static double 
    33form_volume(double r_minor, 
    4         double x_core, 
    5         double thick_rim, 
    6         double thick_face, 
    7         double length) 
     4    double x_core, 
     5    double thick_rim, 
     6    double thick_face, 
     7    double length) 
    88{ 
    99    return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face); 
     
    1212static double 
    1313Iq(double q, 
    14         double r_minor, 
    15         double x_core, 
    16         double thick_rim, 
    17         double thick_face, 
    18         double length, 
    19         double rhoc, 
    20         double rhoh, 
    21         double rhor, 
    22         double rhosolv) 
     14    double r_minor, 
     15    double x_core, 
     16    double thick_rim, 
     17    double thick_face, 
     18    double length, 
     19    double sld_core, 
     20    double sld_face, 
     21    double sld_rim, 
     22    double sld_solvent) 
    2323{ 
    24     double si1,si2,be1,be2; 
    2524     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 
    2625     // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 
    27      //    const double uplim = M_PI_4; 
    2826    const double halfheight = 0.5*length; 
    29     //const double va = 0.0; 
    30     //const double vb = 1.0; 
    31     // inner integral limits 
    32     //const double vaj=0.0; 
    33     //const double vbj=M_PI; 
    34  
    3527    const double r_major = r_minor * x_core; 
    3628    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
    3729    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
    38     const double dr1 = (rhoc-rhoh)   *M_PI*r_minor*r_major*(2.0*halfheight);; 
    39     const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
    40     const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
    41     //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
    42     //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
    43     //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     30    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
     31    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
     32    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     33    const double dr1 = vol1*(sld_core-sld_face); 
     34    const double dr2 = vol2*(sld_rim-sld_solvent); 
     35    const double dr3 = vol3*(sld_face-sld_rim); 
    4436 
    4537    //initialize integral 
     
    4739    for(int i=0;i<76;i++) { 
    4840        //setup inner integral over the ellipsoidal cross-section 
    49         // since we generate these lots of times, why not store them somewhere? 
    50         //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    51         const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
    52         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    53         double inner_sum=0; 
    54         double sinarg1 = q*halfheight*cos_alpha; 
    55         double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 
    56         si1 = sas_sinx_x(sinarg1); 
    57         si2 = sas_sinx_x(sinarg2); 
     41        //const double va = 0.0; 
     42        //const double vb = 1.0; 
     43        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     44        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 
     45        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     46        const double qab = q*sin_theta; 
     47        const double qc = q*cos_theta; 
     48        const double si1 = sas_sinx_x(halfheight*qc); 
     49        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
     50        double inner_sum=0.0; 
    5851        for(int j=0;j<76;j++) { 
    5952            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    60             //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    61             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    62             const double rr = sqrt(r2A - r2B*cos(beta)); 
    63             double besarg1 = q*rr*sin_alpha; 
    64             double besarg2 = q*(rr+thick_rim)*sin_alpha; 
    65             be1 = sas_2J1x_x(besarg1); 
    66             be2 = sas_2J1x_x(besarg2); 
    67             inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 
    68                                               dr2*si2*be2 + 
    69                                               dr3*si2*be1); 
     53            // inner integral limits 
     54            //const double vaj=0.0; 
     55            //const double vbj=M_PI; 
     56            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     57            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 
     58            const double rr = sqrt(r2A - r2B*cos(phi)); 
     59            const double be1 = sas_2J1x_x(rr*qab); 
     60            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
     61            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
     62 
     63            inner_sum += Gauss76Wt[j] * fq * fq; 
    7064        } 
    7165        //now calculate outer integral 
     
    7771 
    7872static double 
    79 Iqxy(double qx, double qy, 
    80           double r_minor, 
    81           double x_core, 
    82           double thick_rim, 
    83           double thick_face, 
    84           double length, 
    85           double rhoc, 
    86           double rhoh, 
    87           double rhor, 
    88           double rhosolv, 
    89           double theta, 
    90           double phi, 
    91           double psi) 
     73Iqxy(double qa, double qb, double qc, 
     74    double r_minor, 
     75    double x_core, 
     76    double thick_rim, 
     77    double thick_face, 
     78    double length, 
     79    double sld_core, 
     80    double sld_face, 
     81    double sld_rim, 
     82    double sld_solvent) 
    9283{ 
    93        // THIS NEEDS TESTING 
    94     double q, xhat, yhat, zhat; 
    95     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    96     const double dr1 = rhoc-rhoh; 
    97     const double dr2 = rhor-rhosolv; 
    98     const double dr3 = rhoh-rhor; 
     84    const double dr1 = sld_core-sld_face; 
     85    const double dr2 = sld_rim-sld_solvent; 
     86    const double dr3 = sld_face-sld_rim; 
    9987    const double r_major = r_minor*x_core; 
    10088    const double halfheight = 0.5*length; 
     
    10492 
    10593    // Compute effective radius in rotated coordinates 
    106     const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat)); 
    107     const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat) 
    108                                    + square((r_minor+thick_rim)*yhat)); 
    109     const double be1 = sas_2J1x_x( q*r_hat ); 
    110     const double be2 = sas_2J1x_x( q*rshell_hat ); 
    111     const double si1 = sas_sinx_x( q*halfheight*zhat ); 
    112     const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat ); 
    113     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1); 
    114     return 1.0e-4 * Aq; 
     94    const double qr_hat = sqrt(square(r_major*qb) + square(r_minor*qa)); 
     95    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qb) 
     96                                   + square((r_minor+thick_rim)*qa)); 
     97    const double be1 = sas_2J1x_x( qr_hat ); 
     98    const double be2 = sas_2J1x_x( qrshell_hat ); 
     99    const double si1 = sas_sinx_x( halfheight*qc ); 
     100    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc ); 
     101    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1; 
     102    return 1.0e-4 * fq*fq; 
    115103} 
    116  
  • sasmodels/models/core_shell_cylinder.c

    r592343f rbecded3  
    1 double form_volume(double radius, double thickness, double length); 
    2 double Iq(double q, double core_sld, double shell_sld, double solvent_sld, 
    3     double radius, double thickness, double length); 
    4 double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld, 
    5     double radius, double thickness, double length, double theta, double phi); 
    6  
    71// vd = volume * delta_rho 
    8 // besarg = q * R * sin(alpha) 
    9 // siarg = q * L/2 * cos(alpha) 
    10 double _cyl(double vd, double besarg, double siarg); 
    11 double _cyl(double vd, double besarg, double siarg) 
     2// besarg = q * R * sin(theta) 
     3// siarg = q * L/2 * cos(theta) 
     4static double _cyl(double vd, double besarg, double siarg) 
    125{ 
    136    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 
    147} 
    158 
    16 double form_volume(double radius, double thickness, double length) 
     9static double 
     10form_volume(double radius, double thickness, double length) 
    1711{ 
    18     return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
     12    return M_PI*square(radius+thickness)*(length+2.0*thickness); 
    1913} 
    2014 
    21 double Iq(double q, 
     15static double 
     16Iq(double q, 
    2217    double core_sld, 
    2318    double shell_sld, 
     
    2823{ 
    2924    // precalculate constants 
    30     const double core_qr = q*radius; 
    31     const double core_qh = q*0.5*length; 
     25    const double core_r = radius; 
     26    const double core_h = 0.5*length; 
    3227    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    33     const double shell_qr = q*(radius + thickness); 
    34     const double shell_qh = q*(0.5*length + thickness); 
     28    const double shell_r = (radius + thickness); 
     29    const double shell_h = (0.5*length + thickness); 
    3530    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3631    double total = 0.0; 
    37     // double lower=0, upper=M_PI_2; 
    3832    for (int i=0; i<76 ;i++) { 
    39         // translate a point in [-1,1] to a point in [lower,upper] 
    40         //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    41         double sn, cn; 
    42         const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    43         SINCOS(alpha, sn, cn); 
    44         const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn) 
    45             + _cyl(shell_vd, shell_qr*sn, shell_qh*cn); 
    46         total += Gauss76Wt[i] * fq * fq * sn; 
     33        // translate a point in [-1,1] to a point in [0, pi/2] 
     34        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     35        double sin_theta, cos_theta; 
     36        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
     37        SINCOS(theta, sin_theta,  cos_theta); 
     38        const double qab = q*sin_theta; 
     39        const double qc = q*cos_theta; 
     40        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     41            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
     42        total += Gauss76Wt[i] * fq * fq * sin_theta; 
    4743    } 
    4844    // translate dx in [-1,1] to dx in [lower,upper] 
     
    5248 
    5349 
    54 double Iqxy(double qx, double qy, 
     50double Iqxy(double qab, double qc, 
    5551    double core_sld, 
    5652    double shell_sld, 
     
    5854    double radius, 
    5955    double thickness, 
    60     double length, 
    61     double theta, 
    62     double phi) 
     56    double length) 
    6357{ 
    64     double q, sin_alpha, cos_alpha; 
    65     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    66  
    67     const double core_qr = q*radius; 
    68     const double core_qh = q*0.5*length; 
     58    const double core_r = radius; 
     59    const double core_h = 0.5*length; 
    6960    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    70     const double shell_qr = q*(radius + thickness); 
    71     const double shell_qh = q*(0.5*length + thickness); 
     61    const double shell_r = (radius + thickness); 
     62    const double shell_h = (0.5*length + thickness); 
    7263    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    7364 
    74     const double fq = _cyl(core_vd, core_qr*sin_alpha, core_qh*cos_alpha) 
    75         + _cyl(shell_vd, shell_qr*sin_alpha, shell_qh*cos_alpha); 
     65    const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     66        + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    7667    return 1.0e-4 * fq * fq; 
    7768} 
  • sasmodels/models/core_shell_ellipsoid.c

    r0a3d9b2 rbecded3  
    1 double form_volume(double radius_equat_core, 
    2                    double polar_core, 
    3                    double equat_shell, 
    4                    double polar_shell); 
    5 double Iq(double q, 
    6           double radius_equat_core, 
    7           double x_core, 
    8           double thick_shell, 
    9           double x_polar_shell, 
    10           double core_sld, 
    11           double shell_sld, 
    12           double solvent_sld); 
    131 
     2// Converted from Igor function gfn4, using the same pattern as ellipsoid 
     3// for evaluating the parts of the integral. 
     4//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN 
     5//                       BY (53) & (58-59) IN CHEN AND 
     6//                       KOTLARCHYK REFERENCE 
     7// 
     8//       <OBLATE ELLIPSOID> 
     9static double 
     10_cs_ellipsoid_kernel(double qab, double qc, 
     11    double equat_core, double polar_core, 
     12    double equat_shell, double polar_shell, 
     13    double sld_core_shell, double sld_shell_solvent) 
     14{ 
     15    const double qr_core = sqrt(square(equat_core*qab) + square(polar_core*qc)); 
     16    const double si_core = sas_3j1x_x(qr_core); 
     17    const double volume_core = M_4PI_3*equat_core*equat_core*polar_core; 
     18    const double fq_core = si_core*volume_core*sld_core_shell; 
    1419 
    15 double Iqxy(double qx, double qy, 
    16           double radius_equat_core, 
    17           double x_core, 
    18           double thick_shell, 
    19           double x_polar_shell, 
    20           double core_sld, 
    21           double shell_sld, 
    22           double solvent_sld, 
    23           double theta, 
    24           double phi); 
     20    const double qr_shell = sqrt(square(equat_shell*qab) + square(polar_shell*qc)); 
     21    const double si_shell = sas_3j1x_x(qr_shell); 
     22    const double volume_shell = M_4PI_3*equat_shell*equat_shell*polar_shell; 
     23    const double fq_shell = si_shell*volume_shell*sld_shell_solvent; 
    2524 
     25    return fq_core + fq_shell; 
     26} 
    2627 
    27 double form_volume(double radius_equat_core, 
    28                    double x_core, 
    29                    double thick_shell, 
    30                    double x_polar_shell) 
     28static double 
     29form_volume(double radius_equat_core, 
     30    double x_core, 
     31    double thick_shell, 
     32    double x_polar_shell) 
    3133{ 
    3234    const double equat_shell = radius_equat_core + thick_shell; 
     
    3739 
    3840static double 
    39 core_shell_ellipsoid_xt_kernel(double q, 
    40           double radius_equat_core, 
    41           double x_core, 
    42           double thick_shell, 
    43           double x_polar_shell, 
    44           double core_sld, 
    45           double shell_sld, 
    46           double solvent_sld) 
     41Iq(double q, 
     42    double radius_equat_core, 
     43    double x_core, 
     44    double thick_shell, 
     45    double x_polar_shell, 
     46    double core_sld, 
     47    double shell_sld, 
     48    double solvent_sld) 
    4749{ 
    48     const double lolim = 0.0; 
    49     const double uplim = 1.0; 
    50  
    51  
    52     const double delpc = core_sld - shell_sld; //core - shell 
    53     const double delps = shell_sld - solvent_sld; //shell - solvent 
    54  
     50    const double sld_core_shell = core_sld - shell_sld; 
     51    const double sld_shell_solvent = shell_sld - solvent_sld; 
    5552 
    5653    const double polar_core = radius_equat_core*x_core; 
     
    5855    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    5956 
    60     double summ = 0.0;   //initialize intergral 
     57    // translate from [-1, 1] => [0, 1] 
     58    const double m = 0.5; 
     59    const double b = 0.5; 
     60    double total = 0.0;     //initialize intergral 
    6161    for(int i=0;i<76;i++) { 
    62         double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim ); 
    63         double yyy = gfn4(zi, radius_equat_core, polar_core, equat_shell, 
    64                           polar_shell, delpc, delps, q); 
    65         summ += Gauss76Wt[i] * yyy; 
     62        const double cos_theta = Gauss76Z[i]*m + b; 
     63        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     64        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, 
     65            radius_equat_core, polar_core, 
     66            equat_shell, polar_shell, 
     67            sld_core_shell, sld_shell_solvent); 
     68        total += Gauss76Wt[i] * fq * fq; 
    6669    } 
    67     summ *= 0.5*(uplim-lolim); 
     70    total *= m; 
    6871 
    6972    // convert to [cm-1] 
    70     return 1.0e-4 * summ; 
     73    return 1.0e-4 * total; 
    7174} 
    7275 
    7376static double 
    74 core_shell_ellipsoid_xt_kernel_2d(double qx, double qy, 
    75           double radius_equat_core, 
    76           double x_core, 
    77           double thick_shell, 
    78           double x_polar_shell, 
    79           double core_sld, 
    80           double shell_sld, 
    81           double solvent_sld, 
    82           double theta, 
    83           double phi) 
     77Iqxy(double qab, double qc, 
     78    double radius_equat_core, 
     79    double x_core, 
     80    double thick_shell, 
     81    double x_polar_shell, 
     82    double core_sld, 
     83    double shell_sld, 
     84    double solvent_sld) 
    8485{ 
    85     double q, sin_alpha, cos_alpha; 
    86     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    87  
    88     const double sldcs = core_sld - shell_sld; 
    89     const double sldss = shell_sld- solvent_sld; 
     86    const double sld_core_shell = core_sld - shell_sld; 
     87    const double sld_shell_solvent = shell_sld - solvent_sld; 
    9088 
    9189    const double polar_core = radius_equat_core*x_core; 
     
    9391    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    9492 
    95     // Call the IGOR library function to get the kernel: 
    96     // MUST use gfn4 not gf2 because of the def of params. 
    97     double answer = gfn4(cos_alpha, 
    98                   radius_equat_core, 
    99                   polar_core, 
    100                   equat_shell, 
    101                   polar_shell, 
    102                   sldcs, 
    103                   sldss, 
    104                   q); 
     93    double fq = _cs_ellipsoid_kernel(qab, qc, 
     94                  radius_equat_core, polar_core, 
     95                  equat_shell, polar_shell, 
     96                  sld_core_shell, sld_shell_solvent); 
    10597 
    10698    //convert to [cm-1] 
    107     answer *= 1.0e-4; 
    108  
    109     return answer; 
     99    return 1.0e-4 * fq * fq; 
    110100} 
    111  
    112 double Iq(double q, 
    113           double radius_equat_core, 
    114           double x_core, 
    115           double thick_shell, 
    116           double x_polar_shell, 
    117           double core_sld, 
    118           double shell_sld, 
    119           double solvent_sld) 
    120 { 
    121     double intensity = core_shell_ellipsoid_xt_kernel(q, 
    122            radius_equat_core, 
    123            x_core, 
    124            thick_shell, 
    125            x_polar_shell, 
    126            core_sld, 
    127            shell_sld, 
    128            solvent_sld); 
    129  
    130     return intensity; 
    131 } 
    132  
    133  
    134 double Iqxy(double qx, double qy, 
    135           double radius_equat_core, 
    136           double x_core, 
    137           double thick_shell, 
    138           double x_polar_shell, 
    139           double core_sld, 
    140           double shell_sld, 
    141           double solvent_sld, 
    142           double theta, 
    143           double phi) 
    144 { 
    145     double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy, 
    146                        radius_equat_core, 
    147                        x_core, 
    148                        thick_shell, 
    149                        x_polar_shell, 
    150                        core_sld, 
    151                        shell_sld, 
    152                        solvent_sld, 
    153                        theta, 
    154                        phi); 
    155  
    156     return intensity; 
    157 } 
  • sasmodels/models/core_shell_ellipsoid.py

    r30b60d2 r8db25bf  
    141141# pylint: enable=bad-whitespace, line-too-long 
    142142 
    143 source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c", 
    144           "core_shell_ellipsoid.c"] 
     143source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    145144 
    146145def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 
  • sasmodels/models/core_shell_parallelepiped.c

    rc69d6d6 r904cd9c  
    1 double form_volume(double length_a, double length_b, double length_c, 
    2                    double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    3 double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 
    4           double solvent_sld, double length_a, double length_b, double length_c, 
    5           double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    6 double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld, 
    7             double crim_sld, double solvent_sld, double length_a, double length_b, 
    8             double length_c, double thick_rim_a, double thick_rim_b, 
    9             double thick_rim_c, double theta, double phi, double psi); 
    10  
    11 double form_volume(double length_a, double length_b, double length_c, 
    12                    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     1static double 
     2form_volume(double length_a, double length_b, double length_c, 
     3    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    134{ 
    145    //return length_a * length_b * length_c; 
     
    1910} 
    2011 
    21 double Iq(double q, 
     12static double 
     13Iq(double q, 
    2214    double core_sld, 
    2315    double arim_sld, 
     
    6153    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    6254    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present 
    63     double Vot = Vin + V1 + V2 + V3; 
    6455 
    6556    // Scale factors (note that drC is not used later) 
     
    10091            const double form_crim = scale11*si1*si2; 
    10192 
    102  
    10393            //  correct FF : sum of square of phase factors 
    10494            inner_total += Gauss76Wt[j] * form * form; 
     
    118108} 
    119109 
    120 double Iqxy(double qx, double qy, 
     110static double 
     111Iqxy(double qa, double qb, double qc, 
    121112    double core_sld, 
    122113    double arim_sld, 
     
    129120    double thick_rim_a, 
    130121    double thick_rim_b, 
    131     double thick_rim_c, 
    132     double theta, 
    133     double phi, 
    134     double psi) 
     122    double thick_rim_c) 
    135123{ 
    136     double q, zhat, yhat, xhat; 
    137     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    138  
    139124    // cspkernel in csparallelepiped recoded here 
    140125    const double dr0 = core_sld-solvent_sld; 
     
    159144    double tc = length_c + 2.0*thick_rim_c; 
    160145    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    161     double siA = sas_sinx_x(0.5*q*length_a*xhat); 
    162     double siB = sas_sinx_x(0.5*q*length_b*yhat); 
    163     double siC = sas_sinx_x(0.5*q*length_c*zhat); 
    164     double siAt = sas_sinx_x(0.5*q*ta*xhat); 
    165     double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    166     double siCt = sas_sinx_x(0.5*q*tc*zhat); 
     146    double siA = sas_sinx_x(0.5*length_a*qa); 
     147    double siB = sas_sinx_x(0.5*length_b*qb); 
     148    double siC = sas_sinx_x(0.5*length_c*qc); 
     149    double siAt = sas_sinx_x(0.5*ta*qa); 
     150    double siBt = sas_sinx_x(0.5*tb*qb); 
     151    double siCt = sas_sinx_x(0.5*tc*qc); 
    167152 
    168153 
  • sasmodels/models/core_shell_parallelepiped.py

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

    r592343f rbecded3  
    1 double form_volume(double radius, double length); 
    2 double fq(double q, double sn, double cn,double radius, double length); 
    3 double orient_avg_1D(double q, double radius, double length); 
    4 double Iq(double q, double sld, double solvent_sld, double radius, double length); 
    5 double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    6     double radius, double length, double theta, double phi); 
    7  
    81#define INVALID(v) (v.radius<0 || v.length<0) 
    92 
    10 double form_volume(double radius, double length) 
     3static double 
     4form_volume(double radius, double length) 
    115{ 
    126    return M_PI*radius*radius*length; 
    137} 
    148 
    15 double fq(double q, double sn, double cn, double radius, double length) 
     9static double 
     10fq(double qab, double qc, double radius, double length) 
    1611{ 
    17     // precompute qr and qh to save time in the loop 
    18     const double qr = q*radius; 
    19     const double qh = q*0.5*length;  
    20     return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn); 
     12    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
    2113} 
    2214 
    23 double orient_avg_1D(double q, double radius, double length) 
     15static double 
     16orient_avg_1D(double q, double radius, double length) 
    2417{ 
    2518    // translate a point in [-1,1] to a point in [0, pi/2] 
    2619    const double zm = M_PI_4; 
    27     const double zb = M_PI_4;  
     20    const double zb = M_PI_4; 
    2821 
    2922    double total = 0.0; 
    3023    for (int i=0; i<76 ;i++) { 
    31         const double alpha = Gauss76Z[i]*zm + zb; 
    32         double sn, cn; // slots to hold sincos function output 
    33         // alpha(theta,phi) the projection of the cylinder on the detector plane 
    34         SINCOS(alpha, sn, cn); 
    35         total += Gauss76Wt[i] * square( fq(q, sn, cn, radius, length) ) * sn; 
     24        const double theta = Gauss76Z[i]*zm + zb; 
     25        double sin_theta, cos_theta; // slots to hold sincos function output 
     26        // theta (theta,phi) the projection of the cylinder on the detector plane 
     27        SINCOS(theta , sin_theta, cos_theta); 
     28        const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
     29        total += Gauss76Wt[i] * form * form * sin_theta; 
    3630    } 
    3731    // translate dx in [-1,1] to dx in [lower,upper] 
     
    3933} 
    4034 
    41 double Iq(double q, 
     35static double 
     36Iq(double q, 
    4237    double sld, 
    4338    double solvent_sld, 
     
    4944} 
    5045 
    51  
    52 double Iqxy(double qx, double qy, 
     46static double 
     47Iqxy(double qab, double qc, 
    5348    double sld, 
    5449    double solvent_sld, 
    5550    double radius, 
    56     double length, 
    57     double theta, 
    58     double phi) 
     51    double length) 
    5952{ 
    60     double q, sin_alpha, cos_alpha; 
    61     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    62     //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha); 
    6353    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    64     const double form = fq(q, sin_alpha, cos_alpha, radius, length); 
     54    const double form = fq(qab, qc, radius, length); 
    6555    return 1.0e-4 * square(s * form); 
    6656} 
  • sasmodels/models/cylinder.py

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

    r3b571ae rbecded3  
    1 double form_volume(double radius_polar, double radius_equatorial); 
    2 double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial); 
    3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    4     double radius_polar, double radius_equatorial, double theta, double phi); 
    5  
    6 double form_volume(double radius_polar, double radius_equatorial) 
     1static double 
     2form_volume(double radius_polar, double radius_equatorial) 
    73{ 
    84    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    95} 
    106 
    11 double Iq(double q, 
     7static  double 
     8Iq(double q, 
    129    double sld, 
    1310    double sld_solvent, 
     
    4138} 
    4239 
    43 double Iqxy(double qx, double qy, 
     40static double 
     41Iqxy(double qab, double qc, 
    4442    double sld, 
    4543    double sld_solvent, 
    4644    double radius_polar, 
    47     double radius_equatorial, 
    48     double theta, 
    49     double phi) 
     45    double radius_equatorial) 
    5046{ 
    51     double q, sin_alpha, cos_alpha; 
    52     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    53     const double r = sqrt(square(radius_equatorial*sin_alpha) 
    54                           + square(radius_polar*cos_alpha)); 
    55     const double f = sas_3j1x_x(q*r); 
     47    const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc)); 
     48    const double f = sas_3j1x_x(qr); 
    5649    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    5750 
    5851    return 1.0e-4 * square(f * s); 
    5952} 
    60  
  • sasmodels/models/ellipsoid.py

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

    r61104c8 r82592da  
    1 double form_volume(double radius_minor, double r_ratio, double length); 
    2 double Iq(double q, double radius_minor, double r_ratio, double length, 
    3           double sld, double solvent_sld); 
    4 double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length, 
    5             double sld, double solvent_sld, double theta, double phi, double psi); 
    6  
    7  
    8 double 
     1static double 
    92form_volume(double radius_minor, double r_ratio, double length) 
    103{ 
     
    125} 
    136 
    14 double 
     7static double 
    158Iq(double q, double radius_minor, double r_ratio, double length, 
    169   double sld, double solvent_sld) 
     
    3528        //const double arg = radius_minor*sin_val; 
    3629        double inner_sum=0; 
    37         for(int j=0;j<20;j++) { 
    38             //20 gauss points for the inner integral 
    39             const double theta = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     30        for(int j=0;j<76;j++) { 
     31            //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017 
     32            const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    4033            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    4134            const double be = sas_2J1x_x(q*r); 
    42             inner_sum += Gauss20Wt[j] * be * be; 
     35            inner_sum += Gauss76Wt[j] * be * be; 
    4336        } 
    4437        //now calculate the value of the inner integral 
     
    6154 
    6255 
    63 double 
    64 Iqxy(double qx, double qy, 
     56static double 
     57Iqxy(double qa, double qb, double qc, 
    6558     double radius_minor, double r_ratio, double length, 
    66      double sld, double solvent_sld, 
    67      double theta, double phi, double psi) 
     59     double sld, double solvent_sld) 
    6860{ 
    69     double q, xhat, yhat, zhat; 
    70     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    71  
    7261    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    7362    // Given:    radius_major = r_ratio * radius_minor 
    74     const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat)); 
    75     const double be = sas_2J1x_x(q*r); 
    76     const double si = sas_sinx_x(q*zhat*0.5*length); 
     63    const double qr = radius_minor*sqrt(square(r_ratio*qb) + square(qa)); 
     64    const double be = sas_2J1x_x(qr); 
     65    const double si = sas_sinx_x(qc*0.5*length); 
    7766    const double Aq = be * si; 
    7867    const double delrho = sld - solvent_sld; 
  • sasmodels/models/elliptical_cylinder.py

    rd9ec8f9 reda8b30  
    11# pylint: disable=line-too-long 
    22r""" 
    3 Definition for 2D (orientated system) 
    4 ------------------------------------- 
    5  
    6 The angles $\theta$ and $\phi$ define the orientation of the axis of the 
    7 cylinder. The angle $\Psi$ is defined as the orientation of the major 
    8 axis of the ellipse with respect to the vector $Q$. A gaussian polydispersity 
    9 can be added to any of the orientation angles, and also for the minor 
    10 radius and the ratio of the ellipse radii. 
    113 
    124.. figure:: img/elliptical_cylinder_geometry.png 
     
    4436 
    4537 
    46 Definition for 1D (no preferred orientation) 
    47 -------------------------------------------- 
    48  
    49 The form factor is averaged over all possible orientation before normalized 
     38For 1D scattering, with no preferred orientation, the form factor is averaged over all possible orientations and normalized 
    5039by the particle volume 
    5140 
     
    5443    P(q) = \text{scale}  <F^2> / V 
    5544 
    56 To provide easy access to the orientation of the elliptical cylinder, we 
    57 define the axis of the cylinder using two angles $\theta$, $\phi$ and $\Psi$ 
    58 (see :ref:`cylinder orientation <cylinder-angle-definition>`). The angle 
    59 $\Psi$ is the rotational angle around its own long_c axis. 
     45For 2d data the orientation of the particle is required, described using a different set  
     46of angles as in the diagrams below, for further details of the calculation and angular  
     47dispersions  see :ref:`orientation` . 
    6048 
    61 All angle parameters are valid and given only for 2D calculation; ie, an 
    62 oriented system. 
    6349 
    6450.. figure:: img/elliptical_cylinder_angle_definition.png 
    6551 
    66     Definition of angles for oriented elliptical cylinder, where axis_ratio is drawn >1, 
    67     and angle $\Psi$ is now a rotation around the axis of the cylinder. 
     52    Note that the angles here are not the same as in the equations for the scattering function. 
     53    Rotation $\theta$, initially in the $xz$ plane, is carried out first, then 
     54    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder. 
     55    The neutron or X-ray beam is along the $z$ axis. 
    6856 
    6957.. figure:: img/elliptical_cylinder_angle_projection.png 
     
    7361 
    7462The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data. 
    75 On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 
    76 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the 
    77 cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) 
    78 The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to 
    79 understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation 
    80 distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.) 
     63 
    8164 
    8265NB: The 2nd virial coefficient of the cylinder is calculated based on the 
  • sasmodels/models/fcc_paracrystal.c

    r50beefe rf728001  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    // Equations from Matsuoka 17-18-19, multiplied by |q| 
     5    const double a1 = ( qa + qb)/2.0; 
     6    const double a2 = ( qa + qc)/2.0; 
     7    const double a3 = ( qb + qc)/2.0; 
    68 
    7 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _FCCeval(double Theta, double Phi, double temp1, double temp3); 
     9    // Matsuoka 23-24-25 
     10    //     Z_k numerator: 1 - exp(a)^2 
     11    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a) 
     12    // Rewriting numerator 
     13    //         => -(exp(2a) - 1) 
     14    //         => -expm1(2a) 
     15    // Rewriting denominator 
     16    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     17    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
     18    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     19    const double exp_arg = exp(arg); 
     20    const double Zq = -cube(expm1(2.0*arg)) 
     21        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     22          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     24 
     25    return Zq; 
     26} 
    927 
    1028 
    11 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    12  
    13         const double Da = d_factor*dnn; 
    14         const double temp1 = q*q*Da*Da; 
    15         const double temp3 = q*dnn; 
    16  
    17         double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    18         return(retVal); 
     29// occupied volume fraction calculated from lattice symmetry and sphere radius 
     30static double 
     31fcc_volume_fraction(double radius, double dnn) 
     32{ 
     33    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    1934} 
    2035 
    21 double _FCCeval(double Theta, double Phi, double temp1, double temp3) { 
    22  
    23         double result; 
    24         double sin_theta,cos_theta,sin_phi,cos_phi; 
    25         SINCOS(Theta, sin_theta, cos_theta); 
    26         SINCOS(Phi, sin_phi, cos_phi); 
    27  
    28         const double temp6 =  sin_theta; 
    29         const double temp7 =  sin_theta*sin_phi + cos_theta; 
    30         const double temp8 = -sin_theta*cos_phi + cos_theta; 
    31         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi; 
    32  
    33         const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
    34         result = cube(1.0-(temp10*temp10))*temp6 
    35             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    36               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    38  
    39         return (result); 
    40 } 
    41  
    42 double form_volume(double radius){ 
     36static double 
     37form_volume(double radius) 
     38{ 
    4339    return sphere_volume(radius); 
    4440} 
    4541 
    4642 
    47 double Iq(double q, double dnn, 
     43static double Iq(double q, double dnn, 
    4844  double d_factor, double radius, 
    49   double sld, double solvent_sld){ 
     45  double sld, double solvent_sld) 
     46{ 
     47    // translate a point in [-1,1] to a point in [0, 2 pi] 
     48    const double phi_m = M_PI; 
     49    const double phi_b = M_PI; 
     50    // translate a point in [-1,1] to a point in [0, pi] 
     51    const double theta_m = M_PI_2; 
     52    const double theta_b = M_PI_2; 
    5053 
    51         //Volume fraction calculated from lattice symmetry and sphere radius 
    52         const double s1 = dnn*sqrt(2.0); 
    53         const double latticescale = 4.0*sphere_volume(radius/s1); 
     54    double outer_sum = 0.0; 
     55    for(int i=0; i<150; i++) { 
     56        double inner_sum = 0.0; 
     57        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     58        double sin_theta, cos_theta; 
     59        SINCOS(theta, sin_theta, cos_theta); 
     60        const double qc = q*cos_theta; 
     61        const double qab = q*sin_theta; 
     62        for(int j=0;j<150;j++) { 
     63            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     64            double sin_phi, cos_phi; 
     65            SINCOS(phi, sin_phi, cos_phi); 
     66            const double qa = qab*cos_phi; 
     67            const double qb = qab*sin_phi; 
     68            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     69            inner_sum += Gauss150Wt[j] * form; 
     70        } 
     71        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     72        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     73    } 
     74    outer_sum *= theta_m; 
     75    const double Zq = outer_sum/(4.0*M_PI); 
     76    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    5477 
    55     const double va = 0.0; 
    56     const double vb = 2.0*M_PI; 
    57     const double vaj = 0.0; 
    58     const double vbj = M_PI; 
    59  
    60     double summ = 0.0; 
    61     double answer = 0.0; 
    62         for(int i=0; i<150; i++) { 
    63                 //setup inner integral over the ellipsoidal cross-section 
    64                 double summj=0.0; 
    65                 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi 
    66                 for(int j=0;j<150;j++) { 
    67                         //20 gauss points for the inner integral 
    68                         double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta 
    69                         double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi); 
    70                         summj += yyy; 
    71                 } 
    72                 //now calculate the value of the inner integral 
    73                 double answer = (vbj-vaj)/2.0*summj; 
    74  
    75                 //now calculate outer integral 
    76                 summ = summ+(Gauss150Wt[i] * answer); 
    77         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    78  
    79         answer = (vb-va)/2.0*summ; 
    80         answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 
    81  
    82     return answer; 
     78    return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     79} 
    8380 
    8481 
     82static double Iqxy(double qa, double qb, double qc, 
     83    double dnn, double d_factor, double radius, 
     84    double sld, double solvent_sld) 
     85{ 
     86    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
     87    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     88    const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     89    return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
    8590} 
    86  
    87 double Iqxy(double qx, double qy, 
    88     double dnn, double d_factor, double radius, 
    89     double sld, double solvent_sld, 
    90     double theta, double phi, double psi) 
    91 { 
    92     double q, zhat, yhat, xhat; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    94  
    95     const double a1 = yhat + xhat; 
    96     const double a2 = xhat + zhat; 
    97     const double a3 = yhat + zhat; 
    98     const double qd = 0.5*q*dnn; 
    99     const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
    100     const double tanh_qd = tanh(arg); 
    101     const double cosh_qd = cosh(arg); 
    102     const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 
    103                     * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 
    104                     * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 
    105  
    106     //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); 
    107  
    108     const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
    109     //the occupied volume of the lattice 
    110     const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    111     return lattice_scale * Fq; 
    112 } 
  • sasmodels/models/fcc_paracrystal.py

    r8f04da4 r1f159bd  
    6464    \end{array} 
    6565 
    66 **NB**: The calculation of $Z(q)$ is a double numerical integral that 
    67 must be carried out with a high density of points to properly capture 
    68 the sharp peaks of the paracrystalline scattering. So be warned that the 
    69 calculation is SLOW. Go get some coffee. Fitting of any experimental data 
    70 must be resolution smeared for any meaningful fit. This makes a triple 
    71 integral. Very, very slow. Go get lunch! 
     66.. note:: 
     67 
     68  The calculation of $Z(q)$ is a double numerical integral that 
     69  must be carried out with a high density of points to properly capture 
     70  the sharp peaks of the paracrystalline scattering. 
     71  So be warned that the calculation is slow. Fitting of any experimental data 
     72  must be resolution smeared for any meaningful fit. This makes a triple integral 
     73  which may be very slow. 
    7274 
    7375The 2D (Anisotropic model) is based on the reference below where $I(q)$ is 
    7476approximated for 1d scattering. Thus the scattering pattern for 2D may not 
    75 be accurate. Note that we are not responsible for any incorrectness of the 
     77be accurate particularly at low $q$. For general details of the calculation 
     78and angular dispersions for oriented particles see :ref:`orientation` . 
     79Note that we are not responsible for any incorrectness of the 
    76802D model computation. 
    7781 
     
    135139 
    136140# april 10 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
     141# TODO: fix the 2d tests 
    137142q = 4.*pi/220. 
    138143tests = [ 
    139144    [{}, [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]], 
    140     [{}, (-0.047, -0.007), 238.103096286], 
    141     [{}, (0.053, 0.063), 0.863609587796], 
     145    #[{}, (-0.047, -0.007), 238.103096286], 
     146    #[{}, (0.053, 0.063), 0.863609587796], 
    142147] 
  • sasmodels/models/hollow_cylinder.c

    r592343f rbecded3  
    1 double form_volume(double radius, double thickness, double length); 
    2 double Iq(double q, double radius, double thickness, double length, double sld, 
    3         double solvent_sld); 
    4 double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld, 
    5         double solvent_sld, double theta, double phi); 
    6  
    71//#define INVALID(v) (v.radius_core >= v.radius) 
    82 
     
    148} 
    159 
    16  
    1710static double 
    18 _hollow_cylinder_kernel(double q, 
    19     double radius, double thickness, double length, double sin_val, double cos_val) 
     11_fq(double qab, double qc, 
     12    double radius, double thickness, double length) 
    2013{ 
    21     const double qs = q*sin_val; 
    22     const double lam1 = sas_2J1x_x((radius+thickness)*qs); 
    23     const double lam2 = sas_2J1x_x(radius*qs); 
     14    const double lam1 = sas_2J1x_x((radius+thickness)*qab); 
     15    const double lam2 = sas_2J1x_x(radius*qab); 
    2416    const double gamma_sq = square(radius/(radius+thickness)); 
    25     //Note: lim_{thickness -> 0} psi = sas_J0(radius*qs) 
    26     //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs) 
    27     const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 
    28     const double t2 = sas_sinx_x(0.5*q*length*cos_val); 
     17    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab) 
     18    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab) 
     19    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00 
     20    const double t2 = sas_sinx_x(0.5*length*qc); 
    2921    return psi*t2; 
    3022} 
    3123 
    32 double 
     24static double 
    3325form_volume(double radius, double thickness, double length) 
    3426{ 
     
    3830 
    3931 
    40 double 
     32static double 
    4133Iq(double q, dou