Changes in / [55d88b4:a1c5758] in sasmodels


Ignore:
Files:
9 added
1 deleted
37 edited

Legend:

Unmodified
Added
Removed

  • 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/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/magnetism/magnetism.rst

    r1f058ea r2c108a3  
    3131to the $x'$ axis, the possible spin states after the sample are then 
    3232 
    33 No spin-flips (+ +) and (- -) 
     33Non spin-flip (+ +) and (- -) 
    3434 
    35 Spin-flips    (+ -) and (- +) 
     35Spin-flip    (+ -) and (- +) 
    3636 
    3737.. figure:: 
    … …  
    4141$\phi$ and $\theta_{up}$, respectively, then, depending on the spin state of the 
    4242neutrons, the scattering length densities, including the nuclear scattering 
    43 length density ($\beta{_N}$) are 
     43length density $(\beta{_N})$ are 
    4444 
    4545.. math:: 
    4646    \beta_{\pm\pm} =  \beta_N \mp D_M M_{\perp x'} 
    47     \text{ when there are no spin-flips} 
     47    \text{ for non spin-flip states} 
    4848 
    4949and 
    … …  
    5151.. math:: 
    5252    \beta_{\pm\mp} =  -D_M (M_{\perp y'} \pm iM_{\perp z'}) 
    53     \text{ when there are} 
     53    \text{ for spin-flip states} 
    5454 
    5555where 
    5656 
    5757.. 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 
     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 
    6363 
    6464Here, $M_{0x}$, $M_{0x}$, $M_{0z}$ are the x, y and z components 
    … …  
    6666 
    6767.. 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 
     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 
    7171 
    7272and the magnetization angles $\theta_M$ and $\phi_M$ are defined in 

  • explore/jitter.py

    • Property mode changed from 100644 to 100755
    r85190c2 raa6989b  
     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 
    511import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
    612import matplotlib.pyplot as plt 
    713from matplotlib.widgets import Slider, CheckButtons 
    814from matplotlib import cm 
    9  
    1015import numpy as np 
    1116from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1217 
    13 def draw_beam(ax): 
     18def draw_beam(ax, view=(0, 0)): 
     19    """ 
     20    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 
     21    """ 
    1422    #ax.plot([0,0],[0,0],[1,-1]) 
    1523    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    … …  
    2230    x = r*np.outer(np.cos(u), np.ones_like(v)) 
    2331    y = r*np.outer(np.sin(u), np.ones_like(v)) 
    24     z = np.outer(np.ones_like(u), v) 
     32    z = 1.3*np.outer(np.ones_like(u), v) 
     33 
     34    theta, phi = view 
     35    shape = x.shape 
     36    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
     37    points = Rz(phi)*Ry(theta)*points 
     38    x, y, z = [v.reshape(shape) for v in points] 
    2539 
    2640    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    2741 
    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 
     42def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0)): 
     43    """ 
     44    Represent jitter as a set of shapes at different orientations. 
     45    """ 
     46    # set max diagonal to 0.95 
     47    scale = 0.95/sqrt(sum(v**2 for v in size)) 
     48    size = tuple(scale*v for v in size) 
     49    draw_shape = draw_parallelepiped 
     50    #draw_shape = draw_ellipsoid 
    3451 
    3552    #np.random.seed(10) 
    … …  
    6481        [ 1,  1,  1], 
    6582    ] 
     83    dtheta, dphi, dpsi = jitter 
    6684    if dtheta == 0: 
    6785        cloud = [v for v in cloud if v[0] == 0] 
    … …  
    7088    if dpsi == 0: 
    7189        cloud = [v for v in cloud if v[2] == 0] 
    72     draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8) 
     90    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
     91    scale = 1/sqrt(3) if dist == 'rectangle' else 1 
    7392    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) 
     93        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
     94        draw_shape(ax, size, view, delta, alpha=0.8) 
    7695    for v in 'xyz': 
    7796        a, b, c = size 
    … …  
    8099        getattr(ax, v+'axis').label.set_text(v) 
    81100 
    82 def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1): 
     101def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
     102    """Draw an ellipsoid.""" 
    83103    a,b,c = size 
    84     theta, phi, psi = view 
    85     dtheta, dphi, dpsi = shimmy 
    86  
    87104    u = np.linspace(0, 2 * np.pi, steps) 
    88105    v = np.linspace(0, np.pi, steps) 
    … …  
    90107    y = b*np.outer(np.sin(u), np.sin(v)) 
    91108    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] 
     109    x, y, z = transform_xyz(view, jitter, x, y, z) 
    98110 
    99111    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    100112 
    101 def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 
     113    draw_labels(ax, view, jitter, [ 
     114         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     115         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     116         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     117         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     118         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     119         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     120    ]) 
     121 
     122def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
     123    """Draw a parallelepiped.""" 
    102124    a,b,c = size 
    103     theta, phi, psi = view 
    104     dtheta, dphi, dpsi = shimmy 
    105  
    106125    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    107126    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    … …  
    118137    ]) 
    119138 
    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] 
     139    x, y, z = transform_xyz(view, jitter, x, y, z) 
    125140    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    126141 
     142    draw_labels(ax, view, jitter, [ 
     143         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     144         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     145         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     146         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     147         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     148         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     149    ]) 
     150 
    127151def draw_sphere(ax, radius=10., steps=100): 
     152    """Draw a sphere""" 
    128153    u = np.linspace(0, 2 * np.pi, steps) 
    129154    v = np.linspace(0, np.pi, steps) 
    … …  
    134159    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    135160 
    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': 
     161def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian'): 
     162    """ 
     163    Draw the dispersion mesh showing the theta-phi orientations at which 
     164    the model will be evaluated. 
     165    """ 
     166    theta, phi, psi = view 
     167    dtheta, dphi, dpsi = jitter 
     168 
     169    if dist == 'gaussian': 
     170        t = np.linspace(-3, 3, n) 
    146171        weights = exp(-0.5*t**2) 
    147     elif dist == 'rect': 
     172    elif dist == 'rectangle': 
     173        # Note: uses sasmodels ridiculous definition of rectangle width 
     174        t = np.linspace(-1, 1, n)*sqrt(3) 
    148175        weights = np.ones_like(t) 
    149176    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  
    163     ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
    164  
    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  
     177        raise ValueError("expected dist to be 'gaussian' or 'rectangle'") 
     178 
     179    # mesh in theta, phi formed by rotating z 
     180    z = np.matrix([[0], [0], [radius]]) 
     181    points = np.hstack([Rx(phi_i)*Ry(theta_i)*z 
     182                        for theta_i in dtheta*t 
     183                        for phi_i in dphi*t]) 
     184    # rotate relative to beam 
     185    points = orient_relative_to_beam(view, points) 
     186 
     187    w = np.outer(weights*cos(radians(dtheta*t)), weights) 
     188 
     189    x, y, z = [np.array(v).flatten() for v in points] 
     190    ax.scatter(x, y, z, c=w.flatten(), marker='o', vmin=0., vmax=1.) 
     191 
     192def draw_labels(ax, view, jitter, text): 
     193    """ 
     194    Draw text at a particular location. 
     195    """ 
     196    labels, locations, orientations = zip(*text) 
     197    px, py, pz = zip(*locations) 
     198    dx, dy, dz = zip(*orientations) 
     199 
     200    px, py, pz = transform_xyz(view, jitter, px, py, pz) 
     201    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz) 
     202 
     203    # TODO: zdir for labels is broken, and labels aren't appearing. 
     204    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 
     205        zdir = np.asarray(zdir).flatten() 
     206        ax.text(p[0], p[1], p[2], label, zdir=zdir) 
     207 
     208# Definition of rotation matrices comes from wikipedia: 
     209#    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations 
    170210def Rx(angle): 
     211    """Construct a matrix to rotate points about *x* by *angle* degrees.""" 
    171212    a = radians(angle) 
    172     R = [[1., 0., 0.], 
    173          [0.,  cos(a), sin(a)], 
    174          [0., -sin(a), cos(a)]] 
     213    R = [[1, 0, 0], 
     214         [0, +cos(a), -sin(a)], 
     215         [0, +sin(a), +cos(a)]] 
    175216    return np.matrix(R) 
    176217 
    177218def Ry(angle): 
     219    """Construct a matrix to rotate points about *y* by *angle* degrees.""" 
    178220    a = radians(angle) 
    179     R = [[cos(a), 0., -sin(a)], 
    180          [0., 1., 0.], 
    181          [sin(a), 0.,  cos(a)]] 
     221    R = [[+cos(a), 0, +sin(a)], 
     222         [0, 1, 0], 
     223         [-sin(a), 0, +cos(a)]] 
    182224    return np.matrix(R) 
    183225 
    184226def Rz(angle): 
     227    """Construct a matrix to rotate points about *z* by *angle* degrees.""" 
    185228    a = radians(angle) 
    186     R = [[cos(a), -sin(a), 0.], 
    187          [sin(a),  cos(a), 0.], 
    188          [0., 0., 1.]] 
     229    R = [[+cos(a), -sin(a), 0], 
     230         [+sin(a), +cos(a), 0], 
     231         [0, 0, 1]] 
    189232    return np.matrix(R) 
    190233 
    191 def main(): 
     234def transform_xyz(view, jitter, x, y, z): 
     235    """ 
     236    Send a set of (x,y,z) points through the jitter and view transforms. 
     237    """ 
     238    x, y, z = [np.asarray(v) for v in (x, y, z)] 
     239    shape = x.shape 
     240    points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     241    points = apply_jitter(jitter, points) 
     242    points = orient_relative_to_beam(view, points) 
     243    x, y, z = [np.array(v).reshape(shape) for v in points] 
     244    return x, y, z 
     245 
     246def apply_jitter(jitter, points): 
     247    """ 
     248    Apply the jitter transform to a set of points. 
     249 
     250    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     251    """ 
     252    dtheta, dphi, dpsi = jitter 
     253    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     254    return points 
     255 
     256def orient_relative_to_beam(view, points): 
     257    """ 
     258    Apply the view transform to a set of points. 
     259 
     260    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     261    """ 
     262    theta, phi, psi = view 
     263    points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     264    return points 
     265 
     266# translate between number of dimension of dispersity and the number of 
     267# points along each dimension. 
     268PD_N_TABLE = { 
     269    (0, 0, 0): (0, 0, 0),     # 0 
     270    (1, 0, 0): (100, 0, 0),   # 100 
     271    (0, 1, 0): (0, 100, 0), 
     272    (0, 0, 1): (0, 0, 100), 
     273    (1, 1, 0): (30, 30, 0),   # 900 
     274    (1, 0, 1): (30, 0, 30), 
     275    (0, 1, 1): (0, 30, 30), 
     276    (1, 1, 1): (15, 15, 15),  # 3375 
     277} 
     278 
     279def clipped_range(data, portion=1.0, mode='central'): 
     280    """ 
     281    Determine range from data. 
     282 
     283    If *portion* is 1, use full range, otherwise use the center of the range 
     284    or the top of the range, depending on whether *mode* is 'central' or 'top'. 
     285    """ 
     286    if portion == 1.0: 
     287        return data.min(), data.max() 
     288    elif mode == 'central': 
     289        data = np.sort(data.flatten()) 
     290        offset = int(portion*len(data)/2 + 0.5) 
     291        return data[offset], data[-offset] 
     292    elif mode == 'top': 
     293        data = np.sort(data.flatten()) 
     294        offset = int(portion*len(data) + 0.5) 
     295        return data[offset], data[-1] 
     296 
     297def draw_scattering(calculator, ax, view, jitter, dist='gaussian'): 
     298    """ 
     299    Plot the scattering for the particular view. 
     300 
     301    *calculator* is returned from :func:`build_model`.  *ax* are the 3D axes 
     302    on which the data will be plotted.  *view* and *jitter* are the current 
     303    orientation and orientation dispersity.  *dist* is one of the sasmodels 
     304    weight distributions. 
     305    """ 
     306    ## Sasmodels use sqrt(3)*width for the rectangle range; scale to the 
     307    ## proper width for comparison. Commented out since now using the 
     308    ## sasmodels definition of width for rectangle. 
     309    #scale = 1/sqrt(3) if dist == 'rectangle' else 1 
     310    scale = 1 
     311 
     312    # add the orientation parameters to the model parameters 
     313    theta, phi, psi = view 
     314    theta_pd, phi_pd, psi_pd = [scale*v for v in jitter] 
     315    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd>0, phi_pd>0, psi_pd>0)] 
     316    ## increase pd_n for testing jitter integration rather than simple viz 
     317    #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)] 
     318 
     319    pars = dict( 
     320        theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n, 
     321        phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n, 
     322        psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n, 
     323    ) 
     324    pars.update(calculator.pars) 
     325 
     326    # compute the pattern 
     327    qx, qy = calculator._data.x_bins, calculator._data.y_bins 
     328    Iqxy = calculator(**pars).reshape(len(qx), len(qy)) 
     329 
     330    # scale it and draw it 
     331    Iqxy = np.log(Iqxy) 
     332    if calculator.limits: 
     333        # use limits from orientation (0,0,0) 
     334        vmin, vmax = calculator.limits 
     335    else: 
     336        vmin, vmax = clipped_range(Iqxy, portion=0.95, mode='top') 
     337    #print("range",(vmin,vmax)) 
     338    #qx, qy = np.meshgrid(qx, qy) 
     339    if 0: 
     340        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
     341        level[level<0] = 0 
     342        colors = plt.get_cmap()(level) 
     343        ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
     344    elif 1: 
     345        ax.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
     346                    levels=np.linspace(vmin, vmax, 24)) 
     347    else: 
     348        ax.pcolormesh(qx, qy, Iqxy) 
     349 
     350def build_model(model_name, n=150, qmax=0.5, **pars): 
     351    """ 
     352    Build a calculator for the given shape. 
     353 
     354    *model_name* is any sasmodels model.  *n* and *qmax* define an n x n mesh 
     355    on which to evaluate the model.  The remaining parameters are stored in 
     356    the returned calculator as *calculator.pars*.  They are used by 
     357    :func:`draw_scattering` to set the non-orientation parameters in the 
     358    calculation. 
     359 
     360    Returns a *calculator* function which takes a dictionary or parameters and 
     361    produces Iqxy.  The Iqxy value needs to be reshaped to an n x n matrix 
     362    for plotting.  See the :class:`sasmodels.direct_model.DirectModel` class 
     363    for details. 
     364    """ 
     365    from sasmodels.core import load_model_info, build_model 
     366    from sasmodels.data import empty_data2D 
     367    from sasmodels.direct_model import DirectModel 
     368 
     369    model_info = load_model_info(model_name) 
     370    model = build_model(model_info) #, dtype='double!') 
     371    q = np.linspace(-qmax, qmax, n) 
     372    data = empty_data2D(q, q) 
     373    calculator = DirectModel(data, model) 
     374 
     375    # stuff the values for non-orientation parameters into the calculator 
     376    calculator.pars = pars.copy() 
     377    calculator.pars.setdefault('backgound', 1e-3) 
     378 
     379    # fix the data limits so that we can see if the pattern fades 
     380    # under rotation or angular dispersion 
     381    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars) 
     382    Iqxy = np.log(Iqxy) 
     383    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top') 
     384    calculator.limits = vmin, vmax+1 
     385 
     386    return calculator 
     387 
     388def select_calculator(model_name, n=150): 
     389    """ 
     390    Create a model calculator for the given shape. 
     391 
     392    *model_name* is one of sphere, cylinder, ellipsoid, triaxial_ellipsoid, 
     393    parallelepiped or bcc_paracrystal. *n* is the number of points to use 
     394    in the q range.  *qmax* is chosen based on model parameters for the 
     395    given model to show something intersting. 
     396 
     397    Returns *calculator* and tuple *size* (a,b,c) giving minor and major 
     398    equitorial axes and polar axis respectively.  See :func:`build_model` 
     399    for details on the returned calculator. 
     400    """ 
     401    a, b, c = 10, 40, 100 
     402    if model_name == 'sphere': 
     403        calculator = build_model('sphere', n=n, radius=c) 
     404        a = b = c 
     405    elif model_name == 'bcc_paracrystal': 
     406        calculator = build_model('bcc_paracrystal', n=n, dnn=c, 
     407                                  d_factor=0.06, radius=40) 
     408        a = b = c 
     409    elif model_name == 'cylinder': 
     410        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c) 
     411        a = b 
     412    elif model_name == 'ellipsoid': 
     413        calculator = build_model('ellipsoid', n=n, qmax=1.0, 
     414                                 radius_polar=c, radius_equatorial=b) 
     415        a = b 
     416    elif model_name == 'triaxial_ellipsoid': 
     417        calculator = build_model('triaxial_ellipsoid', n=n, qmax=0.5, 
     418                                 radius_equat_minor=a, 
     419                                 radius_equat_major=b, 
     420                                 radius_polar=c) 
     421    elif model_name == 'parallelepiped': 
     422        calculator = build_model('parallelepiped', n=n, a=a, b=b, c=c) 
     423    else: 
     424        raise ValueError("unknown model %s"%model_name) 
     425 
     426    return calculator, (a, b, c) 
     427 
     428def main(model_name='parallelepiped'): 
     429    """ 
     430    Show an interactive orientation and jitter demo. 
     431 
     432    *model_name* is one of the models available in :func:`select_model`. 
     433    """ 
     434    # set up calculator 
     435    calculator, size = select_calculator(model_name, n=150) 
     436 
     437    ## uncomment to set an independent the colour range for every view 
     438    ## If left commented, the colour range is fixed for all views 
     439    calculator.limits = None 
     440 
     441    ## use gaussian distribution unless testing integration 
     442    #dist = 'rectangle' 
     443    dist = 'gaussian' 
     444 
     445    ## initial view 
     446    #theta, dtheta = 70., 10. 
     447    #phi, dphi = -45., 3. 
     448    #psi, dpsi = -45., 3. 
     449    theta, phi, psi = 0, 0, 0 
     450    dtheta, dphi, dpsi = 0, 0, 0 
     451 
     452    ## create the plot window 
    192453    #plt.hold(True) 
    193454    plt.set_cmap('gist_earth') 
    … …  
    196457    #ax = plt.subplot(gs[0], projection='3d') 
    197458    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' 
     459    ax.axis('square') 
    206460 
    207461    axcolor = 'lightgoldenrodyellow' 
     462 
     463    ## add control widgets to plot 
    208464    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    209465    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    … …  
    212468    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    213469    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
     470 
    214471    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    215472    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    216473    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  
     474    # Note: using ridiculous definition of rectangle distribution, whose width 
     475    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
     476    # the maximum width to 90. 
     477    dlimit = 30 if dist == 'gaussian' else 90/sqrt(3) 
     478    sdtheta = Slider(axdtheta, 'dTheta', 0, dlimit, valinit=dtheta) 
     479    sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
     480    sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
     481 
     482    ## callback to draw the new view 
    221483    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 
     484        view = stheta.val, sphi.val, spsi.val 
     485        jitter = sdtheta.val, sdphi.val, sdpsi.val 
     486        # set small jitter as 0 if multiple pd dims 
     487        dims = sum(v > 0 for v in jitter) 
     488        limit = [0, 0, 2, 5][dims] 
     489        jitter = [0 if v < limit else v for v in jitter] 
    224490        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) 
     491        draw_beam(ax, (0, 0)) 
     492        draw_jitter(ax, view, jitter, dist=dist, size=size) 
     493        #draw_jitter(ax, view, (0,0,0)) 
     494        draw_mesh(ax, view, jitter, dist=dist) 
     495        draw_scattering(calculator, ax, view, jitter, dist=dist) 
    229496        plt.gcf().canvas.draw() 
    230497 
     498    ## bind control widgets to view updater 
    231499    stheta.on_changed(lambda v: update(v,'theta')) 
    232500    sphi.on_changed(lambda v: update(v, 'phi')) 
    … …  
    236504    sdpsi.on_changed(lambda v: update(v, 'dpsi')) 
    237505 
     506    ## initialize view 
    238507    update(None, 'phi') 
    239508 
     509    ## go interactive 
    240510    plt.show() 
    241511 
    242512if __name__ == "__main__": 
    243     main() 
     513    model_name = sys.argv[1] if len(sys.argv) > 1 else 'parallelepiped' 
     514    main(model_name) 

  • explore/precision.py

    r237c9cf r2a602c7  
    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 r31eea1f  
    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" 
     993<<<<<<< HEAD 
     994        plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
     995        plt.xscale('log' if view == 'log' else linear) 
     996        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best') 
     997======= 
    964998        plot_theory(data, None, resid=err, view=view, use_data=use_data) 
    965999        plt.yscale(errview) 
     1000>>>>>>> master 
    9661001        plt.title("max %s = %.3g"%(errstr, abs(err).max())) 
    9671002        #cbar_title = errstr if errview=="linear" else "log "+errstr 
    … …  
    9951030OPTIONS = [ 
    9961031    # Plotting 
    997     'plot', 'noplot', 
     1032    'plot', 'noplot', 'weights', 
    9981033    'linear', 'log', 'q4', 
    9991034    'rel', 'abs', 
    … …  
    10101045    'demo', 'default',  # TODO: remove demo/default 
    10111046    'nopars', 'pars', 
     1047    'sphere', 'sphere=', # integrate over a sphere in 2d with n points 
    10121048 
    10131049    # Calculation options 
    … …  
    10181054 
    10191055    # Precision options 
    1020     'calc=', 
     1056    'engine=', 
    10211057    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!', 
    10221058    'sasview',  # TODO: remove sasview 3.x support 
    … …  
    11451181        'sets'      : 0, 
    11461182        'engine'    : 'default', 
    1147         'evals'     : '1', 
     1183        'count'     : '1', 
     1184        'show_weights' : False, 
     1185        'sphere'    : 0, 
    11481186    } 
    11491187    for arg in flags: 
    … …  
    11681206        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:] 
    11691207        elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:] 
    1170         elif arg.startswith('-random='):   opts['seed'] = int(arg[8:]) 
    11711208        elif arg.startswith('-title='):    opts['title'] = arg[7:] 
    11721209        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) 
     1210        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
     1211        elif arg.startswith('-neval='):    opts['count'] = arg[7:] 
     1212        elif arg.startswith('-random='): 
     1213            opts['seed'] = int(arg[8:]) 
     1214            opts['sets'] = 0 
     1215        elif arg == '-random': 
     1216            opts['seed'] = np.random.randint(1000000) 
     1217            opts['sets'] = 0 
     1218        elif arg.startswith('-sphere'): 
     1219            opts['sphere'] = int(arg[8:]) if len(arg) > 7 else 150 
     1220            opts['is2d'] = True 
    11761221        elif arg == '-preset':  opts['seed'] = -1 
    11771222        elif arg == '-mono':    opts['mono'] = True 
    … …  
    11961241        elif arg == '-demo':    opts['use_demo'] = True 
    11971242        elif arg == '-default': opts['use_demo'] = False 
     1243        elif arg == '-weights': opts['show_weights'] = True 
    11981244        elif arg == '-html':    opts['html'] = True 
    11991245        elif arg == '-help':    opts['html'] = True 
    … …  
    12321278 
    12331279    if PAR_SPLIT in opts['engine']: 
    1234         engine_types = opts['engine'].split(PAR_SPLIT, 2) 
     1280        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
    12351281        comparison = True 
    12361282    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)] 
     1283        opts['engine'] = [opts['engine']]*2 
     1284 
     1285    if PAR_SPLIT in opts['count']: 
     1286        opts['count'] = [int(k) for k in opts['count'].split(PAR_SPLIT, 2)] 
    12411287        comparison = True 
    12421288    else: 
    1243         evals = [int(opts['evals'])]*2 
     1289        opts['count'] = [int(opts['count'])]*2 
    12441290 
    12451291    if PAR_SPLIT in opts['cutoff']: 
    1246         cutoff = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)] 
     1292        opts['cutoff'] = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)] 
    12471293        comparison = True 
    12481294    else: 
    1249         cutoff = [float(opts['cutoff'])]*2 
    1250  
    1251     base = make_engine(model_info[0], data, engine_types[0], cutoff[0]) 
     1295        opts['cutoff'] = [float(opts['cutoff'])]*2 
     1296 
     1297    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
    12521298    if comparison: 
    1253         comp = make_engine(model_info[1], data, engine_types[1], cutoff[1]) 
     1299        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
    12541300    else: 
    12551301        comp = None 
    … …  
    12601306        'data'      : data, 
    12611307        'name'      : names, 
    1262         'def'       : model_info, 
    1263         'count'     : evals, 
     1308        'info'      : model_info, 
    12641309        'engines'   : [base, comp], 
    12651310        'values'    : values, 
    … …  
    12671312    # pylint: enable=bad-whitespace 
    12681313 
     1314    # Set the integration parameters to the half sphere 
     1315    if opts['sphere'] > 0: 
     1316        set_spherical_integration_parameters(opts, opts['sphere']) 
     1317 
    12691318    return opts 
    12701319 
    1271 def parse_pars(opts): 
    1272     model_info, model_info2 = opts['def'] 
     1320def set_spherical_integration_parameters(opts, steps): 
     1321    """ 
     1322    Set integration parameters for spherical integration over the entire 
     1323    surface in theta-phi coordinates. 
     1324    """ 
     1325    # Set the integration parameters to the half sphere 
     1326    opts['values'].extend([ 
     1327        #'theta=90', 
     1328        'theta_pd=%g'%(90/np.sqrt(3)), 
     1329        'theta_pd_n=%d'%steps, 
     1330        'theta_pd_type=rectangle', 
     1331        #'phi=0', 
     1332        'phi_pd=%g'%(180/np.sqrt(3)), 
     1333        'phi_pd_n=%d'%(2*steps), 
     1334        'phi_pd_type=rectangle', 
     1335        #'background=0', 
     1336    ]) 
     1337    if 'psi' in opts['info'][0].parameters: 
     1338        #opts['values'].append('psi=0') 
     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

    r60335cc 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

    r09141ff re3571cb  
    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 r2bccb5a  
    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 
    … …  
    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    value, 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

    r30b60d2 r6773b02  
    167167import string 
    168168from zlib import crc32 
     169from inspect import currentframe, getframeinfo 
    169170 
    170171import numpy as np  # type: ignore 
    … …  
    370371    """ 
    371372    # Note: need 0xffffffff&val to force an unsigned 32-bit number 
     373    try: 
     374        source = source.encode('utf8') 
     375    except AttributeError: # bytes has no encode attribute in python 3 
     376        pass 
    372377    return "%08X"%(0xffffffff&crc32(source)) 
    373378 
    … …  
    685690    # Load templates and user code 
    686691    kernel_header = load_template('kernel_header.c') 
    687     dll_code = load_template('kernel_iq.c') 
    688     ocl_code = load_template('kernel_iq.cl') 
    689     #ocl_code = load_template('kernel_iq_local.cl') 
     692    kernel_code = load_template('kernel_iq.c') 
    690693    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
     694 
     695    # What kind of 2D model do we need? 
     696    xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str) 
     697               else 'qac' if not partable.is_asymmetric 
     698               else 'qabc') 
    691699 
    692700    # Build initial sources 
    … …  
    713721 
    714722    # Define the parameter table 
    715     # TODO: plug in current line number 
    716     source.append('#line 542 "sasmodels/generate.py"') 
     723    lineno = getframeinfo(currentframe()).lineno + 2 
     724    source.append('#line %d "sasmodels/generate.py"'%lineno) 
     725    #source.append('introduce breakage in generate to test lineno reporting') 
    717726    source.append("#define PARAMETER_TABLE \\") 
    718727    source.append("\\\n".join(p.as_definition() 
    … …  
    730739    source.append(call_volume) 
    731740 
    732     refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 
    733     call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 
    734     if _have_Iqxy(user_code) or isinstance(model_info.Iqxy, str): 
    735         # Call 2D model 
    736         refs = ["_q[2*_i]", "_q[2*_i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 
    737         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 
    738     else: 
    739         # Call 1D model with sqrt(qx^2 + qy^2) 
    740         #warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
    741         # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
    742         pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:] 
    743         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
     741    model_refs = _call_pars("_v.", partable.iq_parameters) 
     742    pars = ",".join(["_q"] + model_refs) 
     743    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     744    if xy_mode == 'qabc': 
     745        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
     746        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
     747        clear_iqxy = "#undef CALL_IQ_ABC" 
     748    elif xy_mode == 'qac': 
     749        pars = ",".join(["_qa", "_qc"] + model_refs) 
     750        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 
     751        clear_iqxy = "#undef CALL_IQ_AC" 
     752    else:  # xy_mode == 'qa' 
     753        pars = ",".join(["_qa"] + model_refs) 
     754        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
     755        clear_iqxy = "#undef CALL_IQ_A" 
    744756 
    745757    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    … …  
    752764    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic) 
    753765    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    754     for k, v in enumerate(magpars[:3]): 
    755         source.append("#define MAGNETIC_PAR%d %d"%(k+1, v)) 
    756766 
    757767    # TODO: allow mixed python/opencl kernels? 
    758768 
    759     ocl = kernels(ocl_code, call_iq, call_iqxy, model_info.name) 
    760     dll = kernels(dll_code, call_iq, call_iqxy, model_info.name) 
     769    ocl = kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
     770    dll = kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    761771    result = { 
    762772        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    … …  
    767777 
    768778 
    769 def kernels(kernel, call_iq, call_iqxy, name): 
     779def kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
    770780    # type: ([str,str], str, str, str) -> List[str] 
    771781    code = kernel[0] 
    … …  
    787797        '#line 1 "%s Iqxy"' % path, 
    788798        code, 
    789         "#undef CALL_IQ", 
     799        clear_iqxy, 
    790800        "#undef KERNEL_NAME", 
    791801    ] 
    … …  
    798808        '#line 1 "%s Imagnetic"' % path, 
    799809        code, 
     810        clear_iqxy, 
    800811        "#undef MAGNETIC", 
    801         "#undef CALL_IQ", 
    802812        "#undef KERNEL_NAME", 
    803813    ] 

  • 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 r2c108a3  
    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 
     16//  NUM_PARS : the number of parameters in the parameter table 
     17//  NUM_VALUES : the number of values to skip at the start of the 
     18//      values array before you get to the dispersity values. 
     19//  PARAMETER_TABLE : list of parameter declarations used to create the 
     20//      ParameterTable type. 
     21//  KERNEL_NAME : model_Iq, model_Iqxy or model_Imagnetic.  This code is 
     22//      included three times, once for each kernel type. 
     23//  MAGNETIC : defined when the magnetic kernel is being instantiated 
     24//  NUM_MAGNETIC : the number of magnetic parameters 
     25//  MAGNETIC_PARS : a comma-separated list of indices to the sld 
     26//      parameters in the parameter table. 
     27//  CALL_VOLUME(table) : call the form volume function 
     28//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
     29//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
     30//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 
     31//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes 
     32//  INVALID(table) : test if the current point is feesible to calculate.  This 
     33//      will be defined in the kernel definition file. 
     34 
    1435#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
    1536#define _PAR_BLOCK_ 
    … …  
    1738typedef struct { 
    1839#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 
     40    int32_t pd_par[MAX_PD];     // id of the nth dispersity variable 
     41    int32_t pd_length[MAX_PD];  // length of the nth dispersity weight vector 
    2142    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector 
    2243    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level 
    … …  
    2546    int32_t num_weights;        // total length of the weights vector 
    2647    int32_t num_active;         // number of non-trivial pd loops 
    27     int32_t theta_par;          // id of spherical correction variable 
     48    int32_t theta_par;          // id of first orientation variable 
    2849} ProblemDetails; 
    2950 
    … …  
    3859#endif // _PAR_BLOCK_ 
    3960 
    40  
    41 #if defined(MAGNETIC) && NUM_MAGNETIC>0 
     61#if defined(MAGNETIC) && NUM_MAGNETIC > 0 
     62// ===== Helper functions for magnetism ===== 
    4263 
    4364// Return value restricted between low and high 
    … …  
    5172//     uu * (sld - m_sigma_x); 
    5273//     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]) 
     74//     ud * (m_sigma_y - 1j*m_sigma_z); 
     75//     du * (m_sigma_y + 1j*m_sigma_z); 
     76// weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag 
     77static void set_spin_weights(double in_spin, double out_spin, double spins[4]) 
    5678{ 
    5779  in_spin = clip(in_spin, 0.0, 1.0); 
    5880  out_spin = clip(out_spin, 0.0, 1.0); 
    5981  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 
     82  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real 
     83  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real 
    6284  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 { 
     85  spins[4] = spins[1]; // du imag 
     86  spins[5] = spins[2]; // ud imag 
     87} 
     88 
     89// Compute the magnetic sld 
     90static double mag_sld( 
     91  const int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 
     92  const double qx, const double qy, 
     93  const double px, const double py, 
     94  const double sld, 
     95  const double mx, const double my, const double mz 
     96) 
     97{ 
     98  if (xs < 4) { 
    6899    const double perp = qy*mx - qx*my; 
    69     return sld + perp*p; 
    70 } 
    71  
    72 #endif // MAGNETIC 
     100    switch (xs) { 
     101      case 0: // uu => sld - D M_perpx 
     102          return sld - px*perp; 
     103      case 1: // ud real => -D M_perpy 
     104          return py*perp; 
     105      case 2: // du real => -D M_perpy 
     106          return py*perp; 
     107      case 3: // dd real => sld + D M_perpx 
     108          return sld + px*perp; 
     109    } 
     110  } else { 
     111    if (xs== 4) { 
     112      return -mz;  // ud imag => -D M_perpz 
     113    } else { // index == 5 
     114      return mz;   // du imag => D M_perpz 
     115    } 
     116  } 
     117} 
     118 
     119 
     120#endif 
     121 
     122// ===== Helper functions for orientation and jitter ===== 
     123 
     124// To change the definition of the angles, run explore/angles.py, which 
     125// uses sympy to generate the equations. 
     126 
     127#if !defined(_QAC_SECTION) && defined(CALL_IQ_AC) 
     128#define _QAC_SECTION 
     129 
     130typedef struct { 
     131    double R31, R32; 
     132} QACRotation; 
     133 
     134// Fill in the rotation matrix R from the view angles (theta, phi) and the 
     135// jitter angles (dtheta, dphi).  This matrix can be applied to all of the 
     136// (qx, qy) points in the image to produce R*[qx,qy]' = [qa,qc]' 
     137static void 
     138qac_rotation( 
     139    QACRotation *rotation, 
     140    double theta, double phi, 
     141    double dtheta, double dphi) 
     142{ 
     143    double sin_theta, cos_theta; 
     144    double sin_phi, cos_phi; 
     145 
     146    // reverse view matrix 
     147    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     148    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     149    const double V11 = cos_phi*cos_theta; 
     150    const double V12 = sin_phi*cos_theta; 
     151    const double V21 = -sin_phi; 
     152    const double V22 = cos_phi; 
     153    const double V31 = sin_theta*cos_phi; 
     154    const double V32 = sin_phi*sin_theta; 
     155 
     156    // reverse jitter matrix 
     157    SINCOS(dtheta*M_PI_180, sin_theta, cos_theta); 
     158    SINCOS(dphi*M_PI_180, sin_phi, cos_phi); 
     159    const double J31 = sin_theta; 
     160    const double J32 = -sin_phi*cos_theta; 
     161    const double J33 = cos_phi*cos_theta; 
     162 
     163    // reverse matrix 
     164    rotation->R31 = J31*V11 + J32*V21 + J33*V31; 
     165    rotation->R32 = J31*V12 + J32*V22 + J33*V32; 
     166} 
     167 
     168// Apply the rotation matrix returned from qac_rotation to the point (qx,qy), 
     169// returning R*[qx,qy]' = [qa,qc]' 
     170static double 
     171qac_apply( 
     172    QACRotation rotation, 
     173    double qx, double qy, 
     174    double *qa_out, double *qc_out) 
     175{ 
     176    const double dqc = rotation.R31*qx + rotation.R32*qy; 
     177    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
     178    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
     179 
     180    *qa_out = dqa; 
     181    *qc_out = dqc; 
     182} 
     183#endif // _QAC_SECTION 
     184 
     185#if !defined(_QABC_SECTION) && defined(CALL_IQ_ABC) 
     186#define _QABC_SECTION 
     187 
     188typedef struct { 
     189    double R11, R12; 
     190    double R21, R22; 
     191    double R31, R32; 
     192} QABCRotation; 
     193 
     194// Fill in the rotation matrix R from the view angles (theta, phi, psi) and the 
     195// jitter angles (dtheta, dphi, dpsi).  This matrix can be applied to all of the 
     196// (qx, qy) points in the image to produce R*[qx,qy]' = [qa,qb,qc]' 
     197static void 
     198qabc_rotation( 
     199    QABCRotation *rotation, 
     200    double theta, double phi, double psi, 
     201    double dtheta, double dphi, double dpsi) 
     202{ 
     203    double sin_theta, cos_theta; 
     204    double sin_phi, cos_phi; 
     205    double sin_psi, cos_psi; 
     206 
     207    // reverse view matrix 
     208    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     209    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     210    SINCOS(psi*M_PI_180, sin_psi, cos_psi); 
     211    const double V11 = -sin_phi*sin_psi + cos_phi*cos_psi*cos_theta; 
     212    const double V12 = sin_phi*cos_psi*cos_theta + sin_psi*cos_phi; 
     213    const double V21 = -sin_phi*cos_psi - sin_psi*cos_phi*cos_theta; 
     214    const double V22 = -sin_phi*sin_psi*cos_theta + cos_phi*cos_psi; 
     215    const double V31 = sin_theta*cos_phi; 
     216    const double V32 = sin_phi*sin_theta; 
     217 
     218    // reverse jitter matrix 
     219    SINCOS(dtheta*M_PI_180, sin_theta, cos_theta); 
     220    SINCOS(dphi*M_PI_180, sin_phi, cos_phi); 
     221    SINCOS(dpsi*M_PI_180, sin_psi, cos_psi); 
     222    const double J11 = cos_psi*cos_theta; 
     223    const double J12 = sin_phi*sin_theta*cos_psi + sin_psi*cos_phi; 
     224    const double J13 = sin_phi*sin_psi - sin_theta*cos_phi*cos_psi; 
     225    const double J21 = -sin_psi*cos_theta; 
     226    const double J22 = -sin_phi*sin_psi*sin_theta + cos_phi*cos_psi; 
     227    const double J23 = sin_phi*cos_psi + sin_psi*sin_theta*cos_phi; 
     228    const double J31 = sin_theta; 
     229    const double J32 = -sin_phi*cos_theta; 
     230    const double J33 = cos_phi*cos_theta; 
     231 
     232    // reverse matrix 
     233    rotation->R11 = J11*V11 + J12*V21 + J13*V31; 
     234    rotation->R12 = J11*V12 + J12*V22 + J13*V32; 
     235    rotation->R21 = J21*V11 + J22*V21 + J23*V31; 
     236    rotation->R22 = J21*V12 + J22*V22 + J23*V32; 
     237    rotation->R31 = J31*V11 + J32*V21 + J33*V31; 
     238    rotation->R32 = J31*V12 + J32*V22 + J33*V32; 
     239} 
     240 
     241// Apply the rotation matrix returned from qabc_rotation to the point (qx,qy), 
     242// returning R*[qx,qy]' = [qa,qb,qc]' 
     243static double 
     244qabc_apply( 
     245    QABCRotation rotation, 
     246    double qx, double qy, 
     247    double *qa_out, double *qb_out, double *qc_out) 
     248{ 
     249    *qa_out = rotation.R11*qx + rotation.R12*qy; 
     250    *qb_out = rotation.R21*qx + rotation.R22*qy; 
     251    *qc_out = rotation.R31*qx + rotation.R32*qy; 
     252} 
     253 
     254#endif // _QABC_SECTION 
     255 
     256 
     257// ==================== KERNEL CODE ======================== 
    73258 
    74259kernel 
    75260void KERNEL_NAME( 
    76261    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 
     262    const int32_t pd_start,     // where we are in the dispersity loop 
     263    const int32_t pd_stop,      // where we are stopping in the dispersity loop 
    79264    global const ProblemDetails *details, 
    80265    global const double *values, 
    81266    global const double *q, // nq q values, with padding to boundary 
    82267    global double *result,  // nq+1 return values, again with padding 
    83     const double cutoff     // cutoff in the polydispersity weight product 
     268    const double cutoff     // cutoff in the dispersity weight product 
    84269    ) 
    85270{ 
     271#ifdef USE_OPENCL 
     272  // who we are and what element we are working with 
     273  const int q_index = get_global_id(0); 
     274  if (q_index >= nq) return; 
     275#else 
     276  // Define q_index here so that debugging statements can be written to work 
     277  // for both OpenCL and DLL using: 
     278  //    if (q_index == 0) {printf(...);} 
     279  int q_index = 0; 
     280#endif 
     281 
    86282  // Storage for the current parameter values.  These will be updated as we 
    87   // walk the polydispersity cube. 
     283  // walk the dispersity mesh. 
    88284  ParameterBlock local_values; 
    89285 
    … …  
    91287  // Location of the sld parameters in the parameter vector. 
    92288  // These parameters are updated with the effective sld due to magnetism. 
    93   #if NUM_MAGNETIC > 3 
    94289  const int32_t slds[] = { MAGNETIC_PARS }; 
    95   #endif 
    96  
    97   // TODO: could precompute these outside of the kernel. 
     290 
    98291  // Interpret polarization cross section. 
    99292  //     up_frac_i = values[NUM_PARS+2]; 
    100293  //     up_frac_f = values[NUM_PARS+3]; 
    101294  //     up_angle = values[NUM_PARS+4]; 
    102   double spins[4]; 
     295  // TODO: could precompute more magnetism parameters before calling the kernel. 
     296  double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
    103297  double cos_mspin, sin_mspin; 
    104   set_spins(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
     298  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
    105299  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
    106300#endif // MAGNETIC 
    … …  
    114308  for (int i=0; i < NUM_PARS; i++) { 
    115309    local_values.vector[i] = values[2+i]; 
    116 //printf("p%d = %g\n",i, local_values.vector[i]); 
     310//if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]); 
    117311  } 
    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; 
     312//if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD); 
     313//if (q_index==0) printf("start:%d  stop:%d\n", pd_start, pd_stop); 
     314 
     315  // If pd_start is zero that means that we are starting a new calculation, 
     316  // and must initialize the result to zero.  Otherwise, we are restarting 
     317  // the calculation from somewhere in the middle of the dispersity mesh, 
     318  // and we update the value rather than reset it. Similarly for the 
     319  // normalization factor, which is stored as the final value in the 
     320  // results vector (one past the number of q values). 
     321  // 
     322  // The code differs slightly between opencl and dll since opencl is only 
     323  // seeing one q value (stored in the variable "this_result") while the dll 
     324  // version must loop over all q. 
     325  #ifdef USE_OPENCL 
     326    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     327    double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
     328  #else // !USE_OPENCL 
     329    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     330    if (pd_start == 0) { 
     331      #ifdef USE_OPENMP 
     332      #pragma omp parallel for 
     333      #endif 
     334      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     335    } 
     336//if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
     337#endif // !USE_OPENCL 
     338 
     339 
     340// ====== macros to set up the parts of the loop ======= 
     341/* 
     342Based on the level of the loop, uses C preprocessor magic to construct 
     343level-specific looping variables, including these from loop level 3: 
     344 
     345  int n3 : length of loop for mesh level 3 
     346  int i3 : current position in the loop for level 3, which is calculated 
     347       from a combination of pd_start, pd_stride[3] and pd_length[3]. 
     348  int p3 : is the index into the parameter table for mesh level 3 
     349  double v3[] : pointer into dispersity array to values for loop 3 
     350  double w3[] : pointer into dispersity array to weights for loop 3 
     351  double weight3 : the product of weights from levels 3 and up, computed 
     352       as weight5*weight4*w3[i3].  Note that we need an outermost 
     353       value weight5 set to 1.0 for this to work properly. 
     354 
     355After expansion, the loop struction will look like the following: 
     356 
     357  // --- PD_INIT(4) --- 
     358  const int n4 = pd_length[4]; 
     359  const int p4 = pd_par[4]; 
     360  global const double *v4 = pd_value + pd_offset[4]; 
     361  global const double *w4 = pd_weight + pd_offset[4]; 
     362  int i4 = (pd_start/pd_stride[4])%n4;  // position in level 4 at pd_start 
     363 
     364  // --- PD_INIT(3) --- 
     365  const int n3 = pd_length[3]; 
     366  ... 
     367  int i3 = (pd_start/pd_stride[3])%n3;  // position in level 3 at pd_start 
     368 
     369  PD_INIT(2) 
     370  PD_INIT(1) 
     371  PD_INIT(0) 
     372 
     373  // --- PD_OUTERMOST_WEIGHT(5) --- 
     374  const double weight5 = 1.0; 
     375 
     376  // --- PD_OPEN(4,5) --- 
     377  while (i4 < n4) { 
     378    parameter[p4] = v4[i4];  // set the value for pd parameter 4 at this mesh point 
     379    const double weight4 = w4[i4] * weight5; 
     380 
     381    // from PD_OPEN(3,4) 
     382    while (i3 < n3) { 
     383      parameter[p3] = v3[i3];  // set the value for pd parameter 3 at this mesh point 
     384      const double weight3 = w3[i3] * weight4; 
     385 
     386      PD_OPEN(3,2) 
     387      PD_OPEN(2,1) 
     388      PD_OPEN(0,1) 
     389 
     390      ... main loop body ... 
     391 
     392      ++step;  // increment counter representing position in dispersity mesh 
     393 
     394      PD_CLOSE(0) 
     395      PD_CLOSE(1) 
     396      PD_CLOSE(2) 
     397 
     398      // --- PD_CLOSE(3) --- 
     399      if (step >= pd_stop) break; 
     400      ++i3; 
     401    } 
     402    i3 = 0; // reset loop counter for next round through the loop 
     403 
     404    // --- PD_CLOSE(4) --- 
     405    if (step >= pd_stop) break; 
     406    ++i4; 
    127407  } 
    128 //printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    129  
     408  i4 = 0; // reset loop counter even though no more rounds through the loop 
     409 
     410*/ 
     411 
     412// Define looping variables 
     413#define PD_INIT(_LOOP) \ 
     414  const int n##_LOOP = details->pd_length[_LOOP]; \ 
     415  const int p##_LOOP = details->pd_par[_LOOP]; \ 
     416  global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \ 
     417  global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \ 
     418  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP; 
     419 
     420// Jump into the middle of the dispersity loop 
     421#define PD_OPEN(_LOOP,_OUTER) \ 
     422  while (i##_LOOP < n##_LOOP) { \ 
     423    local_values.vector[p##_LOOP] = v##_LOOP[i##_LOOP]; \ 
     424    const double weight##_LOOP = w##_LOOP[i##_LOOP] * weight##_OUTER; 
     425 
     426// create the variable "weight#=1.0" where # is the outermost level+1 (=MAX_PD). 
     427#define _PD_OUTERMOST_WEIGHT(_n) const double weight##_n = 1.0; 
     428#define PD_OUTERMOST_WEIGHT(_n) _PD_OUTERMOST_WEIGHT(_n) 
     429 
     430// Close out the loop 
     431#define PD_CLOSE(_LOOP) \ 
     432    if (step >= pd_stop) break; \ 
     433    ++i##_LOOP; \ 
     434  } \ 
     435  i##_LOOP = 0; 
     436 
     437// The main loop body is essentially: 
     438// 
     439//    BUILD_ROTATION      // construct the rotation matrix qxy => qabc 
     440//    for each q 
     441//        FETCH_Q         // set qx,qy from the q input vector 
     442//        APPLY_ROTATION  // convert qx,qy to qa,qb,qc 
     443//        CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...) 
     444// 
     445// Depending on the shape type (radial, axial, triaxial), the variables 
     446// and calling parameters will be slightly different.  These macros 
     447// capture the differences in one spot so the rest of the code is easier 
     448// to read. 
     449#if defined(CALL_IQ) 
     450  // unoriented 1D 
     451  double qk; 
     452  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     453  #define BUILD_ROTATION() do {} while(0) 
     454  #define APPLY_ROTATION() do {} while(0) 
     455  #define CALL_KERNEL() CALL_IQ(qk, local_values.table) 
     456 
     457#elif defined(CALL_IQ_A) 
     458  // unoriented 2D 
     459  double qx, qy; 
     460  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     461  #define BUILD_ROTATION() do {} while(0) 
     462  #define APPLY_ROTATION() do {} while(0) 
     463  #define CALL_KERNEL() CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table) 
     464 
     465#elif defined(CALL_IQ_AC) 
     466  // oriented symmetric 2D 
     467  double qx, qy; 
     468  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     469  double qa, qc; 
     470  QACRotation rotation; 
     471  // Grab the "view" angles (theta, phi) from the initial parameter table. 
     472  // These are separate from the "jitter" angles (dtheta, dphi) which are 
     473  // stored with the dispersity values and copied to the local parameter 
     474  // table as we go through the mesh. 
     475  const double theta = values[details->theta_par+2]; 
     476  const double phi = values[details->theta_par+3]; 
     477  #define BUILD_ROTATION() qac_rotation(&rotation, \ 
     478    theta, phi, local_values.table.theta, local_values.table.phi) 
     479  #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc) 
     480  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table) 
     481 
     482#elif defined(CALL_IQ_ABC) 
     483  // oriented asymmetric 2D 
     484  double qx, qy; 
     485  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     486  double qa, qb, qc; 
     487  QABCRotation rotation; 
     488  // Grab the "view" angles (theta, phi, psi) from the initial parameter table. 
     489  // These are separate from the "jitter" angles (dtheta, dphi, dpsi) which are 
     490  // stored with the dispersity values and copied to the local parameter 
     491  // table as we go through the mesh. 
     492  const double theta = values[details->theta_par+2]; 
     493  const double phi = values[details->theta_par+3]; 
     494  const double psi = values[details->theta_par+4]; 
     495  #define BUILD_ROTATION() qabc_rotation(&rotation, \ 
     496    theta, phi, psi, local_values.table.theta, \ 
     497    local_values.table.phi, local_values.table.psi) 
     498  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
     499  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     500 
     501#endif 
     502 
     503// ====== construct the loops ======= 
     504 
     505// Pointers to the start of the dispersity and weight vectors, if needed. 
    130506#if MAX_PD>0 
    131507  global const double *pd_value = values + NUM_VALUES; 
    … …  
    133509#endif 
    134510 
    135   // Jump into the middle of the polydispersity loop 
     511// The variable "step" is the current position in the dispersity loop. 
     512// It will be incremented each time a new point in the mesh is accumulated, 
     513// and used to test whether we have reached pd_stop. 
     514int step = pd_start; 
     515 
     516// define looping variables 
    136517#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]; 
     518  PD_INIT(4) 
    142519#endif 
    143520#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); 
     521  PD_INIT(3) 
    150522#endif 
    151523#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]; 
     524  PD_INIT(2) 
    157525#endif 
    158526#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]; 
     527  PD_INIT(1) 
    164528#endif 
    165529#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  
     530  PD_INIT(0) 
     531#endif 
     532 
     533// open nested loops 
     534PD_OUTERMOST_WEIGHT(MAX_PD) 
     535#if MAX_PD>4 
     536  PD_OPEN(4,5) 
     537#endif 
     538#if MAX_PD>3 
     539  PD_OPEN(3,4) 
     540#endif 
     541#if MAX_PD>2 
     542  PD_OPEN(2,3) 
     543#endif 
     544#if MAX_PD>1 
     545  PD_OPEN(1,2) 
     546#endif 
    175547#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]; 
     548  PD_OPEN(0,1) 
     549#endif 
     550 
     551 
     552//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");} 
     553 
     554  // ====== loop body ======= 
     555  #ifdef INVALID 
     556  if (!INVALID(local_values.table)) 
     557  #endif 
     558  { 
     559    // Accumulate I(q) 
     560    // Note: weight==0 must always be excluded 
     561    if (weight0 > cutoff) { 
     562      pd_norm += weight0 * CALL_VOLUME(local_values.table); 
     563      BUILD_ROTATION(); 
     564 
     565#ifndef USE_OPENCL 
     566      // DLL needs to explicitly loop over the q values. 
     567      #ifdef USE_OPENMP 
     568      #pragma omp parallel for 
     569      #endif 
     570      for (q_index=0; q_index<nq; q_index++) 
     571#endif // !USE_OPENCL 
     572      { 
     573 
     574        FETCH_Q(); 
     575        APPLY_ROTATION(); 
     576 
     577        // ======= COMPUTE SCATTERING ========== 
     578        #if defined(MAGNETIC) && NUM_MAGNETIC > 0 
     579          // Compute the scattering from the magnetic cross sections. 
     580          double scattering = 0.0; 
    258581          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 
    263582          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); 
     583            // TODO: what is the magnetic scattering at q=0 
     584            const double px = (qy*cos_mspin + qx*sin_mspin)/qsq; 
     585            const double py = (qy*sin_mspin - qx*cos_mspin)/qsq; 
     586 
     587            // loop over uu, ud real, du real, dd, ud imag, du imag 
     588            for (int xs=0; xs<6; xs++) { 
     589              const double xs_weight = spins[xs]; 
     590              if (xs_weight > 1.e-8) { 
     591                // Since the cross section weight is significant, set the slds 
     592                // to the effective slds for this cross section, call the 
     593                // kernel, and add according to weight. 
     594                for (int sk=0; sk<NUM_MAGNETIC; sk++) { 
     595                  const int32_t mag_index = NUM_PARS+5 + 3*sk; 
     596                  const int32_t sld_index = slds[sk]; 
     597                  const double mx = values[mag_index]; 
     598                  const double my = values[mag_index+1]; 
     599                  const double mz = values[mag_index+2]; 
     600                  local_values.vector[sld_index] = 
     601                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz); 
    298602                } 
     603                scattering += xs_weight * CALL_KERNEL(); 
    299604              } 
    300605            } 
    301606          } 
    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); 
    306           result[q_index] += weight * scattering; 
    307         } 
     607        #else  // !MAGNETIC 
     608          const double scattering = CALL_KERNEL(); 
     609        #endif // !MAGNETIC 
     610//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight, weight0); 
     611 
     612        #ifdef USE_OPENCL 
     613          this_result += weight0 * scattering; 
     614        #else // !USE_OPENCL 
     615          result[q_index] += weight0 * scattering; 
     616        #endif // !USE_OPENCL 
    308617      } 
    309618    } 
    310     ++step; 
     619  } 
     620 
     621// close nested loops 
     622++step; 
    311623#if MAX_PD>0 
    312     if (step >= pd_stop) break; 
    313     ++i0; 
    314   } 
    315   i0 = 0; 
     624  PD_CLOSE(0) 
    316625#endif 
    317626#if MAX_PD>1 
    318     if (step >= pd_stop) break; 
    319     ++i1; 
    320   } 
    321   i1 = 0; 
     627  PD_CLOSE(1) 
    322628#endif 
    323629#if MAX_PD>2 
    324     if (step >= pd_stop) break; 
    325     ++i2; 
    326   } 
    327   i2 = 0; 
     630  PD_CLOSE(2) 
    328631#endif 
    329632#if MAX_PD>3 
    330     if (step >= pd_stop) break; 
    331     ++i3; 
    332   } 
    333   i3 = 0; 
     633  PD_CLOSE(3) 
    334634#endif 
    335635#if MAX_PD>4 
    336     if (step >= pd_stop) break; 
    337     ++i4; 
    338   } 
    339   i4 = 0; 
    340 #endif 
    341  
     636  PD_CLOSE(4) 
     637#endif 
     638 
     639// Remember the current result and the updated norm. 
     640#ifdef USE_OPENCL 
     641  result[q_index] = this_result; 
     642  if (q_index == 0) result[nq] = pd_norm; 
     643//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
     644#else // !USE_OPENCL 
     645  result[nq] = pd_norm; 
    342646//printf("res: %g/%g\n", result[0], pd_norm); 
    343   // Remember the updated norm. 
    344   result[nq] = pd_norm; 
    345 } 
     647#endif // !USE_OPENCL 
     648 
     649// clear the macros in preparation for the next kernel 
     650#undef PD_INIT 
     651#undef PD_OPEN 
     652#undef PD_CLOSE 
     653#undef FETCH_Q 
     654#undef BUILD_ROTATION 
     655#undef APPLY_ROTATION 
     656#undef CALL_KERNEL 
     657} 

  • 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 r65314f7  
    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 = [] 

  • sasmodels/modelinfo.py

    r65314f7 re3571cb  
    382382      with vector parameter p sent as p[]. 
    383383 
    384     * *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 
     384    * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 
    385385      function, with vector parameter p sent as p[]. 
    386386 
    … …  
    420420        # type: (List[Parameter]) -> None 
    421421        self.kernel_parameters = parameters 
     422        self._check_angles() 
    422423        self._set_vector_lengths() 
    423424 
    … …  
    438439        self.iq_parameters = [p for p in self.kernel_parameters 
    439440                              if p.type not in ('orientation', 'magnetic')] 
    440         self.iqxy_parameters = [p for p in self.kernel_parameters 
    441                                 if p.type != 'magnetic'] 
     441        # note: orientation no longer sent to Iqxy, so its the same as 
     442        #self.iqxy_parameters = [p for p in self.kernel_parameters 
     443        #                        if p.type != 'magnetic'] 
    442444        self.form_volume_parameters = [p for p in self.kernel_parameters 
    443445                                       if p.type == 'volume'] 
    … …  
    461463        self.has_2d = any(p.type in ('orientation', 'magnetic') 
    462464                          for p in self.kernel_parameters) 
     465        self.is_asymmetric = any(p.name == 'psi' for p in self.kernel_parameters) 
    463466        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    464467                                if p.id.startswith('M0:')] 
    … …  
    467470                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    468471        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
     472 
     473    def _check_angles(self): 
     474        theta = phi = psi = -1 
     475        for k, p in enumerate(self.kernel_parameters): 
     476            if p.name == 'theta': 
     477                theta = k 
     478                if p.type != 'orientation': 
     479                    raise TypeError("theta must be an orientation parameter") 
     480            elif p.name == 'phi': 
     481                phi = k 
     482                if p.type != 'orientation': 
     483                    raise TypeError("phi must be an orientation parameter") 
     484            elif p.name == 'psi': 
     485                psi = k 
     486                if p.type != 'orientation': 
     487                    raise TypeError("psi must be an orientation parameter") 
     488        if theta >= 0 and phi >= 0: 
     489            if phi != theta+1: 
     490                raise TypeError("phi must follow theta") 
     491            if psi >= 0 and psi != phi+1: 
     492                raise TypeError("psi must follow phi") 
     493        elif theta >= 0 or phi >= 0 or psi >= 0: 
     494            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
    469495 
    470496    def __getitem__(self, key): 
    … …  
    476502            raise KeyError("unknown parameter %r"%key) 
    477503        return par 
     504 
     505    def __contains__(self, key): 
     506        for par in self.call_parameters: 
     507            if par.name == key: 
     508                return True 
     509        else: 
     510            return False 
    478511 
    479512    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/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 rbecded3  
    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*qa) + square(r_minor*qb)); 
     95    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 
     96                                   + square((r_minor+thick_rim)*qb)); 
     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 r30b60d2  
    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

    r92dfe0c rbecded3  
    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; 
    15     return length_a * length_b * length_c +  
    16            2.0 * thick_rim_a * length_b * length_c +  
     6    return length_a * length_b * length_c + 
     7           2.0 * thick_rim_a * length_b * length_c + 
    178           2.0 * thick_rim_b * length_a * length_c + 
    189           2.0 * thick_rim_c * length_a * length_b; 
    1910} 
    2011 
    21 double Iq(double q, 
     12static double 
     13Iq(double q, 
    2214    double core_sld, 
    2315    double arim_sld, 
    … …  
    3426    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    3527    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    36      
     28 
    3729    const double mu = 0.5 * q * length_b; 
    38      
     30 
    3931    //calculate volume before rescaling (in original code, but not used) 
    40     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);         
    41     //double vol = length_a * length_b * length_c +  
    42     //       2.0 * thick_rim_a * length_b * length_c +  
     32    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     33    //double vol = length_a * length_b * length_c + 
     34    //       2.0 * thick_rim_a * length_b * length_c + 
    4335    //       2.0 * thick_rim_b * length_a * length_c + 
    4436    //       2.0 * thick_rim_c * length_a * length_b; 
    45      
     37 
    4638    // Scale sides by B 
    4739    const double a_scaled = length_a / length_b; 
    … …  
    10193            //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors 
    10294            // This is the function called by csparallelepiped_analytical_2D_scaled, 
    103             // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    104              
     95            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 
     96 
    10597            //  correct FF : sum of square of phase factors 
    10698            inner_total += Gauss76Wt[j] * form * form; 
    … …  
    118110} 
    119111 
    120 double Iqxy(double qx, double qy, 
     112static double 
     113Iqxy(double qa, double qb, double qc, 
    121114    double core_sld, 
    122115    double arim_sld, 
    … …  
    129122    double thick_rim_a, 
    130123    double thick_rim_b, 
    131     double thick_rim_c, 
    132     double theta, 
    133     double phi, 
    134     double psi) 
     124    double thick_rim_c) 
    135125{ 
    136     double q, zhat, yhat, xhat; 
    137     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    138  
    139126    // cspkernel in csparallelepiped recoded here 
    140127    const double dr0 = core_sld-solvent_sld; 
    … …  
    160147    double tc = length_a + 2.0*thick_rim_c; 
    161148    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162     double siA = sas_sinx_x(0.5*q*length_a*xhat); 
    163     double siB = sas_sinx_x(0.5*q*length_b*yhat); 
    164     double siC = sas_sinx_x(0.5*q*length_c*zhat); 
    165     double siAt = sas_sinx_x(0.5*q*ta*xhat); 
    166     double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    167     double siCt = sas_sinx_x(0.5*q*tc*zhat); 
    168      
     149    double siA = sas_sinx_x(0.5*length_a*qa); 
     150    double siB = sas_sinx_x(0.5*length_b*qb); 
     151    double siC = sas_sinx_x(0.5*length_c*qc); 
     152    double siAt = sas_sinx_x(0.5*ta*qa); 
     153    double siBt = sas_sinx_x(0.5*tb*qb); 
     154    double siCt = sas_sinx_x(0.5*tc*qc); 
     155 
    169156 
    170157    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    … …  
    174161               + drB*siA*(siBt-siB)*siC*V2 
    175162               + drC*siA*siB*(siCt*siCt-siC)*V3); 
    176     
     163 
    177164    return 1.0e-4 * f * f; 
    178165} 

  • 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/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/elliptical_cylinder.c

    r61104c8 rbecded3  
    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) 
    … …  
    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*qa) + square(qb)); 
     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/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/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, double radius, double thickness, double length, 
    4234    double sld, double solvent_sld) 
    4335{ 
    4436    const double lower = 0.0; 
    45     const double upper = 1.0;           //limits of numerical integral 
     37    const double upper = 1.0;        //limits of numerical integral 
    4638 
    47     double summ = 0.0;                  //initialize intergral 
     39    double summ = 0.0;            //initialize intergral 
    4840    for (int i=0;i<76;i++) { 
    49         const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
    50         const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    51         const double inter = _hollow_cylinder_kernel(q, radius, thickness, length, 
    52                                                      sin_val, cos_val); 
    53         summ += Gauss76Wt[i] * inter * inter; 
     41        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     42        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     43        const double form = _fq(q*sin_theta, q*cos_theta, 
     44                                radius, thickness, length); 
     45        summ += Gauss76Wt[i] * form * form; 
    5446    } 
    5547 
    … …  
    5951} 
    6052 
    61 double 
    62 Iqxy(double qx, double qy, 
     53static double 
     54Iqxy(double qab, double qc, 
    6355    double radius, double thickness, double length, 
    64     double sld, double solvent_sld, double theta, double phi) 
     56    double sld, double solvent_sld) 
    6557{ 
    66     double q, sin_alpha, cos_alpha; 
    67     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    68     const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 
    69         sin_alpha, cos_alpha); 
     58    const double form = _fq(qab, qc, radius, thickness, length); 
    7059 
    7160    const double vol = form_volume(radius, thickness, length); 
    72     return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol); 
     61    return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
    7362} 
    74  

  • sasmodels/models/parallelepiped.c

    rd605080 r9b7b23f  
    1 double form_volume(double length_a, double length_b, double length_c); 
    2 double Iq(double q, double sld, double solvent_sld, 
    3     double length_a, double length_b, double length_c); 
    4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    5     double length_a, double length_b, double length_c, 
    6     double theta, double phi, double psi); 
    7  
    8 double form_volume(double length_a, double length_b, double length_c) 
     1static double 
     2form_volume(double length_a, double length_b, double length_c) 
    93{ 
    104    return length_a * length_b * length_c; 
    … …  
    126 
    137 
    14 double Iq(double q, 
     8static double 
     9Iq(double q, 
    1510    double sld, 
    1611    double solvent_sld, 
    … …  
    2015{ 
    2116    const double mu = 0.5 * q * length_b; 
    22      
     17 
    2318    // Scale sides by B 
    2419    const double a_scaled = length_a / length_b; 
    2520    const double c_scaled = length_c / length_b; 
    26          
     21 
    2722    // outer integral (with gauss points), integration limits = 0, 1 
    2823    double outer_total = 0; //initialize integral 
    … …  
    5752 
    5853 
    59 double Iqxy(double qx, double qy, 
     54static double 
     55Iqxy(double qa, double qb, double qc, 
    6056    double sld, 
    6157    double solvent_sld, 
    6258    double length_a, 
    6359    double length_b, 
    64     double length_c, 
    65     double theta, 
    66     double phi, 
    67     double psi) 
     60    double length_c) 
    6861{ 
    69     double q, xhat, yhat, zhat; 
    70     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    71  
    72     const double siA = sas_sinx_x(0.5*length_a*q*xhat); 
    73     const double siB = sas_sinx_x(0.5*length_b*q*yhat); 
    74     const double siC = sas_sinx_x(0.5*length_c*q*zhat); 
     62    const double siA = sas_sinx_x(0.5*length_a*qa); 
     63    const double siB = sas_sinx_x(0.5*length_b*qb); 
     64    const double siC = sas_sinx_x(0.5*length_c*qc); 
    7565    const double V = form_volume(length_a, length_b, length_c); 
    7666    const double drho = (sld - solvent_sld); 

  • sasmodels/models/sc_paracrystal.c

    r50beefe rf728001  
    1 double form_volume(double radius); 
     1static double 
     2sc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    // Equations from Matsuoka 9-10-11, multiplied by |q| 
     5    const double a1 = qa; 
     6    const double a2 = qb; 
     7    const double a3 = qc; 
    28 
    3 double Iq(double q, 
    4           double dnn, 
    5           double d_factor, 
    6           double radius, 
    7           double sphere_sld, 
    8           double solvent_sld); 
     9    // Matsuoka 13-14-15 
     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)); 
    924 
    10 double Iqxy(double qx, double qy, 
    11             double dnn, 
    12             double d_factor, 
    13             double radius, 
    14             double sphere_sld, 
    15             double solvent_sld, 
    16             double theta, 
    17             double phi, 
    18             double psi); 
     25    return Zq; 
     26} 
    1927 
    20 double form_volume(double radius) 
     28// occupied volume fraction calculated from lattice symmetry and sphere radius 
     29static double 
     30sc_volume_fraction(double radius, double dnn) 
     31{ 
     32    return sphere_volume(radius/dnn); 
     33} 
     34 
     35static double 
     36form_volume(double radius) 
    2137{ 
    2238    return sphere_volume(radius); 
    2339} 
    2440 
     41 
    2542static double 
    26 sc_eval(double theta, double phi, double temp3, double temp4, double temp5) 
     43Iq(double q, double dnn, 
     44    double d_factor, double radius, 
     45    double sld, double solvent_sld) 
    2746{ 
    28     double cnt, snt; 
    29     SINCOS(theta, cnt, snt); 
     47    // translate a point in [-1,1] to a point in [0, 2 pi] 
     48    const double phi_m = M_PI_4; 
     49    const double phi_b = M_PI_4; 
     50    // translate a point in [-1,1] to a point in [0, pi] 
     51    const double theta_m = M_PI_4; 
     52    const double theta_b = M_PI_4; 
    3053 
    31     double cnp, snp; 
    32     SINCOS(phi, cnp, snp); 
    3354 
    34         double temp6 = snt; 
    35         double temp7 = -1.0*temp3*snt*cnp; 
    36         double temp8 = temp3*snt*snp; 
    37         double temp9 = temp3*cnt; 
    38         double result = temp6/((1.0-temp4*cos((temp7))+temp5)* 
    39                                (1.0-temp4*cos((temp8))+temp5)* 
    40                                (1.0-temp4*cos((temp9))+temp5)); 
    41         return (result); 
     55    double outer_sum = 0.0; 
     56    for(int i=0; i<150; i++) { 
     57        double inner_sum = 0.0; 
     58        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     59        double sin_theta, cos_theta; 
     60        SINCOS(theta, sin_theta, cos_theta); 
     61        const double qc = q*cos_theta; 
     62        const double qab = q*sin_theta; 
     63        for(int j=0;j<150;j++) { 
     64            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     65            double sin_phi, cos_phi; 
     66            SINCOS(phi, sin_phi, cos_phi); 
     67            const double qa = qab*cos_phi; 
     68            const double qb = qab*sin_phi; 
     69            const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
     70            inner_sum += Gauss150Wt[j] * form; 
     71        } 
     72        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     73        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     74    } 
     75    outer_sum *= theta_m; 
     76    const double Zq = outer_sum/M_PI_2; 
     77    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     78 
     79    return sc_volume_fraction(radius, dnn) * Pq * Zq; 
    4280} 
    4381 
     82 
    4483static double 
    45 sc_integrand(double dnn, double d_factor, double qq, double xx, double yy) 
     84Iqxy(double qa, double qb, double qc, 
     85    double dnn, double d_factor, double radius, 
     86    double sld, double solvent_sld) 
    4687{ 
    47     //Function to calculate integrand values for simple cubic structure 
    48  
    49         double da = d_factor*dnn; 
    50         double temp1 = qq*qq*da*da; 
    51         double temp2 = cube(-expm1(-temp1)); 
    52         double temp3 = qq*dnn; 
    53         double temp4 = 2.0*exp(-0.5*temp1); 
    54         double temp5 = exp(-1.0*temp1); 
    55  
    56         double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2; 
    57  
    58         return(integrand); 
     88    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
     89    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     90    const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor); 
     91    return sc_volume_fraction(radius, dnn) * Pq * Zq; 
    5992} 
    60  
    61 double Iq(double q, 
    62           double dnn, 
    63           double d_factor, 
    64           double radius, 
    65           double sphere_sld, 
    66           double solvent_sld) 
    67 { 
    68         const double va = 0.0; 
    69         const double vb = M_PI_2; //orientation average, outer integral 
    70  
    71     double summ=0.0; 
    72     double answer=0.0; 
    73  
    74         for(int i=0;i<150;i++) { 
    75                 //setup inner integral over the ellipsoidal cross-section 
    76                 double summj=0.0; 
    77                 double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; 
    78                 for(int j=0;j<150;j++) { 
    79                         //150 gauss points for the inner integral 
    80                         double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0; 
    81                         double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij); 
    82                         summj += tmp; 
    83                 } 
    84                 //now calculate the value of the inner integral 
    85                 answer = (vb-va)/2.0*summj; 
    86  
    87                 //now calculate outer integral 
    88                 double tmp = Gauss150Wt[i] * answer; 
    89                 summ += tmp; 
    90         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    91  
    92         answer = (vb-va)/2.0*summ; 
    93  
    94         //Volume fraction calculated from lattice symmetry and sphere radius 
    95         // NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3 
    96         const double latticeScale = sphere_volume(radius/dnn); 
    97  
    98         answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale; 
    99  
    100         return answer; 
    101 } 
    102  
    103 double Iqxy(double qx, double qy, 
    104           double dnn, 
    105           double d_factor, 
    106           double radius, 
    107           double sphere_sld, 
    108           double solvent_sld, 
    109           double theta, 
    110           double phi, 
    111           double psi) 
    112 { 
    113     double q, zhat, yhat, xhat; 
    114     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    115  
    116     const double qd = q*dnn; 
    117     const double arg = 0.5*square(qd*d_factor); 
    118     const double tanh_qd = tanh(arg); 
    119     const double cosh_qd = cosh(arg); 
    120     const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd) 
    121                     * tanh_qd/(1. - cos(qd*yhat)/cosh_qd) 
    122                     * tanh_qd/(1. - cos(qd*xhat)/cosh_qd); 
    123  
    124     const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 
    125     //the occupied volume of the lattice 
    126     const double lattice_scale = sphere_volume(radius/dnn); 
    127     return lattice_scale * Fq; 
    128 } 

  • sasmodels/models/sc_paracrystal.py

    r8f04da4 r8f04da4  
    161161    [{'theta': 10.0, 'phi': 20, 'psi': 30.0}, (0.023, 0.045), 0.0177333171285], 
    162162    ] 
    163  
    164  

  • sasmodels/models/stacked_disks.c

    r19f996b rbecded3  
    1 static double stacked_disks_kernel( 
    2     double q, 
     1static double 
     2stacked_disks_kernel( 
     3    double qab, 
     4    double qc, 
    35    double halfheight, 
    46    double thick_layer, 
    … …  
    911    double layer_sld, 
    1012    double solvent_sld, 
    11     double sin_alpha, 
    12     double cos_alpha, 
    1313    double d) 
    1414 
    … …  
    2020    // zi is the dummy variable for the integration (x in Feigin's notation) 
    2121 
    22     const double besarg1 = q*radius*sin_alpha; 
    23     //const double besarg2 = q*radius*sin_alpha; 
     22    const double besarg1 = radius*qab; 
     23    //const double besarg2 = radius*qab; 
    2424 
    25     const double sinarg1 = q*halfheight*cos_alpha; 
    26     const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 
     25    const double sinarg1 = halfheight*qc; 
     26    const double sinarg2 = (halfheight+thick_layer)*qc; 
    2727 
    2828    const double be1 = sas_2J1x_x(besarg1); 
    … …  
    4343 
    4444    // loop for the structure factor S(q) 
    45     double qd_cos_alpha = q*d*cos_alpha; 
     45    double qd_cos_alpha = d*qc; 
    4646    //d*cos_alpha is the projection of d onto q (in other words the component 
    4747    //of d that is parallel to q. 
    … …  
    6161 
    6262 
    63 static double stacked_disks_1d( 
     63static double 
     64stacked_disks_1d( 
    6465    double q, 
    6566    double thick_core, 
    … …  
    8485        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8586        SINCOS(zi, sin_alpha, cos_alpha); 
    86         double yyy = stacked_disks_kernel(q, 
     87        double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha, 
    8788                           halfheight, 
    8889                           thick_layer, 
    … …  
    9394                           layer_sld, 
    9495                           solvent_sld, 
    95                            sin_alpha, 
    96                            cos_alpha, 
    9796                           d); 
    9897        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    … …  
    105104} 
    106105 
    107 static double form_volume( 
     106static double 
     107form_volume( 
    108108    double thick_core, 
    109109    double thick_layer, 
    … …  
    116116} 
    117117 
    118 static double Iq( 
     118static double 
     119Iq( 
    119120    double q, 
    120121    double thick_core, 
    … …  
    140141 
    141142 
    142 static double Iqxy(double qx, double qy, 
     143static double 
     144Iqxy(double qab, double qc, 
    143145    double thick_core, 
    144146    double thick_layer, 
    … …  
    148150    double core_sld, 
    149151    double layer_sld, 
    150     double solvent_sld, 
    151     double theta, 
    152     double phi) 
     152    double solvent_sld) 
    153153{ 
    154154    int n_stacking = (int)(fp_n_stacking + 0.5); 
    155     double q, sin_alpha, cos_alpha; 
    156     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    157  
    158155    double d = 2.0 * thick_layer + thick_core; 
    159156    double halfheight = 0.5*thick_core; 
    160     double answer = stacked_disks_kernel(q, 
     157    double answer = stacked_disks_kernel(qab, qc, 
    161158                     halfheight, 
    162159                     thick_layer, 
    … …  
    167164                     layer_sld, 
    168165                     solvent_sld, 
    169                      sin_alpha, 
    170                      cos_alpha, 
    171166                     d); 
    172167 

  • sasmodels/models/triaxial_ellipsoid.c

    r68dd6a9 rbecded3  
    1 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar); 
    2 double Iq(double q, double sld, double sld_solvent, 
    3     double radius_equat_minor, double radius_equat_major, double radius_polar); 
    4 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    5     double radius_equat_minor, double radius_equat_major, double radius_polar, double theta, double phi, double psi); 
    6  
    71//#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar) 
    82 
    9  
    10 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     3static double 
     4form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
    115{ 
    126    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
    137} 
    148 
    15 double Iq(double q, 
     9static double 
     10Iq(double q, 
    1611    double sld, 
    1712    double sld_solvent, 
    … …  
    4540    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    4641    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
    47     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    48     return 1.0e-4 * s * s * fqsq; 
     42    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     43    const double drho = (sld - sld_solvent); 
     44    return 1.0e-4 * square(vol*drho) * fqsq; 
    4945} 
    5046 
    51 double Iqxy(double qx, double qy, 
     47static double 
     48Iqxy(double qa, double qb, double qc, 
    5249    double sld, 
    5350    double sld_solvent, 
    5451    double radius_equat_minor, 
    5552    double radius_equat_major, 
    56     double radius_polar, 
    57     double theta, 
    58     double phi, 
    59     double psi) 
     53    double radius_polar) 
    6054{ 
    61     double q, xhat, yhat, zhat; 
    62     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     55    const double qr = sqrt(square(radius_equat_minor*qa) 
     56                           + square(radius_equat_major*qb) 
     57                           + square(radius_polar*qc)); 
     58    const double fq = sas_3j1x_x(qr); 
     59    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     60    const double drho = (sld - sld_solvent); 
    6361 
    64     const double r = sqrt(square(radius_equat_minor*xhat) 
    65                           + square(radius_equat_major*yhat) 
    66                           + square(radius_polar*zhat)); 
    67     const double fq = sas_3j1x_x(q*r); 
    68     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    69  
    70     return 1.0e-4 * square(s * fq); 
     62    return 1.0e-4 * square(vol * drho * fq); 
    7163} 
    72  

  • sasmodels/sasview_model.py

    r9f8ade1 r32f87a5  
    759759        if par.name not in self.params: 
    760760            if par.name == self.multiplicity_info.control: 
    761                 return [self.multiplicity], [1.0] 
     761                return self.multiplicity, [self.multiplicity], [1.0] 
    762762            else: 
    763763                # For hidden parameters use the default value. 
    764                 value = self._model_info.parameters.defaults.get(par.name, np.NaN) 
    765                 return [value], [1.0] 
     764                default = self._model_info.parameters.defaults.get(par.name, np.NaN) 
     765                return [default], [1.0] 
    766766        elif par.polydisperse: 
     767            value = self.params[par.name] 
    767768            dis = self.dispersion[par.name] 
    768769            if dis['type'] == 'array': 
    769                 value, weight = dis['values'], dis['weights'] 
     770                dispersity, weight = dis['values'], dis['weights'] 
    770771            else: 
    771                 value, weight = weights.get_weights( 
     772                dispersity, weight = weights.get_weights( 
    772773                    dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
    773                     self.params[par.name], par.limits, par.relative_pd) 
    774             return value, weight / np.sum(weight) 
    775         else: 
    776             return [self.params[par.name]], [1.0] 
     774                    value, par.limits, par.relative_pd) 
     775            return value, dispersity, weight 
     776        else: 
     777            value = self.params[par.name] 
     778            return value, [value if par.relative_pd else 0.0], [1.0] 
    777779 
    778780def test_cylinder(): 

  • sasmodels/weights.py

    rf1a8811 r3c24ccd  
    5555        """ 
    5656        sigma = self.width * center if relative else self.width 
     57        if not relative: 
     58            # For orientation, the jitter is relative to 0 not the angle 
     59            center = 0 
     60            pass 
    5761        if sigma == 0 or self.npts < 2: 
    5862            if lb <= center <= ub: 
    … …  
    225229    obj = cls(n, width, nsigmas) 
    226230    v, w = obj.get_weights(value, limits[0], limits[1], relative) 
    227     return v, w 
    228  
    229  
    230 def plot_weights(model_info, pairs): 
    231     # type: (ModelInfo, List[Tuple[np.ndarray, np.ndarray]]) -> None 
     231    return v, w/np.sum(w) 
     232 
     233 
     234def plot_weights(model_info, mesh): 
     235    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None 
    232236    """ 
    233237    Plot the weights returned by :func:`get_weights`. 
    234238 
    235     *model_info* is 
    236     :param model_info: 
    237     :param pairs: 
    238     :return: 
     239    *model_info* defines model parameters, etc. 
     240 
     241    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*) 
     242    for each parameter, where (*dispersity*, *weights*) pairs are the 
     243    distributions to be plotted. 
    239244    """ 
    240245    import pylab 
    241246 
    242     if any(len(values)>1 for values, weights in pairs): 
     247    if any(len(dispersity)>1 for value, dispersity, weights in mesh): 
    243248        labels = [p.name for p in model_info.parameters.call_parameters] 
    244         pylab.interactive(True) 
     249        #pylab.interactive(True) 
    245250        pylab.figure() 
    246         for (v,w), s in zip(pairs, labels): 
    247             if len(v) > 1: 
    248                 #print("weights for", s, v, w) 
    249                 pylab.plot(v, w, '-o', label=s) 
     251        for (v,x,w), s in zip(mesh, labels): 
     252            if len(x) > 1: 
     253                pylab.plot(x, w, '-o', label=s) 
    250254        pylab.grid(True) 
    251255        pylab.legend() 
Note: See TracChangeset for help on using the changeset viewer.