Changes in / [9b7fac6:1a8b91c] in sasmodels


Ignore:
Files:
1 added
1 deleted
22 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/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 rd712a0f  
    2020  P(q) = \text{scale} \langle F^* F \rangle / V + \text{background} 
    2121 
    22 where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an 
    23 average over the size distribution. 
     22where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an  
     23average over the size distribution $f(x; \bar x, \sigma)$, giving 
     24 
     25.. math:: 
     26 
     27  P(q) = \frac{\text{scale}}{V} \int_\mathbb{R}  
     28  f(x; \bar x, \sigma) F^2(q, x)\, dx + \text{background} 
    2429 
    2530Each distribution is characterized by a center value $\bar x$ or 
    … …  
    4146with larger values of $N_\sigma$ required for heavier tailed distributions. 
    4247The 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)$ 
     48that $f(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)f(r - 3\sigma_r)$ 
    4449will not contribute much to the average may not hold when particles are large. 
    4550This, too, will require increasing $N_\sigma$. 
    … …  
    6368 
    6469Additional distributions are under consideration. 
     70 
     71.. note:: In 2009 IUPAC decided to introduce the new term 'dispersity' to replace  
     72           the term 'polydispersity' (see `Pure Appl. Chem., (2009), 81(2),  
     73           351-353 <http://media.iupac.org/publications/pac/2009/pdf/8102x0351.pdf>`_  
     74           in order to make the terminology describing distributions of properties  
     75           unambiguous. Throughout the SasView documentation we continue to use the  
     76           term polydispersity because one of the consequences of the IUPAC change is  
     77           that orientational polydispersity would not meet their new criteria (which  
     78           requires dispersity to be dimensionless). 
    6579 
    6680Suggested Applications 

  • doc/guide/plugin.rst

    r7e6bc45e rf796469  
    822822 
    823823        :code:`source = ["lib/Si.c", ...]` 
    824         (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/Si.c>`_) 
     824        (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_Si.c>`_) 
    825825 
    826826    sas_3j1x_x(x): 

  • 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

    r5770493 raf7a97c  
    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        data = load_data(os.path.expanduser(opts['datafile'])) 
     1291    else: 
     1292        # Hack around the fact that make_data doesn't take a pair of resolutions 
     1293        res = opts['res'] 
     1294        opts['res'] = res[0] 
     1295        data0, _ = make_data(opts) 
     1296        if res[0] != res[1]: 
     1297            opts['res'] = res[1] 
     1298            data1, _ = make_data(opts) 
     1299        else: 
     1300            data1 = data0 
     1301        opts['res'] = res 
     1302        data = data0, data1 
     1303 
     1304    base = make_engine(model_info[0], data[0], opts['engine'][0], 
    12191305                       opts['cutoff'][0], opts['ngauss'][0]) 
    12201306    if comparison: 
    1221         comp = make_engine(model_info[1], data, opts['engine'][1], 
     1307        comp = make_engine(model_info[1], data[1], opts['engine'][1], 
    12221308                           opts['cutoff'][1], opts['ngauss'][1]) 
    12231309    else: 
    … …  
    13761462    path = os.path.dirname(info.filename) 
    13771463    url = "file://" + path.replace("\\", "/")[2:] + "/" 
    1378     rst2html.view_html_qtapp(html, url) 
     1464    rst2html.view_html_wxapp(html, url) 
    13791465 
    13801466def explore(opts): 

  • sasmodels/core.py

    r3221de0 r4341dd4  
    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 

  • sasmodels/data.py

    rd86f0fc r1a8c11c  
    3636 
    3737import numpy as np  # type: ignore 
     38from numpy import sqrt, sin, cos, pi 
    3839 
    3940# pylint: disable=unused-import 
    … …  
    301302 
    302303 
    303 def empty_data1D(q, resolution=0.0): 
     304def empty_data1D(q, resolution=0.0, L=0., dL=0.): 
    304305    # type: (np.ndarray, float) -> Data1D 
    305     """ 
     306    r""" 
    306307    Create empty 1D data using the given *q* as the x value. 
    307308 
    308     *resolution* dq/q defaults to 5%. 
     309    rms *resolution* $\Delta q/q$ defaults to 0%.  If wavelength *L* and rms 
     310    wavelength divergence *dL* are defined, then *resolution* defines 
     311    rms $\Delta \theta/\theta$ for the lowest *q*, with $\theta$ derived from 
     312    $q = 4\pi/\lambda \sin(\theta)$. 
    309313    """ 
    310314 
    … …  
    313317    Iq, dIq = None, None 
    314318    q = np.asarray(q) 
    315     data = Data1D(q, Iq, dx=resolution * q, dy=dIq) 
     319    if L != 0 and resolution != 0: 
     320        theta = np.arcsin(q*L/(4*pi)) 
     321        dtheta = theta[0]*resolution 
     322        ## Solving Gaussian error propagation from 
     323        ##   Dq^2 = (dq/dL)^2 DL^2 + (dq/dtheta)^2 Dtheta^2 
     324        ## gives 
     325        ##   (Dq/q)^2 = (DL/L)**2 + (Dtheta/tan(theta))**2 
     326        ## Take the square root and multiply by q, giving 
     327        ##   Dq = (4*pi/L) * sqrt((sin(theta)*dL/L)**2 + (cos(theta)*dtheta)**2) 
     328        dq = (4*pi/L) * sqrt((sin(theta)*dL/L)**2 + (cos(theta)*dtheta)**2) 
     329    else: 
     330        dq = resolution * q 
     331 
     332    data = Data1D(q, Iq, dx=dq, dy=dIq) 
    316333    data.filename = "fake data" 
    317334    return data 
    … …  
    486503            # Note: masks merge, so any masked theory points will stay masked, 
    487504            # and the data mask will be added to it. 
    488             mtheory = masked_array(theory, data.mask.copy()) 
     505            #mtheory = masked_array(theory, data.mask.copy()) 
     506            theory_x = data.x[~data.mask] 
     507            mtheory = masked_array(theory) 
    489508            mtheory[~np.isfinite(mtheory)] = masked 
    490509            if view is 'log': 
    491510                mtheory[mtheory <= 0] = masked 
    492             plt.plot(data.x, scale*mtheory, '-') 
     511            plt.plot(theory_x, scale*mtheory, '-') 
    493512            all_positive = all_positive and (mtheory > 0).all() 
    494513            some_present = some_present or (mtheory.count() > 0) 
    … …  
    526545 
    527546    if use_resid: 
    528         mresid = masked_array(resid, data.mask.copy()) 
     547        theory_x = data.x[~data.mask] 
     548        mresid = masked_array(resid) 
    529549        mresid[~np.isfinite(mresid)] = masked 
    530550        some_present = (mresid.count() > 0) 
    … …  
    532552        if num_plots > 1: 
    533553            plt.subplot(1, num_plots, use_calc + 2) 
    534         plt.plot(data.x, mresid, '.') 
     554        plt.plot(theory_x, mresid, '.') 
    535555        plt.xlabel("$q$/A$^{-1}$") 
    536556        plt.ylabel('residuals') 

  • sasmodels/direct_model.py

    rd18d6dd r1a8c11c  
    250250            qmax = getattr(data, 'qmax', np.inf) 
    251251            accuracy = getattr(data, 'accuracy', 'Low') 
    252             index = ~data.mask & (q >= qmin) & (q <= qmax) 
     252            index = (data.mask == 0) & (q >= qmin) & (q <= qmax) 
    253253            if data.data is not None: 
    254254                index &= ~np.isnan(data.data) 
    … …  
    263263        elif self.data_type == 'Iq': 
    264264            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     265            mask = getattr(data, 'mask', None) 
     266            if mask is not None: 
     267                index &= (mask == 0) 
    265268            if data.y is not None: 
    266269                index &= ~np.isnan(data.y) 

  • 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 r7c35fda  
    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 

  • sasmodels/kerneldll.py

    rbf94e6e r1a3559f  
    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: 
    … …  
    123129    # add openmp support if not running on a mac 
    124130    if sys.platform != "darwin": 
    125         CC.append("-fopenmp") 
     131        # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9) 
     132        # Shut it off for all unix until we can investigate. 
     133        #CC.append("-fopenmp") 
     134        pass 
    126135    def compile_command(source, output): 
    127136        """unix compiler command""" 
    … …  
    158167        return CC + [source, "-o", output, "-lm"] 
    159168 
    160 # Assume the default location of module DLLs is in .sasmodels/compiled_models. 
    161 DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    162  
    163169ALLOW_SINGLE_PRECISION_DLLS = True 
    164170 
    … …  
    197203        return path 
    198204 
    199     return joinpath(DLL_PATH, basename) 
     205    return joinpath(SAS_DLL_PATH, basename) 
    200206 
    201207 
    … …  
    206212    exist yet if it hasn't been compiled. 
    207213    """ 
    208     return os.path.join(DLL_PATH, dll_name(model_info, dtype)) 
     214    return os.path.join(SAS_DLL_PATH, dll_name(model_info, dtype)) 
    209215 
    210216 
    … …  
    225231    models are not allowed as DLLs. 
    226232 
    227     Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 
     233    Set *sasmodels.kerneldll.SAS_DLL_PATH* to the compiled dll output path. 
     234    Alternatively, set the environment variable *SAS_DLL_PATH*. 
    228235    The default is in ~/.sasmodels/compiled_models. 
    229236    """ 
    … …  
    244251    if need_recompile: 
    245252        # Make sure the DLL path exists 
    246         if not os.path.exists(DLL_PATH): 
    247             os.makedirs(DLL_PATH) 
     253        if not os.path.exists(SAS_DLL_PATH): 
     254            os.makedirs(SAS_DLL_PATH) 
    248255        basename = splitext(os.path.basename(dll))[0] + "_" 
    249256        system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 

  • 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 rdc76240  
    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} 

  • 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/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/resolution.py

    r2d81cfe r0b9c6df  
    2020MINIMUM_RESOLUTION = 1e-8 
    2121MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data 
     22PINHOLE_N_SIGMA = 2.5 # From: Barker & Pedersen 1995 JAC 
    2223 
    2324class Resolution(object): 
    … …  
    6566    *q_calc* is the list of points to calculate, or None if this should 
    6667    be estimated from the *q* and *q_width*. 
    67     """ 
    68     def __init__(self, q, q_width, q_calc=None, nsigma=3): 
     68 
     69    *nsigma* is the width of the resolution function.  Should be 2.5. 
     70    See :func:`pinhole_resolution` for details. 
     71    """ 
     72    def __init__(self, q, q_width, q_calc=None, nsigma=PINHOLE_N_SIGMA): 
    6973        #*min_step* is the minimum point spacing to use when computing the 
    7074        #underlying model.  It should be on the order of 
    … …  
    8286 
    8387        # Protect against models which are not defined for very low q.  Limit 
    84         # the smallest q value evaluated (in absolute) to 0.02*min 
     88        # the smallest q value evaluated to 0.02*min.  Note that negative q 
     89        # values are trimmed even for broad resolution.  Although not possible 
     90        # from the geometry, they may appear since we are using a truncated 
     91        # gaussian to represent resolution rather than a skew distribution. 
    8592        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
    86         self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff] 
     93        self.q_calc = self.q_calc[self.q_calc >= cutoff] 
    8794 
    8895        # Build weight matrix from calculated q values 
    8996        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) 
     97            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION), 
     98            nsigma=nsigma) 
    9299 
    93100    def apply(self, theory): 
    … …  
    101108    *q* points at which the data is measured. 
    102109 
    103     *dqx* slit width in qx 
    104  
    105     *dqy* slit height in qy 
     110    *qx_width* slit width in qx 
     111 
     112    *qy_width* slit height in qy 
    106113 
    107114    *q_calc* is the list of points to calculate, or None if this should 
    … …  
    154161 
    155162 
    156 def pinhole_resolution(q_calc, q, q_width): 
    157     """ 
     163def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA): 
     164    r""" 
    158165    Compute the convolution matrix *W* for pinhole resolution 1-D data. 
    159166 
    … …  
    162169    *W*, the resolution smearing can be computed using *dot(W,q)*. 
    163170 
     171    Note that resolution is limited to $\pm 2.5 \sigma$.[1]  The true resolution 
     172    function is a broadened triangle, and does not extend over the entire 
     173    range $(-\infty, +\infty)$.  It is important to impose this limitation 
     174    since some models fall so steeply that the weighted value in gaussian 
     175    tails would otherwise dominate the integral. 
     176 
    164177    *q_calc* must be increasing.  *q_width* must be greater than zero. 
     178 
     179    [1] Barker, J. G., and J. S. Pedersen. 1995. Instrumental Smearing Effects 
     180    in Radially Symmetric Small-Angle Neutron Scattering by Numerical and 
     181    Analytical Methods. Journal of Applied Crystallography 28 (2): 105--14. 
     182    https://doi.org/10.1107/S0021889894010095. 
    165183    """ 
    166184    # The current algorithm is a midpoint rectangle rule.  In the test case, 
    … …  
    170188    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    171189    weights = cdf[1:] - cdf[:-1] 
     190    # Limit q range to +/- 2.5 sigma 
     191    qhigh = q + nsigma*q_width 
     192    #qlow = q - nsigma*q_width  # linear limits 
     193    qlow = q*q/qhigh  # log limits 
     194    weights[q_calc[:, None] < qlow[None, :]] = 0. 
     195    weights[q_calc[:, None] > qhigh[None, :]] = 0. 
    172196    weights /= np.sum(weights, axis=0)[None, :] 
    173197    return weights 
    … …  
    494518 
    495519 
    496 def gaussian(q, q0, dq): 
    497     """ 
    498     Return the Gaussian resolution function. 
     520def gaussian(q, q0, dq, nsigma=2.5): 
     521    """ 
     522    Return the truncated Gaussian resolution function. 
    499523 
    500524    *q0* is the center, *dq* is the width and *q* are the points to evaluate. 
    501525    """ 
    502     return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 
     526    # Calculate the density of the tails; the resulting gaussian needs to be 
     527    # scaled by this amount in ordere to integrate to 1.0 
     528    two_tail_density = 2 * (1 + erf(-nsigma/sqrt(2)))/2 
     529    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density) 
    503530 
    504531 
    … …  
    558585 
    559586 
    560 def romberg_pinhole_1d(q, q_width, form, pars, nsigma=5): 
     587def romberg_pinhole_1d(q, q_width, form, pars, nsigma=2.5): 
    561588    """ 
    562589    Romberg integration for pinhole resolution. 
    … …  
    678705        np.testing.assert_equal(output, self.y) 
    679706 
     707    # TODO: turn pinhole/slit demos into tests 
     708 
     709    @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    680710    def test_pinhole(self): 
    681711        """ 
    … …  
    686716        theory = 12.0-1000.0*resolution.q_calc 
    687717        output = resolution.apply(theory) 
     718        # Note: answer came from output of previous run.  Non-integer 
     719        # values at ends come from the fact that q_calc does not 
     720        # extend beyond q, and so the weights don't balance. 
    688721        answer = [ 
    689             10.44785079, 9.84991299, 8.98101708, 
    690             7.99906585, 6.99998311, 6.00001689, 
    691             5.00093415, 4.01898292, 3.15008701, 2.55214921, 
     722            10.47037734, 9.86925860, 
     723            9., 8., 7., 6., 5., 4., 
     724            3.13074140, 2.52962266, 
    692725            ] 
    693726        np.testing.assert_allclose(output, answer, atol=1e-8) 
    … …  
    732765        self._compare(q, output, answer, 1e-6) 
    733766 
     767    @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    734768    def test_pinhole(self): 
    735769        """ 
    … …  
    746780        self._compare(q, output, answer, 3e-4) 
    747781 
     782    @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    748783    def test_pinhole_romberg(self): 
    749784        """ 
    … …  
    761796        #                     2*np.pi/pars['radius']/200) 
    762797        #tol = 0.001 
    763         ## The default 3 sigma and no extra points gets 1% 
     798        ## The default 2.5 sigma and no extra points gets 1% 
    764799        q_calc = None  # type: np.ndarray 
    765800        tol = 0.01 
    … …  
    10801115 
    10811116    if isinstance(resolution, Slit1D): 
    1082         width, height = resolution.dqx, resolution.dqy 
     1117        width, height = resolution.qx_width, resolution.qy_width 
    10831118        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars) 
    10841119    else: 

  • sasmodels/rst2html.py

    r2d81cfe r1fbadb2  
    5656    # others don't work properly with math_output! 
    5757    if math_output == "mathjax": 
    58         settings = {"math_output": math_output} 
     58        # TODO: this is copied from docs/conf.py; there should be only one 
     59        mathjax_path = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_CHTML" 
     60        settings = {"math_output": math_output + " " + mathjax_path} 
    5961    else: 
    6062        settings = {"math-output": math_output} 

  • sasmodels/sasview_model.py

    r05df1de rd533590  
    686686        self._intermediate_results = getattr(calculator, 'results', None) 
    687687        calculator.release() 
    688         self._model.release() 
     688        #self._model.release() 
    689689        return result 
    690690 
    … …  
    719719 
    720720    def set_dispersion(self, parameter, dispersion): 
    721         # type: (str, weights.Dispersion) -> Dict[str, Any] 
     721        # type: (str, weights.Dispersion) -> None 
    722722        """ 
    723723        Set the dispersion object for a model parameter 
    … …  
    742742            self.dispersion[parameter] = dispersion.get_pars() 
    743743        else: 
    744             raise ValueError("%r is not a dispersity or orientation parameter") 
     744            raise ValueError("%r is not a dispersity or orientation parameter" 
     745                             % parameter) 
    745746 
    746747    def _dispersion_mesh(self): 
Note: See TracChangeset for help on using the changeset viewer.