Changeset 6cff473 in sasmodels


Ignore:
Timestamp:
Apr 4, 2017 11:31:53 AM (8 years ago)
Author:
GitHub <noreply@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
f68e2a5
Parents:
cb0dc22 (diff), a769b54 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Andrew Jackson <andrew.jackson@…> (04/04/17 11:31:53)
git-committer:
GitHub <noreply@…> (04/04/17 11:31:53)
Message:

Merge pull request #29 from SasView?/sascomp_with_data

no ticket: allow sascomp with with q,dq from data file

Files:
1 added
45 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    rf72d70a r6cff473  
    3030 
    3131import sys 
     32import os 
    3233import math 
    3334import datetime 
    … …  
    4041from . import kerneldll 
    4142from . import exception 
    42 from .data import plot_theory, empty_data1D, empty_data2D 
     43from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4344from .direct_model import DirectModel 
    4445from .convert import revert_name, revert_pars, constrain_new_to_old 
    … …  
    8586    -html shows the model docs instead of running the model 
    8687    -title="note" adds note to the plot title, after the model name 
     88    -data="path" uses q, dq from the data file 
    8789 
    8890Any two calculation engines can be selected for comparison: 
    … …  
    747749    comp = opts['engines'][1] if have_comp else None 
    748750    data = opts['data'] 
     751    use_data = have_base ^ have_comp 
    749752 
    750753    # Plot if requested 
    … …  
    763766    if have_base: 
    764767        if have_comp: plt.subplot(131) 
    765         plot_theory(data, base_value, view=view, use_data=False, limits=limits) 
     768        plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
    766769        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    767770        #cbar_title = "log I" 
    … …  
    769772        if have_base: plt.subplot(132) 
    770773        if not opts['is2d'] and have_base: 
    771             plot_theory(data, base_value, view=view, use_data=False, limits=limits) 
    772         plot_theory(data, comp_value, view=view, use_data=False, limits=limits) 
     774            plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
     775        plot_theory(data, comp_value, view=view, use_data=use_data, limits=limits) 
    773776        plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 
    774777        #cbar_title = "log I" 
    … …  
    784787            err[err>cutoff] = cutoff 
    785788        #err,errstr = base/comp,"ratio" 
    786         plot_theory(data, None, resid=err, view=errview, use_data=False) 
     789        plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
    787790        if view == 'linear': 
    788791            plt.xscale('linear') 
    … …  
    834837VALUE_OPTIONS = [ 
    835838    # Note: random is both a name option and a value option 
    836     'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 
     839    'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data', 
    837840    ] 
    838841 
    … …  
    951954        'html'      : False, 
    952955        'title'     : None, 
     956        'data'      : None, 
    953957    } 
    954958    engines = [] 
    … …  
    971975        elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:]) 
    972976        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:]) 
    973         elif arg.startswith('-title'):     opts['title'] = arg[7:] 
     977        elif arg.startswith('-title='):    opts['title'] = arg[7:] 
     978        elif arg.startswith('-data='):     opts['data'] = arg[6:] 
    974979        elif arg == '-random':  opts['seed'] = np.random.randint(1000000) 
    975980        elif arg == '-preset':  opts['seed'] = -1 
    … …  
    11121117 
    11131118    # Create the computational engines 
    1114     data, _ = make_data(opts) 
     1119    if opts['data'] is not None: 
     1120        data = load_data(os.path.expanduser(opts['data'])) 
     1121    else: 
     1122        data, _ = make_data(opts) 
    11151123    if n1: 
    11161124        base = make_engine(model_info, data, engines[0], opts['cutoff']) 

  • sasmodels/data.py

    r40a87fa ra769b54  
    5454    if data is None: 
    5555        raise IOError("Data %r could not be loaded" % filename) 
     56    if hasattr(data, 'x'): 
     57        data.qmin, data.qmax = data.x.min(), data.x.max() 
     58        data.mask = (np.isnan(data.y) if data.y is not None 
     59                     else np.zeros_like(data.x, dtype='bool')) 
    5660    return data 
    5761 
    … …  
    348352    # data, but they already handle the masking and graph markup already, so 
    349353    # do not repeat. 
    350     if hasattr(data, 'lam'): 
     354    if hasattr(data, 'isSesans') and data.isSesans: 
    351355        _plot_result_sesans(data, None, None, use_data=True, limits=limits) 
    352356    elif hasattr(data, 'qx_data'): 
    … …  
    376380    *Iq_calc* is the raw theory values without resolution smearing 
    377381    """ 
    378     if hasattr(data, 'lam'): 
     382    if hasattr(data, 'isSesans') and data.isSesans: 
    379383        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 
    380384    elif hasattr(data, 'qx_data'): 

  • sasmodels/direct_model.py

    rb397165 ra769b54  
    192192 
    193193        # interpret data 
    194         if hasattr(data, 'lam'): 
     194        if hasattr(data, 'isSesans') and data.isSesans: 
    195195            self.data_type = 'sesans' 
    196196        elif hasattr(data, 'qx_data'): 

  • doc/ref/index.rst

    rc34a31f r9f12fbe  
    1010   refs.rst 
    1111   gpu/gpu_computations.rst 
     12   gpu/opencl_installation.rst 
    1213   magnetism/magnetism.rst 
    1314   sesans/sans_to_sesans.rst 

  • example/fit.py

    r1182da5 rf4b36fa  
    2424            % section) 
    2525data = radial_data if section != "tangential" else tan_data 
    26 phi = 0 if section != "tangential" else 90 
     26theta = 89.9 if section != "tangential" else 0 
     27phi = 90 
    2728kernel = load_model(name, dtype="single") 
    2829cutoff = 1e-3 
    … …  
    3031if name == "ellipsoid": 
    3132    model = Model(kernel, 
    32         scale=0.08, 
    33         r_polar=15, r_equatorial=800, 
     33        scale=0.08, background=35, 
     34        radius_polar=15, radius_equatorial=800, 
    3435        sld=.291, sld_solvent=7.105, 
    35         background=0, 
    36         theta=90, phi=phi, 
    37         theta_pd=15, theta_pd_n=40, theta_pd_nsigma=3, 
    38         r_polar_pd=0.222296, r_polar_pd_n=1, r_polar_pd_nsigma=0, 
    39         r_equatorial_pd=.000128, r_equatorial_pd_n=1, r_equatorial_pd_nsigma=0, 
     36        theta=theta, phi=phi, 
     37        theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3, 
    4038        phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3, 
     39        radius_polar_pd=0.222296, radius_polar_pd_n=1, radius_polar_pd_nsigma=0, 
     40        radius_equatorial_pd=.000128, radius_equatorial_pd_n=1, radius_equatorial_pd_nsigma=0, 
    4141        ) 
    4242 
    43  
    4443    # SET THE FITTING PARAMETERS 
    45     model.r_polar.range(15, 1000) 
    46     model.r_equatorial.range(15, 1000) 
    47     model.theta_pd.range(0, 360) 
     44    model.radius_polar.range(15, 1000) 
     45    model.radius_equatorial.range(15, 1000) 
     46    #model.theta.range(0, 90) 
     47    #model.theta_pd.range(0,10) 
     48    model.phi_pd.range(0,20) 
     49    model.phi.range(0, 180) 
    4850    model.background.range(0,1000) 
    4951    model.scale.range(0, 10) 
    5052 
    5153 
    52  
    5354elif name == "lamellar": 
    5455    model = Model(kernel, 
    55         scale=0.08, 
     56        scale=0.08, background=0.003, 
    5657        thickness=19.2946, 
    5758        sld=5.38,sld_sol=7.105, 
    58         background=0.003, 
    5959        thickness_pd= 0.37765, thickness_pd_n=10, thickness_pd_nsigma=3, 
    6060        ) 
    61  
    6261 
    6362    # SET THE FITTING PARAMETERS 
    … …  
    7776        radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3, 
    7877        length_pd=0.493, length_pd_n=10, length_pd_nsigma=3, 
    79         phi_pd=0, phi_pd_n=5, phi_pd_nsigma=3,) 
     78        phi_pd=0, phi_pd_n=5 phi_pd_nsigma=3,) 
    8079        """ 
    8180    pars = dict( 
    8281        scale=.01, background=35, 
    8382        sld=.291, sld_solvent=5.77, 
    84         radius=250, length=178,  
    85         theta=90, phi=phi, 
     83        radius=250, length=178, 
    8684        radius_pd=0.1, radius_pd_n=5, radius_pd_nsigma=3, 
    8785        length_pd=0.1,length_pd_n=5, length_pd_nsigma=3, 
    88         theta_pd=10, theta_pd_n=50, theta_pd_nsigma=3, 
    89         phi_pd=0, phi_pd_n=10, phi_pd_nsigma=3) 
     86        theta=theta, phi=phi, 
     87        theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3, 
     88        phi_pd=10, phi_pd_n=20, phi_pd_nsigma=3) 
    9089    model = Model(kernel, **pars) 
    9190 
    … …  
    9392    model.radius.range(1, 500) 
    9493    model.length.range(1, 5000) 
    95     model.theta.range(-90,100) 
    96     model.theta_pd.range(0, 30) 
    97     model.theta_pd_n = model.theta_pd + 5 
     94    #model.theta.range(0, 90) 
     95    model.phi.range(0, 180) 
     96    model.phi_pd.range(0, 30) 
    9897    model.radius_pd.range(0, 1) 
    99     model.length_pd.range(0, 2) 
     98    model.length_pd.range(0, 1) 
    10099    model.scale.range(0, 10) 
    101100    model.background.range(0, 100) 
    … …  
    104103elif name == "core_shell_cylinder": 
    105104    model = Model(kernel, 
    106         scale= .031, radius=19.5, thickness=30, length=22, 
    107         sld_core=7.105, sld_shell=.291, sdl_solvent=7.105, 
    108         background=0, theta=0, phi=phi, 
    109  
     105        scale= .031, background=0, 
     106        radius=19.5, thickness=30, length=22, 
     107        sld_core=7.105, sld_shell=.291, sld_solvent=7.105, 
    110108        radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3, 
    111109        length_pd=0.26, length_pd_n=10, length_pd_nsigma=3, 
    112110        thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=1, 
    113         theta_pd=1, theta_pd_n=10, theta_pd_nsigma=3, 
    114         phi_pd=0.1, phi_pd_n=1, phi_pd_nsigma=1, 
     111        theta=theta, phi=phi, 
     112        theta_pd=1, theta_pd_n=1, theta_pd_nsigma=3, 
     113        phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3, 
    115114        ) 
    116115 
    117116    # SET THE FITTING PARAMETERS 
    118     #model.radius.range(115, 1000) 
    119     #model.length.range(0, 2500) 
     117    model.radius.range(115, 1000) 
     118    model.length.range(0, 2500) 
    120119    #model.thickness.range(18, 38) 
    121120    #model.thickness_pd.range(0, 1) 
    122121    #model.phi.range(0, 90) 
     122    model.phi_pd.range(0,20) 
    123123    #model.radius_pd.range(0, 1) 
    124124    #model.length_pd.range(0, 1) 
    … …  
    131131elif name == "capped_cylinder": 
    132132    model = Model(kernel, 
    133         scale=.08, radius=20, cap_radius=40, length=400, 
     133        scale=.08, background=35, 
     134        radius=20, cap_radius=40, length=400, 
    134135        sld=1, sld_solvent=6.3, 
    135         background=0, theta=0, phi=phi, 
    136136        radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3, 
    137137        cap_radius_pd=.1, cap_radius_pd_n=5, cap_radius_pd_nsigma=3, 
    138138        length_pd=.1, length_pd_n=1, length_pd_nsigma=0, 
    139         theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0, 
    140         phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
     139        theta=theta, phi=phi, 
     140        theta_pd=0, theta_pd_n=1, theta_pd_nsigma=0, 
     141        phi_pd=10, phi_pd_n=20, phi_pd_nsigma=0, 
    141142        ) 
    142143 
     144    model.radius.range(115, 1000) 
     145    model.length.range(0, 2500) 
     146    #model.thickness.range(18, 38) 
     147    #model.thickness_pd.range(0, 1) 
     148    #model.phi.range(0, 90) 
     149    model.phi_pd.range(0,20) 
     150    #model.radius_pd.range(0, 1) 
     151    #model.length_pd.range(0, 1) 
     152    #model.theta_pd.range(0, 360) 
     153    #model.background.range(0,5) 
    143154    model.scale.range(0, 1) 
    144155 
    … …  
    146157elif name == "triaxial_ellipsoid": 
    147158    model = Model(kernel, 
    148         scale=0.08, req_minor=15, req_major=20, rpolar=500, 
     159        scale=0.08, background=35, 
     160        radius_equat_minor=15, radius_equat_major=20, radius_polar=500, 
    149161        sld=7.105, solvent_sld=.291, 
    150         background=5, theta=0, phi=phi, psi=0, 
     162        radius_equat_minor_pd=.1, radius_equat_minor_pd_n=1, radius_equat_minor_pd_nsigma=0, 
     163        radius_equat_major_pd=.1, radius_equat_major_pd_n=1, radius_equat_major_pd_nsigma=0, 
     164        radius_polar_pd=.1, radius_polar_pd_n=1, radius_polar_pd_nsigma=0, 
     165        theta=theta, phi=phi, psi=0, 
    151166        theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3, 
    152167        phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    153168        psi_pd=30, psi_pd_n=1, psi_pd_nsigma=0, 
    154         req_minor_pd=.1, req_minor_pd_n=1, req_minor_pd_nsigma=0, 
    155         req_major_pd=.1, req_major_pd_n=1, req_major_pd_nsigma=0, 
    156         rpolar_pd=.1, rpolar_pd_n=1, rpolar_pd_nsigma=0, 
    157169        ) 
    158170 
    159171    # SET THE FITTING PARAMETERS 
    160     model.req_minor.range(15, 1000) 
    161     model.req_major.range(15, 1000) 
    162     #model.rpolar.range(15, 1000) 
     172    model.radius_equat_minor.range(15, 1000) 
     173    model.radius_equat_major.range(15, 1000) 
     174    #model.radius_polar.range(15, 1000) 
    163175    #model.background.range(0,1000) 
    164176    #model.theta_pd.range(0, 360) 
    … …  
    173185M = Experiment(data=data, model=model) 
    174186if section == "both": 
    175    tan_model = Model(model.kernel, **model.parameters()) 
     187   tan_model = Model(model.sasmodel, **model.parameters()) 
    176188   tan_model.phi = model.phi - 90 
    177189   tan_model.cutoff = cutoff 

  • sasmodels/compare_many.py

    r424fe00 rf72d70a  
    106106    header = ('\n"Model","%s","Count","%d","Dimension","%s"' 
    107107              % (name, N, "2D" if is_2d else "1D")) 
    108     if not mono: header += ',"Cutoff",%g'%(cutoff,) 
     108    if not mono: 
     109        header += ',"Cutoff",%g'%(cutoff,) 
    109110    print(header) 
    110111 
    … …  
    161162    max_diff = [0] 
    162163    for k in range(N): 
    163         print("%s %d"%(name, k), file=sys.stderr) 
     164        print("Model %s %d"%(name, k+1), file=sys.stderr) 
    164165        seed = np.random.randint(1e6) 
    165166        pars_i = randomize_pars(model_info, pars, seed) 
    166167        constrain_pars(model_info, pars_i) 
    167         constrain_new_to_old(model_info, pars_i) 
     168        if 'sasview' in (base, comp): 
     169            constrain_new_to_old(model_info, pars_i) 
    168170        if mono: 
    169171            pars_i = suppress_pd(pars_i) 
    … …  
    187189    Print the command usage string. 
    188190    """ 
    189     print("usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono) (single|double|quad)") 
     191    print("usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono) (single|double|quad)", 
     192          file=sys.stderr) 
    190193 
    191194 
    … …  
    204207    print("""\ 
    205208 
    206 MODEL is the model name of the model or "all" for all the models 
    207 in alphabetical order. 
     209MODEL is the model name of the model or one of the model types listed in 
     210sasmodels.core.list_models (all, py, c, double, single, opencl, 1d, 2d, 
     211nonmagnetic, magnetic).  Model types can be combined, such as 2d+single. 
    208212 
    209213COUNT is the number of randomly generated parameter sets to try. A value 
    … …  
    220224below the cutoff will be ignored.  Use "mono" for monodisperse models.  The 
    221225choice of polydisperse parameters, and the number of points in the distribution 
    222 is set in compare.py defaults for each model. 
     226is set in compare.py defaults for each model.  Polydispersity is given in the 
     227"demo" attribute of each model. 
    223228 
    224229PRECISION is the floating point precision to use for comparisons.  If two 
    225 precisions are given, then compare one to the other, ignoring sasview. 
     230precisions are given, then compare one to the other.  Precision is one of 
     231fast, single, double for GPU or single!, double!, quad! for DLL.  If no 
     232precision is given, then use single and double! respectively. 
    226233 
    227234Available models: 
    … …  
    233240    Main program. 
    234241    """ 
    235     if len(argv) not in (5, 6): 
     242    if len(argv) not in (3, 4, 5, 6): 
    236243        print_help() 
    237244        return 
    238245 
    239     model = argv[0] 
    240     if not (model in MODELS) and (model != "all"): 
    241         print('Bad model %s.  Use "all" or one of:'%model) 
     246    target = argv[0] 
     247    try: 
     248        model_list = [target] if target in MODELS else core.list_models(target) 
     249    except ValueError: 
     250        print('Bad model %s.  Use model type or one of:' % target, file=sys.stderr) 
    242251        print_models() 
     252        print('model types: all, py, c, double, single, opencl, 1d, 2d, nonmagnetic, magnetic') 
    243253        return 
    244254    try: 
    … …  
    247257        assert argv[2][1] == 'd' 
    248258        Nq = int(argv[2][2:]) 
    249         mono = argv[3] == 'mono' 
     259        mono = len(argv) <= 3 or argv[3] == 'mono' 
    250260        cutoff = float(argv[3]) if not mono else 0 
    251         base = argv[4] 
    252         comp = argv[5] if len(argv) > 5 else "sasview" 
     261        base = argv[4] if len(argv) > 4 else "single" 
     262        comp = argv[5] if len(argv) > 5 else "double!" 
    253263    except Exception: 
    254264        traceback.print_exc() 
    … …  
    258268    data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 
    259269                             'accuracy': 'Low', 'view':'log', 'zero': False}) 
    260     model_list = [model] if model != "all" else MODELS 
    261270    for model in model_list: 
    262271        compare_instance(model, data, index, N=count, mono=mono, 

  • sasmodels/conversion_table.py

    r790b16bb r0c2da4b  
    549549            "radius": "core_radius", 
    550550            "sld_solvent": "core_sld", 
    551             "n_pairs": "n_pairs", 
     551            "n_shells": "n_pairs", 
    552552            "thick_shell": "s_thickness", 
    553553            "sld": "shell_sld", 

  • sasmodels/core.py

    r52e9a45 r5124c969  
    6969        * magnetic: models with an sld 
    7070        * nommagnetic: models without an sld 
    71     """ 
    72     if kind and kind not in KINDS: 
     71 
     72    For multiple conditions, combine with plus.  For example, *c+single+2d* 
     73    would return all oriented models implemented in C which can be computed 
     74    accurately with single precision arithmetic. 
     75    """ 
     76    if kind and any(k not in KINDS for k in kind.split('+')): 
    7377        raise ValueError("kind not in " + ", ".join(KINDS)) 
    7478    files = sorted(glob(joinpath(generate.MODEL_PATH, "[a-zA-Z]*.py"))) 
    7579    available_models = [basename(f)[:-3] for f in files] 
    76     selected = [name for name in available_models if _matches(name, kind)] 
     80    if kind and '+' in kind: 
     81        all_kinds = kind.split('+') 
     82        condition = lambda name: all(_matches(name, k) for k in all_kinds) 
     83    else: 
     84        condition = lambda name: _matches(name, kind) 
     85    selected = [name for name in available_models if condition(name)] 
    7786 
    7887    return selected 

  • sasmodels/model_test.py

    r479d0f3 r598857b  
    9797        is_py = callable(model_info.Iq) 
    9898 
     99        # Some OpenCL drivers seem to be flaky, and are not producing the 
     100        # expected result.  Since we don't have known test values yet for 
     101        # all of our models, we are instead going to compare the results 
     102        # for the 'smoke test' (that is, evaluation at q=0.1 for the default 
     103        # parameters just to see that the model runs to completion) between 
     104        # the OpenCL and the DLL.  To do this, we define a 'stash' which is 
     105        # shared between OpenCL and DLL tests.  This is just a list.  If the 
     106        # list is empty (which it will be when DLL runs, if the DLL runs 
     107        # first), then the results are appended to the list.  If the list 
     108        # is not empty (which it will be when OpenCL runs second), the results 
     109        # are compared to the results stored in the first element of the list. 
     110        # This is a horrible stateful hack which only makes sense because the 
     111        # test suite is thrown away after being run once. 
     112        stash = [] 
     113 
    99114        if is_py:  # kernel implemented in python 
    100115            test_name = "Model: %s, Kernel: python"%model_name 
    … …  
    103118                                 test_method_name, 
    104119                                 platform="dll",  # so that 
    105                                  dtype="double") 
     120                                 dtype="double", 
     121                                 stash=stash) 
    106122            suite.addTest(test) 
    107123        else:   # kernel implemented in C 
     124 
     125            # test using dll if desired 
     126            if 'dll' in loaders or not core.HAVE_OPENCL: 
     127                test_name = "Model: %s, Kernel: dll"%model_name 
     128                test_method_name = "test_%s_dll" % model_info.id 
     129                test = ModelTestCase(test_name, model_info, 
     130                                     test_method_name, 
     131                                     platform="dll", 
     132                                     dtype="double", 
     133                                     stash=stash) 
     134                suite.addTest(test) 
     135 
    108136            # test using opencl if desired and available 
    109137            if 'opencl' in loaders and core.HAVE_OPENCL: 
    … …  
    116144                test = ModelTestCase(test_name, model_info, 
    117145                                     test_method_name, 
    118                                      platform="ocl", dtype=None) 
     146                                     platform="ocl", dtype=None, 
     147                                     stash=stash) 
    119148                #print("defining", test_name) 
    120                 suite.addTest(test) 
    121  
    122             # test using dll if desired 
    123             if 'dll' in loaders or not core.HAVE_OPENCL: 
    124                 test_name = "Model: %s, Kernel: dll"%model_name 
    125                 test_method_name = "test_%s_dll" % model_info.id 
    126                 test = ModelTestCase(test_name, model_info, 
    127                                      test_method_name, 
    128                                      platform="dll", 
    129                                      dtype="double") 
    130149                suite.addTest(test) 
    131150 
    … …  
    144163        """ 
    145164        def __init__(self, test_name, model_info, test_method_name, 
    146                      platform, dtype): 
    147             # type: (str, ModelInfo, str, str, DType) -> None 
     165                     platform, dtype, stash): 
     166            # type: (str, ModelInfo, str, str, DType, List[Any]) -> None 
    148167            self.test_name = test_name 
    149168            self.info = model_info 
    150169            self.platform = platform 
    151170            self.dtype = dtype 
     171            self.stash = stash  # container for the results of the first run 
    152172 
    153173            setattr(self, test_method_name, self.run_all) 
    … …  
    167187                #({}, (0.0, 0.0), None), 
    168188                # test vector form 
    169                 ({}, [0.1]*2, [None]*2), 
     189                ({}, [0.001, 0.01, 0.1], [None]*3), 
    170190                ({}, [(0.1, 0.1)]*2, [None]*2), 
    171191                # test that ER/VR will run if they exist 
    … …  
    174194                ] 
    175195 
    176             tests = self.info.tests 
     196            tests = smoke_tests + self.info.tests 
    177197            try: 
    178198                model = build_model(self.info, dtype=self.dtype, 
    179199                                    platform=self.platform) 
    180                 for test in smoke_tests + tests: 
    181                     self.run_one(model, test) 
    182  
    183                 if not tests and self.platform == "dll": 
    184                     ## Uncomment the following to make forgetting the test 
    185                     ## values an error.  Only do so for the "dll" tests 
    186                     ## to reduce noise from both opencl and dll, and because 
    187                     ## python kernels use platform="dll". 
    188                     #raise Exception("No test cases provided") 
    189                     pass 
     200                results = [self.run_one(model, test) for test in tests] 
     201                if self.stash: 
     202                    for test, target, actual in zip(tests, self.stash[0], results): 
     203                        assert np.all(abs(target-actual) < 5e-5*abs(actual)),\ 
     204                            "GPU/CPU comparison expected %s but got %s for %s"%(target, actual, test[0]) 
     205                else: 
     206                    self.stash.append(results) 
     207 
     208                # Check for missing tests.  Only do so for the "dll" tests 
     209                # to reduce noise from both opencl and dll, and because 
     210                # python kernels use platform="dll". 
     211                if self.platform == "dll": 
     212                    missing = [] 
     213                    ## Uncomment the following to require test cases 
     214                    #missing = self._find_missing_tests() 
     215                    if missing: 
     216                        raise ValueError("Missing tests for "+", ".join(missing)) 
    190217 
    191218            except: 
    192219                annotate_exception(self.test_name) 
    193220                raise 
     221 
     222        def _find_missing_tests(self): 
     223            # type: () -> None 
     224            """make sure there are 1D, 2D, ER and VR tests as appropriate""" 
     225            model_has_VR = callable(self.info.VR) 
     226            model_has_ER = callable(self.info.ER) 
     227            model_has_1D = True 
     228            model_has_2D = any(p.type == 'orientation' 
     229                               for p in self.info.parameters.kernel_parameters) 
     230 
     231            # Lists of tests that have a result that is not None 
     232            single = [test for test in self.info.tests 
     233                      if not isinstance(test[2], list) and test[2] is not None] 
     234            tests_has_VR = any(test[1] == 'VR' for test in single) 
     235            tests_has_ER = any(test[1] == 'ER' for test in single) 
     236            tests_has_1D_single = any(isinstance(test[1], float) for test in single) 
     237            tests_has_2D_single = any(isinstance(test[1], tuple) for test in single) 
     238 
     239            multiple = [test for test in self.info.tests 
     240                        if isinstance(test[2], list) 
     241                            and not all(result is None for result in test[2])] 
     242            tests_has_1D_multiple = any(isinstance(test[1][0], float) 
     243                                        for test in multiple) 
     244            tests_has_2D_multiple = any(isinstance(test[1][0], tuple) 
     245                                        for test in multiple) 
     246 
     247            missing = [] 
     248            if model_has_VR and not tests_has_VR: 
     249                missing.append("VR") 
     250            if model_has_ER and not tests_has_ER: 
     251                missing.append("ER") 
     252            if model_has_1D and not (tests_has_1D_single or tests_has_1D_multiple): 
     253                missing.append("1D") 
     254            if model_has_2D and not (tests_has_2D_single or tests_has_2D_multiple): 
     255                missing.append("2D") 
     256 
     257            return missing 
    194258 
    195259        def run_one(self, model, test): 
    … …  
    207271 
    208272            if x[0] == 'ER': 
    209                 actual = [call_ER(model.info, pars)] 
     273                actual = np.array([call_ER(model.info, pars)]) 
    210274            elif x[0] == 'VR': 
    211                 actual = [call_VR(model.info, pars)] 
     275                actual = np.array([call_VR(model.info, pars)]) 
    212276            elif isinstance(x[0], tuple): 
    213277                qx, qy = zip(*x) 
    … …  
    238302                                    'f(%s); expected:%s; actual:%s' 
    239303                                    % (xi, yi, actual_yi)) 
     304            return actual 
    240305 
    241306    return ModelTestCase 

  • sasmodels/modelinfo.py

    r85fe7f8 r9c44b7b  
    230230    defined as a sublist with the following elements: 
    231231 
    232     *name* is the name that will be used in the call to the kernel 
    233     function and the name that will be displayed to the user.  Names 
     232    *name* is the name that will be displayed to the user.  Names 
    234233    should be lower case, with words separated by underscore.  If 
    235     acronyms are used, the whole acronym should be upper case. 
     234    acronyms are used, the whole acronym should be upper case. For vector 
     235    parameters, the name will be followed by *[len]* where *len* is an 
     236    integer length of the vector, or the name of the parameter which 
     237    controls the length.  The attribute *id* will be created from name 
     238    without the length. 
    236239 
    237240    *units* should be one of *degrees* for angles, *Ang* for lengths, 
    … …  
    603606        # Using the call_parameters table, we already have expanded forms 
    604607        # for each of the vector parameters; put them in a lookup table 
    605         expanded_pars = dict((p.name, p) for p in self.call_parameters) 
     608        # Note: p.id and p.name are currently identical for the call parameters 
     609        expanded_pars = dict((p.id, p) for p in self.call_parameters) 
    606610 
    607611        def append_group(name): 
    … …  
    730734    info.docs = kernel_module.__doc__ 
    731735    info.category = getattr(kernel_module, 'category', None) 
    732     info.single = getattr(kernel_module, 'single', True) 
    733     info.opencl = getattr(kernel_module, 'opencl', True) 
    734736    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
    735737    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    … …  
    745747    info.profile = getattr(kernel_module, 'profile', None) # type: ignore 
    746748    info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore 
     749    # Default single and opencl to True for C models.  Python models have callable Iq. 
     750    info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq)) 
     751    info.single = getattr(kernel_module, 'single', not callable(info.Iq)) 
    747752 
    748753    # multiplicity info 

  • sasmodels/models/core_multi_shell.c

    r925ad6e rc3ccaec  
    88} 
    99 
    10 double 
    11 form_volume(double core_radius, double n, double thickness[]); 
    12 double 
    13 form_volume(double core_radius, double n, double thickness[]) 
     10static double 
     11form_volume(double core_radius, double fp_n, double thickness[]) 
    1412{ 
    1513  double r = core_radius; 
     14  int n = (int)(fp_n+0.5); 
    1615  for (int i=0; i < n; i++) { 
    1716    r += thickness[i]; 
    … …  
    2019} 
    2120 
    22 double 
     21static double 
    2322Iq(double q, double core_sld, double core_radius, 
    24    double solvent_sld, double num_shells, double sld[], double thickness[]); 
    25 double 
    26 Iq(double q, double core_sld, double core_radius, 
    27    double solvent_sld, double num_shells, double sld[], double thickness[]) 
     23   double solvent_sld, double fp_n, double sld[], double thickness[]) 
    2824{ 
    29   const int n = (int)ceil(num_shells); 
     25  const int n = (int)(fp_n+0.5); 
    3026  double f, r, last_sld; 
    3127  r = core_radius; 

  • sasmodels/models/core_multi_shell.py

    r925ad6e r5a0b3d7  
    107107    Returns the SLD profile *r* (Ang), and *rho* (1e-6/Ang^2). 
    108108    """ 
     109    n = int(n+0.5) 
    109110    z = [] 
    110111    rho = [] 
    … …  
    133134def ER(radius, n, thickness): 
    134135    """Effective radius""" 
    135     n = int(n[0])  # n cannot be polydisperse 
     136    n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse 
    136137    return np.sum(thickness[:n], axis=0) + radius 
    137138 

  • sasmodels/models/core_shell_parallelepiped.py

    r797a8e3 rcb0dc22  
    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" 
    7 can be different on all three (pairs) of faces.** 
     6The thickness and the scattering length density of the shell or  
     7"rim" can be different on each (pair) of faces. However at this time 
     8the model does **NOT** actually calculate a c face rim despite the presence of 
     9the parameter. 
     10 
     11.. note:: 
     12   This model was originally ported from NIST IGOR macros. However,t is not 
     13   yet fully understood by the SasView developers and is currently review. 
    814 
    915The form factor is normalized by the particle volume $V$ such that 
    … …  
    3541    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    3642 
    37 **meaning that there are "gaps" at the corners of the solid.** 
     43**meaning that there are "gaps" at the corners of the solid.**  Again note that 
     44$t_C = 0$ currently.  
    3845 
    3946The intensity calculated follows the :ref:`parallelepiped` model, with the 
    … …  
    4148amplitudes of the core and shell, in the same manner as a core-shell model. 
    4249 
    43 .. math:: 
    44  
    45     F_{a}(Q,\alpha,\beta)= 
    46     \Bigg(\frac{sin(Q(L_A+2t_A)/2sin\alpha sin\beta)}{Q(L_A+2t_A)/2sin\alpha 
    47     sin\beta)} 
    48     - \frac{sin(QL_A/2sin\alpha sin\beta)}{QL_A/2sin\alpha sin\beta)} \Bigg) 
    49     + \frac{sin(QL_B/2sin\alpha sin\beta)}{QL_B/2sin\alpha sin\beta)} 
    50     + \frac{sin(QL_C/2sin\alpha sin\beta)}{QL_C/2sin\alpha sin\beta)} 
    5150 
    5251.. note:: 
    … …  
    9392    detector plane. 
    9493 
    95 Validation 
    96 ---------- 
    97  
    98 The model uses the form factor calculations implemented in a c-library provided 
    99 by the NIST Center for Neutron Research (Kline, 2006). 
    100  
    10194References 
    10295---------- 
    … …  
    113106 
    114107* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    115 * **Last Modified by:** Paul Butler **Date:** September 30, 2016 
    116 * **Last Reviewed by:** Miguel Gonzales **Date:** March 21, 2016 
     108* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
     109* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017 
     110* **Currently Under review by:** Paul Butler 
    117111""" 
    118112 

  • sasmodels/models/ellipsoid.c

    r925ad6e r130d4c7  
    44    double radius_polar, double radius_equatorial, double theta, double phi); 
    55 
    6 double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha); 
    7 double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha) 
     6static double 
     7_ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha) 
    88{ 
    99    double ratio = radius_polar/radius_equatorial; 
    10     // Given the following under the radical: 
    11     //     1 + sin^2(T) (v^2 - 1) 
    12     // we can expand to match the form given in Guinier (1955) 
    13     //     = (1 - sin^2(T)) + v^2 sin^2(T) = cos^2(T) + sin^2(T) 
     10    // Using ratio v = Rp/Re, we can expand the following to match the 
     11    // form given in Guinier (1955) 
     12    //     r = Re * sqrt(1 + cos^2(T) (v^2 - 1)) 
     13    //       = Re * sqrt( (1 - cos^2(T)) + v^2 cos^2(T) ) 
     14    //       = Re * sqrt( sin^2(T) + v^2 cos^2(T) ) 
     15    //       = sqrt( Re^2 sin^2(T) + Rp^2 cos^2(T) ) 
     16    // 
    1417    // Instead of using pythagoras we could pass in sin and cos; this may be 
    1518    // slightly better for 2D which has already computed it, but it introduces 
    … …  
    1720    // leave it as is. 
    1821    const double r = radius_equatorial 
    19                      * sqrt(1.0 + sin_alpha*sin_alpha*(ratio*ratio - 1.0)); 
     22                     * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0)); 
    2023    const double f = sas_3j1x_x(q*r); 
    2124 
    … …  
    3942    double total = 0.0; 
    4043    for (int i=0;i<76;i++) { 
    41         //const double sin_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
    42         const double sin_alpha = Gauss76Z[i]*zm + zb; 
    43         total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 
     44        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
     45        const double cos_alpha = Gauss76Z[i]*zm + zb; 
     46        total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
    4447    } 
    4548    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    5962    double q, sin_alpha, cos_alpha; 
    6063    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    61     const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 
     64    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
    6265    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    6366 

  • sasmodels/models/flexible_cylinder.c

    r592343f rc3ccaec  
    11static double 
    2 form_volume(length, kuhn_length, radius) 
     2form_volume(double length, double kuhn_length, double radius) 
    33{ 
    44    return 1.0; 

  • sasmodels/models/fractal.c

    r925ad6e r4788822  
    11#define INVALID(p) (p.fractal_dim < 0.0) 
     2 
     3 static double 
     4 form_volume(double radius) 
     5 { 
     6     return M_4PI_3 * cube(radius); 
     7 } 
    28 
    39static double 
    … …  
    1319 
    1420    //calculate P(q) for the spherical subunits 
    15     const double V = M_4PI_3*cube(radius); 
    16     const double pq = V * square((sld_block-sld_solvent)*sas_3j1x_x(q*radius)); 
     21    const double pq = square(form_volume(radius) * (sld_block-sld_solvent) 
     22                      *sas_3j1x_x(q*radius)); 
    1723 
    1824    // scale to units cm-1 sr-1 (assuming data on absolute scale) 

  • sasmodels/models/fractal.py

    r925ad6e rdf89d77  
    2020.. math:: 
    2121 
    22     P(q)&= F(qR_0)^2 
    23  
    24     F(q)&= \frac{3 (\sin x - x \cos x)}{x^3} 
    25  
    26     V_\text{particle} &= \frac{4}{3}\ \pi R_0 
    27  
     22    P(q)&= F(qR_0)^2 \\ 
     23    F(q)&= \frac{3 (\sin x - x \cos x)}{x^3} \\ 
     24    V_\text{particle} &= \frac{4}{3}\ \pi R_0 \\ 
    2825    S(q) &= 1 + \frac{D_f\  \Gamma\!(D_f-1)}{[1+1/(q \xi)^2\  ]^{(D_f -1)/2}} 
    2926    \frac{\sin[(D_f-1) \tan^{-1}(q \xi) ]}{(q R_0)^{D_f}} 
    … …  
    3229is the fractal dimension, representing the self similarity of the structure.  
    3330Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function  
    34 diverges at $D_f$=0 and $D_f$=1.   
     31diverges at $D_f=0$ and $D_f=1$. 
    3532 
    3633**Polydispersity on the radius is provided for.** 
    … …  
    4744---------- 
    4845 
    49 J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 
     46.. [#] J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 
    5047 
    51 **Author:** NIST IGOR/DANSE **on:** pre 2010 
     48Authorship and Verification 
     49---------------------------- 
    5250 
    53 **Last Modified by:** Paul Butler **on:** March 20, 2016 
    54  
    55 **Last Reviewed by:** Paul Butler **on:** March 20, 2016 
     51* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     52* **Converted to sasmodels by:** Paul Butler **Date:** March 19, 2016 
     53* **Last Modified by:** Paul Butler **Date:** March 12, 2017 
     54* **Last Reviewed by:** Paul Butler **Date:** March 12, 2017 
    5655 
    5756""" 
    … …  
    8483parameters = [["volfraction", "", 0.05, [0.0, 1], "", 
    8584               "volume fraction of blocks"], 
    86               ["radius",    "Ang",  5.0, [0.0, inf], "", 
     85              ["radius",    "Ang",  5.0, [0.0, inf], "volume", 
    8786               "radius of particles"], 
    8887              ["fractal_dim",      "",  2.0, [0.0, 6.0], "", 

  • sasmodels/models/hayter_msa.c

    r4962519 r3fc5d27  
    7070                SofQ=sqhcal(Qdiam, gMSAWave); 
    7171        }else{ 
    72         //SofQ=NaN; 
    73                 SofQ=-1.0; 
     72        SofQ=NAN; 
    7473                //      print "Error Level = ",ierr 
    7574                //      print "Please report HPMSA problem with above error code" 

  • sasmodels/models/lamellar_hg_stack_caille.c

    ra807206 r1c6286d  
    33*/ 
    44 
    5 double Iq(double qval, 
    6       double length_tail, 
    7       double length_head, 
    8       double Nlayers,  
    9       double dd, 
    10       double Cp, 
    11       double tail_sld, 
    12       double head_sld, 
    13       double solvent_sld); 
    14  
    15 double Iq(double qval, 
    16       double length_tail, 
    17       double length_head, 
    18       double Nlayers,  
    19       double dd, 
    20       double Cp, 
    21       double tail_sld, 
    22       double head_sld, 
    23       double solvent_sld) 
     5static double 
     6Iq(double qval, 
     7   double length_tail, 
     8   double length_head, 
     9   double fp_Nlayers, 
     10   double dd, 
     11   double Cp, 
     12   double tail_sld, 
     13   double head_sld, 
     14   double solvent_sld) 
    2415{ 
    25   double NN;   //local variables of coefficient wave 
     16  int Nlayers = (int)(fp_Nlayers+0.5);    //cast to an integer for the loop 
    2617  double inten,Pq,Sq,alpha,temp,t2; 
    2718  //double dQ, dQDefault, t1, t3; 
    28   int ii,NNint; 
    2919  // from wikipedia 0.577215664901532860606512090082402431042159335 
    3020  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015 
    … …  
    3222  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015 
    3323 
    34   NN = trunc(Nlayers);    //be sure that NN is an integer 
    35    
    36   Pq = (head_sld-solvent_sld)*(sin(qval*(length_head+length_tail))-sin(qval*length_tail)) + 
    37               (tail_sld-solvent_sld)*sin(qval*length_tail); 
     24  Pq = (head_sld-solvent_sld)*(sin(qval*(length_head+length_tail))-sin(qval*length_tail)) 
     25       + (tail_sld-solvent_sld)*sin(qval*length_tail); 
    3826  Pq *= Pq; 
    3927  Pq *= 4.0/(qval*qval); 
    4028 
    41   NNint = (int)NN;    //cast to an integer for the loop 
    42   ii=0; 
    4329  Sq = 0.0; 
    44   for(ii=1;ii<=(NNint-1);ii+=1) { 
    45  
    46     //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015 
    47  
     30  for(int ii=1; ii < Nlayers; ii++) { 
    4831    temp = 0.0; 
    4932    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler); 
    … …  
    5235    //t3 = dQ*dQ*dd*dd*ii*ii; 
    5336 
    54     temp = 1.0-ii/NN; 
     37    temp = 1.0-(double)ii/(double)Nlayers; 
    5538    //temp *= cos(dd*qval*ii/(1.0+t1)); 
    5639    temp *= cos(dd*qval*ii); 
    5740    //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii); 
    5841    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 
    59     temp *= exp(-t2/2.0 ); 
     42    temp *= exp(-t2/2.0); 
    6043    //temp /= sqrt(1.0+t1); 
    6144 
    … …  
    7154 
    7255  inten *= 1.0e-04;   // 1/A to 1/cm 
    73   return(inten); 
     56  return inten; 
    7457} 
    7558 

  • sasmodels/models/lamellar_hg_stack_caille.py

    r7c57861 ra57b31d  
    9898    ["length_head", "Ang", 2, [0, inf], "volume", 
    9999     "head thickness"], 
    100     ["Nlayers", "", 30, [0, inf], "", 
     100    ["Nlayers", "", 30, [1, inf], "", 
    101101     "Number of layers"], 
    102102    ["d_spacing", "Ang", 40., [0.0, inf], "volume", 

  • sasmodels/models/lamellar_stack_caille.c

    r0bef47b r1c6286d  
    33*/ 
    44 
    5 double Iq(double qval, 
    6       double del, 
    7       double Nlayers,  
    8       double dd, 
    9           double Cp,  
    10       double sld, 
    11       double solvent_sld); 
    12  
    13 double Iq(double qval, 
    14       double del, 
    15       double Nlayers,  
    16       double dd, 
    17           double Cp,  
    18       double sld, 
    19       double solvent_sld) 
     5static double 
     6Iq(double qval, 
     7   double del, 
     8   double fp_Nlayers, 
     9   double dd, 
     10   double Cp, 
     11   double sld, 
     12   double solvent_sld) 
    2013{ 
    21   double contr,NN;   //local variables of coefficient wave 
     14  int Nlayers = (int)(fp_Nlayers+0.5);    //cast to an integer for the loop 
     15  double contr;   //local variables of coefficient wave 
    2216  double inten,Pq,Sq,alpha,temp,t2; 
    2317  //double dQ, dQDefault, t1, t3; 
    24   int ii,NNint; 
    2518  // from wikipedia 0.577215664901532860606512090082402431042159335 
    2619  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015 
    … …  
    2821  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015 
    2922 
    30   NN = trunc(Nlayers);    //be sure that NN is an integer 
    31    
    3223  contr = sld - solvent_sld; 
    3324 
    3425  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)); 
    3526 
    36   NNint = (int)NN;    //cast to an integer for the loop 
    37   ii=0; 
    3827  Sq = 0.0; 
    39   // the vital "=" in ii<=  added March 2015 
    40   for(ii=1;ii<=(NNint-1);ii+=1) { 
    41  
    42     //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015 
    43  
     28  for (int ii=1; ii < Nlayers; ii++) { 
    4429    temp = 0.0; 
    4530    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler); 
    … …  
    4833    //t3 = dQ*dQ*dd*dd*ii*ii; 
    4934 
    50     temp = 1.0-ii/NN; 
     35    temp = 1.0 - (double)ii / (double)Nlayers; 
    5136    //temp *= cos(dd*qval*ii/(1.0+t1)); 
    5237    temp *= cos(dd*qval*ii); 
    5338    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 
    54     temp *= exp(-t2/2.0 ); 
     39    temp *= exp(-t2/2.0); 
    5540    //temp /= sqrt(1.0+t1); 
    5641 

  • sasmodels/models/lamellar_stack_caille.py

    r7c57861 ra57b31d  
    8888parameters = [ 
    8989    ["thickness",        "Ang",      30.0,  [0, inf],   "volume", "sheet thickness"], 
    90     ["Nlayers",          "",          20,   [0, inf],   "",       "Number of layers"], 
     90    ["Nlayers",          "",          20,   [1, inf],   "",       "Number of layers"], 
    9191    ["d_spacing",        "Ang",      400.,  [0.0,inf],  "volume", "lamellar d-spacing of Caille S(Q)"], 
    9292    ["Caille_parameter", "1/Ang^2",    0.1, [0.0,0.8],  "",       "Caille parameter"], 

  • sasmodels/models/lamellar_stack_paracrystal.c

    r4962519 r5467cd8  
    22 
    33*/ 
    4 double Iq(double qval, 
    5       double th, 
    6       double Nlayers,  
    7           double davg,  
    8           double pd, 
    9       double sld, 
    10       double solvent_sld); 
    11 double paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an); 
    12 double paraCryst_an(double ww, double qval, double davg, long Nlayers); 
     4double paraCryst_sn(double ww, double qval, double davg, int Nlayers, double an); 
     5double paraCryst_an(double ww, double qval, double davg, int Nlayers); 
    136 
    14 double Iq(double qval, 
    15       double th, 
    16       double Nlayers,  
    17           double davg,  
    18           double pd, 
    19       double sld, 
    20       double solvent_sld) 
     7static double 
     8Iq(double qval, 
     9   double th, 
     10   double fp_Nlayers, 
     11   double davg, 
     12   double pd, 
     13   double sld, 
     14   double solvent_sld) 
    2115{ 
    22      
    23         double inten,contr,xn; 
    24         double xi,ww,Pbil,Znq,Snq,an; 
    25         long n1,n2; 
     16        //get the fractional part of Nlayers, to determine the "mixing" of N's 
     17        int n1 = (int)(fp_Nlayers);             //truncate towards zero 
     18        int n2 = n1 + 1; 
     19        const double xn = (double)n2 - fp_Nlayers;      //fractional contribution of n1 
    2620         
    27         contr = sld - solvent_sld; 
    28         //get the fractional part of Nlayers, to determine the "mixing" of N's 
    29          
    30         n1 = (long)trunc(Nlayers);              //rounds towards zero 
    31         n2 = n1 + 1; 
    32         xn = (double)n2 - Nlayers;                      //fractional contribution of n1 
    33          
    34         ww = exp(-qval*qval*pd*pd*davg*davg/2.0); 
     21        const double ww = exp(-0.5*square(qval*pd*davg)); 
    3522 
    3623        //calculate the n1 contribution 
     24        double Znq,Snq,an; 
    3725        an = paraCryst_an(ww,qval,davg,n1); 
    3826        Snq = paraCryst_sn(ww,qval,davg,n1,an); 
    39          
     27 
    4028        Znq = xn*Snq; 
    4129         
    … …  
    5240//      Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg)) 
    5341         
    54         xi = th/2.0;            //use 1/2 the bilayer thickness 
    55         Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi)); 
     42        const double xi = th/2.0;               //use 1/2 the bilayer thickness 
     43        const double Pbil = square(sas_sinx_x(qval*xi)); 
    5644         
    57         inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval); 
    58         inten *= 1.0e-04; 
     45        const double contr = sld - solvent_sld; 
     46        const double inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval); 
    5947//printf("q=%.7e wwm1=%g ww=%.5e an=% 12.5e Snq=% 12.5e Znq=% 12.5e Pbil=% 12.5e\n",qval,wwm1,ww,an,Snq,Znq,Pbil); 
    60         return(inten); 
     48        return 1.0e-4*inten; 
    6149} 
    6250 
    6351// functions for the lamellar paracrystal model 
    6452double 
    65 paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an) { 
     53paraCryst_sn(double ww, double qval, double davg, int Nlayers, double an) { 
    6654         
    6755        double Snq; 
    … …  
    6957        Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) ); 
    7058         
    71         return(Snq); 
     59        return Snq; 
    7260} 
    7361 
    7462double 
    75 paraCryst_an(double ww, double qval, double davg, long Nlayers) { 
    76          
     63paraCryst_an(double ww, double qval, double davg, int Nlayers) { 
    7764        double an; 
    7865         
    … …  
    8269        an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg); 
    8370         
    84         return(an); 
     71        return an; 
    8572} 
    8673 

  • sasmodels/models/lamellar_stack_paracrystal.py

    r7c57861 ra0168e8  
    113113parameters = [["thickness", "Ang", 33.0, [0, inf], "volume", 
    114114               "sheet thickness"], 
    115               ["Nlayers", "", 20, [0, inf], "", 
     115              ["Nlayers", "", 20, [1, inf], "", 
    116116               "Number of layers"], 
    117117              ["d_spacing", "Ang", 250., [0.0, inf], "", 

  • sasmodels/models/linear_pearls.c

    r925ad6e rc3ccaec  
    44            double radius, 
    55            double edge_sep, 
    6             double num_pearls, 
     6            double fp_num_pearls, 
    77            double pearl_sld, 
    88            double solvent_sld); 
    … …  
    1111            double radius, 
    1212            double edge_sep, 
    13             double num_pearls, 
     13            int num_pearls, 
    1414            double pearl_sld, 
    1515            double solvent_sld); 
    1616 
    1717 
    18 double form_volume(double radius, double num_pearls) 
     18double form_volume(double radius, double fp_num_pearls) 
    1919{ 
     20    int num_pearls = (int)(fp_num_pearls + 0.5); 
    2021    // Pearl volume 
    2122    double pearl_vol = M_4PI_3 * cube(radius); 
    … …  
    2728            double radius, 
    2829            double edge_sep, 
    29             double num_pearls, 
     30            int num_pearls, 
    3031            double pearl_sld, 
    3132            double solvent_sld) 
    3233{ 
    33     double n_contrib; 
    3434    //relative sld 
    3535    double contrast_pearl = pearl_sld - solvent_sld; 
    … …  
    4646    double psi = sas_3j1x_x(q * radius); 
    4747 
    48     // N pearls contribution 
    49     int n_max = num_pearls - 1; 
    50     n_contrib = num_pearls; 
    51     for(int num=1; num<=n_max; num++) { 
    52         n_contrib += (2.0*(num_pearls-num)*sas_sinx_x(q*separation*num)); 
     48    // N pearls interaction terms  
     49    double structure_factor = (double)num_pearls; 
     50    for(int num=1; num<num_pearls; num++) { 
     51        structure_factor += 2.0*(num_pearls-num)*sas_sinx_x(q*separation*num); 
    5352    } 
    5453    // form factor for num_pearls 
    55     double form_factor = 1.0e-4 * n_contrib * square(m_s*psi) / tot_vol; 
     54    double form_factor = 1.0e-4 * structure_factor * square(m_s*psi) / tot_vol; 
    5655 
    5756    return form_factor; 
    … …  
    6160            double radius, 
    6261            double edge_sep, 
    63             double num_pearls, 
     62            double fp_num_pearls, 
    6463            double pearl_sld, 
    6564            double solvent_sld) 
    6665{ 
    6766 
     67    int num_pearls = (int)(fp_num_pearls + 0.5); 
    6868        double result = linear_pearls_kernel(q, 
    6969                    radius, 

  • sasmodels/models/linear_pearls.py

    r925ad6e rc3ccaec  
    1616.. math:: 
    1717 
    18     P(Q) = \frac{scale}{V}\left[ m_{p}^2 
    19     \left(N+2\sum_{n-1}^{N-1}(N-n)\frac{sin(qnl)}{qnl}\right) 
    20     \left( 3\frac{sin(qR)-qRcos(qR)}{(qr)^3}\right)^2\right] 
     18    P(Q) = \frac{\text{scale}}{V}\left[ m_{p}^2 
     19    \left(N+2\sum_{n-1}^{N-1}(N-n)\frac{\sin(qnl)}{qnl}\right) 
     20    \left( 3\frac{\sin(qR)-qR\cos(qR)}{(qr)^3}\right)^2\right] 
    2121 
    2222where the mass $m_p$ is $(SLD_{pearl}-SLD_{solvent})*(volume\ of\ N\ pearls)$. 
    … …  
    5656    ["radius",      "Ang",       80.0, [0, inf],     "", "Radius of the pearls"], 
    5757    ["edge_sep",    "Ang",      350.0, [0, inf],     "", "Length of the string segment - surface to surface"], 
    58     ["num_pearls",  "",           3.0, [0, inf],     "", "Number of the pearls"], 
     58    ["num_pearls",  "",           3.0, [1, inf],     "", "Number of the pearls"], 
    5959    ["sld",   "1e-6/Ang^2", 1.0, [-inf, inf],  "sld", "SLD of the pearl spheres"], 
    6060    ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf],  "sld", "SLD of the solvent"], 

  • sasmodels/models/multilayer_vesicle.c

    r925ad6e rec1d4bc  
    1 static 
    2 double multilayer_vesicle_kernel(double q, 
     1static double 
     2form_volume(double radius, 
     3          double thick_shell, 
     4          double thick_solvent, 
     5          double fp_n_shells) 
     6{ 
     7    int n_shells = (int)(fp_n_shells + 0.5); 
     8    double R_N = radius + n_shells*(thick_shell+thick_solvent) - thick_solvent; 
     9    return M_4PI_3*cube(R_N); 
     10} 
     11 
     12static double 
     13multilayer_vesicle_kernel(double q, 
    314          double volfraction, 
    415          double radius, 
    … …  
    718          double sld_solvent, 
    819          double sld, 
    9           double n_pairs) 
     20          int n_shells) 
    1021{ 
    1122    //calculate with a loop, two shells at a time 
    … …  
    2940 
    3041        //do 2 layers at a time 
    31         ii += 1; 
     42        ii++; 
    3243 
    33     } while(ii <= n_pairs-1);  //change to make 0 < n_pairs < 2 correspond to 
     44    } while(ii <= n_shells-1);  //change to make 0 < n_shells < 2 correspond to 
    3445                               //unilamellar vesicles (C. Glinka, 11/24/03) 
    3546 
    36     fval *= volfraction*1.0e-4*fval/voli; 
    37  
    38     return(fval); 
     47    return 1.0e-4*volfraction*fval*fval;  // Volume normalization happens in caller 
    3948} 
    4049 
    41 static 
    42 double Iq(double q, 
     50static double 
     51Iq(double q, 
    4352          double volfraction, 
    4453          double radius, 
    … …  
    4756          double sld_solvent, 
    4857          double sld, 
    49           double n_pairs) 
     58          double fp_n_shells) 
    5059{ 
     60    int n_shells = (int)(fp_n_shells + 0.5); 
    5161    return multilayer_vesicle_kernel(q, 
    5262           volfraction, 
    … …  
    5666           sld_solvent, 
    5767           sld, 
    58            n_pairs); 
     68           n_shells); 
    5969} 
    6070 

  • sasmodels/models/multilayer_vesicle.py

    r925ad6e r5d23de2  
    33---------- 
    44 
    5 This model is a trivial extension of the core_shell_sphere function to include 
    6 *N* shells where the core is filled with solvent and the shells are interleaved 
    7 with layers of solvent. For *N = 1*, this returns the same as the vesicle model, 
    8 except for the normalisation, which here is to outermost volume. 
    9 The shell thicknessess and SLD are constant for all shells as expected for 
    10 a multilayer vesicle. 
     5This model is a trivial extension of the core_shell_sphere function where the 
     6core is filled with solvent and is surrounded by $N$ shells of material 
     7(such as lipids) interleaved with $N - 1$ layers of solvent. For $N = 1$, this 
     8returns the same as the vesicle model, except for the normalisation, which here 
     9is to outermost volume. The shell thicknesses and SLD are constant for all 
     10shells as expected for a multilayer vesicle. 
    1111 
    1212.. figure:: img/multi_shell_geometry.jpg 
    … …  
    1919 
    2020.. math:: 
    21  
    22     P(q) = \frac{\text{scale.volfraction}}{V_t} F^2(q) + \text{background} 
     21    P(q) = \text{scale} \cdot \frac{\phi}{V(R_N)} F^2(q) + \text{background} 
    2322 
    2423where 
    2524 
    2625.. math:: 
    27  
    28      F(q) = (\rho_{shell}-\rho_{solv}) \sum_{i=1}^{n\_pairs} \left[ 
    29      3V(R_i)\frac{\sin(qR_i)-qR_i\cos(qR_i)}{(qR_i)^3} \\ 
    30       - 3V(R_i+t_s)\frac{\sin(q(R_i+t_s))-q(R_i+t_s)\cos(q(R_i+t_s))}{(q(R_i+t_s))^3} 
     26     F(q) = (\rho_\text{shell}-\rho_\text{solv}) \sum_{i=1}^{N} \left[ 
     27     3V(r_i)\frac{\sin(qr_i) - qr_i\cos(qr_i)}{(qr_i)^3} 
     28     - 3V(R_i)\frac{\sin(qR_i) - qR_i\cos(qR_i)}{(qR_i)^3} 
    3129     \right] 
    3230 
     31for 
    3332 
    34 where $R_i = r_c + (i-1)(t_s + t_w)$ 
    35     
    36 where $V_t$ is the volume of the whole particle, $V(R)$ is the volume of a sphere 
    37 of radius $R$, $r_c$ is the radius of the core, $\rho_{shell}$ is the scattering length  
    38 density of a shell, $\rho_{solv}$ is the scattering length density of the solvent. 
     33.. math:: 
    3934 
     35     r_i &= r_c + (i-1)(t_s + t_w) && \text{ solvent radius before shell } i \\ 
     36     R_i &= r_i + t_s && \text{ shell radius for shell } i 
     37 
     38$\phi$ is the volume fraction of particles, $V(r)$ is the volume of a sphere 
     39of radius $r$, $r_c$ is the radius of the core, $t_s$ is the thickness of 
     40the shell, $t_w$ is the thickness of the solvent layer between the shells, 
     41$\rho_\text{shell}$ is the scattering length density of a shell, and 
     42$\rho_\text{solv}$ is the scattering length density of the solvent. 
     43 
     44USAGE NOTES 
     45 
     46* The outer-most shell radius $R_N$ is used as the effective radius 
     47  for $P(Q)$ when $P(Q) * S(Q)$ is applied. 
     48  calculations rather slow. 
     49* The number of shells is always rounded to an integer value as a non interger 
     50  number of layers is not physical. 
     51* Thus Polydispersity should only be applied to number of shells **VERY 
     52  CAREFULLY**.  A possible legitimate use would be for mixed systems in which 
     53  some vesicles have 1 shell, some have 2, etc. A polydispersity on $N$ can be 
     54  used to model the data by using the "array distriubtion" feature. First 
     55  create a file such as *shell_dist.txt* containing the relative portion 
     56  of each vesicle size:: 
     57 
     58    1 20 
     59    2  4 
     60    3  1 
     61 
     62  Turn on polydispersity and select an array distribution for the *n_shells* 
     63  parameter.  Choose the above *shell_dist.txt* file, and the model will be 
     64  computed with 80% 1-shell vesicles, 16% 2-shell vesicles and 4% 
     65  3-shell vesicles. 
     66* This is a highly non-linear, highly oscillatory (especially around the 
     67  q-values that correspond to the repeat distance of the layers), model 
     68  function complicated by the fact that the number of water/shell pairs must 
     69  physically be an integer value, although the optimization treats it as a 
     70  floating point value. Thus it may be that the resolution interpolation is not 
     71  sufficiently fine grained in certain cases. Please report any such occurences 
     72  to the SasView team. Generally, for the best possible experience: 
     73 * Start with the best possible guess 
     74 * Using a priori knowledge, hold as many parameters fixed as possible 
     75 * if N=1, tw (water thickness) must by definition be zero. Both N and tw should 
     76   be fixed during fitting. 
     77 * If N>1, use constraints to keep N > 1 
     78 * Because N only really moves in integer steps, it may get "stuck" if the 
     79   optimizer step size is too small so care should be taken 
     80   If you experience problems with this please contact the SasView team and let 
     81   them know the issue preferably with example data and model which fail to 
     82   converge. 
    4083 
    4184The 2D scattering intensity is the same as 1D, regardless of the orientation 
    … …  
    4689    q = \sqrt{q_x^2 + q_y^2} 
    4790 
    48  
    49 The outer most radius 
    50  
    51 $radius + n\_pairs * thick\_shell + (n\_pairs- 1) * thick\_solvent$ 
    52  
    53 is used for both the volume fraction normalization and for the  
    54 effective radius for *S(Q)* when $P(Q) * S(Q)$ is applied. 
    55  
    5691For information about polarised and magnetic scattering, see 
    5792the :ref:`magnetism` documentation. 
    58  
    59 This code is based on the form factor calculations implemented in the NIST 
    60 Center for Neutron Research provided c-library (Kline, 2006). 
    6193 
    6294References 
    6395---------- 
    6496 
    65 B Cabane, *Small Angle Scattering Methods*, 
    66 in *Surfactant Solutions: New Methods of Investigation*, 
    67 Ch.2, Surfactant Science Series Vol. 22, Ed. R Zana and M Dekker, 
    68 New York, (1987). 
     97.. [#] B Cabane, *Small Angle Scattering Methods*, in *Surfactant Solutions: 
     98   New Methods of Investigation*, Ch.2, Surfactant Science Series Vol. 22, Ed. 
     99   R Zana and M Dekker, New York, (1987). 
    69100 
    70 **Author:** NIST IGOR/DANSE **on:** pre 2010 
     101Authorship and Verification 
     102---------------------------- 
    71103 
    72 **Last Modified by:** Piotr Rozyczko**on:** Feb 24, 2016 
    73  
    74 **Last Reviewed by:** Paul Butler **on:** March 20, 2016 
     104* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     105* **Converted to sasmodels by:** Piotr Rozyczko **Date:** Feb 24, 2016 
     106* **Last Modified by:** Paul Kienzle **Date:** Feb 7, 2017 
     107* **Last Reviewed by:** Paul Butler **Date:** March 12, 2017 
    75108 
    76109""" 
    … …  
    89122    sld_solvent: solvent scattering length density 
    90123    sld: shell scattering length density 
    91     n_pairs:number of "shell plus solvent" layer pairs 
     124    n_shells:number of "shell plus solvent" layer pairs 
    92125    background: incoherent background 
    93126        """ 
    … …  
    98131parameters = [ 
    99132    ["volfraction", "",  0.05, [0.0, 1],  "", "volume fraction of vesicles"], 
    100     ["radius", "Ang", 60.0, [0.0, inf],  "", "radius of solvent filled core"], 
    101     ["thick_shell", "Ang",        10.0, [0.0, inf],  "", "thickness of one shell"], 
    102     ["thick_solvent", "Ang",        10.0, [0.0, inf],  "", "solvent thickness between shells"], 
     133    ["radius", "Ang", 60.0, [0.0, inf],  "volume", "radius of solvent filled core"], 
     134    ["thick_shell", "Ang",        10.0, [0.0, inf],  "volume", "thickness of one shell"], 
     135    ["thick_solvent", "Ang",        10.0, [0.0, inf],  "volume", "solvent thickness between shells"], 
    103136    ["sld_solvent",    "1e-6/Ang^2",  6.4, [-inf, inf], "sld", "solvent scattering length density"], 
    104137    ["sld",   "1e-6/Ang^2",  0.4, [-inf, inf], "sld", "Shell scattering length density"], 
    105     ["n_pairs",     "",            2.0, [1.0, inf],  "", "Number of shell plus solvent layer pairs"], 
     138    ["n_shells",     "",            2.0, [1.0, inf],  "volume", "Number of shell plus solvent layer pairs"], 
    106139    ] 
    107140# pylint: enable=bad-whitespace, line-too-long 
    108141 
     142# TODO: proposed syntax for specifying which parameters can be polydisperse 
     143#polydispersity = ["radius", "thick_shell"] 
     144 
    109145source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 
    110146 
    111 polydispersity = ["radius", "n_pairs"] 
     147def ER(radius, thick_shell, thick_solvent, n_shells): 
     148    n_shells = int(n_shells+0.5) 
     149    return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent 
    112150 
    113151demo = dict(scale=1, background=0, 
    … …  
    118156            sld_solvent=6.4, 
    119157            sld=0.4, 
    120             n_pairs=2.0) 
     158            n_shells=2.0) 
    121159 
    122160tests = [ 
    … …  
    127165      'sld_solvent': 6.4, 
    128166      'sld': 0.4, 
    129       'n_pairs': 2.0, 
     167      'n_shells': 2.0, 
    130168      'scale': 1.0, 
    131169      'background': 0.001, 
    … …  
    138176      'sld_solvent': 6.4, 
    139177      'sld': 0.4, 
    140       'n_pairs': 2.0, 
     178      'n_shells': 2.0, 
    141179      'scale': 1.0, 
    142180      'background': 0.001, 

  • sasmodels/models/onion.c

    r925ad6e rc3ccaec  
    3030 
    3131static double 
    32 form_volume(double radius_core, double n, double thickness[]) 
     32form_volume(double radius_core, double n_shells, double thickness[]) 
    3333{ 
    34   int i; 
     34  int n = (int)(n_shells+0.5); 
    3535  double r = radius_core; 
    36   for (i=0; i < n; i++) { 
     36  for (int i=0; i < n; i++) { 
    3737    r += thickness[i]; 
    3838  } 

  • sasmodels/models/onion.py

    r925ad6e rc3ccaec  
    323323    Returns shape profile with x=radius, y=SLD. 
    324324    """ 
    325  
     325    n_shells = int(n_shells+0.5) 
    326326    total_radius = 1.25*(sum(thickness[:n_shells]) + radius_core + 1) 
    327327    dz = total_radius/400  # 400 points for a smooth plot 
    … …  
    366366    return np.asarray(z), np.asarray(rho) 
    367367 
    368 def ER(radius_core, n, thickness): 
     368def ER(radius_core, n_shells, thickness): 
    369369    """Effective radius""" 
    370     return np.sum(thickness[:int(n[0])], axis=0) + radius_core 
     370    n = int(n_shells[0]+0.5) 
     371    return np.sum(thickness[:n], axis=0) + radius_core 
    371372 
    372373demo = { 

  • sasmodels/models/pearl_necklace.c

    r4b541ac r3f853beb  
    11double form_volume(double radius, double edge_sep, 
    2     double thick_string, double num_pearls); 
     2    double thick_string, double fp_num_pearls); 
    33double Iq(double q, double radius, double edge_sep, 
    4     double thick_string, double num_pearls, double sld,  
     4    double thick_string, double fp_num_pearls, double sld, 
    55    double string_sld, double solvent_sld); 
    66 
    … …  
    99// From Igor library 
    1010static double 
    11 _pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string, 
    12     double num_pearls, double sld_pearl, double sld_string, double sld_solv) 
     11pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string, 
     12    int num_pearls, double sld_pearl, double sld_string, double sld_solv) 
    1313{ 
    1414    // number of string segments 
    15     num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
    16     const double num_strings = num_pearls - 1.0; 
     15    const int num_strings = num_pearls - 1; 
    1716 
    1817    //each masses: contrast * volume 
    … …  
    6968 
    7069double form_volume(double radius, double edge_sep, 
    71     double thick_string, double num_pearls) 
     70    double thick_string, double fp_num_pearls) 
    7271{ 
    73     num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
    74  
    75     const double num_strings = num_pearls - 1.0; 
     72    const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls 
     73    const int num_strings = num_pearls - 1; 
    7674    const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; 
    7775    const double pearl_vol = M_4PI_3 * radius * radius * radius; 
    … …  
    8280 
    8381double Iq(double q, double radius, double edge_sep, 
    84     double thick_string, double num_pearls, double sld,  
     82    double thick_string, double fp_num_pearls, double sld, 
    8583    double string_sld, double solvent_sld) 
    8684{ 
    87     const double form = _pearl_necklace_kernel(q, radius, edge_sep, 
     85    const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls 
     86    const double form = pearl_necklace_kernel(q, radius, edge_sep, 
    8887        thick_string, num_pearls, sld, string_sld, solvent_sld); 
    8988 

  • sasmodels/models/pearl_necklace.py

    r4b541ac r1bd1ea2  
    8282              ["thick_string", "Ang", 2.5, [0, inf], "volume", 
    8383               "Thickness of the chain linkage"], 
    84               ["num_pearls", "none", 3, [0, inf], "volume", 
     84              ["num_pearls", "none", 3, [1, inf], "volume", 
    8585               "Number of pearls in the necklace (must be integer)"], 
    8686              ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", 
    … …  
    100100    Redundant with form_volume. 
    101101    """ 
     102    num_pearls = int(num_pearls + 0.5) 
    102103    number_of_strings = num_pearls - 1.0 
    103104    string_vol = edge_sep * pi * pow((thick_string / 2.0), 2.0) 
    … …  
    111112    Calculation for effective radius. 
    112113    """ 
     114    num_pearls = int(num_pearls + 0.5) 
    113115    tot_vol = volume(radius, edge_sep, thick_string, num_pearls) 
    114116    rad_out = (tot_vol/(4.0/3.0*pi)) ** (1./3.) 

  • sasmodels/models/raspberry.py

    r8e68ea0 rc3ccaec  
    1010    Schematic of the raspberry model 
    1111 
    12 In order to calculate the form factor of the entire complex, the self- 
    13 correlation of the large droplet, the self-correlation of the particles, the 
    14 correlation terms between different particles and the cross terms between large 
    15 droplet and small particles all need to be calculated. 
     12In order to calculate the form factor of the entire complex, the 
     13self-correlation of the large droplet, the self-correlation of the particles, 
     14the correlation terms between different particles and the cross terms between 
     15large droplet and small particles all need to be calculated. 
    1616 
    17 Consider two infinitely thin shells of radii R1 and R2 separated by distance r. 
    18 The general structure of the equation is then the form factor of the two shells 
    19 multiplied by the phase factor that accounts for the separation of their 
    20 centers. 
     17Consider two infinitely thin shells of radii $R_1$ and $R_2$ separated by 
     18distance $r$. The general structure of the equation is then the form factor 
     19of the two shells multiplied by the phase factor that accounts for the 
     20separation of their centers. 
    2121 
    2222.. math:: 

  • sasmodels/models/rpa.c

    r6351bfa reb63414  
    1 double Iq(double q, double case_num, 
     1double Iq(double q, double fp_case_num, 
    22    double N[], double Phi[], double v[], double L[], double b[], 
    33    double Kab, double Kac, double Kad, 
    … …  
    55    ); 
    66 
    7 double Iq(double q, double case_num, 
     7double Iq(double q, double fp_case_num, 
    88    double N[],    // DEGREE OF POLYMERIZATION 
    99    double Phi[],  // VOL FRACTION 
    … …  
    1515    ) 
    1616{ 
    17   int icase = (int)case_num; 
     17  int icase = (int)(fp_case_num+0.5); 
    1818 
    1919  double Nab,Nac,Nad,Nbc,Nbd,Ncd; 
    … …  
    3939  double S31,S32,S33,S34,S41,S42,S43,S44; 
    4040  double Lad,Lbd,Lcd,Nav,Intg; 
    41  
     41   
     42  // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc) 
    4243  //icase was shifted to N-1 from the original code 
    4344  if (icase <= 1){ 
    … …  
    5758    Kab = Kac = Kad = -0.0004; 
    5859  } 
    59  
     60  
     61  // Set volume fraction of component D based on constraint that sum of vol frac =1 
     62  Phi[3]=1.0-Phi[0]-Phi[1]-Phi[2]; 
     63 
     64  //set up values for cross terms in case of block copolymers (1,3,4,6,7,8,9) 
    6065  Nab=sqrt(N[0]*N[1]); 
    6166  Nac=sqrt(N[0]*N[2]); 
    … …  
    7984  Phicd=sqrt(Phi[2]*Phi[3]); 
    8085 
     86  // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk  
    8187  Xa=q*q*b[0]*b[0]*N[0]/6.0; 
    8288  Xb=q*q*b[1]*b[1]*N[1]/6.0; 
    … …  
    8490  Xd=q*q*b[3]*b[3]*N[3]/6.0; 
    8591 
    86   Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); 
    87   S0aa=N[0]*Phi[0]*v[0]*Paa; 
    88   Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); 
     92  //calculate all partial structure factors Pij and normalize n^2 
     93  Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); // free A chain form factor 
     94  S0aa=N[0]*Phi[0]*v[0]*Paa; // Phi * Vp * P(Q)= I(Q0)/delRho^2 
     95  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); //AB diblock (anchored Paa * anchored Pbb) partial form factor 
    8996  S0ab=(Phiab*vab*Nab)*Pab; 
    90   Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 
     97  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor 
    9198  S0ac=(Phiac*vac*Nac)*Pac; 
    92   Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 
     99  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block 
    93100  S0ad=(Phiad*vad*Nad)*Pad; 
    94101 
    95102  S0ba=S0ab; 
    96   Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); 
     103  Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain 
    97104  S0bb=N[1]*Phi[1]*v[1]*Pbb; 
    98   Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 
     105  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock 
    99106  S0bc=(Phibc*vbc*Nbc)*Pbc; 
    100   Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 
     107  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock 
    101108  S0bd=(Phibd*vbd*Nbd)*Pbd; 
    102109 
    103110  S0ca=S0ac; 
    104111  S0cb=S0bc; 
    105   Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); 
     112  Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain 
    106113  S0cc=N[2]*Phi[2]*v[2]*Pcc; 
    107   Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 
     114  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock 
    108115  S0cd=(Phicd*vcd*Ncd)*Pcd; 
    109116 
    … …  
    111118  S0db=S0bd; 
    112119  S0dc=S0cd; 
    113   Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); 
     120  Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain 
    114121  S0dd=N[3]*Phi[3]*v[3]*Pdd; 
     122 
     123  // Reset all unused partial structure factors to 0 (depends on case) 
    115124  //icase was shifted to N-1 from the original code 
    116125  switch(icase){ 
    … …  
    193202  S0dc=S0cd; 
    194203 
     204  // self chi parameter is 0 ... of course 
    195205  Kaa=0.0; 
    196206  Kbb=0.0; 
    … …  
    243253  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 
    244254 
     255  // D is considered the matrix or background component so enters here 
    245256  m=1.0/(S0dd-ZZ); 
    246257 
    … …  
    297308  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 
    298309 
     310  //calculate contrast where L[i] is the scattering length of i and D is the matrix 
     311  //Note that should multiply by Nav to get units of SLD which will become 
     312  // Nav*2 in the next line (SLD^2) but then normalization in that line would 
     313  //need to divide by Nav leaving only Nav or sqrt(Nav) before squaring.  
    299314  Nav=6.022045e+23; 
    300315  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); 
    … …  
    303318 
    304319  Intg=Lad*Lad*S11+Lbd*Lbd*S22+Lcd*Lcd*S33+2.0*Lad*Lbd*S12+2.0*Lbd*Lcd*S23+2.0*Lad*Lcd*S13; 
     320 
     321  //rescale for units of Lij^2 (fm^2 to cm^2) 
     322  Intg *= 1.0e-26;     
    305323 
    306324  return Intg; 

  • sasmodels/models/rpa.py

    r40a87fa r4f9e288  
    11r""" 
     2Definition 
     3---------- 
     4 
    25Calculates the macroscopic scattering intensity for a multi-component 
    36homogeneous mixture of polymers using the Random Phase Approximation. 
    … …  
    2427Case 9: A-B-C-D tetra-block copolymer 
    2528 
    26 **NB: these case numbers are different from those in the NIST SANS package!** 
     29.. note:: 
     30    These case numbers are different from those in the NIST SANS package! 
    2731 
    28 Only one case can be used at any one time. 
     32The models are based on the papers by Akcasu et al. [#Akcasu]_ and by 
     33Hammouda [#Hammouda]_ assuming the polymer follows Gaussian statistics such 
     34that $R_g^2 = n b^2/6$ where $b$ is the statistical segment length and $n$ is 
     35the number of statistical segment lengths. A nice tutorial on how these are 
     36constructed and implemented can be found in chapters 28 and 39 of Boualem 
     37Hammouda's 'SANS Toolbox'[#toolbox]_. 
    2938 
    30 The RPA (mean field) formalism only applies only when the multicomponent 
    31 polymer mixture is in the homogeneous mixed-phase region. 
     39In brief the macroscopic cross sections are derived from the general forms 
     40for homopolymer scattering and the multiblock cross-terms while the inter 
     41polymer cross terms are described in the usual way by the $\chi$ parameter. 
    3242 
    33 **Component D is assumed to be the "background" component (ie, all contrasts 
    34 are calculated with respect to component D).** So the scattering contrast 
    35 for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`. 
     43USAGE NOTES: 
    3644 
    37 Depending on which case is being used, the number of fitting parameters - the 
    38 segment lengths (ba, bb, etc) and $\chi$ parameters (Kab, Kac, etc) - vary. 
    39 The *scale* parameter should be held equal to unity. 
     45* Only one case can be used at any one time. 
     46* The RPA (mean field) formalism only applies only when the multicomponent 
     47  polymer mixture is in the homogeneous mixed-phase region. 
     48* **Component D is assumed to be the "background" component (ie, all contrasts 
     49  are calculated with respect to component D).** So the scattering contrast 
     50  for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`. 
     51* Depending on which case is being used, the number of fitting parameters can  
     52  vary. 
    4053 
    41 The input parameters are the degrees of polymerization, the volume fractions, 
    42 the specific volumes, and the neutron scattering length densities for each 
    43 component. 
     54  .. Note:: 
     55    * In general the degrees of polymerization, the volume 
     56      fractions, the molar volumes, and the neutron scattering lengths for each 
     57      component are obtained from other methods and held fixed while The *scale* 
     58      parameter should be held equal to unity. 
     59    * The variables are normally the segment lengths (b\ :sub:`a`, b\ :sub:`b`, 
     60      etc) and $\chi$ parameters (K\ :sub:`ab`, K\ :sub:`ac`, etc). 
    4461 
    4562 
    … …  
    4764---------- 
    4865 
    49 A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 
     66.. [#Akcasu] A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 
     67   4136. 
     68.. [#Hammouda] B. Hammouda, *Advances in Polymer Science* 106 (1993) 87. 
     69.. [#toolbox] https://www.ncnr.nist.gov/staff/hammouda/the_sans_toolbox.pdf 
     70 
     71Authorship and Verification 
     72---------------------------- 
     73 
     74* **Author:** Boualem Hammouda - NIST IGOR/DANSE **Date:** pre 2010 
     75* **Converted to sasmodels by:** Paul Kienzle **Date:** July 18, 2016 
     76* **Last Modified by:** Paul Butler **Date:** March 12, 2017 
     77* **Last Reviewed by:** Paul Butler **Date:** March 12, 2017 
    5078""" 
    5179 
    … …  
    5381 
    5482name = "rpa" 
    55 title = "Random Phase Approximation - unfinished work in progress" 
     83title = "Random Phase Approximation" 
    5684description = """ 
    5785This formalism applies to multicomponent polymer mixtures in the 
    … …  
    90118    ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    91119    ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 
    92     ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
     120    ["v[4]", "mL/mol", 100.0, [0, inf], "", "molar volume"], 
    93121    ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    94122    ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 
    … …  
    107135 
    108136control = "case_num" 
    109 HIDE_NONE = set() 
    110 HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()) 
     137HIDE_ALL = set("Phi4".split()) 
     138HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()).union(HIDE_ALL) 
    111139HIDE_AB = set("N2 Phi2 v2 L2 b2 K23 K24".split()).union(HIDE_A) 
    112140def hidden(case_num): 
    … …  
    114142    Return a list of parameters to hide depending on the multiplicity parameter. 
    115143    """ 
     144    case_num = int(case_num+0.5) 
    116145    if case_num < 2: 
    117146        return HIDE_AB 
    … …  
    119148        return HIDE_A 
    120149    else: 
    121         return HIDE_NONE 
     150        return HIDE_ALL 
    122151 

  • sasmodels/models/spherical_sld.c

    r925ad6e rc3ccaec  
    11static double form_volume( 
    2     int n_shells, 
     2    double fp_n_shells, 
    33    double thickness[], 
    44    double interface[]) 
    55{ 
     6    int n_shells= (int)(fp_n_shells + 0.5); 
    67    double r = 0.0; 
    78    for (int i=0; i < n_shells; i++) { 
    … …  
    2021        return pow(z, nu); 
    2122    } else if (shape==2) { 
    22         return 1.0 - pow(1. - z, nu); 
     23        return 1.0 - pow(1.0 - z, nu); 
    2324    } else if (shape==3) { 
    2425        return expm1(-nu*z)/expm1(-nu); 
    … …  
    4445static double Iq( 
    4546    double q, 
    46     int n_shells, 
     47    double fp_n_shells, 
    4748    double sld_solvent, 
    4849    double sld[], 
    … …  
    5152    double shape[], 
    5253    double nu[], 
    53     int n_steps) 
     54    double fp_n_steps) 
    5455{ 
    5556    // iteration for # of shells + core + solvent 
     57    int n_shells = (int)(fp_n_shells + 0.5); 
     58    int n_steps = (int)(fp_n_steps + 0.5); 
    5659    double f=0.0; 
    5760    double r=0.0; 

  • sasmodels/models/spherical_sld.py

    r925ad6e rc3ccaec  
    233233    """ 
    234234 
     235    n_shells = int(n_shells + 0.5) 
     236    n_steps = int(n_steps + 0.5) 
    235237    z = [] 
    236238    rho = [] 
    … …  
    240242    rho.append(sld[0]) 
    241243 
    242     for i in range(0, int(n_shells)): 
     244    for i in range(0, n_shells): 
    243245        z_next += thickness[i] 
    244246        z.append(z_next) 
    … …  
    261263def ER(n_shells, thickness, interface): 
    262264    """Effective radius""" 
    263     n_shells = int(n_shells) 
     265    n_shells = int(n_shells + 0.5) 
    264266    total = (np.sum(thickness[:n_shells], axis=1) 
    265267             + np.sum(interface[:n_shells], axis=1)) 

  • sasmodels/models/stacked_disks.c

    r6c3e266 r19f996b  
    1 double form_volume(double thick_core, 
    2                    double thick_layer, 
    3                    double radius, 
    4                    double n_stacking); 
    5  
    6 double Iq(double q, 
    7           double thick_core, 
    8           double thick_layer, 
    9           double radius, 
    10           double n_stacking, 
    11           double sigma_dnn, 
    12           double core_sld, 
    13           double layer_sld, 
    14           double solvent_sld); 
    15  
    16 double Iqxy(double qx, double qy, 
    17           double thick_core, 
    18           double thick_layer, 
    19           double radius, 
    20           double n_stacking, 
    21           double sigma_dnn, 
    22           double core_sld, 
    23           double layer_sld, 
    24           double solvent_sld, 
    25           double theta, 
    26           double phi); 
    27  
    28 static 
    29 double _kernel(double q, 
    30                double radius, 
    31                double core_sld, 
    32                double layer_sld, 
    33                double solvent_sld, 
    34                double halfheight, 
    35                double thick_layer, 
    36                double sin_alpha, 
    37                double cos_alpha, 
    38                double sigma_dnn, 
    39                double d, 
    40                double n_stacking) 
     1static double stacked_disks_kernel( 
     2    double q, 
     3    double halfheight, 
     4    double thick_layer, 
     5    double radius, 
     6    int n_stacking, 
     7    double sigma_dnn, 
     8    double core_sld, 
     9    double layer_sld, 
     10    double solvent_sld, 
     11    double sin_alpha, 
     12    double cos_alpha, 
     13    double d) 
    4114 
    4215{ 
    … …  
    8861 
    8962 
    90 static 
    91 double stacked_disks_kernel(double q, 
    92                             double thick_core, 
    93                             double thick_layer, 
    94                             double radius, 
    95                             double n_stacking, 
    96                             double sigma_dnn, 
    97                             double core_sld, 
    98                             double layer_sld, 
    99                             double solvent_sld) 
     63static double stacked_disks_1d( 
     64    double q, 
     65    double thick_core, 
     66    double thick_layer, 
     67    double radius, 
     68    int n_stacking, 
     69    double sigma_dnn, 
     70    double core_sld, 
     71    double layer_sld, 
     72    double solvent_sld) 
    10073{ 
    10174/*    StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks 
    … …  
    11184        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    11285        SINCOS(zi, sin_alpha, cos_alpha); 
    113         double yyy = _kernel(q, 
     86        double yyy = stacked_disks_kernel(q, 
     87                           halfheight, 
     88                           thick_layer, 
    11489                           radius, 
     90                           n_stacking, 
     91                           sigma_dnn, 
    11592                           core_sld, 
    11693                           layer_sld, 
    11794                           solvent_sld, 
    118                            halfheight, 
    119                            thick_layer, 
    12095                           sin_alpha, 
    12196                           cos_alpha, 
    122                            sigma_dnn, 
    123                            d, 
    124                            n_stacking); 
     97                           d); 
    12598        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    12699    } 
    … …  
    132105} 
    133106 
    134 double form_volume(double thick_core, 
    135                    double thick_layer, 
    136                    double radius, 
    137                    double n_stacking){ 
     107static double form_volume( 
     108    double thick_core, 
     109    double thick_layer, 
     110    double radius, 
     111    double fp_n_stacking) 
     112{ 
     113    int n_stacking = (int)(fp_n_stacking + 0.5); 
    138114    double d = 2.0 * thick_layer + thick_core; 
    139115    return M_PI * radius * radius * d * n_stacking; 
    140116} 
    141117 
    142 double Iq(double q, 
    143           double thick_core, 
    144           double thick_layer, 
    145           double radius, 
    146           double n_stacking, 
    147           double sigma_dnn, 
    148           double core_sld, 
    149           double layer_sld, 
    150           double solvent_sld) 
     118static double Iq( 
     119    double q, 
     120    double thick_core, 
     121    double thick_layer, 
     122    double radius, 
     123    double fp_n_stacking, 
     124    double sigma_dnn, 
     125    double core_sld, 
     126    double layer_sld, 
     127    double solvent_sld) 
    151128{ 
    152     return stacked_disks_kernel(q, 
     129    int n_stacking = (int)(fp_n_stacking + 0.5); 
     130    return stacked_disks_1d(q, 
    153131                    thick_core, 
    154132                    thick_layer, 
    … …  
    162140 
    163141 
    164 double 
    165 Iqxy(double qx, double qy, 
    166      double thick_core, 
    167      double thick_layer, 
    168      double radius, 
    169      double n_stacking, 
    170      double sigma_dnn, 
    171      double core_sld, 
    172      double layer_sld, 
    173      double solvent_sld, 
    174      double theta, 
    175      double phi) 
     142static double Iqxy(double qx, double qy, 
     143    double thick_core, 
     144    double thick_layer, 
     145    double radius, 
     146    double fp_n_stacking, 
     147    double sigma_dnn, 
     148    double core_sld, 
     149    double layer_sld, 
     150    double solvent_sld, 
     151    double theta, 
     152    double phi) 
    176153{ 
     154    int n_stacking = (int)(fp_n_stacking + 0.5); 
    177155    double q, sin_alpha, cos_alpha; 
    178156    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    … …  
    180158    double d = 2.0 * thick_layer + thick_core; 
    181159    double halfheight = 0.5*thick_core; 
    182     double answer = _kernel(q, 
     160    double answer = stacked_disks_kernel(q, 
     161                     halfheight, 
     162                     thick_layer, 
    183163                     radius, 
     164                     n_stacking, 
     165                     sigma_dnn, 
    184166                     core_sld, 
    185167                     layer_sld, 
    186168                     solvent_sld, 
    187                      halfheight, 
    188                      thick_layer, 
    189169                     sin_alpha, 
    190170                     cos_alpha, 
    191                      sigma_dnn, 
    192                      d, 
    193                      n_stacking); 
     171                     d); 
    194172 
    195173    //convert to [cm-1] 

  • sasmodels/models/stacked_disks.py

    rb7e8b94 rc3ccaec  
    126126    ["thick_layer", "Ang",        10.0, [0, inf],    "volume",      "Thickness of layer each side of core"], 
    127127    ["radius",      "Ang",        15.0, [0, inf],    "volume",      "Radius of the stacked disk"], 
    128     ["n_stacking",  "",            1.0, [0, inf],    "volume",      "Number of stacked layer/core/layer disks"], 
     128    ["n_stacking",  "",            1.0, [1, inf],    "volume",      "Number of stacked layer/core/layer disks"], 
    129129    ["sigma_d",     "Ang",         0,   [0, inf],    "",            "Sigma of nearest neighbor spacing"], 
    130130    ["sld_core",    "1e-6/Ang^2",  4,   [-inf, inf], "sld",         "Core scattering length density"], 

  • sasmodels/models/star_polymer.c

    r3a48772 r2586093f  
    33double Iq(double q, double radius2, double arms); 
    44 
    5 static double _mass_fractal_kernel(double q, double radius2, double arms) 
     5static double star_polymer_kernel(double q, double radius2, double arms) 
    66{ 
    77 
    … …  
    2323double Iq(double q, double radius2, double arms) 
    2424{ 
    25     return _mass_fractal_kernel(q, radius2, arms); 
     25    return star_polymer_kernel(q, radius2, arms); 
    2626} 

  • sasmodels/models/unified_power_Rg.py

    r66ca2a6 rc3ccaec  
    9797 
    9898def Iq(q, level, rg, power, B, G): 
    99     ilevel = int(level) 
    100     if ilevel == 0: 
     99    level = int(level + 0.5) 
     100    if level == 0: 
    101101        with errstate(divide='ignore'): 
    102102            return 1./q 
    … …  
    104104    with errstate(divide='ignore', invalid='ignore'): 
    105105        result = np.zeros(q.shape, 'd') 
    106         for i in range(ilevel): 
     106        for i in range(level): 
    107107            exp_now = exp(-(q*rg[i])**2/3.) 
    108108            pow_now = (erf(q*rg[i]/sqrt(6.))**3/q)**power[i] 
    109             if i < ilevel-1: 
     109            if i < level-1: 
    110110                exp_next = exp(-(q*rg[i+1])**2/3.) 
    111111            else: 
    … …  
    113113            result += G[i]*exp_now + B[i]*exp_next*pow_now 
    114114 
    115     result[q == 0] = np.sum(G[:ilevel]) 
     115    result[q == 0] = np.sum(G[:level]) 
    116116    return result 
    117117 

  • sasmodels/product.py

    r9951a86 rf88e248  
    4545    # structure factor calculator.  Structure factors should not 
    4646    # have any magnetic parameters 
    47     assert(s_info.parameters.kernel_parameters[0].id == ER_ID) 
    48     assert(s_info.parameters.kernel_parameters[1].id == VF_ID) 
    49     assert(s_info.parameters.magnetism_index == []) 
     47    if not s_info.parameters.kernel_parameters[0].id == ER_ID: 
     48        raise TypeError("S needs %s as first parameter"%ER_ID) 
     49    if not s_info.parameters.kernel_parameters[1].id == VF_ID: 
     50        raise TypeError("S needs %s as second parameter"%VF_ID) 
     51    if not s_info.parameters.magnetism_index == []: 
     52        raise TypeError("S should not have SLD parameters") 
    5053    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters 
    5154    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
    52     p_set = set(p.id for p in p_pars.call_parameters) 
    53     s_set = set(p.id for p in s_pars.call_parameters) 
    54  
    55     if p_set & s_set: 
    56         # there is some overlap between the parameter names; tag the 
    57         # overlapping S parameters with name_S. 
    58         # Skip the first parameter of s, which is effective radius 
    59         s_list = [(suffix_parameter(par) if par.id in p_set else par) 
    60                   for par in s_pars.kernel_parameters[1:]] 
    61     else: 
    62         # Skip the first parameter of s, which is effective radius 
    63         s_list = s_pars.kernel_parameters[1:] 
     55 
     56    # Create list of parameters for the combined model.  Skip the first 
     57    # parameter of S, which we verified above is effective radius.  If there 
     58    # are any names in P that overlap with those in S, modify the name in S 
     59    # to distinguish it. 
     60    p_set = set(p.id for p in p_pars.kernel_parameters) 
     61    s_list = [(_tag_parameter(par) if par.id in p_set else par) 
     62              for par in s_pars.kernel_parameters[1:]] 
     63    # Check if still a collision after renaming.  This could happen if for 
     64    # example S has volfrac and P has both volfrac and volfrac_S. 
     65    if any(p.id in p_set for p in s_list): 
     66        raise TypeError("name collision: P has P.name and P.name_S while S has S.name") 
     67 
    6468    translate_name = dict((old.id, new.id) for old, new 
    6569                          in zip(s_pars.kernel_parameters[1:], s_list)) 
    6670    demo = {} 
    67     demo.update((k, v) for k, v in p_info.demo.items() 
    68                 if k not in ("background", "scale")) 
     71    demo.update(p_info.demo.items()) 
    6972    demo.update((translate_name[k], v) for k, v in s_info.demo.items() 
    7073                if k not in ("background", "scale") and not k.startswith(ER_ID)) 
    … …  
    9093    # Remember the component info blocks so we can build the model 
    9194    model_info.composition = ('product', [p_info, s_info]) 
    92     model_info.demo = {} 
     95    model_info.demo = demo 
     96 
     97    ## Show the parameter table with the demo values 
     98    #from .compare import get_pars, parlist 
     99    #print("==== %s ====="%model_info.name) 
     100    #values = get_pars(model_info, use_demo=True) 
     101    #print(parlist(model_info, values, is2d=True)) 
    93102    return model_info 
    94103 
    95 def suffix_parameter(par, suffix): 
     104def _tag_parameter(par): 
     105    """ 
     106    Tag the parameter name with _S to indicate that the parameter comes from 
     107    the structure factor parameter set.  This is only necessary if the 
     108    form factor model includes a parameter of the same name as a parameter 
     109    in the structure factor. 
     110    """ 
    96111    par = copy(par) 
    97     par.name = par.name + " S" 
     112    # Protect against a vector parameter in S by appending the vector length 
     113    # to the renamed parameter.  Note: haven't tested this since no existing 
     114    # structure factor models contain vector parameters. 
     115    vector_length = par.name[len(par.id):] 
    98116    par.id = par.id + "_S" 
     117    par.name = par.id + vector_length 
    99118    return par 
    100119 

  • sasmodels/resolution.py

    rb397165 rb32caab  
    1717 
    1818MINIMUM_RESOLUTION = 1e-8 
    19  
    20  
    21 # When extrapolating to -q, what is the minimum positive q relative to q_min 
    22 # that we wish to calculate? 
    23 MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION = 0.01 
     19MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data 
    2420 
    2521class Resolution(object): 
    … …  
    8278        self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma) 
    8379                       if q_calc is None else np.sort(q_calc)) 
     80 
     81        # Protect against models which are not defined for very low q.  Limit 
     82        # the smallest q value evaluated (in absolute) to 0.02*min 
     83        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
     84        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff] 
     85 
     86        # Build weight matrix from calculated q values 
    8487        self.weight_matrix = pinhole_resolution(self.q_calc, self.q, 
    8588                                np.maximum(q_width, MINIMUM_RESOLUTION)) 
     89        self.q_calc = abs(self.q_calc) 
    8690 
    8791    def apply(self, theory): 
    … …  
    123127        self.q_calc = slit_extend_q(q, qx_width, qy_width) \ 
    124128            if q_calc is None else np.sort(q_calc) 
     129 
     130        # Protect against models which are not defined for very low q.  Limit 
     131        # the smallest q value evaluated (in absolute) to 0.02*min 
     132        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
     133        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff] 
     134 
     135        # Build weight matrix from calculated q values 
    125136        self.weight_matrix = \ 
    126137            slit_resolution(self.q_calc, self.q, qx_width, qy_width) 
     138        self.q_calc = abs(self.q_calc) 
    127139 
    128140    def apply(self, theory): 
    … …  
    153165    # neither trapezoid nor Simpson's rule improved the accuracy. 
    154166    edges = bin_edges(q_calc) 
    155     edges[edges < 0.0] = 0.0 # clip edges below zero 
     167    #edges[edges < 0.0] = 0.0 # clip edges below zero 
    156168    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    157169    weights = cdf[1:] - cdf[:-1] 
    … …  
    286298    # The current algorithm is a midpoint rectangle rule. 
    287299    q_edges = bin_edges(q_calc) # Note: requires q > 0 
    288     q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 
     300    #q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 
    289301    weights = np.zeros((len(q), len(q_calc)), 'd') 
    290302 
    … …  
    392404    interval. 
    393405 
    394     if *q_min* is zero or less then *q[0]/10* is used instead. 
     406    Note that extrapolated values may be negative. 
    395407    """ 
    396408    q = np.sort(q) 
    397409    if q_min + 2*MINIMUM_RESOLUTION < q[0]: 
    398         if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    399410        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15 
    400411        q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 
    … …  
    448459        log_delta_q = log(10.) / points_per_decade 
    449460    if q_min < q[0]: 
    450         if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
     461        if q_min < 0: q_min = q[0]*MINIMUM_ABSOLUTE_Q 
    451462        n_low = log_delta_q * (log(q[0])-log(q_min)) 
    452463        q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 

  • sasmodels/sasview_model.py

    rfe8ff99 r749a7d4  
    1616import logging 
    1717from os.path import basename, splitext 
     18import thread 
    1819 
    1920import numpy as np  # type: ignore 
    … …  
    3839except ImportError: 
    3940    pass 
     41 
     42calculation_lock = thread.allocate_lock() 
    4043 
    4144SUPPORT_OLD_STYLE_PLUGINS = True 
    … …  
    605608        to the card for each evaluation. 
    606609        """ 
     610        ## uncomment the following when trying to debug the uncoordinated calls 
     611        ## to calculate_Iq 
     612        #if calculation_lock.locked(): 
     613        #    logging.info("calculation waiting for another thread to complete") 
     614        #    logging.info("\n".join(traceback.format_stack())) 
     615 
     616        with calculation_lock: 
     617            return self._calculate_Iq(qx, qy) 
     618 
     619    def _calculate_Iq(self, qx, qy=None): 
    607620        #core.HAVE_OPENCL = False 
    608621        if self._model is None: 
    … …  
    721734            return [self.params[par.name]], [1.0] 
    722735 
    723 def test_model(): 
     736def test_cylinder(): 
    724737    # type: () -> float 
    725738    """ 
    726     Test that a sasview model (cylinder) can be run. 
     739    Test that the cylinder model runs, returning the value at [0.1,0.1]. 
    727740    """ 
    728741    Cylinder = _make_standard_model('cylinder') 
    … …  
    733746    # type: () -> float 
    734747    """ 
    735     Test that a sasview model (cylinder) can be run. 
     748    Test that 2-D hardsphere model runs and doesn't produce NaN. 
    736749    """ 
    737750    Model = _make_standard_model('hardsphere') 
    … …  
    744757    # type: () -> float 
    745758    """ 
    746     Test that a sasview model (cylinder) can be run. 
     759    Test that the 2-D RPA model runs 
    747760    """ 
    748761    RPA = _make_standard_model('rpa') 
    … …  
    750763    return rpa.evalDistribution([0.1, 0.1]) 
    751764 
     765def test_empty_distribution(): 
     766    # type: () -> None 
     767    """ 
     768    Make sure that sasmodels returns NaN when there are no polydispersity points 
     769    """ 
     770    Cylinder = _make_standard_model('cylinder') 
     771    cylinder = Cylinder() 
     772    cylinder.setParam('radius', -1.0) 
     773    cylinder.setParam('background', 0.) 
     774    Iq = cylinder.evalDistribution(np.asarray([0.1])) 
     775    assert np.isnan(Iq[0]), "empty distribution fails" 
    752776 
    753777def test_model_list(): 
    754778    # type: () -> None 
    755779    """ 
    756     Make sure that all models build as sasview models. 
     780    Make sure that all models build as sasview models 
    757781    """ 
    758782    from .exception import annotate_exception 
    … …  
    781805 
    782806if __name__ == "__main__": 
    783     print("cylinder(0.1,0.1)=%g"%test_model()) 
     807    print("cylinder(0.1,0.1)=%g"%test_cylinder()) 
     808    #test_empty_distribution() 

  • sasmodels/sesans.py

    rb397165 r94d13f1  
    1414import numpy as np  # type: ignore 
    1515from numpy import pi, exp  # type: ignore 
    16 from scipy.special import jv as besselj 
    17 #import direct_model.DataMixin as model 
    18          
    19 def make_q(q_max, Rmax): 
    20     r""" 
    21     Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ 
    22     to $q_max$. This is the integration range of the Hankel transform; bigger range and  
    23     more points makes a better numerical integration. 
    24     Smaller q_min will increase reliable spin echo length range.  
    25     Rmax is the "radius" of the largest expected object and can be set elsewhere. 
    26     q_max is determined by the acceptance angle of the SESANS instrument. 
     16from scipy.special import j0 
     17 
     18class SesansTransform(object): 
    2719    """ 
    28     from sas.sascalc.data_util.nxsunit import Converter 
     20    Spin-Echo SANS transform calculator.  Similar to a resolution function, 
     21    the SesansTransform object takes I(q) for the set of *q_calc* values and 
     22    produces a transformed dataset 
    2923 
    30     q_min = dq = 0.1 * 2*pi / Rmax 
    31     return np.arange(q_min, 
    32                      Converter(q_max[1])(q_max[0], 
    33                                          units="1/A"), 
    34                      dq) 
    35      
    36 def make_all_q(data): 
     24    *SElength* (A) is the set of spin-echo lengths in the measured data. 
     25 
     26    *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin 
     27    echo encoding dimension (for ToF: Q of min(R) and max(lam)). 
     28 
     29    *Rmax* (A) is the maximum size sensitivity; larger radius requires more 
     30    computation time. 
    3731    """ 
    38     Return a $q$ vector suitable for calculating the total scattering cross section for 
    39     calculating the effect of finite acceptance angles on Time of Flight SESANS instruments. 
    40     If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed. 
    41     If the instrument has a rectangular acceptance, 2 all_q vectors are needed. 
    42     If the instrument has a circular acceptance, 1 all_q vector is needed 
    43      
    44     """ 
    45     if not data.has_no_finite_acceptance: 
    46         return [] 
    47     elif data.has_yz_acceptance(data): 
    48         # compute qx, qy 
    49         Qx, Qy = np.meshgrid(qx, qy) 
    50         return [Qx, Qy] 
    51     else: 
    52         # else only need q 
    53         # data.has_z_acceptance 
    54         return [q] 
     32    #: SElength from the data in the original data units; not used by transform 
     33    #: but the GUI uses it, so make sure that it is present. 
     34    q = None  # type: np.ndarray 
    5535 
    56 def transform(data, q_calc, Iq_calc, qmono, Iq_mono): 
    57     """ 
    58     Decides which transform type is to be used, based on the experiment data file contents (header) 
    59     (2016-03-19: currently controlled from parameters script) 
    60     nqmono is the number of q vectors to be used for the detector integration 
    61     """ 
    62     nqmono = len(qmono) 
    63     if nqmono == 0: 
    64         result = call_hankel(data, q_calc, Iq_calc) 
    65     elif nqmono == 1: 
    66         q = qmono[0] 
    67         result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) 
    68     else: 
    69         Qx, Qy = [qmono[0], qmono[1]] 
    70         Qx = np.reshape(Qx, nqx, nqy) 
    71         Qy = np.reshape(Qy, nqx, nqy) 
    72         Iq_mono = np.reshape(Iq_mono, nqx, nqy) 
    73         qx = Qx[0, :] 
    74         qy = Qy[:, 0] 
    75         result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) 
     36    #: q values to calculate when computing transform 
     37    q_calc = None  # type: np.ndarray 
    7638 
    77     return result 
     39    # transform arrays 
     40    _H = None  # type: np.ndarray 
     41    _H0 = None # type: np.ndarray 
    7842 
    79 def call_hankel(data, q_calc, Iq_calc): 
    80     return hankel((data.x, data.x_unit), 
    81                   (data.lam, data.lam_unit), 
    82                   (data.sample.thickness, 
    83                    data.sample.thickness_unit), 
    84                   q_calc, Iq_calc) 
    85    
    86 def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): 
    87     return hankel(data.x, data.lam * 1e-9, 
    88                   data.sample.thickness / 10, 
    89                   q_calc, Iq_calc) 
    90                    
    91 def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): 
    92     return hankel(data.x, data.y, data.lam * 1e-9, 
    93                   data.sample.thickness / 10, 
    94                   q_calc, Iq_calc) 
    95                          
    96 def TotalScatter(model, parameters):  #Work in progress!! 
    97 #    Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) 
    98     allq = np.linspace(0,4*pi/wavelength,1000) 
    99     allIq = 1 
    100     integral = allq*allIq 
    101      
     43    def __init__(self, z, SElength, zaccept, Rmax): 
     44        # type: (np.ndarray, float, float) -> None 
     45        #import logging; logging.info("creating SESANS transform") 
     46        self.q = z 
     47        self._set_hankel(SElength, zaccept, Rmax) 
    10248 
     49    def apply(self, Iq): 
     50        # tye: (np.ndarray) -> np.ndarray 
     51        G0 = np.dot(self._H0, Iq) 
     52        G = np.dot(self._H.T, Iq) 
     53        P = G - G0 
     54        return P 
    10355 
    104 def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still 
    105 #============================================================================== 
    106 #     2D Cosine Transform if "wavelength" is a vector 
    107 #============================================================================== 
    108 #allq is the q-space needed to create the total scattering cross-section 
     56    def _set_hankel(self, SElength, zaccept, Rmax): 
     57        # type: (np.ndarray, float, float) -> None 
     58        # Force float32 arrays, otherwise run into memory problems on some machines 
     59        SElength = np.asarray(SElength, dtype='float32') 
    10960 
    110     Gprime = np.zeros_like(wavelength, 'd') 
    111     s = np.zeros_like(wavelength, 'd') 
    112     sd = np.zeros_like(wavelength, 'd') 
    113     Gprime = np.zeros_like(wavelength, 'd') 
    114     f = np.zeros_like(wavelength, 'd') 
    115     for i, wavelength_i in enumerate(wavelength): 
    116         z = magfield*wavelength_i 
    117         allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
    118         allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
    119         alldq = (allq[1]-allq[0])*1e10 
    120         sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
    121         s[i]=1-exp(-sigma) 
    122         for j, Iqy_j, qy_j in enumerate(qy): 
    123             for k, Iqz_k, qz_k in enumerate(qz): 
    124                 Iq = np.sqrt(Iqy_j^2+Iqz_k^2) 
    125                 q = np.sqrt(qy_j^2 + qz_k^2) 
    126                 Gintegral = Iq*cos(z*Qz_k) 
    127                 Gprime[i] += Gintegral 
    128 #                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] 
    129 #                s[i] += 1-exp(Totalscatter(modelname)*thickness) 
    130 #                For now, work with standard 2-phase scatter 
     61        #Rmax = #value in text box somewhere in FitPage? 
     62        q_max = 2*pi / (SElength[1] - SElength[0]) 
     63        q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1]) 
     64        q = np.arange(q_min, q_max, q_min, dtype='float32') 
     65        dq = q_min 
    13166 
     67        H0 = np.float32(dq/(2*pi)) * q 
    13268 
    133                 sd[i] += Iq 
    134         f[i] = 1-s[i]+sd[i] 
    135         P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] 
     69        repq = np.tile(q, (SElength.size, 1)).T 
     70        repSE = np.tile(SElength, (q.size, 1)) 
     71        H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq 
    13672 
    137  
    138  
    139  
    140 def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): 
    141 #============================================================================== 
    142 #     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS 
    143 #============================================================================== 
    144 #acceptq is the q-space needed to create limited acceptance effect 
    145     SElength= wavelength*magfield 
    146     G = np.zeros_like(SElength, 'd') 
    147     threshold=2*pi*theta/wavelength 
    148     for i, SElength_i in enumerate(SElength): 
    149         allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
    150         allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
    151         alldq = (allq[1]-allq[0])*1e10 
    152         sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
    153         s[i]=1-exp(-sigma) 
    154  
    155         dq = (q[1]-q[0])*1e10 
    156         a = (x<threshold) 
    157         acceptq = a*q 
    158         acceptIq = a*Iq 
    159  
    160         G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) 
    161  
    162 #        G[i]=np.sum(integral) 
    163  
    164     G *= dq*1e10*2*pi 
    165  
    166     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 
    167      
    168 def hankel(SElength, wavelength, thickness, q, Iq): 
    169     r""" 
    170     Compute the expected SESANS polarization for a given SANS pattern. 
    171  
    172     Uses the hankel transform followed by the exponential.  The values for *zz* 
    173     (or spin echo length, or delta), wavelength and sample thickness should 
    174     come from the dataset.  $q$ should be chosen such that the oscillations 
    175     in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). 
    176  
    177     *SElength* [A] is the set of $z$ points at which to compute the 
    178     Hankel transform 
    179  
    180     *wavelength* [m]  is the wavelength of each individual point *zz* 
    181  
    182     *thickness* [cm] is the sample thickness. 
    183  
    184     *q* [A$^{-1}$] is the set of $q$ points at which the model has been 
    185     computed. These should be equally spaced. 
    186  
    187     *I* [cm$^{-1}$] is the value of the SANS model at *q* 
    188     """ 
    189  
    190     from sas.sascalc.data_util.nxsunit import Converter 
    191     wavelength = Converter(wavelength[1])(wavelength[0],"A") 
    192     thickness = Converter(thickness[1])(thickness[0],"A") 
    193     Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters 
    194     SElength = Converter(SElength[1])(SElength[0],"A") 
    195  
    196     G = np.zeros_like(SElength, 'd') 
    197 #============================================================================== 
    198 #     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS 
    199 #============================================================================== 
    200     for i, SElength_i in enumerate(SElength): 
    201         integral = besselj(0, q*SElength_i)*Iq*q 
    202         G[i] = np.sum(integral) 
    203     G0 = np.sum(Iq*q) 
    204  
    205     # [m^-1] step size in q, needed for integration 
    206     dq = (q[1]-q[0]) 
    207  
    208     # integration step, convert q into [m**-1] and 2 pi circle integration 
    209     G *= dq*2*pi 
    210     G0 = np.sum(Iq*q)*dq*2*np.pi 
    211  
    212     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0)) 
    213  
    214     return P 
     73        self.q_calc = q 
     74        self._H, self._H0 = H, H0 
Note: See TracChangeset for help on using the changeset viewer.