Changes in / [a1c5758:55d88b4] in sasmodels


Ignore:
Files:
1 added
9 deleted
37 edited

Legend:

Unmodified
Added
Removed

  • doc/developer/calculator.rst

    r2c108a3 r870a2f4  
    77 
    88This document describes the layer between the form factor kernels and the 
    9 model calculator which implements the dispersity and magnetic SLD 
     9model calculator which implements the polydispersity 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. The C code is defined 
    17 in *kernel_iq.c*, with the minor differences between OpenCL and DLL handled 
    18 by #ifdef statements. 
    19  
    20 The kernel call looks as follows:: 
     16for 1-D, 2-D and 2-D magnetic kernels respectively. 
     17 
     18The C code is defined in *kernel_iq.c* and *kernel_iq.cl* for DLL and OpenCL 
     19respectively.  The kernel call looks as follows:: 
    2120 
    2221  kernel void KERNEL_NAME( 
    2322      int nq,                  // Number of q values in the q vector 
    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 
     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 
    2726      double *values,          // Value and weights vector 
    2827      double *q,               // q or (qx,qy) vector 
    2928      double *result,          // returned I(q), with result[nq] = pd_weight 
    30       double cutoff)           // dispersity weight cutoff 
     29      double cutoff)           // polydispersity weight cutoff 
    3130 
    3231The details for OpenCL and the python loop are slightly different, but these 
    … …  
    3534*nq* indicates the number of q values that will be calculated. 
    3635 
    37 The *pd_start* and *pd_stop* parameters set the range of the dispersity 
    38 loop to compute for the current kernel call.   Give a dispersity 
     36The *pd_start* and *pd_stop* parameters set the range of the polydispersity 
     37loop to compute for the current kernel call.   Give a polydispersity 
    3938calculation with 30 weights for length and 30 weights for radius for example, 
    4039there are a total of 900 calls to the form factor required to compute the 
    … …  
    4342the length index to 3 and the radius index to 10 for a position of 3*30+10=100, 
    4443and could then proceed to position 200.  This allows us to interrupt the 
    45 calculation in the middle of a long dispersity loop without having to 
     44calculation in the middle of a long polydispersity loop without having to 
    4645do special tricks with the C code.  More importantly, it stops the OpenCL 
    4746kernel in a reasonable time; because the GPU is used by the operating 
    … …  
    5049 
    5150The *ProblemDetails* structure is a direct map of the 
    52 :class:`details.CallDetails` buffer.  This indicates which parameters have 
    53 dispersity, and where in the values vector the values and weights can be 
    54 found.  For each parameter with dispersity there is a parameter id, the length 
    55 of the dispersity loop for that parameter, the offset of the parameter 
     51:class:`details.CallDetails` buffer.  This indicates which parameters are 
     52polydisperse, and where in the values vector the values and weights can be 
     53found.  For each polydisperse parameter there is a parameter id, the length 
     54of the polydispersity loop for that parameter, the offset of the parameter 
    5655values in the pd value and pd weight vectors and the 'stride' from one index 
    5756to the next, which is used to translate between the position in the 
    58 dispersity loop and the particular parameter indices.  The *num_eval* 
    59 field is the total size of the dispersity loop.  *num_weights* is the 
     57polydispersity loop and the particular parameter indices.  The *num_eval* 
     58field is the total size of the polydispersity loop.  *num_weights* is the 
    6059number of elements in the pd value and pd weight vectors.  *num_active* is 
    61 the number of non-trivial pd loops (parameters with dispersity should be ordered 
    62 by decreasing pd vector length, with a length of 1 meaning no dispersity). 
     60the number of non-trivial pd loops (polydisperse parameters should be ordered 
     61by decreasing pd vector length, with a length of 1 meaning no polydispersity). 
    6362Oriented objects in 2-D need a cos(theta) spherical correction on the angular 
    6463variation in order to preserve the 'surface area' of the weight distribution. 
    … …  
    7372*(Mx, My, Mz)*.  Sample magnetization is translated from *(M, theta, phi)* 
    7473to *(Mx, My, Mz)* before the kernel is called.   After the fixed values comes 
    75 the pd value vector, with the dispersity values for each parameter 
     74the pd value vector, with the polydispersity values for each parameter 
    7675stacked one after the other.  The order isn't important since the location 
    7776for each parameter is stored in the *pd_offset* field of the *ProblemDetails* 
    … …  
    7978values, the pd weight vector is stored, with the same configuration as the 
    8079pd value vector.  Note that the pd vectors can contain values that are not 
    81 in the dispersity loop; this is used by :class:`mixture.MixtureKernel` 
     80in the polydispersity loop; this is used by :class:`mixture.MixtureKernel` 
    8281to make it easier to call the various component kernels. 
    8382 
    … …  
    8887 
    8988The *results* vector contains one slot for each of the *nq* values, plus 
    90 one extra slot at the end for the weight normalization accumulated across 
    91 all points in the dispersity mesh.  This is required when the dispersity 
    92 loop is broken across several kernel calls. 
     89one extra slot at the end for the current polydisperse normalization.  This 
     90is required when the polydispersity loop is broken across several kernel 
     91calls. 
    9392 
    9493*cutoff* is a importance cutoff so that points which contribute negligibly 
    … …  
    9897 
    9998- USE_OPENCL is defined if running in opencl 
    100 - MAX_PD is the maximum depth of the dispersity loop [model specific] 
     99- MAX_PD is the maximum depth of the polydispersity loop [model specific] 
    101100- NUM_PARS is the number of parameter values in the kernel.  This may be 
    102101  more than the number of parameters if some of the parameters are vector 
    103102  values. 
    104103- NUM_VALUES is the number of fixed values, which defines the offset in the 
    105   value list to the dispersity value and weight vectors. 
     104  value list to the polydisperse value and weight vectors. 
    106105- NUM_MAGNETIC is the number of magnetic SLDs 
    107106- MAGNETIC_PARS is a comma separated list of the magnetic SLDs, indicating 
    108107  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. 
    109111- KERNEL_NAME is the name of the function being declared 
    110112- PARAMETER_TABLE is the declaration of the parameters to the kernel: 
    … …  
    150152    Cylinder2D:: 
    151153 
    152         #define CALL_IQ(q, i, var) Iqxy(qa, qc, \ 
     154        #define CALL_IQ(q, i, var) Iqxy(q[2*i], q[2*i+1], \ 
    153155        var.length, \ 
    154156        var.radius, \ 
    155157        var.sld, \ 
    156         var.sld_solvent) 
     158        var.sld_solvent, \ 
     159        var.theta, \ 
     160        var.phi) 
    157161 
    158162- CALL_VOLUME(var) is similar, but for calling the form volume:: 
    … …  
    178182        #define INVALID(var) constrained(var.p1, var.p2, var.p3) 
    179183 
    180 Our design supports a limited number of dispersity loops, wherein 
    181 we need to cycle through the values of the dispersity, calculate 
     184Our design supports a limited number of polydispersity loops, wherein 
     185we need to cycle through the values of the polydispersity, calculate 
    182186the I(q, p) for each combination of parameters, and perform a normalized 
    183187weighted sum across all the weights.  Parameters may be passed to the 
    184 underlying calculation engine as scalars or vectors, but the dispersity 
     188underlying calculation engine as scalars or vectors, but the polydispersity 
    185189calculator treats the parameter set as one long vector. 
    186190 
    187 Let's assume we have 8 parameters in the model, two of which allow dispersity. 
    188 Since this is a 1-D model the orientation parameters won't be used:: 
     191Let's assume we have 8 parameters in the model, with two polydisperse.  Since 
     192this is a 1-D model the orientation parameters won't be used:: 
    189193 
    190194    0: scale        {scl = constant} 
    … …  
    192196    2: radius       {r = vector of 10pts} 
    193197    3: length       {l = vector of 30pts} 
    194     4: sld          {s1 = constant/(radius**2*length)} 
     198    4: sld          {s = constant/(radius**2*length)} 
    195199    5: sld_solvent  {s2 = constant} 
    196200    6: theta        {not used} 
    … …  
    198202 
    199203This generates the following call to the kernel.  Note that parameters 4 and 
    200 5 are treated as having dispersity even though they don't --- this is because 
     2045 are treated as polydisperse even though they are not --- this is because 
    201205it is harmless to do so and simplifies the looping code:: 
    202206 
    … …  
    206210    NUM_MAGNETIC = 2      // two parameters might be magnetic 
    207211    MAGNETIC_PARS = 4, 5  // they are sld and sld_solvent 
     212    MAGNETIC_PAR0 = 4     // sld index 
     213    MAGNETIC_PAR1 = 5     // solvent index 
    208214 
    209215    details { 
    … …  
    212218        pd_offset = {10, 0, 31, 32}   // *length* starts at index 10 in weights 
    213219        pd_stride = {1, 30, 300, 300} // cumulative product of pd length 
    214         num_eval = 300   // 300 values in the dispersity loop 
     220        num_eval = 300   // 300 values in the polydispersity loop 
    215221        num_weights = 42 // 42 values in the pd vector 
    216222        num_active = 2   // only the first two pd are active 
    … …  
    219225 
    220226    values = { scl, bkg,                                  // universal 
    221                r, l, s1, s2, theta, phi,                  // kernel pars 
     227               r, l, s, s2, theta, phi,                   // kernel pars 
    222228               in spin, out spin, spin angle,             // applied magnetism 
    223                mx s1, my s1, mz s1, mx s2, my s2, mz s2,  // magnetic slds 
     229               mx s, my s, mz s, mx s2, my s2, mz s2,     // magnetic slds 
    224230               r0, .., r9, l0, .., l29, s, s2,            // pd values 
    225231               r0, .., r9, l0, .., l29, s, s2}            // pd weights 
    … …  
    229235    result = {r1, ..., r130, pd_norm, x } 
    230236 
    231 The dispersity parameters are stored in as an array of parameter 
    232 indices, one for each parameter, stored in pd_par[n]. Parameters which do 
    233 not support dispersity do not appear in this array. Each dispersity 
     237The polydisperse parameters are stored in as an array of parameter 
     238indices, one for each polydisperse parameter, stored in pd_par[n]. 
     239Non-polydisperse parameters do not appear in this array. Each polydisperse 
    234240parameter has a weight vector whose length is stored in pd_length[n]. 
    235241The weights are stored in a contiguous vector of weights for all 
    … …  
    237243in pd_offset[n].  The values corresponding to the weights are stored 
    238244together in a separate weights[] vector, with offset stored in 
    239 par_offset[pd_par[n]]. Dispersity parameters should be stored in 
     245par_offset[pd_par[n]]. Polydisperse parameters should be stored in 
    240246decreasing order of length for highest efficiency. 
    241247 
    242 We limit the number of dispersity dimensions to MAX_PD (currently 4), 
    243 though some models may have fewer if they have fewer dispersity 
     248We limit the number of polydisperse dimensions to MAX_PD (currently 4), 
     249though some models may have fewer if they have fewer polydisperse 
    244250parameters.  The main reason for the limit is to reduce code size. 
    245 Each additional dispersity parameter requires a separate dispersity 
    246 loop.  If more than 4 levels of dispersity are needed, then we need to 
    247 switch to a monte carlo importance sampling algorithm with better 
    248 performance for high-dimensional integrals. 
     251Each additional polydisperse parameter requires a separate polydispersity 
     252loop.  If more than 4 levels of polydispersity are needed, then kernel_iq.c 
     253and kernel_iq.cl will need to be extended. 
    249254 
    250255Constraints between parameters are not supported.  Instead users will 
    … …  
    257262theta since the polar coordinates normalization is tied to this parameter. 
    258263 
    259 If there is no dispersity we pretend that we have a disperisty mesh over 
    260 a single parameter with a single point in the distribution, giving 
    261 pd_start=0 and pd_stop=1. 
     264If there is no polydispersity we pretend that it is polydisperisty with one 
     265parameter, pd_start=0 and pd_stop=1.  We may or may not short circuit the 
     266calculation in this case, depending on how much time it saves. 
    262267 
    263268The problem details structure could be allocated and sent in as an integer 
    264269array using the read-only flag.  This would allow us to copy it once per fit 
    265270along with the weights vector, since features such as the number of 
    266 disperity points per pd parameter won't change between function evaluations. 
    267 A new parameter vector must be sent for each I(q) evaluation.  This is 
    268 not currently implemented, and would require some resturcturing of 
    269 the :class:`sasview_model.SasviewModel` interface. 
    270  
    271 The results array will be initialized to zero for dispersity loop 
     271polydisperity elements per pd parameter or the coordinated won't change 
     272between function evaluations.  A new parameter vector must be sent for 
     273each I(q) evaluation.  This is not currently implemented, and would require 
     274some resturcturing of the :class:`sasview_model.SasviewModel` interface. 
     275 
     276The results array will be initialized to zero for polydispersity loop 
    272277entry zero, and preserved between calls to [start, stop] so that the 
    273278results accumulate by the time the loop has completed.  Background and 
    … …  
    290295 
    291296This will require accumulated error for each I(q) value to be preserved 
    292 between kernel calls to implement this fully.  The *kernel_iq.c* code, which 
    293 loops over q for each parameter set in the dispersity loop, will also need 
    294 the accumulation vector. 
     297between kernel calls to implement this fully.  The kernel_iq.c code, which 
     298loops over q for each parameter set in the polydispersity loop, will need 
     299also need the accumalation vector. 

  • doc/genapi.py

    r706f466 r2e66ef5  
    5959    #('alignment', 'GPU data alignment [unused]'), 
    6060    ('bumps_model', 'Bumps interface'), 
     61    ('compare_many', 'Batch compare models on different compute engines'), 
    6162    ('compare', 'Compare models on different compute engines'), 
    62     ('compare_many', 'Batch compare models on different compute engines'), 
    63     ('conversion_table', 'Model conversion table'), 
    6463    ('convert', 'Sasview to sasmodel converter'), 
    6564    ('core', 'Model access'), 
    … …  
    8382    ('sasview_model', 'Sasview interface'), 
    8483    ('sesans', 'SESANS calculation routines'), 
    85     ('special', 'Special functions library'), 
    8684    ('weights', 'Distribution functions'), 
    8785] 

  • doc/guide/magnetism/magnetism.rst

    r2c108a3 r1f058ea  
    3131to the $x'$ axis, the possible spin states after the sample are then 
    3232 
    33 Non spin-flip (+ +) and (- -) 
     33No spin-flips (+ +) and (- -) 
    3434 
    35 Spin-flip    (+ -) and (- +) 
     35Spin-flips    (+ -) 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{ for non spin-flip states} 
     47    \text{ when there are no spin-flips} 
    4848 
    4949and 
    … …  
    5151.. math:: 
    5252    \beta_{\pm\mp} =  -D_M (M_{\perp y'} \pm iM_{\perp z'}) 
    53     \text{ for spin-flip states} 
     53    \text{ when there are} 
    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 100755 to 100644
    raa6989b r85190c2  
    1 #!/usr/bin/env python 
    21""" 
    32Application to explore the difference between sasview 3.x orientation 
    43dispersity and possible replacement algorithms. 
    54""" 
    6 from __future__ import division, print_function 
    7  
    8 import sys, os 
    9 sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
    10  
    115import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
    126import matplotlib.pyplot as plt 
    137from matplotlib.widgets import Slider, CheckButtons 
    148from matplotlib import cm 
     9 
    1510import numpy as np 
    1611from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1712 
    18 def 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     """ 
     13def draw_beam(ax): 
    2214    #ax.plot([0,0],[0,0],[1,-1]) 
    2315    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    … …  
    3022    x = r*np.outer(np.cos(u), np.ones_like(v)) 
    3123    y = r*np.outer(np.sin(u), np.ones_like(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] 
     24    z = np.outer(np.ones_like(u), v) 
    3925 
    4026    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    4127 
    42 def 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 
     28def 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 
    5134 
    5235    #np.random.seed(10) 
    … …  
    8164        [ 1,  1,  1], 
    8265    ] 
    83     dtheta, dphi, dpsi = jitter 
    8466    if dtheta == 0: 
    8567        cloud = [v for v in cloud if v[0] == 0] 
    … …  
    8870    if dpsi == 0: 
    8971        cloud = [v for v in cloud if v[2] == 0] 
    90     draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    91     scale = 1/sqrt(3) if dist == 'rectangle' else 1 
     72    draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8) 
    9273    for point in cloud: 
    93         delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
    94         draw_shape(ax, size, view, delta, alpha=0.8) 
     74        shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]] 
     75        draw_shape(ax, size, view, shimmy, alpha=0.8) 
    9576    for v in 'xyz': 
    9677        a, b, c = size 
    … …  
    9980        getattr(ax, v+'axis').label.set_text(v) 
    10081 
    101 def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
    102     """Draw an ellipsoid.""" 
     82def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1): 
    10383    a,b,c = size 
     84    theta, phi, psi = view 
     85    dtheta, dphi, dpsi = shimmy 
     86 
    10487    u = np.linspace(0, 2 * np.pi, steps) 
    10588    v = np.linspace(0, np.pi, steps) 
    … …  
    10790    y = b*np.outer(np.sin(u), np.sin(v)) 
    10891    z = c*np.outer(np.ones_like(u), np.cos(v)) 
    109     x, y, z = transform_xyz(view, jitter, x, y, z) 
     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] 
    11098 
    11199    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    112100 
    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  
    122 def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
    123     """Draw a parallelepiped.""" 
     101def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 
    124102    a,b,c = size 
     103    theta, phi, psi = view 
     104    dtheta, dphi, dpsi = shimmy 
     105 
    125106    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    126107    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    … …  
    137118    ]) 
    138119 
    139     x, y, z = transform_xyz(view, jitter, x, y, z) 
     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] 
    140125    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    141126 
    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  
    151127def draw_sphere(ax, radius=10., steps=100): 
    152     """Draw a sphere""" 
    153128    u = np.linspace(0, 2 * np.pi, steps) 
    154129    v = np.linspace(0, np.pi, steps) 
    … …  
    159134    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    160135 
    161 def 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) 
     136def 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': 
    171146        weights = exp(-0.5*t**2) 
    172     elif dist == 'rectangle': 
    173         # Note: uses sasmodels ridiculous definition of rectangle width 
    174         t = np.linspace(-1, 1, n)*sqrt(3) 
     147    elif dist == 'rect': 
    175148        weights = np.ones_like(t) 
    176149    else: 
    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  
    192 def 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 
     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 
     165def 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 
    210170def Rx(angle): 
    211     """Construct a matrix to rotate points about *x* by *angle* degrees.""" 
    212171    a = radians(angle) 
    213     R = [[1, 0, 0], 
    214          [0, +cos(a), -sin(a)], 
    215          [0, +sin(a), +cos(a)]] 
     172    R = [[1., 0., 0.], 
     173         [0.,  cos(a), sin(a)], 
     174         [0., -sin(a), cos(a)]] 
    216175    return np.matrix(R) 
    217176 
    218177def Ry(angle): 
    219     """Construct a matrix to rotate points about *y* by *angle* degrees.""" 
    220178    a = radians(angle) 
    221     R = [[+cos(a), 0, +sin(a)], 
    222          [0, 1, 0], 
    223          [-sin(a), 0, +cos(a)]] 
     179    R = [[cos(a), 0., -sin(a)], 
     180         [0., 1., 0.], 
     181         [sin(a), 0.,  cos(a)]] 
    224182    return np.matrix(R) 
    225183 
    226184def Rz(angle): 
    227     """Construct a matrix to rotate points about *z* by *angle* degrees.""" 
    228185    a = radians(angle) 
    229     R = [[+cos(a), -sin(a), 0], 
    230          [+sin(a), +cos(a), 0], 
    231          [0, 0, 1]] 
     186    R = [[cos(a), -sin(a), 0.], 
     187         [sin(a),  cos(a), 0.], 
     188         [0., 0., 1.]] 
    232189    return np.matrix(R) 
    233190 
    234 def 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  
    246 def 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  
    256 def 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. 
    268 PD_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  
    279 def 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  
    297 def 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  
    350 def 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  
    388 def 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  
    428 def 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 
     191def main(): 
    453192    #plt.hold(True) 
    454193    plt.set_cmap('gist_earth') 
    … …  
    457196    #ax = plt.subplot(gs[0], projection='3d') 
    458197    ax = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d') 
    459     ax.axis('square') 
     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' 
    460206 
    461207    axcolor = 'lightgoldenrodyellow' 
    462  
    463     ## add control widgets to plot 
    464208    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    465209    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    … …  
    468212    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    469213    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
    470  
    471214    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    472215    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    473216    axdpsi= plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
    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 
     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 
    483221    def update(val, axis=None): 
    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] 
     222        theta, phi, psi = stheta.val, sphi.val, spsi.val 
     223        dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val 
    490224        ax.cla() 
    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) 
     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) 
    496229        plt.gcf().canvas.draw() 
    497230 
    498     ## bind control widgets to view updater 
    499231    stheta.on_changed(lambda v: update(v,'theta')) 
    500232    sphi.on_changed(lambda v: update(v, 'phi')) 
    … …  
    504236    sdpsi.on_changed(lambda v: update(v, 'dpsi')) 
    505237 
    506     ## initialize view 
    507238    update(None, 'phi') 
    508239 
    509     ## go interactive 
    510240    plt.show() 
    511241 
    512242if __name__ == "__main__": 
    513     model_name = sys.argv[1] if len(sys.argv) > 1 else 'parallelepiped' 
    514     main(model_name) 
     243    main() 

  • explore/precision.py

    r2a602c7 r237c9cf  
    11#!/usr/bin/env python 
    22r""" 
    3 Show numerical precision of various expressions. 
    4  
    5 Evaluates the same function(s) in single and double precision and compares 
    6 the results to 500 digit mpmath evaluation of the same function. 
    7  
    8 Note: a quick way to generation C and python code for taylor series 
    9 expansions 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. 
     3Show numerical precision of $2 J_1(x)/x$. 
    304""" 
    315from __future__ import division, print_function 
    … …  
    310284    np_function=scipy.special.erfc, 
    311285    ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]), 
    312     limits=(-5., 5.), 
    313 ) 
    314 add_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"), 
    319286    limits=(-5., 5.), 
    320287) 
    … …  
    481448) 
    482449 
    483 replacement_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 """ 
    507 add_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  
    514450# Alternate versions of 3 j1(x)/x, for posterity 
    515451def taylor_3j1x_x(x): 

  • sasmodels/compare.py

    r31eea1f rced5bd2  
    4242from . import exception 
    4343from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    44 from .direct_model import DirectModel, get_mesh 
     44from .direct_model import DirectModel 
    4545from .convert import revert_name, revert_pars, constrain_new_to_old 
    4646from .generate import FLOAT_RE 
    47 from .weights import plot_weights 
    4847 
    4948try: 
    … …  
    8483    -pars/-nopars* prints the parameter set or not 
    8584    -default/-demo* use demo vs default parameters 
    86     -sphere[=150] set up spherical integration over theta/phi using n points 
    8785 
    8886    === calculation options === 
    89     -mono*/-poly force monodisperse or allow polydisperse random parameters 
     87    -mono*/-poly force monodisperse or allow polydisperse demo parameters 
    9088    -cutoff=1e-5* cutoff value for including a point in polydispersity 
    9189    -magnetic/-nonmagnetic* suppress magnetism 
    … …  
    9492 
    9593    === precision options === 
    96     -engine=default uses the default calcution precision 
     94    -calc=default uses the default calcution precision 
    9795    -single/-double/-half/-fast sets an OpenCL calculation engine 
    9896    -single!/-double!/-quad! sets an OpenMP calculation engine 
    … …  
    105103    -abs/-rel* plot relative or absolute error 
    106104    -title="note" adds note to the plot title, after the model name 
    107     -weights shows weights plots for the polydisperse parameters 
    108105 
    109106    === output options === 
    … …  
    114111vary from 64-bit to 128-bit, with 80-bit floats being common (1e-19 precision). 
    115112On unix and mac you may need single quotes around the DLL computation 
    116 engines, such as -engine='single!,double!' since !, is treated as a history 
     113engines, such as -calc='single!,double!' since !, is treated as a history 
    117114expansion request in the shell. 
    118115 
    … …  
    126123 
    127124    # compare single and double precision calculation for a barbell 
    128     sascomp barbell -engine=single,double 
     125    sascomp barbell -calc=single,double 
    129126 
    130127    # generate 10 random lorentz models, with seed=27 
    … …  
    135132 
    136133    # model timing test requires multiple evals to perform the estimate 
    137     sascomp pringle -engine=single,double -timing=100,100 -noplot 
     134    sascomp pringle -calc=single,double -timing=100,100 -noplot 
    138135""" 
    139136 
    … …  
    281278        # orientation in [-180,180], orientation pd in [0,45] 
    282279        if p.endswith('_pd'): 
    283             return 0., 180. 
     280            return 0., 45. 
    284281        else: 
    285282            return -180., 180. 
    … …  
    436433 
    437434 
    438 def 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  
    448435def constrain_pars(model_info, pars): 
    449436    # type: (ModelInfo, ParameterSet) -> None 
    … …  
    508495    Format the parameter list for printing. 
    509496    """ 
    510     is2d = True 
    511497    lines = [] 
    512498    parameters = model_info.parameters 
    … …  
    829815 
    830816 
    831 def compare(opts, limits=None, maxdim=np.inf): 
     817def compare(opts, limits=None): 
    832818    # type: (Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float] 
    833819    """ 
    … …  
    840826    as the values are adjusted, making it easier to see the effects of the 
    841827    parameters. 
    842  
    843     *maxdim* is the maximum value for any parameter with units of Angstrom. 
    844     """ 
     828    """ 
     829    limits = np.Inf, -np.Inf 
    845830    for k in range(opts['sets']): 
    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) 
     831        opts['pars'] = parse_pars(opts) 
    852832        if opts['pars'] is None: 
    853833            return 
    … …  
    855835        if opts['plot']: 
    856836            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)) 
    863837    if opts['plot']: 
    864838        import matplotlib.pyplot as plt 
    865839        plt.show() 
    866     return limits 
    867840 
    868841def run_models(opts, verbose=False): 
    … …  
    939912 
    940913 
    941 def plot_models(opts, result, limits=None, setnum=0): 
     914def plot_models(opts, result, limits=(np.Inf, -np.Inf), setnum=0): 
    942915    # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float] 
    943     import matplotlib.pyplot as plt 
    944  
    945916    base_value, comp_value = result['base_value'], result['comp_value'] 
    946917    base_time, comp_time = result['base_time'], result['comp_time'] 
    … …  
    954925    # Plot if requested 
    955926    view = opts['view'] 
    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 
     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 
    965936 
    966937    if have_base: 
    … …  
    991962            err[err > cutoff] = cutoff 
    992963        #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 ======= 
    998964        plot_theory(data, None, resid=err, view=view, use_data=use_data) 
    999965        plt.yscale(errview) 
    1000 >>>>>>> master 
    1001966        plt.title("max %s = %.3g"%(errstr, abs(err).max())) 
    1002967        #cbar_title = errstr if errview=="linear" else "log "+errstr 
    … …  
    1030995OPTIONS = [ 
    1031996    # Plotting 
    1032     'plot', 'noplot', 'weights', 
     997    'plot', 'noplot', 
    1033998    'linear', 'log', 'q4', 
    1034999    'rel', 'abs', 
    … …  
    10451010    'demo', 'default',  # TODO: remove demo/default 
    10461011    'nopars', 'pars', 
    1047     'sphere', 'sphere=', # integrate over a sphere in 2d with n points 
    10481012 
    10491013    # Calculation options 
    … …  
    10541018 
    10551019    # Precision options 
    1056     'engine=', 
     1020    'calc=', 
    10571021    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!', 
    10581022    'sasview',  # TODO: remove sasview 3.x support 
    … …  
    11811145        'sets'      : 0, 
    11821146        'engine'    : 'default', 
    1183         'count'     : '1', 
    1184         'show_weights' : False, 
    1185         'sphere'    : 0, 
     1147        'evals'     : '1', 
    11861148    } 
    11871149    for arg in flags: 
    … …  
    12061168        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:] 
    12071169        elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:] 
     1170        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:]) 
    12081171        elif arg.startswith('-title='):    opts['title'] = arg[7:] 
    12091172        elif arg.startswith('-data='):     opts['datafile'] = arg[6:] 
    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 
     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) 
    12211176        elif arg == '-preset':  opts['seed'] = -1 
    12221177        elif arg == '-mono':    opts['mono'] = True 
    … …  
    12411196        elif arg == '-demo':    opts['use_demo'] = True 
    12421197        elif arg == '-default': opts['use_demo'] = False 
    1243         elif arg == '-weights': opts['show_weights'] = True 
    12441198        elif arg == '-html':    opts['html'] = True 
    12451199        elif arg == '-help':    opts['html'] = True 
    … …  
    12781232 
    12791233    if PAR_SPLIT in opts['engine']: 
    1280         opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
     1234        engine_types = opts['engine'].split(PAR_SPLIT, 2) 
    12811235        comparison = True 
    12821236    else: 
    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)] 
     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)] 
    12871241        comparison = True 
    12881242    else: 
    1289         opts['count'] = [int(opts['count'])]*2 
     1243        evals = [int(opts['evals'])]*2 
    12901244 
    12911245    if PAR_SPLIT in opts['cutoff']: 
    1292         opts['cutoff'] = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)] 
     1246        cutoff = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)] 
    12931247        comparison = True 
    12941248    else: 
    1295         opts['cutoff'] = [float(opts['cutoff'])]*2 
    1296  
    1297     base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
     1249        cutoff = [float(opts['cutoff'])]*2 
     1250 
     1251    base = make_engine(model_info[0], data, engine_types[0], cutoff[0]) 
    12981252    if comparison: 
    1299         comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
     1253        comp = make_engine(model_info[1], data, engine_types[1], cutoff[1]) 
    13001254    else: 
    13011255        comp = None 
    … …  
    13061260        'data'      : data, 
    13071261        'name'      : names, 
    1308         'info'      : model_info, 
     1262        'def'       : model_info, 
     1263        'count'     : evals, 
    13091264        'engines'   : [base, comp], 
    13101265        'values'    : values, 
    … …  
    13121267    # pylint: enable=bad-whitespace 
    13131268 
    1314     # Set the integration parameters to the half sphere 
    1315     if opts['sphere'] > 0: 
    1316         set_spherical_integration_parameters(opts, opts['sphere']) 
    1317  
    13181269    return opts 
    13191270 
    1320 def 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  
    1341 def parse_pars(opts, maxdim=np.inf): 
    1342     model_info, model_info2 = opts['info'] 
     1271def parse_pars(opts): 
     1272    model_info, model_info2 = opts['def'] 
    13431273 
    13441274    # Get demo parameters from model definition, or use default parameters 
    … …  
    13511281    if opts['seed'] > -1: 
    13521282        pars = randomize_pars(model_info, pars) 
    1353         limit_dimensions(model_info, pars, maxdim) 
    13541283        if model_info != model_info2: 
    13551284            pars2 = randomize_pars(model_info2, pars2) 
    1356             limit_dimensions(model_info, pars2, maxdim) 
    13571285            # Share values for parameters with the same name 
    13581286            for k, v in pars.items(): 
    … …  
    14311359    from . import rst2html 
    14321360 
    1433     info = opts['info'][0] 
     1361    info = opts['def'][0] 
    14341362    html = make_html(info) 
    14351363    path = os.path.dirname(info.filename) 
    … …  
    14821410        opts['pars'] = list(opts['pars']) 
    14831411        p1, p2 = opts['pars'] 
    1484         m1, m2 = opts['info'] 
     1412        m1, m2 = opts['def'] 
    14851413        self.fix_p2 = m1 != m2 or p1 != p2 
    14861414        model_info = m1 
    … …  
    15011429        self.starting_values = dict((k, v.value) for k, v in pars.items()) 
    15021430        self.pd_types = pd_types 
    1503         self.limits = None 
     1431        self.limits = np.Inf, -np.Inf 
    15041432 
    15051433    def revert_values(self): 

  • sasmodels/core.py

    r9e771a3 r60335cc  
    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 
    275     guaranteed by the OpenCL standard.  'default' will choose the appropriate 
    276     default for the model and platform. 
     274    fast native functions rather than those with the precision level guaranteed 
     275    by the OpenCL standard. 
    277276 
    278277    Platform preference can be specfied ("ocl" vs "dll"), with the default 

  • sasmodels/data.py

    re3571cb r09141ff  
    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') and not getattr(data, 'radial', False): 
     365    elif hasattr(data, 'qx_data'): 
    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') and not getattr(data, 'radial', False): 
     393    elif hasattr(data, 'qx_data'): 
    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 
    431427 
    432428    use_data = use_data and data.y is not None 
    … …  
    655651    ymin, ymax = min(data.qy_data), max(data.qy_data) 
    656652    if view == 'log': 
    657         vmin_scaled, vmax_scaled= np.log10(vmin), np.log10(vmax) 
    658     else: 
    659         vmin_scaled, vmax_scaled = vmin, vmax 
     653        vmin, vmax = np.log10(vmin), np.log10(vmax) 
    660654    plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)), 
    661655               interpolation='nearest', aspect=1, origin='lower', 
    662                extent=[xmin, xmax, ymin, ymax], 
    663                vmin=vmin_scaled, vmax=vmax_scaled) 
     656               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
    664657    plt.xlabel("$q_x$/A$^{-1}$") 
    665658    plt.ylabel("$q_y$/A$^{-1}$") 

  • sasmodels/details.py

    r2bccb5a rccd5f01  
    1515 
    1616import numpy as np  # type: ignore 
    17 from numpy import pi, cos, sin, radians 
     17from numpy import pi, cos, sin 
    1818 
    1919try: 
    … …  
    2929 
    3030try: 
    31     from typing import List, Tuple, Sequence 
     31    from typing import List 
    3232except ImportError: 
    3333    pass 
    3434else: 
    3535    from .modelinfo import ModelInfo 
    36     from .modelinfo import ParameterTable 
    3736 
    3837 
    … …  
    5453    coordinates, the total circumference decreases as latitude varies from 
    5554    pi r^2 at the equator to 0 at the pole, and the weight associated 
    56     with a range of latitude values needs to be scaled by this circumference. 
     55    with a range of phi values needs to be scaled by this circumference. 
    5756    This scale factor needs to be updated each time the theta value 
    5857    changes.  *theta_par* indicates which of the values in the parameter 
    … …  
    219218 
    220219ZEROS = tuple([0.]*31) 
    221 def make_kernel_args(kernel, mesh): 
     220def make_kernel_args(kernel, pairs): 
    222221    # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool] 
    223222    """ 
    224     Converts (value, dispersity, weight) for each parameter into kernel pars. 
     223    Converts (value, weight) pairs into parameters for the kernel call. 
    225224 
    226225    Returns a CallDetails object indicating the polydispersity, a data object 
    … …  
    231230    npars = kernel.info.parameters.npars 
    232231    nvalues = kernel.info.parameters.nvalues 
    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) 
     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 ((),()) 
    237234    length = np.array([len(w) for w in weights]) 
    238235    offset = np.cumsum(np.hstack((0, length))) 
    239236    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
    240237    # Pad value array to a 32 value boundaryd 
    241     data_len = nvalues + 2*sum(len(v) for v in dispersity) 
     238    data_len = nvalues + 2*sum(len(v) for v in values) 
    242239    extra = (32 - data_len%32)%32 
    243     data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra]) 
     240    data = np.hstack((scalars,) + values + weights + ZEROS[:extra]) 
    244241    data = data.astype(kernel.dtype) 
    245242    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
    … …  
    247244    return call_details, data, is_magnetic 
    248245 
    249 def 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  
    279246 
    280247def convert_magnetism(parameters, values): 
    281     # type: (ParameterTable, Sequence[np.ndarray]) -> bool 
    282248    """ 
    283249    Convert magnetism values from polar to rectangular coordinates. 
    … …  
    289255    scale = mag[:,0] 
    290256    if np.any(scale): 
    291         theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
     257        theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. 
    292258        cos_theta = cos(theta) 
    293259        mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
    … …  
    299265 
    300266 
    301 def dispersion_mesh(model_info, mesh): 
     267def dispersion_mesh(model_info, pars): 
    302268    # type: (ModelInfo) -> Tuple[List[np.ndarray], List[np.ndarray]] 
    303269    """ 
    304270    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. 
    315271 
    316272    Returns [p1,p2,...],w where pj is a vector of values for parameter j 
    … …  
    318274    parameter set in the vector. 
    319275    """ 
    320     value, dispersity, weight = zip(*mesh) 
     276    value, weight = zip(*pars) 
    321277    #weight = [w if len(w)>0 else [1.] for w in weight] 
    322278    weight = np.vstack([v.flatten() for v in meshgrid(*weight)]) 
    323279    weight = np.prod(weight, axis=0) 
    324     dispersity = [v.flatten() for v in meshgrid(*dispersity)] 
     280    value = [v.flatten() for v in meshgrid(*value)] 
    325281    lengths = [par.length for par in model_info.parameters.kernel_parameters 
    326282               if par.type == 'volume'] 
    … …  
    329285        offset = 0 
    330286        for n in lengths: 
    331             pars.append(np.vstack(dispersity[offset:offset+n]) 
    332                         if n > 1 else dispersity[offset]) 
     287            pars.append(np.vstack(value[offset:offset+n]) 
     288                        if n > 1 else value[offset]) 
    333289            offset += n 
    334         dispersity = pars 
    335     return dispersity, weight 
     290        value = pars 
     291    return value, weight 

  • sasmodels/direct_model.py

    r9e771a3 rd1ff3a5  
    5555    *mono* is True if polydispersity should be set to none on all parameters. 
    5656    """ 
    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) 
     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) 
    6073    #print("values:", values) 
    6174    return calculator(call_details, values, cutoff, is_magnetic) 
     75 
    6276 
    6377def call_ER(model_info, pars): 
    … …  
    115129    return x, y, model_info.profile_axes 
    116130 
    117 def 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  
    142 def _get_par_weights(parameter, values, active=True): 
    143     # type: (Parameter, Dict[str, float]) -> Tuple[float, np.ndarray, np.ndarray] 
     131 
     132def get_weights(parameter, values): 
     133    # type: (Parameter, Dict[str, float]) -> Tuple[np.ndarray, np.ndarray] 
    144134    """ 
    145135    Generate the distribution for parameter *name* given the parameter values 
    … …  
    150140    """ 
    151141    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') 
    152145    npts = values.get(parameter.name+'_pd_n', 0) 
    153146    width = values.get(parameter.name+'_pd', 0.0) 
    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  
    169 def _vol_pars(model_info, values): 
     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 
     155def _vol_pars(model_info, pars): 
    170156    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray] 
    171     vol_pars = [_get_par_weights(p, values) 
     157    vol_pars = [get_weights(p, pars) 
    172158                for p in model_info.parameters.call_parameters 
    173159                if p.type == 'volume'] 
    174160    #import pylab; pylab.plot(vol_pars[0][0],vol_pars[0][1]); pylab.show() 
    175     dispersity, weight = dispersion_mesh(model_info, vol_pars) 
    176     return dispersity, weight 
     161    value, weight = dispersion_mesh(model_info, vol_pars) 
     162    return value, weight 
    177163 
    178164 

  • sasmodels/generate.py

    r6773b02 r30b60d2  
    167167import string 
    168168from zlib import crc32 
    169 from inspect import currentframe, getframeinfo 
    170169 
    171170import numpy as np  # type: ignore 
    … …  
    371370    """ 
    372371    # 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 
    377372    return "%08X"%(0xffffffff&crc32(source)) 
    378373 
    … …  
    690685    # Load templates and user code 
    691686    kernel_header = load_template('kernel_header.c') 
    692     kernel_code = load_template('kernel_iq.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') 
    693690    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') 
    699691 
    700692    # Build initial sources 
    … …  
    721713 
    722714    # Define the parameter table 
    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') 
     715    # TODO: plug in current line number 
     716    source.append('#line 542 "sasmodels/generate.py"') 
    726717    source.append("#define PARAMETER_TABLE \\") 
    727718    source.append("\\\n".join(p.as_definition() 
    … …  
    739730    source.append(call_volume) 
    740731 
    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" 
     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)) 
    756744 
    757745    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    … …  
    764752    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic) 
    765753    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)) 
    766756 
    767757    # TODO: allow mixed python/opencl kernels? 
    768758 
    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) 
     759    ocl = kernels(ocl_code, call_iq, call_iqxy, model_info.name) 
     760    dll = kernels(dll_code, call_iq, call_iqxy, model_info.name) 
    771761    result = { 
    772762        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    … …  
    777767 
    778768 
    779 def kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
     769def kernels(kernel, call_iq, call_iqxy, name): 
    780770    # type: ([str,str], str, str, str) -> List[str] 
    781771    code = kernel[0] 
    … …  
    797787        '#line 1 "%s Iqxy"' % path, 
    798788        code, 
    799         clear_iqxy, 
     789        "#undef CALL_IQ", 
    800790        "#undef KERNEL_NAME", 
    801791    ] 
    … …  
    808798        '#line 1 "%s Imagnetic"' % path, 
    809799        code, 
    810         clear_iqxy, 
    811800        "#undef MAGNETIC", 
     801        "#undef CALL_IQ", 
    812802        "#undef KERNEL_NAME", 
    813803    ] 

  • sasmodels/kernel_header.c

    r8698a0d rbb4b509  
    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). 
    9289   static SAS_DOUBLE expm1(SAS_DOUBLE x_in) { 
    9390      double x = (double)x_in;  // go back to float for single precision kernels 
    … …  
    150147inline double cube(double x) { return x*x*x; } 
    151148inline 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

    r2c108a3 rbde38b5  
     1 
    12/* 
    23    ########################################################## 
    … …  
    1112*/ 
    1213 
    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  
    3514#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
    3615#define _PAR_BLOCK_ 
    … …  
    3817typedef struct { 
    3918#if MAX_PD > 0 
    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 
     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 
    4221    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector 
    4322    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level 
    … …  
    4625    int32_t num_weights;        // total length of the weights vector 
    4726    int32_t num_active;         // number of non-trivial pd loops 
    48     int32_t theta_par;          // id of first orientation variable 
     27    int32_t theta_par;          // id of spherical correction variable 
    4928} ProblemDetails; 
    5029 
    … …  
    5938#endif // _PAR_BLOCK_ 
    6039 
    61 #if defined(MAGNETIC) && NUM_MAGNETIC > 0 
    62 // ===== Helper functions for magnetism ===== 
     40 
     41#if defined(MAGNETIC) && NUM_MAGNETIC>0 
    6342 
    6443// Return value restricted between low and high 
    … …  
    7251//     uu * (sld - m_sigma_x); 
    7352//     dd * (sld + m_sigma_x); 
    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 
    77 static void set_spin_weights(double in_spin, double out_spin, double spins[4]) 
     53//     ud * (m_sigma_y + 1j*m_sigma_z); 
     54//     du * (m_sigma_y - 1j*m_sigma_z); 
     55static void set_spins(double in_spin, double out_spin, double spins[4]) 
    7856{ 
    7957  in_spin = clip(in_spin, 0.0, 1.0); 
    8058  out_spin = clip(out_spin, 0.0, 1.0); 
    8159  spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd 
    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 
     60  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du 
     61  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud 
    8462  spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu 
    85   spins[4] = spins[1]; // du imag 
    86   spins[5] = spins[2]; // ud imag 
    8763} 
    8864 
    89 // Compute the magnetic sld 
    90 static 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 ) 
     65static double mag_sld(double qx, double qy, double p, 
     66                       double mx, double my, double sld) 
    9767{ 
    98   if (xs < 4) { 
    9968    const double perp = qy*mx - qx*my; 
    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   } 
     69    return sld + perp*p; 
    11770} 
    11871 
    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  
    130 typedef 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]' 
    137 static void 
    138 qac_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]' 
    170 static double 
    171 qac_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  
    188 typedef 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]' 
    197 static void 
    198 qabc_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]' 
    243 static double 
    244 qabc_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 ======================== 
     72#endif // MAGNETIC 
    25873 
    25974kernel 
    26075void KERNEL_NAME( 
    26176    int32_t nq,                 // number of q values 
    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 
     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 
    26479    global const ProblemDetails *details, 
    26580    global const double *values, 
    26681    global const double *q, // nq q values, with padding to boundary 
    26782    global double *result,  // nq+1 return values, again with padding 
    268     const double cutoff     // cutoff in the dispersity weight product 
     83    const double cutoff     // cutoff in the polydispersity weight product 
    26984    ) 
    27085{ 
    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  
    28286  // Storage for the current parameter values.  These will be updated as we 
    283   // walk the dispersity mesh. 
     87  // walk the polydispersity cube. 
    28488  ParameterBlock local_values; 
    28589 
    … …  
    28791  // Location of the sld parameters in the parameter vector. 
    28892  // These parameters are updated with the effective sld due to magnetism. 
     93  #if NUM_MAGNETIC > 3 
    28994  const int32_t slds[] = { MAGNETIC_PARS }; 
    290  
     95  #endif 
     96 
     97  // TODO: could precompute these outside of the kernel. 
    29198  // Interpret polarization cross section. 
    29299  //     up_frac_i = values[NUM_PARS+2]; 
    293100  //     up_frac_f = values[NUM_PARS+3]; 
    294101  //     up_angle = values[NUM_PARS+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 
     102  double spins[4]; 
    297103  double cos_mspin, sin_mspin; 
    298   set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
     104  set_spins(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
    299105  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
    300106#endif // MAGNETIC 
    … …  
    308114  for (int i=0; i < NUM_PARS; i++) { 
    309115    local_values.vector[i] = values[2+i]; 
    310 //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]); 
    311   } 
    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 /* 
    342 Based on the level of the loop, uses C preprocessor magic to construct 
    343 level-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  
    355 After 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; 
    407   } 
    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. 
     116//printf("p%d = %g\n",i, local_values.vector[i]); 
     117  } 
     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; 
     127  } 
     128//printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
     129 
    506130#if MAX_PD>0 
    507131  global const double *pd_value = values + NUM_VALUES; 
    … …  
    509133#endif 
    510134 
    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. 
    514 int step = pd_start; 
    515  
    516 // define looping variables 
     135  // Jump into the middle of the polydispersity loop 
    517136#if MAX_PD>4 
    518   PD_INIT(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]; 
    519142#endif 
    520143#if MAX_PD>3 
    521   PD_INIT(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); 
    522150#endif 
    523151#if MAX_PD>2 
    524   PD_INIT(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]; 
    525157#endif 
    526158#if MAX_PD>1 
    527   PD_INIT(1) 
    528 #endif 
    529 #if MAX_PD>0 
    530   PD_INIT(0) 
    531 #endif 
    532  
    533 // open nested loops 
    534 PD_OUTERMOST_WEIGHT(MAX_PD) 
     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]; 
     164#endif 
     165#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 
     175#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 
    535188#if MAX_PD>4 
    536   PD_OPEN(4,5) 
     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; 
    537196#endif 
    538197#if MAX_PD>3 
    539   PD_OPEN(3,4) 
     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; 
    540204#endif 
    541205#if MAX_PD>2 
    542   PD_OPEN(2,3) 
     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; 
    543212#endif 
    544213#if MAX_PD>1 
    545   PD_OPEN(1,2) 
    546 #endif 
    547 #if MAX_PD>0 
    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. 
     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]; 
     258          const double qsq = qx*qx + qy*qy; 
     259 
     260          // Constant across orientation, polydispersity for given qx, qy 
    580261          double scattering = 0.0; 
    581           const double qsq = qx*qx + qy*qy; 
     262          // TODO: what is the magnetic scattering at q=0 
    582263          if (qsq > 1.e-16) { 
    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); 
     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); 
    602298                } 
    603                 scattering += xs_weight * CALL_KERNEL(); 
    604299              } 
    605300            } 
    606301          } 
    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 
     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        } 
    617308      } 
    618309    } 
    619   } 
    620  
    621 // close nested loops 
    622 ++step; 
    623 #if MAX_PD>0 
    624   PD_CLOSE(0) 
     310    ++step; 
     311#if MAX_PD>0 
     312    if (step >= pd_stop) break; 
     313    ++i0; 
     314  } 
     315  i0 = 0; 
    625316#endif 
    626317#if MAX_PD>1 
    627   PD_CLOSE(1) 
     318    if (step >= pd_stop) break; 
     319    ++i1; 
     320  } 
     321  i1 = 0; 
    628322#endif 
    629323#if MAX_PD>2 
    630   PD_CLOSE(2) 
     324    if (step >= pd_stop) break; 
     325    ++i2; 
     326  } 
     327  i2 = 0; 
    631328#endif 
    632329#if MAX_PD>3 
    633   PD_CLOSE(3) 
     330    if (step >= pd_stop) break; 
     331    ++i3; 
     332  } 
     333  i3 = 0; 
    634334#endif 
    635335#if MAX_PD>4 
    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 
     336    if (step >= pd_stop) break; 
     337    ++i4; 
     338  } 
     339  i4 = 0; 
     340#endif 
     341 
     342//printf("res: %g/%g\n", result[0], pd_norm); 
     343  // Remember the updated norm. 
    645344  result[nq] = pd_norm; 
    646 //printf("res: %g/%g\n", result[0], pd_norm); 
    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 
    657345} 

  • sasmodels/kernelpy.py

    r8698a0d r3b32f8d  
    113113 
    114114        partable = model_info.parameters 
    115         #kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 
    116         #                     else partable.iq_parameters) 
    117         kernel_parameters = partable.iq_parameters 
     115        kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 
     116                             else partable.iq_parameters) 
    118117        volume_parameters = partable.form_volume_parameters 
    119118 
    … …  
    202201 
    203202    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) 
    208209    p0_length = call_details.pd_length[0] 
    209210    p0_index = p0_length 
    … …  
    222223            parameters[pd_par] = pd_value[pd_offset+pd_index] 
    223224            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) 
    224228            p0_index = loop_index%p0_length 
    225229 
    226230        weight = partial_weight * pd_weight[p0_offset + p0_index] 
    227231        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) 
    228235        p0_index += 1 
    229236        if weight > cutoff: 
    … …  
    237244 
    238245            # update value and norm 
     246            weight *= spherical_correction 
    239247            total += weight * Iq 
    240248            pd_norm += weight * form_volume() 

  • sasmodels/model_test.py

    r65314f7 r65314f7  
    8686    suite = unittest.TestSuite() 
    8787 
    88     if models[0] in core.KINDS: 
     88    if models[0] == 'all': 
    8989        skip = models[1:] 
    90         models = list_models(models[0]) 
     90        models = list_models() 
    9191    else: 
    9292        skip = [] 

  • sasmodels/modelinfo.py

    re3571cb r65314f7  
    382382      with vector parameter p sent as p[]. 
    383383 
    384     * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 
     384    * *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() 
    423422        self._set_vector_lengths() 
    424423 
    … …  
    439438        self.iq_parameters = [p for p in self.kernel_parameters 
    440439                              if p.type not in ('orientation', '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'] 
     440        self.iqxy_parameters = [p for p in self.kernel_parameters 
     441                                if p.type != 'magnetic'] 
    444442        self.form_volume_parameters = [p for p in self.kernel_parameters 
    445443                                       if p.type == 'volume'] 
    … …  
    463461        self.has_2d = any(p.type in ('orientation', 'magnetic') 
    464462                          for p in self.kernel_parameters) 
    465         self.is_asymmetric = any(p.name == 'psi' for p in self.kernel_parameters) 
    466463        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    467464                                if p.id.startswith('M0:')] 
    … …  
    470467                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    471468        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") 
    495469 
    496470    def __getitem__(self, key): 
    … …  
    502476            raise KeyError("unknown parameter %r"%key) 
    503477        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 
    511478 
    512479    def _set_vector_lengths(self): 

  • sasmodels/models/barbell.c

    rbecded3 r592343f  
     1double form_volume(double radius_bell, double radius, double length); 
     2double Iq(double q, double sld, double solvent_sld, 
     3        double radius_bell, double radius, double length); 
     4double Iqxy(double qx, double qy, double sld, double solvent_sld, 
     5        double radius_bell, double radius, double length, 
     6        double theta, double phi); 
     7 
    18#define INVALID(v) (v.radius_bell < v.radius) 
    29 
    310//barbell kernel - same as dumbell 
    411static double 
    5 _bell_kernel(double qab, double qc, double h, double radius_bell, 
    6              double half_length) 
     12_bell_kernel(double q, double h, double radius_bell, 
     13             double half_length, double sin_alpha, double cos_alpha) 
    714{ 
    815    // translate a point in [-1,1] to a point in [lower,upper] 
    … …  
    1926    //    m = q R cos(alpha) 
    2027    //    b = q(L/2-h) cos(alpha) 
    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) 
     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) 
    2431    double total = 0.0; 
    2532    for (int i = 0; i < 76; i++){ 
    2633        const double t = Gauss76Z[i]*zm + zb; 
    2734        const double radical = 1.0 - t*t; 
    28         const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
     35        const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
    2936        const double Fq = cos(m*t + b) * radical * bj; 
    3037        total += Gauss76Wt[i] * Fq; 
    … …  
    3744 
    3845static double 
    39 _fq(double qab, double qc, double h, 
    40     double radius_bell, double radius, double half_length) 
     46_fq(double q, double h, 
     47    double radius_bell, double radius, double half_length, 
     48    double sin_alpha, double cos_alpha) 
    4149{ 
    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); 
     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); 
    4553    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    4654    const double Aq = bell_fq + cyl_fq; 
    … …  
    4856} 
    4957 
    50 static double 
    51 form_volume(double radius_bell, 
    52     double radius, 
    53     double length) 
     58 
     59double form_volume(double radius_bell, 
     60        double radius, 
     61        double length) 
    5462{ 
    5563    // bell radius should never be less than radius when this is called 
    … …  
    6270} 
    6371 
    64 static double 
    65 Iq(double q, double sld, double solvent_sld, 
    66     double radius_bell, double radius, double length) 
     72double Iq(double q, double sld, double solvent_sld, 
     73          double radius_bell, double radius, double length) 
    6774{ 
    6875    const double h = -sqrt(radius_bell*radius_bell - radius*radius); 
    … …  
    7784        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    7885        SINCOS(alpha, sin_alpha, cos_alpha); 
    79         const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
     86        const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 
    8087        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8188    } 
    … …  
    8996 
    9097 
    91 static double 
    92 Iqxy(double qab, double qc, 
    93     double sld, double solvent_sld, 
    94     double radius_bell, double radius, double length) 
     98double Iqxy(double qx, double qy, 
     99        double sld, double solvent_sld, 
     100        double radius_bell, double radius, double length, 
     101        double theta, double phi) 
    95102{ 
     103    double q, sin_alpha, cos_alpha; 
     104    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     105 
    96106    const double h = -sqrt(square(radius_bell) - square(radius)); 
    97     const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length); 
     107    const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 
    98108 
    99109    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/bcc_paracrystal.c

    rea60e08 r50beefe  
    1 static double 
    2 bcc_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; 
     1double form_volume(double radius); 
     2double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
     3double 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); 
    86 
    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)); 
     7double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
     8double _BCCeval(double Theta, double Phi, double temp1, double temp3); 
     9double _sphereform(double q, double radius, double sld, double solvent_sld); 
    2510 
    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 
    5211 
    53     return Zq; 
     12double _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); 
    5420} 
    5521 
     22double _BCCeval(double Theta, double Phi, double temp1, double temp3) { 
    5623 
    57 // occupied volume fraction calculated from lattice symmetry and sphere radius 
    58 static double 
    59 bcc_volume_fraction(double radius, double dnn) 
    60 { 
    61     return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
     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); 
    6241} 
    6342 
    64 static double 
    65 form_volume(double radius) 
    66 { 
     43double form_volume(double radius){ 
    6744    return sphere_volume(radius); 
    6845} 
    6946 
    7047 
    71 static 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; 
     48double Iq(double q, double dnn, 
     49  double d_factor, double radius, 
     50  double sld, double solvent_sld){ 
    8151 
    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; 
     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; 
    10684} 
    10785 
    10886 
    109 static double Iqxy(double qa, double qb, double qc, 
     87double Iqxy(double qx, double qy, 
    11088    double dnn, double d_factor, double radius, 
    111     double sld, double solvent_sld) 
     89    double sld, double solvent_sld, 
     90    double theta, double phi, double psi) 
    11291{ 
    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; 
     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; 
    117111} 

  • sasmodels/models/capped_cylinder.c

    rbecded3 r592343f  
     1double form_volume(double radius, double radius_cap, double length); 
     2double Iq(double q, double sld, double solvent_sld, 
     3    double radius, double radius_cap, double length); 
     4double Iqxy(double qx, double qy, double sld, double solvent_sld, 
     5    double radius, double radius_cap, double length, double theta, double phi); 
     6 
    17#define INVALID(v) (v.radius_cap < v.radius) 
    28 
    … …  
    814//   radius_cap is the radius of the lens 
    915//   length is the cylinder length, or the separation between the lens halves 
    10 //   theta is the angle of the cylinder wrt q. 
     16//   alpha is the angle of the cylinder wrt q. 
    1117static double 
    12 _cap_kernel(double qab, double qc, double h, double radius_cap, 
    13     double half_length) 
     18_cap_kernel(double q, double h, double radius_cap, 
     19                      double half_length, double sin_alpha, double cos_alpha) 
    1420{ 
    1521    // translate a point in [-1,1] to a point in [lower,upper] 
    … …  
    2026 
    2127    // cos term in integral is: 
    22     //    cos (q (R t - h + L/2) cos(theta)) 
     28    //    cos (q (R t - h + L/2) cos(alpha)) 
    2329    // so turn it into: 
    2430    //    cos (m t + b) 
    2531    // where: 
    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) 
     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) 
    3137    double total = 0.0; 
    3238    for (int i=0; i<76 ;i++) { 
    3339        const double t = Gauss76Z[i]*zm + zb; 
    3440        const double radical = 1.0 - t*t; 
    35         const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
     41        const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
    3642        const double Fq = cos(m*t + b) * radical * bj; 
    3743        total += Gauss76Wt[i] * Fq; 
    … …  
    4450 
    4551static double 
    46 _fq(double qab, double qc, double h, double radius_cap, double radius, double half_length) 
     52_fq(double q, double h, double radius_cap, double radius, double half_length, 
     53    double sin_alpha, double cos_alpha) 
    4754{ 
    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); 
     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); 
    5158    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5259    const double Aq = cap_Fq + cyl_Fq; 
    … …  
    5461} 
    5562 
    56 static double 
    57 form_volume(double radius, double radius_cap, double length) 
     63double form_volume(double radius, double radius_cap, double length) 
    5864{ 
    5965    // cap radius should never be less than radius when this is called 
    … …  
    8490} 
    8591 
    86 static double 
    87 Iq(double q, double sld, double solvent_sld, 
    88     double radius, double radius_cap, double length) 
     92double Iq(double q, double sld, double solvent_sld, 
     93          double radius, double radius_cap, double length) 
    8994{ 
    9095    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    … …  
    96101    double total = 0.0; 
    97102    for (int i=0; i<76 ;i++) { 
    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; 
     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; 
    106110    } 
    107111    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    114118 
    115119 
    116 static double 
    117 Iqxy(double qab, double qc, 
     120double Iqxy(double qx, double qy, 
    118121    double sld, double solvent_sld, double radius, 
    119     double radius_cap, double length) 
     122    double radius_cap, double length, 
     123    double theta, double phi) 
    120124{ 
     125    double q, sin_alpha, cos_alpha; 
     126    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     127 
    121128    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    122     const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length); 
     129    const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 
    123130 
    124131    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/core_shell_bicelle.c

    rbecded3 rb260926  
    1 static double 
    2 form_volume(double radius, double thick_rim, double thick_face, double length) 
     1double form_volume(double radius, double thick_rim, double thick_face, double length); 
     2double 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 
     13double 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 
     26double form_volume(double radius, double thick_rim, double thick_face, double length) 
    327{ 
    4     return M_PI*square(radius+thick_rim)*(length+2.0*thick_face); 
     28    return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face); 
    529} 
    630 
    731static double 
    8 bicelle_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) 
     32bicelle_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) 
    1843{ 
    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); 
     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); 
    2550 
    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); 
     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); 
    3055 
    3156    const double t = vol1*dr1*si1*be1 + 
    … …  
    3358                     vol3*dr3*si2*be1; 
    3459 
    35     return t; 
     60    const double retval = t*t; 
     61 
     62    return retval; 
     63 
    3664} 
    3765 
    3866static double 
    39 Iq(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) 
     67bicelle_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) 
    4876{ 
    4977    // set up the integration end points 
    … …  
    5179    const double halflength = 0.5*length; 
    5280 
    53     double total = 0.0; 
     81    double summ = 0.0; 
    5482    for(int i=0;i<N_POINTS_76;i++) { 
    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; 
     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; 
    6190    } 
    6291 
    6392    // calculate value of integral to return 
    64     double answer = total*uplim; 
     93    double answer = uplim*summ; 
     94    return answer; 
     95} 
     96 
     97static double 
     98bicelle_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) 
     109{ 
     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, 
     114                           0.5*length, core_sld, face_sld, rim_sld, 
     115                           solvent_sld, sin_alpha, cos_alpha); 
    65116    return 1.0e-4*answer; 
    66117} 
    67118 
    68 static double 
    69 Iqxy(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) 
     119double 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) 
    78128{ 
    79     double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face, 
    80                            0.5*length, core_sld, face_sld, rim_sld, 
    81                            solvent_sld); 
    82     return 1.0e-4*fq*fq; 
     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; 
    83132} 
     133 
     134 
     135double 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

    rbecded3 rdedcf34  
    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 sld_core, 
    20     double sld_face, 
    21     double sld_rim, 
    22     double sld_solvent) 
     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) 
    2323{ 
     24    double si1,si2,be1,be2; 
    2425     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 
    2526     // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 
     27     //    const double uplim = M_PI_4; 
    2628    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 
    2735    const double r_major = r_minor * x_core; 
    2836    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
    2937    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
    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); 
     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); 
    3644 
    3745    //initialize integral 
    … …  
    3947    for(int i=0;i<76;i++) { 
    4048        //setup inner integral over the ellipsoidal cross-section 
    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; 
     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); 
    5158        for(int j=0;j<76;j++) { 
    5259            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    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; 
     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); 
    6470        } 
    6571        //now calculate outer integral 
    … …  
    7177 
    7278static double 
    73 Iqxy(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) 
     79Iqxy(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) 
    8392{ 
    84     const double dr1 = sld_core-sld_face; 
    85     const double dr2 = sld_rim-sld_solvent; 
    86     const double dr3 = sld_face-sld_rim; 
     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; 
    8799    const double r_major = r_minor*x_core; 
    88100    const double halfheight = 0.5*length; 
    … …  
    92104 
    93105    // Compute effective radius in rotated coordinates 
    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; 
     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; 
    103115} 
     116 

  • sasmodels/models/core_shell_cylinder.c

    rbecded3 r592343f  
     1double form_volume(double radius, double thickness, double length); 
     2double Iq(double q, double core_sld, double shell_sld, double solvent_sld, 
     3    double radius, double thickness, double length); 
     4double 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 
    17// vd = volume * delta_rho 
    2 // besarg = q * R * sin(theta) 
    3 // siarg = q * L/2 * cos(theta) 
    4 static double _cyl(double vd, double besarg, double siarg) 
     8// besarg = q * R * sin(alpha) 
     9// siarg = q * L/2 * cos(alpha) 
     10double _cyl(double vd, double besarg, double siarg); 
     11double _cyl(double vd, double besarg, double siarg) 
    512{ 
    613    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 
    714} 
    815 
    9 static double 
    10 form_volume(double radius, double thickness, double length) 
     16double form_volume(double radius, double thickness, double length) 
    1117{ 
    12     return M_PI*square(radius+thickness)*(length+2.0*thickness); 
     18    return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
    1319} 
    1420 
    15 static double 
    16 Iq(double q, 
     21double Iq(double q, 
    1722    double core_sld, 
    1823    double shell_sld, 
    … …  
    2328{ 
    2429    // precalculate constants 
    25     const double core_r = radius; 
    26     const double core_h = 0.5*length; 
     30    const double core_qr = q*radius; 
     31    const double core_qh = q*0.5*length; 
    2732    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    28     const double shell_r = (radius + thickness); 
    29     const double shell_h = (0.5*length + thickness); 
     33    const double shell_qr = q*(radius + thickness); 
     34    const double shell_qh = q*(0.5*length + thickness); 
    3035    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3136    double total = 0.0; 
     37    // double lower=0, upper=M_PI_2; 
    3238    for (int i=0; i<76 ;i++) { 
    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; 
     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; 
    4347    } 
    4448    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4852 
    4953 
    50 double Iqxy(double qab, double qc, 
     54double Iqxy(double qx, double qy, 
    5155    double core_sld, 
    5256    double shell_sld, 
    … …  
    5458    double radius, 
    5559    double thickness, 
    56     double length) 
     60    double length, 
     61    double theta, 
     62    double phi) 
    5763{ 
    58     const double core_r = radius; 
    59     const double core_h = 0.5*length; 
     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; 
    6069    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    61     const double shell_r = (radius + thickness); 
    62     const double shell_h = (0.5*length + thickness); 
     70    const double shell_qr = q*(radius + thickness); 
     71    const double shell_qh = q*(0.5*length + thickness); 
    6372    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    6473 
    65     const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
    66         + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
     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); 
    6776    return 1.0e-4 * fq * fq; 
    6877} 

  • sasmodels/models/core_shell_ellipsoid.c

    rbecded3 r0a3d9b2  
     1double form_volume(double radius_equat_core, 
     2                   double polar_core, 
     3                   double equat_shell, 
     4                   double polar_shell); 
     5double 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); 
    113 
    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> 
    9 static 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; 
    1914 
    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; 
     15double 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); 
    2425 
    25     return fq_core + fq_shell; 
    26 } 
    2726 
    28 static double 
    29 form_volume(double radius_equat_core, 
    30     double x_core, 
    31     double thick_shell, 
    32     double x_polar_shell) 
     27double form_volume(double radius_equat_core, 
     28                   double x_core, 
     29                   double thick_shell, 
     30                   double x_polar_shell) 
    3331{ 
    3432    const double equat_shell = radius_equat_core + thick_shell; 
    … …  
    3937 
    4038static double 
    41 Iq(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) 
     39core_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) 
    4947{ 
    50     const double sld_core_shell = core_sld - shell_sld; 
    51     const double sld_shell_solvent = shell_sld - solvent_sld; 
     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 
    5255 
    5356    const double polar_core = radius_equat_core*x_core; 
    … …  
    5558    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    5659 
    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 
     60    double summ = 0.0;   //initialize intergral 
    6161    for(int i=0;i<76;i++) { 
    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; 
     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; 
    6966    } 
    70     total *= m; 
     67    summ *= 0.5*(uplim-lolim); 
    7168 
    7269    // convert to [cm-1] 
    73     return 1.0e-4 * total; 
     70    return 1.0e-4 * summ; 
    7471} 
    7572 
    7673static double 
    77 Iqxy(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) 
     74core_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) 
    8584{ 
    86     const double sld_core_shell = core_sld - shell_sld; 
    87     const double sld_shell_solvent = shell_sld - solvent_sld; 
     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; 
    8890 
    8991    const double polar_core = radius_equat_core*x_core; 
    … …  
    9193    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    9294 
    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); 
     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); 
    97105 
    98106    //convert to [cm-1] 
    99     return 1.0e-4 * fq * fq; 
     107    answer *= 1.0e-4; 
     108 
     109    return answer; 
    100110} 
     111 
     112double 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 
     134double 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/gauss76.c", "core_shell_ellipsoid.c"] 
     143source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c", 
     144          "core_shell_ellipsoid.c"] 
    144145 
    145146def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 

  • sasmodels/models/core_shell_parallelepiped.c

    rbecded3 r92dfe0c  
    1 static double 
    2 form_volume(double length_a, double length_b, double length_c, 
    3     double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     1double form_volume(double length_a, double length_b, double length_c,  
     2                   double thick_rim_a, double thick_rim_b, double thick_rim_c); 
     3double 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); 
     6double 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 
     11double form_volume(double length_a, double length_b, double length_c,  
     12                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    413{ 
    514    //return length_a * length_b * length_c; 
    6     return length_a * length_b * length_c + 
    7            2.0 * thick_rim_a * length_b * length_c + 
     15    return length_a * length_b * length_c +  
     16           2.0 * thick_rim_a * length_b * length_c +  
    817           2.0 * thick_rim_b * length_a * length_c + 
    918           2.0 * thick_rim_c * length_a * length_b; 
    1019} 
    1120 
    12 static double 
    13 Iq(double q, 
     21double Iq(double q, 
    1422    double core_sld, 
    1523    double arim_sld, 
    … …  
    2634    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    2735    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    28  
     36     
    2937    const double mu = 0.5 * q * length_b; 
    30  
     38     
    3139    //calculate volume before rescaling (in original code, but not used) 
    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 + 
     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 +  
    3543    //       2.0 * thick_rim_b * length_a * length_c + 
    3644    //       2.0 * thick_rim_c * length_a * length_b; 
    37  
     45     
    3846    // Scale sides by B 
    3947    const double a_scaled = length_a / length_b; 
    … …  
    93101            //   ( 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 
    94102            // This is the function called by csparallelepiped_analytical_2D_scaled, 
    95             // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 
    96  
     103            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
     104             
    97105            //  correct FF : sum of square of phase factors 
    98106            inner_total += Gauss76Wt[j] * form * form; 
    … …  
    110118} 
    111119 
    112 static double 
    113 Iqxy(double qa, double qb, double qc, 
     120double Iqxy(double qx, double qy, 
    114121    double core_sld, 
    115122    double arim_sld, 
    … …  
    122129    double thick_rim_a, 
    123130    double thick_rim_b, 
    124     double thick_rim_c) 
     131    double thick_rim_c, 
     132    double theta, 
     133    double phi, 
     134    double psi) 
    125135{ 
     136    double q, zhat, yhat, xhat; 
     137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     138 
    126139    // cspkernel in csparallelepiped recoded here 
    127140    const double dr0 = core_sld-solvent_sld; 
    … …  
    147160    double tc = length_a + 2.0*thick_rim_c; 
    148161    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    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  
     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     
    156169 
    157170    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    … …  
    161174               + drB*siA*(siBt-siB)*siC*V2 
    162175               + drC*siA*siB*(siCt*siCt-siC)*V3); 
    163  
     176    
    164177    return 1.0e-4 * f * f; 
    165178} 

  • sasmodels/models/cylinder.c

    rbecded3 r592343f  
     1double form_volume(double radius, double length); 
     2double fq(double q, double sn, double cn,double radius, double length); 
     3double orient_avg_1D(double q, double radius, double length); 
     4double Iq(double q, double sld, double solvent_sld, double radius, double length); 
     5double Iqxy(double qx, double qy, double sld, double solvent_sld, 
     6    double radius, double length, double theta, double phi); 
     7 
    18#define INVALID(v) (v.radius<0 || v.length<0) 
    29 
    3 static double 
    4 form_volume(double radius, double length) 
     10double form_volume(double radius, double length) 
    511{ 
    612    return M_PI*radius*radius*length; 
    713} 
    814 
    9 static double 
    10 fq(double qab, double qc, double radius, double length) 
     15double fq(double q, double sn, double cn, double radius, double length) 
    1116{ 
    12     return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
     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); 
    1321} 
    1422 
    15 static double 
    16 orient_avg_1D(double q, double radius, double length) 
     23double orient_avg_1D(double q, double radius, double length) 
    1724{ 
    1825    // translate a point in [-1,1] to a point in [0, pi/2] 
    1926    const double zm = M_PI_4; 
    20     const double zb = M_PI_4; 
     27    const double zb = M_PI_4;  
    2128 
    2229    double total = 0.0; 
    2330    for (int i=0; i<76 ;i++) { 
    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; 
     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; 
    3036    } 
    3137    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    3339} 
    3440 
    35 static double 
    36 Iq(double q, 
     41double Iq(double q, 
    3742    double sld, 
    3843    double solvent_sld, 
    … …  
    4449} 
    4550 
    46 static double 
    47 Iqxy(double qab, double qc, 
     51 
     52double Iqxy(double qx, double qy, 
    4853    double sld, 
    4954    double solvent_sld, 
    5055    double radius, 
    51     double length) 
     56    double length, 
     57    double theta, 
     58    double phi) 
    5259{ 
     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); 
    5363    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    54     const double form = fq(qab, qc, radius, length); 
     64    const double form = fq(q, sin_alpha, cos_alpha, radius, length); 
    5565    return 1.0e-4 * square(s * form); 
    5666} 

  • sasmodels/models/ellipsoid.c

    rbecded3 r3b571ae  
    1 static double 
    2 form_volume(double radius_polar, double radius_equatorial) 
     1double form_volume(double radius_polar, double radius_equatorial); 
     2double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial); 
     3double Iqxy(double qx, double qy, double sld, double sld_solvent, 
     4    double radius_polar, double radius_equatorial, double theta, double phi); 
     5 
     6double form_volume(double radius_polar, double radius_equatorial) 
    37{ 
    48    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    59} 
    610 
    7 static  double 
    8 Iq(double q, 
     11double Iq(double q, 
    912    double sld, 
    1013    double sld_solvent, 
    … …  
    3841} 
    3942 
    40 static double 
    41 Iqxy(double qab, double qc, 
     43double Iqxy(double qx, double qy, 
    4244    double sld, 
    4345    double sld_solvent, 
    4446    double radius_polar, 
    45     double radius_equatorial) 
     47    double radius_equatorial, 
     48    double theta, 
     49    double phi) 
    4650{ 
    47     const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc)); 
    48     const double f = sas_3j1x_x(qr); 
     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); 
    4956    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    5057 
    5158    return 1.0e-4 * square(f * s); 
    5259} 
     60 

  • sasmodels/models/elliptical_cylinder.c

    rbecded3 r61104c8  
    1 static double 
     1double form_volume(double radius_minor, double r_ratio, double length); 
     2double Iq(double q, double radius_minor, double r_ratio, double length, 
     3          double sld, double solvent_sld); 
     4double 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 
     8double 
    29form_volume(double radius_minor, double r_ratio, double length) 
    310{ 
    … …  
    512} 
    613 
    7 static double 
     14double 
    815Iq(double q, double radius_minor, double r_ratio, double length, 
    916   double sld, double solvent_sld) 
    … …  
    5461 
    5562 
    56 static double 
    57 Iqxy(double qa, double qb, double qc, 
     63double 
     64Iqxy(double qx, double qy, 
    5865     double radius_minor, double r_ratio, double length, 
    59      double sld, double solvent_sld) 
     66     double sld, double solvent_sld, 
     67     double theta, double phi, double psi) 
    6068{ 
     69    double q, xhat, yhat, zhat; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     71 
    6172    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    6273    // Given:    radius_major = r_ratio * radius_minor 
    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); 
     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); 
    6677    const double Aq = be * si; 
    6778    const double delrho = sld - solvent_sld; 

  • sasmodels/models/fcc_paracrystal.c

    rf728001 r50beefe  
    1 static double 
    2 fcc_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; 
     1double form_volume(double radius); 
     2double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
     3double 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); 
    86 
    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)); 
     7double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
     8double _FCCeval(double Theta, double Phi, double temp1, double temp3); 
    249 
    25     return Zq; 
     10 
     11double _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); 
    2619} 
    2720 
     21double _FCCeval(double Theta, double Phi, double temp1, double temp3) { 
    2822 
    29 // occupied volume fraction calculated from lattice symmetry and sphere radius 
    30 static double 
    31 fcc_volume_fraction(double radius, double dnn) 
    32 { 
    33     return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
     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); 
    3440} 
    3541 
    36 static double 
    37 form_volume(double radius) 
    38 { 
     42double form_volume(double radius){ 
    3943    return sphere_volume(radius); 
    4044} 
    4145 
    4246 
    43 static double Iq(double q, double dnn, 
     47double Iq(double q, double dnn, 
    4448  double d_factor, double radius, 
    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; 
     49  double sld, double solvent_sld){ 
    5350 
    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); 
     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); 
    7754 
    78     return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     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; 
     83 
     84 
    7985} 
    8086 
     87double 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); 
    8194 
    82 static 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; 
     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; 
    90112} 

  • sasmodels/models/hollow_cylinder.c

    rbecded3 r592343f  
     1double form_volume(double radius, double thickness, double length); 
     2double Iq(double q, double radius, double thickness, double length, double sld, 
     3        double solvent_sld); 
     4double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld, 
     5        double solvent_sld, double theta, double phi); 
     6 
    17//#define INVALID(v) (v.radius_core >= v.radius) 
    28 
    … …  
    814} 
    915 
     16 
    1017static double 
    11 _fq(double qab, double qc, 
    12     double radius, double thickness, double length) 
     18_hollow_cylinder_kernel(double q, 
     19    double radius, double thickness, double length, double sin_val, double cos_val) 
    1320{ 
    14     const double lam1 = sas_2J1x_x((radius+thickness)*qab); 
    15     const double lam2 = sas_2J1x_x(radius*qab); 
     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); 
    1624    const double gamma_sq = square(radius/(radius+thickness)); 
    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); 
     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); 
    2129    return psi*t2; 
    2230} 
    2331 
    24 static double 
     32double 
    2533form_volume(double radius, double thickness, double length) 
    2634{ 
    … …  
    3038 
    3139 
    32 static double 
     40double 
    3341Iq(double q, double radius, double thickness, double length, 
    3442    double sld, double solvent_sld) 
    3543{ 
    3644    const double lower = 0.0; 
    37     const double upper = 1.0;        //limits of numerical integral 
     45    const double upper = 1.0;           //limits of numerical integral 
    3846 
    39     double summ = 0.0;            //initialize intergral 
     47    double summ = 0.0;                  //initialize intergral 
    4048    for (int i=0;i<76;i++) { 
    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; 
     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; 
    4654    } 
    4755 
    … …  
    5159} 
    5260 
    53 static double 
    54 Iqxy(double qab, double qc, 
     61double 
     62Iqxy(double qx, double qy, 
    5563    double radius, double thickness, double length, 
    56     double sld, double solvent_sld) 
     64    double sld, double solvent_sld, double theta, double phi) 
    5765{ 
    58     const double form = _fq(qab, qc, radius, thickness, length); 
     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); 
    5970 
    6071    const double vol = form_volume(radius, thickness, length); 
    61     return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
     72    return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol); 
    6273} 
     74 

  • sasmodels/models/parallelepiped.c

    r9b7b23f rd605080  
    1 static double 
    2 form_volume(double length_a, double length_b, double length_c) 
     1double form_volume(double length_a, double length_b, double length_c); 
     2double Iq(double q, double sld, double solvent_sld, 
     3    double length_a, double length_b, double length_c); 
     4double 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 
     8double form_volume(double length_a, double length_b, double length_c) 
    39{ 
    410    return length_a * length_b * length_c; 
    … …  
    612 
    713 
    8 static double 
    9 Iq(double q, 
     14double Iq(double q, 
    1015    double sld, 
    1116    double solvent_sld, 
    … …  
    1520{ 
    1621    const double mu = 0.5 * q * length_b; 
    17  
     22     
    1823    // Scale sides by B 
    1924    const double a_scaled = length_a / length_b; 
    2025    const double c_scaled = length_c / length_b; 
    21  
     26         
    2227    // outer integral (with gauss points), integration limits = 0, 1 
    2328    double outer_total = 0; //initialize integral 
    … …  
    5257 
    5358 
    54 static double 
    55 Iqxy(double qa, double qb, double qc, 
     59double Iqxy(double qx, double qy, 
    5660    double sld, 
    5761    double solvent_sld, 
    5862    double length_a, 
    5963    double length_b, 
    60     double length_c) 
     64    double length_c, 
     65    double theta, 
     66    double phi, 
     67    double psi) 
    6168{ 
    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); 
     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); 
    6575    const double V = form_volume(length_a, length_b, length_c); 
    6676    const double drho = (sld - solvent_sld); 

  • sasmodels/models/sc_paracrystal.c

    rf728001 r50beefe  
    1 static double 
    2 sc_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; 
     1double form_volume(double radius); 
    82 
    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)); 
     3double Iq(double q, 
     4          double dnn, 
     5          double d_factor, 
     6          double radius, 
     7          double sphere_sld, 
     8          double solvent_sld); 
    249 
    25     return Zq; 
    26 } 
     10double 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); 
    2719 
    28 // occupied volume fraction calculated from lattice symmetry and sphere radius 
    29 static double 
    30 sc_volume_fraction(double radius, double dnn) 
    31 { 
    32     return sphere_volume(radius/dnn); 
    33 } 
    34  
    35 static double 
    36 form_volume(double radius) 
     20double form_volume(double radius) 
    3721{ 
    3822    return sphere_volume(radius); 
    3923} 
    4024 
     25static double 
     26sc_eval(double theta, double phi, double temp3, double temp4, double temp5) 
     27{ 
     28    double cnt, snt; 
     29    SINCOS(theta, cnt, snt); 
     30 
     31    double cnp, snp; 
     32    SINCOS(phi, cnp, snp); 
     33 
     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); 
     42} 
    4143 
    4244static double 
    43 Iq(double q, double dnn, 
    44     double d_factor, double radius, 
    45     double sld, double solvent_sld) 
     45sc_integrand(double dnn, double d_factor, double qq, double xx, double yy) 
    4646{ 
    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; 
     47    //Function to calculate integrand values for simple cubic structure 
    5348 
     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); 
    5455 
    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); 
     56        double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2; 
    7857 
    79     return sc_volume_fraction(radius, dnn) * Pq * Zq; 
     58        return(integrand); 
    8059} 
    8160 
     61double 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 
    8270 
    83 static double 
    84 Iqxy(double qa, double qb, double qc, 
    85     double dnn, double d_factor, double radius, 
    86     double sld, double solvent_sld) 
     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 
     103double 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) 
    87112{ 
    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; 
     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; 
    92128} 

  • 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

    rbecded3 r19f996b  
    1 static double 
    2 stacked_disks_kernel( 
    3     double qab, 
    4     double qc, 
     1static double stacked_disks_kernel( 
     2    double q, 
    53    double halfheight, 
    64    double thick_layer, 
    … …  
    119    double layer_sld, 
    1210    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 = radius*qab; 
    23     //const double besarg2 = radius*qab; 
     22    const double besarg1 = q*radius*sin_alpha; 
     23    //const double besarg2 = q*radius*sin_alpha; 
    2424 
    25     const double sinarg1 = halfheight*qc; 
    26     const double sinarg2 = (halfheight+thick_layer)*qc; 
     25    const double sinarg1 = q*halfheight*cos_alpha; 
     26    const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 
    2727 
    2828    const double be1 = sas_2J1x_x(besarg1); 
    … …  
    4343 
    4444    // loop for the structure factor S(q) 
    45     double qd_cos_alpha = d*qc; 
     45    double qd_cos_alpha = q*d*cos_alpha; 
    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 
    64 stacked_disks_1d( 
     63static double stacked_disks_1d( 
    6564    double q, 
    6665    double thick_core, 
    … …  
    8584        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8685        SINCOS(zi, sin_alpha, cos_alpha); 
    87         double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha, 
     86        double yyy = stacked_disks_kernel(q, 
    8887                           halfheight, 
    8988                           thick_layer, 
    … …  
    9493                           layer_sld, 
    9594                           solvent_sld, 
     95                           sin_alpha, 
     96                           cos_alpha, 
    9697                           d); 
    9798        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    … …  
    104105} 
    105106 
    106 static double 
    107 form_volume( 
     107static double form_volume( 
    108108    double thick_core, 
    109109    double thick_layer, 
    … …  
    116116} 
    117117 
    118 static double 
    119 Iq( 
     118static double Iq( 
    120119    double q, 
    121120    double thick_core, 
    … …  
    141140 
    142141 
    143 static double 
    144 Iqxy(double qab, double qc, 
     142static double Iqxy(double qx, double qy, 
    145143    double thick_core, 
    146144    double thick_layer, 
    … …  
    150148    double core_sld, 
    151149    double layer_sld, 
    152     double solvent_sld) 
     150    double solvent_sld, 
     151    double theta, 
     152    double phi) 
    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 
    155158    double d = 2.0 * thick_layer + thick_core; 
    156159    double halfheight = 0.5*thick_core; 
    157     double answer = stacked_disks_kernel(qab, qc, 
     160    double answer = stacked_disks_kernel(q, 
    158161                     halfheight, 
    159162                     thick_layer, 
    … …  
    164167                     layer_sld, 
    165168                     solvent_sld, 
     169                     sin_alpha, 
     170                     cos_alpha, 
    166171                     d); 
    167172 

  • sasmodels/models/triaxial_ellipsoid.c

    rbecded3 r68dd6a9  
     1double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar); 
     2double Iq(double q, double sld, double sld_solvent, 
     3    double radius_equat_minor, double radius_equat_major, double radius_polar); 
     4double 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 
    17//#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar) 
    28 
    3 static double 
    4 form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     9 
     10double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
    511{ 
    612    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
    713} 
    814 
    9 static double 
    10 Iq(double q, 
     15double Iq(double q, 
    1116    double sld, 
    1217    double sld_solvent, 
    … …  
    4045    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    4146    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
    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; 
     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; 
    4549} 
    4650 
    47 static double 
    48 Iqxy(double qa, double qb, double qc, 
     51double Iqxy(double qx, double qy, 
    4952    double sld, 
    5053    double sld_solvent, 
    5154    double radius_equat_minor, 
    5255    double radius_equat_major, 
    53     double radius_polar) 
     56    double radius_polar, 
     57    double theta, 
     58    double phi, 
     59    double psi) 
    5460{ 
    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); 
     61    double q, xhat, yhat, zhat; 
     62    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    6163 
    62     return 1.0e-4 * square(vol * drho * fq); 
     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); 
    6371} 
     72 

  • sasmodels/sasview_model.py

    r32f87a5 r9f8ade1  
    759759        if par.name not in self.params: 
    760760            if par.name == self.multiplicity_info.control: 
    761                 return self.multiplicity, [self.multiplicity], [1.0] 
     761                return [self.multiplicity], [1.0] 
    762762            else: 
    763763                # For hidden parameters use the default value. 
    764                 default = self._model_info.parameters.defaults.get(par.name, np.NaN) 
    765                 return [default], [1.0] 
     764                value = self._model_info.parameters.defaults.get(par.name, np.NaN) 
     765                return [value], [1.0] 
    766766        elif par.polydisperse: 
    767             value = self.params[par.name] 
    768767            dis = self.dispersion[par.name] 
    769768            if dis['type'] == 'array': 
    770                 dispersity, weight = dis['values'], dis['weights'] 
     769                value, weight = dis['values'], dis['weights'] 
    771770            else: 
    772                 dispersity, weight = weights.get_weights( 
     771                value, weight = weights.get_weights( 
    773772                    dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
    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] 
     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] 
    779777 
    780778def test_cylinder(): 

  • sasmodels/weights.py

    r3c24ccd rf1a8811  
    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 
    6157        if sigma == 0 or self.npts < 2: 
    6258            if lb <= center <= ub: 
    … …  
    229225    obj = cls(n, width, nsigmas) 
    230226    v, w = obj.get_weights(value, limits[0], limits[1], relative) 
    231     return v, w/np.sum(w) 
    232  
    233  
    234 def plot_weights(model_info, mesh): 
    235     # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None 
     227    return v, w 
     228 
     229 
     230def plot_weights(model_info, pairs): 
     231    # type: (ModelInfo, List[Tuple[np.ndarray, np.ndarray]]) -> None 
    236232    """ 
    237233    Plot the weights returned by :func:`get_weights`. 
    238234 
    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. 
     235    *model_info* is 
     236    :param model_info: 
     237    :param pairs: 
     238    :return: 
    244239    """ 
    245240    import pylab 
    246241 
    247     if any(len(dispersity)>1 for value, dispersity, weights in mesh): 
     242    if any(len(values)>1 for values, weights in pairs): 
    248243        labels = [p.name for p in model_info.parameters.call_parameters] 
    249         #pylab.interactive(True) 
     244        pylab.interactive(True) 
    250245        pylab.figure() 
    251         for (v,x,w), s in zip(mesh, labels): 
    252             if len(x) > 1: 
    253                 pylab.plot(x, w, '-o', label=s) 
     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) 
    254250        pylab.grid(True) 
    255251        pylab.legend() 
Note: See TracChangeset for help on using the changeset viewer.