Changeset 797a8e3 in sasmodels
- Timestamp:
- Nov 26, 2016 6:28:36 PM (8 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 98ce141
- Parents:
- 5810f00 (diff), f80f334 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - Files:
-
- 1 added
- 77 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified README.rst ¶
r84bc3c1 re30d645 46 46 lamellar.py is an example of a single file model with embedded C code. 47 47 48 Magnetism hasn't been implemented yet. We may want a separate Imagnetic49 calculator with the extra parameters and calculations. We should50 return all desired spin states together so we can share the work of51 computing the form factors for the different magnetic contrasts. This52 will mean extending the data handler to support multiple cross sections53 in the same data set.54 55 48 |TravisStatus|_ 56 49 -
TabularUnified sasmodels/compare.py ¶
ra0d75ce r8c65a33 33 33 import datetime 34 34 import traceback 35 import re 35 36 36 37 import numpy as np # type: ignore … … 38 39 from . import core 39 40 from . import kerneldll 40 from . import weights41 from . import exception 41 42 from .data import plot_theory, empty_data1D, empty_data2D 42 43 from .direct_model import DirectModel 43 44 from .convert import revert_name, revert_pars, constrain_new_to_old 45 from .generate import FLOAT_RE 44 46 45 47 try: 46 48 from typing import Optional, Dict, Any, Callable, Tuple 47 except :49 except Exception: 48 50 pass 49 51 else: … … 58 60 sasmodels rewrite. 59 61 60 model is the name of the modelto compare (see below).62 model or model1,model2 are the names of the models to compare (see below). 61 63 N1 is the number of times to run sasmodels (default=1). 62 64 N2 is the number times to run sasview (default=1). … … 71 73 -preset*/-random[=seed] preset or random parameters 72 74 -mono/-poly* force monodisperse/polydisperse 75 -magnetic/-nonmagnetic* suppress magnetism 73 76 -cutoff=1e-5* cutoff value for including a point in polydispersity 74 77 -pars/-nopars* prints the parameter set or not … … 81 84 -default/-demo* use demo vs default parameters 82 85 -html shows the model docs instead of running the model 83 -title=" graph title" adds a title to the plot86 -title="note" adds note to the plot title, after the model name 84 87 85 88 Any two calculation engines can be selected for comparison: … … 89 92 -sasview sets the sasview calculation engine 90 93 91 The default is -single - sasview. Note that the interpretation of quad94 The default is -single -double. Note that the interpretation of quad 92 95 precision depends on architecture, and may vary from 64-bit to 128-bit, 93 96 with 80-bit floats being common (1e-19 precision). 94 97 95 98 Key=value pairs allow you to set specific values for the model parameters. 99 Key=value1,value2 to compare different values of the same parameter. 100 value can be an expression including other parameters 96 101 """ 97 102 … … 108 113 109 114 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True 115 116 # list of math functions for use in evaluating parameters 117 MATH = dict((k,getattr(math, k)) for k in dir(math) if not k.startswith('_')) 110 118 111 119 # CRUFT python 2.6 … … 314 322 name = name.split('*')[0] 315 323 316 if name == 'capped_cylinder' and pars['radius_cap'] < pars['radius']: 317 pars['radius'], pars['radius_cap'] = pars['radius_cap'], pars['radius'] 318 if name == 'barbell' and pars['radius_bell'] < pars['radius']: 319 pars['radius'], pars['radius_bell'] = pars['radius_bell'], pars['radius'] 320 321 # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff) 322 if name == 'guinier': 324 if name == 'barbell': 325 if pars['radius_bell'] < pars['radius']: 326 pars['radius'], pars['radius_bell'] = pars['radius_bell'], pars['radius'] 327 328 elif name == 'capped_cylinder': 329 if pars['radius_cap'] < pars['radius']: 330 pars['radius'], pars['radius_cap'] = pars['radius_cap'], pars['radius'] 331 332 elif name == 'guinier': 333 # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff) 323 334 #q_max = 0.2 # mid q maximum 324 335 q_max = 1.0 # high q maximum … … 326 337 pars['rg'] = min(pars['rg'], rg_max) 327 338 328 if name == 'rpa': 339 elif name == 'pearl_necklace': 340 if pars['radius'] < pars['thick_string']: 341 pars['radius'], pars['thick_string'] = pars['thick_string'], pars['radius'] 342 pars['num_pearls'] = math.ceil(pars['num_pearls']) 343 pass 344 345 elif name == 'rpa': 329 346 # Make sure phi sums to 1.0 330 347 if pars['case_num'] < 2: … … 337 354 pars['Phi'+c] /= total 338 355 356 elif name == 'stacked_disks': 357 pars['n_stacking'] = math.ceil(pars['n_stacking']) 358 339 359 def parlist(model_info, pars, is2d): 340 360 # type: (ModelInfo, ParameterSet, bool) -> str … … 344 364 lines = [] 345 365 parameters = model_info.parameters 366 magnetic = False 346 367 for p in parameters.user_parameters(pars, is2d): 368 if any(p.id.startswith(x) for x in ('M0:', 'mtheta:', 'mphi:')): 369 continue 370 if p.id.startswith('up:') and not magnetic: 371 continue 347 372 fields = dict( 348 373 value=pars.get(p.id, p.default), … … 352 377 pdtype=pars.get(p.id+"_pd_type", 'gaussian'), 353 378 relative_pd=p.relative_pd, 379 M0=pars.get('M0:'+p.id, 0.), 380 mphi=pars.get('mphi:'+p.id, 0.), 381 mtheta=pars.get('mtheta:'+p.id, 0.), 354 382 ) 355 383 lines.append(_format_par(p.name, **fields)) 384 magnetic = magnetic or fields['M0'] != 0. 356 385 return "\n".join(lines) 357 386 … … 359 388 360 389 def _format_par(name, value=0., pd=0., n=0, nsigma=3., pdtype='gaussian', 361 relative_pd=False ):390 relative_pd=False, M0=0., mphi=0., mtheta=0.): 362 391 # type: (str, float, float, int, float, str) -> str 363 392 line = "%s: %g"%(name, value) … … 367 396 line += " +/- %g (%d points in [-%g,%g] sigma %s)"\ 368 397 % (pd, n, nsigma, nsigma, pdtype) 398 if M0 != 0.: 399 line += " M0:%.3f mphi:%.1f mtheta:%.1f" % (M0, mphi, mtheta) 369 400 return line 370 401 … … 380 411 for p in pars: 381 412 if p.endswith("_pd_n"): pars[p] = 0 413 return pars 414 415 def suppress_magnetism(pars): 416 # type: (ParameterSet) -> ParameterSet 417 """ 418 Suppress theta_pd for now until the normalization is resolved. 419 420 May also suppress complete polydispersity of the model to test 421 models more quickly. 422 """ 423 pars = pars.copy() 424 for p in pars: 425 if p.startswith("M0:"): pars[p] = 0 382 426 return pars 383 427 … … 466 510 if k.endswith('.type'): 467 511 par = k[:-5] 512 if v == 'gaussian': continue 468 513 cls = dispersers[v if v != 'rectangle' else 'rectangula'] 469 514 handle = cls() 470 515 model[0].disperser_handles[par] = handle 471 model[0].set_dispersion(par, handle) 516 try: 517 model[0].set_dispersion(par, handle) 518 except Exception: 519 exception.annotate_exception("while setting %s to %r" 520 %(par, v)) 521 raise 522 472 523 473 524 #print("sasview pars",oldpars) … … 605 656 parameters. 606 657 """ 607 n_base, n_comp = opts['n1'], opts['n2'] 608 pars = opts['pars'] 658 result = run_models(opts, verbose=True) 659 if opts['plot']: # Note: never called from explore 660 plot_models(opts, result, limits=limits) 661 662 def run_models(opts, verbose=False): 663 # type: (Dict[str, Any]) -> Dict[str, Any] 664 665 n_base, n_comp = opts['count'] 666 pars, pars2 = opts['pars'] 609 667 data = opts['data'] 610 668 … … 612 670 base = opts['engines'][0] if n_base else None 613 671 comp = opts['engines'][1] if n_comp else None 672 614 673 base_time = comp_time = None 615 674 base_value = comp_value = resid = relerr = None … … 620 679 base_raw, base_time = time_calculation(base, pars, n_base) 621 680 base_value = np.ma.masked_invalid(base_raw) 622 print("%s t=%.2f ms, intensity=%.0f" 623 % (base.engine, base_time, base_value.sum())) 681 if verbose: 682 print("%s t=%.2f ms, intensity=%.0f" 683 % (base.engine, base_time, base_value.sum())) 624 684 _show_invalid(data, base_value) 625 685 except ImportError: … … 630 690 if n_comp > 0: 631 691 try: 632 comp_raw, comp_time = time_calculation(comp, pars , n_comp)692 comp_raw, comp_time = time_calculation(comp, pars2, n_comp) 633 693 comp_value = np.ma.masked_invalid(comp_raw) 634 print("%s t=%.2f ms, intensity=%.0f" 635 % (comp.engine, comp_time, comp_value.sum())) 694 if verbose: 695 print("%s t=%.2f ms, intensity=%.0f" 696 % (comp.engine, comp_time, comp_value.sum())) 636 697 _show_invalid(data, comp_value) 637 698 except ImportError: … … 643 704 resid = (base_value - comp_value) 644 705 relerr = resid/np.where(comp_value != 0., abs(comp_value), 1.0) 645 _print_stats("|%s-%s|" 646 % (base.engine, comp.engine) + (" "*(3+len(comp.engine))), 647 resid) 648 _print_stats("|(%s-%s)/%s|" 649 % (base.engine, comp.engine, comp.engine), 650 relerr) 706 if verbose: 707 _print_stats("|%s-%s|" 708 % (base.engine, comp.engine) + (" "*(3+len(comp.engine))), 709 resid) 710 _print_stats("|(%s-%s)/%s|" 711 % (base.engine, comp.engine, comp.engine), 712 relerr) 713 714 return dict(base_value=base_value, comp_value=comp_value, 715 base_time=base_time, comp_time=comp_time, 716 resid=resid, relerr=relerr) 717 718 719 def _print_stats(label, err): 720 # type: (str, np.ma.ndarray) -> None 721 # work with trimmed data, not the full set 722 sorted_err = np.sort(abs(err.compressed())) 723 if len(sorted_err) == 0.: 724 print(label + " no valid values") 725 return 726 727 p50 = int((len(sorted_err)-1)*0.50) 728 p98 = int((len(sorted_err)-1)*0.98) 729 data = [ 730 "max:%.3e"%sorted_err[-1], 731 "median:%.3e"%sorted_err[p50], 732 "98%%:%.3e"%sorted_err[p98], 733 "rms:%.3e"%np.sqrt(np.mean(sorted_err**2)), 734 "zero-offset:%+.3e"%np.mean(sorted_err), 735 ] 736 print(label+" "+" ".join(data)) 737 738 739 def plot_models(opts, result, limits=None): 740 # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float] 741 base_value, comp_value= result['base_value'], result['comp_value'] 742 base_time, comp_time = result['base_time'], result['comp_time'] 743 resid, relerr = result['resid'], result['relerr'] 744 745 have_base, have_comp = (base_value is not None), (comp_value is not None) 746 base = opts['engines'][0] if have_base else None 747 comp = opts['engines'][1] if have_comp else None 748 data = opts['data'] 651 749 652 750 # Plot if requested 653 if not opts['plot'] and not opts['explore']: return654 751 view = opts['view'] 655 752 import matplotlib.pyplot as plt 656 753 if limits is None: 657 754 vmin, vmax = np.Inf, -np.Inf 658 if n_base > 0:755 if have_base: 659 756 vmin = min(vmin, base_value.min()) 660 757 vmax = max(vmax, base_value.max()) 661 if n_comp > 0:758 if have_comp: 662 759 vmin = min(vmin, comp_value.min()) 663 760 vmax = max(vmax, comp_value.max()) 664 761 limits = vmin, vmax 665 762 666 if n_base > 0:667 if n_comp > 0: plt.subplot(131)763 if have_base: 764 if have_comp: plt.subplot(131) 668 765 plot_theory(data, base_value, view=view, use_data=False, limits=limits) 669 766 plt.title("%s t=%.2f ms"%(base.engine, base_time)) 670 767 #cbar_title = "log I" 671 if n_comp > 0: 672 if n_base > 0: plt.subplot(132) 768 if have_comp: 769 if have_base: plt.subplot(132) 770 if not opts['is2d'] and have_base: 771 plot_theory(data, base_value, view=view, use_data=False, limits=limits) 673 772 plot_theory(data, comp_value, view=view, use_data=False, limits=limits) 674 773 plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 675 774 #cbar_title = "log I" 676 if n_comp > 0 and n_base > 0:775 if have_base and have_comp: 677 776 plt.subplot(133) 678 777 if not opts['rel_err']: … … 680 779 else: 681 780 err, errstr, errview = abs(relerr), "rel err", "log" 781 if 0: # 95% cutoff 782 sorted = np.sort(err.flatten()) 783 cutoff = sorted[int(sorted.size*0.95)] 784 err[err>cutoff] = cutoff 682 785 #err,errstr = base/comp,"ratio" 683 786 plot_theory(data, None, resid=err, view=errview, use_data=False) … … 691 794 fig = plt.gcf() 692 795 extra_title = ' '+opts['title'] if opts['title'] else '' 693 fig.suptitle( opts['name']+ extra_title)694 695 if n_comp > 0 and n_base > 0 and '-hist' in opts:796 fig.suptitle(":".join(opts['name']) + extra_title) 797 798 if have_base and have_comp and opts['show_hist']: 696 799 plt.figure() 697 800 v = relerr … … 708 811 return limits 709 812 710 def _print_stats(label, err):711 # type: (str, np.ma.ndarray) -> None712 # work with trimmed data, not the full set713 sorted_err = np.sort(abs(err.compressed()))714 p50 = int((len(sorted_err)-1)*0.50)715 p98 = int((len(sorted_err)-1)*0.98)716 data = [717 "max:%.3e"%sorted_err[-1],718 "median:%.3e"%sorted_err[p50],719 "98%%:%.3e"%sorted_err[p98],720 "rms:%.3e"%np.sqrt(np.mean(sorted_err**2)),721 "zero-offset:%+.3e"%np.mean(sorted_err),722 ]723 print(label+" "+" ".join(data))724 813 725 814 … … 735 824 'preset', 'random', 736 825 'poly', 'mono', 826 'magnetic', 'nonmagnetic', 737 827 'nopars', 'pars', 738 828 'rel', 'abs', … … 791 881 return pars 792 882 793 883 INTEGER_RE = re.compile("^[+-]?[1-9][0-9]*$") 884 def isnumber(str): 885 match = FLOAT_RE.match(str) 886 isfloat = (match and not str[match.end():]) 887 return isfloat or INTEGER_RE.match(str) 888 889 # For distinguishing pairs of models for comparison 890 # key-value pair separator = 891 # shell characters | & ; <> $ % ' " \ # ` 892 # model and parameter names _ 893 # parameter expressions - + * / . ( ) 894 # path characters including tilde expansion and windows drive ~ / : 895 # not sure about brackets [] {} 896 # maybe one of the following @ ? ^ ! , 897 MODEL_SPLIT = ',' 794 898 def parse_opts(argv): 795 899 # type: (List[str]) -> Dict[str, Any] … … 813 917 print("expected parameters: model N1 N2") 814 918 815 name = positional_args[0]816 try:817 model_info = core.load_model_info(name)818 except ImportError as exc:819 print(str(exc))820 print("Could not find model; use one of:\n " + models)821 return None822 823 919 invalid = [o[1:] for o in flags 824 920 if o[1:] not in NAME_OPTIONS … … 828 924 return None 829 925 926 name = positional_args[0] 927 n1 = int(positional_args[1]) if len(positional_args) > 1 else 1 928 n2 = int(positional_args[2]) if len(positional_args) > 2 else 1 830 929 831 930 # pylint: disable=bad-whitespace … … 842 941 'seed' : -1, # default to preset 843 942 'mono' : False, 943 # Default to magnetic a magnetic moment is set on the command line 944 'magnetic' : False, 844 945 'show_pars' : False, 845 946 'show_hist' : False, … … 875 976 elif arg == '-mono': opts['mono'] = True 876 977 elif arg == '-poly': opts['mono'] = False 978 elif arg == '-magnetic': opts['magnetic'] = True 979 elif arg == '-nonmagnetic': opts['magnetic'] = False 877 980 elif arg == '-pars': opts['show_pars'] = True 878 981 elif arg == '-nopars': opts['show_pars'] = False … … 895 998 # pylint: enable=bad-whitespace 896 999 1000 if MODEL_SPLIT in name: 1001 name, name2 = name.split(MODEL_SPLIT, 2) 1002 else: 1003 name2 = name 1004 try: 1005 model_info = core.load_model_info(name) 1006 model_info2 = core.load_model_info(name2) if name2 != name else model_info 1007 except ImportError as exc: 1008 print(str(exc)) 1009 print("Could not find model; use one of:\n " + models) 1010 return None 1011 1012 # Get demo parameters from model definition, or use default parameters 1013 # if model does not define demo parameters 1014 pars = get_pars(model_info, opts['use_demo']) 1015 pars2 = get_pars(model_info2, opts['use_demo']) 1016 pars2.update((k, v) for k, v in pars.items() if k in pars2) 1017 # randomize parameters 1018 #pars.update(set_pars) # set value before random to control range 1019 if opts['seed'] > -1: 1020 pars = randomize_pars(model_info, pars, seed=opts['seed']) 1021 if model_info != model_info2: 1022 pars2 = randomize_pars(model_info2, pars2, seed=opts['seed']) 1023 # Share values for parameters with the same name 1024 for k, v in pars.items(): 1025 if k in pars2: 1026 pars2[k] = v 1027 else: 1028 pars2 = pars.copy() 1029 constrain_pars(model_info, pars) 1030 constrain_pars(model_info2, pars2) 1031 print("Randomize using -random=%i"%opts['seed']) 1032 if opts['mono']: 1033 pars = suppress_pd(pars) 1034 pars2 = suppress_pd(pars2) 1035 if not opts['magnetic']: 1036 pars = suppress_magnetism(pars) 1037 pars2 = suppress_magnetism(pars2) 1038 1039 # Fill in parameters given on the command line 1040 presets = {} 1041 presets2 = {} 1042 for arg in values: 1043 k, v = arg.split('=', 1) 1044 if k not in pars and k not in pars2: 1045 # extract base name without polydispersity info 1046 s = set(p.split('_pd')[0] for p in pars) 1047 print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 1048 return None 1049 v1, v2 = v.split(MODEL_SPLIT, 2) if MODEL_SPLIT in v else (v,v) 1050 if v1 and k in pars: 1051 presets[k] = float(v1) if isnumber(v1) else v1 1052 if v2 and k in pars2: 1053 presets2[k] = float(v2) if isnumber(v2) else v2 1054 1055 # If pd given on the command line, default pd_n to 35 1056 for k, v in list(presets.items()): 1057 if k.endswith('_pd'): 1058 presets.setdefault(k+'_n', 35.) 1059 for k, v in list(presets2.items()): 1060 if k.endswith('_pd'): 1061 presets2.setdefault(k+'_n', 35.) 1062 1063 # Evaluate preset parameter expressions 1064 context = MATH.copy() 1065 context.update(pars) 1066 context.update((k,v) for k,v in presets.items() if isinstance(v, float)) 1067 for k, v in presets.items(): 1068 if not isinstance(v, float) and not k.endswith('_type'): 1069 presets[k] = eval(v, context) 1070 context.update(presets) 1071 context.update((k,v) for k,v in presets2.items() if isinstance(v, float)) 1072 for k, v in presets2.items(): 1073 if not isinstance(v, float) and not k.endswith('_type'): 1074 presets2[k] = eval(v, context) 1075 1076 # update parameters with presets 1077 pars.update(presets) # set value after random to control value 1078 pars2.update(presets2) # set value after random to control value 1079 #import pprint; pprint.pprint(model_info) 1080 1081 same_model = name == name2 and pars == pars 897 1082 if len(engines) == 0: 898 engines.extend(['single', 'double']) 1083 if same_model: 1084 engines.extend(['single', 'double']) 1085 else: 1086 engines.extend(['single', 'single']) 899 1087 elif len(engines) == 1: 900 if engines[0][0] == 'double': 1088 if not same_model: 1089 engines.append(engines[0]) 1090 elif engines[0] == 'double': 901 1091 engines.append('single') 902 1092 else: … … 905 1095 del engines[2:] 906 1096 907 n1 = int(positional_args[1]) if len(positional_args) > 1 else 1908 n2 = int(positional_args[2]) if len(positional_args) > 2 else 1909 1097 use_sasview = any(engine == 'sasview' and count > 0 910 1098 for engine, count in zip(engines, [n1, n2])) 911 912 # Get demo parameters from model definition, or use default parameters913 # if model does not define demo parameters914 pars = get_pars(model_info, opts['use_demo'])915 916 917 # Fill in parameters given on the command line918 presets = {}919 for arg in values:920 k, v = arg.split('=', 1)921 if k not in pars:922 # extract base name without polydispersity info923 s = set(p.split('_pd')[0] for p in pars)924 print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))925 return None926 presets[k] = float(v) if not k.endswith('type') else v927 928 # randomize parameters929 #pars.update(set_pars) # set value before random to control range930 if opts['seed'] > -1:931 pars = randomize_pars(model_info, pars, seed=opts['seed'])932 print("Randomize using -random=%i"%opts['seed'])933 if opts['mono']:934 pars = suppress_pd(pars)935 pars.update(presets) # set value after random to control value936 #import pprint; pprint.pprint(model_info)937 constrain_pars(model_info, pars)938 1099 if use_sasview: 939 1100 constrain_new_to_old(model_info, pars) 1101 constrain_new_to_old(model_info2, pars2) 1102 940 1103 if opts['show_pars']: 941 print(str(parlist(model_info, pars, opts['is2d']))) 1104 if not same_model: 1105 print("==== %s ====="%model_info.name) 1106 print(str(parlist(model_info, pars, opts['is2d']))) 1107 print("==== %s ====="%model_info2.name) 1108 print(str(parlist(model_info2, pars2, opts['is2d']))) 1109 else: 1110 print(str(parlist(model_info, pars, opts['is2d']))) 942 1111 943 1112 # Create the computational engines … … 948 1117 base = None 949 1118 if n2: 950 comp = make_engine(model_info , data, engines[1], opts['cutoff'])1119 comp = make_engine(model_info2, data, engines[1], opts['cutoff']) 951 1120 else: 952 1121 comp = None … … 955 1124 # Remember it all 956 1125 opts.update({ 957 'name' : name,958 'def' : model_info,959 'n1' : n1,960 'n2' : n2,961 'presets' : presets,962 'pars' : pars,963 1126 'data' : data, 1127 'name' : [name, name2], 1128 'def' : [model_info, model_info2], 1129 'count' : [n1, n2], 1130 'presets' : [presets, presets2], 1131 'pars' : [pars, pars2], 964 1132 'engines' : [base, comp], 965 1133 }) … … 976 1144 from .generate import view_html_from_info 977 1145 app = wx.App() if wx.GetApp() is None else None 978 view_html_from_info(opts['def'] )1146 view_html_from_info(opts['def'][0]) 979 1147 if app: app.MainLoop() 980 1148 … … 988 1156 from bumps.names import FitProblem # type: ignore 989 1157 from bumps.gui.app_frame import AppFrame # type: ignore 1158 from bumps.gui import signal 990 1159 991 1160 is_mac = "cocoa" in wx.version() 992 1161 # Create an app if not running embedded 993 1162 app = wx.App() if wx.GetApp() is None else None 994 problem = FitProblem(Explore(opts)) 1163 model = Explore(opts) 1164 problem = FitProblem(model) 995 1165 frame = AppFrame(parent=None, title="explore", size=(1000,700)) 996 1166 if not is_mac: frame.Show() … … 998 1168 frame.panel.Layout() 999 1169 frame.panel.aui.Split(0, wx.TOP) 1170 def reset_parameters(event): 1171 model.revert_values() 1172 signal.update_parameters(problem) 1173 frame.Bind(wx.EVT_TOOL, reset_parameters, frame.ToolBar.GetToolByPos(1)) 1000 1174 if is_mac: frame.Show() 1001 1175 # If running withing an app, start the main loop … … 1015 1189 config_matplotlib() 1016 1190 self.opts = opts 1017 model_info = opts['def'] 1018 pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 1191 p1, p2 = opts['pars'] 1192 m1, m2 = opts['def'] 1193 self.fix_p2 = m1 != m2 or p1 != p2 1194 model_info = m1 1195 pars, pd_types = bumps_model.create_parameters(model_info, **p1) 1019 1196 # Initialize parameter ranges, fixing the 2D parameters for 1D data. 1020 1197 if not opts['is2d']: 1021 for p in model_info.parameters.user_parameters( is2d=False):1198 for p in model_info.parameters.user_parameters({}, is2d=False): 1022 1199 for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 1023 1200 k = p.name+ext … … 1030 1207 1031 1208 self.pars = pars 1209 self.starting_values = dict((k, v.value) for k, v in pars.items()) 1032 1210 self.pd_types = pd_types 1033 1211 self.limits = None 1212 1213 def revert_values(self): 1214 for k, v in self.starting_values.items(): 1215 self.pars[k].value = v 1216 1217 def model_update(self): 1218 pass 1034 1219 1035 1220 def numpoints(self): … … 1062 1247 pars = dict((k, v.value) for k, v in self.pars.items()) 1063 1248 pars.update(self.pd_types) 1064 self.opts['pars'] = pars 1065 limits = compare(self.opts, limits=self.limits) 1249 self.opts['pars'][0] = pars 1250 if not self.fix_p2: 1251 self.opts['pars'][1] = pars 1252 result = run_models(self.opts) 1253 limits = plot_models(self.opts, result, limits=self.limits) 1066 1254 if self.limits is None: 1067 1255 vmin, vmax = limits 1068 1256 self.limits = vmax*1e-7, 1.3*vmax 1257 import pylab; pylab.clf() 1258 plot_models(self.opts, result, limits=self.limits) 1069 1259 1070 1260 -
TabularUnified sasmodels/conversion_table.py ¶
r5f1acda rfa6d6fc 143 143 } 144 144 ], 145 "core_shell_ellipsoid ": [145 "core_shell_ellipsoid:1": [ 146 146 "CoreShellEllipsoidModel", 147 147 { 148 "sld_core": "sld_core", 149 "sld_shell": "sld_shell", 150 "sld_solvent": "sld_solvent", 151 "radius_equat_core": "equat_core", 152 "x_core": "polar_core", 153 "thick_shell": "equat_shell", 154 "x_polar_shell": "polar_shell", 155 "theta": "axis_theta", 148 156 "phi": "axis_phi", 149 "sld_core": "sld_core",150 "polar_shell": "polar_shell",151 "sld_solvent": "sld_solvent",152 "equat_shell": "equat_shell",153 "equat_core": "equat_core",154 "theta": "axis_theta",155 "polar_core": "polar_core",156 "sld_shell": "sld_shell"157 157 } 158 158 ], … … 160 160 "CoreShellEllipsoidXTModel", 161 161 { 162 "phi": "axis_phi",163 162 "sld_core": "sld_core", 163 "sld_shell": "sld_shell", 164 "sld_solvent": "sld_solvent", 165 "radius_equat_core": "equat_core", 166 "thick_shell": "T_shell", 164 167 "x_core": "X_core", 165 "sld_solvent": "sld_solvent",166 "thick_shell": "T_shell",167 168 "x_polar_shell": "XpolarShell", 168 169 "theta": "axis_theta", 169 " sld_shell": "sld_shell"170 "phi": "axis_phi", 170 171 } 171 172 ], … … 173 174 "CSParallelepipedModel", 174 175 { 176 "sld_core": "sld_pcore", 177 "sld_a": "sld_rimA", 178 "sld_b": "sld_rimB", 179 "sld_c": "sld_rimC", 180 "sld_solvent": "sld_solv", 181 "length_a": "shortA", 182 "length_b": "midB", 183 "length_c": "longC", 184 "thick_rim_a": "rimA", 185 "thick_rim_c": "rimC", 186 "thick_rim_b": "rimB", 187 "theta": "parallel_theta", 175 188 "phi": "parallel_phi", 176 189 "psi": "parallel_psi", 177 "sld_core": "sld_pcore",178 "sld_c": "sld_rimC",179 "sld_b": "sld_rimB",180 "sld_solvent": "sld_solv",181 "length_a": "shortA",182 "sld_a": "sld_rimA",183 "length_b": "midB",184 "thick_rimc": "rimC",185 "theta": "parallel_theta",186 "thick_rim_a": "rimA",187 "length_c": "longC",188 "thick_rim_b": "rimB"189 190 } 190 191 ], … … 240 241 "DABModel", 241 242 { 242 " length": "length"243 "cor_length": "length" 243 244 } 244 245 ], … … 257 258 "EllipticalCylinderModel", 258 259 { 260 "axis_ratio": "r_ratio", 261 "radius_minor": "r_minor", 262 "sld": "sldCyl", 263 "sld_solvent": "sldSolv", 264 "theta": "cyl_theta", 259 265 "phi": "cyl_phi", 260 266 "psi": "cyl_psi", 261 "theta": "cyl_theta",262 "sld": "sldCyl",263 "axis_ratio": "r_ratio",264 "sld_solvent": "sldSolv"265 267 } 266 268 ], … … 299 301 "FractalCoreShellModel", 300 302 { 303 "sld_core": "core_sld", 301 304 "sld_shell": "shell_sld", 302 305 "sld_solvent": "solvent_sld", 303 "sld_core": "core_sld" 306 "radius": "radius", 307 "thickness": "thickness", 308 "fractal_dim": "frac_dim", 309 "cor_length": "cor_length", 310 "volfraction": "volfraction", 304 311 } 305 312 ], … … 326 333 "PeakGaussModel", 327 334 { 328 "sigma": "B" 335 "peak_pos": "q0", 336 "sigma": "B", 329 337 } 330 338 ], … … 334 342 "rg": "radius", 335 343 "lorentz_scale": "lScale", 336 "fractal_dim": "FractalExp", 344 "guinier_scale": "gScale", 345 "fractal_dim": "scale", 337 346 "cor_length": "zeta", 338 "guinier_scale": "gScale"339 347 } 340 348 ], … … 350 358 "s": "dim", 351 359 "rg": "rg", 352 " m": "m",360 "porod_exp": "m", 353 361 "scale": "scale", 354 362 "background": "background" … … 358 366 "HardsphereStructure", 359 367 { 360 " radius_effective_pd": "effect_radius_pd",368 "scale": "scale_factor", 361 369 "radius_effective": "effect_radius", 362 "radius_effective_pd_n": "effect_radius_pd_n"363 370 } 364 371 ], … … 366 373 "HayterMSAStructure", 367 374 { 368 "salt_concentration": "saltconc", 369 "radius_effective_pd": "effect_radius_pd", 375 "scale": "scale_factor", 370 376 "radius_effective": "effect_radius", 371 "radius_effective_pd_n": "effect_radius_pd_n" 377 "volfraction": "volfraction", 378 "charge": "charge", 379 "temperature": "temperature", 380 "concentration_salt": "saltconc", 381 "dielectconst": "dielectconst", 372 382 } 373 383 ], … … 375 385 "HollowCylinderModel", 376 386 { 387 "sld": "sldCyl", 388 "sld_solvent": "sldSolv", 389 "radius": "core_radius", 390 "thickness": "radius", 391 "length": "length", 392 "theta": "axis_theta", 377 393 "phi": "axis_phi", 378 "scale": "scale",379 "radius_core": "core_radius",380 "sld_solvent": "sldSolv",381 "length": "length",382 "radius": "radius",383 "background": "background",384 "sld": "sldCyl",385 "theta": "axis_theta"386 394 } 387 395 ], … … 389 397 "RectangularHollowPrismModel", 390 398 { 399 "sld": "sldPipe", 400 "sld_solvent": "sldSolv", 401 "length_a": "short_side", 391 402 "b2a_ratio": "b2a_ratio", 392 "length_a": "short_side", 393 "sld": "sldPipe", 394 "length_c": "c2a_ratio", 395 "sld_solvent": "sldSolv", 396 "thickness": "thickness" 403 "c2a_ratio": "c2a_ratio", 404 "thickness": "thickness", 397 405 } 398 406 ], … … 401 409 { 402 410 "sld": "sldPipe", 411 "sld_solvent": "sldSolv", 412 "length_a": "short_side", 403 413 "b2a_ratio": "b2a_ratio", 404 "length_a": "short_side", 405 "length_c": "c2a_ratio", 406 "sld_solvent": "sldSolv" 414 "c2a_ratio": "c2a_ratio", 407 415 } 408 416 ], … … 428 436 "LamellarPSHGModel", 429 437 { 438 "sld": "sld_tail", 439 "sld_head": "sld_head", 440 "sld_solvent": "sld_solvent", 441 "length_tail": "deltaT", 442 "length_head": "deltaH", 443 "d_spacing": "spacing", 430 444 "Caille_parameter": "caille", 431 445 "Nlayers": "n_plates", 432 "sld_head": "sld_head",433 "length_tail": "deltaT",434 "length_head": "deltaH",435 "sld": "sld_tail",436 "sld_solvent": "sld_solvent"437 446 } 438 447 ], … … 441 450 { 442 451 "sld": "sld_bi", 452 "sld_solvent": "sld_sol", 453 "thickness": "delta", 454 "d_spacing": "spacing", 443 455 "Caille_parameter": "caille", 444 456 "Nlayers": "N_plates", 445 "sld_solvent": "sld_sol",446 "thickness": "delta"447 457 } 448 458 ], … … 451 461 { 452 462 "sld": "sld_layer", 463 "sld_solvent": "sld_solvent", 464 "thickness": "thickness", 465 "d_spacing": "spacing", 453 466 "sigma_d": "pd_spacing", 454 " sld_solvent": "sld_solvent"467 "Nlayers": "Nlayers", 455 468 } 456 469 ], … … 500 513 { 501 514 "rg": "rg", 502 " scale": "scale",503 "background": "background" 515 "i_zero": "scale", 516 "background": "background", 504 517 } 505 518 ], … … 576 589 "rg": "rg", 577 590 "polydispersity": "poly_m", 591 "i_zero": "scale", 592 "background": "background", 593 } 594 ], 595 "polymer_excl_volume": [ 596 "PolymerExclVolume", 597 { 598 "rg": "rg", 599 "scale": "scale", 600 "background": "background", 601 "porod_exp": "m" 602 } 603 ], 604 "polymer_micelle": [ 605 "MicelleSphCoreModel", 606 { 607 "sld_corona": "rho_corona", 608 "sld_solvent": "rho_solv", 609 "sld_core": "rho_core", 610 "ndensity": "ndensity", 611 "v_core": "v_core", 612 "v_corona": "v_corona", 613 "radius_core": "radius_core", 614 "rg": "radius_gyr", 615 "d_penetration": "d_penetration", 616 "n_aggreg": "n_aggreg", 617 } 618 ], 619 "porod": [ 620 "PorodModel", 621 { 578 622 "scale": "scale", 579 623 "background": "background" 580 624 } 581 625 ], 582 "polymer_excl_volume": [583 "PolymerExclVolume",584 {585 "rg": "rg",586 "scale": "scale",587 "background": "background",588 "porod_exp": "m"589 }590 ],591 "polymer_micelle": [592 "MicelleSphCoreModel",593 {594 "sld_corona": "rho_corona",595 "sld_solvent": "rho_solv",596 "sld_core": "rho_core"597 }598 ],599 "porod": [600 "PorodModel",601 {602 "scale": "scale",603 "background": "background"604 }605 ],606 626 "power_law": [ 607 627 "PowerLawAbsModel", … … 616 636 { 617 637 "scale": "scale", 618 "s olvent_sld": "sld_solvent",638 "sld_solvent": "sld_solvent", 619 639 "thickness": "thickness", 620 640 "beta": "beta", … … 643 663 { 644 664 "sld": "sldPipe", 665 "length_a": "short_side", 645 666 "b2a_ratio": "b2a_ratio", 646 "length_a": "short_side", 647 "length_c": "c2a_ratio", 667 "c2a_ratio": "c2a_ratio", 648 668 "sld_solvent": "sldSolv" 649 669 } … … 716 736 "SquareWellStructure", 717 737 { 718 " radius_effective_pd": "effect_radius_pd",738 "scale": "scale_factor", 719 739 "radius_effective": "effect_radius", 720 "radius_effective_pd_n": "effect_radius_pd_n" 740 "wellwidth": "wellwidth", 741 "welldepth": "welldepth", 721 742 } 722 743 ], … … 729 750 "theta": "axis_theta", 730 751 "sld_solvent": "solvent_sld", 731 "n_stacking": "n_stacking" 752 "n_stacking": "n_stacking", 753 "thick_layer": "layer_thick", 754 "thick_core": "core_thick", 732 755 } 733 756 ], … … 742 765 "StickyHSStructure", 743 766 { 744 " radius_effective_pd": "effect_radius_pd",767 "scale": "scale_factor", 745 768 "radius_effective": "effect_radius", 746 "radius_effective_pd_n": "effect_radius_pd_n"747 769 } 748 770 ], … … 758 780 "TeubnerStreyModel", 759 781 { 760 "a2": "scale" 782 # Note: parameters are completely rewritten in convert.py 783 "volfraction_a": "volfraction_a", 784 "sld_a": "sld_a", 785 "sld_b": "sld_b", 786 "d": "d", 787 "xi": "xi", 761 788 } 762 789 ], -
TabularUnified sasmodels/convert.py ¶
r51241113 rf80f334 2 2 Convert models to and from sasview. 3 3 """ 4 from __future__ import print_function 5 6 from os.path import join as joinpath, abspath, dirname4 from __future__ import print_function, division 5 6 import re 7 7 import math 8 8 import warnings 9 9 10 10 from .conversion_table import CONVERSION_TABLE 11 from .core import load_model_info 11 12 12 13 # List of models which SasView versions don't contain the explicit 'scale' argument. … … 17 18 'two_lorentzian', 18 19 "two_power_law", 19 'gel_fit',20 20 'gauss_lorentz_gel', 21 21 'be_polyelectrolyte', … … 49 49 # Convert new style names for polydispersity info to old style names 50 50 PD_DOT = [ 51 ("", ""),52 51 ("_pd", ".width"), 53 52 ("_pd_n", ".npts"), … … 56 55 ] 57 56 58 def _convert_pars(pars, mapping): 59 """ 60 Rename the parameters and any associated polydispersity attributes. 61 """ 62 newpars = pars.copy() 63 for new, old in mapping.items(): 64 if old == new: continue 65 for underscore, dot in PD_DOT: 66 if old+dot in newpars: 67 if new is not None: 68 newpars[new+underscore] = pars[old+dot] 69 del newpars[old+dot] 70 return newpars 71 72 def convert_model(name, pars): 73 """ 74 Convert model from old style parameter names to new style. 75 """ 76 _, _ = name, pars # lint 77 raise NotImplementedError 78 # need to load all new models in order to determine old=>new 79 # model name mapping 80 81 def _unscale(par, scale): 57 def _rescale(par, scale): 82 58 return [pk*scale for pk in par] if isinstance(par, list) else par*scale 83 59 84 def _is_sld(modelinfo, id): 60 def _is_sld(model_info, id): 61 """ 62 Return True if parameter is a magnetic magnitude or SLD parameter. 63 """ 85 64 if id.startswith('M0:'): 86 65 return True 87 if (id.endswith('_pd') or id.endswith('_pd_n') or id.endswith('_pd_nsigma') 88 or id.endswith('_pd_width') or id.endswith('_pd_type')): 66 if id.startswith('volfraction') or id.startswith('radius_effective'): 89 67 return False 90 for p in modelinfo.parameters.call_parameters: 68 if '_pd' in id or '.' in id: 69 return False 70 for p in model_info.parameters.call_parameters: 91 71 if p.id == id: 92 72 return p.type == 'sld' 93 73 # check through kernel parameters in case it is a named as a vector 94 for p in model info.parameters.kernel_parameters:74 for p in model_info.parameters.kernel_parameters: 95 75 if p.id == id: 96 76 return p.type == 'sld' 97 77 raise ValueError("unknown parameter %r in conversion"%id) 98 78 99 def _ unscale_sld(modelinfo, pars):100 """ 101 rescale all sld parameters in the new model definition by 1e6so the79 def _rescale_sld(model_info, pars, scale): 80 """ 81 rescale all sld parameters in the new model definition by *scale* so the 102 82 numbers are nicer. Relies on the fact that all sld parameters in the 103 new model definition end with sld. 104 """ 105 return dict((id, (_unscale(v, 1e-6) if _is_sld(modelinfo, id) else v)) 83 new model definition end with sld. For backward conversion use 84 *scale=1e-6*. For forward conversion use *scale=1e6*. 85 """ 86 return dict((id, (_rescale(v, scale) if _is_sld(model_info, id) else v)) 106 87 for id, v in pars.items()) 107 88 108 def _remove_pd(pars, key, name):109 """110 Remove polydispersity from the parameter list.111 112 Note: operates in place113 """114 # Bumps style parameter names115 width = pars.pop(key+".width", 0.0)116 n_points = pars.pop(key+".npts", 0)117 if width != 0.0 and n_points != 0:118 warnings.warn("parameter %s not polydisperse in sasview %s"%(key, name))119 pars.pop(key+".nsigmas", None)120 pars.pop(key+".type", None)121 return pars122 123 def _revert_pars(pars, mapping):124 """125 Rename the parameters and any associated polydispersity attributes.126 """127 newpars = pars.copy()128 129 for new, old in mapping.items():130 for underscore, dot in PD_DOT:131 if old and old+underscore == new+dot:132 continue133 if new+underscore in newpars:134 if old is not None:135 newpars[old+dot] = pars[new+underscore]136 del newpars[new+underscore]137 for k in list(newpars.keys()):138 for underscore, dot in PD_DOT[1:]: # skip "" => ""139 if k.endswith(underscore):140 newpars[k[:-len(underscore)]+dot] = newpars[k]141 del newpars[k]142 return newpars143 144 def revert_name(model_info):145 oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}])146 return oldname147 89 148 90 def _get_translation_table(model_info): … … 162 104 return translation 163 105 106 # ========= FORWARD CONVERSION sasview 3.x => sasmodels =========== 107 def _dot_pd_to_underscore_pd(par): 108 if par.endswith(".width"): 109 return par[:-6]+"_pd" 110 elif par.endswith(".type"): 111 return par[:-5]+"_pd_type" 112 elif par.endswith(".nsigmas"): 113 return par[:-8]+"_pd_nsigma" 114 elif par.endswith(".npts"): 115 return par[:-5]+"_pd_n" 116 else: 117 return par 118 119 def _pd_to_underscores(pars): 120 return dict((_dot_pd_to_underscore_pd(k), v) for k, v in pars.items()) 121 122 def _convert_name(conv_dict, pars): 123 """ 124 Renames parameter values (upper, lower, etc) to v4.0 names 125 :param conv_dict: conversion dictionary mapping new name : old name 126 :param pars: parameters to convert 127 :return: 128 """ 129 new_pars = {} 130 i = 0 131 j = 0 132 for key_par, value_par in pars.iteritems(): 133 j += 1 134 for key_conv, value_conv in conv_dict.iteritems(): 135 if re.search(value_conv, key_par): 136 new_pars[key_par.replace(value_conv, key_conv)] = value_par 137 i += 1 138 break 139 elif re.search("background", key_par) or re.search("scale", key_par): 140 new_pars[key_par] = value_par 141 i += 1 142 break 143 if i != j: 144 new_pars[key_par] = value_par 145 i += 1 146 return new_pars 147 148 def _convert_pars(pars, mapping): 149 """ 150 Rename the parameters and any associated polydispersity attributes. 151 """ 152 newpars = pars.copy() 153 for new, old in mapping.items(): 154 if old == new: continue 155 if old is None: continue 156 for underscore, dot in PD_DOT: 157 source = old+dot 158 if source in newpars: 159 if new is not None: 160 target = new+dot 161 else: 162 target = None 163 if source != target: 164 if target: 165 newpars[target] = pars[old+dot] 166 del newpars[source] 167 return newpars 168 169 170 def _conversion_target(model_name): 171 """ 172 Find the sasmodel name which translates into the sasview name. 173 174 Note: *CoreShellEllipsoidModel* translates into *core_shell_ellipsoid:1*. 175 This is necessary since there is only one variant in sasmodels for the 176 two variants in sasview. 177 """ 178 for sasmodels_name, [sasview_name, _] in CONVERSION_TABLE.items(): 179 if sasview_name == model_name: 180 return sasmodels_name 181 return None 182 183 184 def _hand_convert(name, oldpars): 185 if name == 'core_shell_parallelepiped': 186 # Make sure pd on rim parameters defaults to zero 187 # ... probably not necessary. 188 oldpars['rimA.width'] = 0.0 189 oldpars['rimB.width'] = 0.0 190 oldpars['rimC.width'] = 0.0 191 elif name == 'core_shell_ellipsoid:1': 192 # Reverse translation (from new to old), from core_shell_ellipsoid.c 193 # equat_shell = equat_core + thick_shell 194 # polar_core = equat_core * x_core 195 # polar_shell = equat_core * x_core + thick_shell*x_polar_shell 196 # Forward translation (from old to new), inverting reverse translation: 197 # thick_shell = equat_shell - equat_core 198 # x_core = polar_core / equat_core 199 # x_polar_shell = (polar_shell - polar_core)/(equat_shell - equat_core) 200 # Auto translation (old <=> new) happens after hand_convert 201 # equat_shell <=> thick_shell 202 # polar_core <=> x_core 203 # polar_shell <=> x_polar_shell 204 # So... 205 equat_core, equat_shell = oldpars['equat_core'], oldpars['equat_shell'] 206 polar_core, polar_shell = oldpars['polar_core'], oldpars['polar_shell'] 207 oldpars['equat_shell'] = equat_shell - equat_core 208 oldpars['polar_core'] = polar_core / equat_core 209 oldpars['polar_shell'] = (polar_shell-polar_core)/(equat_shell-equat_core) 210 elif name == 'hollow_cylinder': 211 # now uses radius and thickness 212 thickness = oldpars['radius'] - oldpars['core_radius'] 213 oldpars['radius'] = thickness 214 if 'radius.width' in oldpars: 215 pd = oldpars['radius.width']*oldpars['radius']/thickness 216 oldpars['radius.width'] = pd 217 elif name == 'pearl_necklace': 218 pass 219 #_remove_pd(oldpars, 'num_pearls', name) 220 #_remove_pd(oldpars, 'thick_string', name) 221 elif name == 'polymer_micelle': 222 if 'ndensity' in oldpars: 223 oldpars['ndensity'] /= 1e15 224 if 'ndensity.lower' in oldpars: 225 oldpars['ndensity.lower'] /= 1e15 226 if 'ndensity.upper' in oldpars: 227 oldpars['ndensity.upper'] /= 1e15 228 elif name == 'rpa': 229 # convert scattering lengths from femtometers to centimeters 230 for p in "L1", "L2", "L3", "L4": 231 if p in oldpars: 232 oldpars[p] /= 1e-13 233 if p + ".lower" in oldpars: 234 oldpars[p + ".lower"] /= 1e-13 235 if p + ".upper" in oldpars: 236 oldpars[p + ".upper"] /= 1e-13 237 elif name == 'spherical_sld': 238 oldpars["CONTROL"] += 1 239 elif name == 'teubner_strey': 240 # basically undoing the entire Teubner-Strey calculations here. 241 # drho = (sld_a - sld_b) 242 # k = 2.0*math.pi*xi/d 243 # a2 = (1.0 + k**2)**2 244 # c1 = 2.0 * xi**2 * (1.0 - k**2) 245 # c2 = xi**4 246 # prefactor = 8.0*math.pi*phi*(1.0-phi)*drho**2*c2/xi 247 # scale = 1e-4*prefactor 248 # oldpars['scale'] = a2/scale 249 # oldpars['c1'] = c1/scale 250 # oldpars['c2'] = c2/scale 251 252 # need xi, d, sld_a, sld_b, phi=volfraction_a 253 # assume contrast is 1.0e-6, scale=1, background=0 254 sld_a, sld_b = 1.0, 0. 255 drho = sld_a - sld_b 256 257 # find xi 258 p_scale = oldpars['scale'] 259 p_c1 = oldpars['c1'] 260 p_c2= oldpars['c2'] 261 i_1 = 0.5*p_c1/p_c2 262 i_2 = math.sqrt(math.fabs(p_scale/p_c2)) 263 i_3 = 2/(i_1 + i_2) 264 xi = math.sqrt(math.fabs(i_3)) 265 266 # find d from xi 267 k = math.sqrt(math.fabs(1 - 0.5*p_c1/p_c2*xi**2)) 268 d = 2*math.pi*xi/k 269 270 # solve quadratic phi (1-phi) = xi/(1e-4 8 pi drho^2 c2) 271 # favour volume fraction in [0, 0.5] 272 c = xi / (1e-4 * 8.0 * math.pi * drho**2 * p_c2) 273 phi = 0.5 - math.sqrt(0.25 - c) 274 275 # scale sld_a by 1e-6 because the translator will scale it back 276 oldpars.update(volfraction_a=phi, xi=xi, d=d, sld_a=sld_a*1e-6, 277 sld_b=sld_b, scale=1.0) 278 oldpars.pop('c1') 279 oldpars.pop('c2') 280 281 return oldpars 282 283 def convert_model(name, pars, use_underscore=False): 284 """ 285 Convert model from old style parameter names to new style. 286 """ 287 newname = _conversion_target(name) 288 if newname is None: 289 return name, pars 290 if ':' in newname: # core_shell_ellipsoid:1 291 model_info = load_model_info(newname[:-2]) 292 # Know that the table exists and isn't multiplicity so grab it directly 293 # Can't use _get_translation_table since that will return the 'bare' 294 # version. 295 translation = CONVERSION_TABLE[newname][1] 296 else: 297 model_info = load_model_info(newname) 298 translation = _get_translation_table(model_info) 299 newpars = _hand_convert(newname, pars.copy()) 300 newpars = _convert_name(translation, newpars) 301 newpars = _convert_pars(newpars, translation) 302 if not model_info.structure_factor: 303 newpars = _rescale_sld(model_info, newpars, 1e6) 304 newpars.setdefault('scale', 1.0) 305 newpars.setdefault('background', 0.0) 306 if use_underscore: 307 newpars = _pd_to_underscores(newpars) 308 return newname, newpars 309 310 311 # ========= BACKWARD CONVERSION sasmodels => sasview 3.x =========== 312 313 def _revert_pars(pars, mapping): 314 """ 315 Rename the parameters and any associated polydispersity attributes. 316 """ 317 newpars = pars.copy() 318 319 for new, old in mapping.items(): 320 for underscore, dot in PD_DOT: 321 if old and old+underscore == new+dot: 322 continue 323 if new+underscore in newpars: 324 if old is not None: 325 newpars[old+dot] = pars[new+underscore] 326 del newpars[new+underscore] 327 for k in list(newpars.keys()): 328 for underscore, dot in PD_DOT[1:]: # skip "" => "" 329 if k.endswith(underscore): 330 newpars[k[:-len(underscore)]+dot] = newpars[k] 331 del newpars[k] 332 return newpars 333 334 def revert_name(model_info): 335 oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}]) 336 return oldname 337 338 def _remove_pd(pars, key, name): 339 """ 340 Remove polydispersity from the parameter list. 341 342 Note: operates in place 343 """ 344 # Bumps style parameter names 345 width = pars.pop(key+".width", 0.0) 346 n_points = pars.pop(key+".npts", 0) 347 if width != 0.0 and n_points != 0: 348 warnings.warn("parameter %s not polydisperse in sasview %s"%(key, name)) 349 pars.pop(key+".nsigmas", None) 350 pars.pop(key+".type", None) 351 return pars 352 164 353 def _trim_vectors(model_info, pars, oldpars): 165 354 _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}]) … … 180 369 composition_type, parts = model_info.composition 181 370 if composition_type == 'product': 182 translation = {'scale':'scale_factor'}183 translation.update(_get_translation_table(parts[0]))371 translation = _get_translation_table(parts[0]) 372 # structure factor models include scale:scale_factor mapping 184 373 translation.update(_get_translation_table(parts[1])) 185 374 else: … … 187 376 else: 188 377 translation = _get_translation_table(model_info) 189 oldpars = _revert_pars(_ unscale_sld(model_info, pars), translation)378 oldpars = _revert_pars(_rescale_sld(model_info, pars, 1e-6), translation) 190 379 oldpars = _trim_vectors(model_info, pars, oldpars) 191 380 … … 215 404 if name in MODELS_WITHOUT_VOLFRACTION: 216 405 del oldpars['volfraction'] 217 if name == 'stacked_disks':218 _remove_pd(oldpars, 'n_stacking', name)219 elif name == 'pearl_necklace':220 _remove_pd(oldpars, 'num_pearls', name)221 _remove_pd(oldpars, 'thick_string', name)222 elif name == 'core_shell_parallelepiped':223 _remove_pd(oldpars, 'rimA', name)224 _remove_pd(oldpars, 'rimB', name)225 _remove_pd(oldpars, 'rimC', name)226 elif name == 'polymer_micelle':227 if 'ndensity' in oldpars:228 oldpars['ndensity'] *= 1e15229 elif name == 'spherical_sld':230 oldpars["CONTROL"] -= 1231 # remove polydispersity from shells232 for k in range(1, 11):233 _remove_pd(oldpars, 'thick_flat'+str(k), 'thickness')234 _remove_pd(oldpars, 'thick_inter'+str(k), 'interface')235 # remove extra shells236 for k in range(int(pars['n_shells']), 11):237 oldpars.pop('sld_flat'+str(k), 0)238 oldpars.pop('thick_flat'+str(k), 0)239 oldpars.pop('thick_inter'+str(k), 0)240 oldpars.pop('func_inter'+str(k), 0)241 oldpars.pop('nu_inter'+str(k), 0)242 406 elif name == 'core_multi_shell': 243 407 # kill extra shells … … 252 416 elif name == 'core_shell_parallelepiped': 253 417 _remove_pd(oldpars, 'rimA', name) 254 elif name in ['mono_gauss_coil', 'poly_gauss_coil']: 255 del oldpars['i_zero'] 418 _remove_pd(oldpars, 'rimB', name) 419 _remove_pd(oldpars, 'rimC', name) 420 elif name == 'hollow_cylinder': 421 # now uses radius and thickness 422 thickness = oldpars['core_radius'] 423 oldpars['radius'] += thickness 424 oldpars['radius.width'] *= thickness/oldpars['radius'] 425 #elif name in ['mono_gauss_coil', 'poly_gauss_coil']: 426 # del oldpars['i_zero'] 256 427 elif name == 'onion': 257 428 oldpars.pop('n_shells', None) 429 elif name == 'pearl_necklace': 430 _remove_pd(oldpars, 'num_pearls', name) 431 _remove_pd(oldpars, 'thick_string', name) 432 elif name == 'polymer_micelle': 433 if 'ndensity' in oldpars: 434 oldpars['ndensity'] *= 1e15 258 435 elif name == 'rpa': 259 436 # convert scattering lengths from femtometers to centimeters … … 272 449 for k in "Kab,Kac,Kad".split(','): 273 450 oldpars.pop(k, None) 451 elif name == 'spherical_sld': 452 oldpars["CONTROL"] -= 1 453 # remove polydispersity from shells 454 for k in range(1, 11): 455 _remove_pd(oldpars, 'thick_flat'+str(k), 'thickness') 456 _remove_pd(oldpars, 'thick_inter'+str(k), 'interface') 457 # remove extra shells 458 for k in range(int(pars['n_shells']), 11): 459 oldpars.pop('sld_flat'+str(k), 0) 460 oldpars.pop('thick_flat'+str(k), 0) 461 oldpars.pop('thick_inter'+str(k), 0) 462 oldpars.pop('func_inter'+str(k), 0) 463 oldpars.pop('nu_inter'+str(k), 0) 464 elif name == 'stacked_disks': 465 _remove_pd(oldpars, 'n_stacking', name) 466 elif name == 'teubner_strey': 467 # basically redoing the entire Teubner-Strey calculations here. 468 volfraction = oldpars.pop('volfraction_a') 469 xi = oldpars.pop('xi') 470 d = oldpars.pop('d') 471 sld_a = oldpars.pop('sld_a') 472 sld_b = oldpars.pop('sld_b') 473 drho = 1e6*(sld_a - sld_b) # conversion autoscaled these 474 k = 2.0*math.pi*xi/d 475 a2 = (1.0 + k**2)**2 476 c1 = 2.0 * xi**2 * (1.0 - k**2) 477 c2 = xi**4 478 prefactor = 8.0*math.pi*volfraction*(1.0-volfraction)*drho**2*c2/xi 479 scale = 1e-4*prefactor 480 oldpars['scale'] = a2/scale 481 oldpars['c1'] = c1/scale 482 oldpars['c2'] = c2/scale 274 483 275 484 #print("convert from",list(sorted(pars))) … … 315 524 if name in MODELS_WITHOUT_VOLFRACTION: 316 525 pars['volfraction'] = 1 317 if name == 'pearl_necklace': 318 pars['string_thickness_pd_n'] = 0 319 pars['number_of_pearls_pd_n'] = 0 526 if name == 'core_multi_shell': 527 pars['n'] = min(math.ceil(pars['n']), 4) 528 elif name == 'gel_fit': 529 pars['scale'] = 1 320 530 elif name == 'line': 321 531 pars['scale'] = 1 322 532 pars['background'] = 0 533 elif name == 'mono_gauss_coil': 534 pars['scale'] = 1 535 elif name == 'onion': 536 pars['n_shells'] = math.ceil(pars['n_shells']) 537 elif name == 'pearl_necklace': 538 pars['string_thickness_pd_n'] = 0 539 pars['number_of_pearls_pd_n'] = 0 540 elif name == 'poly_gauss_coil': 541 pars['scale'] = 1 323 542 elif name == 'rpa': 324 543 pars['case_num'] = int(pars['case_num']) 325 elif name == 'mono_gauss_coil':326 pars['i_zero'] = 1327 elif name == 'poly_gauss_coil':328 pars['i_zero'] = 1329 elif name == 'core_multi_shell':330 pars['n'] = min(math.ceil(pars['n']), 4)331 elif name == 'onion':332 pars['n_shells'] = math.ceil(pars['n_shells'])333 544 elif name == 'spherical_sld': 334 545 pars['n_shells'] = math.ceil(pars['n_shells']) … … 339 550 pars['thickness%d_pd_n'%k] = 0 340 551 pars['interface%d_pd_n'%k] = 0 341 552 elif name == 'teubner_strey': 553 pars['scale'] = 1 554 if pars['volfraction_a'] > 0.5: 555 pars['volfraction_a'] = 1.0 - pars['volfraction_a'] 556 elif name == 'unified_power_Rg': 557 pars['level'] = int(pars['level']) 558 559 def _check_one(name, seed=None): 560 """ 561 Generate a random set of parameters for *name*, and check that they can 562 be converted back to SasView 3.x and forward again to sasmodels. Raises 563 an error if the parameters are changed. 564 """ 565 from . import compare 566 567 model_info = load_model_info(name) 568 569 old_name = revert_name(model_info) 570 if old_name is None: 571 return 572 573 pars = compare.get_pars(model_info, use_demo=False) 574 pars = compare.randomize_pars(model_info, pars, seed=seed) 575 if name == "teubner_strey": 576 # T-S model is underconstrained, so fix the assumptions. 577 pars['sld_a'], pars['sld_b'] = 1.0, 0.0 578 compare.constrain_pars(model_info, pars) 579 constrain_new_to_old(model_info, pars) 580 old_pars = revert_pars(model_info, pars) 581 new_name, new_pars = convert_model(old_name, old_pars, use_underscore=True) 582 if 1: 583 print("==== %s in ====="%name) 584 print(str(compare.parlist(model_info, pars, True))) 585 print("==== %s ====="%old_name) 586 for k, v in sorted(old_pars.items()): 587 print(k, v) 588 print("==== %s out ====="%new_name) 589 print(str(compare.parlist(model_info, new_pars, True))) 590 assert name==new_name, "%r != %r"%(name, new_name) 591 for k, v in new_pars.items(): 592 assert k in pars, "%s: %r appeared from conversion"%(name, k) 593 if isinstance(v, float): 594 assert abs(v-pars[k])<=abs(1e-12*v), "%s: %r %s != %s"%(name, k, v, pars[k]) 595 else: 596 assert v == pars[k], "%s: %r %s != %s"%(name, k, v, pars[k]) 597 for k, v in pars.items(): 598 assert k in pars, "%s: %r not converted"%(name, k) 599 600 def test_backward_forward(): 601 from .core import list_models 602 for name in list_models('all'): 603 L = lambda: _check_one(name, seed=1) 604 L.description = name 605 yield L -
TabularUnified sasmodels/core.py ¶
r1875f4e r52e9a45 53 53 "nonmagnetic", "magnetic") 54 54 def list_models(kind=None): 55 # type: ( ) -> List[str]55 # type: (str) -> List[str] 56 56 """ 57 57 Return the list of available models on the model path. -
TabularUnified sasmodels/generate.py ¶
r234c532 r7fcdc9f 166 166 import re 167 167 import string 168 from zlib import crc32 168 169 169 170 import numpy as np # type: ignore … … 356 357 return newest 357 358 359 def tag_source(source): 360 # type: (str) -> str 361 """ 362 Return a unique tag for the source code. 363 """ 364 # Note: need 0xffffffff&val to force an unsigned 32-bit number 365 return "%08X"%(0xffffffff&crc32(source)) 358 366 359 367 def convert_type(source, dtype): -
TabularUnified sasmodels/kernel_header.c ¶
re1d6983 r218cdbc 148 148 inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 149 149 150 #if 0 151 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 152 SINCOS(phi*M_PI_180, sn, cn); \ 153 q = sqrt(qx*qx + qy*qy); \ 154 cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * cos(theta*M_PI_180)); \ 155 sn = sqrt(1 - cn*cn); \ 156 } while (0) 157 #else 158 // SasView 3.x definition of orientation 159 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 160 SINCOS(theta*M_PI_180, sn, cn); \ 161 q = sqrt(qx*qx + qy*qy);\ 162 cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \ 163 sn = sqrt(1 - cn*cn); \ 164 } while (0) 165 #endif 166 167 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 168 q = sqrt(qx*qx + qy*qy); \ 169 const double qxhat = qx/q; \ 170 const double qyhat = qy/q; \ 171 double sin_theta, cos_theta; \ 172 double sin_phi, cos_phi; \ 173 double sin_psi, cos_psi; \ 174 SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 175 SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 176 SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 177 cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ 178 cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ 179 cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 180 } while (0) -
TabularUnified sasmodels/kernelcl.py ¶
r6e5b2a7 r7fcdc9f 294 294 # so we don't really need to cache things for ourselves. I'll do so 295 295 # anyway just to save some data munging time. 296 key = "%s-%s%s"%(name, dtype, ("-fast" if fast else "")) 296 tag = generate.tag_source(source) 297 key = "%s-%s-%s%s"%(name, dtype, tag, ("-fast" if fast else "")) 297 298 # Check timestamp on program 298 299 program, program_timestamp = self.compiled.get(key, (None, np.inf)) -
TabularUnified sasmodels/modelinfo.py ¶
rca55aa0 r85fe7f8 548 548 return full_list 549 549 550 def user_parameters(self, pars ={}, is2d=True):550 def user_parameters(self, pars, is2d=True): 551 551 # type: (Dict[str, float], bool) -> List[Parameter] 552 552 """ -
TabularUnified sasmodels/models/barbell.c ¶
r0d6e865 r3a48772 39 39 // translate dx in [-1,1] to dx in [lower,upper] 40 40 const double integral = total*zm; 41 const double bell_ Fq = 2*M_PI*cube(radius_bell)*integral;42 return bell_ Fq;41 const double bell_fq = 2.0*M_PI*cube(radius_bell)*integral; 42 return bell_fq; 43 43 } 44 45 static double 46 _fq(double q, double h, 47 double radius_bell, double radius, double half_length, 48 double sin_alpha, double cos_alpha) 49 { 50 const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 51 const double bj = sas_J1c(q*radius*sin_alpha); 52 const double si = sinc(q*half_length*cos_alpha); 53 const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 54 const double Aq = bell_fq + cyl_fq; 55 return Aq; 56 } 57 44 58 45 59 double form_volume(double radius_bell, … … 47 61 double length) 48 62 { 49 50 63 // bell radius should never be less than radius when this is called 51 64 const double hdist = sqrt(square(radius_bell) - square(radius)); … … 71 84 double sin_alpha, cos_alpha; // slots to hold sincos function output 72 85 SINCOS(alpha, sin_alpha, cos_alpha); 73 74 const double bell_Fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 75 const double bj = sas_J1c(q*radius*sin_alpha); 76 const double si = sinc(q*half_length*cos_alpha); 77 const double cyl_Fq = M_PI*radius*radius*length*bj*si; 78 const double Aq = bell_Fq + cyl_Fq; 86 const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 79 87 total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 80 88 } … … 93 101 double theta, double phi) 94 102 { 95 // Compute angle alpha between q and the cylinder axis 96 double sn, cn; // slots to hold sincos function output 97 SINCOS(phi*M_PI_180, sn, cn); 98 const double q = sqrt(qx*qx + qy*qy); 99 const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 100 const double alpha = acos(cos_val); // rod angle relative to q 103 double q, sin_alpha, cos_alpha; 104 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 101 105 102 106 const double h = -sqrt(square(radius_bell) - square(radius)); 103 const double half_length = 0.5*length; 104 105 double sin_alpha, cos_alpha; // slots to hold sincos function output 106 SINCOS(alpha, sin_alpha, cos_alpha); 107 const double bell_Fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 108 const double bj = sas_J1c(q*radius*sin_alpha); 109 const double si = sinc(q*half_length*cos_alpha); 110 const double cyl_Fq = M_PI*radius*radius*length*bj*si; 111 const double Aq = cyl_Fq + bell_Fq; 107 const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 112 108 113 109 // Multiply by contrast^2 and convert to cm-1 -
TabularUnified sasmodels/models/bcc_paracrystal.c ¶
r0bef47b r4962519 31 31 const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 32 32 const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 33 33 34 const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 34 result = pow(1.0-(temp10*temp10),3)*temp635 result = cube(1.0 - (temp10*temp10))*temp6 35 36 / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 36 37 * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) … … 81 82 82 83 return answer; 83 84 85 84 } 86 85 87 86 88 double Iqxy(double qx, double qy, double dnn, 89 double d_factor, double radius,double sld, double solvent_sld, 90 double theta, double phi, double psi){ 87 double Iqxy(double qx, double qy, 88 double dnn, double d_factor, double radius, 89 double sld, double solvent_sld, 90 double theta, double phi, double psi) 91 { 92 double q, cos_a1, cos_a2, cos_a3; 93 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 91 94 92 double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z, 93 //double q_z; 94 double cos_val_b3, cos_val_b2, cos_val_b1; 95 double a1_dot_q, a2_dot_q,a3_dot_q; 96 double answer; 97 double Zq, Fkq, Fkq_2; 95 const double a1 = +cos_a3 - cos_a1 + cos_a2; 96 const double a2 = +cos_a3 + cos_a1 - cos_a2; 97 const double a3 = -cos_a3 + cos_a1 + cos_a2; 98 98 99 //convert to q and make scaled values 100 double q = sqrt(qx*qx+qy*qy); 101 double q_x = qx/q; 102 double q_y = qy/q; 99 const double qd = 0.5*q*dnn; 100 const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 101 const double tanh_qd = tanh(arg); 102 const double cosh_qd = cosh(arg); 103 const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 104 * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 105 * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 103 106 104 //convert angle degree to radian 105 theta = theta * M_PI_180; 106 phi = phi * M_PI_180; 107 psi = psi * M_PI_180; 108 109 const double Da = d_factor*dnn; 110 const double s1 = dnn/sqrt(0.75); 111 112 113 //the occupied volume of the lattice 114 const double latticescale = 2.0*sphere_volume(radius/s1); 115 // q vector 116 //q_z = 0.0; // for SANS; assuming qz is negligible 117 /// Angles here are respect to detector coordinate 118 /// instead of against q coordinate(PRB 36(46), 3(6), 1754(3854)) 119 // b3 axis orientation 120 b3_x = cos(theta) * cos(phi); 121 b3_y = sin(theta); 122 //b3_z = -cos(theta) * sin(phi); 123 cos_val_b3 = b3_x*q_x + b3_y*q_y;// + b3_z*q_z; 124 125 //alpha = acos(cos_val_b3); 126 // b1 axis orientation 127 b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 128 b1_y = sin(psi)*cos(theta); 129 cos_val_b1 = b1_x*q_x + b1_y*q_y; 130 // b2 axis orientation 131 b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 132 b2_y = cos(theta)*cos(psi); 133 cos_val_b2 = b2_x*q_x + b2_y*q_y; 134 135 // The following test should always pass 136 if (fabs(cos_val_b3)>1.0) { 137 //printf("bcc_ana_2D: Unexpected error: cos()>1\n"); 138 cos_val_b3 = 1.0; 139 } 140 if (fabs(cos_val_b2)>1.0) { 141 //printf("bcc_ana_2D: Unexpected error: cos()>1\n"); 142 cos_val_b2 = 1.0; 143 } 144 if (fabs(cos_val_b1)>1.0) { 145 //printf("bcc_ana_2D: Unexpected error: cos()>1\n"); 146 cos_val_b1 = 1.0; 147 } 148 // Compute the angle btw vector q and the a3 axis 149 a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3); 150 151 // a1 axis 152 a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1); 153 154 // a2 axis 155 a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2); 156 157 158 // Get Fkq and Fkq_2 159 Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q)); 160 Fkq_2 = Fkq*Fkq; 161 // Call Zq=Z1*Z2*Z3 162 Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2); 163 Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2); 164 Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2); 165 166 // Use SphereForm directly from libigor 167 answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale; 168 169 return answer; 170 } 107 const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 108 //the occupied volume of the lattice 109 const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 110 return lattice_scale * Fq; 111 } -
TabularUnified sasmodels/models/binary_hard_sphere.c ¶
re481a39 r4f79d94 7 7 ); 8 8 9 double Iqxy(double qx, double qy,10 double lg_radius, double sm_radius,11 double lg_vol_frac, double sm_vol_frac,12 double lg_sld, double sm_sld, double solvent_sld13 );14 15 9 void calculate_psfs(double qval, 16 10 double r2, double nf2, … … 55 49 // /* do form factor calculations */ 56 50 57 v1 = 4.0*M_PI/3.0*r1*r1*r1;58 v2 = 4.0*M_PI/3.0*r2*r2*r2;51 v1 = M_4PI_3*r1*r1*r1; 52 v2 = M_4PI_3*r2*r2*r2; 59 53 60 54 n1 = phi1/v1; … … 64 58 qr2 = r2*q; 65 59 66 //if (qr1 == 0){67 //sc1 = 1.0/3.0;68 //}else{69 //sc1 = (sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;70 //}71 //if (qr2 == 0){72 //sc2 = 1.0/3.0;73 //}else{74 //sc2 = (sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;75 //}76 60 sc1 = sph_j1c(qr1); 77 61 sc2 = sph_j1c(qr2); 78 b1 = r1*r1*r1*(rho1-rhos)* 4.0/3.0*M_PI*sc1;79 b2 = r2*r2*r2*(rho2-rhos)* 4.0/3.0*M_PI*sc2;62 b1 = r1*r1*r1*(rho1-rhos)*M_4PI_3*sc1; 63 b2 = r2*r2*r2*(rho2-rhos)*M_4PI_3*sc2; 80 64 inten = n1*b1*b1*psf11; 81 65 inten += sqrt(n1*n2)*2.0*b1*b2*psf12; … … 88 72 } 89 73 90 91 double Iqxy(double qx, double qy,92 double lg_radius, double sm_radius,93 double lg_vol_frac, double sm_vol_frac,94 double lg_sld, double sm_sld, double solvent_sld)95 96 {97 double q = sqrt(qx*qx + qy*qy);98 return Iq(q,99 lg_radius, sm_radius,100 lg_vol_frac, sm_vol_frac,101 lg_sld, sm_sld, solvent_sld);102 }103 74 104 75 void calculate_psfs(double qval, -
TabularUnified sasmodels/models/capped_cylinder.c ¶
r0d6e865 r3a48772 45 45 // translate dx in [-1,1] to dx in [lower,upper] 46 46 const double integral = total*zm; 47 const double cap_Fq = 2 *M_PI*cube(radius_cap)*integral;47 const double cap_Fq = 2.0*M_PI*cube(radius_cap)*integral; 48 48 return cap_Fq; 49 } 50 51 static double 52 _fq(double q, double h, double radius_cap, double radius, double half_length, 53 double sin_alpha, double cos_alpha) 54 { 55 const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 56 const double bj = sas_J1c(q*radius*sin_alpha); 57 const double si = sinc(q*half_length*cos_alpha); 58 const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 59 const double Aq = cap_Fq + cyl_Fq; 60 return Aq; 49 61 } 50 62 … … 93 105 SINCOS(alpha, sin_alpha, cos_alpha); 94 106 95 const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 96 const double bj = sas_J1c(q*radius*sin_alpha); 97 const double si = sinc(q*half_length*cos_alpha); 98 const double cyl_Fq = M_PI*radius*radius*length*bj*si; 99 const double Aq = cap_Fq + cyl_Fq; 100 total += Gauss76Wt[i] * Aq * Aq * sin_alpha; // sin_alpha for spherical coord integration 107 const double Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha); 108 // sin_alpha for spherical coord integration 109 total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 101 110 } 102 111 // translate dx in [-1,1] to dx in [lower,upper] … … 114 123 double theta, double phi) 115 124 { 116 // Compute angle alpha between q and the cylinder axis 117 double sn, cn; 118 SINCOS(phi*M_PI_180, sn, cn); 119 const double q = sqrt(qx*qx+qy*qy); 120 const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 121 const double alpha = acos(cos_val); // rod angle relative to q 125 double q, sin_alpha, cos_alpha; 126 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 122 127 123 128 const double h = sqrt(radius_cap*radius_cap - radius*radius); 124 const double half_length = 0.5*length; 125 126 double sin_alpha, cos_alpha; // slots to hold sincos function output 127 SINCOS(alpha, sin_alpha, cos_alpha); 128 const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 129 const double bj = sas_J1c(q*radius*sin_alpha); 130 const double si = sinc(q*half_length*cos_alpha); 131 const double cyl_Fq = M_PI*radius*radius*length*bj*si; 132 const double Aq = cap_Fq + cyl_Fq; 129 const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 133 130 134 131 // Multiply by contrast^2 and convert to cm-1 -
TabularUnified sasmodels/models/core_shell_bicelle.c ¶
r0d6e865 r5bddd89 26 26 double form_volume(double radius, double thick_rim, double thick_face, double length) 27 27 { 28 return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2 *thick_face);28 return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face); 29 29 } 30 30 … … 39 39 double rhor, 40 40 double rhosolv, 41 double dum) 41 double sin_alpha, 42 double cos_alpha) 42 43 { 43 44 double si1,si2,be1,be2; … … 49 50 const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick); 50 51 const double vol3 = M_PI*rad*rad*2.0*(length+facthick); 51 double sn,cn; 52 SINCOS(dum, sn, cn); 53 double besarg1 = qq*rad*sn; 54 double besarg2 = qq*(rad+radthick)*sn; 55 double sinarg1 = qq*length*cn; 56 double sinarg2 = qq*(length+facthick)*cn; 52 double besarg1 = qq*rad*sin_alpha; 53 double besarg2 = qq*(rad+radthick)*sin_alpha; 54 double sinarg1 = qq*length*cos_alpha; 55 double sinarg2 = qq*(length+facthick)*cos_alpha; 57 56 58 57 be1 = sas_J1c(besarg1); … … 65 64 vol3*dr3*si2*be1; 66 65 67 const double retval = t*t*s n;66 const double retval = t*t*sin_alpha; 68 67 69 return (retval);68 return retval; 70 69 71 70 } … … 83 82 { 84 83 // set up the integration end points 85 const double uplim = M_PI /4;86 const double halfheight = length/2.0;84 const double uplim = M_PI_4; 85 const double halfheight = 0.5*length; 87 86 88 87 double summ = 0.0; 89 88 for(int i=0;i<N_POINTS_76;i++) { 90 double zi = (Gauss76Z[i] + 1.0)*uplim; 89 double alpha = (Gauss76Z[i] + 1.0)*uplim; 90 double sin_alpha, cos_alpha; // slots to hold sincos function output 91 SINCOS(alpha, sin_alpha, cos_alpha); 91 92 double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 92 halfheight, rhoc, rhoh, rhor,rhosolv, zi); 93 halfheight, rhoc, rhoh, rhor, rhosolv, 94 sin_alpha, cos_alpha); 93 95 summ += yyy; 94 96 } … … 96 98 // calculate value of integral to return 97 99 double answer = uplim*summ; 98 return (answer);100 return answer; 99 101 } 100 102 101 103 static double 102 bicelle_kernel_2d(double q , double q_x, double q_y,104 bicelle_kernel_2d(double qx, double qy, 103 105 double radius, 104 106 double thick_rim, … … 112 114 double phi) 113 115 { 114 //convert angle degree to radian 115 theta *= M_PI_180; 116 phi *= M_PI_180; 116 double q, sin_alpha, cos_alpha; 117 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 117 118 118 // Cylinder orientation119 const double cyl_x = sin(theta) * cos(phi);120 const double cyl_y = sin(theta) * sin(phi);121 122 // Compute the angle btw vector q and the axis of the cylinder123 const double cos_val = cyl_x*q_x + cyl_y*q_y;124 const double alpha = acos( cos_val );125 126 // Get the kernel127 119 double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 128 length/2.0, core_sld, face_sld, rim_sld,129 solvent_sld, alpha) / fabs(sin(alpha));120 0.5*length, core_sld, face_sld, rim_sld, 121 solvent_sld, sin_alpha, cos_alpha) / fabs(sin_alpha); 130 122 131 123 answer *= 1.0e-4; … … 162 154 double phi) 163 155 { 164 double q; 165 q = sqrt(qx*qx+qy*qy); 166 double intensity = bicelle_kernel_2d(q, qx/q, qy/q, 156 double intensity = bicelle_kernel_2d(qx, qy, 167 157 radius, 168 158 thick_rim, -
TabularUnified sasmodels/models/core_shell_cylinder.c ¶
r0d6e865 r9aa4881 5 5 double radius, double thickness, double length, double theta, double phi); 6 6 7 // twovd = 2 *volume * delta_rho7 // vd = volume * delta_rho 8 8 // besarg = q * R * sin(alpha) 9 9 // siarg = q * L/2 * cos(alpha) 10 double _cyl(double twovd, double besarg, double siarg);11 double _cyl(double twovd, double besarg, double siarg)10 double _cyl(double vd, double besarg, double siarg); 11 double _cyl(double vd, double besarg, double siarg) 12 12 { 13 const double bj = (besarg == 0.0 ? 0.5 : 0.5*sas_J1c(besarg)); 14 const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 15 return twovd*si*bj; 13 return vd * sinc(siarg) * sas_J1c(besarg); 16 14 } 17 15 18 16 double form_volume(double radius, double thickness, double length) 19 17 { 20 return M_PI*(radius+thickness)*(radius+thickness)*(length+2 *thickness);18 return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 21 19 } 22 20 … … 32 30 const double core_qr = q*radius; 33 31 const double core_qh = q*0.5*length; 34 const double core_twovd = 2.0 * form_volume(radius,0,length) 35 * (core_sld-shell_sld); 32 const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 36 33 const double shell_qr = q*(radius + thickness); 37 34 const double shell_qh = q*(0.5*length + thickness); 38 const double shell_twovd = 2.0 * form_volume(radius,thickness,length) 39 * (shell_sld-solvent_sld); 35 const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 40 36 double total = 0.0; 41 37 // double lower=0, upper=M_PI_2; … … 46 42 const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 47 43 SINCOS(alpha, sn, cn); 48 const double fq = _cyl(core_ twovd, core_qr*sn, core_qh*cn)49 + _cyl(shell_ twovd, shell_qr*sn, shell_qh*cn);44 const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn) 45 + _cyl(shell_vd, shell_qr*sn, shell_qh*cn); 50 46 total += Gauss76Wt[i] * fq * fq * sn; 51 47 } … … 66 62 double phi) 67 63 { 68 double sn, cn; // slots to hold sincos function output 69 70 // Compute angle alpha between q and the cylinder axis 71 SINCOS(phi*M_PI_180, sn, cn); 72 // # The following correction factor exists in sasview, but it can't be 73 // # right, so we are leaving it out for now. 74 // const double correction = fabs(cn)*M_PI_2; 75 const double q = sqrt(qx*qx+qy*qy); 76 const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 77 const double alpha = acos(cos_val); 64 double q, sin_alpha, cos_alpha; 65 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 78 66 79 67 const double core_qr = q*radius; 80 68 const double core_qh = q*0.5*length; 81 const double core_twovd = 2.0 * form_volume(radius,0,length) 82 * (core_sld-shell_sld); 69 const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 83 70 const double shell_qr = q*(radius + thickness); 84 71 const double shell_qh = q*(0.5*length + thickness); 85 const double shell_twovd = 2.0 * form_volume(radius,thickness,length) 86 * (shell_sld-solvent_sld); 72 const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 87 73 88 SINCOS(alpha, sn, cn); 89 const double fq = _cyl(core_twovd, core_qr*sn, core_qh*cn) 90 + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn); 74 const double fq = _cyl(core_vd, core_qr*sin_alpha, core_qh*cos_alpha) 75 + _cyl(shell_vd, shell_qr*sin_alpha, shell_qh*cos_alpha); 91 76 return 1.0e-4 * fq * fq; 92 77 } -
TabularUnified sasmodels/models/core_shell_ellipsoid.c ¶
r0d6e865 r0a3d9b2 32 32 const double equat_shell = radius_equat_core + thick_shell; 33 33 const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 34 double vol = 4.0*M_PI/3.0*equat_shell*equat_shell*polar_shell;34 double vol = M_4PI_3*equat_shell*equat_shell*polar_shell; 35 35 return vol; 36 36 } … … 49 49 const double uplim = 1.0; 50 50 51 double summ = 0.0; //initialize intergral52 51 53 52 const double delpc = core_sld - shell_sld; //core - shell … … 59 58 const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 60 59 61 for(int i=0;i<N_POINTS_76;i++) { 62 double zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 63 double yyy = Gauss76Wt[i] * gfn4(zi, 64 radius_equat_core, 65 polar_core, 66 equat_shell, 67 polar_shell, 68 delpc, 69 delps, 70 q); 71 summ += yyy; 60 double summ = 0.0; //initialize intergral 61 for(int i=0;i<76;i++) { 62 double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim ); 63 double yyy = gfn4(zi, radius_equat_core, polar_core, equat_shell, 64 polar_shell, delpc, delps, q); 65 summ += Gauss76Wt[i] * yyy; 72 66 } 67 summ *= 0.5*(uplim-lolim); 73 68 74 double answer = (uplim-lolim)/2.0*summ; 75 //convert to [cm-1] 76 answer *= 1.0e-4; 77 78 return answer; 69 // convert to [cm-1] 70 return 1.0e-4 * summ; 79 71 } 80 72 81 73 static double 82 core_shell_ellipsoid_xt_kernel_2d(double q , double q_x, double q_y,74 core_shell_ellipsoid_xt_kernel_2d(double qx, double qy, 83 75 double radius_equat_core, 84 76 double x_core, … … 91 83 double phi) 92 84 { 93 //convert angle degree to radian 94 theta = theta * M_PI_180; 95 phi = phi * M_PI_180; 96 97 // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis. 98 const double cyl_x = sin(theta) * cos(phi); 99 const double cyl_y = sin(theta) * sin(phi); 85 double q, sin_alpha, cos_alpha; 86 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 100 87 101 88 const double sldcs = core_sld - shell_sld; 102 89 const double sldss = shell_sld- solvent_sld; 103 104 // Compute the angle btw vector q and the105 // axis of the cylinder106 const double cos_val = cyl_x*q_x + cyl_y*q_y;107 90 108 91 const double polar_core = radius_equat_core*x_core; … … 112 95 // Call the IGOR library function to get the kernel: 113 96 // MUST use gfn4 not gf2 because of the def of params. 114 double answer = gfn4(cos_ val,97 double answer = gfn4(cos_alpha, 115 98 radius_equat_core, 116 99 polar_core, … … 160 143 double phi) 161 144 { 162 double q; 163 q = sqrt(qx*qx+qy*qy); 164 double intensity = core_shell_ellipsoid_xt_kernel_2d(q, qx/q, qy/q, 145 double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy, 165 146 radius_equat_core, 166 147 x_core, -
TabularUnified sasmodels/models/core_shell_ellipsoid.py ¶
ref5a314 r73e08ae 125 125 # ["name", "units", default, [lower, upper], "type", "description"], 126 126 parameters = [ 127 ["radius_equat_core","Ang", 20, [0, inf], 127 ["radius_equat_core","Ang", 20, [0, inf], "volume", "Equatorial radius of core"], 128 128 ["x_core", "None", 3, [0, inf], "volume", "axial ratio of core, X = r_polar/r_equatorial"], 129 129 ["thick_shell", "Ang", 30, [0, inf], "volume", "thickness of shell at equator"], -
TabularUnified sasmodels/models/core_shell_parallelepiped.c ¶
r2222134 r14838a3 35 35 // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 36 36 37 double t1, t2, t3, t4, tmp, answer; 38 double mu = q * length_b; 37 const double mu = 0.5 * q * length_b; 39 38 40 39 //calculate volume before rescaling (in original code, but not used) 41 //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 40 //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 42 41 //double vol = length_a * length_b * length_c + 43 42 // 2.0 * thick_rim_a * length_b * length_c + … … 46 45 47 46 // Scale sides by B 48 double a_scaled = length_a / length_b;49 double c_scaled = length_c / length_b;47 const double a_scaled = length_a / length_b; 48 const double c_scaled = length_c / length_b; 50 49 51 // DelRho values (note that drC is not used later) 52 double dr0 = core_sld-solvent_sld; 53 double drA = arim_sld-solvent_sld; 54 double drB = brim_sld-solvent_sld; 55 //double drC = crim_sld-solvent_sld; 50 // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) 51 // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c, 52 // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 53 // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!) 54 double ta = (a_scaled + 2.0*thick_rim_a)/length_b; 55 double tb = (a_scaled + 2.0*thick_rim_b)/length_b; 56 56 57 //Order of integration 58 int nordi=76; 59 int nordj=76; 60 61 // outer integral (with nordi gauss points), integration limits = 0, 1 62 double summ = 0; //initialize integral 57 double Vin = length_a * length_b * length_c; 58 //double Vot = (length_a * length_b * length_c + 59 // 2.0 * thick_rim_a * length_b * length_c + 60 // 2.0 * length_a * thick_rim_b * length_c + 61 // 2.0 * length_a * length_b * thick_rim_c); 62 double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 63 double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 63 64 64 for( int i=0; i<nordi; i++) { 65 66 // inner integral (with nordj gauss points), integration limits = 0, 1 67 68 double summj = 0.0; 69 double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 70 71 for(int j=0; j<nordj; j++) { 65 // Scale factors (note that drC is not used later) 66 const double drho0 = (core_sld-solvent_sld); 67 const double drhoA = (arim_sld-solvent_sld); 68 const double drhoB = (brim_sld-solvent_sld); 69 //const double drC_Vot = (crim_sld-solvent_sld)*Vot; 72 70 73 double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 74 double mudum = mu * sqrt(1.0-sigma*sigma); 71 // Precompute scale factors for combining cross terms from the shape 72 const double scale23 = drhoA*V1; 73 const double scale14 = drhoB*V2; 74 const double scale12 = drho0*Vin - scale23 - scale14; 75 75 76 double Vin = length_a * length_b * length_c; 77 //double Vot = (length_a * length_b * length_c + 78 // 2.0 * thick_rim_a * length_b * length_c + 79 // 2.0 * length_a * thick_rim_b * length_c + 80 // 2.0 * length_a * length_b * thick_rim_c); 81 double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 82 double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 83 84 // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) 85 // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c, 86 // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 87 // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!) 88 double ta = (a_scaled+2.0*thick_rim_a)/length_b; 89 double tb = (a_scaled+2.0*thick_rim_b)/length_b; 90 91 double arg1 = (0.5*mudum*a_scaled) * sin(0.5*M_PI*uu); 92 double arg2 = (0.5*mudum) * cos(0.5*M_PI*uu); 93 double arg3= (0.5*mudum*ta) * sin(0.5*M_PI*uu); 94 double arg4= (0.5*mudum*tb) * cos(0.5*M_PI*uu); 76 // outer integral (with gauss points), integration limits = 0, 1 77 double outer_total = 0; //initialize integral 95 78 96 if(arg1==0.0){ 97 t1 = 1.0; 98 } else { 99 t1 = (sin(arg1)/arg1); //defn for CSPP model sin(arg1)/arg1 test: (sin(arg1)/arg1)*(sin(arg1)/arg1) 100 } 101 if(arg2==0.0){ 102 t2 = 1.0; 103 } else { 104 t2 = (sin(arg2)/arg2); //defn for CSPP model sin(arg2)/arg2 test: (sin(arg2)/arg2)*(sin(arg2)/arg2) 105 } 106 if(arg3==0.0){ 107 t3 = 1.0; 108 } else { 109 t3 = sin(arg3)/arg3; 110 } 111 if(arg4==0.0){ 112 t4 = 1.0; 113 } else { 114 t4 = sin(arg4)/arg4; 115 } 116 79 for( int i=0; i<76; i++) { 80 double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 81 double mu_proj = mu * sqrt(1.0-sigma*sigma); 82 83 // inner integral (with gauss points), integration limits = 0, 1 84 double inner_total = 0.0; 85 for(int j=0; j<76; j++) { 86 const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 87 double sin_uu, cos_uu; 88 SINCOS(M_PI_2*uu, sin_uu, cos_uu); 89 const double si1 = sinc(mu_proj * sin_uu * a_scaled); 90 const double si2 = sinc(mu_proj * cos_uu); 91 const double si3 = sinc(mu_proj * sin_uu * ta); 92 const double si4 = sinc(mu_proj * cos_uu * tb); 93 117 94 // Expression in libCylinder.c (neither drC nor Vot are used) 118 tmp =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 ); // correct FF : square of sum of phase factors 119 95 const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 96 120 97 // To note also that in csparallelepiped.cpp, there is a function called 121 98 // cspkernel, which seems to make more sense and has the following comment: … … 126 103 // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 127 104 128 summj += Gauss76Wt[j] * tmp; 105 // correct FF : sum of square of phase factors 106 inner_total += Gauss76Wt[j] * form * form; 107 } 108 inner_total *= 0.5; 129 109 130 } 131 132 // value of the inner integral 133 answer = 0.5 * summj; 110 // now sum up the outer integral 111 const double si = sinc(mu * c_scaled * sigma); 112 outer_total += Gauss76Wt[i] * inner_total * si * si; 113 } 114 outer_total *= 0.5; 134 115 135 // finish the outer integral 136 double arg = 0.5 * mu* c_scaled *sigma; 137 if ( arg == 0.0 ) { 138 answer *= 1.0; 139 } else { 140 answer *= sin(arg)*sin(arg)/arg/arg; 141 } 142 143 // now sum up the outer integral 144 summ += Gauss76Wt[i] * answer; 145 146 } 147 148 answer = 0.5 * summ; 149 150 //convert from [1e-12 A-1] to [cm-1] 151 answer *= 1.0e-4; 152 153 return answer; 116 //convert from [1e-12 A-1] to [cm-1] 117 return 1.0e-4 * outer_total; 154 118 } 155 119 … … 170 134 double psi) 171 135 { 172 double tmp1, tmp2, tmp3, tmpt1, tmpt2, tmpt3; 136 double q, cos_val_a, cos_val_b, cos_val_c; 137 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 173 138 174 double q = sqrt(qx*qx+qy*qy); 175 double qx_scaled = qx/q; 176 double qy_scaled = qy/q; 139 // cspkernel in csparallelepiped recoded here 140 const double dr0 = core_sld-solvent_sld; 141 const double drA = arim_sld-solvent_sld; 142 const double drB = brim_sld-solvent_sld; 143 const double drC = crim_sld-solvent_sld; 177 144 178 // Convert angles given in degrees to radians 179 theta *= M_PI_180; 180 phi *= M_PI_180; 181 psi *= M_PI_180; 182 183 // Parallelepiped c axis orientation 184 double cparallel_x = cos(theta) * cos(phi); 185 double cparallel_y = sin(theta); 186 187 // Compute angle between q and parallelepiped axis 188 double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz; 189 190 // Parallelepiped a axis orientation 191 double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 192 double parallel_y = sin(psi)*cos(theta); 193 double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled; 194 195 // Parallelepiped b axis orientation 196 double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 197 double bparallel_y = cos(theta)*cos(psi); 198 double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled; 199 200 // The following tests should always pass 201 if (fabs(cos_val_c)>1.0) { 202 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 203 cos_val_c = 1.0; 204 } 205 if (fabs(cos_val_a)>1.0) { 206 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 207 cos_val_a = 1.0; 208 } 209 if (fabs(cos_val_b)>1.0) { 210 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 211 cos_val_b = 1.0; 212 } 213 214 // cspkernel in csparallelepiped recoded here 215 double dr0 = core_sld-solvent_sld; 216 double drA = arim_sld-solvent_sld; 217 double drB = brim_sld-solvent_sld; 218 double drC = crim_sld-solvent_sld; 219 double Vin = length_a * length_b * length_c; 145 double Vin = length_a * length_b * length_c; 146 double V1 = 2.0 * thick_rim_a * length_b * length_c; // incorrect V1(aa*bb*cc+2*ta*bb*cc) 147 double V2 = 2.0 * length_a * thick_rim_b * length_c; // incorrect V2(aa*bb*cc+2*aa*tb*cc) 148 double V3 = 2.0 * length_a * length_b * thick_rim_c; 220 149 // As for the 1D case, Vot is not used 221 150 //double Vot = (length_a * length_b * length_c + 222 151 // 2.0 * thick_rim_a * length_b * length_c + 223 152 // 2.0 * length_a * thick_rim_b * length_c + 224 153 // 2.0 * length_a * length_b * thick_rim_c); 225 double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 226 double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 227 double V3 = (2.0 * length_a * length_b * thick_rim_c); 154 228 155 // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 229 156 // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, 230 157 // but for the moment I let it like this until understanding better the code. 231 158 double ta = length_a + 2.0*thick_rim_a; 232 159 double tb = length_a + 2.0*thick_rim_b; 233 160 double tc = length_a + 2.0*thick_rim_c; 234 161 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 235 double argA = 0.5*q*length_a*cos_val_a;236 double argB = 0.5*q*length_b*cos_val_b;237 double argC = 0.5*q*length_c*cos_val_c;238 double argtA = 0.5*q*ta*cos_val_a;239 double argtB = 0.5*q*tb*cos_val_b;240 double argtC = 0.5*q*tc*cos_val_c;162 double siA = sinc(0.5*q*length_a*cos_val_a); 163 double siB = sinc(0.5*q*length_b*cos_val_b); 164 double siC = sinc(0.5*q*length_c*cos_val_c); 165 double siAt = sinc(0.5*q*ta*cos_val_a); 166 double siBt = sinc(0.5*q*tb*cos_val_b); 167 double siCt = sinc(0.5*q*tc*cos_val_c); 241 168 242 if(argA==0.0) {243 tmp1 = 1.0;244 } else {245 tmp1 = sin(argA)/argA;246 }247 if (argB==0.0) {248 tmp2 = 1.0;249 } else {250 tmp2 = sin(argB)/argB;251 }252 if (argC==0.0) {253 tmp3 = 1.0;254 } else {255 tmp3 = sin(argC)/argC;256 }257 if(argtA==0.0) {258 tmpt1 = 1.0;259 } else {260 tmpt1 = sin(argtA)/argtA;261 }262 if (argtB==0.0) {263 tmpt2 = 1.0;264 } else {265 tmpt2 = sin(argtB)/argtB;266 }267 if (argtC==0.0) {268 tmpt3 = 1.0;269 } else {270 tmpt3 = sin(argtC)*sin(argtC)/argtC/argtC;271 }272 169 273 170 // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 274 171 // in the 1D code, but should be checked! 275 double f = ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1 + 276 drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3); 172 double f = ( dr0*siA*siB*siC*Vin 173 + drA*(siAt-siA)*siB*siC*V1 174 + drB*siA*(siBt-siB)*siC*V2 175 + drC*siA*siB*(siCt*siCt-siC)*V3); 277 176 278 177 return 1.0e-4 * f * f; -
TabularUnified sasmodels/models/core_shell_parallelepiped.py ¶
r5810f00 r797a8e3 178 178 # parameters for demo 179 179 demo = dict(scale=1, background=0.0, 180 sld_core=1e-6, sld_a=2e-6, sld_b=4e-6, 181 sld_c=2e-6, sld_solvent=6e-6, 180 sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6, 182 181 length_a=35, length_b=75, length_c=400, 183 182 thick_rim_a=10, thick_rim_b=10, thick_rim_c=10, … … 191 190 theta_pd=10, theta_pd_n=1, 192 191 phi_pd=10, phi_pd_n=1, 193 psi_pd=10, psi_pd_n=1 0)192 psi_pd=10, psi_pd_n=1) 194 193 195 194 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) -
TabularUnified sasmodels/models/core_shell_sphere.c ¶
r2c74c11 r3a48772 16 16 double form_volume(double radius, double thickness) 17 17 { 18 return 4.0 * M_PI / 3.0 * pow((radius + thickness), 3);18 return M_4PI_3 * cube(radius + thickness); 19 19 } -
TabularUnified sasmodels/models/core_shell_sphere.py ¶
r9a4811a r4962519 95 95 """ 96 96 return (1, 1) 97 whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3)98 core = 4.0 * pi / 3.0 * radius * radius * radius97 whole = 4.0/3.0 * pi * (radius + thickness)**3 98 core = 4.0/3.0 * pi * radius**3 99 99 return whole, whole - core 100 100 -
TabularUnified sasmodels/models/cylinder.c ¶
r9cc7fca rb829b16 18 18 const double qr = q*radius; 19 19 const double qh = q*0.5*length; 20 return sas_J1c(qr*sn) * sinc(qh*cn);20 return sas_J1c(qr*sn) * sinc(qh*cn); 21 21 } 22 22 … … 58 58 double phi) 59 59 { 60 double sn, cn; // slots to hold sincos function output 61 62 // Compute angle alpha between q and the cylinder axis 63 SINCOS(phi*M_PI_180, sn, cn); 64 const double q = sqrt(qx*qx + qy*qy); 65 const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 66 67 const double alpha = acos(cos_val); 68 69 SINCOS(alpha, sn, cn); 60 double q, sin_alpha, cos_alpha; 61 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 62 //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha); 70 63 const double s = (sld-solvent_sld) * form_volume(radius, length); 71 return 1.0e-4 * square(s * fq(q, sn, cn, radius, length)); 64 const double form = fq(q, sin_alpha, cos_alpha, radius, length); 65 return 1.0e-4 * square(s * form); 72 66 } -
TabularUnified sasmodels/models/cylinder.py ¶
r0d6e865 r4cdd0cc 111 111 112 112 # [ "name", "units", default, [lower, upper], "type", "description"], 113 parameters = [["sld", " 4e-6/Ang^2", 4, [-inf, inf], "sld",113 parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", 114 114 "Cylinder scattering length density"], 115 115 ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld", -
TabularUnified sasmodels/models/dab.py ¶
ra807206 r4962519 59 59 60 60 Iq = """ 61 double numerator = pow(cor_length, 3);62 double denominator = pow(1 + pow(q*cor_length,2), 2);61 double numerator = cube(cor_length); 62 double denominator = square(1 + square(q*cor_length)); 63 63 64 64 return numerator / denominator ; -
TabularUnified sasmodels/models/ellipsoid.c ¶
r0d6e865 r73e08ae 8 8 { 9 9 double ratio = radius_polar/radius_equatorial; 10 const double u = q*radius_equatorial*sqrt(1.0 11 + sin_alpha*sin_alpha*(ratio*ratio - 1.0)); 12 const double f = sph_j1c(u); 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) 14 // Instead of using pythagoras we could pass in sin and cos; this may be 15 // slightly better for 2D which has already computed it, but it introduces 16 // an extra sqrt and square for 1-D not required by the current form, so 17 // leave it as is. 18 const double r = radius_equatorial 19 * sqrt(1.0 + sin_alpha*sin_alpha*(ratio*ratio - 1.0)); 20 const double f = sph_j1c(q*r); 13 21 14 22 return f*f; … … 49 57 double phi) 50 58 { 51 double sn, cn; 52 53 const double q = sqrt(qx*qx + qy*qy); 54 SINCOS(phi*M_PI_180, sn, cn); 55 const double cos_alpha = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 56 const double alpha = acos(cos_alpha); 57 SINCOS(alpha, sn, cn); 58 const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sn); 59 double q, sin_alpha, cos_alpha; 60 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 61 const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 59 62 const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 60 63 -
TabularUnified sasmodels/models/elliptical_cylinder.c ¶
ra807206 r251f54b 6 6 7 7 8 double _elliptical_cylinder_kernel(double q, double radius_minor, double r_ratio, double theta); 9 10 double _elliptical_cylinder_kernel(double q, double radius_minor, double r_ratio, double theta) 11 { 12 // This is the function LAMBDA1^2 in Feigin's notation 13 // q is the q-value for the calculation (1/A) 14 // radius_minor is the transformed radius"a" in Feigin's notation 15 // r_ratio is the ratio (major radius)/(minor radius) of the Ellipsoid [=] --- 16 // theta is the dummy variable of the integration 17 18 double retval,arg; 19 20 arg = q*radius_minor*sqrt((1.0+r_ratio*r_ratio)/2+(1.0-r_ratio*r_ratio)*cos(theta)/2); 21 //retval = 2.0*J1(arg)/arg; 22 retval = sas_J1c(arg); 23 return retval*retval ; 24 } 25 26 27 double form_volume(double radius_minor, double r_ratio, double length) 8 double 9 form_volume(double radius_minor, double r_ratio, double length) 28 10 { 29 11 return M_PI * radius_minor * radius_minor * r_ratio * length; 30 12 } 31 13 32 double Iq(double q, double radius_minor, double r_ratio, double length, 33 double sld, double solvent_sld) { 14 double 15 Iq(double q, double radius_minor, double r_ratio, double length, 16 double sld, double solvent_sld) 17 { 18 // orientational average limits 19 const double va = 0.0; 20 const double vb = 1.0; 21 // inner integral limits 22 const double vaj=0.0; 23 const double vbj=M_PI; 34 24 35 const int nordi=76; //order of integration 36 const int nordj=20; 37 double va,vb; //upper and lower integration limits 38 double summ,zi,yyy,answer; //running tally of integration 39 double summj,vaj,vbj,zij,arg,si; //for the inner integration 40 41 // orientational average limits 42 va = 0.0; 43 vb = 1.0; 44 // inner integral limits 45 vaj=0.0; 46 vbj=M_PI; 25 const double radius_major = r_ratio * radius_minor; 26 const double rA = 0.5*(square(radius_major) + square(radius_minor)); 27 const double rB = 0.5*(square(radius_major) - square(radius_minor)); 47 28 48 29 //initialize integral 49 summ = 0.0; 50 51 const double delrho = sld - solvent_sld; 52 53 for(int i=0;i<nordi;i++) { 30 double outer_sum = 0.0; 31 for(int i=0;i<76;i++) { 54 32 //setup inner integral over the ellipsoidal cross-section 55 summj=0; 56 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy 57 arg = radius_minor*sqrt(1.0-zi*zi); 58 for(int j=0;j<nordj;j++) { 33 const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 34 const double sin_val = sqrt(1.0 - cos_val*cos_val); 35 //const double arg = radius_minor*sin_val; 36 double inner_sum=0; 37 for(int j=0;j<20;j++) { 59 38 //20 gauss points for the inner integral 60 zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy 61 yyy = Gauss20Wt[j] * _elliptical_cylinder_kernel(q, arg, r_ratio, zij); 62 summj += yyy; 39 const double theta = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 40 const double r = sin_val*sqrt(rA - rB*cos(theta)); 41 const double be = sas_J1c(q*r); 42 inner_sum += Gauss20Wt[j] * be * be; 63 43 } 64 44 //now calculate the value of the inner integral 65 answer = (vbj-vaj)/2.0*summj;45 inner_sum *= 0.5*(vbj-vaj); 66 46 67 47 //now calculate outer integral 68 arg = q*length*zi/2.0; 69 si = square(sinc(arg)); 70 yyy = Gauss76Wt[i] * answer * si; 71 summ += yyy; 48 const double si = sinc(q*0.5*length*cos_val); 49 outer_sum += Gauss76Wt[i] * inner_sum * si * si; 72 50 } 51 outer_sum *= 0.5*(vb-va); 73 52 74 53 //divide integral by Pi 75 answer = (vb-va)/2.0*summ/M_PI; 76 // Multiply by contrast^2 77 answer *= delrho*delrho; 54 const double form = outer_sum/M_PI; 78 55 56 // scale by contrast and volume, and convert to to 1/cm units 79 57 const double vol = form_volume(radius_minor, r_ratio, length); 80 return answer*vol*vol*1.0e-4; 58 const double delrho = sld - solvent_sld; 59 return 1.0e-4*square(delrho*vol)*form; 81 60 } 82 61 83 62 84 double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length,85 double sld, double solvent_sld, double theta, double phi, double psi) { 86 const double _theta = theta * M_PI / 180.0;87 const double _phi = phi * M_PI / 180.0;88 const double _psi = psi * M_PI / 180.0;89 const double q = sqrt(qx*qx+qy*qy); 90 const double q_x = qx/q;91 const double q_y = qy/q;63 double 64 Iqxy(double qx, double qy, 65 double radius_minor, double r_ratio, double length, 66 double sld, double solvent_sld, 67 double theta, double phi, double psi) 68 { 69 double q, cos_val, cos_mu, cos_nu; 70 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val, cos_mu, cos_nu); 92 71 93 //Cylinder orientation 94 double cyl_x = cos(_theta) * cos(_phi); 95 double cyl_y = sin(_theta); 96 97 //cyl_z = -cos(_theta) * sin(_phi); 98 99 // q vector 100 //q_z = 0; 101 102 // Note: cos(alpha) = 0 and 1 will get an 103 // undefined value from CylKernel 104 //alpha = acos( cos_val ); 105 106 //ellipse orientation: 107 // the elliptical corss section was transformed and projected 108 // into the detector plane already through sin(alpha)and furthermore psi remains as same 109 // on the detector plane. 110 // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt 111 // the wave vector q. 112 113 //x- y- component on the detector plane. 114 const double ella_x = -cos(_phi)*sin(_psi) * sin(_theta)+sin(_phi)*cos(_psi); 115 const double ella_y = sin(_psi)*cos(_theta); 116 const double ellb_x = -sin(_theta)*cos(_psi)*cos(_phi)-sin(_psi)*sin(_phi); 117 const double ellb_y = cos(_theta)*cos(_psi); 118 119 // Compute the angle btw vector q and the 120 // axis of the cylinder 121 double cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; 122 123 // calculate the axis of the ellipse wrt q-coord. 124 double cos_nu = ella_x*q_x + ella_y*q_y; 125 double cos_mu = ellb_x*q_x + ellb_y*q_y; 126 127 // The following test should always pass 128 if (fabs(cos_val)>1.0) { 129 //printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n"); 130 cos_val = 1.0; 131 } 132 if (fabs(cos_nu)>1.0) { 133 //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n"); 134 cos_nu = 1.0; 135 } 136 if (fabs(cos_mu)>1.0) { 137 //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n"); 138 cos_mu = 1.0; 139 } 140 141 const double r_major = r_ratio * radius_minor; 142 const double qr = q*sqrt( r_major*r_major*cos_nu*cos_nu + radius_minor*radius_minor*cos_mu*cos_mu ); 143 const double qL = q*length*cos_val/2.0; 144 145 double Be; 146 if (qr==0){ 147 Be = 0.5; 148 }else{ 149 //Be = NR_BessJ1(qr)/qr; 150 Be = 0.5*sas_J1c(qr); 151 } 152 153 double Si; 154 if (qL==0){ 155 Si = 1.0; 156 }else{ 157 Si = sin(qL)/qL; 158 } 159 160 const double k = 2.0 * Be * Si; 72 // Compute: r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 73 // Given: radius_major = r_ratio * radius_minor 74 const double r = radius_minor*sqrt(square(r_ratio*cos_nu) + cos_mu*cos_mu); 75 const double be = sas_J1c(q*r); 76 const double si = sinc(q*0.5*length*cos_val); 77 const double Aq = be * si; 78 const double delrho = sld - solvent_sld; 161 79 const double vol = form_volume(radius_minor, r_ratio, length); 162 return (sld - solvent_sld) * (sld - solvent_sld) * k * k *vol*vol*1.0e-4;80 return 1.0e-4 * square(delrho * vol * Aq); 163 81 } -
TabularUnified sasmodels/models/fcc_paracrystal.c ¶
r0bef47b r4962519 32 32 33 33 const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 34 result = pow((1.0-(temp10*temp10)),3)*temp634 result = cube(1.0-(temp10*temp10))*temp6 35 35 / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 36 36 * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) … … 85 85 } 86 86 87 double Iqxy(double qx, double qy, 88 double dnn, double d_factor, double radius, 89 double sld, double solvent_sld, 90 double theta, double phi, double psi) 91 { 92 double q, cos_a1, cos_a2, cos_a3; 93 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 87 94 88 double Iqxy(double qx, double qy, double dnn, 89 double d_factor, double radius,double sld, double solvent_sld, 90 double theta, double phi, double psi){ 95 const double a1 = cos_a2 + cos_a3; 96 const double a2 = cos_a3 + cos_a1; 97 const double a3 = cos_a2 + cos_a1; 98 const double qd = 0.5*q*dnn; 99 const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 100 const double tanh_qd = tanh(arg); 101 const double cosh_qd = cosh(arg); 102 const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 103 * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 104 * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 91 105 92 double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z, 93 // double q_z; 94 double cos_val_b3, cos_val_b2, cos_val_b1; 95 double a1_dot_q, a2_dot_q,a3_dot_q; 96 double answer; 97 double Zq, Fkq, Fkq_2; 106 //if (isnan(Zq)) printf("q:(%g,%g) qd: %g a1: %g a2: %g a3: %g arg: %g\n", qx, qy, qd, a1, a2, a3, arg); 98 107 99 //convert to q and make scaled values 100 double q = sqrt(qx*qx+qy*qy); 101 double q_x = qx/q; 102 double q_y = qy/q; 103 104 //convert angle degree to radian 105 theta = theta * M_PI_180; 106 phi = phi * M_PI_180; 107 psi = psi * M_PI_180; 108 109 const double Da = d_factor*dnn; 110 const double s1 = dnn/sqrt(0.75); 111 112 113 //the occupied volume of the lattice 114 const double latticescale = 2.0*sphere_volume(radius/s1); 115 // q vector 116 // q_z = 0.0; // for SANS; assuming qz is negligible 117 /// Angles here are respect to detector coordinate 118 /// instead of against q coordinate(PRB 36(46), 3(6), 1754(3854)) 119 // b3 axis orientation 120 b3_x = cos(theta) * cos(phi); 121 b3_y = sin(theta); 122 //b3_z = -cos(theta) * sin(phi); 123 cos_val_b3 = b3_x*q_x + b3_y*q_y;// + b3_z*q_z; 124 125 //alpha = acos(cos_val_b3); 126 // b1 axis orientation 127 b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 128 b1_y = sin(psi)*cos(theta); 129 cos_val_b1 = b1_x*q_x + b1_y*q_y; 130 // b2 axis orientation 131 b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 132 b2_y = cos(theta)*cos(psi); 133 cos_val_b2 = b2_x*q_x + b2_y*q_y; 134 135 // The following test should always pass 136 if (fabs(cos_val_b3)>1.0) { 137 //printf("FCC_ana_2D: Unexpected error: cos()>1\n"); 138 cos_val_b3 = 1.0; 139 } 140 if (fabs(cos_val_b2)>1.0) { 141 //printf("FCC_ana_2D: Unexpected error: cos()>1\n"); 142 cos_val_b2 = 1.0; 143 } 144 if (fabs(cos_val_b1)>1.0) { 145 //printf("FCC_ana_2D: Unexpected error: cos()>1\n"); 146 cos_val_b1 = 1.0; 147 } 148 // Compute the angle btw vector q and the a3 axis 149 a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3); 150 151 // a1 axis 152 a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1); 153 154 // a2 axis 155 a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2); 156 157 158 // Get Fkq and Fkq_2 159 Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q)); 160 Fkq_2 = Fkq*Fkq; 161 // Call Zq=Z1*Z2*Z3 162 Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2); 163 Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2); 164 Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2); 165 166 // Use SphereForm directly from libigor 167 answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale; 168 169 return answer; 170 } 108 const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 109 //the occupied volume of the lattice 110 const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 111 return lattice_scale * Fq; 112 } -
TabularUnified sasmodels/models/flexible_cylinder_elliptical.c ¶
rb66d38e r92ce163 2 2 double Iq(double q, double length, double kuhn_length, double radius, 3 3 double axis_ratio, double sld, double solvent_sld); 4 double Iqxy(double qx, double qy, double length, double kuhn_length,5 double radius, double axis_ratio, double sld, double solvent_sld);6 4 double flexible_cylinder_ex_kernel(double q, double length, double kuhn_length, 7 5 double radius, double axis_ratio, double sld, … … 17 15 elliptical_crosssection(double q, double a, double b) 18 16 { 19 double sum m=0.0;17 double sum=0.0; 20 18 21 19 for(int i=0;i<N_POINTS_76;i++) { 22 double zi = ( Gauss76Z[i] + 1.0 )*M_PI/4.0;20 const double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4; 23 21 double sn, cn; 24 22 SINCOS(zi, sn, cn); 25 double arg = q*sqrt(a*a*sn*sn+b*b*cn*cn); 26 double yyy = pow((double)sas_J1c(arg),2); 27 yyy *= Gauss76Wt[i]; 28 summ += yyy; 23 const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn); 24 const double yyy = sas_J1c(arg); 25 sum += Gauss76Wt[i] * yyy * yyy; 29 26 } 30 31 summ /= 2.0; 32 return(summ); 27 sum *= 0.5; 28 return(sum); 33 29 34 30 } … … 77 73 } 78 74 79 double Iqxy(double qx, double qy,80 double length,81 double kuhn_length,82 double radius,83 double axis_ratio,84 double sld,85 double solvent_sld)86 {87 double q;88 q = sqrt(qx*qx+qy*qy);89 double result = flexible_cylinder_ex_kernel(q,90 length,91 kuhn_length,92 radius,93 axis_ratio,94 sld,95 solvent_sld);96 97 return result;98 } -
TabularUnified sasmodels/models/fractal.c ¶
r2c74c11 r217590b 1 double Iq(double q, 2 double volfraction, 3 double radius, 4 double fractal_dim, 5 double cor_length, 6 double sld_block, 7 double sld_solvent); 1 #define INVALID(p) (p.fractal_dim < 0.0) 8 2 9 double Iq(double q, 10 double volfraction, 11 double radius, 12 double fractal_dim, 13 double cor_length, 14 double sld_block, 15 double sld_solvent) 3 static double 4 Iq(double q, 5 double volfraction, 6 double radius, 7 double fractal_dim, 8 double cor_length, 9 double sld_block, 10 double sld_solvent) 16 11 { 17 double qr,r0,Df,corr,phi,sldp,sldm; 18 double pq,sq,inten; 19 20 // Actively check the argument - needed for mass fractal - is it needie 21 //here? 22 if (fractal_dim <= 0.0){ 23 return 0.0; 24 } 25 26 phi = volfraction; // volume fraction of building block spheres... 27 r0 = radius; // radius of building block 28 Df = fractal_dim; // fractal dimension 29 corr = cor_length; // correlation length of fractal-like aggregates 30 sldp = sld_block; // SLD of building block 31 sldm = sld_solvent; // SLD of matrix or solution 32 33 qr=q*r0; 34 12 const double sq = fractal_sq(q, radius, fractal_dim, cor_length); 13 35 14 //calculate P(q) for the spherical subunits 36 pq = phi*M_4PI_3*r0*r0*r0*(sldp-sldm)*(sldp-sldm)*sph_j1c(qr)*sph_j1c(qr); 37 38 //calculate S(q) 39 sq = Df*sas_gamma(Df-1.0)*sin((Df-1.0)*atan(q*corr)); 40 sq /= pow(qr,Df) * pow((1.0 + 1.0/(q*corr)/(q*corr)),((Df-1.0)/2.0)); 41 sq += 1.0; 42 43 //combine, scale to units cm-1 sr-1 (assuming data on absolute scale) 44 //and return 45 inten = pq*sq; 46 // convert I(1/A) to (1/cm) 47 inten *= 1.0e8; 48 //convert rho^2 in 10^-6 1/A to 1/A 49 inten *= 1.0e-12; 50 51 52 return(inten); 15 const double V = M_4PI_3*cube(radius); 16 const double pq = V * square((sld_block-sld_solvent)*sph_j1c(q*radius)); 17 18 // scale to units cm-1 sr-1 (assuming data on absolute scale) 19 // convert I(1/A) to (1/cm) => 1e8 * I(q) 20 // convert rho^2 in 10^-6 1/A to 1/A => 1e-12 * I(q) 21 // combined: 1e-4 * I(q) 22 23 return 1.e-4 * volfraction * sq * pq; 53 24 } 54 25 -
TabularUnified sasmodels/models/fractal.py ¶
r785cbec rfef353f 30 30 31 31 where $\xi$ is the correlation length representing the cluster size and $D_f$ 32 is the fractal dimension, representing the self similarity of the structure. 32 is the fractal dimension, representing the self similarity of the structure. 33 Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function 34 diverges at $D_f$=0 and $D_f$=1. 33 35 34 36 **Polydispersity on the radius is provided for.** … … 95 97 # pylint: enable=bad-whitespace, line-too-long 96 98 97 source = ["lib/sph_j1c.c", "lib/sas_gamma.c", " fractal.c"]99 source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"] 98 100 99 101 demo = dict(volfraction=0.05, -
TabularUnified sasmodels/models/fractal_core_shell.c ¶
ra807206 r217590b 1 double form_volume(double radius, double thickness); 2 3 double Iq(double q, 4 double radius, 5 double thickness, 6 double core_sld, 7 double shell_sld, 8 double solvent_sld, 9 double volfraction, 10 double fractal_dim, 11 double cor_length); 12 13 double form_volume(double radius, double thickness) 1 static double 2 form_volume(double radius, double thickness) 14 3 { 15 return 4.0 * M_PI / 3.0 * pow((radius + thickness), 3);4 return M_4PI_3 * cube(radius + thickness); 16 5 } 17 6 18 double Iq(double q, 19 double radius, 20 double thickness, 21 double core_sld, 22 double shell_sld, 23 double solvent_sld, 24 double volfraction, 25 double fractal_dim, 26 double cor_length) { 7 static double 8 Iq(double q, 9 double radius, 10 double thickness, 11 double core_sld, 12 double shell_sld, 13 double solvent_sld, 14 double volfraction, 15 double fractal_dim, 16 double cor_length) 17 { 18 const double sq = fractal_sq(q, radius, fractal_dim, cor_length); 19 const double pq = core_shell_kernel(q, radius, thickness, 20 core_sld, shell_sld, solvent_sld); 27 21 28 29 double intensity = core_shell_kernel(q, 30 radius, 31 thickness, 32 core_sld, 33 shell_sld, 34 solvent_sld); 35 //calculate S(q) 36 double frac_1 = fractal_dim-1.0; 37 double qr = q*radius; 38 39 double t1 = fractal_dim*sas_gamma(frac_1)*sin(frac_1*atan(q*cor_length)); 40 double t2 = (1.0 + 1.0/(q*cor_length)/(q*cor_length)); 41 double t3 = pow(qr, fractal_dim) * pow(t2, (frac_1/2.0)); 42 double sq = t1/t3; 43 sq += 1.0; 44 45 return sq*intensity*volfraction; 22 // Note: core_shell_kernel already performs the 1e-4 unit conversion 23 return volfraction * sq * pq; 46 24 } 47 25 -
TabularUnified sasmodels/models/fractal_core_shell.py ¶
ra807206 r217590b 63 63 # ["name", "units", default, [lower, upper], "type","description"], 64 64 parameters = [ 65 ["radius", "Ang", 60.0, [0 , inf],"volume", "Sphere core radius"],66 ["thickness", "Ang", 10.0, [0 , inf],"volume", "Sphere shell thickness"],65 ["radius", "Ang", 60.0, [0.0, inf], "volume", "Sphere core radius"], 66 ["thickness", "Ang", 10.0, [0.0, inf], "volume", "Sphere shell thickness"], 67 67 ["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", "Sphere core scattering length density"], 68 68 ["sld_shell", "1e-6/Ang^2", 2.0, [-inf, inf], "sld", "Sphere shell scattering length density"], 69 69 ["sld_solvent", "1e-6/Ang^2", 3.0, [-inf, inf], "sld", "Solvent scattering length density"], 70 ["volfraction", "", 1.0, [0, inf], "", "Volume fraction of building block spheres"], 71 ["fractal_dim", "", 2.0, [-inf, inf], "", "Fractal dimension"], 72 ["cor_length", "Ang", 100.0, [0, inf], "", "Correlation length of fractal-like aggregates"]] 70 ["volfraction", "", 1.0, [0.0, inf], "", "Volume fraction of building block spheres"], 71 ["fractal_dim", "", 2.0, [0.0, 6.0], "", "Fractal dimension"], 72 ["cor_length", "Ang", 100.0, [0.0, inf], "", "Correlation length of fractal-like aggregates"], 73 ] 73 74 # pylint: enable=bad-whitespace, line-too-long 74 75 75 source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c", "fractal_core_shell.c"] 76 source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c", 77 "lib/fractal_sq.c", "fractal_core_shell.c"] 76 78 77 79 demo = dict(scale=0.05, … … 100 102 @param thickness: shell thickness 101 103 """ 102 whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3)103 core = 4.0 * pi / 3.0 * radius * radius * radius104 whole = 4.0/3.0 * pi * (radius + thickness)**3 105 core = 4.0/3.0 * pi * radius**3 104 106 return whole, whole-core 105 107 -
TabularUnified sasmodels/models/fuzzy_sphere.py ¶
r9a4811a r3a48772 86 86 # This should perhaps be volume normalized? 87 87 form_volume = """ 88 return 1.333333333333333*M_PI*radius*radius*radius;88 return M_4PI_3*cube(radius); 89 89 """ 90 90 -
TabularUnified sasmodels/models/hayter_msa.c ¶
r0bef47b r4962519 23 23 double SIdiam, diam, Kappa, cs, IonSt; 24 24 double Perm, Beta; 25 double pi,charge;25 double charge; 26 26 int ierr; 27 27 28 pi = M_PI;29 30 28 diam=2*radius_effective; //in A 31 29 … … 38 36 charge=zz*Elcharge; //in Coulomb (C) 39 37 SIdiam = diam*1.0E-10; //in m 40 Vp= 4.0*pi/3.0*(SIdiam/2.0)*(SIdiam/2.0)*(SIdiam/2.0); //in m^338 Vp=M_4PI_3*cube(SIdiam/2.0); //in m^3 41 39 cs=csalt*6.022E23*1.0E3; //# salt molecules/m^3 42 40 … … 50 48 Kappa=sqrt(2*Beta*IonSt/Perm); //Kappa calc from Ionic strength 51 49 // Kappa=2/SIdiam // Use to compare with HP paper 52 gMSAWave[5]=Beta*charge*charge/( pi*Perm*SIdiam*pow((2.0+Kappa*SIdiam),2));50 gMSAWave[5]=Beta*charge*charge/(M_PI*Perm*SIdiam*square(2.0+Kappa*SIdiam)); 53 51 54 52 // Finally set up dimensionless parameters -
TabularUnified sasmodels/models/hollow_cylinder.c ¶
r0d6e865 rf8f0991 1 1 double form_volume(double radius, double thickness, double length); 2 3 2 double Iq(double q, double radius, double thickness, double length, double sld, 4 3 double solvent_sld); … … 9 8 10 9 // From Igor library 11 static double hollow_cylinder_scaling(12 10 static double 11 _hollow_cylinder_scaling(double integrand, double delrho, double volume) 13 12 { 14 double answer; 15 // Multiply by contrast^2 16 answer = integrand*delrho*delrho; 17 18 //normalize by cylinder volume 19 answer *= volume*volume; 20 21 //convert to [cm-1] 22 answer *= 1.0e-4; 23 24 return answer; 13 return 1.0e-4 * square(volume * delrho) * integrand; 25 14 } 26 15 27 16 28 static double _hollow_cylinder_kernel( 29 double q, double radius, double thickness, double length, double dum) 17 static double 18 _hollow_cylinder_kernel(double q, 19 double radius, double thickness, double length, double sin_val, double cos_val) 30 20 { 31 const double qs = q*s qrt(1.0-dum*dum);21 const double qs = q*sin_val; 32 22 const double lam1 = sas_J1c((radius+thickness)*qs); 33 23 const double lam2 = sas_J1c(radius*qs); 34 24 const double gamma_sq = square(radius/(radius+thickness)); 35 //Note: lim_{r -> r_c} psi = J0(radius_core*qs) 25 //Note: lim_{thickness -> 0} psi = J0(radius*qs) 26 //Note: lim_{radius -> 0} psi = sas_J1c(thickness*qs) 36 27 const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 37 const double t2 = sinc(q*length*dum/2.0); 38 return square(psi*t2); 28 const double t2 = sinc(0.5*q*length*cos_val); 29 return psi*t2; 30 } 31 32 double 33 form_volume(double radius, double thickness, double length) 34 { 35 double v_shell = M_PI*length*(square(radius+thickness) - radius*radius); 36 return v_shell; 39 37 } 40 38 41 39 42 static double hollow_cylinder_analytical_2D_scaled( 43 double q, double q_x, double q_y, double radius, double thickness,44 double length, double sld, double solvent_sld, double theta, double phi)40 double 41 Iq(double q, double radius, double thickness, double length, 42 double sld, double solvent_sld) 45 43 { 46 double cyl_x, cyl_y; //, cyl_z 47 //double q_z; 48 double vol, cos_val, delrho; 49 double answer; 50 //convert angle degree to radian 51 theta = theta * M_PI_180; 52 phi = phi * M_PI_180; 53 delrho = solvent_sld - sld; 44 const double lower = 0.0; 45 const double upper = 1.0; //limits of numerical integral 54 46 55 // Cylinder orientation 56 cyl_x = sin(theta) * cos(phi); 57 cyl_y = sin(theta) * sin(phi); 58 //cyl_z = -cos(theta) * sin(phi); 47 double summ = 0.0; //initialize intergral 48 for (int i=0;i<76;i++) { 49 const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 50 const double sin_val = sqrt(1.0 - cos_val*cos_val); 51 const double inter = _hollow_cylinder_kernel(q, radius, thickness, length, 52 sin_val, cos_val); 53 summ += Gauss76Wt[i] * inter * inter; 54 } 59 55 60 // q vector 61 //q_z = 0; 62 63 // Compute the angle btw vector q and the 64 // axis of the cylinder 65 cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; 66 67 answer = _hollow_cylinder_kernel(q, radius, thickness, length, cos_val); 68 69 vol = form_volume(radius, thickness, length); 70 answer = hollow_cylinder_scaling(answer, delrho, vol); 71 72 return answer; 56 const double Aq = 0.5*summ*(upper-lower); 57 const double volume = form_volume(radius, thickness, length); 58 return _hollow_cylinder_scaling(Aq, solvent_sld - sld, volume); 73 59 } 74 60 61 double 62 Iqxy(double qx, double qy, 63 double radius, double thickness, double length, 64 double sld, double solvent_sld, double theta, double phi) 65 { 66 double q, sin_alpha, cos_alpha; 67 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 68 const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 69 sin_alpha, cos_alpha); 75 70 76 double form_volume(double radius, double thickness, double length) 77 { 78 double v_shell = M_PI*length*((radius+thickness)*(radius+thickness)-radius*radius); 79 return(v_shell); 71 const double vol = form_volume(radius, thickness, length); 72 return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol); 80 73 } 81 74 82 83 double Iq(double q, double radius, double thickness, double length,84 double sld, double solvent_sld)85 {86 int i;87 double lower,upper,zi, inter; //upper and lower integration limits88 double summ,answer,delrho; //running tally of integration89 double norm,volume; //final calculation variables90 91 lower = 0.0;92 upper = 1.0; //limits of numerical integral93 94 summ = 0.0; //initialize intergral95 for (i=0;i<76;i++) {96 zi = ( Gauss76Z[i] * (upper-lower) + lower + upper )/2.0;97 inter = Gauss76Wt[i] * _hollow_cylinder_kernel(q, radius, thickness, length, zi);98 summ += inter;99 }100 101 norm = summ*(upper-lower)/2.0;102 volume = form_volume(radius, thickness, length);103 delrho = solvent_sld - sld;104 answer = hollow_cylinder_scaling(norm, delrho, volume);105 106 return(answer);107 }108 109 110 double Iqxy(double qx, double qy, double radius, double thickness,111 double length, double sld, double solvent_sld, double theta, double phi)112 {113 const double q = sqrt(qx*qx+qy*qy);114 return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, thickness, length, sld, solvent_sld, theta, phi);115 } -
TabularUnified sasmodels/models/hollow_rectangular_prism.c ¶
ra807206 r6f676fb 2 2 double Iq(double q, double sld, double solvent_sld, double length_a, 3 3 double b2a_ratio, double c2a_ratio, double thickness); 4 double Iqxy(double qx, double qy, double sld, double solvent_sld,5 double length_a, double b2a_ratio, double c2a_ratio, double thickness);6 4 7 5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 8 6 { 9 double b_side= length_a * b2a_ratio;10 double c_side= length_a * c2a_ratio;7 double length_b = length_a * b2a_ratio; 8 double length_c = length_a * c2a_ratio; 11 9 double a_core = length_a - 2.0*thickness; 12 double b_core = b_side- 2.0*thickness;13 double c_core = c_side- 2.0*thickness;10 double b_core = length_b - 2.0*thickness; 11 double c_core = length_c - 2.0*thickness; 14 12 double vol_core = a_core * b_core * c_core; 15 double vol_total = length_a * b_side * c_side;13 double vol_total = length_a * length_b * length_c; 16 14 double vol_shell = vol_total - vol_core; 17 15 return vol_shell; … … 26 24 double thickness) 27 25 { 28 double termA1, termA2, termB1, termB2, termC1, termC2; 26 const double length_b = length_a * b2a_ratio; 27 const double length_c = length_a * c2a_ratio; 28 const double a_half = 0.5 * length_a; 29 const double b_half = 0.5 * length_b; 30 const double c_half = 0.5 * length_c; 31 const double vol_total = length_a * length_b * length_c; 32 const double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness); 33 34 //Integration limits to use in Gaussian quadrature 35 const double v1a = 0.0; 36 const double v1b = M_PI_2; //theta integration limits 37 const double v2a = 0.0; 38 const double v2b = M_PI_2; //phi integration limits 29 39 30 double b_side = length_a * b2a_ratio; 31 double c_side = length_a * c2a_ratio; 32 double a_half = 0.5 * length_a; 33 double b_half = 0.5 * b_side; 34 double c_half = 0.5 * c_side; 40 double outer_sum = 0.0; 41 for(int i=0; i<76; i++) { 35 42 36 //Integration limits to use in Gaussian quadrature 37 double v1a = 0.0; 38 double v1b = 0.5 * M_PI; //theta integration limits 39 double v2a = 0.0; 40 double v2b = 0.5 * M_PI; //phi integration limits 41 42 //Order of integration 43 int nordi=76; 44 int nordj=76; 43 const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 44 double sin_theta, cos_theta; 45 SINCOS(theta, sin_theta, cos_theta); 45 46 46 double sumi = 0.0; 47 48 for(int i=0; i<nordi; i++) { 47 const double termC1 = sinc(q * c_half * cos(theta)); 48 const double termC2 = sinc(q * (c_half-thickness)*cos(theta)); 49 49 50 double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 50 double inner_sum = 0.0; 51 for(int j=0; j<76; j++) { 51 52 52 double arg = q * c_half * cos(theta); 53 if (fabs(arg) > 1.e-16) {termC1 = sin(arg)/arg;} else {termC1 = 1.0;} 54 arg = q * (c_half-thickness)*cos(theta); 55 if (fabs(arg) > 1.e-16) {termC2 = sin(arg)/arg;} else {termC2 = 1.0;} 56 57 double sumj = 0.0; 58 59 for(int j=0; j<nordj; j++) { 60 61 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 53 const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 54 double sin_phi, cos_phi; 55 SINCOS(phi, sin_phi, cos_phi); 62 56 63 57 // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0 64 58 65 arg = q * a_half * sin(theta) * sin(phi); 66 if (fabs(arg) > 1.e-16) {termA1 = sin(arg)/arg;} else {termA1 = 1.0;} 67 arg = q * (a_half-thickness) * sin(theta) * sin(phi); 68 if (fabs(arg) > 1.e-16) {termA2 = sin(arg)/arg;} else {termA2 = 1.0;} 59 const double termA1 = sinc(q * a_half * sin_theta * sin_phi); 60 const double termA2 = sinc(q * (a_half-thickness) * sin_theta * sin_phi); 69 61 70 arg = q * b_half * sin(theta) * cos(phi); 71 if (fabs(arg) > 1.e-16) {termB1 = sin(arg)/arg;} else {termB1 = 1.0;} 72 arg = q * (b_half-thickness) * sin(theta) * cos(phi); 73 if (fabs(arg) > 1.e-16) {termB2 = sin(arg)/arg;} else {termB2 = 1.0;} 62 const double termB1 = sinc(q * b_half * sin_theta * cos_phi); 63 const double termB2 = sinc(q * (b_half-thickness) * sin_theta * cos_phi); 74 64 75 double AP1 = (length_a*b_side*c_side) * termA1 * termB1 * termC1; 76 double AP2 = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness) * termA2 * termB2 * termC2; 77 double AP = AP1 - AP2; 65 const double AP1 = vol_total * termA1 * termB1 * termC1; 66 const double AP2 = vol_core * termA2 * termB2 * termC2; 78 67 79 sumj += Gauss76Wt[j] * (AP*AP); 68 inner_sum += Gauss76Wt[j] * square(AP1-AP2); 69 } 70 inner_sum *= 0.5 * (v2b-v2a); 80 71 81 } 82 83 sumj = 0.5 * (v2b-v2a) * sumj; 84 sumi += Gauss76Wt[i] * sumj * sin(theta); 85 72 outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 86 73 } 87 88 double answer = 0.5*(v1b-v1a)*sumi; 74 outer_sum *= 0.5*(v1b-v1a); 89 75 90 76 // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 91 77 // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 92 answer *= (2.0/M_PI);78 const double form = outer_sum/M_PI_2; 93 79 94 80 // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 95 answer *= (sld-solvent_sld)*(sld-solvent_sld);81 const double delrho = sld - solvent_sld; 96 82 97 83 // Convert from [1e-12 A-1] to [cm-1] 98 answer *= 1.0e-4; 99 100 return answer; 101 84 return 1.0e-4 * delrho * delrho * form; 102 85 } 103 104 double Iqxy(double qx, double qy,105 double sld,106 double solvent_sld,107 double length_a,108 double b2a_ratio,109 double c2a_ratio,110 double thickness)111 {112 double q = sqrt(qx*qx + qy*qy);113 double intensity = Iq(q, sld, solvent_sld, length_a, b2a_ratio, c2a_ratio, thickness);114 return intensity;115 } -
TabularUnified sasmodels/models/hollow_rectangular_prism.py ¶
ra807206 rab2aea8 13 13 the difference of the amplitudes of two massive parallelepipeds 14 14 differing in their outermost dimensions in each direction by the 15 same length increment :math:`2\Delta`(Nayuk, 2012).15 same length increment $2\Delta$ (Nayuk, 2012). 16 16 17 17 As in the case of the massive parallelepiped model (:ref:`rectangular-prism`), … … 58 58 59 59 .. math:: 60 I(q) = \text{scale} \times V \times (\rho_ {\text{p}} -61 \rho_ {\text{solvent}})^2 \times P(q) + \text{background}60 I(q) = \text{scale} \times V \times (\rho_\text{p} - 61 \rho_\text{solvent})^2 \times P(q) + \text{background} 62 62 63 where $\rho_ {\text{p}}$ is the scattering length of the parallelepiped,64 $\rho_ {\text{solvent}}$ is the scattering length of the solvent,63 where $\rho_\text{p}$ is the scattering length of the parallelepiped, 64 $\rho_\text{solvent}$ is the scattering length of the solvent, 65 65 and (if the data are in absolute units) *scale* represents the volume fraction 66 66 (which is unitless). … … 148 148 # parameters for demo 149 149 demo = dict(scale=1, background=0, 150 sld=6.3 e-6, sld_solvent=1.0e-6,150 sld=6.3, sld_solvent=1.0, 151 151 length_a=35, b2a_ratio=1, c2a_ratio=1, thickness=1, 152 152 length_a_pd=0.1, length_a_pd_n=10, -
TabularUnified sasmodels/models/hollow_rectangular_prism_thin_walls.c ¶
ra807206 rab2aea8 2 2 double Iq(double q, double sld, double solvent_sld, double length_a, 3 3 double b2a_ratio, double c2a_ratio); 4 double Iqxy(double qx, double qy, double sld, double solvent_sld,5 double length_a, double b2a_ratio, double c2a_ratio);6 4 7 5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 8 6 { 9 double b_side= length_a * b2a_ratio;10 double c_side= length_a * c2a_ratio;11 double vol_shell = 2.0 * (length_a* b_side + length_a*c_side + b_side*c_side);7 double length_b = length_a * b2a_ratio; 8 double length_c = length_a * c2a_ratio; 9 double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 12 10 return vol_shell; 13 11 } … … 20 18 double c2a_ratio) 21 19 { 22 double b_side= length_a * b2a_ratio;23 double c_side= length_a * c2a_ratio;24 double a_half = 0.5 * length_a;25 double b_half = 0.5 * b_side;26 double c_half = 0.5 * c_side;20 const double length_b = length_a * b2a_ratio; 21 const double length_c = length_a * c2a_ratio; 22 const double a_half = 0.5 * length_a; 23 const double b_half = 0.5 * length_b; 24 const double c_half = 0.5 * length_c; 27 25 28 26 //Integration limits to use in Gaussian quadrature 29 double v1a = 0.0;30 double v1b = 0.5 * M_PI; //theta integration limits31 double v2a = 0.0;32 double v2b = 0.5 * M_PI; //phi integration limits27 const double v1a = 0.0; 28 const double v1b = M_PI_2; //theta integration limits 29 const double v2a = 0.0; 30 const double v2b = M_PI_2; //phi integration limits 33 31 34 //Order of integration35 int nordi=76;36 int nordj=76;32 double outer_sum = 0.0; 33 for(int i=0; i<76; i++) { 34 const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 37 35 38 double sumi = 0.0; 39 40 for(int i=0; i<nordi; i++) { 36 double sin_theta, cos_theta; 37 double sin_c, cos_c; 38 SINCOS(theta, sin_theta, cos_theta); 39 SINCOS(q*c_half*cos_theta, sin_c, cos_c); 41 40 42 double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );43 44 41 // To check potential problems if denominator goes to zero here !!! 45 double termAL_theta = 8.0*cos(q*c_half*cos(theta)) / (q*q*sin(theta)*sin(theta));46 double termAT_theta = 8.0*sin(q*c_half*cos(theta)) / (q*q*sin(theta)*cos(theta));42 const double termAL_theta = 8.0 * cos_c / (q*q*sin_theta*sin_theta); 43 const double termAT_theta = 8.0 * sin_c / (q*q*sin_theta*cos_theta); 47 44 48 double sumj= 0.0;49 50 for(int j=0; j<nordj; j++) { 45 double inner_sum = 0.0; 46 for(int j=0; j<76; j++) { 47 const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 51 48 52 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 53 49 double sin_phi, cos_phi; 50 double sin_a, cos_a; 51 double sin_b, cos_b; 52 SINCOS(phi, sin_phi, cos_phi); 53 SINCOS(q*a_half*sin_theta*sin_phi, sin_a, cos_a); 54 SINCOS(q*b_half*sin_theta*cos_phi, sin_b, cos_b); 55 54 56 // Amplitude AL from eqn. (7c) 55 double AL = termAL_theta * sin(q*a_half*sin(theta)*sin(phi)) *56 sin(q*b_half*sin(theta)*cos(phi)) / (sin(phi)*cos(phi));57 const double AL = termAL_theta 58 * sin_a*sin_b / (sin_phi*cos_phi); 57 59 58 60 // Amplitude AT from eqn. (9) 59 double AT = termAT_theta * ( (cos(q*a_half*sin(theta)*sin(phi))*sin(q*b_half*sin(theta)*cos(phi))/cos(phi))60 + (cos(q*b_half*sin(theta)*cos(phi))*sin(q*a_half*sin(theta)*sin(phi))/sin(phi)));61 const double AT = termAT_theta 62 * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 61 63 62 sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT); 64 inner_sum += Gauss76Wt[j] * square(AL+AT); 65 } 63 66 64 } 65 66 sumj = 0.5 * (v2b-v2a) * sumj; 67 sumi += Gauss76Wt[i] * sumj * sin(theta); 68 67 inner_sum *= 0.5 * (v2b-v2a); 68 outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 69 69 } 70 70 71 double answer = 0.5*(v1b-v1a)*sumi;71 outer_sum *= 0.5*(v1b-v1a); 72 72 73 73 // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 74 74 // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 75 answer *= (2.0/M_PI);75 double answer = outer_sum/M_PI_2; 76 76 77 77 // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 78 answer *= (sld-solvent_sld)*(sld-solvent_sld);78 answer *= square(sld-solvent_sld); 79 79 80 80 // Convert from [1e-12 A-1] to [cm-1] … … 82 82 83 83 return answer; 84 85 84 } 86 87 double Iqxy(double qx, double qy,88 double sld,89 double solvent_sld,90 double length_a,91 double b2a_ratio,92 double c2a_ratio)93 {94 double q = sqrt(qx*qx + qy*qy);95 double intensity = Iq(q, sld, solvent_sld, length_a, b2a_ratio, c2a_ratio);96 return intensity;97 } -
TabularUnified sasmodels/models/hollow_rectangular_prism_thin_walls.py ¶
ra807206 rab2aea8 3 3 r""" 4 4 5 This model provides the form factor, *P(q)*, for a hollow rectangular5 This model provides the form factor, $P(q)$, for a hollow rectangular 6 6 prism with infinitely thin walls. It computes only the 1D scattering, not the 2D. 7 7 … … 14 14 15 15 Assuming a hollow parallelepiped with infinitely thin walls, edge lengths 16 :math:`A \le B \le C`and presenting an orientation with respect to the17 scattering vector given by |theta| and |phi|, where |theta|is the angle18 between the *z* axis and the longest axis of the parallelepiped *C*, and19 |phi| is the angle between the scattering vector (lying in the *xy*plane)20 and the *y*axis, the form factor is given by16 $A \le B \le C$ and presenting an orientation with respect to the 17 scattering vector given by $\theta$ and $\phi$, where $\theta$ is the angle 18 between the $z$ axis and the longest axis of the parallelepiped $C$, and 19 $\phi$ is the angle between the scattering vector (lying in the $xy$ plane) 20 and the $y$ axis, the form factor is given by 21 21 22 22 .. math:: 23 P(q) = \frac{1}{V^2} \frac{2}{\pi} \int_0^{\frac{\pi}{2}} 24 \int_0^{\frac{\pi}{2}} [A_L(q)+A_T(q)]^2 \sin\theta d\theta d\phi 23 24 P(q) = \frac{1}{V^2} \frac{2}{\pi} \int_0^{\frac{\pi}{2}} 25 \int_0^{\frac{\pi}{2}} [A_L(q)+A_T(q)]^2 \sin\theta\,d\theta\,d\phi 25 26 26 27 where 27 28 28 29 .. math:: 29 V = 2AB + 2AC + 2BC30 30 31 .. math:: 32 A_L(q) = 8 \times \frac{ \sin \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)33 \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr)34 \cos \bigl( q \frac{C}{2} \cos\theta \bigr) }35 {q^2 \, \sin^2\theta \, \sin\phi \cos\phi}36 37 .. math:: 38 A_T(q) = A_F(q) \times \frac{2 \, \sin \bigl( q \frac{C}{2} \cos\theta \bigr)}{q \,\cos\theta}31 V &= 2AB + 2AC + 2BC \\ 32 A_L(q) &= 8 \times \frac{ 33 \sin \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right) 34 \sin \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) 35 \cos \left( \tfrac{1}{2} q C \cos\theta \right) 36 }{q^2 \, \sin^2\theta \, \sin\phi \cos\phi} \\ 37 A_T(q) &= A_F(q) \times 38 \frac{2\,\sin \left( \tfrac{1}{2} q C \cos\theta \right)}{q\,\cos\theta} 39 39 40 40 and 41 41 42 42 .. math:: 43 A_F(q) = 4 \frac{ \cos \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr) 44 \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) } 43 44 A_F(q) = 4 \frac{ \cos \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right) 45 \sin \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) } 45 46 {q \, \cos\phi \, \sin\theta} + 46 4 \frac{ \sin \ bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)47 \cos \ bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) }47 4 \frac{ \sin \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right) 48 \cos \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) } 48 49 {q \, \sin\phi \, \sin\theta} 49 50 … … 51 52 52 53 .. math:: 53 I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} - \rho_{\mbox{solvent}})^2 \times P(q)54 54 55 where *V* is the volume of the rectangular prism, :math:`\rho_{\mbox{p}}` 56 is the scattering length of the parallelepiped, :math:`\rho_{\mbox{solvent}}` 55 I(q) = \text{scale} \times V \times (\rho_\text{p} - \rho_\text{solvent})^2 \times P(q) 56 57 where $V$ is the volume of the rectangular prism, $\rho_\text{p}$ 58 is the scattering length of the parallelepiped, $\rho_\text{solvent}$ 57 59 is the scattering length of the solvent, and (if the data are in absolute 58 60 units) *scale* represents the volume fraction (which is unitless). … … 127 129 # parameters for demo 128 130 demo = dict(scale=1, background=0, 129 sld=6.3 e-6, sld_solvent=1.0e-6,131 sld=6.3, sld_solvent=1.0, 130 132 length_a=35, b2a_ratio=1, c2a_ratio=1, 131 133 length_a_pd=0.1, length_a_pd_n=10, -
TabularUnified sasmodels/models/lamellar_stack_paracrystal.c ¶
r0bef47b r4962519 67 67 double Snq; 68 68 69 Snq = an/( (double)Nlayers* pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) );69 Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) ); 70 70 71 71 return(Snq); -
TabularUnified sasmodels/models/lib/Si.c ¶
r0278e3f rf719764 3 3 double Si(double x) 4 4 { 5 double out; 5 if (x >= M_PI*6.2/4.0) { 6 const double xxinv = 1./(x*x); 7 // Explicitly writing factorial values triples the speed of the calculation 8 const double out_cos = (((-720.*xxinv + 24.)*xxinv - 2.)*xxinv + 1.)/x; 9 const double out_sin = (((-5040.*xxinv + 120.)*xxinv - 6.)*xxinv + 1)*xxinv; 6 10 7 if (x >= M_PI*6.2/4.0){8 double out_sin = 0.0;9 double out_cos = 0.0;10 out = M_PI/2.0;11 11 double sin_x, cos_x; 12 SINCOS(x, sin_x, cos_x); 13 return M_PI_2 - cos_x*out_cos - sin_x*out_sin; 14 } else { 15 const double xx = x*x; 12 16 // Explicitly writing factorial values triples the speed of the calculation 13 out_cos = 1./x - 2./pow(x,3) + 24./pow(x,5) - 720./pow(x,7);14 out_sin = 1./pow(x,2) - 6./pow(x,4) + 120./pow(x,6) - 5040./pow(x,8);15 16 out -= cos(x) * out_cos;17 out -= sin(x) * out_sin;18 return out;17 return (((((-1./439084800.*xx 18 + 1./3265920.)*xx 19 - 1./35280.)*xx 20 + 1./600.)*xx 21 - 1./18.)*xx 22 + 1.)*x; 19 23 } 20 21 // Explicitly writing factorial values triples the speed of the calculation22 out = x - pow(x, 3)/18. + pow(x,5)/600. - pow(x,7)/35280. + pow(x,9)/3265920. - pow(x,11)/439084800.;23 24 return out;25 24 } -
TabularUnified sasmodels/models/lib/core_shell.c ¶
r7d4b2ae r3a48772 17 17 const double core_contrast = core_sld - shell_sld; 18 18 const double core_bes = sph_j1c(core_qr); 19 const double core_volume = 4.0 * M_PI / 3.0 * radius * radius * radius;19 const double core_volume = M_4PI_3 * cube(radius); 20 20 double f = core_volume * core_bes * core_contrast; 21 21 … … 24 24 const double shell_contrast = shell_sld - solvent_sld; 25 25 const double shell_bes = sph_j1c(shell_qr); 26 const double shell_volume = 4.0 * M_PI / 3.0 * pow((radius + thickness), 3);26 const double shell_volume = M_4PI_3 * cube(radius + thickness); 27 27 f += shell_volume * shell_bes * shell_contrast; 28 28 return f * f * 1.0e-4; -
TabularUnified sasmodels/models/lib/gfn.c ¶
r177c1a1 r3a48772 12 12 { 13 13 // local variables 14 const double pi43=4.0/3.0*M_PI;15 14 const double aa = crmaj; 16 15 const double bb = crmin; … … 20 19 //const double siq = (uq == 0.0 ? 1.0 : 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq); 21 20 const double siq = sph_j1c(uq); 22 const double vc = pi43*aa*aa*bb;21 const double vc = M_4PI_3*aa*aa*bb; 23 22 const double gfnc = siq*vc*delpc; 24 23 25 24 const double ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx)); 26 25 const double ut= sqrt(ut2)*qq; 27 const double vt = pi43*trmaj*trmaj*trmin;26 const double vt = M_4PI_3*trmaj*trmaj*trmin; 28 27 //const double sit = (ut == 0.0 ? 1.0 : 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut); 29 28 const double sit = sph_j1c(ut); … … 33 32 const double result = tgfn*tgfn; 34 33 35 return (result);34 return result; 36 35 } -
TabularUnified sasmodels/models/lib/sas_gamma.c ¶
r1596de3 r97f9b46 5 5 We use gamma definition Gamma(t + 1) = t * Gamma(t) to compute 6 6 to function for values lower than 1.0. Namely Gamma(t) = 1/t * Gamma(t + 1) 7 For t < 0, we use Gamma(t) = pi / ( Gamma(1 - t) * sin(pi * t) ) 7 8 */ 8 9 … … 137 138 138 139 139 inline double sas_gamma( double x) { return tgamma(x+1)/x; } 140 inline double sas_gamma(double x) 141 { 142 // Note: the builtin tgamma can give slow and unreliable results for x<1. 143 // The following transform extends it to zero and to negative values. 144 // It should return NaN for zero and negative integers but doesn't. 145 // The accuracy is okay but not wonderful for negative numbers, maybe 146 // one or two digits lost in the calculation. If higher accuracy is 147 // needed, you could test the following loop: 148 // double norm = 1.; 149 // while (x<1.) { norm*=x; x+=1.; } 150 // return tgamma(x)/norm; 151 return (x<0. ? M_PI/tgamma(1.-x)/sin(M_PI*x) : tgamma(x+1)/x); 152 } -
TabularUnified sasmodels/models/lib/wrc_cyl.c ¶
rba32cdd r92ce163 14 14 //return t; 15 15 16 return pow( (1.0 + (x/3.12)*(x/3.12) + 17 (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 16 return pow(1.0+square(x/3.12)+cube(x/8.67), 0.176/3.0); 18 17 } 19 18 … … 23 22 { 24 23 const double r = b/L; 25 return (L*b/6.0) * 26 (1.0 - r*1.5 + 1.5*r*r - 0.75*r*r*r*(1.0 - exp(-2.0/r))); 24 return (L*b/6.0) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r)))); 25 27 26 } 28 27 … … 41 40 } 42 41 43 static inlinedouble42 static double 44 43 sech_WR(double x) 45 44 { … … 51 50 { 52 51 double C; 53 const double onehalf = 1.0/2.0;54 52 55 53 if( L/b > 10.0) { … … 86 84 87 85 const double t2 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 88 Rg02/b2))*((1.0 + onehalf*(((-1.0) -86 Rg02/b2))*((1.0 + 0.5*(((-1.0) - 89 87 tanh((-C4 + Rgb/C5))))))); 90 88 … … 112 110 113 111 const double t9 = (2.0*b4*(((2.0*q0*Rg2)/b - 114 (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + onehalf*(((-1.0) -112 (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + 0.5*(((-1.0) - 115 113 tanh(((-C4) + Rgb)/C5)))))); 116 114 117 115 const double t10 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 118 Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) +116 Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 119 117 Rgb)/C5)))))); 120 118 … … 131 129 132 130 const double t14 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 133 Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) +131 Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 134 132 Rgb)/C5)))))); 135 133 … … 140 138 141 139 double yy = (pow(q0,p1)*(((-((b*M_PI)/(L*q0))) +t1/L +t2/(q04*Rg22) + 142 onehalf*t3*t4)) + (t5*((pow(q0,(p1 - p2))*140 0.5*t3*t4)) + (t5*((pow(q0,(p1 - p2))* 143 141 (((-pow(q0,(-p1)))*(((b2*M_PI)/(L*q02) +t6/L +t7/(2.0*C5) - 144 t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + onehalf*t11*t12)) -142 t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 0.5*t11*t12)) - 145 143 b*p1*pow(q0,((-1.0) - p1))*(((-((b*M_PI)/(L*q0))) + t13/L + 146 t14/(q04*Rg22) + onehalf*t15*((1.0 + tanh(((-C4) +144 t14/(q04*Rg22) + 0.5*t15*((1.0 + tanh(((-C4) + 147 145 Rgb)/C5))))))))))); 148 146 … … 154 152 { 155 153 double C; 156 const double onehalf = 1.0/2.0;157 154 158 155 if( L/b > 10.0) { … … 201 198 const double t5 = (2.0*b4*(((2.0*q0*Rg2)/b - 202 199 (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))* 203 ((1.0 + onehalf*(((-1.0) - tanh(((-C4) +200 ((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 204 201 Rgb)/C5))))))/(q04*Rg22); 205 202 206 203 const double t6 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 207 Rg02/b2))*((1.0 + onehalf*(((-1) - tanh(((-C4) +204 Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 208 205 Rgb)/C5))))))/(q05*Rg22); 209 206 … … 219 216 220 217 const double t10 = (2.0*b4*(((-1) + pow((double)M_E,(-(Rg02/b2))) + 221 Rg02/b2))*((1.0 + onehalf*(((-1) - tanh(((-C4) +218 Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 222 219 Rgb)/C5))))))/(q04*Rg22); 223 220 224 221 const double yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*M_PI)/(L*q02) + 225 t2 + t3 - t4 + t5 - t6 + onehalf*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*222 t2 + t3 - t4 + t5 - t6 + 0.5*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))* 226 223 (((-((b*M_PI)/(L*q0))) + t9 + t10 + 227 onehalf*((C3*pow(((Rgb)),((-3.0)/miu)) +224 0.5*((C3*pow(((Rgb)),((-3.0)/miu)) + 228 225 C2*pow(((Rgb)),((-2.0)/miu)) + 229 226 C1*pow(((Rgb)),((-1.0)/miu))))* -
TabularUnified sasmodels/models/linear_pearls.c ¶
r2c74c11 r4962519 19 19 { 20 20 // Pearl volume 21 double pearl_vol = 4.0 /3.0 * M_PI * pow(radius, 3.0);21 double pearl_vol = M_4PI_3 * cube(radius); 22 22 // Return total volume 23 23 return num_pearls * pearl_vol;; … … 35 35 double contrast_pearl = pearl_sld - solvent_sld; 36 36 //each volume 37 double pearl_vol = 4.0 /3.0 * M_PI * pow(radius, 3.0);37 double pearl_vol = M_4PI_3 * cube(radius); 38 38 //total volume 39 39 double tot_vol = num_pearls * pearl_vol; … … 43 43 double separation = edge_sep + 2.0 * radius; 44 44 45 double x=q*radius;46 47 // Try Taylor on x*xos(x)48 // double out_cos = x - pow(x,3)/2 + pow(x,5)/24 - pow(x,7)/720 + pow(x,9)/40320;49 // psi -= x*out_cos;50 51 45 //sine functions of a pearl 52 double psi = sin(q * radius); 53 psi -= x * cos(x); 54 psi /= pow((q * radius), 3.0); 46 double psi = sph_j1c(q * radius); 55 47 56 48 // N pearls contribution … … 61 53 } 62 54 // form factor for num_pearls 63 double form_factor = n_contrib; 64 form_factor *= pow((m_s*psi*3.0), 2.0); 65 form_factor /= (tot_vol * 1.0e4); 55 double form_factor = 1.0e-4 * n_contrib * square(m_s*psi) / tot_vol; 66 56 67 57 return form_factor; -
TabularUnified sasmodels/models/linear_pearls.py ¶
r40a87fa r4962519 63 63 single = False 64 64 65 source = ["li near_pearls.c"]65 source = ["lib/sph_j1c.c", "linear_pearls.c"] 66 66 67 67 demo = dict(scale=1.0, background=0.0, -
TabularUnified sasmodels/models/mass_fractal.c ¶
ra807206 r6d96b66 1 double form_volume(double radius); 1 #define INVALID(p) (p.fractal_dim_mass < 1.0) 2 2 3 double Iq(double q, 4 double radius, 5 double fractal_dim_mass, 6 double cutoff_length); 3 static double 4 Iq(double q, double radius, double fractal_dim_mass, double cutoff_length) 5 { 6 //calculate P(q) 7 const double pq = square(sph_j1c(q*radius)); 7 8 8 static double _mass_fractal_kernel(double q, 9 double radius, 10 double fractal_dim_mass, 11 double cutoff_length) 12 { 13 // Actively check the argument. 14 if (fractal_dim_mass <= 1.0){ 15 return 0.0; 9 //calculate S(q) 10 // S(q) = gamma(D-1) sin((D-1)atan(q c))/q c^(D-1) (1+(q c)^2)^(-(D-1)/2) 11 // lim D->1 [gamma(D-1) sin((D-1) a)] = a 12 // lim q->0 [sin((D-1) atan(q c))/q] = (D-1) c 13 // lim q,D->0,1 [gamma(D-1) sin((D-1)atan(q c))/q] = c 14 double sq; 15 if (q > 0. && fractal_dim_mass > 1.) { 16 const double Dm1 = fractal_dim_mass - 1.0; 17 const double t1 = sas_gamma(Dm1)*sin(Dm1*atan(q*cutoff_length)); 18 const double t2 = pow(cutoff_length, Dm1); 19 const double t3 = pow(1.0 + square(q*cutoff_length), -0.5*Dm1); 20 sq = t1 * t2 * t3 / q; 21 } else if (q > 0.) { 22 sq = atan(q*cutoff_length)/q; 23 } else if (fractal_dim_mass > 1.) { 24 const double D = fractal_dim_mass; 25 sq = pow(cutoff_length, D) * sas_gamma(D); 26 } else { 27 sq = cutoff_length; 16 28 } 17 29 18 //calculate P(q) 19 double pq = sph_j1c(q*radius); 20 pq = pq*pq; 21 22 //calculate S(q) 23 double mmo = fractal_dim_mass-1.0; 24 double sq = sas_gamma(mmo)*sin((mmo)*atan(q*cutoff_length)); 25 sq *= pow(cutoff_length, mmo); 26 sq /= pow((1.0 + (q*cutoff_length)*(q*cutoff_length)),(mmo/2.0)); 27 sq /= q; 28 29 //combine and return 30 double result = pq * sq; 31 32 return result; 30 return pq * sq; 33 31 } 34 double form_volume(double radius){35 36 return 1.333333333333333*M_PI*radius*radius*radius;37 }38 39 double Iq(double q,40 double radius,41 double fractal_dim_mass,42 double cutoff_length)43 {44 return _mass_fractal_kernel(q,45 radius,46 fractal_dim_mass,47 cutoff_length);48 } -
TabularUnified sasmodels/models/mass_fractal.py ¶
ra807206 r6d96b66 78 78 79 79 # pylint: disable=bad-whitespace, line-too-long 80 # ["name", "units", default, [lower, upper], "type","description"], 81 parameters = [["radius", "Ang", 10.0, [0.0, inf], "", "Particle radius"], 82 ["fractal_dim_mass", "", 1.9, [1.0, 6.0], "", "Mass fractal dimension"], 83 ["cutoff_length", "Ang", 100.0, [0.0, inf], "", "Cut-off length"], 84 ] 80 # ["name", "units", default, [lower, upper], "type","description"], 81 parameters = [ 82 ["radius", "Ang", 10.0, [0.0, inf], "", "Particle radius"], 83 ["fractal_dim_mass", "", 1.9, [1.0, 6.0], "", "Mass fractal dimension"], 84 ["cutoff_length", "Ang", 100.0, [0.0, inf], "", "Cut-off length"], 85 ] 85 86 # pylint: enable=bad-whitespace, line-too-long 86 87 -
TabularUnified sasmodels/models/mass_surface_fractal.c ¶
r30fbe2e r6d96b66 1 double form_volume(double radius); 2 3 double Iq(double q, 4 double fractal_dim_mass, 5 double fractal_dim_surf, 6 double rg_cluster, 7 double rg_primary); 8 9 static double _mass_surface_fractal_kernel(double q, 1 #define INVALID(p) (p.fractal_dim_mass + p.fractal_dim_surf > 6.) 2 static double 3 Iq(double q, 10 4 double fractal_dim_mass, 11 5 double fractal_dim_surf, … … 14 8 { 15 9 //computation 16 double tot_dim = 6.0 - fractal_dim_surf - fractal_dim_mass; 17 fractal_dim_mass /= 2.0; 18 tot_dim /= 2.0; 10 const double Dm = 0.5*fractal_dim_mass; 11 const double Dt = 0.5*(6.0 - (fractal_dim_mass + fractal_dim_surf)); 19 12 20 double rc_norm = rg_cluster * rg_cluster / (3.0 * fractal_dim_mass); 21 double rp_norm = rg_primary * rg_primary / (3.0 * tot_dim); 13 const double t1 = Dm==0. ? 1.0 : pow(1.0 + square(q*rg_cluster)/(3.0*Dm), -Dm); 14 const double t2 = Dt==0. ? 1.0 : pow(1.0 + square(q*rg_primary)/(3.0*Dt), -Dt); 15 const double form = t1*t2; 22 16 23 //x for P 24 double x_val1 = 1.0 + q * q * rc_norm; 25 double x_val2 = 1.0 + q * q * rp_norm; 26 27 double inv_form = pow(x_val1, fractal_dim_mass) * pow(x_val2, tot_dim); 28 29 //another singular 30 if (inv_form == 0.0) return 0.0; 31 32 double form_factor = 1.0; 33 form_factor /= inv_form; 34 35 return (form_factor); 17 return form; 36 18 } 37 double form_volume(double radius){38 39 return 1.333333333333333*M_PI*radius*radius*radius;40 }41 42 double Iq(double q,43 double fractal_dim_mass,44 double fractal_dim_surf,45 double rg_cluster,46 double rg_primary)47 {48 return _mass_surface_fractal_kernel(q,49 fractal_dim_mass,50 fractal_dim_surf,51 rg_cluster,52 rg_primary);53 } -
TabularUnified sasmodels/models/mass_surface_fractal.py ¶
r30fbe2e r6d96b66 78 78 79 79 # pylint: disable=bad-whitespace, line-too-long 80 # ["name", "units", default, [lower, upper], "type","description"], 81 parameters = [["fractal_dim_mass", "", 1.8, [1e-16, 6.0], "", 82 "Mass fractal dimension"], 83 ["fractal_dim_surf", "", 2.3, [1e-16, 6.0], "", 84 "Surface fractal dimension"], 85 ["rg_cluster", "Ang", 86.7, [0.0, inf], "", 86 "Cluster radius of gyration"], 87 ["rg_primary", "Ang", 4000., [0.0, inf], "", 88 "Primary particle radius of gyration"], 89 ] 80 # ["name", "units", default, [lower, upper], "type","description"], 81 parameters = [ 82 ["fractal_dim_mass", "", 1.8, [0.0, 6.0], "", "Mass fractal dimension"], 83 ["fractal_dim_surf", "", 2.3, [0.0, 6.0], "", "Surface fractal dimension"], 84 ["rg_cluster", "Ang", 86.7, [0.0, inf], "", "Cluster radius of gyration"], 85 ["rg_primary", "Ang", 4000., [0.0, inf], "", "Primary particle radius of gyration"], 86 ] 90 87 # pylint: enable=bad-whitespace, line-too-long 91 88 -
TabularUnified sasmodels/models/multilayer_vesicle.c ¶
r2c74c11 r3a48772 19 19 20 20 // layer 1 21 voli = 4.0*M_PI/3.0*ri*ri*ri;21 voli = M_4PI_3*ri*ri*ri; 22 22 fval += voli*sldi*sph_j1c(ri*q); 23 23 … … 25 25 26 26 // layer 2 27 voli = 4.0*M_PI/3.0*ri*ri*ri;27 voli = M_4PI_3*ri*ri*ri; 28 28 fval -= voli*sldi*sph_j1c(ri*q); 29 29 -
TabularUnified sasmodels/models/onion.c ¶
rd119f34 r9762341 30 30 31 31 static double 32 form_volume(double core_radius, double n, double thickness[])32 form_volume(double radius_core, double n, double thickness[]) 33 33 { 34 34 int i; 35 double r = core_radius;35 double r = radius_core; 36 36 for (i=0; i < n; i++) { 37 37 r += thickness[i]; … … 41 41 42 42 static double 43 Iq(double q, double sld_core, double core_radius, double sld_solvent,43 Iq(double q, double sld_core, double radius_core, double sld_solvent, 44 44 double n_shells, double sld_in[], double sld_out[], double thickness[], 45 45 double A[]) 46 46 { 47 47 int n = (int)(n_shells+0.5); 48 double r_out = core_radius;48 double r_out = radius_core; 49 49 double f = f_exp(q, r_out, sld_core, 0.0, 0.0, 0.0, 0.0); 50 50 for (int i=0; i < n; i++){ -
TabularUnified sasmodels/models/onion.py ¶
r7b68dc5 r9762341 366 366 return np.asarray(z), np.asarray(rho) 367 367 368 def ER( core_radius, n, thickness):368 def ER(radius_core, n, thickness): 369 369 """Effective radius""" 370 return np.sum(thickness[:int(n[0])], axis=0) + core_radius370 return np.sum(thickness[:int(n[0])], axis=0) + radius_core 371 371 372 372 demo = { 373 373 "sld_solvent": 2.2, 374 374 "sld_core": 1.0, 375 " core_radius": 100,375 "radius_core": 100, 376 376 "n_shells": 4, 377 377 "sld_in": [0.5, 1.5, 0.9, 2.0], … … 381 381 # Could also specify them individually as 382 382 # "A1": 0, "A2": -1, "A3": 1e-4, "A4": 1, 383 #" core_radius_pd_n": 10,384 #" core_radius_pd": 0.4,383 #"radius_core_pd_n": 10, 384 #"radius_core_pd": 0.4, 385 385 #"thickness4_pd_n": 10, 386 386 #"thickness4_pd": 0.4, -
TabularUnified sasmodels/models/parallelepiped.c ¶
ra807206 r14838a3 1 1 double form_volume(double length_a, double length_b, double length_c); 2 double Iq(double q, double sld, double solvent_sld, double length_a, double length_b, double length_c); 2 double Iq(double q, double sld, double solvent_sld, 3 double length_a, double length_b, double length_c); 3 4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 4 double length_a, double length_b, double length_c, double theta, double phi, double psi); 5 6 // From Igor library 7 double _pkernel(double a, double b,double c, double ala, double alb, double alc); 8 double _pkernel(double a, double b,double c, double ala, double alb, double alc){ 9 double argA,argB,argC,tmp1,tmp2,tmp3; 10 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 11 argA = 0.5*a*ala; 12 argB = 0.5*b*alb; 13 argC = 0.5*c*alc; 14 if(argA==0.0) { 15 tmp1 = 1.0; 16 } else { 17 tmp1 = sin(argA)*sin(argA)/argA/argA; 18 } 19 if (argB==0.0) { 20 tmp2 = 1.0; 21 } else { 22 tmp2 = sin(argB)*sin(argB)/argB/argB; 23 } 24 if (argC==0.0) { 25 tmp3 = 1.0; 26 } else { 27 tmp3 = sin(argC)*sin(argC)/argC/argC; 28 } 29 return (tmp1*tmp2*tmp3); 30 31 } 32 5 double length_a, double length_b, double length_c, 6 double theta, double phi, double psi); 33 7 34 8 double form_volume(double length_a, double length_b, double length_c) … … 45 19 double length_c) 46 20 { 47 double tmp1, tmp2; 48 49 double mu = q * length_b; 21 const double mu = 0.5 * q * length_b; 50 22 51 23 // Scale sides by B 52 double a_scaled = length_a / length_b;53 double c_scaled = length_c / length_b;24 const double a_scaled = length_a / length_b; 25 const double c_scaled = length_c / length_b; 54 26 55 //Order of integration 56 int nordi=76; 57 int nordj=76; 27 // outer integral (with gauss points), integration limits = 0, 1 28 double outer_total = 0; //initialize integral 58 29 59 // outer integral (with nordi gauss points), integration limits = 0, 1 60 double summ = 0; //initialize integral 30 for( int i=0; i<76; i++) { 31 const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 32 const double mu_proj = mu * sqrt(1.0-sigma*sigma); 61 33 62 for( int i=0; i<nordi; i++) { 63 64 // inner integral (with nordj gauss points), integration limits = 0, 1 65 66 double summj = 0.0; 67 double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 68 69 for(int j=0; j<nordj; j++) { 34 // inner integral (with gauss points), integration limits = 0, 1 35 // corresponding to angles from 0 to pi/2. 36 double inner_total = 0.0; 37 for(int j=0; j<76; j++) { 38 const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 39 double sin_uu, cos_uu; 40 SINCOS(M_PI_2*uu, sin_uu, cos_uu); 41 const double si1 = sinc(mu_proj * sin_uu * a_scaled); 42 const double si2 = sinc(mu_proj * cos_uu); 43 inner_total += Gauss76Wt[j] * square(si1 * si2); 44 } 45 inner_total *= 0.5; 70 46 71 double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 72 double mudum = mu * sqrt(1.0-sigma*sigma); 73 double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu); 74 double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu); 75 if(arg1==0.0) { 76 tmp1 = 1.0; 77 } else { 78 tmp1 = sin(arg1)*sin(arg1)/arg1/arg1; 79 } 80 if (arg2==0.0) { 81 tmp2 = 1.0; 82 } else { 83 tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; 84 } 47 const double si = sinc(mu * c_scaled * sigma); 48 outer_total += Gauss76Wt[i] * inner_total * si * si; 49 } 50 outer_total *= 0.5; 85 51 86 summj += Gauss76Wt[j] * tmp1 * tmp2; 87 } 88 89 // value of the inner integral 90 double answer = 0.5 * summj; 91 92 double arg = 0.5 * mu * c_scaled * sigma; 93 if ( arg == 0.0 ) { 94 answer *= 1.0; 95 } else { 96 answer *= sin(arg)*sin(arg)/arg/arg; 97 } 98 99 // sum of outer integral 100 summ += Gauss76Wt[i] * answer; 101 102 } 103 104 const double vd = (sld-solvent_sld) * form_volume(length_a, length_b, length_c); 105 106 // convert from [1e-12 A-1] to [cm-1] and 0.5 factor for outer integral 107 return 1.0e-4 * 0.5 * vd * vd * summ; 108 52 // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1] 53 const double V = form_volume(length_a, length_b, length_c); 54 const double drho = (sld-solvent_sld); 55 return 1.0e-4 * square(drho * V) * outer_total; 109 56 } 110 57 … … 120 67 double psi) 121 68 { 122 double q = sqrt(qx*qx+qy*qy); 123 double qx_scaled = qx/q; 124 double qy_scaled = qy/q; 69 double q, cos_val_a, cos_val_b, cos_val_c; 70 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 125 71 126 // Convert angles given in degrees to radians 127 theta *= M_PI_180; 128 phi *= M_PI_180; 129 psi *= M_PI_180; 130 131 // Parallelepiped c axis orientation 132 double cparallel_x = cos(theta) * cos(phi); 133 double cparallel_y = sin(theta); 134 135 // Compute angle between q and parallelepiped axis 136 double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz; 137 138 // Parallelepiped a axis orientation 139 double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 140 double parallel_y = sin(psi)*cos(theta); 141 double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled; 142 143 // Parallelepiped b axis orientation 144 double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 145 double bparallel_y = cos(theta)*cos(psi); 146 double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled; 147 148 // The following tests should always pass 149 if (fabs(cos_val_c)>1.0) { 150 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 151 cos_val_c = 1.0; 152 } 153 if (fabs(cos_val_a)>1.0) { 154 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 155 cos_val_a = 1.0; 156 } 157 if (fabs(cos_val_b)>1.0) { 158 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 159 cos_val_b = 1.0; 160 } 161 162 // Call the IGOR library function to get the kernel 163 double form = _pkernel( q*length_a, q*length_b, q*length_c, cos_val_a, cos_val_b, cos_val_c); 164 165 // Multiply by contrast^2 166 const double vd = (sld - solvent_sld) * form_volume(length_a, length_b, length_c); 167 return 1.0e-4 * vd * vd * form; 72 const double siA = sinc(0.5*q*length_a*cos_val_a); 73 const double siB = sinc(0.5*q*length_b*cos_val_b); 74 const double siC = sinc(0.5*q*length_c*cos_val_c); 75 const double V = form_volume(length_a, length_b, length_c); 76 const double drho = (sld - solvent_sld); 77 const double form = V * drho * siA * siB * siC; 78 // Square and convert from [1e-12 A-1] to [cm-1] 79 return 1.0e-4 * form * form; 168 80 } -
TabularUnified sasmodels/models/parallelepiped.py ¶
r416f5c7 red0827a 221 221 # parameters for demo 222 222 demo = dict(scale=1, background=0, 223 sld=6.3 e-6, sld_solvent=1.0e-6,223 sld=6.3, sld_solvent=1.0, 224 224 length_a=35, length_b=75, length_c=400, 225 225 theta=45, phi=30, psi=15, -
TabularUnified sasmodels/models/pearl_necklace.c ¶
ra807206 r2126131 1 double _pearl_necklace_kernel(double q, double radius, double edge_sep,2 double thick_string, double num_pearls, double sld_pearl,3 double sld_string, double sld_solv);4 1 double form_volume(double radius, double edge_sep, 5 double thick_string, double num_pearls); 6 2 double thick_string, double num_pearls); 7 3 double Iq(double q, double radius, double edge_sep, 8 9 4 double thick_string, double num_pearls, double sld, 5 double string_sld, double solvent_sld); 10 6 11 7 #define INVALID(v) (v.thick_string >= v.radius || v.num_pearls <= 0) 12 8 13 9 // From Igor library 14 double _pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string, 15 double num_pearls, double sld_pearl, double sld_string, double sld_solv) 10 static double 11 _pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string, 12 double num_pearls, double sld_pearl, double sld_string, double sld_solv) 16 13 { 17 //relative slds 18 double contrast_pearl = sld_pearl - sld_solv; 19 double contrast_string = sld_string - sld_solv; 20 21 // number of string segments 22 num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 23 double num_strings = num_pearls - 1.0; 24 25 //Pi 26 double pi = 4.0*atan(1.0); 27 28 // center to center distance between the neighboring pearls 29 double A_s = edge_sep + 2.0 * radius; 30 31 // Repeated Calculations 32 double sincasq = sinc(q*A_s); 33 double oneminussinc = 1 - sincasq; 34 double q_r = q * radius; 35 double q_edge = q * edge_sep; 36 37 // each volume 38 double string_vol = edge_sep * pi * thick_string * thick_string / 4.0; 39 double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius; 14 // number of string segments 15 num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 16 const double num_strings = num_pearls - 1.0; 40 17 41 //total volume 42 double tot_vol; 43 //each masses 44 double m_r= contrast_string * string_vol; 45 double m_s= contrast_pearl * pearl_vol; 46 double psi, gamma, beta; 47 //form factors 48 double sss, srr, srs; //cross 49 double srr_1, srr_2, srr_3; 50 double form_factor; 51 tot_vol = num_strings * string_vol; 52 tot_vol += num_pearls * pearl_vol; 18 //each masses: contrast * volume 19 const double contrast_pearl = sld_pearl - sld_solv; 20 const double contrast_string = sld_string - sld_solv; 21 const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; 22 const double pearl_vol = M_4PI_3 * radius * radius * radius; 23 const double m_string = contrast_string * string_vol; 24 const double m_pearl = contrast_pearl * pearl_vol; 53 25 54 //sine functions of a pearl 55 psi = sin(q_r); 56 psi -= q_r * cos(q_r); 57 psi *= 3.0; 58 psi /= q_r * q_r * q_r; 26 // center to center distance between the neighboring pearls 27 const double A_s = edge_sep + 2.0 * radius; 59 28 60 // Note take only 20 terms in Si series: 10 terms may be enough though. 61 gamma = Si(q_edge); 62 gamma /= (q_edge); 63 beta = Si(q * (A_s - radius)); 64 beta -= Si(q_r); 65 beta /= q_edge; 29 //sine functions of a pearl 30 // Note: lim_(q->0) Si(q*a)/(q*b) = a/b 31 // So therefore: 32 // beta = q==0. ? 1.0 : (Si(q*(A_s-radius)) - Si(q*radius))/q_edge; 33 // gamma = q==0. ? 1.0 : Si(q_edge)/q_edge; 34 // But there is a 1/(1-sinc) term below which blows up so don't bother 35 const double q_edge = q * edge_sep; 36 const double beta = (Si(q*(A_s-radius)) - Si(q*radius)) / q_edge; 37 const double gamma = Si(q_edge) / q_edge; 38 const double psi = sph_j1c(q*radius); 66 39 67 // form factor for num_pearls 68 sss = 1.0 - pow(sincasq, num_pearls); 69 sss /= oneminussinc * oneminussinc; 70 sss *= -sincasq; 71 sss -= num_pearls / 2.0; 72 sss += num_pearls / oneminussinc; 73 sss *= 2.0 * m_s * psi * m_s * psi; 40 // Precomputed sinc terms 41 const double si = sinc(q*A_s); 42 const double omsi = 1.0 - si; 43 const double pow_si = pow(si, num_pearls); 74 44 75 // form factor for num_strings (like thin rods) 76 srr_1 = -sinc(q_edge/2.0) * sinc(q_edge/2.0); 45 // form factor for num_pearls 46 const double sss = 2.0*square(m_pearl*psi) * ( 47 - si * (1.0 - pow_si) / (omsi*omsi) 48 + num_pearls / omsi 49 - 0.5 * num_pearls 50 ); 77 51 78 srr_1 += 2.0 * gamma; 79 srr_1 *= num_strings; 80 srr_2 = 2.0/oneminussinc; 81 srr_2 *= num_strings; 82 srr_2 *= beta * beta; 83 srr_3 = 1.0 - pow(sincasq, num_strings); 84 srr_3 /= oneminussinc * oneminussinc; 85 srr_3 *= beta * beta; 86 srr_3 *= -2.0; 52 // form factor for num_strings (like thin rods) 53 const double srr = m_string * m_string * ( 54 - 2.0 * (1.0 - pow_si/si)*beta*beta / (omsi*omsi) 55 + 2.0 * num_strings*beta*beta / omsi 56 + num_strings * (2.0*gamma - square(sinc(q_edge/2.0))) 57 ); 87 58 88 // total srr 89 srr = srr_1 + srr_2 + srr_3; 90 srr *= m_r * m_r; 59 // form factor for correlations 60 const double srs = 4.0 * m_string * m_pearl * beta * psi * ( 61 - si * (1.0 - pow_si/si) / (omsi*omsi) 62 + num_strings / omsi 63 ); 91 64 92 // form factor for correlations 93 srs = 1.0; 94 srs -= pow(sincasq, num_strings); 95 srs /= oneminussinc * oneminussinc; 96 srs *= -sincasq; 97 srs += num_strings/oneminussinc; 98 srs *= 4.0; 99 srs *= (m_r * m_s * beta * psi); 65 const double form = sss + srr + srs; 100 66 101 form_factor = sss + srr + srs; 102 form_factor /= (tot_vol * 1.0e4); // norm by volume and A^-1 to cm^-1 103 104 return (form_factor); 67 return 1.0e-4 * form; 105 68 } 106 69 107 70 double form_volume(double radius, double edge_sep, 108 71 double thick_string, double num_pearls) 109 72 { 110 double total_vol; 73 num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 111 74 112 double pi = 4.0*atan(1.0); 113 double number_of_strings = num_pearls - 1.0; 114 115 double string_vol = edge_sep * pi * thick_string * thick_string / 4.0; 116 double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius; 75 const double num_strings = num_pearls - 1.0; 76 const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; 77 const double pearl_vol = M_4PI_3 * radius * radius * radius; 78 const double volume = num_strings*string_vol + num_pearls*pearl_vol; 117 79 118 total_vol = number_of_strings * string_vol; 119 total_vol += num_pearls * pearl_vol; 120 121 return(total_vol); 80 return volume; 122 81 } 123 82 124 83 double Iq(double q, double radius, double edge_sep, 125 126 84 double thick_string, double num_pearls, double sld, 85 double string_sld, double solvent_sld) 127 86 { 128 double value, tot_vol; 129 130 value = _pearl_necklace_kernel(q, radius, edge_sep, thick_string, 131 num_pearls, sld, string_sld, solvent_sld); 132 tot_vol = form_volume(radius, edge_sep, thick_string, num_pearls); 87 const double form = _pearl_necklace_kernel(q, radius, edge_sep, 88 thick_string, num_pearls, sld, string_sld, solvent_sld); 133 89 134 return value*tot_vol;90 return form; 135 91 } -
TabularUnified sasmodels/models/pearl_necklace.py ¶
ra807206 r2126131 92 92 ] 93 93 94 source = ["lib/Si.c", " pearl_necklace.c"]94 source = ["lib/Si.c", "lib/sph_j1c.c", "pearl_necklace.c"] 95 95 single = False # use double precision unless told otherwise 96 96 … … 112 112 """ 113 113 tot_vol = volume(radius, edge_sep, thick_string, num_pearls) 114 rad_out = pow((3.0*tot_vol/4.0/pi), 0.33333)114 rad_out = (tot_vol/(4.0/3.0*pi)) ** (1./3.) 115 115 return rad_out 116 116 -
TabularUnified sasmodels/models/poly_gauss_coil.py ¶
ra807206 rf0afad2 59 59 """ 60 60 61 from numpy import inf, exp, power 61 import numpy as np 62 from numpy import inf, expm1, power 62 63 63 64 name = "poly_gauss_coil" … … 83 84 # pylint: disable = missing-docstring 84 85 u = polydispersity - 1.0 85 z = (q*rg)**2 / (1.0 + 2.0*u) 86 z = q**2 * (rg**2 / (1.0 + 2.0*u)) 87 86 88 # need to trap the case of the polydispersity being 1 (ie, monodisperse!) 87 89 if polydispersity == 1.0: 88 inten = i_zero * 2.0 * (exp(-z) + z - 1.0) 90 result = 2.0 * (expm1(-z) + z) 91 index = q != 0. 92 result[index] /= z[index]**2 93 result[~index] = 1.0 89 94 else: 90 inten = i_zero * 2.0 * (power(1.0 + u*z, -1.0/u) + z - 1.0) / (1.0 + u) 91 index = q != 0. 92 inten[~index] = i_zero 93 inten[index] /= z[index]**2 94 return inten 95 # Taylor series around z=0 of (2*(1+uz)^(-1/u) + z - 1) / (z^2(u+1)) 96 p = [ 97 #(-1 - 20*u - 155*u**2 - 580*u**3 - 1044*u**4 - 720*u**5) / 2520., 98 #( 1 + 14*u + 71*u**2 + 154*u**3 + 120*u**4) / 360., 99 #(-1 - 9*u - 26*u**2 - 24*u**3) / 60., 100 ( 1 + 5*u + 6*u**2) / 12., 101 (-1 - 2*u) / 3., 102 ( 1 ), 103 ] 104 result = 2.0 * (power(1.0 + u*z, -1.0/u) + z - 1.0) / (1.0 + u) 105 index = z > 1e-4 106 result[index] /= z[index]**2 107 result[~index] = np.polyval(p, z[~index]) 108 return i_zero * result 95 109 Iq.vectorized = True # Iq accepts an array of q values 96 110 -
TabularUnified sasmodels/models/polymer_micelle.c ¶
ra807206 rc3ebc71 32 32 // Self-correlation term of the core 33 33 const double bes_core = sph_j1c(q*radius_core); 34 const double term1 = n_aggreg*n_aggreg*beta_core*beta_core*bes_core*bes_core;34 const double term1 = square(n_aggreg*beta_core*bes_core); 35 35 36 36 // Self-correlation term of the chains 37 const double qrg2 = q*rg*q*rg;37 const double qrg2 = square(q*rg); 38 38 const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(expm1(-qrg2)+qrg2)/(qrg2*qrg2); 39 39 const double term2 = n_aggreg * beta_corona * beta_corona * debye_chain; … … 42 42 const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2; 43 43 const double bes_corona = sinc(q*(radius_core + d_penetration * rg)); 44 const double term3 = 2 * n_aggreg * n_aggreg * beta_core * beta_corona *44 const double term3 = 2.0 * n_aggreg * n_aggreg * beta_core * beta_corona * 45 45 bes_core * chain_ampl * bes_corona; 46 46 47 47 // Interference cross-term between chains 48 const double term4 = n_aggreg * (n_aggreg - 1.0) * beta_corona * beta_corona *49 chain_ampl * chain_ampl * bes_corona * bes_corona;48 const double term4 = n_aggreg * (n_aggreg - 1.0) 49 * square(beta_corona * chain_ampl * bes_corona); 50 50 51 51 // I(q)_micelle : Sum of 4 terms computed above -
TabularUnified sasmodels/models/porod.py ¶
r40a87fa r4962519 41 41 """ 42 42 with errstate(divide='ignore'): 43 return power(q, -4)43 return q**-4 44 44 45 45 Iq.vectorized = True # Iq accepts an array of q values -
TabularUnified sasmodels/models/rectangular_prism.c ¶
ra807206 rab2aea8 2 2 double Iq(double q, double sld, double solvent_sld, double length_a, 3 3 double b2a_ratio, double c2a_ratio); 4 double Iqxy(double qx, double qy, double sld, double solvent_sld,5 double length_a, double b2a_ratio, double c2a_ratio);6 4 7 5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio) … … 17 15 double c2a_ratio) 18 16 { 19 double termA, termB, termC; 20 21 double b_side = length_a * b2a_ratio; 22 double c_side = length_a * c2a_ratio; 23 double volume = length_a * b_side * c_side; 24 double a_half = 0.5 * length_a; 25 double b_half = 0.5 * b_side; 26 double c_half = 0.5 * c_side; 17 const double length_b = length_a * b2a_ratio; 18 const double length_c = length_a * c2a_ratio; 19 const double a_half = 0.5 * length_a; 20 const double b_half = 0.5 * length_b; 21 const double c_half = 0.5 * length_c; 27 22 28 23 //Integration limits to use in Gaussian quadrature 29 double v1a = 0.0;30 double v1b = 0.5 * M_PI; //theta integration limits31 double v2a = 0.0;32 double v2b = 0.5 * M_PI; //phi integration limits24 const double v1a = 0.0; 25 const double v1b = M_PI_2; //theta integration limits 26 const double v2a = 0.0; 27 const double v2b = M_PI_2; //phi integration limits 33 28 34 //Order of integration 35 int nordi=76; 36 int nordj=76; 29 double outer_sum = 0.0; 30 for(int i=0; i<76; i++) { 31 const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 32 double sin_theta, cos_theta; 33 SINCOS(theta, sin_theta, cos_theta); 37 34 38 double sumi = 0.0; 39 40 for(int i=0; i<nordi; i++) { 35 const double termC = sinc(q * c_half * cos_theta); 41 36 42 double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 37 double inner_sum = 0.0; 38 for(int j=0; j<76; j++) { 39 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 40 double sin_phi, cos_phi; 41 SINCOS(phi, sin_phi, cos_phi); 43 42 44 double arg = q * c_half * cos(theta); 45 if (fabs(arg) > 1.e-16) {termC = sin(arg)/arg;} else {termC = 1.0;} 46 47 double sumj = 0.0; 48 49 for(int j=0; j<nordj; j++) { 50 51 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 52 53 // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 54 55 arg = q * a_half * sin(theta) * sin(phi); 56 if (fabs(arg) > 1.e-16) {termA = sin(arg)/arg;} else {termA = 1.0;} 57 58 arg = q * b_half * sin(theta) * cos(phi); 59 if (fabs(arg) > 1.e-16) {termB = sin(arg)/arg;} else {termB = 1.0;} 60 61 double AP = termA * termB * termC; 62 63 sumj += Gauss76Wt[j] * (AP*AP); 64 65 } 66 67 sumj = 0.5 * (v2b-v2a) * sumj; 68 sumi += Gauss76Wt[i] * sumj * sin(theta); 69 43 // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 44 const double termA = sinc(q * a_half * sin_theta * sin_phi); 45 const double termB = sinc(q * b_half * sin_theta * cos_phi); 46 const double AP = termA * termB * termC; 47 inner_sum += Gauss76Wt[j] * AP * AP; 48 } 49 inner_sum = 0.5 * (v2b-v2a) * inner_sum; 50 outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 70 51 } 71 52 72 double answer = 0.5*(v1b-v1a)* sumi;53 double answer = 0.5*(v1b-v1a)*outer_sum; 73 54 74 55 // Normalize by Pi (Eqn. 16). … … 77 58 // The factor 2 appears because the theta integral has been defined between 78 59 // 0 and pi/2, instead of 0 to pi. 79 answer *= (2.0/M_PI); //Form factor P(q)60 answer /= M_PI_2; //Form factor P(q) 80 61 81 62 // Multiply by contrast^2 and volume^2 82 answer *= (sld-solvent_sld)*(sld-solvent_sld)*volume*volume; 63 const double volume = length_a * length_b * length_c; 64 answer *= square((sld-solvent_sld)*volume); 83 65 84 66 // Convert from [1e-12 A-1] to [cm-1] … … 86 68 87 69 return answer; 88 89 70 } 90 91 double Iqxy(double qx, double qy,92 double sld,93 double solvent_sld,94 double length_a,95 double b2a_ratio,96 double c2a_ratio)97 {98 double q = sqrt(qx*qx + qy*qy);99 double intensity = Iq(q, sld, solvent_sld, length_a, b2a_ratio, c2a_ratio);100 return intensity;101 } -
TabularUnified sasmodels/models/rectangular_prism.py ¶
ra807206 rab2aea8 3 3 r""" 4 4 5 This model provides the form factor, *P(q)*, for a rectangular prism.5 This model provides the form factor, $P(q)$, for a rectangular prism. 6 6 7 7 Note that this model is almost totally equivalent to the existing 8 8 :ref:`parallelepiped` model. 9 9 The only difference is that the way the relevant 10 parameters are defined here ( *a*, *b/a*, *c/a* instead of *a*, *b*, *c*)10 parameters are defined here ($a$, $b/a$, $c/a$ instead of $a$, $b$, $c$) 11 11 which allows use of polydispersity with this model while keeping the shape of 12 the prism (e.g. setting *b/a* = 1 and *c/a* = 1and applying polydispersity12 the prism (e.g. setting $b/a = 1$ and $c/a = 1$ and applying polydispersity 13 13 to *a* will generate a distribution of cubes of different sizes). 14 14 Note also that, contrary to :ref:`parallelepiped`, it does not compute … … 24 24 Note also that the angle definitions used in the code and the present 25 25 documentation correspond to those used in (Nayuk, 2012) (see Fig. 1 of 26 that reference), with |theta| corresponding to |alpha|in that paper,26 that reference), with $\theta$ corresponding to $\alpha$ in that paper, 27 27 and not to the usual convention used for example in the 28 28 :ref:`parallelepiped` model. As the present model does not compute … … 30 30 31 31 In this model the scattering from a massive parallelepiped with an 32 orientation with respect to the scattering vector given by |theta|33 and |phi|32 orientation with respect to the scattering vector given by $\theta$ 33 and $\phi$ 34 34 35 35 .. math:: 36 A_P\,(q) = \frac{\sin \bigl( q \frac{C}{2} \cos\theta \bigr)}{\left( q \frac{C}{2}37 \cos\theta \right)} \, \times \, \frac{\sin \bigl( q \frac{A}{2} \sin\theta \sin\phi38 \bigr)}{\left( q \frac{A}{2} \sin\theta \sin\phi \right)} \, \times \, \frac{\sin \bigl(39 q \frac{B}{2} \sin\theta \cos\phi \bigr)}{\left( q \frac{B}{2} \sin\theta \cos\phi \right)}40 36 41 where *A*, *B* and *C* are the sides of the parallelepiped and must fulfill 42 :math:`A \le B \le C`, |theta| is the angle between the *z* axis and the 43 longest axis of the parallelepiped *C*, and |phi| is the angle between the 44 scattering vector (lying in the *xy* plane) and the *y* axis. 37 A_P\,(q) = 38 \frac{\sin \left( \tfrac{1}{2}qC \cos\theta \right) }{\tfrac{1}{2} qC \cos\theta} 39 \,\times\, 40 \frac{\sin \left( \tfrac{1}{2}qA \cos\theta \right) }{\tfrac{1}{2} qA \cos\theta} 41 \,\times\ , 42 \frac{\sin \left( \tfrac{1}{2}qB \cos\theta \right) }{\tfrac{1}{2} qB \cos\theta} 43 44 where $A$, $B$ and $C$ are the sides of the parallelepiped and must fulfill 45 $A \le B \le C$, $\theta$ is the angle between the $z$ axis and the 46 longest axis of the parallelepiped $C$, and $\phi$ is the angle between the 47 scattering vector (lying in the $xy$ plane) and the $y$ axis. 45 48 46 49 The normalized form factor in 1D is obtained averaging over all possible … … 48 51 49 52 .. math:: 50 P(q) = \frac{2}{\pi} \ times \, \int_0^{\frac{\pi}{2}} \,53 P(q) = \frac{2}{\pi} \int_0^{\frac{\pi}{2}} \, 51 54 \int_0^{\frac{\pi}{2}} A_P^2(q) \, \sin\theta \, d\theta \, d\phi 52 55 … … 54 57 55 58 .. math:: 56 I(q) = \ mbox{scale} \times V \times (\rho_{\mbox{p}} -57 \rho_ {\mbox{solvent}})^2 \times P(q)59 I(q) = \text{scale} \times V \times (\rho_\text{p} - 60 \rho_\text{solvent})^2 \times P(q) 58 61 59 where *V* is the volume of the rectangular prism, :math:`\rho_{\mbox{p}}`60 is the scattering length of the parallelepiped, :math:`\rho_{\mbox{solvent}}`62 where $V$ is the volume of the rectangular prism, $\rho_\text{p}$ 63 is the scattering length of the parallelepiped, $\rho_\text{solvent}$ 61 64 is the scattering length of the solvent, and (if the data are in absolute 62 65 units) *scale* represents the volume fraction (which is unitless). … … 125 128 # parameters for demo 126 129 demo = dict(scale=1, background=0, 127 sld=6.3 e-6, sld_solvent=1.0e-6,130 sld=6.3, sld_solvent=1.0, 128 131 length_a=35, b2a_ratio=1, c2a_ratio=1, 129 132 length_a_pd=0.1, length_a_pd_n=10, -
TabularUnified sasmodels/models/sc_paracrystal.c ¶
r0bef47b r4962519 49 49 double da = d_factor*dnn; 50 50 double temp1 = qq*qq*da*da; 51 double temp2 = pow( 1.0-exp(-1.0*temp1) ,3);51 double temp2 = cube(-expm1(-temp1)); 52 52 double temp3 = qq*dnn; 53 53 double temp4 = 2.0*exp(-0.5*temp1); 54 54 double temp5 = exp(-1.0*temp1); 55 55 56 double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5); 57 integrand *= 2.0/M_PI; 56 double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2; 58 57 59 58 return(integrand); 60 59 } 61 60 62 static 63 double sc_crystal_kernel(double q, 61 double Iq(double q, 64 62 double dnn, 65 63 double d_factor, … … 69 67 { 70 68 const double va = 0.0; 71 const double vb = M_PI /2.0; //orientation average, outer integral69 const double vb = M_PI_2; //orientation average, outer integral 72 70 73 71 double summ=0.0; … … 103 101 } 104 102 105 static 106 double sc_crystal_kernel_2d(double q, double q_x, double q_y, 103 double Iqxy(double qx, double qy, 107 104 double dnn, 108 105 double d_factor, … … 114 111 double psi) 115 112 { 116 //convert angle degree to radian 117 theta = theta * M_PI_180; 118 phi = phi * M_PI_180; 119 psi = psi * M_PI_180; 113 double q, cos_a1, cos_a2, cos_a3; 114 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 120 115 121 const double qda_2 = pow(q*d_factor*dnn,2.0); 116 const double qd = q*dnn; 117 const double arg = 0.5*square(qd*d_factor); 118 const double tanh_qd = tanh(arg); 119 const double cosh_qd = cosh(arg); 120 const double Zq = tanh_qd/(1. - cos(qd*cos_a1)/cosh_qd) 121 * tanh_qd/(1. - cos(qd*cos_a2)/cosh_qd) 122 * tanh_qd/(1. - cos(qd*cos_a3)/cosh_qd); 122 123 123 double snt, cnt; 124 SINCOS(theta, snt, cnt); 125 126 double snp, cnp; 127 SINCOS(phi, snp, cnp); 128 129 double sns, cns; 130 SINCOS(psi, sns, cns); 131 132 /// Angles here are respect to detector coordinate instead of against 133 // q coordinate(PRB 36, 3, 1754) 134 // a3 axis orientation 135 136 const double a3_x = cnt * cnp; 137 const double a3_y = snt; 138 139 // Compute the angle btw vector q and the a3 axis 140 double cos_val_a3 = a3_x*q_x + a3_y*q_y; 141 142 // a1 axis orientation 143 const double a1_x = -cnp*sns * snt+snp*cns; 144 const double a1_y = sns*cnt; 145 146 double cos_val_a1 = a1_x*q_x + a1_y*q_y; 147 148 // a2 axis orientation 149 const double a2_x = -snt*cns*cnp-sns*snp; 150 const double a2_y = cnt*cns; 151 152 // a2 axis 153 const double cos_val_a2 = a2_x*q_x + a2_y*q_y; 154 155 // The following test should always pass 156 if (fabs(cos_val_a3)>1.0) { 157 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 158 cos_val_a3 = 1.0; 159 } 160 if (fabs(cos_val_a1)>1.0) { 161 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 162 cos_val_a1 = 1.0; 163 } 164 if (fabs(cos_val_a2)>1.0) { 165 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 166 cos_val_a3 = 1.0; 167 } 168 169 const double a3_dot_q = dnn*q*cos_val_a3; 170 const double a1_dot_q = dnn*q*cos_val_a1; 171 const double a2_dot_q = dnn*q*cos_val_a2; 172 173 // Call Zq=Z1*Z2*Z3 174 double Zq = (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a1_dot_q)+exp(-qda_2)); 175 Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a2_dot_q)+exp(-qda_2)); 176 Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a3_dot_q)+exp(-qda_2)); 177 178 // Use SphereForm directly from libigor 179 double answer = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 180 181 //consider scales 182 const double latticeScale = sphere_volume(radius/dnn); 183 answer *= latticeScale; 184 185 return answer; 124 const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 125 //the occupied volume of the lattice 126 const double lattice_scale = sphere_volume(radius/dnn); 127 return lattice_scale * Fq; 186 128 } 187 188 double Iq(double q,189 double dnn,190 double d_factor,191 double radius,192 double sphere_sld,193 double solvent_sld)194 {195 return sc_crystal_kernel(q,196 dnn,197 d_factor,198 radius,199 sphere_sld,200 solvent_sld);201 }202 203 // Iqxy is never called since no orientation or magnetic parameters.204 double Iqxy(double qx, double qy,205 double dnn,206 double d_factor,207 double radius,208 double sphere_sld,209 double solvent_sld,210 double theta,211 double phi,212 double psi)213 {214 double q = sqrt(qx*qx + qy*qy);215 216 217 return sc_crystal_kernel_2d(q, qx/q, qy/q,218 dnn,219 d_factor,220 radius,221 sphere_sld,222 solvent_sld,223 theta,224 phi,225 psi);226 227 }228 -
TabularUnified sasmodels/models/stacked_disks.c ¶
r0d6e865 r3ac4e1b 14 14 double solvent_sld); 15 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 16 28 static 17 double _kernel(double q q,29 double _kernel(double q, 18 30 double radius, 19 31 double core_sld, … … 22 34 double halfheight, 23 35 double thick_layer, 24 double zi, 36 double sin_alpha, 37 double cos_alpha, 25 38 double sigma_dnn, 26 39 double d, … … 28 41 29 42 { 30 // qq is the q-value for the calculation (1/A)31 32 33 34 43 // q is the q-value for the calculation (1/A) 44 // radius is the core radius of the cylinder (A) 45 // *_sld are the respective SLD's 46 // halfheight is the *Half* CORE-LENGTH of the cylinder = L (A) 47 // zi is the dummy variable for the integration (x in Feigin's notation) 35 48 36 const double besarg1 = qq*radius*sin(zi);37 const double besarg2 = qq*radius*sin(zi);49 const double besarg1 = q*radius*sin_alpha; 50 //const double besarg2 = q*radius*sin_alpha; 38 51 39 const double sinarg1 = qq*halfheight*cos(zi);40 const double sinarg2 = qq*(halfheight+thick_layer)*cos(zi);52 const double sinarg1 = q*halfheight*cos_alpha; 53 const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 41 54 42 const double be1 = sas_J1c(besarg1); 43 const double be2 = sas_J1c(besarg2); 44 const double si1 = sin(sinarg1)/sinarg1; 45 const double si2 = sin(sinarg2)/sinarg2; 55 const double be1 = sas_J1c(besarg1); 56 //const double be2 = sas_J1c(besarg2); 57 const double be2 = be1; 58 const double si1 = sinc(sinarg1); 59 const double si2 = sinc(sinarg2); 46 60 47 const double dr1 = (core_sld-solvent_sld);48 const double dr2 = (layer_sld-solvent_sld);49 50 const double totald=2.0*(thick_layer+halfheight);61 const double dr1 = core_sld - solvent_sld; 62 const double dr2 = layer_sld - solvent_sld; 63 const double area = M_PI*radius*radius; 64 const double totald = 2.0*(thick_layer + halfheight); 51 65 52 const double t1 = area*(2.0*halfheight)*dr1*(si1)*(be1);53 const double t2 = area*dr2*(totald*si2-2.0*halfheight*si1)*(be2);66 const double t1 = area * (2.0*halfheight) * dr1 * si1 * be1; 67 const double t2 = area * dr2 * (totald*si2 - 2.0*halfheight*si1) * be2; 54 68 69 double pq = square(t1 + t2); 55 70 56 double retval =((t1+t2)*(t1+t2))*sin(zi); 71 // loop for the structure factor S(q) 72 double qd_cos_alpha = q*d*cos_alpha; 73 double debye_arg = -0.5*square(qd_cos_alpha*sigma_dnn); 74 double sq=0.0; 75 for (int kk=1; kk<n_stacking; kk++) { 76 sq += (n_stacking-kk) * cos(qd_cos_alpha*kk) * exp(debye_arg*kk); 77 } 78 // end of loop for S(q) 79 sq = 1.0 + 2.0*sq/n_stacking; 57 80 58 // loop for the structure facture S(q) 59 double sqq=0.0; 60 for(int kk=1;kk<n_stacking;kk+=1) { 61 double dexpt=qq*cos(zi)*qq*cos(zi)*d*d*sigma_dnn*sigma_dnn*kk/2.0; 62 sqq=sqq+(n_stacking-kk)*cos(qq*cos(zi)*d*kk)*exp(-1.*dexpt); 63 } 64 65 // end of loop for S(q) 66 sqq=1.0+2.0*sqq/n_stacking; 67 68 retval *= sqq; 69 70 return(retval); 81 return pq * sq; 71 82 } 72 83 … … 83 94 double solvent_sld) 84 95 { 85 /* 96 /* StackedDiscsX : calculates the form factor of a stacked "tactoid" of core shell disks 86 97 like clay platelets that are not exfoliated 87 98 */ 88 double summ = 0.0;//initialize integral99 double summ = 0.0; //initialize integral 89 100 90 double d=2.0*thick_layer+thick_core;91 double halfheight = thick_core/2.0;101 double d = 2.0*thick_layer+thick_core; 102 double halfheight = 0.5*thick_core; 92 103 93 for(int i=0;i<N_POINTS_76;i++) { 94 double zi = (Gauss76Z[i] + 1.0)*M_PI/4.0; 95 double yyy = Gauss76Wt[i] * 96 _kernel(q, 97 radius, 98 core_sld, 99 layer_sld, 100 solvent_sld, 101 halfheight, 102 thick_layer, 103 zi, 104 sigma_dnn, 105 d, 106 n_stacking); 107 summ += yyy; 108 } 104 for(int i=0; i<N_POINTS_76; i++) { 105 double zi = (Gauss76Z[i] + 1.0)*M_PI_4; 106 double sin_alpha, cos_alpha; // slots to hold sincos function output 107 SINCOS(zi, sin_alpha, cos_alpha); 108 double yyy = _kernel(q, 109 radius, 110 core_sld, 111 layer_sld, 112 solvent_sld, 113 halfheight, 114 thick_layer, 115 sin_alpha, 116 cos_alpha, 117 sigma_dnn, 118 d, 119 n_stacking); 120 summ += Gauss76Wt[i] * yyy * sin_alpha; 121 } 109 122 110 double answer = M_PI/4.0*summ;123 double answer = M_PI_4*summ; 111 124 112 //Convert to [cm-1] 113 answer *= 1.0e-4; 114 115 return answer; 116 } 117 118 static double stacked_disks_kernel_2d(double q, double q_x, double q_y, 119 double thick_core, 120 double thick_layer, 121 double radius, 122 double n_stacking, 123 double sigma_dnn, 124 double core_sld, 125 double layer_sld, 126 double solvent_sld, 127 double theta, 128 double phi) 129 { 130 131 double ct, st, cp, sp; 132 133 //convert angle degree to radian 134 theta = theta * M_PI/180.0; 135 phi = phi * M_PI/180.0; 136 137 SINCOS(theta, st, ct); 138 SINCOS(phi, sp, cp); 139 140 // silence compiler warnings about unused variable 141 (void) sp; 142 143 // parallelepiped orientation 144 const double cyl_x = st * cp; 145 const double cyl_y = st * sp; 146 147 // Compute the angle btw vector q and the 148 // axis of the parallelepiped 149 const double cos_val = cyl_x*q_x + cyl_y*q_y; 150 151 // Note: cos(alpha) = 0 and 1 will get an 152 // undefined value from Stackdisc_kern 153 double alpha = acos( cos_val ); 154 155 // Call the IGOR library function to get the kernel 156 double d = 2 * thick_layer + thick_core; 157 double halfheight = thick_core/2.0; 158 double answer = _kernel(q, 159 radius, 160 core_sld, 161 layer_sld, 162 solvent_sld, 163 halfheight, 164 thick_layer, 165 alpha, 166 sigma_dnn, 167 d, 168 n_stacking); 169 170 answer /= sin(alpha); 171 //convert to [cm-1] 172 answer *= 1.0e-4; 173 174 return answer; 125 //Convert to [cm-1] 126 return 1.0e-4*answer; 175 127 } 176 128 … … 179 131 double radius, 180 132 double n_stacking){ 181 double d = 2 * thick_layer + thick_core;182 return acos(-1.0)* radius * radius * d * n_stacking;133 double d = 2.0 * thick_layer + thick_core; 134 return M_PI * radius * radius * d * n_stacking; 183 135 } 184 136 … … 203 155 solvent_sld); 204 156 } 157 158 159 double 160 Iqxy(double qx, double qy, 161 double thick_core, 162 double thick_layer, 163 double radius, 164 double n_stacking, 165 double sigma_dnn, 166 double core_sld, 167 double layer_sld, 168 double solvent_sld, 169 double theta, 170 double phi) 171 { 172 double q, sin_alpha, cos_alpha; 173 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 174 175 double d = 2.0 * thick_layer + thick_core; 176 double halfheight = 0.5*thick_core; 177 double answer = _kernel(q, 178 radius, 179 core_sld, 180 layer_sld, 181 solvent_sld, 182 halfheight, 183 thick_layer, 184 sin_alpha, 185 cos_alpha, 186 sigma_dnn, 187 d, 188 n_stacking); 189 190 //convert to [cm-1] 191 answer *= 1.0e-4; 192 193 return answer; 194 } 195 -
TabularUnified sasmodels/models/star_polymer.c ¶
r2c74c11 r3a48772 6 6 { 7 7 8 double u_2 = radius2 * pow(q,2);8 double u_2 = radius2 * q * q; 9 9 double v = u_2 * arms / (3.0 * arms - 2.0); 10 10 11 double term1 = v - 1.0 + exp(-v);12 double term2 = ((arms - 1.0)/2.0) * pow((1.0 - exp(-v)),2.0);11 double term1 = v + expm1(-v); 12 double term2 = ((arms - 1.0)/2.0) * square(expm1(-v)); 13 13 14 return (2.0 * (term1 + term2)) / (arms * pow(v,2.0));14 return (2.0 * (term1 + term2)) / (arms * v * v); 15 15 16 16 } -
TabularUnified sasmodels/models/surface_fractal.c ¶
ra807206 rb716cc6 1 // Don't need invalid test since fractal_dim_surf is not polydisperse 2 // #define INVALID(v) (v.fractal_dim_surf <= 1.0 || v.fractal_dim_surf >= 3.0) 3 1 4 double form_volume(double radius); 2 5 … … 11 14 double cutoff_length) 12 15 { 13 double pq, sq, mmo, result; 16 // calculate P(q) 17 const double pq = square(sph_j1c(q*radius)); 14 18 15 //Replaced the original formula with Taylor expansion near zero. 16 //pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2); 19 // calculate S(q) 20 // Note: lim q->0 S(q) = -gamma(mmo) cutoff_length^mmo (mmo cutoff_length) 21 // however, the surface fractal formula is invalid outside the range 22 const double mmo = 5.0 - fractal_dim_surf; 23 const double sq = sas_gamma(mmo) * pow(cutoff_length, mmo) 24 * pow(1.0 + square(q*cutoff_length), -0.5*mmo) 25 * sin(-mmo * atan(q*cutoff_length)) / q; 17 26 18 pq = sph_j1c(q*radius); 19 pq = pq*pq; 27 // Empirically determined that the results are valid within this range. 28 // Above 1/r, the form starts to oscillate; below 29 //const double result = (q > 5./(3-fractal_dim_surf)/cutoff_length) && q < 1./radius 30 // ? pq * sq : 0.); 20 31 21 //calculate S(q) 22 mmo = 5.0 - fractal_dim_surf; 23 sq = sas_gamma(mmo)*sin(-(mmo)*atan(q*cutoff_length)); 24 sq *= pow(cutoff_length, mmo); 25 sq /= pow((1.0 + (q*cutoff_length)*(q*cutoff_length)),(mmo/2.0)); 26 sq /= q; 32 double result = pq * sq; 27 33 28 //combine and return 29 result = pq * sq; 30 31 return result; 34 // exclude negative results 35 return result > 0. ? result : 0.; 32 36 } 33 double form_volume(double radius) {34 35 return 1.333333333333333*M_PI*radius*radius*radius;37 double form_volume(double radius) 38 { 39 return M_4PI_3*cube(radius); 36 40 } 37 41 -
TabularUnified sasmodels/models/surface_fractal.py ¶
ra807206 r5c94f41 10 10 .. math:: 11 11 12 I(q) = scale \times P(q)S(q) + background 13 14 .. math:: 15 16 P(q) = F(qR)^2 17 18 .. math:: 19 20 F(x) = \frac{3\left[sin(x)-xcos(x)\right]}{x^3} 21 22 .. math:: 23 24 S(q) = \frac{\Gamma(5-D_S)\zeta^{5-D_S}}{\left[1+(q\zeta)^2 25 \right]^{(5-D_S)/2}} 26 \frac{sin\left[(D_S - 5) tan^{-1}(q\zeta) \right]}{q} 27 28 .. math:: 29 30 scale = scale\_factor \times NV^2(\rho_{particle} - \rho_{solvent})^2 31 32 .. math:: 33 34 V = \frac{4}{3}\pi R^3 12 I(q) &= \text{scale} \times P(q)S(q) + \text{background} \\ 13 P(q) &= F(qR)^2 \\ 14 F(x) &= \frac{3\left[\sin(x)-x\cos(x)\right]}{x^3} \\ 15 S(q) &= \Gamma(5-D_S)\xi^{\,5-D_S}\left[1+(q\xi)^2 \right]^{-(5-D_S)/2} 16 \sin\left[-(5-D_S) \tan^{-1}(q\xi) \right] q^{-1} \\ 17 \text{scale} &= \text{scale_factor}\, N V^2(\rho_\text{particle} - \rho_\text{solvent})^2 \\ 18 V &= \frac{4}{3}\pi R^3 35 19 36 20 where $R$ is the radius of the building block, $D_S$ is the **surface** fractal 37 dimension, | \zeta\| is the cut-off length, $\rho_{solvent}$ is the scattering38 length density of the solvent ,39 and $\rho_{particle}$ is the scatteringlength density of particles.21 dimension, $\xi$ is the cut-off length, $\rho_\text{solvent}$ is the scattering 22 length density of the solvent and $\rho_\text{particle}$ is the scattering 23 length density of particles. 40 24 41 25 .. note:: 42 The surface fractal dimension $D_s$ is only valid if $1<surface\_dim<3$. 43 It is also only valid over a limited $q$ range (see the reference for 44 details) 26 27 The surface fractal dimension is only valid if $1<D_S<3$. The result is 28 only valid over a limited $q$ range, $\tfrac{5}{3-D_S}\xi^{\,-1} < q < R^{-1}$. 29 See the reference for details. 45 30 46 31 … … 89 74 source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "surface_fractal.c"] 90 75 91 demo = dict(scale=1, background= 0,76 demo = dict(scale=1, background=1e-5, 92 77 radius=10, fractal_dim_surf=2.0, cutoff_length=500) 93 78 -
TabularUnified sasmodels/models/teubner_strey.py ¶
rb3f2a24 r8393c74 69 69 B Jakobs, T Sottmann, R Strey, and I Grillo, *J. Chem. Phys.*, 115 (2001), 580 70 70 """ 71 from __future__ import division 71 72 72 73 import numpy as np 73 from numpy import inf, power,pi74 from numpy import inf, pi 74 75 75 76 name = "teubner_strey" … … 89 90 ] 90 91 91 def Iq(q, volfraction , sld, sld_solvent,d,xi):92 def Iq(q, volfraction_a, sld_a, sld_b, d, xi): 92 93 """SAS form""" 93 drho 2 = (sld-sld_solvent)*(sld-sld_solvent)94 drho = sld_a - sld_b 94 95 k = 2.0*pi*xi/d 95 a2 = power(1.0+power(k,2.0),2.0) 96 c1 = -2.0*xi*xi*power(k,2.0)+2*xi*xi 97 c2 = power(xi,4.0) 98 prefactor = 8.0*pi*volfraction*(1.0-volfraction)*drho2*c2/xi 99 #k2 = (2.0*pi/d)*(2.0*pi/d) 100 #xi2 = 1/(xi*xi) 101 #q2 = q*q 102 #result = prefactor/((xi2+k2)*(xi2+k2)+2.0*(xi2-k2)*q2+q2*q2) 96 a2 = (1.0 + k**2)**2 97 c1 = 2.0*xi**2 * (1.0 - k**2) 98 c2 = xi**4 99 prefactor = 8.0*pi * volfraction_a*(1.0 - volfraction_a) * drho**2 * c2/xi 103 100 return 1.0e-4*prefactor / np.polyval([c2, c1, a2], q**2) 104 101 -
TabularUnified sasmodels/models/triaxial_ellipsoid.c ¶
ra807206 r3a48772 10 10 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 11 11 { 12 return 1.333333333333333*M_PI*radius_equat_minor*radius_equat_major*radius_polar;12 return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 13 13 } 14 14 … … 58 58 double psi) 59 59 { 60 double stheta, ctheta; 61 double sphi, cphi; 62 double spsi, cpsi; 60 double q, calpha, cmu, cnu; 61 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu); 63 62 64 const double q = sqrt(qx*qx + qy*qy);65 const double qxhat = qx/q;66 const double qyhat = qy/q;67 SINCOS(theta*M_PI_180, stheta, ctheta);68 SINCOS(phi*M_PI_180, sphi, cphi);69 SINCOS(psi*M_PI_180, spsi, cpsi);70 const double calpha = ctheta*cphi*qxhat + stheta*qyhat;71 const double cnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat;72 const double cmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat;73 63 const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu 74 64 + radius_equat_major*radius_equat_major*cmu*cmu -
TabularUnified sasmodels/models/vesicle.c ¶
r2c74c11 r3a48772 8 8 { 9 9 //note that for the vesicle model, the volume is ONLY the shell volume 10 double volume; 11 volume =4.*M_PI*(radius+thickness)*(radius+thickness)*(radius+thickness)/3; 12 volume -=4.*M_PI*radius*radius*radius/3.; 13 return volume; 10 return M_4PI_3*(cube(radius+thickness) - cube(radius)); 14 11 } 15 12 … … 32 29 // core first, then add in shell 33 30 contrast = sld_solvent-sld; 34 vol = 4.0*M_PI/3.0*radius*radius*radius;35 f = vol *sph_j1c(q*radius)*contrast;31 vol = M_4PI_3*cube(radius); 32 f = vol * sph_j1c(q*radius) * contrast; 36 33 37 34 //now the shell. No volume normalization as this is done by the caller 38 35 contrast = sld-sld_solvent; 39 vol = 4.0*M_PI/3.0*(radius+thickness)*(radius+thickness)*(radius+thickness);40 f += vol *sph_j1c(q*(radius+thickness))*contrast;36 vol = M_4PI_3*cube(radius+thickness); 37 f += vol * sph_j1c(q*(radius+thickness)) * contrast; 41 38 42 39 //rescale to [cm-1]. 43 f2 = volfraction *f*f*1.0e-4;40 f2 = volfraction * f*f*1.0e-4; 44 41 45 return (f2);42 return f2; 46 43 } -
TabularUnified sasmodels/models/vesicle.py ¶
re77872e r3a48772 116 116 ''' 117 117 118 whole = 4. * pi * (radius + thickness) ** 3. / 3.119 core = 4. * pi * radius ** 3. / 3.118 whole = 4./3. * pi * (radius + thickness)**3 119 core = 4./3. * pi * radius**3 120 120 return whole, whole - core 121 121 -
TabularUnified sasmodels/rst2html.py ¶
rb217c71 r0890871 155 155 assert replace_dollar(u"a (again $in parens$) a") == u"a (again :math:`in parens`) a" 156 156 157 def view_rst_app(filename): 158 import wx # type: ignore 159 app = wx.App() 160 view_rst(filename) 161 app.MainLoop() 162 163 157 164 if __name__ == "__main__": 158 test_dollar() 165 import sys 166 view_rst_app(sys.argv[1]) 167 -
TabularUnified sasmodels/sasview_model.py ¶
ra80e64c r8977226 50 50 """ 51 51 import sys 52 import sas 52 import sas # needed in order to set sas.models 53 53 import sas.sascalc.fit 54 54 sys.modules['sas.models'] = sas.sascalc.fit … … 58 58 from sasmodels.conversion_table import CONVERSION_TABLE 59 59 for new_name, conversion in CONVERSION_TABLE.items(): 60 # CoreShellEllipsoidModel => core_shell_ellipsoid:1 61 new_name = new_name.split(':')[0] 60 62 old_name = conversion[0] 61 63 module_attrs = {old_name: find_model(new_name)} … … 228 230 magnetic_params = [] 229 231 fixed = [] 230 for p in model_info.parameters.user_parameters( ):232 for p in model_info.parameters.user_parameters({}, is2d=True): 231 233 if p.type == 'orientation': 232 234 orientation_params.append(p.name) … … 352 354 self.dispersion = collections.OrderedDict() 353 355 self.details = {} 354 for p in self._model_info.parameters.user_parameters( ):356 for p in self._model_info.parameters.user_parameters({}, is2d=True): 355 357 if p.name in hidden: 356 358 continue -
TabularUnified sasmodels/models/be_polyelectrolyte.py ¶
r5df888c rbf9de53 67 67 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 68 68 * **Last Modified by:** Paul Kienzle **Date:** July 24, 2016 69 * **Last Reviewed by:** Piotr rozyczko **Date:** January 27, 2016 69 * **Last Reviewed by:** Paul Butler and Richard Heenan **Date:** 70 October 07, 2016 70 71 """ 71 72 -
TabularUnified sasmodels/models/core_multi_shell.py ¶
rb0c4271 r2d73a53 13 13 14 14 The 2D scattering intensity is the same as $P(q)$ above, regardless of the 15 orientation of the $ q$ vector which is defined as15 orientation of the $\vec q$ vector which is defined as 16 16 17 17 .. math:: … … 29 29 30 30 Our model uses the form factor calculations implemented in a c-library provided 31 by the NIST Center for Neutron Research (Kline, 2006) .31 by the NIST Center for Neutron Research (Kline, 2006) [#kline]_. 32 32 33 33 References … … 35 35 36 36 .. [#] See the :ref:`core-shell-sphere` model documentation. 37 .. [#] L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 38 Plenum Press, New York, 1987. 37 .. [#kline] S R Kline, *J Appl. Cryst.*, 39 (2006) 895 38 .. [#] L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 39 Neutron Scattering*, Plenum Press, New York, 1987. 39 40 40 41 Authorship and Verification … … 43 44 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 44 45 * **Last Modified by:** Paul Kienzle **Date:** September 12, 2016 45 * **Last Reviewed by:** Under Review **Date:** as of October 5, 201646 * **Last Reviewed by:** Paul Kienzle **Date:** September 12, 2016 46 47 """ 47 48 -
TabularUnified sasmodels/models/core_shell_cylinder.py ¶
r40a87fa r755ecc2 1 # core shell cylinder model2 # Note: model title and parameter table are inserted automatically3 1 r""" 4 The form factor is normalized by the particle volume.5 6 2 Definition 7 3 ---------- 8 4 9 5 The output of the 2D scattering intensity function for oriented core-shell 10 cylinders is given by (Kline, 2006) 6 cylinders is given by (Kline, 2006 [#kline]_). The form factor is normalized 7 by the particle volume. 11 8 12 9 .. math:: … … 61 58 The $\theta$ and $\phi$ parameters are not used for the 1D output. 62 59 63 Validation64 ----------65 66 Validation of our code was done by comparing the output of the 1D model to67 the output of the software provided by the NIST (Kline, 2006).68 69 Averaging over a distribution of orientation is done by evaluating the70 equation above. Since we have no other software to compare the71 implementation of the intensity for fully oriented cylinders, we72 compared the result of averaging our 2D output using a uniform73 distribution $p(\theta,\phi) = 1.0$.74 75 60 Reference 76 61 --------- 77 see, for example, Ian Livsey J. Chem. Soc., Faraday Trans. 2, 1987,83, 1445-145278 62 79 2016/03/18 - Description reviewed by RKH 63 .. [#] see, for example, Ian Livsey J. Chem. Soc., Faraday Trans. 2, 1987,83, 64 1445-1452 65 .. [#kline] S R Kline, *J Appl. Cryst.*, 39 (2006) 895 66 67 Authorship and Verification 68 ---------------------------- 69 70 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 71 * **Last Modified by:** Paul Kienzle **Date:** Aug 8, 2016 72 * **Last Reviewed by:** Richard Heenan **Date:** March 18, 2016 80 73 """ 81 74 … … 158 151 theta_pd=15, theta_pd_n=45, 159 152 phi_pd=15, phi_pd_n=1) 160 # ADDED by: RKH ON: 18Mar2016 renamed sld's etc 153
Note: See TracChangeset
for help on using the changeset viewer.