Changes in / [37c7d5e:ae32bb8] in sasmodels
- Files:
-
- 34 edited
Legend:
- Unmodified
- Added
- Removed
-
example/fit.py
rf4b36fa r1182da5 24 24 % section) 25 25 data = radial_data if section != "tangential" else tan_data 26 theta = 89.9 if section != "tangential" else 0 27 phi = 90 26 phi = 0 if section != "tangential" else 90 28 27 kernel = load_model(name, dtype="single") 29 28 cutoff = 1e-3 … … 31 30 if name == "ellipsoid": 32 31 model = Model(kernel, 33 scale=0.08, background=35,34 r adius_polar=15, radius_equatorial=800,32 scale=0.08, 33 r_polar=15, r_equatorial=800, 35 34 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, 38 40 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,41 41 ) 42 42 43 43 44 # 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) 50 48 model.background.range(0,1000) 51 49 model.scale.range(0, 10) 52 50 53 51 52 54 53 elif name == "lamellar": 55 54 model = Model(kernel, 56 scale=0.08, background=0.003,55 scale=0.08, 57 56 thickness=19.2946, 58 57 sld=5.38,sld_sol=7.105, 58 background=0.003, 59 59 thickness_pd= 0.37765, thickness_pd_n=10, thickness_pd_nsigma=3, 60 60 ) 61 61 62 62 63 # SET THE FITTING PARAMETERS … … 76 77 radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3, 77 78 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,) 79 80 """ 80 81 pars = dict( 81 82 scale=.01, background=35, 82 83 sld=.291, sld_solvent=5.77, 83 radius=250, length=178, 84 radius=250, length=178, 85 theta=90, phi=phi, 84 86 radius_pd=0.1, radius_pd_n=5, radius_pd_nsigma=3, 85 87 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) 89 90 model = Model(kernel, **pars) 90 91 … … 92 93 model.radius.range(1, 500) 93 94 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 97 98 model.radius_pd.range(0, 1) 98 model.length_pd.range(0, 1)99 model.length_pd.range(0, 2) 99 100 model.scale.range(0, 10) 100 101 model.background.range(0, 100) … … 103 104 elif name == "core_shell_cylinder": 104 105 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 108 110 radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3, 109 111 length_pd=0.26, length_pd_n=10, length_pd_nsigma=3, 110 112 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, 114 115 ) 115 116 116 117 # 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) 119 120 #model.thickness.range(18, 38) 120 121 #model.thickness_pd.range(0, 1) 121 122 #model.phi.range(0, 90) 122 model.phi_pd.range(0,20)123 123 #model.radius_pd.range(0, 1) 124 124 #model.length_pd.range(0, 1) … … 131 131 elif name == "capped_cylinder": 132 132 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, 135 134 sld=1, sld_solvent=6.3, 135 background=0, theta=0, phi=phi, 136 136 radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3, 137 137 cap_radius_pd=.1, cap_radius_pd_n=5, cap_radius_pd_nsigma=3, 138 138 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, 142 141 ) 143 142 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)154 143 model.scale.range(0, 1) 155 144 … … 157 146 elif name == "triaxial_ellipsoid": 158 147 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, 161 149 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, 166 151 theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3, 167 152 phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 168 153 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, 169 157 ) 170 158 171 159 # SET THE FITTING PARAMETERS 172 model.r adius_equat_minor.range(15, 1000)173 model.r adius_equat_major.range(15, 1000)174 #model.r adius_polar.range(15, 1000)160 model.req_minor.range(15, 1000) 161 model.req_major.range(15, 1000) 162 #model.rpolar.range(15, 1000) 175 163 #model.background.range(0,1000) 176 164 #model.theta_pd.range(0, 360) … … 185 173 M = Experiment(data=data, model=model) 186 174 if section == "both": 187 tan_model = Model(model. sasmodel, **model.parameters())175 tan_model = Model(model.kernel, **model.parameters()) 188 176 tan_model.phi = model.phi - 90 189 177 tan_model.cutoff = cutoff -
sasmodels/compare.py
rf72d70a rfe25eda 322 322 name = name.split('*')[0] 323 323 324 # Suppress magnetism for python models (not yet implemented)325 if callable(model_info.Iq):326 pars.update(suppress_magnetism(pars))327 328 324 if name == 'barbell': 329 325 if pars['radius_bell'] < pars['radius']: … … 344 340 if pars['radius'] < pars['thick_string']: 345 341 pars['radius'], pars['thick_string'] = pars['thick_string'], pars['radius'] 342 pars['num_pearls'] = math.ceil(pars['num_pearls']) 346 343 pass 347 344 … … 356 353 for c in '1234': 357 354 pars['Phi'+c] /= total 355 356 elif name == 'stacked_disks': 357 pars['n_stacking'] = math.ceil(pars['n_stacking']) 358 358 359 359 def parlist(model_info, pars, is2d): -
sasmodels/compare_many.py
rf72d70a r424fe00 106 106 header = ('\n"Model","%s","Count","%d","Dimension","%s"' 107 107 % (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,) 110 109 print(header) 111 110 … … 162 161 max_diff = [0] 163 162 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) 165 164 seed = np.random.randint(1e6) 166 165 pars_i = randomize_pars(model_info, pars, seed) 167 166 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) 170 168 if mono: 171 169 pars_i = suppress_pd(pars_i) … … 189 187 Print the command usage string. 190 188 """ 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)") 193 190 194 191 … … 207 204 print("""\ 208 205 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. 206 MODEL is the model name of the model or "all" for all the models 207 in alphabetical order. 212 208 213 209 COUNT is the number of randomly generated parameter sets to try. A value … … 224 220 below the cutoff will be ignored. Use "mono" for monodisperse models. The 225 221 choice 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. 222 is set in compare.py defaults for each model. 228 223 229 224 PRECISION 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. 225 precisions are given, then compare one to the other, ignoring sasview. 233 226 234 227 Available models: … … 240 233 Main program. 241 234 """ 242 if len(argv) not in ( 3, 4,5, 6):235 if len(argv) not in (5, 6): 243 236 print_help() 244 237 return 245 238 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) 251 242 print_models() 252 print('model types: all, py, c, double, single, opencl, 1d, 2d, nonmagnetic, magnetic')253 243 return 254 244 try: … … 257 247 assert argv[2][1] == 'd' 258 248 Nq = int(argv[2][2:]) 259 mono = len(argv) <= 3 orargv[3] == 'mono'249 mono = argv[3] == 'mono' 260 250 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" 263 253 except Exception: 264 254 traceback.print_exc() … … 268 258 data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 269 259 'accuracy': 'Low', 'view':'log', 'zero': False}) 260 model_list = [model] if model != "all" else MODELS 270 261 for model in model_list: 271 262 compare_instance(model, data, index, N=count, mono=mono, -
sasmodels/conversion_table.py
rd3e3f756 rbb584b3 549 549 "radius": "core_radius", 550 550 "sld_solvent": "core_sld", 551 "n_ shells": "n_pairs",551 "n_pairs": "n_pairs", 552 552 "thick_shell": "s_thickness", 553 553 "sld": "shell_sld", -
sasmodels/core.py
r5124c969 r52e9a45 69 69 * magnetic: models with an sld 70 70 * 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: 77 73 raise ValueError("kind not in " + ", ".join(KINDS)) 78 74 files = sorted(glob(joinpath(generate.MODEL_PATH, "[a-zA-Z]*.py"))) 79 75 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)] 86 77 87 78 return selected -
sasmodels/model_test.py
r598857b r1a6cd57 97 97 is_py = callable(model_info.Iq) 98 98 99 # Some OpenCL drivers seem to be flaky, and are not producing the100 # expected result. Since we don't have known test values yet for101 # all of our models, we are instead going to compare the results102 # for the 'smoke test' (that is, evaluation at q=0.1 for the default103 # parameters just to see that the model runs to completion) between104 # the OpenCL and the DLL. To do this, we define a 'stash' which is105 # shared between OpenCL and DLL tests. This is just a list. If the106 # list is empty (which it will be when DLL runs, if the DLL runs107 # first), then the results are appended to the list. If the list108 # is not empty (which it will be when OpenCL runs second), the results109 # 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 the111 # test suite is thrown away after being run once.112 stash = []113 114 99 if is_py: # kernel implemented in python 115 100 test_name = "Model: %s, Kernel: python"%model_name … … 118 103 test_method_name, 119 104 platform="dll", # so that 120 dtype="double", 121 stash=stash) 105 dtype="double") 122 106 suite.addTest(test) 123 107 else: # kernel implemented in C 124 125 # test using dll if desired126 if 'dll' in loaders or not core.HAVE_OPENCL:127 test_name = "Model: %s, Kernel: dll"%model_name128 test_method_name = "test_%s_dll" % model_info.id129 test = ModelTestCase(test_name, model_info,130 test_method_name,131 platform="dll",132 dtype="double",133 stash=stash)134 suite.addTest(test)135 136 108 # test using opencl if desired and available 137 109 if 'opencl' in loaders and core.HAVE_OPENCL: … … 144 116 test = ModelTestCase(test_name, model_info, 145 117 test_method_name, 146 platform="ocl", dtype=None, 147 stash=stash) 118 platform="ocl", dtype=None) 148 119 #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") 149 130 suite.addTest(test) 150 131 … … 163 144 """ 164 145 def __init__(self, test_name, model_info, test_method_name, 165 platform, dtype , stash):166 # type: (str, ModelInfo, str, str, DType , List[Any]) -> None146 platform, dtype): 147 # type: (str, ModelInfo, str, str, DType) -> None 167 148 self.test_name = test_name 168 149 self.info = model_info 169 150 self.platform = platform 170 151 self.dtype = dtype 171 self.stash = stash # container for the results of the first run172 152 173 153 setattr(self, test_method_name, self.run_all) … … 187 167 #({}, (0.0, 0.0), None), 188 168 # test vector form 189 ({}, [0. 001, 0.01, 0.1], [None]*3),169 ({}, [0.1]*2, [None]*2), 190 170 ({}, [(0.1, 0.1)]*2, [None]*2), 191 171 # test that ER/VR will run if they exist … … 194 174 ] 195 175 196 tests = s moke_tests + self.info.tests176 tests = self.info.tests 197 177 try: 198 178 model = build_model(self.info, dtype=self.dtype, 199 179 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 217 190 218 191 except: 219 192 annotate_exception(self.test_name) 220 193 raise 221 222 def _find_missing_tests(self):223 # type: () -> None224 """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 = True228 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 None232 single = [test for test in self.info.tests233 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.tests240 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 missing258 194 259 195 def run_one(self, model, test): … … 271 207 272 208 if x[0] == 'ER': 273 actual = np.array([call_ER(model.info, pars)])209 actual = [call_ER(model.info, pars)] 274 210 elif x[0] == 'VR': 275 actual = np.array([call_VR(model.info, pars)])211 actual = [call_VR(model.info, pars)] 276 212 elif isinstance(x[0], tuple): 277 213 qx, qy = zip(*x) … … 302 238 'f(%s); expected:%s; actual:%s' 303 239 % (xi, yi, actual_yi)) 304 return actual305 240 306 241 return ModelTestCase -
sasmodels/modelinfo.py
rf88e248 r85fe7f8 230 230 defined as a sublist with the following elements: 231 231 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 233 234 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. 239 236 240 237 *units* should be one of *degrees* for angles, *Ang* for lengths, … … 606 603 # Using the call_parameters table, we already have expanded forms 607 604 # 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) 610 606 611 607 def append_group(name): … … 734 730 info.docs = kernel_module.__doc__ 735 731 info.category = getattr(kernel_module, 'category', None) 732 info.single = getattr(kernel_module, 'single', True) 733 info.opencl = getattr(kernel_module, 'opencl', True) 736 734 info.structure_factor = getattr(kernel_module, 'structure_factor', False) 737 735 info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) … … 747 745 info.profile = getattr(kernel_module, 'profile', None) # type: ignore 748 746 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))752 747 753 748 # multiplicity info -
sasmodels/models/core_multi_shell.c
r925ad6e r925ad6e 8 8 } 9 9 10 static double 11 form_volume(double core_radius, double fp_n, double thickness[]) 10 double 11 form_volume(double core_radius, double n, double thickness[]); 12 double 13 form_volume(double core_radius, double n, double thickness[]) 12 14 { 13 15 double r = core_radius; 14 int n = (int)(fp_n+0.5);15 16 for (int i=0; i < n; i++) { 16 17 r += thickness[i]; … … 19 20 } 20 21 21 staticdouble22 double 22 23 Iq(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[]); 25 double 26 Iq(double q, double core_sld, double core_radius, 27 double solvent_sld, double num_shells, double sld[], double thickness[]) 24 28 { 25 const int n = (int) (fp_n+0.5);29 const int n = (int)ceil(num_shells); 26 30 double f, r, last_sld; 27 31 r = core_radius; -
sasmodels/models/core_multi_shell.py
r5a0b3d7 r925ad6e 107 107 Returns the SLD profile *r* (Ang), and *rho* (1e-6/Ang^2). 108 108 """ 109 n = int(n+0.5)110 109 z = [] 111 110 rho = [] … … 134 133 def ER(radius, n, thickness): 135 134 """Effective radius""" 136 n = int(n[0] +0.5) # n is a control parameter and is notpolydisperse135 n = int(n[0]) # n cannot be polydisperse 137 136 return np.sum(thickness[:n], axis=0) + radius 138 137 -
sasmodels/models/flexible_cylinder.c
r592343f r592343f 1 1 static double 2 form_volume( double length, double kuhn_length, doubleradius)2 form_volume(length, kuhn_length, radius) 3 3 { 4 4 return 1.0; -
sasmodels/models/lamellar_hg_stack_caille.c
r1c6286d ra807206 3 3 */ 4 4 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) 5 double Iq(double qval, 6 double length_tail, 7 double length_head, 8 double Nlayers, 9 double dd, 10 double Cp, 11 double tail_sld, 12 double head_sld, 13 double solvent_sld); 14 15 double Iq(double qval, 16 double length_tail, 17 double length_head, 18 double Nlayers, 19 double dd, 20 double Cp, 21 double tail_sld, 22 double head_sld, 23 double solvent_sld) 15 24 { 16 int Nlayers = (int)(fp_Nlayers+0.5); //cast to an integer for the loop25 double NN; //local variables of coefficient wave 17 26 double inten,Pq,Sq,alpha,temp,t2; 18 27 //double dQ, dQDefault, t1, t3; 28 int ii,NNint; 19 29 // from wikipedia 0.577215664901532860606512090082402431042159335 20 30 const double Euler = 0.577215664901533; // Euler's constant, increased sig figs for new models Feb 2015 … … 22 32 //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015 23 33 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); 26 38 Pq *= Pq; 27 39 Pq *= 4.0/(qval*qval); 28 40 41 NNint = (int)NN; //cast to an integer for the loop 42 ii=0; 29 43 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 31 48 temp = 0.0; 32 49 alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler); … … 35 52 //t3 = dQ*dQ*dd*dd*ii*ii; 36 53 37 temp = 1.0- (double)ii/(double)Nlayers;54 temp = 1.0-ii/NN; 38 55 //temp *= cos(dd*qval*ii/(1.0+t1)); 39 56 temp *= cos(dd*qval*ii); 40 57 //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii); 41 58 //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 42 temp *= exp(-t2/2.0 );59 temp *= exp(-t2/2.0 ); 43 60 //temp /= sqrt(1.0+t1); 44 61 … … 54 71 55 72 inten *= 1.0e-04; // 1/A to 1/cm 56 return inten;73 return(inten); 57 74 } 58 75 -
sasmodels/models/lamellar_hg_stack_caille.py
ra57b31d r7c57861 98 98 ["length_head", "Ang", 2, [0, inf], "volume", 99 99 "head thickness"], 100 ["Nlayers", "", 30, [ 1, inf], "",100 ["Nlayers", "", 30, [0, inf], "", 101 101 "Number of layers"], 102 102 ["d_spacing", "Ang", 40., [0.0, inf], "volume", -
sasmodels/models/lamellar_stack_caille.c
r1c6286d r0bef47b 3 3 */ 4 4 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) 5 double Iq(double qval, 6 double del, 7 double Nlayers, 8 double dd, 9 double Cp, 10 double sld, 11 double solvent_sld); 12 13 double Iq(double qval, 14 double del, 15 double Nlayers, 16 double dd, 17 double Cp, 18 double sld, 19 double solvent_sld) 13 20 { 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 16 22 double inten,Pq,Sq,alpha,temp,t2; 17 23 //double dQ, dQDefault, t1, t3; 24 int ii,NNint; 18 25 // from wikipedia 0.577215664901532860606512090082402431042159335 19 26 const double Euler = 0.577215664901533; // Euler's constant, increased sig figs for new models Feb 2015 … … 21 28 //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015 22 29 30 NN = trunc(Nlayers); //be sure that NN is an integer 31 23 32 contr = sld - solvent_sld; 24 33 25 34 Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)); 26 35 36 NNint = (int)NN; //cast to an integer for the loop 37 ii=0; 27 38 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 29 44 temp = 0.0; 30 45 alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler); … … 33 48 //t3 = dQ*dQ*dd*dd*ii*ii; 34 49 35 temp = 1.0 - (double)ii / (double)Nlayers;50 temp = 1.0-ii/NN; 36 51 //temp *= cos(dd*qval*ii/(1.0+t1)); 37 52 temp *= cos(dd*qval*ii); 38 53 //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 39 temp *= exp(-t2/2.0 );54 temp *= exp(-t2/2.0 ); 40 55 //temp /= sqrt(1.0+t1); 41 56 -
sasmodels/models/lamellar_stack_caille.py
ra57b31d r7c57861 88 88 parameters = [ 89 89 ["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"], 91 91 ["d_spacing", "Ang", 400., [0.0,inf], "volume", "lamellar d-spacing of Caille S(Q)"], 92 92 ["Caille_parameter", "1/Ang^2", 0.1, [0.0,0.8], "", "Caille parameter"], -
sasmodels/models/lamellar_stack_paracrystal.c
r5467cd8 r4962519 2 2 3 3 */ 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); 4 double Iq(double qval, 5 double th, 6 double Nlayers, 7 double davg, 8 double pd, 9 double sld, 10 double solvent_sld); 11 double paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an); 12 double paraCryst_an(double ww, double qval, double davg, long Nlayers); 6 13 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) 14 double Iq(double qval, 15 double th, 16 double Nlayers, 17 double davg, 18 double pd, 19 double sld, 20 double solvent_sld) 15 21 { 22 23 double inten,contr,xn; 24 double xi,ww,Pbil,Znq,Snq,an; 25 long n1,n2; 26 27 contr = sld - solvent_sld; 16 28 //get the fractional part of Nlayers, to determine the "mixing" of N's 17 int n1 = (int)(fp_Nlayers); //truncate towards zero18 int n2 = n1 + 1;19 const double xn = (double)n2 - fp_Nlayers; //fractional contribution of n120 29 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); 22 35 23 36 //calculate the n1 contribution 24 double Znq,Snq,an;25 37 an = paraCryst_an(ww,qval,davg,n1); 26 38 Snq = paraCryst_sn(ww,qval,davg,n1,an); 27 39 28 40 Znq = xn*Snq; 29 41 … … 40 52 // Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg)) 41 53 42 const doublexi = th/2.0; //use 1/2 the bilayer thickness43 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)); 44 56 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; 47 59 //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); 49 61 } 50 62 51 63 // functions for the lamellar paracrystal model 52 64 double 53 paraCryst_sn(double ww, double qval, double davg, intNlayers, double an) {65 paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an) { 54 66 55 67 double Snq; … … 57 69 Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) ); 58 70 59 return Snq;71 return(Snq); 60 72 } 61 73 62 74 double 63 paraCryst_an(double ww, double qval, double davg, int Nlayers) { 75 paraCryst_an(double ww, double qval, double davg, long Nlayers) { 76 64 77 double an; 65 78 … … 69 82 an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg); 70 83 71 return an;84 return(an); 72 85 } 73 86 -
sasmodels/models/lamellar_stack_paracrystal.py
ra0168e8 r7c57861 113 113 parameters = [["thickness", "Ang", 33.0, [0, inf], "volume", 114 114 "sheet thickness"], 115 ["Nlayers", "", 20, [ 1, inf], "",115 ["Nlayers", "", 20, [0, inf], "", 116 116 "Number of layers"], 117 117 ["d_spacing", "Ang", 250., [0.0, inf], "", -
sasmodels/models/linear_pearls.c
r925ad6e r925ad6e 4 4 double radius, 5 5 double edge_sep, 6 double fp_num_pearls,6 double num_pearls, 7 7 double pearl_sld, 8 8 double solvent_sld); … … 11 11 double radius, 12 12 double edge_sep, 13 intnum_pearls,13 double num_pearls, 14 14 double pearl_sld, 15 15 double solvent_sld); 16 16 17 17 18 double form_volume(double radius, double fp_num_pearls)18 double form_volume(double radius, double num_pearls) 19 19 { 20 int num_pearls = (int)(fp_num_pearls + 0.5);21 20 // Pearl volume 22 21 double pearl_vol = M_4PI_3 * cube(radius); … … 28 27 double radius, 29 28 double edge_sep, 30 intnum_pearls,29 double num_pearls, 31 30 double pearl_sld, 32 31 double solvent_sld) 33 32 { 33 double n_contrib; 34 34 //relative sld 35 35 double contrast_pearl = pearl_sld - solvent_sld; … … 46 46 double psi = sas_3j1x_x(q * radius); 47 47 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)); 52 53 } 53 54 // 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; 55 56 56 57 return form_factor; … … 60 61 double radius, 61 62 double edge_sep, 62 double fp_num_pearls,63 double num_pearls, 63 64 double pearl_sld, 64 65 double solvent_sld) 65 66 { 66 67 67 int num_pearls = (int)(fp_num_pearls + 0.5);68 68 double result = linear_pearls_kernel(q, 69 69 radius, -
sasmodels/models/linear_pearls.py
r925ad6e r925ad6e 16 16 .. math:: 17 17 18 P(Q) = \frac{ \text{scale}}{V}\left[ m_{p}^219 \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] 21 21 22 22 where the mass $m_p$ is $(SLD_{pearl}-SLD_{solvent})*(volume\ of\ N\ pearls)$. … … 56 56 ["radius", "Ang", 80.0, [0, inf], "", "Radius of the pearls"], 57 57 ["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"], 59 59 ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", "SLD of the pearl spheres"], 60 60 ["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, 1 static 2 double multilayer_vesicle_kernel(double q, 14 3 double volfraction, 15 4 double radius, … … 18 7 double sld_solvent, 19 8 double sld, 20 int n_shells)9 double n_pairs) 21 10 { 22 11 //calculate with a loop, two shells at a time … … 40 29 41 30 //do 2 layers at a time 42 ii ++;31 ii += 1; 43 32 44 } while(ii <= n_ shells-1); //change to make 0 < n_shells < 2 correspond to33 } while(ii <= n_pairs-1); //change to make 0 < n_pairs < 2 correspond to 45 34 //unilamellar vesicles (C. Glinka, 11/24/03) 46 35 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); 48 39 } 49 40 50 static double51 Iq(double q,41 static 42 double Iq(double q, 52 43 double volfraction, 53 44 double radius, … … 56 47 double sld_solvent, 57 48 double sld, 58 double fp_n_shells)49 double n_pairs) 59 50 { 60 int n_shells = (int)(fp_n_shells + 0.5);61 51 return multilayer_vesicle_kernel(q, 62 52 volfraction, … … 66 56 sld_solvent, 67 57 sld, 68 n_ shells);58 n_pairs); 69 59 } 70 60 -
sasmodels/models/multilayer_vesicle.py
r5d23de2 r925ad6e 3 3 ---------- 4 4 5 This model is a trivial extension of the core_shell_sphere function where the6 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 fora multilayer vesicle.5 This model is a trivial extension of the core_shell_sphere function to include 6 *N* shells where the core is filled with solvent and the shells are interleaved 7 with layers of solvent. For *N = 1*, this returns the same as the vesicle model, 8 except for the normalisation, which here is to outermost volume. 9 The shell thicknessess and SLD are constant for all shells as expected for 10 a multilayer vesicle. 11 11 12 12 .. figure:: img/multi_shell_geometry.jpg … … 19 19 20 20 .. 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} 22 23 23 24 where 24 25 25 26 .. 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} 29 31 \right] 30 32 31 for32 33 33 .. math:: 34 where $R_i = r_c + (i-1)(t_s + t_w)$ 35 36 where $V_t$ is the volume of the whole particle, $V(R)$ is the volume of a sphere 37 of radius $R$, $r_c$ is the radius of the core, $\rho_{shell}$ is the scattering length 38 density of a shell, $\rho_{solv}$ is the scattering length density of the solvent. 34 39 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 } i37 38 $\phi$ is the volume fraction of particles, $V(r)$ is the volume of a sphere39 of radius $r$, $r_c$ is the radius of the core, $t_s$ is the thickness of40 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, and42 $\rho_\text{solv}$ is the scattering length density of the solvent.43 44 USAGE NOTES45 46 * The outer-most shell radius $R_N$ is used as the effective radius47 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 interger50 number of layers is not physical.51 * Thus Polydispersity should only be applied to number of shells **VERY52 CAREFULLY**. A possible legitimate use would be for mixed systems in which53 some vesicles have 1 shell, some have 2, etc. A polydispersity on $N$ can be54 used to model the data by using the "array distriubtion" feature. First55 create a file such as *shell_dist.txt* containing the relative portion56 of each vesicle size::57 58 1 2059 2 460 3 161 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 be64 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 the67 q-values that correspond to the repeat distance of the layers), model68 function complicated by the fact that the number of water/shell pairs must69 physically be an integer value, although the optimization treats it as a70 floating point value. Thus it may be that the resolution interpolation is not71 sufficiently fine grained in certain cases. Please report any such occurences72 to the SasView team. Generally, for the best possible experience:73 * Start with the best possible guess74 * Using a priori knowledge, hold as many parameters fixed as possible75 * if N=1, tw (water thickness) must by definition be zero. Both N and tw should76 be fixed during fitting.77 * If N>1, use constraints to keep N > 178 * Because N only really moves in integer steps, it may get "stuck" if the79 optimizer step size is too small so care should be taken80 If you experience problems with this please contact the SasView team and let81 them know the issue preferably with example data and model which fail to82 converge.83 40 84 41 The 2D scattering intensity is the same as 1D, regardless of the orientation … … 89 46 q = \sqrt{q_x^2 + q_y^2} 90 47 48 49 The outer most radius 50 51 $radius + n\_pairs * thick\_shell + (n\_pairs- 1) * thick\_solvent$ 52 53 is used for both the volume fraction normalization and for the 54 effective radius for *S(Q)* when $P(Q) * S(Q)$ is applied. 55 91 56 For information about polarised and magnetic scattering, see 92 57 the :ref:`magnetism` documentation. 58 59 This code is based on the form factor calculations implemented in the NIST 60 Center for Neutron Research provided c-library (Kline, 2006). 93 61 94 62 References 95 63 ---------- 96 64 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). 65 B Cabane, *Small Angle Scattering Methods*, 66 in *Surfactant Solutions: New Methods of Investigation*, 67 Ch.2, Surfactant Science Series Vol. 22, Ed. R Zana and M Dekker, 68 New York, (1987). 100 69 101 Authorship and Verification 102 ---------------------------- 70 **Author:** NIST IGOR/DANSE **on:** pre 2010 103 71 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 108 75 109 76 """ … … 122 89 sld_solvent: solvent scattering length density 123 90 sld: shell scattering length density 124 n_ shells:number of "shell plus solvent" layer pairs91 n_pairs:number of "shell plus solvent" layer pairs 125 92 background: incoherent background 126 93 """ … … 131 98 parameters = [ 132 99 ["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"], 136 103 ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "sld", "solvent scattering length density"], 137 104 ["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"], 139 106 ] 140 107 # pylint: enable=bad-whitespace, line-too-long 141 108 142 # TODO: proposed syntax for specifying which parameters can be polydisperse143 #polydispersity = ["radius", "thick_shell"]144 145 109 source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 146 110 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 111 polydispersity = ["radius", "n_pairs"] 150 112 151 113 demo = dict(scale=1, background=0, … … 156 118 sld_solvent=6.4, 157 119 sld=0.4, 158 n_ shells=2.0)120 n_pairs=2.0) 159 121 160 122 tests = [ … … 165 127 'sld_solvent': 6.4, 166 128 'sld': 0.4, 167 'n_ shells': 2.0,129 'n_pairs': 2.0, 168 130 'scale': 1.0, 169 131 'background': 0.001, … … 176 138 'sld_solvent': 6.4, 177 139 'sld': 0.4, 178 'n_ shells': 2.0,140 'n_pairs': 2.0, 179 141 'scale': 1.0, 180 142 'background': 0.001, -
sasmodels/models/onion.c
r925ad6e r925ad6e 30 30 31 31 static double 32 form_volume(double radius_core, double n _shells, double thickness[])32 form_volume(double radius_core, double n, double thickness[]) 33 33 { 34 int n = (int)(n_shells+0.5);34 int i; 35 35 double r = radius_core; 36 for (i nt i=0; i < n; i++) {36 for (i=0; i < n; i++) { 37 37 r += thickness[i]; 38 38 } -
sasmodels/models/onion.py
r925ad6e r925ad6e 323 323 Returns shape profile with x=radius, y=SLD. 324 324 """ 325 n_shells = int(n_shells+0.5) 325 326 326 total_radius = 1.25*(sum(thickness[:n_shells]) + radius_core + 1) 327 327 dz = total_radius/400 # 400 points for a smooth plot … … 366 366 return np.asarray(z), np.asarray(rho) 367 367 368 def ER(radius_core, n _shells, thickness):368 def ER(radius_core, n, thickness): 369 369 """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 372 371 373 372 demo = { -
sasmodels/models/pearl_necklace.c
r3f853beb r4b541ac 1 1 double form_volume(double radius, double edge_sep, 2 double thick_string, double fp_num_pearls);2 double thick_string, double num_pearls); 3 3 double 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, 5 5 double string_sld, double solvent_sld); 6 6 … … 9 9 // From Igor library 10 10 static double 11 pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string,12 intnum_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) 13 13 { 14 14 // 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; 16 17 17 18 //each masses: contrast * volume … … 68 69 69 70 double form_volume(double radius, double edge_sep, 70 double thick_string, double fp_num_pearls)71 double thick_string, double num_pearls) 71 72 { 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; 74 76 const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; 75 77 const double pearl_vol = M_4PI_3 * radius * radius * radius; … … 80 82 81 83 double 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, 83 85 double string_sld, double solvent_sld) 84 86 { 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, 87 88 thick_string, num_pearls, sld, string_sld, solvent_sld); 88 89 -
sasmodels/models/pearl_necklace.py
r1bd1ea2 r4b541ac 82 82 ["thick_string", "Ang", 2.5, [0, inf], "volume", 83 83 "Thickness of the chain linkage"], 84 ["num_pearls", "none", 3, [ 1, inf], "volume",84 ["num_pearls", "none", 3, [0, inf], "volume", 85 85 "Number of pearls in the necklace (must be integer)"], 86 86 ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", … … 100 100 Redundant with form_volume. 101 101 """ 102 num_pearls = int(num_pearls + 0.5)103 102 number_of_strings = num_pearls - 1.0 104 103 string_vol = edge_sep * pi * pow((thick_string / 2.0), 2.0) … … 112 111 Calculation for effective radius. 113 112 """ 114 num_pearls = int(num_pearls + 0.5)115 113 tot_vol = volume(radius, edge_sep, thick_string, num_pearls) 116 114 rad_out = (tot_vol/(4.0/3.0*pi)) ** (1./3.) -
sasmodels/models/raspberry.py
r8e68ea0 r8e68ea0 10 10 Schematic of the raspberry model 11 11 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 largedroplet and small particles all need to be calculated.12 In order to calculate the form factor of the entire complex, the self- 13 correlation of the large droplet, the self-correlation of the particles, the 14 correlation terms between different particles and the cross terms between large 15 droplet and small particles all need to be calculated. 16 16 17 Consider two infinitely thin shells of radii $R_1$ and $R_2$ separated by18 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 theircenters.17 Consider two infinitely thin shells of radii R1 and R2 separated by distance r. 18 The general structure of the equation is then the form factor of the two shells 19 multiplied by the phase factor that accounts for the separation of their 20 centers. 21 21 22 22 .. math:: -
sasmodels/models/rpa.c
r20c856a r6351bfa 1 double Iq(double q, double fp_case_num,1 double Iq(double q, double case_num, 2 2 double N[], double Phi[], double v[], double L[], double b[], 3 3 double Kab, double Kac, double Kad, … … 5 5 ); 6 6 7 double Iq(double q, double fp_case_num,7 double Iq(double q, double case_num, 8 8 double N[], // DEGREE OF POLYMERIZATION 9 9 double Phi[], // VOL FRACTION … … 15 15 ) 16 16 { 17 int icase = (int) (fp_case_num+0.5);17 int icase = (int)case_num; 18 18 19 19 double Nab,Nac,Nad,Nbc,Nbd,Ncd; … … 39 39 double S31,S32,S33,S34,S41,S42,S43,S44; 40 40 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 43 42 //icase was shifted to N-1 from the original code 44 43 if (icase <= 1){ … … 58 57 Kab = Kac = Kad = -0.0004; 59 58 } 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 65 60 Nab=sqrt(N[0]*N[1]); 66 61 Nac=sqrt(N[0]*N[2]); … … 84 79 Phicd=sqrt(Phi[2]*Phi[3]); 85 80 86 // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk87 81 Xa=q*q*b[0]*b[0]*N[0]/6.0; 88 82 Xb=q*q*b[1]*b[1]*N[1]/6.0; … … 90 84 Xd=q*q*b[3]*b[3]*N[3]/6.0; 91 85 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); 96 89 S0ab=(Phiab*vab*Nab)*Pab; 97 Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor90 Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 98 91 S0ac=(Phiac*vac*Nac)*Pac; 99 Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block92 Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 100 93 S0ad=(Phiad*vad*Nad)*Pad; 101 94 102 95 S0ba=S0ab; 103 Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain96 Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); 104 97 S0bb=N[1]*Phi[1]*v[1]*Pbb; 105 Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock98 Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 106 99 S0bc=(Phibc*vbc*Nbc)*Pbc; 107 Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock100 Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 108 101 S0bd=(Phibd*vbd*Nbd)*Pbd; 109 102 110 103 S0ca=S0ac; 111 104 S0cb=S0bc; 112 Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain105 Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); 113 106 S0cc=N[2]*Phi[2]*v[2]*Pcc; 114 Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock107 Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 115 108 S0cd=(Phicd*vcd*Ncd)*Pcd; 116 109 … … 118 111 S0db=S0bd; 119 112 S0dc=S0cd; 120 Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain113 Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); 121 114 S0dd=N[3]*Phi[3]*v[3]*Pdd; 122 123 // Reset all unused partial structure factors to 0 (depends on case)124 115 //icase was shifted to N-1 from the original code 125 116 switch(icase){ … … 202 193 S0dc=S0cd; 203 194 204 // self chi parameter is 0 ... of course205 195 Kaa=0.0; 206 196 Kbb=0.0; … … 253 243 ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 254 244 255 // D is considered the matrix or background component so enters here256 245 m=1.0/(S0dd-ZZ); 257 246 … … 308 297 S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 309 298 310 //calculate contrast where L[i] is the scattering length of i and D is the matrix311 //Note that should multiply by Nav to get units of SLD which will become312 // Nav*2 in the next line (SLD^2) but then normalization in that line would313 //need to divide by Nav leaving only Nav or sqrt(Nav) before squaring.314 299 Nav=6.022045e+23; 315 300 Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); … … 318 303 319 304 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;323 305 324 306 return Intg; -
sasmodels/models/rpa.py
r4f9e288 r40a87fa 1 1 r""" 2 Definition3 ----------4 5 2 Calculates the macroscopic scattering intensity for a multi-component 6 3 homogeneous mixture of polymers using the Random Phase Approximation. … … 27 24 Case 9: A-B-C-D tetra-block copolymer 28 25 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!** 31 27 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]_. 28 Only one case can be used at any one time. 38 29 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. 30 The RPA (mean field) formalism only applies only when the multicomponent 31 polymer mixture is in the homogeneous mixed-phase region. 42 32 43 USAGE NOTES: 33 **Component D is assumed to be the "background" component (ie, all contrasts 34 are calculated with respect to component D).** So the scattering contrast 35 for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`. 44 36 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. 37 Depending on which case is being used, the number of fitting parameters - the 38 segment lengths (ba, bb, etc) and $\chi$ parameters (Kab, Kac, etc) - vary. 39 The *scale* parameter should be held equal to unity. 53 40 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). 41 The input parameters are the degrees of polymerization, the volume fractions, 42 the specific volumes, and the neutron scattering length densities for each 43 component. 61 44 62 45 … … 64 47 ---------- 65 48 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 49 A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 78 50 """ 79 51 … … 81 53 82 54 name = "rpa" 83 title = "Random Phase Approximation "55 title = "Random Phase Approximation - unfinished work in progress" 84 56 description = """ 85 57 This formalism applies to multicomponent polymer mixtures in the … … 118 90 ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 119 91 ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 120 ["v[4]", "mL/mol", 100.0, [0, inf], "", " molarvolume"],92 ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 121 93 ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 122 94 ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], … … 135 107 136 108 control = "case_num" 137 HIDE_ ALL = set("Phi4".split())138 HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()) .union(HIDE_ALL)109 HIDE_NONE = set() 110 HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()) 139 111 HIDE_AB = set("N2 Phi2 v2 L2 b2 K23 K24".split()).union(HIDE_A) 140 112 def hidden(case_num): … … 142 114 Return a list of parameters to hide depending on the multiplicity parameter. 143 115 """ 144 case_num = int(case_num+0.5)145 116 if case_num < 2: 146 117 return HIDE_AB … … 148 119 return HIDE_A 149 120 else: 150 return HIDE_ ALL121 return HIDE_NONE 151 122 -
sasmodels/models/spherical_sld.c
r925ad6e r925ad6e 1 1 static double form_volume( 2 double fp_n_shells,2 int n_shells, 3 3 double thickness[], 4 4 double interface[]) 5 5 { 6 int n_shells= (int)(fp_n_shells + 0.5);7 6 double r = 0.0; 8 7 for (int i=0; i < n_shells; i++) { … … 21 20 return pow(z, nu); 22 21 } else if (shape==2) { 23 return 1.0 - pow(1. 0- z, nu);22 return 1.0 - pow(1. - z, nu); 24 23 } else if (shape==3) { 25 24 return expm1(-nu*z)/expm1(-nu); … … 45 44 static double Iq( 46 45 double q, 47 double fp_n_shells,46 int n_shells, 48 47 double sld_solvent, 49 48 double sld[], … … 52 51 double shape[], 53 52 double nu[], 54 double fp_n_steps)53 int n_steps) 55 54 { 56 55 // 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);59 56 double f=0.0; 60 57 double r=0.0; -
sasmodels/models/spherical_sld.py
r925ad6e r925ad6e 233 233 """ 234 234 235 n_shells = int(n_shells + 0.5)236 n_steps = int(n_steps + 0.5)237 235 z = [] 238 236 rho = [] … … 242 240 rho.append(sld[0]) 243 241 244 for i in range(0, n_shells):242 for i in range(0, int(n_shells)): 245 243 z_next += thickness[i] 246 244 z.append(z_next) … … 263 261 def ER(n_shells, thickness, interface): 264 262 """Effective radius""" 265 n_shells = int(n_shells + 0.5)263 n_shells = int(n_shells) 266 264 total = (np.sum(thickness[:n_shells], axis=1) 267 265 + 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) 1 double form_volume(double thick_core, 2 double thick_layer, 3 double radius, 4 double n_stacking); 5 6 double Iq(double q, 7 double thick_core, 8 double thick_layer, 9 double radius, 10 double n_stacking, 11 double sigma_dnn, 12 double core_sld, 13 double layer_sld, 14 double solvent_sld); 15 16 double Iqxy(double qx, double qy, 17 double thick_core, 18 double thick_layer, 19 double radius, 20 double n_stacking, 21 double sigma_dnn, 22 double core_sld, 23 double layer_sld, 24 double solvent_sld, 25 double theta, 26 double phi); 27 28 static 29 double _kernel(double q, 30 double radius, 31 double core_sld, 32 double layer_sld, 33 double solvent_sld, 34 double halfheight, 35 double thick_layer, 36 double sin_alpha, 37 double cos_alpha, 38 double sigma_dnn, 39 double d, 40 double n_stacking) 14 41 15 42 { … … 61 88 62 89 63 static double stacked_disks_1d(64 65 double thick_core,66 double thick_layer,67 double radius,68 intn_stacking,69 double sigma_dnn,70 double core_sld,71 double layer_sld,72 double solvent_sld)90 static 91 double stacked_disks_kernel(double q, 92 double thick_core, 93 double thick_layer, 94 double radius, 95 double n_stacking, 96 double sigma_dnn, 97 double core_sld, 98 double layer_sld, 99 double solvent_sld) 73 100 { 74 101 /* StackedDiscsX : calculates the form factor of a stacked "tactoid" of core shell disks … … 84 111 double sin_alpha, cos_alpha; // slots to hold sincos function output 85 112 SINCOS(zi, sin_alpha, cos_alpha); 86 double yyy = stacked_disks_kernel(q, 87 halfheight, 88 thick_layer, 113 double yyy = _kernel(q, 89 114 radius, 90 n_stacking,91 sigma_dnn,92 115 core_sld, 93 116 layer_sld, 94 117 solvent_sld, 118 halfheight, 119 thick_layer, 95 120 sin_alpha, 96 121 cos_alpha, 97 d); 122 sigma_dnn, 123 d, 124 n_stacking); 98 125 summ += Gauss76Wt[i] * yyy * sin_alpha; 99 126 } … … 105 132 } 106 133 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); 134 double form_volume(double thick_core, 135 double thick_layer, 136 double radius, 137 double n_stacking){ 114 138 double d = 2.0 * thick_layer + thick_core; 115 139 return M_PI * radius * radius * d * n_stacking; 116 140 } 117 141 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, 142 double Iq(double q, 143 double thick_core, 144 double thick_layer, 145 double radius, 146 double n_stacking, 147 double sigma_dnn, 148 double core_sld, 149 double layer_sld, 150 double solvent_sld) 151 { 152 return stacked_disks_kernel(q, 131 153 thick_core, 132 154 thick_layer, … … 140 162 141 163 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); 164 double 165 Iqxy(double qx, double qy, 166 double thick_core, 167 double thick_layer, 168 double radius, 169 double n_stacking, 170 double sigma_dnn, 171 double core_sld, 172 double layer_sld, 173 double solvent_sld, 174 double theta, 175 double phi) 176 { 155 177 double q, sin_alpha, cos_alpha; 156 178 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); … … 158 180 double d = 2.0 * thick_layer + thick_core; 159 181 double halfheight = 0.5*thick_core; 160 double answer = stacked_disks_kernel(q, 161 halfheight, 162 thick_layer, 182 double answer = _kernel(q, 163 183 radius, 164 n_stacking,165 sigma_dnn,166 184 core_sld, 167 185 layer_sld, 168 186 solvent_sld, 187 halfheight, 188 thick_layer, 169 189 sin_alpha, 170 190 cos_alpha, 171 d); 191 sigma_dnn, 192 d, 193 n_stacking); 172 194 173 195 //convert to [cm-1] -
sasmodels/models/stacked_disks.py
rb7e8b94 rb7e8b94 126 126 ["thick_layer", "Ang", 10.0, [0, inf], "volume", "Thickness of layer each side of core"], 127 127 ["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"], 129 129 ["sigma_d", "Ang", 0, [0, inf], "", "Sigma of nearest neighbor spacing"], 130 130 ["sld_core", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Core scattering length density"], -
sasmodels/models/star_polymer.c
r2586093f r3a48772 3 3 double Iq(double q, double radius2, double arms); 4 4 5 static double star_polymer_kernel(double q, double radius2, double arms)5 static double _mass_fractal_kernel(double q, double radius2, double arms) 6 6 { 7 7 … … 23 23 double Iq(double q, double radius2, double arms) 24 24 { 25 return star_polymer_kernel(q, radius2, arms);25 return _mass_fractal_kernel(q, radius2, arms); 26 26 } -
sasmodels/models/unified_power_Rg.py
r66ca2a6 r66ca2a6 97 97 98 98 def 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: 101 101 with errstate(divide='ignore'): 102 102 return 1./q … … 104 104 with errstate(divide='ignore', invalid='ignore'): 105 105 result = np.zeros(q.shape, 'd') 106 for i in range( level):106 for i in range(ilevel): 107 107 exp_now = exp(-(q*rg[i])**2/3.) 108 108 pow_now = (erf(q*rg[i]/sqrt(6.))**3/q)**power[i] 109 if i < level-1:109 if i < ilevel-1: 110 110 exp_next = exp(-(q*rg[i+1])**2/3.) 111 111 else: … … 113 113 result += G[i]*exp_now + B[i]*exp_next*pow_now 114 114 115 result[q == 0] = np.sum(G[: level])115 result[q == 0] = np.sum(G[:ilevel]) 116 116 return result 117 117 -
sasmodels/product.py
rf88e248 r9951a86 45 45 # structure factor calculator. Structure factors should not 46 46 # 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 == []) 53 50 p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters 54 51 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:] 68 64 translate_name = dict((old.id, new.id) for old, new 69 65 in zip(s_pars.kernel_parameters[1:], s_list)) 70 66 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")) 72 69 demo.update((translate_name[k], v) for k, v in s_info.demo.items() 73 70 if k not in ("background", "scale") and not k.startswith(ER_ID)) … … 93 90 # Remember the component info blocks so we can build the model 94 91 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 = {} 102 93 return model_info 103 94 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 """ 95 def suffix_parameter(par, suffix): 111 96 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" 116 98 par.id = par.id + "_S" 117 par.name = par.id + vector_length118 99 return par 119 100
Note: See TracChangeset
for help on using the changeset viewer.