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


Ignore:
Files:
1 added
1 deleted
22 edited

Legend:

Unmodified
Added
Removed

  • .travis.yml

    r5c36bf1 r335271e  
    66    env: 
    77    - PY=2.7 
    8     #- DEPLOY=True 
     8    - DEPLOY=True 
    99  - os: linux 
    1010    env: 

  • doc/guide/magnetism/magnetism.rst

    rbefe905 r4f5afc9  
    3939 
    4040.. math:: 
    41     -- &= (1-u_i)(1-u_f) \\ 
    42     -+ &= (1-u_i)(u_f) \\ 
    43     +- &= (u_i)(1-u_f) \\ 
    44     ++ &= (u_i)(u_f) 
     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} 
    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

    rd712a0f r29afc50  
    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 $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} 
     22where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an 
     23average over the size distribution. 
    2924 
    3025Each distribution is characterized by a center value $\bar x$ or 
    … …  
    4641with larger values of $N_\sigma$ required for heavier tailed distributions. 
    4742The scattering in general falls rapidly with $qr$ so the usual assumption 
    48 that $f(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)f(r - 3\sigma_r)$ 
     43that $G(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)G(r - 3\sigma_r)$ 
    4944will not contribute much to the average may not hold when particles are large. 
    5045This, too, will require increasing $N_\sigma$. 
    … …  
    6863 
    6964Additional 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). 
    7965 
    8066Suggested Applications 

  • doc/guide/plugin.rst

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

  • example/model_ellipsoid_hayter_msa.py

    r9a99993 r8a5f021  
    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

    raf7a97c r5770493  
    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 
    110109 
    111110    === output options === 
    112111    -edit starts the parameter explorer 
    113112    -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 
    122113 
    123114The interpretation of quad precision depends on architecture, and may 
    … …  
    654645 
    655646def make_data(opts): 
    656     # type: (Dict[str, Any], float) -> Tuple[Data, np.ndarray] 
     647    # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray] 
    657648    """ 
    658649    Generate an empty dataset, used with the model to set Q points 
    … …  
    676667        if opts['zero']: 
    677668            q = np.hstack((0, q)) 
    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) 
     669        data = empty_data1D(q, resolution=res) 
    682670        index = slice(None, None) 
    683671    return data, index 
    … …  
    784772            dim = base._kernel.dim 
    785773            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() 
    803774    if opts['plot']: 
    804775        import matplotlib.pyplot as plt 
    … …  
    806777    return limits 
    807778 
    808 def 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  
    837779def run_models(opts, verbose=False): 
    838780    # type: (Dict[str, Any]) -> Dict[str, Any] 
    … …  
    844786    base_n, comp_n = opts['count'] 
    845787    base_pars, comp_pars = opts['pars'] 
    846     base_data, comp_data = opts['data'] 
     788    data = opts['data'] 
    847789 
    848790    comparison = comp is not None 
    … …  
    858800            print("%s t=%.2f ms, intensity=%.0f" 
    859801                  % (base.engine, base_time, base_value.sum())) 
    860         _show_invalid(base_data, base_value) 
     802        _show_invalid(data, base_value) 
    861803    except ImportError: 
    862804        traceback.print_exc() 
    … …  
    870812                print("%s t=%.2f ms, intensity=%.0f" 
    871813                      % (comp.engine, comp_time, comp_value.sum())) 
    872             _show_invalid(base_data, comp_value) 
     814            _show_invalid(data, comp_value) 
    873815        except ImportError: 
    874816            traceback.print_exc() 
    … …  
    924866    have_base, have_comp = (base_value is not None), (comp_value is not None) 
    925867    base, comp = opts['engines'] 
    926     base_data, comp_data = opts['data'] 
     868    data = opts['data'] 
    927869    use_data = (opts['datafile'] is not None) and (have_base ^ have_comp) 
    928870 
    929871    # Plot if requested 
    930872    view = opts['view'] 
    931     #view = 'log' 
    932873    if limits is None: 
    933874        vmin, vmax = np.inf, -np.inf 
    … …  
    943884        if have_comp: 
    944885            plt.subplot(131) 
    945         plot_theory(base_data, base_value, view=view, use_data=use_data, limits=limits) 
     886        plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
    946887        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    947888        #cbar_title = "log I" 
    … …  
    950891            plt.subplot(132) 
    951892        if not opts['is2d'] and have_base: 
    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) 
     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) 
    954895        plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 
    955896        #cbar_title = "log I" 
    … …  
    967908            err[err > cutoff] = cutoff 
    968909        #err,errstr = base/comp,"ratio" 
    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) 
     910        plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
    973911        plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear') 
    974912        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best') 
    … …  
    1004942OPTIONS = [ 
    1005943    # Plotting 
    1006     'plot', 'noplot', 
    1007     'weights', 'profile', 
     944    'plot', 'noplot', 'weights', 
    1008945    'linear', 'log', 'q4', 
    1009946    'rel', 'abs', 
    … …  
    11211058 
    11221059    invalid = [o[1:] for o in flags 
    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'))] 
     1060               if o[1:] not in NAME_OPTIONS 
     1061               and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] 
    11261062    if invalid: 
    11271063        print("Invalid options: %s"%(", ".join(invalid))) 
    … …  
    11391075        'qmax'      : 0.05, 
    11401076        'nq'        : 128, 
    1141         'res'       : '0.0', 
     1077        'res'       : 0.0, 
    11421078        'noise'     : 0.0, 
    11431079        'accuracy'  : 'Low', 
    … …  
    11601096        'count'     : '1', 
    11611097        'show_weights' : False, 
    1162         'show_profile' : False, 
    11631098        'sphere'    : 0, 
    11641099        'ngauss'    : '0', 
    … …  
    11801115        elif arg.startswith('-q='): 
    11811116            opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')] 
    1182         elif arg.startswith('-res='):      opts['res'] = arg[5:] 
     1117        elif arg.startswith('-res='):      opts['res'] = float(arg[5:]) 
    11831118        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:]) 
    11841119        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:]) 
    … …  
    12211156        elif arg == '-default': opts['use_demo'] = False 
    12221157        elif arg == '-weights': opts['show_weights'] = True 
    1223         elif arg == '-profile': opts['show_profile'] = True 
    12241158        elif arg == '-html':    opts['html'] = True 
    12251159        elif arg == '-help':    opts['html'] = True 
    1226         elif arg.startswith('-D'): 
    1227             var, val = arg[2:].split('=') 
    1228             os.environ[var] = val 
    12291160    # pylint: enable=bad-whitespace,C0321 
    12301161 
    … …  
    12421173    if opts['qmin'] is None: 
    12431174        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) 
    12441179 
    12451180    comparison = any(PAR_SPLIT in v for v in values) 
    … …  
    12811216        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12821217 
    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], 
     1218    base = make_engine(model_info[0], data, opts['engine'][0], 
    13051219                       opts['cutoff'][0], opts['ngauss'][0]) 
    13061220    if comparison: 
    1307         comp = make_engine(model_info[1], data[1], opts['engine'][1], 
     1221        comp = make_engine(model_info[1], data, opts['engine'][1], 
    13081222                           opts['cutoff'][1], opts['ngauss'][1]) 
    13091223    else: 
    … …  
    14621376    path = os.path.dirname(info.filename) 
    14631377    url = "file://" + path.replace("\\", "/")[2:] + "/" 
    1464     rst2html.view_html_wxapp(html, url) 
     1378    rst2html.view_html_qtapp(html, url) 
    14651379 
    14661380def explore(opts): 

  • sasmodels/core.py

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

    r1a8c11c rd86f0fc  
    3636 
    3737import numpy as np  # type: ignore 
    38 from numpy import sqrt, sin, cos, pi 
    3938 
    4039# pylint: disable=unused-import 
    … …  
    302301 
    303302 
    304 def empty_data1D(q, resolution=0.0, L=0., dL=0.): 
     303def empty_data1D(q, resolution=0.0): 
    305304    # type: (np.ndarray, float) -> Data1D 
    306     r""" 
     305    """ 
    307306    Create empty 1D data using the given *q* as the x value. 
    308307 
    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)$. 
     308    *resolution* dq/q defaults to 5%. 
    313309    """ 
    314310 
    … …  
    317313    Iq, dIq = None, None 
    318314    q = np.asarray(q) 
    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) 
     315    data = Data1D(q, Iq, dx=resolution * q, dy=dIq) 
    333316    data.filename = "fake data" 
    334317    return data 
    … …  
    503486            # Note: masks merge, so any masked theory points will stay masked, 
    504487            # and the data mask will be added to it. 
    505             #mtheory = masked_array(theory, data.mask.copy()) 
    506             theory_x = data.x[~data.mask] 
    507             mtheory = masked_array(theory) 
     488            mtheory = masked_array(theory, data.mask.copy()) 
    508489            mtheory[~np.isfinite(mtheory)] = masked 
    509490            if view is 'log': 
    510491                mtheory[mtheory <= 0] = masked 
    511             plt.plot(theory_x, scale*mtheory, '-') 
     492            plt.plot(data.x, scale*mtheory, '-') 
    512493            all_positive = all_positive and (mtheory > 0).all() 
    513494            some_present = some_present or (mtheory.count() > 0) 
    … …  
    545526 
    546527    if use_resid: 
    547         theory_x = data.x[~data.mask] 
    548         mresid = masked_array(resid) 
     528        mresid = masked_array(resid, data.mask.copy()) 
    549529        mresid[~np.isfinite(mresid)] = masked 
    550530        some_present = (mresid.count() > 0) 
    … …  
    552532        if num_plots > 1: 
    553533            plt.subplot(1, num_plots, use_calc + 2) 
    554         plt.plot(theory_x, mresid, '.') 
     534        plt.plot(data.x, mresid, '.') 
    555535        plt.xlabel("$q$/A$^{-1}$") 
    556536        plt.ylabel('residuals') 

  • sasmodels/direct_model.py

    r1a8c11c rd18d6dd  
    250250            qmax = getattr(data, 'qmax', np.inf) 
    251251            accuracy = getattr(data, 'accuracy', 'Low') 
    252             index = (data.mask == 0) & (q >= qmin) & (q <= qmax) 
     252            index = ~data.mask & (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) 
    268265            if data.y is not None: 
    269266                index &= ~np.isnan(data.y) 

  • sasmodels/jitter.py

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

  • sasmodels/kernel_iq.c

    r7c35fda rdc6f601  
    8383  in_spin = clip(in_spin, 0.0, 1.0); 
    8484  out_spin = clip(out_spin, 0.0, 1.0); 
    85   // Previous version of this function took the square root of the weights, 
    86   // under the assumption that  
    87   // 
     85  // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that 
    8886  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
    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 
     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 
    10092  weight[4] = weight[1]; // du.imag 
    10193  weight[5] = weight[2]; // ud.imag 

  • sasmodels/kerneldll.py

    r1a3559f rbf94e6e  
    9999    pass 
    100100# pylint: enable=unused-import 
    101  
    102 if "SAS_DLL_PATH" in os.environ: 
    103     SAS_DLL_PATH = os.environ["SAS_DLL_PATH"] 
    104 else: 
    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") 
    107101 
    108102if "SAS_COMPILER" in os.environ: 
    … …  
    129123    # add openmp support if not running on a mac 
    130124    if sys.platform != "darwin": 
    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 
     125        CC.append("-fopenmp") 
    135126    def compile_command(source, output): 
    136127        """unix compiler command""" 
    … …  
    167158        return CC + [source, "-o", output, "-lm"] 
    168159 
     160# Assume the default location of module DLLs is in .sasmodels/compiled_models. 
     161DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
     162 
    169163ALLOW_SINGLE_PRECISION_DLLS = True 
    170164 
    … …  
    203197        return path 
    204198 
    205     return joinpath(SAS_DLL_PATH, basename) 
     199    return joinpath(DLL_PATH, basename) 
    206200 
    207201 
    … …  
    212206    exist yet if it hasn't been compiled. 
    213207    """ 
    214     return os.path.join(SAS_DLL_PATH, dll_name(model_info, dtype)) 
     208    return os.path.join(DLL_PATH, dll_name(model_info, dtype)) 
    215209 
    216210 
    … …  
    231225    models are not allowed as DLLs. 
    232226 
    233     Set *sasmodels.kerneldll.SAS_DLL_PATH* to the compiled dll output path. 
    234     Alternatively, set the environment variable *SAS_DLL_PATH*. 
     227    Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 
    235228    The default is in ~/.sasmodels/compiled_models. 
    236229    """ 
    … …  
    251244    if need_recompile: 
    252245        # Make sure the DLL path exists 
    253         if not os.path.exists(SAS_DLL_PATH): 
    254             os.makedirs(SAS_DLL_PATH) 
     246        if not os.path.exists(DLL_PATH): 
     247            os.makedirs(DLL_PATH) 
    255248        basename = splitext(os.path.basename(dll))[0] + "_" 
    256249        system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 

  • sasmodels/models/core_shell_parallelepiped.c

    rdbf1a60 re077231  
    5959 
    6060    // outer integral (with gauss points), integration limits = 0, 1 
    61     // substitute d_cos_alpha for sin_alpha d_alpha 
    6261    double outer_sum = 0; //initialize integral 
    6362    for( int i=0; i<GAUSS_N; i++) { 
    6463        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    6564        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
     65 
     66        // inner integral (with gauss points), integration limits = 0, pi/2 
    6667        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
    6768        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) 
    7169        double inner_sum = 0.0; 
    7270        for(int j=0; j<GAUSS_N; j++) { 
    73             const double u = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     71            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    7472            double sin_beta, cos_beta; 
    75             SINCOS(M_PI_2*u, sin_beta, cos_beta); 
     73            SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
    7674            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
    7775            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
    … …  
    9391            inner_sum += GAUSS_W[j] * f * f; 
    9492        } 
    95         // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    9693        inner_sum *= 0.5; 
    9794        // now sum up the outer integral 
    9895        outer_sum += GAUSS_W[i] * inner_sum; 
    9996    } 
    100     // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    10197    outer_sum *= 0.5; 
    10298 

  • sasmodels/models/core_shell_parallelepiped.py

    rf89ec96 r97be877  
    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 "rim" can be 
    7 different on each (pair) of faces. The three dimensions of the core of the 
    8 parallelepiped (strictly here a cuboid) may be given in *any* size order as 
    9 long as the particles are randomly oriented (i.e. take on all possible 
    10 orientations see notes on 2D below). To avoid multiple fit solutions, 
    11 especially with Monte-Carlo fit methods, it may be advisable to restrict their 
    12 ranges. There may be a number of closely similar "best fits", so some trial and 
    13 error, or fixing of some dimensions at expected values, may help. 
     6The thickness and the scattering length density of the shell or 
     7"rim" can be different on each (pair) of faces. 
    148 
    159The form factor is normalized by the particle volume $V$ such that 
    … …  
    1711.. math:: 
    1812 
    19     I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle 
    20     + \text{background} 
     13    I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background} 
    2114 
    2215where $\langle \ldots \rangle$ is an average over all possible orientations 
    23 of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be 
    24 pulled out of the form factor term due to the multiple slds in the model. 
    25  
    26 The core of the solid is defined by the dimensions $A$, $B$, $C$ here shown 
    27 such 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  
     16of the rectangular solid. 
     17 
     18The function calculated is the form factor of the rectangular solid below. 
     19The 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 
    3423 
    3524There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 
    3625(on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$ 
    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  
    46 The total volume of the solid is thus given as 
     26$(=t_C)$ faces. The projection in the $AB$ plane is then 
     27 
     28.. image:: img/core_shell_parallelepiped_projection.jpg 
     29 
     30The volume of the solid is 
    4731 
    4832.. math:: 
    4933 
    5034    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
     35 
     36**meaning that there are "gaps" at the corners of the solid.** 
    5137 
    5238The intensity calculated follows the :ref:`parallelepiped` model, with the 
    5339core-shell intensity being calculated as the square of the sum of the 
    54 amplitudes of the core and the slabs on the edges. The scattering amplitude is 
    55 computed for a particular orientation of the core-shell parallelepiped with 
    56 respect to the scattering vector and then averaged over all possible 
    57 orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis 
    58 of the parallelepiped, and $\beta$ is the angle between the projection of the 
    59 particle 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  
    66 and 
    67  
    68 .. math:: 
    69  
    70     F(q,\alpha,\beta) 
     40amplitudes of the core and the slabs on the edges. 
     41 
     42the scattering amplitude is computed for a particular orientation of the 
     43core-shell parallelepiped with respect to the scattering vector and then 
     44averaged over all possible orientations, where $\alpha$ is the angle between 
     45the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
     46the angle between projection of the particle in the $xy$ detector plane 
     47and the $y$ axis. 
     48 
     49.. math:: 
     50 
     51    F(Q) 
    7152    &= (\rho_\text{core}-\rho_\text{solvent}) 
    7253       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
    7354    &+ (\rho_\text{A}-\rho_\text{solvent}) 
    74         \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\ 
     55        \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
    7556    &+ (\rho_\text{B}-\rho_\text{solvent}) 
    7657        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
    … …  
    8263.. math:: 
    8364 
    84     S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L} 
     65    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
    8566 
    8667and 
    … …  
    8869.. math:: 
    8970 
    90     Q_A &= q \sin\alpha \sin\beta \\ 
    91     Q_B &= q \sin\alpha \cos\beta \\ 
    92     Q_C &= q \cos\alpha 
     71    Q_A &= \sin\alpha \sin\beta \\ 
     72    Q_B &= \sin\alpha \cos\beta \\ 
     73    Q_C &= \cos\alpha 
    9374 
    9475 
    9576where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
    96 are the scattering lengths of the parallelepiped core, and the rectangular 
     77are the scattering length of the parallelepiped core, and the rectangular 
    9778slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
    9879is the scattering length of the solvent. 
    9980 
    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  
    10981FITTING NOTES 
    11082~~~~~~~~~~~~~ 
    11183 
    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. 
     84If the scale is set equal to the particle volume fraction, $\phi$, the returned 
     85value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 
     86**no interparticle interference effects are included in this calculation.** 
     87 
     88There are many parameters in this model. Hold as many fixed as possible with 
     89known values, or you will certainly end up at a solution that is unphysical. 
     90 
     91The returned value is in units of |cm^-1|, on absolute scale. 
     92 
     93NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
     94based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
     95and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
     96to give an oblate or prolate particle, to give an effective radius, 
     97for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     98 
     99For 2d data the orientation of the particle is required, described using 
     100angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
     101details of the calculation and angular dispersions see :ref:`orientation`. 
     102The 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 
     105For 2d, constraints must be applied during fitting to ensure that the 
     106inequality $A < B < C$ is not violated, and hence the correct definition 
     107of angles is preserved. The calculation will not report an error, 
     108but the results may be not correct. 
    143109 
    144110.. figure:: img/parallelepiped_angle_definition.png 
    145111 
    146112    Definition of the angles for oriented core-shell parallelepipeds. 
    147     Note that rotation $\theta$, initially in the $x-z$ plane, is carried 
     113    Note that rotation $\theta$, initially in the $xz$ plane, is carried 
    148114    out first, then rotation $\phi$ about the $z$ axis, finally rotation 
    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. 
     115    $\Psi$ is now around the axis of the cylinder. The neutron or X-ray 
     116    beam is along the $z$ axis. 
    151117 
    152118.. figure:: img/parallelepiped_angle_projection.png 
    … …  
    154120    Examples of the angles for oriented core-shell parallelepipeds against the 
    155121    detector plane. 
    156  
    157  
    158 Validation 
    159 ---------- 
    160  
    161 Cross-checked against hollow rectangular prism and rectangular prism for equal 
    162 thickness overlapping sides, and by Monte Carlo sampling of points within the 
    163 shape for non-uniform, non-overlapping sides. 
    164  
    165122 
    166123References 
    … …  
    178135 
    179136* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    180 * **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016 
     137* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
    181138* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 
    182 * **Last Reviewed by:** Paul Butler **Date:** May 24, 2018 - documentation 
    183   updated 
     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. 
    184142""" 
    185143 

  • sasmodels/models/core_shell_sphere.py

    rdc76240 r2d81cfe  
    2121.. math:: 
    2222 
    23     F(q) = \frac{3}{V_s}\left[ 
     23    F^2(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

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

  • sasmodels/models/mass_surface_fractal.py

    r7994359 r2d81cfe  
    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     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      
     41    $(surface\_dim + mass\_dim ) < 6$ . 
     42 
    4943 
    5044References 
    5145---------- 
    5246 
    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) 
     47P Schmidt, *J Appl. Cryst.*, 24 (1991) 414-435 Equation(19) 
    5648 
    57 Authorship 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 
     49A J Hurd, D W Schaefer, J E Martin, *Phys. Rev. A*, 
     5035 (1987) 2361-2364 Equation(2) 
    6251""" 
    6352 
    … …  
    7867        rg_primary    =  rg 
    7968        background   =  background 
     69        Ref: Schmidt, J Appl Cryst, eq(19), (1991), 24, 414-435 
    8070        Hurd, Schaefer, Martin, Phys Rev A, eq(2),(1987),35, 2361-2364 
    8171        Note that 0 < Ds< 6 and 0 < Dm < 6. 
    … …  
    8878    ["fractal_dim_mass", "",      1.8, [0.0, 6.0], "", "Mass fractal dimension"], 
    8979    ["fractal_dim_surf", "",      2.3, [0.0, 6.0], "", "Surface fractal dimension"], 
    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"], 
     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"], 
    9282] 
    9383# pylint: enable=bad-whitespace, line-too-long 
    … …  
    117107            fractal_dim_mass=1.8, 
    118108            fractal_dim_surf=2.3, 
    119             rg_cluster=4000.0, 
    120             rg_primary=86.7) 
     109            rg_cluster=86.7, 
     110            rg_primary=4000.0) 
    121111 
    122112tests = [ 
    123113 
    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, 
     114    # Accuracy tests based on content in test/utest_other_models.py 
     115    [{'fractal_dim_mass':      1.8, 
    126116      'fractal_dim_surf':   2.3, 
    127117      'rg_cluster':   86.7, 
    … …  
    133123    [{'fractal_dim_mass':      3.3, 
    134124      'fractal_dim_surf':   1.0, 
    135       'rg_cluster': 4000.0, 
    136       'rg_primary':   90.0, 
    137      }, 0.001, 0.0932516614456], 
     125      'rg_cluster':   90.0, 
     126      'rg_primary': 4000.0, 
     127     }, 0.001, 0.18562699016], 
    138128 
    139129    [{'fractal_dim_mass':      1.3, 
    140       'fractal_dim_surf':   2.0, 
    141       'rg_cluster': 2000.0, 
    142       'rg_primary':   90.0, 
     130      'fractal_dim_surf':   1.0, 
     131      'rg_cluster':   90.0, 
     132      'rg_primary': 2000.0, 
    143133      'background':    0.8, 
    144      }, 0.001, 1.28296431786], 
     134     }, 0.001, 1.16539753641], 
    145135 
    146136    [{'fractal_dim_mass':      2.3, 
    147       'fractal_dim_surf':   3.1, 
    148       'rg_cluster':  1000.0, 
    149       'rg_primary':  30.0, 
     137      'fractal_dim_surf':   1.0, 
     138      'rg_cluster':   90.0, 
     139      'rg_primary': 1000.0, 
    150140      'scale':        10.0, 
    151141      'background':    0.0, 
    152      }, 0.051, 0.00333804044899], 
     142     }, 0.051, 0.000169548800377], 
    153143    ] 

  • sasmodels/models/parallelepiped.c

    rdbf1a60 r108e70e  
    3838            inner_total += GAUSS_W[j] * square(si1 * si2); 
    3939        } 
    40         // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    4140        inner_total *= 0.5; 
    4241 
    … …  
    4443        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4544    } 
    46     // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    4745    outer_total *= 0.5; 
    4846 

  • sasmodels/models/parallelepiped.py

    rf89ec96 ref07e95  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
     4The form factor is normalized by the particle volume. 
     5For information about polarised and magnetic scattering, see 
     6the :ref:`magnetism` documentation. 
     7 
    48Definition 
    59---------- 
    610 
    7 This model calculates the scattering from a rectangular solid 
    8 (:numref:`parallelepiped-image`). 
    9 If you need to apply polydispersity, see also :ref:`rectangular-prism`. For 
    10 information about polarised and magnetic scattering, see 
    11 the :ref:`magnetism` documentation. 
     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`. 
    1214 
    1315.. _parallelepiped-image: 
    … …  
    1921 
    2022The three dimensions of the parallelepiped (strictly here a cuboid) may be 
    21 given in *any* size order as long as the particles are randomly oriented (i.e. 
    22 take on all possible orientations see notes on 2D below). To avoid multiple fit 
    23 solutions, especially with Monte-Carlo fit methods, it may be advisable to 
    24 restrict their ranges. There may be a number of closely similar "best fits", so 
    25 some trial and error, or fixing of some dimensions at expected values, may 
    26 help. 
    27  
    28 The form factor is normalized by the particle volume and the 1D scattering 
    29 intensity $I(q)$ is then calculated as: 
     23given in *any* size order. To avoid multiple fit solutions, especially 
     24with Monte-Carlo fit methods, it may be advisable to restrict their ranges. 
     25There may be a number of closely similar "best fits", so some trial and 
     26error, or fixing of some dimensions at expected values, may help. 
     27 
     28The 1D scattering intensity $I(q)$ is calculated as: 
    3029 
    3130.. Comment by Miguel Gonzalez: 
    … …  
    4039 
    4140    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 
    42            \left< P(q, \alpha, \beta) \right> + \text{background} 
     41           \left< P(q, \alpha) \right> + \text{background} 
    4342 
    4443where the volume $V = A B C$, the contrast is defined as 
    45 $\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$ 
    46 is the form factor corresponding to a parallelepiped oriented 
    47 at 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 
    49 and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all 
    50 orientations. 
     44$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, 
     45$P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented 
     46at an angle $\alpha$ (angle between the long axis C and $\vec q$), 
     47and the averaging $\left<\ldots\right>$ is applied over all orientations. 
    5148 
    5249Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
    53 form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_) 
     50form factor is given by (Mittelbach and Porod, 1961) 
    5451 
    5552.. math:: 
    … …  
    6966    \mu &= qB 
    7067 
    71 where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been 
    72 applied. 
    73  
    74 For **oriented** particles, the 2D scattering intensity, $I(q_x, q_y)$, is 
    75 given 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) 
     68The scattering intensity per unit volume is returned in units of |cm^-1|. 
     69 
     70NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
     71the averaged effective radius, after appropriately sorting the three 
     72dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and 
     73length $(= C)$ values, and used as the effective radius for 
     74$S(q)$ when $P(q) \cdot S(q)$ is applied. 
     75 
     76For 2d data the orientation of the particle is required, described using 
     77angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     78of 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 
     89The angle $\Psi$ is the rotational angle around the $C$ axis. 
     90For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis 
     91oriented parallel to the y-axis of the detector with $A$ along the x-axis. 
     92For 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, 
     94before doing a final rotation of $\Psi$ degrees around the resulting $C$ axis 
     95of 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 
     108On introducing "Orientational Distribution" in the angles, "distribution of 
     109theta" and "distribution of phi" parameters will appear. These are actually 
     110rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, 
     111perpendicular to the $a$ x $c$ and $b$ x $c$ faces. (When $\theta = \phi = 0$ 
     112these are parallel to the $Y$ and $X$ axes of the instrument.) The third 
     113orientation distribution, in $\psi$, is about the $c$ axis of the particle, 
     114perpendicular to the $a$ x $b$ face. Some experimentation may be required to 
     115understand the 2d patterns fully as discussed in :ref:`orientation` . 
     116 
     117For a given orientation of the parallelepiped, the 2D form factor is 
     118calculated 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 
     126with 
     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 
     134and 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) 
    80139            + \text{background} 
    81140 
    … …  
    89148   with scale being the volume fraction. 
    90149 
    91 Where $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  
    102 with 
    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  
    111 FITTING 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  
    168150 
    169151Validation 
    … …  
    174156angles. 
    175157 
     158 
    176159References 
    177160---------- 
    178161 
    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 
     162P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 
     163 
     164R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
    182165 
    183166Authorship and Verification 
    … …  
    186169* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    187170* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017 
    188 * **Last Reviewed by:**  Miguel Gonzales and Paul Butler **Date:** May 24, 
    189   2018 - documentation updated 
     171* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017 
    190172""" 
    191173 

  • sasmodels/resolution.py

    r0b9c6df r2d81cfe  
    2020MINIMUM_RESOLUTION = 1e-8 
    2121MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data 
    22 PINHOLE_N_SIGMA = 2.5 # From: Barker & Pedersen 1995 JAC 
    2322 
    2423class Resolution(object): 
    … …  
    6665    *q_calc* is the list of points to calculate, or None if this should 
    6766    be estimated from the *q* and *q_width*. 
    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): 
     67    """ 
     68    def __init__(self, q, q_width, q_calc=None, nsigma=3): 
    7369        #*min_step* is the minimum point spacing to use when computing the 
    7470        #underlying model.  It should be on the order of 
    … …  
    8682 
    8783        # Protect against models which are not defined for very low q.  Limit 
    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. 
     84        # the smallest q value evaluated (in absolute) to 0.02*min 
    9285        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
    93         self.q_calc = self.q_calc[self.q_calc >= cutoff] 
     86        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff] 
    9487 
    9588        # Build weight matrix from calculated q values 
    9689        self.weight_matrix = pinhole_resolution( 
    97             self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION), 
    98             nsigma=nsigma) 
     90            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) 
     91        self.q_calc = abs(self.q_calc) 
    9992 
    10093    def apply(self, theory): 
    … …  
    108101    *q* points at which the data is measured. 
    109102 
    110     *qx_width* slit width in qx 
    111  
    112     *qy_width* slit height in qy 
     103    *dqx* slit width in qx 
     104 
     105    *dqy* slit height in qy 
    113106 
    114107    *q_calc* is the list of points to calculate, or None if this should 
    … …  
    161154 
    162155 
    163 def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA): 
    164     r""" 
     156def pinhole_resolution(q_calc, q, q_width): 
     157    """ 
    165158    Compute the convolution matrix *W* for pinhole resolution 1-D data. 
    166159 
    … …  
    169162    *W*, the resolution smearing can be computed using *dot(W,q)*. 
    170163 
    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  
    177164    *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. 
    183165    """ 
    184166    # The current algorithm is a midpoint rectangle rule.  In the test case, 
    … …  
    188170    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    189171    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. 
    196172    weights /= np.sum(weights, axis=0)[None, :] 
    197173    return weights 
    … …  
    518494 
    519495 
    520 def gaussian(q, q0, dq, nsigma=2.5): 
    521     """ 
    522     Return the truncated Gaussian resolution function. 
     496def gaussian(q, q0, dq): 
     497    """ 
     498    Return the Gaussian resolution function. 
    523499 
    524500    *q0* is the center, *dq* is the width and *q* are the points to evaluate. 
    525501    """ 
    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) 
     502    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 
    530503 
    531504 
    … …  
    585558 
    586559 
    587 def romberg_pinhole_1d(q, q_width, form, pars, nsigma=2.5): 
     560def romberg_pinhole_1d(q, q_width, form, pars, nsigma=5): 
    588561    """ 
    589562    Romberg integration for pinhole resolution. 
    … …  
    705678        np.testing.assert_equal(output, self.y) 
    706679 
    707     # TODO: turn pinhole/slit demos into tests 
    708  
    709     @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    710680    def test_pinhole(self): 
    711681        """ 
    … …  
    716686        theory = 12.0-1000.0*resolution.q_calc 
    717687        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. 
    721688        answer = [ 
    722             10.47037734, 9.86925860, 
    723             9., 8., 7., 6., 5., 4., 
    724             3.13074140, 2.52962266, 
     689            10.44785079, 9.84991299, 8.98101708, 
     690            7.99906585, 6.99998311, 6.00001689, 
     691            5.00093415, 4.01898292, 3.15008701, 2.55214921, 
    725692            ] 
    726693        np.testing.assert_allclose(output, answer, atol=1e-8) 
    … …  
    765732        self._compare(q, output, answer, 1e-6) 
    766733 
    767     @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    768734    def test_pinhole(self): 
    769735        """ 
    … …  
    780746        self._compare(q, output, answer, 3e-4) 
    781747 
    782     @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    783748    def test_pinhole_romberg(self): 
    784749        """ 
    … …  
    796761        #                     2*np.pi/pars['radius']/200) 
    797762        #tol = 0.001 
    798         ## The default 2.5 sigma and no extra points gets 1% 
     763        ## The default 3 sigma and no extra points gets 1% 
    799764        q_calc = None  # type: np.ndarray 
    800765        tol = 0.01 
    … …  
    11151080 
    11161081    if isinstance(resolution, Slit1D): 
    1117         width, height = resolution.qx_width, resolution.qy_width 
     1082        width, height = resolution.dqx, resolution.dqy 
    11181083        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars) 
    11191084    else: 

  • sasmodels/rst2html.py

    r1fbadb2 r2d81cfe  
    5656    # others don't work properly with math_output! 
    5757    if math_output == "mathjax": 
    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} 
     58        settings = {"math_output": math_output} 
    6159    else: 
    6260        settings = {"math-output": math_output} 

  • sasmodels/sasview_model.py

    rd533590 r05df1de  
    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) -> None 
     721        # type: (str, weights.Dispersion) -> Dict[str, Any] 
    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" 
    745                              % parameter) 
     744            raise ValueError("%r is not a dispersity or orientation parameter") 
    746745 
    747746    def _dispersion_mesh(self): 
Note: See TracChangeset for help on using the changeset viewer.