Changeset d2b939c in sasmodels


Ignore:
Timestamp:
Feb 7, 2017 2:34:46 PM (3 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
3a45c2c
Parents:
b8ddf2e (diff), fe8ff99 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-794

Files:
4 added
1 deleted
107 edited
1 moved

Legend:

Unmodified
Added
Removed
  • README.rst

    r84bc3c1 re30d645  
    4646lamellar.py is an example of a single file model with embedded C code. 
    4747 
    48 Magnetism hasn't been implemented yet.  We may want a separate Imagnetic 
    49 calculator with the extra parameters and calculations.   We should 
    50 return all desired spin states together so we can share the work of 
    51 computing the form factors for the different magnetic contrasts.  This 
    52 will mean extending the data handler to support multiple cross sections 
    53 in the same data set. 
    54  
    5548|TravisStatus|_ 
    5649 
  • sasmodels/compare.py

    ra0d75ce rfe25eda  
    3333import datetime 
    3434import traceback 
     35import re 
    3536 
    3637import numpy as np  # type: ignore 
     
    3839from . import core 
    3940from . import kerneldll 
    40 from . import weights 
     41from . import exception 
    4142from .data import plot_theory, empty_data1D, empty_data2D 
    4243from .direct_model import DirectModel 
    4344from .convert import revert_name, revert_pars, constrain_new_to_old 
     45from .generate import FLOAT_RE 
    4446 
    4547try: 
    4648    from typing import Optional, Dict, Any, Callable, Tuple 
    47 except: 
     49except Exception: 
    4850    pass 
    4951else: 
     
    5860sasmodels rewrite. 
    5961 
    60 model is the name of the model to compare (see below). 
     62model or model1,model2 are the names of the models to compare (see below). 
    6163N1 is the number of times to run sasmodels (default=1). 
    6264N2 is the number times to run sasview (default=1). 
     
    7173    -preset*/-random[=seed] preset or random parameters 
    7274    -mono/-poly* force monodisperse/polydisperse 
     75    -magnetic/-nonmagnetic* suppress magnetism 
    7376    -cutoff=1e-5* cutoff value for including a point in polydispersity 
    7477    -pars/-nopars* prints the parameter set or not 
     
    8184    -default/-demo* use demo vs default parameters 
    8285    -html shows the model docs instead of running the model 
    83     -title="graph title" adds a title to the plot 
     86    -title="note" adds note to the plot title, after the model name 
    8487 
    8588Any two calculation engines can be selected for comparison: 
     
    8992    -sasview sets the sasview calculation engine 
    9093 
    91 The default is -single -sasview.  Note that the interpretation of quad 
     94The default is -single -double.  Note that the interpretation of quad 
    9295precision depends on architecture, and may vary from 64-bit to 128-bit, 
    9396with 80-bit floats being common (1e-19 precision). 
    9497 
    9598Key=value pairs allow you to set specific values for the model parameters. 
     99Key=value1,value2 to compare different values of the same parameter. 
     100value can be an expression including other parameters 
    96101""" 
    97102 
     
    108113 
    109114kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True 
     115 
     116# list of math functions for use in evaluating parameters 
     117MATH = dict((k,getattr(math, k)) for k in dir(math) if not k.startswith('_')) 
    110118 
    111119# CRUFT python 2.6 
     
    314322        name = name.split('*')[0] 
    315323 
    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) 
    323334        #q_max = 0.2  # mid q maximum 
    324335        q_max = 1.0  # high q maximum 
     
    326337        pars['rg'] = min(pars['rg'], rg_max) 
    327338 
    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': 
    329346        # Make sure phi sums to 1.0 
    330347        if pars['case_num'] < 2: 
     
    337354            pars['Phi'+c] /= total 
    338355 
     356    elif name == 'stacked_disks': 
     357        pars['n_stacking'] = math.ceil(pars['n_stacking']) 
     358 
    339359def parlist(model_info, pars, is2d): 
    340360    # type: (ModelInfo, ParameterSet, bool) -> str 
     
    344364    lines = [] 
    345365    parameters = model_info.parameters 
     366    magnetic = False 
    346367    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 
    347372        fields = dict( 
    348373            value=pars.get(p.id, p.default), 
     
    352377            pdtype=pars.get(p.id+"_pd_type", 'gaussian'), 
    353378            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.), 
    354382        ) 
    355383        lines.append(_format_par(p.name, **fields)) 
     384        magnetic = magnetic or fields['M0'] != 0. 
    356385    return "\n".join(lines) 
    357386 
     
    359388 
    360389def _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.): 
    362391    # type: (str, float, float, int, float, str) -> str 
    363392    line = "%s: %g"%(name, value) 
     
    367396        line += " +/- %g  (%d points in [-%g,%g] sigma %s)"\ 
    368397                % (pd, n, nsigma, nsigma, pdtype) 
     398    if M0 != 0.: 
     399        line += "  M0:%.3f  mphi:%.1f  mtheta:%.1f" % (M0, mphi, mtheta) 
    369400    return line 
    370401 
     
    380411    for p in pars: 
    381412        if p.endswith("_pd_n"): pars[p] = 0 
     413    return pars 
     414 
     415def 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 
    382426    return pars 
    383427 
     
    466510            if k.endswith('.type'): 
    467511                par = k[:-5] 
     512                if v == 'gaussian': continue 
    468513                cls = dispersers[v if v != 'rectangle' else 'rectangula'] 
    469514                handle = cls() 
    470515                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 
    472523 
    473524        #print("sasview pars",oldpars) 
     
    605656    parameters. 
    606657    """ 
    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 
     662def 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'] 
    609667    data = opts['data'] 
    610668 
     
    612670    base = opts['engines'][0] if n_base else None 
    613671    comp = opts['engines'][1] if n_comp else None 
     672 
    614673    base_time = comp_time = None 
    615674    base_value = comp_value = resid = relerr = None 
     
    620679            base_raw, base_time = time_calculation(base, pars, n_base) 
    621680            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())) 
    624684            _show_invalid(data, base_value) 
    625685        except ImportError: 
     
    630690    if n_comp > 0: 
    631691        try: 
    632             comp_raw, comp_time = time_calculation(comp, pars, n_comp) 
     692            comp_raw, comp_time = time_calculation(comp, pars2, n_comp) 
    633693            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())) 
    636697            _show_invalid(data, comp_value) 
    637698        except ImportError: 
     
    643704        resid = (base_value - comp_value) 
    644705        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 
     719def _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 
     739def 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'] 
    651749 
    652750    # Plot if requested 
    653     if not opts['plot'] and not opts['explore']: return 
    654751    view = opts['view'] 
    655752    import matplotlib.pyplot as plt 
    656753    if limits is None: 
    657754        vmin, vmax = np.Inf, -np.Inf 
    658         if n_base > 0: 
     755        if have_base: 
    659756            vmin = min(vmin, base_value.min()) 
    660757            vmax = max(vmax, base_value.max()) 
    661         if n_comp > 0: 
     758        if have_comp: 
    662759            vmin = min(vmin, comp_value.min()) 
    663760            vmax = max(vmax, comp_value.max()) 
    664761        limits = vmin, vmax 
    665762 
    666     if n_base > 0: 
    667         if n_comp > 0: plt.subplot(131) 
     763    if have_base: 
     764        if have_comp: plt.subplot(131) 
    668765        plot_theory(data, base_value, view=view, use_data=False, limits=limits) 
    669766        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    670767        #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) 
    673772        plot_theory(data, comp_value, view=view, use_data=False, limits=limits) 
    674773        plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 
    675774        #cbar_title = "log I" 
    676     if n_comp > 0 and n_base > 0: 
     775    if have_base and have_comp: 
    677776        plt.subplot(133) 
    678777        if not opts['rel_err']: 
     
    680779        else: 
    681780            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 
    682785        #err,errstr = base/comp,"ratio" 
    683786        plot_theory(data, None, resid=err, view=errview, use_data=False) 
     
    691794    fig = plt.gcf() 
    692795    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']: 
    696799        plt.figure() 
    697800        v = relerr 
     
    708811    return limits 
    709812 
    710 def _print_stats(label, err): 
    711     # type: (str, np.ma.ndarray) -> None 
    712     # work with trimmed data, not the full set 
    713     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)) 
    724813 
    725814 
     
    735824    'preset', 'random', 
    736825    'poly', 'mono', 
     826    'magnetic', 'nonmagnetic', 
    737827    'nopars', 'pars', 
    738828    'rel', 'abs', 
     
    791881    return pars 
    792882 
    793  
     883INTEGER_RE = re.compile("^[+-]?[1-9][0-9]*$") 
     884def 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 @ ? ^ ! , 
     897MODEL_SPLIT = ',' 
    794898def parse_opts(argv): 
    795899    # type: (List[str]) -> Dict[str, Any] 
     
    813917        print("expected parameters: model N1 N2") 
    814918 
    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 None 
    822  
    823919    invalid = [o[1:] for o in flags 
    824920               if o[1:] not in NAME_OPTIONS 
     
    828924        return None 
    829925 
     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 
    830929 
    831930    # pylint: disable=bad-whitespace 
     
    842941        'seed'      : -1,  # default to preset 
    843942        'mono'      : False, 
     943        # Default to magnetic a magnetic moment is set on the command line 
     944        'magnetic'  : False, 
    844945        'show_pars' : False, 
    845946        'show_hist' : False, 
     
    875976        elif arg == '-mono':    opts['mono'] = True 
    876977        elif arg == '-poly':    opts['mono'] = False 
     978        elif arg == '-magnetic':       opts['magnetic'] = True 
     979        elif arg == '-nonmagnetic':    opts['magnetic'] = False 
    877980        elif arg == '-pars':    opts['show_pars'] = True 
    878981        elif arg == '-nopars':  opts['show_pars'] = False 
     
    895998    # pylint: enable=bad-whitespace 
    896999 
     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['np'] = np 
     1066    context.update(pars) 
     1067    context.update((k,v) for k,v in presets.items() if isinstance(v, float)) 
     1068    for k, v in presets.items(): 
     1069        if not isinstance(v, float) and not k.endswith('_type'): 
     1070            presets[k] = eval(v, context) 
     1071    context.update(presets) 
     1072    context.update((k,v) for k,v in presets2.items() if isinstance(v, float)) 
     1073    for k, v in presets2.items(): 
     1074        if not isinstance(v, float) and not k.endswith('_type'): 
     1075            presets2[k] = eval(v, context) 
     1076 
     1077    # update parameters with presets 
     1078    pars.update(presets)  # set value after random to control value 
     1079    pars2.update(presets2)  # set value after random to control value 
     1080    #import pprint; pprint.pprint(model_info) 
     1081 
     1082    same_model = name == name2 and pars == pars 
    8971083    if len(engines) == 0: 
    898         engines.extend(['single', 'double']) 
     1084        if same_model: 
     1085            engines.extend(['single', 'double']) 
     1086        else: 
     1087            engines.extend(['single', 'single']) 
    8991088    elif len(engines) == 1: 
    900         if engines[0][0] == 'double': 
     1089        if not same_model: 
     1090            engines.append(engines[0]) 
     1091        elif engines[0] == 'double': 
    9011092            engines.append('single') 
    9021093        else: 
     
    9051096        del engines[2:] 
    9061097 
    907     n1 = int(positional_args[1]) if len(positional_args) > 1 else 1 
    908     n2 = int(positional_args[2]) if len(positional_args) > 2 else 1 
    9091098    use_sasview = any(engine == 'sasview' and count > 0 
    9101099                      for engine, count in zip(engines, [n1, n2])) 
    911  
    912     # Get demo parameters from model definition, or use default parameters 
    913     # if model does not define demo parameters 
    914     pars = get_pars(model_info, opts['use_demo']) 
    915  
    916  
    917     # Fill in parameters given on the command line 
    918     presets = {} 
    919     for arg in values: 
    920         k, v = arg.split('=', 1) 
    921         if k not in pars: 
    922             # extract base name without polydispersity info 
    923             s = set(p.split('_pd')[0] for p in pars) 
    924             print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 
    925             return None 
    926         presets[k] = float(v) if not k.endswith('type') else v 
    927  
    928     # randomize parameters 
    929     #pars.update(set_pars)  # set value before random to control range 
    930     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 value 
    936     #import pprint; pprint.pprint(model_info) 
    937     constrain_pars(model_info, pars) 
    9381100    if use_sasview: 
    9391101        constrain_new_to_old(model_info, pars) 
     1102        constrain_new_to_old(model_info2, pars2) 
     1103 
    9401104    if opts['show_pars']: 
    941         print(str(parlist(model_info, pars, opts['is2d']))) 
     1105        if not same_model: 
     1106            print("==== %s ====="%model_info.name) 
     1107            print(str(parlist(model_info, pars, opts['is2d']))) 
     1108            print("==== %s ====="%model_info2.name) 
     1109            print(str(parlist(model_info2, pars2, opts['is2d']))) 
     1110        else: 
     1111            print(str(parlist(model_info, pars, opts['is2d']))) 
    9421112 
    9431113    # Create the computational engines 
     
    9481118        base = None 
    9491119    if n2: 
    950         comp = make_engine(model_info, data, engines[1], opts['cutoff']) 
     1120        comp = make_engine(model_info2, data, engines[1], opts['cutoff']) 
    9511121    else: 
    9521122        comp = None 
     
    9551125    # Remember it all 
    9561126    opts.update({ 
    957         'name'      : name, 
    958         'def'       : model_info, 
    959         'n1'        : n1, 
    960         'n2'        : n2, 
    961         'presets'   : presets, 
    962         'pars'      : pars, 
    9631127        'data'      : data, 
     1128        'name'      : [name, name2], 
     1129        'def'       : [model_info, model_info2], 
     1130        'count'     : [n1, n2], 
     1131        'presets'   : [presets, presets2], 
     1132        'pars'      : [pars, pars2], 
    9641133        'engines'   : [base, comp], 
    9651134    }) 
     
    9761145    from .generate import view_html_from_info 
    9771146    app = wx.App() if wx.GetApp() is None else None 
    978     view_html_from_info(opts['def']) 
     1147    view_html_from_info(opts['def'][0]) 
    9791148    if app: app.MainLoop() 
    9801149 
     
    9881157    from bumps.names import FitProblem  # type: ignore 
    9891158    from bumps.gui.app_frame import AppFrame  # type: ignore 
     1159    from bumps.gui import signal 
    9901160 
    9911161    is_mac = "cocoa" in wx.version() 
    9921162    # Create an app if not running embedded 
    9931163    app = wx.App() if wx.GetApp() is None else None 
    994     problem = FitProblem(Explore(opts)) 
     1164    model = Explore(opts) 
     1165    problem = FitProblem(model) 
    9951166    frame = AppFrame(parent=None, title="explore", size=(1000,700)) 
    9961167    if not is_mac: frame.Show() 
     
    9981169    frame.panel.Layout() 
    9991170    frame.panel.aui.Split(0, wx.TOP) 
     1171    def reset_parameters(event): 
     1172        model.revert_values() 
     1173        signal.update_parameters(problem) 
     1174    frame.Bind(wx.EVT_TOOL, reset_parameters, frame.ToolBar.GetToolByPos(1)) 
    10001175    if is_mac: frame.Show() 
    10011176    # If running withing an app, start the main loop 
     
    10151190        config_matplotlib() 
    10161191        self.opts = opts 
    1017         model_info = opts['def'] 
    1018         pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 
     1192        p1, p2 = opts['pars'] 
     1193        m1, m2 = opts['def'] 
     1194        self.fix_p2 = m1 != m2 or p1 != p2 
     1195        model_info = m1 
     1196        pars, pd_types = bumps_model.create_parameters(model_info, **p1) 
    10191197        # Initialize parameter ranges, fixing the 2D parameters for 1D data. 
    10201198        if not opts['is2d']: 
    1021             for p in model_info.parameters.user_parameters(is2d=False): 
     1199            for p in model_info.parameters.user_parameters({}, is2d=False): 
    10221200                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 
    10231201                    k = p.name+ext 
     
    10301208 
    10311209        self.pars = pars 
     1210        self.starting_values = dict((k, v.value) for k, v in pars.items()) 
    10321211        self.pd_types = pd_types 
    10331212        self.limits = None 
     1213 
     1214    def revert_values(self): 
     1215        for k, v in self.starting_values.items(): 
     1216            self.pars[k].value = v 
     1217 
     1218    def model_update(self): 
     1219        pass 
    10341220 
    10351221    def numpoints(self): 
     
    10621248        pars = dict((k, v.value) for k, v in self.pars.items()) 
    10631249        pars.update(self.pd_types) 
    1064         self.opts['pars'] = pars 
    1065         limits = compare(self.opts, limits=self.limits) 
     1250        self.opts['pars'][0] = pars 
     1251        if not self.fix_p2: 
     1252            self.opts['pars'][1] = pars 
     1253        result = run_models(self.opts) 
     1254        limits = plot_models(self.opts, result, limits=self.limits) 
    10661255        if self.limits is None: 
    10671256            vmin, vmax = limits 
    10681257            self.limits = vmax*1e-7, 1.3*vmax 
     1258            import pylab; pylab.clf() 
     1259            plot_models(self.opts, result, limits=self.limits) 
    10691260 
    10701261 
  • sasmodels/conversion_table.py

    r5f1acda r790b16bb  
    55names for each parameter in sasmodels.  This is used by :mod:`convert` to 
    66determine the equivalent parameter set when comparing a sasmodels model to 
    7 the models defined in SasView 3.1. 
     7the models defined in previous versions of SasView and sasmodels. This is now 
     8versioned based on the version number of SasView. 
     9 
     10When any sasmodels parameter or model name is changed, this must be modified to 
     11account for that. 
     12 
     13Usage: 
     14<old_Sasview_version> : { 
     15    <new_model_name> : [ 
     16        <old_model_name> , 
     17        { 
     18            <new_param_name_1> : <old_param_name_1>, 
     19            ... 
     20            <new_param_name_n> : <old_param_name_n> 
     21        } 
     22    ] 
     23} 
     24 
     25Any future parameter and model name changes can and should be given in this 
     26table for future compatibility. 
    827""" 
    928 
    10  
    1129CONVERSION_TABLE = { 
     30    (3,1,2) : { 
    1231    "adsorbed_layer": [ 
    1332        "Core2ndMomentModel", 
     
    143162        } 
    144163    ], 
    145     "core_shell_ellipsoid": [ 
     164    "core_shell_ellipsoid:1": [ 
    146165        "CoreShellEllipsoidModel", 
    147166        { 
     167            "sld_core": "sld_core", 
     168            "sld_shell": "sld_shell", 
     169            "sld_solvent": "sld_solvent", 
     170            "radius_equat_core": "equat_core", 
     171            "x_core": "polar_core", 
     172            "thick_shell": "equat_shell", 
     173            "x_polar_shell": "polar_shell", 
     174            "theta": "axis_theta", 
    148175            "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" 
    157176        } 
    158177    ], 
     
    160179        "CoreShellEllipsoidXTModel", 
    161180        { 
    162             "phi": "axis_phi", 
    163181            "sld_core": "sld_core", 
     182            "sld_shell": "sld_shell", 
     183            "sld_solvent": "sld_solvent", 
     184            "radius_equat_core": "equat_core", 
     185            "thick_shell": "T_shell", 
    164186            "x_core": "X_core", 
    165             "sld_solvent": "sld_solvent", 
    166             "thick_shell": "T_shell", 
    167187            "x_polar_shell": "XpolarShell", 
    168188            "theta": "axis_theta", 
    169             "sld_shell": "sld_shell" 
     189            "phi": "axis_phi", 
    170190        } 
    171191    ], 
     
    173193        "CSParallelepipedModel", 
    174194        { 
     195            "sld_core": "sld_pcore", 
     196            "sld_a": "sld_rimA", 
     197            "sld_b": "sld_rimB", 
     198            "sld_c": "sld_rimC", 
     199            "sld_solvent": "sld_solv", 
     200            "length_a": "shortA", 
     201            "length_b": "midB", 
     202            "length_c": "longC", 
     203            "thick_rim_a": "rimA", 
     204            "thick_rim_c": "rimC", 
     205            "thick_rim_b": "rimB", 
     206            "theta": "parallel_theta", 
    175207            "phi": "parallel_phi", 
    176208            "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" 
    189209        } 
    190210    ], 
     
    210230    ], 
    211231    "correlation_length": [ 
    212         "CorrLengthModel", 
     232        "CorrLength", 
    213233        { 
    214234            "porod_scale": "scale_p", 
     
    217237            "lorentz_exp": "exponent_l", 
    218238            "cor_length": "length_l" 
    219         } 
     239        }, 
     240        "CorrLengthModel" 
    220241    ], 
    221242    "cylinder": [ 
     
    240261        "DABModel", 
    241262        { 
    242             "length": "length" 
     263            "cor_length": "length" 
    243264        } 
    244265    ], 
     
    257278        "EllipticalCylinderModel", 
    258279        { 
     280            "axis_ratio": "r_ratio", 
     281            "radius_minor": "r_minor", 
     282            "sld": "sldCyl", 
     283            "sld_solvent": "sldSolv", 
     284            "theta": "cyl_theta", 
    259285            "phi": "cyl_phi", 
    260286            "psi": "cyl_psi", 
    261             "theta": "cyl_theta", 
    262             "sld": "sldCyl", 
    263             "axis_ratio": "r_ratio", 
    264             "sld_solvent": "sldSolv" 
    265287        } 
    266288    ], 
     
    297319    ], 
    298320    "fractal_core_shell": [ 
    299         "FractalCoreShellModel", 
    300         { 
     321        "FractalCoreShell", 
     322        { 
     323            "sld_core": "core_sld", 
    301324            "sld_shell": "shell_sld", 
    302325            "sld_solvent": "solvent_sld", 
    303             "sld_core": "core_sld" 
    304         } 
     326            "radius": "radius", 
     327            "thickness": "thickness", 
     328            "fractal_dim": "frac_dim", 
     329            "cor_length": "cor_length", 
     330            "volfraction": "volfraction", 
     331        }, 
     332        "FractalCoreShellModel" 
    305333    ], 
    306334    "fuzzy_sphere": [ 
     
    314342    ], 
    315343    "gauss_lorentz_gel": [ 
    316         "GaussLorentzGelModel", 
     344        "GaussLorentzGel", 
    317345        { 
    318346            "gauss_scale": "scale_g", 
     
    321349            "background": "background", 
    322350            "lorentz_scale": "scale_l" 
    323         } 
     351        }, 
     352        "GaussLorentzGelModel" 
    324353    ], 
    325354    "gaussian_peak": [ 
     355        "Peak Gauss Model", 
     356        { 
     357            "peak_pos": "q0", 
     358            "sigma": "B", 
     359        }, 
    326360        "PeakGaussModel", 
    327         { 
    328             "sigma": "B" 
    329         } 
    330361    ], 
    331362    "gel_fit": [ 
     
    334365            "rg": "radius", 
    335366            "lorentz_scale": "lScale", 
     367            "guinier_scale": "gScale", 
    336368            "fractal_dim": "FractalExp", 
    337369            "cor_length": "zeta", 
    338             "guinier_scale": "gScale" 
    339370        } 
    340371    ], 
    341372    "guinier": [ 
     373        "Guinier", 
     374        { 
     375            "rg": "rg" 
     376        }, 
    342377        "GuinierModel", 
    343         { 
    344             "rg": "rg" 
    345         } 
    346378    ], 
    347379    "guinier_porod": [ 
    348         "GuinierPorodModel", 
     380        "GuinierPorod", 
    349381        { 
    350382            "s": "dim", 
    351383            "rg": "rg", 
    352             "m": "m", 
     384            "porod_exp": "m", 
    353385            "scale": "scale", 
    354386            "background": "background" 
    355         } 
     387        }, 
     388        "GuinierPorodModel", 
    356389    ], 
    357390    "hardsphere": [ 
    358391        "HardsphereStructure", 
    359392        { 
    360             "radius_effective_pd": "effect_radius_pd", 
     393            "scale": "scale_factor", 
    361394            "radius_effective": "effect_radius", 
    362             "radius_effective_pd_n": "effect_radius_pd_n" 
    363395        } 
    364396    ], 
     
    366398        "HayterMSAStructure", 
    367399        { 
    368             "salt_concentration": "saltconc", 
    369             "radius_effective_pd": "effect_radius_pd", 
     400            "scale": "scale_factor", 
    370401            "radius_effective": "effect_radius", 
    371             "radius_effective_pd_n": "effect_radius_pd_n" 
     402            "volfraction": "volfraction", 
     403            "charge": "charge", 
     404            "temperature": "temperature", 
     405            "concentration_salt": "saltconc", 
     406            "dielectconst": "dielectconst", 
    372407        } 
    373408    ], 
     
    375410        "HollowCylinderModel", 
    376411        { 
     412            "sld": "sldCyl", 
     413            "sld_solvent": "sldSolv", 
     414            "radius": "core_radius", 
     415            "thickness": "radius", 
     416            "length": "length", 
     417            "theta": "axis_theta", 
    377418            "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" 
    386419        } 
    387420    ], 
     
    389422        "RectangularHollowPrismModel", 
    390423        { 
     424            "sld": "sldPipe", 
     425            "sld_solvent": "sldSolv", 
     426            "length_a": "short_side", 
    391427            "b2a_ratio": "b2a_ratio", 
    392             "length_a": "short_side", 
    393             "sld": "sldPipe", 
    394             "length_c": "c2a_ratio", 
    395             "sld_solvent": "sldSolv", 
    396             "thickness": "thickness" 
     428            "c2a_ratio": "c2a_ratio", 
     429            "thickness": "thickness", 
    397430        } 
    398431    ], 
     
    401434        { 
    402435            "sld": "sldPipe", 
     436            "sld_solvent": "sldSolv", 
     437            "length_a": "short_side", 
    403438            "b2a_ratio": "b2a_ratio", 
    404             "length_a": "short_side", 
    405             "length_c": "c2a_ratio", 
    406             "sld_solvent": "sldSolv" 
     439            "c2a_ratio": "c2a_ratio", 
    407440        } 
    408441    ], 
     
    428461        "LamellarPSHGModel", 
    429462        { 
     463            "sld": "sld_tail", 
     464            "sld_head": "sld_head", 
     465            "sld_solvent": "sld_solvent", 
     466            "length_tail": "deltaT", 
     467            "length_head": "deltaH", 
     468            "d_spacing": "spacing", 
    430469            "Caille_parameter": "caille", 
    431470            "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" 
    437471        } 
    438472    ], 
     
    441475        { 
    442476            "sld": "sld_bi", 
     477            "sld_solvent": "sld_sol", 
     478            "thickness": "delta", 
     479            "d_spacing": "spacing", 
    443480            "Caille_parameter": "caille", 
    444             "Nlayers": "N_plates", 
    445             "sld_solvent": "sld_sol", 
    446             "thickness": "delta" 
     481            "Nlayers": "n_plates", 
    447482        } 
    448483    ], 
     
    451486        { 
    452487            "sld": "sld_layer", 
     488            "sld_solvent": "sld_solvent", 
     489            "thickness": "thickness", 
     490            "d_spacing": "spacing", 
    453491            "sigma_d": "pd_spacing", 
    454             "sld_solvent": "sld_solvent" 
     492            "Nlayers": "Nlayers", 
    455493        } 
    456494    ], 
     
    473511    ], 
    474512    "lorentz": [ 
     513        "Lorentz", 
     514        { 
     515            "cor_length": "length" 
     516        }, 
    475517        "LorentzModel", 
    476         { 
    477             "cor_length": "length" 
    478         } 
    479518    ], 
    480519    "mass_fractal": [ 
     
    497536    ], 
    498537    "mono_gauss_coil": [ 
     538        "Debye", 
     539        { 
     540            "rg": "rg", 
     541            "i_zero": "scale", 
     542            "background": "background", 
     543        }, 
    499544        "DebyeModel", 
    500         { 
    501             "rg": "rg", 
    502             "scale": "scale", 
    503             "background": "background" 
    504         } 
    505545    ], 
    506546    "multilayer_vesicle": [ 
     
    551591    ], 
    552592    "peak_lorentz": [ 
    553         "PeakLorentzModel", 
     593        "Peak Lorentz Model", 
    554594        { 
    555595            "peak_pos": "q0", 
    556596            "peak_hwhm": "B" 
    557         } 
     597        }, 
     598        "PeakLorentzModel", 
    558599    ], 
    559600    "pearl_necklace": [ 
     
    576617            "rg": "rg", 
    577618            "polydispersity": "poly_m", 
     619            "i_zero": "scale", 
     620            "background": "background", 
     621        } 
     622    ], 
     623    "polymer_excl_volume": [ 
     624        "PolymerExclVolume", 
     625        { 
     626            "rg": "rg", 
     627            "scale": "scale", 
     628            "background": "background", 
     629            "porod_exp": "m" 
     630        } 
     631    ], 
     632    "polymer_micelle": [ 
     633        "MicelleSphCoreModel", 
     634        { 
     635            "sld_corona": "rho_corona", 
     636            "sld_solvent": "rho_solv", 
     637            "sld_core": "rho_core", 
     638            "ndensity": "ndensity", 
     639            "v_core": "v_core", 
     640            "v_corona": "v_corona", 
     641            "radius_core": "radius_core", 
     642            "rg": "radius_gyr", 
     643            "d_penetration": "d_penetration", 
     644            "n_aggreg": "n_aggreg", 
     645        } 
     646    ], 
     647    "porod": [ 
     648        "PorodModel", 
     649        { 
    578650            "scale": "scale", 
    579651            "background": "background" 
    580652        } 
    581653    ], 
    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     ], 
    606654    "power_law": [ 
    607655        "PowerLawAbsModel", 
     
    616664        { 
    617665            "scale": "scale", 
    618             "solvent_sld": "sld_solvent", 
     666            "sld_solvent": "sld_solvent", 
    619667            "thickness": "thickness", 
    620668            "beta": "beta", 
     
    643691        { 
    644692            "sld": "sldPipe", 
     693            "length_a": "short_side", 
    645694            "b2a_ratio": "b2a_ratio", 
    646             "length_a": "short_side", 
    647             "length_c": "c2a_ratio", 
     695            "c2a_ratio": "c2a_ratio", 
    648696            "sld_solvent": "sldSolv" 
    649697        } 
     
    716764        "SquareWellStructure", 
    717765        { 
    718             "radius_effective_pd": "effect_radius_pd", 
     766            "scale": "scale_factor", 
    719767            "radius_effective": "effect_radius", 
    720             "radius_effective_pd_n": "effect_radius_pd_n" 
     768            "wellwidth": "wellwidth", 
     769            "welldepth": "welldepth", 
    721770        } 
    722771    ], 
     
    729778            "theta": "axis_theta", 
    730779            "sld_solvent": "solvent_sld", 
    731             "n_stacking": "n_stacking" 
     780            "n_stacking": "n_stacking", 
     781            "thick_layer": "layer_thick", 
     782            "thick_core": "core_thick", 
    732783        } 
    733784    ], 
     
    742793        "StickyHSStructure", 
    743794        { 
    744             "radius_effective_pd": "effect_radius_pd", 
     795            "scale": "scale_factor", 
    745796            "radius_effective": "effect_radius", 
    746             "radius_effective_pd_n": "effect_radius_pd_n" 
    747797        } 
    748798    ], 
     
    758808        "TeubnerStreyModel", 
    759809        { 
    760             "a2": "scale" 
     810            # Note: parameters are completely rewritten in convert.py 
     811            "volfraction_a": "volfraction_a", 
     812            "sld_a": "sld_a", 
     813            "sld_b": "sld_b", 
     814            "d": "d", 
     815            "xi": "xi", 
    761816        } 
    762817    ], 
     
    775830    ], 
    776831    "two_lorentzian": [ 
    777         "TwoLorentzianModel", 
     832        "TwoLorentzian", 
    778833        { 
    779834            "lorentz_scale_1": "scale_1", 
     
    784839            "lorentz_length_1": "length_1", 
    785840            "background": "background" 
    786         } 
     841        }, 
     842        "TwoLorentzianModel", 
    787843    ], 
    788844    "two_power_law": [ 
    789         "TwoPowerLawModel", 
     845        "TwoPowerLaw", 
    790846        { 
    791847            "coefficent_1": "coef_A", 
     
    794850            "background": "background", 
    795851            "crossover": "qc" 
    796         } 
     852        }, 
     853        "TwoPowerLawModel", 
    797854    ], 
    798855    "unified_power_Rg": [ 
     856        "UnifiedPowerRg", 
     857        dict(((field_new+str(index), field_old+str(index)) 
     858              for field_new, field_old in [("rg", "Rg"), 
     859                                           ("power", "power"), 
     860                                           ("G", "G"), 
     861                                           ("B", "B"),] 
     862              for index in range(11)), 
     863             **{ 
     864                   "background": "background", 
     865                   "scale": "scale", 
     866               }), 
    799867        "UnifiedPowerRgModel", 
    800         { 
    801         } 
    802868    ], 
    803869    "vesicle": [ 
     
    808874        } 
    809875    ] 
     876    } 
    810877} 
  • sasmodels/convert.py

    r51241113 r07c8d46  
    22Convert models to and from sasview. 
    33""" 
    4 from __future__ import print_function 
    5  
    6 from os.path import join as joinpath, abspath, dirname 
     4from __future__ import print_function, division 
     5 
     6import re 
    77import math 
    88import warnings 
    99 
    1010from .conversion_table import CONVERSION_TABLE 
     11from .core import load_model_info 
    1112 
    1213# List of models which SasView versions don't contain the explicit 'scale' argument. 
     
    1718    'two_lorentzian', 
    1819    "two_power_law", 
    19     'gel_fit', 
    2020    'gauss_lorentz_gel', 
    2121    'be_polyelectrolyte', 
     
    4949# Convert new style names for polydispersity info to old style names 
    5050PD_DOT = [ 
    51     ("", ""), 
    5251    ("_pd", ".width"), 
    5352    ("_pd_n", ".npts"), 
    5453    ("_pd_nsigma", ".nsigmas"), 
    5554    ("_pd_type", ".type"), 
     55    (".lower", ".lower"), 
     56    (".upper", ".upper"), 
     57    (".fittable", ".fittable"), 
     58    (".std", ".std"), 
     59    (".units", ".units"), 
     60    ("", "") 
    5661    ] 
     62 
     63def _rescale(par, scale): 
     64    return [pk*scale for pk in par] if isinstance(par, list) else par*scale 
     65 
     66def _is_sld(model_info, id): 
     67    """ 
     68    Return True if parameter is a magnetic magnitude or SLD parameter. 
     69    """ 
     70    if id.startswith('M0:'): 
     71        return True 
     72    if '_pd' in id or '.' in id: 
     73        return False 
     74    for p in model_info.parameters.call_parameters: 
     75        if p.id == id: 
     76            return p.type == 'sld' 
     77    # check through kernel parameters in case it is a named as a vector 
     78    for p in model_info.parameters.kernel_parameters: 
     79        if p.id == id: 
     80            return p.type == 'sld' 
     81    return False 
     82 
     83def _rescale_sld(model_info, pars, scale): 
     84    """ 
     85    rescale all sld parameters in the new model definition by *scale* so the 
     86    numbers are nicer.  Relies on the fact that all sld parameters in the 
     87    new model definition end with sld.  For backward conversion use 
     88    *scale=1e-6*.  For forward conversion use *scale=1e6*. 
     89    """ 
     90    return dict((id, (_rescale(v, scale) if _is_sld(model_info, id) else v)) 
     91                for id, v in pars.items()) 
     92 
     93 
     94def _get_translation_table(model_info, version=(3,1,2)): 
     95    conv_param = CONVERSION_TABLE.get(version, {}).get(model_info.id, [None, {}]) 
     96    translation = conv_param[1].copy() 
     97    for p in model_info.parameters.kernel_parameters: 
     98        if p.length > 1: 
     99            newid = p.id 
     100            oldid = translation.get(p.id, p.id) 
     101            translation.pop(newid, None) 
     102            for k in range(1, p.length+1): 
     103                if newid+str(k) not in translation: 
     104                    translation[newid+str(k)] = oldid+str(k) 
     105    # Remove control parameter from the result 
     106    if model_info.control: 
     107        translation[model_info.control] = "CONTROL" 
     108    return translation 
     109 
     110# ========= FORWARD CONVERSION sasview 3.x => sasmodels =========== 
     111def _dot_pd_to_underscore_pd(par): 
     112    if par.endswith(".width"): 
     113        return par[:-6]+"_pd" 
     114    elif par.endswith(".type"): 
     115        return par[:-5]+"_pd_type" 
     116    elif par.endswith(".nsigmas"): 
     117        return par[:-8]+"_pd_nsigma" 
     118    elif par.endswith(".npts"): 
     119        return par[:-5]+"_pd_n" 
     120    else: 
     121        return par 
     122 
     123def _pd_to_underscores(pars): 
     124    return dict((_dot_pd_to_underscore_pd(k), v) for k, v in pars.items()) 
    57125 
    58126def _convert_pars(pars, mapping): 
     
    63131    for new, old in mapping.items(): 
    64132        if old == new: continue 
     133        if old is None: continue 
    65134        for underscore, dot in PD_DOT: 
    66             if old+dot in newpars: 
     135            source = old+dot 
     136            if source in newpars: 
    67137                if new is not None: 
    68                     newpars[new+underscore] = pars[old+dot] 
    69                 del newpars[old+dot] 
     138                    target = new+dot 
     139                else: 
     140                    target = None 
     141                if source != target: 
     142                    if target: 
     143                        newpars[target] = pars[old+dot] 
     144                    del newpars[source] 
    70145    return newpars 
    71146 
    72 def convert_model(name, pars): 
     147def _conversion_target(model_name, version=(3,1,2)): 
     148    """ 
     149    Find the sasmodel name which translates into the sasview name. 
     150 
     151    Note: *CoreShellEllipsoidModel* translates into *core_shell_ellipsoid:1*. 
     152    This is necessary since there is only one variant in sasmodels for the 
     153    two variants in sasview. 
     154    """ 
     155    for sasmodels_name, sasview_dict in \ 
     156            CONVERSION_TABLE.get(version, {}).items(): 
     157        if sasview_dict[0] == model_name: 
     158            return sasmodels_name 
     159    return None 
     160 
     161def _hand_convert(name, oldpars, version=(3,1,2)): 
     162    if version == (3,1,2): 
     163        oldpars = _hand_convert_3_1_2_to_4_1(name, oldpars) 
     164    return oldpars 
     165 
     166def _hand_convert_3_1_2_to_4_1(name, oldpars): 
     167    if name == 'core_shell_parallelepiped': 
     168        # Make sure pd on rim parameters defaults to zero 
     169        # ... probably not necessary. 
     170        oldpars['rimA.width'] = 0.0 
     171        oldpars['rimB.width'] = 0.0 
     172        oldpars['rimC.width'] = 0.0 
     173    elif name == 'core_shell_ellipsoid:1': 
     174        # Reverse translation (from new to old), from core_shell_ellipsoid.c 
     175        #    equat_shell = equat_core + thick_shell 
     176        #    polar_core = equat_core * x_core 
     177        #    polar_shell = equat_core * x_core + thick_shell*x_polar_shell 
     178        # Forward translation (from old to new), inverting reverse translation: 
     179        #    thick_shell = equat_shell - equat_core 
     180        #    x_core = polar_core / equat_core 
     181        #    x_polar_shell = (polar_shell - polar_core)/(equat_shell - equat_core) 
     182        # Auto translation (old <=> new) happens after hand_convert 
     183        #    equat_shell <=> thick_shell 
     184        #    polar_core <=> x_core 
     185        #    polar_shell <=> x_polar_shell 
     186        # So... 
     187        equat_core, equat_shell = oldpars['equat_core'], oldpars['equat_shell'] 
     188        polar_core, polar_shell = oldpars['polar_core'], oldpars['polar_shell'] 
     189        oldpars['equat_shell'] = equat_shell - equat_core 
     190        oldpars['polar_core'] = polar_core / equat_core 
     191        oldpars['polar_shell'] = (polar_shell-polar_core)/(equat_shell-equat_core) 
     192    elif name == 'hollow_cylinder': 
     193        # now uses radius and thickness 
     194        thickness = oldpars['radius'] - oldpars['core_radius'] 
     195        oldpars['radius'] = thickness 
     196        if 'radius.width' in oldpars: 
     197            pd = oldpars['radius.width']*oldpars['radius']/thickness 
     198            oldpars['radius.width'] = pd 
     199    elif name == 'multilayer_vesicle': 
     200        if 'scale' in oldpars: 
     201            oldpars['volfraction'] = oldpars['scale'] 
     202            oldpars['scale'] = 1.0 
     203        if 'scale.lower' in oldpars: 
     204            oldpars['volfraction.lower'] = oldpars['scale.lower'] 
     205        if 'scale.upper' in oldpars: 
     206            oldpars['volfraction.upper'] = oldpars['scale.upper'] 
     207        if 'scale.fittable' in oldpars: 
     208            oldpars['volfraction.fittable'] = oldpars['scale.fittable'] 
     209        if 'scale.std' in oldpars: 
     210            oldpars['volfraction.std'] = oldpars['scale.std'] 
     211        if 'scale.units' in oldpars: 
     212            oldpars['volfraction.units'] = oldpars['scale.units'] 
     213    elif name == 'pearl_necklace': 
     214        pass 
     215        #_remove_pd(oldpars, 'num_pearls', name) 
     216        #_remove_pd(oldpars, 'thick_string', name) 
     217    elif name == 'polymer_micelle': 
     218        if 'ndensity' in oldpars: 
     219            oldpars['ndensity'] /= 1e15 
     220        if 'ndensity.lower' in oldpars: 
     221            oldpars['ndensity.lower'] /= 1e15 
     222        if 'ndensity.upper' in oldpars: 
     223            oldpars['ndensity.upper'] /= 1e15 
     224    elif name == 'rpa': 
     225        # convert scattering lengths from femtometers to centimeters 
     226        for p in "L1", "L2", "L3", "L4": 
     227            if p in oldpars: 
     228                oldpars[p] /= 1e-13 
     229            if p + ".lower" in oldpars: 
     230                oldpars[p + ".lower"] /= 1e-13 
     231            if p + ".upper" in oldpars: 
     232                oldpars[p + ".upper"] /= 1e-13 
     233    elif name == 'spherical_sld': 
     234        j = 0 
     235        while "func_inter" + str(j) in oldpars: 
     236            name = "func_inter" + str(j) 
     237            new_name = "shape" + str(j + 1) 
     238            if oldpars[name] == 'Erf(|nu|*z)': 
     239                oldpars[new_name] = int(0) 
     240            elif oldpars[name] == 'RPower(z^|nu|)': 
     241                oldpars[new_name] = int(1) 
     242            elif oldpars[name] == 'LPower(z^|nu|)': 
     243                oldpars[new_name] = int(2) 
     244            elif oldpars[name] == 'RExp(-|nu|*z)': 
     245                oldpars[new_name] = int(3) 
     246            elif oldpars[name] == 'LExp(-|nu|*z)': 
     247                oldpars[new_name] = int(4) 
     248            else: 
     249                oldpars[new_name] = int(0) 
     250            oldpars.pop(name) 
     251            oldpars['n_shells'] = str(j + 1) 
     252            j += 1 
     253    elif name == 'teubner_strey': 
     254        # basically undoing the entire Teubner-Strey calculations here. 
     255        #    drho = (sld_a - sld_b) 
     256        #    k = 2.0*math.pi*xi/d 
     257        #    a2 = (1.0 + k**2)**2 
     258        #    c1 = 2.0 * xi**2 * (1.0 - k**2) 
     259        #    c2 = xi**4 
     260        #    prefactor = 8.0*math.pi*phi*(1.0-phi)*drho**2*c2/xi 
     261        #    scale = 1e-4*prefactor 
     262        #    oldpars['scale'] = a2/scale 
     263        #    oldpars['c1'] = c1/scale 
     264        #    oldpars['c2'] = c2/scale 
     265 
     266        # need xi, d, sld_a, sld_b, phi=volfraction_a 
     267        # assume contrast is 1.0e-6, scale=1, background=0 
     268        sld_a, sld_b = 1.0, 0. 
     269        drho = sld_a - sld_b 
     270 
     271        # find xi 
     272        p_scale = oldpars['scale'] 
     273        p_c1 = oldpars['c1'] 
     274        p_c2= oldpars['c2'] 
     275        i_1 = 0.5*p_c1/p_c2 
     276        i_2 = math.sqrt(math.fabs(p_scale/p_c2)) 
     277        i_3 = 2/(i_1 + i_2) 
     278        xi = math.sqrt(math.fabs(i_3)) 
     279 
     280        # find d from xi 
     281        k = math.sqrt(math.fabs(1 - 0.5*p_c1/p_c2*xi**2)) 
     282        d = 2*math.pi*xi/k 
     283 
     284        # solve quadratic phi (1-phi) = xi/(1e-4 8 pi drho^2 c2) 
     285        # favour volume fraction in [0, 0.5] 
     286        c = xi / (1e-4 * 8.0 * math.pi * drho**2 * p_c2) 
     287        phi = 0.5 - math.sqrt(0.25 - c) 
     288 
     289        # scale sld_a by 1e-6 because the translator will scale it back 
     290        oldpars.update(volfraction_a=phi, xi=xi, d=d, sld_a=sld_a*1e-6, 
     291                       sld_b=sld_b, scale=1.0) 
     292        oldpars.pop('c1') 
     293        oldpars.pop('c2') 
     294 
     295    return oldpars 
     296 
     297def convert_model(name, pars, use_underscore=False, model_version=(3,1,2)): 
    73298    """ 
    74299    Convert model from old style parameter names to new style. 
    75300    """ 
    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): 
    82     return [pk*scale for pk in par] if isinstance(par, list) else par*scale 
    83  
    84 def _is_sld(modelinfo, id): 
    85     if id.startswith('M0:'): 
    86         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')): 
    89         return False 
    90     for p in modelinfo.parameters.call_parameters: 
    91         if p.id == id: 
    92             return p.type == 'sld' 
    93     # check through kernel parameters in case it is a named as a vector 
    94     for p in modelinfo.parameters.kernel_parameters: 
    95         if p.id == id: 
    96             return p.type == 'sld' 
    97     raise ValueError("unknown parameter %r in conversion"%id) 
    98  
    99 def _unscale_sld(modelinfo, pars): 
    100     """ 
    101     rescale all sld parameters in the new model definition by 1e6 so the 
    102     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)) 
    106                 for id, v in pars.items()) 
    107  
    108 def _remove_pd(pars, key, name): 
    109     """ 
    110     Remove polydispersity from the parameter list. 
    111  
    112     Note: operates in place 
    113     """ 
    114     # Bumps style parameter names 
    115     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 pars 
     301    newpars = pars 
     302    keys = sorted(CONVERSION_TABLE.keys()) 
     303    for i, version in enumerate(keys): 
     304        # Don't allow indices outside list 
     305        next_i = i + 1 
     306        if next_i == len(keys): 
     307            next_i = i 
     308        # If the save state is from a later version, skip the check 
     309        if model_version <= keys[next_i]: 
     310            newname = _conversion_target(name, version) 
     311        else: 
     312            newname = None 
     313        # If no conversion is found, move on 
     314        if newname is None: 
     315            newname = name 
     316            continue 
     317        if ':' in newname:   # core_shell_ellipsoid:1 
     318            model_info = load_model_info(newname[:-2]) 
     319            # Know the table exists and isn't multiplicity so grab it directly 
     320            # Can't use _get_translation_table since that will return the 'bare' 
     321            # version. 
     322            translation = CONVERSION_TABLE.get(version, {})[newname][1] 
     323        else: 
     324            model_info = load_model_info(newname) 
     325            translation = _get_translation_table(model_info, version) 
     326        newpars = _hand_convert(newname, newpars, version) 
     327        newpars = _convert_pars(newpars, translation) 
     328        # TODO: Still not convinced this is the best check 
     329        if not model_info.structure_factor and version == (3,1,2): 
     330            newpars = _rescale_sld(model_info, newpars, 1e6) 
     331        newpars.setdefault('scale', 1.0) 
     332        newpars.setdefault('background', 0.0) 
     333        if use_underscore: 
     334            newpars = _pd_to_underscores(newpars) 
     335        name = newname 
     336    return newname, newpars 
     337 
     338# ========= BACKWARD CONVERSION sasmodels => sasview 3.x =========== 
    122339 
    123340def _revert_pars(pars, mapping): 
     
    146363    return oldname 
    147364 
    148 def _get_translation_table(model_info): 
    149     _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    150     translation = translation.copy() 
    151     for p in model_info.parameters.kernel_parameters: 
    152         if p.length > 1: 
    153             newid = p.id 
    154             oldid = translation.get(p.id, p.id) 
    155             translation.pop(newid, None) 
    156             for k in range(1, p.length+1): 
    157                 if newid+str(k) not in translation: 
    158                     translation[newid+str(k)] = oldid+str(k) 
    159     # Remove control parameter from the result 
    160     if model_info.control: 
    161         translation[model_info.control] = "CONTROL" 
    162     return translation 
     365def _remove_pd(pars, key, name): 
     366    """ 
     367    Remove polydispersity from the parameter list. 
     368 
     369    Note: operates in place 
     370    """ 
     371    # Bumps style parameter names 
     372    width = pars.pop(key+".width", 0.0) 
     373    n_points = pars.pop(key+".npts", 0) 
     374    if width != 0.0 and n_points != 0: 
     375        warnings.warn("parameter %s not polydisperse in sasview %s"%(key, name)) 
     376    pars.pop(key+".nsigmas", None) 
     377    pars.pop(key+".type", None) 
     378    return pars 
    163379 
    164380def _trim_vectors(model_info, pars, oldpars): 
     
    180396        composition_type, parts = model_info.composition 
    181397        if composition_type == 'product': 
    182             translation = {'scale':'scale_factor'} 
    183             translation.update(_get_translation_table(parts[0])) 
     398            translation = _get_translation_table(parts[0]) 
     399            # structure factor models include scale:scale_factor mapping 
    184400            translation.update(_get_translation_table(parts[1])) 
    185401        else: 
     
    187403    else: 
    188404        translation = _get_translation_table(model_info) 
    189     oldpars = _revert_pars(_unscale_sld(model_info, pars), translation) 
     405    oldpars = _revert_pars(_rescale_sld(model_info, pars, 1e-6), translation) 
    190406    oldpars = _trim_vectors(model_info, pars, oldpars) 
    191407 
     
    215431        if name in MODELS_WITHOUT_VOLFRACTION: 
    216432            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'] *= 1e15 
    229         elif name == 'spherical_sld': 
    230             oldpars["CONTROL"] -= 1 
    231             # remove polydispersity from shells 
    232             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 shells 
    236             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) 
    242433        elif name == 'core_multi_shell': 
    243434            # kill extra shells 
     
    252443        elif name == 'core_shell_parallelepiped': 
    253444            _remove_pd(oldpars, 'rimA', name) 
    254         elif name in ['mono_gauss_coil', 'poly_gauss_coil']: 
    255             del oldpars['i_zero'] 
     445            _remove_pd(oldpars, 'rimB', name) 
     446            _remove_pd(oldpars, 'rimC', name) 
     447        elif name == 'hollow_cylinder': 
     448            # now uses radius and thickness 
     449            thickness = oldpars['core_radius'] 
     450            oldpars['radius'] += thickness 
     451            oldpars['radius.width'] *= thickness/oldpars['radius'] 
     452        #elif name in ['mono_gauss_coil', 'poly_gauss_coil']: 
     453        #    del oldpars['i_zero'] 
    256454        elif name == 'onion': 
    257455            oldpars.pop('n_shells', None) 
     456        elif name == 'pearl_necklace': 
     457            _remove_pd(oldpars, 'num_pearls', name) 
     458            _remove_pd(oldpars, 'thick_string', name) 
     459        elif name == 'polymer_micelle': 
     460            if 'ndensity' in oldpars: 
     461                oldpars['ndensity'] *= 1e15 
    258462        elif name == 'rpa': 
    259463            # convert scattering lengths from femtometers to centimeters 
     
    272476                for k in "Kab,Kac,Kad".split(','): 
    273477                    oldpars.pop(k, None) 
     478        elif name == 'spherical_sld': 
     479            oldpars["CONTROL"] -= 1 
     480            # remove polydispersity from shells 
     481            for k in range(1, 11): 
     482                _remove_pd(oldpars, 'thick_flat'+str(k), 'thickness') 
     483                _remove_pd(oldpars, 'thick_inter'+str(k), 'interface') 
     484            # remove extra shells 
     485            for k in range(int(pars['n_shells']), 11): 
     486                oldpars.pop('sld_flat'+str(k), 0) 
     487                oldpars.pop('thick_flat'+str(k), 0) 
     488                oldpars.pop('thick_inter'+str(k), 0) 
     489                oldpars.pop('func_inter'+str(k), 0) 
     490                oldpars.pop('nu_inter'+str(k), 0) 
     491        elif name == 'stacked_disks': 
     492            _remove_pd(oldpars, 'n_stacking', name) 
     493        elif name == 'teubner_strey': 
     494            # basically redoing the entire Teubner-Strey calculations here. 
     495            volfraction = oldpars.pop('volfraction_a') 
     496            xi = oldpars.pop('xi') 
     497            d = oldpars.pop('d') 
     498            sld_a = oldpars.pop('sld_a') 
     499            sld_b = oldpars.pop('sld_b') 
     500            drho = 1e6*(sld_a - sld_b)  # conversion autoscaled these 
     501            k = 2.0*math.pi*xi/d 
     502            a2 = (1.0 + k**2)**2 
     503            c1 = 2.0 * xi**2 * (1.0 - k**2) 
     504            c2 = xi**4 
     505            prefactor = 8.0*math.pi*volfraction*(1.0-volfraction)*drho**2*c2/xi 
     506            scale = 1e-4*prefactor 
     507            oldpars['scale'] = a2/scale 
     508            oldpars['c1'] = c1/scale 
     509            oldpars['c2'] = c2/scale 
    274510 
    275511    #print("convert from",list(sorted(pars))) 
     
    315551        if name in MODELS_WITHOUT_VOLFRACTION: 
    316552            pars['volfraction'] = 1 
    317         if name == 'pearl_necklace': 
    318             pars['string_thickness_pd_n'] = 0 
    319             pars['number_of_pearls_pd_n'] = 0 
     553        if name == 'core_multi_shell': 
     554            pars['n'] = min(math.ceil(pars['n']), 4) 
     555        elif name == 'gel_fit': 
     556            pars['scale'] = 1 
    320557        elif name == 'line': 
    321558            pars['scale'] = 1 
    322559            pars['background'] = 0 
     560        elif name == 'mono_gauss_coil': 
     561            pars['scale'] = 1 
     562        elif name == 'onion': 
     563            pars['n_shells'] = math.ceil(pars['n_shells']) 
     564        elif name == 'pearl_necklace': 
     565            pars['string_thickness_pd_n'] = 0 
     566            pars['number_of_pearls_pd_n'] = 0 
     567        elif name == 'poly_gauss_coil': 
     568            pars['scale'] = 1 
    323569        elif name == 'rpa': 
    324570            pars['case_num'] = int(pars['case_num']) 
    325         elif name == 'mono_gauss_coil': 
    326             pars['i_zero'] = 1 
    327         elif name == 'poly_gauss_coil': 
    328             pars['i_zero'] = 1 
    329         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']) 
    333571        elif name == 'spherical_sld': 
    334572            pars['n_shells'] = math.ceil(pars['n_shells']) 
     
    339577                pars['thickness%d_pd_n'%k] = 0 
    340578                pars['interface%d_pd_n'%k] = 0 
    341  
     579        elif name == 'teubner_strey': 
     580            pars['scale'] = 1 
     581            if pars['volfraction_a'] > 0.5: 
     582                pars['volfraction_a'] = 1.0 - pars['volfraction_a'] 
     583        elif name == 'unified_power_Rg': 
     584            pars['level'] = int(pars['level']) 
     585 
     586def _check_one(name, seed=None): 
     587    """ 
     588    Generate a random set of parameters for *name*, and check that they can 
     589    be converted back to SasView 3.x and forward again to sasmodels.  Raises 
     590    an error if the parameters are changed. 
     591    """ 
     592    from . import compare 
     593 
     594    model_info = load_model_info(name) 
     595 
     596    old_name = revert_name(model_info) 
     597    if old_name is None: 
     598        return 
     599 
     600    pars = compare.get_pars(model_info, use_demo=False) 
     601    pars = compare.randomize_pars(model_info, pars, seed=seed) 
     602    if name == "teubner_strey": 
     603        # T-S model is underconstrained, so fix the assumptions. 
     604        pars['sld_a'], pars['sld_b'] = 1.0, 0.0 
     605    compare.constrain_pars(model_info, pars) 
     606    constrain_new_to_old(model_info, pars) 
     607    old_pars = revert_pars(model_info, pars) 
     608    new_name, new_pars = convert_model(old_name, old_pars, use_underscore=True) 
     609    if 1: 
     610        print("==== %s in ====="%name) 
     611        print(str(compare.parlist(model_info, pars, True))) 
     612        print("==== %s ====="%old_name) 
     613        for k, v in sorted(old_pars.items()): 
     614            print(k, v) 
     615        print("==== %s out ====="%new_name) 
     616        print(str(compare.parlist(model_info, new_pars, True))) 
     617    assert name==new_name, "%r != %r"%(name, new_name) 
     618    for k, v in new_pars.items(): 
     619        assert k in pars, "%s: %r appeared from conversion"%(name, k) 
     620        if isinstance(v, float): 
     621            assert abs(v-pars[k])<=abs(1e-12*v), "%s: %r  %s != %s"%(name, k, v, pars[k]) 
     622        else: 
     623            assert v == pars[k], "%s: %r  %s != %s"%(name, k, v, pars[k]) 
     624    for k, v in pars.items(): 
     625        assert k in pars, "%s: %r not converted"%(name, k) 
     626 
     627def test_backward_forward(): 
     628    from .core import list_models 
     629    for name in list_models('all'): 
     630        L = lambda: _check_one(name, seed=1) 
     631        L.description = name 
     632        yield L 
  • sasmodels/core.py

    r1875f4e r52e9a45  
    5353         "nonmagnetic", "magnetic") 
    5454def list_models(kind=None): 
    55     # type: () -> List[str] 
     55    # type: (str) -> List[str] 
    5656    """ 
    5757    Return the list of available models on the model path. 
  • sasmodels/details.py

    r6dc78e4 rccd5f01  
    230230    npars = kernel.info.parameters.npars 
    231231    nvalues = kernel.info.parameters.nvalues 
    232     scalars = [v[0][0] for v in pairs] 
     232    scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs] 
    233233    values, weights = zip(*pairs[2:npars+2]) if npars else ((),()) 
    234234    length = np.array([len(w) for w in weights]) 
  • sasmodels/generate.py

    r234c532 r1e7b0db0  
    139139    square(x) = x*x 
    140140    cube(x) = x*x*x 
    141     sinc(x) = sin(x)/x, with sin(0)/0 -> 1 
     141    sas_sinx_x(x) = sin(x)/x, with sin(0)/0 -> 1 
    142142    all double precision constants must include the decimal point 
    143143    all double declarations may be converted to half, float, or long double 
     
    166166import re 
    167167import string 
     168from zlib import crc32 
    168169 
    169170import numpy as np  # type: ignore 
     
    356357    return newest 
    357358 
     359def 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)) 
    358366 
    359367def convert_type(source, dtype): 
  • sasmodels/kernel_header.c

    re1d6983 rdaeef4c  
    146146inline double square(double x) { return x*x; } 
    147147inline double cube(double x) { return x*x*x; } 
    148 inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
     148inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    149149 
     150#if 1 
     151//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
     152#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
     153    SINCOS(phi*M_PI_180, sn, cn); \ 
     154    q = sqrt(qx*qx + qy*qy); \ 
     155    cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
     156    sn = sqrt(1 - cn*cn); \ 
     157    } while (0) 
     158#else 
     159// SasView 3.x definition of orientation 
     160#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
     161    SINCOS(theta*M_PI_180, sn, cn); \ 
     162    q = sqrt(qx*qx + qy*qy);\ 
     163    cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \ 
     164    sn = sqrt(1 - cn*cn); \ 
     165    } while (0) 
     166#endif 
     167 
     168#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
     169    q = sqrt(qx*qx + qy*qy); \ 
     170    const double qxhat = qx/q; \ 
     171    const double qyhat = qy/q; \ 
     172    double sin_theta, cos_theta; \ 
     173    double sin_phi, cos_phi; \ 
     174    double sin_psi, cos_psi; \ 
     175    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
     176    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
     177    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
     178    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ 
     179    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ 
     180    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
     181    } while (0) 
  • sasmodels/kernel_template.c

    r0d6e865 r1e7b0db0  
    133133inline double square(double x) { return x*x; } 
    134134inline double cube(double x) { return x*x*x; } 
    135 inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
     135inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    136136 
    137137 
  • sasmodels/kernelcl.py

    rb8ddf2e rd2b939c  
    294294        # so we don't really need to cache things for ourselves.  I'll do so 
    295295        # 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 "")) 
    297298        # Check timestamp on program 
    298299        program, program_timestamp = self.compiled.get(key, (None, np.inf)) 
  • sasmodels/modelinfo.py

    rca55aa0 r85fe7f8  
    548548        return full_list 
    549549 
    550     def user_parameters(self, pars={}, is2d=True): 
     550    def user_parameters(self, pars, is2d=True): 
    551551        # type: (Dict[str, float], bool) -> List[Parameter] 
    552552        """ 
  • sasmodels/models/barbell.c

    r0d6e865 r592343f  
    3333        const double t = Gauss76Z[i]*zm + zb; 
    3434        const double radical = 1.0 - t*t; 
    35         const double bj = sas_J1c(qrst*sqrt(radical)); 
     35        const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
    3636        const double Fq = cos(m*t + b) * radical * bj; 
    3737        total += Gauss76Wt[i] * Fq; 
     
    3939    // translate dx in [-1,1] to dx in [lower,upper] 
    4040    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; 
    4343} 
     44 
     45static 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_2J1x_x(q*radius*sin_alpha); 
     52    const double si = sas_sinx_x(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 
    4458 
    4559double form_volume(double radius_bell, 
     
    4761        double length) 
    4862{ 
    49  
    5063    // bell radius should never be less than radius when this is called 
    5164    const double hdist = sqrt(square(radius_bell) - square(radius)); 
     
    7184        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    7285        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); 
    7987        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8088    } 
     
    93101        double theta, double phi) 
    94102{ 
    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); 
    101105 
    102106    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); 
    112108 
    113109    // Multiply by contrast^2 and convert to cm-1 
  • sasmodels/models/barbell.py

    r0d6e865 rfcb33e4  
    2020.. math:: 
    2121 
    22     I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q)\right> 
     22    I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q,\alpha).sin(\alpha)\right> 
    2323 
    24 where the amplitude $A(q)$ is given as 
     24where the amplitude $A(q,\alpha)$ with the rod axis at angle $\alpha$ to $q$ is given as 
    2525 
    2626.. math:: 
    2727 
    2828    A(q) =&\ \pi r^2L 
    29         \frac{\sin\left(\tfrac12 qL\cos\theta\right)} 
    30              {\tfrac12 qL\cos\theta} 
    31         \frac{2 J_1(qr\sin\theta)}{qr\sin\theta} \\ 
     29        \frac{\sin\left(\tfrac12 qL\cos\alpha\right)} 
     30             {\tfrac12 qL\cos\alpha} 
     31        \frac{2 J_1(qr\sin\alpha)}{qr\sin\alpha} \\ 
    3232        &\ + 4 \pi R^3 \int_{-h/R}^1 dt 
    33         \cos\left[ q\cos\theta 
     33        \cos\left[ q\cos\alpha 
    3434            \left(Rt + h + {\tfrac12} L\right)\right] 
    3535        \times (1-t^2) 
    36         \frac{J_1\left[qR\sin\theta \left(1-t^2\right)^{1/2}\right]} 
    37              {qR\sin\theta \left(1-t^2\right)^{1/2}} 
     36        \frac{J_1\left[qR\sin\alpha \left(1-t^2\right)^{1/2}\right]} 
     37             {qR\sin\alpha \left(1-t^2\right)^{1/2}} 
    3838 
    3939The $\left<\ldots\right>$ brackets denote an average of the structure over 
    40 all orientations. $\left<A^2(q)\right>$ is then the form factor, $P(q)$. 
     40all orientations. $\left<A^2(q,\alpha)\right>$ is then the form factor, $P(q)$. 
    4141The scale factor is equivalent to the volume fraction of cylinders, each of 
    4242volume, $V$. Contrast $\Delta\rho$ is the difference of scattering length 
     
    8585* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    8686* **Last Modified by:** Paul Butler **Date:** March 20, 2016 
    87 * **Last Reviewed by:** Paul Butler **Date:** March 20, 2016 
     87* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
    8888""" 
    8989from numpy import inf 
  • sasmodels/models/bcc_paracrystal.c

    r0bef47b r4962519  
    3131        const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 
    3232        const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 
     33 
    3334        const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 
    34         result = pow(1.0-(temp10*temp10),3)*temp6 
     35        result = cube(1.0 - (temp10*temp10))*temp6 
    3536            / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    3637              * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
     
    8182 
    8283    return answer; 
    83  
    84  
    8584} 
    8685 
    8786 
    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){ 
     87double 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); 
    9194 
    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; 
    9898 
    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); 
    103106 
    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} 
  • sasmodels/models/bcc_paracrystal.py

    rb0c4271 r925ad6e  
    128128# pylint: enable=bad-whitespace, line-too-long 
    129129 
    130 source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal.c"] 
     130source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal.c"] 
    131131 
    132132# parameters for demo 
  • sasmodels/models/be_polyelectrolyte.py

    r5df888c rbf9de53  
    6767* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    6868* **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 
    7071""" 
    7172 
  • sasmodels/models/binary_hard_sphere.c

    re481a39 r925ad6e  
    77    ); 
    88     
    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_sld 
    13     ); 
    14  
    159void calculate_psfs(double qval, 
    1610    double r2, double nf2, 
     
    5549    // /* do form factor calculations  */ 
    5650     
    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; 
    5953     
    6054    n1 = phi1/v1; 
     
    6458    qr2 = r2*q; 
    6559 
    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     sc1 = sph_j1c(qr1); 
    77     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; 
     60    sc1 = sas_3j1x_x(qr1); 
     61    sc2 = sas_3j1x_x(qr2); 
     62    b1 = r1*r1*r1*(rho1-rhos)*M_4PI_3*sc1; 
     63    b2 = r2*r2*r2*(rho2-rhos)*M_4PI_3*sc2; 
    8064    inten = n1*b1*b1*psf11; 
    8165    inten += sqrt(n1*n2)*2.0*b1*b2*psf12; 
     
    8872} 
    8973 
    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 } 
    10374 
    10475void calculate_psfs(double qval, 
  • sasmodels/models/binary_hard_sphere.py

    rb0c4271 r925ad6e  
    108108             ] 
    109109 
    110 source = ["lib/sph_j1c.c", "binary_hard_sphere.c"] 
     110source = ["lib/sas_3j1x_x.c", "binary_hard_sphere.c"] 
    111111 
    112112# parameters for demo and documentation 
  • sasmodels/models/capped_cylinder.c

    r0d6e865 r592343f  
    3939        const double t = Gauss76Z[i]*zm + zb; 
    4040        const double radical = 1.0 - t*t; 
    41         const double bj = sas_J1c(qrst*sqrt(radical)); 
     41        const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
    4242        const double Fq = cos(m*t + b) * radical * bj; 
    4343        total += Gauss76Wt[i] * Fq; 
     
    4545    // translate dx in [-1,1] to dx in [lower,upper] 
    4646    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; 
    4848    return cap_Fq; 
     49} 
     50 
     51static 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_2J1x_x(q*radius*sin_alpha); 
     57    const double si = sas_sinx_x(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; 
    4961} 
    5062 
     
    93105        SINCOS(alpha, sin_alpha, cos_alpha); 
    94106 
    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; 
    101110    } 
    102111    // translate dx in [-1,1] to dx in [lower,upper] 
     
    114123    double theta, double phi) 
    115124{ 
    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); 
    122127 
    123128    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); 
    133130 
    134131    // Multiply by contrast^2 and convert to cm-1 
  • sasmodels/models/capped_cylinder.py

    r0d6e865 rfcb33e4  
    2121.. math:: 
    2222 
    23     I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q)\right> 
     23    I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q,\alpha).sin(\alpha)\right> 
    2424 
    25 where the amplitude $A(q)$ is given as 
     25where the amplitude $A(q,\alpha)$ with the rod axis at angle $\alpha$ to $q$ is given as 
    2626 
    2727.. math:: 
    2828 
    2929    A(q) =&\ \pi r^2L 
    30         \frac{\sin\left(\tfrac12 qL\cos\theta\right)} 
    31             {\tfrac12 qL\cos\theta} 
    32         \frac{2 J_1(qr\sin\theta)}{qr\sin\theta} \\ 
     30        \frac{\sin\left(\tfrac12 qL\cos\alpha\right)} 
     31            {\tfrac12 qL\cos\alpha} 
     32        \frac{2 J_1(qr\sin\alpha)}{qr\sin\alpha} \\ 
    3333        &\ + 4 \pi R^3 \int_{-h/R}^1 dt 
    34         \cos\left[ q\cos\theta 
     34        \cos\left[ q\cos\alpha 
    3535            \left(Rt + h + {\tfrac12} L\right)\right] 
    3636        \times (1-t^2) 
    37         \frac{J_1\left[qR\sin\theta \left(1-t^2\right)^{1/2}\right]} 
    38              {qR\sin\theta \left(1-t^2\right)^{1/2}} 
     37        \frac{J_1\left[qR\sin\alpha \left(1-t^2\right)^{1/2}\right]} 
     38             {qR\sin\alpha \left(1-t^2\right)^{1/2}} 
    3939 
    4040The $\left<\ldots\right>$ brackets denote an average of the structure over 
     
    8888* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    8989* **Last Modified by:** Paul Butler **Date:** September 30, 2016 
    90 * **Last Reviewed by:** Richard Heenan **Date:** March 19, 2016 
     90* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
     91 
    9192""" 
    9293from numpy import inf 
  • sasmodels/models/core_multi_shell.c

    rc5ac2b2 r925ad6e  
    33f_constant(double q, double r, double sld) 
    44{ 
    5   const double bes = sph_j1c(q * r); 
     5  const double bes = sas_3j1x_x(q * r); 
    66  const double vol = M_4PI_3 * cube(r); 
    77  return sld * vol * bes; 
     
    3333  f = 0.; 
    3434  for (int i=0; i<n; i++) { 
    35     f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sph_j1c(q*r); 
     35    f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sas_3j1x_x(q*r); 
    3636    last_sld = sld[i]; 
    3737    r += thickness[i]; 
    3838  } 
    39   f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sph_j1c(q*r); 
     39  f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sas_3j1x_x(q*r); 
    4040  return f * f * 1.0e-4; 
    4141} 
  • sasmodels/models/core_multi_shell.py

    rb0c4271 r925ad6e  
    1313 
    1414The 2D scattering intensity is the same as $P(q)$ above, regardless of the 
    15 orientation of the $q$ vector which is defined as 
     15orientation of the $\vec q$ vector which is defined as 
    1616 
    1717.. math:: 
     
    2929 
    3030Our model uses the form factor calculations implemented in a c-library provided 
    31 by the NIST Center for Neutron Research (Kline, 2006). 
     31by the NIST Center for Neutron Research (Kline, 2006) [#kline]_. 
    3232 
    3333References 
     
    3535 
    3636.. [#] 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. 
    3940 
    4041Authorship and Verification 
     
    4344* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    4445* **Last Modified by:** Paul Kienzle **Date:** September 12, 2016 
    45 * **Last Reviewed by:** Under Review **Date:** as of October 5, 2016 
     46* **Last Reviewed by:** Paul Kienzle **Date:** September 12, 2016 
    4647""" 
    4748 
     
    100101             ] 
    101102 
    102 source = ["lib/sph_j1c.c", "core_multi_shell.c"] 
     103source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 
    103104 
    104105def profile(sld_core, radius, sld_solvent, n, sld, thickness): 
  • sasmodels/models/core_shell_bicelle.c

    r0d6e865 r592343f  
    2626double form_volume(double radius, double thick_rim, double thick_face, double length) 
    2727{ 
    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); 
    2929} 
    3030 
     
    3939              double rhor, 
    4040              double rhosolv, 
    41               double dum) 
     41              double sin_alpha, 
     42              double cos_alpha) 
    4243{ 
    4344    double si1,si2,be1,be2; 
     
    4950    const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick); 
    5051    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; 
    5756 
    58     be1 = sas_J1c(besarg1); 
    59     be2 = sas_J1c(besarg2); 
    60     si1 = sinc(sinarg1); 
    61     si2 = sinc(sinarg2); 
     57    be1 = sas_2J1x_x(besarg1); 
     58    be2 = sas_2J1x_x(besarg2); 
     59    si1 = sas_sinx_x(sinarg1); 
     60    si2 = sas_sinx_x(sinarg2); 
    6261 
    6362    const double t = vol1*dr1*si1*be1 + 
     
    6564                     vol3*dr3*si2*be1; 
    6665 
    67     const double retval = t*t*sn; 
     66    const double retval = t*t*sin_alpha; 
    6867 
    69     return(retval); 
     68    return retval; 
    7069 
    7170} 
     
    8382{ 
    8483    // 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; 
    8786 
    8887    double summ = 0.0; 
    8988    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); 
    9192        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); 
    9395        summ += yyy; 
    9496    } 
     
    9698    // calculate value of integral to return 
    9799    double answer = uplim*summ; 
    98     return(answer); 
     100    return answer; 
    99101} 
    100102 
    101103static double 
    102 bicelle_kernel_2d(double q, double q_x, double q_y, 
     104bicelle_kernel_2d(double qx, double qy, 
    103105          double radius, 
    104106          double thick_rim, 
     
    112114          double phi) 
    113115{ 
    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); 
    117118 
    118     // Cylinder orientation 
    119     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 cylinder 
    123     const double cos_val = cyl_x*q_x + cyl_y*q_y; 
    124     const double alpha = acos( cos_val ); 
    125  
    126     // Get the kernel 
    127119    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); 
    130122 
    131123    answer *= 1.0e-4; 
     
    162154          double phi) 
    163155{ 
    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, 
    167157                      radius, 
    168158                      thick_rim, 
  • sasmodels/models/core_shell_bicelle.py

    ra23639a r3b9a526  
    4141 
    4242    I(Q,\alpha) = \frac{\text{scale}}{V_t} \cdot 
    43         F(Q,\alpha)^2 + \text{background} 
     43        F(Q,\alpha)^2.sin(\alpha) + \text{background} 
    4444 
    4545where 
     
    8585* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    8686* **Last Modified by:** Paul Butler **Date:** September 30, 2016 
    87 * **Last Reviewed by:** Richard Heenan **Date:** October 5, 2016 
     87* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
    8888""" 
    8989 
     
    141141# pylint: enable=bad-whitespace, line-too-long 
    142142 
    143 source = ["lib/Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
     143source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
    144144          "core_shell_bicelle.c"] 
    145145 
     
    156156            phi=0) 
    157157 
    158 qx, qy = 0.4 * cos(90), 0.5 * sin(0) 
     158#qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 
  • sasmodels/models/core_shell_cylinder.c

    r0d6e865 r592343f  
    55    double radius, double thickness, double length, double theta, double phi); 
    66 
    7 // twovd = 2 * volume * delta_rho 
     7// vd = volume * delta_rho 
    88// besarg = q * R * sin(alpha) 
    99// siarg = q * L/2 * cos(alpha) 
    10 double _cyl(double twovd, double besarg, double siarg); 
    11 double _cyl(double twovd, double besarg, double siarg) 
     10double _cyl(double vd, double besarg, double siarg); 
     11double _cyl(double vd, double besarg, double siarg) 
    1212{ 
    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 * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 
    1614} 
    1715 
    1816double form_volume(double radius, double thickness, double length) 
    1917{ 
    20     return M_PI*(radius+thickness)*(radius+thickness)*(length+2*thickness); 
     18    return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
    2119} 
    2220 
     
    3230    const double core_qr = q*radius; 
    3331    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); 
    3633    const double shell_qr = q*(radius + thickness); 
    3734    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); 
    4036    double total = 0.0; 
    4137    // double lower=0, upper=M_PI_2; 
     
    4642        const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    4743        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); 
    5046        total += Gauss76Wt[i] * fq * fq * sn; 
    5147    } 
     
    6662    double phi) 
    6763{ 
    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); 
    7866 
    7967    const double core_qr = q*radius; 
    8068    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); 
    8370    const double shell_qr = q*(radius + thickness); 
    8471    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); 
    8773 
    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); 
    9176    return 1.0e-4 * fq * fq; 
    9277} 
  • sasmodels/models/core_shell_cylinder.py

    r40a87fa rfcb33e4  
    1 # core shell cylinder model 
    2 # Note: model title and parameter table are inserted automatically 
    31r""" 
    4 The form factor is normalized by the particle volume. 
    5  
    62Definition 
    73---------- 
    84 
    95The output of the 2D scattering intensity function for oriented core-shell 
    10 cylinders is given by (Kline, 2006) 
     6cylinders is given by (Kline, 2006 [#kline]_). The form factor is normalized 
     7by the particle volume. 
    118 
    129.. math:: 
    1310 
    14     I(q,\alpha) = \frac{\text{scale}}{V_s} F^2(q) + \text{background} 
     11    I(q,\alpha) = \frac{\text{scale}}{V_s} F^2(q,\alpha).sin(\alpha) + \text{background} 
    1512 
    1613where 
     
    1815.. math:: 
    1916 
    20     F(q) = &\ (\rho_c - \rho_s) V_c 
     17    F(q,\alpha) = &\ (\rho_c - \rho_s) V_c 
    2118           \frac{\sin \left( q \tfrac12 L\cos\alpha \right)} 
    2219                {q \tfrac12 L\cos\alpha} 
     
    6158The $\theta$ and $\phi$ parameters are not used for the 1D output. 
    6259 
    63 Validation 
    64 ---------- 
    65  
    66 Validation of our code was done by comparing the output of the 1D model to 
    67 the output of the software provided by the NIST (Kline, 2006). 
    68  
    69 Averaging over a distribution of orientation is done by evaluating the 
    70 equation above. Since we have no other software to compare the 
    71 implementation of the intensity for fully oriented cylinders, we 
    72 compared the result of averaging our 2D output using a uniform 
    73 distribution $p(\theta,\phi) = 1.0$. 
    74  
    7560Reference 
    7661--------- 
    77 see, for example, Ian Livsey  J. Chem. Soc., Faraday Trans. 2, 1987,83, 1445-1452 
    7862 
    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 
     67Authorship 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 
    8073""" 
    8174 
     
    158151            theta_pd=15, theta_pd_n=45, 
    159152            phi_pd=15, phi_pd_n=1) 
    160 # ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc 
     153 
  • sasmodels/models/core_shell_ellipsoid.c

    r0d6e865 r0a3d9b2  
    3232    const double equat_shell = radius_equat_core + thick_shell; 
    3333    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; 
    3535    return vol; 
    3636} 
     
    4949    const double uplim = 1.0; 
    5050 
    51     double summ = 0.0;   //initialize intergral 
    5251 
    5352    const double delpc = core_sld - shell_sld; //core - shell 
     
    5958    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    6059 
    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; 
    7266    } 
     67    summ *= 0.5*(uplim-lolim); 
    7368 
    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; 
    7971} 
    8072 
    8173static double 
    82 core_shell_ellipsoid_xt_kernel_2d(double q, double q_x, double q_y, 
     74core_shell_ellipsoid_xt_kernel_2d(double qx, double qy, 
    8375          double radius_equat_core, 
    8476          double x_core, 
     
    9183          double phi) 
    9284{ 
    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); 
    10087 
    10188    const double sldcs = core_sld - shell_sld; 
    10289    const double sldss = shell_sld- solvent_sld; 
    103  
    104     // Compute the angle btw vector q and the 
    105     // axis of the cylinder 
    106     const double cos_val = cyl_x*q_x + cyl_y*q_y; 
    10790 
    10891    const double polar_core = radius_equat_core*x_core; 
     
    11295    // Call the IGOR library function to get the kernel: 
    11396    // MUST use gfn4 not gf2 because of the def of params. 
    114     double answer = gfn4(cos_val, 
     97    double answer = gfn4(cos_alpha, 
    11598                  radius_equat_core, 
    11699                  polar_core, 
     
    160143          double phi) 
    161144{ 
    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, 
    165146                       radius_equat_core, 
    166147                       x_core, 
  • sasmodels/models/core_shell_ellipsoid.py

    ref5a314 rdaeef4c  
    125125#             ["name", "units", default, [lower, upper], "type", "description"], 
    126126parameters = [ 
    127     ["radius_equat_core","Ang",     20,   [0, inf],    "volume",      "Equatorial radius of core"], 
     127    ["radius_equat_core","Ang",     20,   [0, inf],   "volume",      "Equatorial radius of core"], 
    128128    ["x_core",        "None",       3,   [0, inf],    "volume",      "axial ratio of core, X = r_polar/r_equatorial"], 
    129129    ["thick_shell",   "Ang",       30,   [0, inf],    "volume",      "thickness of shell at equator"], 
     
    137137# pylint: enable=bad-whitespace, line-too-long 
    138138 
    139 source = ["lib/sph_j1c.c", "lib/gfn.c", "lib/gauss76.c", 
     139source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c", 
    140140          "core_shell_ellipsoid.c"] 
    141141 
     
    162162 
    163163q = 0.1 
    164 phi = pi/6 
    165 qx = q*cos(phi) 
    166 qy = q*sin(phi) 
    167 # After redefinition of angles find new reasonable values for unit test 
    168 #tests = [ 
    169 #    # Accuracy tests based on content in test/utest_coreshellellipsoidXTmodel.py 
    170 #    [{'radius_equat_core': 200.0, 
    171 #      'x_core': 0.1, 
    172 #      'thick_shell': 50.0, 
    173 #      'x_polar_shell': 0.2, 
    174 #      'sld_core': 2.0, 
    175 #      'sld_shell': 1.0, 
    176 #      'sld_solvent': 6.3, 
    177 #      'background': 0.001, 
    178 #      'scale': 1.0, 
    179 #     }, 1.0, 0.00189402], 
     164# tests had in old coords theta=0, phi=0; new coords theta=90, phi=0 
     165qx = q*cos(pi/6.0) 
     166qy = q*sin(pi/6.0) 
     167# 11Jan2017 RKH sorted tests after redefinition of angles 
     168tests = [ 
     169     # Accuracy tests based on content in test/utest_coreshellellipsoidXTmodel.py 
     170    [{'radius_equat_core': 200.0, 
     171      'x_core': 0.1, 
     172      'thick_shell': 50.0, 
     173      'x_polar_shell': 0.2, 
     174      'sld_core': 2.0, 
     175      'sld_shell': 1.0, 
     176      'sld_solvent': 6.3, 
     177      'background': 0.001, 
     178      'scale': 1.0, 
     179     }, 1.0, 0.00189402], 
    180180 
    181181    # Additional tests with larger range of parameters 
    182 #    [{'background': 0.01}, 0.1, 11.6915], 
    183  
    184 #    [{'radius_equat_core': 20.0, 
    185 #      'x_core': 200.0, 
    186 #      'thick_shell': 54.0, 
    187 #      'x_polar_shell': 3.0, 
    188 #      'sld_core': 20.0, 
    189 #      'sld_shell': 10.0, 
    190 #      'sld_solvent': 6.0, 
    191 #      'background': 0.0, 
    192 #      'scale': 1.0, 
    193 #     }, 0.01, 8688.53], 
    194  
    195 #   [{'background': 0.001}, (0.4, 0.5), 0.00690673], 
    196  
    197 #   [{'radius_equat_core': 20.0, 
    198 #      'x_core': 200.0, 
    199 #      'thick_shell': 54.0, 
    200 #      'x_polar_shell': 3.0, 
    201 #      'sld_core': 20.0, 
    202 #      'sld_shell': 10.0, 
    203 #      'sld_solvent': 6.0, 
    204 #      'background': 0.01, 
    205 #      'scale': 0.01, 
    206 #     }, (qx, qy), 0.0100002], 
    207 #    ] 
     182    [{'background': 0.01}, 0.1, 11.6915], 
     183 
     184    [{'radius_equat_core': 20.0, 
     185      'x_core': 200.0, 
     186      'thick_shell': 54.0, 
     187      'x_polar_shell': 3.0, 
     188      'sld_core': 20.0, 
     189      'sld_shell': 10.0, 
     190      'sld_solvent': 6.0, 
     191      'background': 0.0, 
     192      'scale': 1.0, 
     193     }, 0.01, 8688.53], 
     194 
     195   # 2D tests 
     196   [{'background': 0.001, 
     197     'theta': 90.0, 
     198     'phi': 0.0, 
     199     }, (0.4, 0.5), 0.00690673], 
     200 
     201   [{'radius_equat_core': 20.0, 
     202      'x_core': 200.0, 
     203      'thick_shell': 54.0, 
     204      'x_polar_shell': 3.0, 
     205      'sld_core': 20.0, 
     206      'sld_shell': 10.0, 
     207      'sld_solvent': 6.0, 
     208      'background': 0.01, 
     209      'scale': 0.01, 
     210      'theta': 90.0, 
     211      'phi': 0.0, 
     212     }, (qx, qy), 0.01000025], 
     213    ] 
  • sasmodels/models/core_shell_parallelepiped.c

    r2222134 r1e7b0db0  
    3535    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    3636     
    37     double t1, t2, t3, t4, tmp, answer;    
    38     double mu = q * length_b; 
     37    const double mu = 0.5 * q * length_b; 
    3938     
    4039    //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);         
    4241    //double vol = length_a * length_b * length_c +  
    4342    //       2.0 * thick_rim_a * length_b * length_c +  
     
    4645     
    4746    // 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; 
    5049 
    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; 
    5656 
    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) 
    6364 
    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; 
    7270 
    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; 
    7575 
    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 
    9578 
    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 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
     90            const double si2 = sas_sinx_x(mu_proj * cos_uu); 
     91            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); 
     92            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 
     93 
    11794            // 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 
    12097            // To note also that in csparallelepiped.cpp, there is a function called 
    12198            // cspkernel, which seems to make more sense and has the following comment: 
     
    126103            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    127104             
    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; 
    129109 
    130         } 
    131                  
    132         // value of the inner integral 
    133         answer = 0.5 * summj; 
     110        // now sum up the outer integral 
     111        const double si = sas_sinx_x(mu * c_scaled * sigma); 
     112        outer_total += Gauss76Wt[i] * inner_total * si * si; 
     113    } 
     114    outer_total *= 0.5; 
    134115 
    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; 
    154118} 
    155119 
     
    170134    double psi) 
    171135{ 
    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); 
    173138 
    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; 
    177144 
    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; 
    220149    // As for the 1D case, Vot is not used 
    221         //double Vot = (length_a * length_b * length_c + 
     150    //double Vot = (length_a * length_b * length_c + 
    222151    //              2.0 * thick_rim_a * length_b * length_c + 
    223152    //              2.0 * length_a * thick_rim_b * length_c + 
    224153    //              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 
    228155    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    229156    // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, 
    230157    // but for the moment I let it like this until understanding better the code. 
    231         double ta = length_a + 2.0*thick_rim_a; 
     158    double ta = length_a + 2.0*thick_rim_a; 
    232159    double tb = length_a + 2.0*thick_rim_b; 
    233160    double tc = length_a + 2.0*thick_rim_c; 
    234161    //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 = sas_sinx_x(0.5*q*length_a*cos_val_a); 
     163    double siB = sas_sinx_x(0.5*q*length_b*cos_val_b); 
     164    double siC = sas_sinx_x(0.5*q*length_c*cos_val_c); 
     165    double siAt = sas_sinx_x(0.5*q*ta*cos_val_a); 
     166    double siBt = sas_sinx_x(0.5*q*tb*cos_val_b); 
     167    double siCt = sas_sinx_x(0.5*q*tc*cos_val_c); 
    241168     
    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     } 
    272169 
    273170    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    274171    // 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); 
    277176    
    278177    return 1.0e-4 * f * f; 
  • sasmodels/models/core_shell_parallelepiped.py

    r2222134 r797a8e3  
    1 # core_shell_parallelepiped model 
    2 # Note: model title and parameter table are inserted automatically 
    31r""" 
     2Definition 
     3---------- 
     4 
    45Calculates the form factor for a rectangular solid with a core-shell structure. 
    56**The thickness and the scattering length density of the shell or "rim" 
     
    1516of the rectangular solid. 
    1617 
    17 An instrument resolution smeared version of the model is also provided. 
    18  
    19  
    20 Definition 
    21 ---------- 
    2218 
    2319The function calculated is the form factor of the rectangular solid below. 
     
    4137**meaning that there are "gaps" at the corners of the solid.** 
    4238 
    43 The intensity calculated follows the :ref:`parallelepiped` model, with the core-shell 
    44 intensity being calculated as the square of the sum of the amplitudes of the 
    45 core and shell, in the same manner as a core-shell model. 
    46  
    47 **For the calculation of the form factor to be valid, the sides of the solid 
    48 MUST be chosen such that** $A < B < C$. 
    49 **If this inequality is not satisfied, the model will not report an error, 
    50 and the calculation will not be correct.** 
     39The intensity calculated follows the :ref:`parallelepiped` model, with the 
     40core-shell intensity being calculated as the square of the sum of the 
     41amplitudes of the core and shell, in the same manner as a core-shell model. 
     42 
     43.. math:: 
     44 
     45    F_{a}(Q,\alpha,\beta)= 
     46    \Bigg(\frac{sin(Q(L_A+2t_A)/2sin\alpha sin\beta)}{Q(L_A+2t_A)/2sin\alpha 
     47    sin\beta)} 
     48    - \frac{sin(QL_A/2sin\alpha sin\beta)}{QL_A/2sin\alpha sin\beta)} \Bigg) 
     49    + \frac{sin(QL_B/2sin\alpha sin\beta)}{QL_B/2sin\alpha sin\beta)} 
     50    + \frac{sin(QL_C/2sin\alpha sin\beta)}{QL_C/2sin\alpha sin\beta)} 
     51 
     52.. note:: 
     53 
     54    For the calculation of the form factor to be valid, the sides of the solid 
     55    MUST be chosen such that** $A < B < C$. 
     56    If this inequality is not satisfied, the model will not report an error, 
     57    but the calculation will not be correct and thus the result wrong. 
    5158 
    5259FITTING NOTES 
    5360If the scale is set equal to the particle volume fraction, |phi|, the returned 
    5461value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    55 However, **no interparticle interference effects are included in this calculation.** 
     62However, **no interparticle interference effects are included in this 
     63calculation.** 
    5664 
    5765There are many parameters in this model. Hold as many fixed as possible with 
     
    6876and length $(C+2t_C)$ values, and used as the effective radius 
    6977for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    70  
    71 .. Comment by Miguel Gonzalez: 
    72    The later seems to contradict the previous statement that interparticle interference 
    73    effects are not included. 
    7478 
    7579To provide easy access to the orientation of the parallelepiped, we define the 
     
    98102---------- 
    99103 
    100 P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 
    101 Equations (1), (13-14). (in German) 
    102  
     104.. [#] P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 
     105    Equations (1), (13-14). (in German) 
     106.. [#] D Singh (2009). *Small angle scattering studies of self assembly in 
     107   lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available 
     108   from Proquest <http://search.proquest.com/docview/304915826?accountid 
     109   =26379>`_ 
     110 
     111Authorship and Verification 
     112---------------------------- 
     113 
     114* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     115* **Last Modified by:** Paul Butler **Date:** September 30, 2016 
     116* **Last Reviewed by:** Miguel Gonzales **Date:** March 21, 2016 
    103117""" 
    104118 
     
    164178# parameters for demo 
    165179demo = dict(scale=1, background=0.0, 
    166             sld_core=1e-6, sld_a=2e-6, sld_b=4e-6, 
    167             sld_c=2e-6, sld_solvent=6e-6, 
     180            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6, 
    168181            length_a=35, length_b=75, length_c=400, 
    169182            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10, 
     
    177190            theta_pd=10, theta_pd_n=1, 
    178191            phi_pd=10, phi_pd_n=1, 
    179             psi_pd=10, psi_pd_n=10) 
     192            psi_pd=10, psi_pd_n=1) 
    180193 
    181194qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
  • sasmodels/models/core_shell_sphere.c

    r2c74c11 r3a48772  
    1616double form_volume(double radius, double thickness) 
    1717{ 
    18     return 4.0 * M_PI / 3.0 * pow((radius + thickness), 3); 
     18    return M_4PI_3 * cube(radius + thickness); 
    1919} 
  • sasmodels/models/core_shell_sphere.py

    r9a4811a r925ad6e  
    7575# pylint: enable=bad-whitespace, line-too-long 
    7676 
    77 source = ["lib/sph_j1c.c", "lib/core_shell.c", "core_shell_sphere.c"] 
     77source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] 
    7878 
    7979demo = dict(scale=1, background=0, radius=60, thickness=10, 
     
    9595    """ 
    9696    return (1, 1) 
    97     whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3) 
    98     core = 4.0 * pi / 3.0 * radius * radius * radius 
     97    whole = 4.0/3.0 * pi * (radius + thickness)**3 
     98    core = 4.0/3.0 * pi * radius**3 
    9999    return whole, whole - core 
    100100 
  • sasmodels/models/cylinder.c

    r9cc7fca r592343f  
    1818    const double qr = q*radius; 
    1919    const double qh = q*0.5*length;  
    20     return  sas_J1c(qr*sn) * sinc(qh*cn) ; 
     20    return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn); 
    2121} 
    2222 
     
    5858    double phi) 
    5959{ 
    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); 
    7063    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); 
    7266} 
  • sasmodels/models/cylinder.py

    r4cdd0cc rb7e8b94  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    4 The form factor is normalized by the particle volume V = \piR^2L. 
     4 
    55For information about polarised and magnetic scattering, see 
    66the :ref:`magnetism` documentation. 
     
    1414.. math:: 
    1515 
    16     P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha) + \text{background} 
     16    P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha).sin(\alpha) + \text{background} 
    1717 
    1818where 
     
    2525           \frac{J_1 \left(q R \sin \alpha\right)}{q R \sin \alpha} 
    2626 
    27 and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V$ 
     27and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V =\pi R^2L$ 
    2828is the volume of the cylinder, $L$ is the length of the cylinder, $R$ is the 
    2929radius of the cylinder, and $\Delta\rho$ (contrast) is the scattering length 
     
    3535.. math:: 
    3636 
    37     F^2(q)=\int_{0}^{\pi/2}{F^2(q,\theta)\sin(\theta)d\theta} 
     37    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}=\int_{0}^{1}{F^2(q,u)du} 
    3838 
    3939 
    40 To provide easy access to the orientation of the cylinder, we define the 
    41 axis of the cylinder using two angles $\theta$ and $\phi$. Those angles 
     40Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with  
     41$sin(\alpha)=\sqrt{1-u^2}$.  
     42 
     43The output of the 1D scattering intensity function for randomly oriented 
     44cylinders is thus given by 
     45 
     46.. math:: 
     47 
     48    P(q) = \frac{\text{scale}}{V} 
     49        \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background} 
     50 
     51 
     52NB: The 2nd virial coefficient of the cylinder is calculated based on the 
     53radius and length values, and used as the effective radius for $S(q)$ 
     54when $P(q) \cdot S(q)$ is applied. 
     55 
     56For oriented cylinders, we define the direction of the 
     57axis of the cylinder using two angles $\theta$ (note this is not the 
     58same as the scattering angle used in q) and $\phi$. Those angles 
    4259are defined in :numref:`cylinder-angle-definition` . 
    4360 
     
    4865    Definition of the angles for oriented cylinders. 
    4966 
    50  
    51 NB: The 2nd virial coefficient of the cylinder is calculated based on the 
    52 radius and length values, and used as the effective radius for $S(q)$ 
    53 when $P(q) \cdot S(q)$ is applied. 
    54  
    55 The output of the 1D scattering intensity function for randomly oriented 
    56 cylinders is then given by 
    57  
    58 .. math:: 
    59  
    60     P(q) = \frac{\text{scale}}{V} 
    61         \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background} 
    62  
    63 The $\theta$ and $\phi$ parameters are not used for the 1D output. 
     67The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data. 
    6468 
    6569Validation 
     
    7478 
    7579    P(q) = \int_0^{\pi/2} d\phi 
    76         \int_0^\pi p(\alpha) P_0(q,\alpha) \sin \alpha\ d\alpha 
     80        \int_0^\pi p(\theta) P_0(q,\theta) \sin \theta\ d\theta 
    7781 
    7882 
    7983where $p(\theta,\phi) = 1$ is the probability distribution for the orientation 
    80 and $P_0(q,\alpha)$ is the scattering intensity for the fully oriented 
     84and $P_0(q,\theta)$ is the scattering intensity for the fully oriented 
    8185system, and then comparing to the 1D result. 
    8286 
     
    145149 
    146150qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    147 # After redefinition of angles, find new tests values  
    148 #tests = [[{}, 0.2, 0.042761386790780453], 
    149 #         [{}, [0.2], [0.042761386790780453]], 
    150 #         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 
    151 #         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 
    152 #        ] 
     151# After redefinition of angles, find new tests values.  Was 10 10 in old coords 
     152tests = [[{}, 0.2, 0.042761386790780453], 
     153        [{}, [0.2], [0.042761386790780453]], 
     154#  new coords     
     155        [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], 
     156        [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], 
     157# old coords   [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 
     158#              [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 
     159        ] 
    153160del qx, qy  # not necessary to delete, but cleaner 
    154161# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc 
  • sasmodels/models/dab.py

    ra807206 r4962519  
    5959 
    6060Iq = """ 
    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)); 
    6363     
    6464    return numerator / denominator ; 
  • sasmodels/models/ellipsoid.c

    r0d6e865 r925ad6e  
    88{ 
    99    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 = sas_3j1x_x(q*r); 
    1321 
    1422    return f*f; 
     
    4957    double phi) 
    5058{ 
    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); 
    5962    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    6063 
  • sasmodels/models/ellipsoid.py

    r0d6e865 r925ad6e  
    135135             ] 
    136136 
    137 source = ["lib/sph_j1c.c", "lib/gauss76.c", "ellipsoid.c"] 
     137source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
    138138 
    139139def ER(radius_polar, radius_equatorial): 
  • sasmodels/models/elliptical_cylinder.c

    ra807206 r592343f  
    66 
    77 
    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) 
     8double 
     9form_volume(double radius_minor, double r_ratio, double length) 
    2810{ 
    2911    return M_PI * radius_minor * radius_minor * r_ratio * length; 
    3012} 
    3113 
    32 double Iq(double q, double radius_minor, double r_ratio, double length, 
    33           double sld, double solvent_sld) { 
     14double 
     15Iq(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; 
    3424 
    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)); 
    4728 
    4829    //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++) { 
    5432        //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++) { 
    5938            //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_2J1x_x(q*r); 
     42            inner_sum += Gauss20Wt[j] * be * be; 
    6343        } 
    6444        //now calculate the value of the inner integral 
    65         answer = (vbj-vaj)/2.0*summj; 
     45        inner_sum *= 0.5*(vbj-vaj); 
    6646 
    6747        //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 = sas_sinx_x(q*0.5*length*cos_val); 
     49        outer_sum += Gauss76Wt[i] * inner_sum * si * si; 
    7250    } 
     51    outer_sum *= 0.5*(vb-va); 
    7352 
    7453    //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; 
    7855 
     56    // scale by contrast and volume, and convert to to 1/cm units 
    7957    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; 
    8160} 
    8261 
    8362 
    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; 
     63double 
     64Iqxy(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); 
    9271 
    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_2J1x_x(q*r); 
     76    const double si = sas_sinx_x(q*0.5*length*cos_val); 
     77    const double Aq = be * si; 
     78    const double delrho = sld - solvent_sld; 
    16179    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); 
    16381} 
  • sasmodels/models/elliptical_cylinder.py

    ra807206 rfcb33e4  
    2020 
    2121    I(\vec q)=\frac{1}{V_\text{cyl}}\int{d\psi}\int{d\phi}\int{ 
    22         p(\theta,\phi,\psi)F^2(\vec q,\alpha,\psi)\sin(\theta)d\theta} 
     22        p(\theta,\phi,\psi)F^2(\vec q,\alpha,\psi)\sin(\alpha)d\alpha} 
    2323 
    2424with the functions 
     
    2626.. math:: 
    2727 
    28     F(\vec q,\alpha,\psi) = 2\frac{J_1(a)\sin(b)}{ab} 
     28    F(q,\alpha,\psi) = 2\frac{J_1(a)\sin(b)}{ab} 
    2929 
    3030where 
     
    3232.. math:: 
    3333 
    34     a &= \vec q\sin(\alpha)\left[ 
    35         r^2_\text{major}\sin^2(\psi)+r^2_\text{minor}\cos(\psi) \right]^{1/2} 
     34    a = qr'\sin(\alpha) 
     35     
     36    b = q\frac{L}{2}\cos(\alpha) 
     37     
     38    r'=\frac{r_{minor}}{\sqrt{2}}\sqrt{(1+\nu^{2}) + (1-\nu^{2})cos(\psi)} 
    3639 
    37     b &= \vec q\frac{L}{2}\cos(\alpha) 
    3840 
    39 and the angle $\Psi$ is defined as the orientation of the major axis of the 
     41and the angle $\psi$ is defined as the orientation of the major axis of the 
    4042ellipse with respect to the vector $\vec q$. The angle $\alpha$ is the angle 
    4143between the axis of the cylinder and $\vec q$. 
     
    9597 
    9698L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 
    97 Neutron Scattering*, Plenum, New York, (1987) 
     99Neutron Scattering*, Plenum, New York, (1987) [see table 3.4] 
     100 
     101Authorship and Verification 
     102---------------------------- 
     103 
     104* **Author:** 
     105* **Last Modified by:**  
     106* **Last Reviewed by:**  Richard Heenan - corrected equation in docs **Date:** December 21, 2016 
     107 
    98108""" 
    99109 
  • sasmodels/models/fcc_paracrystal.c

    r0bef47b r4962519  
    3232 
    3333        const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
    34         result = pow((1.0-(temp10*temp10)),3)*temp6 
     34        result = cube(1.0-(temp10*temp10))*temp6 
    3535            / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    3636              * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
     
    8585} 
    8686 
     87double 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); 
    8794 
    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); 
    91105 
    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); 
    98107 
    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} 
  • sasmodels/models/fcc_paracrystal.py

    r0bef47b r925ad6e  
    116116# pylint: enable=bad-whitespace, line-too-long 
    117117 
    118 source = ["lib/sph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal.c"] 
     118source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal.c"] 
    119119 
    120120# parameters for demo 
  • sasmodels/models/flexible_cylinder.c

    r4937980 r592343f  
    1414{ 
    1515    const double contrast = sld - solvent_sld; 
    16     const double cross_section = sas_J1c(q*radius); 
     16    const double cross_section = sas_2J1x_x(q*radius); 
    1717    const double volume = M_PI*radius*radius*length; 
    1818    const double flex = Sk_WR(q, length, kuhn_length); 
  • sasmodels/models/flexible_cylinder_elliptical.c

    rb66d38e r592343f  
    22double Iq(double q, double length, double kuhn_length, double radius, 
    33          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); 
    64double flexible_cylinder_ex_kernel(double q, double length, double kuhn_length, 
    75                                double radius, double axis_ratio, double sld, 
     
    1715elliptical_crosssection(double q, double a, double b) 
    1816{ 
    19     double summ=0.0; 
     17    double sum=0.0; 
    2018 
    2119    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; 
    2321        double sn, cn; 
    2422        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_2J1x_x(arg); 
     25        sum += Gauss76Wt[i] * yyy * yyy; 
    2926    } 
    30  
    31     summ /= 2.0; 
    32     return(summ); 
     27    sum *= 0.5; 
     28    return(sum); 
    3329 
    3430} 
     
    7773} 
    7874 
    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 } 
  • sasmodels/models/fractal.c

    r2c74c11 r925ad6e  
    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) 
    82 
    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) 
     3static double 
     4Iq(double q, 
     5   double volfraction, 
     6   double radius, 
     7   double fractal_dim, 
     8   double cor_length, 
     9   double sld_block, 
     10   double sld_solvent) 
    1611{ 
    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 
    3514    //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)*sas_3j1x_x(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; 
    5324} 
    5425 
  • sasmodels/models/fractal.py

    rd1cfa86 r925ad6e  
    9797# pylint: enable=bad-whitespace, line-too-long 
    9898 
    99 source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "fractal.c"] 
     99source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"] 
    100100 
    101101demo = dict(volfraction=0.05, 
  • sasmodels/models/fractal_core_shell.c

    ra807206 rbdd08df  
    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) 
     1static double 
     2form_volume(double radius, double thickness) 
    143{ 
    15     return 4.0 * M_PI / 3.0 * pow((radius + thickness), 3); 
     4    return M_4PI_3 * cube(radius + thickness); 
    165} 
    176 
    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) { 
     7static double 
     8Iq(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    //The radius for the building block of the core shell particle that is  
     19    //needed by the Teixeira fractal S(q) is the radius of the whole particle. 
     20    const double cs_radius = radius + thickness; 
     21    const double sq = fractal_sq(q, cs_radius, fractal_dim, cor_length); 
     22    const double pq = core_shell_kernel(q, radius, thickness, 
     23                                        core_sld, shell_sld, solvent_sld); 
    2724 
    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; 
     25    // Note: core_shell_kernel already performs the 1e-4 unit conversion 
     26    return volfraction * sq * pq; 
    4627} 
    4728 
  • sasmodels/models/fractal_core_shell.py

    ra807206 r925ad6e  
    11r""" 
     2Definition 
     3---------- 
    24Calculates the scattering from a fractal structure with a primary building 
    35block of core-shell spheres, as opposed to just homogeneous spheres in 
    4 the fractal model. 
    5 This model could find use for aggregates of coated particles, or aggregates 
    6 of vesicles. 
    7  
    8 Definition 
    9 ---------- 
     6the fractal model. It is an extension of the well known Teixeira\ [#teixeira]_ 
     7fractal model replacing the $P(q)$ of a solid sphere with that of a core-shell 
     8sphere. This model could find use for aggregates of coated particles, or 
     9aggregates of vesicles for example. 
    1010 
    1111.. math:: 
    1212 
    13     I(q) = \text{background} + P(q)S(q) 
     13    I(q) = P(q)S(q) + \text{background} 
    1414 
    15 The form factor $P(q)$ is that from core_shell model with $bkg$ = 0 
    16  
     15Where $P(q)$ is the core-shell form factor and $S(q)$ is the 
     16Teixeira\ [#teixeira]_ fractal structure factor both of which are given again 
     17below: 
    1718 
    1819.. math:: 
    1920 
    20     P(q)=\frac{scale}{V_s}\left[3V_c(\rho_c-\rho_s) 
     21    P(q) &= \frac{\phi}{V_s}\left[3V_c(\rho_c-\rho_s) 
    2122    \frac{\sin(qr_c)-qr_c\cos(qr_c)}{(qr_c)^3}+ 
    2223    3V_s(\rho_s-\rho_{solv}) 
    2324    \frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3}\right]^2 
    2425 
     26    S(q) &= 1 + \frac{D_f\ \Gamma\!(D_f-1)}{[1+1/(q\xi)^2]^{(D_f-1)/2}}  
     27    \frac{\sin[(D_f-1)\tan^{-1}(q\xi)]}{(qr_s)^{D_f}} 
    2528 
    26 while the fractal structure factor $S(q)$ is 
     29where $\phi$ is the volume fraction of particles, $V_s$ is the volume of the 
     30whole particle, $V_c$ is the volume of the core, $\rho_c$, $\rho_s$, and 
     31$\rho_{solv}$ are the scattering length densities of the core, shell, and 
     32solvent respectively, $r_c$ and $r_s$ are the radius of the core and the radius 
     33of the whole particle respectively, $D_f$ is the fractal dimension, and |xi| the 
     34correlation length. 
     35  
     36Polydispersity of radius and thickness are also provided for. 
    2737 
    28 .. math:: 
    29  
    30     S(q) = \frac{D_f\Gamma(D_f-1)\sin((D_f-1)\tan^{-1}(q\xi))} 
    31     {(qr_c)^{D_f}\left(1+\frac{1}{q^2\xi ^2} \right)^{\frac{D_f-1}{2}}} 
    32  
    33 where $D_f$ = fractal_dim, |xi| = cor_length, $r_c$ = (core) radius, and 
    34 $scale$ = volume fraction. 
    35  
    36 The fractal structure is as documented in the fractal model. 
    37 Polydispersity of radius and thickness is provided for. 
    38  
    39 For 2D data: The 2D scattering intensity is calculated in the same way as 1D, 
    40 where the $q$ vector is defined as 
     38This model does not allow for anisotropy and thus the 2D scattering intensity 
     39is calculated in the same way as 1D, where the $q$ vector is defined as 
    4140 
    4241.. math:: 
     
    4443    q = \sqrt{q_x^2 + q_y^2} 
    4544 
     45Our model is derived from the form factor calculations implemented in IGOR 
     46macros by the NIST Center for Neutron Research\ [#Kline]_ 
     47 
    4648References 
    4749---------- 
    4850 
    49 See the core_shell and fractal model descriptions 
     51.. [#teixeira] J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 
     52.. [#Kline]  S R Kline, *J Appl. Cryst.*, 39 (2006) 895 
     53 
     54Authorship and Verification 
     55---------------------------- 
     56 
     57* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     58* **Last Modified by:** Paul Butler and Paul Kienzle **on:** November 27, 2016 
     59* **Last Reviewed by:** Paul Butler and Paul Kienzle **on:** November 27, 2016 
    5060 
    5161""" 
     
    5464 
    5565name = "fractal_core_shell" 
    56 title = "" 
    57 description = """ 
     66title = "Scattering from a fractal structure formed from core shell spheres" 
     67description = """\ 
     68    Model for fractal aggregates of core-shell primary particles. It is based on 
     69    the Teixeira model for the S(q) of a fractal * P(q) for a core-shell sphere 
    5870 
    59 """ 
     71    radius =  the radius of the core 
     72    thickness = thickness of the shell 
     73    thick_layer = thickness of a layer 
     74    sld_core = the SLD of the core 
     75    sld_shell = the SLD of the shell 
     76    sld_solvent = the SLD of the solvent 
     77    volfraction = volume fraction of core-shell particles 
     78    fractal_dim = fractal dimension 
     79    cor_length = correlation length of the fractal like aggretates 
     80    """ 
    6081category = "shape-independent" 
    6182 
     
    6384#   ["name", "units", default, [lower, upper], "type","description"], 
    6485parameters = [ 
    65     ["radius",      "Ang",        60.0, [0, inf],    "volume", "Sphere core radius"], 
    66     ["thickness",   "Ang",        10.0, [0, inf],    "volume", "Sphere shell thickness"], 
     86    ["radius",      "Ang",        60.0, [0.0, inf],  "volume", "Sphere core radius"], 
     87    ["thickness",   "Ang",        10.0, [0.0, inf],  "volume", "Sphere shell thickness"], 
    6788    ["sld_core",    "1e-6/Ang^2", 1.0,  [-inf, inf], "sld",    "Sphere core scattering length density"], 
    6889    ["sld_shell",   "1e-6/Ang^2", 2.0,  [-inf, inf], "sld",    "Sphere shell scattering length density"], 
    6990    ["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"]] 
     91    ["volfraction", "",           1.0,  [0.0, inf],  "",       "Volume fraction of building block spheres"], 
     92    ["fractal_dim",    "",        2.0,  [0.0, 6.0],  "",       "Fractal dimension"], 
     93    ["cor_length",  "Ang",      100.0,  [0.0, inf],  "",       "Correlation length of fractal-like aggregates"], 
     94] 
    7395# pylint: enable=bad-whitespace, line-too-long 
    7496 
    75 source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c", "fractal_core_shell.c"] 
     97source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "lib/core_shell.c", 
     98          "lib/fractal_sq.c", "fractal_core_shell.c"] 
    7699 
    77100demo = dict(scale=0.05, 
     
    100123        @param thickness: shell thickness 
    101124    """ 
    102     whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3) 
    103     core = 4.0 * pi / 3.0 * radius * radius * radius 
     125    whole = 4.0/3.0 * pi * (radius + thickness)**3 
     126    core = 4.0/3.0 * pi * radius**3 
    104127    return whole, whole-core 
    105128 
    106129tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
    107          [{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704], 
     130         [{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704]] 
    108131 
    109          # The SasView test result was 0.00169, with a background of 0.001 
    110          [{'radius': 60.0, 
    111            'thickness': 10.0, 
    112            'sld_core': 1.0, 
    113            'sld_shell': 2.0, 
    114            'sld_solvent': 3.0, 
    115            'background': 0.0 
    116           }, 0.4, 0.00070126]] 
     132#         # The SasView test result was 0.00169, with a background of 0.001 
     133#         # They are however wrong as we now know.  IGOR might be a more 
     134#         # appropriate source. Otherwise will just have to assume this is now 
     135#         # correct and self generate a correct answer for the future. Until we 
     136#         # figure it out leave the tests commented out 
     137#         [{'radius': 60.0, 
     138#           'thickness': 10.0, 
     139#           'sld_core': 1.0, 
     140#           'sld_shell': 2.0, 
     141#           'sld_solvent': 3.0, 
     142#           'background': 0.0 
     143#          }, 0.015211, 692.84]] 
  • sasmodels/models/fuzzy_sphere.py

    r9a4811a r925ad6e  
    8181# pylint: enable=bad-whitespace,line-too-long 
    8282 
    83 source = ["lib/sph_j1c.c"] 
     83source = ["lib/sas_3j1x_x.c"] 
    8484 
    8585# No volume normalization despite having a volume parameter 
    8686# This should perhaps be volume normalized? 
    8787form_volume = """ 
    88     return 1.333333333333333*M_PI*radius*radius*radius; 
     88    return M_4PI_3*cube(radius); 
    8989    """ 
    9090 
    9191Iq = """ 
    9292    const double qr = q*radius; 
    93     const double bes = sph_j1c(qr); 
     93    const double bes = sas_3j1x_x(qr); 
    9494    const double qf = q*fuzziness; 
    9595    const double fq = bes * (sld - sld_solvent) * form_volume(radius) * exp(-0.5*qf*qf); 
  • sasmodels/models/guinier_porod.py

    ra807206 rcdcebf1  
    44and dimensionality of scattering objects, including asymmetric objects 
    55such as rods or platelets, and shapes intermediate between spheres 
    6 and rods or between rods and platelets. 
     6and rods or between rods and platelets, and overcomes some of the  
     7deficiencies of the (Beaucage) Unified_Power_Rg model (see Hammouda, 2010). 
    78 
    89Definition 
     
    5960--------- 
    6061 
    61 A Guinier, G Fournet, *Small-Angle Scattering of X-Rays*, 
    62 John Wiley and Sons, New York, (1955) 
     62B Hammouda, *A new Guinier-Porod model, J. Appl. Cryst.*, (2010), 43, 716-719 
    6363 
    64 O Glatter, O Kratky, *Small-Angle X-Ray Scattering*, Academic Press (1982) 
    65 Check out Chapter 4 on Data Treatment, pages 155-156. 
     64B Hammouda, *Analysis of the Beaucage model, J. Appl. Cryst.*, (2010), 43, 1474-1478 
     65 
    6666""" 
    6767 
  • sasmodels/models/hayter_msa.c

    r0bef47b r4962519  
    2323        double SIdiam, diam, Kappa, cs, IonSt; 
    2424        double  Perm, Beta; 
    25         double pi, charge; 
     25        double charge; 
    2626        int ierr; 
    2727         
    28         pi = M_PI; 
    29  
    3028        diam=2*radius_effective;                //in A 
    3129 
     
    3836        charge=zz*Elcharge;             //in Coulomb (C) 
    3937        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^3 
     38        Vp=M_4PI_3*cube(SIdiam/2.0);    //in m^3 
    4139        cs=csalt*6.022E23*1.0E3;        //# salt molecules/m^3 
    4240         
     
    5048        Kappa=sqrt(2*Beta*IonSt/Perm);     //Kappa calc from Ionic strength 
    5149                                                                           //   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)); 
    5351         
    5452        //         Finally set up dimensionless parameters  
  • sasmodels/models/hollow_cylinder.c

    r0d6e865 r592343f  
    11double form_volume(double radius, double thickness, double length); 
    2  
    32double Iq(double q, double radius, double thickness, double length, double sld, 
    43        double solvent_sld); 
     
    98 
    109// From Igor library 
    11 static double hollow_cylinder_scaling( 
    12     double integrand, double delrho, double volume) 
     10static double 
     11_hollow_cylinder_scaling(double integrand, double delrho, double volume) 
    1312{ 
    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; 
    2514} 
    2615 
    2716 
    28 static double _hollow_cylinder_kernel( 
    29     double q, double radius, double thickness, double length, double dum) 
     17static double 
     18_hollow_cylinder_kernel(double q, 
     19    double radius, double thickness, double length, double sin_val, double cos_val) 
    3020{ 
    31     const double qs = q*sqrt(1.0-dum*dum); 
    32     const double lam1 = sas_J1c((radius+thickness)*qs); 
    33     const double lam2 = sas_J1c(radius*qs); 
     21    const double qs = q*sin_val; 
     22    const double lam1 = sas_2J1x_x((radius+thickness)*qs); 
     23    const double lam2 = sas_2J1x_x(radius*qs); 
    3424    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 = sas_J0(radius*qs) 
     26    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs) 
    3627    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 = sas_sinx_x(0.5*q*length*cos_val); 
     29    return psi*t2; 
     30} 
     31 
     32double 
     33form_volume(double radius, double thickness, double length) 
     34{ 
     35    double v_shell = M_PI*length*(square(radius+thickness) - radius*radius); 
     36    return v_shell; 
    3937} 
    4038 
    4139 
    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) 
     40double 
     41Iq(double q, double radius, double thickness, double length, 
     42    double sld, double solvent_sld) 
    4543{ 
    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 
    5446 
    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    } 
    5955 
    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); 
    7359} 
    7460 
     61double 
     62Iqxy(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); 
    7570 
    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); 
    8073} 
    8174 
    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 limits 
    88     double summ,answer,delrho;                  //running tally of integration 
    89     double norm,volume; //final calculation variables 
    90  
    91     lower = 0.0; 
    92     upper = 1.0;                //limits of numerical integral 
    93  
    94     summ = 0.0;                 //initialize intergral 
    95     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 } 
  • sasmodels/models/hollow_rectangular_prism.c

    ra807206 r1e7b0db0  
    22double Iq(double q, double sld, double solvent_sld, double length_a,  
    33          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); 
    64 
    75double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    86{ 
    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; 
    119    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; 
    1412    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; 
    1614    double vol_shell = vol_total - vol_core; 
    1715    return vol_shell; 
     
    2624    double thickness) 
    2725{ 
    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 
    2939     
    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++) { 
    3542 
    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); 
    4546 
    46     double sumi = 0.0; 
    47      
    48     for(int i=0; i<nordi; i++) { 
     47        const double termC1 = sas_sinx_x(q * c_half * cos(theta)); 
     48        const double termC2 = sas_sinx_x(q * (c_half-thickness)*cos(theta)); 
    4949 
    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++) { 
    5152 
    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); 
    6256 
    6357            // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0 
    6458 
    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 = sas_sinx_x(q * a_half * sin_theta * sin_phi); 
     60            const double termA2 = sas_sinx_x(q * (a_half-thickness) * sin_theta * sin_phi); 
    6961 
    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 = sas_sinx_x(q * b_half * sin_theta * cos_phi); 
     63            const double termB2 = sas_sinx_x(q * (b_half-thickness) * sin_theta * cos_phi); 
    7464 
    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; 
    7867 
    79                 sumj += Gauss76Wt[j] * (AP*AP); 
     68            inner_sum += Gauss76Wt[j] * square(AP1-AP2); 
     69        } 
     70        inner_sum *= 0.5 * (v2b-v2a); 
    8071 
    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); 
    8673    } 
    87  
    88     double answer = 0.5*(v1b-v1a)*sumi; 
     74    outer_sum *= 0.5*(v1b-v1a); 
    8975 
    9076    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 
    9177    // 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; 
    9379 
    9480    // 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; 
    9682 
    9783    // 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; 
    10285} 
    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 } 
  • sasmodels/models/hollow_rectangular_prism.py

    ra807206 rab2aea8  
    1313the difference of the amplitudes of two massive parallelepipeds 
    1414differing in their outermost dimensions in each direction by the 
    15 same length increment :math:`2\Delta` (Nayuk, 2012). 
     15same length increment $2\Delta$ (Nayuk, 2012). 
    1616 
    1717As in the case of the massive parallelepiped model (:ref:`rectangular-prism`), 
     
    5858 
    5959.. 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} 
    6262 
    63 where $\rho_{\text{p}}$ is the scattering length of the parallelepiped, 
    64 $\rho_{\text{solvent}}$ is the scattering length of the solvent, 
     63where $\rho_\text{p}$ is the scattering length of the parallelepiped, 
     64$\rho_\text{solvent}$ is the scattering length of the solvent, 
    6565and (if the data are in absolute units) *scale* represents the volume fraction 
    6666(which is unitless). 
     
    148148# parameters for demo 
    149149demo = dict(scale=1, background=0, 
    150             sld=6.3e-6, sld_solvent=1.0e-6, 
     150            sld=6.3, sld_solvent=1.0, 
    151151            length_a=35, b2a_ratio=1, c2a_ratio=1, thickness=1, 
    152152            length_a_pd=0.1, length_a_pd_n=10, 
  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    ra807206 rab2aea8  
    22double Iq(double q, double sld, double solvent_sld, double length_a,  
    33          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); 
    64 
    75double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    86{ 
    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); 
    1210    return vol_shell; 
    1311} 
     
    2018    double c2a_ratio) 
    2119{ 
    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; 
    2725 
    2826   //Integration limits to use in Gaussian quadrature 
    29     double v1a = 0.0; 
    30     double v1b = 0.5 * M_PI;  //theta integration limits 
    31     double v2a = 0.0; 
    32     double v2b = 0.5 * M_PI;  //phi integration limits 
     27    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 
    3331     
    34     //Order of integration 
    35     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 ); 
    3735 
    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); 
    4140 
    42             double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );  
    43          
    4441        // 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); 
    4744 
    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 ); 
    5148 
    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 
    5456            // 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); 
    5759 
    5860            // 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 ); 
    6163 
    62             sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT); 
     64            inner_sum += Gauss76Wt[j] * square(AL+AT); 
     65        } 
    6366 
    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; 
    6969    } 
    7070 
    71     double answer = 0.5*(v1b-v1a)*sumi; 
     71    outer_sum *= 0.5*(v1b-v1a); 
    7272 
    7373    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 
    7474    // 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; 
    7676 
    7777    // 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); 
    7979 
    8080    // Convert from [1e-12 A-1] to [cm-1] 
     
    8282 
    8383    return answer; 
    84      
    8584} 
    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 } 
  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    ra807206 rab2aea8  
    33r""" 
    44 
    5 This model provides the form factor, *P(q)*, for a hollow rectangular 
     5This model provides the form factor, $P(q)$, for a hollow rectangular 
    66prism with infinitely thin walls. It computes only the 1D scattering, not the 2D. 
    77 
     
    1414 
    1515Assuming a hollow parallelepiped with infinitely thin walls, edge lengths 
    16 :math:`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 
     16$A \le B \le C$ and presenting an orientation with respect to the 
     17scattering vector given by $\theta$ and $\phi$, where $\theta$ is the angle 
     18between 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) 
     20and the $y$ axis, the form factor is given by 
    2121 
    2222.. 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 
    2526 
    2627where 
    2728 
    2829.. math:: 
    29   V = 2AB + 2AC + 2BC 
    3030 
    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} 
    3939 
    4040and 
    4141 
    4242.. 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) } 
    4546                     {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) } 
    4849                     {q \, \sin\phi \, \sin\theta} 
    4950 
     
    5152 
    5253.. math:: 
    53   I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} - \rho_{\mbox{solvent}})^2 \times P(q) 
    5454