Changeset cdcd80f in sasmodels


Ignore:
Timestamp:
Sep 25, 2018 9:58:34 AM (4 weeks ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
ticket-1084
Parents:
bfe4245 (diff), 2015f02 (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-1084

Files:
1 added
1 deleted
36 edited

Legend:

Unmodified
Added
Removed
  • .travis.yml

    r335271e r5c36bf1  
    66    env: 
    77    - PY=2.7 
    8     - DEPLOY=True 
     8    #- DEPLOY=True 
    99  - os: linux 
    1010    env: 
  • doc/guide/gpu_setup.rst

    r59485a4 r63602b1  
    139139the compiler. 
    140140 
    141 On Windows, set *SASCOMPILER=tinycc* for the tinycc compiler, 
    142 *SASCOMPILER=msvc* for the Microsoft Visual C compiler, 
    143 or *SASCOMPILER=mingw* for the MinGW compiler. If TinyCC is available 
     141On Windows, set *SAS_COMPILER=tinycc* for the tinycc compiler, 
     142*SAS_COMPILER=msvc* for the Microsoft Visual C compiler, 
     143or *SAS_COMPILER=mingw* for the MinGW compiler. If TinyCC is available 
    144144on the python path (it is provided with SasView), that will be the 
    145145default. If you want one of the other compilers, be sure to have it 
  • doc/guide/magnetism/magnetism.rst

    r4f5afc9 rbefe905  
    3939 
    4040.. math:: 
    41     -- &= ((1-u_i)(1-u_f))^{1/4} \\ 
    42     -+ &= ((1-u_i)(u_f))^{1/4} \\ 
    43     +- &= ((u_i)(1-u_f))^{1/4} \\ 
    44     ++ &= ((u_i)(u_f))^{1/4} 
     41    -- &= (1-u_i)(1-u_f) \\ 
     42    -+ &= (1-u_i)(u_f) \\ 
     43    +- &= (u_i)(1-u_f) \\ 
     44    ++ &= (u_i)(u_f) 
    4545 
    4646Ideally the experiment would measure the pure spin states independently and 
     
    104104| 2015-05-02 Steve King 
    105105| 2017-11-15 Paul Kienzle 
     106| 2018-06-02 Adam Washington 
  • doc/guide/pd/polydispersity.rst

    r29afc50 rd089a00  
    88.. _polydispersityhelp: 
    99 
    10 Polydispersity Distributions 
    11 ---------------------------- 
    12  
    13 With some models in sasmodels we can calculate the average intensity for a 
    14 population of particles that exhibit size and/or orientational 
    15 polydispersity. The resultant intensity is normalized by the average 
    16 particle volume such that 
     10Polydispersity & Orientational Distributions 
     11-------------------------------------------- 
     12 
     13For some models we can calculate the average intensity for a population of  
     14particles that possess size and/or orientational (ie, angular) distributions.  
     15In SasView we call the former *polydispersity* but use the parameter *PD* to  
     16parameterise both. In other words, the meaning of *PD* in a model depends on  
     17the actual parameter it is being applied too. 
     18 
     19The resultant intensity is then normalized by the average particle volume such  
     20that 
    1721 
    1822.. math:: 
     
    2024  P(q) = \text{scale} \langle F^* F \rangle / V + \text{background} 
    2125 
    22 where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an 
    23 average over the size distribution. 
     26where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an  
     27average over the distribution $f(x; \bar x, \sigma)$, giving 
     28 
     29.. math:: 
     30 
     31  P(q) = \frac{\text{scale}}{V} \int_\mathbb{R}  
     32  f(x; \bar x, \sigma) F^2(q, x)\, dx + \text{background} 
    2433 
    2534Each distribution is characterized by a center value $\bar x$ or 
    2635$x_\text{med}$, a width parameter $\sigma$ (note this is *not necessarily* 
    27 the standard deviation, so read the description carefully), the number of 
    28 sigmas $N_\sigma$ to include from the tails of the distribution, and the 
    29 number of points used to compute the average. The center of the distribution 
    30 is set by the value of the model parameter. The meaning of a polydispersity  
    31 parameter *PD* (not to be confused with a molecular weight distributions  
    32 in polymer science) in a model depends on the type of parameter it is being  
    33 applied too. 
     36the standard deviation, so read the description of the distribution carefully),  
     37the number of sigmas $N_\sigma$ to include from the tails of the distribution,  
     38and the number of points used to compute the average. The center of the  
     39distribution is set by the value of the model parameter. 
    3440 
    3541The distribution width applied to *volume* (ie, shape-describing) parameters  
    3642is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$.  
    37 However, the distribution width applied to *orientation* (ie, angle-describing)  
    38 parameters is just $\sigma = \mathrm{PD}$. 
     43However, the distribution width applied to *orientation* parameters is just  
     44$\sigma = \mathrm{PD}$. 
    3945 
    4046$N_\sigma$ determines how far into the tails to evaluate the distribution, 
    4147with larger values of $N_\sigma$ required for heavier tailed distributions. 
    4248The scattering in general falls rapidly with $qr$ so the usual assumption 
    43 that $G(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)G(r - 3\sigma_r)$ 
     49that $f(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)f(r - 3\sigma_r)$ 
    4450will not contribute much to the average may not hold when particles are large. 
    4551This, too, will require increasing $N_\sigma$. 
    4652 
    4753Users should note that the averaging computation is very intensive. Applying 
    48 polydispersion to multiple parameters at the same time or increasing the 
    49 number of points in the distribution will require patience! However, the 
    50 calculations are generally more robust with more data points or more angles. 
     54polydispersion and/or orientational distributions to multiple parameters at  
     55the same time, or increasing the number of points in the distribution, will  
     56require patience! However, the calculations are generally more robust with  
     57more data points or more angles. 
    5158 
    5259The following distribution functions are provided: 
     
    6471Additional distributions are under consideration. 
    6572 
     73**Beware: when the Polydispersity & Orientational Distribution panel in SasView is** 
     74**first opened, the default distribution for all parameters is the Gaussian Distribution.** 
     75**This may not be suitable. See Suggested Applications below.** 
     76 
     77.. note:: In 2009 IUPAC decided to introduce the new term 'dispersity' to replace  
     78           the term 'polydispersity' (see `Pure Appl. Chem., (2009), 81(2),  
     79           351-353 <http://media.iupac.org/publications/pac/2009/pdf/8102x0351.pdf>`_  
     80           in order to make the terminology describing distributions of chemical  
     81           properties unambiguous. However, these terms are unrelated to the  
     82           proportional size distributions and orientational distributions used in  
     83           SasView models. 
     84 
    6685Suggested Applications 
    6786^^^^^^^^^^^^^^^^^^^^^^ 
    6887 
    69 If applying polydispersion to parameters describing particle sizes, use 
     88If applying polydispersion to parameters describing particle sizes, consider using 
    7089the Lognormal or Schulz distributions. 
    7190 
    7291If applying polydispersion to parameters describing interfacial thicknesses 
    73 or angular orientations, use the Gaussian or Boltzmann distributions. 
     92or angular orientations, consider using the Gaussian or Boltzmann distributions. 
    7493 
    7594If applying polydispersion to parameters describing angles, use the Uniform  
     
    318337^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
    319338 
    320 Many commercial Dynamic Light Scattering (DLS) instruments produce a size 
    321 polydispersity parameter, sometimes even given the symbol $p$\ ! This 
    322 parameter is defined as the relative standard deviation coefficient of 
    323 variation of the size distribution and is NOT the same as the polydispersity 
    324 parameters in the Lognormal and Schulz distributions above (though they all 
    325 related) except when the DLS polydispersity parameter is <0.13. 
    326  
    327 .. math:: 
    328  
    329     p_{DLS} = \sqrt(\nu / \bar x^2) 
    330  
    331 where $\nu$ is the variance of the distribution and $\bar x$ is the mean 
    332 value of $x$. 
     339Several measures of polydispersity abound in Dynamic Light Scattering (DLS) and  
     340it should not be assumed that any of the following can be simply equated with  
     341the polydispersity *PD* parameter used in SasView. 
     342 
     343The dimensionless **Polydispersity Index (PI)** is a measure of the width of the  
     344distribution of autocorrelation function decay rates (*not* the distribution of  
     345particle sizes itself, though the two are inversely related) and is defined by  
     346ISO 22412:2017 as 
     347 
     348.. math:: 
     349 
     350    PI = \mu_{2} / \bar \Gamma^2 
     351 
     352where $\mu_\text{2}$ is the second cumulant, and $\bar \Gamma^2$ is the  
     353intensity-weighted average value, of the distribution of decay rates. 
     354 
     355*If the distribution of decay rates is Gaussian* then 
     356 
     357.. math:: 
     358 
     359    PI = \sigma^2 / 2\bar \Gamma^2 
     360 
     361where $\sigma$ is the standard deviation, allowing a **Relative Polydispersity (RP)**  
     362to be defined as 
     363 
     364.. math:: 
     365 
     366    RP = \sigma / \bar \Gamma = \sqrt{2 \cdot PI} 
     367 
     368PI values smaller than 0.05 indicate a highly monodisperse system. Values  
     369greater than 0.7 indicate significant polydispersity. 
     370 
     371The **size polydispersity P-parameter** is defined as the relative standard  
     372deviation coefficient of variation   
     373 
     374.. math:: 
     375 
     376    P = \sqrt\nu / \bar R 
     377 
     378where $\nu$ is the variance of the distribution and $\bar R$ is the mean 
     379value of $R$. Here, the product $P \bar R$ is *equal* to the standard  
     380deviation of the Lognormal distribution. 
     381 
     382P values smaller than 0.13 indicate a monodisperse system. 
    333383 
    334384For more information see: 
    335 S King, C Washington & R Heenan, *Phys Chem Chem Phys*, (2005), 7, 143 
     385 
     386`ISO 22412:2017, International Standards Organisation (2017) <https://www.iso.org/standard/65410.html>`_. 
     387 
     388`Polydispersity: What does it mean for DLS and Chromatography <http://www.materials-talks.com/blog/2014/10/23/polydispersity-what-does-it-mean-for-dls-and-chromatography/>`_. 
     389 
     390`Dynamic Light Scattering: Common Terms Defined, Whitepaper WP111214. Malvern Instruments (2011) <http://www.biophysics.bioc.cam.ac.uk/wp-content/uploads/2011/02/DLS_Terms_defined_Malvern.pdf>`_. 
     391 
     392S King, C Washington & R Heenan, *Phys Chem Chem Phys*, (2005), 7, 143. 
     393 
     394T Allen, in *Particle Size Measurement*, 4th Edition, Chapman & Hall, London (1990). 
    336395 
    337396.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     
    343402| 2018-03-20 Steve King 
    344403| 2018-04-04 Steve King 
     404| 2018-08-09 Steve King 
  • doc/guide/plugin.rst

    r7e6bc45e r2015f02  
    423423calculations, but instead rely on numerical integration to compute the 
    424424appropriately smeared pattern. 
     425 
     426Each .py file also contains a function:: 
     427 
     428        def random(): 
     429        ... 
     430         
     431This function provides a model-specific random parameter set which shows model  
     432features in the USANS to SANS range.  For example, core-shell sphere sets the  
     433outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q  
     434value for the transition from flat to falling.  It then uses a beta distribution  
     435to set the percentage of the shape which is shell, giving a preference for very  
     436thin or very thick shells (but never 0% or 100%).  Using `-sets=10` in sascomp  
     437should show a reasonable variety of curves over the default sascomp q range.   
     438The parameter set is returned as a dictionary of `{parameter: value, ...}`.   
     439Any model parameters not included in the dictionary will default according to  
     440the code in the `_randomize_one()` function from sasmodels/compare.py. 
    425441 
    426442Python Models 
     
    822838 
    823839        :code:`source = ["lib/Si.c", ...]` 
    824         (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/Si.c>`_) 
     840        (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_Si.c>`_) 
    825841 
    826842    sas_3j1x_x(x): 
  • doc/guide/scripting.rst

    r4aa5dce rbd7630d  
    1010The key functions are :func:`sasmodels.core.load_model` for loading the 
    1111model definition and compiling the kernel and 
    12 :func:`sasmodels.data.load_data` for calling sasview to load the data. Need 
    13 the data because that defines the resolution function and the q values to 
    14 evaluate. If there is no data, then use :func:`sasmodels.data.empty_data1D` 
    15 or :func:`sasmodels.data.empty_data2D` to create some data with a given $q$. 
    16  
    17 Using sasmodels through bumps 
    18 ============================= 
    19  
    20 With the data and the model, you can wrap it in a *bumps* model with 
     12:func:`sasmodels.data.load_data` for calling sasview to load the data. 
     13 
     14Preparing data 
     15============== 
     16 
     17Usually you will load data via the sasview loader, with the 
     18:func:`sasmodels.data.load_data` function.  For example:: 
     19 
     20    from sasmodels.data import load_data 
     21    data = load_data("sasmodels/example/093191_201.dat") 
     22 
     23You may want to apply a data mask, such a beam stop, and trim high $q$:: 
     24 
     25    from sasmodels.data import set_beam_stop 
     26    set_beam_stop(data, qmin, qmax) 
     27 
     28The :func:`sasmodels.data.set_beam_stop` method simply sets the *mask* 
     29attribute for the data. 
     30 
     31The data defines the resolution function and the q values to evaluate, so 
     32even if you simulating experiments prior to making measurements, you still 
     33need a data object for reference. Use :func:`sasmodels.data.empty_data1D` 
     34or :func:`sasmodels.data.empty_data2D` to create a container with a 
     35given $q$ and $\Delta q/q$.  For example:: 
     36 
     37    import numpy as np 
     38    from sasmodels.data import empty_data1D 
     39 
     40    # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5% 
     41    q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120) 
     42    data = empty_data1D(q, resolution=0.05) 
     43 
     44To use a more realistic model of resolution, or to load data from a file 
     45format not understood by SasView, you can use :class:`sasmodels.data.Data1D` 
     46or :class:`sasmodels.data.Data2D` directly.  The 1D data uses 
     47*x*, *y*, *dx* and *dy* for $x = q$ and $y = I(q)$, and 2D data uses 
     48*x*, *y*, *z*, *dx*, *dy*, *dz* for $x, y = qx, qy$ and $z = I(qx, qy)$. 
     49[Note: internally, the Data2D object uses SasView conventions, 
     50*qx_data*, *qy_data*, *data*, *dqx_data*, *dqy_data*, and *err_data*.] 
     51 
     52For USANS data, use 1D data, but set *dxl* and *dxw* attributes to 
     53indicate slit resolution:: 
     54 
     55    data.dxl = 0.117 
     56 
     57See :func:`sasmodels.resolution.slit_resolution` for details. 
     58 
     59SESANS data is more complicated; if your SESANS format is not supported by 
     60SasView you need to define a number of attributes beyond *x*, *y*.  For 
     61example:: 
     62 
     63    SElength = np.linspace(0, 2400, 61) # [A] 
     64    data = np.ones_like(SElength) 
     65    err_data = np.ones_like(SElength)*0.03 
     66 
     67    class Source: 
     68        wavelength = 6 # [A] 
     69        wavelength_unit = "A" 
     70    class Sample: 
     71        zacceptance = 0.1 # [A^-1] 
     72        thickness = 0.2 # [cm] 
     73 
     74    class SESANSData1D: 
     75        #q_zmax = 0.23 # [A^-1] 
     76        lam = 0.2 # [nm] 
     77        x = SElength 
     78        y = data 
     79        dy = err_data 
     80        sample = Sample() 
     81    data = SESANSData1D() 
     82 
     83    x, y = ... # create or load sesans 
     84    data = smd.Data 
     85 
     86The *data* module defines various data plotters as well. 
     87 
     88Using sasmodels directly 
     89======================== 
     90 
     91Once you have a computational kernel and a data object, you can evaluate 
     92the model for various parameters using 
     93:class:`sasmodels.direct_model.DirectModel`.  The resulting object *f* 
     94will be callable as *f(par=value, ...)*, returning the $I(q)$ for the $q$ 
     95values in the data.  For example:: 
     96 
     97    import numpy as np 
     98    from sasmodels.data import empty_data1D 
     99    from sasmodels.core import load_model 
     100    from sasmodels.direct_model import DirectModel 
     101 
     102    # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5% 
     103    q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120) 
     104    data = empty_data1D(q, resolution=0.05) 
     105    kernel = load_model("ellipsoid) 
     106    f = DirectModel(data, kernel) 
     107    Iq = f(radius_polar=100) 
     108 
     109Polydispersity information is set with special parameter names: 
     110 
     111    * *par_pd* for polydispersity width, $\Delta p/p$, 
     112    * *par_pd_n* for the number of points in the distribution, 
     113    * *par_pd_type* for the distribution type (as a string), and 
     114    * *par_pd_nsigmas* for the limits of the distribution. 
     115 
     116Using sasmodels through the bumps optimizer 
     117=========================================== 
     118 
     119Like DirectModel, you can wrap data and a kernel in a *bumps* model with 
    21120class:`sasmodels.bumps_model.Model` and create an 
    22 class:`sasmodels.bump_model.Experiment` that you can fit with the *bumps* 
     121class:`sasmodels.bumps_model.Experiment` that you can fit with the *bumps* 
    23122interface. Here is an example from the *example* directory such as 
    24123*example/model.py*:: 
     
    75174    SasViewCom bumps.cli example/model.py --preview 
    76175 
    77 Using sasmodels directly 
    78 ======================== 
    79  
    80 Bumps has a notion of parameter boxes in which you can set and retrieve 
    81 values.  Instead of using bumps, you can create a directly callable function 
    82 with :class:`sasmodels.direct_model.DirectModel`.  The resulting object *f* 
    83 will be callable as *f(par=value, ...)*, returning the $I(q)$ for the $q$ 
    84 values in the data.  Polydisperse parameters use the same naming conventions 
    85 as in the bumps model, with e.g., radius_pd being the polydispersity associated 
    86 with radius. 
     176Calling the computation kernel 
     177============================== 
    87178 
    88179Getting a simple function that you can call on a set of q values and return 
  • doc/rst_prolog

    r30b60d2 r2c12061  
    99.. |Ang^-3| replace:: |Ang|\ :sup:`-3` 
    1010.. |Ang^-4| replace:: |Ang|\ :sup:`-4` 
     11.. |nm^-1| replace:: nm\ :sup:`-1` 
    1112.. |cm^-1| replace:: cm\ :sup:`-1` 
    1213.. |cm^2| replace:: cm\ :sup:`2` 
  • example/model_ellipsoid_hayter_msa.py

    r8a5f021 r9a99993  
    1616 
    1717# DEFINE THE MODEL 
    18 kernel = load_model('ellipsoid*hayter_msa') 
     18kernel = load_model('ellipsoid@hayter_msa') 
    1919 
    2020pars = dict(scale=6.4, background=0.06, sld=0.33, sld_solvent=2.15, radius_polar=14.0, 
  • sasmodels/compare.py

    r1fbadb2 rbd7630d  
    107107    -title="note" adds note to the plot title, after the model name 
    108108    -weights shows weights plots for the polydisperse parameters 
     109    -profile shows the sld profile if the model has a plottable sld profile 
    109110 
    110111    === output options === 
    111112    -edit starts the parameter explorer 
    112113    -help/-html shows the model docs instead of running the model 
     114 
     115    === environment variables === 
     116    -DSAS_MODELPATH=path sets directory containing custom models 
     117    -DSAS_OPENCL=vendor:device|none sets the target OpenCL device 
     118    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only) 
     119    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler 
     120    -DSAS_OPENMP=1 turns on OpenMP for the DLLs 
     121    -DSAS_DLL_PATH=path sets the path to the compiled modules 
    113122 
    114123The interpretation of quad precision depends on architecture, and may 
     
    645654 
    646655def make_data(opts): 
    647     # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray] 
     656    # type: (Dict[str, Any], float) -> Tuple[Data, np.ndarray] 
    648657    """ 
    649658    Generate an empty dataset, used with the model to set Q points 
     
    667676        if opts['zero']: 
    668677            q = np.hstack((0, q)) 
    669         data = empty_data1D(q, resolution=res) 
     678        # TODO: provide command line control of lambda and Delta lambda/lambda 
     679        #L, dLoL = 5, 0.14/np.sqrt(6)  # wavelength and 14% triangular FWHM 
     680        L, dLoL = 0, 0 
     681        data = empty_data1D(q, resolution=res, L=L, dL=L*dLoL) 
    670682        index = slice(None, None) 
    671683    return data, index 
     
    772784            dim = base._kernel.dim 
    773785            plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim)) 
     786        if opts['show_profile']: 
     787            import pylab 
     788            base, comp = opts['engines'] 
     789            base_pars, comp_pars = opts['pars'] 
     790            have_base = base._kernel.info.profile is not None 
     791            have_comp = ( 
     792                comp is not None 
     793                and comp._kernel.info.profile is not None 
     794                and base_pars != comp_pars 
     795            ) 
     796            if have_base or have_comp: 
     797                pylab.figure() 
     798            if have_base: 
     799                plot_profile(base._kernel.info, **base_pars) 
     800            if have_comp: 
     801                plot_profile(comp._kernel.info, label='comp', **comp_pars) 
     802                pylab.legend() 
    774803    if opts['plot']: 
    775804        import matplotlib.pyplot as plt 
     
    777806    return limits 
    778807 
     808def plot_profile(model_info, label='base', **args): 
     809    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None 
     810    """ 
     811    Plot the profile returned by the model profile method. 
     812 
     813    *model_info* defines model parameters, etc. 
     814 
     815    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*) 
     816    for each parameter, where (*dispersity*, *weights*) pairs are the 
     817    distributions to be plotted. 
     818    """ 
     819    import pylab 
     820 
     821    args = dict((k, v) for k, v in args.items() 
     822                if "_pd" not in k 
     823                   and ":" not in k 
     824                   and k not in ("background", "scale", "theta", "phi", "psi")) 
     825    args = args.copy() 
     826 
     827    args.pop('scale', 1.) 
     828    args.pop('background', 0.) 
     829    z, rho = model_info.profile(**args) 
     830    #pylab.interactive(True) 
     831    pylab.plot(z, rho, '-', label=label) 
     832    pylab.grid(True) 
     833    #pylab.show() 
     834 
     835 
     836 
    779837def run_models(opts, verbose=False): 
    780838    # type: (Dict[str, Any]) -> Dict[str, Any] 
     
    786844    base_n, comp_n = opts['count'] 
    787845    base_pars, comp_pars = opts['pars'] 
    788     data = opts['data'] 
     846    base_data, comp_data = opts['data'] 
    789847 
    790848    comparison = comp is not None 
     
    800858            print("%s t=%.2f ms, intensity=%.0f" 
    801859                  % (base.engine, base_time, base_value.sum())) 
    802         _show_invalid(data, base_value) 
     860        _show_invalid(base_data, base_value) 
    803861    except ImportError: 
    804862        traceback.print_exc() 
     
    812870                print("%s t=%.2f ms, intensity=%.0f" 
    813871                      % (comp.engine, comp_time, comp_value.sum())) 
    814             _show_invalid(data, comp_value) 
     872            _show_invalid(base_data, comp_value) 
    815873        except ImportError: 
    816874            traceback.print_exc() 
     
    866924    have_base, have_comp = (base_value is not None), (comp_value is not None) 
    867925    base, comp = opts['engines'] 
    868     data = opts['data'] 
     926    base_data, comp_data = opts['data'] 
    869927    use_data = (opts['datafile'] is not None) and (have_base ^ have_comp) 
    870928 
    871929    # Plot if requested 
    872930    view = opts['view'] 
     931    #view = 'log' 
    873932    if limits is None: 
    874933        vmin, vmax = np.inf, -np.inf 
     
    884943        if have_comp: 
    885944            plt.subplot(131) 
    886         plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
     945        plot_theory(base_data, base_value, view=view, use_data=use_data, limits=limits) 
    887946        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    888947        #cbar_title = "log I" 
     
    891950            plt.subplot(132) 
    892951        if not opts['is2d'] and have_base: 
    893             plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
    894         plot_theory(data, comp_value, view=view, use_data=use_data, limits=limits) 
     952            plot_theory(comp_data, base_value, view=view, use_data=use_data, limits=limits) 
     953        plot_theory(comp_data, comp_value, view=view, use_data=use_data, limits=limits) 
    895954        plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 
    896955        #cbar_title = "log I" 
     
    908967            err[err > cutoff] = cutoff 
    909968        #err,errstr = base/comp,"ratio" 
    910         plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
     969        # Note: base_data only since base and comp have same q values (though 
     970        # perhaps different resolution), and we are plotting the difference 
     971        # at each q 
     972        plot_theory(base_data, None, resid=err, view=errview, use_data=use_data) 
    911973        plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear') 
    912974        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best') 
     
    9421004OPTIONS = [ 
    9431005    # Plotting 
    944     'plot', 'noplot', 'weights', 
     1006    'plot', 'noplot', 
     1007    'weights', 'profile', 
    9451008    'linear', 'log', 'q4', 
    9461009    'rel', 'abs', 
     
    10581121 
    10591122    invalid = [o[1:] for o in flags 
    1060                if o[1:] not in NAME_OPTIONS 
    1061                and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] 
     1123               if not (o[1:] in NAME_OPTIONS 
     1124                       or any(o.startswith('-%s='%t) for t in VALUE_OPTIONS) 
     1125                       or o.startswith('-D'))] 
    10621126    if invalid: 
    10631127        print("Invalid options: %s"%(", ".join(invalid))) 
     
    10751139        'qmax'      : 0.05, 
    10761140        'nq'        : 128, 
    1077         'res'       : 0.0, 
     1141        'res'       : '0.0', 
    10781142        'noise'     : 0.0, 
    10791143        'accuracy'  : 'Low', 
     
    10961160        'count'     : '1', 
    10971161        'show_weights' : False, 
     1162        'show_profile' : False, 
    10981163        'sphere'    : 0, 
    10991164        'ngauss'    : '0', 
     
    11151180        elif arg.startswith('-q='): 
    11161181            opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')] 
    1117         elif arg.startswith('-res='):      opts['res'] = float(arg[5:]) 
     1182        elif arg.startswith('-res='):      opts['res'] = arg[5:] 
    11181183        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:]) 
    11191184        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:]) 
     
    11561221        elif arg == '-default': opts['use_demo'] = False 
    11571222        elif arg == '-weights': opts['show_weights'] = True 
     1223        elif arg == '-profile': opts['show_profile'] = True 
    11581224        elif arg == '-html':    opts['html'] = True 
    11591225        elif arg == '-help':    opts['html'] = True 
     1226        elif arg.startswith('-D'): 
     1227            var, val = arg[2:].split('=') 
     1228            os.environ[var] = val 
    11601229    # pylint: enable=bad-whitespace,C0321 
    11611230 
     
    11731242    if opts['qmin'] is None: 
    11741243        opts['qmin'] = 0.001*opts['qmax'] 
    1175     if opts['datafile'] is not None: 
    1176         data = load_data(os.path.expanduser(opts['datafile'])) 
    1177     else: 
    1178         data, _ = make_data(opts) 
    11791244 
    11801245    comparison = any(PAR_SPLIT in v for v in values) 
     
    12161281        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12171282 
    1218     base = make_engine(model_info[0], data, opts['engine'][0], 
     1283    if PAR_SPLIT in opts['res']: 
     1284        opts['res'] = [float(k) for k in opts['res'].split(PAR_SPLIT, 2)] 
     1285        comparison = True 
     1286    else: 
     1287        opts['res'] = [float(opts['res'])]*2 
     1288 
     1289    if opts['datafile'] is not None: 
     1290        data0 = load_data(os.path.expanduser(opts['datafile'])) 
     1291        data = data0, data0 
     1292    else: 
     1293        # Hack around the fact that make_data doesn't take a pair of resolutions 
     1294        res = opts['res'] 
     1295        opts['res'] = res[0] 
     1296        data0, _ = make_data(opts) 
     1297        if res[0] != res[1]: 
     1298            opts['res'] = res[1] 
     1299            data1, _ = make_data(opts) 
     1300        else: 
     1301            data1 = data0 
     1302        opts['res'] = res 
     1303        data = data0, data1 
     1304 
     1305    base = make_engine(model_info[0], data[0], opts['engine'][0], 
    12191306                       opts['cutoff'][0], opts['ngauss'][0]) 
    12201307    if comparison: 
    1221         comp = make_engine(model_info[1], data, opts['engine'][1], 
     1308        comp = make_engine(model_info[1], data[1], opts['engine'][1], 
    12221309                           opts['cutoff'][1], opts['ngauss'][1]) 
    12231310    else: 
  • sasmodels/core.py

    r3221de0 r2dcd6e7  
    3737if CUSTOM_MODEL_PATH == "": 
    3838    CUSTOM_MODEL_PATH = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models") 
    39     if not os.path.isdir(CUSTOM_MODEL_PATH): 
    40         os.makedirs(CUSTOM_MODEL_PATH) 
     39    #if not os.path.isdir(CUSTOM_MODEL_PATH): 
     40    #    os.makedirs(CUSTOM_MODEL_PATH) 
    4141 
    4242# TODO: refactor composite model support 
     
    233233        if not callable(model_info.Iq): 
    234234            source = generate.make_source(model_info)['dll'] 
    235             old_path = kerneldll.DLL_PATH 
     235            old_path = kerneldll.SAS_DLL_PATH 
    236236            try: 
    237                 kerneldll.DLL_PATH = path 
     237                kerneldll.SAS_DLL_PATH = path 
    238238                dll = kerneldll.make_dll(source, model_info, dtype=numpy_dtype) 
    239239            finally: 
    240                 kerneldll.DLL_PATH = old_path 
     240                kerneldll.SAS_DLL_PATH = old_path 
    241241            compiled_dlls.append(dll) 
    242242    return compiled_dlls 
  • sasmodels/data.py

    rd86f0fc rbd7630d  
    3636 
    3737import numpy as np  # type: ignore 
     38from numpy import sqrt, sin, cos, pi 
    3839 
    3940# pylint: disable=unused-import 
     
    182183    *x* is spin echo length and *y* is polarization (P/P0). 
    183184    """ 
     185    isSesans = True 
    184186    def __init__(self, **kw): 
    185187        Data1D.__init__(self, **kw) 
     
    300302        self.wavelength_unit = "A" 
    301303 
    302  
    303 def empty_data1D(q, resolution=0.0): 
     304class Sample(object): 
     305    """ 
     306    Sample attributes. 
     307    """ 
     308    def __init__(self): 
     309        # type: () -> None 
     310        pass 
     311 
     312def empty_data1D(q, resolution=0.0, L=0., dL=0.): 
    304313    # type: (np.ndarray, float) -> Data1D 
    305     """ 
     314    r""" 
    306315    Create empty 1D data using the given *q* as the x value. 
    307316 
    308     *resolution* dq/q defaults to 5%. 
     317    rms *resolution* $\Delta q/q$ defaults to 0%.  If wavelength *L* and rms 
     318    wavelength divergence *dL* are defined, then *resolution* defines 
     319    rms $\Delta \theta/\theta$ for the lowest *q*, with $\theta$ derived from 
     320    $q = 4\pi/\lambda \sin(\theta)$. 
    309321    """ 
    310322 
     
    313325    Iq, dIq = None, None 
    314326    q = np.asarray(q) 
    315     data = Data1D(q, Iq, dx=resolution * q, dy=dIq) 
     327    if L != 0 and resolution != 0: 
     328        theta = np.arcsin(q*L/(4*pi)) 
     329        dtheta = theta[0]*resolution 
     330        ## Solving Gaussian error propagation from 
     331        ##   Dq^2 = (dq/dL)^2 DL^2 + (dq/dtheta)^2 Dtheta^2 
     332        ## gives 
     333        ##   (Dq/q)^2 = (DL/L)**2 + (Dtheta/tan(theta))**2 
     334        ## Take the square root and multiply by q, giving 
     335        ##   Dq = (4*pi/L) * sqrt((sin(theta)*dL/L)**2 + (cos(theta)*dtheta)**2) 
     336        dq = (4*pi/L) * sqrt((sin(theta)*dL/L)**2 + (cos(theta)*dtheta)**2) 
     337    else: 
     338        dq = resolution * q 
     339 
     340    data = Data1D(q, Iq, dx=dq, dy=dIq) 
    316341    data.filename = "fake data" 
    317342    return data 
     
    486511            # Note: masks merge, so any masked theory points will stay masked, 
    487512            # and the data mask will be added to it. 
    488             mtheory = masked_array(theory, data.mask.copy()) 
     513            #mtheory = masked_array(theory, data.mask.copy()) 
     514            theory_x = data.x[data.mask == 0] 
     515            mtheory = masked_array(theory) 
    489516            mtheory[~np.isfinite(mtheory)] = masked 
    490517            if view is 'log': 
    491518                mtheory[mtheory <= 0] = masked 
    492             plt.plot(data.x, scale*mtheory, '-') 
     519            plt.plot(theory_x, scale*mtheory, '-') 
    493520            all_positive = all_positive and (mtheory > 0).all() 
    494521            some_present = some_present or (mtheory.count() > 0) 
     
    526553 
    527554    if use_resid: 
    528         mresid = masked_array(resid, data.mask.copy()) 
     555        theory_x = data.x[data.mask == 0] 
     556        mresid = masked_array(resid) 
    529557        mresid[~np.isfinite(mresid)] = masked 
    530558        some_present = (mresid.count() > 0) 
     
    532560        if num_plots > 1: 
    533561            plt.subplot(1, num_plots, use_calc + 2) 
    534         plt.plot(data.x, mresid, '.') 
     562        plt.plot(theory_x, mresid, '.') 
    535563        plt.xlabel("$q$/A$^{-1}$") 
    536564        plt.ylabel('residuals') 
  • sasmodels/direct_model.py

    rd18d6dd r7b9e4dd  
    3131from . import resolution2d 
    3232from .details import make_kernel_args, dispersion_mesh 
     33from .modelinfo import DEFAULT_BACKGROUND 
    3334 
    3435# pylint: disable=unused-import 
     
    250251            qmax = getattr(data, 'qmax', np.inf) 
    251252            accuracy = getattr(data, 'accuracy', 'Low') 
    252             index = ~data.mask & (q >= qmin) & (q <= qmax) 
     253            index = (data.mask == 0) & (q >= qmin) & (q <= qmax) 
    253254            if data.data is not None: 
    254255                index &= ~np.isnan(data.data) 
     
    263264        elif self.data_type == 'Iq': 
    264265            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     266            mask = getattr(data, 'mask', None) 
     267            if mask is not None: 
     268                index &= (mask == 0) 
    265269            if data.y is not None: 
    266270                index &= ~np.isnan(data.y) 
     
    346350 
    347351        # Need to pull background out of resolution for multiple scattering 
    348         background = pars.get('background', 0.) 
     352        background = pars.get('background', DEFAULT_BACKGROUND) 
    349353        pars = pars.copy() 
    350354        pars['background'] = 0. 
  • sasmodels/generate.py

    rd86f0fc r6e45516  
    965965    docs = model_info.docs if model_info.docs is not None else "" 
    966966    docs = convert_section_titles_to_boldface(docs) 
    967     pars = make_partable(model_info.parameters.COMMON 
    968                          + model_info.parameters.kernel_parameters) 
     967    if model_info.structure_factor: 
     968        pars = model_info.parameters.kernel_parameters 
     969    else: 
     970        pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters 
     971    partable = make_partable(pars) 
    969972    subst = dict(id=model_info.id.replace('_', '-'), 
    970973                 name=model_info.name, 
    971974                 title=model_info.title, 
    972                  parameters=pars, 
     975                 parameters=partable, 
    973976                 returns=Sq_units if model_info.structure_factor else Iq_units, 
    974977                 docs=docs) 
  • sasmodels/jitter.py

    rb3703f5 r1198f90  
    774774        # set small jitter as 0 if multiple pd dims 
    775775        dims = sum(v > 0 for v in jitter) 
    776         limit = [0, 0, 0.5, 5][dims] 
     776        limit = [0, 0.5, 5][dims] 
    777777        jitter = [0 if v < limit else v for v in jitter] 
    778778        axes.cla() 
  • sasmodels/kernel_iq.c

    rdc6f601 r70530778  
    8383  in_spin = clip(in_spin, 0.0, 1.0); 
    8484  out_spin = clip(out_spin, 0.0, 1.0); 
    85   // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that 
     85  // Previous version of this function took the square root of the weights, 
     86  // under the assumption that 
     87  // 
    8688  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
    87   // which is likely to be the case for simple models. 
    88   weight[0] = sqrt((1.0-in_spin) * (1.0-out_spin)); // dd 
    89   weight[1] = sqrt((1.0-in_spin) * out_spin);       // du.real 
    90   weight[2] = sqrt(in_spin * (1.0-out_spin));       // ud.real 
    91   weight[3] = sqrt(in_spin * out_spin);             // uu 
     89  // 
     90  // However, since the weights are applied to the final intensity and 
     91  // are not interned inside the I(q) function, we want the full 
     92  // weight and not the square root.  Any function using 
     93  // set_spin_weights as part of calculating an amplitude will need to 
     94  // manually take that square root, but there is currently no such 
     95  // function. 
     96  weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd 
     97  weight[1] = (1.0-in_spin) * out_spin;       // du 
     98  weight[2] = in_spin * (1.0-out_spin);       // ud 
     99  weight[3] = in_spin * out_spin;             // uu 
    92100  weight[4] = weight[1]; // du.imag 
    93101  weight[5] = weight[2]; // ud.imag 
     
    180188    QACRotation *rotation, 
    181189    double qx, double qy, 
    182     double *qa_out, double *qc_out) 
     190    double *qab_out, double *qc_out) 
    183191{ 
     192    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    184193    const double dqc = rotation->R31*qx + rotation->R32*qy; 
    185     // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    186     const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
    187  
    188     *qa_out = dqa; 
     194    const double dqab_sq = -dqc*dqc + qx*qx + qy*qy; 
     195    //*qab_out = sqrt(fabs(dqab_sq)); 
     196    *qab_out = dqab_sq > 0.0 ? sqrt(dqab_sq) : 0.0; 
    189197    *qc_out = dqc; 
    190198} 
  • sasmodels/kerneldll.py

    rbfe4245 rcdcd80f  
    9999    pass 
    100100# pylint: enable=unused-import 
     101 
     102if "SAS_DLL_PATH" in os.environ: 
     103    SAS_DLL_PATH = os.environ["SAS_DLL_PATH"] 
     104else: 
     105    # Assume the default location of module DLLs is in .sasmodels/compiled_models. 
     106    SAS_DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    101107 
    102108if "SAS_COMPILER" in os.environ: 
     
    166172        return CC + [source, "-o", output, "-lm"] 
    167173 
    168 # Assume the default location of module DLLs is in .sasmodels/compiled_models. 
    169 DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    170  
    171174ALLOW_SINGLE_PRECISION_DLLS = True 
    172175 
     
    205208        return path 
    206209 
    207     return joinpath(DLL_PATH, basename) 
     210    return joinpath(SAS_DLL_PATH, basename) 
    208211 
    209212 
     
    214217    exist yet if it hasn't been compiled. 
    215218    """ 
    216     return os.path.join(DLL_PATH, dll_name(model_info, dtype)) 
     219    return os.path.join(SAS_DLL_PATH, dll_name(model_info, dtype)) 
    217220 
    218221 
     
    233236    models are not allowed as DLLs. 
    234237 
    235     Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 
     238    Set *sasmodels.kerneldll.SAS_DLL_PATH* to the compiled dll output path. 
     239    Alternatively, set the environment variable *SAS_DLL_PATH*. 
    236240    The default is in ~/.sasmodels/compiled_models. 
    237241    """ 
     
    252256    if need_recompile: 
    253257        # Make sure the DLL path exists 
    254         if not os.path.exists(DLL_PATH): 
    255             os.makedirs(DLL_PATH) 
     258        if not os.path.exists(SAS_DLL_PATH): 
     259            os.makedirs(SAS_DLL_PATH) 
    256260        basename = splitext(os.path.basename(dll))[0] + "_" 
    257261        system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 
  • sasmodels/model_test.py

    • Property mode changed from 100644 to 100755
    r3221de0 r012cd34  
    376376        stream.writeln(traceback.format_exc()) 
    377377        return 
    378     # Run the test suite 
    379     suite.run(result) 
    380  
    381     # Print the failures and errors 
    382     for _, tb in result.errors: 
    383         stream.writeln(tb) 
    384     for _, tb in result.failures: 
    385         stream.writeln(tb) 
    386378 
    387379    # Warn if there are no user defined tests. 
     
    393385    # iterator since we don't have direct access to the list of tests in the 
    394386    # test suite. 
     387    # In Qt5 suite.run() will clear all tests in the suite after running 
     388    # with no way of retaining them for the test below, so let's check 
     389    # for user tests before running the suite. 
    395390    for test in suite: 
    396391        if not test.info.tests: 
     
    399394    else: 
    400395        stream.writeln("Note: no test suite created --- this should never happen") 
     396 
     397    # Run the test suite 
     398    suite.run(result) 
     399 
     400    # Print the failures and errors 
     401    for _, tb in result.errors: 
     402        stream.writeln(tb) 
     403    for _, tb in result.failures: 
     404        stream.writeln(tb) 
    401405 
    402406    output = stream.getvalue() 
  • sasmodels/modelinfo.py

    r95498a3 r7b9e4dd  
    4545# Note that scale and background cannot be coordinated parameters whose value 
    4646# depends on the some polydisperse parameter with the current implementation 
     47DEFAULT_BACKGROUND = 1e-3 
    4748COMMON_PARAMETERS = [ 
    4849    ("scale", "", 1, (0.0, np.inf), "", "Source intensity"), 
    49     ("background", "1/cm", 1e-3, (-np.inf, np.inf), "", "Source background"), 
     50    ("background", "1/cm", DEFAULT_BACKGROUND, (-np.inf, np.inf), "", "Source background"), 
    5051] 
    5152assert (len(COMMON_PARAMETERS) == 2 
     
    589590                Parameter('up:frac_f', '', 0., [0., 1.], 
    590591                          'magnetic', 'fraction of spin up final'), 
    591                 Parameter('up:angle', 'degress', 0., [0., 360.], 
     592                Parameter('up:angle', 'degrees', 0., [0., 360.], 
    592593                          'magnetic', 'spin up angle'), 
    593594            ]) 
  • sasmodels/models/_spherepy.py

    r108e70e rca4444f  
    11r""" 
    22For information about polarised and magnetic scattering, see 
    3 the :doc:`magnetic help <../sasgui/perspectives/fitting/mag_help>` documentation. 
     3the :ref:`magnetism` documentation. 
    44 
    55Definition 
  • sasmodels/models/core_shell_cylinder.py

    r2d81cfe re31b19a  
    55The output of the 2D scattering intensity function for oriented core-shell 
    66cylinders is given by (Kline, 2006 [#kline]_). The form factor is normalized 
    7 by the particle volume. 
     7by the particle volume. Note that in this model the shell envelops the entire 
     8core so that besides a "sleeve" around the core, the shell also provides two 
     9flat end caps of thickness = shell thickness. In other words the length of the 
     10total cyclinder is the length of the core cylinder plus twice the thickness of 
     11the shell. If no end caps are desired one should use the 
     12:ref:`core-shell-bicelle` and set the thickness of the end caps (in this case 
     13the "thick_face") to zero. 
    814 
    915.. math:: 
     
    3339 
    3440and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, 
    35 $V_s$ is the volume of the outer shell (i.e. the total volume, including 
    36 the shell), $V_c$ is the volume of the core, $L$ is the length of the core, 
     41$V_s$ is the total volume (i.e. including both the core and the outer shell), 
     42$V_c$ is the volume of the core, $L$ is the length of the core, 
    3743$R$ is the radius of the core, $T$ is the thickness of the shell, $\rho_c$ 
    3844is the scattering length density of the core, $\rho_s$ is the scattering 
     
    135141    return 0.5 * (ddd) ** (1. / 3.) 
    136142 
    137 def VR(radius, thickness, length): 
    138     """ 
    139     Returns volume ratio 
    140     """ 
    141     whole = pi * (radius + thickness) ** 2 * (length + 2 * thickness) 
    142     core = pi * radius ** 2 * length 
    143     return whole, whole - core 
    144  
    145143def random(): 
    146144    outer_radius = 10**np.random.uniform(1, 4.7) 
  • sasmodels/models/core_shell_parallelepiped.c

    re077231 rdbf1a60  
    5959 
    6060    // outer integral (with gauss points), integration limits = 0, 1 
     61    // substitute d_cos_alpha for sin_alpha d_alpha 
    6162    double outer_sum = 0; //initialize integral 
    6263    for( int i=0; i<GAUSS_N; i++) { 
    6364        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    6465        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    65  
    66         // inner integral (with gauss points), integration limits = 0, pi/2 
    6766        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
    6867        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
     68 
     69        // inner integral (with gauss points), integration limits = 0, 1 
     70        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta) 
    6971        double inner_sum = 0.0; 
    7072        for(int j=0; j<GAUSS_N; j++) { 
    71             const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     73            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    7274            double sin_beta, cos_beta; 
    73             SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
     75            SINCOS(M_PI_2*u, sin_beta, cos_beta); 
    7476            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
    7577            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
     
    9193            inner_sum += GAUSS_W[j] * f * f; 
    9294        } 
     95        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    9396        inner_sum *= 0.5; 
    9497        // now sum up the outer integral 
    9598        outer_sum += GAUSS_W[i] * inner_sum; 
    9699    } 
     100    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    97101    outer_sum *= 0.5; 
    98102 
  • sasmodels/models/core_shell_parallelepiped.py

    r97be877 rf89ec96  
    44 
    55Calculates the form factor for a rectangular solid with a core-shell structure. 
    6 The thickness and the scattering length density of the shell or 
    7 "rim" can be different on each (pair) of faces. 
     6The thickness and the scattering length density of the shell or "rim" can be 
     7different on each (pair) of faces. The three dimensions of the core of the 
     8parallelepiped (strictly here a cuboid) may be given in *any* size order as 
     9long as the particles are randomly oriented (i.e. take on all possible 
     10orientations see notes on 2D below). To avoid multiple fit solutions, 
     11especially with Monte-Carlo fit methods, it may be advisable to restrict their 
     12ranges. There may be a number of closely similar "best fits", so some trial and 
     13error, or fixing of some dimensions at expected values, may help. 
    814 
    915The form factor is normalized by the particle volume $V$ such that 
     
    1117.. math:: 
    1218 
    13     I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background} 
     19    I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle 
     20    + \text{background} 
    1421 
    1522where $\langle \ldots \rangle$ is an average over all possible orientations 
    16 of the rectangular solid. 
    17  
    18 The function calculated is the form factor of the rectangular solid below. 
    19 The core of the solid is defined by the dimensions $A$, $B$, $C$ such that 
    20 $A < B < C$. 
    21  
    22 .. image:: img/core_shell_parallelepiped_geometry.jpg 
     23of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be 
     24pulled out of the form factor term due to the multiple slds in the model. 
     25 
     26The core of the solid is defined by the dimensions $A$, $B$, $C$ here shown 
     27such that $A < B < C$. 
     28 
     29.. figure:: img/parallelepiped_geometry.jpg 
     30 
     31   Core of the core shell parallelepiped with the corresponding definition 
     32   of sides. 
     33 
    2334 
    2435There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 
    2536(on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$ 
    26 $(=t_C)$ faces. The projection in the $AB$ plane is then 
    27  
    28 .. image:: img/core_shell_parallelepiped_projection.jpg 
    29  
    30 The volume of the solid is 
     37$(=t_C)$ faces. The projection in the $AB$ plane is 
     38 
     39.. figure:: img/core_shell_parallelepiped_projection.jpg 
     40 
     41   AB cut through the core-shell parallelipiped showing the cross secion of 
     42   four of the six shell slabs. As can be seen, this model leaves **"gaps"** 
     43   at the corners of the solid. 
     44 
     45 
     46The total volume of the solid is thus given as 
    3147 
    3248.. math:: 
    3349 
    3450    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    35  
    36 **meaning that there are "gaps" at the corners of the solid.** 
    3751 
    3852The intensity calculated follows the :ref:`parallelepiped` model, with the 
    3953core-shell intensity being calculated as the square of the sum of the 
    40 amplitudes of the core and the slabs on the edges. 
    41  
    42 the scattering amplitude is computed for a particular orientation of the 
    43 core-shell parallelepiped with respect to the scattering vector and then 
    44 averaged over all possible orientations, where $\alpha$ is the angle between 
    45 the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
    46 the angle between projection of the particle in the $xy$ detector plane 
    47 and the $y$ axis. 
    48  
    49 .. math:: 
    50  
    51     F(Q) 
     54amplitudes of the core and the slabs on the edges. The scattering amplitude is 
     55computed for a particular orientation of the core-shell parallelepiped with 
     56respect to the scattering vector and then averaged over all possible 
     57orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis 
     58of the parallelepiped, and $\beta$ is the angle between the projection of the 
     59particle in the $xy$ detector plane and the $y$ axis. 
     60 
     61.. math:: 
     62 
     63    P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha 
     64    \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta} 
     65 
     66and 
     67 
     68.. math:: 
     69 
     70    F(q,\alpha,\beta) 
    5271    &= (\rho_\text{core}-\rho_\text{solvent}) 
    5372       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
    5473    &+ (\rho_\text{A}-\rho_\text{solvent}) 
    55         \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
     74        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\ 
    5675    &+ (\rho_\text{B}-\rho_\text{solvent}) 
    5776        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
     
    6382.. math:: 
    6483 
    65     S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
     84    S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L} 
    6685 
    6786and 
     
    6988.. math:: 
    7089 
    71     Q_A &= \sin\alpha \sin\beta \\ 
    72     Q_B &= \sin\alpha \cos\beta \\ 
    73     Q_C &= \cos\alpha 
     90    Q_A &= q \sin\alpha \sin\beta \\ 
     91    Q_B &= q \sin\alpha \cos\beta \\ 
     92    Q_C &= q \cos\alpha 
    7493 
    7594 
    7695where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
    77 are the scattering length of the parallelepiped core, and the rectangular 
     96are the scattering lengths of the parallelepiped core, and the rectangular 
    7897slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
    7998is the scattering length of the solvent. 
    8099 
     100.. note::  
     101 
     102   the code actually implements two substitutions: $d(cos\alpha)$ is 
     103   substituted for -$sin\alpha \ d\alpha$ (note that in the 
     104   :ref:`parallelepiped` code this is explicitly implemented with 
     105   $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that 
     106   $du = \pi/2 \ d\beta$.  Thus both integrals go from 0 to 1 rather than 0 
     107   to $\pi/2$. 
     108 
    81109FITTING NOTES 
    82110~~~~~~~~~~~~~ 
    83111 
    84 If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    85 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 
    86 **no interparticle interference effects are included in this calculation.** 
    87  
    88 There are many parameters in this model. Hold as many fixed as possible with 
    89 known values, or you will certainly end up at a solution that is unphysical. 
    90  
    91 The returned value is in units of |cm^-1|, on absolute scale. 
    92  
    93 NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    94 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    95 and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
    96 to give an oblate or prolate particle, to give an effective radius, 
    97 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    98  
    99 For 2d data the orientation of the particle is required, described using 
    100 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
    101 details of the calculation and angular dispersions see :ref:`orientation`. 
    102 The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    103 $\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    104  
    105 For 2d, constraints must be applied during fitting to ensure that the 
    106 inequality $A < B < C$ is not violated, and hence the correct definition 
    107 of angles is preserved. The calculation will not report an error, 
    108 but the results may be not correct. 
     112#. There are many parameters in this model. Hold as many fixed as possible with 
     113   known values, or you will certainly end up at a solution that is unphysical. 
     114 
     115#. The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
     116   based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
     117   and length $(C+2t_C)$ values, after appropriately sorting the three 
     118   dimensions to give an oblate or prolate particle, to give an effective radius 
     119   for $S(q)$ when $P(q) * S(q)$ is applied. 
     120 
     121#. For 2d data the orientation of the particle is required, described using 
     122   angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, where $\theta$ 
     123   and $\phi$ define the orientation of the director in the laboratry reference 
     124   frame of the beam direction ($z$) and detector plane ($x-y$ plane), while 
     125   the angle $\Psi$ is effectively the rotational angle around the particle 
     126   $C$ axis. For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the 
     127   $B$ axis oriented parallel to the y-axis of the detector with $A$ along 
     128   the x-axis. For other $\theta$, $\phi$ values, the order of rotations 
     129   matters. In particular, the parallelepiped must first be rotated $\theta$ 
     130   degrees in the $x-z$ plane before rotating $\phi$ degrees around the $z$ 
     131   axis (in the $x-y$ plane). Applying orientational distribution to the 
     132   particle orientation (i.e  `jitter` to one or more of these angles) can get 
     133   more confusing as `jitter` is defined **NOT** with respect to the laboratory 
     134   frame but the particle reference frame. It is thus highly recmmended to 
     135   read :ref:`orientation` for further details of the calculation and angular 
     136   dispersions. 
     137 
     138.. note:: For 2d, constraints must be applied during fitting to ensure that the 
     139   order of sides chosen is not altered, and hence that the correct definition 
     140   of angles is preserved. For the default choice shown here, that means 
     141   ensuring that the inequality $A < B < C$ is not violated,  The calculation 
     142   will not report an error, but the results may be not correct. 
    109143 
    110144.. figure:: img/parallelepiped_angle_definition.png 
    111145 
    112146    Definition of the angles for oriented core-shell parallelepipeds. 
    113     Note that rotation $\theta$, initially in the $xz$ plane, is carried 
     147    Note that rotation $\theta$, initially in the $x-z$ plane, is carried 
    114148    out first, then rotation $\phi$ about the $z$ axis, finally rotation 
    115     $\Psi$ is now around the axis of the cylinder. The neutron or X-ray 
    116     beam is along the $z$ axis. 
     149    $\Psi$ is now around the $C$ axis of the particle. The neutron or X-ray 
     150    beam is along the $z$ axis and the detecotr defines the $x-y$ plane. 
    117151 
    118152.. figure:: img/parallelepiped_angle_projection.png 
     
    120154    Examples of the angles for oriented core-shell parallelepipeds against the 
    121155    detector plane. 
     156 
     157 
     158Validation 
     159---------- 
     160 
     161Cross-checked against hollow rectangular prism and rectangular prism for equal 
     162thickness overlapping sides, and by Monte Carlo sampling of points within the 
     163shape for non-uniform, non-overlapping sides. 
     164 
    122165 
    123166References 
     
    135178 
    136179* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    137 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
     180* **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016 
    138181* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 
    139 * Cross-checked against hollow rectangular prism and rectangular prism for 
    140   equal thickness overlapping sides, and by Monte Carlo sampling of points 
    141   within the shape for non-uniform, non-overlapping sides. 
     182* **Last Reviewed by:** Paul Butler **Date:** May 24, 2018 - documentation 
     183  updated 
    142184""" 
    143185 
  • sasmodels/models/core_shell_sphere.py

    r2d81cfe rda1c8d1  
    2121.. math:: 
    2222 
    23     F^2(q) = \frac{3}{V_s}\left[ 
     23    F(q) = \frac{3}{V_s}\left[ 
    2424       V_c(\rho_c-\rho_s)\frac{\sin(qr_c)-qr_c\cos(qr_c)}{(qr_c)^3} + 
    2525       V_s(\rho_s-\rho_\text{solv})\frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3} 
     
    8989    return radius + thickness 
    9090 
    91 def VR(radius, thickness): 
    92     """ 
    93         Volume ratio 
    94         @param radius: core radius 
    95         @param thickness: shell thickness 
    96     """ 
    97     return (1, 1) 
    98     whole = 4.0/3.0 * pi * (radius + thickness)**3 
    99     core = 4.0/3.0 * pi * radius**3 
    100     return whole, whole - core 
    101  
    10291def random(): 
    10392    outer_radius = 10**np.random.uniform(1.3, 4.3) 
     
    114103tests = [ 
    115104    [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
    116     # TODO: VR test suppressed until we sort out new product model 
    117     # and determine what to do with volume ratio. 
    118     #[{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704], 
    119105 
    120106    # The SasView test result was 0.00169, with a background of 0.001 
  • sasmodels/models/ellipsoid.py

    r2d81cfe r0168844  
    125125import numpy as np 
    126126from numpy import inf, sin, cos, pi 
     127 
     128try: 
     129    from numpy import cbrt 
     130except ImportError: 
     131    def cbrt(x): return x ** (1.0/3.0) 
    127132 
    128133name = "ellipsoid" 
     
    170175    idx = radius_polar < radius_equatorial 
    171176    ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2 
    172     idx = radius_polar == radius_equatorial 
    173     ee[idx] = 2 * radius_polar[idx] 
    174     valid = (radius_polar * radius_equatorial != 0) 
     177    valid = (radius_polar * radius_equatorial != 0) & (radius_polar != radius_equatorial) 
    175178    bd = 1.0 - ee[valid] 
    176179    e1 = np.sqrt(ee[valid]) 
     
    179182    b2 = 1.0 + bd / 2 / e1 * np.log(bL) 
    180183    delta = 0.75 * b1 * b2 
    181  
    182     ddd = np.zeros_like(radius_polar) 
    183     ddd[valid] = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial ** 2 
    184     return 0.5 * ddd ** (1.0 / 3.0) 
     184    ddd = 2.0 * (delta + 1.0) * (radius_polar * radius_equatorial**2)[valid] 
     185 
     186    r = np.zeros_like(radius_polar) 
     187    r[valid] = 0.5 * cbrt(ddd) 
     188    idx = radius_polar == radius_equatorial 
     189    r[idx] = radius_polar[idx] 
     190    return r 
    185191 
    186192def random(): 
  • sasmodels/models/fractal_core_shell.py

    ref07e95 reb3eb38  
    8888    ["sld_shell",   "1e-6/Ang^2", 2.0,  [-inf, inf], "sld",    "Sphere shell scattering length density"], 
    8989    ["sld_solvent", "1e-6/Ang^2", 3.0,  [-inf, inf], "sld",    "Solvent scattering length density"], 
    90     ["volfraction", "",           1.0,  [0.0, inf],  "",       "Volume fraction of building block spheres"], 
     90    ["volfraction", "",           0.05,  [0.0, inf],  "",       "Volume fraction of building block spheres"], 
    9191    ["fractal_dim",    "",        2.0,  [0.0, 6.0],  "",       "Fractal dimension"], 
    9292    ["cor_length",  "Ang",      100.0,  [0.0, inf],  "",       "Correlation length of fractal-like aggregates"], 
     
    134134    return radius + thickness 
    135135 
    136 def VR(radius, thickness): 
    137     """ 
    138         Volume ratio 
    139         @param radius: core radius 
    140         @param thickness: shell thickness 
    141     """ 
    142     whole = 4.0/3.0 * pi * (radius + thickness)**3 
    143     core = 4.0/3.0 * pi * radius**3 
    144     return whole, whole-core 
     136tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
    145137 
    146 tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
    147          [{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704]] 
    148  
    149 #         # The SasView test result was 0.00169, with a background of 0.001 
    150 #         # They are however wrong as we now know.  IGOR might be a more 
    151 #         # appropriate source. Otherwise will just have to assume this is now 
    152 #         # correct and self generate a correct answer for the future. Until we 
    153 #         # figure it out leave the tests commented out 
    154 #         [{'radius': 60.0, 
    155 #           'thickness': 10.0, 
    156 #           'sld_core': 1.0, 
    157 #           'sld_shell': 2.0, 
    158 #           'sld_solvent': 3.0, 
    159 #           'background': 0.0 
    160 #          }, 0.015211, 692.84]] 
     138#         # At some point the SasView 3.x test result was deemed incorrect. The 
     139          #following tests were verified against NIST IGOR macros ver 7.850. 
     140          #NOTE: NIST macros do only provide for a polydispers core (no option 
     141          #for a poly shell or for a monodisperse core.  The results seemed 
     142          #extremely sensitive to the core PD, varying non monotonically all 
     143          #the way to a PD of 1e-6. From 1e-6 to 1e-9 no changes in the 
     144          #results were observed and the values below were taken using PD=1e-9. 
     145          #Non-monotonically = I(0.001)=188 to 140 to 177 back to 160 etc. 
     146         [{'radius': 20.0, 
     147           'thickness': 5.0, 
     148           'sld_core': 3.5, 
     149           'sld_shell': 1.0, 
     150           'sld_solvent': 6.35, 
     151           'volfraction': 0.05, 
     152           'background': 0.0}, 
     153           [0.001,0.00291,0.0107944,0.029923,0.100726,0.476304], 
     154           [177.146,165.151,84.1596,20.1466,1.40906,0.00622666]]] 
  • sasmodels/models/guinier.py

    r2d81cfe rc9fc873  
    77.. math:: 
    88 
    9     I(q) = \text{scale} \cdot \exp{\left[ \frac{-Q^2R_g^2}{3} \right]} 
     9    I(q) = \text{scale} \cdot \exp{\left[ \frac{-Q^2 R_g^2 }{3} \right]} 
    1010            + \text{background} 
    1111 
     
    1919 
    2020.. math:: q = \sqrt{q_x^2 + q_y^2} 
     21 
     22In scattering, the radius of gyration $R_g$ quantifies the objects's 
     23distribution of SLD (not mass density, as in mechanics) from the objects's 
     24SLD centre of mass. It is defined by 
     25 
     26.. math:: R_g^2 = \frac{\sum_i\rho_i\left(r_i-r_0\right)^2}{\sum_i\rho_i} 
     27 
     28where $r_0$ denotes the object's SLD centre of mass and $\rho_i$ is the SLD at 
     29a point $i$. 
     30 
     31Notice that $R_g^2$ may be negative (since SLD can be negative), which happens 
     32when a form factor $P(Q)$ is increasing with $Q$ rather than decreasing. This 
     33can occur for core/shell particles, hollow particles, or for composite 
     34particles with domains of different SLDs in a solvent with an SLD close to the 
     35average match point. (Alternatively, this might be regarded as there being an 
     36internal inter-domain "structure factor" within a single particle which gives 
     37rise to a peak in the scattering). 
     38 
     39To specify a negative value of $R_g^2$ in SasView, simply give $R_g$ a negative 
     40value ($R_g^2$ will be evaluated as $R_g |R_g|$). Note that the physical radius  
     41of gyration, of the exterior of the particle, will still be large and positive.  
     42It is only the apparent size from the small $Q$ data that will give a small or  
     43negative value of $R_g^2$. 
    2144 
    2245References 
     
    4265 
    4366#             ["name", "units", default, [lower, upper], "type","description"], 
    44 parameters = [["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"]] 
     67parameters = [["rg", "Ang", 60.0, [-inf, inf], "", "Radius of Gyration"]] 
    4568 
    4669Iq = """ 
    47     double exponent = rg*rg*q*q/3.0; 
     70    double exponent = fabs(rg)*rg*q*q/3.0; 
    4871    double value = exp(-exponent); 
    4972    return value; 
     
    6689 
    6790# parameters for demo 
    68 demo = dict(scale=1.0, rg=60.0) 
     91demo = dict(scale=1.0,  background=0.001, rg=60.0 ) 
    6992 
    7093# parameters for unit tests 
  • sasmodels/models/hollow_cylinder.py

    r2d81cfe r455aaa1  
    11r""" 
     2Definition 
     3---------- 
     4 
    25This model provides the form factor, $P(q)$, for a monodisperse hollow right 
    3 angle circular cylinder (rigid tube) where the form factor is normalized by the 
    4 volume of the tube (i.e. not by the external volume). 
     6angle circular cylinder (rigid tube) where the The inside and outside of the 
     7hollow cylinder are assumed to have the same SLD and the form factor is thus 
     8normalized by the volume of the tube (i.e. not by the total cylinder volume). 
    59 
    610.. math:: 
     
    812    P(q) = \text{scale} \left<F^2\right>/V_\text{shell} + \text{background} 
    913 
    10 where the averaging $\left<\ldots\right>$ is applied only for the 1D calculation. 
     14where the averaging $\left<\ldots\right>$ is applied only for the 1D 
     15calculation. If Intensity is given on an absolute scale, the scale factor here 
     16is the volume fraction of the shell.  This differs from 
     17the :ref:`core-shell-cylinder` in that, in that case, scale is the volume 
     18fraction of the entire cylinder (core+shell). The application might be for a 
     19bilayer which wraps into a hollow tube and the volume fraction of material is 
     20all in the shell, whereas the :ref:`core-shell-cylinder` model might be used for 
     21a cylindrical micelle where the tails in the core have a different SLD than the 
     22headgroups (in the shell) and the volume fraction of material comes fromm the 
     23whole cyclinder.  NOTE: the hollow_cylinder represents a tube whereas the 
     24core_shell_cylinder includes a shell layer covering the ends (end caps) as well. 
    1125 
    12 The inside and outside of the hollow cylinder are assumed have the same SLD. 
    13  
    14 Definition 
    15 ---------- 
    1626 
    1727The 1D scattering intensity is calculated in the following way (Guinier, 1955) 
     
    4858---------- 
    4959 
    50 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 
    51 Neutron Scattering*, Plenum Press, New York, (1987) 
     60.. [#] L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 
     61   Neutron Scattering*, Plenum Press, New York, (1987) 
    5262 
    5363Authorship and Verification 
     
    5565 
    5666* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    57 * **Last Modified by:** Richard Heenan **Date:** October 06, 2016 
    58    (reparametrised to use thickness, not outer radius) 
    59 * **Last Reviewed by:** Richard Heenan **Date:** October 06, 2016 
     67* **Last Modified by:** Paul Butler **Date:** September 06, 2018 
     68   (corrected VR calculation) 
     69* **Last Reviewed by:** Paul Butler **Date:** September 06, 2018 
    6070""" 
    6171 
     
    120130    vol_total = pi*router*router*length 
    121131    vol_shell = vol_total - vol_core 
    122     return vol_shell, vol_total 
     132    return vol_total, vol_shell 
    123133 
    124134def random(): 
     
    151161tests = [ 
    152162    [{}, 0.00005, 1764.926], 
    153     [{}, 'VR', 1.8], 
     163    [{}, 'VR', 0.55555556], 
    154164    [{}, 0.001, 1756.76], 
    155165    [{}, (qx, qy), 2.36885476192], 
  • sasmodels/models/hollow_rectangular_prism.py

    r0e55afe r455aaa1  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    4  
    5 This model provides the form factor, $P(q)$, for a hollow rectangular 
    6 parallelepiped with a wall of thickness $\Delta$. 
    7  
    8  
    94Definition 
    105---------- 
    116 
    12 The 1D scattering intensity for this model is calculated by forming 
    13 the difference of the amplitudes of two massive parallelepipeds 
    14 differing in their outermost dimensions in each direction by the 
    15 same length increment $2\Delta$ (Nayuk, 2012). 
     7This model provides the form factor, $P(q)$, for a hollow rectangular 
     8parallelepiped with a wall of thickness $\Delta$. The 1D scattering intensity 
     9for this model is calculated by forming the difference of the amplitudes of two 
     10massive parallelepipeds differing in their outermost dimensions in each 
     11direction by the same length increment $2\Delta$ (\ [#Nayuk2012]_ Nayuk, 2012). 
    1612 
    1713As in the case of the massive parallelepiped model (:ref:`rectangular-prism`), 
     
    6157  \rho_\text{solvent})^2 \times P(q) + \text{background} 
    6258 
    63 where $\rho_\text{p}$ is the scattering length of the parallelepiped, 
    64 $\rho_\text{solvent}$ is the scattering length of the solvent, 
     59where $\rho_\text{p}$ is the scattering length density of the parallelepiped, 
     60$\rho_\text{solvent}$ is the scattering length density of the solvent, 
    6561and (if the data are in absolute units) *scale* represents the volume fraction 
    66 (which is unitless). 
     62(which is unitless) of the rectangular shell of material (i.e. not including 
     63the volume of the solvent filled core). 
    6764 
    6865For 2d data the orientation of the particle is required, described using 
     
    7370 
    7471For 2d, constraints must be applied during fitting to ensure that the inequality 
    75 $A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, 
    76 but the results may be not correct. 
     72$A < B < C$ is not violated, and hence the correct definition of angles is 
     73preserved. The calculation will not report an error if the inequality is *not* 
     74preserved, but the results may be not correct. 
    7775 
    7876.. figure:: img/parallelepiped_angle_definition.png 
     
    9997---------- 
    10098 
    101 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     99.. [#Nayuk2012] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     100 
     101 
     102Authorship and Verification 
     103---------------------------- 
     104 
     105* **Author:** Miguel Gonzales **Date:** February 26, 2016 
     106* **Last Modified by:** Paul Kienzle **Date:** December 14, 2017 
     107* **Last Reviewed by:** Paul Butler **Date:** September 06, 2018 
    102108""" 
    103109 
  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    r2d81cfe r6e7d7b6  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
     4Definition 
     5---------- 
     6 
    47 
    58This model provides the form factor, $P(q)$, for a hollow rectangular 
    69prism with infinitely thin walls. It computes only the 1D scattering, not the 2D. 
    7  
    8  
    9 Definition 
    10 ---------- 
    11  
    1210The 1D scattering intensity for this model is calculated according to the 
    13 equations given by Nayuk and Huber (Nayuk, 2012). 
     11equations given by Nayuk and Huber\ [#Nayuk2012]_. 
    1412 
    1513Assuming a hollow parallelepiped with infinitely thin walls, edge lengths 
     
    5553  I(q) = \text{scale} \times V \times (\rho_\text{p} - \rho_\text{solvent})^2 \times P(q) 
    5654 
    57 where $V$ is the volume of the rectangular prism, $\rho_\text{p}$ 
    58 is the scattering length of the parallelepiped, $\rho_\text{solvent}$ 
    59 is the scattering length of the solvent, and (if the data are in absolute 
    60 units) *scale* represents the volume fraction (which is unitless). 
     55where $V$ is the surface area of the rectangular prism, $\rho_\text{p}$ 
     56is the scattering length density of the parallelepiped, $\rho_\text{solvent}$ 
     57is the scattering length density of the solvent, and (if the data are in 
     58absolute units) *scale* is related to the total surface area. 
    6159 
    6260**The 2D scattering intensity is not computed by this model.** 
     
    6765 
    6866Validation of the code was conducted  by qualitatively comparing the output 
    69 of the 1D model to the curves shown in (Nayuk, 2012). 
     67of the 1D model to the curves shown in (Nayuk, 2012\ [#Nayuk2012]_). 
    7068 
    7169 
     
    7371---------- 
    7472 
    75 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     73.. [#Nayuk2012] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     74 
     75 
     76Authorship and Verification 
     77---------------------------- 
     78 
     79* **Author:** Miguel Gonzales **Date:** February 26, 2016 
     80* **Last Modified by:** Paul Kienzle **Date:** October 15, 2016 
     81* **Last Reviewed by:** Paul Butler **Date:** September 07, 2018 
    7682""" 
    7783 
  • sasmodels/models/mass_surface_fractal.py

    r2d81cfe r7994359  
    3939    The surface ( $D_s$ ) and mass ( $D_m$ ) fractal dimensions are only 
    4040    valid if $0 < surface\_dim < 6$ , $0 < mass\_dim < 6$ , and 
    41     $(surface\_dim + mass\_dim ) < 6$ . 
    42  
     41    $(surface\_dim + mass\_dim ) < 6$ .  
     42    Older versions of sasview may have the default primary particle radius 
     43    larger than the cluster radius, this was an error, also present in the  
     44    Schmidt review paper below. The primary particle should be the smaller  
     45    as described in the original Hurd et.al. who also point out that  
     46    polydispersity in the primary particle sizes may affect their  
     47    apparent surface fractal dimension. 
     48     
    4349 
    4450References 
    4551---------- 
    4652 
    47 P Schmidt, *J Appl. Cryst.*, 24 (1991) 414-435 Equation(19) 
     53.. [#] P Schmidt, *J Appl. Cryst.*, 24 (1991) 414-435 Equation(19) 
     54.. [#] A J Hurd, D W Schaefer, J E Martin, *Phys. Rev. A*, 
     55   35 (1987) 2361-2364 Equation(2) 
    4856 
    49 A J Hurd, D W Schaefer, J E Martin, *Phys. Rev. A*, 
    50 35 (1987) 2361-2364 Equation(2) 
     57Authorship and Verification 
     58---------------------------- 
     59 
     60* **Converted to sasmodels by:** Piotr Rozyczko **Date:** Jan 20, 2016 
     61* **Last Reviewed by:** Richard Heenan **Date:** May 30, 2018 
    5162""" 
    5263 
     
    6778        rg_primary    =  rg 
    6879        background   =  background 
    69         Ref: Schmidt, J Appl Cryst, eq(19), (1991), 24, 414-435 
    7080        Hurd, Schaefer, Martin, Phys Rev A, eq(2),(1987),35, 2361-2364 
    7181        Note that 0 < Ds< 6 and 0 < Dm < 6. 
     
    7888    ["fractal_dim_mass", "",      1.8, [0.0, 6.0], "", "Mass fractal dimension"], 
    7989    ["fractal_dim_surf", "",      2.3, [0.0, 6.0], "", "Surface fractal dimension"], 
    80     ["rg_cluster",       "Ang",  86.7, [0.0, inf], "", "Cluster radius of gyration"], 
    81     ["rg_primary",       "Ang", 4000., [0.0, inf], "", "Primary particle radius of gyration"], 
     90    ["rg_cluster",       "Ang", 4000., [0.0, inf], "", "Cluster radius of gyration"], 
     91    ["rg_primary",       "Ang",  86.7, [0.0, inf], "", "Primary particle radius of gyration"], 
    8292] 
    8393# pylint: enable=bad-whitespace, line-too-long 
     
    107117            fractal_dim_mass=1.8, 
    108118            fractal_dim_surf=2.3, 
    109             rg_cluster=86.7, 
    110             rg_primary=4000.0) 
     119            rg_cluster=4000.0, 
     120            rg_primary=86.7) 
    111121 
    112122tests = [ 
    113123 
    114     # Accuracy tests based on content in test/utest_other_models.py 
    115     [{'fractal_dim_mass':      1.8, 
     124    # Accuracy tests based on content in test/utest_other_models.py  All except first, changed so rg_cluster is the larger, RKH 30 May 2018 
     125    [{'fractal_dim_mass':   1.8, 
    116126      'fractal_dim_surf':   2.3, 
    117127      'rg_cluster':   86.7, 
     
    123133    [{'fractal_dim_mass':      3.3, 
    124134      'fractal_dim_surf':   1.0, 
    125       'rg_cluster':   90.0, 
    126       'rg_primary': 4000.0, 
    127      }, 0.001, 0.18562699016], 
     135      'rg_cluster': 4000.0, 
     136      'rg_primary':   90.0, 
     137     }, 0.001, 0.0932516614456], 
    128138 
    129139    [{'fractal_dim_mass':      1.3, 
    130       'fractal_dim_surf':   1.0, 
    131       'rg_cluster':   90.0, 
    132       'rg_primary': 2000.0, 
     140      'fractal_dim_surf':   2.0, 
     141      'rg_cluster': 2000.0, 
     142      'rg_primary':   90.0, 
    133143      'background':    0.8, 
    134      }, 0.001, 1.16539753641], 
     144     }, 0.001, 1.28296431786], 
    135145 
    136146    [{'fractal_dim_mass':      2.3, 
    137       'fractal_dim_surf':   1.0, 
    138       'rg_cluster':   90.0, 
    139       'rg_primary': 1000.0, 
     147      'fractal_dim_surf':   3.1, 
     148      'rg_cluster':  1000.0, 
     149      'rg_primary':  30.0, 
    140150      'scale':        10.0, 
    141151      'background':    0.0, 
    142      }, 0.051, 0.000169548800377], 
     152     }, 0.051, 0.00333804044899], 
    143153    ] 
  • sasmodels/models/parallelepiped.c

    r108e70e rdbf1a60  
    3838            inner_total += GAUSS_W[j] * square(si1 * si2); 
    3939        } 
     40        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    4041        inner_total *= 0.5; 
    4142 
     
    4344        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4445    } 
     46    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    4547    outer_total *= 0.5; 
    4648 
  • sasmodels/models/parallelepiped.py

    ref07e95 rf89ec96  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    4 The form factor is normalized by the particle volume. 
    5 For information about polarised and magnetic scattering, see 
    6 the :ref:`magnetism` documentation. 
    7  
    84Definition 
    95---------- 
    106 
    11  This model calculates the scattering from a rectangular parallelepiped 
    12  (\:numref:`parallelepiped-image`\). 
    13  If you need to apply polydispersity, see also :ref:`rectangular-prism`. 
     7This model calculates the scattering from a rectangular solid 
     8(:numref:`parallelepiped-image`). 
     9If you need to apply polydispersity, see also :ref:`rectangular-prism`. For 
     10information about polarised and magnetic scattering, see 
     11the :ref:`magnetism` documentation. 
    1412 
    1513.. _parallelepiped-image: 
     
    2119 
    2220The three dimensions of the parallelepiped (strictly here a cuboid) may be 
    23 given in *any* size order. To avoid multiple fit solutions, especially 
    24 with Monte-Carlo fit methods, it may be advisable to restrict their ranges. 
    25 There may be a number of closely similar "best fits", so some trial and 
    26 error, or fixing of some dimensions at expected values, may help. 
    27  
    28 The 1D scattering intensity $I(q)$ is calculated as: 
     21given in *any* size order as long as the particles are randomly oriented (i.e. 
     22take on all possible orientations see notes on 2D below). To avoid multiple fit 
     23solutions, especially with Monte-Carlo fit methods, it may be advisable to 
     24restrict their ranges. There may be a number of closely similar "best fits", so 
     25some trial and error, or fixing of some dimensions at expected values, may 
     26help. 
     27 
     28The form factor is normalized by the particle volume and the 1D scattering 
     29intensity $I(q)$ is then calculated as: 
    2930 
    3031.. Comment by Miguel Gonzalez: 
     
    3940 
    4041    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 
    41            \left< P(q, \alpha) \right> + \text{background} 
     42           \left< P(q, \alpha, \beta) \right> + \text{background} 
    4243 
    4344where the volume $V = A B C$, the contrast is defined as 
    44 $\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, 
    45 $P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented 
    46 at an angle $\alpha$ (angle between the long axis C and $\vec q$), 
    47 and the averaging $\left<\ldots\right>$ is applied over all orientations. 
     45$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$ 
     46is the form factor corresponding to a parallelepiped oriented 
     47at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$ 
     48(the angle between the projection of the particle in the $xy$ detector plane 
     49and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all 
     50orientations. 
    4851 
    4952Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
    50 form factor is given by (Mittelbach and Porod, 1961) 
     53form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_) 
    5154 
    5255.. math:: 
     
    6669    \mu &= qB 
    6770 
    68 The scattering intensity per unit volume is returned in units of |cm^-1|. 
    69  
    70 NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    71 the averaged effective radius, after appropriately sorting the three 
    72 dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and 
    73 length $(= C)$ values, and used as the effective radius for 
    74 $S(q)$ when $P(q) \cdot S(q)$ is applied. 
    75  
    76 For 2d data the orientation of the particle is required, described using 
    77 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
    78 of the calculation and angular dispersions see :ref:`orientation` . 
    79  
    80 .. Comment by Miguel Gonzalez: 
    81    The following text has been commented because I think there are two 
    82    mistakes. Psi is the rotational angle around C (but I cannot understand 
    83    what it means against the q plane) and psi=0 corresponds to a||x and b||y. 
    84  
    85    The angle $\Psi$ is the rotational angle around the $C$ axis against 
    86    the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel 
    87    to the $x$-axis of the detector. 
    88  
    89 The angle $\Psi$ is the rotational angle around the $C$ axis. 
    90 For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis 
    91 oriented parallel to the y-axis of the detector with $A$ along the x-axis. 
    92 For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated 
    93 $\theta$ degrees in the $z-x$ plane and then $\phi$ degrees around the $z$ axis, 
    94 before doing a final rotation of $\Psi$ degrees around the resulting $C$ axis 
    95 of the particle to obtain the final orientation of the parallelepiped. 
    96  
    97 .. _parallelepiped-orientation: 
    98  
    99 .. figure:: img/parallelepiped_angle_definition.png 
    100  
    101     Definition of the angles for oriented parallelepiped, shown with $A<B<C$. 
    102  
    103 .. figure:: img/parallelepiped_angle_projection.png 
    104  
    105     Examples of the angles for an oriented parallelepiped against the 
    106     detector plane. 
    107  
    108 On introducing "Orientational Distribution" in the angles, "distribution of 
    109 theta" and "distribution of phi" parameters will appear. These are actually 
    110 rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, 
    111 perpendicular to the $a$ x $c$ and $b$ x $c$ faces. (When $\theta = \phi = 0$ 
    112 these are parallel to the $Y$ and $X$ axes of the instrument.) The third 
    113 orientation distribution, in $\psi$, is about the $c$ axis of the particle, 
    114 perpendicular to the $a$ x $b$ face. Some experimentation may be required to 
    115 understand the 2d patterns fully as discussed in :ref:`orientation` . 
    116  
    117 For a given orientation of the parallelepiped, the 2D form factor is 
    118 calculated as 
    119  
    120 .. math:: 
    121  
    122     P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2 
    123                   \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2 
    124                   \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2 
    125  
    126 with 
    127  
    128 .. math:: 
    129  
    130     \cos\alpha &= \hat A \cdot \hat q, \\ 
    131     \cos\beta  &= \hat B \cdot \hat q, \\ 
    132     \cos\gamma &= \hat C \cdot \hat q 
    133  
    134 and the scattering intensity as: 
    135  
    136 .. math:: 
    137  
    138     I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y) 
     71where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been 
     72applied. 
     73 
     74For **oriented** particles, the 2D scattering intensity, $I(q_x, q_y)$, is 
     75given as: 
     76 
     77.. math:: 
     78 
     79    I(q_x, q_y) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 P(q_x, q_y) 
    13980            + \text{background} 
    14081 
     
    14889   with scale being the volume fraction. 
    14990 
     91Where $P(q_x, q_y)$ for a given orientation of the form factor is calculated as 
     92 
     93.. math:: 
     94 
     95    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1} 
     96                   {2}qA\cos\alpha)}\right]^2 
     97                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1} 
     98                   {2}qB\cos\beta)}\right]^2 
     99                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1} 
     100                   {2}qC\cos\gamma)}\right]^2 
     101 
     102with 
     103 
     104.. math:: 
     105 
     106    \cos\alpha &= \hat A \cdot \hat q, \\ 
     107    \cos\beta  &= \hat B \cdot \hat q, \\ 
     108    \cos\gamma &= \hat C \cdot \hat q 
     109 
     110 
     111FITTING NOTES 
     112~~~~~~~~~~~~~ 
     113 
     114#. The 2nd virial coefficient of the parallelepiped is calculated based on 
     115   the averaged effective radius, after appropriately sorting the three 
     116   dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and 
     117   length $(= C)$ values, and used as the effective radius for 
     118   $S(q)$ when $P(q) \cdot S(q)$ is applied. 
     119 
     120#. For 2d data the orientation of the particle is required, described using 
     121   angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, where $\theta$ 
     122   and $\phi$ define the orientation of the director in the laboratry reference 
     123   frame of the beam direction ($z$) and detector plane ($x-y$ plane), while 
     124   the angle $\Psi$ is effectively the rotational angle around the particle 
     125   $C$ axis. For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the 
     126   $B$ axis oriented parallel to the y-axis of the detector with $A$ along 
     127   the x-axis. For other $\theta$, $\phi$ values, the order of rotations 
     128   matters. In particular, the parallelepiped must first be rotated $\theta$ 
     129   degrees in the $x-z$ plane before rotating $\phi$ degrees around the $z$ 
     130   axis (in the $x-y$ plane). Applying orientational distribution to the 
     131   particle orientation (i.e  `jitter` to one or more of these angles) can get 
     132   more confusing as `jitter` is defined **NOT** with respect to the laboratory 
     133   frame but the particle reference frame. It is thus highly recmmended to 
     134   read :ref:`orientation` for further details of the calculation and angular 
     135   dispersions. 
     136 
     137.. note:: For 2d, constraints must be applied during fitting to ensure that the 
     138   order of sides chosen is not altered, and hence that the correct definition 
     139   of angles is preserved. For the default choice shown here, that means 
     140   ensuring that the inequality $A < B < C$ is not violated,  The calculation 
     141   will not report an error, but the results may be not correct. 
     142    
     143.. _parallelepiped-orientation: 
     144 
     145.. figure:: img/parallelepiped_angle_definition.png 
     146 
     147    Definition of the angles for oriented parallelepiped, shown with $A<B<C$. 
     148 
     149.. figure:: img/parallelepiped_angle_projection.png 
     150 
     151    Examples of the angles for an oriented parallelepiped against the 
     152    detector plane. 
     153 
     154.. Comment by Paul Butler 
     155   I am commenting this section out as we are trying to minimize the amount of 
     156   oritentational detail here and encourage the user to go to the full 
     157   orientation documentation so that changes can be made in just one place. 
     158   below is the commented paragrah: 
     159   On introducing "Orientational Distribution" in the angles, "distribution of 
     160   theta" and "distribution of phi" parameters will appear. These are actually 
     161   rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, 
     162   perpendicular to the $a$ x $c$ and $b$ x $c$ faces. (When $\theta = \phi = 0$ 
     163   these are parallel to the $Y$ and $X$ axes of the instrument.) The third 
     164   orientation distribution, in $\psi$, is about the $c$ axis of the particle, 
     165   perpendicular to the $a$ x $b$ face. Some experimentation may be required to 
     166   understand the 2d patterns fully as discussed in :ref:`orientation` . 
     167 
    150168 
    151169Validation 
     
    156174angles. 
    157175 
    158  
    159176References 
    160177---------- 
    161178 
    162 P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 
    163  
    164 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     179.. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*, 
     180   14 (1961) 185-211 
     181.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
    165182 
    166183Authorship and Verification 
     
    169186* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    170187* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017 
    171 * **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017 
     188* **Last Reviewed by:**  Miguel Gonzales and Paul Butler **Date:** May 24, 
     189  2018 - documentation updated 
    172190""" 
    173191 
  • sasmodels/models/spherical_sld.py

    r2d81cfe r5601947  
    11r""" 
     2Definition 
     3---------- 
     4 
    25Similarly to the onion, this model provides the form factor, $P(q)$, for 
    36a multi-shell sphere, where the interface between the each neighboring 
     
    1619interface. The form factor is normalized by the total volume of the sphere. 
    1720 
    18 Interface shapes are as follows:: 
     21Interface shapes are as follows: 
    1922 
    2023    0: erf($\nu z$) 
     24     
    2125    1: Rpow($z^\nu$) 
     26     
    2227    2: Lpow($z^\nu$) 
     28     
    2329    3: Rexp($-\nu z$) 
     30     
    2431    4: Lexp($-\nu z$) 
    25  
    26 Definition 
    27 ---------- 
    2832 
    2933The form factor $P(q)$ in 1D is calculated by: 
     
    174178    when $P(Q) * S(Q)$ is applied. 
    175179 
     180 
    176181References 
    177182---------- 
    178 L A Feigin and D I Svergun, Structure Analysis by Small-Angle X-Ray 
    179 and Neutron Scattering, Plenum Press, New York, (1987) 
     183 
     184.. [#] L A Feigin and D I Svergun, Structure Analysis by Small-Angle X-Ray 
     185   and Neutron Scattering, Plenum Press, New York, (1987) 
     186 
     187 
     188Authorship and Verification 
     189---------------------------- 
     190 
     191* **Author:** Jae-Hie Cho **Date:** Nov 1, 2010 
     192* **Last Modified by:** Paul Kienzle **Date:** Dec 20, 2016 
     193* **Last Reviewed by:** Paul Butler **Date:** September 8, 2018 
    180194""" 
    181195 
  • sasmodels/models/spinodal.py

    ref07e95 r475ff58  
    33---------- 
    44 
    5 This model calculates the SAS signal of a phase separating solution 
    6 under spinodal decomposition. The scattering intensity $I(q)$ is calculated as 
     5This model calculates the SAS signal of a phase separating system  
     6undergoing spinodal decomposition. The scattering intensity $I(q)$ is calculated  
     7as  
    78 
    89.. math:: 
    910    I(q) = I_{max}\frac{(1+\gamma/2)x^2}{\gamma/2+x^{2+\gamma}}+B 
    1011 
    11 where $x=q/q_0$ and $B$ is a flat background. The characteristic structure 
    12 length scales with the correlation peak at $q_0$. The exponent $\gamma$ is 
    13 equal to $d+1$ with d the dimensionality of the off-critical concentration 
    14 mixtures. A transition to $\gamma=2d$ is seen near the percolation threshold 
    15 into the critical concentration regime. 
     12where $x=q/q_0$, $q_0$ is the peak position, $I_{max}$ is the intensity  
     13at $q_0$ (parameterised as the $scale$ parameter), and $B$ is a flat  
     14background. The spinodal wavelength is given by $2\pi/q_0$.  
     15 
     16The exponent $\gamma$ is equal to $d+1$ for off-critical concentration  
     17mixtures (smooth interfaces) and $2d$ for critical concentration mixtures  
     18(entangled interfaces), where $d$ is the dimensionality (ie, 1, 2, 3) of the  
     19system. Thus 2 <= $\gamma$ <= 6. A transition from $\gamma=d+1$ to $\gamma=2d$  
     20is expected near the percolation threshold.  
     21 
     22As this function tends to zero as $q$ tends to zero, in practice it may be  
     23necessary to combine it with another function describing the low-angle  
     24scattering, or to simply omit the low-angle scattering from the fit. 
    1625 
    1726References 
     
    2231Physica A 123,497 (1984). 
    2332 
    24 Authorship and Verification 
    25 ---------------------------- 
     33Revision History 
     34---------------- 
    2635 
    27 * **Author:** Dirk Honecker **Date:** Oct 7, 2016 
     36* **Author:**  Dirk Honecker **Date:** Oct 7, 2016 
     37* **Revised:** Steve King    **Date:** Sep 7, 2018 
    2838""" 
    2939 
     
    3444title = "Spinodal decomposition model" 
    3545description = """\ 
    36       I(q) = scale ((1+gamma/2)x^2)/(gamma/2+x^(2+gamma))+background 
     46      I(q) = Imax ((1+gamma/2)x^2)/(gamma/2+x^(2+gamma)) + background 
    3747 
    3848      List of default parameters: 
    39       scale = scaling 
    40       gamma = exponent 
    41       x = q/q_0 
     49       
     50      Imax = correlation peak intensity at q_0 
     51      background = incoherent background 
     52      gamma = exponent (see model documentation) 
    4253      q_0 = correlation peak position [1/A] 
    43       background = Incoherent background""" 
     54      x = q/q_0""" 
     55       
    4456category = "shape-independent" 
    4557 
  • sasmodels/models/vesicle.py

    ref07e95 rb477605  
    33---------- 
    44 
    5 The 1D scattering intensity is calculated in the following way (Guinier, 1955) 
     5This model provides the form factor, *P(q)*, for an unilamellar vesicle and is 
     6effectively identical to the hollow sphere reparameterized to be 
     7more intuitive for a vesicle and normalizing the form factor by the volume of 
     8the shell. The 1D scattering intensity is calculated in the following way 
     9(Guinier,1955\ [#Guinier1955]_) 
    610 
    711.. math:: 
     
    5357---------- 
    5458 
    55 A Guinier and G. Fournet, *Small-Angle Scattering of X-Rays*, John Wiley and 
    56 Sons, New York, (1955) 
     59.. [#Guinier1955] A Guinier and G. Fournet, *Small-Angle Scattering of X-Rays*, John Wiley and 
     60   Sons, New York, (1955) 
     61 
     62 
     63Authorship and Verification 
     64---------------------------- 
    5765 
    5866* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    5967* **Last Modified by:** Paul Butler **Date:** March 20, 2016 
    60 * **Last Reviewed by:** Paul Butler **Date:** March 20, 2016 
     68* **Last Reviewed by:** Paul Butler **Date:** September 7, 2018 
    6169""" 
    6270 
     
    6573 
    6674name = "vesicle" 
    67 title = "This model provides the form factor, *P(q)*, for an unilamellar \ 
    68     vesicle. This is model is effectively identical to the hollow sphere \ 
    69     reparameterized to be more intuitive for a vesicle and normalizing the \ 
    70     form factor by the volume of the shell." 
     75title = "Vesicle model representing a hollow sphere" 
    7176description = """ 
    7277    Model parameters: 
  • sasmodels/resolution.py

    r2d81cfe r9e7837a  
    2020MINIMUM_RESOLUTION = 1e-8 
    2121MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data 
     22# According to (Barker & Pedersen 1995 JAC), 2.5 sigma is a good limit. 
     23# According to simulations with github.com:scattering/sansresolution.git 
     24# it is better to use asymmetric bounds (2.5, 3.0) 
     25PINHOLE_N_SIGMA = (2.5, 3.0) 
    2226 
    2327class Resolution(object): 
     
    6569    *q_calc* is the list of points to calculate, or None if this should 
    6670    be estimated from the *q* and *q_width*. 
    67     """ 
    68     def __init__(self, q, q_width, q_calc=None, nsigma=3): 
     71 
     72    *nsigma* is the width of the resolution function.  Should be 2.5. 
     73    See :func:`pinhole_resolution` for details. 
     74    """ 
     75    def __init__(self, q, q_width, q_calc=None, nsigma=PINHOLE_N_SIGMA): 
    6976        #*min_step* is the minimum point spacing to use when computing the 
    7077        #underlying model.  It should be on the order of 
     
    8289 
    8390        # Protect against models which are not defined for very low q.  Limit 
    84         # the smallest q value evaluated (in absolute) to 0.02*min 
    85         cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
    86         self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff] 
     91        # the smallest q value evaluated to 0.02*min.  Note that negative q 
     92        # values are trimmed even for broad resolution.  Although not possible 
     93        # from the geometry, they may appear since we are using a truncated 
     94        # gaussian to represent resolution rather than a skew distribution. 
     95        #cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
     96        #self.q_calc = self.q_calc[self.q_calc >= cutoff] 
    8797 
    8898        # Build weight matrix from calculated q values 
    8999        self.weight_matrix = pinhole_resolution( 
    90             self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) 
    91         self.q_calc = abs(self.q_calc) 
     100            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION), 
     101            nsigma=nsigma) 
    92102 
    93103    def apply(self, theory): 
     
    101111    *q* points at which the data is measured. 
    102112 
    103     *dqx* slit width in qx 
    104  
    105     *dqy* slit height in qy 
     113    *qx_width* slit width in qx 
     114 
     115    *qy_width* slit height in qy 
    106116 
    107117    *q_calc* is the list of points to calculate, or None if this should 
     
    154164 
    155165 
    156 def pinhole_resolution(q_calc, q, q_width): 
    157     """ 
     166def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA): 
     167    r""" 
    158168    Compute the convolution matrix *W* for pinhole resolution 1-D data. 
    159169 
     
    162172    *W*, the resolution smearing can be computed using *dot(W,q)*. 
    163173 
     174    Note that resolution is limited to $\pm 2.5 \sigma$.[1]  The true resolution 
     175    function is a broadened triangle, and does not extend over the entire 
     176    range $(-\infty, +\infty)$.  It is important to impose this limitation 
     177    since some models fall so steeply that the weighted value in gaussian 
     178    tails would otherwise dominate the integral. 
     179 
    164180    *q_calc* must be increasing.  *q_width* must be greater than zero. 
     181 
     182    [1] Barker, J. G., and J. S. Pedersen. 1995. Instrumental Smearing Effects 
     183    in Radially Symmetric Small-Angle Neutron Scattering by Numerical and 
     184    Analytical Methods. Journal of Applied Crystallography 28 (2): 105--14. 
     185    https://doi.org/10.1107/S0021889894010095. 
    165186    """ 
    166187    # The current algorithm is a midpoint rectangle rule.  In the test case, 
     
    170191    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    171192    weights = cdf[1:] - cdf[:-1] 
     193    # Limit q range to (-2.5,+3) sigma 
     194    try: 
     195        nsigma_low, nsigma_high = nsigma 
     196    except TypeError: 
     197        nsigma_low = nsigma_high = nsigma 
     198    qhigh = q + nsigma_high*q_width 
     199    qlow = q - nsigma_low*q_width  # linear limits 
     200    ##qlow = q*q/qhigh  # log limits 
     201    weights[q_calc[:, None] < qlow[None, :]] = 0. 
     202    weights[q_calc[:, None] > qhigh[None, :]] = 0. 
    172203    weights /= np.sum(weights, axis=0)[None, :] 
    173204    return weights 
     
    341372 
    342373 
    343 def pinhole_extend_q(q, q_width, nsigma=3): 
     374def pinhole_extend_q(q, q_width, nsigma=PINHOLE_N_SIGMA): 
    344375    """ 
    345376    Given *q* and *q_width*, find a set of sampling points *q_calc* so 
     
    347378    function. 
    348379    """ 
    349     q_min, q_max = np.min(q - nsigma*q_width), np.max(q + nsigma*q_width) 
     380    try: 
     381        nsigma_low, nsigma_high = nsigma 
     382    except TypeError: 
     383        nsigma_low = nsigma_high = nsigma 
     384    q_min, q_max = np.min(q - nsigma_low*q_width), np.max(q + nsigma_high*q_width) 
    350385    return linear_extrapolation(q, q_min, q_max) 
    351386 
     
    494529 
    495530 
    496 def gaussian(q, q0, dq): 
    497     """ 
    498     Return the Gaussian resolution function. 
     531def gaussian(q, q0, dq, nsigma=2.5): 
     532    """ 
     533    Return the truncated Gaussian resolution function. 
    499534 
    500535    *q0* is the center, *dq* is the width and *q* are the points to evaluate. 
    501536    """ 
    502     return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 
     537    # Calculate the density of the tails; the resulting gaussian needs to be 
     538    # scaled by this amount in ordere to integrate to 1.0 
     539    two_tail_density = 2 * (1 + erf(-nsigma/sqrt(2)))/2 
     540    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density) 
    503541 
    504542 
     
    558596 
    559597 
    560 def romberg_pinhole_1d(q, q_width, form, pars, nsigma=5): 
     598def romberg_pinhole_1d(q, q_width, form, pars, nsigma=2.5): 
    561599    """ 
    562600    Romberg integration for pinhole resolution. 
     
    678716        np.testing.assert_equal(output, self.y) 
    679717 
     718    # TODO: turn pinhole/slit demos into tests 
     719 
     720    @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    680721    def test_pinhole(self): 
    681722        """ 
     
    686727        theory = 12.0-1000.0*resolution.q_calc 
    687728        output = resolution.apply(theory) 
     729        # Note: answer came from output of previous run.  Non-integer 
     730        # values at ends come from the fact that q_calc does not 
     731        # extend beyond q, and so the weights don't balance. 
    688732        answer = [ 
    689             10.44785079, 9.84991299, 8.98101708, 
    690             7.99906585, 6.99998311, 6.00001689, 
    691             5.00093415, 4.01898292, 3.15008701, 2.55214921, 
     733            10.47037734, 9.86925860, 
     734            9., 8., 7., 6., 5., 4., 
     735            3.13074140, 2.52962266, 
    692736            ] 
    693737        np.testing.assert_allclose(output, answer, atol=1e-8) 
     
    732776        self._compare(q, output, answer, 1e-6) 
    733777 
     778    @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    734779    def test_pinhole(self): 
    735780        """ 
     
    746791        self._compare(q, output, answer, 3e-4) 
    747792 
     793    @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    748794    def test_pinhole_romberg(self): 
    749795        """ 
     
    761807        #                     2*np.pi/pars['radius']/200) 
    762808        #tol = 0.001 
    763         ## The default 3 sigma and no extra points gets 1% 
     809        ## The default 2.5 sigma and no extra points gets 1% 
    764810        q_calc = None  # type: np.ndarray 
    765811        tol = 0.01 
     
    10801126 
    10811127    if isinstance(resolution, Slit1D): 
    1082         width, height = resolution.dqx, resolution.dqy 
     1128        width, height = resolution.qx_width, resolution.qy_width 
    10831129        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars) 
    10841130    else: 
Note: See TracChangeset for help on using the changeset viewer.