Changes in / [b9c7379:a08c47f] in sasmodels


Ignore:
Files:
8 added
6 deleted
59 edited

Legend:

Unmodified
Added
Removed

  • README.rst

    re30d645 r2a64722  
    1010is available. 
    1111 
    12 Example 
     12Install 
    1313------- 
     14 
     15The easiest way to use sasmodels is from `SasView <http://www.sasview.org/>`_. 
     16 
     17You can also install sasmodels as a standalone package in python. Use 
     18`miniconda <https://docs.conda.io/en/latest/miniconda.html>`_ 
     19or `anaconda <https://www.anaconda.com/>`_ 
     20to create a python environment with the sasmodels dependencies:: 
     21 
     22    $ conda create -n sasmodels -c conda-forge numpy scipy matplotlib pyopencl 
     23 
     24The option ``-n sasmodels`` names the environment sasmodels, and the option 
     25``-c conda-forge`` selects the conda-forge package channel because pyopencl 
     26is not part of the base anaconda distribution. 
     27 
     28Activate the environment and install sasmodels:: 
     29 
     30    $ conda activate sasmodels 
     31    (sasmodels) $ pip install sasmodels 
     32 
     33Install `bumps <https://github.com/bumps/bumps>`_ if you want to use it to fit 
     34your data:: 
     35 
     36    (sasmodels) $ pip install bumps 
     37 
     38Usage 
     39----- 
     40 
     41Check that the works:: 
     42 
     43    (sasmodels) $ python -m sasmodels.compare cylinder 
     44 
     45To show the orientation explorer:: 
     46 
     47    (sasmodels) $ python -m sasmodels.jitter 
     48 
     49Documentation is available online as part of the SasView 
     50`fitting perspective <http://www.sasview.org/docs/index.html>`_ 
     51as well as separate pages for 
     52`individual models <http://www.sasview.org/docs/user/sasgui/perspectives/fitting/models/index.html>`_. 
     53Programming details for sasmodels are available in the 
     54`developer documentation <http://www.sasview.org/docs/dev/dev.html>`_. 
     55 
     56 
     57Fitting Example 
     58--------------- 
    1459 
    1560The example directory contains a radial+tangential data set for an oriented 
    1661rod-like shape. 
    1762 
    18 The data is loaded by sas.dataloader from the sasview package, so sasview 
    19 is needed to run the example. 
     63To load the example data, you will need the SAS data loader from the sasview 
     64package. This is not yet available on PyPI, so you will need a copy of the 
     65SasView source code to run it.  Create a directory somewhere to hold the 
     66sasview and sasmodels source code, which we will refer to as $SOURCE. 
    2067 
    21 To run the example, you need sasview, sasmodels and bumps.  Assuming these 
    22 repositories are installed side by side, change to the sasmodels/example 
    23 directory and enter:: 
     68Use the following to install sasview, and the sasmodels examples:: 
    2469 
    25     PYTHONPATH=..:../../sasview/src ../../bumps/run.py fit.py \ 
    26         cylinder --preview 
     70    (sasmodels) $ cd $SOURCE 
     71    (sasmodels) $ conda install git 
     72    (sasmodels) $ git clone https://github.com/sasview/sasview.git 
     73    (sasmodels) $ git clone https://github.com/sasview/sasmodels.git 
    2774 
    28 See bumps documentation for instructions on running the fit.  With the 
    29 python packages installed, e.g., into a virtual environment, then the 
    30 python path need not be set, and the command would be:: 
     75Set the path to the sasview source on your python path within the sasmodels 
     76environment.  On Windows, this will be:: 
    3177 
    32     bumps fit.py cylinder --preview 
     78    (sasmodels)> set PYTHONPATH="$SOURCE\sasview\src" 
     79    (sasmodels)> cd $SOURCE/sasmodels/example 
     80    (sasmodels)> python -m bumps.cli fit.py cylinder --preview 
     81 
     82On Mac/Linux with the standard shell this will be:: 
     83 
     84    (sasmodels) $ export PYTHONPATH="$SOURCE/sasview/src" 
     85    (sasmodels) $ cd $SOURCE/sasmodels/example 
     86    (sasmodels) $ bumps fit.py cylinder --preview 
    3387 
    3488The fit.py model accepts up to two arguments.  The first argument is the 
    … …  
    3892both radial and tangential simultaneously, use the word "both". 
    3993 
    40 Notes 
    41 ----- 
    42  
    43 cylinder.c + cylinder.py is the cylinder model with renamed variables and 
    44 sld scaled by 1e6 so the numbers are nicer.  The model name is "cylinder" 
    45  
    46 lamellar.py is an example of a single file model with embedded C code. 
     94See `bumps documentation <https://bumps.readthedocs.io/>`_ for detailed 
     95instructions on running the fit. 
    4796 
    4897|TravisStatus|_ 

  • doc/guide/plugin.rst

    r81751c2 r81751c2  
    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:: 

  • doc/guide/scripting.rst

    rbd7630d r23df833  
    188188python kernel.  Once the kernel is in hand, we can then marshal a set of 
    189189parameters into a :class:`sasmodels.details.CallDetails` object and ship it to 
    190 the kernel using the :func:`sansmodels.direct_model.call_kernel` function.  An 
    191 example should help, *example/cylinder_eval.py*:: 
    192  
    193     from numpy import logspace 
     190the kernel using the :func:`sansmodels.direct_model.call_kernel` function.  To 
     191accesses the underlying $<F(q)>$ and $<F^2(q)>$, use 
     192:func:`sasmodels.direct_model.call_Fq` instead. 
     193 
     194The following example should 
     195help, *example/cylinder_eval.py*:: 
     196 
     197    from numpy import logspace, sqrt 
    194198    from matplotlib import pyplot as plt 
    195199    from sasmodels.core import load_model 
    196     from sasmodels.direct_model import call_kernel 
     200    from sasmodels.direct_model import call_kernel, call_Fq 
    197201 
    198202    model = load_model('cylinder') 
    199203    q = logspace(-3, -1, 200) 
    200204    kernel = model.make_kernel([q]) 
    201     Iq = call_kernel(kernel, dict(radius=200.)) 
    202     plt.loglog(q, Iq) 
     205    pars = {'radius': 200, 'radius_pd': 0.1, 'scale': 2} 
     206    Iq = call_kernel(kernel, pars) 
     207    F, Fsq, Reff, V, Vratio = call_Fq(kernel, pars) 
     208 
     209    plt.loglog(q, Iq, label='2 I(q)') 
     210    plt.loglog(q, F**2/V, label='<F(q)>^2/V') 
     211    plt.loglog(q, Fsq/V, label='<F^2(q)>/V') 
     212    plt.xlabel('q (1/A)') 
     213    plt.ylabel('I(q) (1/cm)') 
     214    plt.title('Cylinder with radius 200.') 
     215    plt.legend() 
    203216    plt.show() 
    204217 
    205 On windows, this can be called from the cmd prompt using sasview as:: 
     218.. figure:: direct_call.png 
     219 
     220    Comparison between $I(q)$, $<F(q)>$ and $<F^2(q)>$ for cylinder model. 
     221 
     222This compares $I(q)$ with $<F(q)>$ and $<F^2(q)>$ for a cylinder 
     223with *radius=200 +/- 20* and *scale=2*. Note that *call_Fq* does not 
     224include scale and background, nor does it normalize by the average volume. 
     225The definition of $F = \rho V \hat F$ scaled by the contrast and 
     226volume, compared to the canonical cylinder $\hat F$, with $\hat F(0) = 1$. 
     227Integrating over polydispersity and orientation, the returned values are 
     228$\sum_{r,w\in N(r_o, r_o/10)} \sum_\theta w F(q,r_o,\theta)\sin\theta$ and 
     229$\sum_{r,w\in N(r_o, r_o/10)} \sum_\theta w F^2(q,r_o,\theta)\sin\theta$. 
     230 
     231On windows, this example can be called from the cmd prompt using sasview as 
     232as the python interpreter:: 
    206233 
    207234    SasViewCom example/cylinder_eval.py 

  • example/cylinder_eval.py

    r2e66ef5 r23df833  
    33""" 
    44 
    5 from numpy import logspace 
     5from numpy import logspace, sqrt 
    66from matplotlib import pyplot as plt 
    77from sasmodels.core import load_model 
    8 from sasmodels.direct_model import call_kernel 
     8from sasmodels.direct_model import call_kernel, call_Fq 
    99 
    1010model = load_model('cylinder') 
    1111q = logspace(-3, -1, 200) 
    1212kernel = model.make_kernel([q]) 
    13 Iq = call_kernel(kernel, dict(radius=200.)) 
    14 plt.loglog(q, Iq) 
     13pars = {'radius': 200, 'radius_pd': 0.1, 'scale': 2} 
     14Iq = call_kernel(kernel, pars) 
     15F, Fsq, Reff, V, Vratio = call_Fq(kernel, pars) 
     16plt.loglog(q, Iq, label='2 I(q)') 
     17plt.loglog(q, F**2/V, label='<F(q)>^2/V') 
     18plt.loglog(q, Fsq/V, label='<F^2(q)>/V') 
    1519plt.xlabel('q (1/A)') 
    16 plt.ylabel('I(q)') 
     20plt.ylabel('I(q) (1/cm)') 
    1721plt.title('Cylinder with radius 200.') 
     22plt.legend() 
    1823plt.show() 

  • explore/beta/sasfit_compare.py

    r2a12351b r119073a  
    505505    } 
    506506 
    507     Q, IQ = load_sasfit(data_file('richard_test.txt')) 
    508     Q, IQSD = load_sasfit(data_file('richard_test2.txt')) 
    509     Q, IQBD = load_sasfit(data_file('richard_test3.txt')) 
     507    Q, IQ = load_sasfit(data_file('sasfit_sphere_schulz_IQD.txt')) 
     508    Q, IQSD = load_sasfit(data_file('sasfit_sphere_schulz_IQSD.txt')) 
     509    Q, IQBD = load_sasfit(data_file('sasfit_sphere_schulz_IQBD.txt')) 
    510510    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
    511511    actual = sphere_r(Q, norm="sasfit", **pars) 
    … …  
    526526    } 
    527527 
    528     Q, IQ = load_sasfit(data_file('richard_test4.txt')) 
    529     Q, IQSD = load_sasfit(data_file('richard_test5.txt')) 
    530     Q, IQBD = load_sasfit(data_file('richard_test6.txt')) 
     528    Q, IQ = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQD.txt')) 
     529    Q, IQSD = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQSD.txt')) 
     530    Q, IQBD = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQBD.txt')) 
    531531    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
    532532    actual = ellipsoid_pe(Q, norm="sasfit", **pars) 

  • explore/precision.py

    ree60aa7 rcd28947  
    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 the gammainc and gammaincc functions. 
     210# Create a 2-D variant of the precision test if we need to handle other two 
     211# parameter functions. 
     212A = 1 
     213def parse_extra_pars(): 
     214    """ 
     215    Parse the command line looking for the second parameter "A=..." for the 
     216    gammainc/gammaincc functions. 
     217    """ 
     218    global A 
     219 
     220    A_str = str(A) 
     221    pop = [] 
     222    for k, v in enumerate(sys.argv[1:]): 
     223        if v.startswith("A="): 
     224            A_str = v[2:] 
     225            pop.append(k+1) 
     226    if pop: 
     227        sys.argv = [v for k, v in enumerate(sys.argv) if k not in pop] 
     228        A = float(A_str) 
     229 
     230parse_extra_pars() 
     231 
    199232 
    200233# =============== FUNCTION DEFINITIONS ================ 
    … …  
    297330    ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]), 
    298331    limits=(-3.1, 10), 
     332) 
     333add_function( 
     334    name="gammaln(x)", 
     335    mp_function=mp.loggamma, 
     336    np_function=scipy.special.gammaln, 
     337    ocl_function=make_ocl("return sas_gammaln(q);", "sas_gammaln", ["lib/sas_gammainc.c"]), 
     338    #ocl_function=make_ocl("return lgamma(q);", "sas_gammaln"), 
     339) 
     340add_function( 
     341    # Note: "a" is given as A=... on the command line via parse_extra_pars 
     342    name="gammainc(x)", 
     343    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a), 
     344    np_function=lambda x, a=A: scipy.special.gammainc(a, x), 
     345    ocl_function=make_ocl("return sas_gammainc(%.15g,q);"%A, "sas_gammainc", ["lib/sas_gammainc.c"]), 
     346) 
     347add_function( 
     348    # Note: "a" is given as A=... on the command line via parse_extra_pars 
     349    name="gammaincc(x)", 
     350    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a), 
     351    np_function=lambda x, a=A: scipy.special.gammaincc(a, x), 
     352    ocl_function=make_ocl("return sas_gammaincc(%.15g,q);"%A, "sas_gammaincc", ["lib/sas_gammainc.c"]), 
    299353) 
    300354add_function( 
    … …  
    465519lanczos_gamma = """\ 
    466520    const double coeff[] = { 
    467             76.18009172947146,     -86.50532032941677, 
    468             24.01409824083091,     -1.231739572450155, 
     521            76.18009172947146, -86.50532032941677, 
     522            24.01409824083091, -1.231739572450155, 
    469523            0.1208650973866179e-2,-0.5395239384953e-5 
    470524            }; 
    … …  
    477531""" 
    478532add_function( 
    479     name="log gamma(x)", 
     533    name="loggamma(x)", 
    480534    mp_function=mp.loggamma, 
    481535    np_function=scipy.special.gammaln, 
    … …  
    601655 
    602656ALL_FUNCTIONS = set(FUNCTIONS.keys()) 
    603 ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet 
     657ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead 
    604658ALL_FUNCTIONS.discard("3j1/x:taylor") 
    605659ALL_FUNCTIONS.discard("3j1/x:trig") 
    … …  
    617671    -r indicates that the relative error should be plotted (default), 
    618672    -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: 
     673        log indicates log stepping in [10^-3, 10^5] (default) 
     674        logq indicates log stepping in [10^-4, 10^1] 
     675        linear indicates linear stepping in [1, 1000] 
     676        zoom indicates linear stepping in [1000, 1010] 
     677        neg indicates linear stepping in [-100.1, 100.1] 
     678        start:stop:n[:stepping] indicates an n-step plot in [start, stop] 
     679            or [10^start, 10^stop] if stepping is "log" (default n=400) 
     680Some functions (notably gammainc/gammaincc) have an additional parameter A 
     681which can be set from the command line as A=value.  Default is A=1. 
     682 
     683Name is one of: 
    625684    """+names) 
    626685    sys.exit(1) 

  • sasmodels/compare.py

    r4de14584 r4de14584  
    11521152        'rel_err'   : True, 
    11531153        'explore'   : False, 
    1154         'use_demo'  : True, 
     1154        'use_demo'  : False, 
    11551155        'zero'      : False, 
    11561156        'html'      : False, 

  • sasmodels/direct_model.py

    r304c775 r5024a56  
    3232from .details import make_kernel_args, dispersion_mesh 
    3333from .modelinfo import DEFAULT_BACKGROUND 
     34from .product import RADIUS_MODE_ID 
    3435 
    3536# pylint: disable=unused-import 
    … …  
    6768    # type: (Kernel, ParameterSet, float, bool) -> np.ndarray 
    6869    """ 
    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)) 
     70    Like :func:`call_kernel`, but returning F, F^2, R_eff, V_shell, V_form/V_shell. 
     71 
     72    For solid objects V_shell is equal to V_form and the volume ratio is 1. 
     73 
     74    Use parameter *radius_effective_mode* to select the effective radius 
     75    calculation. 
     76    """ 
     77    R_eff_type = int(pars.pop(RADIUS_MODE_ID, 1.0)) 
    7278    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
    7379    #print("pars", list(zip(*mesh))[0]) 
    … …  
    7682    return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type) 
    7783 
    78 def call_profile(model_info, **pars): 
    79     # type: (ModelInfo, ...) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]] 
     84def call_profile(model_info, pars=None): 
     85    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]] 
    8086    """ 
    8187    Returns the profile *x, y, (xlabel, ylabel)* representing the model. 
    8288    """ 
     89    if pars is None: 
     90        pars = {} 
    8391    args = {} 
    8492    for p in model_info.parameters.kernel_parameters: 
    … …  
    324332 
    325333        # Need to pull background out of resolution for multiple scattering 
    326         background = pars.get('background', DEFAULT_BACKGROUND) 
     334        default_background = self._model.info.parameters.common_parameters[1].default 
     335        background = pars.get('background', default_background) 
    327336        pars = pars.copy() 
    328337        pars['background'] = 0. 
    … …  
    377386        Generate a plottable profile. 
    378387        """ 
    379         return call_profile(self.model.info, **pars) 
     388        return call_profile(self.model.info, pars) 
    380389 
    381390def main(): 

  • sasmodels/generate.py

    r39a06c9 rcd28947  
    703703    """ 
    704704    for code in source: 
    705         m = _FQ_PATTERN.search(code) 
    706         if m is not None: 
     705        if _FQ_PATTERN.search(code) is not None: 
    707706            return True 
    708707    return False 
    … …  
    712711    # type: (List[str]) -> bool 
    713712    """ 
    714     Return True if C source defines "void Fq(". 
     713    Return True if C source defines "double shell_volume(". 
    715714    """ 
    716715    for code in source: 
    717         m = _SHELL_VOLUME_PATTERN.search(code) 
    718         if m is not None: 
     716        if _SHELL_VOLUME_PATTERN.search(code) is not None: 
    719717            return True 
    720718    return False 
    … …  
    10081006        pars = model_info.parameters.kernel_parameters 
    10091007    else: 
    1010         pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters 
     1008        pars = (model_info.parameters.common_parameters 
     1009                + model_info.parameters.kernel_parameters) 
    10111010    partable = make_partable(pars) 
    10121011    subst = dict(id=model_info.id.replace('_', '-'), 

  • sasmodels/jitter.py

    r1198f90 r7d97437  
    1515    pass 
    1616 
     17import matplotlib as mpl 
    1718import matplotlib.pyplot as plt 
    1819from matplotlib.widgets import Slider 
    … …  
    746747        pass 
    747748 
    748     axcolor = 'lightgoldenrodyellow' 
     749    # CRUFT: use axisbg instead of facecolor for matplotlib<2 
     750    facecolor_prop = 'facecolor' if mpl.__version__ > '2' else 'axisbg' 
     751    props = {facecolor_prop: 'lightgoldenrodyellow'} 
    749752 
    750753    ## add control widgets to plot 
    751     axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    752     axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    753     axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) 
     754    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], **props) 
     755    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], **props) 
     756    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], **props) 
    754757    stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) 
    755758    sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) 
    756759    spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) 
    757760 
    758     axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    759     axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    760     axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
     761    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], **props) 
     762    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props) 
     763    axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props) 
    761764    # Note: using ridiculous definition of rectangle distribution, whose width 
    762765    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 

  • sasmodels/kernel.py

    re44432d rcd28947  
    133133        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
    134134        total_weight = self.result[nout*self.q_input.nq + 0] 
     135        # Note: total_weight = sum(weight > cutoff), with cutoff >= 0, so it 
     136        # is okay to test directly against zero.  If weight is zero then I(q), 
     137        # etc. must also be zero. 
    135138        if total_weight == 0.: 
    136139            total_weight = 1. 

  • sasmodels/kernelcl.py

    rf872fd1 rf872fd1  
    5858import time 
    5959 
     60try: 
     61    from time import perf_counter as clock 
     62except ImportError: # CRUFT: python < 3.3 
     63    import sys 
     64    if sys.platform.count("darwin") > 0: 
     65        from time import time as clock 
     66    else: 
     67        from time import clock 
     68 
    6069import numpy as np  # type: ignore 
    61  
    6270 
    6371# Attempt to setup opencl. This may fail if the pyopencl package is not 
    … …  
    611619        #Call kernel and retrieve results 
    612620        wait_for = None 
    613         last_nap = time.clock() 
     621        last_nap = clock() 
    614622        step = 1000000//self.q_input.nq + 1 
    615623        for start in range(0, call_details.num_eval, step): 
    … …  
    622630                # Allow other processes to run 
    623631                wait_for[0].wait() 
    624                 current_time = time.clock() 
     632                current_time = clock() 
    625633                if current_time - last_nap > 0.5: 
    626634                    time.sleep(0.001) 

  • sasmodels/kernelpy.py

    re44432d re44432d  
    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

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

  • sasmodels/modelinfo.py

    r39a06c9 r39a06c9  
    404404      parameters counted as n individual parameters p1, p2, ... 
    405405 
     406    * *common_parameters* is the list of common parameters, with a unique 
     407      copy for each model so that structure factors can have a default 
     408      background of 0.0. 
     409 
    406410    * *call_parameters* is the complete list of parameters to the kernel, 
    407411      including scale and background, with vector parameters recorded as 
    … …  
    416420    parameters don't use vector notation, and instead use p1, p2, ... 
    417421    """ 
    418     # scale and background are implicit parameters 
    419     COMMON = [Parameter(*p) for p in COMMON_PARAMETERS] 
    420  
    421422    def __init__(self, parameters): 
    422423        # type: (List[Parameter]) -> None 
     424 
     425        # scale and background are implicit parameters 
     426        # Need them to be unique to each model in case they have different 
     427        # properties, such as default=0.0 for structure factor backgrounds. 
     428        self.common_parameters = [Parameter(*p) for p in COMMON_PARAMETERS] 
     429 
    423430        self.kernel_parameters = parameters 
    424431        self._set_vector_lengths() 
    … …  
    468475                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    469476        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
     477 
     478    def set_zero_background(self): 
     479        """ 
     480        Set the default background to zero for this model.  This is done for 
     481        structure factor models. 
     482        """ 
     483        # type: () -> None 
     484        # Make sure background is the second common parameter. 
     485        assert self.common_parameters[1].id == "background" 
     486        self.common_parameters[1].default = 0.0 
     487        self.defaults = self._get_defaults() 
    470488 
    471489    def check_angles(self): 
    … …  
    567585    def _get_call_parameters(self): 
    568586        # type: () -> List[Parameter] 
    569         full_list = self.COMMON[:] 
     587        full_list = self.common_parameters[:] 
    570588        for p in self.kernel_parameters: 
    571589            if p.length == 1: 
    … …  
    670688 
    671689        # Gather the user parameters in order 
    672         result = control + self.COMMON 
     690        result = control + self.common_parameters 
    673691        for p in self.kernel_parameters: 
    674692            if not is2d and p.type in ('orientation', 'magnetic'): 
    … …  
    770788 
    771789    info = ModelInfo() 
     790 
     791    # Build the parameter table 
    772792    #print("make parameter table", kernel_module.parameters) 
    773793    parameters = make_parameter_table(getattr(kernel_module, 'parameters', [])) 
     794 
     795    # background defaults to zero for structure factor models 
     796    structure_factor = getattr(kernel_module, 'structure_factor', False) 
     797    if structure_factor: 
     798        parameters.set_zero_background() 
     799 
     800    # TODO: remove demo parameters 
     801    # The plots in the docs are generated from the model default values. 
     802    # Sascomp set parameters from the command line, and so doesn't need 
     803    # demo values for testing. 
    774804    demo = expand_pars(parameters, getattr(kernel_module, 'demo', None)) 
     805 
    775806    filename = abspath(kernel_module.__file__).replace('.pyc', '.py') 
    776807    kernel_id = splitext(basename(filename))[0] 

  • sasmodels/models/barbell.c

    rd42dd4a r99658f6  
    6363 
    6464static double 
     65radius_from_excluded_volume(double radius_bell, double radius, double length) 
     66{ 
     67    const double hdist = sqrt(square(radius_bell) - square(radius)); 
     68    const double length_tot = length + 2.0*(hdist+ radius); 
     69    return 0.5*cbrt(0.75*radius_bell*(2.0*radius_bell*length_tot + (radius_bell + length_tot)*(M_PI*radius_bell + length_tot))); 
     70} 
     71 
     72static double 
    6573radius_from_volume(double radius_bell, double radius, double length) 
    6674{ 
    … …  
    8189    switch (mode) { 
    8290    default: 
    83     case 1: // equivalent sphere 
     91    case 1: // equivalent cylinder excluded volume 
     92        return radius_from_excluded_volume(radius_bell, radius , length); 
     93    case 2: // equivalent volume sphere 
    8494        return radius_from_volume(radius_bell, radius , length); 
    85     case 2: // radius 
     95    case 3: // radius 
    8696        return radius; 
    87     case 3: // half length 
     97    case 4: // half length 
    8898        return 0.5*length; 
    89     case 4: // half total length 
     99    case 5: // half total length 
    90100        return radius_from_totallength(radius_bell,radius,length); 
    91101    } 

  • sasmodels/models/barbell.py

    ree60aa7 r99658f6  
    7979.. [#] H Kaya and N R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda 
    8080   and errata) 
     81L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    8182 
    8283Authorship and Verification 
    … …  
    116117 
    117118source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] 
    118 have_Fq = True 
     119have_Fq = True  
    119120effective_radius_type = [ 
    120     "equivalent sphere", "radius", "half length", "half total length", 
     121    "equivalent cylinder excluded volume","equivalent volume sphere", "radius", "half length", "half total length", 
    121122    ] 
    122123 

  • sasmodels/models/capped_cylinder.c

    rd42dd4a r99658f6  
    8585 
    8686static double 
     87radius_from_excluded_volume(double radius, double radius_cap, double length) 
     88{ 
     89    const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); 
     90    const double length_tot = length + 2.0*hc; 
     91    return 0.5*cbrt(0.75*radius*(2.0*radius*length_tot + (radius + length_tot)*(M_PI*radius + length_tot))); 
     92} 
     93 
     94static double 
    8795radius_from_volume(double radius, double radius_cap, double length) 
    8896{ 
    … …  
    103111    switch (mode) { 
    104112    default: 
    105     case 1: // equivalent sphere 
     113    case 1: // equivalent cylinder excluded volume 
     114        return radius_from_excluded_volume(radius, radius_cap, length); 
     115    case 2: // equivalent volume sphere 
    106116        return radius_from_volume(radius, radius_cap, length); 
    107     case 2: // radius 
     117    case 3: // radius 
    108118        return radius; 
    109     case 3: // half length 
     119    case 4: // half length 
    110120        return 0.5*length; 
    111     case 4: // half total length 
     121    case 5: // half total length 
    112122        return radius_from_totallength(radius, radius_cap,length); 
    113123    } 

  • sasmodels/models/capped_cylinder.py

    ree60aa7 r99658f6  
    8282.. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda 
    8383   and errata) 
     84L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    8485 
    8586Authorship and Verification 
    … …  
    138139have_Fq = True 
    139140effective_radius_type = [ 
    140     "equivalent sphere", "radius", "half length", "half total length", 
     141    "equivalent cylinder excluded volume", "equivalent volume sphere", "radius", "half length", "half total length", 
    141142    ] 
    142143 

  • sasmodels/models/core_shell_bicelle.c

    rd42dd4a r99658f6  
    3737 
    3838static double 
     39radius_from_excluded_volume(double radius, double thick_rim, double thick_face, double length) 
     40{ 
     41    const double radius_tot = radius + thick_rim; 
     42    const double length_tot = length + 2.0*thick_face; 
     43    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length_tot + (radius_tot + length_tot)*(M_PI*radius_tot + length_tot))); 
     44} 
     45 
     46static double 
    3947radius_from_volume(double radius, double thick_rim, double thick_face, double length) 
    4048{ 
    … …  
    5664    switch (mode) { 
    5765    default: 
    58     case 1: // equivalent sphere 
     66    case 1: // equivalent cylinder excluded volume 
     67        return radius_from_excluded_volume(radius, thick_rim, thick_face, length); 
     68    case 2: // equivalent sphere 
    5969        return radius_from_volume(radius, thick_rim, thick_face, length); 
    60     case 2: // outer rim radius 
     70    case 3: // outer rim radius 
    6171        return radius + thick_rim; 
    62     case 3: // half outer thickness 
     72    case 4: // half outer thickness 
    6373        return 0.5*length + thick_face; 
    64     case 4: // half diagonal 
     74    case 5: // half diagonal 
    6575        return radius_from_diagonal(radius,thick_rim,thick_face,length); 
    6676    } 

  • sasmodels/models/core_shell_bicelle.py

    ree60aa7 r99658f6  
    8989   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    9090   =26379>`_ 
     91    
     92   L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    9193 
    9294Authorship and Verification 
    … …  
    156158have_Fq = True 
    157159effective_radius_type = [ 
    158     "equivalent sphere", "outer rim radius", 
     160    "excluded volume","equivalent volume sphere", "outer rim radius", 
    159161    "half outer thickness", "half diagonal", 
    160162    ] 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    rd42dd4a r99658f6  
    88{ 
    99    return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face); 
     10} 
     11 
     12static double 
     13radius_from_excluded_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     14{ 
     15    const double r_equiv     = sqrt((r_minor + thick_rim)*(r_minor*x_core + thick_rim)); 
     16    const double length_tot  = length + 2.0*thick_face; 
     17    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_tot + (r_equiv + length_tot)*(M_PI*r_equiv + length_tot))); 
    1018} 
    1119 
    … …  
    3139    switch (mode) { 
    3240    default: 
    33     case 1: // equivalent sphere 
     41    case 1: // equivalent cylinder excluded volume 
     42        return radius_from_excluded_volume(r_minor, x_core, thick_rim, thick_face, length); 
     43    case 2: // equivalent volume sphere 
    3444        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
    35     case 2: // outer rim average radius 
     45    case 3: // outer rim average radius 
    3646        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
    37     case 3: // outer rim min radius 
     47    case 4: // outer rim min radius 
    3848        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
    39     case 4: // outer max radius 
     49    case 5: // outer max radius 
    4050        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
    41     case 5: // half outer thickness 
     51    case 6: // half outer thickness 
    4252        return 0.5*length + thick_face; 
    43     case 6: // half diagonal 
     53    case 7: // half diagonal 
    4454        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
    4555    } 

  • sasmodels/models/core_shell_bicelle_elliptical.py

    r304c775 r99658f6  
    100100 
    101101.. [#] 
     102L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    102103 
    103104Authorship and Verification 
    … …  
    148149have_Fq = True 
    149150effective_radius_type = [ 
    150     "equivalent sphere", "outer rim average radius", "outer rim min radius", 
     151    "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius", 
    151152    "outer max radius", "half outer thickness", "half diagonal", 
    152153    ] 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    rd42dd4a r99658f6  
    99    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 
    1010                 square(r_minor)*x_core*2.0*thick_face  ); 
     11} 
     12 
     13static double 
     14radius_from_excluded_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     15{ 
     16    const double r_equiv     = sqrt((r_minor + thick_rim)*(r_minor*x_core + thick_rim)); 
     17    const double length_tot  = length + 2.0*thick_face; 
     18    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_tot + (r_equiv + length_tot)*(M_PI*r_equiv + length_tot))); 
    1119} 
    1220 
    … …  
    3240    switch (mode) { 
    3341    default: 
    34     case 1: // equivalent sphere 
     42    case 1: // equivalent cylinder excluded volume 
     43        return radius_from_excluded_volume(r_minor, x_core, thick_rim, thick_face, length); 
     44    case 2: // equivalent sphere 
    3545        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
    36     case 2: // outer rim average radius 
     46    case 3: // outer rim average radius 
    3747        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
    38     case 3: // outer rim min radius 
     48    case 4: // outer rim min radius 
    3949        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
    40     case 4: // outer max radius 
     50    case 5: // outer max radius 
    4151        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
    42     case 5: // half outer thickness 
     52    case 6: // half outer thickness 
    4353        return 0.5*length + thick_face; 
    44     case 6: // half diagonal (this ignores the missing "corners", so may give unexpected answer if thick_face 
     54    case 7: // half diagonal (this ignores the missing "corners", so may give unexpected answer if thick_face 
    4555            //  or thick_rim is extremely large) 
    4656        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    r304c775 r99658f6  
    112112 
    113113.. [#] 
     114L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    114115 
    115116Authorship and Verification 
    … …  
    161162have_Fq = True 
    162163effective_radius_type = [ 
    163     "equivalent sphere", "outer rim average radius", "outer rim min radius", 
     164    "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius", 
    164165    "outer max radius", "half outer thickness", "half diagonal", 
    165166    ] 

  • sasmodels/models/core_shell_cylinder.c

    rd42dd4a r99658f6  
    1111{ 
    1212    return M_PI*square(radius+thickness)*(length+2.0*thickness); 
     13} 
     14 
     15static double 
     16radius_from_excluded_volume(double radius, double thickness, double length) 
     17{ 
     18    const double radius_tot = radius + thickness; 
     19    const double length_tot = length + 2.0*thickness; 
     20    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length_tot + (radius_tot + length_tot)*(M_PI*radius_tot + length_tot))); 
    1321} 
    1422 
    … …  
    3341    switch (mode) { 
    3442    default: 
    35     case 1: // equivalent sphere 
     43    case 1: //cylinder excluded volume 
     44        return radius_from_excluded_volume(radius, thickness, length); 
     45    case 2: // equivalent volume sphere 
    3646        return radius_from_volume(radius, thickness, length); 
    37     case 2: // outer radius 
     47    case 3: // outer radius 
    3848        return radius + thickness; 
    39     case 3: // half outer length 
     49    case 4: // half outer length 
    4050        return 0.5*length + thickness; 
    41     case 4: // half min outer length 
     51    case 5: // half min outer length 
    4252        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness); 
    43     case 5: // half max outer length 
     53    case 6: // half max outer length 
    4454        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness); 
    45     case 6: // half outer diagonal 
     55    case 7: // half outer diagonal 
    4656        return radius_from_diagonal(radius,thickness,length); 
    4757    } 

  • sasmodels/models/core_shell_cylinder.py

    ree60aa7 r99658f6  
    7070   1445-1452 
    7171.. [#kline] S R Kline, *J Appl. Cryst.*, 39 (2006) 895 
     72L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    7273 
    7374Authorship and Verification 
    … …  
    133134have_Fq = True 
    134135effective_radius_type = [ 
    135     "equivalent sphere", "outer radius", "half outer length", 
    136     "half min outer dimension", "half max outer dimension", 
    137     "half outer diagonal", 
     136    "excluded volume", "equivalent volume sphere", "outer radius", "half outer length", 
     137    "half min outer dimension", "half max outer dimension", "half outer diagonal", 
    138138    ] 
    139139 

  • sasmodels/models/core_shell_ellipsoid.c

    rd42dd4a r99658f6  
    7575    switch (mode) { 
    7676    default: 
    77     case 1: // equivalent sphere 
     77    case 1: // average outer curvature 
     78        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     79    case 2: // equivalent volume sphere 
    7880        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); 
    8181    case 3: // min outer radius 
    8282        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot); 

  • sasmodels/models/core_shell_ellipsoid.py

    ree60aa7 r99658f6  
    147147have_Fq = True 
    148148effective_radius_type = [ 
    149     "equivalent sphere", "average outer curvature", 
     149    "average outer curvature", "equivalent volume sphere",      
    150150    "min outer radius", "max outer radius", 
    151151    ] 

  • sasmodels/models/core_shell_parallelepiped.c

    rd42dd4a r99658f6  
    2828 
    2929static double 
     30radius_from_excluded_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    double r_equiv, length; 
     34    double lengths[3] = {length_a+thick_rim_a, length_b+thick_rim_b, length_c+thick_rim_c}; 
     35    double lengthmax = fmax(lengths[0],fmax(lengths[1],lengths[2])); 
     36    double length_1 = lengthmax; 
     37    double length_2 = lengthmax; 
     38    double length_3 = lengthmax; 
     39 
     40    for(int ilen=0; ilen<3; ilen++) { 
     41        if (lengths[ilen] < length_1) { 
     42            length_2 = length_1; 
     43            length_1 = lengths[ilen]; 
     44            } else { 
     45                if (lengths[ilen] < length_2) { 
     46                        length_2 = lengths[ilen]; 
     47                } 
     48            } 
     49    } 
     50    if(length_2-length_1 > length_3-length_2) { 
     51        r_equiv = sqrt(length_2*length_3/M_PI); 
     52        length  = length_1; 
     53    } else  { 
     54        r_equiv = sqrt(length_1*length_2/M_PI); 
     55        length  = length_3; 
     56    } 
     57 
     58    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length))); 
     59} 
     60 
     61static double 
    3062radius_from_volume(double length_a, double length_b, double length_c, 
    3163                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    … …  
    4880    switch (mode) { 
    4981    default: 
    50     case 1: // equivalent sphere 
     82    case 1: // equivalent cylinder excluded volume 
     83        return radius_from_excluded_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     84    case 2: // equivalent volume sphere 
    5185        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 
     86    case 3: // half outer length a 
    5387        return 0.5 * length_a + thick_rim_a; 
    54     case 3: // half outer length b 
     88    case 4: // half outer length b 
    5589        return 0.5 * length_b + thick_rim_b; 
    56     case 4: // half outer length c 
     90    case 5: // half outer length c 
    5791        return 0.5 * length_c + thick_rim_c; 
    58     case 5: // equivalent circular cross-section 
     92    case 6: // equivalent circular cross-section 
    5993        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b); 
    60     case 6: // half outer ab diagonal 
     94    case 7: // half outer ab diagonal 
    6195        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 
     96    case 8: // half outer diagonal 
    6397        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)); 
    6498    } 

  • sasmodels/models/core_shell_parallelepiped.py

    ree60aa7 r99658f6  
    173173   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    174174   =26379>`_ 
     175L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    175176 
    176177Authorship and Verification 
    … …  
    228229have_Fq = True 
    229230effective_radius_type = [ 
    230     "equivalent sphere", 
     231    "equivalent cylinder excluded volume",  
     232    "equivalent volume sphere", 
    231233    "half outer length_a", "half outer length_b", "half outer length_c", 
    232234    "equivalent circular cross-section", 

  • sasmodels/models/cylinder.c

    rd42dd4a r99658f6  
    1111{ 
    1212    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
     13} 
     14 
     15static double 
     16radius_from_excluded_volume(double radius, double length) 
     17{ 
     18    return 0.5*cbrt(0.75*radius*(2.0*radius*length + (radius + length)*(M_PI*radius + length))); 
    1319} 
    1420 
    … …  
    3137    default: 
    3238    case 1: 
     39        return radius_from_excluded_volume(radius, length); 
     40    case 2: 
    3341        return radius_from_volume(radius, length); 
    34     case 2: 
     42    case 3: 
    3543        return radius; 
    36     case 3: 
     44    case 4: 
    3745        return 0.5*length; 
    38     case 4: 
     46    case 5: 
    3947        return (radius < 0.5*length ? radius : 0.5*length); 
    40     case 5: 
     48    case 6: 
    4149        return (radius > 0.5*length ? radius : 0.5*length); 
    42     case 6: 
     50    case 7: 
    4351        return radius_from_diagonal(radius,length); 
    4452    } 

  • sasmodels/models/cylinder.py

    r304c775 r99658f6  
    9898J. S. Pedersen, Adv. Colloid Interface Sci. 70, 171-210 (1997). 
    9999G. Fournet, Bull. Soc. Fr. Mineral. Cristallogr. 74, 39-113 (1951). 
     100L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    100101""" 
    101102 
    … …  
    140141have_Fq = True 
    141142effective_radius_type = [ 
    142     "equivalent sphere", "radius", 
     143    "excluded volume", "equivalent volume sphere", "radius", 
    143144    "half length", "half min dimension", "half max dimension", "half diagonal", 
    144145    ] 
    … …  
    185186radius, length = parameters[2][2], parameters[3][2] 
    186187tests.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), 
     188    ({'radius_effective_mode': 0}, 0.1, None, None, 0., pi*radius**2*length, 1.0),    
     189    ({'radius_effective_mode': 1}, 0.1, None, None, 0.5*(0.75*radius*(2.0*radius*length + (radius + length)*(pi*radius + length)))**(1./3.), None, None),     
     190    ({'radius_effective_mode': 2}, 0.1, None, None, (0.75*radius**2*length)**(1./3.), None, None), 
     191    ({'radius_effective_mode': 3}, 0.1, None, None, radius, None, None), 
     192    ({'radius_effective_mode': 4}, 0.1, None, None, length/2., None, None), 
     193    ({'radius_effective_mode': 5}, 0.1, None, None, min(radius, length/2.), None, None), 
     194    ({'radius_effective_mode': 6}, 0.1, None, None, max(radius, length/2.), None, None), 
     195    ({'radius_effective_mode': 7}, 0.1, None, None, np.sqrt(4*radius**2 + length**2)/2., None, None), 
    194196]) 
    195197del radius, length 

  • sasmodels/models/ellipsoid.c

    rd42dd4a r99658f6  
    3636    switch (mode) { 
    3737    default: 
    38     case 1: // equivalent sphere 
     38    case 1: // average curvature 
     39        return radius_from_curvature(radius_polar, radius_equatorial); 
     40    case 2: // equivalent volume sphere 
    3941        return radius_from_volume(radius_polar, radius_equatorial); 
    40     case 2: // average curvature 
    41         return radius_from_curvature(radius_polar, radius_equatorial); 
    4242    case 3: // min radius 
    4343        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial); 

  • sasmodels/models/ellipsoid.py

    ree60aa7 r99658f6  
    170170have_Fq = True 
    171171effective_radius_type = [ 
    172     "equivalent sphere", "average curvature", "min radius", "max radius", 
     172    "average curvature", "equivalent volume sphere", "min radius", "max radius", 
    173173    ] 
    174174 

  • sasmodels/models/elliptical_cylinder.c

    rd42dd4a r99658f6  
    33{ 
    44    return M_PI * radius_minor * radius_minor * r_ratio * length; 
     5} 
     6 
     7static double 
     8radius_from_excluded_volume(double radius_minor, double r_ratio, double length) 
     9{ 
     10    const double r_equiv = sqrt(radius_minor*radius_minor*r_ratio); 
     11    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length))); 
    512} 
    613 
    … …  
    3845    switch (mode) { 
    3946    default: 
    40     case 1: // equivalent sphere 
     47    case 1: // equivalent cylinder excluded volume 
     48        return radius_from_excluded_volume(radius_minor, r_ratio, length); 
     49    case 2: // equivalent volume sphere 
    4150        return radius_from_volume(radius_minor, r_ratio, length); 
    42     case 2: // average radius 
     51    case 3: // average radius 
    4352        return 0.5*radius_minor*(1.0 + r_ratio); 
    44     case 3: // min radius 
     53    case 4: // min radius 
    4554        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
    46     case 4: // max radius 
     55    case 5: // max radius 
    4756        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
    48     case 5: // equivalent circular cross-section 
     57    case 6: // equivalent circular cross-section 
    4958        return sqrt(radius_minor*radius_minor*r_ratio); 
    50     case 6: // half length 
     59    case 7: // half length 
    5160        return 0.5*length; 
    52     case 7: // half min dimension 
     61    case 8: // half min dimension 
    5362        return radius_from_min_dimension(radius_minor,r_ratio,0.5*length); 
    54     case 8: // half max dimension 
     63    case 9: // half max dimension 
    5564        return radius_from_max_dimension(radius_minor,r_ratio,0.5*length); 
    56     case 9: // half diagonal 
     65    case 10: // half diagonal 
    5766        return radius_from_diagonal(radius_minor,r_ratio,length); 
    5867    } 

  • sasmodels/models/elliptical_cylinder.py

    ree60aa7 r99658f6  
    8888L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 
    8989Neutron Scattering*, Plenum, New York, (1987) [see table 3.4] 
     90L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    9091 
    9192Authorship and Verification 
    … …  
    124125have_Fq = True 
    125126effective_radius_type = [ 
    126     "equivalent sphere", "average radius", "min radius", "max radius", 
     127    "equivalent cylinder excluded volume", "equivalent volume sphere", "average radius", "min radius", "max radius", 
    127128    "equivalent circular cross-section", 
    128129    "half length", "half min dimension", "half max dimension", "half diagonal", 

  • sasmodels/models/hardsphere.py

    r304c775 r304c775  
    162162    return pars 
    163163 
    164 demo = dict(radius_effective=200, volfraction=0.2, 
    165             radius_effective_pd=0.1, radius_effective_pd_n=40) 
    166164# Q=0.001 is in the Taylor series, low Q part, so add Q=0.1, 
    167165# assuming double precision sasview is correct 

  • sasmodels/models/hollow_cylinder.c

    rd42dd4a r99658f6  
    1111{ 
    1212    return M_PI*length*square(radius+thickness); 
     13} 
     14 
     15static double 
     16radius_from_excluded_volume(double radius, double thickness, double length) 
     17{ 
     18    const double radius_tot = radius + thickness; 
     19    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length + (radius_tot + length)*(M_PI*radius_tot + length))); 
    1320} 
    1421 
    … …  
    3138    switch (mode) { 
    3239    default: 
    33     case 1: // equivalent sphere 
     40    case 1: // excluded volume 
     41        return radius_from_excluded_volume(radius, thickness, length); 
     42    case 2: // equivalent volume sphere 
    3443        return radius_from_volume(radius, thickness, length); 
    35     case 2: // outer radius 
     44    case 3: // outer radius 
    3645        return radius + thickness; 
    37     case 3: // half length 
     46    case 4: // half length 
    3847        return 0.5*length; 
    39     case 4: // half outer min dimension 
     48    case 5: // half outer min dimension 
    4049        return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length); 
    41     case 5: // half outer max dimension 
     50    case 6: // half outer max dimension 
    4251        return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length); 
    43     case 6: // half outer diagonal 
     52    case 7: // half outer diagonal 
    4453        return radius_from_diagonal(radius,thickness,length); 
    4554    } 

  • sasmodels/models/hollow_cylinder.py

    r304c775 r99658f6  
    6060.. [#] L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 
    6161   Neutron Scattering*, Plenum Press, New York, (1987) 
     62L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    6263 
    6364Authorship and Verification 
    … …  
    102103have_Fq = True 
    103104effective_radius_type = [ 
    104     "equivalent sphere", "outer radius", "half length", 
     105    "excluded volume", "equivalent outer volume sphere", "outer radius", "half length", 
    105106    "half outer min dimension", "half outer max dimension", 
    106107    "half outer diagonal", 
    … …  
    140141    [{}, 0.00005, 1764.926], 
    141142    [{}, 0.1, None, None, 
    142      (3./4*(radius+thickness)**2*length)**(1./3),  # R_eff from volume 
     143     0.5*(0.75*(radius+thickness)*(2.0*(radius+thickness)*length + ((radius+thickness) + length)*(pi*(radius+thickness) + length)))**(1./3.),  # R_eff from excluded volume 
    143144     pi*((radius+thickness)**2-radius**2)*length,  # shell volume 
    144145     (radius+thickness)**2/((radius+thickness)**2 - radius**2), # form:shell ratio 

  • sasmodels/models/hollow_rectangular_prism.c

    rd42dd4a r99658f6  
    2222 
    2323static double 
     24radius_from_excluded_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     25{ 
     26    const double r_equiv = sqrt(length_a*length_a*b2a_ratio/M_PI); 
     27    const double length_c = length_a*c2a_ratio; 
     28    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_c + (r_equiv + length_c)*(M_PI*r_equiv + length_c))); 
     29} 
     30 
     31static double 
    2432effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    2533// NOTE length_a is external dimension! 
    … …  
    2735    switch (mode) { 
    2836    default: 
    29     case 1: // equivalent sphere 
     37    case 1: // equivalent cylinder excluded volume 
     38        return radius_from_excluded_volume(length_a, b2a_ratio, c2a_ratio); 
     39    case 2: // equivalent outer volume sphere 
    3040        return cbrt(cube(length_a)*b2a_ratio*c2a_ratio/M_4PI_3); 
    31     case 2: // half length_a 
     41    case 3: // half length_a 
    3242        return 0.5 * length_a; 
    33     case 3: // half length_b 
     43    case 4: // half length_b 
    3444        return 0.5 * length_a*b2a_ratio; 
    35     case 4: // half length_c 
     45    case 5: // half length_c 
    3646        return 0.5 * length_a*c2a_ratio; 
    37     case 5: // equivalent outer circular cross-section 
     47    case 6: // equivalent outer circular cross-section 
    3848        return length_a*sqrt(b2a_ratio/M_PI); 
    39     case 6: // half ab diagonal 
     49    case 7: // half ab diagonal 
    4050        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
    41     case 7: // half diagonal 
     51    case 8: // half diagonal 
    4252        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
    4353    } 

  • sasmodels/models/hollow_rectangular_prism.py

    ree60aa7 r99658f6  
    9898 
    9999.. [#Nayuk2012] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     100L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    100101 
    101102 
    … …  
    150151have_Fq = True 
    151152effective_radius_type = [ 
    152     "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     153    "equivalent cylinder excluded volume", "equivalent outer volume sphere",  
     154    "half length_a", "half length_b", "half length_c", 
    153155    "equivalent outer circular cross-section", 
    154156    "half ab diagonal", "half diagonal", 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    rd42dd4a r99658f6  
    1818 
    1919static double 
     20radius_from_excluded_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     21{ 
     22    const double r_equiv = sqrt(length_a*length_a*b2a_ratio/M_PI); 
     23    const double length_c = length_a*c2a_ratio; 
     24    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_c + (r_equiv + length_c)*(M_PI*r_equiv + length_c))); 
     25} 
     26 
     27static double 
    2028effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio) 
    2129{ 
    2230    switch (mode) { 
    2331    default: 
    24     case 1: // equivalent sphere 
     32    case 1: // equivalent cylinder excluded volume 
     33        return radius_from_excluded_volume(length_a, b2a_ratio, c2a_ratio); 
     34    case 2: // equivalent outer volume sphere 
    2535        return cbrt(cube(length_a)*b2a_ratio*c2a_ratio/M_4PI_3); 
    26     case 2: // half length_a 
     36    case 3: // half length_a 
    2737        return 0.5 * length_a; 
    28     case 3: // half length_b 
     38    case 4: // half length_b 
    2939        return 0.5 * length_a*b2a_ratio; 
    30     case 4: // half length_c 
     40    case 5: // half length_c 
    3141        return 0.5 * length_a*c2a_ratio; 
    32     case 5: // equivalent outer circular cross-section 
     42    case 6: // equivalent outer circular cross-section 
    3343        return length_a*sqrt(b2a_ratio/M_PI); 
    34     case 6: // half ab diagonal 
     44    case 7: // half ab diagonal 
    3545        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
    36     case 7: // half diagonal 
     46    case 8: // half diagonal 
    3747        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
    3848    } 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    ree60aa7 r99658f6  
    7272 
    7373.. [#Nayuk2012] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     74L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    7475 
    7576 
    … …  
    110111have_Fq = True 
    111112effective_radius_type = [ 
    112     "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     113    "equivalent cylinder excluded volume", "equivalent outer volume sphere",  
     114    "half length_a", "half length_b", "half length_c", 
    113115    "equivalent outer circular cross-section", 
    114116    "half ab diagonal", "half diagonal", 

  • sasmodels/models/parallelepiped.c

    rd42dd4a r99658f6  
    66 
    77static double 
     8radius_from_excluded_volume(double length_a, double length_b, double length_c) 
     9{ 
     10    double r_equiv, length; 
     11    double lengths[3] = {length_a, length_b, length_c}; 
     12    double lengthmax = fmax(lengths[0],fmax(lengths[1],lengths[2])); 
     13    double length_1 = lengthmax; 
     14    double length_2 = lengthmax; 
     15    double length_3 = lengthmax; 
     16 
     17    for(int ilen=0; ilen<3; ilen++) { 
     18        if (lengths[ilen] < length_1) { 
     19            length_2 = length_1; 
     20            length_1 = lengths[ilen]; 
     21            } else { 
     22                if (lengths[ilen] < length_2) { 
     23                        length_2 = lengths[ilen]; 
     24                } 
     25            } 
     26    } 
     27    if(length_2-length_1 > length_3-length_2) { 
     28        r_equiv = sqrt(length_2*length_3/M_PI); 
     29        length  = length_1; 
     30    } else  { 
     31        r_equiv = sqrt(length_1*length_2/M_PI); 
     32        length  = length_3; 
     33    } 
     34 
     35    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length))); 
     36} 
     37 
     38static double 
    839effective_radius(int mode, double length_a, double length_b, double length_c) 
    940{ 
    1041    switch (mode) { 
    1142    default: 
    12     case 1: // equivalent sphere 
     43    case 1: // equivalent cylinder excluded volume 
     44        return radius_from_excluded_volume(length_a,length_b,length_c); 
     45    case 2: // equivalent volume sphere 
    1346        return cbrt(length_a*length_b*length_c/M_4PI_3); 
    14     case 2: // half length_a 
     47    case 3: // half length_a 
    1548        return 0.5 * length_a; 
    16     case 3: // half length_b 
     49    case 4: // half length_b 
    1750        return 0.5 * length_b; 
    18     case 4: // half length_c 
     51    case 5: // half length_c 
    1952        return 0.5 * length_c; 
    20     case 5: // equivalent circular cross-section 
     53    case 6: // equivalent circular cross-section 
    2154        return sqrt(length_a*length_b/M_PI); 
    22     case 6: // half ab diagonal 
     55    case 7: // half ab diagonal 
    2356        return 0.5*sqrt(length_a*length_a + length_b*length_b); 
    24     case 7: // half diagonal 
     57    case 8: // half diagonal 
    2558        return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c); 
    2659    } 

  • sasmodels/models/parallelepiped.py

    ree60aa7 r99658f6  
    180180   14 (1961) 185-211 
    181181.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     182L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    182183 
    183184Authorship and Verification 
    … …  
    232233have_Fq = True 
    233234effective_radius_type = [ 
    234     "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     235    "equivalent cylinder excluded volume", "equivalent volume sphere",  
     236    "half length_a", "half length_b", "half length_c", 
    235237    "equivalent circular cross-section", "half ab diagonal", "half diagonal", 
    236238    ] 

  • sasmodels/models/pearl_necklace.c

    r3f853beb r9b5fd42  
    4040    const double si = sas_sinx_x(q*A_s); 
    4141    const double omsi = 1.0 - si; 
    42     const double pow_si = pow(si, num_pearls); 
     42    const double pow_si = pown(si, num_pearls); 
    4343 
    4444    // form factor for num_pearls 
    … …  
    6767} 
    6868 
    69 double form_volume(double radius, double edge_sep, 
    70     double thick_string, double fp_num_pearls) 
     69double form_volume(double radius, double edge_sep, double thick_string, double fp_num_pearls) 
    7170{ 
    7271    const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls 
    … …  
    7776 
    7877    return volume; 
     78} 
     79 
     80static double 
     81radius_from_volume(double radius, double edge_sep, double thick_string, double fp_num_pearls) 
     82{ 
     83    const double vol_tot = form_volume(radius, edge_sep, thick_string, fp_num_pearls); 
     84    return cbrt(vol_tot/M_4PI_3); 
     85} 
     86 
     87static double 
     88effective_radius(int mode, double radius, double edge_sep, double thick_string, double fp_num_pearls) 
     89{ 
     90    switch (mode) { 
     91    default: 
     92    case 1: 
     93        return radius_from_volume(radius, edge_sep, thick_string, fp_num_pearls); 
     94    } 
    7995} 
    8096 

  • sasmodels/models/pearl_necklace.py

    r2cc8aa2 rcf3d0ce  
    5353R Schweins and K Huber, *Particle Scattering Factor of Pearl Necklace Chains*, 
    5454*Macromol. Symp.* 211 (2004) 25-42 2004 
     55L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    5556""" 
    5657 
    … …  
    9596source = ["lib/sas_Si.c", "lib/sas_3j1x_x.c", "pearl_necklace.c"] 
    9697single = False  # use double precision unless told otherwise 
    97  
    98 def volume(radius, edge_sep, thick_string, num_pearls): 
    99     """ 
    100     Calculates the total particle volume of the necklace. 
    101     Redundant with form_volume. 
    102     """ 
    103     num_pearls = int(num_pearls + 0.5) 
    104     number_of_strings = num_pearls - 1.0 
    105     string_vol = edge_sep * pi * pow((thick_string / 2.0), 2.0) 
    106     pearl_vol = 4.0 /3.0 * pi * pow(radius, 3.0) 
    107     total_vol = number_of_strings * string_vol 
    108     total_vol += num_pearls * pearl_vol 
    109     return total_vol 
    110  
    111 def ER(radius, edge_sep, thick_string, num_pearls): 
    112     """ 
    113     Calculation for effective radius. 
    114     """ 
    115     num_pearls = int(num_pearls + 0.5) 
    116     tot_vol = volume(radius, edge_sep, thick_string, num_pearls) 
    117     rad_out = (tot_vol/(4.0/3.0*pi)) ** (1./3.) 
    118     return rad_out 
    119  
     98effective_radius_type = [ 
     99    "equivalent volume sphere",  
     100    ] 
     101     
    120102def random(): 
    121103    radius = 10**np.random.uniform(1, 3) # 1 - 1000 

  • sasmodels/models/pringle.c

    rd42dd4a r99658f6  
    105105 
    106106static double 
     107radius_from_excluded_volume(double radius, double thickness) 
     108{ 
     109    return 0.5*cbrt(0.75*radius*(2.0*radius*thickness + (radius + thickness)*(M_PI*radius + thickness))); 
     110} 
     111 
     112static double 
    107113effective_radius(int mode, double radius, double thickness, double alpha, double beta) 
    108114{ 
    109115    switch (mode) { 
    110116    default: 
    111     case 1: // equivalent sphere 
     117    case 1: // equivalent cylinder excluded volume 
     118        return radius_from_excluded_volume(radius, thickness); 
     119    case 2: // equivalent volume sphere 
    112120        return cbrt(M_PI*radius*radius*thickness/M_4PI_3); 
    113     case 2: // radius 
     121    case 3: // radius 
    114122        return radius; 
    115123    } 

  • sasmodels/models/pringle.py

    ree60aa7 r99658f6  
    4242Karen Edler, Universtiy of Bath, Private Communication. 2012. 
    4343Derivation by Stefan Alexandru Rautu. 
     44L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949). 
    4445 
    4546* **Author:** Andrew Jackson **Date:** 2008 
    … …  
    7475source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", 
    7576          "lib/sas_JN.c", "lib/gauss76.c", "pringle.c"] 
    76 effective_radius_type = ["equivalent sphere", "radius"] 
     77effective_radius_type = ["equivalent cylinder excluded volume", "equivalent volume sphere", "radius"] 
    7778 
    7879def random(): 

  • sasmodels/models/rectangular_prism.c

    rd42dd4a r99658f6  
    66 
    77static double 
     8radius_from_excluded_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     9{ 
     10    double const r_equiv   = sqrt(length_a*length_a*b2a_ratio/M_PI); 
     11    double const length_c  = c2a_ratio*length_a; 
     12    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_c + (r_equiv + length_c)*(M_PI*r_equiv + length_c))); 
     13} 
     14 
     15static double 
    816effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio) 
    917{ 
    1018    switch (mode) { 
    1119    default: 
    12     case 1: // equivalent sphere 
     20    case 1: // equivalent cylinder excluded volume 
     21        return radius_from_excluded_volume(length_a,b2a_ratio,c2a_ratio); 
     22    case 2: // equivalent volume sphere 
    1323        return cbrt(cube(length_a)*b2a_ratio*c2a_ratio/M_4PI_3); 
    14     case 2: // half length_a 
     24    case 3: // half length_a 
    1525        return 0.5 * length_a; 
    16     case 3: // half length_b 
     26    case 4: // half length_b 
    1727        return 0.5 * length_a*b2a_ratio; 
    18     case 4: // half length_c 
     28    case 5: // half length_c 
    1929        return 0.5 * length_a*c2a_ratio; 
    20     case 5: // equivalent circular cross-section 
     30    case 6: // equivalent circular cross-section 
    2131        return length_a*sqrt(b2a_ratio/M_PI); 
    22     case 6: // half ab diagonal 
     32    case 7: // half ab diagonal 
    2333        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
    24     case 7: // half diagonal 
     34    case 8: // half diagonal 
    2535        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
    2636    } 

  • sasmodels/models/rectangular_prism.py

    ree60aa7 r99658f6  
    9999 
    100100R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     101 
     102L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
     103 
    101104""" 
    102105 
    … …  
    137140have_Fq = True 
    138141effective_radius_type = [ 
    139     "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     142    "equivalent cylinder excluded volume", "equivalent volume sphere",  
     143    "half length_a", "half length_b", "half length_c", 
    140144    "equivalent circular cross-section", "half ab diagonal", "half diagonal", 
    141145    ] 

  • sasmodels/models/triaxial_ellipsoid.c

    rd42dd4a r99658f6  
    55{ 
    66    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
     7} 
     8 
     9static double 
     10radius_from_curvature(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     11{ 
     12    // Trivial cases 
     13    if (radius_equat_minor == radius_equat_major == radius_polar) return radius_polar; 
     14    if (radius_equat_minor * radius_equat_major * radius_polar == 0.)  return 0.; 
     15 
     16 
     17    double r_equat_equiv, r_polar_equiv; 
     18    double radii[3] = {radius_equat_minor, radius_equat_major, radius_polar}; 
     19    double radmax = fmax(radii[0],fmax(radii[1],radii[2])); 
     20 
     21    double radius_1 = radmax; 
     22    double radius_2 = radmax; 
     23    double radius_3 = radmax; 
     24 
     25    for(int irad=0; irad<3; irad++) { 
     26        if (radii[irad] < radius_1) { 
     27            radius_3 = radius_2; 
     28            radius_2 = radius_1; 
     29            radius_1 = radii[irad]; 
     30            } else { 
     31                if (radii[irad] < radius_2) { 
     32                        radius_2 = radii[irad]; 
     33                } 
     34            } 
     35    } 
     36    if(radius_2-radius_1 > radius_3-radius_2) { 
     37        r_equat_equiv = sqrt(radius_2*radius_3); 
     38        r_polar_equiv = radius_1; 
     39    } else  { 
     40        r_equat_equiv = sqrt(radius_1*radius_2); 
     41        r_polar_equiv = radius_3; 
     42    } 
     43 
     44    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     45    const double ratio = (r_polar_equiv < r_equat_equiv 
     46                          ? r_polar_equiv / r_equat_equiv 
     47                          : r_equat_equiv / r_polar_equiv); 
     48    const double e1 = sqrt(1.0 - ratio*ratio); 
     49    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     50    const double bL = (1.0 + e1) / (1.0 - e1); 
     51    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     52    const double delta = 0.75 * b1 * b2; 
     53    const double ddd = 2.0 * (delta + 1.0) * r_polar_equiv * r_equat_equiv * r_equat_equiv; 
     54    return 0.5 * cbrt(ddd); 
    755} 
    856 
    … …  
    3280    switch (mode) { 
    3381    default: 
    34     case 1: // equivalent sphere 
     82    case 1: // equivalent biaxial ellipsoid average curvature 
     83        return radius_from_curvature(radius_equat_minor,radius_equat_major, radius_polar); 
     84    case 2: // equivalent volume sphere 
    3585        return radius_from_volume(radius_equat_minor,radius_equat_major, radius_polar); 
    36     case 2: // min radius 
     86    case 3: // min radius 
    3787        return radius_from_min_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
    38     case 3: // max radius 
     88    case 4: // max radius 
    3989        return radius_from_max_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
    4090    } 

  • sasmodels/models/triaxial_ellipsoid.py

    ree60aa7 r99658f6  
    158158source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"] 
    159159have_Fq = True 
    160 effective_radius_type = ["equivalent sphere", "min radius", "max radius"] 
     160effective_radius_type = ["equivalent biaxial ellipsoid average curvature", "equivalent volume sphere", "min radius", "max radius"] 
    161161 
    162162def random(): 

  • sasmodels/product.py

    r39a06c9 r99658f6  
    3737#] 
    3838 
     39STRUCTURE_MODE_ID = "structure_factor_mode" 
     40RADIUS_MODE_ID = "radius_effective_mode" 
    3941RADIUS_ID = "radius_effective" 
    4042VOLFRAC_ID = "volfraction" 
    … …  
    4345    if p_info.have_Fq: 
    4446        par = parse_parameter( 
    45                 "structure_factor_mode", 
     47                STRUCTURE_MODE_ID, 
    4648                "", 
    4749                0, 
    … …  
    5254    if p_info.effective_radius_type is not None: 
    5355        par = parse_parameter( 
    54                 "radius_effective_mode", 
     56                RADIUS_MODE_ID, 
    5557                "", 
    56                 0, 
     58                1, 
    5759                [["unconstrained"] + p_info.effective_radius_type], 
    5860                "", 

  • sasmodels/resolution.py

    r9e7837a re2592f0  
    445445    q = np.sort(q) 
    446446    if q_min + 2*MINIMUM_RESOLUTION < q[0]: 
    447         n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15 
     447        n_low = int(np.ceil((q[0]-q_min) / (q[1]-q[0]))) if q[1] > q[0] else 15 
    448448        q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 
    449449    else: 
    450450        q_low = [] 
    451451    if q_max - 2*MINIMUM_RESOLUTION > q[-1]: 
    452         n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15 
     452        n_high = int(np.ceil((q_max-q[-1]) / (q[-1]-q[-2]))) if q[-1] > q[-2] else 15 
    453453        q_high = np.linspace(q[-1], q_max, n_high+1)[1:] 
    454454    else: 
    … …  
    499499            q_min = q[0]*MINIMUM_ABSOLUTE_Q 
    500500        n_low = log_delta_q * (log(q[0])-log(q_min)) 
    501         q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 
     501        q_low = np.logspace(log10(q_min), log10(q[0]), int(np.ceil(n_low))+1)[:-1] 
    502502    else: 
    503503        q_low = [] 
    504504    if q_max > q[-1]: 
    505505        n_high = log_delta_q * (log(q_max)-log(q[-1])) 
    506         q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:] 
     506        q_high = np.logspace(log10(q[-1]), log10(q_max), int(np.ceil(n_high))+1)[1:] 
    507507    else: 
    508508        q_high = [] 

  • sasmodels/sasview_model.py

    r39a06c9 r3a1afed  
    382382            hidden.add('scale') 
    383383            hidden.add('background') 
    384             self._model_info.parameters.defaults['background'] = 0. 
    385384 
    386385        # Update the parameter lists to exclude any hidden parameters 
    … …  
    695694            return self._calculate_Iq(qx, qy) 
    696695 
    697     def _calculate_Iq(self, qx, qy=None, Fq=False, effective_radius_type=1): 
     696    def _calculate_Iq(self, qx, qy=None): 
    698697        if self._model is None: 
    699698            self._model = core.build_model(self._model_info) 
    … …  
    715714        #print("values", values) 
    716715        #print("is_mag", is_magnetic) 
    717         if Fq: 
    718             result = calculator.Fq(call_details, values, cutoff=self.cutoff, 
    719                                    magnetic=is_magnetic, 
    720                                    effective_radius_type=effective_radius_type) 
    721716        result = calculator(call_details, values, cutoff=self.cutoff, 
    722717                            magnetic=is_magnetic) 
    … …  
    736731        Calculate the effective radius for P(q)*S(q) 
    737732 
     733        *mode* is the R_eff type, which defaults to 1 to match the ER 
     734        calculation for sasview models from version 3.x. 
     735 
    738736        :return: the value of the effective radius 
    739737        """ 
    740         Fq = self._calculate_Iq([0.1], True, mode) 
    741         return Fq[2] 
     738        # ER and VR are only needed for old multiplication models, based on 
     739        # sas.sascalc.fit.MultiplicationModel.  Fail for now.  If we want to 
     740        # continue supporting them then add some test cases so that the code 
     741        # is exercised.  We can access ER/VR using the kernel Fq function by 
     742        # extending _calculate_Iq so that it calls: 
     743        #    if er_mode > 0: 
     744        #        res = calculator.Fq(call_details, values, cutoff=self.cutoff, 
     745        #                            magnetic=False, effective_radius_type=mode) 
     746        #        R_eff, form_shell_ratio = res[2], res[4] 
     747        #        return R_eff, form_shell_ratio 
     748        # Then use the following in calculate_ER: 
     749        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=mode) 
     750        #    return ER 
     751        # Similarly, for calculate_VR: 
     752        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=1) 
     753        #    return VR 
     754        # Obviously a combined calculate_ER_VR method would be better, but 
     755        # we only need them to support very old models, so ignore the 2x 
     756        # performance hit. 
     757        raise NotImplementedError("ER function is no longer available.") 
    742758 
    743759    def calculate_VR(self): 
    … …  
    748764        :return: the value of the form:shell volume ratio 
    749765        """ 
    750         Fq = self._calculate_Iq([0.1], True, mode) 
    751         return Fq[4] 
     766        # See comments in calculate_ER. 
     767        raise NotImplementedError("VR function is no longer available.") 
    752768 
    753769    def set_dispersion(self, parameter, dispersion): 
    … …  
    819835            return value, [value], [1.0] 
    820836 
     837    @classmethod 
     838    def runTests(cls): 
     839        """ 
     840        Run any tests built into the model and captures the test output. 
     841 
     842        Returns success flag and output 
     843        """ 
     844        from .model_test import check_model 
     845        return check_model(cls._model_info) 
     846 
    821847def test_cylinder(): 
    822848    # type: () -> float 
    … …  
    849875    P = _make_standard_model('cylinder')() 
    850876    model = MultiplicationModel(P, S) 
    851     model.setParam('radius_effective_mode', 1.0) 
     877    model.setParam(product.RADIUS_MODE_ID, 1.0) 
    852878    value = model.evalDistribution([0.1, 0.1]) 
    853879    if np.isnan(value): 
    … …  
    904930    CylinderModel().evalDistribution([0.1, 0.1]) 
    905931 
     932def test_structure_factor_background(): 
     933    # type: () -> None 
     934    """ 
     935    Check that sasview model and direct model match, with background=0. 
     936    """ 
     937    from .data import empty_data1D 
     938    from .core import load_model_info, build_model 
     939    from .direct_model import DirectModel 
     940 
     941    model_name = "hardsphere" 
     942    q = [0.0] 
     943 
     944    sasview_model = _make_standard_model(model_name)() 
     945    sasview_value = sasview_model.evalDistribution(np.array(q))[0] 
     946 
     947    data = empty_data1D(q) 
     948    model_info = load_model_info(model_name) 
     949    model = build_model(model_info) 
     950    direct_model = DirectModel(data, model) 
     951    direct_value_zero_background = direct_model(background=0.0) 
     952 
     953    assert sasview_value == direct_value_zero_background 
     954 
     955    # Additionally check that direct value background defaults to zero 
     956    direct_value_default = direct_model() 
     957    assert sasview_value == direct_value_default 
     958 
     959 
    906960def magnetic_demo(): 
    907961    Model = _make_standard_model('sphere') 
    … …  
    924978    #print("rpa:", test_rpa()) 
    925979    #test_empty_distribution() 
     980    #test_structure_factor_background() 

  • 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 

  • sasmodels/weights.py

    r3d58247 re2592f0  
    2323    default = dict(npts=35, width=0, nsigmas=3) 
    2424    def __init__(self, npts=None, width=None, nsigmas=None): 
    25         self.npts = self.default['npts'] if npts is None else npts 
     25        self.npts = self.default['npts'] if npts is None else int(npts) 
    2626        self.width = self.default['width'] if width is None else width 
    2727        self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas 
Note: See TracChangeset for help on using the changeset viewer.