Changeset 62d7601 in sasmodels


Ignore:
Timestamp:
Sep 28, 2017 4:24:09 PM (5 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
706f466
Parents:
bf09f55 (diff), 505d0ad (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-776-orientation

Files:
9 added
35 edited

Legend:

Unmodified
Added
Removed
  • doc/genmodel.py

    r745b7bb rb866abf  
    33import sys, os, math, re 
    44import numpy as np 
     5import matplotlib 
     6matplotlib.use('Agg') 
    57import matplotlib.pyplot as plt 
    68sys.path.insert(0, os.path.abspath('..')) 
  • doc/guide/index.rst

    r2e66ef5 rc0d7ab3  
    99   intro.rst 
    1010   install.rst 
     11   gpu_setup.rst 
    1112   pd/polydispersity.rst 
    1213   resolution.rst 
  • doc/guide/install.rst

    rf8a2baa rc0d7ab3  
    5353This will allow you to edit the files in the package and have the changes 
    5454show up immediately in python the next time you load your program. 
    55  
    56 OpenCL Installation 
    57 ******************* 
    58 *Warning! GPU devices do not in general offer the same level of memory 
    59 protection as CPU devices. If your code attempts to write outside allocated 
    60 memory buffers unpredicatable behaviour may result (eg, your video display 
    61 may freeze, or your system may crash, etc). Do not install OpenCL drivers 
    62 without first checking for known issues (eg, some computer manufacturers 
    63 install modified graphics drivers so replacing these may not be a good 
    64 idea!). If in doubt, seek advice from an IT professional before proceeding 
    65 further.* 
    66  
    67 Check if you have OpenCL already installed 
    68 ========================================== 
    69  
    70 **Windows** 
    71  
    72 The following instructions are based on 
    73 http://web.engr.oregonstate.edu/~mjb/cs475/DoIHaveOpenCL.pdf 
    74  
    75 * Go to: Start -> Control Panel -> System & Security -> Administrative Tools 
    76 * Double Click on Computer Managment 
    77 * Click on Device Manager 
    78 * Click open Display Adapters 
    79 * Right-click on available adapter and select Properties 
    80 * Click on Driver 
    81 * Go to Driver Details 
    82 * Scroll down and see if OpenCL is installed (look for OpenCL*.dll files) 
    83  
    84 **Mac OSX** 
    85  
    86 For OS X operating systems higher than 10.6 OpenCL is shipped along with 
    87 the system. 
    88  
    89 However, OpenCL has had a rocky history on Macs. Apple provide a useful 
    90 compatibility table at https://support.apple.com/en-us/HT202823 
    91  
    92  
    93 Installation 
    94 ============ 
    95  
    96 **Windows** 
    97  
    98 Depending on the graphic card in your system, drivers 
    99 can be obtained from different sources: 
    100  
    101 * NVIDIA: https://developer.nvidia.com/opencl 
    102 * AMD: http://developer.amd.com/tools-and-sdks/opencl-zone/ 
    103  
    104  
    105 **Mac OSX** 
    106  
    107 N/A 
    108  
    109 You cannot download OpenCL driver updates for your Mac. They are packaged 
    110 with the normal quarterly OS X updates from Apple. 
    111  
    112  
    113 .. note:: 
    114     Intel provides OpenCL drivers for Intel processors at 
    115     https://software.intel.com/en-us/articles/opencl-drivers 
    116     These can sometimes make use of special vector instructions across multiple 
    117     processors, so it is worth installing if the GPU does not support double 
    118     precision. You can install this driver alongside the GPU driver for NVIDIA 
    119     or AMD. 
    120  
    121  
    122 GPU Selection 
    123 ************* 
    124  
    125 sasmodels evaluations can run on your graphics card (GPU) or they can run 
    126 on the processor (CPU). In general, calculations performed on the GPU will run faster. 
    127  
    128 To run on the GPU, your computer must have OpenCL drivers installed. 
    129 For information about OpenCL installation see this 
    130 :ref:`opencl-installation` guidance. 
    131  
    132 Where the model is evaluated is a little bit complicated. 
    133 If the model has the line *single=False* then it requires double precision. 
    134 If the GPU is single precision only, then it will try running via OpenCL 
    135 on the CPU.  If the OpenCL driver is not available for the CPU then 
    136 it will run as a normal program on the CPU. 
    137  
    138 For models with a large number of parameters or with a lot of code, 
    139 the GPU may be too small to run the program effectively. 
    140 In this case, you should try simplifying the model, maybe breaking it 
    141 into several different modules so that you don't need *IF* statements in your code. 
    142 If it is still too big, you can set *opencl=False* in the model file and 
    143 the model will only run as a normal program on the CPU. 
    144 This will not usually be necessary. 
    145  
    146 Device Selection 
    147 ================ 
    148 If you have multiple GPU devices you can tell SasView which device to use. 
    149 By default, SasView looks for one GPU and one CPU device 
    150 from available OpenCL platforms. 
    151  
    152 SasView prefers AMD or NVIDIA drivers for GPU, and prefers Intel or 
    153 Apple drivers for CPU. Both GPU and CPU are included on the assumption that CPU 
    154 is always available and supports double precision. 
    155  
    156 The device order is important: GPU is checked before CPU on the assumption that 
    157 it will be faster. By examining ~/sasview.log you can see which device SasView 
    158 chose to run the model. 
    159  
    160 **If you don't want to use OpenCL, you can set** *SAS_OPENCL=None* 
    161 **in your environment settings, and it will only use normal programs.** 
    162  
    163 If you want to use one of the other devices, you can run the following 
    164 from the python console in SasView:: 
    165  
    166     import pyopencl as cl 
    167     cl.create_some_context() 
    168  
    169 This will provide a menu of different OpenCL drivers available. 
    170 When one is selected, it will say "set PYOPENCL_CTX=..." 
    171 Use that value as the value of *SAS_OPENCL*. 
    172  
    173 Compiler Selection 
    174 ================== 
    175 For models run as normal programs, you may need to specify a compiler. 
    176 This is done using the SAS_COMPILER environment variable. 
    177 Set it to *tinycc* for the tinycc compiler, *msvc* for the 
    178 Microsoft Visual C compiler, or *mingw* for the MinGW compiler. 
    179 TinyCC is provided with SasView so that is the default. 
    180 If you want one of the other compilers, be sure to have it available 
    181 in your *PATH* so SasView can find it! 
    182  
    183  
    184 *Document History* 
    185  
    186 | 2017-05-17 Paul Kienzle 
  • doc/guide/magnetism/magnetism.rst

    r64eecf7 r59485a4  
    8787    The values of the 'Up_frac_i' and 'Up_frac_f' must be in the range 0 to 1. 
    8888 
    89 .. note:: 
    90     This help document was last changed by Steve King, 02May2015 
     89*Document History* 
    9190 
    92 * Document History * 
    93  
     91| 2015-05-02 Steve King 
    9492| 2017-05-08 Paul Kienzle 
  • doc/guide/plugin.rst

    r30b60d2 r3048ec6  
    7878     at the start of the model documentation and as a tooltip in the SasView GUI 
    7979 
    80 - a **short discription**: 
     80- a **short description**: 
    8181   - this is the **description** string in the *.py* file 
    8282   - this is a medium length description which appears when you click 
     
    233233 
    234234**name = "mymodel"** defines the name of the model that is shown to the user. 
    235 If it is not provided, it will use the name of the model file, with '_' 
    236 replaced by spaces and the parts capitalized.  So *adsorbed_layer.py* will 
    237 become *Adsorbed Layer*.  The predefined models all use the name of the 
    238 model file as the name of the model, so the default may be changed. 
     235If it is not provided it will use the name of the model file. The name must 
     236be a valid variable name, starting with a letter and contains only letters, 
     237numbers or underscore.  Spaces, dashes, and other symbols are not permitted. 
    239238 
    240239**title = "short description"** is short description of the model which 
     
    298297- **"name"** is the name of the parameter shown on the FitPage. 
    299298 
     299  - the name must be a valid variable name, starting with a letter and 
     300    containing only letters, numbers and underscore. 
     301 
    300302  - parameter names should follow the mathematical convention; e.g., 
    301303    *radius_core* not *core_radius*, or *sld_solvent* not *solvent_sld*. 
     
    366368    dispersion. 
    367369 
     370Some models will have integer parameters, such as number of pearls in the 
     371pearl necklace model, or number of shells in the multi-layer vesicle model. 
     372The optimizers in BUMPS treat all parameters as floating point numbers which 
     373can take arbitrary values, even for integer parameters, so your model should 
     374round the incoming parameter value to the nearest integer inside your model 
     375you should round to the nearest integer.  In C code, you can do this using:: 
     376 
     377    static double 
     378    Iq(double q, ..., double fp_n, ...) 
     379    { 
     380        int n = (int)(fp_n + 0.5); 
     381        ... 
     382    } 
     383 
     384in python:: 
     385 
     386    def Iq(q, ..., n, ...): 
     387        n = int(n+0.5) 
     388        ... 
     389 
     390Derivative based optimizers such as Levenberg-Marquardt will not work 
     391for integer parameters since the partial derivative is always zero, but 
     392the remaining optimizers (DREAM, differential evolution, Nelder-Mead simplex) 
     393will still function. 
    368394 
    369395Model Computation 
     
    391417     \approx \frac{\sum_{i=1}^n G(x_i) I(q,x_i)}{\sum_{i=1}^n G(x_i) V(x_i)} 
    392418 
    393 That is, the indivdual models do not need to include polydispersity 
     419That is, the individual models do not need to include polydispersity 
    394420calculations, but instead rely on numerical integration to compute the 
    395421appropriately smeared pattern.   Angular dispersion values over polar angle 
     
    411437The parameters *par1, par2, ...* are the list of non-orientation parameters 
    412438to the model in the order that they appear in the parameter table. 
    413 **Note that the autogenerated model file uses** *x* **rather than** *q*. 
     439**Note that the auto-generated model file uses** *x* **rather than** *q*. 
    414440 
    415441The *.py* file should import trigonometric and exponential functions from 
     
    935961across multiple parameters can be very slow). 
    936962 
    937 If your model has 2D orientational calculation, then you should also 
     963If your model has 2D orientation calculation, then you should also 
    938964test with:: 
    939965 
  • sasmodels/conversion_table.py

    r0c2da4b r505d0ad  
    1111account for that. 
    1212 
    13 Usage: 
    14 <old_Sasview_version> : { 
    15     <new_model_name> : [ 
    16         <old_model_name> , 
    17         { 
    18             <new_param_name_1> : <old_param_name_1>, 
    19             ... 
    20             <new_param_name_n> : <old_param_name_n> 
    21         } 
    22     ] 
    23 } 
     13Usage:: 
     14 
     15    <old_Sasview_version> : { 
     16        <new_model_name> : [ 
     17            <old_model_name> , 
     18            { 
     19                <new_param_name_1> : <old_param_name_1>, 
     20                ... 
     21                <new_param_name_n> : <old_param_name_n> 
     22            } 
     23        ] 
     24    } 
    2425 
    2526Any future parameter and model name changes can and should be given in this 
  • sasmodels/custom/__init__.py

    r3c852d4 r0f48f1e  
    3131        if fullname in sys.modules: 
    3232            del sys.modules[fullname] 
     33        if path.endswith(".py") and os.path.exists(path) and os.path.exists(path+"c"): 
     34            # remove automatic pyc file before loading a py file 
     35            os.unlink(path+"c") 
    3336        module = imp.load_source(fullname, os.path.expanduser(path)) 
    34         #os.unlink(path+"c")  # remove the automatic pyc file 
    3537        return module 
    3638 
  • sasmodels/kernelpy.py

    r9bc4882 r62d7601  
    3434        self.info = model_info 
    3535        self.dtype = np.dtype('d') 
     36        logging.info("load python model " + self.info.name) 
    3637 
    3738    def make_kernel(self, q_vectors): 
    38         logging.info("creating python kernel " + self.info.name) 
    3939        q_input = PyInput(q_vectors, dtype=F64) 
    4040        kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 
  • sasmodels/sasview_model.py

    r9644b5a r9f8ade1  
    1616import logging 
    1717from os.path import basename, splitext, abspath, getmtime 
    18 import thread 
     18try: 
     19    import _thread as thread 
     20except ImportError: 
     21    import thread 
    1922 
    2023import numpy as np  # type: ignore 
  • explore/jitter.py

    r85190c2 rd4c33d6  
    33dispersity and possible replacement algorithms. 
    44""" 
     5import sys 
     6 
    57import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
    68import matplotlib.pyplot as plt 
    79from matplotlib.widgets import Slider, CheckButtons 
    810from matplotlib import cm 
    9  
    1011import numpy as np 
    1112from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1213 
    13 def draw_beam(ax): 
     14def draw_beam(ax, view=(0, 0)): 
    1415    #ax.plot([0,0],[0,0],[1,-1]) 
    1516    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
     
    2223    x = r*np.outer(np.cos(u), np.ones_like(v)) 
    2324    y = r*np.outer(np.sin(u), np.ones_like(v)) 
    24     z = np.outer(np.ones_like(u), v) 
     25    z = 1.3*np.outer(np.ones_like(u), v) 
     26 
     27    theta, phi = view 
     28    shape = x.shape 
     29    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
     30    points = Rz(phi)*Ry(theta)*points 
     31    x, y, z = [v.reshape(shape) for v in points] 
    2532 
    2633    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    2734 
    28 def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi): 
    29     size=[0.1, 0.4, 1.0] 
    30     view=[theta, phi, psi] 
    31     shimmy=[0,0,0] 
    32     #draw_shape = draw_parallelepiped 
    33     draw_shape = draw_ellipsoid 
     35def draw_jitter(ax, view, jitter): 
     36    size = [0.1, 0.4, 1.0] 
     37    draw_shape = draw_parallelepiped 
     38    #draw_shape = draw_ellipsoid 
    3439 
    3540    #np.random.seed(10) 
     
    6469        [ 1,  1,  1], 
    6570    ] 
     71    dtheta, dphi, dpsi = jitter 
    6672    if dtheta == 0: 
    6773        cloud = [v for v in cloud if v[0] == 0] 
     
    7076    if dpsi == 0: 
    7177        cloud = [v for v in cloud if v[2] == 0] 
    72     draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8) 
     78    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    7379    for point in cloud: 
    74         shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]] 
    75         draw_shape(ax, size, view, shimmy, alpha=0.8) 
     80        delta = [dtheta*point[0], dphi*point[1], dpsi*point[2]] 
     81        draw_shape(ax, size, view, delta, alpha=0.8) 
    7682    for v in 'xyz': 
    7783        a, b, c = size 
     
    8086        getattr(ax, v+'axis').label.set_text(v) 
    8187 
    82 def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1): 
     88def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
    8389    a,b,c = size 
    84     theta, phi, psi = view 
    85     dtheta, dphi, dpsi = shimmy 
    86  
    8790    u = np.linspace(0, 2 * np.pi, steps) 
    8891    v = np.linspace(0, np.pi, steps) 
     
    9093    y = b*np.outer(np.sin(u), np.sin(v)) 
    9194    z = c*np.outer(np.ones_like(u), np.cos(v)) 
    92  
    93     shape = x.shape 
    94     points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
    95     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    96     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    97     x,y,z = [v.reshape(shape) for v in points] 
     95    x, y, z = transform_xyz(view, jitter, x, y, z) 
    9896 
    9997    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    10098 
    101 def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 
     99    draw_labels(ax, view, jitter, [ 
     100         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     101         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     102         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     103         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     104         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     105         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     106    ]) 
     107 
     108def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
    102109    a,b,c = size 
    103     theta, phi, psi = view 
    104     dtheta, dphi, dpsi = shimmy 
    105  
    106110    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    107111    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
     
    118122    ]) 
    119123 
    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] 
     124    x, y, z = transform_xyz(view, jitter, x, y, z) 
    125125    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    126126 
    127 def draw_sphere(ax, radius=10., steps=100): 
    128     u = np.linspace(0, 2 * np.pi, steps) 
    129     v = np.linspace(0, np.pi, steps) 
    130  
    131     x = radius * np.outer(np.cos(u), np.sin(v)) 
    132     y = radius * np.outer(np.sin(u), np.sin(v)) 
    133     z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
    134     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    135  
    136 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    137     theta_center = radians(theta) 
    138     phi_center = radians(phi) 
    139     flow_center = radians(flow) 
    140     dtheta = radians(dtheta) 
    141     dphi = radians(dphi) 
    142  
    143     # 10 point 3-sigma gaussian weights 
    144     t = np.linspace(-3., 3., 11) 
     127    draw_labels(ax, view, jitter, [ 
     128         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     129         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     130         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     131         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     132         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     133         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     134    ]) 
     135 
     136def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gauss'): 
     137    theta, phi, psi = view 
     138    dtheta, dphi, dpsi = jitter 
    145139    if dist == 'gauss': 
     140        t = np.linspace(-3, 3, n) 
    146141        weights = exp(-0.5*t**2) 
    147142    elif dist == 'rect': 
     143        t = np.linspace(0, 1, n) 
    148144        weights = np.ones_like(t) 
    149145    else: 
    150146        raise ValueError("expected dist to be 'gauss' or 'rect'") 
    151     theta = dtheta*t 
    152     phi = dphi*t 
    153  
    154     x = radius * np.outer(cos(phi), cos(theta)) 
    155     y = radius * np.outer(sin(phi), cos(theta)) 
    156     z = radius * np.outer(np.ones_like(phi), sin(theta)) 
    157     #w = np.outer(weights, weights*abs(cos(dtheta*t))) 
    158     w = np.outer(weights, weights*abs(cos(theta))) 
    159  
    160     x, y, z, w = [v.flatten() for v in (x,y,z,w)] 
    161     x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center) 
    162  
    163     ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
    164  
    165 def rotate(x, y, z, phi, theta, psi): 
    166     R = Rz(phi)*Ry(theta)*Rz(psi) 
    167     p = np.vstack([x,y,z]) 
    168     return R*p 
     147 
     148    # mesh in theta, phi formed by rotating z 
     149    z = np.matrix([[0], [0], [radius]]) 
     150    points = np.hstack([Rx(phi_i)*Ry(theta_i)*z 
     151                        for theta_i in dtheta*t 
     152                        for phi_i in dphi*t]) 
     153    # rotate relative to beam 
     154    points = orient_relative_to_beam(view, points) 
     155 
     156    w = np.outer(weights, weights) 
     157 
     158    x, y, z = [np.array(v).flatten() for v in points] 
     159    ax.scatter(x, y, z, c=w.flatten(), marker='o', vmin=0., vmax=1.) 
    169160 
    170161def Rx(angle): 
     
    188179         [0., 0., 1.]] 
    189180    return np.matrix(R) 
     181 
     182def transform_xyz(view, jitter, x, y, z): 
     183    x, y, z = [np.asarray(v) for v in (x, y, z)] 
     184    shape = x.shape 
     185    points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     186    points = apply_jitter(jitter, points) 
     187    points = orient_relative_to_beam(view, points) 
     188    x, y, z = [np.array(v).reshape(shape) for v in points] 
     189    return x, y, z 
     190 
     191def apply_jitter(jitter, points): 
     192    dtheta, dphi, dpsi = jitter 
     193    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     194    return points 
     195 
     196def orient_relative_to_beam(view, points): 
     197    theta, phi, psi = view 
     198    points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     199    return points 
     200 
     201def draw_labels(ax, view, jitter, text): 
     202    labels, locations, orientations = zip(*text) 
     203    px, py, pz = zip(*locations) 
     204    dx, dy, dz = zip(*orientations) 
     205 
     206    px, py, pz = transform_xyz(view, jitter, px, py, pz) 
     207    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz) 
     208 
     209    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 
     210        zdir = np.asarray(zdir).flatten() 
     211        ax.text(p[0], p[1], p[2], label, zdir=zdir) 
     212 
     213def draw_sphere(ax, radius=10., steps=100): 
     214    u = np.linspace(0, 2 * np.pi, steps) 
     215    v = np.linspace(0, np.pi, steps) 
     216 
     217    x = radius * np.outer(np.cos(u), np.sin(v)) 
     218    y = radius * np.outer(np.sin(u), np.sin(v)) 
     219    z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
     220    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    190221 
    191222def main(): 
     
    206237 
    207238    axcolor = 'lightgoldenrodyellow' 
     239 
    208240    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    209241    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
     
    212244    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    213245    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
     246 
    214247    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    215248    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
     
    217250    sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta) 
    218251    sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi) 
    219     sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi) 
     252    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dpsi) 
    220253 
    221254    def update(val, axis=None): 
    222         theta, phi, psi = stheta.val, sphi.val, spsi.val 
    223         dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val 
     255        view = stheta.val, sphi.val, spsi.val 
     256        jitter = sdtheta.val, sdphi.val, sdpsi.val 
    224257        ax.cla() 
    225         draw_beam(ax) 
    226         draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi) 
    227         #if not axis.startswith('d'): 
    228         #    ax.view_init(elev=theta, azim=phi) 
     258        draw_beam(ax, (0, 0)) 
     259        draw_jitter(ax, view, jitter) 
     260        #draw_jitter(ax, view, (0,0,0)) 
     261        draw_mesh(ax, view, jitter) 
    229262        plt.gcf().canvas.draw() 
    230263 
  • sasmodels/details.py

    rccd5f01 rf39759c  
    1515 
    1616import numpy as np  # type: ignore 
    17 from numpy import pi, cos, sin 
     17from numpy import pi, cos, sin, radians 
    1818 
    1919try: 
     
    2929 
    3030try: 
    31     from typing import List 
     31    from typing import List, Tuple, Sequence 
    3232except ImportError: 
    3333    pass 
    3434else: 
    3535    from .modelinfo import ModelInfo 
     36    from .modelinfo import ParameterTable 
    3637 
    3738 
     
    5354    coordinates, the total circumference decreases as latitude varies from 
    5455    pi r^2 at the equator to 0 at the pole, and the weight associated 
    55     with a range of phi values needs to be scaled by this circumference. 
     56    with a range of latitude values needs to be scaled by this circumference. 
    5657    This scale factor needs to be updated each time the theta value 
    5758    changes.  *theta_par* indicates which of the values in the parameter 
     
    231232    nvalues = kernel.info.parameters.nvalues 
    232233    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 ((),()) 
     234    # skipping scale and background when building values and weights 
     235    values, weights = zip(*pairs[2:npars+2]) if npars else ((), ()) 
     236    weights = correct_theta_weights(kernel.info.parameters, values, weights) 
    234237    length = np.array([len(w) for w in weights]) 
    235238    offset = np.cumsum(np.hstack((0, length))) 
     
    244247    return call_details, data, is_magnetic 
    245248 
     249def correct_theta_weights(parameters, values, 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.  Avoid 
     254    evaluation strictly at the pole, which would otherwise send the weight to 
     255    zero. 
     256 
     257    Note: values and weights do not include scale and background 
     258    """ 
     259    # TODO: document code, explaining why scale and background are skipped 
     260    # given that we don't have scale and background in the list, we 
     261    # should be calling the variables something other than values and weights 
     262    # Apparently the parameters.theta_offset similarly skips scale and 
     263    # and background, so the indexing works out. 
     264    if parameters.theta_offset >= 0: 
     265        index = parameters.theta_offset 
     266        theta = values[index] 
     267        theta_weight = np.minimum(abs(cos(radians(theta))), 1e-6) 
     268        # copy the weights list so we can update it 
     269        weights = list(weights) 
     270        weights[index] = theta_weight*np.asarray(weights[index]) 
     271        weights = tuple(weights) 
     272    return weights 
     273 
    246274 
    247275def convert_magnetism(parameters, values): 
     276    # type: (ParameterTable, Sequence[np.ndarray]) -> bool 
    248277    """ 
    249278    Convert magnetism values from polar to rectangular coordinates. 
     
    255284    scale = mag[:,0] 
    256285    if np.any(scale): 
    257         theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. 
     286        theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
    258287        cos_theta = cos(theta) 
    259288        mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
     
    269298    """ 
    270299    Create a mesh grid of dispersion parameters and weights. 
     300 
     301    pars is a list of pairs (values, weights), where the values are the 
     302    individual parameter values at which to evaluate the polydispersity 
     303    distribution and weights are the weights associated with each value. 
     304 
     305    Only the volume parameters should be included in this list.  Orientation 
     306    parameters do not affect the calculation of effective radius or volume 
     307    ratio. 
    271308 
    272309    Returns [p1,p2,...],w where pj is a vector of values for parameter j 
  • sasmodels/kernel_header.c

    rbb4b509 r73cbc5b  
    164164    SINCOS(phi*M_PI_180, sn, cn); \ 
    165165    q = sqrt(qx*qx + qy*qy); \ 
    166     cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
     166    cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
    167167    sn = sqrt(1 - cn*cn); \ 
    168168    } while (0) 
  • sasmodels/kernel_iq.c

    rbde38b5 rd4c33d6  
    2525    int32_t num_weights;        // total length of the weights vector 
    2626    int32_t num_active;         // number of non-trivial pd loops 
    27     int32_t theta_par;          // id of spherical correction variable 
     27    int32_t theta_par;          // id of spherical correction variable (not used) 
    2828} ProblemDetails; 
    2929 
     
    173173 
    174174 
    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  
    186175  int step = pd_start; 
    187176 
     
    220209#endif 
    221210#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   } 
    225211  while(i0 < n0) { 
    226212    local_values.vector[p0] = v0[i0]; 
    227213    double weight0 = w0[i0] * weight1; 
    228214//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     } 
    232215#else 
    233216    const double weight0 = 1.0; 
     
    244227      // Note: weight==0 must always be excluded 
    245228      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); 
     229        pd_norm += weight0 * CALL_VOLUME(local_values.table); 
    250230 
    251231        #ifdef USE_OPENMP 
     
    304284#endif // !MAGNETIC 
    305285//printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0); 
    306           result[q_index] += weight * scattering; 
     286          result[q_index] += weight0 * scattering; 
    307287        } 
    308288      } 
  • sasmodels/kernel_iq.cl

    rbde38b5 rd4c33d6  
    2525    int32_t num_weights;        // total length of the weights vector 
    2626    int32_t num_active;         // number of non-trivial pd loops 
    27     int32_t theta_par;          // id of spherical correction variable 
     27    int32_t theta_par;          // id of spherical correction variable (not used) 
    2828} ProblemDetails; 
    2929 
     
    169169 
    170170 
    171 #if MAX_PD>0 
    172   const int theta_par = details->theta_par; 
    173   const bool fast_theta = (theta_par == p0); 
    174   const bool slow_theta = (theta_par >= 0 && !fast_theta); 
    175   double spherical_correction = 1.0; 
    176 #else 
    177   // Note: if not polydisperse the weights cancel and we don't need the 
    178   // spherical correction. 
    179   const double spherical_correction = 1.0; 
    180 #endif 
    181  
    182171  int step = pd_start; 
    183172 
     
    217206#endif 
    218207#if MAX_PD>0 
    219   if (slow_theta) { // Theta is not in inner loop 
    220     spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6); 
    221   } 
    222208  while(i0 < n0) { 
    223209    local_values.vector[p0] = v0[i0]; 
    224210    double weight0 = w0[i0] * weight1; 
    225211//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0); 
    226     if (fast_theta) { // Theta is in inner loop 
    227       spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6); 
    228     } 
    229212#else 
    230213    const double weight0 = 1.0; 
     
    232215 
    233216//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"); } 
    234 //if (q_index == 0) printf("sphcor: %g\n", spherical_correction); 
    235217 
    236218    #ifdef INVALID 
     
    241223      // Note: weight==0 must always be excluded 
    242224      if (weight0 > cutoff) { 
    243         // spherical correction is set at a minimum of 1e-6, otherwise there 
    244         // would be problems looking at models with theta=90. 
    245         const double weight = weight0 * spherical_correction; 
    246         pd_norm += weight * CALL_VOLUME(local_values.table); 
     225        pd_norm += weight0 * CALL_VOLUME(local_values.table); 
    247226 
    248227#if defined(MAGNETIC) && NUM_MAGNETIC > 0 
     
    296275        const double scattering = CALL_IQ(q, q_index, local_values.table); 
    297276#endif // !MAGNETIC 
    298         this_result += weight * scattering; 
     277        this_result += weight0 * scattering; 
    299278      } 
    300279    } 
  • sasmodels/model_test.py

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

    r592343f r2a0b2b1  
    1010//barbell kernel - same as dumbell 
    1111static double 
    12 _bell_kernel(double q, double h, double radius_bell, 
    13              double half_length, double sin_alpha, double cos_alpha) 
     12_bell_kernel(double qab, double qc, double h, double radius_bell, 
     13             double half_length) 
    1414{ 
    1515    // translate a point in [-1,1] to a point in [lower,upper] 
     
    2626    //    m = q R cos(alpha) 
    2727    //    b = q(L/2-h) cos(alpha) 
    28     const double m = q*radius_bell*cos_alpha; // cos argument slope 
    29     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    30     const double qrst = q*radius_bell*sin_alpha; // Q*R*sin(theta) 
     28    const double m = radius_bell*qc; // cos argument slope 
     29    const double b = (half_length-h)*qc; // cos argument intercept 
     30    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    3232    for (int i = 0; i < 76; i++){ 
    3333        const double t = Gauss76Z[i]*zm + zb; 
    3434        const double radical = 1.0 - t*t; 
    35         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     35        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    3636        const double Fq = cos(m*t + b) * radical * bj; 
    3737        total += Gauss76Wt[i] * Fq; 
     
    4444 
    4545static double 
    46 _fq(double q, double h, 
    47     double radius_bell, double radius, double half_length, 
    48     double sin_alpha, double cos_alpha) 
     46_fq(double qab, double qc, double h, 
     47    double radius_bell, double radius, double half_length) 
    4948{ 
    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); 
     49    const double bell_fq = _bell_kernel(qab, qc, h, radius_bell, half_length); 
     50    const double bj = sas_2J1x_x(radius*qab); 
     51    const double si = sas_sinx_x(half_length*qc); 
    5352    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5453    const double Aq = bell_fq + cyl_fq; 
     
    8483        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8584        SINCOS(alpha, sin_alpha, cos_alpha); 
    86         const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 
     85        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    8786        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8887    } 
     
    103102    double q, sin_alpha, cos_alpha; 
    104103    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     104    const double qab = q*sin_alpha; 
     105    const double qc = q*cos_alpha; 
    105106 
    106107    const double h = -sqrt(square(radius_bell) - square(radius)); 
    107     const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 
     108    const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length); 
    108109 
    109110    // Multiply by contrast^2 and convert to cm-1 
  • sasmodels/models/bcc_paracrystal.c

    r50beefe r7e0b281  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4#if 0  // Equations as written in Matsuoka 
     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; 
     8#else 
     9    const double a1 = (+qa + qb - qc)/2.0; 
     10    const double a2 = (+qa - qb + qc)/2.0; 
     11    const double a3 = (-qa + qb + qc)/2.0; 
     12#endif 
    613 
    7 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _BCCeval(double Theta, double Phi, double temp1, double temp3); 
    9 double _sphereform(double q, double radius, double sld, double solvent_sld); 
     14#if 1 
     15    // Numerator: (1 - exp(a)^2)^3 
     16    //         => (-(exp(2a) - 1))^3 
     17    //         => -expm1(2a)^3 
     18    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 
     19    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     20    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 
     21    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     22    const double exp_arg = exp(arg); 
     23    const double Zq = -cube(expm1(2.0*arg)) 
     24        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     25          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     26          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     27#else 
     28    // Alternate form, which perhaps is more approachable 
     29    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     30    const double sinh_qd = sinh(arg); 
     31    const double cosh_qd = cosh(arg); 
     32    const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1)) 
     33                    * sinh_qd/(cosh_qd - cos(dnn*a2)) 
     34                    * sinh_qd/(cosh_qd - cos(dnn*a3)); 
     35#endif 
     36 
     37    return Zq; 
     38} 
    1039 
    1140 
    12 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    13  
    14         const double Da = d_factor*dnn; 
    15         const double temp1 = q*q*Da*Da; 
    16         const double temp3 = q*dnn; 
    17  
    18         double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    19         return(retVal); 
     41// occupied volume fraction calculated from lattice symmetry and sphere radius 
     42static double 
     43bcc_volume_fraction(double radius, double dnn) 
     44{ 
     45    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
    2046} 
    2147 
    22 double _BCCeval(double Theta, double Phi, double temp1, double temp3) { 
    23  
    24         double result; 
    25         double sin_theta,cos_theta,sin_phi,cos_phi; 
    26         SINCOS(Theta, sin_theta, cos_theta); 
    27         SINCOS(Phi, sin_phi, cos_phi); 
    28  
    29         const double temp6 =  sin_theta; 
    30         const double temp7 =  sin_theta*cos_phi + sin_theta*sin_phi + cos_theta; 
    31         const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 
    32         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 
    33  
    34         const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 
    35         result = cube(1.0 - (temp10*temp10))*temp6 
    36             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    38               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    39  
    40         return (result); 
    41 } 
    42  
    43 double form_volume(double radius){ 
     48static double 
     49form_volume(double radius) 
     50{ 
    4451    return sphere_volume(radius); 
    4552} 
    4653 
    4754 
    48 double Iq(double q, double dnn, 
     55static double Iq(double q, double dnn, 
    4956  double d_factor, double radius, 
    50   double sld, double solvent_sld){ 
     57  double sld, double solvent_sld) 
     58{ 
     59    // translate a point in [-1,1] to a point in [0, 2 pi] 
     60    const double phi_m = M_PI; 
     61    const double phi_b = M_PI; 
     62    // translate a point in [-1,1] to a point in [0, pi] 
     63    const double theta_m = M_PI_2; 
     64    const double theta_b = M_PI_2; 
    5165 
    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; 
     66    double outer_sum = 0.0; 
     67    for(int i=0; i<150; i++) { 
     68        double inner_sum = 0.0; 
     69        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     70        double sin_theta, cos_theta; 
     71        SINCOS(theta, sin_theta, cos_theta); 
     72        const double qc = q*cos_theta; 
     73        const double qab = q*sin_theta; 
     74        for(int j=0;j<150;j++) { 
     75            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     76            double sin_phi, cos_phi; 
     77            SINCOS(phi, sin_phi, cos_phi); 
     78            const double qa = qab*cos_phi; 
     79            const double qb = qab*sin_phi; 
     80            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     81            inner_sum += Gauss150Wt[j] * form; 
     82        } 
     83        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     84        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     85    } 
     86    outer_sum *= theta_m; 
     87    const double Zq = outer_sum/(4.0*M_PI); 
     88    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     89    return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
    8490} 
    8591 
    8692 
    87 double Iqxy(double qx, double qy, 
     93static double Iqxy(double qx, double qy, 
    8894    double dnn, double d_factor, double radius, 
    8995    double sld, double solvent_sld, 
     
    9298    double q, zhat, yhat, xhat; 
    9399    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     100    const double qa = q*xhat; 
     101    const double qb = q*yhat; 
     102    const double qc = q*zhat; 
    94103 
    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; 
     104    q = sqrt(qa*qa + qb*qb + qc*qc); 
     105    const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     106    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     107    return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
    111108} 
  • sasmodels/models/capped_cylinder.c

    r592343f r2a0b2b1  
    1414//   radius_cap is the radius of the lens 
    1515//   length is the cylinder length, or the separation between the lens halves 
    16 //   alpha is the angle of the cylinder wrt q. 
     16//   theta is the angle of the cylinder wrt q. 
    1717static double 
    18 _cap_kernel(double q, double h, double radius_cap, 
    19                       double half_length, double sin_alpha, double cos_alpha) 
     18_cap_kernel(double qab, double qc, double h, double radius_cap, 
     19            double half_length) 
    2020{ 
    2121    // translate a point in [-1,1] to a point in [lower,upper] 
     
    2626 
    2727    // cos term in integral is: 
    28     //    cos (q (R t - h + L/2) cos(alpha)) 
     28    //    cos (q (R t - h + L/2) cos(theta)) 
    2929    // so turn it into: 
    3030    //    cos (m t + b) 
    3131    // where: 
    32     //    m = q R cos(alpha) 
    33     //    b = q(L/2-h) cos(alpha) 
    34     const double m = q*radius_cap*cos_alpha; // cos argument slope 
    35     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    36     const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta) 
     32    //    m = q R cos(theta) 
     33    //    b = q(L/2-h) cos(theta) 
     34    const double m = radius_cap*qc; // cos argument slope 
     35    const double b = (half_length-h)*qc; // cos argument intercept 
     36    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3737    double total = 0.0; 
    3838    for (int i=0; i<76 ;i++) { 
    3939        const double t = Gauss76Z[i]*zm + zb; 
    4040        const double radical = 1.0 - t*t; 
    41         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     41        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    4242        const double Fq = cos(m*t + b) * radical * bj; 
    4343        total += Gauss76Wt[i] * Fq; 
     
    5050 
    5151static double 
    52 _fq(double q, double h, double radius_cap, double radius, double half_length, 
    53     double sin_alpha, double cos_alpha) 
     52_fq(double qab, double qc, double h, double radius_cap, double radius, double half_length) 
    5453{ 
    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); 
     54    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length); 
     55    const double bj = sas_2J1x_x(radius*qab); 
     56    const double si = sas_sinx_x(half_length*qc); 
    5857    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5958    const double Aq = cap_Fq + cyl_Fq; 
     
    101100    double total = 0.0; 
    102101    for (int i=0; i<76 ;i++) { 
    103         const double alpha= Gauss76Z[i]*zm + zb; 
    104         double sin_alpha, cos_alpha; // slots to hold sincos function output 
    105         SINCOS(alpha, sin_alpha, cos_alpha); 
    106  
    107         const double Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha); 
    108         // sin_alpha for spherical coord integration 
    109         total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
     102        const double theta = Gauss76Z[i]*zm + zb; 
     103        double sin_theta, cos_theta; // slots to hold sincos function output 
     104        SINCOS(theta, sin_theta, cos_theta); 
     105        const double qab = q*sin_theta; 
     106        const double qc = q*cos_theta; 
     107        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
     108        // scale by sin_theta for spherical coord integration 
     109        total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
    110110    } 
    111111    // translate dx in [-1,1] to dx in [lower,upper] 
     
    125125    double q, sin_alpha, cos_alpha; 
    126126    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     127    const double qab = q*sin_alpha; 
     128    const double qc = q*cos_alpha; 
    127129 
    128130    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    129     const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 
     131    const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length); 
    130132 
    131133    // Multiply by contrast^2 and convert to cm-1 
  • sasmodels/models/core_shell_bicelle.c

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

    rdedcf34 r2a0b2b1  
    1717        double thick_face, 
    1818        double length, 
    19         double rhoc, 
    20         double rhoh, 
    21         double rhor, 
    22         double rhosolv) 
     19        double sld_core, 
     20        double sld_face, 
     21        double sld_rim, 
     22        double sld_solvent) 
    2323{ 
    24     double si1,si2,be1,be2; 
    2524     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 
    2625     // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 
    27      //    const double uplim = M_PI_4; 
    2826    const double halfheight = 0.5*length; 
    29     //const double va = 0.0; 
    30     //const double vb = 1.0; 
    31     // inner integral limits 
    32     //const double vaj=0.0; 
    33     //const double vbj=M_PI; 
    34  
    3527    const double r_major = r_minor * x_core; 
    3628    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
    3729    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
    38     const double dr1 = (rhoc-rhoh)   *M_PI*r_minor*r_major*(2.0*halfheight);; 
    39     const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
    40     const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
    41     //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
    42     //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
    43     //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     30    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
     31    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
     32    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     33    const double dr1 = vol1*(sld_core-sld_face); 
     34    const double dr2 = vol2*(sld_rim-sld_solvent); 
     35    const double dr3 = vol3*(sld_face-sld_rim); 
    4436 
    4537    //initialize integral 
     
    4739    for(int i=0;i<76;i++) { 
    4840        //setup inner integral over the ellipsoidal cross-section 
    49         // since we generate these lots of times, why not store them somewhere? 
    50         //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    51         const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
    52         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    53         double inner_sum=0; 
    54         double sinarg1 = q*halfheight*cos_alpha; 
    55         double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 
    56         si1 = sas_sinx_x(sinarg1); 
    57         si2 = sas_sinx_x(sinarg2); 
     41        //const double va = 0.0; 
     42        //const double vb = 1.0; 
     43        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     44        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 
     45        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     46        const double qab = q*sin_theta; 
     47        const double qc = q*cos_theta; 
     48        const double si1 = sas_sinx_x(halfheight*qc); 
     49        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
     50        double inner_sum=0.0; 
    5851        for(int j=0;j<76;j++) { 
    5952            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    60             //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    61             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    62             const double rr = sqrt(r2A - r2B*cos(beta)); 
    63             double besarg1 = q*rr*sin_alpha; 
    64             double besarg2 = q*(rr+thick_rim)*sin_alpha; 
    65             be1 = sas_2J1x_x(besarg1); 
    66             be2 = sas_2J1x_x(besarg2); 
    67             inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 
    68                                               dr2*si2*be2 + 
    69                                               dr3*si2*be1); 
     53            // inner integral limits 
     54            //const double vaj=0.0; 
     55            //const double vbj=M_PI; 
     56            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     57            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 
     58            const double rr = sqrt(r2A - r2B*cos(phi)); 
     59            const double be1 = sas_2J1x_x(rr*qab); 
     60            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
     61            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
     62 
     63            inner_sum += Gauss76Wt[j] * fq * fq; 
    7064        } 
    7165        //now calculate outer integral 
     
    8377          double thick_face, 
    8478          double length, 
    85           double rhoc, 
    86           double rhoh, 
    87           double rhor, 
    88           double rhosolv, 
     79          double sld_core, 
     80          double sld_face, 
     81          double sld_rim, 
     82          double sld_solvent, 
    8983          double theta, 
    9084          double phi, 
    9185          double psi) 
    9286{ 
    93        // THIS NEEDS TESTING 
    9487    double q, xhat, yhat, zhat; 
    9588    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; 
     89    const double qa = q*xhat; 
     90    const double qb = q*yhat; 
     91    const double qc = q*zhat; 
     92 
     93    const double dr1 = sld_core-sld_face; 
     94    const double dr2 = sld_rim-sld_solvent; 
     95    const double dr3 = sld_face-sld_rim; 
    9996    const double r_major = r_minor*x_core; 
    10097    const double halfheight = 0.5*length; 
     
    104101 
    105102    // Compute effective radius in rotated coordinates 
    106     const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat)); 
    107     const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat) 
    108                                    + square((r_minor+thick_rim)*yhat)); 
    109     const double be1 = sas_2J1x_x( q*r_hat ); 
    110     const double be2 = sas_2J1x_x( q*rshell_hat ); 
    111     const double si1 = sas_sinx_x( q*halfheight*zhat ); 
    112     const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat ); 
    113     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1); 
    114     return 1.0e-4 * Aq; 
     103    const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb)); 
     104    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 
     105                                   + square((r_minor+thick_rim)*qb)); 
     106    const double be1 = sas_2J1x_x( qr_hat ); 
     107    const double be2 = sas_2J1x_x( qrshell_hat ); 
     108    const double si1 = sas_sinx_x( halfheight*qc ); 
     109    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc ); 
     110    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1; 
     111    return 1.0e-4 * fq*fq; 
    115112} 
    116  
  • sasmodels/models/core_shell_cylinder.c

    r592343f r2a0b2b1  
    1 double form_volume(double radius, double thickness, double length); 
    2 double Iq(double q, double core_sld, double shell_sld, double solvent_sld, 
    3     double radius, double thickness, double length); 
    4 double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld, 
    5     double radius, double thickness, double length, double theta, double phi); 
    6  
    71// vd = volume * delta_rho 
    8 // besarg = q * R * sin(alpha) 
    9 // siarg = q * L/2 * cos(alpha) 
    10 double _cyl(double vd, double besarg, double siarg); 
    11 double _cyl(double vd, double besarg, double siarg) 
     2// besarg = q * R * sin(theta) 
     3// siarg = q * L/2 * cos(theta) 
     4static double _cyl(double vd, double besarg, double siarg) 
    125{ 
    136    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 
    147} 
    158 
    16 double form_volume(double radius, double thickness, double length) 
     9static double 
     10form_volume(double radius, double thickness, double length) 
    1711{ 
    18     return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
     12    return M_PI*square(radius+thickness)*(length+2.0*thickness); 
    1913} 
    2014 
    21 double Iq(double q, 
     15static double 
     16Iq(double q, 
    2217    double core_sld, 
    2318    double shell_sld, 
     
    2823{ 
    2924    // precalculate constants 
    30     const double core_qr = q*radius; 
    31     const double core_qh = q*0.5*length; 
     25    const double core_r = radius; 
     26    const double core_h = 0.5*length; 
    3227    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    33     const double shell_qr = q*(radius + thickness); 
    34     const double shell_qh = q*(0.5*length + thickness); 
     28    const double shell_r = (radius + thickness); 
     29    const double shell_h = (0.5*length + thickness); 
    3530    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3631    double total = 0.0; 
    37     // double lower=0, upper=M_PI_2; 
    3832    for (int i=0; i<76 ;i++) { 
    39         // translate a point in [-1,1] to a point in [lower,upper] 
    40         //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    41         double sn, cn; 
    42         const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    43         SINCOS(alpha, sn, cn); 
    44         const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn) 
    45             + _cyl(shell_vd, shell_qr*sn, shell_qh*cn); 
    46         total += Gauss76Wt[i] * fq * fq * sn; 
     33        // translate a point in [-1,1] to a point in [0, pi/2] 
     34        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     35        double sin_theta, cos_theta; 
     36        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
     37        SINCOS(theta, sin_theta,  cos_theta); 
     38        const double qab = q*sin_theta; 
     39        const double qc = q*cos_theta; 
     40        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     41            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
     42        total += Gauss76Wt[i] * fq * fq * sin_theta; 
    4743    } 
    4844    // translate dx in [-1,1] to dx in [lower,upper] 
     
    6460    double q, sin_alpha, cos_alpha; 
    6561    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     62    const double qab = q*sin_alpha; 
     63    const double qc = q*cos_alpha; 
    6664 
    67     const double core_qr = q*radius; 
    68     const double core_qh = q*0.5*length; 
     65    const double core_r = radius; 
     66    const double core_h = 0.5*length; 
    6967    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    70     const double shell_qr = q*(radius + thickness); 
    71     const double shell_qh = q*(0.5*length + thickness); 
     68    const double shell_r = (radius + thickness); 
     69    const double shell_h = (0.5*length + thickness); 
    7270    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    7371 
    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); 
     72    const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     73        + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    7674    return 1.0e-4 * fq * fq; 
    7775} 
  • sasmodels/models/core_shell_ellipsoid.c

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

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

    r92dfe0c r2a0b2b1  
    1 double form_volume(double length_a, double length_b, double length_c,  
     1double form_volume(double length_a, double length_b, double length_c, 
    22                   double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    33double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 
     
    99            double thick_rim_c, double theta, double phi, double psi); 
    1010 
    11 double form_volume(double length_a, double length_b, double length_c,  
     11double form_volume(double length_a, double length_b, double length_c, 
    1212                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    1313{ 
    1414    //return length_a * length_b * length_c; 
    15     return length_a * length_b * length_c +  
    16            2.0 * thick_rim_a * length_b * length_c +  
     15    return length_a * length_b * length_c + 
     16           2.0 * thick_rim_a * length_b * length_c + 
    1717           2.0 * thick_rim_b * length_a * length_c + 
    1818           2.0 * thick_rim_c * length_a * length_b; 
     
    3434    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    3535    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    36      
     36 
    3737    const double mu = 0.5 * q * length_b; 
    38      
     38 
    3939    //calculate volume before rescaling (in original code, but not used) 
    40     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);         
    41     //double vol = length_a * length_b * length_c +  
    42     //       2.0 * thick_rim_a * length_b * length_c +  
     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 + 
    4343    //       2.0 * thick_rim_b * length_a * length_c + 
    4444    //       2.0 * thick_rim_c * length_a * length_b; 
    45      
     45 
    4646    // Scale sides by B 
    4747    const double a_scaled = length_a / length_b; 
     
    101101            //   ( 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 
    102102            // This is the function called by csparallelepiped_analytical_2D_scaled, 
    103             // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    104              
     103            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 
     104 
    105105            //  correct FF : sum of square of phase factors 
    106106            inner_total += Gauss76Wt[j] * form * form; 
     
    136136    double q, zhat, yhat, xhat; 
    137137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     138    const double qa = q*xhat; 
     139    const double qb = q*yhat; 
     140    const double qc = q*zhat; 
    138141 
    139142    // cspkernel in csparallelepiped recoded here 
     
    160163    double tc = length_a + 2.0*thick_rim_c; 
    161164    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162     double siA = sas_sinx_x(0.5*q*length_a*xhat); 
    163     double siB = sas_sinx_x(0.5*q*length_b*yhat); 
    164     double siC = sas_sinx_x(0.5*q*length_c*zhat); 
    165     double siAt = sas_sinx_x(0.5*q*ta*xhat); 
    166     double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    167     double siCt = sas_sinx_x(0.5*q*tc*zhat); 
    168      
     165    double siA = sas_sinx_x(0.5*length_a*qa); 
     166    double siB = sas_sinx_x(0.5*length_b*qb); 
     167    double siC = sas_sinx_x(0.5*length_c*qc); 
     168    double siAt = sas_sinx_x(0.5*ta*qa); 
     169    double siBt = sas_sinx_x(0.5*tb*qb); 
     170    double siCt = sas_sinx_x(0.5*tc*qc); 
     171 
    169172 
    170173    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
     
    174177               + drB*siA*(siBt-siB)*siC*V2 
    175178               + drC*siA*siB*(siCt*siCt-siC)*V3); 
    176     
     179 
    177180    return 1.0e-4 * f * f; 
    178181} 
  • sasmodels/models/cylinder.c

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

    r3b571ae r2a0b2b1  
    1 double form_volume(double radius_polar, double radius_equatorial); 
    2 double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial); 
    3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    4     double radius_polar, double radius_equatorial, double theta, double phi); 
    5  
    6 double form_volume(double radius_polar, double radius_equatorial) 
     1static double 
     2form_volume(double radius_polar, double radius_equatorial) 
    73{ 
    84    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    95} 
    106 
    11 double Iq(double q, 
     7static  double 
     8Iq(double q, 
    129    double sld, 
    1310    double sld_solvent, 
     
    4138} 
    4239 
    43 double Iqxy(double qx, double qy, 
     40static double 
     41Iqxy(double qx, double qy, 
    4442    double sld, 
    4543    double sld_solvent, 
     
    5149    double q, sin_alpha, cos_alpha; 
    5250    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); 
     51    const double qab = q*sin_alpha; 
     52    const double qc = q*cos_alpha; 
     53 
     54    const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc)); 
     55    const double f = sas_3j1x_x(qr); 
    5656    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    5757 
    5858    return 1.0e-4 * square(f * s); 
    5959} 
    60  
  • sasmodels/models/elliptical_cylinder.c

    r61104c8 r2a0b2b1  
    1 double form_volume(double radius_minor, double r_ratio, double length); 
    2 double Iq(double q, double radius_minor, double r_ratio, double length, 
    3           double sld, double solvent_sld); 
    4 double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length, 
    5             double sld, double solvent_sld, double theta, double phi, double psi); 
    6  
    7  
    8 double 
     1static double 
    92form_volume(double radius_minor, double r_ratio, double length) 
    103{ 
     
    125} 
    136 
    14 double 
     7static double 
    158Iq(double q, double radius_minor, double r_ratio, double length, 
    169   double sld, double solvent_sld) 
     
    6154 
    6255 
    63 double 
     56static double 
    6457Iqxy(double qx, double qy, 
    6558     double radius_minor, double r_ratio, double length, 
     
    6962    double q, xhat, yhat, zhat; 
    7063    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     64    const double qa = q*xhat; 
     65    const double qb = q*yhat; 
     66    const double qc = q*zhat; 
    7167 
    7268    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    7369    // Given:    radius_major = r_ratio * radius_minor 
    74     const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat)); 
    75     const double be = sas_2J1x_x(q*r); 
    76     const double si = sas_sinx_x(q*zhat*0.5*length); 
     70    const double qr = radius_minor*sqrt(square(r_ratio*qa) + square(qb)); 
     71    const double be = sas_2J1x_x(qr); 
     72    const double si = sas_sinx_x(qc*0.5*length); 
    7773    const double Aq = be * si; 
    7874    const double delrho = sld - solvent_sld; 
  • sasmodels/models/fcc_paracrystal.c

    r50beefe r7e0b281  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4#if 0  // Equations as written in Matsuoka 
     5    const double a1 = ( qa + qb)/2.0; 
     6    const double a2 = (-qa + qc)/2.0; 
     7    const double a3 = (-qa + qb)/2.0; 
     8#else 
     9    const double a1 = ( qa + qb)/2.0; 
     10    const double a2 = ( qa + qc)/2.0; 
     11    const double a3 = ( qb + qc)/2.0; 
     12#endif 
    613 
    7 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _FCCeval(double Theta, double Phi, double temp1, double temp3); 
     14    // Numerator: (1 - exp(a)^2)^3 
     15    //         => (-(exp(2a) - 1))^3 
     16    //         => -expm1(2a)^3 
     17    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 
     18    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     19    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 
     20    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     21    const double exp_arg = exp(arg); 
     22    const double Zq = -cube(expm1(2.0*arg)) 
     23        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     24          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     25          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     26 
     27    return Zq; 
     28} 
    929 
    1030 
    11 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    12  
    13         const double Da = d_factor*dnn; 
    14         const double temp1 = q*q*Da*Da; 
    15         const double temp3 = q*dnn; 
    16  
    17         double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    18         return(retVal); 
     31// occupied volume fraction calculated from lattice symmetry and sphere radius 
     32static double 
     33fcc_volume_fraction(double radius, double dnn) 
     34{ 
     35    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    1936} 
    2037 
    21 double _FCCeval(double Theta, double Phi, double temp1, double temp3) { 
    22  
    23         double result; 
    24         double sin_theta,cos_theta,sin_phi,cos_phi; 
    25         SINCOS(Theta, sin_theta, cos_theta); 
    26         SINCOS(Phi, sin_phi, cos_phi); 
    27  
    28         const double temp6 =  sin_theta; 
    29         const double temp7 =  sin_theta*sin_phi + cos_theta; 
    30         const double temp8 = -sin_theta*cos_phi + cos_theta; 
    31         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi; 
    32  
    33         const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
    34         result = cube(1.0-(temp10*temp10))*temp6 
    35             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    36               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    38  
    39         return (result); 
    40 } 
    41  
    42 double form_volume(double radius){ 
     38static double 
     39form_volume(double radius) 
     40{ 
    4341    return sphere_volume(radius); 
    4442} 
    4543 
    4644 
    47 double Iq(double q, double dnn, 
     45static double Iq(double q, double dnn, 
    4846  double d_factor, double radius, 
    49   double sld, double solvent_sld){ 
     47  double sld, double solvent_sld) 
     48{ 
     49    // translate a point in [-1,1] to a point in [0, 2 pi] 
     50    const double phi_m = M_PI; 
     51    const double phi_b = M_PI; 
     52    // translate a point in [-1,1] to a point in [0, pi] 
     53    const double theta_m = M_PI_2; 
     54    const double theta_b = M_PI_2; 
    5055 
    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); 
     56    double outer_sum = 0.0; 
     57    for(int i=0; i<150; i++) { 
     58        double inner_sum = 0.0; 
     59        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     60        double sin_theta, cos_theta; 
     61        SINCOS(theta, sin_theta, cos_theta); 
     62        const double qc = q*cos_theta; 
     63        const double qab = q*sin_theta; 
     64        for(int j=0;j<150;j++) { 
     65            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     66            double sin_phi, cos_phi; 
     67            SINCOS(phi, sin_phi, cos_phi); 
     68            const double qa = qab*cos_phi; 
     69            const double qb = qab*sin_phi; 
     70            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     71            inner_sum += Gauss150Wt[j] * form; 
     72        } 
     73        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     74        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     75    } 
     76    outer_sum *= theta_m; 
     77    const double Zq = outer_sum/(4.0*M_PI); 
     78    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    5479 
    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; 
     80    return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     81} 
    8382 
    8483 
    85 } 
    86  
    87 double Iqxy(double qx, double qy, 
     84static double Iqxy(double qx, double qy, 
    8885    double dnn, double d_factor, double radius, 
    8986    double sld, double solvent_sld, 
     
    9289    double q, zhat, yhat, xhat; 
    9390    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     91    const double qa = q*xhat; 
     92    const double qb = q*yhat; 
     93    const double qc = q*zhat; 
    9494 
    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; 
     95    q = sqrt(qa*qa + qb*qb + qc*qc); 
     96    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     97    const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     98    return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
    11299} 
  • sasmodels/models/hollow_cylinder.c

    r592343f r2a0b2b1  
    11double form_volume(double radius, double thickness, double length); 
    22double Iq(double q, double radius, double thickness, double length, double sld, 
    3         double solvent_sld); 
     3    double solvent_sld); 
    44double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld, 
    5         double solvent_sld, double theta, double phi); 
     5    double solvent_sld, double theta, double phi); 
    66 
    77//#define INVALID(v) (v.radius_core >= v.radius) 
     
    1414} 
    1515 
    16  
    1716static double 
    18 _hollow_cylinder_kernel(double q, 
    19     double radius, double thickness, double length, double sin_val, double cos_val) 
     17_fq(double qab, double qc, 
     18    double radius, double thickness, double length) 
    2019{ 
    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); 
     20    const double lam1 = sas_2J1x_x((radius+thickness)*qab); 
     21    const double lam2 = sas_2J1x_x(radius*qab); 
    2422    const double gamma_sq = square(radius/(radius+thickness)); 
    25     //Note: lim_{thickness -> 0} psi = sas_J0(radius*qs) 
    26     //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs) 
    27     const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 
    28     const double t2 = sas_sinx_x(0.5*q*length*cos_val); 
     23    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab) 
     24    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab) 
     25    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00 
     26    const double t2 = sas_sinx_x(0.5*length*qc); 
    2927    return psi*t2; 
    3028} 
     
    4341{ 
    4442    const double lower = 0.0; 
    45     const double upper = 1.0;           //limits of numerical integral 
     43    const double upper = 1.0;        //limits of numerical integral 
    4644 
    47     double summ = 0.0;                  //initialize intergral 
     45    double summ = 0.0;            //initialize intergral 
    4846    for (int i=0;i<76;i++) { 
    49         const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
    50         const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    51         const double inter = _hollow_cylinder_kernel(q, radius, thickness, length, 
    52                                                      sin_val, cos_val); 
    53         summ += Gauss76Wt[i] * inter * inter; 
     47        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     48        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     49        const double form = _fq(q*sin_theta, q*cos_theta, 
     50                                radius, thickness, length); 
     51        summ += Gauss76Wt[i] * form * form; 
    5452    } 
    5553 
     
    6664    double q, sin_alpha, cos_alpha; 
    6765    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); 
     66    const double qab = q*sin_alpha; 
     67    const double qc = q*cos_alpha; 
     68 
     69    const double form = _fq(qab, qc, radius, thickness, length); 
    7070 
    7171    const double vol = form_volume(radius, thickness, length); 
    72     return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol); 
     72    return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
    7373} 
    74  
  • sasmodels/models/parallelepiped.c

    rd605080 r2a0b2b1  
    2020{ 
    2121    const double mu = 0.5 * q * length_b; 
    22      
     22 
    2323    // Scale sides by B 
    2424    const double a_scaled = length_a / length_b; 
    2525    const double c_scaled = length_c / length_b; 
    26          
     26 
    2727    // outer integral (with gauss points), integration limits = 0, 1 
    2828    double outer_total = 0; //initialize integral 
     
    6969    double q, xhat, yhat, zhat; 
    7070    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     71    const double qa = q*xhat; 
     72    const double qb = q*yhat; 
     73    const double qc = q*zhat; 
    7174 
    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); 
     75    const double siA = sas_sinx_x(0.5*length_a*qa); 
     76    const double siB = sas_sinx_x(0.5*length_b*qb); 
     77    const double siC = sas_sinx_x(0.5*length_c*qc); 
    7578    const double V = form_volume(length_a, length_b, length_c); 
    7679    const double drho = (sld - solvent_sld); 
  • sasmodels/models/sc_paracrystal.c

    r50beefe r7e0b281  
    1 double form_volume(double radius); 
     1static double 
     2sc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    const double a1 = qa; 
     5    const double a2 = qb; 
     6    const double a3 = qc; 
    27 
    3 double Iq(double q, 
    4           double dnn, 
    5           double d_factor, 
    6           double radius, 
    7           double sphere_sld, 
    8           double solvent_sld); 
     8    // Numerator: (1 - exp(a)^2)^3 
     9    //         => (-(exp(2a) - 1))^3 
     10    //         => -expm1(2a)^3 
     11    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 
     12    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     13    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 
     14    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     15    const double exp_arg = exp(arg); 
     16    const double Zq = -cube(expm1(2.0*arg)) 
     17        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     18          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     19          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
    920 
    10 double Iqxy(double qx, double qy, 
    11             double dnn, 
    12             double d_factor, 
    13             double radius, 
    14             double sphere_sld, 
    15             double solvent_sld, 
    16             double theta, 
    17             double phi, 
    18             double psi); 
     21    return Zq; 
     22} 
    1923 
    20 double form_volume(double radius) 
     24// occupied volume fraction calculated from lattice symmetry and sphere radius 
     25static double 
     26sc_volume_fraction(double radius, double dnn) 
     27{ 
     28    return sphere_volume(radius/dnn); 
     29} 
     30 
     31static double 
     32form_volume(double radius) 
    2133{ 
    2234    return sphere_volume(radius); 
    2335} 
    2436 
    25 static double 
    26 sc_eval(double theta, double phi, double temp3, double temp4, double temp5) 
     37 
     38static double Iq(double q, double dnn, 
     39  double d_factor, double radius, 
     40  double sld, double solvent_sld) 
    2741{ 
    28     double cnt, snt; 
    29     SINCOS(theta, cnt, snt); 
     42    // translate a point in [-1,1] to a point in [0, 2 pi] 
     43    const double phi_m = M_PI_4; 
     44    const double phi_b = M_PI_4; 
     45    // translate a point in [-1,1] to a point in [0, pi] 
     46    const double theta_m = M_PI_4; 
     47    const double theta_b = M_PI_4; 
    3048 
    31     double cnp, snp; 
    32     SINCOS(phi, cnp, snp); 
    3349 
    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); 
     50    double outer_sum = 0.0; 
     51    for(int i=0; i<150; i++) { 
     52        double inner_sum = 0.0; 
     53        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     54        double sin_theta, cos_theta; 
     55        SINCOS(theta, sin_theta, cos_theta); 
     56        const double qc = q*cos_theta; 
     57        const double qab = q*sin_theta; 
     58        for(int j=0;j<150;j++) { 
     59            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     60            double sin_phi, cos_phi; 
     61            SINCOS(phi, sin_phi, cos_phi); 
     62            const double qa = qab*cos_phi; 
     63            const double qb = qab*sin_phi; 
     64            const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
     65            inner_sum += Gauss150Wt[j] * form; 
     66        } 
     67        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     68        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     69    } 
     70    outer_sum *= theta_m; 
     71    const double Zq = outer_sum/M_PI_2; 
     72    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     73 
     74    return sc_volume_fraction(radius, dnn) * Pq * Zq; 
    4275} 
    4376 
    44 static double 
    45 sc_integrand(double dnn, double d_factor, double qq, double xx, double yy) 
    46 { 
    47     //Function to calculate integrand values for simple cubic structure 
    4877 
    49         double da = d_factor*dnn; 
    50         double temp1 = qq*qq*da*da; 
    51         double temp2 = cube(-expm1(-temp1)); 
    52         double temp3 = qq*dnn; 
    53         double temp4 = 2.0*exp(-0.5*temp1); 
    54         double temp5 = exp(-1.0*temp1); 
    55  
    56         double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2; 
    57  
    58         return(integrand); 
    59 } 
    60  
    61 double Iq(double q, 
    62           double dnn, 
    63           double d_factor, 
    64           double radius, 
    65           double sphere_sld, 
    66           double solvent_sld) 
    67 { 
    68         const double va = 0.0; 
    69         const double vb = M_PI_2; //orientation average, outer integral 
    70  
    71     double summ=0.0; 
    72     double answer=0.0; 
    73  
    74         for(int i=0;i<150;i++) { 
    75                 //setup inner integral over the ellipsoidal cross-section 
    76                 double summj=0.0; 
    77                 double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; 
    78                 for(int j=0;j<150;j++) { 
    79                         //150 gauss points for the inner integral 
    80                         double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0; 
    81                         double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij); 
    82                         summj += tmp; 
    83                 } 
    84                 //now calculate the value of the inner integral 
    85                 answer = (vb-va)/2.0*summj; 
    86  
    87                 //now calculate outer integral 
    88                 double tmp = Gauss150Wt[i] * answer; 
    89                 summ += tmp; 
    90         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    91  
    92         answer = (vb-va)/2.0*summ; 
    93  
    94         //Volume fraction calculated from lattice symmetry and sphere radius 
    95         // NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3 
    96         const double latticeScale = sphere_volume(radius/dnn); 
    97  
    98         answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale; 
    99  
    100         return answer; 
    101 } 
    102  
    103 double Iqxy(double qx, double qy, 
    104           double dnn, 
    105           double d_factor, 
    106           double radius, 
    107           double sphere_sld, 
    108           double solvent_sld, 
    109           double theta, 
    110           double phi, 
    111           double psi) 
     78static double Iqxy(double qx, double qy, 
     79    double dnn, double d_factor, double radius, 
     80    double sld, double solvent_sld, 
     81    double theta, double phi, double psi) 
    11282{ 
    11383    double q, zhat, yhat, xhat; 
    11484    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     85    const double qa = q*xhat; 
     86    const double qb = q*yhat; 
     87    const double qc = q*zhat; 
    11588 
    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; 
     89    q = sqrt(qa*qa + qb*qb + qc*qc); 
     90    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     91    const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor); 
     92    return sc_volume_fraction(radius, dnn) * Pq * Zq; 
    12893} 
  • sasmodels/models/sc_paracrystal.py

    r8f04da4 r9bc4882  
    161161    [{'theta': 10.0, 'phi': 20, 'psi': 30.0}, (0.023, 0.045), 0.0177333171285], 
    162162    ] 
    163  
    164  
  • sasmodels/models/stacked_disks.c

    r19f996b rb34fc77  
    11static double stacked_disks_kernel( 
    2     double q, 
     2    double qab, 
     3    double qc, 
    34    double halfheight, 
    45    double thick_layer, 
     
    910    double layer_sld, 
    1011    double solvent_sld, 
    11     double sin_alpha, 
    12     double cos_alpha, 
    1312    double d) 
    1413 
     
    2019    // zi is the dummy variable for the integration (x in Feigin's notation) 
    2120 
    22     const double besarg1 = q*radius*sin_alpha; 
    23     //const double besarg2 = q*radius*sin_alpha; 
     21    const double besarg1 = radius*qab; 
     22    //const double besarg2 = radius*qab; 
    2423 
    25     const double sinarg1 = q*halfheight*cos_alpha; 
    26     const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 
     24    const double sinarg1 = halfheight*qc; 
     25    const double sinarg2 = (halfheight+thick_layer)*qc; 
    2726 
    2827    const double be1 = sas_2J1x_x(besarg1); 
     
    4342 
    4443    // loop for the structure factor S(q) 
    45     double qd_cos_alpha = q*d*cos_alpha; 
     44    double qd_cos_alpha = d*qc; 
    4645    //d*cos_alpha is the projection of d onto q (in other words the component 
    4746    //of d that is parallel to q. 
     
    8483        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8584        SINCOS(zi, sin_alpha, cos_alpha); 
    86         double yyy = stacked_disks_kernel(q, 
     85        double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha, 
    8786                           halfheight, 
    8887                           thick_layer, 
     
    9392                           layer_sld, 
    9493                           solvent_sld, 
    95                            sin_alpha, 
    96                            cos_alpha, 
    9794                           d); 
    9895        summ += Gauss76Wt[i] * yyy * sin_alpha; 
     
    152149    double phi) 
    153150{ 
    154     int n_stacking = (int)(fp_n_stacking + 0.5); 
    155151    double q, sin_alpha, cos_alpha; 
    156152    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     153    const double qab = q*sin_alpha; 
     154    const double qc = q*cos_alpha; 
    157155 
     156    int n_stacking = (int)(fp_n_stacking + 0.5); 
    158157    double d = 2.0 * thick_layer + thick_core; 
    159158    double halfheight = 0.5*thick_core; 
    160     double answer = stacked_disks_kernel(q, 
     159    double answer = stacked_disks_kernel(qab, qc, 
    161160                     halfheight, 
    162161                     thick_layer, 
     
    167166                     layer_sld, 
    168167                     solvent_sld, 
    169                      sin_alpha, 
    170                      cos_alpha, 
    171168                     d); 
    172169 
  • sasmodels/models/triaxial_ellipsoid.c

    r68dd6a9 r2a0b2b1  
    1 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar); 
    2 double Iq(double q, double sld, double sld_solvent, 
    3     double radius_equat_minor, double radius_equat_major, double radius_polar); 
    4 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    5     double radius_equat_minor, double radius_equat_major, double radius_polar, double theta, double phi, double psi); 
    6  
    71//#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar) 
    82 
    9  
    10 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     3static double 
     4form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
    115{ 
    126    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
    137} 
    148 
    15 double Iq(double q, 
     9static double 
     10Iq(double q, 
    1611    double sld, 
    1712    double sld_solvent, 
     
    4540    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    4641    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
    47     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    48     return 1.0e-4 * s * s * fqsq; 
     42    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     43    const double drho = (sld - sld_solvent); 
     44    return 1.0e-4 * square(vol*drho) * fqsq; 
    4945} 
    5046 
    51 double Iqxy(double qx, double qy, 
     47static double 
     48Iqxy(double qx, double qy, 
    5249    double sld, 
    5350    double sld_solvent, 
     
    6158    double q, xhat, yhat, zhat; 
    6259    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     60    const double qa = q*xhat; 
     61    const double qb = q*yhat; 
     62    const double qc = q*zhat; 
    6363 
    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); 
     64    const double qr = sqrt(square(radius_equat_minor*qa) 
     65                           + square(radius_equat_major*qb) 
     66                           + square(radius_polar*qc)); 
     67    const double fq = sas_3j1x_x(qr); 
     68    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     69    const double drho = (sld - sld_solvent); 
    6970 
    70     return 1.0e-4 * square(s * fq); 
     71    return 1.0e-4 * square(vol * drho * fq); 
    7172} 
    72  
  • sasmodels/weights.py

    r41e7f2e rd4c33d6  
    5555        """ 
    5656        sigma = self.width * center if relative else self.width 
     57        if not relative: 
     58            # For orientation, the jitter is relative to 0 not the angle 
     59            #center = 0 
     60            pass 
    5761        if sigma == 0 or self.npts < 2: 
    5862            if lb <= center <= ub: 
Note: See TracChangeset for help on using the changeset viewer.