Changeset aa8c6e0 in sasmodels


Ignore:
Timestamp:
Oct 30, 2018 10:43:35 AM (23 months ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
765d025, 23df833
Parents:
153f8f6 (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 beta_approx

Files:
50 added
106 edited

Legend:

Unmodified
Added
Removed
  • doc/guide/plugin.rst

    r81751c2 raa8c6e0  
    744744    erf, erfc, tgamma, lgamma:  **do not use** 
    745745        Special functions that should be part of the standard, but are missing 
    746         or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma 
    747         instead (see below). Note: lgamma(x) has not yet been tested. 
     746        or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma 
     747        and sas_lgamma instead (see below). 
    748748 
    749749Some non-standard constants and functions are also provided: 
     
    812812        Gamma function sas_gamma\ $(x) = \Gamma(x)$. 
    813813 
    814         The standard math function, tgamma(x) is unstable for $x < 1$ 
     814        The standard math function, tgamma(x), is unstable for $x < 1$ 
    815815        on some platforms. 
    816816 
    817817        :code:`source = ["lib/sas_gamma.c", ...]` 
    818818        (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_) 
     819 
     820    sas_gammaln(x): 
     821        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. 
     822 
     823        The standard math function, lgamma(x), is incorrect for single 
     824        precision on some platforms. 
     825 
     826        :code:`source = ["lib/sas_gammainc.c", ...]` 
     827        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
     828 
     829    sas_gammainc(a, x), sas_gammaincc(a, x): 
     830        Incomplete gamma function 
     831        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     832        and complementary incomplete gamma function 
     833        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     834 
     835        :code:`source = ["lib/sas_gammainc.c", ...]` 
     836        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
    819837 
    820838    sas_erf(x), sas_erfc(x): 
     
    854872        If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively. 
    855873 
     874        Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100]. 
     875 
    856876        The standard math function jn(n, x) is not available on all platforms. 
    857877 
     
    862882        Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$. 
    863883 
     884        Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100]. 
     885 
    864886        This function uses Taylor series for small and large arguments: 
    865887 
    866         For large arguments, 
     888        For large arguments use the following Taylor series, 
    867889 
    868890        .. math:: 
     
    872894             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 
    873895 
    874         For small arguments, 
     896        For small arguments , 
    875897 
    876898        .. math:: 
  • explore/precision.py

    ree60aa7 raa8c6e0  
    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( 
     
    465511lanczos_gamma = """\ 
    466512    const double coeff[] = { 
    467             76.18009172947146,     -86.50532032941677, 
    468             24.01409824083091,     -1.231739572450155, 
     513            76.18009172947146, -86.50532032941677, 
     514            24.01409824083091, -1.231739572450155, 
    469515            0.1208650973866179e-2,-0.5395239384953e-5 
    470516            }; 
     
    477523""" 
    478524add_function( 
    479     name="log gamma(x)", 
     525    name="loggamma(x)", 
    480526    mp_function=mp.loggamma, 
    481527    np_function=scipy.special.gammaln, 
     
    601647 
    602648ALL_FUNCTIONS = set(FUNCTIONS.keys()) 
    603 ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet 
     649ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead 
    604650ALL_FUNCTIONS.discard("3j1/x:taylor") 
    605651ALL_FUNCTIONS.discard("3j1/x:trig") 
     
    617663    -r indicates that the relative error should be plotted (default), 
    618664    -x<range> indicates the steps in x, where <range> is one of the following 
    619       log indicates log stepping in [10^-3, 10^5] (default) 
    620       logq indicates log stepping in [10^-4, 10^1] 
    621       linear indicates linear stepping in [1, 1000] 
    622       zoom indicates linear stepping in [1000, 1010] 
    623       neg indicates linear stepping in [-100.1, 100.1] 
    624 and name is "all" or one of: 
     665        log indicates log stepping in [10^-3, 10^5] (default) 
     666        logq indicates log stepping in [10^-4, 10^1] 
     667        linear indicates linear stepping in [1, 1000] 
     668        zoom indicates linear stepping in [1000, 1010] 
     669        neg indicates linear stepping in [-100.1, 100.1] 
     670        start:stop:n[:stepping] indicates an n-step plot in [start, stop] 
     671            or [10^start, 10^stop] if stepping is "log" (default n=400) 
     672Some functions (notably gammainc/gammaincc) have an additional parameter A 
     673which can be set from the command line as A=value.  Default is A=1. 
     674 
     675Name is one of: 
    625676    """+names) 
    626677    sys.exit(1) 
  • sasmodels/kernelpy.py

    re44432d raa8c6e0  
    4242        self.info = model_info 
    4343        self.dtype = np.dtype('d') 
     44        logger.info("make python model " + self.info.name) 
    4445 
    4546    def make_kernel(self, q_vectors): 
  • sasmodels/model_test.py

    r81751c2 raa8c6e0  
    4747import sys 
    4848import unittest 
     49import traceback 
    4950 
    5051try: 
     
    7677# pylint: enable=unused-import 
    7778 
    78  
    7979def make_suite(loaders, models): 
    8080    # type: (List[str], List[str]) -> unittest.TestSuite 
     
    8787    *models* is the list of models to test, or *["all"]* to test all models. 
    8888    """ 
    89     ModelTestCase = _hide_model_case_from_nose() 
    9089    suite = unittest.TestSuite() 
    9190 
     
    9695        skip = [] 
    9796    for model_name in models: 
    98         if model_name in skip: 
    99             continue 
    100         model_info = load_model_info(model_name) 
    101  
    102         #print('------') 
    103         #print('found tests in', model_name) 
    104         #print('------') 
    105  
    106         # if ispy then use the dll loader to call pykernel 
    107         # don't try to call cl kernel since it will not be 
    108         # available in some environmentes. 
    109         is_py = callable(model_info.Iq) 
    110  
    111         # Some OpenCL drivers seem to be flaky, and are not producing the 
    112         # expected result.  Since we don't have known test values yet for 
    113         # all of our models, we are instead going to compare the results 
    114         # for the 'smoke test' (that is, evaluation at q=0.1 for the default 
    115         # parameters just to see that the model runs to completion) between 
    116         # the OpenCL and the DLL.  To do this, we define a 'stash' which is 
    117         # shared between OpenCL and DLL tests.  This is just a list.  If the 
    118         # list is empty (which it will be when DLL runs, if the DLL runs 
    119         # first), then the results are appended to the list.  If the list 
    120         # is not empty (which it will be when OpenCL runs second), the results 
    121         # are compared to the results stored in the first element of the list. 
    122         # This is a horrible stateful hack which only makes sense because the 
    123         # test suite is thrown away after being run once. 
    124         stash = [] 
    125  
    126         if is_py:  # kernel implemented in python 
    127             test_name = "%s-python"%model_name 
    128             test_method_name = "test_%s_python" % model_info.id 
     97        if model_name not in skip: 
     98            model_info = load_model_info(model_name) 
     99            _add_model_to_suite(loaders, suite, model_info) 
     100 
     101    return suite 
     102 
     103def _add_model_to_suite(loaders, suite, model_info): 
     104    ModelTestCase = _hide_model_case_from_nose() 
     105 
     106    #print('------') 
     107    #print('found tests in', model_name) 
     108    #print('------') 
     109 
     110    # if ispy then use the dll loader to call pykernel 
     111    # don't try to call cl kernel since it will not be 
     112    # available in some environmentes. 
     113    is_py = callable(model_info.Iq) 
     114 
     115    # Some OpenCL drivers seem to be flaky, and are not producing the 
     116    # expected result.  Since we don't have known test values yet for 
     117    # all of our models, we are instead going to compare the results 
     118    # for the 'smoke test' (that is, evaluation at q=0.1 for the default 
     119    # parameters just to see that the model runs to completion) between 
     120    # the OpenCL and the DLL.  To do this, we define a 'stash' which is 
     121    # shared between OpenCL and DLL tests.  This is just a list.  If the 
     122    # list is empty (which it will be when DLL runs, if the DLL runs 
     123    # first), then the results are appended to the list.  If the list 
     124    # is not empty (which it will be when OpenCL runs second), the results 
     125    # are compared to the results stored in the first element of the list. 
     126    # This is a horrible stateful hack which only makes sense because the 
     127    # test suite is thrown away after being run once. 
     128    stash = [] 
     129 
     130    if is_py:  # kernel implemented in python 
     131        test_name = "%s-python"%model_info.name 
     132        test_method_name = "test_%s_python" % model_info.id 
     133        test = ModelTestCase(test_name, model_info, 
     134                                test_method_name, 
     135                                platform="dll",  # so that 
     136                                dtype="double", 
     137                                stash=stash) 
     138        suite.addTest(test) 
     139    else:   # kernel implemented in C 
     140 
     141        # test using dll if desired 
     142        if 'dll' in loaders or not use_opencl(): 
     143            test_name = "%s-dll"%model_info.name 
     144            test_method_name = "test_%s_dll" % model_info.id 
    129145            test = ModelTestCase(test_name, model_info, 
    130                                  test_method_name, 
    131                                  platform="dll",  # so that 
    132                                  dtype="double", 
    133                                  stash=stash) 
     146                                    test_method_name, 
     147                                    platform="dll", 
     148                                    dtype="double", 
     149                                    stash=stash) 
    134150            suite.addTest(test) 
    135         else:   # kernel implemented in C 
    136  
    137             # test using dll if desired 
    138             if 'dll' in loaders: 
    139                 test_name = "%s-dll"%model_name 
    140                 test_method_name = "test_%s_dll" % model_info.id 
    141                 test = ModelTestCase(test_name, model_info, 
    142                                      test_method_name, 
    143                                      platform="dll", 
    144                                      dtype="double", 
    145                                      stash=stash) 
    146                 suite.addTest(test) 
    147  
    148             # test using opencl if desired and available 
    149             if 'opencl' in loaders and use_opencl(): 
    150                 test_name = "%s-opencl"%model_name 
    151                 test_method_name = "test_%s_opencl" % model_info.id 
    152                 # Using dtype=None so that the models that are only 
    153                 # correct for double precision are not tested using 
    154                 # single precision.  The choice is determined by the 
    155                 # presence of *single=False* in the model file. 
    156                 test = ModelTestCase(test_name, model_info, 
    157                                      test_method_name, 
    158                                      platform="ocl", dtype=None, 
    159                                      stash=stash) 
    160                 #print("defining", test_name) 
    161                 suite.addTest(test) 
    162  
    163             # test using cuda if desired and available 
    164             if 'cuda' in loaders and use_cuda(): 
    165                 test_name = "%s-cuda"%model_name 
    166                 test_method_name = "test_%s_cuda" % model_info.id 
    167                 # Using dtype=None so that the models that are only 
    168                 # correct for double precision are not tested using 
    169                 # single precision.  The choice is determined by the 
    170                 # presence of *single=False* in the model file. 
    171                 test = ModelTestCase(test_name, model_info, 
    172                                      test_method_name, 
    173                                      platform="cuda", dtype=None, 
    174                                      stash=stash) 
    175                 #print("defining", test_name) 
    176                 suite.addTest(test) 
    177  
    178     return suite 
     151 
     152        # test using opencl if desired and available 
     153        if 'opencl' in loaders and use_opencl(): 
     154            test_name = "%s-opencl"%model_info.name 
     155            test_method_name = "test_%s_opencl" % model_info.id 
     156            # Using dtype=None so that the models that are only 
     157            # correct for double precision are not tested using 
     158            # single precision.  The choice is determined by the 
     159            # presence of *single=False* in the model file. 
     160            test = ModelTestCase(test_name, model_info, 
     161                                    test_method_name, 
     162                                    platform="ocl", dtype=None, 
     163                                    stash=stash) 
     164            #print("defining", test_name) 
     165            suite.addTest(test) 
     166 
     167        # test using cuda if desired and available 
     168        if 'cuda' in loaders and use_cuda(): 
     169            test_name = "%s-cuda"%model_name 
     170            test_method_name = "test_%s_cuda" % model_info.id 
     171            # Using dtype=None so that the models that are only 
     172            # correct for double precision are not tested using 
     173            # single precision.  The choice is determined by the 
     174            # presence of *single=False* in the model file. 
     175            test = ModelTestCase(test_name, model_info, 
     176                                    test_method_name, 
     177                                    platform="cuda", dtype=None, 
     178                                    stash=stash) 
     179            #print("defining", test_name) 
     180            suite.addTest(test) 
     181 
    179182 
    180183def _hide_model_case_from_nose(): 
     
    403406    return abs(target-actual)/shift < 1.5*10**-digits 
    404407 
    405 def run_one(model): 
    406     # type: (str) -> str 
    407     """ 
    408     Run the tests for a single model, printing the results to stdout. 
    409  
    410     *model* can by a python file, which is handy for checking user defined 
    411     plugin models. 
     408# CRUFT: old interface; should be deprecated and removed 
     409def run_one(model_name): 
     410    # msg = "use check_model(model_info) rather than run_one(model_name)" 
     411    # warnings.warn(msg, category=DeprecationWarning, stacklevel=2) 
     412    try: 
     413        model_info = load_model_info(model_name) 
     414    except Exception: 
     415        output = traceback.format_exc() 
     416        return output 
     417 
     418    success, output = check_model(model_info) 
     419    return output 
     420 
     421def check_model(model_info): 
     422    # type: (ModelInfo) -> str 
     423    """ 
     424    Run the tests for a single model, capturing the output. 
     425 
     426    Returns success status and the output string. 
    412427    """ 
    413428    # Note that running main() directly did not work from within the 
     
    423438 
    424439    # Build a test suite containing just the model 
    425     loader = 'opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll' 
    426     models = [model] 
    427     try: 
    428         suite = make_suite([loader], models) 
    429     except Exception: 
    430         import traceback 
    431         stream.writeln(traceback.format_exc()) 
    432         return 
     440    loaders = ['opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll'] 
     441    suite = unittest.TestSuite() 
     442    _add_model_to_suite(loaders, suite, model_info) 
    433443 
    434444    # Warn if there are no user defined tests. 
     
    445455    for test in suite: 
    446456        if not test.info.tests: 
    447             stream.writeln("Note: %s has no user defined tests."%model) 
     457            stream.writeln("Note: %s has no user defined tests."%model_info.name) 
    448458        break 
    449459    else: 
     
    461471    output = stream.getvalue() 
    462472    stream.close() 
    463     return output 
     473    return result.wasSuccessful(), output 
    464474 
    465475 
  • sasmodels/sasview_model.py

    r39a06c9 raa8c6e0  
    819819            return value, [value], [1.0] 
    820820 
     821    @classmethod 
     822    def runTests(cls): 
     823        """ 
     824        Run any tests built into the model and captures the test output. 
     825 
     826        Returns success flag and output 
     827        """ 
     828        from .model_test import check_model 
     829        return check_model(cls._model_info) 
     830 
    821831def test_cylinder(): 
    822832    # type: () -> float 
  • 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 
  • doc/guide/gpu_setup.rst

    r63602b1 r8b31efa  
    9494Device Selection 
    9595================ 
     96**OpenCL drivers** 
     97 
    9698If you have multiple GPU devices you can tell the program which device to use. 
    9799By default, the program looks for one GPU and one CPU device from available 
     
    104106was used to run the model. 
    105107 
    106 **If you don't want to use OpenCL, you can set** *SAS_OPENCL=None* 
    107 **in your environment settings, and it will only use normal programs.** 
    108  
    109 If you want to use one of the other devices, you can run the following 
     108If you want to use a specific driver and devices, you can run the following 
    110109from the python console:: 
    111110 
     
    115114This will provide a menu of different OpenCL drivers available. 
    116115When one is selected, it will say "set PYOPENCL_CTX=..." 
    117 Use that value as the value of *SAS_OPENCL*. 
     116Use that value as the value of *SAS_OPENCL=driver:device*. 
     117 
     118To use the default OpenCL device (rather than CUDA or None), 
     119set *SAS_OPENCL=opencl*. 
     120 
     121In batch queues, you may need to set *XDG_CACHE_HOME=~/.cache*  
     122(Linux only) to a different directory, depending on how the filesystem  
     123is configured.  You should also set *SAS_DLL_PATH* for CPU-only modules. 
     124 
     125    -DSAS_MODELPATH=path sets directory containing custom models 
     126    -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU device 
     127    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only) 
     128    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler 
     129    -DSAS_OPENMP=1 turns on OpenMP for the DLLs 
     130    -DSAS_DLL_PATH=path sets the path to the compiled modules 
     131 
     132 
     133**CUDA drivers** 
     134 
     135If OpenCL drivers are not available on your system, but NVidia CUDA 
     136drivers are available, then set *SAS_OPENCL=cuda* or 
     137*SAS_OPENCL=cuda:n* for a particular device number *n*.  If no device 
     138number is specified, then the CUDA drivers looks for look for 
     139*CUDA_DEVICE=n* or a file ~/.cuda-device containing n for the device number. 
     140 
     141In batch queues, the SLURM command *sbatch --gres=gpu:1 ...* will set 
     142*CUDA_VISIBLE_DEVICES=n*, which ought to set the correct device 
     143number for *SAS_OPENCL=cuda*.  If not, then set 
     144*CUDA_DEVICE=$CUDA_VISIBLE_DEVICES* within the batch script.  You may 
     145need to set the CUDA cache directory to a folder accessible across the 
     146cluster with *PYCUDA_CACHE_DIR* (or *PYCUDA_DISABLE_CACHE* to disable 
     147caching), and you may need to set environment specific compiler flags 
     148with *PYCUDA_DEFAULT_NVCC_FLAGS*.  You should also set *SAS_DLL_PATH*  
     149for CPU-only modules. 
     150 
     151**No GPU support** 
     152 
     153If you don't want to use OpenCL or CUDA, you can set *SAS_OPENCL=None* 
     154in your environment settings, and it will only use normal programs. 
     155 
     156In batch queues, you may need to set *SAS_DLL_PATH* to a directory 
     157accessible on the compute node. 
     158 
    118159 
    119160Device Testing 
     
    154195*Document History* 
    155196 
    156 | 2017-09-27 Paul Kienzle 
     197| 2018-10-15 Paul Kienzle 
  • sasmodels/compare.py

    r610ef23 r07646b6  
    4141from . import kerneldll 
    4242from . import kernelcl 
     43from . import kernelcuda 
    4344from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4445from .direct_model import DirectModel, get_mesh 
     
    115116    === environment variables === 
    116117    -DSAS_MODELPATH=path sets directory containing custom models 
    117     -DSAS_OPENCL=vendor:device|none sets the target OpenCL device 
     118    -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU device 
    118119    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only) 
    119120    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler 
     
    352353 
    353354    # If it is a list of choices, pick one at random with equal probability 
    354     # In practice, the model specific random generator will override. 
    355355    par = model_info.parameters[name] 
    356     if len(par.limits) > 2:  # choice list 
    357         return np.random.randint(len(par.limits)) 
     356    if par.choices:  # choice list 
     357        return np.random.randint(len(par.choices)) 
    358358 
    359359    # If it is a fixed range, pick from it with equal probability. 
     
    725725        set_integration_size(model_info, ngauss) 
    726726 
    727     if dtype != "default" and not dtype.endswith('!') and not kernelcl.use_opencl(): 
     727    if (dtype != "default" and not dtype.endswith('!')  
     728            and not (kernelcl.use_opencl() or kernelcuda.use_cuda())): 
    728729        raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR) 
    729730 
  • sasmodels/core.py

    r2dcd6e7 r07646b6  
    2121from . import mixture 
    2222from . import kernelpy 
     23from . import kernelcuda 
    2324from . import kernelcl 
    2425from . import kerneldll 
     
    205206 
    206207    numpy_dtype, fast, platform = parse_dtype(model_info, dtype, platform) 
    207  
    208208    source = generate.make_source(model_info) 
    209209    if platform == "dll": 
    210210        #print("building dll", numpy_dtype) 
    211211        return kerneldll.load_dll(source['dll'], model_info, numpy_dtype) 
     212    elif platform == "cuda": 
     213        return kernelcuda.GpuModel(source, model_info, numpy_dtype, fast=fast) 
    212214    else: 
    213215        #print("building ocl", numpy_dtype) 
     
    245247    # type: (ModelInfo, str, str) -> (np.dtype, bool, str) 
    246248    """ 
    247     Interpret dtype string, returning np.dtype and fast flag. 
     249    Interpret dtype string, returning np.dtype, fast flag and platform. 
    248250 
    249251    Possible types include 'half', 'single', 'double' and 'quad'.  If the 
     
    253255    default for the model and platform. 
    254256 
    255     Platform preference can be specfied ("ocl" vs "dll"), with the default 
    256     being OpenCL if it is availabe.  If the dtype name ends with '!' then 
    257     platform is forced to be DLL rather than OpenCL. 
     257    Platform preference can be specfied ("ocl", "cuda", "dll"), with the 
     258    default being OpenCL or CUDA if available, otherwise DLL.  If the dtype 
     259    name ends with '!' then platform is forced to be DLL rather than GPU. 
     260    The default platform is set by the environment variable SAS_OPENCL, 
     261    SAS_OPENCL=driver:device for OpenCL, SAS_OPENCL=cuda:device for CUDA 
     262    or SAS_OPENCL=none for DLL. 
    258263 
    259264    This routine ignores the preferences within the model definition.  This 
     
    265270    # Assign default platform, overriding ocl with dll if OpenCL is unavailable 
    266271    # If opencl=False OpenCL is switched off 
    267  
    268272    if platform is None: 
    269273        platform = "ocl" 
    270     if not kernelcl.use_opencl() or not model_info.opencl: 
    271         platform = "dll" 
    272274 
    273275    # Check if type indicates dll regardless of which platform is given 
     
    275277        platform = "dll" 
    276278        dtype = dtype[:-1] 
     279 
     280    # Make sure model allows opencl/gpu 
     281    if not model_info.opencl: 
     282        platform = "dll" 
     283 
     284    # Make sure opencl is available, or fallback to cuda then to dll 
     285    if platform == "ocl" and not kernelcl.use_opencl(): 
     286        platform = "cuda" if kernelcuda.use_cuda() else "dll" 
    277287 
    278288    # Convert special type names "half", "fast", and "quad" 
     
    285295        dtype = "float16" 
    286296 
    287     # Convert dtype string to numpy dtype. 
     297    # Convert dtype string to numpy dtype.  Use single precision for GPU 
     298    # if model allows it, otherwise use double precision. 
    288299    if dtype is None or dtype == "default": 
    289         numpy_dtype = (generate.F32 if platform == "ocl" and model_info.single 
     300        numpy_dtype = (generate.F32 if model_info.single and platform in ("ocl", "cuda") 
    290301                       else generate.F64) 
    291302    else: 
    292303        numpy_dtype = np.dtype(dtype) 
    293304 
    294     # Make sure that the type is supported by opencl, otherwise use dll 
     305    # Make sure that the type is supported by GPU, otherwise use dll 
    295306    if platform == "ocl": 
    296307        env = kernelcl.environment() 
    297         if not env.has_type(numpy_dtype): 
    298             platform = "dll" 
    299             if dtype is None: 
    300                 numpy_dtype = generate.F64 
     308    elif platform == "cuda": 
     309        env = kernelcuda.environment() 
     310    else: 
     311        env = None 
     312    if env is not None and not env.has_type(numpy_dtype): 
     313        platform = "dll" 
     314        if dtype is None: 
     315            numpy_dtype = generate.F64 
    301316 
    302317    return numpy_dtype, fast, platform 
     
    365380    model = load_model("cylinder@hardsphere*sphere") 
    366381    actual = [p.name for p in model.info.parameters.kernel_parameters] 
    367     target = ("sld sld_solvent radius length theta phi volfraction" 
     382    target = ("sld sld_solvent radius length theta phi" 
     383              " radius_effective volfraction " 
     384              " structure_factor_mode radius_effective_mode" 
    368385              " A_sld A_sld_solvent A_radius").split() 
    369386    assert target == actual, "%s != %s"%(target, actual) 
  • sasmodels/data.py

    rbd7630d rc88f983  
    337337    else: 
    338338        dq = resolution * q 
    339  
    340339    data = Data1D(q, Iq, dx=dq, dy=dIq) 
    341340    data.filename = "fake data" 
  • sasmodels/direct_model.py

    r7b9e4dd r304c775  
    6464    return calculator(call_details, values, cutoff, is_magnetic) 
    6565 
    66 def call_ER(model_info, pars): 
    67     # type: (ModelInfo, ParameterSet) -> float 
    68     """ 
    69     Call the model ER function using *values*. 
    70  
    71     *model_info* is either *model.info* if you have a loaded model, 
    72     or *kernel.info* if you have a model kernel prepared for evaluation. 
    73     """ 
    74     if model_info.ER is None: 
    75         return 1.0 
    76     elif not model_info.parameters.form_volume_parameters: 
    77         # handle the case where ER is provided but model is not polydisperse 
    78         return model_info.ER() 
    79     else: 
    80         value, weight = _vol_pars(model_info, pars) 
    81         individual_radii = model_info.ER(*value) 
    82         return np.sum(weight*individual_radii) / np.sum(weight) 
    83  
    84  
    85 def call_VR(model_info, pars): 
    86     # type: (ModelInfo, ParameterSet) -> float 
    87     """ 
    88     Call the model VR function using *pars*. 
    89  
    90     *model_info* is either *model.info* if you have a loaded model, 
    91     or *kernel.info* if you have a model kernel prepared for evaluation. 
    92     """ 
    93     if model_info.VR is None: 
    94         return 1.0 
    95     elif not model_info.parameters.form_volume_parameters: 
    96         # handle the case where ER is provided but model is not polydisperse 
    97         return model_info.VR() 
    98     else: 
    99         value, weight = _vol_pars(model_info, pars) 
    100         whole, part = model_info.VR(*value) 
    101         return np.sum(weight*part)/np.sum(weight*whole) 
    102  
     66def call_Fq(calculator, pars, cutoff=0., mono=False): 
     67    # type: (Kernel, ParameterSet, float, bool) -> np.ndarray 
     68    """ 
     69    Like :func:`call_kernel`, but returning F, F^2, R_eff, V, V_form/V_shell. 
     70    """ 
     71    R_eff_type = int(pars.pop('radius_effective_type', 1.0)) 
     72    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
     73    #print("pars", list(zip(*mesh))[0]) 
     74    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
     75    #print("values:", values) 
     76    return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type) 
    10377 
    10478def call_profile(model_info, **pars): 
     
    419393    model_name = sys.argv[1] 
    420394    call = sys.argv[2].upper() 
    421     if call != "ER_VR": 
    422         try: 
    423             values = [float(v) for v in call.split(',')] 
    424         except ValueError: 
    425             values = [] 
    426         if len(values) == 1: 
    427             q, = values 
    428             data = empty_data1D([q]) 
    429         elif len(values) == 2: 
    430             qx, qy = values 
    431             data = empty_data2D([qx], [qy]) 
    432         else: 
    433             print("use q or qx,qy or ER or VR") 
    434             sys.exit(1) 
     395    try: 
     396        values = [float(v) for v in call.split(',')] 
     397    except ValueError: 
     398        values = [] 
     399    if len(values) == 1: 
     400        q, = values 
     401        data = empty_data1D([q]) 
     402    elif len(values) == 2: 
     403        qx, qy = values 
     404        data = empty_data2D([qx], [qy]) 
    435405    else: 
    436         data = empty_data1D([0.001])  # Data not used in ER/VR 
     406        print("use q or qx,qy") 
     407        sys.exit(1) 
    437408 
    438409    model_info = load_model_info(model_name) 
     
    442413                for pair in sys.argv[3:] 
    443414                for k, v in [pair.split('=')]) 
    444     if call == "ER_VR": 
    445         ER = call_ER(model_info, pars) 
    446         VR = call_VR(model_info, pars) 
    447         print(ER, VR) 
    448     else: 
    449         Iq = calculator(**pars) 
    450         print(Iq[0]) 
     415    Iq = calculator(**pars) 
     416    print(Iq[0]) 
    451417 
    452418if __name__ == "__main__": 
  • sasmodels/generate.py

    r6e45516 r39a06c9  
    2424    dimension, or 1.0 if no volume normalization is required. 
    2525 
    26     *ER(p1, p2, ...)* returns the effective radius of the form with 
    27     particular dimensions. 
    28  
    29     *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms. 
     26    *shell_volume(p1, p2, ...)* returns the volume of the shell for forms 
     27    which are hollow. 
     28 
     29    *effective_radius(mode, p1, p2, ...)* returns the effective radius of 
     30    the form with particular dimensions.  Mode determines the type of 
     31    effective radius returned, with mode=1 for equivalent volume. 
    3032 
    3133    #define INVALID(v) (expr)  returns False if v.parameter is invalid 
     
    7274background.  This will work correctly even when polydispersity is off. 
    7375 
    74 *ER* and *VR* are python functions which operate on parameter vectors. 
    75 The constructor code will generate the necessary vectors for computing 
    76 them with the desired polydispersity. 
    7776The kernel module must set variables defining the kernel meta data: 
    7877 
     
    106105    create the OpenCL kernel functions.  The files defining the functions 
    107106    need to be listed before the files which use the functions. 
    108  
    109     *ER* is a python function defining the effective radius.  If it is 
    110     not present, the effective radius is 0. 
    111  
    112     *VR* is a python function defining the volume ratio.  If it is not 
    113     present, the volume ratio is 1. 
    114107 
    115108    *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing 
     
    671664 
    672665 
    673 # type in IQXY pattern could be single, float, double, long double, ... 
    674 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
     666_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]", 
    675667                           flags=re.MULTILINE) 
    676668def find_xy_mode(source): 
     
    701693    return 'qa' 
    702694 
     695# Note: not presently used.  Need to know whether Fq is available before 
     696# trying to compile the source.  Leave the code here in case we decide that 
     697# define have_Fq for each form factor is too tedious and error prone. 
     698_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE) 
     699def contains_Fq(source): 
     700    # type: (List[str]) -> bool 
     701    """ 
     702    Return True if C source defines "void Fq(". 
     703    """ 
     704    for code in source: 
     705        m = _FQ_PATTERN.search(code) 
     706        if m is not None: 
     707            return True 
     708    return False 
     709 
     710_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE) 
     711def contains_shell_volume(source): 
     712    # type: (List[str]) -> bool 
     713    """ 
     714    Return True if C source defines "void Fq(". 
     715    """ 
     716    for code in source: 
     717        m = _SHELL_VOLUME_PATTERN.search(code) 
     718        if m is not None: 
     719            return True 
     720    return False 
    703721 
    704722def _add_source(source, code, path, lineno=1): 
     
    730748    # dispersion.  Need to be careful that necessary parameters are available 
    731749    # for computing volume even if we allow non-disperse volume parameters. 
    732  
    733750    partable = model_info.parameters 
    734751 
     
    743760    for path, code in user_code: 
    744761        _add_source(source, code, path) 
    745  
    746762    if model_info.c_code: 
    747763        _add_source(source, model_info.c_code, model_info.filename, 
     
    751767    q, qx, qy, qab, qa, qb, qc \ 
    752768        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
     769 
    753770    # Generate form_volume function, etc. from body only 
    754771    if isinstance(model_info.form_volume, str): 
    755772        pars = partable.form_volume_parameters 
    756773        source.append(_gen_fn(model_info, 'form_volume', pars)) 
     774    if isinstance(model_info.shell_volume, str): 
     775        pars = partable.form_volume_parameters 
     776        source.append(_gen_fn(model_info, 'shell_volume', pars)) 
    757777    if isinstance(model_info.Iq, str): 
    758778        pars = [q] + partable.iq_parameters 
     
    767787        pars = [qa, qb, qc] + partable.iq_parameters 
    768788        source.append(_gen_fn(model_info, 'Iqabc', pars)) 
     789 
     790    # Check for shell_volume in source 
     791    is_hollow = contains_shell_volume(source) 
    769792 
    770793    # What kind of 2D model do we need?  Is it consistent with the parameters? 
     
    789812    source.append("\\\n".join(p.as_definition() 
    790813                              for p in partable.kernel_parameters)) 
    791  
    792814    # Define the function calls 
     815    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0" 
    793816    if partable.form_volume_parameters: 
    794817        refs = _call_pars("_v.", partable.form_volume_parameters) 
    795         call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs)) 
     818        if is_hollow: 
     819            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2) 
     820        else: 
     821            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs)) 
     822        if model_info.effective_radius_type: 
     823            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs)) 
    796824    else: 
    797825        # Model doesn't have volume.  We could make the kernel run a little 
    798826        # faster by not using/transferring the volume normalizations, but 
    799827        # the ifdef's reduce readability more than is worthwhile. 
    800         call_volume = "#define CALL_VOLUME(v) 1.0" 
     828        call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)" 
    801829    source.append(call_volume) 
    802  
     830    source.append(call_effective_radius) 
    803831    model_refs = _call_pars("_v.", partable.iq_parameters) 
    804     pars = ",".join(["_q"] + model_refs) 
    805     call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     832 
     833    if model_info.have_Fq: 
     834        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs) 
     835        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars 
     836        clear_iq = "#undef CALL_FQ" 
     837    else: 
     838        pars = ",".join(["_q"] + model_refs) 
     839        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     840        clear_iq = "#undef CALL_IQ" 
    806841    if xy_mode == 'qabc': 
    807842        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
     
    812847        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    813848        clear_iqxy = "#undef CALL_IQ_AC" 
    814     elif xy_mode == 'qa': 
     849    elif xy_mode == 'qa' and not model_info.have_Fq: 
    815850        pars = ",".join(["_qa"] + model_refs) 
    816851        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    817852        clear_iqxy = "#undef CALL_IQ_A" 
     853    elif xy_mode == 'qa' and model_info.have_Fq: 
     854        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs) 
     855        # Note: uses rare C construction (expr1, expr2) which computes 
     856        # expr1 then expr2 and evaluates to expr2.  This allows us to 
     857        # leave it looking like a function even though it is returning 
     858        # its values by reference. 
     859        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars 
     860        clear_iqxy = "#undef CALL_FQ_A" 
    818861    elif xy_mode == 'qxy': 
    819862        orientation_refs = _call_pars("_v.", partable.orientation_parameters) 
     
    831874    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    832875               if p.type == 'sld'] 
    833  
    834876    # Fill in definitions for numbers of parameters 
    835877    source.append("#define MAX_PD %s"%partable.max_pd) 
     
    839881    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    840882    source.append("#define PROJECTION %d"%PROJECTION) 
    841  
    842883    # TODO: allow mixed python/opencl kernels? 
    843  
    844     ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    845     dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
     884    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
     885    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
     886 
    846887    result = { 
    847888        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    848889        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    849890    } 
    850  
    851891    return result 
    852892 
    853893 
    854 def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
     894def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name): 
    855895    # type: ([str,str], str, str, str) -> List[str] 
    856896    code = kernel[0] 
     
    862902        '#line 1 "%s Iq"' % path, 
    863903        code, 
    864         "#undef CALL_IQ", 
     904        clear_iq, 
    865905        "#undef KERNEL_NAME", 
    866906        ] 
  • sasmodels/kernel.py

    r2d81cfe re44432d  
    4141    dtype = None  # type: np.dtype 
    4242 
    43     def __call__(self, call_details, values, cutoff, magnetic): 
     43    def Iq(self, call_details, values, cutoff, magnetic): 
    4444        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    45         raise NotImplementedError("need to implement __call__") 
     45        r""" 
     46        Returns I(q) from the polydisperse average scattering. 
     47 
     48        .. math:: 
     49 
     50            I(q) = \text{scale} \cdot P(q) + \text{background} 
     51 
     52        With the correct choice of model and contrast, setting *scale* to 
     53        the volume fraction $V_f$ of particles should match the measured 
     54        absolute scattering.  Some models (e.g., vesicle) have volume fraction 
     55        built into the model, and do not need an additional scale. 
     56        """ 
     57        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic, 
     58                              effective_radius_type=0) 
     59        combined_scale = values[0]/shell_volume 
     60        background = values[1] 
     61        return combined_scale*F2 + background 
     62    __call__ = Iq 
     63 
     64    def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0): 
     65        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     66        r""" 
     67        Returns <F(q)>, <F(q)^2>, effective radius, shell volume and 
     68        form:shell volume ratio.  The <F(q)> term may be None if the 
     69        form factor does not support direct computation of $F(q)$ 
     70 
     71        $P(q) = <F^2(q)>/<V>$ is used for structure factor calculations, 
     72 
     73        .. math:: 
     74 
     75            I(q) = \text{scale} \cdot P(q) \cdot S(q) + \text{background} 
     76 
     77        For the beta approximation, this becomes 
     78 
     79        .. math:: 
     80 
     81            I(q) = \text{scale} * P (1 + <F>^2/<F^2> (S - 1)) + \text{background} 
     82                 = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background} 
     83 
     84        $<F(q)>$ and $<F^2(q)>$ are averaged by polydispersity in shape 
     85        and orientation, with each configuration $x_k$ having form factor 
     86        $F(q, x_k)$, weight $w_k$ and volume $V_k$.  The result is: 
     87 
     88        .. math:: 
     89 
     90            P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} 
     91 
     92        The form factor itself is scaled by volume and contrast to compute the 
     93        total scattering.  This is then squared, and the volume weighted 
     94        F^2 is then normalized by volume F.  For a given density, the number 
     95        of scattering centers is assumed to scale linearly with volume.  Later 
     96        scaling the resulting $P(q)$ by the volume fraction of particles 
     97        gives the total scattering on an absolute scale. Most models 
     98        incorporate the volume fraction into the overall scale parameter.  An 
     99        exception is vesicle, which includes the volume fraction parameter in 
     100        the model itself, scaling $F$ by $\surd V_f$ so that the math for 
     101        the beta approximation works out. 
     102 
     103        By scaling $P(q)$ by total weight $\sum w_k$, there is no need to make 
     104        sure that the polydisperisity distributions normalize to one.  In 
     105        particular, any distibution values $x_k$ outside the valid domain 
     106        of $F$ will not be included, and the distribution will be implicitly 
     107        truncated.  This is controlled by the parameter limits defined in the 
     108        model (which truncate the distribution before calling the kernel) as 
     109        well as any region excluded using the *INVALID* macro defined within 
     110        the model itself. 
     111 
     112        The volume used in the polydispersity calculation is the form volume 
     113        for solid objects or the shell volume for hollow objects.  Shell 
     114        volume should be used within $F$ so that the normalizing scale 
     115        represents the volume fraction of the shell rather than the entire 
     116        form.  This corresponds to the volume fraction of shell-forming 
     117        material added to the solvent. 
     118 
     119        The calculation of $S$ requires the effective radius and the 
     120        volume fraction of the particles.  The model can have several 
     121        different ways to compute effective radius, with the 
     122        *effective_radius_type* parameter used to select amongst them.  The 
     123        volume fraction of particles should be determined from the total 
     124        volume fraction of the form, not just the shell volume fraction. 
     125        This makes a difference for hollow shapes, which need to scale 
     126        the volume fraction by the returned volume ratio when computing $S$. 
     127        For solid objects, the shell volume is set to the form volume so 
     128        this scale factor evaluates to one and so can be used for both 
     129        hollow and solid shapes. 
     130        """ 
     131        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     132        #print("returned",self.q_input.q, self.result) 
     133        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     134        total_weight = self.result[nout*self.q_input.nq + 0] 
     135        if total_weight == 0.: 
     136            total_weight = 1. 
     137        # Note: shell_volume == form_volume for solid objects 
     138        form_volume = self.result[nout*self.q_input.nq + 1]/total_weight 
     139        shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight 
     140        effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight 
     141        if shell_volume == 0.: 
     142            shell_volume = 1. 
     143        F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None 
     144        F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight 
     145        return F1, F2, effective_radius, shell_volume, form_volume/shell_volume 
    46146 
    47147    def release(self): 
    48148        # type: () -> None 
    49149        pass 
     150 
     151    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     152        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     153        """ 
     154        Call the kernel.  Subclasses defining kernels for particular execution 
     155        engines need to provide an implementation for this. 
     156        """ 
     157        raise NotImplementedError() 
  • sasmodels/kernel_header.c

    r108e70e r07646b6  
    11#ifdef __OPENCL_VERSION__ 
    22# define USE_OPENCL 
     3#elif defined(__CUDACC__) 
     4# define USE_CUDA 
    35#elif defined(_OPENMP) 
    46# define USE_OPENMP 
    57#endif 
     8 
     9// Use SAS_DOUBLE to force the use of double even for float kernels 
     10#define SAS_DOUBLE dou ## ble 
    611 
    712// If opencl is not available, then we are compiling a C function 
    813// Note: if using a C++ compiler, then define kernel as extern "C" 
    914#ifdef USE_OPENCL 
     15 
     16   #define USE_GPU 
     17   #define pglobal global 
     18   #define pconstant constant 
     19 
    1020   typedef int int32_t; 
    11 #  if defined(USE_SINCOS) 
    12 #    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar) 
    13 #  else 
    14 #    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) 
    15 #  endif 
     21 
     22   #if defined(USE_SINCOS) 
     23   #  define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar) 
     24   #else 
     25   #  define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) 
     26   #endif 
    1627   // Intel CPU on Mac gives strange values for erf(); on the verified 
    1728   // platforms (intel, nvidia, amd), the cephes erf() is significantly 
     
    2435   #  define erfcf erfc 
    2536   #endif 
    26 #else // !USE_OPENCL 
    27 // Use SAS_DOUBLE to force the use of double even for float kernels 
    28 #  define SAS_DOUBLE dou ## ble 
    29 #  ifdef __cplusplus 
     37 
     38#elif defined(USE_CUDA) 
     39 
     40   #define USE_GPU 
     41   #define local __shared__ 
     42   #define pglobal 
     43   #define constant __constant__ 
     44   #define pconstant const 
     45   #define kernel extern "C" __global__ 
     46 
     47   // OpenCL powr(a,b) = C99 pow(a,b), b >= 0 
     48   // OpenCL pown(a,b) = C99 pow(a,b), b integer 
     49   #define powr(a,b) pow(a,b) 
     50   #define pown(a,b) pow(a,b) 
     51   //typedef int int32_t; 
     52   #if defined(USE_SINCOS) 
     53   #  define SINCOS(angle,svar,cvar) sincos(angle,&svar,&cvar) 
     54   #else 
     55   #  define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) 
     56   #endif 
     57 
     58#else // !USE_OPENCL && !USE_CUDA 
     59 
     60   #define local 
     61   #define pglobal 
     62   #define constant const 
     63   #define pconstant const 
     64 
     65   #ifdef __cplusplus 
    3066      #include <cstdio> 
    3167      #include <cmath> 
     
    5187     #endif 
    5288     inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); } 
    53 else // !__cplusplus 
     89   #else // !__cplusplus 
    5490     #include <inttypes.h>  // C99 guarantees that int32_t types is here 
    5591     #include <stdio.h> 
     
    68104         #define NEED_EXPM1 
    69105         #define NEED_TGAMMA 
     106         #define NEED_CBRT 
    70107         // expf missing from windows? 
    71108         #define expf exp 
     
    76113     #define kernel 
    77114     #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) 
    78 #  endif  // !__cplusplus 
    79 #  define global 
    80 #  define local 
    81 #  define constant const 
    82 // OpenCL powr(a,b) = C99 pow(a,b), b >= 0 
    83 // OpenCL pown(a,b) = C99 pow(a,b), b integer 
    84 #  define powr(a,b) pow(a,b) 
    85 #  define pown(a,b) pow(a,b) 
     115   #endif  // !__cplusplus 
     116   // OpenCL powr(a,b) = C99 pow(a,b), b >= 0 
     117   // OpenCL pown(a,b) = C99 pow(a,b), b integer 
     118   #define powr(a,b) pow(a,b) 
     119   #define pown(a,b) pow(a,b) 
     120 
    86121#endif // !USE_OPENCL 
     122 
     123#if defined(NEED_CBRT) 
     124   #define cbrt(_x) pow(_x, 0.33333333333333333333333) 
     125#endif 
    87126 
    88127#if defined(NEED_EXPM1) 
  • sasmodels/kernel_iq.c

    r70530778 r12f4c19  
    2626//  MAGNETIC_PARS : a comma-separated list of indices to the sld 
    2727//      parameters in the parameter table. 
    28 //  CALL_VOLUME(table) : call the form volume function 
     28//  CALL_VOLUME(form, shell, table) : assign form and shell values 
     29//  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function 
    2930//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
    3031//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
     32//  CALL_FQ(q, F1, F2, table) : call the Fq function for 1D calcs. 
     33//  CALL_FQ_A(q, F1, F2, table) : call the Iq function with |q| for 2D data. 
    3134//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 
    3235//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes 
     
    3639//  PROJECTION : equirectangular=1, sinusoidal=2 
    3740//      see explore/jitter.py for definitions. 
     41 
    3842 
    3943#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
     
    270274#endif // _QABC_SECTION 
    271275 
    272  
    273276// ==================== KERNEL CODE ======================== 
    274  
    275277kernel 
    276278void KERNEL_NAME( 
    277     int32_t nq,                 // number of q values 
    278     const int32_t pd_start,     // where we are in the dispersity loop 
    279     const int32_t pd_stop,      // where we are stopping in the dispersity loop 
    280     global const ProblemDetails *details, 
    281     global const double *values, 
    282     global const double *q, // nq q values, with padding to boundary 
    283     global double *result,  // nq+1 return values, again with padding 
    284     const double cutoff     // cutoff in the dispersity weight product 
     279    int32_t nq,                   // number of q values 
     280    const int32_t pd_start,       // where we are in the dispersity loop 
     281    const int32_t pd_stop,        // where we are stopping in the dispersity loop 
     282    pglobal const ProblemDetails *details, 
     283    pglobal const double *values, // parameter values and distributions 
     284    pglobal const double *q,      // nq q values, with padding to boundary 
     285    pglobal double *result,       // nq+1 return values, again with padding 
     286    const double cutoff,          // cutoff in the dispersity weight product 
     287    int32_t effective_radius_type // which effective radius to compute 
    285288    ) 
    286289{ 
    287 #ifdef USE_OPENCL 
     290#if defined(USE_GPU) 
    288291  // who we are and what element we are working with 
     292  #if defined(USE_OPENCL) 
    289293  const int q_index = get_global_id(0); 
     294  #else // USE_CUDA 
     295  const int q_index = threadIdx.x + blockIdx.x * blockDim.x; 
     296  #endif 
    290297  if (q_index >= nq) return; 
    291298#else 
     
    338345  // 
    339346  // The code differs slightly between opencl and dll since opencl is only 
    340   // seeing one q value (stored in the variable "this_result") while the dll 
     347  // seeing one q value (stored in the variable "this_F2") while the dll 
    341348  // version must loop over all q. 
    342   #ifdef USE_OPENCL 
    343     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    344     double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
    345   #else // !USE_OPENCL 
    346     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     349  #if defined(CALL_FQ) 
     350    double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]); 
     351    double weighted_form = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     352    double weighted_shell = (pd_start == 0 ? 0.0 : result[2*nq+2]); 
     353    double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+3]); 
     354  #else 
     355    double weight_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     356    double weighted_form = (pd_start == 0 ? 0.0 : result[nq+1]); 
     357    double weighted_shell = (pd_start == 0 ? 0.0 : result[nq+2]); 
     358    double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+3]); 
     359  #endif 
     360  #if defined(USE_GPU) 
     361    #if defined(CALL_FQ) 
     362      double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]); 
     363      double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]); 
     364    #else 
     365      double this_F2 = (pd_start == 0 ? 0.0 : result[q_index]); 
     366    #endif 
     367  #else // !USE_GPU 
    347368    if (pd_start == 0) { 
    348369      #ifdef USE_OPENMP 
    349370      #pragma omp parallel for 
    350371      #endif 
    351       for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     372      #if defined(CALL_FQ) 
     373          // 2*nq for F^2,F pairs 
     374          for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0; 
     375      #else 
     376          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     377      #endif 
    352378    } 
    353379    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    354 #endif // !USE_OPENCL 
     380#endif // !USE_GPU 
    355381 
    356382 
     
    375401  const int n4 = pd_length[4]; 
    376402  const int p4 = pd_par[4]; 
    377   global const double *v4 = pd_value + pd_offset[4]; 
    378   global const double *w4 = pd_weight + pd_offset[4]; 
     403  pglobal const double *v4 = pd_value + pd_offset[4]; 
     404  pglobal const double *w4 = pd_weight + pd_offset[4]; 
    379405  int i4 = (pd_start/pd_stride[4])%n4;  // position in level 4 at pd_start 
    380406 
     
    411437          FETCH_Q         // set qx,qy from the q input vector 
    412438          APPLY_ROTATION  // convert qx,qy to qa,qb,qc 
    413           CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...) 
     439          CALL_KERNEL     // F2 = Iqxy(qa, qb, qc, p1, p2, ...) 
    414440 
    415441      ++step;  // increment counter representing position in dispersity mesh 
     
    442468// inner loop and defines the macros that use them. 
    443469 
    444 #if defined(CALL_IQ) 
    445   // unoriented 1D 
     470 
     471#if defined(CALL_FQ) // COMPUTE_F1_F2 is true 
     472  // unoriented 1D returning <F> and <F^2> 
     473  // Note that F1 and F2 are returned from CALL_FQ by reference, and the 
     474  // user of the CALL_KERNEL macro below is assuming that F1 and F2 are defined. 
     475  double qk; 
     476  double F1, F2; 
     477  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     478  #define BUILD_ROTATION() do {} while(0) 
     479  #define APPLY_ROTATION() do {} while(0) 
     480  #define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table) 
     481 
     482#elif defined(CALL_FQ_A) 
     483  // unoriented 2D return <F> and <F^2> 
     484  // Note that the CALL_FQ_A macro is computing _F1_slot and _F2_slot by 
     485  // reference then returning _F2_slot.  We are calling them _F1_slot and 
     486  // _F2_slot here so they don't conflict with _F1 and _F2 in the macro 
     487  // expansion, or with the use of F2 = CALL_KERNEL() when it is used below. 
     488  double qx, qy; 
     489  double _F1_slot, _F2_slot; 
     490  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     491  #define BUILD_ROTATION() do {} while(0) 
     492  #define APPLY_ROTATION() do {} while(0) 
     493  #define CALL_KERNEL() CALL_FQ_A(sqrt(qx*qx+qy*qy),_F1_slot,_F2_slot,local_values.table) 
     494 
     495#elif defined(CALL_IQ) 
     496  // unoriented 1D return <F^2> 
    446497  double qk; 
    447498  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
    448499  #define BUILD_ROTATION() do {} while(0) 
    449500  #define APPLY_ROTATION() do {} while(0) 
    450   #define CALL_KERNEL() CALL_IQ(qk, local_values.table) 
     501  #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 
    451502 
    452503#elif defined(CALL_IQ_A) 
     
    482533  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
    483534  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     535 
    484536#elif defined(CALL_IQ_XY) 
    485537  // direct call to qx,qy calculator 
     
    495547// logic in the IQ_AC and IQ_ABC branches.  This will become more important 
    496548// if we implement more projections, or more complicated projections. 
    497 #if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation 
     549#if defined(CALL_IQ) || defined(CALL_IQ_A) || defined(CALL_FQ) || defined(CALL_FQ_A) 
     550  // no orientation 
    498551  #define APPLY_PROJECTION() const double weight=weight0 
    499552#elif defined(CALL_IQ_XY) // pass orientation to the model 
     
    562615  const int n##_LOOP = details->pd_length[_LOOP]; \ 
    563616  const int p##_LOOP = details->pd_par[_LOOP]; \ 
    564   global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \ 
    565   global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \ 
     617  pglobal const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \ 
     618  pglobal const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \ 
    566619  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP; 
    567620 
     
    587640// Pointers to the start of the dispersity and weight vectors, if needed. 
    588641#if MAX_PD>0 
    589   global const double *pd_value = values + NUM_VALUES; 
    590   global const double *pd_weight = pd_value + details->num_weights; 
     642  pglobal const double *pd_value = values + NUM_VALUES; 
     643  pglobal const double *pd_weight = pd_value + details->num_weights; 
    591644#endif 
    592645 
     
    645698    // Note: weight==0 must always be excluded 
    646699    if (weight > cutoff) { 
    647       pd_norm += weight * CALL_VOLUME(local_values.table); 
     700      double form, shell; 
     701      CALL_VOLUME(form, shell, local_values.table); 
     702      weight_norm += weight; 
     703      weighted_form += weight * form; 
     704      weighted_shell += weight * shell; 
     705      if (effective_radius_type != 0) { 
     706        weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); 
     707      } 
    648708      BUILD_ROTATION(); 
    649709 
    650 #ifndef USE_OPENCL 
     710#if !defined(USE_GPU) 
    651711      // DLL needs to explicitly loop over the q values. 
    652712      #ifdef USE_OPENMP 
     
    654714      #endif 
    655715      for (q_index=0; q_index<nq; q_index++) 
    656 #endif // !USE_OPENCL 
     716#endif // !USE_GPU 
    657717      { 
    658718 
     
    663723        #if defined(MAGNETIC) && NUM_MAGNETIC > 0 
    664724          // Compute the scattering from the magnetic cross sections. 
    665           double scattering = 0.0; 
     725          double F2 = 0.0; 
    666726          const double qsq = qx*qx + qy*qy; 
    667727          if (qsq > 1.e-16) { 
     
    688748//  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz); 
    689749                } 
    690                 scattering += xs_weight * CALL_KERNEL(); 
     750                F2 += xs_weight * CALL_KERNEL(); 
    691751              } 
    692752            } 
    693753          } 
    694754        #else  // !MAGNETIC 
    695           const double scattering = CALL_KERNEL(); 
     755          #if defined(CALL_FQ) 
     756            CALL_KERNEL(); // sets F1 and F2 by reference 
     757          #else 
     758            const double F2 = CALL_KERNEL(); 
     759          #endif 
    696760        #endif // !MAGNETIC 
    697 //printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 
    698  
    699         #ifdef USE_OPENCL 
    700           this_result += weight * scattering; 
     761//printf("q_index:%d %g %g %g %g\n", q_index, F2, weight0); 
     762 
     763        #if defined(USE_GPU) 
     764          #if defined(CALL_FQ) 
     765            this_F2 += weight * F2; 
     766            this_F1 += weight * F1; 
     767          #else 
     768            this_F2 += weight * F2; 
     769          #endif 
    701770        #else // !USE_OPENCL 
    702           result[q_index] += weight * scattering; 
     771          #if defined(CALL_FQ) 
     772            result[2*q_index+0] += weight * F2; 
     773            result[2*q_index+1] += weight * F1; 
     774          #else 
     775            result[q_index] += weight * F2; 
     776          #endif 
    703777        #endif // !USE_OPENCL 
    704778      } 
    705779    } 
    706780  } 
    707  
    708781// close nested loops 
    709782++step; 
     
    724797#endif 
    725798 
    726 // Remember the current result and the updated norm. 
    727 #ifdef USE_OPENCL 
    728   result[q_index] = this_result; 
    729   if (q_index == 0) result[nq] = pd_norm; 
    730 //if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
    731 #else // !USE_OPENCL 
    732   result[nq] = pd_norm; 
    733 //printf("res: %g/%g\n", result[0], pd_norm); 
    734 #endif // !USE_OPENCL 
     799// Remember the results and the updated norm. 
     800#if defined(USE_GPU) 
     801  #if defined(CALL_FQ) 
     802  result[2*q_index+0] = this_F2; 
     803  result[2*q_index+1] = this_F1; 
     804  #else 
     805  result[q_index] = this_F2; 
     806  #endif 
     807  if (q_index == 0) 
     808#endif 
     809  { 
     810#if defined(CALL_FQ) 
     811    result[2*nq] = weight_norm; 
     812    result[2*nq+1] = weighted_form; 
     813    result[2*nq+2] = weighted_shell; 
     814    result[2*nq+3] = weighted_radius; 
     815#else 
     816    result[nq] = weight_norm; 
     817    result[nq+1] = weighted_form; 
     818    result[nq+2] = weighted_shell; 
     819    result[nq+3] = weighted_radius; 
     820#endif 
     821  } 
    735822 
    736823// ** clear the macros in preparation for the next kernel ** 
  • sasmodels/kernelcl.py

    rd86f0fc rf872fd1  
    11""" 
    22GPU driver for C kernels 
     3 
     4TODO: docs are out of date 
    35 
    46There should be a single GPU environment running on the system.  This 
     
    5961 
    6062 
    61 # Attempt to setup opencl. This may fail if the opencl package is not 
     63# Attempt to setup opencl. This may fail if the pyopencl package is not 
    6264# installed or if it is installed but there are no devices available. 
    6365try: 
     
    7476 
    7577from . import generate 
     78from .generate import F32, F64 
    7679from .kernel import KernelModel, Kernel 
    7780 
     
    131134 
    132135def use_opencl(): 
    133     return HAVE_OPENCL and os.environ.get("SAS_OPENCL", "").lower() != "none" 
     136    sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower() 
     137    return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda") 
    134138 
    135139ENV = None 
     
    162166    Return true if device supports the requested precision. 
    163167    """ 
    164     if dtype == generate.F32: 
     168    if dtype == F32: 
    165169        return True 
    166     elif dtype == generate.F64: 
     170    elif dtype == F64: 
    167171        return "cl_khr_fp64" in device.extensions 
    168     elif dtype == generate.F16: 
    169         return "cl_khr_fp16" in device.extensions 
    170172    else: 
     173        # Not supporting F16 type since it isn't accurate enough 
    171174        return False 
    172175 
     
    179182        cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE, 
    180183        queue.device) 
    181  
    182 def _stretch_input(vector, dtype, extra=1e-3, boundary=32): 
    183     # type: (np.ndarray, np.dtype, float, int) -> np.ndarray 
    184     """ 
    185     Stretch an input vector to the correct boundary. 
    186  
    187     Performance on the kernels can drop by a factor of two or more if the 
    188     number of values to compute does not fall on a nice power of two 
    189     boundary.   The trailing additional vector elements are given a 
    190     value of *extra*, and so f(*extra*) will be computed for each of 
    191     them.  The returned array will thus be a subset of the computed array. 
    192  
    193     *boundary* should be a power of 2 which is at least 32 for good 
    194     performance on current platforms (as of Jan 2015).  It should 
    195     probably be the max of get_warp(kernel,queue) and 
    196     device.min_data_type_align_size//4. 
    197     """ 
    198     remainder = vector.size % boundary 
    199     if remainder != 0: 
    200         size = vector.size + (boundary - remainder) 
    201         vector = np.hstack((vector, [extra] * (size - vector.size))) 
    202     return np.ascontiguousarray(vector, dtype=dtype) 
    203  
    204184 
    205185def compile_model(context, source, dtype, fast=False): 
     
    239219    """ 
    240220    GPU context, with possibly many devices, and one queue per device. 
     221 
     222    Because the environment can be reset during a live program (e.g., if the 
     223    user changes the active GPU device in the GUI), everything associated 
     224    with the device context must be cached in the environment and recreated 
     225    if the environment changes.  The *cache* attribute is a simple dictionary 
     226    which holds keys and references to objects, such as compiled kernels and 
     227    allocated buffers.  The running program should check in the cache for 
     228    long lived objects and create them if they are not there.  The program 
     229    should not hold onto cached objects, but instead only keep them active 
     230    for the duration of a function call.  When the environment is destroyed 
     231    then the *release* method for each active cache item is called before 
     232    the environment is freed.  This means that each cl buffer should be 
     233    in its own cache entry. 
    241234    """ 
    242235    def __init__(self): 
    243236        # type: () -> None 
    244237        # find gpu context 
    245         #self.context = cl.create_some_context() 
    246  
    247         self.context = None 
    248         if 'SAS_OPENCL' in os.environ: 
    249             #Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context 
    250             os.environ["PYOPENCL_CTX"] = os.environ["SAS_OPENCL"] 
    251         if 'PYOPENCL_CTX' in os.environ: 
    252             self._create_some_context() 
    253  
    254         if not self.context: 
    255             self.context = _get_default_context() 
     238        context_list = _create_some_context() 
     239 
     240        # Find a context for F32 and for F64 (maybe the same one). 
     241        # F16 isn't good enough. 
     242        self.context = {} 
     243        for dtype in (F32, F64): 
     244            for context in context_list: 
     245                if has_type(context.devices[0], dtype): 
     246                    self.context[dtype] = context 
     247                    break 
     248            else: 
     249                self.context[dtype] = None 
     250 
     251        # Build a queue for each context 
     252        self.queue = {} 
     253        context = self.context[F32] 
     254        self.queue[F32] = cl.CommandQueue(context, context.devices[0]) 
     255        if self.context[F64] == self.context[F32]: 
     256            self.queue[F64] = self.queue[F32] 
     257        else: 
     258            context = self.context[F64] 
     259            self.queue[F64] = cl.CommandQueue(context, context.devices[0]) 
    256260 
    257261        # Byte boundary for data alignment 
    258         #self.data_boundary = max(d.min_data_type_align_size 
    259         #                         for d in self.context.devices) 
    260         self.queues = [cl.CommandQueue(context, context.devices[0]) 
    261                        for context in self.context] 
     262        #self.data_boundary = max(context.devices[0].min_data_type_align_size 
     263        #                         for context in self.context.values()) 
     264 
     265        # Cache for compiled programs, and for items in context 
    262266        self.compiled = {} 
     267        self.cache = {} 
    263268 
    264269    def has_type(self, dtype): 
     
    267272        Return True if all devices support a given type. 
    268273        """ 
    269         return any(has_type(d, dtype) 
    270                    for context in self.context 
    271                    for d in context.devices) 
    272  
    273     def get_queue(self, dtype): 
    274         # type: (np.dtype) -> cl.CommandQueue 
    275         """ 
    276         Return a command queue for the kernels of type dtype. 
    277         """ 
    278         for context, queue in zip(self.context, self.queues): 
    279             if all(has_type(d, dtype) for d in context.devices): 
    280                 return queue 
    281  
    282     def get_context(self, dtype): 
    283         # type: (np.dtype) -> cl.Context 
    284         """ 
    285         Return a OpenCL context for the kernels of type dtype. 
    286         """ 
    287         for context in self.context: 
    288             if all(has_type(d, dtype) for d in context.devices): 
    289                 return context 
    290  
    291     def _create_some_context(self): 
    292         # type: () -> cl.Context 
    293         """ 
    294         Protected call to cl.create_some_context without interactivity.  Use 
    295         this if SAS_OPENCL is set in the environment.  Sets the *context* 
    296         attribute. 
    297         """ 
    298         try: 
    299             self.context = [cl.create_some_context(interactive=False)] 
    300         except Exception as exc: 
    301             warnings.warn(str(exc)) 
    302             warnings.warn("pyopencl.create_some_context() failed") 
    303             warnings.warn("the environment variable 'SAS_OPENCL' might not be set correctly") 
     274        return self.context.get(dtype, None) is not None 
    304275 
    305276    def compile_program(self, name, source, dtype, fast, timestamp): 
     
    318289            del self.compiled[key] 
    319290        if key not in self.compiled: 
    320             context = self.get_context(dtype) 
     291            context = self.context[dtype] 
    321292            logging.info("building %s for OpenCL %s", key, 
    322293                         context.devices[0].name.strip()) 
    323             program = compile_model(self.get_context(dtype), 
     294            program = compile_model(self.context[dtype], 
    324295                                    str(source), dtype, fast) 
    325296            self.compiled[key] = (program, timestamp) 
    326297        return program 
     298 
     299    def free_buffer(self, key): 
     300        if key in self.cache: 
     301            self.cache[key].release() 
     302            del self.cache[key] 
     303 
     304    def __del__(self): 
     305        for v in self.cache.values(): 
     306            release = getattr(v, 'release', lambda: None) 
     307            release() 
     308        self.cache = {} 
     309 
     310_CURRENT_ID = 0 
     311def unique_id(): 
     312    global _CURRENT_ID 
     313    _CURRENT_ID += 1 
     314    return _CURRENT_ID 
     315 
     316def _create_some_context(): 
     317    # type: () -> cl.Context 
     318    """ 
     319    Protected call to cl.create_some_context without interactivity. 
     320 
     321    Uses SAS_OPENCL or PYOPENCL_CTX if they are set in the environment, 
     322    otherwise scans for the most appropriate device using 
     323    :func:`_get_default_context`.  Ignore *SAS_OPENCL=OpenCL*, which 
     324    indicates that an OpenCL device should be used without specifying 
     325    which one (and not a CUDA device, or no GPU). 
     326    """ 
     327    # Assume we do not get here if SAS_OPENCL is None or CUDA 
     328    sas_opencl = os.environ.get('SAS_OPENCL', 'opencl') 
     329    if sas_opencl.lower() != 'opencl': 
     330        # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context 
     331        os.environ["PYOPENCL_CTX"] = sas_opencl 
     332 
     333    if 'PYOPENCL_CTX' in os.environ: 
     334        try: 
     335            return [cl.create_some_context(interactive=False)] 
     336        except Exception as exc: 
     337            warnings.warn(str(exc)) 
     338            warnings.warn("pyopencl.create_some_context() failed") 
     339            warnings.warn("the environment variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set correctly") 
     340 
     341    return _get_default_context() 
    327342 
    328343def _get_default_context(): 
     
    404419        self.dtype = dtype 
    405420        self.fast = fast 
    406         self.program = None # delay program creation 
    407         self._kernels = None 
     421        self.timestamp = generate.ocl_timestamp(self.info) 
     422        self._cache_key = unique_id() 
    408423 
    409424    def __getstate__(self): 
     
    414429        # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None 
    415430        self.info, self.source, self.dtype, self.fast = state 
    416         self.program = None 
    417431 
    418432    def make_kernel(self, q_vectors): 
    419433        # type: (List[np.ndarray]) -> "GpuKernel" 
    420         if self.program is None: 
    421             compile_program = environment().compile_program 
    422             timestamp = generate.ocl_timestamp(self.info) 
    423             self.program = compile_program( 
     434        return GpuKernel(self, q_vectors) 
     435 
     436    @property 
     437    def Iq(self): 
     438        return self._fetch_kernel('Iq') 
     439 
     440    def fetch_kernel(self, name): 
     441        # type: (str) -> cl.Kernel 
     442        """ 
     443        Fetch the kernel from the environment by name, compiling it if it 
     444        does not already exist. 
     445        """ 
     446        gpu = environment() 
     447        key = self._cache_key 
     448        if key not in gpu.cache: 
     449            program = gpu.compile_program( 
    424450                self.info.name, 
    425451                self.source['opencl'], 
    426452                self.dtype, 
    427453                self.fast, 
    428                 timestamp) 
     454                self.timestamp) 
    429455            variants = ['Iq', 'Iqxy', 'Imagnetic'] 
    430456            names = [generate.kernel_name(self.info, k) for k in variants] 
    431             kernels = [getattr(self.program, k) for k in names] 
    432             self._kernels = dict((k, v) for k, v in zip(variants, kernels)) 
    433         is_2d = len(q_vectors) == 2 
    434         if is_2d: 
    435             kernel = [self._kernels['Iqxy'], self._kernels['Imagnetic']] 
     457            kernels = [getattr(program, k) for k in names] 
     458            data = dict((k, v) for k, v in zip(variants, kernels)) 
     459            # keep a handle to program so GC doesn't collect 
     460            data['program'] = program 
     461            gpu.cache[key] = data 
    436462        else: 
    437             kernel = [self._kernels['Iq']]*2 
    438         return GpuKernel(kernel, self.dtype, self.info, q_vectors) 
    439  
    440     def release(self): 
    441         # type: () -> None 
    442         """ 
    443         Free the resources associated with the model. 
    444         """ 
    445         if self.program is not None: 
    446             self.program = None 
    447  
    448     def __del__(self): 
    449         # type: () -> None 
    450         self.release() 
     463            data = gpu.cache[key] 
     464        return data[name] 
    451465 
    452466# TODO: check that we don't need a destructor for buffers which go out of scope 
     
    473487        # type: (List[np.ndarray], np.dtype) -> None 
    474488        # TODO: do we ever need double precision q? 
    475         env = environment() 
    476489        self.nq = q_vectors[0].size 
    477490        self.dtype = np.dtype(dtype) 
     
    482495        # architectures tested so far. 
    483496        if self.is_2d: 
    484             # Note: 16 rather than 15 because result is 1 longer than input. 
    485             width = ((self.nq+16)//16)*16 
     497            width = ((self.nq+15)//16)*16 
    486498            self.q = np.empty((width, 2), dtype=dtype) 
    487499            self.q[:self.nq, 0] = q_vectors[0] 
    488500            self.q[:self.nq, 1] = q_vectors[1] 
    489501        else: 
    490             # Note: 32 rather than 31 because result is 1 longer than input. 
    491             width = ((self.nq+32)//32)*32 
     502            width = ((self.nq+31)//32)*32 
    492503            self.q = np.empty(width, dtype=dtype) 
    493504            self.q[:self.nq] = q_vectors[0] 
    494505        self.global_size = [self.q.shape[0]] 
    495         context = env.get_context(self.dtype) 
    496         #print("creating inputs of size", self.global_size) 
    497         self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    498                              hostbuf=self.q) 
     506        self._cache_key = unique_id() 
     507 
     508    @property 
     509    def q_b(self): 
     510        """Lazy creation of q buffer so it can survive context reset""" 
     511        env = environment() 
     512        key = self._cache_key 
     513        if key not in env.cache: 
     514            context = env.context[self.dtype] 
     515            #print("creating inputs of size", self.global_size) 
     516            buffer = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     517                               hostbuf=self.q) 
     518            env.cache[key] = buffer 
     519        return env.cache[key] 
    499520 
    500521    def release(self): 
    501522        # type: () -> None 
    502523        """ 
    503         Free the memory. 
    504         """ 
    505         if self.q_b is not None: 
    506             self.q_b.release() 
    507             self.q_b = None 
     524        Free the buffer associated with the q value 
     525        """ 
     526        environment().free_buffer(id(self)) 
    508527 
    509528    def __del__(self): 
     
    515534    Callable SAS kernel. 
    516535 
    517     *kernel* is the GpuKernel object to call 
    518  
    519     *model_info* is the module information 
    520  
    521     *q_vectors* is the q vectors at which the kernel should be evaluated 
     536    *model* is the GpuModel object to call 
     537 
     538    The following attributes are defined: 
     539 
     540    *info* is the module information 
    522541 
    523542    *dtype* is the kernel precision 
     543 
     544    *dim* is '1d' or '2d' 
     545 
     546    *result* is a vector to contain the results of the call 
    524547 
    525548    The resulting call method takes the *pars*, a list of values for 
     
    531554    Call :meth:`release` when done with the kernel instance. 
    532555    """ 
    533     def __init__(self, kernel, dtype, model_info, q_vectors): 
     556    def __init__(self, model, q_vectors): 
    534557        # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None 
    535         q_input = GpuInput(q_vectors, dtype) 
    536         self.kernel = kernel 
    537         self.info = model_info 
    538         self.dtype = dtype 
    539         self.dim = '2d' if q_input.is_2d else '1d' 
    540         # plus three for the normalization values 
    541         self.result = np.empty(q_input.nq+1, dtype) 
    542  
    543         # Inputs and outputs for each kernel call 
    544         # Note: res may be shorter than res_b if global_size != nq 
     558        dtype = model.dtype 
     559        self.q_input = GpuInput(q_vectors, dtype) 
     560        self._model = model 
     561        # F16 isn't sufficient, so don't support it 
     562        self._as_dtype = np.float64 if dtype == generate.F64 else np.float32 
     563        self._cache_key = unique_id() 
     564 
     565        # attributes accessed from the outside 
     566        self.dim = '2d' if self.q_input.is_2d else '1d' 
     567        self.info = model.info 
     568        self.dtype = model.dtype 
     569 
     570        # holding place for the returned value 
     571        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     572        extra_q = 4  # total weight, form volume, shell volume and R_eff 
     573        self.result = np.empty(self.q_input.nq*nout+extra_q, dtype) 
     574 
     575    @property 
     576    def _result_b(self): 
     577        """Lazy creation of result buffer so it can survive context reset""" 
    545578        env = environment() 
    546         self.queue = env.get_queue(dtype) 
    547  
    548         self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    549                                   q_input.global_size[0] * dtype.itemsize) 
    550         self.q_input = q_input # allocated by GpuInput above 
    551  
    552         self._need_release = [self.result_b, self.q_input] 
    553         self.real = (np.float32 if dtype == generate.F32 
    554                      else np.float64 if dtype == generate.F64 
    555                      else np.float16 if dtype == generate.F16 
    556                      else np.float32)  # will never get here, so use np.float32 
    557  
    558     def __call__(self, call_details, values, cutoff, magnetic): 
     579        key = self._cache_key 
     580        if key not in env.cache: 
     581            context = env.context[self.dtype] 
     582            width = ((self.result.size+31)//32)*32 * self.dtype.itemsize 
     583            buffer = cl.Buffer(context, mf.READ_WRITE, width) 
     584            env.cache[key] = buffer 
     585        return env.cache[key] 
     586 
     587    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    559588        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    560         context = self.queue.context 
    561         # Arrange data transfer to card 
     589        env = environment() 
     590        queue = env.queue[self._model.dtype] 
     591        context = queue.context 
     592 
     593        # Arrange data transfer to/from card 
     594        q_b = self.q_input.q_b 
     595        result_b = self._result_b 
    562596        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    563597                              hostbuf=call_details.buffer) 
     
    565599                             hostbuf=values) 
    566600 
    567         kernel = self.kernel[1 if magnetic else 0] 
    568         args = [ 
     601        name = 'Iq' if self.dim == '1d' else 'Imagnetic' if magnetic else 'Iqxy' 
     602        kernel = self._model.fetch_kernel(name) 
     603        kernel_args = [ 
    569604            np.uint32(self.q_input.nq), None, None, 
    570             details_b, values_b, self.q_input.q_b, self.result_b, 
    571             self.real(cutoff), 
     605            details_b, values_b, q_b, result_b, 
     606            self._as_dtype(cutoff), 
     607            np.uint32(effective_radius_type), 
    572608        ] 
    573609        #print("Calling OpenCL") 
    574610        #call_details.show(values) 
    575         # Call kernel and retrieve results 
     611        #Call kernel and retrieve results 
    576612        wait_for = None 
    577613        last_nap = time.clock() 
     
    580616            stop = min(start + step, call_details.num_eval) 
    581617            #print("queuing",start,stop) 
    582             args[1:3] = [np.int32(start), np.int32(stop)] 
    583             wait_for = [kernel(self.queue, self.q_input.global_size, None, 
    584                                *args, wait_for=wait_for)] 
     618            kernel_args[1:3] = [np.int32(start), np.int32(stop)] 
     619            wait_for = [kernel(queue, self.q_input.global_size, None, 
     620                               *kernel_args, wait_for=wait_for)] 
    585621            if stop < call_details.num_eval: 
    586622                # Allow other processes to run 
     
    588624                current_time = time.clock() 
    589625                if current_time - last_nap > 0.5: 
    590                     time.sleep(0.05) 
     626                    time.sleep(0.001) 
    591627                    last_nap = current_time 
    592         cl.enqueue_copy(self.queue, self.result, self.result_b) 
     628        cl.enqueue_copy(queue, self.result, result_b, wait_for=wait_for) 
    593629        #print("result", self.result) 
    594630 
     
    598634                v.release() 
    599635 
    600         pd_norm = self.result[self.q_input.nq] 
    601         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    602         background = values[1] 
    603         #print("scale",scale,values[0],self.result[self.q_input.nq],background) 
    604         return scale*self.result[:self.q_input.nq] + background 
    605         # return self.result[:self.q_input.nq] 
    606  
    607636    def release(self): 
    608637        # type: () -> None 
     
    610639        Release resources associated with the kernel. 
    611640        """ 
    612         for v in self._need_release: 
    613             v.release() 
    614         self._need_release = [] 
     641        environment().free_buffer(id(self)) 
     642        self.q_input.release() 
    615643 
    616644    def __del__(self): 
  • sasmodels/kerneldll.py

    r1a3559f re44432d  
    9999    pass 
    100100# pylint: enable=unused-import 
     101 
     102# Compiler output is a byte stream that needs to be decode in python 3 
     103decode = (lambda s: s) if sys.version_info[0] < 3 else (lambda s: s.decode('utf8')) 
    101104 
    102105if "SAS_DLL_PATH" in os.environ: 
     
    184187        subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT) 
    185188    except subprocess.CalledProcessError as exc: 
    186         raise RuntimeError("compile failed.\n%s\n%s"%(command_str, exc.output)) 
     189        output = decode(exc.output) 
     190        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, output)) 
    187191    if not os.path.exists(output): 
    188192        raise RuntimeError("compile failed.  File is in %r"%source) 
     
    315319 
    316320        # int, int, int, int*, double*, double*, double*, double*, double 
    317         argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type] 
     321        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type, ct.c_int32] 
    318322        names = [generate.kernel_name(self.info, variant) 
    319323                 for variant in ("Iq", "Iqxy", "Imagnetic")] 
     
    375379    def __init__(self, kernel, model_info, q_input): 
    376380        # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 
     381        #,model_info,q_input) 
    377382        self.kernel = kernel 
    378383        self.info = model_info 
     
    380385        self.dtype = q_input.dtype 
    381386        self.dim = '2d' if q_input.is_2d else '1d' 
    382         self.result = np.empty(q_input.nq+1, q_input.dtype) 
     387        # leave room for f1/f2 results in case we need to compute beta for 1d models 
     388        nout = 2 if self.info.have_Fq else 1 
     389        # +4 for total weight, shell volume, effective radius, form volume 
     390        self.result = np.empty(q_input.nq*nout + 4, self.dtype) 
    383391        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    384392                     else np.float64 if self.q_input.dtype == generate.F64 
    385393                     else np.float128) 
    386394 
    387     def __call__(self, call_details, values, cutoff, magnetic): 
    388         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    389  
     395    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     396        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    390397        kernel = self.kernel[1 if magnetic else 0] 
    391398        args = [ 
     
    394401            None, # pd_stop pd_stride[MAX_PD] 
    395402            call_details.buffer.ctypes.data, # problem 
    396             values.ctypes.data,  #pars 
    397             self.q_input.q.ctypes.data, #q 
     403            values.ctypes.data,  # pars 
     404            self.q_input.q.ctypes.data, # q 
    398405            self.result.ctypes.data,   # results 
    399406            self.real(cutoff), # cutoff 
     407            effective_radius_type, # cutoff 
    400408        ] 
    401409        #print("Calling DLL") 
     
    407415            kernel(*args) # type: ignore 
    408416 
    409         #print("returned",self.q_input.q, self.result) 
    410         pd_norm = self.result[self.q_input.nq] 
    411         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    412         background = values[1] 
    413         #print("scale",scale,background) 
    414         return scale*self.result[:self.q_input.nq] + background 
    415  
    416417    def release(self): 
    417418        # type: () -> None 
  • sasmodels/mixture.py

    recf895e r39a06c9  
    143143    #model_info.tests = [] 
    144144    #model_info.source = [] 
    145     # Iq, Iqxy, form_volume, ER, VR and sesans 
    146145    # Remember the component info blocks so we can build the model 
    147146    model_info.composition = ('mixture', parts) 
  • sasmodels/modelinfo.py

    rbd547d0 r39a06c9  
    164164    parameter.length = length 
    165165    parameter.length_control = control 
    166  
    167166    return parameter 
    168167 
     
    265264    will have a magnitude and a direction, which may be different from 
    266265    other sld parameters. The volume parameters are used for calls 
    267     to form_volume within the kernel (required for volume normalization) 
    268     and for calls to ER and VR for effective radius and volume ratio 
    269     respectively. 
     266    to form_volume within the kernel (required for volume normalization), 
     267    to shell_volume (for hollow shapes), and to effective_radius (for 
     268    structure factor interactions) respectively. 
    270269 
    271270    *description* is a short description of the parameter.  This will 
     
    424423        self.kernel_parameters = parameters 
    425424        self._set_vector_lengths() 
    426  
    427425        self.npars = sum(p.length for p in self.kernel_parameters) 
    428426        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
     
    431429        if self.nmagnetic: 
    432430            self.nvalues += 3 + 3*self.nmagnetic 
    433  
    434431        self.call_parameters = self._get_call_parameters() 
    435432        self.defaults = self._get_defaults() 
     
    722719 
    723720#: Set of variables defined in the model that might contain C code 
    724 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
     721C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code'] 
    725722 
    726723def _find_source_lines(model_info, kernel_module): 
     
    771768        # Custom sum/multi models 
    772769        return kernel_module.model_info 
     770 
    773771    info = ModelInfo() 
    774772    #print("make parameter table", kernel_module.parameters) 
     
    792790    info.category = getattr(kernel_module, 'category', None) 
    793791    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
     792    # TODO: find Fq by inspection 
     793    info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None) 
     794    info.have_Fq = getattr(kernel_module, 'have_Fq', False) 
    794795    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    795796    # Note: custom.load_custom_kernel_module assumes the C sources are defined 
     
    797798    info.source = getattr(kernel_module, 'source', []) 
    798799    info.c_code = getattr(kernel_module, 'c_code', None) 
     800    info.effective_radius = getattr(kernel_module, 'effective_radius', None) 
    799801    # TODO: check the structure of the tests 
    800802    info.tests = getattr(kernel_module, 'tests', []) 
    801     info.ER = getattr(kernel_module, 'ER', None) # type: ignore 
    802     info.VR = getattr(kernel_module, 'VR', None) # type: ignore 
    803803    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 
     804    info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore 
    804805    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    805806    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 
     
    825826    info.lineno = {} 
    826827    _find_source_lines(info, kernel_module) 
    827  
    828828    return info 
    829829 
     
    918918    #: provided in the file. 
    919919    structure_factor = None # type: bool 
     920    #: True if the model defines an Fq function with signature 
     921    #: void Fq(double q, double *F1, double *F2, ...) 
     922    have_Fq = False 
     923    #: List of options for computing the effective radius of the shape, 
     924    #: or None if the model is not usable as a form factor model. 
     925    effective_radius_type = None   # type: List[str] 
    920926    #: List of C source files used to define the model.  The source files 
    921927    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 
    922928    #: model defines orientation parameters. Files containing the most basic 
    923929    #: functions must appear first in the list, followed by the files that 
    924     #: use those functions.  Form factors are indicated by providing 
    925     #: an :attr:`ER` function. 
     930    #: use those functions. 
    926931    source = None           # type: List[str] 
    927     #: The set of tests that must pass.  The format of the tests is described 
    928     #: in :mod:`model_test`. 
    929     tests = None            # type: List[TestCondition] 
    930     #: Returns the effective radius of the model given its volume parameters. 
    931     #: The presence of *ER* indicates that the model is a form factor model 
    932     #: that may be used together with a structure factor to form an implicit 
    933     #: multiplication model. 
    934     #: 
    935     #: The parameters to the *ER* function must be marked with type *volume*. 
    936     #: in the parameter table.  They will appear in the same order as they 
    937     #: do in the table.  The values passed to *ER* will be vectors, with one 
    938     #: value for each polydispersity condition.  For example, if the model 
    939     #: is polydisperse over both length and radius, then both length and 
    940     #: radius will have the same number of values in the vector, with one 
    941     #: value for each *length X radius*.  If only *radius* is polydisperse, 
    942     #: then the value for *length* will be repeated once for each value of 
    943     #: *radius*.  The *ER* function should return one effective radius for 
    944     #: each parameter set.  Multiplicity parameters will be received as 
    945     #: arrays, with one row per polydispersity condition. 
    946     ER = None               # type: Optional[Callable[[np.ndarray], np.ndarray]] 
    947     #: Returns the occupied volume and the total volume for each parameter set. 
    948     #: See :attr:`ER` for details on the parameters. 
    949     VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]] 
    950     #: Arbitrary C code containing supporting functions, etc., to be inserted 
    951     #: after everything in source.  This can include Iq and Iqxy functions with 
    952     #: the full function signature, including all parameters. 
    953     c_code = None 
     932    #: inline source code, added after all elements of source 
     933    c_code = None           # type: Optional[str] 
    954934    #: Returns the form volume for python-based models.  Form volume is needed 
    955935    #: for volume normalization in the polydispersity integral.  If no 
     
    959939    #: C code, either defined as a string, or in the sources. 
    960940    form_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
     941    #: Returns the shell volume for python-based models.  Form volume and 
     942    #: shell volume are needed for volume normalization in the polydispersity 
     943    #: integral and structure interactions for hollow shapes.  If no 
     944    #: parameters are *volume* parameters, then shell volume is not needed. 
     945    #: For C-based models, (with :attr:`sources` defined, or with :attr:`Iq` 
     946    #: defined using a string containing C code), shell_volume must also be 
     947    #: C code, either defined as a string, or in the sources. 
     948    shell_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
    961949    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined 
    962950    #: by the parameter table.  *Iq* can be defined as a python function, or 
     
    993981    #: Line numbers for symbols defining C code 
    994982    lineno = None           # type: Dict[str, int] 
     983    #: The set of tests that must pass.  The format of the tests is described 
     984    #: in :mod:`model_test`. 
     985    tests = None            # type: List[TestCondition] 
    995986 
    996987    def __init__(self): 
  • sasmodels/models/_spherepy.py

    rca4444f r304c775  
    7171    return 1.333333333333333 * pi * radius ** 3 
    7272 
     73def effective_radius(mode, radius): 
     74    return radius 
     75 
    7376def Iq(q, sld, sld_solvent, radius): 
    7477    #print "q",q 
     
    105108sesans.vectorized = True  # sesans accepts an array of z values 
    106109 
    107 def ER(radius): 
    108     return radius 
    109  
    110 # VR defaults to 1.0 
    111  
    112110demo = dict(scale=1, background=0, 
    113111            sld=6, sld_solvent=1, 
  • sasmodels/models/barbell.c

    r108e70e rd42dd4a  
    6363 
    6464static double 
    65 Iq(double q, double sld, double solvent_sld, 
     65radius_from_volume(double radius_bell, double radius, double length) 
     66{ 
     67    const double vol_barbell = form_volume(radius_bell,radius,length); 
     68    return cbrt(vol_barbell/M_4PI_3); 
     69} 
     70 
     71static double 
     72radius_from_totallength(double radius_bell, double radius, double length) 
     73{ 
     74    const double hdist = sqrt(square(radius_bell) - square(radius)); 
     75    return 0.5*length + hdist + radius_bell; 
     76} 
     77 
     78static double 
     79effective_radius(int mode, double radius_bell, double radius, double length) 
     80{ 
     81    switch (mode) { 
     82    default: 
     83    case 1: // equivalent sphere 
     84        return radius_from_volume(radius_bell, radius , length); 
     85    case 2: // radius 
     86        return radius; 
     87    case 3: // half length 
     88        return 0.5*length; 
     89    case 4: // half total length 
     90        return radius_from_totallength(radius_bell,radius,length); 
     91    } 
     92} 
     93 
     94static void 
     95Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 
    6696    double radius_bell, double radius, double length) 
    6797{ 
     
    72102    const double zm = M_PI_4; 
    73103    const double zb = M_PI_4; 
    74     double total = 0.0; 
     104    double total_F1 = 0.0; 
     105    double total_F2 = 0.0; 
    75106    for (int i = 0; i < GAUSS_N; i++){ 
    76107        const double alpha= GAUSS_Z[i]*zm + zb; 
     
    78109        SINCOS(alpha, sin_alpha, cos_alpha); 
    79110        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    80         total += GAUSS_W[i] * Aq * Aq * sin_alpha; 
     111        total_F1 += GAUSS_W[i] * Aq * sin_alpha; 
     112        total_F2 += GAUSS_W[i] * Aq * Aq * sin_alpha; 
    81113    } 
    82114    // translate dx in [-1,1] to dx in [lower,upper] 
    83     const double form = total*zm; 
     115    const double form_avg = total_F1*zm; 
     116    const double form_squared_avg = total_F2*zm; 
    84117 
    85118    //Contrast 
    86119    const double s = (sld - solvent_sld); 
    87     return 1.0e-4 * s * s * form; 
     120    *F1 = 1.0e-2 * s * form_avg; 
     121    *F2 = 1.0e-4 * s * s * form_squared_avg; 
    88122} 
    89  
    90123 
    91124static double 
  • sasmodels/models/barbell.py

    r2d81cfe ree60aa7  
    116116 
    117117source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] 
     118have_Fq = True 
     119effective_radius_type = [ 
     120    "equivalent sphere", "radius", "half length", "half total length", 
     121    ] 
    118122 
    119123def random(): 
  • sasmodels/models/binary_hard_sphere.c

    r925ad6e r71b751d  
    11double form_volume(void); 
    22 
    3 double Iq(double q,  
     3double Iq(double q, 
    44    double lg_radius, double sm_radius, 
    55    double lg_vol_frac, double sm_vol_frac, 
    66    double lg_sld, double sm_sld, double solvent_sld 
    77    ); 
    8      
     8 
    99void calculate_psfs(double qval, 
    1010    double r2, double nf2, 
     
    2222    double lg_vol_frac, double sm_vol_frac, 
    2323    double lg_sld, double sm_sld, double solvent_sld) 
    24      
    2524{ 
    2625    double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten;       //my local names 
     
    2827    double phi1,phi2,phr,a3; 
    2928    double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2; 
    30      
     29 
    3130    r2 = lg_radius; 
    3231    r1 = sm_radius; 
     
    3635    rho1 = sm_sld; 
    3736    rhos = solvent_sld; 
    38      
    39      
     37 
     38 
    4039    phi = phi1 + phi2; 
    4140    aa = r1/r2; 
     
    4645    // calculate the PSF's here 
    4746    calculate_psfs(q,r2,nf2,aa,phi,&psf11,&psf22,&psf12); 
    48      
     47 
    4948    // /* do form factor calculations  */ 
    50      
     49 
    5150    v1 = M_4PI_3*r1*r1*r1; 
    5251    v2 = M_4PI_3*r2*r2*r2; 
    53      
     52 
    5453    n1 = phi1/v1; 
    5554    n2 = phi2/v2; 
    56      
     55 
    5756    qr1 = r1*q; 
    5857    qr2 = r2*q; 
     
    6867    inten *= 1.0e8; 
    6968    ///*convert rho^2 in 10^-6A to A*/ 
    70     inten *= 1.0e-12;     
     69    inten *= 1.0e-12; 
    7170    return(inten); 
    7271} 
     
    7776    double aa, double phi, 
    7877    double *s11, double *s22, double *s12) 
    79  
    8078{ 
    8179    //  variable qval,r2,nf2,aa,phi,&s11,&s22,&s12 
    82      
     80 
    8381    //   calculate constant terms 
    8482    double s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4; 
     
    8785    double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13; 
    8886    double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1; 
    89      
     87 
    9088    s2 = 2.0*r2; 
    9189//    s1 = aa*s2;  why is this never used?  check original paper? 
     
    108106    gm1=(v1*a1+a3*v2*a2)*.5; 
    109107    gm12=2.*gm1*(1.-aa)/aa; 
    110     //c   
     108    //c 
    111109    //c   calculate the direct correlation functions and print results 
    112110    //c 
    113111    //  do 20 j=1,npts 
    114      
     112 
    115113    yy=qval*s2; 
    116114    //c   calculate direct correlation functions 
     
    123121    t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3; 
    124122    f11=24.*v1*(t1+t2+t3)/ay3; 
    125      
     123 
    126124    //c ------c22 
    127125    y2=yy*yy; 
     
    131129    tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3; 
    132130    f22=24.*v2*(tt1+tt2+tt3)/y3; 
    133      
     131 
    134132    //c   -----c12 
    135133    yl=.5*yy*(1.-aa); 
     
    151149    ttt4=a1*(t41+t42)/y1; 
    152150    f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3); 
    153      
     151 
    154152    c11=f11; 
    155153    c22=f22; 
    156154    c12=f12; 
    157155    *s11=1./(1.+c11-(c12)*c12/(1.+c22)); 
    158     *s22=1./(1.+c22-(c12)*c12/(1.+c11));  
    159     *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));    
    160      
     156    *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
     157    *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12)); 
     158 
    161159    return; 
    162160} 
    163  
  • sasmodels/models/capped_cylinder.c

    r108e70e rd42dd4a  
    8585 
    8686static double 
    87 Iq(double q, double sld, double solvent_sld, 
     87radius_from_volume(double radius, double radius_cap, double length) 
     88{ 
     89    const double vol_cappedcyl = form_volume(radius,radius_cap,length); 
     90    return cbrt(vol_cappedcyl/M_4PI_3); 
     91} 
     92 
     93static double 
     94radius_from_totallength(double radius, double radius_cap, double length) 
     95{ 
     96    const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); 
     97    return 0.5*length + hc; 
     98} 
     99 
     100static double 
     101effective_radius(int mode, double radius, double radius_cap, double length) 
     102{ 
     103    switch (mode) { 
     104    default: 
     105    case 1: // equivalent sphere 
     106        return radius_from_volume(radius, radius_cap, length); 
     107    case 2: // radius 
     108        return radius; 
     109    case 3: // half length 
     110        return 0.5*length; 
     111    case 4: // half total length 
     112        return radius_from_totallength(radius, radius_cap,length); 
     113    } 
     114} 
     115 
     116static void 
     117Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 
    88118    double radius, double radius_cap, double length) 
    89119{ 
     
    94124    const double zm = M_PI_4; 
    95125    const double zb = M_PI_4; 
    96     double total = 0.0; 
     126    double total_F1 = 0.0; 
     127    double total_F2 = 0.0; 
    97128    for (int i=0; i<GAUSS_N ;i++) { 
    98129        const double theta = GAUSS_Z[i]*zm + zb; 
     
    103134        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
    104135        // scale by sin_theta for spherical coord integration 
    105         total += GAUSS_W[i] * Aq * Aq * sin_theta; 
     136        total_F1 += GAUSS_W[i] * Aq * sin_theta; 
     137        total_F2 += GAUSS_W[i] * Aq * Aq * sin_theta; 
    106138    } 
    107139    // translate dx in [-1,1] to dx in [lower,upper] 
    108     const double form = total * zm; 
     140    const double form_avg = total_F1 * zm; 
     141    const double form_squared_avg = total_F2 * zm; 
    109142 
    110143    // Contrast 
    111144    const double s = (sld - solvent_sld); 
    112     return 1.0e-4 * s * s * form; 
     145    *F1 = 1.0e-2 * s * form_avg; 
     146    *F2 = 1.0e-4 * s * s * form_squared_avg; 
    113147} 
    114148 
  • sasmodels/models/capped_cylinder.py

    r2d81cfe ree60aa7  
    136136 
    137137source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"] 
     138have_Fq = True 
     139effective_radius_type = [ 
     140    "equivalent sphere", "radius", "half length", "half total length", 
     141    ] 
    138142 
    139143def random(): 
  • sasmodels/models/core_multi_shell.c

    rc3ccaec rd42dd4a  
    99 
    1010static double 
    11 form_volume(double core_radius, double fp_n, double thickness[]) 
     11outer_radius(double core_radius, double fp_n, double thickness[]) 
    1212{ 
    1313  double r = core_radius; 
     
    1616    r += thickness[i]; 
    1717  } 
    18   return M_4PI_3 * cube(r); 
     18  return r; 
    1919} 
    2020 
    2121static double 
    22 Iq(double q, double core_sld, double core_radius, 
     22form_volume(double core_radius, double fp_n, double thickness[]) 
     23{ 
     24  return M_4PI_3 * cube(outer_radius(core_radius, fp_n, thickness)); 
     25} 
     26 
     27static double 
     28effective_radius(int mode, double core_radius, double fp_n, double thickness[]) 
     29{ 
     30  switch (mode) { 
     31  default: 
     32  case 1: // outer radius 
     33    return outer_radius(core_radius, fp_n, thickness); 
     34  case 2: // core radius 
     35    return core_radius; 
     36  } 
     37} 
     38 
     39static void 
     40Fq(double q, double *F1, double *F2, double core_sld, double core_radius, 
    2341   double solvent_sld, double fp_n, double sld[], double thickness[]) 
    2442{ 
     
    3452  } 
    3553  f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sas_3j1x_x(q*r); 
    36   return f * f * 1.0e-4; 
     54  *F1 = 1e-2 * f; 
     55  *F2 = 1e-4 * f * f; 
    3756} 
  • sasmodels/models/core_multi_shell.py

    r2d81cfe ree60aa7  
    9999 
    100100source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 
     101have_Fq = True 
     102effective_radius_type = ["outer radius", "core radius"] 
    101103 
    102104def random(): 
     
    143145    return np.asarray(z), np.asarray(rho) 
    144146 
    145 def ER(radius, n, thickness): 
    146     """Effective radius""" 
    147     n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse 
    148     return np.sum(thickness[:n], axis=0) + radius 
    149  
    150147demo = dict(sld_core=6.4, 
    151148            radius=60, 
  • sasmodels/models/core_shell_bicelle.c

    r108e70e rd42dd4a  
    3737 
    3838static double 
    39 Iq(double q, 
     39radius_from_volume(double radius, double thick_rim, double thick_face, double length) 
     40{ 
     41    const double volume_bicelle = form_volume(radius,thick_rim,thick_face,length); 
     42    return cbrt(volume_bicelle/M_4PI_3); 
     43} 
     44 
     45static double 
     46radius_from_diagonal(double radius, double thick_rim, double thick_face, double length) 
     47{ 
     48    const double radius_tot = radius + thick_rim; 
     49    const double length_tot = length + 2.0*thick_face; 
     50    return sqrt(radius_tot*radius_tot + 0.25*length_tot*length_tot); 
     51} 
     52 
     53static double 
     54effective_radius(int mode, double radius, double thick_rim, double thick_face, double length) 
     55{ 
     56    switch (mode) { 
     57    default: 
     58    case 1: // equivalent sphere 
     59        return radius_from_volume(radius, thick_rim, thick_face, length); 
     60    case 2: // outer rim radius 
     61        return radius + thick_rim; 
     62    case 3: // half outer thickness 
     63        return 0.5*length + thick_face; 
     64    case 4: // half diagonal 
     65        return radius_from_diagonal(radius,thick_rim,thick_face,length); 
     66    } 
     67} 
     68 
     69static void 
     70Fq(double q, 
     71    double *F1, 
     72    double *F2, 
    4073    double radius, 
    4174    double thick_radius, 
     
    5184    const double halflength = 0.5*length; 
    5285 
    53     double total = 0.0; 
     86    double total_F1 = 0.0; 
     87    double total_F2 = 0.0; 
    5488    for(int i=0;i<GAUSS_N;i++) { 
    5589        double theta = (GAUSS_Z[i] + 1.0)*uplim; 
    5690        double sin_theta, cos_theta; // slots to hold sincos function output 
    5791        SINCOS(theta, sin_theta, cos_theta); 
    58         double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
     92        double form = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
    5993                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
    60         total += GAUSS_W[i]*fq*fq*sin_theta; 
     94        total_F1 += GAUSS_W[i]*form*sin_theta; 
     95        total_F2 += GAUSS_W[i]*form*form*sin_theta; 
    6196    } 
     97    // Correct for integration range 
     98    total_F1 *= uplim; 
     99    total_F2 *= uplim; 
    62100 
    63     // calculate value of integral to return 
    64     double answer = total*uplim; 
    65     return 1.0e-4*answer; 
     101    *F1 = 1.0e-2*total_F1; 
     102    *F2 = 1.0e-4*total_F2; 
    66103} 
    67104 
  • sasmodels/models/core_shell_bicelle.py

    r2d81cfe ree60aa7  
    154154source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
    155155          "core_shell_bicelle.c"] 
     156have_Fq = True 
     157effective_radius_type = [ 
     158    "equivalent sphere", "outer rim radius", 
     159    "half outer thickness", "half diagonal", 
     160    ] 
    156161 
    157162def random(): 
  • sasmodels/models/core_shell_bicelle_elliptical.c

    r108e70e rd42dd4a  
    1111 
    1212static double 
    13 Iq(double q, 
     13radius_from_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     14{ 
     15    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length); 
     16    return cbrt(volume_bicelle/M_4PI_3); 
     17} 
     18 
     19static double 
     20radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     21{ 
     22    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor); 
     23    const double radius_max_tot = radius_max + thick_rim; 
     24    const double length_tot = length + 2.0*thick_face; 
     25    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot); 
     26} 
     27 
     28static double 
     29effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     30{ 
     31    switch (mode) { 
     32    default: 
     33    case 1: // equivalent sphere 
     34        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
     35    case 2: // outer rim average radius 
     36        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
     37    case 3: // outer rim min radius 
     38        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     39    case 4: // outer max radius 
     40        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     41    case 5: // half outer thickness 
     42        return 0.5*length + thick_face; 
     43    case 6: // half diagonal 
     44        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
     45    } 
     46} 
     47 
     48static void 
     49Fq(double q, 
     50    double *F1, 
     51    double *F2, 
    1452    double r_minor, 
    1553    double x_core, 
     
    3674 
    3775    //initialize integral 
    38     double outer_sum = 0.0; 
     76    double outer_total_F1 = 0.0; 
     77    double outer_total_F2 = 0.0; 
    3978    for(int i=0;i<GAUSS_N;i++) { 
    4079        //setup inner integral over the ellipsoidal cross-section 
    41         //const double va = 0.0; 
    42         //const double vb = 1.0; 
    4380        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    4481        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 
     
    4885        const double si1 = sas_sinx_x(halfheight*qc); 
    4986        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
    50         double inner_sum=0.0; 
     87        double inner_total_F1 = 0; 
     88        double inner_total_F2 = 0; 
    5189        for(int j=0;j<GAUSS_N;j++) { 
    5290            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    53             // inner integral limits 
    54             //const double vaj=0.0; 
    55             //const double vbj=M_PI; 
    56             //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    57             const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    58             const double rr = sqrt(r2A - r2B*cos(phi)); 
     91            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     92            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
     93            const double rr = sqrt(r2A - r2B*cos(beta)); 
    5994            const double be1 = sas_2J1x_x(rr*qab); 
    6095            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
    61             const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
     96            const double f = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
    6297 
    63             inner_sum += GAUSS_W[j] * fq * fq; 
     98            inner_total_F1 += GAUSS_W[j] * f; 
     99            inner_total_F2 += GAUSS_W[j] * f * f; 
    64100        } 
    65101        //now calculate outer integral 
    66         outer_sum += GAUSS_W[i] * inner_sum; 
     102        outer_total_F1 += GAUSS_W[i] * inner_total_F1; 
     103        outer_total_F2 += GAUSS_W[i] * inner_total_F2; 
    67104    } 
     105    // now complete change of integration variables (1-0)/(1-(-1))= 0.5 
     106    outer_total_F1 *= 0.25; 
     107    outer_total_F2 *= 0.25; 
    68108 
    69     return outer_sum*2.5e-05; 
     109    //convert from [1e-12 A-1] to [cm-1] 
     110    *F1 = 1e-2*outer_total_F1; 
     111    *F2 = 1e-4*outer_total_F2; 
    70112} 
    71113 
  • sasmodels/models/core_shell_bicelle_elliptical.py

    rfc3ae1b r304c775  
    146146source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
    147147          "core_shell_bicelle_elliptical.c"] 
     148have_Fq = True 
     149effective_radius_type = [ 
     150    "equivalent sphere", "outer rim average radius", "outer rim min radius", 
     151    "outer max radius", "half outer thickness", "half diagonal", 
     152    ] 
    148153 
    149154def random(): 
     
    179184 
    180185tests = [ 
    181     [{'radius': 30.0, 'x_core': 3.0, 
    182       'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'ER', 1], 
    183     [{'radius': 30.0, 'x_core': 3.0, 
    184       'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'VR', 1], 
     186    #[{'radius': 30.0, 'x_core': 3.0, 
     187    #  'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'ER', 1], 
     188    #[{'radius': 30.0, 'x_core': 3.0, 
     189    #  'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'VR', 1], 
    185190 
    186191    [{'radius': 30.0, 'x_core': 3.0, 
  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r108e70e rd42dd4a  
    1212 
    1313static double 
    14 Iq(double q, 
     14radius_from_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     15{ 
     16    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length); 
     17    return cbrt(volume_bicelle/M_4PI_3); 
     18} 
     19 
     20static double 
     21radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     22{ 
     23    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor); 
     24    const double radius_max_tot = radius_max + thick_rim; 
     25    const double length_tot = length + 2.0*thick_face; 
     26    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot); 
     27} 
     28 
     29static double 
     30effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     31{ 
     32    switch (mode) { 
     33    default: 
     34    case 1: // equivalent sphere 
     35        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
     36    case 2: // outer rim average radius 
     37        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
     38    case 3: // outer rim min radius 
     39        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     40    case 4: // outer max radius 
     41        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     42    case 5: // half outer thickness 
     43        return 0.5*length + thick_face; 
     44    case 6: // half diagonal (this ignores the missing "corners", so may give unexpected answer if thick_face 
     45            //  or thick_rim is extremely large) 
     46        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
     47    } 
     48} 
     49 
     50static void 
     51Fq(double q, 
     52        double *F1, 
     53        double *F2, 
    1554        double r_minor, 
    1655        double x_core, 
     
    2463        double sigma) 
    2564{ 
    26     double si1,si2,be1,be2; 
    2765     // core_shell_bicelle_elliptical_belt, RKH 5th Oct 2017, core_shell_bicelle_elliptical 
    2866     // tested briefly against limiting cases of cylinder, hollow cylinder & elliptical cylinder models 
     
    3876    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
    3977    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
     78    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
     79    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*halfheight; 
     80    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
    4081    // dr1,2,3 are now for Vcore, Vcore+rim, Vcore+face, 
    41     const double dr1 = (-rhor - rhoh + rhoc + rhosolv) *M_PI*r_minor*r_major* 
    42           2.0*halfheight; 
    43     const double dr2 = (rhor-rhosolv) *M_PI*(r_minor+thick_rim)*( 
    44          r_major+thick_rim)* 2.0*halfheight; 
    45     const double dr3 = (rhoh-rhosolv) *M_PI*r_minor*r_major* 
    46          2.0*(halfheight+thick_face); 
     82    const double dr1 = vol1*(-rhor - rhoh + rhoc + rhosolv); 
     83    const double dr2 = vol2*(rhor-rhosolv); 
     84    const double dr3 = vol3*(rhoh-rhosolv); 
     85 
    4786    //initialize integral 
    48     double outer_sum = 0.0; 
     87    double outer_total_F1 = 0.0; 
     88    double outer_total_F2 = 0.0; 
    4989    for(int i=0;i<GAUSS_N;i++) { 
    5090        //setup inner integral over the ellipsoidal cross-section 
    5191        // since we generate these lots of times, why not store them somewhere? 
    52         //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    53         const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0; 
    54         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    55         double inner_sum=0; 
    56         double sinarg1 = q*halfheight*cos_alpha; 
    57         double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 
    58         si1 = sas_sinx_x(sinarg1); 
    59         si2 = sas_sinx_x(sinarg2); 
     92        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     93        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 
     94        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     95        const double qab = q*sin_theta; 
     96        const double qc = q*cos_theta; 
     97        const double si1 = sas_sinx_x(halfheight*qc); 
     98        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
     99        double inner_total_F1 = 0; 
     100        double inner_total_F2 = 0; 
    60101        for(int j=0;j<GAUSS_N;j++) { 
    61102            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
     
    63104            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    64105            const double rr = sqrt(r2A - r2B*cos(beta)); 
    65             double besarg1 = q*rr*sin_alpha; 
    66             double besarg2 = q*(rr+thick_rim)*sin_alpha; 
    67             be1 = sas_2J1x_x(besarg1); 
    68             be2 = sas_2J1x_x(besarg2); 
    69             inner_sum += GAUSS_W[j] *square(dr1*si1*be1 + 
    70                                               dr2*si1*be2 + 
    71                                               dr3*si2*be1); 
     106            const double be1 = sas_2J1x_x(rr*qab); 
     107            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
     108            const double f = dr1*si1*be1 + dr2*si1*be2 + dr3*si2*be1; 
     109 
     110            inner_total_F1 += GAUSS_W[j] * f; 
     111            inner_total_F2 += GAUSS_W[j] * f * f; 
    72112        } 
    73113        //now calculate outer integral 
    74         outer_sum += GAUSS_W[i] * inner_sum; 
     114        outer_total_F1 += GAUSS_W[i] * inner_total_F1; 
     115        outer_total_F2 += GAUSS_W[i] * inner_total_F2; 
    75116    } 
     117    // now complete change of integration variables (1-0)/(1-(-1))= 0.5 
     118    outer_total_F1 *= 0.25; 
     119    outer_total_F2 *= 0.25; 
    76120 
    77     return outer_sum*2.5e-05*exp(-0.5*square(q*sigma)); 
     121    //convert from [1e-12 A-1] to [cm-1] 
     122    *F1 = 1e-2*outer_total_F1*exp(-0.25*square(q*sigma)); 
     123    *F2 = 1e-4*outer_total_F2*exp(-0.5*square(q*sigma)); 
    78124} 
    79125 
     
    111157    const double si1 = sas_sinx_x( halfheight*qc ); 
    112158    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc ); 
    113     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1); 
    114     return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma)); 
     159    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1; 
     160    const double atten = exp(-0.5*(qa*qa + qb*qb + qc*qc)*(sigma*sigma)); 
     161    return 1.0e-4 * fq*fq * atten; 
    115162} 
  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    rfc3ae1b r304c775  
    159159source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
    160160          "core_shell_bicelle_elliptical_belt_rough.c"] 
     161have_Fq = True 
     162effective_radius_type = [ 
     163    "equivalent sphere", "outer rim average radius", "outer rim min radius", 
     164    "outer max radius", "half outer thickness", "half diagonal", 
     165    ] 
    161166 
    162167demo = dict(scale=1, background=0, 
     
    181186 
    182187tests = [ 
    183     [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'ER', 1], 
    184     [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'VR', 1], 
     188    #[{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'ER', 1], 
     189    #[{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'VR', 1], 
    185190 
    186191    [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0, 
  • sasmodels/models/core_shell_cylinder.c

    rd86f0fc rd42dd4a  
    1414 
    1515static double 
    16 Iq(double q, 
     16radius_from_volume(double radius, double thickness, double length) 
     17{ 
     18    const double volume_outer_cyl = form_volume(radius,thickness,length); 
     19    return cbrt(volume_outer_cyl/M_4PI_3); 
     20} 
     21 
     22static double 
     23radius_from_diagonal(double radius, double thickness, double length) 
     24{ 
     25    const double radius_outer = radius + thickness; 
     26    const double length_outer = length + 2.0*thickness; 
     27    return sqrt(radius_outer*radius_outer + 0.25*length_outer*length_outer); 
     28} 
     29 
     30static double 
     31effective_radius(int mode, double radius, double thickness, double length) 
     32{ 
     33    switch (mode) { 
     34    default: 
     35    case 1: // equivalent sphere 
     36        return radius_from_volume(radius, thickness, length); 
     37    case 2: // outer radius 
     38        return radius + thickness; 
     39    case 3: // half outer length 
     40        return 0.5*length + thickness; 
     41    case 4: // half min outer length 
     42        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness); 
     43    case 5: // half max outer length 
     44        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness); 
     45    case 6: // half outer diagonal 
     46        return radius_from_diagonal(radius,thickness,length); 
     47    } 
     48} 
     49 
     50static void 
     51Fq(double q, 
     52    double *F1, 
     53    double *F2, 
    1754    double core_sld, 
    1855    double shell_sld, 
     
    2966    const double shell_h = (0.5*length + thickness); 
    3067    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    31     double total = 0.0; 
     68    double total_F1 = 0.0; 
     69    double total_F2 = 0.0; 
    3270    for (int i=0; i<GAUSS_N ;i++) { 
    3371        // translate a point in [-1,1] to a point in [0, pi/2] 
     
    4078        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
    4179            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    42         total += GAUSS_W[i] * fq * fq * sin_theta; 
     80        total_F1 += GAUSS_W[i] * fq * sin_theta; 
     81        total_F2 += GAUSS_W[i] * fq * fq * sin_theta; 
    4382    } 
    4483    // translate dx in [-1,1] to dx in [lower,upper] 
    4584    //const double form = (upper-lower)/2.0*total; 
    46     return 1.0e-4 * total * M_PI_4; 
     85    *F1 = 1.0e-2 * total_F1 * M_PI_4; 
     86    *F2 = 1.0e-4 * total_F2 * M_PI_4; 
    4787} 
    48  
    4988 
    5089static double 
  • sasmodels/models/core_shell_cylinder.py

    re31b19a ree60aa7  
    131131 
    132132source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] 
    133  
    134 def ER(radius, thickness, length): 
    135     """ 
    136     Returns the effective radius used in the S*P calculation 
    137     """ 
    138     radius = radius + thickness 
    139     length = length + 2 * thickness 
    140     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    141     return 0.5 * (ddd) ** (1. / 3.) 
     133have_Fq = True 
     134effective_radius_type = [ 
     135    "equivalent sphere", "outer radius", "half outer length", 
     136    "half min outer dimension", "half max outer dimension", 
     137    "half outer diagonal", 
     138    ] 
    142139 
    143140def random(): 
  • sasmodels/models/core_shell_ellipsoid.c

    r108e70e rd42dd4a  
    3939 
    4040static double 
    41 Iq(double q, 
     41radius_from_volume(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     42{ 
     43    const double volume_ellipsoid = form_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     44    return cbrt(volume_ellipsoid/M_4PI_3); 
     45} 
     46 
     47static double 
     48radius_from_curvature(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     49{ 
     50    // Trivial cases 
     51    if (1.0 == x_core && 1.0 == x_polar_shell) return radius_equat_core + thick_shell; 
     52    if ((radius_equat_core + thick_shell)*(radius_equat_core*x_core + thick_shell*x_polar_shell) == 0.)  return 0.; 
     53 
     54    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     55    const double radius_equat_tot = radius_equat_core + thick_shell; 
     56    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 
     57    const double ratio = (radius_polar_tot < radius_equat_tot 
     58                          ? radius_polar_tot / radius_equat_tot 
     59                          : radius_equat_tot / radius_polar_tot); 
     60    const double e1 = sqrt(1.0 - ratio*ratio); 
     61    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     62    const double bL = (1.0 + e1) / (1.0 - e1); 
     63    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     64    const double delta = 0.75 * b1 * b2; 
     65    const double ddd = 2.0 * (delta + 1.0) * radius_polar_tot * radius_equat_tot * radius_equat_tot; 
     66    return 0.5 * cbrt(ddd); 
     67} 
     68 
     69static double 
     70effective_radius(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     71{ 
     72    const double radius_equat_tot = radius_equat_core + thick_shell; 
     73    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 
     74 
     75    switch (mode) { 
     76    default: 
     77    case 1: // equivalent sphere 
     78        return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     79    case 2: // average outer curvature 
     80        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     81    case 3: // min outer radius 
     82        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
     83    case 4: // max outer radius 
     84        return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
     85    } 
     86} 
     87 
     88static void 
     89Fq(double q, 
     90    double *F1, 
     91    double *F2, 
    4292    double radius_equat_core, 
    4393    double x_core, 
     
    58108    const double m = 0.5; 
    59109    const double b = 0.5; 
    60     double total = 0.0;     //initialize intergral 
     110    double total_F1 = 0.0;     //initialize intergral 
     111    double total_F2 = 0.0;     //initialize intergral 
    61112    for(int i=0;i<GAUSS_N;i++) { 
    62113        const double cos_theta = GAUSS_Z[i]*m + b; 
     
    66117            equat_shell, polar_shell, 
    67118            sld_core_shell, sld_shell_solvent); 
    68         total += GAUSS_W[i] * fq * fq; 
     119        total_F1 += GAUSS_W[i] * fq; 
     120        total_F2 += GAUSS_W[i] * fq * fq; 
    69121    } 
    70     total *= m; 
     122    total_F1 *= m; 
     123    total_F2 *= m; 
    71124 
    72125    // convert to [cm-1] 
    73     return 1.0e-4 * total; 
     126    *F1 = 1.0e-2 * total_F1; 
     127    *F2 = 1.0e-4 * total_F2; 
    74128} 
     129 
    75130 
    76131static double 
  • sasmodels/models/core_shell_ellipsoid.py

    r2d81cfe ree60aa7  
    145145 
    146146source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    147  
    148 def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 
    149     """ 
    150         Returns the effective radius used in the S*P calculation 
    151     """ 
    152     from .ellipsoid import ER as ellipsoid_ER 
    153     polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell 
    154     equat_outer = radius_equat_core + thick_shell 
    155     return ellipsoid_ER(polar_outer, equat_outer) 
     147have_Fq = True 
     148effective_radius_type = [ 
     149    "equivalent sphere", "average outer curvature", 
     150    "min outer radius", "max outer radius", 
     151    ] 
    156152 
    157153def random(): 
  • sasmodels/models/core_shell_parallelepiped.c

    rdbf1a60 rd42dd4a  
    2828 
    2929static double 
    30 Iq(double q, 
     30radius_from_volume(double length_a, double length_b, double length_c, 
     31                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     32{ 
     33    const double volume = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     34    return cbrt(volume/M_4PI_3); 
     35} 
     36 
     37static double 
     38radius_from_crosssection(double length_a, double length_b, double thick_rim_a, double thick_rim_b) 
     39{ 
     40    const double area_xsec_paral = length_a*length_b + 2.0*thick_rim_a*length_b + 2.0*thick_rim_b*length_a; 
     41    return sqrt(area_xsec_paral/M_PI); 
     42} 
     43 
     44static double 
     45effective_radius(int mode, double length_a, double length_b, double length_c, 
     46                 double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     47{ 
     48    switch (mode) { 
     49    default: 
     50    case 1: // equivalent sphere 
     51        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     52    case 2: // half outer length a 
     53        return 0.5 * length_a + thick_rim_a; 
     54    case 3: // half outer length b 
     55        return 0.5 * length_b + thick_rim_b; 
     56    case 4: // half outer length c 
     57        return 0.5 * length_c + thick_rim_c; 
     58    case 5: // equivalent circular cross-section 
     59        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b); 
     60    case 6: // half outer ab diagonal 
     61        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b)); 
     62    case 7: // half outer diagonal 
     63        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b) + square(length_c+ 2.0*thick_rim_c)); 
     64    } 
     65} 
     66 
     67static void 
     68Fq(double q, 
     69    double *F1, 
     70    double *F2, 
    3171    double core_sld, 
    3272    double arim_sld, 
     
    60100    // outer integral (with gauss points), integration limits = 0, 1 
    61101    // substitute d_cos_alpha for sin_alpha d_alpha 
    62     double outer_sum = 0; //initialize integral 
     102    double outer_sum_F1 = 0; //initialize integral 
     103    double outer_sum_F2 = 0; //initialize integral 
    63104    for( int i=0; i<GAUSS_N; i++) { 
    64105        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
     
    69110        // inner integral (with gauss points), integration limits = 0, 1 
    70111        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta) 
    71         double inner_sum = 0.0; 
     112        double inner_sum_F1 = 0.0; 
     113        double inner_sum_F2 = 0.0; 
    72114        for(int j=0; j<GAUSS_N; j++) { 
    73115            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     
    91133#endif 
    92134 
    93             inner_sum += GAUSS_W[j] * f * f; 
     135            inner_sum_F1 += GAUSS_W[j] * f; 
     136            inner_sum_F2 += GAUSS_W[j] * f * f; 
    94137        } 
    95138        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    96         inner_sum *= 0.5; 
    97         // now sum up the outer integral 
    98         outer_sum += GAUSS_W[i] * inner_sum; 
     139        // and sum up the outer integral 
     140        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * 0.5; 
     141        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * 0.5; 
    99142    } 
    100143    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    101     outer_sum *= 0.5; 
     144    outer_sum_F1 *= 0.5; 
     145    outer_sum_F2 *= 0.5; 
    102146 
    103147    //convert from [1e-12 A-1] to [cm-1] 
    104     return 1.0e-4 * outer_sum; 
     148    *F1 = 1.0e-2 * outer_sum_F1; 
     149    *F2 = 1.0e-4 * outer_sum_F2; 
    105150} 
    106151 
  • sasmodels/models/core_shell_parallelepiped.py

    rf89ec96 ree60aa7  
    9898is the scattering length of the solvent. 
    9999 
    100 .. note::  
     100.. note:: 
    101101 
    102102   the code actually implements two substitutions: $d(cos\alpha)$ is 
     
    226226 
    227227source = ["lib/gauss76.c", "core_shell_parallelepiped.c"] 
    228  
    229  
    230 def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c): 
    231     """ 
    232         Return equivalent radius (ER) 
    233     """ 
    234     from .parallelepiped import ER as ER_p 
    235  
    236     a = length_a + 2*thick_rim_a 
    237     b = length_b + 2*thick_rim_b 
    238     c = length_c + 2*thick_rim_c 
    239     return ER_p(a, b, c) 
    240  
    241 # VR defaults to 1.0 
     228have_Fq = True 
     229effective_radius_type = [ 
     230    "equivalent sphere", 
     231    "half outer length_a", "half outer length_b", "half outer length_c", 
     232    "equivalent circular cross-section", 
     233    "half outer ab diagonal", "half outer diagonal", 
     234    ] 
    242235 
    243236def random(): 
  • sasmodels/models/core_shell_sphere.c

    r3a48772 rd42dd4a  
    1 double form_volume(double radius, double thickness); 
    2 double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld); 
     1static double 
     2form_volume(double radius, double thickness) 
     3{ 
     4    return M_4PI_3 * cube(radius + thickness); 
     5} 
    36 
    4 double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) { 
     7static double 
     8effective_radius(int mode, double radius, double thickness) 
     9{ 
     10    switch (mode) { 
     11    default: 
     12    case 1: // outer radius 
     13        return radius + thickness; 
     14    case 2: // core radius 
     15        return radius; 
     16    } 
     17} 
    518 
    6  
    7     double intensity = core_shell_kernel(q, 
     19static void 
     20Fq(double q, double *F1, double *F2, double radius, 
     21   double thickness, double core_sld, double shell_sld, double solvent_sld) { 
     22    double form = core_shell_fq(q, 
    823                              radius, 
    924                              thickness, 
     
    1126                              shell_sld, 
    1227                              solvent_sld); 
    13     return intensity; 
     28    *F1 = 1.0e-2*form; 
     29    *F2 = 1.0e-4*form*form; 
    1430} 
    15  
    16 double form_volume(double radius, double thickness) 
    17 { 
    18     return M_4PI_3 * cube(radius + thickness); 
    19 } 
  • sasmodels/models/core_shell_sphere.py

    rda1c8d1 r304c775  
    5858title = "Form factor for a monodisperse spherical particle with particle with a core-shell structure." 
    5959description = """ 
    60     F^2(q) = 3/V_s [V_c (sld_core-sld_shell) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
    61                    + V_s (sld_shell-sld_solvent) (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3] 
     60    F(q) = [V_c (sld_core-sld_shell) 3 (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
     61            + V_s (sld_shell-sld_solvent) 3 (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3] 
    6262 
    6363            V_s: Volume of the sphere shell 
     
    7777 
    7878source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] 
     79have_Fq = True 
     80effective_radius_type = ["outer radius", "core radius"] 
    7981 
    8082demo = dict(scale=1, background=0, radius=60, thickness=10, 
    8183            sld_core=1.0, sld_shell=2.0, sld_solvent=0.0) 
    82  
    83 def ER(radius, thickness): 
    84     """ 
    85         Equivalent radius 
    86         @param radius: core radius 
    87         @param thickness: shell thickness 
    88     """ 
    89     return radius + thickness 
    9084 
    9185def random(): 
     
    10296 
    10397tests = [ 
    104     [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
     98    [{'radius': 20.0, 'thickness': 10.0}, 0.1, None, None, 30.0, 4.*pi/3*30**3, 1.0], 
    10599 
    106100    # The SasView test result was 0.00169, with a background of 0.001 
  • sasmodels/models/cylinder.c

    r108e70e rd42dd4a  
    88 
    99static double 
    10 fq(double qab, double qc, double radius, double length) 
     10_fq(double qab, double qc, double radius, double length) 
    1111{ 
    1212    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
     
    1414 
    1515static double 
    16 orient_avg_1D(double q, double radius, double length) 
     16radius_from_volume(double radius, double length) 
     17{ 
     18    return cbrt(M_PI*radius*radius*length/M_4PI_3); 
     19} 
     20 
     21static double 
     22radius_from_diagonal(double radius, double length) 
     23{ 
     24    return sqrt(radius*radius + 0.25*length*length); 
     25} 
     26 
     27static double 
     28effective_radius(int mode, double radius, double length) 
     29{ 
     30    switch (mode) { 
     31    default: 
     32    case 1: 
     33        return radius_from_volume(radius, length); 
     34    case 2: 
     35        return radius; 
     36    case 3: 
     37        return 0.5*length; 
     38    case 4: 
     39        return (radius < 0.5*length ? radius : 0.5*length); 
     40    case 5: 
     41        return (radius > 0.5*length ? radius : 0.5*length); 
     42    case 6: 
     43        return radius_from_diagonal(radius,length); 
     44    } 
     45} 
     46 
     47static void 
     48Fq(double q, 
     49    double *F1, 
     50    double *F2, 
     51    double sld, 
     52    double solvent_sld, 
     53    double radius, 
     54    double length) 
    1755{ 
    1856    // translate a point in [-1,1] to a point in [0, pi/2] 
     
    2058    const double zb = M_PI_4; 
    2159 
    22     double total = 0.0; 
     60    double total_F1 = 0.0; 
     61    double total_F2 = 0.0; 
    2362    for (int i=0; i<GAUSS_N ;i++) { 
    2463        const double theta = GAUSS_Z[i]*zm + zb; 
     
    2665        // theta (theta,phi) the projection of the cylinder on the detector plane 
    2766        SINCOS(theta , sin_theta, cos_theta); 
    28         const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
    29         total += GAUSS_W[i] * form * form * sin_theta; 
     67        const double form = _fq(q*sin_theta, q*cos_theta, radius, length); 
     68        total_F1 += GAUSS_W[i] * form * sin_theta; 
     69        total_F2 += GAUSS_W[i] * form * form * sin_theta; 
    3070    } 
    3171    // translate dx in [-1,1] to dx in [lower,upper] 
    32     return total*zm; 
     72    total_F1 *= zm; 
     73    total_F2 *= zm; 
     74    const double s = (sld - solvent_sld) * form_volume(radius, length); 
     75    *F1 = 1e-2 * s * total_F1; 
     76    *F2 = 1e-4 * s * s * total_F2; 
    3377} 
    3478 
    35 static double 
    36 Iq(double q, 
    37     double sld, 
    38     double solvent_sld, 
    39     double radius, 
    40     double length) 
    41 { 
    42     const double s = (sld - solvent_sld) * form_volume(radius, length); 
    43     return 1.0e-4 * s * s * orient_avg_1D(q, radius, length); 
    44 } 
     79 
    4580 
    4681static double 
     
    5186    double length) 
    5287{ 
     88    const double form = _fq(qab, qc, radius, length); 
    5389    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    54     const double form = fq(qab, qc, radius, length); 
    5590    return 1.0e-4 * square(s * form); 
    5691} 
     92 
  • sasmodels/models/cylinder.py

    r2d81cfe r304c775  
    138138 
    139139source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 
    140  
    141 def ER(radius, length): 
    142     """ 
    143         Return equivalent radius (ER) 
    144     """ 
    145     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    146     return 0.5 * (ddd) ** (1. / 3.) 
     140have_Fq = True 
     141effective_radius_type = [ 
     142    "equivalent sphere", "radius", 
     143    "half length", "half min dimension", "half max dimension", "half diagonal", 
     144    ] 
    147145 
    148146def random(): 
     
    169167            phi_pd=10, phi_pd_n=5) 
    170168 
     169# pylint: disable=bad-whitespace, line-too-long 
    171170qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    172171# After redefinition of angles, find new tests values.  Was 10 10 in old coords 
     
    182181] 
    183182del qx, qy  # not necessary to delete, but cleaner 
     183 
     184# Default radius and length 
     185radius, length = parameters[2][2], parameters[3][2] 
     186tests.extend([ 
     187    ({'radius_effective_type': 0}, 0.1, None, None, 0., pi*radius**2*length, 1.0), 
     188    ({'radius_effective_type': 1}, 0.1, None, None, (0.75*radius**2*length)**(1./3.), None, None), 
     189    ({'radius_effective_type': 2}, 0.1, None, None, radius, None, None), 
     190    ({'radius_effective_type': 3}, 0.1, None, None, length/2., None, None), 
     191    ({'radius_effective_type': 4}, 0.1, None, None, min(radius, length/2.), None, None), 
     192    ({'radius_effective_type': 5}, 0.1, None, None, max(radius, length/2.), None, None), 
     193    ({'radius_effective_type': 6}, 0.1, None, None, np.sqrt(4*radius**2 + length**2)/2., None, None), 
     194]) 
     195del radius, length 
     196# pylint: enable=bad-whitespace, line-too-long 
     197 
    184198# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc 
  • sasmodels/models/dab.py

    r2d81cfe r304c775  
    7474    return pars 
    7575 
    76 # ER defaults to 1.0 
    77  
    78 # VR defaults to 1.0 
    79  
    8076demo = dict(scale=1, background=0, cor_length=50) 
  • sasmodels/models/ellipsoid.c

    r108e70e rd42dd4a  
    55} 
    66 
    7 static  double 
    8 Iq(double q, 
     7static double 
     8radius_from_volume(double radius_polar, double radius_equatorial) 
     9{ 
     10    return cbrt(radius_polar*radius_equatorial*radius_equatorial); 
     11} 
     12 
     13static double 
     14radius_from_curvature(double radius_polar, double radius_equatorial) 
     15{ 
     16    // Trivial cases 
     17    if (radius_polar == radius_equatorial) return radius_polar; 
     18    if (radius_polar * radius_equatorial == 0.)  return 0.; 
     19 
     20    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     21    const double ratio = (radius_polar < radius_equatorial 
     22                          ? radius_polar / radius_equatorial 
     23                          : radius_equatorial / radius_polar); 
     24    const double e1 = sqrt(1.0 - ratio*ratio); 
     25    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     26    const double bL = (1.0 + e1) / (1.0 - e1); 
     27    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     28    const double delta = 0.75 * b1 * b2; 
     29    const double ddd = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial * radius_equatorial; 
     30    return 0.5 * cbrt(ddd); 
     31} 
     32 
     33static double 
     34effective_radius(int mode, double radius_polar, double radius_equatorial) 
     35{ 
     36    switch (mode) { 
     37    default: 
     38    case 1: // equivalent sphere 
     39        return radius_from_volume(radius_polar, radius_equatorial); 
     40    case 2: // average curvature 
     41        return radius_from_curvature(radius_polar, radius_equatorial); 
     42    case 3: // min radius 
     43        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial); 
     44    case 4: // max radius 
     45        return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial); 
     46    } 
     47} 
     48 
     49 
     50static void 
     51Fq(double q, 
     52    double *F1, 
     53    double *F2, 
    954    double sld, 
    1055    double sld_solvent, 
     
    2065    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
    2166    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
    22  
    2367    // translate a point in [-1,1] to a point in [0, 1] 
    2468    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    2569    const double zm = 0.5; 
    2670    const double zb = 0.5; 
    27     double total = 0.0; 
     71    double total_F2 = 0.0; 
     72    double total_F1 = 0.0; 
    2873    for (int i=0;i<GAUSS_N;i++) { 
    2974        const double u = GAUSS_Z[i]*zm + zb; 
    3075        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
    3176        const double f = sas_3j1x_x(q*r); 
    32         total += GAUSS_W[i] * f * f; 
     77        total_F2 += GAUSS_W[i] * f * f; 
     78        total_F1 += GAUSS_W[i] * f; 
    3379    } 
    3480    // translate dx in [-1,1] to dx in [lower,upper] 
    35     const double form = total*zm; 
     81    total_F1 *= zm; 
     82    total_F2 *= zm; 
    3683    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    37     return 1.0e-4 * s * s * form; 
     84    *F1 = 1e-2 * s * total_F1; 
     85    *F2 = 1e-4 * s * s * total_F2; 
    3886} 
    3987 
  • sasmodels/models/ellipsoid.py

    r0168844 ree60aa7  
    166166             ] 
    167167 
     168 
    168169source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
    169  
    170 def ER(radius_polar, radius_equatorial): 
    171     # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
    172     ee = np.empty_like(radius_polar) 
    173     idx = radius_polar > radius_equatorial 
    174     ee[idx] = (radius_polar[idx] ** 2 - radius_equatorial[idx] ** 2) / radius_polar[idx] ** 2 
    175     idx = radius_polar < radius_equatorial 
    176     ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2 
    177     valid = (radius_polar * radius_equatorial != 0) & (radius_polar != radius_equatorial) 
    178     bd = 1.0 - ee[valid] 
    179     e1 = np.sqrt(ee[valid]) 
    180     b1 = 1.0 + np.arcsin(e1) / (e1 * np.sqrt(bd)) 
    181     bL = (1.0 + e1) / (1.0 - e1) 
    182     b2 = 1.0 + bd / 2 / e1 * np.log(bL) 
    183     delta = 0.75 * b1 * b2 
    184     ddd = 2.0 * (delta + 1.0) * (radius_polar * radius_equatorial**2)[valid] 
    185  
    186     r = np.zeros_like(radius_polar) 
    187     r[valid] = 0.5 * cbrt(ddd) 
    188     idx = radius_polar == radius_equatorial 
    189     r[idx] = radius_polar[idx] 
    190     return r 
     170have_Fq = True 
     171effective_radius_type = [ 
     172    "equivalent sphere", "average curvature", "min radius", "max radius", 
     173    ] 
    191174 
    192175def random(): 
  • sasmodels/models/elliptical_cylinder.c

    r108e70e rd42dd4a  
    66 
    77static double 
    8 Iq(double q, double radius_minor, double r_ratio, double length, 
     8radius_from_volume(double radius_minor, double r_ratio, double length) 
     9{ 
     10    const double volume_ellcyl = form_volume(radius_minor,r_ratio,length); 
     11    return cbrt(volume_ellcyl/M_4PI_3); 
     12} 
     13 
     14static double 
     15radius_from_min_dimension(double radius_minor, double r_ratio, double hlength) 
     16{ 
     17    const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
     18    return (rad_min < hlength ? rad_min : hlength); 
     19} 
     20 
     21static double 
     22radius_from_max_dimension(double radius_minor, double r_ratio, double hlength) 
     23{ 
     24    const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
     25    return (rad_max > hlength ? rad_max : hlength); 
     26} 
     27 
     28static double 
     29radius_from_diagonal(double radius_minor, double r_ratio, double length) 
     30{ 
     31    const double radius_max = (r_ratio > 1.0 ? radius_minor*r_ratio : radius_minor); 
     32    return sqrt(radius_max*radius_max + 0.25*length*length); 
     33} 
     34 
     35static double 
     36effective_radius(int mode, double radius_minor, double r_ratio, double length) 
     37{ 
     38    switch (mode) { 
     39    default: 
     40    case 1: // equivalent sphere 
     41        return radius_from_volume(radius_minor, r_ratio, length); 
     42    case 2: // average radius 
     43        return 0.5*radius_minor*(1.0 + r_ratio); 
     44    case 3: // min radius 
     45        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
     46    case 4: // max radius 
     47        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
     48    case 5: // equivalent circular cross-section 
     49        return sqrt(radius_minor*radius_minor*r_ratio); 
     50    case 6: // half length 
     51        return 0.5*length; 
     52    case 7: // half min dimension 
     53        return radius_from_min_dimension(radius_minor,r_ratio,0.5*length); 
     54    case 8: // half max dimension 
     55        return radius_from_max_dimension(radius_minor,r_ratio,0.5*length); 
     56    case 9: // half diagonal 
     57        return radius_from_diagonal(radius_minor,r_ratio,length); 
     58    } 
     59} 
     60 
     61static void 
     62Fq(double q, double *F1, double *F2, double radius_minor, double r_ratio, double length, 
    963   double sld, double solvent_sld) 
    1064{ 
     
    2175 
    2276    //initialize integral 
    23     double outer_sum = 0.0; 
     77    double outer_sum_F1 = 0.0; 
     78    double outer_sum_F2 = 0.0; 
    2479    for(int i=0;i<GAUSS_N;i++) { 
    2580        //setup inner integral over the ellipsoidal cross-section 
     
    2782        const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    2883        //const double arg = radius_minor*sin_val; 
    29         double inner_sum=0; 
     84        double inner_sum_F1 = 0.0; 
     85        double inner_sum_F2 = 0.0; 
    3086        for(int j=0;j<GAUSS_N;j++) { 
    3187            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    3288            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3389            const double be = sas_2J1x_x(q*r); 
    34             inner_sum += GAUSS_W[j] * be * be; 
     90            inner_sum_F1 += GAUSS_W[j] * be; 
     91            inner_sum_F2 += GAUSS_W[j] * be * be; 
    3592        } 
    3693        //now calculate the value of the inner integral 
    37         inner_sum *= 0.5*(vbj-vaj); 
     94        inner_sum_F1 *= 0.5*(vbj-vaj); 
     95        inner_sum_F2 *= 0.5*(vbj-vaj); 
    3896 
    3997        //now calculate outer integral 
    4098        const double si = sas_sinx_x(q*0.5*length*cos_val); 
    41         outer_sum += GAUSS_W[i] * inner_sum * si * si; 
     99        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * si; 
     100        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * si * si; 
    42101    } 
    43     outer_sum *= 0.5*(vb-va); 
    44  
    45     //divide integral by Pi 
    46     const double form = outer_sum/M_PI; 
     102    // correct limits and divide integral by pi 
     103    outer_sum_F1 *= 0.5*(vb-va)/M_PI; 
     104    outer_sum_F2 *= 0.5*(vb-va)/M_PI; 
    47105 
    48106    // scale by contrast and volume, and convert to to 1/cm units 
    49     const double vol = form_volume(radius_minor, r_ratio, length); 
    50     const double delrho = sld - solvent_sld; 
    51     return 1.0e-4*square(delrho*vol)*form; 
     107    const double volume = form_volume(radius_minor, r_ratio, length); 
     108    const double contrast = sld - solvent_sld; 
     109    const double s = contrast*volume; 
     110    *F1 = 1.0e-2*s*outer_sum_F1; 
     111    *F2 = 1.0e-4*s*s*outer_sum_F2; 
    52112} 
    53113 
     
    63123    const double be = sas_2J1x_x(qr); 
    64124    const double si = sas_sinx_x(qc*0.5*length); 
    65     const double Aq = be * si; 
    66     const double delrho = sld - solvent_sld; 
    67     const double vol = form_volume(radius_minor, r_ratio, length); 
    68     return 1.0e-4 * square(delrho * vol * Aq); 
     125    const double fq = be * si; 
     126    const double contrast = sld - solvent_sld; 
     127    const double volume = form_volume(radius_minor, r_ratio, length); 
     128    return 1.0e-4 * square(contrast * volume * fq); 
    69129} 
  • sasmodels/models/elliptical_cylinder.py

    ra261a83 ree60aa7  
    122122 
    123123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
     124have_Fq = True 
     125effective_radius_type = [ 
     126    "equivalent sphere", "average radius", "min radius", "max radius", 
     127    "equivalent circular cross-section", 
     128    "half length", "half min dimension", "half max dimension", "half diagonal", 
     129    ] 
    124130 
    125131demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 
    126132            sld=4.0, sld_solvent=1.0, theta=10.0, phi=20, psi=30, 
    127133            theta_pd=10, phi_pd=2, psi_pd=3) 
    128  
    129 def ER(radius_minor, axis_ratio, length): 
    130     """ 
    131         Equivalent radius 
    132         @param radius_minor: Ellipse minor radius 
    133         @param axis_ratio: Ratio of major radius over minor radius 
    134         @param length: Length of the cylinder 
    135     """ 
    136     radius = sqrt(radius_minor * radius_minor * axis_ratio) 
    137     ddd = 0.75 * radius * (2 * radius * length 
    138                            + (length + radius) * (length + pi * radius)) 
    139     return 0.5 * (ddd) ** (1. / 3.) 
    140134 
    141135def random(): 
     
    161155 
    162156tests = [ 
    163     [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], 
    164     [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], 
     157#    [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], 
     158#    [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], 
    165159 
    166160    # The SasView test result was 0.00169, with a background of 0.001 
  • sasmodels/models/fcc_paracrystal.c

    r108e70e r71b751d  
    7979} 
    8080 
    81  
    8281static double Iqabc(double qa, double qb, double qc, 
    8382    double dnn, double d_factor, double radius, 
  • sasmodels/models/flexible_cylinder_elliptical.c

    r74768cb r71b751d  
    7272    return result; 
    7373} 
    74  
  • sasmodels/models/fractal.c

    r4788822 r71b751d  
    1919 
    2020    //calculate P(q) for the spherical subunits 
    21     const double pq = square(form_volume(radius) * (sld_block-sld_solvent) 
    22                       *sas_3j1x_x(q*radius)); 
     21    const double fq = form_volume(radius) * (sld_block-sld_solvent) 
     22                      *sas_3j1x_x(q*radius); 
    2323 
    2424    // scale to units cm-1 sr-1 (assuming data on absolute scale) 
     
    2727    //    combined: 1e-4 * I(q) 
    2828 
    29     return 1.e-4 * volfraction * sq * pq; 
     29    return 1.e-4 * volfraction * sq * fq * fq; 
    3030} 
    31  
  • sasmodels/models/fractal_core_shell.c

    rbdd08df rbad3093  
    1616   double cor_length) 
    1717{ 
    18     //The radius for the building block of the core shell particle that is  
     18    //The radius for the building block of the core shell particle that is 
    1919    //needed by the Teixeira fractal S(q) is the radius of the whole particle. 
    2020    const double cs_radius = radius + thickness; 
    2121    const double sq = fractal_sq(q, cs_radius, fractal_dim, cor_length); 
    22     const double pq = core_shell_kernel(q, radius, thickness, 
    23                                         core_sld, shell_sld, solvent_sld); 
     22    const double fq = core_shell_fq(q, radius, thickness, 
     23                                    core_sld, shell_sld, solvent_sld); 
    2424 
    25     // Note: core_shell_kernel already performs the 1e-4 unit conversion 
    26     return volfraction * sq * pq; 
     25    return 1.0e-4 * volfraction * sq * fq * fq; 
    2726} 
    28  
  • sasmodels/models/fractal_core_shell.py

    reb3eb38 r2cc8aa2  
    134134    return radius + thickness 
    135135 
    136 tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
    137  
     136#tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
     137tests = [ 
    138138#         # At some point the SasView 3.x test result was deemed incorrect. The 
    139139          #following tests were verified against NIST IGOR macros ver 7.850. 
  • sasmodels/models/fuzzy_sphere.py

    r2d81cfe ree60aa7  
    7878              ["sld_solvent", "1e-6/Ang^2",  3, [-inf, inf], "sld",    "Solvent scattering length density"], 
    7979              ["radius",      "Ang",        60, [0, inf],    "volume", "Sphere radius"], 
    80               ["fuzziness",   "Ang",        10, [0, inf],    "",       "std deviation of Gaussian convolution for interface (must be << radius)"], 
     80              ["fuzziness",   "Ang",        10, [0, inf],    "volume",       "std deviation of Gaussian convolution for interface (must be << radius)"], 
    8181             ] 
    8282# pylint: enable=bad-whitespace,line-too-long 
    8383 
    84 source = ["lib/sas_3j1x_x.c"] 
    85  
    86 # No volume normalization despite having a volume parameter 
    87 # This should perhaps be volume normalized? 
    88 form_volume = """ 
    89     return M_4PI_3*cube(radius); 
    90     """ 
    91  
    92 Iq = """ 
    93     const double qr = q*radius; 
    94     const double bes = sas_3j1x_x(qr); 
    95     const double qf = q*fuzziness; 
    96     const double fq = bes * (sld - sld_solvent) * form_volume(radius) * exp(-0.5*qf*qf); 
    97     return 1.0e-4*fq*fq; 
    98     """ 
    99  
    100 def ER(radius): 
    101     """ 
    102     Return radius 
    103     """ 
    104     return radius 
    105  
    106 # VR defaults to 1.0 
     84source = ["lib/sas_3j1x_x.c","fuzzy_sphere.c"] 
     85have_Fq = True 
     86effective_radius_type = ["radius", "radius + fuzziness"] 
    10787 
    10888def random(): 
  • sasmodels/models/gaussian_peak.py

    r2d81cfe r304c775  
    5050    """ 
    5151 
    52 # VR defaults to 1.0 
    53  
    5452def random(): 
    5553    peak_pos = 10**np.random.uniform(-3, -1) 
  • sasmodels/models/hardsphere.py

    r2d81cfe r304c775  
    6262 
    6363#             ["name", "units", default, [lower, upper], "type","description"], 
    64 parameters = [["radius_effective", "Ang", 50.0, [0, inf], "volume", 
     64parameters = [["radius_effective", "Ang", 50.0, [0, inf], "", 
    6565               "effective radius of hard sphere"], 
    6666              ["volfraction", "", 0.2, [0, 0.74], "", 
    6767               "volume fraction of hard spheres"], 
    6868             ] 
    69  
    70 # No volume normalization despite having a volume parameter 
    71 # This should perhaps be volume normalized? 
    72 form_volume = """ 
    73     return 1.0; 
    74     """ 
    7569 
    7670Iq = r""" 
     
    168162    return pars 
    169163 
    170 # ER defaults to 0.0 
    171 # VR defaults to 1.0 
    172  
    173164demo = dict(radius_effective=200, volfraction=0.2, 
    174165            radius_effective_pd=0.1, radius_effective_pd_n=40) 
  • sasmodels/models/hayter_msa.py

    r2d81cfe r304c775  
    8989    return 1.0; 
    9090    """ 
    91 # ER defaults to 0.0 
    92 # VR defaults to 1.0 
    9391 
    9492def random(): 
  • sasmodels/models/hollow_cylinder.c

    r108e70e rd42dd4a  
    11//#define INVALID(v) (v.radius_core >= v.radius) 
    22 
    3 // From Igor library 
    43static double 
    5 _hollow_cylinder_scaling(double integrand, double delrho, double volume) 
     4shell_volume(double radius, double thickness, double length) 
    65{ 
    7     return 1.0e-4 * square(volume * delrho) * integrand; 
     6    return M_PI*length*(square(radius+thickness) - radius*radius); 
     7} 
     8 
     9static double 
     10form_volume(double radius, double thickness, double length) 
     11{ 
     12    return M_PI*length*square(radius+thickness); 
     13} 
     14 
     15static double 
     16radius_from_volume(double radius, double thickness, double length) 
     17{ 
     18    const double volume_outer_cyl = M_PI*square(radius + thickness)*length; 
     19    return cbrt(volume_outer_cyl/M_4PI_3); 
     20} 
     21 
     22static double 
     23radius_from_diagonal(double radius, double thickness, double length) 
     24{ 
     25    return sqrt(square(radius + thickness) + 0.25*square(length)); 
     26} 
     27 
     28static double 
     29effective_radius(int mode, double radius, double thickness, double length) 
     30{ 
     31    switch (mode) { 
     32    default: 
     33    case 1: // equivalent sphere 
     34        return radius_from_volume(radius, thickness, length); 
     35    case 2: // outer radius 
     36        return radius + thickness; 
     37    case 3: // half length 
     38        return 0.5*length; 
     39    case 4: // half outer min dimension 
     40        return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length); 
     41    case 5: // half outer max dimension 
     42        return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length); 
     43    case 6: // half outer diagonal 
     44        return radius_from_diagonal(radius,thickness,length); 
     45    } 
    846} 
    947 
     
    2260} 
    2361 
    24 static double 
    25 form_volume(double radius, double thickness, double length) 
    26 { 
    27     double v_shell = M_PI*length*(square(radius+thickness) - radius*radius); 
    28     return v_shell; 
    29 } 
    30  
    31  
    32 static double 
    33 Iq(double q, double radius, double thickness, double length, 
     62static void 
     63Fq(double q, double *F1, double *F2, double radius, double thickness, double length, 
    3464    double sld, double solvent_sld) 
    3565{ 
     
    3767    const double upper = 1.0;        //limits of numerical integral 
    3868 
    39     double summ = 0.0;            //initialize intergral 
     69    double total_F1 = 0.0;            //initialize intergral 
     70    double total_F2 = 0.0; 
    4071    for (int i=0;i<GAUSS_N;i++) { 
    4172        const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper ); 
     
    4374        const double form = _fq(q*sin_theta, q*cos_theta, 
    4475                                radius, thickness, length); 
    45         summ += GAUSS_W[i] * form * form; 
     76        total_F1 += GAUSS_W[i] * form; 
     77        total_F2 += GAUSS_W[i] * form * form; 
    4678    } 
     79    total_F1 *= 0.5*(upper-lower); 
     80    total_F2 *= 0.5*(upper-lower); 
     81    const double s = (sld - solvent_sld) * shell_volume(radius, thickness, length); 
     82    *F1 = 1e-2 * s * total_F1; 
     83    *F2 = 1e-4 * s*s * total_F2; 
     84} 
    4785 
    48     const double Aq = 0.5*summ*(upper-lower); 
    49     const double volume = form_volume(radius, thickness, length); 
    50     return _hollow_cylinder_scaling(Aq, solvent_sld - sld, volume); 
    51 } 
    5286 
    5387static double 
     
    5791{ 
    5892    const double form = _fq(qab, qc, radius, thickness, length); 
    59  
    60     const double vol = form_volume(radius, thickness, length); 
    61     return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
     93    const double s = (sld - solvent_sld) * shell_volume(radius, thickness, length); 
     94    return 1.0e-4*square(s * form); 
    6295} 
  • sasmodels/models/hollow_cylinder.py

    r455aaa1 r304c775  
    6969* **Last Reviewed by:** Paul Butler **Date:** September 06, 2018 
    7070""" 
     71from __future__ import division 
    7172 
    7273import numpy as np 
     
    99100 
    100101source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
    101  
    102 # pylint: disable=W0613 
    103 def ER(radius, thickness, length): 
    104     """ 
    105     :param radius:      Cylinder core radius 
    106     :param thickness:   Cylinder wall thickness 
    107     :param length:      Cylinder length 
    108     :return:            Effective radius 
    109     """ 
    110     router = radius + thickness 
    111     if router == 0 or length == 0: 
    112         return 0.0 
    113     len1 = router 
    114     len2 = length/2.0 
    115     term1 = len1*len1*2.0*len2/2.0 
    116     term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0) 
    117     ddd = 3.0*term1*term2 
    118     diam = pow(ddd, (1.0/3.0)) 
    119     return diam 
    120  
    121 def VR(radius, thickness, length): 
    122     """ 
    123     :param radius:      Cylinder radius 
    124     :param thickness:   Cylinder wall thickness 
    125     :param length:      Cylinder length 
    126     :return:            Volf ratio for P(q)*S(q) 
    127     """ 
    128     router = radius + thickness 
    129     vol_core = pi*radius*radius*length 
    130     vol_total = pi*router*router*length 
    131     vol_shell = vol_total - vol_core 
    132     return vol_total, vol_shell 
     102have_Fq = True 
     103effective_radius_type = [ 
     104    "equivalent sphere", "outer radius", "half length", 
     105    "half outer min dimension", "half outer max dimension", 
     106    "half outer diagonal", 
     107    ] 
    133108 
    134109def random(): 
     
    158133qx = q*cos(pi/6.0) 
    159134qy = q*sin(pi/6.0) 
     135radius = parameters[0][2] 
     136thickness = parameters[1][2] 
     137length = parameters[2][2] 
    160138# Parameters for unit tests 
    161139tests = [ 
    162140    [{}, 0.00005, 1764.926], 
    163     [{}, 'VR', 0.55555556], 
     141    [{}, 0.1, None, None, 
     142     (3./4*(radius+thickness)**2*length)**(1./3),  # R_eff from volume 
     143     pi*((radius+thickness)**2-radius**2)*length,  # shell volume 
     144     (radius+thickness)**2/((radius+thickness)**2 - radius**2), # form:shell ratio 
     145    ], 
    164146    [{}, 0.001, 1756.76], 
    165147    [{}, (qx, qy), 2.36885476192], 
  • sasmodels/models/hollow_rectangular_prism.c

    rd86f0fc rd42dd4a  
     1static double 
     2shell_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     3{ 
     4    const double length_b = length_a * b2a_ratio; 
     5    const double length_c = length_a * c2a_ratio; 
     6    const double form_volume = length_a * length_b * length_c; 
     7    const double a_core = length_a - 2.0*thickness; 
     8    const double b_core = length_b - 2.0*thickness; 
     9    const double c_core = length_c - 2.0*thickness; 
     10    const double core_volume = a_core * b_core * c_core; 
     11    return form_volume - core_volume; 
     12} 
     13 
    114static double 
    215form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    316{ 
    4     double length_b = length_a * b2a_ratio; 
    5     double length_c = length_a * c2a_ratio; 
    6     double a_core = length_a - 2.0*thickness; 
    7     double b_core = length_b - 2.0*thickness; 
    8     double c_core = length_c - 2.0*thickness; 
    9     double vol_core = a_core * b_core * c_core; 
    10     double vol_total = length_a * length_b * length_c; 
    11     double vol_shell = vol_total - vol_core; 
    12     return vol_shell; 
     17    const double length_b = length_a * b2a_ratio; 
     18    const double length_c = length_a * c2a_ratio; 
     19    const double form_volume = length_a * length_b * length_c; 
     20    return form_volume; 
    1321} 
    1422 
    1523static double 
    16 Iq(double q, 
     24effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     25// NOTE length_a is external dimension! 
     26{ 
     27    switch (mode) { 
     28    default: 
     29    case 1: // equivalent sphere 
     30        return cbrt(cube(length_a)*b2a_ratio*c2a_ratio/M_4PI_3); 
     31    case 2: // half length_a 
     32        return 0.5 * length_a; 
     33    case 3: // half length_b 
     34        return 0.5 * length_a*b2a_ratio; 
     35    case 4: // half length_c 
     36        return 0.5 * length_a*c2a_ratio; 
     37    case 5: // equivalent outer circular cross-section 
     38        return length_a*sqrt(b2a_ratio/M_PI); 
     39    case 6: // half ab diagonal 
     40        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     41    case 7: // half diagonal 
     42        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     43    } 
     44} 
     45 
     46static void 
     47Fq(double q, 
     48    double *F1, 
     49    double *F2, 
    1750    double sld, 
    1851    double solvent_sld, 
     
    3669    const double v2b = M_PI_2;  //phi integration limits 
    3770 
    38     double outer_sum = 0.0; 
     71    double outer_sum_F1 = 0.0; 
     72    double outer_sum_F2 = 0.0; 
    3973    for(int i=0; i<GAUSS_N; i++) { 
    4074 
     
    4680        const double termC2 = sas_sinx_x(q * (c_half-thickness)*cos(theta)); 
    4781 
    48         double inner_sum = 0.0; 
     82        double inner_sum_F1 = 0.0; 
     83        double inner_sum_F2 = 0.0; 
    4984        for(int j=0; j<GAUSS_N; j++) { 
    5085 
     
    6499            const double AP2 = vol_core * termA2 * termB2 * termC2; 
    65100 
    66             inner_sum += GAUSS_W[j] * square(AP1-AP2); 
     101            inner_sum_F1 += GAUSS_W[j] * (AP1-AP2); 
     102            inner_sum_F2 += GAUSS_W[j] * square(AP1-AP2); 
    67103        } 
    68         inner_sum *= 0.5 * (v2b-v2a); 
     104        inner_sum_F1 *= 0.5 * (v2b-v2a); 
     105        inner_sum_F2 *= 0.5 * (v2b-v2a); 
    69106 
    70         outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 
     107        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * sin(theta); 
     108        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * sin(theta); 
    71109    } 
    72     outer_sum *= 0.5*(v1b-v1a); 
     110    outer_sum_F1 *= 0.5*(v1b-v1a); 
     111    outer_sum_F2 *= 0.5*(v1b-v1a); 
    73112 
    74113    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 
    75114    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 
    76     const double form = outer_sum/M_PI_2; 
     115    const double form_avg = outer_sum_F1/M_PI_2; 
     116    const double form_squared_avg = outer_sum_F2/M_PI_2; 
    77117 
    78118    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
    79     const double delrho = sld - solvent_sld; 
     119    const double contrast = sld - solvent_sld; 
    80120 
    81121    // Convert from [1e-12 A-1] to [cm-1] 
    82     return 1.0e-4 * delrho * delrho * form; 
     122    *F1 = 1.0e-2 * contrast * form_avg; 
     123    *F2 = 1.0e-4 * contrast * contrast * form_squared_avg; 
    83124} 
    84125 
  • sasmodels/models/hollow_rectangular_prism.py

    r455aaa1 ree60aa7  
    132132               "Solvent scattering length density"], 
    133133              ["length_a", "Ang", 35, [0, inf], "volume", 
    134                "Shorter side of the parallelepiped"], 
     134               "Shortest, external, size of the parallelepiped"], 
    135135              ["b2a_ratio", "Ang", 1, [0, inf], "volume", 
    136136               "Ratio sides b/a"], 
     
    148148 
    149149source = ["lib/gauss76.c", "hollow_rectangular_prism.c"] 
    150  
    151 def ER(length_a, b2a_ratio, c2a_ratio, thickness): 
    152     """ 
    153     Return equivalent radius (ER) 
    154     thickness parameter not used 
    155     """ 
    156     b_side = length_a * b2a_ratio 
    157     c_side = length_a * c2a_ratio 
    158  
    159     # surface average radius (rough approximation) 
    160     surf_rad = sqrt(length_a * b_side / pi) 
    161  
    162     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    163     return 0.5 * (ddd) ** (1. / 3.) 
    164  
    165 def VR(length_a, b2a_ratio, c2a_ratio, thickness): 
    166     """ 
    167     Return shell volume and total volume 
    168     """ 
    169     b_side = length_a * b2a_ratio 
    170     c_side = length_a * c2a_ratio 
    171     a_core = length_a - 2.0*thickness 
    172     b_core = b_side - 2.0*thickness 
    173     c_core = c_side - 2.0*thickness 
    174     vol_core = a_core * b_core * c_core 
    175     vol_total = length_a * b_side * c_side 
    176     vol_shell = vol_total - vol_core 
    177     return vol_total, vol_shell 
    178  
     150have_Fq = True 
     151effective_radius_type = [ 
     152    "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     153    "equivalent outer circular cross-section", 
     154    "half ab diagonal", "half diagonal", 
     155    ] 
    179156 
    180157def random(): 
  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    rd86f0fc rd42dd4a  
     1static double 
     2shell_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     3{ 
     4    const double length_b = length_a * b2a_ratio; 
     5    const double length_c = length_a * c2a_ratio; 
     6    const double shell_volume = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
     7    return shell_volume; 
     8} 
     9 
    110static double 
    211form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    312{ 
    4     double length_b = length_a * b2a_ratio; 
    5     double length_c = length_a * c2a_ratio; 
    6     double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
    7     return vol_shell; 
     13    const double length_b = length_a * b2a_ratio; 
     14    const double length_c = length_a * c2a_ratio; 
     15    const double form_volume = length_a * length_b * length_c; 
     16    return form_volume; 
    817} 
    918 
    1019static double 
    11 Iq(double q, 
     20effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio) 
     21{ 
     22    switch (mode) { 
     23    default: 
     24    case 1: // equivalent sphere 
     25        return cbrt(cube(length_a)*b2a_ratio*c2a_ratio/M_4PI_3); 
     26    case 2: // half length_a 
     27        return 0.5 * length_a; 
     28    case 3: // half length_b 
     29        return 0.5 * length_a*b2a_ratio; 
     30    case 4: // half length_c 
     31        return 0.5 * length_a*c2a_ratio; 
     32    case 5: // equivalent outer circular cross-section 
     33        return length_a*sqrt(b2a_ratio/M_PI); 
     34    case 6: // half ab diagonal 
     35        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     36    case 7: // half diagonal 
     37        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     38    } 
     39} 
     40 
     41static void 
     42Fq(double q, 
     43    double *F1, 
     44    double *F2, 
    1245    double sld, 
    1346    double solvent_sld, 
     
    2861    const double v2b = M_PI_2;  //phi integration limits 
    2962 
    30     double outer_sum = 0.0; 
     63    double outer_sum_F1 = 0.0; 
     64    double outer_sum_F2 = 0.0; 
    3165    for(int i=0; i<GAUSS_N; i++) { 
    3266        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     
    4175        const double termAT_theta = 8.0 * sin_c / (q*q*sin_theta*cos_theta); 
    4276 
    43         double inner_sum = 0.0; 
     77        double inner_sum_F1 = 0.0; 
     78        double inner_sum_F2 = 0.0; 
    4479        for(int j=0; j<GAUSS_N; j++) { 
    4580            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     
    6095                * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 
    6196 
    62             inner_sum += GAUSS_W[j] * square(AL+AT); 
     97            inner_sum_F1 += GAUSS_W[j] * (AL+AT); 
     98            inner_sum_F2 += GAUSS_W[j] * square(AL+AT); 
    6399        } 
    64100 
    65         inner_sum *= 0.5 * (v2b-v2a); 
    66         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     101        inner_sum_F1 *= 0.5 * (v2b-v2a); 
     102        inner_sum_F2 *= 0.5 * (v2b-v2a); 
     103        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * sin_theta; 
     104        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * sin_theta; 
    67105    } 
    68106 
    69     outer_sum *= 0.5*(v1b-v1a); 
     107    outer_sum_F1 *= 0.5*(v1b-v1a); 
     108    outer_sum_F2 *= 0.5*(v1b-v1a); 
    70109 
    71110    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 
    72111    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 
    73     double answer = outer_sum/M_PI_2; 
     112    const double form_avg = outer_sum_F1/M_PI_2; 
     113    const double form_squared_avg = outer_sum_F2/M_PI_2; 
    74114 
    75115    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
    76     answer *= square(sld-solvent_sld); 
     116    const double contrast = sld - solvent_sld; 
    77117 
    78118    // Convert from [1e-12 A-1] to [cm-1] 
    79     answer *= 1.0e-4; 
    80  
    81     return answer; 
     119    *F1 = 1e-2 * contrast * form_avg; 
     120    *F2 = 1e-4 * contrast * contrast * form_squared_avg; 
    82121} 
  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    r6e7d7b6 ree60aa7  
    108108 
    109109source = ["lib/gauss76.c", "hollow_rectangular_prism_thin_walls.c"] 
    110  
    111 def ER(length_a, b2a_ratio, c2a_ratio): 
    112     """ 
    113         Return equivalent radius (ER) 
    114     """ 
    115     b_side = length_a * b2a_ratio 
    116     c_side = length_a * c2a_ratio 
    117  
    118     # surface average radius (rough approximation) 
    119     surf_rad = sqrt(length_a * b_side / pi) 
    120  
    121     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    122     return 0.5 * (ddd) ** (1. / 3.) 
    123  
    124 def VR(length_a, b2a_ratio, c2a_ratio): 
    125     """ 
    126         Return shell volume and total volume 
    127     """ 
    128     b_side = length_a * b2a_ratio 
    129     c_side = length_a * c2a_ratio 
    130     vol_total = length_a * b_side * c_side 
    131     vol_shell = 2.0 * (length_a*b_side + length_a*c_side + b_side*c_side) 
    132     return vol_shell, vol_total 
     110have_Fq = True 
     111effective_radius_type = [ 
     112    "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     113    "equivalent outer circular cross-section", 
     114    "half ab diagonal", "half diagonal", 
     115    ] 
    133116 
    134117 
  • sasmodels/models/lamellar.py

    r2d81cfe r304c775  
    9595    return pars 
    9696 
    97 # ER defaults to 0.0 
    98 # VR defaults to 1.0 
    99  
    10097demo = dict(scale=1, background=0, 
    10198            sld=6, sld_solvent=1, 
  • sasmodels/models/lamellar_hg.py

    r2d81cfe r304c775  
    109109    return pars 
    110110 
    111 # ER defaults to 0.0 
    112 # VR defaults to 1.0 
    113  
    114111demo = dict(scale=1, background=0, 
    115112            length_tail=15, length_head=10, 
  • sasmodels/models/lamellar_hg_stack_caille.c

    r1c6286d r71b751d  
    5656  return inten; 
    5757} 
    58  
  • sasmodels/models/lamellar_hg_stack_caille.py

    r2d81cfe r304c775  
    122122    """ 
    123123 
    124 # ER defaults to 0.0 
    125 # VR defaults to 1.0 
    126  
    127124def random(): 
    128125    total_thickness = 10**np.random.uniform(2, 4.7) 
  • sasmodels/models/lamellar_stack_caille.c

    r1c6286d r71b751d  
    5252  return(inten); 
    5353} 
    54  
  • sasmodels/models/lamellar_stack_caille.py

    r2d81cfe r304c775  
    120120    """ 
    121121 
    122 # ER defaults to 0.0 
    123 # VR defaults to 1.0 
    124  
    125122demo = dict(scale=1, background=0, 
    126123            thickness=67., Nlayers=3.75, d_spacing=200., 
  • sasmodels/models/lamellar_stack_paracrystal.c

    r5467cd8 r71b751d  
    1818        int n2 = n1 + 1; 
    1919        const double xn = (double)n2 - fp_Nlayers;      //fractional contribution of n1 
    20          
     20 
    2121        const double ww = exp(-0.5*square(qval*pd*davg)); 
    2222 
     
    2727 
    2828        Znq = xn*Snq; 
    29          
     29 
    3030        //calculate the n2 contribution 
    3131        an = paraCryst_an(ww,qval,davg,n2); 
     
    3333 
    3434        Znq += (1.0-xn)*Snq; 
    35          
     35 
    3636        //and the independent contribution 
    3737        Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg)); 
    38          
     38 
    3939        //the limit when Nlayers approaches infinity 
    4040//      Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg)) 
    41          
     41 
    4242        const double xi = th/2.0;               //use 1/2 the bilayer thickness 
    4343        const double Pbil = square(sas_sinx_x(qval*xi)); 
    44          
     44 
    4545        const double contr = sld - solvent_sld; 
    4646        const double inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval); 
     
    5252double 
    5353paraCryst_sn(double ww, double qval, double davg, int Nlayers, double an) { 
    54          
     54 
    5555        double Snq; 
    5656 
    5757        Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) ); 
    58          
     58 
    5959        return Snq; 
    6060} 
     
    6363paraCryst_an(double ww, double qval, double davg, int Nlayers) { 
    6464        double an; 
    65          
     65 
    6666        an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg); 
    6767        an -= 4.0*pow(ww,(Nlayers+2))*cos((double)Nlayers*qval*davg); 
    6868        an += 2.0*pow(ww,(Nlayers+3))*cos((double)(Nlayers-1)*qval*davg); 
    6969        an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg); 
    70          
     70 
    7171        return an; 
    7272} 
  • sasmodels/models/lamellar_stack_paracrystal.py

    r2d81cfe r304c775  
    147147    return pars 
    148148 
    149 # ER defaults to 0.0 
    150 # VR defaults to 1.0 
    151  
    152149demo = dict(scale=1, background=0, 
    153150            thickness=33, Nlayers=20, d_spacing=250, sigma_d=0.2, 
  • sasmodels/models/lib/core_shell.c

    r925ad6e r71b751d  
    77********************************************************************/ 
    88static 
    9 double core_shell_kernel(double q, 
     9double core_shell_fq(double q, 
    1010                         double radius, 
    1111                         double thickness, 
    1212                         double core_sld, 
    1313                         double shell_sld, 
    14                          double solvent_sld) { 
     14                         double solvent_sld) 
     15{ 
    1516    // Core first, then add in shell 
    1617    const double core_qr = q * radius; 
     
    2627    const double shell_volume = M_4PI_3 * cube(radius + thickness); 
    2728    f += shell_volume * shell_bes * shell_contrast; 
    28     return f * f * 1.0e-4; 
     29    return f; 
    2930} 
     31 
     32// Deprecated function: use core_shell_fq instead. 
     33static 
     34double core_shell_kernel(double q, 
     35                         double radius, 
     36                         double thickness, 
     37                         double core_sld, 
     38                         double shell_sld, 
     39                         double solvent_sld) 
     40{ 
     41    const double fq = core_shell_fq(q, radius, thickness, core_sld, shell_sld, solvent_sld); 
     42    return 1.0e-4 * fq*fq; 
     43} 
  • sasmodels/models/lib/gauss76.c

    r99b84ec r74e9b5f  
    1111 
    1212// Gaussians 
    13 constant double Gauss76Wt[76]={ 
     13constant double Gauss76Wt[76] = { 
    1414        .00126779163408536,             //0 
    1515        .00294910295364247, 
     
    9090}; 
    9191 
    92 constant double Gauss76Z[76]={ 
     92constant double Gauss76Z[76] = { 
    9393        -.999505948362153,              //0 
    9494        -.997397786355355, 
  • sasmodels/models/lib/polevl.c

    r447e9aa r74e9b5f  
    5151*/ 
    5252 
    53 double polevl( double x, constant double *coef, int N ); 
    54 double polevl( double x, constant double *coef, int N ) 
     53static 
     54double polevl( double x, pconstant double *coef, int N ) 
    5555{ 
    5656 
     
    7272 */ 
    7373 
    74 double p1evl( double x, constant double *coef, int N ); 
    75 double p1evl( double x, constant double *coef, int N ) 
     74static 
     75double p1evl( double x, pconstant double *coef, int N ) 
    7676{ 
    7777    int i=0; 
  • sasmodels/models/lib/sas_J1.c

    r5181ccc r74e9b5f  
    4242#if FLOAT_SIZE>4 
    4343//Cephes double pression function 
    44 double cephes_j1(double x); 
    4544 
    4645constant double RPJ1[8] = { 
     
    106105    0.0 }; 
    107106 
     107static 
    108108double cephes_j1(double x) 
    109109{ 
     
    155155#else 
    156156//Single precission version of cephes 
    157 float cephes_j1f(float x); 
    158  
    159157constant float JPJ1[8] = { 
    160158    -4.878788132172128E-009, 
     
    190188    }; 
    191189 
     190static 
    192191float cephes_j1f(float xx) 
    193192{ 
     
    240239 
    241240//Finally J1c function that equals 2*J1(x)/x 
    242 double sas_2J1x_x(double x); 
     241static 
    243242double sas_2J1x_x(double x) 
    244243{ 
  • sasmodels/models/lib/sphere_form.c

    r925ad6e r01c8d9e  
    1313    return 1.0e-4*square(contrast * fq); 
    1414} 
     15 
  • sasmodels/models/linear_pearls.c

    rc3ccaec r71b751d  
    4646    double psi = sas_3j1x_x(q * radius); 
    4747 
    48     // N pearls interaction terms  
     48    // N pearls interaction terms 
    4949    double structure_factor = (double)num_pearls; 
    5050    for(int num=1; num<num_pearls; num++) { 
  • sasmodels/models/mass_surface_fractal.py

    r7994359 r71b751d  
    3838 
    3939    The surface ( $D_s$ ) and mass ( $D_m$ ) fractal dimensions are only 
    40     valid if $0 < surface\_dim < 6$ , $0 < mass\_dim < 6$ , and 
    41     $(surface\_dim + mass\_dim ) < 6$ .  
     40    valid if $0 < surface\_dim < 6$, $0 < mass\_dim < 6$, and 
     41    $(surface\_dim + mass\_dim ) < 6$. 
    4242    Older versions of sasview may have the default primary particle radius 
    43     larger than the cluster radius, this was an error, also present in the  
    44     Schmidt review paper below. The primary particle should be the smaller  
    45     as described in the original Hurd et.al. who also point out that  
    46     polydispersity in the primary particle sizes may affect their  
     43    larger than the cluster radius, this was an error, also present in the 
     44    Schmidt review paper below. The primary particle should be the smaller 
     45    as described in the original Hurd, et al., who also point out that 
     46    polydispersity in the primary particle sizes may affect their 
    4747    apparent surface fractal dimension. 
    48      
     48 
    4949 
    5050References 
  • sasmodels/models/mono_gauss_coil.py

    r2d81cfe