Changes in / [6cff473:a769b54] in sasmodels


Ignore:
Files:
1 deleted
43 edited

Legend:

Unmodified
Added
Removed

  • doc/ref/index.rst

    r9f12fbe rc34a31f  
    1010   refs.rst 
    1111   gpu/gpu_computations.rst 
    12    gpu/opencl_installation.rst 
    1312   magnetism/magnetism.rst 
    1413   sesans/sans_to_sesans.rst 

  • example/fit.py

    rf4b36fa r1182da5  
    2424            % section) 
    2525data = radial_data if section != "tangential" else tan_data 
    26 theta = 89.9 if section != "tangential" else 0 
    27 phi = 90 
     26phi = 0 if section != "tangential" else 90 
    2827kernel = load_model(name, dtype="single") 
    2928cutoff = 1e-3 
    … …  
    3130if name == "ellipsoid": 
    3231    model = Model(kernel, 
    33         scale=0.08, background=35, 
    34         radius_polar=15, radius_equatorial=800, 
     32        scale=0.08, 
     33        r_polar=15, r_equatorial=800, 
    3534        sld=.291, sld_solvent=7.105, 
    36         theta=theta, phi=phi, 
    37         theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3, 
     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, 
    3840        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 
    4344    # SET THE FITTING PARAMETERS 
    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) 
     45    model.r_polar.range(15, 1000) 
     46    model.r_equatorial.range(15, 1000) 
     47    model.theta_pd.range(0, 360) 
    5048    model.background.range(0,1000) 
    5149    model.scale.range(0, 10) 
    5250 
    5351 
     52 
    5453elif name == "lamellar": 
    5554    model = Model(kernel, 
    56         scale=0.08, background=0.003, 
     55        scale=0.08, 
    5756        thickness=19.2946, 
    5857        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 
    6162 
    6263    # SET THE FITTING PARAMETERS 
    … …  
    7677        radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3, 
    7778        length_pd=0.493, length_pd_n=10, length_pd_nsigma=3, 
    78         phi_pd=0, phi_pd_n=5 phi_pd_nsigma=3,) 
     79        phi_pd=0, phi_pd_n=5, phi_pd_nsigma=3,) 
    7980        """ 
    8081    pars = dict( 
    8182        scale=.01, background=35, 
    8283        sld=.291, sld_solvent=5.77, 
    83         radius=250, length=178, 
     84        radius=250, length=178,  
     85        theta=90, phi=phi, 
    8486        radius_pd=0.1, radius_pd_n=5, radius_pd_nsigma=3, 
    8587        length_pd=0.1,length_pd_n=5, length_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) 
     88        theta_pd=10, theta_pd_n=50, theta_pd_nsigma=3, 
     89        phi_pd=0, phi_pd_n=10, phi_pd_nsigma=3) 
    8990    model = Model(kernel, **pars) 
    9091 
    … …  
    9293    model.radius.range(1, 500) 
    9394    model.length.range(1, 5000) 
    94     #model.theta.range(0, 90) 
    95     model.phi.range(0, 180) 
    96     model.phi_pd.range(0, 30) 
     95    model.theta.range(-90,100) 
     96    model.theta_pd.range(0, 30) 
     97    model.theta_pd_n = model.theta_pd + 5 
    9798    model.radius_pd.range(0, 1) 
    98     model.length_pd.range(0, 1) 
     99    model.length_pd.range(0, 2) 
    99100    model.scale.range(0, 10) 
    100101    model.background.range(0, 100) 
    … …  
    103104elif name == "core_shell_cylinder": 
    104105    model = Model(kernel, 
    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, 
     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 
    108110        radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3, 
    109111        length_pd=0.26, length_pd_n=10, length_pd_nsigma=3, 
    110112        thickness_pd=1, thickness_pd_n=1, thickness_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, 
     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, 
    114115        ) 
    115116 
    116117    # SET THE FITTING PARAMETERS 
    117     model.radius.range(115, 1000) 
    118     model.length.range(0, 2500) 
     118    #model.radius.range(115, 1000) 
     119    #model.length.range(0, 2500) 
    119120    #model.thickness.range(18, 38) 
    120121    #model.thickness_pd.range(0, 1) 
    121122    #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, background=35, 
    134         radius=20, cap_radius=40, length=400, 
     133        scale=.08, radius=20, cap_radius=40, length=400, 
    135134        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=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, 
     139        theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0, 
     140        phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    142141        ) 
    143142 
    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) 
    154143    model.scale.range(0, 1) 
    155144 
    … …  
    157146elif name == "triaxial_ellipsoid": 
    158147    model = Model(kernel, 
    159         scale=0.08, background=35, 
    160         radius_equat_minor=15, radius_equat_major=20, radius_polar=500, 
     148        scale=0.08, req_minor=15, req_major=20, rpolar=500, 
    161149        sld=7.105, solvent_sld=.291, 
    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, 
     150        background=5, theta=0, phi=phi, psi=0, 
    166151        theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3, 
    167152        phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    168153        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, 
    169157        ) 
    170158 
    171159    # SET THE FITTING PARAMETERS 
    172     model.radius_equat_minor.range(15, 1000) 
    173     model.radius_equat_major.range(15, 1000) 
    174     #model.radius_polar.range(15, 1000) 
     160    model.req_minor.range(15, 1000) 
     161    model.req_major.range(15, 1000) 
     162    #model.rpolar.range(15, 1000) 
    175163    #model.background.range(0,1000) 
    176164    #model.theta_pd.range(0, 360) 
    … …  
    185173M = Experiment(data=data, model=model) 
    186174if section == "both": 
    187    tan_model = Model(model.sasmodel, **model.parameters()) 
     175   tan_model = Model(model.kernel, **model.parameters()) 
    188176   tan_model.phi = model.phi - 90 
    189177   tan_model.cutoff = cutoff 

  • sasmodels/compare.py

    rf72d70a ra769b54  
    324324        name = name.split('*')[0] 
    325325 
    326     # Suppress magnetism for python models (not yet implemented) 
    327     if callable(model_info.Iq): 
    328         pars.update(suppress_magnetism(pars)) 
    329  
    330326    if name == 'barbell': 
    331327        if pars['radius_bell'] < pars['radius']: 
    … …  
    346342        if pars['radius'] < pars['thick_string']: 
    347343            pars['radius'], pars['thick_string'] = pars['thick_string'], pars['radius'] 
     344        pars['num_pearls'] = math.ceil(pars['num_pearls']) 
    348345        pass 
    349346 
    … …  
    358355        for c in '1234': 
    359356            pars['Phi'+c] /= total 
     357 
     358    elif name == 'stacked_disks': 
     359        pars['n_stacking'] = math.ceil(pars['n_stacking']) 
    360360 
    361361def parlist(model_info, pars, is2d): 

  • sasmodels/compare_many.py

    rf72d70a r424fe00  
    106106    header = ('\n"Model","%s","Count","%d","Dimension","%s"' 
    107107              % (name, N, "2D" if is_2d else "1D")) 
    108     if not mono: 
    109         header += ',"Cutoff",%g'%(cutoff,) 
     108    if not mono: header += ',"Cutoff",%g'%(cutoff,) 
    110109    print(header) 
    111110 
    … …  
    162161    max_diff = [0] 
    163162    for k in range(N): 
    164         print("Model %s %d"%(name, k+1), file=sys.stderr) 
     163        print("%s %d"%(name, k), file=sys.stderr) 
    165164        seed = np.random.randint(1e6) 
    166165        pars_i = randomize_pars(model_info, pars, seed) 
    167166        constrain_pars(model_info, pars_i) 
    168         if 'sasview' in (base, comp): 
    169             constrain_new_to_old(model_info, pars_i) 
     167        constrain_new_to_old(model_info, pars_i) 
    170168        if mono: 
    171169            pars_i = suppress_pd(pars_i) 
    … …  
    189187    Print the command usage string. 
    190188    """ 
    191     print("usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono) (single|double|quad)", 
    192           file=sys.stderr) 
     189    print("usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono) (single|double|quad)") 
    193190 
    194191 
    … …  
    207204    print("""\ 
    208205 
    209 MODEL is the model name of the model or one of the model types listed in 
    210 sasmodels.core.list_models (all, py, c, double, single, opencl, 1d, 2d, 
    211 nonmagnetic, magnetic).  Model types can be combined, such as 2d+single. 
     206MODEL is the model name of the model or "all" for all the models 
     207in alphabetical order. 
    212208 
    213209COUNT is the number of randomly generated parameter sets to try. A value 
    … …  
    224220below the cutoff will be ignored.  Use "mono" for monodisperse models.  The 
    225221choice of polydisperse parameters, and the number of points in the distribution 
    226 is set in compare.py defaults for each model.  Polydispersity is given in the 
    227 "demo" attribute of each model. 
     222is set in compare.py defaults for each model. 
    228223 
    229224PRECISION is the floating point precision to use for comparisons.  If two 
    230 precisions are given, then compare one to the other.  Precision is one of 
    231 fast, single, double for GPU or single!, double!, quad! for DLL.  If no 
    232 precision is given, then use single and double! respectively. 
     225precisions are given, then compare one to the other, ignoring sasview. 
    233226 
    234227Available models: 
    … …  
    240233    Main program. 
    241234    """ 
    242     if len(argv) not in (3, 4, 5, 6): 
     235    if len(argv) not in (5, 6): 
    243236        print_help() 
    244237        return 
    245238 
    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) 
     239    model = argv[0] 
     240    if not (model in MODELS) and (model != "all"): 
     241        print('Bad model %s.  Use "all" or one of:'%model) 
    251242        print_models() 
    252         print('model types: all, py, c, double, single, opencl, 1d, 2d, nonmagnetic, magnetic') 
    253243        return 
    254244    try: 
    … …  
    257247        assert argv[2][1] == 'd' 
    258248        Nq = int(argv[2][2:]) 
    259         mono = len(argv) <= 3 or argv[3] == 'mono' 
     249        mono = argv[3] == 'mono' 
    260250        cutoff = float(argv[3]) if not mono else 0 
    261         base = argv[4] if len(argv) > 4 else "single" 
    262         comp = argv[5] if len(argv) > 5 else "double!" 
     251        base = argv[4] 
     252        comp = argv[5] if len(argv) > 5 else "sasview" 
    263253    except Exception: 
    264254        traceback.print_exc() 
    … …  
    268258    data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 
    269259                             'accuracy': 'Low', 'view':'log', 'zero': False}) 
     260    model_list = [model] if model != "all" else MODELS 
    270261    for model in model_list: 
    271262        compare_instance(model, data, index, N=count, mono=mono, 

  • sasmodels/conversion_table.py

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

  • sasmodels/core.py

    r5124c969 r52e9a45  
    6969        * magnetic: models with an sld 
    7070        * nommagnetic: models without an sld 
    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('+')): 
     71    """ 
     72    if kind and kind not in KINDS: 
    7773        raise ValueError("kind not in " + ", ".join(KINDS)) 
    7874    files = sorted(glob(joinpath(generate.MODEL_PATH, "[a-zA-Z]*.py"))) 
    7975    available_models = [basename(f)[:-3] for f in files] 
    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)] 
     76    selected = [name for name in available_models if _matches(name, kind)] 
    8677 
    8778    return selected 

  • sasmodels/model_test.py

    r598857b r1a6cd57  
    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  
    11499        if is_py:  # kernel implemented in python 
    115100            test_name = "Model: %s, Kernel: python"%model_name 
    … …  
    118103                                 test_method_name, 
    119104                                 platform="dll",  # so that 
    120                                  dtype="double", 
    121                                  stash=stash) 
     105                                 dtype="double") 
    122106            suite.addTest(test) 
    123107        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  
    136108            # test using opencl if desired and available 
    137109            if 'opencl' in loaders and core.HAVE_OPENCL: 
    … …  
    144116                test = ModelTestCase(test_name, model_info, 
    145117                                     test_method_name, 
    146                                      platform="ocl", dtype=None, 
    147                                      stash=stash) 
     118                                     platform="ocl", dtype=None) 
    148119                #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") 
    149130                suite.addTest(test) 
    150131 
    … …  
    163144        """ 
    164145        def __init__(self, test_name, model_info, test_method_name, 
    165                      platform, dtype, stash): 
    166             # type: (str, ModelInfo, str, str, DType, List[Any]) -> None 
     146                     platform, dtype): 
     147            # type: (str, ModelInfo, str, str, DType) -> None 
    167148            self.test_name = test_name 
    168149            self.info = model_info 
    169150            self.platform = platform 
    170151            self.dtype = dtype 
    171             self.stash = stash  # container for the results of the first run 
    172152 
    173153            setattr(self, test_method_name, self.run_all) 
    … …  
    187167                #({}, (0.0, 0.0), None), 
    188168                # test vector form 
    189                 ({}, [0.001, 0.01, 0.1], [None]*3), 
     169                ({}, [0.1]*2, [None]*2), 
    190170                ({}, [(0.1, 0.1)]*2, [None]*2), 
    191171                # test that ER/VR will run if they exist 
    … …  
    194174                ] 
    195175 
    196             tests = smoke_tests + self.info.tests 
     176            tests = self.info.tests 
    197177            try: 
    198178                model = build_model(self.info, dtype=self.dtype, 
    199179                                    platform=self.platform) 
    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)) 
     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 
    217190 
    218191            except: 
    219192                annotate_exception(self.test_name) 
    220193                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 
    258194 
    259195        def run_one(self, model, test): 
    … …  
    271207 
    272208            if x[0] == 'ER': 
    273                 actual = np.array([call_ER(model.info, pars)]) 
     209                actual = [call_ER(model.info, pars)] 
    274210            elif x[0] == 'VR': 
    275                 actual = np.array([call_VR(model.info, pars)]) 
     211                actual = [call_VR(model.info, pars)] 
    276212            elif isinstance(x[0], tuple): 
    277213                qx, qy = zip(*x) 
    … …  
    302238                                    'f(%s); expected:%s; actual:%s' 
    303239                                    % (xi, yi, actual_yi)) 
    304             return actual 
    305240 
    306241    return ModelTestCase 

  • sasmodels/modelinfo.py

    rf88e248 r85fe7f8  
    230230    defined as a sublist with the following elements: 
    231231 
    232     *name* is the name that will be displayed to the user.  Names 
     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 
    233234    should be lower case, with words separated by underscore.  If 
    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. 
     235    acronyms are used, the whole acronym should be upper case. 
    239236 
    240237    *units* should be one of *degrees* for angles, *Ang* for lengths, 
    … …  
    606603        # Using the call_parameters table, we already have expanded forms 
    607604        # for each of the vector parameters; put them in a lookup table 
    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) 
     605        expanded_pars = dict((p.name, p) for p in self.call_parameters) 
    610606 
    611607        def append_group(name): 
    … …  
    734730    info.docs = kernel_module.__doc__ 
    735731    info.category = getattr(kernel_module, 'category', None) 
     732    info.single = getattr(kernel_module, 'single', True) 
     733    info.opencl = getattr(kernel_module, 'opencl', True) 
    736734    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
    737735    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    … …  
    747745    info.profile = getattr(kernel_module, 'profile', None) # type: ignore 
    748746    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)) 
    752747 
    753748    # multiplicity info 

  • sasmodels/models/core_multi_shell.c

    r925ad6e r925ad6e  
    88} 
    99 
    10 static double 
    11 form_volume(double core_radius, double fp_n, double thickness[]) 
     10double 
     11form_volume(double core_radius, double n, double thickness[]); 
     12double 
     13form_volume(double core_radius, double n, double thickness[]) 
    1214{ 
    1315  double r = core_radius; 
    14   int n = (int)(fp_n+0.5); 
    1516  for (int i=0; i < n; i++) { 
    1617    r += thickness[i]; 
    … …  
    1920} 
    2021 
    21 static double 
     22double 
    2223Iq(double q, double core_sld, double core_radius, 
    23    double solvent_sld, double fp_n, double sld[], double thickness[]) 
     24   double solvent_sld, double num_shells, double sld[], double thickness[]); 
     25double 
     26Iq(double q, double core_sld, double core_radius, 
     27   double solvent_sld, double num_shells, double sld[], double thickness[]) 
    2428{ 
    25   const int n = (int)(fp_n+0.5); 
     29  const int n = (int)ceil(num_shells); 
    2630  double f, r, last_sld; 
    2731  r = core_radius; 

  • sasmodels/models/core_multi_shell.py

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

  • sasmodels/models/core_shell_parallelepiped.py

    rcb0dc22 r5810f00  
    44 
    55Calculates the form factor for a rectangular solid with a core-shell structure. 
    6 The thickness and the scattering length density of the shell or  
    7 "rim" can be different on each (pair) of faces. However at this time 
    8 the model does **NOT** actually calculate a c face rim despite the presence of 
    9 the 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. 
     6**The thickness and the scattering length density of the shell or "rim" 
     7can be different on all three (pairs) of faces.** 
    148 
    159The form factor is normalized by the particle volume $V$ such that 
    … …  
    4135    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    4236 
    43 **meaning that there are "gaps" at the corners of the solid.**  Again note that 
    44 $t_C = 0$ currently.  
     37**meaning that there are "gaps" at the corners of the solid.** 
    4538 
    4639The intensity calculated follows the :ref:`parallelepiped` model, with the 
    … …  
    4841amplitudes of the core and shell, in the same manner as a core-shell model. 
    4942 
     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)} 
    5051 
    5152.. note:: 
    … …  
    9293    detector plane. 
    9394 
     95Validation 
     96---------- 
     97 
     98The model uses the form factor calculations implemented in a c-library provided 
     99by the NIST Center for Neutron Research (Kline, 2006). 
     100 
    94101References 
    95102---------- 
    … …  
    106113 
    107114* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    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 
     115* **Last Modified by:** Paul Butler **Date:** September 30, 2016 
     116* **Last Reviewed by:** Miguel Gonzales **Date:** March 21, 2016 
    111117""" 
    112118 

  • sasmodels/models/ellipsoid.c

    r130d4c7 r925ad6e  
    44    double radius_polar, double radius_equatorial, double theta, double phi); 
    55 
    6 static double 
    7 _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha) 
     6double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha); 
     7double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha) 
    88{ 
    99    double ratio = radius_polar/radius_equatorial; 
    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     // 
     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) 
    1714    // Instead of using pythagoras we could pass in sin and cos; this may be 
    1815    // slightly better for 2D which has already computed it, but it introduces 
    … …  
    2017    // leave it as is. 
    2118    const double r = radius_equatorial 
    22                      * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0)); 
     19                     * sqrt(1.0 + sin_alpha*sin_alpha*(ratio*ratio - 1.0)); 
    2320    const double f = sas_3j1x_x(q*r); 
    2421 
    … …  
    4239    double total = 0.0; 
    4340    for (int i=0;i<76;i++) { 
    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); 
     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); 
    4744    } 
    4845    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    6259    double q, sin_alpha, cos_alpha; 
    6360    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    64     const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
     61    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 
    6562    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    6663 

  • sasmodels/models/flexible_cylinder.c

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

  • sasmodels/models/fractal.c

    r4788822 r925ad6e  
    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  } 
    82 
    93static double 
    … …  
    1913 
    2014    //calculate P(q) for the spherical subunits 
    21     const double pq = square(form_volume(radius) * (sld_block-sld_solvent) 
    22                       *sas_3j1x_x(q*radius)); 
     15    const double V = M_4PI_3*cube(radius); 
     16    const double pq = V * square((sld_block-sld_solvent)*sas_3j1x_x(q*radius)); 
    2317 
    2418    // scale to units cm-1 sr-1 (assuming data on absolute scale) 

  • sasmodels/models/fractal.py

    rdf89d77 r925ad6e  
    2020.. math:: 
    2121 
    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 \\ 
     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 
    2528    S(q) &= 1 + \frac{D_f\  \Gamma\!(D_f-1)}{[1+1/(q \xi)^2\  ]^{(D_f -1)/2}} 
    2629    \frac{\sin[(D_f-1) \tan^{-1}(q \xi) ]}{(q R_0)^{D_f}} 
    … …  
    2932is the fractal dimension, representing the self similarity of the structure.  
    3033Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function  
    31 diverges at $D_f=0$ and $D_f=1$. 
     34diverges at $D_f$=0 and $D_f$=1.   
    3235 
    3336**Polydispersity on the radius is provided for.** 
    … …  
    4447---------- 
    4548 
    46 .. [#] J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 
     49J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 
    4750 
    48 Authorship and Verification 
    49 ---------------------------- 
     51**Author:** NIST IGOR/DANSE **on:** pre 2010 
    5052 
    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 
     53**Last Modified by:** Paul Butler **on:** March 20, 2016 
     54 
     55**Last Reviewed by:** Paul Butler **on:** March 20, 2016 
    5556 
    5657""" 
    … …  
    8384parameters = [["volfraction", "", 0.05, [0.0, 1], "", 
    8485               "volume fraction of blocks"], 
    85               ["radius",    "Ang",  5.0, [0.0, inf], "volume", 
     86              ["radius",    "Ang",  5.0, [0.0, inf], "", 
    8687               "radius of particles"], 
    8788              ["fractal_dim",      "",  2.0, [0.0, 6.0], "", 

  • sasmodels/models/hayter_msa.c

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

  • sasmodels/models/lamellar_hg_stack_caille.c

    r1c6286d ra807206  
    33*/ 
    44 
    5 static double 
    6 Iq(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) 
     5double 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 
     15double 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) 
    1524{ 
    16   int Nlayers = (int)(fp_Nlayers+0.5);    //cast to an integer for the loop 
     25  double NN;   //local variables of coefficient wave 
    1726  double inten,Pq,Sq,alpha,temp,t2; 
    1827  //double dQ, dQDefault, t1, t3; 
     28  int ii,NNint; 
    1929  // from wikipedia 0.577215664901532860606512090082402431042159335 
    2030  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015 
    … …  
    2232  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015 
    2333 
    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); 
     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); 
    2638  Pq *= Pq; 
    2739  Pq *= 4.0/(qval*qval); 
    2840 
     41  NNint = (int)NN;    //cast to an integer for the loop 
     42  ii=0; 
    2943  Sq = 0.0; 
    30   for(int ii=1; ii < Nlayers; ii++) { 
     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 
    3148    temp = 0.0; 
    3249    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler); 
    … …  
    3552    //t3 = dQ*dQ*dd*dd*ii*ii; 
    3653 
    37     temp = 1.0-(double)ii/(double)Nlayers; 
     54    temp = 1.0-ii/NN; 
    3855    //temp *= cos(dd*qval*ii/(1.0+t1)); 
    3956    temp *= cos(dd*qval*ii); 
    4057    //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii); 
    4158    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 
    42     temp *= exp(-t2/2.0); 
     59    temp *= exp(-t2/2.0 ); 
    4360    //temp /= sqrt(1.0+t1); 
    4461 
    … …  
    5471 
    5572  inten *= 1.0e-04;   // 1/A to 1/cm 
    56   return inten; 
     73  return(inten); 
    5774} 
    5875 

  • sasmodels/models/lamellar_hg_stack_caille.py

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

  • sasmodels/models/lamellar_stack_caille.c

    r1c6286d r0bef47b  
    33*/ 
    44 
    5 static double 
    6 Iq(double qval, 
    7    double del, 
    8    double fp_Nlayers, 
    9    double dd, 
    10    double Cp, 
    11    double sld, 
    12    double solvent_sld) 
     5double Iq(double qval, 
     6      double del, 
     7      double Nlayers,  
     8      double dd, 
     9          double Cp,  
     10      double sld, 
     11      double solvent_sld); 
     12 
     13double Iq(double qval, 
     14      double del, 
     15      double Nlayers,  
     16      double dd, 
     17          double Cp,  
     18      double sld, 
     19      double solvent_sld) 
    1320{ 
    14   int Nlayers = (int)(fp_Nlayers+0.5);    //cast to an integer for the loop 
    15   double contr;   //local variables of coefficient wave 
     21  double contr,NN;   //local variables of coefficient wave 
    1622  double inten,Pq,Sq,alpha,temp,t2; 
    1723  //double dQ, dQDefault, t1, t3; 
     24  int ii,NNint; 
    1825  // from wikipedia 0.577215664901532860606512090082402431042159335 
    1926  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015 
    … …  
    2128  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015 
    2229 
     30  NN = trunc(Nlayers);    //be sure that NN is an integer 
     31   
    2332  contr = sld - solvent_sld; 
    2433 
    2534  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)); 
    2635 
     36  NNint = (int)NN;    //cast to an integer for the loop 
     37  ii=0; 
    2738  Sq = 0.0; 
    28   for (int ii=1; ii < Nlayers; ii++) { 
     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 
    2944    temp = 0.0; 
    3045    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler); 
    … …  
    3348    //t3 = dQ*dQ*dd*dd*ii*ii; 
    3449 
    35     temp = 1.0 - (double)ii / (double)Nlayers; 
     50    temp = 1.0-ii/NN; 
    3651    //temp *= cos(dd*qval*ii/(1.0+t1)); 
    3752    temp *= cos(dd*qval*ii); 
    3853    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 
    39     temp *= exp(-t2/2.0); 
     54    temp *= exp(-t2/2.0 ); 
    4055    //temp /= sqrt(1.0+t1); 
    4156 

  • sasmodels/models/lamellar_stack_caille.py

    ra57b31d r7c57861  
    8888parameters = [ 
    8989    ["thickness",        "Ang",      30.0,  [0, inf],   "volume", "sheet thickness"], 
    90     ["Nlayers",          "",          20,   [1, inf],   "",       "Number of layers"], 
     90    ["Nlayers",          "",          20,   [0, 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

    r5467cd8 r4962519  
    22 
    33*/ 
    4 double paraCryst_sn(double ww, double qval, double davg, int Nlayers, double an); 
    5 double paraCryst_an(double ww, double qval, double davg, int Nlayers); 
     4double Iq(double qval, 
     5      double th, 
     6      double Nlayers,  
     7          double davg,  
     8          double pd, 
     9      double sld, 
     10      double solvent_sld); 
     11double paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an); 
     12double paraCryst_an(double ww, double qval, double davg, long Nlayers); 
    613 
    7 static double 
    8 Iq(double qval, 
    9    double th, 
    10    double fp_Nlayers, 
    11    double davg, 
    12    double pd, 
    13    double sld, 
    14    double solvent_sld) 
     14double Iq(double qval, 
     15      double th, 
     16      double Nlayers,  
     17          double davg,  
     18          double pd, 
     19      double sld, 
     20      double solvent_sld) 
    1521{ 
     22     
     23        double inten,contr,xn; 
     24        double xi,ww,Pbil,Znq,Snq,an; 
     25        long n1,n2; 
     26         
     27        contr = sld - solvent_sld; 
    1628        //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 
    2029         
    21         const double ww = exp(-0.5*square(qval*pd*davg)); 
     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); 
    2235 
    2336        //calculate the n1 contribution 
    24         double Znq,Snq,an; 
    2537        an = paraCryst_an(ww,qval,davg,n1); 
    2638        Snq = paraCryst_sn(ww,qval,davg,n1,an); 
    27  
     39         
    2840        Znq = xn*Snq; 
    2941         
    … …  
    4052//      Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg)) 
    4153         
    42         const double xi = th/2.0;               //use 1/2 the bilayer thickness 
    43         const double Pbil = square(sas_sinx_x(qval*xi)); 
     54        xi = th/2.0;            //use 1/2 the bilayer thickness 
     55        Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi)); 
    4456         
    45         const double contr = sld - solvent_sld; 
    46         const double inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval); 
     57        inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval); 
     58        inten *= 1.0e-04; 
    4759//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); 
    48         return 1.0e-4*inten; 
     60        return(inten); 
    4961} 
    5062 
    5163// functions for the lamellar paracrystal model 
    5264double 
    53 paraCryst_sn(double ww, double qval, double davg, int Nlayers, double an) { 
     65paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an) { 
    5466         
    5567        double Snq; 
    … …  
    5769        Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) ); 
    5870         
    59         return Snq; 
     71        return(Snq); 
    6072} 
    6173 
    6274double 
    63 paraCryst_an(double ww, double qval, double davg, int Nlayers) { 
     75paraCryst_an(double ww, double qval, double davg, long Nlayers) { 
     76         
    6477        double an; 
    6578         
    … …  
    6982        an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg); 
    7083         
    71         return an; 
     84        return(an); 
    7285} 
    7386 

  • sasmodels/models/lamellar_stack_paracrystal.py

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

  • sasmodels/models/linear_pearls.c

    r925ad6e r925ad6e  
    44            double radius, 
    55            double edge_sep, 
    6             double fp_num_pearls, 
     6            double num_pearls, 
    77            double pearl_sld, 
    88            double solvent_sld); 
    … …  
    1111            double radius, 
    1212            double edge_sep, 
    13             int num_pearls, 
     13            double num_pearls, 
    1414            double pearl_sld, 
    1515            double solvent_sld); 
    1616 
    1717 
    18 double form_volume(double radius, double fp_num_pearls) 
     18double form_volume(double radius, double num_pearls) 
    1919{ 
    20     int num_pearls = (int)(fp_num_pearls + 0.5); 
    2120    // Pearl volume 
    2221    double pearl_vol = M_4PI_3 * cube(radius); 
    … …  
    2827            double radius, 
    2928            double edge_sep, 
    30             int num_pearls, 
     29            double num_pearls, 
    3130            double pearl_sld, 
    3231            double solvent_sld) 
    3332{ 
     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 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); 
     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)); 
    5253    } 
    5354    // form factor for num_pearls 
    54     double form_factor = 1.0e-4 * structure_factor * square(m_s*psi) / tot_vol; 
     55    double form_factor = 1.0e-4 * n_contrib * square(m_s*psi) / tot_vol; 
    5556 
    5657    return form_factor; 
    … …  
    6061            double radius, 
    6162            double edge_sep, 
    62             double fp_num_pearls, 
     63            double num_pearls, 
    6364            double pearl_sld, 
    6465            double solvent_sld) 
    6566{ 
    6667 
    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 r925ad6e  
    1616.. math:: 
    1717 
    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] 
     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] 
    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, [1, inf],     "", "Number of the pearls"], 
     58    ["num_pearls",  "",           3.0, [0, 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

    rec1d4bc r925ad6e  
    1 static double 
    2 form_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  
    12 static double 
    13 multilayer_vesicle_kernel(double q, 
     1static 
     2double multilayer_vesicle_kernel(double q, 
    143          double volfraction, 
    154          double radius, 
    … …  
    187          double sld_solvent, 
    198          double sld, 
    20           int n_shells) 
     9          double n_pairs) 
    2110{ 
    2211    //calculate with a loop, two shells at a time 
    … …  
    4029 
    4130        //do 2 layers at a time 
    42         ii++; 
     31        ii += 1; 
    4332 
    44     } while(ii <= n_shells-1);  //change to make 0 < n_shells < 2 correspond to 
     33    } while(ii <= n_pairs-1);  //change to make 0 < n_pairs < 2 correspond to 
    4534                               //unilamellar vesicles (C. Glinka, 11/24/03) 
    4635 
    47     return 1.0e-4*volfraction*fval*fval;  // Volume normalization happens in caller 
     36    fval *= volfraction*1.0e-4*fval/voli; 
     37 
     38    return(fval); 
    4839} 
    4940 
    50 static double 
    51 Iq(double q, 
     41static 
     42double Iq(double q, 
    5243          double volfraction, 
    5344          double radius, 
    … …  
    5647          double sld_solvent, 
    5748          double sld, 
    58           double fp_n_shells) 
     49          double n_pairs) 
    5950{ 
    60     int n_shells = (int)(fp_n_shells + 0.5); 
    6151    return multilayer_vesicle_kernel(q, 
    6252           volfraction, 
    … …  
    6656           sld_solvent, 
    6757           sld, 
    68            n_shells); 
     58           n_pairs); 
    6959} 
    7060 

  • sasmodels/models/multilayer_vesicle.py

    r5d23de2 r925ad6e  
    33---------- 
    44 
    5 This model is a trivial extension of the core_shell_sphere function where the 
    6 core 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 
    8 returns the same as the vesicle model, except for the normalisation, which here 
    9 is to outermost volume. The shell thicknesses and SLD are constant for all 
    10 shells as expected for a multilayer vesicle. 
     5This 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 
     7with layers of solvent. For *N = 1*, this returns the same as the vesicle model, 
     8except for the normalisation, which here is to outermost volume. 
     9The shell thicknessess and SLD are constant for all shells as expected for 
     10a multilayer vesicle. 
    1111 
    1212.. figure:: img/multi_shell_geometry.jpg 
    … …  
    1919 
    2020.. math:: 
    21     P(q) = \text{scale} \cdot \frac{\phi}{V(R_N)} F^2(q) + \text{background} 
     21 
     22    P(q) = \frac{\text{scale.volfraction}}{V_t} F^2(q) + \text{background} 
    2223 
    2324where 
    2425 
    2526.. math:: 
    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} 
     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} 
    2931     \right] 
    3032 
    31 for 
    3233 
    33 .. math:: 
     34where $R_i = r_c + (i-1)(t_s + t_w)$ 
     35    
     36where $V_t$ is the volume of the whole particle, $V(R)$ is the volume of a sphere 
     37of radius $R$, $r_c$ is the radius of the core, $\rho_{shell}$ is the scattering length  
     38density of a shell, $\rho_{solv}$ is the scattering length density of the solvent. 
    3439 
    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 
    39 of radius $r$, $r_c$ is the radius of the core, $t_s$ is the thickness of 
    40 the 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  
    44 USAGE 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. 
    8340 
    8441The 2D scattering intensity is the same as 1D, regardless of the orientation 
    … …  
    8946    q = \sqrt{q_x^2 + q_y^2} 
    9047 
     48 
     49The outer most radius 
     50 
     51$radius + n\_pairs * thick\_shell + (n\_pairs- 1) * thick\_solvent$ 
     52 
     53is used for both the volume fraction normalization and for the  
     54effective radius for *S(Q)* when $P(Q) * S(Q)$ is applied. 
     55 
    9156For information about polarised and magnetic scattering, see 
    9257the :ref:`magnetism` documentation. 
     58 
     59This code is based on the form factor calculations implemented in the NIST 
     60Center for Neutron Research provided c-library (Kline, 2006). 
    9361 
    9462References 
    9563---------- 
    9664 
    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). 
     65B Cabane, *Small Angle Scattering Methods*, 
     66in *Surfactant Solutions: New Methods of Investigation*, 
     67Ch.2, Surfactant Science Series Vol. 22, Ed. R Zana and M Dekker, 
     68New York, (1987). 
    10069 
    101 Authorship and Verification 
    102 ---------------------------- 
     70**Author:** NIST IGOR/DANSE **on:** pre 2010 
    10371 
    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 
     72**Last Modified by:** Piotr Rozyczko**on:** Feb 24, 2016 
     73 
     74**Last Reviewed by:** Paul Butler **on:** March 20, 2016 
    10875 
    10976""" 
    … …  
    12289    sld_solvent: solvent scattering length density 
    12390    sld: shell scattering length density 
    124     n_shells:number of "shell plus solvent" layer pairs 
     91    n_pairs:number of "shell plus solvent" layer pairs 
    12592    background: incoherent background 
    12693        """ 
    … …  
    13198parameters = [ 
    13299    ["volfraction", "",  0.05, [0.0, 1],  "", "volume fraction of vesicles"], 
    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"], 
     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"], 
    136103    ["sld_solvent",    "1e-6/Ang^2",  6.4, [-inf, inf], "sld", "solvent scattering length density"], 
    137104    ["sld",   "1e-6/Ang^2",  0.4, [-inf, inf], "sld", "Shell scattering length density"], 
    138     ["n_shells",     "",            2.0, [1.0, inf],  "volume", "Number of shell plus solvent layer pairs"], 
     105    ["n_pairs",     "",            2.0, [1.0, inf],  "", "Number of shell plus solvent layer pairs"], 
    139106    ] 
    140107# pylint: enable=bad-whitespace, line-too-long 
    141108 
    142 # TODO: proposed syntax for specifying which parameters can be polydisperse 
    143 #polydispersity = ["radius", "thick_shell"] 
    144  
    145109source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 
    146110 
    147 def 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 
     111polydispersity = ["radius", "n_pairs"] 
    150112 
    151113demo = dict(scale=1, background=0, 
    … …  
    156118            sld_solvent=6.4, 
    157119            sld=0.4, 
    158             n_shells=2.0) 
     120            n_pairs=2.0) 
    159121 
    160122tests = [ 
    … …  
    165127      'sld_solvent': 6.4, 
    166128      'sld': 0.4, 
    167       'n_shells': 2.0, 
     129      'n_pairs': 2.0, 
    168130      'scale': 1.0, 
    169131      'background': 0.001, 
    … …  
    176138      'sld_solvent': 6.4, 
    177139      'sld': 0.4, 
    178       'n_shells': 2.0, 
     140      'n_pairs': 2.0, 
    179141      'scale': 1.0, 
    180142      'background': 0.001, 

  • sasmodels/models/onion.c

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

  • sasmodels/models/onion.py

    r925ad6e r925ad6e  
    323323    Returns shape profile with x=radius, y=SLD. 
    324324    """ 
    325     n_shells = int(n_shells+0.5) 
     325 
    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_shells, thickness): 
     368def ER(radius_core, n, thickness): 
    369369    """Effective radius""" 
    370     n = int(n_shells[0]+0.5) 
    371     return np.sum(thickness[:n], axis=0) + radius_core 
     370    return np.sum(thickness[:int(n[0])], axis=0) + radius_core 
    372371 
    373372demo = { 

  • sasmodels/models/pearl_necklace.c

    r3f853beb r4b541ac  
    11double form_volume(double radius, double edge_sep, 
    2     double thick_string, double fp_num_pearls); 
     2    double thick_string, double num_pearls); 
    33double Iq(double q, double radius, double edge_sep, 
    4     double thick_string, double fp_num_pearls, double sld, 
     4    double thick_string, double 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     int num_pearls, double sld_pearl, double sld_string, double sld_solv) 
     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) 
    1313{ 
    1414    // number of string segments 
    15     const int num_strings = num_pearls - 1; 
     15    num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
     16    const double num_strings = num_pearls - 1.0; 
    1617 
    1718    //each masses: contrast * volume 
    … …  
    6869 
    6970double form_volume(double radius, double edge_sep, 
    70     double thick_string, double fp_num_pearls) 
     71    double thick_string, double num_pearls) 
    7172{ 
    72     const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls 
    73     const int num_strings = num_pearls - 1; 
     73    num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
     74 
     75    const double num_strings = num_pearls - 1.0; 
    7476    const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; 
    7577    const double pearl_vol = M_4PI_3 * radius * radius * radius; 
    … …  
    8082 
    8183double Iq(double q, double radius, double edge_sep, 
    82     double thick_string, double fp_num_pearls, double sld, 
     84    double thick_string, double num_pearls, double sld,  
    8385    double string_sld, double solvent_sld) 
    8486{ 
    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, 
     87    const double form = _pearl_necklace_kernel(q, radius, edge_sep, 
    8788        thick_string, num_pearls, sld, string_sld, solvent_sld); 
    8889 

  • sasmodels/models/pearl_necklace.py

    r1bd1ea2 r4b541ac  
    8282              ["thick_string", "Ang", 2.5, [0, inf], "volume", 
    8383               "Thickness of the chain linkage"], 
    84               ["num_pearls", "none", 3, [1, inf], "volume", 
     84              ["num_pearls", "none", 3, [0, 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) 
    103102    number_of_strings = num_pearls - 1.0 
    104103    string_vol = edge_sep * pi * pow((thick_string / 2.0), 2.0) 
    … …  
    112111    Calculation for effective radius. 
    113112    """ 
    114     num_pearls = int(num_pearls + 0.5) 
    115113    tot_vol = volume(radius, edge_sep, thick_string, num_pearls) 
    116114    rad_out = (tot_vol/(4.0/3.0*pi)) ** (1./3.) 

  • sasmodels/models/raspberry.py

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

  • sasmodels/models/rpa.c

    r20c856a r6351bfa  
    1 double Iq(double q, double fp_case_num, 
     1double Iq(double q, double 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 fp_case_num, 
     7double Iq(double q, double case_num, 
    88    double N[],    // DEGREE OF POLYMERIZATION 
    99    double Phi[],  // VOL FRACTION 
    … …  
    1515    ) 
    1616{ 
    17   int icase = (int)(fp_case_num+0.5); 
     17  int icase = (int)case_num; 
    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    
    42   // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc) 
     41 
    4342  //icase was shifted to N-1 from the original code 
    4443  if (icase <= 1){ 
    … …  
    5857    Kab = Kac = Kad = -0.0004; 
    5958  } 
    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) 
     59 
    6560  Nab=sqrt(N[0]*N[1]); 
    6661  Nac=sqrt(N[0]*N[2]); 
    … …  
    8479  Phicd=sqrt(Phi[2]*Phi[3]); 
    8580 
    86   // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk  
    8781  Xa=q*q*b[0]*b[0]*N[0]/6.0; 
    8882  Xb=q*q*b[1]*b[1]*N[1]/6.0; 
    … …  
    9084  Xd=q*q*b[3]*b[3]*N[3]/6.0; 
    9185 
    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 
     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); 
    9689  S0ab=(Phiab*vab*Nab)*Pab; 
    97   Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor 
     90  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 
    9891  S0ac=(Phiac*vac*Nac)*Pac; 
    99   Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block 
     92  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 
    10093  S0ad=(Phiad*vad*Nad)*Pad; 
    10194 
    10295  S0ba=S0ab; 
    103   Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain 
     96  Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); 
    10497  S0bb=N[1]*Phi[1]*v[1]*Pbb; 
    105   Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock 
     98  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 
    10699  S0bc=(Phibc*vbc*Nbc)*Pbc; 
    107   Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock 
     100  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 
    108101  S0bd=(Phibd*vbd*Nbd)*Pbd; 
    109102 
    110103  S0ca=S0ac; 
    111104  S0cb=S0bc; 
    112   Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain 
     105  Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); 
    113106  S0cc=N[2]*Phi[2]*v[2]*Pcc; 
    114   Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock 
     107  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 
    115108  S0cd=(Phicd*vcd*Ncd)*Pcd; 
    116109 
    … …  
    118111  S0db=S0bd; 
    119112  S0dc=S0cd; 
    120   Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain 
     113  Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); 
    121114  S0dd=N[3]*Phi[3]*v[3]*Pdd; 
    122  
    123   // Reset all unused partial structure factors to 0 (depends on case) 
    124115  //icase was shifted to N-1 from the original code 
    125116  switch(icase){ 
    … …  
    202193  S0dc=S0cd; 
    203194 
    204   // self chi parameter is 0 ... of course 
    205195  Kaa=0.0; 
    206196  Kbb=0.0; 
    … …  
    253243  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 
    254244 
    255   // D is considered the matrix or background component so enters here 
    256245  m=1.0/(S0dd-ZZ); 
    257246 
    … …  
    308297  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 
    309298 
    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.  
    314299  Nav=6.022045e+23; 
    315300  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); 
    … …  
    318303 
    319304  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;     
    323305 
    324306  return Intg; 

  • sasmodels/models/rpa.py

    r4f9e288 r40a87fa  
    11r""" 
    2 Definition 
    3 ---------- 
    4  
    52Calculates the macroscopic scattering intensity for a multi-component 
    63homogeneous mixture of polymers using the Random Phase Approximation. 
    … …  
    2724Case 9: A-B-C-D tetra-block copolymer 
    2825 
    29 .. note:: 
    30     These case numbers are different from those in the NIST SANS package! 
     26**NB: these case numbers are different from those in the NIST SANS package!** 
    3127 
    32 The models are based on the papers by Akcasu et al. [#Akcasu]_ and by 
    33 Hammouda [#Hammouda]_ assuming the polymer follows Gaussian statistics such 
    34 that $R_g^2 = n b^2/6$ where $b$ is the statistical segment length and $n$ is 
    35 the number of statistical segment lengths. A nice tutorial on how these are 
    36 constructed and implemented can be found in chapters 28 and 39 of Boualem 
    37 Hammouda's 'SANS Toolbox'[#toolbox]_. 
     28Only one case can be used at any one time. 
    3829 
    39 In brief the macroscopic cross sections are derived from the general forms 
    40 for homopolymer scattering and the multiblock cross-terms while the inter 
    41 polymer cross terms are described in the usual way by the $\chi$ parameter. 
     30The RPA (mean field) formalism only applies only when the multicomponent 
     31polymer mixture is in the homogeneous mixed-phase region. 
    4232 
    43 USAGE NOTES: 
     33**Component D is assumed to be the "background" component (ie, all contrasts 
     34are calculated with respect to component D).** So the scattering contrast 
     35for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`. 
    4436 
    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. 
     37Depending on which case is being used, the number of fitting parameters - the 
     38segment lengths (ba, bb, etc) and $\chi$ parameters (Kab, Kac, etc) - vary. 
     39The *scale* parameter should be held equal to unity. 
    5340 
    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). 
     41The input parameters are the degrees of polymerization, the volume fractions, 
     42the specific volumes, and the neutron scattering length densities for each 
     43component. 
    6144 
    6245 
    … …  
    6447---------- 
    6548 
    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  
    71 Authorship 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 
     49A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 
    7850""" 
    7951 
    … …  
    8153 
    8254name = "rpa" 
    83 title = "Random Phase Approximation" 
     55title = "Random Phase Approximation - unfinished work in progress" 
    8456description = """ 
    8557This formalism applies to multicomponent polymer mixtures in the 
    … …  
    11890    ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    11991    ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 
    120     ["v[4]", "mL/mol", 100.0, [0, inf], "", "molar volume"], 
     92    ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    12193    ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    12294    ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 
    … …  
    135107 
    136108control = "case_num" 
    137 HIDE_ALL = set("Phi4".split()) 
    138 HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()).union(HIDE_ALL) 
     109HIDE_NONE = set() 
     110HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()) 
    139111HIDE_AB = set("N2 Phi2 v2 L2 b2 K23 K24".split()).union(HIDE_A) 
    140112def hidden(case_num): 
    … …  
    142114    Return a list of parameters to hide depending on the multiplicity parameter. 
    143115    """ 
    144     case_num = int(case_num+0.5) 
    145116    if case_num < 2: 
    146117        return HIDE_AB 
    … …  
    148119        return HIDE_A 
    149120    else: 
    150         return HIDE_ALL 
     121        return HIDE_NONE 
    151122 

  • sasmodels/models/spherical_sld.c

    r925ad6e r925ad6e  
    11static double form_volume( 
    2     double fp_n_shells, 
     2    int n_shells, 
    33    double thickness[], 
    44    double interface[]) 
    55{ 
    6     int n_shells= (int)(fp_n_shells + 0.5); 
    76    double r = 0.0; 
    87    for (int i=0; i < n_shells; i++) { 
    … …  
    2120        return pow(z, nu); 
    2221    } else if (shape==2) { 
    23         return 1.0 - pow(1.0 - z, nu); 
     22        return 1.0 - pow(1. - z, nu); 
    2423    } else if (shape==3) { 
    2524        return expm1(-nu*z)/expm1(-nu); 
    … …  
    4544static double Iq( 
    4645    double q, 
    47     double fp_n_shells, 
     46    int n_shells, 
    4847    double sld_solvent, 
    4948    double sld[], 
    … …  
    5251    double shape[], 
    5352    double nu[], 
    54     double fp_n_steps) 
     53    int n_steps) 
    5554{ 
    5655    // 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); 
    5956    double f=0.0; 
    6057    double r=0.0; 

  • sasmodels/models/spherical_sld.py

    r925ad6e r925ad6e  
    233233    """ 
    234234 
    235     n_shells = int(n_shells + 0.5) 
    236     n_steps = int(n_steps + 0.5) 
    237235    z = [] 
    238236    rho = [] 
    … …  
    242240    rho.append(sld[0]) 
    243241 
    244     for i in range(0, n_shells): 
     242    for i in range(0, int(n_shells)): 
    245243        z_next += thickness[i] 
    246244        z.append(z_next) 
    … …  
    263261def ER(n_shells, thickness, interface): 
    264262    """Effective radius""" 
    265     n_shells = int(n_shells + 0.5) 
     263    n_shells = int(n_shells) 
    266264    total = (np.sum(thickness[:n_shells], axis=1) 
    267265             + np.sum(interface[:n_shells], axis=1)) 

  • sasmodels/models/stacked_disks.c

    r19f996b r6c3e266  
    1 static 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) 
     1double form_volume(double thick_core, 
     2                   double thick_layer, 
     3                   double radius, 
     4                   double n_stacking); 
     5 
     6double 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 
     16double 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 
     28static 
     29double _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) 
    1441 
    1542{ 
    … …  
    6188 
    6289 
    63 static 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) 
     90static 
     91double 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) 
    73100{ 
    74101/*    StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks 
    … …  
    84111        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    85112        SINCOS(zi, sin_alpha, cos_alpha); 
    86         double yyy = stacked_disks_kernel(q, 
    87                            halfheight, 
    88                            thick_layer, 
     113        double yyy = _kernel(q, 
    89114                           radius, 
    90                            n_stacking, 
    91                            sigma_dnn, 
    92115                           core_sld, 
    93116                           layer_sld, 
    94117                           solvent_sld, 
     118                           halfheight, 
     119                           thick_layer, 
    95120                           sin_alpha, 
    96121                           cos_alpha, 
    97                            d); 
     122                           sigma_dnn, 
     123                           d, 
     124                           n_stacking); 
    98125        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    99126    } 
    … …  
    105132} 
    106133 
    107 static 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); 
     134double form_volume(double thick_core, 
     135                   double thick_layer, 
     136                   double radius, 
     137                   double n_stacking){ 
    114138    double d = 2.0 * thick_layer + thick_core; 
    115139    return M_PI * radius * radius * d * n_stacking; 
    116140} 
    117141 
    118 static 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) 
    128 { 
    129     int n_stacking = (int)(fp_n_stacking + 0.5); 
    130     return stacked_disks_1d(q, 
     142double 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) 
     151{ 
     152    return stacked_disks_kernel(q, 
    131153                    thick_core, 
    132154                    thick_layer, 
    … …  
    140162 
    141163 
    142 static 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) 
    153 { 
    154     int n_stacking = (int)(fp_n_stacking + 0.5); 
     164double 
     165Iqxy(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) 
     176{ 
    155177    double q, sin_alpha, cos_alpha; 
    156178    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    … …  
    158180    double d = 2.0 * thick_layer + thick_core; 
    159181    double halfheight = 0.5*thick_core; 
    160     double answer = stacked_disks_kernel(q, 
    161                      halfheight, 
    162                      thick_layer, 
     182    double answer = _kernel(q, 
    163183                     radius, 
    164                      n_stacking, 
    165                      sigma_dnn, 
    166184                     core_sld, 
    167185                     layer_sld, 
    168186                     solvent_sld, 
     187                     halfheight, 
     188                     thick_layer, 
    169189                     sin_alpha, 
    170190                     cos_alpha, 
    171                      d); 
     191                     sigma_dnn, 
     192                     d, 
     193                     n_stacking); 
    172194 
    173195    //convert to [cm-1] 

  • sasmodels/models/stacked_disks.py

    rb7e8b94 rb7e8b94  
    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, [1, inf],    "volume",      "Number of stacked layer/core/layer disks"], 
     128    ["n_stacking",  "",            1.0, [0, 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

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

  • sasmodels/models/unified_power_Rg.py

    r66ca2a6 r66ca2a6  
    9797 
    9898def Iq(q, level, rg, power, B, G): 
    99     level = int(level + 0.5) 
    100     if level == 0: 
     99    ilevel = int(level) 
     100    if ilevel == 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(level): 
     106        for i in range(ilevel): 
    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 < level-1: 
     109            if i < ilevel-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[:level]) 
     115    result[q == 0] = np.sum(G[:ilevel]) 
    116116    return result 
    117117 

  • sasmodels/product.py

    rf88e248 r9951a86  
    4545    # structure factor calculator.  Structure factors should not 
    4646    # have any magnetic parameters 
    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") 
     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 == []) 
    5350    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters 
    5451    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
    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  
     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:] 
    6864    translate_name = dict((old.id, new.id) for old, new 
    6965                          in zip(s_pars.kernel_parameters[1:], s_list)) 
    7066    demo = {} 
    71     demo.update(p_info.demo.items()) 
     67    demo.update((k, v) for k, v in p_info.demo.items() 
     68                if k not in ("background", "scale")) 
    7269    demo.update((translate_name[k], v) for k, v in s_info.demo.items() 
    7370                if k not in ("background", "scale") and not k.startswith(ER_ID)) 
    … …  
    9390    # Remember the component info blocks so we can build the model 
    9491    model_info.composition = ('product', [p_info, s_info]) 
    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)) 
     92    model_info.demo = {} 
    10293    return model_info 
    10394 
    104 def _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     """ 
     95def suffix_parameter(par, suffix): 
    11196    par = copy(par) 
    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):] 
     97    par.name = par.name + " S" 
    11698    par.id = par.id + "_S" 
    117     par.name = par.id + vector_length 
    11899    return par 
    119100 

  • sasmodels/resolution.py

    rb32caab rb397165  
    1717 
    1818MINIMUM_RESOLUTION = 1e-8 
    19 MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data 
     19 
     20 
     21# When extrapolating to -q, what is the minimum positive q relative to q_min 
     22# that we wish to calculate? 
     23MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION = 0.01 
    2024 
    2125class Resolution(object): 
    … …  
    7882        self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma) 
    7983                       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 
    8784        self.weight_matrix = pinhole_resolution(self.q_calc, self.q, 
    8885                                np.maximum(q_width, MINIMUM_RESOLUTION)) 
    89         self.q_calc = abs(self.q_calc) 
    9086 
    9187    def apply(self, theory): 
    … …  
    127123        self.q_calc = slit_extend_q(q, qx_width, qy_width) \ 
    128124            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 
    136125        self.weight_matrix = \ 
    137126            slit_resolution(self.q_calc, self.q, qx_width, qy_width) 
    138         self.q_calc = abs(self.q_calc) 
    139127 
    140128    def apply(self, theory): 
    … …  
    165153    # neither trapezoid nor Simpson's rule improved the accuracy. 
    166154    edges = bin_edges(q_calc) 
    167     #edges[edges < 0.0] = 0.0 # clip edges below zero 
     155    edges[edges < 0.0] = 0.0 # clip edges below zero 
    168156    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    169157    weights = cdf[1:] - cdf[:-1] 
    … …  
    298286    # The current algorithm is a midpoint rectangle rule. 
    299287    q_edges = bin_edges(q_calc) # Note: requires q > 0 
    300     #q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 
     288    q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 
    301289    weights = np.zeros((len(q), len(q_calc)), 'd') 
    302290 
    … …  
    404392    interval. 
    405393 
    406     Note that extrapolated values may be negative. 
     394    if *q_min* is zero or less then *q[0]/10* is used instead. 
    407395    """ 
    408396    q = np.sort(q) 
    409397    if q_min + 2*MINIMUM_RESOLUTION < q[0]: 
     398        if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    410399        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15 
    411400        q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 
    … …  
    459448        log_delta_q = log(10.) / points_per_decade 
    460449    if q_min < q[0]: 
    461         if q_min < 0: q_min = q[0]*MINIMUM_ABSOLUTE_Q 
     450        if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    462451        n_low = log_delta_q * (log(q[0])-log(q_min)) 
    463452        q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 

  • sasmodels/sasview_model.py

    r749a7d4 r64614ad  
    1616import logging 
    1717from os.path import basename, splitext 
    18 import thread 
    1918 
    2019import numpy as np  # type: ignore 
    … …  
    3938except ImportError: 
    4039    pass 
    41  
    42 calculation_lock = thread.allocate_lock() 
    4340 
    4441SUPPORT_OLD_STYLE_PLUGINS = True 
    … …  
    608605        to the card for each evaluation. 
    609606        """ 
    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): 
    620607        #core.HAVE_OPENCL = False 
    621608        if self._model is None: 
    … …  
    734721            return [self.params[par.name]], [1.0] 
    735722 
    736 def test_cylinder(): 
     723def test_model(): 
    737724    # type: () -> float 
    738725    """ 
    739     Test that the cylinder model runs, returning the value at [0.1,0.1]. 
     726    Test that a sasview model (cylinder) can be run. 
    740727    """ 
    741728    Cylinder = _make_standard_model('cylinder') 
    … …  
    746733    # type: () -> float 
    747734    """ 
    748     Test that 2-D hardsphere model runs and doesn't produce NaN. 
     735    Test that a sasview model (cylinder) can be run. 
    749736    """ 
    750737    Model = _make_standard_model('hardsphere') 
    … …  
    757744    # type: () -> float 
    758745    """ 
    759     Test that the 2-D RPA model runs 
     746    Test that a sasview model (cylinder) can be run. 
    760747    """ 
    761748    RPA = _make_standard_model('rpa') 
    … …  
    763750    return rpa.evalDistribution([0.1, 0.1]) 
    764751 
    765 def 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" 
    776752 
    777753def test_model_list(): 
    778754    # type: () -> None 
    779755    """ 
    780     Make sure that all models build as sasview models 
     756    Make sure that all models build as sasview models. 
    781757    """ 
    782758    from .exception import annotate_exception 
    … …  
    805781 
    806782if __name__ == "__main__": 
    807     print("cylinder(0.1,0.1)=%g"%test_cylinder()) 
    808     #test_empty_distribution() 
     783    print("cylinder(0.1,0.1)=%g"%test_model()) 

  • sasmodels/sesans.py

    r94d13f1 rb397165  
    1414import numpy as np  # type: ignore 
    1515from numpy import pi, exp  # type: ignore 
    16 from scipy.special import j0 
    17  
    18 class SesansTransform(object): 
    19     """ 
    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 
    23  
    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. 
    31     """ 
    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 
    35  
    36     #: q values to calculate when computing transform 
    37     q_calc = None  # type: np.ndarray 
    38  
    39     # transform arrays 
    40     _H = None  # type: np.ndarray 
    41     _H0 = None # type: np.ndarray 
    42  
    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) 
    48  
    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 
    55  
    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') 
    60  
    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 
    66  
    67         H0 = np.float32(dq/(2*pi)) * q 
    68  
    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 
    72  
    73         self.q_calc = q 
    74         self._H, self._H0 = H, H0 
     16from scipy.special import jv as besselj 
     17#import direct_model.DataMixin as model 
     18         
     19def 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. 
     27    """ 
     28    from sas.sascalc.data_util.nxsunit import Converter 
     29 
     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     
     36def make_all_q(data): 
     37    """ 
     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] 
     55 
     56def 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) 
     76 
     77    return result 
     78 
     79def 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   
     86def 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                   
     91def 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                         
     96def 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     
     102 
     103 
     104def 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 
     109 
     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 
     131 
     132 
     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] 
     136 
     137 
     138 
     139 
     140def 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     
     168def 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 
Note: See TracChangeset for help on using the changeset viewer.