Changeset 37c7d5e in sasmodels


Ignore:
Timestamp:
Mar 13, 2017 4:19:34 AM (8 years ago)
Author:
wojciech
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
352494a
Parents:
ae32bb8 (diff), 5d23de2 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' of ​https://github.com/SasView/sasmodels

Files:
35 edited

Legend:

Unmodified
Added
Removed

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

    rfe25eda rf72d70a  
    322322        name = name.split('*')[0] 
    323323 
     324    # Suppress magnetism for python models (not yet implemented) 
     325    if callable(model_info.Iq): 
     326        pars.update(suppress_magnetism(pars)) 
     327 
    324328    if name == 'barbell': 
    325329        if pars['radius_bell'] < pars['radius']: 
    … …  
    340344        if pars['radius'] < pars['thick_string']: 
    341345            pars['radius'], pars['thick_string'] = pars['thick_string'], pars['radius'] 
    342         pars['num_pearls'] = math.ceil(pars['num_pearls']) 
    343346        pass 
    344347 
    … …  
    353356        for c in '1234': 
    354357            pars['Phi'+c] /= total 
    355  
    356     elif name == 'stacked_disks': 
    357         pars['n_stacking'] = math.ceil(pars['n_stacking']) 
    358358 
    359359def parlist(model_info, pars, is2d): 

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

  • doc/ref/gpu/gpu_computations.rst

    r3f5a566 r7e74ed5  
    3131from available OpenCL platforms. 
    3232 
     33OpenCL devices can be set from OpenCL options dialog in Fitting menu or as 
     34enviromental variables. 
     35 
     36**If you don't want to use OpenCL, you can select "No OpenCL" option from** 
     37**GUI dialog or set *SAS_OPENCL=None* in your environment settings** 
     38**This will only use normal programs.** 
     39 
    3340SasView prefers AMD or NVIDIA drivers for GPU, and prefers Intel or 
    3441Apple drivers for CPU. Both GPU and CPU are included on the assumption that CPU  
    … …  
    3946chose to run the model. 
    4047 
    41 **If you don't want to use OpenCL, you can set** *SAS_OPENCL=None* 
    42 **in your environment settings, and it will only use normal programs.** 
    43  
    44 If you want to use one of the other devices, you can run the following 
    45 from the python console in SasView:: 
     48If you want to use one of the other devices, you can select it from OpenCL 
     49options dialog (accessible from Fitting menu) or run the following from 
     50the python console in SasView:: 
    4651 
    4752    import pyopencl as cl 
Note: See TracChangeset for help on using the changeset viewer.