Changeset 6245e42 in sasmodels


Ignore:
Timestamp:
Nov 18, 2016 9:21:25 AM (5 years ago)
Author:
GitHub <noreply@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
f80f334
Parents:
0e0e867 (diff), f8f5a37 (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.
git-author:
Andrew Jackson <andrew.jackson@…> (11/18/16 09:21:25)
git-committer:
GitHub <noreply@…> (11/18/16 09:21:25)
Message:

Merge pull request #12 from SasView?/ticket776

Ticket776

Files:
1 added
73 edited

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 r8c65a33  
    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.update(pars) 
     1066    context.update((k,v) for k,v in presets.items() if isinstance(v, float)) 
     1067    for k, v in presets.items(): 
     1068        if not isinstance(v, float) and not k.endswith('_type'): 
     1069            presets[k] = eval(v, context) 
     1070    context.update(presets) 
     1071    context.update((k,v) for k,v in presets2.items() if isinstance(v, float)) 
     1072    for k, v in presets2.items(): 
     1073        if not isinstance(v, float) and not k.endswith('_type'): 
     1074            presets2[k] = eval(v, context) 
     1075 
     1076    # update parameters with presets 
     1077    pars.update(presets)  # set value after random to control value 
     1078    pars2.update(presets2)  # set value after random to control value 
     1079    #import pprint; pprint.pprint(model_info) 
     1080 
     1081    same_model = name == name2 and pars == pars 
    8971082    if len(engines) == 0: 
    898         engines.extend(['single', 'double']) 
     1083        if same_model: 
     1084            engines.extend(['single', 'double']) 
     1085        else: 
     1086            engines.extend(['single', 'single']) 
    8991087    elif len(engines) == 1: 
    900         if engines[0][0] == 'double': 
     1088        if not same_model: 
     1089            engines.append(engines[0]) 
     1090        elif engines[0] == 'double': 
    9011091            engines.append('single') 
    9021092        else: 
     
    9051095        del engines[2:] 
    9061096 
    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 
    9091097    use_sasview = any(engine == 'sasview' and count > 0 
    9101098                      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) 
    9381099    if use_sasview: 
    9391100        constrain_new_to_old(model_info, pars) 
     1101        constrain_new_to_old(model_info2, pars2) 
     1102 
    9401103    if opts['show_pars']: 
    941         print(str(parlist(model_info, pars, opts['is2d']))) 
     1104        if not same_model: 
     1105            print("==== %s ====="%model_info.name) 
     1106            print(str(parlist(model_info, pars, opts['is2d']))) 
     1107            print("==== %s ====="%model_info2.name) 
     1108            print(str(parlist(model_info2, pars2, opts['is2d']))) 
     1109        else: 
     1110            print(str(parlist(model_info, pars, opts['is2d']))) 
    9421111 
    9431112    # Create the computational engines 
     
    9481117        base = None 
    9491118    if n2: 
    950         comp = make_engine(model_info, data, engines[1], opts['cutoff']) 
     1119        comp = make_engine(model_info2, data, engines[1], opts['cutoff']) 
    9511120    else: 
    9521121        comp = None 
     
    9551124    # Remember it all 
    9561125    opts.update({ 
    957         'name'      : name, 
    958         'def'       : model_info, 
    959         'n1'        : n1, 
    960         'n2'        : n2, 
    961         'presets'   : presets, 
    962         'pars'      : pars, 
    9631126        'data'      : data, 
     1127        'name'      : [name, name2], 
     1128        'def'       : [model_info, model_info2], 
     1129        'count'     : [n1, n2], 
     1130        'presets'   : [presets, presets2], 
     1131        'pars'      : [pars, pars2], 
    9641132        'engines'   : [base, comp], 
    9651133    }) 
     
    9761144    from .generate import view_html_from_info 
    9771145    app = wx.App() if wx.GetApp() is None else None 
    978     view_html_from_info(opts['def']) 
     1146    view_html_from_info(opts['def'][0]) 
    9791147    if app: app.MainLoop() 
    9801148 
     
    9881156    from bumps.names import FitProblem  # type: ignore 
    9891157    from bumps.gui.app_frame import AppFrame  # type: ignore 
     1158    from bumps.gui import signal 
    9901159 
    9911160    is_mac = "cocoa" in wx.version() 
    9921161    # Create an app if not running embedded 
    9931162    app = wx.App() if wx.GetApp() is None else None 
    994     problem = FitProblem(Explore(opts)) 
     1163    model = Explore(opts) 
     1164    problem = FitProblem(model) 
    9951165    frame = AppFrame(parent=None, title="explore", size=(1000,700)) 
    9961166    if not is_mac: frame.Show() 
     
    9981168    frame.panel.Layout() 
    9991169    frame.panel.aui.Split(0, wx.TOP) 
     1170    def reset_parameters(event): 
     1171        model.revert_values() 
     1172        signal.update_parameters(problem) 
     1173    frame.Bind(wx.EVT_TOOL, reset_parameters, frame.ToolBar.GetToolByPos(1)) 
    10001174    if is_mac: frame.Show() 
    10011175    # If running withing an app, start the main loop 
     
    10151189        config_matplotlib() 
    10161190        self.opts = opts 
    1017         model_info = opts['def'] 
    1018         pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 
     1191        p1, p2 = opts['pars'] 
     1192        m1, m2 = opts['def'] 
     1193        self.fix_p2 = m1 != m2 or p1 != p2 
     1194        model_info = m1 
     1195        pars, pd_types = bumps_model.create_parameters(model_info, **p1) 
    10191196        # Initialize parameter ranges, fixing the 2D parameters for 1D data. 
    10201197        if not opts['is2d']: 
    1021             for p in model_info.parameters.user_parameters(is2d=False): 
     1198            for p in model_info.parameters.user_parameters({}, is2d=False): 
    10221199                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 
    10231200                    k = p.name+ext 
     
    10301207 
    10311208        self.pars = pars 
     1209        self.starting_values = dict((k, v.value) for k, v in pars.items()) 
    10321210        self.pd_types = pd_types 
    10331211        self.limits = None 
     1212 
     1213    def revert_values(self): 
     1214        for k, v in self.starting_values.items(): 
     1215            self.pars[k].value = v 
     1216 
     1217    def model_update(self): 
     1218        pass 
    10341219 
    10351220    def numpoints(self): 
     
    10621247        pars = dict((k, v.value) for k, v in self.pars.items()) 
    10631248        pars.update(self.pd_types) 
    1064         self.opts['pars'] = pars 
    1065         limits = compare(self.opts, limits=self.limits) 
     1249        self.opts['pars'][0] = pars 
     1250        if not self.fix_p2: 
     1251            self.opts['pars'][1] = pars 
     1252        result = run_models(self.opts) 
     1253        limits = plot_models(self.opts, result, limits=self.limits) 
    10661254        if self.limits is None: 
    10671255            vmin, vmax = limits 
    10681256            self.limits = vmax*1e-7, 1.3*vmax 
     1257            import pylab; pylab.clf() 
     1258            plot_models(self.opts, result, limits=self.limits) 
    10691259 
    10701260 
  • sasmodels/conversion_table.py

    r5f1acda rfa6d6fc  
    143143        } 
    144144    ], 
    145     "core_shell_ellipsoid": [ 
     145    "core_shell_ellipsoid:1": [ 
    146146        "CoreShellEllipsoidModel", 
    147147        { 
     148            "sld_core": "sld_core", 
     149            "sld_shell": "sld_shell", 
     150            "sld_solvent": "sld_solvent", 
     151            "radius_equat_core": "equat_core", 
     152            "x_core": "polar_core", 
     153            "thick_shell": "equat_shell", 
     154            "x_polar_shell": "polar_shell", 
     155            "theta": "axis_theta", 
    148156            "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" 
    157157        } 
    158158    ], 
     
    160160        "CoreShellEllipsoidXTModel", 
    161161        { 
    162             "phi": "axis_phi", 
    163162            "sld_core": "sld_core", 
     163            "sld_shell": "sld_shell", 
     164            "sld_solvent": "sld_solvent", 
     165            "radius_equat_core": "equat_core", 
     166            "thick_shell": "T_shell", 
    164167            "x_core": "X_core", 
    165             "sld_solvent": "sld_solvent", 
    166             "thick_shell": "T_shell", 
    167168            "x_polar_shell": "XpolarShell", 
    168169            "theta": "axis_theta", 
    169             "sld_shell": "sld_shell" 
     170            "phi": "axis_phi", 
    170171        } 
    171172    ], 
     
    173174        "CSParallelepipedModel", 
    174175        { 
     176            "sld_core": "sld_pcore", 
     177            "sld_a": "sld_rimA", 
     178            "sld_b": "sld_rimB", 
     179            "sld_c": "sld_rimC", 
     180            "sld_solvent": "sld_solv", 
     181            "length_a": "shortA", 
     182            "length_b": "midB", 
     183            "length_c": "longC", 
     184            "thick_rim_a": "rimA", 
     185            "thick_rim_c": "rimC", 
     186            "thick_rim_b": "rimB", 
     187            "theta": "parallel_theta", 
    175188            "phi": "parallel_phi", 
    176189            "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" 
    189190        } 
    190191    ], 
     
    240241        "DABModel", 
    241242        { 
    242             "length": "length" 
     243            "cor_length": "length" 
    243244        } 
    244245    ], 
     
    257258        "EllipticalCylinderModel", 
    258259        { 
     260            "axis_ratio": "r_ratio", 
     261            "radius_minor": "r_minor", 
     262            "sld": "sldCyl", 
     263            "sld_solvent": "sldSolv", 
     264            "theta": "cyl_theta", 
    259265            "phi": "cyl_phi", 
    260266            "psi": "cyl_psi", 
    261             "theta": "cyl_theta", 
    262             "sld": "sldCyl", 
    263             "axis_ratio": "r_ratio", 
    264             "sld_solvent": "sldSolv" 
    265267        } 
    266268    ], 
     
    299301        "FractalCoreShellModel", 
    300302        { 
     303            "sld_core": "core_sld", 
    301304            "sld_shell": "shell_sld", 
    302305            "sld_solvent": "solvent_sld", 
    303             "sld_core": "core_sld" 
     306            "radius": "radius", 
     307            "thickness": "thickness", 
     308            "fractal_dim": "frac_dim", 
     309            "cor_length": "cor_length", 
     310            "volfraction": "volfraction", 
    304311        } 
    305312    ], 
     
    326333        "PeakGaussModel", 
    327334        { 
    328             "sigma": "B" 
     335            "peak_pos": "q0", 
     336            "sigma": "B", 
    329337        } 
    330338    ], 
     
    334342            "rg": "radius", 
    335343            "lorentz_scale": "lScale", 
    336             "fractal_dim": "FractalExp", 
     344            "guinier_scale": "gScale", 
     345            "fractal_dim": "scale", 
    337346            "cor_length": "zeta", 
    338             "guinier_scale": "gScale" 
    339347        } 
    340348    ], 
     
    350358            "s": "dim", 
    351359            "rg": "rg", 
    352             "m": "m", 
     360            "porod_exp": "m", 
    353361            "scale": "scale", 
    354362            "background": "background" 
     
    358366        "HardsphereStructure", 
    359367        { 
    360             "radius_effective_pd": "effect_radius_pd", 
     368            "scale": "scale_factor", 
    361369            "radius_effective": "effect_radius", 
    362             "radius_effective_pd_n": "effect_radius_pd_n" 
    363370        } 
    364371    ], 
     
    366373        "HayterMSAStructure", 
    367374        { 
    368             "salt_concentration": "saltconc", 
    369             "radius_effective_pd": "effect_radius_pd", 
     375            "scale": "scale_factor", 
    370376            "radius_effective": "effect_radius", 
    371             "radius_effective_pd_n": "effect_radius_pd_n" 
     377            "volfraction": "volfraction", 
     378            "charge": "charge", 
     379            "temperature": "temperature", 
     380            "concentration_salt": "saltconc", 
     381            "dielectconst": "dielectconst", 
    372382        } 
    373383    ], 
     
    375385        "HollowCylinderModel", 
    376386        { 
     387            "sld": "sldCyl", 
     388            "sld_solvent": "sldSolv", 
     389            "radius": "core_radius", 
     390            "thickness": "radius", 
     391            "length": "length", 
     392            "theta": "axis_theta", 
    377393            "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" 
    386394        } 
    387395    ], 
     
    389397        "RectangularHollowPrismModel", 
    390398        { 
     399            "sld": "sldPipe", 
     400            "sld_solvent": "sldSolv", 
     401            "length_a": "short_side", 
    391402            "b2a_ratio": "b2a_ratio", 
    392             "length_a": "short_side", 
    393             "sld": "sldPipe", 
    394             "length_c": "c2a_ratio", 
    395             "sld_solvent": "sldSolv", 
    396             "thickness": "thickness" 
     403            "c2a_ratio": "c2a_ratio", 
     404            "thickness": "thickness", 
    397405        } 
    398406    ], 
     
    401409        { 
    402410            "sld": "sldPipe", 
     411            "sld_solvent": "sldSolv", 
     412            "length_a": "short_side", 
    403413            "b2a_ratio": "b2a_ratio", 
    404             "length_a": "short_side", 
    405             "length_c": "c2a_ratio", 
    406             "sld_solvent": "sldSolv" 
     414            "c2a_ratio": "c2a_ratio", 
    407415        } 
    408416    ], 
     
    428436        "LamellarPSHGModel", 
    429437        { 
     438            "sld": "sld_tail", 
     439            "sld_head": "sld_head", 
     440            "sld_solvent": "sld_solvent", 
     441            "length_tail": "deltaT", 
     442            "length_head": "deltaH", 
     443            "d_spacing": "spacing", 
    430444            "Caille_parameter": "caille", 
    431445            "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" 
    437446        } 
    438447    ], 
     
    441450        { 
    442451            "sld": "sld_bi", 
     452            "sld_solvent": "sld_sol", 
     453            "thickness": "delta", 
     454            "d_spacing": "spacing", 
    443455            "Caille_parameter": "caille", 
    444456            "Nlayers": "N_plates", 
    445             "sld_solvent": "sld_sol", 
    446             "thickness": "delta" 
    447457        } 
    448458    ], 
     
    451461        { 
    452462            "sld": "sld_layer", 
     463            "sld_solvent": "sld_solvent", 
     464            "thickness": "thickness", 
     465            "d_spacing": "spacing", 
    453466            "sigma_d": "pd_spacing", 
    454             "sld_solvent": "sld_solvent" 
     467            "Nlayers": "Nlayers", 
    455468        } 
    456469    ], 
     
    500513        { 
    501514            "rg": "rg", 
    502             "scale": "scale", 
    503             "background": "background" 
     515            "i_zero": "scale", 
     516            "background": "background", 
    504517        } 
    505518    ], 
     
    576589            "rg": "rg", 
    577590            "polydispersity": "poly_m", 
     591            "i_zero": "scale", 
     592            "background": "background", 
     593        } 
     594    ], 
     595    "polymer_excl_volume": [ 
     596        "PolymerExclVolume", 
     597        { 
     598            "rg": "rg", 
     599            "scale": "scale", 
     600            "background": "background", 
     601            "porod_exp": "m" 
     602        } 
     603    ], 
     604    "polymer_micelle": [ 
     605        "MicelleSphCoreModel", 
     606        { 
     607            "sld_corona": "rho_corona", 
     608            "sld_solvent": "rho_solv", 
     609            "sld_core": "rho_core", 
     610            "ndensity": "ndensity", 
     611            "v_core": "v_core", 
     612            "v_corona": "v_corona", 
     613            "radius_core": "radius_core", 
     614            "rg": "radius_gyr", 
     615            "d_penetration": "d_penetration", 
     616            "n_aggreg": "n_aggreg", 
     617        } 
     618    ], 
     619    "porod": [ 
     620        "PorodModel", 
     621        { 
    578622            "scale": "scale", 
    579623            "background": "background" 
    580624        } 
    581625    ], 
    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     ], 
    606626    "power_law": [ 
    607627        "PowerLawAbsModel", 
     
    616636        { 
    617637            "scale": "scale", 
    618             "solvent_sld": "sld_solvent", 
     638            "sld_solvent": "sld_solvent", 
    619639            "thickness": "thickness", 
    620640            "beta": "beta", 
     
    643663        { 
    644664            "sld": "sldPipe", 
     665            "length_a": "short_side", 
    645666            "b2a_ratio": "b2a_ratio", 
    646             "length_a": "short_side", 
    647             "length_c": "c2a_ratio", 
     667            "c2a_ratio": "c2a_ratio", 
    648668            "sld_solvent": "sldSolv" 
    649669        } 
     
    716736        "SquareWellStructure", 
    717737        { 
    718             "radius_effective_pd": "effect_radius_pd", 
     738            "scale": "scale_factor", 
    719739            "radius_effective": "effect_radius", 
    720             "radius_effective_pd_n": "effect_radius_pd_n" 
     740            "wellwidth": "wellwidth", 
     741            "welldepth": "welldepth", 
    721742        } 
    722743    ], 
     
    729750            "theta": "axis_theta", 
    730751            "sld_solvent": "solvent_sld", 
    731             "n_stacking": "n_stacking" 
     752            "n_stacking": "n_stacking", 
     753            "thick_layer": "layer_thick", 
     754            "thick_core": "core_thick", 
    732755        } 
    733756    ], 
     
    742765        "StickyHSStructure", 
    743766        { 
    744             "radius_effective_pd": "effect_radius_pd", 
     767            "scale": "scale_factor", 
    745768            "radius_effective": "effect_radius", 
    746             "radius_effective_pd_n": "effect_radius_pd_n" 
    747769        } 
    748770    ], 
     
    758780        "TeubnerStreyModel", 
    759781        { 
    760             "a2": "scale" 
     782            # Note: parameters are completely rewritten in convert.py 
     783            "volfraction_a": "volfraction_a", 
     784            "sld_a": "sld_a", 
     785            "sld_b": "sld_b", 
     786            "d": "d", 
     787            "xi": "xi", 
    761788        } 
    762789    ], 
  • sasmodels/convert.py

    r51241113 rf8f5a37  
    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"), 
     
    5655    ] 
    5756 
    58 def _convert_pars(pars, mapping): 
    59     """ 
    60     Rename the parameters and any associated polydispersity attributes. 
    61     """ 
    62     newpars = pars.copy() 
    63     for new, old in mapping.items(): 
    64         if old == new: continue 
    65         for underscore, dot in PD_DOT: 
    66             if old+dot in newpars: 
    67                 if new is not None: 
    68                     newpars[new+underscore] = pars[old+dot] 
    69                 del newpars[old+dot] 
    70     return newpars 
    71  
    72 def convert_model(name, pars): 
    73     """ 
    74     Convert model from old style parameter names to new style. 
    75     """ 
    76     _, _ = name, pars # lint 
    77     raise NotImplementedError 
    78     # need to load all new models in order to determine old=>new 
    79     # model name mapping 
    80  
    81 def _unscale(par, scale): 
     57def _rescale(par, scale): 
    8258    return [pk*scale for pk in par] if isinstance(par, list) else par*scale 
    8359 
    84 def _is_sld(modelinfo, id): 
     60def _is_sld(model_info, id): 
     61    """ 
     62    Return True if parameter is a magnetic magnitude or SLD parameter. 
     63    """ 
    8564    if id.startswith('M0:'): 
    8665        return True 
    87     if (id.endswith('_pd') or id.endswith('_pd_n') or id.endswith('_pd_nsigma') 
    88             or id.endswith('_pd_width') or id.endswith('_pd_type')): 
     66    if id.startswith('volfraction') or id.startswith('radius_effective'): 
    8967        return False 
    90     for p in modelinfo.parameters.call_parameters: 
     68    if '_pd' in id or '.' in id: 
     69        return False 
     70    for p in model_info.parameters.call_parameters: 
    9171        if p.id == id: 
    9272            return p.type == 'sld' 
    9373    # check through kernel parameters in case it is a named as a vector 
    94     for p in modelinfo.parameters.kernel_parameters: 
     74    for p in model_info.parameters.kernel_parameters: 
    9575        if p.id == id: 
    9676            return p.type == 'sld' 
    9777    raise ValueError("unknown parameter %r in conversion"%id) 
    9878 
    99 def _unscale_sld(modelinfo, pars): 
    100     """ 
    101     rescale all sld parameters in the new model definition by 1e6 so the 
     79def _rescale_sld(model_info, pars, scale): 
     80    """ 
     81    rescale all sld parameters in the new model definition by *scale* so the 
    10282    numbers are nicer.  Relies on the fact that all sld parameters in the 
    103     new model definition end with sld. 
    104     """ 
    105     return dict((id, (_unscale(v, 1e-6) if _is_sld(modelinfo, id) else v)) 
     83    new model definition end with sld.  For backward conversion use 
     84    *scale=1e-6*.  For forward conversion use *scale=1e6*. 
     85    """ 
     86    return dict((id, (_rescale(v, scale) if _is_sld(model_info, id) else v)) 
    10687                for id, v in pars.items()) 
    10788 
    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 
    122  
    123 def _revert_pars(pars, mapping): 
    124     """ 
    125     Rename the parameters and any associated polydispersity attributes. 
    126     """ 
    127     newpars = pars.copy() 
    128  
    129     for new, old in mapping.items(): 
    130         for underscore, dot in PD_DOT: 
    131             if old and old+underscore == new+dot: 
    132                 continue 
    133             if new+underscore in newpars: 
    134                 if old is not None: 
    135                     newpars[old+dot] = pars[new+underscore] 
    136                 del newpars[new+underscore] 
    137     for k in list(newpars.keys()): 
    138         for underscore, dot in PD_DOT[1:]:  # skip "" => "" 
    139             if k.endswith(underscore): 
    140                 newpars[k[:-len(underscore)]+dot] = newpars[k] 
    141                 del newpars[k] 
    142     return newpars 
    143  
    144 def revert_name(model_info): 
    145     oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    146     return oldname 
    14789 
    14890def _get_translation_table(model_info): 
     
    162104    return translation 
    163105 
     106# ========= FORWARD CONVERSION sasview 3.x => sasmodels =========== 
     107def _dot_pd_to_underscore_pd(par): 
     108    if par.endswith(".width"): 
     109        return par[:-6]+"_pd" 
     110    elif par.endswith(".type"): 
     111        return par[:-5]+"_pd_type" 
     112    elif par.endswith(".nsigmas"): 
     113        return par[:-8]+"_pd_nsigma" 
     114    elif par.endswith(".npts"): 
     115        return par[:-5]+"_pd_n" 
     116    else: 
     117        return par 
     118 
     119def _pd_to_underscores(pars): 
     120    return dict((_dot_pd_to_underscore_pd(k), v) for k, v in pars.items()) 
     121 
     122def _convert_name(conv_dict, pars): 
     123    """ 
     124    Renames parameter values (upper, lower, etc) to v4.0 names 
     125    :param conv_dict: conversion dictionary mapping new name : old name 
     126    :param pars: parameters to convert 
     127    :return: 
     128    """ 
     129    new_pars = {} 
     130    i = 0 
     131    j = 0 
     132    for key_par, value_par in pars.iteritems(): 
     133        j += 1 
     134        for key_conv, value_conv in conv_dict.iteritems(): 
     135            if re.search(value_conv, key_par): 
     136                new_pars[key_par.replace(value_conv, key_conv)] = value_par 
     137                i += 1 
     138                break 
     139            elif re.search("background", key_par) or re.search("scale", key_par): 
     140                new_pars[key_par] = value_par 
     141                i += 1 
     142                break 
     143        if i != j: 
     144            new_pars[key_par] = value_par 
     145            i += 1 
     146    return new_pars 
     147 
     148def _convert_pars(pars, mapping): 
     149    """ 
     150    Rename the parameters and any associated polydispersity attributes. 
     151    """ 
     152    newpars = pars.copy() 
     153    for new, old in mapping.items(): 
     154        if old == new: continue 
     155        if old is None: continue 
     156        for underscore, dot in PD_DOT: 
     157            source = old+dot 
     158            if source in newpars: 
     159                if new is not None: 
     160                    target = new+dot 
     161                else: 
     162                    target = None 
     163                if source != target: 
     164                    if target: 
     165                        newpars[target] = pars[old+dot] 
     166                    del newpars[source] 
     167    return newpars 
     168 
     169 
     170def _conversion_target(model_name): 
     171    """ 
     172    Find the sasmodel name which translates into the sasview name. 
     173 
     174    Note: *CoreShellEllipsoidModel* translates into *core_shell_ellipsoid:1*. 
     175    This is necessary since there is only one variant in sasmodels for the 
     176    two variants in sasview. 
     177    """ 
     178    for sasmodels_name, [sasview_name, _] in CONVERSION_TABLE.items(): 
     179        if sasview_name == model_name: 
     180            return sasmodels_name 
     181    return None 
     182 
     183 
     184def _hand_convert(name, oldpars): 
     185    if name == 'core_shell_parallelepiped': 
     186        # Make sure pd on rim parameters defaults to zero 
     187        # ... probably not necessary. 
     188        oldpars['rimA.width'] = 0.0 
     189        oldpars['rimB.width'] = 0.0 
     190        oldpars['rimC.width'] = 0.0 
     191    elif name == 'core_shell_ellipsoid:1': 
     192        # Reverse translation (from new to old), from core_shell_ellipsoid.c 
     193        #    equat_shell = equat_core + thick_shell 
     194        #    polar_core = equat_core * x_core 
     195        #    polar_shell = equat_core * x_core + thick_shell*x_polar_shell 
     196        # Forward translation (from old to new), inverting reverse translation: 
     197        #    thick_shell = equat_shell - equat_core 
     198        #    x_core = polar_core / equat_core 
     199        #    x_polar_shell = (polar_shell - polar_core)/(equat_shell - equat_core) 
     200        # Auto translation (old <=> new) happens after hand_convert 
     201        #    equat_shell <=> thick_shell 
     202        #    polar_core <=> x_core 
     203        #    polar_shell <=> x_polar_shell 
     204        # So... 
     205        equat_core, equat_shell = oldpars['equat_core'], oldpars['equat_shell'] 
     206        polar_core, polar_shell = oldpars['polar_core'], oldpars['polar_shell'] 
     207        oldpars['equat_shell'] = equat_shell - equat_core 
     208        oldpars['polar_core'] = polar_core / equat_core 
     209        oldpars['polar_shell'] = (polar_shell-polar_core)/(equat_shell-equat_core) 
     210    elif name == 'hollow_cylinder': 
     211        # now uses radius and thickness 
     212        thickness = oldpars['radius'] - oldpars['core_radius'] 
     213        oldpars['radius'] = thickness 
     214        if 'radius.width' in oldpars: 
     215            pd = oldpars['radius.width']*oldpars['radius']/thickness 
     216            oldpars['radius.width'] = pd 
     217    elif name == 'pearl_necklace': 
     218        pass 
     219        #_remove_pd(oldpars, 'num_pearls', name) 
     220        #_remove_pd(oldpars, 'thick_string', name) 
     221    elif name == 'polymer_micelle': 
     222        if 'ndensity' in oldpars: 
     223            oldpars['ndensity'] /= 1e15 
     224        if 'ndensity.lower' in oldpars: 
     225            oldpars['ndensity.lower'] /= 1e15 
     226        if 'ndensity.upper' in oldpars: 
     227            oldpars['ndensity.upper'] /= 1e15 
     228    elif name == 'rpa': 
     229        # convert scattering lengths from femtometers to centimeters 
     230        for p in "L1", "L2", "L3", "L4": 
     231            if p in oldpars: 
     232                oldpars[p] /= 1e-13 
     233            if p + ".lower" in oldpars: 
     234                oldpars[p + ".lower"] /= 1e-13 
     235            if p + ".upper" in oldpars: 
     236                oldpars[p + ".upper"] /= 1e-13 
     237    elif name == 'spherical_sld': 
     238        oldpars["CONTROL"] += 1 
     239    elif name == 'teubner_strey': 
     240        # basically undoing the entire Teubner-Strey calculations here. 
     241        #    drho = (sld_a - sld_b) 
     242        #    k = 2.0*math.pi*xi/d 
     243        #    a2 = (1.0 + k**2)**2 
     244        #    c1 = 2.0 * xi**2 * (1.0 - k**2) 
     245        #    c2 = xi**4 
     246        #    prefactor = 8.0*math.pi*phi*(1.0-phi)*drho**2*c2/xi 
     247        #    scale = 1e-4*prefactor 
     248        #    oldpars['scale'] = a2/scale 
     249        #    oldpars['c1'] = c1/scale 
     250        #    oldpars['c2'] = c2/scale 
     251 
     252        # need xi, d, sld_a, sld_b, phi=volfraction_a 
     253        # assume contrast is 1.0e-6, scale=1, background=0 
     254        sld_a, sld_b = 1.0, 0. 
     255        drho = sld_a - sld_b 
     256 
     257        # find xi 
     258        p_scale = oldpars['scale'] 
     259        p_c1 = oldpars['c1'] 
     260        p_c2= oldpars['c2'] 
     261        i_1 = 0.5*p_c1/p_c2 
     262        i_2 = math.sqrt(math.fabs(p_scale/p_c2)) 
     263        i_3 = 2/(i_1 + i_2) 
     264        xi = math.sqrt(math.fabs(i_3)) 
     265 
     266        # find d from xi 
     267        k = math.sqrt(math.fabs(1 - 0.5*p_c1/p_c2*xi**2)) 
     268        d = 2*math.pi*xi/k 
     269 
     270        # solve quadratic phi (1-phi) = xi/(1e-4 8 pi drho^2 c2) 
     271        # favour volume fraction in [0, 0.5] 
     272        c = xi / (1e-4 * 8.0 * math.pi * drho**2 * p_c2) 
     273        phi = 0.5 - math.sqrt(0.25 - c) 
     274 
     275        # scale sld_a by 1e-6 because the translator will scale it back 
     276        oldpars.update(volfraction_a=phi, xi=xi, d=d, sld_a=sld_a*1e-6, 
     277                       sld_b=sld_b, scale=1.0) 
     278        oldpars.pop('c1') 
     279        oldpars.pop('c2') 
     280 
     281    return oldpars 
     282 
     283def convert_model(name, pars, use_underscore=False): 
     284    """ 
     285    Convert model from old style parameter names to new style. 
     286    """ 
     287    newname = _conversion_target(name) 
     288    if newname is None: 
     289        return name, pars 
     290    if ':' in newname:   # core_shell_ellipsoid:1 
     291        model_info = load_model_info(newname[:-2]) 
     292        # Know that the table exists and isn't multiplicity so grab it directly 
     293        # Can't use _get_translation_table since that will return the 'bare' 
     294        # version. 
     295        translation = CONVERSION_TABLE[newname] 
     296    else: 
     297        model_info = load_model_info(newname) 
     298        translation = _get_translation_table(model_info) 
     299    newpars = _hand_convert(newname, pars.copy()) 
     300    newpars = _convert_name(translation, newpars) 
     301    newpars = _convert_pars(newpars, translation) 
     302    if not model_info.structure_factor: 
     303        newpars = _rescale_sld(model_info, newpars, 1e6) 
     304    newpars.setdefault('scale', 1.0) 
     305    newpars.setdefault('background', 0.0) 
     306    if use_underscore: 
     307        newpars = _pd_to_underscores(newpars) 
     308    return newname, newpars 
     309 
     310 
     311# ========= BACKWARD CONVERSION sasmodels => sasview 3.x =========== 
     312 
     313def _revert_pars(pars, mapping): 
     314    """ 
     315    Rename the parameters and any associated polydispersity attributes. 
     316    """ 
     317    newpars = pars.copy() 
     318 
     319    for new, old in mapping.items(): 
     320        for underscore, dot in PD_DOT: 
     321            if old and old+underscore == new+dot: 
     322                continue 
     323            if new+underscore in newpars: 
     324                if old is not None: 
     325                    newpars[old+dot] = pars[new+underscore] 
     326                del newpars[new+underscore] 
     327    for k in list(newpars.keys()): 
     328        for underscore, dot in PD_DOT[1:]:  # skip "" => "" 
     329            if k.endswith(underscore): 
     330                newpars[k[:-len(underscore)]+dot] = newpars[k] 
     331                del newpars[k] 
     332    return newpars 
     333 
     334def revert_name(model_info): 
     335    oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
     336    return oldname 
     337 
     338def _remove_pd(pars, key, name): 
     339    """ 
     340    Remove polydispersity from the parameter list. 
     341 
     342    Note: operates in place 
     343    """ 
     344    # Bumps style parameter names 
     345    width = pars.pop(key+".width", 0.0) 
     346    n_points = pars.pop(key+".npts", 0) 
     347    if width != 0.0 and n_points != 0: 
     348        warnings.warn("parameter %s not polydisperse in sasview %s"%(key, name)) 
     349    pars.pop(key+".nsigmas", None) 
     350    pars.pop(key+".type", None) 
     351    return pars 
     352 
    164353def _trim_vectors(model_info, pars, oldpars): 
    165354    _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
     
    180369        composition_type, parts = model_info.composition 
    181370        if composition_type == 'product': 
    182             translation = {'scale':'scale_factor'} 
    183             translation.update(_get_translation_table(parts[0])) 
     371            translation = _get_translation_table(parts[0]) 
     372            # structure factor models include scale:scale_factor mapping 
    184373            translation.update(_get_translation_table(parts[1])) 
    185374        else: 
     
    187376    else: 
    188377        translation = _get_translation_table(model_info) 
    189     oldpars = _revert_pars(_unscale_sld(model_info, pars), translation) 
     378    oldpars = _revert_pars(_rescale_sld(model_info, pars, 1e-6), translation) 
    190379    oldpars = _trim_vectors(model_info, pars, oldpars) 
    191380 
     
    215404        if name in MODELS_WITHOUT_VOLFRACTION: 
    216405            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) 
    242406        elif name == 'core_multi_shell': 
    243407            # kill extra shells 
     
    252416        elif name == 'core_shell_parallelepiped': 
    253417            _remove_pd(oldpars, 'rimA', name) 
    254         elif name in ['mono_gauss_coil', 'poly_gauss_coil']: 
    255             del oldpars['i_zero'] 
     418            _remove_pd(oldpars, 'rimB', name) 
     419            _remove_pd(oldpars, 'rimC', name) 
     420        elif name == 'hollow_cylinder': 
     421            # now uses radius and thickness 
     422            thickness = oldpars['core_radius'] 
     423            oldpars['radius'] += thickness 
     424            oldpars['radius.width'] *= thickness/oldpars['radius'] 
     425        #elif name in ['mono_gauss_coil', 'poly_gauss_coil']: 
     426        #    del oldpars['i_zero'] 
    256427        elif name == 'onion': 
    257428            oldpars.pop('n_shells', None) 
     429        elif name == 'pearl_necklace': 
     430            _remove_pd(oldpars, 'num_pearls', name) 
     431            _remove_pd(oldpars, 'thick_string', name) 
     432        elif name == 'polymer_micelle': 
     433            if 'ndensity' in oldpars: 
     434                oldpars['ndensity'] *= 1e15 
    258435        elif name == 'rpa': 
    259436            # convert scattering lengths from femtometers to centimeters 
     
    272449                for k in "Kab,Kac,Kad".split(','): 
    273450                    oldpars.pop(k, None) 
     451        elif name == 'spherical_sld': 
     452            oldpars["CONTROL"] -= 1 
     453            # remove polydispersity from shells 
     454            for k in range(1, 11): 
     455                _remove_pd(oldpars, 'thick_flat'+str(k), 'thickness') 
     456                _remove_pd(oldpars, 'thick_inter'+str(k), 'interface') 
     457            # remove extra shells 
     458            for k in range(int(pars['n_shells']), 11): 
     459                oldpars.pop('sld_flat'+str(k), 0) 
     460                oldpars.pop('thick_flat'+str(k), 0) 
     461                oldpars.pop('thick_inter'+str(k), 0) 
     462                oldpars.pop('func_inter'+str(k), 0) 
     463                oldpars.pop('nu_inter'+str(k), 0) 
     464        elif name == 'stacked_disks': 
     465            _remove_pd(oldpars, 'n_stacking', name) 
     466        elif name == 'teubner_strey': 
     467            # basically redoing the entire Teubner-Strey calculations here. 
     468            volfraction = oldpars.pop('volfraction_a') 
     469            xi = oldpars.pop('xi') 
     470            d = oldpars.pop('d') 
     471            sld_a = oldpars.pop('sld_a') 
     472            sld_b = oldpars.pop('sld_b') 
     473            drho = 1e6*(sld_a - sld_b)  # conversion autoscaled these 
     474            k = 2.0*math.pi*xi/d 
     475            a2 = (1.0 + k**2)**2 
     476            c1 = 2.0 * xi**2 * (1.0 - k**2) 
     477            c2 = xi**4 
     478            prefactor = 8.0*math.pi*volfraction*(1.0-volfraction)*drho**2*c2/xi 
     479            scale = 1e-4*prefactor 
     480            oldpars['scale'] = a2/scale 
     481            oldpars['c1'] = c1/scale 
     482            oldpars['c2'] = c2/scale 
    274483 
    275484    #print("convert from",list(sorted(pars))) 
     
    315524        if name in MODELS_WITHOUT_VOLFRACTION: 
    316525            pars['volfraction'] = 1 
    317         if name == 'pearl_necklace': 
    318             pars['string_thickness_pd_n'] = 0 
    319             pars['number_of_pearls_pd_n'] = 0 
     526        if name == 'core_multi_shell': 
     527            pars['n'] = min(math.ceil(pars['n']), 4) 
     528        elif name == 'gel_fit': 
     529            pars['scale'] = 1 
    320530        elif name == 'line': 
    321531            pars['scale'] = 1 
    322532            pars['background'] = 0 
     533        elif name == 'mono_gauss_coil': 
     534            pars['scale'] = 1 
     535        elif name == 'onion': 
     536            pars['n_shells'] = math.ceil(pars['n_shells']) 
     537        elif name == 'pearl_necklace': 
     538            pars['string_thickness_pd_n'] = 0 
     539            pars['number_of_pearls_pd_n'] = 0 
     540        elif name == 'poly_gauss_coil': 
     541            pars['scale'] = 1 
    323542        elif name == 'rpa': 
    324543            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']) 
    333544        elif name == 'spherical_sld': 
    334545            pars['n_shells'] = math.ceil(pars['n_shells']) 
     
    339550                pars['thickness%d_pd_n'%k] = 0 
    340551                pars['interface%d_pd_n'%k] = 0 
    341  
     552        elif name == 'teubner_strey': 
     553            pars['scale'] = 1 
     554            if pars['volfraction_a'] > 0.5: 
     555                pars['volfraction_a'] = 1.0 - pars['volfraction_a'] 
     556        elif name == 'unified_power_Rg': 
     557            pars['level'] = int(pars['level']) 
     558 
     559def _check_one(name, seed=None): 
     560    """ 
     561    Generate a random set of parameters for *name*, and check that they can 
     562    be converted back to SasView 3.x and forward again to sasmodels.  Raises 
     563    an error if the parameters are changed. 
     564    """ 
     565    from . import compare 
     566 
     567    model_info = load_model_info(name) 
     568 
     569    old_name = revert_name(model_info) 
     570    if old_name is None: 
     571        return 
     572 
     573    pars = compare.get_pars(model_info, use_demo=False) 
     574    pars = compare.randomize_pars(model_info, pars, seed=seed) 
     575    if name == "teubner_strey": 
     576        # T-S model is underconstrained, so fix the assumptions. 
     577        pars['sld_a'], pars['sld_b'] = 1.0, 0.0 
     578    compare.constrain_pars(model_info, pars) 
     579    constrain_new_to_old(model_info, pars) 
     580    old_pars = revert_pars(model_info, pars) 
     581    new_name, new_pars = convert_model(old_name, old_pars, use_underscore=True) 
     582    if 1: 
     583        print("==== %s in ====="%name) 
     584        print(str(compare.parlist(model_info, pars, True))) 
     585        print("==== %s ====="%old_name) 
     586        for k, v in sorted(old_pars.items()): 
     587            print(k, v) 
     588        print("==== %s out ====="%new_name) 
     589        print(str(compare.parlist(model_info, new_pars, True))) 
     590    assert name==new_name, "%r != %r"%(name, new_name) 
     591    for k, v in new_pars.items(): 
     592        assert k in pars, "%s: %r appeared from conversion"%(name, k) 
     593        if isinstance(v, float): 
     594            assert abs(v-pars[k])<=abs(1e-12*v), "%s: %r  %s != %s"%(name, k, v, pars[k]) 
     595        else: 
     596            assert v == pars[k], "%s: %r  %s != %s"%(name, k, v, pars[k]) 
     597    for k, v in pars.items(): 
     598        assert k in pars, "%s: %r not converted"%(name, k) 
     599 
     600def test_backward_forward(): 
     601    from .core import list_models 
     602    for name in list_models('all'): 
     603        L = lambda: _check_one(name, seed=1) 
     604        L.description = name 
     605        yield L 
  • 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/generate.py

    r234c532 r7fcdc9f  
    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 r218cdbc  
    148148inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    149149 
     150#if 0 
     151#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
     152    SINCOS(phi*M_PI_180, sn, cn); \ 
     153    q = sqrt(qx*qx + qy*qy); \ 
     154    cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * cos(theta*M_PI_180));  \ 
     155    sn = sqrt(1 - cn*cn); \ 
     156    } while (0) 
     157#else 
     158// SasView 3.x definition of orientation 
     159#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
     160    SINCOS(theta*M_PI_180, sn, cn); \ 
     161    q = sqrt(qx*qx + qy*qy);\ 
     162    cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \ 
     163    sn = sqrt(1 - cn*cn); \ 
     164    } while (0) 
     165#endif 
     166 
     167#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
     168    q = sqrt(qx*qx + qy*qy); \ 
     169    const double qxhat = qx/q; \ 
     170    const double qyhat = qy/q; \ 
     171    double sin_theta, cos_theta; \ 
     172    double sin_phi, cos_phi; \ 
     173    double sin_psi, cos_psi; \ 
     174    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
     175    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
     176    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
     177    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ 
     178    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ 
     179    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
     180    } while (0) 
  • sasmodels/kernelcl.py

    r6e5b2a7 r7fcdc9f  
    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 r3a48772  
    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_J1c(q*radius*sin_alpha); 
     52    const double si = sinc(q*half_length*cos_alpha); 
     53    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
     54    const double Aq = bell_fq + cyl_fq; 
     55    return Aq; 
     56} 
     57 
    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/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/binary_hard_sphere.c

    re481a39 r4f79d94  
    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     //} 
    7660    sc1 = sph_j1c(qr1); 
    7761    sc2 = sph_j1c(qr2); 
    78     b1 = r1*r1*r1*(rho1-rhos)*4.0/3.0*M_PI*sc1; 
    79     b2 = r2*r2*r2*(rho2-rhos)*4.0/3.0*M_PI*sc2; 
     62    b1 = r1*r1*r1*(rho1-rhos)*M_4PI_3*sc1; 
     63    b2 = r2*r2*r2*(rho2-rhos)*M_4PI_3*sc2; 
    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/capped_cylinder.c

    r0d6e865 r3a48772  
    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_J1c(q*radius*sin_alpha); 
     57    const double si = sinc(q*half_length*cos_alpha); 
     58    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
     59    const double Aq = cap_Fq + cyl_Fq; 
     60    return Aq; 
    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/core_shell_bicelle.c

    r0d6e865 r5bddd89  
    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 
    5857    be1 = sas_J1c(besarg1); 
     
    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_cylinder.c

    r0d6e865 r9aa4881  
    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 * sinc(siarg) * sas_J1c(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_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 r73e08ae  
    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"], 
  • sasmodels/models/core_shell_parallelepiped.c

    r2222134 r14838a3  
    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 = sinc(mu_proj * sin_uu * a_scaled); 
     90            const double si2 = sinc(mu_proj * cos_uu); 
     91            const double si3 = sinc(mu_proj * sin_uu * ta); 
     92            const double si4 = sinc(mu_proj * cos_uu * tb); 
     93 
    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 = sinc(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 = sinc(0.5*q*length_a*cos_val_a); 
     163    double siB = sinc(0.5*q*length_b*cos_val_b); 
     164    double siC = sinc(0.5*q*length_c*cos_val_c); 
     165    double siAt = sinc(0.5*q*ta*cos_val_a); 
     166    double siBt = sinc(0.5*q*tb*cos_val_b); 
     167    double siCt = sinc(0.5*q*tc*cos_val_c); 
    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 r14838a3  
    164164# parameters for demo 
    165165demo = 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, 
     166            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6, 
    168167            length_a=35, length_b=75, length_c=400, 
    169168            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10, 
     
    177176            theta_pd=10, theta_pd_n=1, 
    178177            phi_pd=10, phi_pd_n=1, 
    179             psi_pd=10, psi_pd_n=10) 
     178            psi_pd=10, psi_pd_n=1) 
    180179 
    181180qx, 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 r4962519  
    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 rb829b16  
    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_J1c(qr*sn) * sinc(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/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 r73e08ae  
    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 = sph_j1c(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/elliptical_cylinder.c

    ra807206 r251f54b  
    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_J1c(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 = sinc(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_J1c(q*r); 
     76    const double si = sinc(q*0.5*length*cos_val); 
     77    const double Aq = be * si; 
     78    const double delrho = sld - solvent_sld; 
    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/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/flexible_cylinder_elliptical.c

    rb66d38e r92ce163  
    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_J1c(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 r217590b  
    1 double Iq(double q, 
    2           double volfraction, 
    3           double radius, 
    4           double fractal_dim, 
    5           double cor_length, 
    6           double sld_block, 
    7           double sld_solvent); 
     1#define INVALID(p) (p.fractal_dim < 0.0) 
    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)*sph_j1c(q*radius)); 
     17 
     18    // scale to units cm-1 sr-1 (assuming data on absolute scale) 
     19    //    convert I(1/A) to (1/cm)  => 1e8 * I(q) 
     20    //    convert rho^2 in 10^-6 1/A to 1/A  => 1e-12 * I(q) 
     21    //    combined: 1e-4 * I(q) 
     22 
     23    return 1.e-4 * volfraction * sq * pq; 
    5324} 
    5425 
  • sasmodels/models/fractal.py

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

    ra807206 r217590b  
    1 double form_volume(double radius, double thickness); 
    2  
    3 double Iq(double q, 
    4           double radius, 
    5           double thickness, 
    6           double core_sld, 
    7           double shell_sld, 
    8           double solvent_sld, 
    9           double volfraction, 
    10           double fractal_dim, 
    11           double cor_length); 
    12  
    13 double form_volume(double radius, double thickness) 
     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    const double sq = fractal_sq(q, radius, fractal_dim, cor_length); 
     19    const double pq = core_shell_kernel(q, radius, thickness, 
     20                                        core_sld, shell_sld, solvent_sld); 
    2721 
    28  
    29    double intensity = core_shell_kernel(q, 
    30                               radius, 
    31                               thickness, 
    32                               core_sld, 
    33                               shell_sld, 
    34                               solvent_sld); 
    35     //calculate S(q) 
    36     double frac_1 = fractal_dim-1.0; 
    37     double qr = q*radius; 
    38  
    39     double t1 = fractal_dim*sas_gamma(frac_1)*sin(frac_1*atan(q*cor_length)); 
    40     double t2 = (1.0 + 1.0/(q*cor_length)/(q*cor_length)); 
    41     double t3 = pow(qr, fractal_dim) * pow(t2, (frac_1/2.0)); 
    42     double sq = t1/t3; 
    43     sq += 1.0; 
    44  
    45     return sq*intensity*volfraction; 
     22    // Note: core_shell_kernel already performs the 1e-4 unit conversion 
     23    return volfraction * sq * pq; 
    4624} 
    4725 
  • sasmodels/models/fractal_core_shell.py

    ra807206 r217590b  
    6363#   ["name", "units", default, [lower, upper], "type","description"], 
    6464parameters = [ 
    65     ["radius",      "Ang",        60.0, [0, inf],    "volume", "Sphere core radius"], 
    66     ["thickness",   "Ang",        10.0, [0, inf],    "volume", "Sphere shell thickness"], 
     65    ["radius",      "Ang",        60.0, [0.0, inf],  "volume", "Sphere core radius"], 
     66    ["thickness",   "Ang",        10.0, [0.0, inf],  "volume", "Sphere shell thickness"], 
    6767    ["sld_core",    "1e-6/Ang^2", 1.0,  [-inf, inf], "sld",    "Sphere core scattering length density"], 
    6868    ["sld_shell",   "1e-6/Ang^2", 2.0,  [-inf, inf], "sld",    "Sphere shell scattering length density"], 
    6969    ["sld_solvent", "1e-6/Ang^2", 3.0,  [-inf, inf], "sld",    "Solvent scattering length density"], 
    70     ["volfraction", "",           1.0,  [0, inf],    "",       "Volume fraction of building block spheres"], 
    71     ["fractal_dim",    "",           2.0,  [-inf, inf], "",       "Fractal dimension"], 
    72     ["cor_length",  "Ang",      100.0,  [0, inf],    "",       "Correlation length of fractal-like aggregates"]] 
     70    ["volfraction", "",           1.0,  [0.0, inf],  "",       "Volume fraction of building block spheres"], 
     71    ["fractal_dim",    "",        2.0,  [0.0, 6.0],  "",       "Fractal dimension"], 
     72    ["cor_length",  "Ang",      100.0,  [0.0, inf],  "",       "Correlation length of fractal-like aggregates"], 
     73] 
    7374# pylint: enable=bad-whitespace, line-too-long 
    7475 
    75 source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c", "fractal_core_shell.c"] 
     76source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c", 
     77          "lib/fractal_sq.c", "fractal_core_shell.c"] 
    7678 
    7779demo = dict(scale=0.05, 
     
    100102        @param thickness: shell thickness 
    101103    """ 
    102     whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3) 
    103     core = 4.0 * pi / 3.0 * radius * radius * radius 
     104    whole = 4.0/3.0 * pi * (radius + thickness)**3 
     105    core = 4.0/3.0 * pi * radius**3 
    104106    return whole, whole-core 
    105107 
  • sasmodels/models/fuzzy_sphere.py

    r9a4811a r3a48772  
    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 
  • 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 rf8f0991  
    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); 
     21    const double qs = q*sin_val; 
    3222    const double lam1 = sas_J1c((radius+thickness)*qs); 
    3323    const double lam2 = sas_J1c(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 = J0(radius*qs) 
     26    //Note: lim_{radius -> 0} psi = sas_J1c(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 = sinc(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 r6f676fb  
    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 = sinc(q * c_half * cos(theta)); 
     48        const double termC2 = sinc(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 = sinc(q * a_half * sin_theta * sin_phi); 
     60            const double termA2 = sinc(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 = sinc(q * b_half * sin_theta * cos_phi); 
     63            const double termB2 = sinc(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 
    55 where *V* is the volume of the rectangular prism, :math:`\rho_{\mbox{p}}` 
    56 is the scattering length of the parallelepiped, :math:`\rho_{\mbox{solvent}}` 
     55  I(q) = \text{scale} \times V \times (\rho_\text{p} - \rho_\text{solvent})^2 \times P(q) 
     56 
     57where $V$ is the volume of the rectangular prism, $\rho_\text{p}$ 
     58is the scattering length of the parallelepiped, $\rho_\text{solvent}$ 
    5759is the scattering length of the solvent, and (if the data are in absolute 
    5860units) *scale* represents the volume fraction (which is unitless). 
     
    127129# parameters for demo 
    128130demo = dict(scale=1, background=0, 
    129             sld=6.3e-6, sld_solvent=1.0e-6, 
     131            sld=6.3, sld_solvent=1.0, 
    130132            length_a=35, b2a_ratio=1, c2a_ratio=1, 
    131133            length_a_pd=0.1, length_a_pd_n=10, 
  • sasmodels/models/lamellar_stack_paracrystal.c

    r0bef47b r4962519  
    6767        double Snq; 
    6868 
    69         Snq = an/( (double)Nlayers*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) ); 
     69        Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) ); 
    7070         
    7171        return(Snq); 
  • sasmodels/models/lib/Si.c

    r0278e3f rf719764  
    33double Si(double x) 
    44{ 
    5     double out; 
     5    if (x >= M_PI*6.2/4.0) { 
     6        const double xxinv = 1./(x*x); 
     7        // Explicitly writing factorial values triples the speed of the calculation 
     8        const double out_cos = (((-720.*xxinv + 24.)*xxinv - 2.)*xxinv + 1.)/x; 
     9        const double out_sin = (((-5040.*xxinv + 120.)*xxinv - 6.)*xxinv + 1)*xxinv; 
    610 
    7     if (x >= M_PI*6.2/4.0){ 
    8         double out_sin = 0.0; 
    9         double out_cos = 0.0; 
    10         out = M_PI/2.0; 
    11  
     11        double sin_x, cos_x; 
     12        SINCOS(x, sin_x, cos_x); 
     13        return M_PI_2 - cos_x*out_cos - sin_x*out_sin; 
     14    } else { 
     15        const double xx = x*x; 
    1216        // Explicitly writing factorial values triples the speed of the calculation 
    13         out_cos = 1./x - 2./pow(x,3) + 24./pow(x,5) - 720./pow(x,7); 
    14         out_sin = 1./pow(x,2) - 6./pow(x,4) + 120./pow(x,6) - 5040./pow(x,8); 
    15  
    16         out -= cos(x) * out_cos; 
    17         out -= sin(x) * out_sin; 
    18         return out; 
     17        return (((((-1./439084800.*xx 
     18            + 1./3265920.)*xx 
     19            - 1./35280.)*xx 
     20            + 1./600.)*xx 
     21            - 1./18.)*xx 
     22            + 1.)*x; 
    1923    } 
    20  
    21     // Explicitly writing factorial values triples the speed of the calculation 
    22     out = x - pow(x, 3)/18. + pow(x,5)/600. - pow(x,7)/35280. + pow(x,9)/3265920. - pow(x,11)/439084800.; 
    23  
    24     return out; 
    2524} 
  • sasmodels/models/lib/core_shell.c

    r7d4b2ae r3a48772  
    1717    const double core_contrast = core_sld - shell_sld; 
    1818    const double core_bes = sph_j1c(core_qr); 
    19     const double core_volume = 4.0 * M_PI / 3.0 * radius * radius * radius; 
     19    const double core_volume = M_4PI_3 * cube(radius); 
    2020    double f = core_volume * core_bes * core_contrast; 
    2121 
     
    2424    const double shell_contrast = shell_sld - solvent_sld; 
    2525    const double shell_bes = sph_j1c(shell_qr); 
    26     const double shell_volume = 4.0 * M_PI / 3.0 * pow((radius + thickness), 3); 
     26    const double shell_volume = M_4PI_3 * cube(radius + thickness); 
    2727    f += shell_volume * shell_bes * shell_contrast; 
    2828    return f * f * 1.0e-4; 
  • sasmodels/models/lib/gfn.c

    r177c1a1 r3a48772  
    1212{ 
    1313    // local variables 
    14     const double pi43=4.0/3.0*M_PI; 
    1514    const double aa = crmaj; 
    1615    const double bb = crmin; 
     
    2019    //const double siq = (uq == 0.0 ? 1.0 : 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq); 
    2120    const double siq = sph_j1c(uq); 
    22     const double vc = pi43*aa*aa*bb; 
     21    const double vc = M_4PI_3*aa*aa*bb; 
    2322    const double gfnc = siq*vc*delpc; 
    2423 
    2524    const double ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx)); 
    2625    const double ut= sqrt(ut2)*qq; 
    27     const double vt = pi43*trmaj*trmaj*trmin; 
     26    const double vt = M_4PI_3*trmaj*trmaj*trmin; 
    2827    //const double sit = (ut == 0.0 ? 1.0 : 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut); 
    2928    const double sit = sph_j1c(ut); 
     
    3332    const double result = tgfn*tgfn; 
    3433 
    35     return (result); 
     34    return result; 
    3635} 
  • sasmodels/models/lib/sas_gamma.c

    r1596de3 r97f9b46  
    55We use gamma definition Gamma(t + 1) = t * Gamma(t) to compute 
    66to function for values lower than 1.0. Namely Gamma(t) = 1/t * Gamma(t + 1) 
     7For t < 0, we use Gamma(t) = pi / ( Gamma(1 - t) * sin(pi * t) ) 
    78*/ 
    89 
     
    137138 
    138139 
    139 inline double sas_gamma( double x) { return tgamma(x+1)/x; } 
     140inline double sas_gamma(double x) 
     141{ 
     142    // Note: the builtin tgamma can give slow and unreliable results for x<1. 
     143    // The following transform extends it to zero and to negative values. 
     144    // It should return NaN for zero and negative integers but doesn't. 
     145    // The accuracy is okay but not wonderful for negative numbers, maybe 
     146    // one or two digits lost in the calculation. If higher accuracy is 
     147    // needed, you could test the following loop: 
     148    //    double norm = 1.; 
     149    //    while (x<1.) { norm*=x; x+=1.; } 
     150    //    return tgamma(x)/norm; 
     151    return (x<0. ? M_PI/tgamma(1.-x)/sin(M_PI*x) : tgamma(x+1)/x); 
     152} 
  • sasmodels/models/lib/wrc_cyl.c

    rba32cdd r92ce163  
    1414    //return t; 
    1515 
    16     return pow( (1.0 + (x/3.12)*(x/3.12) + 
    17          (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 
     16    return pow(1.0+square(x/3.12)+cube(x/8.67), 0.176/3.0); 
    1817} 
    1918 
     
    2322{ 
    2423    const double r = b/L; 
    25     return (L*b/6.0) * 
    26            (1.0 - r*1.5  + 1.5*r*r - 0.75*r*r*r*(1.0 - exp(-2.0/r))); 
     24    return (L*b/6.0) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r)))); 
     25 
    2726} 
    2827 
     
    4140} 
    4241 
    43 static inline double 
     42static double 
    4443sech_WR(double x) 
    4544{ 
     
    5150{ 
    5251    double C; 
    53     const double onehalf = 1.0/2.0; 
    5452 
    5553    if( L/b > 10.0) { 
     
    8684 
    8785    const double t2 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    88          Rg02/b2))*((1.0 + onehalf*(((-1.0) - 
     86         Rg02/b2))*((1.0 + 0.5*(((-1.0) - 
    8987         tanh((-C4 + Rgb/C5))))))); 
    9088 
     
    112110 
    113111    const double t9 = (2.0*b4*(((2.0*q0*Rg2)/b - 
    114          (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + onehalf*(((-1.0) - 
     112         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + 0.5*(((-1.0) - 
    115113         tanh(((-C4) + Rgb)/C5)))))); 
    116114 
    117115    const double t10 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    118          Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) + 
     116         Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    119117         Rgb)/C5)))))); 
    120118 
     
    131129 
    132130    const double t14 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    133           Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) + 
     131          Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    134132          Rgb)/C5)))))); 
    135133 
     
    140138 
    141139    double yy = (pow(q0,p1)*(((-((b*M_PI)/(L*q0))) +t1/L +t2/(q04*Rg22) + 
    142         onehalf*t3*t4)) + (t5*((pow(q0,(p1 - p2))* 
     140        0.5*t3*t4)) + (t5*((pow(q0,(p1 - p2))* 
    143141        (((-pow(q0,(-p1)))*(((b2*M_PI)/(L*q02) +t6/L +t7/(2.0*C5) - 
    144         t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + onehalf*t11*t12)) - 
     142        t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 0.5*t11*t12)) - 
    145143        b*p1*pow(q0,((-1.0) - p1))*(((-((b*M_PI)/(L*q0))) + t13/L + 
    146         t14/(q04*Rg22) + onehalf*t15*((1.0 + tanh(((-C4) + 
     144        t14/(q04*Rg22) + 0.5*t15*((1.0 + tanh(((-C4) + 
    147145        Rgb)/C5))))))))))); 
    148146 
     
    154152{ 
    155153    double C; 
    156     const double onehalf = 1.0/2.0; 
    157154 
    158155    if( L/b > 10.0) { 
     
    201198    const double t5 = (2.0*b4*(((2.0*q0*Rg2)/b - 
    202199         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))* 
    203          ((1.0 + onehalf*(((-1.0) - tanh(((-C4) + 
     200         ((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    204201         Rgb)/C5))))))/(q04*Rg22); 
    205202 
    206203    const double t6 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    207          Rg02/b2))*((1.0 + onehalf*(((-1) - tanh(((-C4) + 
     204         Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 
    208205         Rgb)/C5))))))/(q05*Rg22); 
    209206 
     
    219216 
    220217    const double t10 = (2.0*b4*(((-1) + pow((double)M_E,(-(Rg02/b2))) + 
    221           Rg02/b2))*((1.0 + onehalf*(((-1) - tanh(((-C4) + 
     218          Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 
    222219          Rgb)/C5))))))/(q04*Rg22); 
    223220 
    224221    const double yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*M_PI)/(L*q02) + 
    225          t2 + t3 - t4 + t5 - t6 + onehalf*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))* 
     222         t2 + t3 - t4 + t5 - t6 + 0.5*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))* 
    226223         (((-((b*M_PI)/(L*q0))) + t9 + t10 + 
    227          onehalf*((C3*pow(((Rgb)),((-3.0)/miu)) + 
     224         0.5*((C3*pow(((Rgb)),((-3.0)/miu)) + 
    228225         C2*pow(((Rgb)),((-2.0)/miu)) + 
    229226         C1*pow(((Rgb)),((-1.0)/miu))))* 
  • sasmodels/models/linear_pearls.c

    r2c74c11 r4962519  
    1919{ 
    2020    // Pearl volume 
    21     double pearl_vol = 4.0 /3.0 * M_PI * pow(radius, 3.0); 
     21    double pearl_vol = M_4PI_3 * cube(radius); 
    2222    // Return total volume 
    2323    return num_pearls * pearl_vol;; 
     
    3535    double contrast_pearl = pearl_sld - solvent_sld; 
    3636    //each volume 
    37     double pearl_vol = 4.0 /3.0 * M_PI * pow(radius, 3.0); 
     37    double pearl_vol = M_4PI_3 * cube(radius); 
    3838    //total volume 
    3939    double tot_vol = num_pearls * pearl_vol; 
     
    4343    double separation = edge_sep + 2.0 * radius; 
    4444 
    45     double x=q*radius; 
    46  
    47     // Try Taylor on x*xos(x) 
    48         // double out_cos = x - pow(x,3)/2 + pow(x,5)/24 - pow(x,7)/720 + pow(x,9)/40320; 
    49     // psi -= x*out_cos; 
    50  
    5145    //sine functions of a pearl 
    52     double psi = sin(q * radius); 
    53     psi -= x * cos(x); 
    54     psi /= pow((q * radius), 3.0); 
     46    double psi = sph_j1c(q * radius); 
    5547 
    5648    // N pearls contribution 
     
    6153    } 
    6254    // form factor for num_pearls 
    63     double form_factor = n_contrib; 
    64     form_factor *= pow((m_s*psi*3.0), 2.0); 
    65     form_factor /= (tot_vol * 1.0e4); 
     55    double form_factor = 1.0e-4 * n_contrib * square(m_s*psi) / tot_vol; 
    6656 
    6757    return form_factor; 
  • sasmodels/models/linear_pearls.py

    r40a87fa r4962519  
    6363single = False 
    6464 
    65 source = ["linear_pearls.c"] 
     65source = ["lib/sph_j1c.c", "linear_pearls.c"] 
    6666 
    6767demo = dict(scale=1.0, background=0.0, 
  • sasmodels/models/mass_fractal.c

    ra807206 r6d96b66  
    1 double form_volume(double radius); 
     1#define INVALID(p) (p.fractal_dim_mass < 1.0) 
    22 
    3 double Iq(double q, 
    4           double radius, 
    5           double fractal_dim_mass, 
    6           double cutoff_length); 
     3static double 
     4Iq(double q, double radius, double fractal_dim_mass, double cutoff_length) 
     5{ 
     6    //calculate P(q) 
     7    const double pq = square(sph_j1c(q*radius)); 
    78 
    8 static double _mass_fractal_kernel(double q, 
    9           double radius, 
    10           double fractal_dim_mass, 
    11           double cutoff_length) 
    12 { 
    13     // Actively check the argument. 
    14     if (fractal_dim_mass <= 1.0){ 
    15        return 0.0; 
     9    //calculate S(q) 
     10    // S(q) = gamma(D-1) sin((D-1)atan(q c))/q c^(D-1) (1+(q c)^2)^(-(D-1)/2) 
     11    // lim D->1 [gamma(D-1) sin((D-1) a)] = a 
     12    // lim q->0 [sin((D-1) atan(q c))/q] = (D-1) c 
     13    // lim q,D->0,1 [gamma(D-1) sin((D-1)atan(q c))/q] = c 
     14    double sq; 
     15    if (q > 0. && fractal_dim_mass > 1.) { 
     16        const double Dm1 = fractal_dim_mass - 1.0; 
     17        const double t1 = sas_gamma(Dm1)*sin(Dm1*atan(q*cutoff_length)); 
     18        const double t2 = pow(cutoff_length, Dm1); 
     19        const double t3 = pow(1.0 + square(q*cutoff_length), -0.5*Dm1); 
     20        sq = t1 * t2 * t3 / q; 
     21    } else if (q > 0.) { 
     22        sq = atan(q*cutoff_length)/q; 
     23    } else if (fractal_dim_mass > 1.) { 
     24        const double D = fractal_dim_mass; 
     25        sq = pow(cutoff_length, D) * sas_gamma(D); 
     26    } else { 
     27        sq = cutoff_length; 
    1628    } 
    1729 
    18     //calculate P(q) 
    19     double pq = sph_j1c(q*radius); 
    20     pq = pq*pq; 
    21  
    22     //calculate S(q) 
    23     double mmo = fractal_dim_mass-1.0; 
    24     double sq = sas_gamma(mmo)*sin((mmo)*atan(q*cutoff_length)); 
    25     sq *= pow(cutoff_length, mmo); 
    26     sq /= pow((1.0 + (q*cutoff_length)*(q*cutoff_length)),(mmo/2.0)); 
    27     sq /= q; 
    28  
    29     //combine and return 
    30     double result = pq * sq; 
    31  
    32     return result; 
     30    return pq * sq; 
    3331} 
    34 double form_volume(double radius){ 
    35  
    36     return 1.333333333333333*M_PI*radius*radius*radius; 
    37 } 
    38  
    39 double Iq(double q, 
    40           double radius, 
    41           double fractal_dim_mass, 
    42           double cutoff_length) 
    43 { 
    44     return _mass_fractal_kernel(q, 
    45            radius, 
    46            fractal_dim_mass, 
    47            cutoff_length); 
    48 } 
  • sasmodels/models/mass_fractal.py

    ra807206 r6d96b66  
    7878 
    7979# pylint: disable=bad-whitespace, line-too-long 
    80 #             ["name", "units", default, [lower, upper], "type","description"], 
    81 parameters = [["radius",        "Ang",  10.0, [0.0, inf], "", "Particle radius"], 
    82               ["fractal_dim_mass",      "",      1.9, [1.0, 6.0], "", "Mass fractal dimension"], 
    83               ["cutoff_length", "Ang", 100.0, [0.0, inf], "", "Cut-off length"], 
    84              ] 
     80#   ["name", "units", default, [lower, upper], "type","description"], 
     81parameters = [ 
     82    ["radius",           "Ang",  10.0, [0.0, inf], "", "Particle radius"], 
     83    ["fractal_dim_mass", "",      1.9, [1.0, 6.0], "", "Mass fractal dimension"], 
     84    ["cutoff_length",    "Ang", 100.0, [0.0, inf], "", "Cut-off length"], 
     85] 
    8586# pylint: enable=bad-whitespace, line-too-long 
    8687 
  • sasmodels/models/mass_surface_fractal.c

    r30fbe2e r6d96b66  
    1 double form_volume(double radius); 
    2  
    3 double Iq(double q, 
    4           double fractal_dim_mass, 
    5           double fractal_dim_surf, 
    6           double rg_cluster, 
    7           double rg_primary); 
    8  
    9 static double _mass_surface_fractal_kernel(double q, 
     1#define INVALID(p) (p.fractal_dim_mass + p.fractal_dim_surf > 6.) 
     2static double 
     3Iq(double q, 
    104          double fractal_dim_mass, 
    115          double fractal_dim_surf, 
     
    148{ 
    159     //computation 
    16     double tot_dim = 6.0 - fractal_dim_surf - fractal_dim_mass; 
    17     fractal_dim_mass /= 2.0; 
    18     tot_dim /= 2.0; 
     10    const double Dm = 0.5*fractal_dim_mass; 
     11    const double Dt = 0.5*(6.0 - (fractal_dim_mass + fractal_dim_surf)); 
    1912 
    20     double rc_norm = rg_cluster * rg_cluster / (3.0 * fractal_dim_mass); 
    21     double rp_norm = rg_primary * rg_primary / (3.0 * tot_dim); 
     13    const double t1 = Dm==0. ? 1.0 : pow(1.0 + square(q*rg_cluster)/(3.0*Dm), -Dm); 
     14    const double t2 = Dt==0. ? 1.0 : pow(1.0 + square(q*rg_primary)/(3.0*Dt), -Dt); 
     15    const double form = t1*t2; 
    2216 
    23     //x for P 
    24     double x_val1 = 1.0 +  q * q * rc_norm; 
    25     double x_val2 = 1.0 +  q * q * rp_norm; 
    26  
    27     double inv_form = pow(x_val1, fractal_dim_mass) * pow(x_val2, tot_dim); 
    28  
    29     //another singular 
    30     if (inv_form == 0.0) return 0.0; 
    31  
    32     double form_factor = 1.0; 
    33     form_factor /= inv_form; 
    34  
    35     return (form_factor); 
     17    return form; 
    3618} 
    37 double form_volume(double radius){ 
    38  
    39     return 1.333333333333333*M_PI*radius*radius*radius; 
    40 } 
    41  
    42 double Iq(double q, 
    43           double fractal_dim_mass, 
    44           double fractal_dim_surf, 
    45           double rg_cluster, 
    46           double rg_primary) 
    47 { 
    48     return _mass_surface_fractal_kernel(q, 
    49             fractal_dim_mass, 
    50             fractal_dim_surf, 
    51             rg_cluster, 
    52             rg_primary); 
    53 } 
  • sasmodels/models/mass_surface_fractal.py

    r30fbe2e r6d96b66  
    7878 
    7979# pylint: disable=bad-whitespace, line-too-long 
    80 #             ["name", "units", default, [lower, upper], "type","description"], 
    81 parameters = [["fractal_dim_mass",      "",    1.8, [1e-16, 6.0], "", 
    82                "Mass fractal dimension"], 
    83               ["fractal_dim_surf",   "",    2.3, [1e-16, 6.0], "", 
    84                "Surface fractal dimension"], 
    85               ["rg_cluster", "Ang",   86.7, [0.0, inf], "", 
    86                "Cluster radius of gyration"], 
    87               ["rg_primary", "Ang", 4000.,  [0.0, inf], "", 
    88                "Primary particle radius of gyration"], 
    89              ] 
     80#   ["name", "units", default, [lower, upper], "type","description"], 
     81parameters = [ 
     82    ["fractal_dim_mass", "",      1.8, [0.0, 6.0], "", "Mass fractal dimension"], 
     83    ["fractal_dim_surf", "",      2.3, [0.0, 6.0], "", "Surface fractal dimension"], 
     84    ["rg_cluster",       "Ang",  86.7, [0.0, inf], "", "Cluster radius of gyration"], 
     85    ["rg_primary",       "Ang", 4000., [0.0, inf], "", "Primary particle radius of gyration"], 
     86] 
    9087# pylint: enable=bad-whitespace, line-too-long 
    9188 
  • sasmodels/models/multilayer_vesicle.c

    r2c74c11 r3a48772  
    1919 
    2020        // layer 1 
    21         voli = 4.0*M_PI/3.0*ri*ri*ri; 
     21        voli = M_4PI_3*ri*ri*ri; 
    2222        fval += voli*sldi*sph_j1c(ri*q); 
    2323 
     
    2525 
    2626        // layer 2 
    27         voli = 4.0*M_PI/3.0*ri*ri*ri; 
     27        voli = M_4PI_3*ri*ri*ri; 
    2828        fval -= voli*sldi*sph_j1c(ri*q); 
    2929 
  • sasmodels/models/onion.c

    rd119f34 r9762341  
    3030 
    3131static double 
    32 form_volume(double core_radius, double n, double thickness[]) 
     32form_volume(double radius_core, double n, double thickness[]) 
    3333{ 
    3434  int i; 
    35   double r = core_radius; 
     35  double r = radius_core; 
    3636  for (i=0; i < n; i++) { 
    3737    r += thickness[i]; 
     
    4141 
    4242static double 
    43 Iq(double q, double sld_core, double core_radius, double sld_solvent, 
     43Iq(double q, double sld_core, double radius_core, double sld_solvent, 
    4444    double n_shells, double sld_in[], double sld_out[], double thickness[], 
    4545    double A[]) 
    4646{ 
    4747  int n = (int)(n_shells+0.5); 
    48   double r_out = core_radius; 
     48  double r_out = radius_core; 
    4949  double f = f_exp(q, r_out, sld_core, 0.0, 0.0, 0.0, 0.0); 
    5050  for (int i=0; i < n; i++){ 
  • sasmodels/models/onion.py

    r7b68dc5 r9762341  
    366366    return np.asarray(z), np.asarray(rho) 
    367367 
    368 def ER(core_radius, n, thickness): 
     368def ER(radius_core, n, thickness): 
    369369    """Effective radius""" 
    370     return np.sum(thickness[:int(n[0])], axis=0) + core_radius 
     370    return np.sum(thickness[:int(n[0])], axis=0) + radius_core 
    371371 
    372372demo = { 
    373373    "sld_solvent": 2.2, 
    374374    "sld_core": 1.0, 
    375     "core_radius": 100, 
     375    "radius_core": 100, 
    376376    "n_shells": 4, 
    377377    "sld_in": [0.5, 1.5, 0.9, 2.0], 
     
    381381    # Could also specify them individually as 
    382382    # "A1": 0, "A2": -1, "A3": 1e-4, "A4": 1, 
    383     #"core_radius_pd_n": 10, 
    384     #"core_radius_pd": 0.4, 
     383    #"radius_core_pd_n": 10, 
     384    #"radius_core_pd": 0.4, 
    385385    #"thickness4_pd_n": 10, 
    386386    #"thickness4_pd": 0.4, 
  • sasmodels/models/parallelepiped.c

    ra807206 r14838a3  
    11double form_volume(double length_a, double length_b, double length_c); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, double length_b, double length_c); 
     2double Iq(double q, double sld, double solvent_sld, 
     3    double length_a, double length_b, double length_c); 
    34double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    4     double length_a, double length_b, double length_c, double theta, double phi, double psi); 
    5  
    6 // From Igor library 
    7 double _pkernel(double a, double b,double c, double ala, double alb, double alc); 
    8 double _pkernel(double a, double b,double c, double ala, double alb, double alc){ 
    9     double argA,argB,argC,tmp1,tmp2,tmp3; 
    10     //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    11     argA = 0.5*a*ala; 
    12     argB = 0.5*b*alb; 
    13     argC = 0.5*c*alc; 
    14     if(argA==0.0) { 
    15         tmp1 = 1.0; 
    16     } else { 
    17         tmp1 = sin(argA)*sin(argA)/argA/argA; 
    18     } 
    19     if (argB==0.0) { 
    20         tmp2 = 1.0; 
    21     } else { 
    22         tmp2 = sin(argB)*sin(argB)/argB/argB; 
    23     } 
    24     if (argC==0.0) { 
    25         tmp3 = 1.0; 
    26     } else { 
    27         tmp3 = sin(argC)*sin(argC)/argC/argC; 
    28     } 
    29     return (tmp1*tmp2*tmp3); 
    30  
    31 } 
    32  
     5    double length_a, double length_b, double length_c, 
     6    double theta, double phi, double psi); 
    337 
    348double form_volume(double length_a, double length_b, double length_c) 
     
    4519    double length_c) 
    4620{ 
    47     double tmp1, tmp2; 
    48      
    49     double mu = q * length_b; 
     21    const double mu = 0.5 * q * length_b; 
    5022     
    5123    // Scale sides by B 
    52     double a_scaled = length_a / length_b; 
    53     double c_scaled = length_c / length_b; 
     24    const double a_scaled = length_a / length_b; 
     25    const double c_scaled = length_c / length_b; 
    5426         
    55     //Order of integration 
    56     int nordi=76;                                
    57     int nordj=76; 
     27    // outer integral (with gauss points), integration limits = 0, 1 
     28    double outer_total = 0; //initialize integral 
    5829 
    59     // outer integral (with nordi gauss points), integration limits = 0, 1 
    60     double summ = 0; //initialize integral 
     30    for( int i=0; i<76; i++) { 
     31        const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     32        const double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    6133 
    62     for( int i=0; i<nordi; i++) { 
    63                  
    64         // inner integral (with nordj gauss points), integration limits = 0, 1 
    65          
    66         double summj = 0.0; 
    67             double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );          
    68                  
    69             for(int j=0; j<nordj; j++) { 
     34        // inner integral (with gauss points), integration limits = 0, 1 
     35        // corresponding to angles from 0 to pi/2. 
     36        double inner_total = 0.0; 
     37        for(int j=0; j<76; j++) { 
     38            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     39            double sin_uu, cos_uu; 
     40            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
     41            const double si1 = sinc(mu_proj * sin_uu * a_scaled); 
     42            const double si2 = sinc(mu_proj * cos_uu); 
     43            inner_total += Gauss76Wt[j] * square(si1 * si2); 
     44        } 
     45        inner_total *= 0.5; 
    7046 
    71             double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    72             double mudum = mu * sqrt(1.0-sigma*sigma); 
    73                 double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu); 
    74                 double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu); 
    75             if(arg1==0.0) { 
    76                 tmp1 = 1.0; 
    77             } else { 
    78                 tmp1 = sin(arg1)*sin(arg1)/arg1/arg1; 
    79             } 
    80             if (arg2==0.0) { 
    81                 tmp2 = 1.0; 
    82             } else { 
    83                 tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; 
    84             } 
     47        const double si = sinc(mu * c_scaled * sigma); 
     48        outer_total += Gauss76Wt[i] * inner_total * si * si; 
     49    } 
     50    outer_total *= 0.5; 
    8551 
    86             summj += Gauss76Wt[j] * tmp1 * tmp2; 
    87         } 
    88                  
    89         // value of the inner integral 
    90         double answer = 0.5 * summj; 
    91  
    92         double arg = 0.5 * mu * c_scaled * sigma; 
    93         if ( arg == 0.0 ) { 
    94             answer *= 1.0; 
    95         } else { 
    96             answer *= sin(arg)*sin(arg)/arg/arg; 
    97         } 
    98                  
    99             // sum of outer integral 
    100         summ += Gauss76Wt[i] * answer; 
    101          
    102     }    
    103     
    104     const double vd = (sld-solvent_sld) * form_volume(length_a, length_b, length_c); 
    105      
    106     // convert from [1e-12 A-1] to [cm-1] and 0.5 factor for outer integral 
    107     return 1.0e-4 * 0.5 * vd * vd * summ; 
    108      
     52    // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1] 
     53    const double V = form_volume(length_a, length_b, length_c); 
     54    const double drho = (sld-solvent_sld); 
     55    return 1.0e-4 * square(drho * V) * outer_total; 
    10956} 
    11057 
     
    12067    double psi) 
    12168{ 
    122     double q = sqrt(qx*qx+qy*qy); 
    123     double qx_scaled = qx/q; 
    124     double qy_scaled = qy/q; 
     69    double q, cos_val_a, cos_val_b, cos_val_c; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 
    12571 
    126     // Convert angles given in degrees to radians 
    127     theta *= M_PI_180; 
    128     phi   *= M_PI_180; 
    129     psi   *= M_PI_180; 
    130      
    131     // Parallelepiped c axis orientation 
    132     double cparallel_x = cos(theta) * cos(phi); 
    133     double cparallel_y = sin(theta); 
    134      
    135     // Compute angle between q and parallelepiped axis 
    136     double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz; 
    137  
    138     // Parallelepiped a axis orientation 
    139     double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 
    140     double parallel_y = sin(psi)*cos(theta); 
    141     double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled; 
    142  
    143     // Parallelepiped b axis orientation 
    144     double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 
    145     double bparallel_y = cos(theta)*cos(psi); 
    146     double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled; 
    147  
    148     // The following tests should always pass 
    149     if (fabs(cos_val_c)>1.0) { 
    150       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    151       cos_val_c = 1.0; 
    152     } 
    153     if (fabs(cos_val_a)>1.0) { 
    154       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    155       cos_val_a = 1.0; 
    156     } 
    157     if (fabs(cos_val_b)>1.0) { 
    158       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    159       cos_val_b = 1.0; 
    160     } 
    161      
    162     // Call the IGOR library function to get the kernel 
    163     double form = _pkernel( q*length_a, q*length_b, q*length_c, cos_val_a, cos_val_b, cos_val_c); 
    164    
    165     // Multiply by contrast^2 
    166     const double vd = (sld - solvent_sld) * form_volume(length_a, length_b, length_c); 
    167     return 1.0e-4 * vd * vd * form; 
     72    const double siA = sinc(0.5*q*length_a*cos_val_a); 
     73    const double siB = sinc(0.5*q*length_b*cos_val_b); 
     74    const double siC = sinc(0.5*q*length_c*cos_val_c); 
     75    const double V = form_volume(length_a, length_b, length_c); 
     76    const double drho = (sld - solvent_sld); 
     77    const double form = V * drho * siA * siB * siC; 
     78    // Square and convert from [1e-12 A-1] to [cm-1] 
     79    return 1.0e-4 * form * form; 
    16880} 
  • sasmodels/models/parallelepiped.py

    r416f5c7 red0827a  
    221221# parameters for demo 
    222222demo = dict(scale=1, background=0, 
    223             sld=6.3e-6, sld_solvent=1.0e-6, 
     223            sld=6.3, sld_solvent=1.0, 
    224224            length_a=35, length_b=75, length_c=400, 
    225225            theta=45, phi=30, psi=15, 
  • sasmodels/models/pearl_necklace.c

    ra807206 r2126131  
    1 double _pearl_necklace_kernel(double q, double radius, double edge_sep, 
    2         double thick_string, double num_pearls, double sld_pearl, 
    3         double sld_string, double sld_solv); 
    41double form_volume(double radius, double edge_sep, 
    5         double thick_string, double num_pearls); 
    6  
     2    double thick_string, double num_pearls); 
    73double Iq(double q, double radius, double edge_sep, 
    8         double thick_string, double num_pearls, double sld,  
    9         double string_sld, double solvent_sld); 
     4    double thick_string, double num_pearls, double sld,  
     5    double string_sld, double solvent_sld); 
    106 
    117#define INVALID(v) (v.thick_string >= v.radius || v.num_pearls <= 0) 
    128 
    139// From Igor library 
    14 double _pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string, 
    15         double num_pearls, double sld_pearl, double sld_string, double sld_solv) 
     10static double 
     11_pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string, 
     12    double num_pearls, double sld_pearl, double sld_string, double sld_solv) 
    1613{ 
    17         //relative slds 
    18         double contrast_pearl = sld_pearl - sld_solv; 
    19         double contrast_string = sld_string - sld_solv; 
    20          
    21         // number of string segments 
    22         num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
    23         double num_strings = num_pearls - 1.0; 
    24          
    25         //Pi 
    26         double pi = 4.0*atan(1.0); 
    27          
    28         // center to center distance between the neighboring pearls 
    29         double A_s = edge_sep + 2.0 * radius; 
    30          
    31         // Repeated Calculations 
    32         double sincasq = sinc(q*A_s); 
    33         double oneminussinc = 1 - sincasq; 
    34         double q_r = q * radius; 
    35         double q_edge = q * edge_sep; 
    36          
    37         // each volume 
    38         double string_vol = edge_sep * pi * thick_string * thick_string / 4.0; 
    39         double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius; 
     14    // number of string segments 
     15    num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
     16    const double num_strings = num_pearls - 1.0; 
    4017 
    41         //total volume 
    42         double tot_vol; 
    43         //each masses 
    44         double m_r= contrast_string * string_vol; 
    45         double m_s= contrast_pearl * pearl_vol; 
    46         double psi, gamma, beta; 
    47         //form factors 
    48         double sss, srr, srs; //cross 
    49         double srr_1, srr_2, srr_3; 
    50         double form_factor; 
    51         tot_vol = num_strings * string_vol; 
    52         tot_vol += num_pearls * pearl_vol; 
     18    //each masses: contrast * volume 
     19    const double contrast_pearl = sld_pearl - sld_solv; 
     20    const double contrast_string = sld_string - sld_solv; 
     21    const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; 
     22    const double pearl_vol = M_4PI_3 * radius * radius * radius; 
     23    const double m_string = contrast_string * string_vol; 
     24    const double m_pearl = contrast_pearl * pearl_vol; 
    5325 
    54         //sine functions of a pearl 
    55         psi = sin(q_r); 
    56         psi -= q_r * cos(q_r); 
    57         psi *= 3.0; 
    58         psi /= q_r * q_r * q_r; 
     26    // center to center distance between the neighboring pearls 
     27    const double A_s = edge_sep + 2.0 * radius; 
    5928 
    60         // Note take only 20 terms in Si series: 10 terms may be enough though. 
    61         gamma = Si(q_edge); 
    62         gamma /= (q_edge); 
    63         beta = Si(q * (A_s - radius)); 
    64         beta -= Si(q_r); 
    65         beta /= q_edge; 
     29    //sine functions of a pearl 
     30    // Note: lim_(q->0) Si(q*a)/(q*b) = a/b 
     31    // So therefore: 
     32    //    beta = q==0. ? 1.0 : (Si(q*(A_s-radius)) - Si(q*radius))/q_edge; 
     33    //    gamma = q==0. ? 1.0 : Si(q_edge)/q_edge; 
     34    // But there is a 1/(1-sinc) term below which blows up so don't bother 
     35    const double q_edge = q * edge_sep; 
     36    const double beta = (Si(q*(A_s-radius)) - Si(q*radius)) / q_edge; 
     37    const double gamma = Si(q_edge) / q_edge; 
     38    const double psi = sph_j1c(q*radius); 
    6639 
    67         // form factor for num_pearls 
    68         sss = 1.0 - pow(sincasq, num_pearls); 
    69         sss /= oneminussinc * oneminussinc; 
    70         sss *= -sincasq; 
    71         sss -= num_pearls / 2.0; 
    72         sss += num_pearls / oneminussinc; 
    73         sss *= 2.0 * m_s * psi * m_s * psi; 
     40    // Precomputed sinc terms 
     41    const double si = sinc(q*A_s); 
     42    const double omsi = 1.0 - si; 
     43    const double pow_si = pow(si, num_pearls); 
    7444 
    75         // form factor for num_strings (like thin rods) 
    76         srr_1 = -sinc(q_edge/2.0) * sinc(q_edge/2.0); 
     45    // form factor for num_pearls 
     46    const double sss = 2.0*square(m_pearl*psi) * ( 
     47        - si * (1.0 - pow_si) / (omsi*omsi) 
     48        + num_pearls / omsi 
     49        - 0.5 * num_pearls 
     50        ); 
    7751 
    78         srr_1 += 2.0 * gamma; 
    79         srr_1 *= num_strings; 
    80         srr_2 = 2.0/oneminussinc; 
    81         srr_2 *= num_strings; 
    82         srr_2 *= beta * beta; 
    83         srr_3 = 1.0 - pow(sincasq, num_strings); 
    84         srr_3 /= oneminussinc * oneminussinc; 
    85         srr_3 *= beta * beta; 
    86         srr_3 *= -2.0; 
     52    // form factor for num_strings (like thin rods) 
     53    const double srr = m_string * m_string * ( 
     54        - 2.0 * (1.0 - pow_si/si)*beta*beta / (omsi*omsi) 
     55        + 2.0 * num_strings*beta*beta / omsi 
     56        + num_strings * (2.0*gamma - square(sinc(q_edge/2.0))) 
     57        ); 
    8758 
    88         // total srr 
    89         srr = srr_1 + srr_2 + srr_3; 
    90         srr *= m_r * m_r; 
     59    // form factor for correlations 
     60    const double srs = 4.0 * m_string * m_pearl * beta * psi * ( 
     61        - si * (1.0 - pow_si/si) / (omsi*omsi) 
     62        + num_strings / omsi 
     63        ); 
    9164 
    92         // form factor for correlations 
    93         srs = 1.0; 
    94         srs -= pow(sincasq, num_strings); 
    95         srs /= oneminussinc * oneminussinc; 
    96         srs *= -sincasq; 
    97         srs += num_strings/oneminussinc; 
    98         srs *= 4.0; 
    99         srs *= (m_r * m_s * beta * psi); 
     65    const double form = sss + srr + srs; 
    10066 
    101         form_factor = sss + srr + srs; 
    102         form_factor /= (tot_vol * 1.0e4); // norm by volume and A^-1 to cm^-1 
    103  
    104         return (form_factor); 
     67    return 1.0e-4 * form; 
    10568} 
    10669 
    10770double form_volume(double radius, double edge_sep, 
    108         double thick_string, double num_pearls) 
     71    double thick_string, double num_pearls) 
    10972{ 
    110         double total_vol; 
     73    num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
    11174 
    112         double pi = 4.0*atan(1.0); 
    113         double number_of_strings = num_pearls - 1.0; 
    114          
    115         double string_vol = edge_sep * pi * thick_string * thick_string / 4.0; 
    116         double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius; 
     75    const double num_strings = num_pearls - 1.0; 
     76    const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; 
     77    const double pearl_vol = M_4PI_3 * radius * radius * radius; 
     78    const double volume = num_strings*string_vol + num_pearls*pearl_vol; 
    11779 
    118         total_vol = number_of_strings * string_vol; 
    119         total_vol += num_pearls * pearl_vol; 
    120  
    121         return(total_vol); 
     80    return volume; 
    12281} 
    12382 
    12483double Iq(double q, double radius, double edge_sep, 
    125         double thick_string, double num_pearls, double sld,  
    126         double string_sld, double solvent_sld) 
     84    double thick_string, double num_pearls, double sld,  
     85    double string_sld, double solvent_sld) 
    12786{ 
    128         double value, tot_vol; 
    129          
    130         value = _pearl_necklace_kernel(q, radius, edge_sep, thick_string, 
    131                 num_pearls, sld, string_sld, solvent_sld); 
    132         tot_vol = form_volume(radius, edge_sep, thick_string, num_pearls); 
     87    const double form = _pearl_necklace_kernel(q, radius, edge_sep, 
     88        thick_string, num_pearls, sld, string_sld, solvent_sld); 
    13389 
    134         return value*tot_vol; 
     90    return form; 
    13591} 
  • sasmodels/models/pearl_necklace.py

    ra807206 r2126131  
    9292             ] 
    9393 
    94 source = ["lib/Si.c", "pearl_necklace.c"] 
     94source = ["lib/Si.c", "lib/sph_j1c.c", "pearl_necklace.c"] 
    9595single = False  # use double precision unless told otherwise 
    9696 
     
    112112    """ 
    113113    tot_vol = volume(radius, edge_sep, thick_string, num_pearls) 
    114     rad_out = pow((3.0*tot_vol/4.0/pi), 0.33333) 
     114    rad_out = (tot_vol/(4.0/3.0*pi)) ** (1./3.) 
    115115    return rad_out 
    116116 
  • sasmodels/models/poly_gauss_coil.py

    ra807206 rf0afad2  
    5959""" 
    6060 
    61 from numpy import inf, exp, power 
     61import numpy as np 
     62from numpy import inf, expm1, power 
    6263 
    6364name = "poly_gauss_coil" 
     
    8384    # pylint: disable = missing-docstring 
    8485    u = polydispersity - 1.0 
    85     z = (q*rg)**2 / (1.0 + 2.0*u) 
     86    z = q**2 * (rg**2 / (1.0 + 2.0*u)) 
     87 
    8688    # need to trap the case of the polydispersity being 1 (ie, monodisperse!) 
    8789    if polydispersity == 1.0: 
    88         inten = i_zero * 2.0 * (exp(-z) + z - 1.0) 
     90        result = 2.0 * (expm1(-z) + z) 
     91        index = q != 0. 
     92        result[index] /= z[index]**2 
     93        result[~index] = 1.0 
    8994    else: 
    90         inten = i_zero * 2.0 * (power(1.0 + u*z, -1.0/u) + z - 1.0) / (1.0 + u) 
    91     index = q != 0. 
    92     inten[~index] = i_zero 
    93     inten[index] /= z[index]**2 
    94     return inten 
     95        # Taylor series around z=0 of (2*(1+uz)^(-1/u) + z - 1) / (z^2(u+1)) 
     96        p = [ 
     97            #(-1 - 20*u - 155*u**2 - 580*u**3 - 1044*u**4 - 720*u**5) / 2520., 
     98            #( 1 + 14*u + 71*u**2 + 154*u**3 + 120*u**4) / 360., 
     99            #(-1 - 9*u - 26*u**2 - 24*u**3) / 60., 
     100            ( 1 + 5*u + 6*u**2) / 12., 
     101            (-1 - 2*u) / 3., 
     102            ( 1 ), 
     103            ] 
     104        result = 2.0 * (power(1.0 + u*z, -1.0/u) + z - 1.0) / (1.0 + u) 
     105        index = z > 1e-4 
     106        result[index] /= z[index]**2 
     107        result[~index] = np.polyval(p, z[~index]) 
     108    return i_zero * result 
    95109Iq.vectorized = True  # Iq accepts an array of q values 
    96110 
  • sasmodels/models/polymer_micelle.c

    ra807206 rc3ebc71  
    3232    // Self-correlation term of the core 
    3333    const double bes_core = sph_j1c(q*radius_core); 
    34     const double term1 = n_aggreg*n_aggreg*beta_core*beta_core*bes_core*bes_core; 
     34    const double term1 = square(n_aggreg*beta_core*bes_core); 
    3535 
    3636    // Self-correlation term of the chains 
    37     const double qrg2 = q*rg*q*rg; 
     37    const double qrg2 = square(q*rg); 
    3838    const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(expm1(-qrg2)+qrg2)/(qrg2*qrg2); 
    3939    const double term2 = n_aggreg * beta_corona * beta_corona * debye_chain; 
     
    4242    const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2; 
    4343    const double bes_corona = sinc(q*(radius_core + d_penetration * rg)); 
    44     const double term3 = 2 * n_aggreg * n_aggreg * beta_core * beta_corona * 
     44    const double term3 = 2.0 * n_aggreg * n_aggreg * beta_core * beta_corona * 
    4545                 bes_core * chain_ampl * bes_corona; 
    4646 
    4747    // Interference cross-term between chains 
    48     const double term4 = n_aggreg * (n_aggreg - 1.0) * beta_corona * beta_corona * 
    49                  chain_ampl * chain_ampl * bes_corona * bes_corona; 
     48    const double term4 = n_aggreg * (n_aggreg - 1.0) 
     49                 * square(beta_corona * chain_ampl * bes_corona); 
    5050 
    5151    // I(q)_micelle : Sum of 4 terms computed above 
  • sasmodels/models/porod.py

    r40a87fa r4962519  
    4141    """ 
    4242    with errstate(divide='ignore'): 
    43         return power(q, -4) 
     43        return q**-4 
    4444 
    4545Iq.vectorized = True  # Iq accepts an array of q values 
  • sasmodels/models/rectangular_prism.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) 
     
    1715    double c2a_ratio) 
    1816{ 
    19     double termA, termB, termC; 
    20      
    21     double b_side = length_a * b2a_ratio; 
    22     double c_side = length_a * c2a_ratio; 
    23     double volume = length_a * b_side * c_side; 
    24     double a_half = 0.5 * length_a; 
    25     double b_half = 0.5 * b_side; 
    26     double c_half = 0.5 * c_side; 
     17    const double length_b = length_a * b2a_ratio; 
     18    const double length_c = length_a * c2a_ratio; 
     19    const double a_half = 0.5 * length_a; 
     20    const double b_half = 0.5 * length_b; 
     21    const double c_half = 0.5 * length_c; 
    2722 
    2823   //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 
     24    const double v1a = 0.0; 
     25    const double v1b = M_PI_2;  //theta integration limits 
     26    const double v2a = 0.0; 
     27    const double v2b = M_PI_2;  //phi integration limits 
    3328     
    34     //Order of integration 
    35     int nordi=76;                                
    36     int nordj=76; 
     29    double outer_sum = 0.0; 
     30    for(int i=0; i<76; i++) { 
     31        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     32        double sin_theta, cos_theta; 
     33        SINCOS(theta, sin_theta, cos_theta); 
    3734 
    38     double sumi = 0.0; 
    39      
    40     for(int i=0; i<nordi; i++) { 
     35        const double termC = sinc(q * c_half * cos_theta); 
    4136 
    42             double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );  
     37        double inner_sum = 0.0; 
     38        for(int j=0; j<76; j++) { 
     39            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     40            double sin_phi, cos_phi; 
     41            SINCOS(phi, sin_phi, cos_phi); 
    4342 
    44             double arg = q * c_half * cos(theta); 
    45             if (fabs(arg) > 1.e-16) {termC = sin(arg)/arg;} else {termC = 1.0;}   
    46  
    47             double sumj = 0.0; 
    48          
    49             for(int j=0; j<nordj; j++) { 
    50  
    51             double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );  
    52  
    53                 // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 
    54  
    55                 arg = q * a_half * sin(theta) * sin(phi);  
    56                 if (fabs(arg) > 1.e-16) {termA = sin(arg)/arg;} else {termA = 1.0;} 
    57                 
    58                 arg = q * b_half * sin(theta) * cos(phi);  
    59                 if (fabs(arg) > 1.e-16) {termB = sin(arg)/arg;} else {termB = 1.0;}        
    60                 
    61                 double AP = termA * termB * termC;   
    62  
    63                 sumj += Gauss76Wt[j] * (AP*AP); 
    64  
    65             } 
    66  
    67