Changeset b3f4831 in sasmodels


Ignore:
Timestamp:
Oct 30, 2018 11:07:41 AM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
ticket_1156
Children:
cc8b183
Parents:
5778279 (diff), c6084f1 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket_1156

Files:
1 added
22 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/magnetism/magnetism.rst

    rbefe905 rdf87acf  
    8989 
    9090===========   ================================================================ 
    91  M0:sld       $D_M M_0$ 
    92  mtheta:sld   $\theta_M$ 
    93  mphi:sld     $\phi_M$ 
    94  up:angle     $\theta_\mathrm{up}$ 
    95  up:frac_i    $u_i$ = (spin up)/(spin up + spin down) *before* the sample 
    96  up:frac_f    $u_f$ = (spin up)/(spin up + spin down) *after* the sample 
     91 sld_M0       $D_M M_0$ 
     92 sld_mtheta   $\theta_M$ 
     93 sld_mphi     $\phi_M$ 
     94 up_frac_i    $u_i$ = (spin up)/(spin up + spin down) *before* the sample 
     95 up_frac_f    $u_f$ = (spin up)/(spin up + spin down) *after* the sample 
     96 up_angle     $\theta_\mathrm{up}$ 
    9797===========   ================================================================ 
    9898 
    9999.. note:: 
    100     The values of the 'up:frac_i' and 'up:frac_f' must be in the range 0 to 1. 
     100    The values of the 'up_frac_i' and 'up_frac_f' must be in the range 0 to 1. 
    101101 
    102102*Document History* 

  • doc/guide/plugin.rst

    rf796469 r57c609b  
    423423calculations, but instead rely on numerical integration to compute the 
    424424appropriately smeared pattern. 
     425 
     426Each .py file also contains a function:: 
     427 
     428        def random(): 
     429        ... 
     430 
     431This function provides a model-specific random parameter set which shows model 
     432features in the USANS to SANS range.  For example, core-shell sphere sets the 
     433outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q 
     434value for the transition from flat to falling.  It then uses a beta distribution 
     435to set the percentage of the shape which is shell, giving a preference for very 
     436thin or very thick shells (but never 0% or 100%).  Using `-sets=10` in sascomp 
     437should show a reasonable variety of curves over the default sascomp q range. 
     438The parameter set is returned as a dictionary of `{parameter: value, ...}`. 
     439Any model parameters not included in the dictionary will default according to 
     440the code in the `_randomize_one()` function from sasmodels/compare.py. 
    425441 
    426442Python Models 
    … …  
    685701    erf, erfc, tgamma, lgamma:  **do not use** 
    686702        Special functions that should be part of the standard, but are missing 
    687         or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma 
    688         instead (see below). Note: lgamma(x) has not yet been tested. 
     703        or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma 
     704        and sas_lgamma instead (see below). 
    689705 
    690706Some non-standard constants and functions are also provided: 
    … …  
    753769        Gamma function sas_gamma\ $(x) = \Gamma(x)$. 
    754770 
    755         The standard math function, tgamma(x) is unstable for $x < 1$ 
     771        The standard math function, tgamma(x), is unstable for $x < 1$ 
    756772        on some platforms. 
    757773 
    758774        :code:`source = ["lib/sas_gamma.c", ...]` 
    759775        (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_) 
     776 
     777    sas_gammaln(x): 
     778        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. 
     779 
     780        The standard math function, lgamma(x), is incorrect for single 
     781        precision on some platforms. 
     782 
     783        :code:`source = ["lib/sas_gammainc.c", ...]` 
     784        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
     785 
     786    sas_gammainc(a, x), sas_gammaincc(a, x): 
     787        Incomplete gamma function 
     788        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     789        and complementary incomplete gamma function 
     790        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     791 
     792        :code:`source = ["lib/sas_gammainc.c", ...]` 
     793        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
    760794 
    761795    sas_erf(x), sas_erfc(x): 
    … …  
    795829        If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively. 
    796830 
     831        Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100]. 
     832 
    797833        The standard math function jn(n, x) is not available on all platforms. 
    798834 
    … …  
    803839        Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$. 
    804840 
     841        Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100]. 
     842 
    805843        This function uses Taylor series for small and large arguments: 
    806844 
    807         For large arguments, 
     845        For large arguments use the following Taylor series, 
    808846 
    809847        .. math:: 
    … …  
    813851             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 
    814852 
    815         For small arguments, 
     853        For small arguments , 
    816854 
    817855        .. math:: 

  • explore/precision.py

    r2a7e20e rfba9ca0  
    9595            neg:    [-100,100] 
    9696 
     97        For arbitrary range use "start:stop:steps:scale" where scale is 
     98        one of log, lin, or linear. 
     99 
    97100        *diff* is "relative", "absolute" or "none" 
    98101 
    … …  
    102105        linear = not xrange.startswith("log") 
    103106        if xrange == "zoom": 
    104             lin_min, lin_max, lin_steps = 1000, 1010, 2000 
     107            start, stop, steps = 1000, 1010, 2000 
    105108        elif xrange == "neg": 
    106             lin_min, lin_max, lin_steps = -100.1, 100.1, 2000 
     109            start, stop, steps = -100.1, 100.1, 2000 
    107110        elif xrange == "linear": 
    108             lin_min, lin_max, lin_steps = 1, 1000, 2000 
    109             lin_min, lin_max, lin_steps = 0.001, 2, 2000 
     111            start, stop, steps = 1, 1000, 2000 
     112            start, stop, steps = 0.001, 2, 2000 
    110113        elif xrange == "log": 
    111             log_min, log_max, log_steps = -3, 5, 400 
     114            start, stop, steps = -3, 5, 400 
    112115        elif xrange == "logq": 
    113             log_min, log_max, log_steps = -4, 1, 400 
     116            start, stop, steps = -4, 1, 400 
     117        elif ':' in xrange: 
     118            parts = xrange.split(':') 
     119            linear = parts[3] != "log" if len(parts) == 4 else True 
     120            steps = int(parts[2]) if len(parts) > 2 else 400 
     121            start = float(parts[0]) 
     122            stop = float(parts[1]) 
     123 
    114124        else: 
    115125            raise ValueError("unknown range "+xrange) 
    … …  
    121131            # value to x in the given precision. 
    122132            if linear: 
    123                 lin_min = max(lin_min, self.limits[0]) 
    124                 lin_max = min(lin_max, self.limits[1]) 
    125                 qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single') 
    126                 #qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='double') 
     133                start = max(start, self.limits[0]) 
     134                stop = min(stop, self.limits[1]) 
     135                qrf = np.linspace(start, stop, steps, dtype='single') 
     136                #qrf = np.linspace(start, stop, steps, dtype='double') 
    127137                qr = [mp.mpf(float(v)) for v in qrf] 
    128                 #qr = mp.linspace(lin_min, lin_max, lin_steps) 
     138                #qr = mp.linspace(start, stop, steps) 
    129139            else: 
    130                 log_min = np.log10(max(10**log_min, self.limits[0])) 
    131                 log_max = np.log10(min(10**log_max, self.limits[1])) 
    132                 qrf = np.logspace(log_min, log_max, log_steps, dtype='single') 
    133                 #qrf = np.logspace(log_min, log_max, log_steps, dtype='double') 
     140                start = np.log10(max(10**start, self.limits[0])) 
     141                stop = np.log10(min(10**stop, self.limits[1])) 
     142                qrf = np.logspace(start, stop, steps, dtype='single') 
     143                #qrf = np.logspace(start, stop, steps, dtype='double') 
    134144                qr = [mp.mpf(float(v)) for v in qrf] 
    135                 #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)] 
     145                #qr = [10**v for v in mp.linspace(start, stop, steps)] 
    136146 
    137147        target = self.call_mpmath(qr, bits=500) 
    … …  
    176186    """ 
    177187    if diff == "relative": 
    178         err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd') 
     188        err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd') 
    179189        #err = np.clip(err, 0, 1) 
    180190        pylab.loglog(x, err, '-', label=label) 
    … …  
    197207    return model_info 
    198208 
     209# Hack to allow second parameter A in two parameter functions 
     210A = 1 
     211def parse_extra_pars(): 
     212    global A 
     213 
     214    A_str = str(A) 
     215    pop = [] 
     216    for k, v in enumerate(sys.argv[1:]): 
     217        if v.startswith("A="): 
     218            A_str = v[2:] 
     219            pop.append(k+1) 
     220    if pop: 
     221        sys.argv = [v for k, v in enumerate(sys.argv) if k not in pop] 
     222        A = float(A_str) 
     223 
     224parse_extra_pars() 
     225 
    199226 
    200227# =============== FUNCTION DEFINITIONS ================ 
    … …  
    297324    ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]), 
    298325    limits=(-3.1, 10), 
     326) 
     327add_function( 
     328    name="gammaln(x)", 
     329    mp_function=mp.loggamma, 
     330    np_function=scipy.special.gammaln, 
     331    ocl_function=make_ocl("return sas_gammaln(q);", "sas_gammaln", ["lib/sas_gammainc.c"]), 
     332    #ocl_function=make_ocl("return lgamma(q);", "sas_gammaln"), 
     333) 
     334add_function( 
     335    name="gammainc(x)", 
     336    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a), 
     337    np_function=lambda x, a=A: scipy.special.gammainc(a, x), 
     338    ocl_function=make_ocl("return sas_gammainc(%.15g,q);"%A, "sas_gammainc", ["lib/sas_gammainc.c"]), 
     339) 
     340add_function( 
     341    name="gammaincc(x)", 
     342    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a), 
     343    np_function=lambda x, a=A: scipy.special.gammaincc(a, x), 
     344    ocl_function=make_ocl("return sas_gammaincc(%.15g,q);"%A, "sas_gammaincc", ["lib/sas_gammainc.c"]), 
    299345) 
    300346add_function( 
    … …  
    463509lanczos_gamma = """\ 
    464510    const double coeff[] = { 
    465             76.18009172947146,     -86.50532032941677, 
    466             24.01409824083091,     -1.231739572450155, 
     511            76.18009172947146, -86.50532032941677, 
     512            24.01409824083091, -1.231739572450155, 
    467513            0.1208650973866179e-2,-0.5395239384953e-5 
    468514            }; 
    … …  
    475521""" 
    476522add_function( 
    477     name="log gamma(x)", 
     523    name="loggamma(x)", 
    478524    mp_function=mp.loggamma, 
    479525    np_function=scipy.special.gammaln, 
    … …  
    599645 
    600646ALL_FUNCTIONS = set(FUNCTIONS.keys()) 
    601 ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet 
     647ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead 
    602648ALL_FUNCTIONS.discard("3j1/x:taylor") 
    603649ALL_FUNCTIONS.discard("3j1/x:trig") 
    … …  
    615661    -r indicates that the relative error should be plotted (default), 
    616662    -x<range> indicates the steps in x, where <range> is one of the following 
    617       log indicates log stepping in [10^-3, 10^5] (default) 
    618       logq indicates log stepping in [10^-4, 10^1] 
    619       linear indicates linear stepping in [1, 1000] 
    620       zoom indicates linear stepping in [1000, 1010] 
    621       neg indicates linear stepping in [-100.1, 100.1] 
    622 and name is "all" or one of: 
     663        log indicates log stepping in [10^-3, 10^5] (default) 
     664        logq indicates log stepping in [10^-4, 10^1] 
     665        linear indicates linear stepping in [1, 1000] 
     666        zoom indicates linear stepping in [1000, 1010] 
     667        neg indicates linear stepping in [-100.1, 100.1] 
     668        start:stop:n[:stepping] indicates an n-step plot in [start, stop] 
     669            or [10^start, 10^stop] if stepping is "log" (default n=400) 
     670Some functions (notably gammainc/gammaincc) have an additional parameter A 
     671which can be set from the command line as A=value.  Default is A=1. 
     672 
     673Name is one of: 
    623674    """+names) 
    624675    sys.exit(1) 

  • sasmodels/__init__.py

    re65c3ba ra1ec908  
    1414defining new models. 
    1515""" 
    16 __version__ = "0.97" 
     16__version__ = "0.98" 
    1717 
    1818def data_files(): 

  • sasmodels/compare.py

    rbd7630d r610ef23  
    368368 
    369369    # Limit magnetic SLDs to a smaller range, from zero to iron=5/A^2 
    370     if par.name.startswith('M0:'): 
     370    if par.name.endswith('_M0'): 
    371371        return np.random.uniform(0, 5) 
    372372 
    … …  
    538538    magnetic_pars = [] 
    539539    for p in parameters.user_parameters(pars, is2d): 
    540         if any(p.id.startswith(x) for x in ('M0:', 'mtheta:', 'mphi:')): 
     540        if any(p.id.endswith(x) for x in ('_M0', '_mtheta', '_mphi')): 
    541541            continue 
    542542        if p.id.startswith('up:'): 
    … …  
    550550            pdtype=pars.get(p.id+"_pd_type", 'gaussian'), 
    551551            relative_pd=p.relative_pd, 
    552             M0=pars.get('M0:'+p.id, 0.), 
    553             mphi=pars.get('mphi:'+p.id, 0.), 
    554             mtheta=pars.get('mtheta:'+p.id, 0.), 
     552            M0=pars.get(p.id+'_M0', 0.), 
     553            mphi=pars.get(p.id+'_mphi', 0.), 
     554            mtheta=pars.get(p.id+'_mtheta', 0.), 
    555555        ) 
    556556        lines.append(_format_par(p.name, **fields)) 
    … …  
    618618    if suppress: 
    619619        for p in pars: 
    620             if p.startswith("M0:"): 
     620            if p.endswith("_M0"): 
    621621                pars[p] = 0 
    622622    else: 
    … …  
    624624        first_mag = None 
    625625        for p in pars: 
    626             if p.startswith("M0:"): 
     626            if p.endswith("_M0"): 
    627627                any_mag |= (pars[p] != 0) 
    628628                if first_mag is None: 

  • sasmodels/convert.py

    r78f8308 rb3f4831  
    167167    if version < (4, 2, 0): 
    168168        oldpars = _hand_convert_4_1_to_4_2(name, oldpars) 
     169        oldpars = _rename_magnetic_pars(oldpars) 
    169170    return oldpars 
    170171 
    … …  
    174175        oldpars['lattice_distortion'] = oldpars.pop('d_factor') 
    175176    return oldpars 
     177 
     178def _rename_magnetic_pars(pars): 
     179    """ 
     180    Change from M0:par to par_M0, etc. 
     181    """ 
     182    keys = list(pars.items()) 
     183    for k in keys: 
     184        if k.startswith('M0:'): 
     185            pars[k[3:]+'_M0'] = pars.pop(k) 
     186        elif k.startswith('mtheta:'): 
     187            pars[k[7:]+'_mtheta'] = pars.pop(k) 
     188        elif k.startswith('mphi:'): 
     189            pars[k[5:]+'_mphi'] = pars.pop(k) 
     190        elif k.startswith('up:'): 
     191            pars['up_'+k[3:]] = pars.pop(k) 
     192    return pars 
    176193 
    177194def _hand_convert_3_1_2_to_4_1(name, oldpars): 

  • sasmodels/custom/__init__.py

    r0f48f1e rd321747  
    1212import sys 
    1313import os 
    14 from os.path import basename, splitext 
     14from os.path import basename, splitext, join as joinpath, exists, dirname 
    1515 
    1616try: 
    … …  
    1818    from importlib.util import spec_from_file_location, module_from_spec  # type: ignore 
    1919    def load_module_from_path(fullname, path): 
     20        # type: (str, str) -> "module" 
    2021        """load module from *path* as *fullname*""" 
    2122        spec = spec_from_file_location(fullname, os.path.expanduser(path)) 
    … …  
    2728    import imp 
    2829    def load_module_from_path(fullname, path): 
     30        # type: (str, str) -> "module" 
    2931        """load module from *path* as *fullname*""" 
    3032        # Clear out old definitions, if any 
    … …  
    3739        return module 
    3840 
     41_MODULE_CACHE = {} # type: Dict[str, Tuple("module", int)] 
     42_MODULE_DEPENDS = {} # type: Dict[str, List[str]] 
     43_MODULE_DEPENDS_STACK = [] # type: List[str] 
    3944def load_custom_kernel_module(path): 
     45    # type: str -> "module" 
    4046    """load SAS kernel from *path* as *sasmodels.custom.modelname*""" 
    4147    # Pull off the last .ext if it exists; there may be others 
    4248    name = basename(splitext(path)[0]) 
    43     # Placing the model in the 'sasmodels.custom' name space. 
    44     kernel_module = load_module_from_path('sasmodels.custom.'+name, 
    45                                           os.path.expanduser(path)) 
    46     return kernel_module 
     49    path = os.path.expanduser(path) 
     50 
     51    # Reload module if necessary. 
     52    if need_reload(path): 
     53        # Assume the module file is the only dependency 
     54        _MODULE_DEPENDS[path] = set([path]) 
     55 
     56        # Load the module while pushing it onto the dependency stack.  If 
     57        # this triggers any submodules, then they will add their dependencies 
     58        # to this module as the "working_on" parent.  Pop the stack when the 
     59        # module is loaded. 
     60        _MODULE_DEPENDS_STACK.append(path) 
     61        module = load_module_from_path('sasmodels.custom.'+name, path) 
     62        _MODULE_DEPENDS_STACK.pop() 
     63 
     64        # Include external C code in the dependencies.  We are looking 
     65        # for module.source and assuming that it is a list of C source files 
     66        # relative to the module itself.  Any files that do not exist, 
     67        # such as those in the standard libraries, will be ignored. 
     68        # TODO: look in builtin module path for standard c sources 
     69        # TODO: share code with generate.model_sources 
     70        c_sources = getattr(module, 'source', None) 
     71        if isinstance(c_sources, (list, tuple)): 
     72            _MODULE_DEPENDS[path].update(_find_sources(path, c_sources)) 
     73 
     74        # Cache the module, and tag it with the newest timestamp 
     75        timestamp = max(os.path.getmtime(f) for f in _MODULE_DEPENDS[path]) 
     76        _MODULE_CACHE[path] = module, timestamp 
     77 
     78        #print("loading", os.path.basename(path), _MODULE_CACHE[path][1], 
     79        #    [os.path.basename(p) for p in _MODULE_DEPENDS[path]]) 
     80 
     81    # Add path and all its dependence to the parent module, if there is one. 
     82    if _MODULE_DEPENDS_STACK: 
     83        working_on = _MODULE_DEPENDS_STACK[-1] 
     84        _MODULE_DEPENDS[working_on].update(_MODULE_DEPENDS[path]) 
     85 
     86    return _MODULE_CACHE[path][0] 
     87 
     88def need_reload(path): 
     89    # type: str -> bool 
     90    """ 
     91    Return True if any path dependencies have a timestamp newer than the time 
     92    when the path was most recently loaded. 
     93    """ 
     94    # TODO: fails if a dependency has a modification time in the future 
     95    # If the newest dependency has a time stamp in the future, then this 
     96    # will be recorded as the cached time.  When a second dependency 
     97    # is updated to the current time stamp, it will still be considered 
     98    # older than the current build and the reload will not be triggered. 
     99    # Could instead treat all future times as 0 here and in the code above 
     100    # which records the newest timestamp.  This will force a reload when 
     101    # the future time is reached, but other than that should perform 
     102    # correctly.  Probably not worth the extra code... 
     103    _, cache_time = _MODULE_CACHE.get(path, (None, -1)) 
     104    depends = _MODULE_DEPENDS.get(path, [path]) 
     105    #print("reload", any(cache_time < os.path.getmtime(p) for p in depends)) 
     106    #for f in depends: print(">>>  ", f, os.path.getmtime(f)) 
     107    return any(cache_time < os.path.getmtime(p) for p in depends) 
     108 
     109def _find_sources(path, source_list): 
     110    # type: (str, List[str]) -> List[str] 
     111    """ 
     112    Return a list of the sources relative to base file; ignore any that 
     113    are not found. 
     114    """ 
     115    root = dirname(path) 
     116    found = [] 
     117    for source_name in source_list: 
     118        source_path = joinpath(root, source_name) 
     119        if exists(source_path): 
     120            found.append(source_path) 
     121    return found 

  • sasmodels/kernelpy.py

    r108e70e r12eec1e  
    3737        self.info = model_info 
    3838        self.dtype = np.dtype('d') 
    39         logger.info("load python model " + self.info.name) 
     39        logger.info("make python model " + self.info.name) 
    4040 
    4141    def make_kernel(self, q_vectors): 

  • sasmodels/model_test.py

    r012cd34 r12eec1e  
    4747import sys 
    4848import unittest 
     49import traceback 
    4950 
    5051try: 
    … …  
    7475# pylint: enable=unused-import 
    7576 
    76  
    7777def make_suite(loaders, models): 
    7878    # type: (List[str], List[str]) -> unittest.TestSuite 
    … …  
    8686    *models* is the list of models to test, or *["all"]* to test all models. 
    8787    """ 
    88     ModelTestCase = _hide_model_case_from_nose() 
    8988    suite = unittest.TestSuite() 
    9089 
    … …  
    9594        skip = [] 
    9695    for model_name in models: 
    97         if model_name in skip: 
    98             continue 
    99         model_info = load_model_info(model_name) 
    100  
    101         #print('------') 
    102         #print('found tests in', model_name) 
    103         #print('------') 
    104  
    105         # if ispy then use the dll loader to call pykernel 
    106         # don't try to call cl kernel since it will not be 
    107         # available in some environmentes. 
    108         is_py = callable(model_info.Iq) 
    109  
    110         # Some OpenCL drivers seem to be flaky, and are not producing the 
    111         # expected result.  Since we don't have known test values yet for 
    112         # all of our models, we are instead going to compare the results 
    113         # for the 'smoke test' (that is, evaluation at q=0.1 for the default 
    114         # parameters just to see that the model runs to completion) between 
    115         # the OpenCL and the DLL.  To do this, we define a 'stash' which is 
    116         # shared between OpenCL and DLL tests.  This is just a list.  If the 
    117         # list is empty (which it will be when DLL runs, if the DLL runs 
    118         # first), then the results are appended to the list.  If the list 
    119         # is not empty (which it will be when OpenCL runs second), the results 
    120         # are compared to the results stored in the first element of the list. 
    121         # This is a horrible stateful hack which only makes sense because the 
    122         # test suite is thrown away after being run once. 
    123         stash = [] 
    124  
    125         if is_py:  # kernel implemented in python 
    126             test_name = "%s-python"%model_name 
    127             test_method_name = "test_%s_python" % model_info.id 
     96        if model_name not in skip: 
     97            model_info = load_model_info(model_name) 
     98            _add_model_to_suite(loaders, suite, model_info) 
     99 
     100    return suite 
     101 
     102def _add_model_to_suite(loaders, suite, model_info): 
     103    ModelTestCase = _hide_model_case_from_nose() 
     104 
     105    #print('------') 
     106    #print('found tests in', model_name) 
     107    #print('------') 
     108 
     109    # if ispy then use the dll loader to call pykernel 
     110    # don't try to call cl kernel since it will not be 
     111    # available in some environmentes. 
     112    is_py = callable(model_info.Iq) 
     113 
     114    # Some OpenCL drivers seem to be flaky, and are not producing the 
     115    # expected result.  Since we don't have known test values yet for 
     116    # all of our models, we are instead going to compare the results 
     117    # for the 'smoke test' (that is, evaluation at q=0.1 for the default 
     118    # parameters just to see that the model runs to completion) between 
     119    # the OpenCL and the DLL.  To do this, we define a 'stash' which is 
     120    # shared between OpenCL and DLL tests.  This is just a list.  If the 
     121    # list is empty (which it will be when DLL runs, if the DLL runs 
     122    # first), then the results are appended to the list.  If the list 
     123    # is not empty (which it will be when OpenCL runs second), the results 
     124    # are compared to the results stored in the first element of the list. 
     125    # This is a horrible stateful hack which only makes sense because the 
     126    # test suite is thrown away after being run once. 
     127    stash = [] 
     128 
     129    if is_py:  # kernel implemented in python 
     130        test_name = "%s-python"%model_info.name 
     131        test_method_name = "test_%s_python" % model_info.id 
     132        test = ModelTestCase(test_name, model_info, 
     133                                test_method_name, 
     134                                platform="dll",  # so that 
     135                                dtype="double", 
     136                                stash=stash) 
     137        suite.addTest(test) 
     138    else:   # kernel implemented in C 
     139 
     140        # test using dll if desired 
     141        if 'dll' in loaders or not use_opencl(): 
     142            test_name = "%s-dll"%model_info.name 
     143            test_method_name = "test_%s_dll" % model_info.id 
    128144            test = ModelTestCase(test_name, model_info, 
    129                                  test_method_name, 
    130                                  platform="dll",  # so that 
    131                                  dtype="double", 
    132                                  stash=stash) 
     145                                    test_method_name, 
     146                                    platform="dll", 
     147                                    dtype="double", 
     148                                    stash=stash) 
    133149            suite.addTest(test) 
    134         else:   # kernel implemented in C 
    135  
    136             # test using dll if desired 
    137             if 'dll' in loaders or not use_opencl(): 
    138                 test_name = "%s-dll"%model_name 
    139                 test_method_name = "test_%s_dll" % model_info.id 
    140                 test = ModelTestCase(test_name, model_info, 
    141                                      test_method_name, 
    142                                      platform="dll", 
    143                                      dtype="double", 
    144                                      stash=stash) 
    145                 suite.addTest(test) 
    146  
    147             # test using opencl if desired and available 
    148             if 'opencl' in loaders and use_opencl(): 
    149                 test_name = "%s-opencl"%model_name 
    150                 test_method_name = "test_%s_opencl" % model_info.id 
    151                 # Using dtype=None so that the models that are only 
    152                 # correct for double precision are not tested using 
    153                 # single precision.  The choice is determined by the 
    154                 # presence of *single=False* in the model file. 
    155                 test = ModelTestCase(test_name, model_info, 
    156                                      test_method_name, 
    157                                      platform="ocl", dtype=None, 
    158                                      stash=stash) 
    159                 #print("defining", test_name) 
    160                 suite.addTest(test) 
    161  
    162     return suite 
     150 
     151        # test using opencl if desired and available 
     152        if 'opencl' in loaders and use_opencl(): 
     153            test_name = "%s-opencl"%model_info.name 
     154            test_method_name = "test_%s_opencl" % model_info.id 
     155            # Using dtype=None so that the models that are only 
     156            # correct for double precision are not tested using 
     157            # single precision.  The choice is determined by the 
     158            # presence of *single=False* in the model file. 
     159            test = ModelTestCase(test_name, model_info, 
     160                                    test_method_name, 
     161                                    platform="ocl", dtype=None, 
     162                                    stash=stash) 
     163            #print("defining", test_name) 
     164            suite.addTest(test) 
     165 
    163166 
    164167def _hide_model_case_from_nose(): 
    … …  
    348351    return abs(target-actual)/shift < 1.5*10**-digits 
    349352 
    350 def run_one(model): 
    351     # type: (str) -> str 
    352     """ 
    353     Run the tests for a single model, printing the results to stdout. 
    354  
    355     *model* can by a python file, which is handy for checking user defined 
    356     plugin models. 
     353# CRUFT: old interface; should be deprecated and removed 
     354def run_one(model_name): 
     355    # msg = "use check_model(model_info) rather than run_one(model_name)" 
     356    # warnings.warn(msg, category=DeprecationWarning, stacklevel=2) 
     357    try: 
     358        model_info = load_model_info(model_name) 
     359    except Exception: 
     360        output = traceback.format_exc() 
     361        return output 
     362 
     363    success, output = check_model(model_info) 
     364    return output 
     365 
     366def check_model(model_info): 
     367    # type: (ModelInfo) -> str 
     368    """ 
     369    Run the tests for a single model, capturing the output. 
     370 
     371    Returns success status and the output string. 
    357372    """ 
    358373    # Note that running main() directly did not work from within the 
    … …  
    369384    # Build a test suite containing just the model 
    370385    loaders = ['opencl'] if use_opencl() else ['dll'] 
    371     models = [model] 
    372     try: 
    373         suite = make_suite(loaders, models) 
    374     except Exception: 
    375         import traceback 
    376         stream.writeln(traceback.format_exc()) 
    377         return 
     386    suite = unittest.TestSuite() 
     387    _add_model_to_suite(loaders, suite, model_info) 
    378388 
    379389    # Warn if there are no user defined tests. 
    … …  
    390400    for test in suite: 
    391401        if not test.info.tests: 
    392             stream.writeln("Note: %s has no user defined tests."%model) 
     402            stream.writeln("Note: %s has no user defined tests."%model_info.name) 
    393403        break 
    394404    else: 
    … …  
    406416    output = stream.getvalue() 
    407417    stream.close() 
    408     return output 
     418    return result.wasSuccessful(), output 
    409419 
    410420 

  • sasmodels/modelinfo.py

    r7b9e4dd rbd547d0  
    466466        self.is_asymmetric = any(p.name == 'psi' for p in self.kernel_parameters) 
    467467        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    468                                 if p.id.startswith('M0:')] 
     468                                if p.id.endswith('_M0')] 
    469469 
    470470        self.pd_1d = set(p.name for p in self.call_parameters 
    … …  
    586586        if self.nmagnetic > 0: 
    587587            full_list.extend([ 
    588                 Parameter('up:frac_i', '', 0., [0., 1.], 
     588                Parameter('up_frac_i', '', 0., [0., 1.], 
    589589                          'magnetic', 'fraction of spin up incident'), 
    590                 Parameter('up:frac_f', '', 0., [0., 1.], 
     590                Parameter('up_frac_f', '', 0., [0., 1.], 
    591591                          'magnetic', 'fraction of spin up final'), 
    592                 Parameter('up:angle', 'degrees', 0., [0., 360.], 
     592                Parameter('up_angle', 'degrees', 0., [0., 360.], 
    593593                          'magnetic', 'spin up angle'), 
    594594            ]) 
    … …  
    596596            for p in slds: 
    597597                full_list.extend([ 
    598                     Parameter('M0:'+p.id, '1e-6/Ang^2', 0., [-np.inf, np.inf], 
     598                    Parameter(p.id+'_M0', '1e-6/Ang^2', 0., [-np.inf, np.inf], 
    599599                              'magnetic', 'magnetic amplitude for '+p.description), 
    600                     Parameter('mtheta:'+p.id, 'degrees', 0., [-90., 90.], 
     600                    Parameter(p.id+'_mtheta', 'degrees', 0., [-90., 90.], 
    601601                              'magnetic', 'magnetic latitude for '+p.description), 
    602                     Parameter('mphi:'+p.id, 'degrees', 0., [-180., 180.], 
     602                    Parameter(p.id+'_mphi', 'degrees', 0., [-180., 180.], 
    603603                              'magnetic', 'magnetic longitude for '+p.description), 
    604604                ]) 
    … …  
    640640 
    641641        Parameters marked as sld will automatically have a set of associated 
    642         magnetic parameters (m0:p, mtheta:p, mphi:p), as well as polarization 
    643         information (up:theta, up:frac_i, up:frac_f). 
     642        magnetic parameters (p_M0, p_mtheta, p_mphi), as well as polarization 
     643        information (up_theta, up_frac_i, up_frac_f). 
    644644        """ 
    645645        # control parameters go first 
    … …  
    668668            result.append(expanded_pars[name]) 
    669669            if is2d: 
    670                 for tag in 'M0:', 'mtheta:', 'mphi:': 
    671                     if tag+name in expanded_pars: 
    672                         result.append(expanded_pars[tag+name]) 
     670                for tag in '_M0', '_mtheta', '_mphi': 
     671                    if name+tag in expanded_pars: 
     672                        result.append(expanded_pars[name+tag]) 
    673673 
    674674        # Gather the user parameters in order 
    … …  
    703703                append_group(p.id) 
    704704 
    705         if is2d and 'up:angle' in expanded_pars: 
     705        if is2d and 'up_angle' in expanded_pars: 
    706706            result.extend([ 
    707                 expanded_pars['up:frac_i'], 
    708                 expanded_pars['up:frac_f'], 
    709                 expanded_pars['up:angle'], 
     707                expanded_pars['up_frac_i'], 
     708                expanded_pars['up_frac_f'], 
     709                expanded_pars['up_angle'], 
    710710            ]) 
    711711 
    … …  
    793793    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
    794794    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
     795    # Note: custom.load_custom_kernel_module assumes the C sources are defined 
     796    # by this attribute. 
    795797    info.source = getattr(kernel_module, 'source', []) 
    796798    info.c_code = getattr(kernel_module, 'c_code', None) 
    … …  
    10141016                         for k in range(control+1, p.length+1) 
    10151017                         if p.length > 1) 
     1018            for p in self.parameters.kernel_parameters: 
     1019                if p.length > 1 and p.type == "sld": 
     1020                    for k in range(control+1, p.length+1): 
     1021                        base = p.id+str(k) 
     1022                        hidden.update((base+"_M0", base+"_mtheta", base+"_mphi")) 
    10161023        return hidden 

  • sasmodels/models/bcc_paracrystal.py

    re7e9231 rb3f4831  
    11r""" 
     2.. warning:: This model and this model description are under review following  
     3             concerns raised by SasView users. If you need to use this model,  
     4             please email help@sasview.org for the latest situation. *The  
     5             SasView Developers. September 2018.* 
     6 
    27Definition 
    38---------- 
    … …  
    1318 
    1419    I(q) = \frac{\text{scale}}{V_p} V_\text{lattice} P(q) Z(q) 
    15  
    1620 
    1721where *scale* is the volume fraction of spheres, $V_p$ is the volume of the 
    … …  
    97101 
    98102Authorship and Verification 
    99 ---------------------------- 
     103--------------------------- 
    100104 
    101105* **Author:** NIST IGOR/DANSE **Date:** pre 2010 

  • sasmodels/models/be_polyelectrolyte.py

    ref07e95 rca77fc1  
    11r""" 
     2.. note:: Please read the Validation section below. 
     3 
    24Definition 
    35---------- 
    … …  
    1113 
    1214    I(q) = K\frac{q^2+k^2}{4\pi L_b\alpha ^2} 
    13     \frac{1}{1+r_{0}^2(q^2+k^2)(q^2-12hC_a/b^2)} + background 
     15    \frac{1}{1+r_{0}^4(q^2+k^2)(q^2-12hC_a/b^2)} + background 
    1416 
    1517    k^2 = 4\pi L_b(2C_s + \alpha C_a) 
    1618 
    17     r_{0}^2 = \frac{1}{\alpha \sqrt{C_a} \left( b/\sqrt{48\pi L_b}\right)} 
     19    r_{0}^2 = \frac{b}{\alpha \sqrt{C_a 48\pi L_b}} 
    1820 
    1921where 
    2022 
    2123$K$ is the contrast factor for the polymer which is defined differently than in 
    22 other models and is given in barns where $1 barn = 10^{-24} cm^2$.  $K$ is 
     24other models and is given in barns where 1 $barn = 10^{-24}$ $cm^2$.  $K$ is 
    2325defined as: 
    2426 
    … …  
    2931    a = b_p - (v_p/v_s) b_s 
    3032 
    31 where $b_p$ and $b_s$ are sum of the scattering lengths of the atoms 
    32 constituting the monomer of the polymer and the sum of the scattering lengths 
    33 of the atoms constituting the solvent molecules respectively, and $v_p$ and 
    34 $v_s$ are the partial molar volume of the polymer and the solvent respectively 
    35  
    36 $L_b$ is the Bjerrum length(|Ang|) - **Note:** This parameter needs to be 
    37 kept constant for a given solvent and temperature! 
    38  
    39 $h$ is the virial parameter (|Ang^3|/mol) - **Note:** See [#Borue]_ for the 
    40 correct interpretation of this parameter.  It incorporates second and third 
    41 virial coefficients and can be Negative. 
    42  
    43 $b$ is the monomer length(|Ang|), $C_s$ is the concentration of monovalent 
    44 salt(mol/L), $\alpha$ is the ionization degree (ionization degree : ratio of 
    45 charged monomers  to total number of monomers), $C_a$ is the polymer molar 
    46 concentration(mol/L), and $background$ is the incoherent background. 
     33where: 
     34 
     35- $b_p$ and $b_s$ are **sum of the scattering lengths of the atoms** 
     36  constituting the polymer monomer and the solvent molecules, respectively. 
     37 
     38- $v_p$ and $v_s$ are the partial molar volume of the polymer and the  
     39  solvent, respectively. 
     40 
     41- $L_b$ is the Bjerrum length (|Ang|) - **Note:** This parameter needs to be 
     42  kept constant for a given solvent and temperature! 
     43 
     44- $h$ is the virial parameter (|Ang^3|) - **Note:** See [#Borue]_ for the 
     45  correct interpretation of this parameter.  It incorporates second and third 
     46  virial coefficients and can be *negative*. 
     47 
     48- $b$ is the monomer length (|Ang|). 
     49 
     50- $C_s$ is the concentration of monovalent salt(1/|Ang^3| - internally converted from mol/L). 
     51 
     52- $\alpha$ is the degree of ionization (the ratio of charged monomers to the total  
     53  number of monomers) 
     54 
     55- $C_a$ is the polymer molar concentration (1/|Ang^3| - internally converted from mol/L) 
     56 
     57- $background$ is the incoherent background. 
    4758 
    4859For 2D data the scattering intensity is calculated in the same way as 1D, 
    … …  
    5263 
    5364    q = \sqrt{q_x^2 + q_y^2} 
     65 
     66Validation 
     67---------- 
     68 
     69As of the last revision, this code is believed to be correct.  However it 
     70needs further validation and should be used with caution at this time.  The 
     71history of this code goes back to a 1998 implementation. It was recently noted 
     72that in that implementation, while both the polymer concentration and salt 
     73concentration were converted from experimental units of mol/L to more 
     74dimensionally useful units of 1/|Ang^3|, only the converted version of the 
     75polymer concentration was actually being used in the calculation while the 
     76unconverted salt concentration (still in apparent units of mol/L) was being  
     77used.  This was carried through to Sasmodels as used for SasView 4.1 (though  
     78the line of code converting the salt concentration to the new units was removed  
     79somewhere along the line). Simple dimensional analysis of the calculation shows  
     80that the converted salt concentration should be used, which the original code  
     81suggests was the intention, so this has now been corrected (for SasView 4.2).  
     82Once better validation has been performed this note will be removed. 
    5483 
    5584References 
    … …  
    6695 
    6796* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    68 * **Last Modified by:** Paul Kienzle **Date:** July 24, 2016 
    69 * **Last Reviewed by:** Paul Butler and Richard Heenan **Date:** October 07, 2016 
     97* **Last Modified by:** Paul Butler **Date:** September 25, 2018 
     98* **Last Reviewed by:** Paul Butler **Date:** September 25, 2018 
    7099""" 
    71100 
    … …  
    92121    ["contrast_factor",       "barns",   10.0,  [-inf, inf], "", "Contrast factor of the polymer"], 
    93122    ["bjerrum_length",        "Ang",      7.1,  [0, inf],    "", "Bjerrum length"], 
    94     ["virial_param",          "Ang^3/mol", 12.0,  [-inf, inf], "", "Virial parameter"], 
     123    ["virial_param",          "Ang^3", 12.0,  [-inf, inf], "", "Virial parameter"], 
    95124    ["monomer_length",        "Ang",     10.0,  [0, inf],    "", "Monomer length"], 
    96125    ["salt_concentration",    "mol/L",    0.0,  [-inf, inf], "", "Concentration of monovalent salt"], 
    … …  
    102131 
    103132def Iq(q, 
    104        contrast_factor=10.0, 
    105        bjerrum_length=7.1, 
    106        virial_param=12.0, 
    107        monomer_length=10.0, 
    108        salt_concentration=0.0, 
    109        ionization_degree=0.05, 
    110        polymer_concentration=0.7): 
     133       contrast_factor, 
     134       bjerrum_length, 
     135       virial_param, 
     136       monomer_length, 
     137       salt_concentration, 
     138       ionization_degree, 
     139       polymer_concentration): 
    111140    """ 
    112     :param q:                     Input q-value 
    113     :param contrast_factor:       Contrast factor of the polymer 
    114     :param bjerrum_length:        Bjerrum length 
    115     :param virial_param:          Virial parameter 
    116     :param monomer_length:        Monomer length 
    117     :param salt_concentration:    Concentration of monovalent salt 
    118     :param ionization_degree:     Degree of ionization 
    119     :param polymer_concentration: Polymer molar concentration 
    120     :return:                      1-D intensity 
     141    :params: see parameter table 
     142    :return: 1-D form factor for polyelectrolytes in low salt 
     143     
     144    parameter names, units, default values, and behavior (volume, sld etc) are 
     145    defined in the parameter table.  The concentrations are converted from 
     146    experimental mol/L to dimensionaly useful 1/A3 in first two lines 
    121147    """ 
    122148 
    123     concentration = polymer_concentration * 6.022136e-4 
    124  
    125     k_square = 4.0 * pi * bjerrum_length * (2*salt_concentration + 
    126                                             ionization_degree * concentration) 
    127  
    128     r0_square = 1.0/ionization_degree/sqrt(concentration) * \ 
     149    concentration_pol = polymer_concentration * 6.022136e-4 
     150    concentration_salt = salt_concentration * 6.022136e-4 
     151 
     152    k_square = 4.0 * pi * bjerrum_length * (2*concentration_salt + 
     153                                            ionization_degree * concentration_pol) 
     154 
     155    r0_square = 1.0/ionization_degree/sqrt(concentration_pol) * \ 
    129156                (monomer_length/sqrt((48.0*pi*bjerrum_length))) 
    130157 
    … …  
    133160 
    134161    term2 = 1.0 + r0_square**2 * (q**2 + k_square) * \ 
    135         (q**2 - (12.0 * virial_param * concentration/(monomer_length**2))) 
     162        (q**2 - (12.0 * virial_param * concentration_pol/(monomer_length**2))) 
    136163 
    137164    return term1/term2 
    … …  
    174201 
    175202    # Accuracy tests based on content in test/utest_other_models.py 
     203    # Note that these should some day be validated beyond this self validation 
     204    # (circular reasoning). -- i.e. the "good value," at least for those with 
     205    # non zero salt concentrations, were obtained by running the current 
     206    # model in SasView and copying the appropriate result here. 
     207    #    PDB -- Sep 26, 2018 
    176208    [{'contrast_factor':       10.0, 
    177209      'bjerrum_length':         7.1, 
    … …  
    184216     }, 0.001, 0.0948379], 
    185217 
    186     # Additional tests with larger range of parameters 
    187218    [{'contrast_factor':       10.0, 
    188219      'bjerrum_length':       100.0, 
    189220      'virial_param':           3.0, 
    190       'monomer_length':         1.0, 
    191       'salt_concentration':    10.0, 
    192       'ionization_degree':      2.0, 
    193       'polymer_concentration': 10.0, 
     221      'monomer_length':         5.0, 
     222      'salt_concentration':     1.0, 
     223      'ionization_degree':      0.1, 
     224      'polymer_concentration':  1.0, 
    194225      'background':             0.0, 
    195      }, 0.1, -3.75693800588], 
     226     }, 0.1, 0.253469484], 
    196227 
    197228    [{'contrast_factor':       10.0, 
    198229      'bjerrum_length':       100.0, 
    199230      'virial_param':           3.0, 
    200       'monomer_length':         1.0, 
    201       'salt_concentration':    10.0, 
    202       'ionization_degree':      2.0, 
    203       'polymer_concentration': 10.0, 
    204       'background':           100.0 
    205      }, 5.0, 100.029142149], 
     231      'monomer_length':         5.0, 
     232      'salt_concentration':     1.0, 
     233      'ionization_degree':      0.1, 
     234      'polymer_concentration':  1.0, 
     235      'background':             1.0, 
     236     }, 0.05, 1.738358122], 
    206237 
    207238    [{'contrast_factor':     100.0, 
    208239      'bjerrum_length':       10.0, 
    209       'virial_param':        180.0, 
    210       'monomer_length':        1.0, 
     240      'virial_param':         12.0, 
     241      'monomer_length':       10.0, 
    211242      'salt_concentration':    0.1, 
    212243      'ionization_degree':     0.5, 
    213244      'polymer_concentration': 0.1, 
    214       'background':             0.0, 
    215      }, 200., 1.80664667511e-06], 
     245      'background':           0.01, 
     246     }, 0.5, 0.012881893], 
    216247    ] 

  • sasmodels/models/fcc_paracrystal.py

    re7e9231 rb3f4831  
    33#note - calculation requires double precision 
    44r""" 
     5.. warning:: This model and this model description are under review following 
     6             concerns raised by SasView users. If you need to use this model, 
     7             please email help@sasview.org for the latest situation. *The 
     8             SasView Developers. September 2018.* 
     9 
    510Definition 
    611---------- 
    … …  
    1015negligible, and the size of the paracrystal is infinitely large. 
    1116Paracrystalline distortion is assumed to be isotropic and characterized by 
    12 a Gaussian distribution.  
     17a Gaussian distribution. 
    1318 
    1419The scattering intensity $I(q)$ is calculated as 

  • sasmodels/models/sc_paracrystal.py

    r0b906ea rb3f4831  
    11r""" 
     2.. warning:: This model and this model description are under review following 
     3             concerns raised by SasView users. If you need to use this model, 
     4             please email help@sasview.org for the latest situation. *The 
     5             SasView Developers. September 2018.* 
     6 
    27Definition 
    38---------- 
    … …  
    914by a Gaussian distribution. 
    1015 
    11 he scattering intensity $I(q)$ is calculated as 
     16The scattering intensity $I(q)$ is calculated as 
    1217 
    1318.. math:: 
    … …  
    2126 
    2227Equation (16) of the 1987 reference\ [#CIT1987]_ is used to calculate $Z(q)$, 
    23 using equations (13)-(15) from the 1987 paper\ [#CIT1987]_ for Z1, Z2, and Z3. 
     28using equations (13)-(15) from the 1987 paper\ [#CIT1987]_ for $Z1$, $Z2$, and 
     29$Z3$. 
    2430 
    2531The lattice correction (the occupied volume of the lattice) for a simple cubic 

  • sasmodels/models/spinodal.py

    r475ff58 r93fe8a1  
    1212where $x=q/q_0$, $q_0$ is the peak position, $I_{max}$ is the intensity  
    1313at $q_0$ (parameterised as the $scale$ parameter), and $B$ is a flat  
    14 background. The spinodal wavelength is given by $2\pi/q_0$.  
     14background. The spinodal wavelength, $\Lambda$, is given by $2\pi/q_0$.  
     15 
     16The definition of $I_{max}$ in the literature varies. Hashimoto *et al* (1991)  
     17define it as  
     18 
     19.. math:: 
     20    I_{max} = \Lambda^3\Delta\rho^2 
     21     
     22whereas Meier & Strobl (1987) give  
     23 
     24.. math:: 
     25    I_{max} = V_z\Delta\rho^2 
     26     
     27where $V_z$ is the volume per monomer unit. 
    1528 
    1629The exponent $\gamma$ is equal to $d+1$ for off-critical concentration  
    … …  
    2841 
    2942H. Furukawa. Dynamics-scaling theory for phase-separating unmixing mixtures: 
    30 Growth rates of droplets and scaling properties of autocorrelation functions. 
    31 Physica A 123,497 (1984). 
     43Growth rates of droplets and scaling properties of autocorrelation functions.  
     44Physica A 123, 497 (1984). 
     45 
     46H. Meier & G. Strobl. Small-Angle X-ray Scattering Study of Spinodal  
     47Decomposition in Polystyrene/Poly(styrene-co-bromostyrene) Blends.  
     48Macromolecules 20, 649-654 (1987). 
     49 
     50T. Hashimoto, M. Takenaka & H. Jinnai. Scattering Studies of Self-Assembling  
     51Processes of Polymer Blends in Spinodal Decomposition.  
     52J. Appl. Cryst. 24, 457-466 (1991). 
    3253 
    3354Revision History 
    … …  
    3556 
    3657* **Author:**  Dirk Honecker **Date:** Oct 7, 2016 
    37 * **Revised:** Steve King    **Date:** Sep 7, 2018 
     58* **Revised:** Steve King    **Date:** Oct 25, 2018 
    3859""" 
    3960 

  • sasmodels/sasview_model.py

    rd533590 rce1eed5  
    6262#: set of defined models (standard and custom) 
    6363MODELS = {}  # type: Dict[str, SasviewModelType] 
     64# TODO: remove unused MODEL_BY_PATH cache once sasview no longer references it 
    6465#: custom model {path: model} mapping so we can check timestamps 
    6566MODEL_BY_PATH = {}  # type: Dict[str, SasviewModelType] 
     67#: Track modules that we have loaded so we can determine whether the model 
     68#: has changed since we last reloaded. 
     69_CACHED_MODULE = {}  # type: Dict[str, "module"] 
    6670 
    6771def find_model(modelname): 
    … …  
    106110    Load a custom model given the model path. 
    107111    """ 
    108     model = MODEL_BY_PATH.get(path, None) 
    109     if model is not None and model.timestamp == getmtime(path): 
    110         #logger.info("Model already loaded %s", path) 
    111         return model 
    112  
    113112    #logger.info("Loading model %s", path) 
     113 
     114    # Load the kernel module.  This may already be cached by the loader, so 
     115    # only requires checking the timestamps of the dependents. 
    114116    kernel_module = custom.load_custom_kernel_module(path) 
    115     if hasattr(kernel_module, 'Model'): 
    116         model = kernel_module.Model 
     117 
     118    # Check if the module has changed since we last looked. 
     119    reloaded = kernel_module != _CACHED_MODULE.get(path, None) 
     120    _CACHED_MODULE[path] = kernel_module 
     121 
     122    # Turn the module into a model.  We need to do this in even if the 
     123    # model has already been loaded so that we can determine the model 
     124    # name and retrieve it from the MODELS cache. 
     125    model = getattr(kernel_module, 'Model', None) 
     126    if model is not None: 
    117127        # Old style models do not set the name in the class attributes, so 
    118128        # set it here; this name will be overridden when the object is created 
    … …  
    127137        model_info = modelinfo.make_model_info(kernel_module) 
    128138        model = make_model_from_info(model_info) 
    129     model.timestamp = getmtime(path) 
    130139 
    131140    # If a model name already exists and we are loading a different model, 
    … …  
    143152                    _previous_name, model.name, model.filename) 
    144153 
    145     MODELS[model.name] = model 
    146     MODEL_BY_PATH[path] = model 
    147     return model 
     154    # Only update the model if the module has changed 
     155    if reloaded or model.name not in MODELS: 
     156        MODELS[model.name] = model 
     157 
     158    return MODELS[model.name] 
    148159 
    149160 
    … …  
    372383            hidden.add('background') 
    373384            self._model_info.parameters.defaults['background'] = 0. 
     385 
     386        # Update the parameter lists to exclude any hidden parameters 
     387        self.magnetic_params = tuple(pname for pname in self.magnetic_params 
     388                                     if pname not in hidden) 
     389        self.orientation_params = tuple(pname for pname in self.orientation_params 
     390                                        if pname not in hidden) 
    374391 
    375392        self._persistency_dict = {} 
    … …  
    786803            return value, [value], [1.0] 
    787804 
     805    @classmethod 
     806    def runTests(cls): 
     807        """ 
     808        Run any tests built into the model and captures the test output. 
     809 
     810        Returns success flag and output 
     811        """ 
     812        from .model_test import check_model 
     813        return check_model(cls._model_info) 
     814 
    788815def test_cylinder(): 
    789816    # type: () -> float 
    … …  
    873900    Model = _make_standard_model('sphere') 
    874901    model = Model() 
    875     model.setParam('M0:sld', 8) 
     902    model.setParam('sld_M0', 8) 
    876903    q = np.linspace(-0.35, 0.35, 500) 
    877904    qx, qy = np.meshgrid(q, q) 

  • sasmodels/special.py

    rdf69efa rfba9ca0  
    113113        The standard math function, tgamma(x) is unstable for $x < 1$ 
    114114        on some platforms. 
     115 
     116    sas_gammaln(x): 
     117        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. 
     118 
     119        The standard math function, lgamma(x), is incorrect for single 
     120        precision on some platforms. 
     121 
     122    sas_gammainc(a, x), sas_gammaincc(a, x): 
     123        Incomplete gamma function 
     124        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     125        and complementary incomplete gamma function 
     126        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
    115127 
    116128    sas_erf(x), sas_erfc(x): 
    … …  
    207219from numpy import pi, nan, inf 
    208220from scipy.special import gamma as sas_gamma 
     221from scipy.special import gammaln as sas_gammaln 
     222from scipy.special import gammainc as sas_gammainc 
     223from scipy.special import gammaincc as sas_gammaincc 
    209224from scipy.special import erf as sas_erf 
    210225from scipy.special import erfc as sas_erfc 

  • setup.py

    r1f991d6 r783e76f  
    2929                return version[1:-1] 
    3030    raise RuntimeError("Could not read version from %s/__init__.py"%package) 
     31 
     32install_requires = ['numpy', 'scipy'] 
     33 
     34if sys.platform=='win32' or sys.platform=='cygwin': 
     35    install_requires.append('tinycc') 
    3136 
    3237setup( 
    … …  
    6166        'sasmodels': ['*.c', '*.cl'], 
    6267    }, 
    63     install_requires=[ 
    64     ], 
     68    install_requires=install_requires, 
    6569    extras_require={ 
     70        'full': ['docutils', 'bumps', 'matplotlib'], 
     71        'server': ['bumps'], 
    6672        'OpenCL': ["pyopencl"], 
    67         'Bumps': ["bumps"], 
    68         'TinyCC': ["tinycc"], 
    6973    }, 
    7074    build_requires=['setuptools'], 

  • explore/realspace.py

    r362a66f r5778279  
    99import numpy as np 
    1010from numpy import pi, radians, sin, cos, sqrt 
    11 from numpy.random import poisson, uniform, randn, rand 
     11from numpy.random import poisson, uniform, randn, rand, randint 
    1212from numpy.polynomial.legendre import leggauss 
    1313from scipy.integrate import simps 
    … …  
    7878 
    7979 
     80I3 = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]]) 
     81 
    8082class Shape: 
    81     rotation = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]]) 
     83    rotation = I3 
    8284    center = np.array([0., 0., 0.])[:, None] 
    8385    r_max = None 
     86    lattice_size = np.array((1, 1, 1)) 
     87    lattice_spacing = np.array((1., 1., 1.)) 
     88    lattice_distortion = 0.0 
     89    lattice_rotation = 0.0 
     90    lattice_type = "" 
    8491 
    8592    def volume(self): 
    … …  
    96103 
    97104    def rotate(self, theta, phi, psi): 
    98         self.rotation = rotation(theta, phi, psi) * self.rotation 
     105        if theta != 0. or phi != 0. or psi != 0.: 
     106            self.rotation = rotation(theta, phi, psi) * self.rotation 
    99107        return self 
    100108 
    … …  
    103111        return self 
    104112 
     113    def lattice(self, size=(1, 1, 1), spacing=(2, 2, 2), type="sc", 
     114                distortion=0.0, rotation=0.0): 
     115        self.lattice_size = np.asarray(size, 'i') 
     116        self.lattice_spacing = np.asarray(spacing, 'd') 
     117        self.lattice_type = type 
     118        self.lattice_distortion = distortion 
     119        self.lattice_rotation = rotation 
     120 
    105121    def _adjust(self, points): 
    106         points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 
     122        if self.rotation is I3: 
     123            points = points.T + self.center 
     124        else: 
     125            points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 
     126        if self.lattice_type: 
     127            points = self._apply_lattice(points) 
    107128        return points.T 
    108129 
    109     def r_bins(self, q, over_sampling=1, r_step=0.): 
    110         r_max = min(2 * pi / q[0], self.r_max) 
     130    def r_bins(self, q, over_sampling=10, r_step=0.): 
     131        if self.lattice_type: 
     132            r_max = np.sqrt(np.sum(self.lattice_size*self.lattice_spacing*self.dims)**2)/2 
     133        else: 
     134            r_max = self.r_max 
     135        #r_max = min(2 * pi / q[0], r_max) 
    111136        if r_step == 0.: 
    112137            r_step = 2 * pi / q[-1] / over_sampling 
    113138        #r_step = 0.01 
    114139        return np.arange(r_step, r_max, r_step) 
     140 
     141    def _apply_lattice(self, points): 
     142        """Spread points to different lattice positions""" 
     143        size = self.lattice_size 
     144        spacing = self.lattice_spacing 
     145        shuffle = self.lattice_distortion 
     146        rotate = self.lattice_rotation 
     147        lattice = self.lattice_type 
     148 
     149        if rotate != 0: 
     150            # To vectorize the rotations we will need to unwrap the matrix multiply 
     151            raise NotImplementedError("don't handle rotations yet") 
     152 
     153        # Determine the number of lattice points in the lattice 
     154        shapes_per_cell = 2 if lattice == "bcc" else 4 if lattice == "fcc" else 1 
     155        number_of_lattice_points = np.prod(size) * shapes_per_cell 
     156 
     157        # For each point in the original shape, figure out which lattice point 
     158        # to translate it to.  This is both cell index (i*ny*nz + j*nz  + k) as 
     159        # well as the point in the cell (corner, body center or face center). 
     160        nsamples = points.shape[1] 
     161        lattice_point = randint(number_of_lattice_points, size=nsamples) 
     162 
     163        # Translate the cell index into the i,j,k coordinates of the senter 
     164        cell_index = lattice_point // shapes_per_cell 
     165        center = np.vstack((cell_index//(size[1]*size[2]), 
     166                            (cell_index%(size[1]*size[2]))//size[2], 
     167                            cell_index%size[2])) 
     168        center = np.asarray(center, dtype='d') 
     169        if lattice == "bcc": 
     170            center[:, lattice_point % shapes_per_cell == 1] += [[0.5], [0.5], [0.5]] 
     171        elif lattice == "fcc": 
     172            center[:, lattice_point % shapes_per_cell == 1] += [[0.0], [0.5], [0.5]] 
     173            center[:, lattice_point % shapes_per_cell == 2] += [[0.5], [0.0], [0.5]] 
     174            center[:, lattice_point % shapes_per_cell == 3] += [[0.5], [0.5], [0.0]] 
     175 
     176        # Each lattice point has its own displacement from the ideal position. 
     177        # Not checking that shapes do not overlap if displacement is too large. 
     178        offset = shuffle*(randn(3, number_of_lattice_points) if shuffle < 0.3 
     179                          else rand(3, number_of_lattice_points)) 
     180        center += offset[:, cell_index] 
     181 
     182        # Each lattice point has its own rotation.  Rotate the point prior to 
     183        # applying any displacement. 
     184        # rotation = rotate*(randn(size=(shapes, 3)) if shuffle < 30 else rand(size=(nsamples, 3))) 
     185        # for k in shapes: points[k] = rotation[k]*points[k] 
     186        points += center*(np.array([spacing])*np.array(self.dims)).T 
     187        return points 
    115188 
    116189class Composite(Shape): 
    … …  
    669742    Iq = 100 * np.ones_like(qx) 
    670743    data = Data2D(x=qx, y=qy, z=Iq, dx=None, dy=None, dz=np.sqrt(Iq)) 
    671     data.x_bins = qx[0,:] 
    672     data.y_bins = qy[:,0] 
     744    data.x_bins = qx[0, :] 
     745    data.y_bins = qy[:, 0] 
    673746    data.filename = "fake data" 
    674747 
    … …  
    695768    return shape, fn, fn_xy 
    696769 
    697 def build_sphere(radius=125, rho=2): 
     770DEFAULT_SPHERE_RADIUS = 125 
     771DEFAULT_SPHERE_CONTRAST = 2 
     772def build_sphere(radius=DEFAULT_SPHERE_RADIUS, rho=DEFAULT_SPHERE_CONTRAST): 
    698773    shape = TriaxialEllipsoid(radius, radius, radius, rho) 
    699774    fn = lambda q: sphere_Iq(q, radius)*rho**2 
    … …  
    751826    return shape, fn, fn_xy 
    752827 
    753 def build_cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
    754                   shuffle=0, rotate=0): 
     828def build_sc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     829                        shuffle=0, rotate=0): 
    755830    a, b, c = shape.dims 
    756     shapes = [copy(shape) 
     831    corners= [copy(shape) 
    757832              .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
    758833                     (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
    … …  
    762837              for iy in range(ny) 
    763838              for iz in range(nz)] 
    764     lattice = Composite(shapes) 
     839    lattice = Composite(corners) 
    765840    return lattice 
    766841 
     842def build_bcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     843                      shuffle=0, rotate=0): 
     844    a, b, c = shape.dims 
     845    corners = [copy(shape) 
     846               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     847                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     848                      (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     849               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     850               for ix in range(nx) 
     851               for iy in range(ny) 
     852               for iz in range(nz)] 
     853    centers = [copy(shape) 
     854               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     855                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     856                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     857               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     858               for ix in range(nx) 
     859               for iy in range(ny) 
     860               for iz in range(nz)] 
     861    lattice = Composite(corners + centers) 
     862    return lattice 
     863 
     864def build_fcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     865                      shuffle=0, rotate=0): 
     866    a, b, c = shape.dims 
     867    corners = [copy(shape) 
     868               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     869                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     870                      (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     871               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     872               for ix in range(nx) 
     873               for iy in range(ny) 
     874               for iz in range(nz)] 
     875    faces_a = [copy(shape) 
     876               .shift((ix+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     877                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     878                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     879               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     880               for ix in range(nx) 
     881               for iy in range(ny) 
     882               for iz in range(nz)] 
     883    faces_b = [copy(shape) 
     884               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     885                      (iy+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     886                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     887               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     888               for ix in range(nx) 
     889               for iy in range(ny) 
     890               for iz in range(nz)] 
     891    faces_c = [copy(shape) 
     892               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     893                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     894                      (iz+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     895               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     896               for ix in range(nx) 
     897               for iy in range(ny) 
     898               for iz in range(nz)] 
     899    lattice = Composite(corners + faces_a + faces_b + faces_c) 
     900    return lattice 
    767901 
    768902SHAPE_FUNCTIONS = OrderedDict([ 
    … …  
    775909]) 
    776910SHAPES = list(SHAPE_FUNCTIONS.keys()) 
     911LATTICE_FUNCTIONS = OrderedDict([ 
     912    ("sc", build_sc_lattice), 
     913    ("bcc", build_bcc_lattice), 
     914    ("fcc", build_fcc_lattice), 
     915]) 
     916LATTICE_TYPES = list(LATTICE_FUNCTIONS.keys()) 
    777917 
    778918def check_shape(title, shape, fn=None, show_points=False, 
    … …  
    783923    r = shape.r_bins(q, r_step=r_step) 
    784924    sampling_density = samples / shape.volume 
     925    print("sampling points") 
    785926    rho, points = shape.sample(sampling_density) 
     927    print("calculating Pr") 
    786928    t0 = time.time() 
    787929    Pr = calc_Pr(r, rho-rho_solvent, points) 
    … …  
    792934    import pylab 
    793935    if show_points: 
    794          plot_points(rho, points); pylab.figure() 
     936        plot_points(rho, points); pylab.figure() 
    795937    plot_calc(r, Pr, q, Iq, theory=theory, title=title) 
    796938    pylab.gcf().canvas.set_window_title(title) 
    … …  
    806948    Qx, Qy = np.meshgrid(qx, qy) 
    807949    sampling_density = samples / shape.volume 
     950    print("sampling points") 
    808951    t0 = time.time() 
    809952    rho, points = shape.sample(sampling_density) 
    … …  
    844987                        help='lattice size') 
    845988    parser.add_argument('-z', '--spacing', type=str, default='2,2,2', 
    846                         help='lattice spacing') 
     989                        help='lattice spacing (relative to shape)') 
     990    parser.add_argument('-t', '--type', choices=LATTICE_TYPES, 
     991                        default=LATTICE_TYPES[0], 
     992                        help='lattice type') 
    847993    parser.add_argument('-r', '--rotate', type=float, default=0., 
    848994                        help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise") 
    … …  
    8581004    nx, ny, nz = [int(v) for v in opts.lattice.split(',')] 
    8591005    dx, dy, dz = [float(v) for v in opts.spacing.split(',')] 
    860     shuffle, rotate = opts.shuffle, opts.rotate 
     1006    distortion, rotation = opts.shuffle, opts.rotate 
    8611007    shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 
    862     if nx > 1 or ny > 1 or nz > 1: 
    863         shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 
     1008    view = tuple(float(v) for v in opts.view.split(',')) 
     1009    # If comparing a sphere in a cubic lattice, compare against the 
     1010    # corresponding paracrystalline model. 
     1011    if opts.shape == "sphere" and dx == dy == dz and nx*ny*nz > 1: 
     1012        radius = pars.get('radius', DEFAULT_SPHERE_RADIUS) 
     1013        model_name = opts.type + "_paracrystal" 
     1014        model_pars = { 
     1015            "scale": 1., 
     1016            "background": 0., 
     1017            "lattice_spacing": 2*radius*dx, 
     1018            "lattice_distortion": distortion, 
     1019            "radius": radius, 
     1020            "sld": pars.get('rho', DEFAULT_SPHERE_CONTRAST), 
     1021            "sld_solvent": 0., 
     1022            "theta": view[0], 
     1023            "phi": view[1], 
     1024            "psi": view[2], 
     1025        } 
     1026        fn, fn_xy = wrap_sasmodel(model_name, **model_pars) 
     1027    if nx*ny*nz > 1: 
     1028        if rotation != 0: 
     1029            print("building %s lattice"%opts.type) 
     1030            build_lattice = LATTICE_FUNCTIONS[opts.type] 
     1031            shape = build_lattice(shape, nx, ny, nz, dx, dy, dz, 
     1032                                  distortion, rotation) 
     1033        else: 
     1034            shape.lattice(size=(nx, ny, nz), spacing=(dx, dy, dz), 
     1035                          type=opts.type, 
     1036                          rotation=rotation, distortion=distortion) 
     1037 
    8641038    title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 
    8651039    if opts.dim == 1: 
    … …  
    8671041                    mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 
    8681042    else: 
    869         view = tuple(float(v) for v in opts.view.split(',')) 
    8701043        check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot, 
    8711044                       mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 

  • sasmodels/models/bcc_paracrystal.c

    r108e70e r642046e  
    11static double 
    2 bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     2bcc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion) 
    33{ 
    44    // Equations from Matsuoka 26-27-28, multiplied by |q| 
    … …  
    1717    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
    1818    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
    19     const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     19    const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 
    2020    const double exp_arg = exp(arg); 
    2121    const double Zq = -cube(expm1(2.0*arg)) 
    22         / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
    23           * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
    24           * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     22        / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0) 
     24          * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0)); 
    2525 
    2626#elif 0 
    … …  
    3636    //            = tanh(-a) / [1 - cos(d a_k)/cosh(-a)] 
    3737    // 
    38     const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     38    const double arg = 0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 
    3939    const double sinh_qd = sinh(arg); 
    4040    const double cosh_qd = cosh(arg); 
    41     const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1)) 
    42                     * sinh_qd/(cosh_qd - cos(dnn*a2)) 
    43                     * sinh_qd/(cosh_qd - cos(dnn*a3)); 
     41    const double Zq = sinh_qd/(cosh_qd - cos(lattice_spacing*a1)) 
     42                    * sinh_qd/(cosh_qd - cos(lattice_spacing*a2)) 
     43                    * sinh_qd/(cosh_qd - cos(lattice_spacing*a3)); 
    4444#else 
    45     const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     45    const double arg = 0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 
    4646    const double tanh_qd = tanh(arg); 
    4747    const double cosh_qd = cosh(arg); 
    48     const double Zq = tanh_qd/(1.0 - cos(dnn*a1)/cosh_qd) 
    49                     * tanh_qd/(1.0 - cos(dnn*a2)/cosh_qd) 
    50                     * tanh_qd/(1.0 - cos(dnn*a3)/cosh_qd); 
     48    const double Zq = tanh_qd/(1.0 - cos(lattice_spacing*a1)/cosh_qd) 
     49                    * tanh_qd/(1.0 - cos(lattice_spacing*a2)/cosh_qd) 
     50                    * tanh_qd/(1.0 - cos(lattice_spacing*a3)/cosh_qd); 
    5151#endif 
    5252 
    … …  
    5757// occupied volume fraction calculated from lattice symmetry and sphere radius 
    5858static double 
    59 bcc_volume_fraction(double radius, double dnn) 
     59bcc_volume_fraction(double radius, double lattice_spacing) 
    6060{ 
    61     return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
     61    return 2.0*sphere_volume(radius/lattice_spacing); 
    6262} 
    6363 
    … …  
    6969 
    7070 
    71 static double Iq(double q, double dnn, 
    72     double d_factor, double radius, 
     71static double Iq(double q, double lattice_spacing, 
     72    double lattice_distortion, double radius, 
    7373    double sld, double solvent_sld) 
    7474{ 
    … …  
    9494            const double qa = qab*cos_phi; 
    9595            const double qb = qab*sin_phi; 
    96             const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     96            const double form = bcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
    9797            inner_sum += GAUSS_W[j] * form; 
    9898        } 
    … …  
    103103    const double Zq = outer_sum/(4.0*M_PI); 
    104104    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    105     return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
     105    return bcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    106106} 
    107107 
    108108 
    109109static double Iqabc(double qa, double qb, double qc, 
    110     double dnn, double d_factor, double radius, 
     110    double lattice_spacing, double lattice_distortion, double radius, 
    111111    double sld, double solvent_sld) 
    112112{ 
    113113    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
    114     const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     114    const double Zq = bcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
    115115    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    116     return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
     116    return bcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    117117} 

  • sasmodels/models/fcc_paracrystal.c

    r108e70e r642046e  
    11static double 
    2 fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     2fcc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion) 
    33{ 
    44    // Equations from Matsuoka 17-18-19, multiplied by |q| 
    … …  
    1616    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
    1717    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
    18     const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     18    const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 
    1919    const double exp_arg = exp(arg); 
    2020    const double Zq = -cube(expm1(2.0*arg)) 
    21         / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
    22           * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
    23           * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     21        / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0) 
     22          * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0)); 
    2424 
    2525    return Zq; 
    … …  
    2929// occupied volume fraction calculated from lattice symmetry and sphere radius 
    3030static double 
    31 fcc_volume_fraction(double radius, double dnn) 
     31fcc_volume_fraction(double radius, double lattice_spacing) 
    3232{ 
    33     return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
     33    return 4.0*sphere_volume(radius/lattice_spacing); 
    3434} 
    3535 
    … …  
    4141 
    4242 
    43 static double Iq(double q, double dnn, 
    44   double d_factor, double radius, 
     43static double Iq(double q, double lattice_spacing, 
     44  double lattice_distortion, double radius, 
    4545  double sld, double solvent_sld) 
    4646{ 
    … …  
    6666            const double qa = qab*cos_phi; 
    6767            const double qb = qab*sin_phi; 
    68             const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     68            const double form = fcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
    6969            inner_sum += GAUSS_W[j] * form; 
    7070        } 
    … …  
    7676    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    7777 
    78     return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     78    return fcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    7979} 
    8080 
    8181 
    8282static double Iqabc(double qa, double qb, double qc, 
    83     double dnn, double d_factor, double radius, 
     83    double lattice_spacing, double lattice_distortion, double radius, 
    8484    double sld, double solvent_sld) 
    8585{ 
    8686    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
    8787    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    88     const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor); 
    89     return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     88    const double Zq = fcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
     89    return fcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    9090} 

  • sasmodels/models/sc_paracrystal.c

    r108e70e r6530963  
    11static double 
    2 sc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     2sc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion) 
    33{ 
    44    // Equations from Matsuoka 9-10-11, multiplied by |q| 
    … …  
    1616    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
    1717    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
    18     const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     18    const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 
    1919    const double exp_arg = exp(arg); 
    2020    const double Zq = -cube(expm1(2.0*arg)) 
    21         / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
    22           * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
    23           * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     21        / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0) 
     22          * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0)); 
    2424 
    2525    return Zq; 
    … …  
    2828// occupied volume fraction calculated from lattice symmetry and sphere radius 
    2929static double 
    30 sc_volume_fraction(double radius, double dnn) 
     30sc_volume_fraction(double radius, double lattice_spacing) 
    3131{ 
    32     return sphere_volume(radius/dnn); 
     32    return sphere_volume(radius/lattice_spacing); 
    3333} 
    3434 
    … …  
    4141 
    4242static double 
    43 Iq(double q, double dnn, 
    44     double d_factor, double radius, 
     43Iq(double q, double lattice_spacing, 
     44    double lattice_distortion, double radius, 
    4545    double sld, double solvent_sld) 
    4646{ 
    … …  
    6767            const double qa = qab*cos_phi; 
    6868            const double qb = qab*sin_phi; 
    69             const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
     69            const double form = sc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
    7070            inner_sum += GAUSS_W[j] * form; 
    7171        } 
    … …  
    7777    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    7878 
    79     return sc_volume_fraction(radius, dnn) * Pq * Zq; 
     79    return sc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    8080} 
    8181 
    … …  
    8383static double 
    8484Iqabc(double qa, double qb, double qc, 
    85     double dnn, double d_factor, double radius, 
     85    double lattice_spacing, double lattice_distortion, double radius, 
    8686    double sld, double solvent_sld) 
    8787{ 
    8888    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
    8989    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    90     const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor); 
    91     return sc_volume_fraction(radius, dnn) * Pq * Zq; 
     90    const double Zq = sc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
     91    return sc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    9292} 
Note: See TracChangeset for help on using the changeset viewer.