Changeset 2773c66 in sasmodels


Ignore:
Timestamp:
Sep 8, 2018 10:29:05 AM (7 years ago)
Author:
Torin Cooper-Bennun <torin.cooper-bennun@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
ba51e00
Parents:
b763f9d (diff), c88f983 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'beta_approx' into beta_approx_new_R_eff

Files:
3 added
70 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/scripting.rst

    r4aa5dce rbd7630d  
    1010The key functions are :func:`sasmodels.core.load_model` for loading the 
    1111model definition and compiling the kernel and 
    12 :func:`sasmodels.data.load_data` for calling sasview to load the data. Need 
    13 the data because that defines the resolution function and the q values to 
    14 evaluate. If there is no data, then use :func:`sasmodels.data.empty_data1D` 
    15 or :func:`sasmodels.data.empty_data2D` to create some data with a given $q$. 
    16  
    17 Using sasmodels through bumps 
    18 ============================= 
    19  
    20 With the data and the model, you can wrap it in a *bumps* model with 
     12:func:`sasmodels.data.load_data` for calling sasview to load the data. 
     13 
     14Preparing data 
     15============== 
     16 
     17Usually you will load data via the sasview loader, with the 
     18:func:`sasmodels.data.load_data` function.  For example:: 
     19 
     20    from sasmodels.data import load_data 
     21    data = load_data("sasmodels/example/093191_201.dat") 
     22 
     23You may want to apply a data mask, such a beam stop, and trim high $q$:: 
     24 
     25    from sasmodels.data import set_beam_stop 
     26    set_beam_stop(data, qmin, qmax) 
     27 
     28The :func:`sasmodels.data.set_beam_stop` method simply sets the *mask* 
     29attribute for the data. 
     30 
     31The data defines the resolution function and the q values to evaluate, so 
     32even if you simulating experiments prior to making measurements, you still 
     33need a data object for reference. Use :func:`sasmodels.data.empty_data1D` 
     34or :func:`sasmodels.data.empty_data2D` to create a container with a 
     35given $q$ and $\Delta q/q$.  For example:: 
     36 
     37    import numpy as np 
     38    from sasmodels.data import empty_data1D 
     39 
     40    # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5% 
     41    q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120) 
     42    data = empty_data1D(q, resolution=0.05) 
     43 
     44To use a more realistic model of resolution, or to load data from a file 
     45format not understood by SasView, you can use :class:`sasmodels.data.Data1D` 
     46or :class:`sasmodels.data.Data2D` directly.  The 1D data uses 
     47*x*, *y*, *dx* and *dy* for $x = q$ and $y = I(q)$, and 2D data uses 
     48*x*, *y*, *z*, *dx*, *dy*, *dz* for $x, y = qx, qy$ and $z = I(qx, qy)$. 
     49[Note: internally, the Data2D object uses SasView conventions, 
     50*qx_data*, *qy_data*, *data*, *dqx_data*, *dqy_data*, and *err_data*.] 
     51 
     52For USANS data, use 1D data, but set *dxl* and *dxw* attributes to 
     53indicate slit resolution:: 
     54 
     55    data.dxl = 0.117 
     56 
     57See :func:`sasmodels.resolution.slit_resolution` for details. 
     58 
     59SESANS data is more complicated; if your SESANS format is not supported by 
     60SasView you need to define a number of attributes beyond *x*, *y*.  For 
     61example:: 
     62 
     63    SElength = np.linspace(0, 2400, 61) # [A] 
     64    data = np.ones_like(SElength) 
     65    err_data = np.ones_like(SElength)*0.03 
     66 
     67    class Source: 
     68        wavelength = 6 # [A] 
     69        wavelength_unit = "A" 
     70    class Sample: 
     71        zacceptance = 0.1 # [A^-1] 
     72        thickness = 0.2 # [cm] 
     73 
     74    class SESANSData1D: 
     75        #q_zmax = 0.23 # [A^-1] 
     76        lam = 0.2 # [nm] 
     77        x = SElength 
     78        y = data 
     79        dy = err_data 
     80        sample = Sample() 
     81    data = SESANSData1D() 
     82 
     83    x, y = ... # create or load sesans 
     84    data = smd.Data 
     85 
     86The *data* module defines various data plotters as well. 
     87 
     88Using sasmodels directly 
     89======================== 
     90 
     91Once you have a computational kernel and a data object, you can evaluate 
     92the model for various parameters using 
     93:class:`sasmodels.direct_model.DirectModel`.  The resulting object *f* 
     94will be callable as *f(par=value, ...)*, returning the $I(q)$ for the $q$ 
     95values in the data.  For example:: 
     96 
     97    import numpy as np 
     98    from sasmodels.data import empty_data1D 
     99    from sasmodels.core import load_model 
     100    from sasmodels.direct_model import DirectModel 
     101 
     102    # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5% 
     103    q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120) 
     104    data = empty_data1D(q, resolution=0.05) 
     105    kernel = load_model("ellipsoid) 
     106    f = DirectModel(data, kernel) 
     107    Iq = f(radius_polar=100) 
     108 
     109Polydispersity information is set with special parameter names: 
     110 
     111    * *par_pd* for polydispersity width, $\Delta p/p$, 
     112    * *par_pd_n* for the number of points in the distribution, 
     113    * *par_pd_type* for the distribution type (as a string), and 
     114    * *par_pd_nsigmas* for the limits of the distribution. 
     115 
     116Using sasmodels through the bumps optimizer 
     117=========================================== 
     118 
     119Like DirectModel, you can wrap data and a kernel in a *bumps* model with 
    21120class:`sasmodels.bumps_model.Model` and create an 
    22 class:`sasmodels.bump_model.Experiment` that you can fit with the *bumps* 
     121class:`sasmodels.bumps_model.Experiment` that you can fit with the *bumps* 
    23122interface. Here is an example from the *example* directory such as 
    24123*example/model.py*:: 
    … …  
    75174    SasViewCom bumps.cli example/model.py --preview 
    76175 
    77 Using sasmodels directly 
    78 ======================== 
    79  
    80 Bumps has a notion of parameter boxes in which you can set and retrieve 
    81 values.  Instead of using bumps, you can create a directly callable function 
    82 with :class:`sasmodels.direct_model.DirectModel`.  The resulting object *f* 
    83 will be callable as *f(par=value, ...)*, returning the $I(q)$ for the $q$ 
    84 values in the data.  Polydisperse parameters use the same naming conventions 
    85 as in the bumps model, with e.g., radius_pd being the polydispersity associated 
    86 with radius. 
     176Calling the computation kernel 
     177============================== 
    87178 
    88179Getting a simple function that you can call on a set of q values and return 

  • sasmodels/compare.py

    raf7a97c rbd7630d  
    12881288 
    12891289    if opts['datafile'] is not None: 
    1290         data = load_data(os.path.expanduser(opts['datafile'])) 
     1290        data0 = load_data(os.path.expanduser(opts['datafile'])) 
     1291        data = data0, data0 
    12911292    else: 
    12921293        # Hack around the fact that make_data doesn't take a pair of resolutions 

  • sasmodels/data.py

    r7e923c2 rc88f983  
    183183    *x* is spin echo length and *y* is polarization (P/P0). 
    184184    """ 
     185    isSesans = True 
    185186    def __init__(self, **kw): 
    186187        Data1D.__init__(self, **kw) 
    … …  
    301302        self.wavelength_unit = "A" 
    302303 
     304class Sample(object): 
     305    """ 
     306    Sample attributes. 
     307    """ 
     308    def __init__(self): 
     309        # type: () -> None 
     310        pass 
    303311 
    304312def empty_data1D(q, resolution=0.0, L=0., dL=0.): 

  • sasmodels/direct_model.py

    r1a8c11c r7b9e4dd  
    3131from . import resolution2d 
    3232from .details import make_kernel_args, dispersion_mesh 
     33from .modelinfo import DEFAULT_BACKGROUND 
    3334 
    3435# pylint: disable=unused-import 
    … …  
    349350 
    350351        # Need to pull background out of resolution for multiple scattering 
    351         background = pars.get('background', 0.) 
     352        background = pars.get('background', DEFAULT_BACKGROUND) 
    352353        pars = pars.copy() 
    353354        pars['background'] = 0. 

  • sasmodels/generate.py

    r5399809 r2773c66  
    991991    docs = model_info.docs if model_info.docs is not None else "" 
    992992    docs = convert_section_titles_to_boldface(docs) 
    993     pars = make_partable(model_info.parameters.COMMON 
    994                          + model_info.parameters.kernel_parameters) 
     993    if model_info.structure_factor: 
     994        pars = model_info.parameters.kernel_parameters 
     995    else: 
     996        pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters 
     997    partable = make_partable(pars) 
    995998    subst = dict(id=model_info.id.replace('_', '-'), 
    996999                 name=model_info.name, 
    9971000                 title=model_info.title, 
    998                  parameters=pars, 
     1001                 parameters=partable, 
    9991002                 returns=Sq_units if model_info.structure_factor else Iq_units, 
    10001003                 docs=docs) 

  • sasmodels/kernel_iq.c

    rc57ee9e r2773c66  
    192192    QACRotation *rotation, 
    193193    double qx, double qy, 
    194     double *qa_out, double *qc_out) 
     194    double *qab_out, double *qc_out) 
    195195{ 
     196    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    196197    const double dqc = rotation->R31*qx + rotation->R32*qy; 
    197     // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    198     const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
    199  
    200     *qa_out = dqa; 
     198    const double dqab_sq = -dqc*dqc + qx*qx + qy*qy; 
     199    //*qab_out = sqrt(fabs(dqab_sq)); 
     200    *qab_out = dqab_sq > 0.0 ? sqrt(dqab_sq) : 0.0; 
    201201    *qc_out = dqc; 
    202202} 

  • sasmodels/modelinfo.py

    r5399809 r2773c66  
    4545# Note that scale and background cannot be coordinated parameters whose value 
    4646# depends on the some polydisperse parameter with the current implementation 
     47DEFAULT_BACKGROUND = 1e-3 
    4748COMMON_PARAMETERS = [ 
    4849    ("scale", "", 1, (0.0, np.inf), "", "Source intensity"), 
    49     ("background", "1/cm", 1e-3, (-np.inf, np.inf), "", "Source background"), 
     50    ("background", "1/cm", DEFAULT_BACKGROUND, (-np.inf, np.inf), "", "Source background"), 
    5051] 
    5152assert (len(COMMON_PARAMETERS) == 2 

  • sasmodels/models/ellipsoid.py

    rd277229 r2773c66  
    125125import numpy as np 
    126126from numpy import inf, sin, cos, pi 
     127 
     128try: 
     129    from numpy import cbrt 
     130except ImportError: 
     131    def cbrt(x): return x ** (1.0/3.0) 
    127132 
    128133name = "ellipsoid" 

  • sasmodels/models/spinodal.py

    r71b751d rc88f983  
    33---------- 
    44 
    5 This model calculates the SAS signal of a phase separating solution 
    6 under spinodal decomposition. The scattering intensity $I(q)$ is calculated as 
     5This model calculates the SAS signal of a phase separating system  
     6undergoing spinodal decomposition. The scattering intensity $I(q)$ is calculated  
     7as  
    78 
    89.. math:: 
    910    I(q) = I_{max}\frac{(1+\gamma/2)x^2}{\gamma/2+x^{2+\gamma}}+B 
    1011 
    11 where $x=q/q_0$ and $B$ is a flat background. The characteristic structure 
    12 length scales with the correlation peak at $q_0$. The exponent $\gamma$ is 
    13 equal to $d+1$ with d the dimensionality of the off-critical concentration 
    14 mixtures. A transition to $\gamma=2d$ is seen near the percolation threshold 
    15 into the critical concentration regime. 
     12where $x=q/q_0$, $q_0$ is the peak position, $I_{max}$ is the intensity  
     13at $q_0$ (parameterised as the $scale$ parameter), and $B$ is a flat  
     14background. The spinodal wavelength is given by $2\pi/q_0$.  
     15 
     16The exponent $\gamma$ is equal to $d+1$ for off-critical concentration  
     17mixtures (smooth interfaces) and $2d$ for critical concentration mixtures  
     18(entangled interfaces), where $d$ is the dimensionality (ie, 1, 2, 3) of the  
     19system. Thus 2 <= $\gamma$ <= 6. A transition from $\gamma=d+1$ to $\gamma=2d$  
     20is expected near the percolation threshold.  
     21 
     22As this function tends to zero as $q$ tends to zero, in practice it may be  
     23necessary to combine it with another function describing the low-angle  
     24scattering, or to simply omit the low-angle scattering from the fit. 
    1625 
    1726References 
    … …  
    2231Physica A 123,497 (1984). 
    2332 
    24 Authorship and Verification 
    25 ---------------------------- 
     33Revision History 
     34---------------- 
    2635 
    27 * **Author:** Dirk Honecker **Date:** Oct 7, 2016 
     36* **Author:**  Dirk Honecker **Date:** Oct 7, 2016 
     37* **Revised:** Steve King    **Date:** Sep 7, 2018 
    2838""" 
    2939 
    … …  
    3444title = "Spinodal decomposition model" 
    3545description = """\ 
    36       I(q) = scale ((1+gamma/2)x^2)/(gamma/2+x^(2+gamma))+background 
     46      I(q) = Imax ((1+gamma/2)x^2)/(gamma/2+x^(2+gamma)) + background 
    3747 
    3848      List of default parameters: 
    39       scale = scaling 
    40       gamma = exponent 
    41       x = q/q_0 
     49       
     50      Imax = correlation peak intensity at q_0 
     51      background = incoherent background 
     52      gamma = exponent (see model documentation) 
    4253      q_0 = correlation peak position [1/A] 
    43       background = Incoherent background""" 
     54      x = q/q_0""" 
     55       
    4456category = "shape-independent" 
    4557 

  • sasmodels/product.py

    r33d7be3 r2773c66  
    4242    pars = [] 
    4343    if p_info.have_Fq: 
    44         par = parse_parameter("structure_factor_mode", "", 0, [["P*S","P*(1+beta*(S-1))"]], "", 
    45                         "Structure factor calculation") 
     44        par = parse_parameter( 
     45                "structure_factor_mode", 
     46                "", 
     47                0, 
     48                [["P*S","P*(1+beta*(S-1))"]], 
     49                "", 
     50                "Structure factor calculation") 
    4651        pars.append(par) 
    4752    if p_info.effective_radius_type is not None: 
    48         par = parse_parameter("radius_effective_mode", "", 0, [["unconstrained"] + p_info.effective_radius_type], "", 
    49                         "Effective radius calculation") 
     53        par = parse_parameter( 
     54                "radius_effective_mode", 
     55                "", 
     56                0, 
     57                [["unconstrained"] + p_info.effective_radius_type], 
     58                "", 
     59                "Effective radius calculation") 
    5060        pars.append(par) 
    5161    return pars 
    … …  
    169179    )) 
    170180 
    171 def _intermediates_beta(F1, F2, S, scale, bg, effective_radius): 
    172     # type: (np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, float) -> OrderedDict[str, Union[np.ndarray, float]] 
     181def _intermediates_beta(F1,              # type: np.ndarray 
     182                        F2,              # type: np.ndarray 
     183                        S,               # type: np.ndarray 
     184                        scale,           # type: np.ndarray 
     185                        bg,              # type: np.ndarray 
     186                        effective_radius # type: float 
     187                        ): 
     188    # type: (...) -> OrderedDict[str, Union[np.ndarray, float]] 
    173189    """ 
    174190    Returns intermediate results for beta approximation-enabled product. The result may be an array or a float. 

  • sasmodels/sasview_model.py

    raa44a6a r2773c66  
    543543            # pylint: disable=unpacking-non-sequence 
    544544            q, phi = x 
    545             return self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)])[0] 
    546         else: 
    547             return self.calculate_Iq([x])[0] 
     545            result, _ = self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)]) 
     546            return result[0] 
     547        else: 
     548            result, _ = self.calculate_Iq([x]) 
     549            return result[0] 
    548550 
    549551 
    … …  
    560562        """ 
    561563        if isinstance(x, (list, tuple)): 
    562             return self.calculate_Iq([x[0]], [x[1]])[0] 
    563         else: 
    564             return self.calculate_Iq([x])[0] 
     564            result, _ = self.calculate_Iq([x[0]], [x[1]]) 
     565            return result[0] 
     566        else: 
     567            result, _ = self.calculate_Iq([x]) 
     568            return result[0] 
    565569 
    566570    def evalDistribution(self, qdist): 
    … …  
    596600            # Check whether we have a list of ndarrays [qx,qy] 
    597601            qx, qy = qdist 
    598             return self.calculate_Iq(qx, qy) 
     602            result, _ = self.calculate_Iq(qx, qy) 
     603            return result 
    599604 
    600605        elif isinstance(qdist, np.ndarray): 
    601606            # We have a simple 1D distribution of q-values 
    602             return self.calculate_Iq(qdist) 
     607            result, _ = self.calculate_Iq(qdist) 
     608            return result 
    603609 
    604610        else: 
    … …  
    640646        if composition and composition[0] == 'product': # only P*S for now 
    641647            with calculation_lock: 
    642                 self._calculate_Iq(qx) 
    643                 # for compatibility with sasview 4.3 
    644                 results = self._intermediate_results() 
    645                 return results["P(Q)"], results["S(Q)"] 
     648                _, lazy_results = self._calculate_Iq(qx) 
     649                # for compatibility with sasview 4.x 
     650                results = lazy_results() 
     651                pq_data = results.get("P(Q)") 
     652                sq_data = results.get("S(Q)") 
     653                return pq_data, sq_data 
    646654        else: 
    647655            return None 
    648656 
    649     def calculate_Iq(self, qx, qy=None): 
    650         # type: (Sequence[float], Optional[Sequence[float]]) -> np.ndarray 
     657    def calculate_Iq(self, 
     658                     qx,     # type: Sequence[float] 
     659                     qy=None # type: Optional[Sequence[float]] 
     660                     ): 
     661        # type: (...) -> Tuple[np.ndarray, Callable[[], collections.OrderedDict[str, np.ndarray]]] 
    651662        """ 
    652663        Calculate Iq for one set of q with the current parameters. 
    … …  
    656667        This should NOT be used for fitting since it copies the *q* vectors 
    657668        to the card for each evaluation. 
     669 
     670        The returned tuple contains the scattering intensity followed by a 
     671        callable which returns a dictionary of intermediate data from 
     672        ProductKernel. 
    658673        """ 
    659674        ## uncomment the following when trying to debug the uncoordinated calls 
    … …  
    688703        result = calculator(call_details, values, cutoff=self.cutoff, 
    689704                            magnetic=is_magnetic) 
     705        lazy_results = getattr(calculator, 'results', 
     706                               lambda: collections.OrderedDict()) 
    690707        #print("result", result) 
    691         self._intermediate_results = getattr(calculator, 'results', None) 
     708 
    692709        calculator.release() 
    693710        #self._model.release() 
    694         return result 
     711 
     712        return result, lazy_results 
    695713 
    696714    def calculate_ER(self): 
    … …  
    881899    q = np.linspace(-0.35, 0.35, 500) 
    882900    qx, qy = np.meshgrid(q, q) 
    883     result = model.calculate_Iq(qx.flatten(), qy.flatten()) 
     901    result, _ = model.calculate_Iq(qx.flatten(), qy.flatten()) 
    884902    result = result.reshape(qx.shape) 
    885903 

  • explore/beta/sasfit_compare.py

    r7b0abf8 r2a12351b  
    208208    F1, F2 = np.zeros_like(q), np.zeros_like(q) 
    209209    for k, Rpk in enumerate(Rp_val): 
     210        print("ellipsoid cycle", k, "of", len(Rp_val)) 
    210211        for i, Rei in enumerate(Re_val): 
    211212            theory = ellipsoid_theta(q, Rpk, Rei, sld, sld_solvent) 
    … …  
    343344    #Sq = Iq/Pq 
    344345    #Iq = None#= Sq = None 
    345     r = I._kernel.results 
    346     return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r[1], Ibeta=Iq) 
     346    r = dict(I._kernel.results()) 
     347    return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r["S_eff(Q)"], Ibeta=Iq) 
    347348 
    348349def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'): 
    … …  
    401402        radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, 
    402403        background=0, 
    403         radius_polar_pd=.1, radius_polar_pd_type=pd_type, 
    404         radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type, 
     404        radius_polar_pd=0.1, radius_polar_pd_type=pd_type, 
     405        radius_equatorial_pd=0.1, radius_equatorial_pd_type=pd_type, 
    405406        volfraction=0.15, 
    406407        radius_effective=270.7543927018, 
    407408        #radius_effective=12.59921049894873, 
    408409        ) 
    409     target = sasmodels_theory(q, model, beta_mode=1, **pars) 
     410    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
    410411    actual = ellipsoid_pe(q, norm='sasview', **pars) 
    411412    title = " ".join(("sasmodels", model, pd_type)) 

  • sasmodels/core.py

    r71b751d r5399809  
    363363    model = load_model("cylinder@hardsphere*sphere") 
    364364    actual = [p.name for p in model.info.parameters.kernel_parameters] 
    365     target = ("sld sld_solvent radius length theta phi volfraction" 
    366               " beta_mode A_sld A_sld_solvent A_radius").split() 
     365    target = ("sld sld_solvent radius length theta phi" 
     366              " radius_effective volfraction structure_factor_mode" 
     367              " A_sld A_sld_solvent A_radius").split() 
    367368    assert target == actual, "%s != %s"%(target, actual) 
    368369 

  • sasmodels/kernel.py

    r2d81cfe r5399809  
    4141    dtype = None  # type: np.dtype 
    4242 
    43     def __call__(self, call_details, values, cutoff, magnetic): 
     43    def Iq(self, call_details, values, cutoff, magnetic): 
    4444        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    45         raise NotImplementedError("need to implement __call__") 
     45        Pq, Reff = self.Pq_Reff(call_details, values, cutoff, magnetic, effective_radius_type=0) 
     46        return Pq 
     47    __call__ = Iq 
     48 
     49    def Pq_Reff(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     50        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     51        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     52        #print("returned",self.q_input.q, self.result) 
     53        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     54        total_weight = self.result[nout*self.q_input.nq + 0] 
     55        if total_weight == 0.: 
     56            total_weight = 1. 
     57        weighted_volume = self.result[nout*self.q_input.nq + 1] 
     58        weighted_radius = self.result[nout*self.q_input.nq + 2] 
     59        effective_radius = weighted_radius/total_weight 
     60        # compute I = scale*P + background 
     61        #           = scale*(sum(w*F^2)/sum w)/(sum (w*V)/sum w) + background 
     62        #           = scale/sum (w*V) * sum(w*F^2) + background 
     63        F2 = self.result[0:nout*self.q_input.nq:nout] 
     64        scale = values[0]/(weighted_volume if weighted_volume != 0.0 else 1.0) 
     65        background = values[1] 
     66        Pq = scale*F2 + background 
     67        #print("scale",scale,background) 
     68        return Pq, effective_radius 
     69 
     70    def beta(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     71        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     72        if self.dim == '2d': 
     73            raise NotImplementedError("beta not yet supported for 2D") 
     74        if not self.info.have_Fq: 
     75            raise NotImplementedError("beta not yet supported for "+self.info.id) 
     76        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     77        total_weight = self.result[2*self.q_input.nq + 0] 
     78        if total_weight == 0.: 
     79            total_weight = 1. 
     80        weighted_volume = self.result[2*self.q_input.nq + 1] 
     81        weighted_radius = self.result[2*self.q_input.nq + 2] 
     82        volume_average = weighted_volume/total_weight 
     83        effective_radius = weighted_radius/total_weight 
     84        F2 = self.result[0:2*self.q_input.nq:2]/total_weight 
     85        F1 = self.result[1:2*self.q_input.nq:2]/total_weight 
     86        return F1, F2, volume_average, effective_radius 
    4687 
    4788    def release(self): 
    4889        # type: () -> None 
    4990        pass 
     91 
     92    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     93        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     94        raise NotImpmentedError() 

  • sasmodels/kernel_header.c

    r108e70e r296c52b  
    6868         #define NEED_EXPM1 
    6969         #define NEED_TGAMMA 
     70         #define NEED_CBRT 
    7071         // expf missing from windows? 
    7172         #define expf exp 
    … …  
    8586#  define pown(a,b) pow(a,b) 
    8687#endif // !USE_OPENCL 
     88 
     89#if defined(NEED_CBRT) 
     90   #define cbrt(_x) pow(_x, 0.33333333333333333333333) 
     91#endif 
    8792 
    8893#if defined(NEED_EXPM1) 

  • sasmodels/kernelcl.py

    rc036ddb r5399809  
    481481        # at this point, so instead using 32, which is good on the set of 
    482482        # architectures tested so far. 
     483        extra_q = 3  # total weight, weighted volume and weighted radius 
    483484        if self.is_2d: 
    484             # Note: 16 rather than 15 because result is 1 longer than input. 
    485             width = ((self.nq+16)//16)*16 
     485            width = ((self.nq+15+extra_q)//16)*16 
    486486            self.q = np.empty((width, 2), dtype=dtype) 
    487487            self.q[:self.nq, 0] = q_vectors[0] 
    488488            self.q[:self.nq, 1] = q_vectors[1] 
    489489        else: 
    490             # Note: 32 rather than 31 because result is 1 longer than input. 
    491             width = ((self.nq+32)//32)*32 
     490            width = ((self.nq+31+extra_q)//32)*32 
    492491            self.q = np.empty(width, dtype=dtype) 
    493492            self.q[:self.nq] = q_vectors[0] 
    … …  
    539538        self.dim = '2d' if q_input.is_2d else '1d' 
    540539        # leave room for f1/f2 results in case we need to compute beta for 1d models 
    541         num_returns = 1 if self.dim == '2d' else 2  # 
    542         # plus 1 for the normalization value 
    543         self.result = np.empty((q_input.nq+1)*num_returns, dtype) 
     540        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     541        # plus 3 weight, volume, radius 
     542        self.result = np.empty(q_input.nq*nout + 3, self.dtype) 
    544543 
    545544        # Inputs and outputs for each kernel call 
    … …  
    549548 
    550549        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    551                                   q_input.global_size[0] * num_returns * dtype.itemsize) 
     550                                  q_input.global_size[0] * nout * dtype.itemsize) 
    552551        self.q_input = q_input # allocated by GpuInput above 
    553552 
    … …  
    558557                     else np.float32)  # will never get here, so use np.float32 
    559558 
    560     def Iq(self, call_details, values, cutoff, magnetic): 
    561         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    562         self._call_kernel(call_details, values, cutoff, magnetic) 
    563         #print("returned",self.q_input.q, self.result) 
    564         pd_norm = self.result[self.q_input.nq] 
    565         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    566         background = values[1] 
    567         #print("scale",scale,background) 
    568         return scale*self.result[:self.q_input.nq] + background 
    569     __call__ = Iq 
    570  
    571     def beta(self, call_details, values, cutoff, magnetic): 
    572         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    573         if self.dim == '2d': 
    574             raise NotImplementedError("beta not yet supported for 2D") 
    575         self._call_kernel(call_details, values, cutoff, magnetic) 
    576         w_norm = self.result[2*self.q_input.nq + 1] 
    577         pd_norm = self.result[self.q_input.nq] 
    578         if w_norm == 0.: 
    579             w_norm = 1. 
    580         F2 = self.result[:self.q_input.nq]/w_norm 
    581         F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm 
    582         volume_avg = pd_norm/w_norm 
    583         return F1, F2, volume_avg 
    584  
    585     def _call_kernel(self, call_details, values, cutoff, magnetic): 
     559    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    586560        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    587561        context = self.queue.context 
    … …  
    597571            details_b, values_b, self.q_input.q_b, self.result_b, 
    598572            self.real(cutoff), 
     573            np.uint32(effective_radius_type), 
    599574        ] 
    600575        #print("Calling OpenCL") 

  • sasmodels/kerneldll.py

    r2f8cbb9 r6e7ba14  
    315315 
    316316        # int, int, int, int*, double*, double*, double*, double*, double 
    317         argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type] 
     317        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type, ct.c_int32] 
    318318        names = [generate.kernel_name(self.info, variant) 
    319319                 for variant in ("Iq", "Iqxy", "Imagnetic")] 
    … …  
    382382        self.dim = '2d' if q_input.is_2d else '1d' 
    383383        # leave room for f1/f2 results in case we need to compute beta for 1d models 
    384         num_returns = 1 if self.dim == '2d' else 2  # 
    385         # plus 1 for the normalization value 
    386         self.result = np.empty((q_input.nq+1)*num_returns, self.dtype) 
     384        nout = 2 if self.info.have_Fq else 1 
     385        # plus 3 weight, volume, radius 
     386        self.result = np.empty(q_input.nq*nout + 3, self.dtype) 
    387387        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    388388                     else np.float64 if self.q_input.dtype == generate.F64 
    389389                     else np.float128) 
    390390 
    391     def Iq(self, call_details, values, cutoff, magnetic): 
    392         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    393         self._call_kernel(call_details, values, cutoff, magnetic) 
    394         #print("returned",self.q_input.q, self.result) 
    395         pd_norm = self.result[self.q_input.nq] 
    396         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    397         background = values[1] 
    398         #print("scale",scale,background) 
    399         return scale*self.result[:self.q_input.nq] + background 
    400     __call__ = Iq 
    401  
    402     def beta(self, call_details, values, cutoff, magnetic): 
    403         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    404         if self.dim == '2d': 
    405             raise NotImplementedError("beta not yet supported for 2D") 
    406         self._call_kernel(call_details, values, cutoff, magnetic) 
    407         w_norm = self.result[2*self.q_input.nq + 1] 
    408         pd_norm = self.result[self.q_input.nq] 
    409         if w_norm == 0.: 
    410             w_norm = 1. 
    411         F2 = self.result[:self.q_input.nq]/w_norm 
    412         F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm 
    413         volume_avg = pd_norm/w_norm 
    414         return F1, F2, volume_avg 
    415  
    416     def _call_kernel(self, call_details, values, cutoff, magnetic): 
     391    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     392        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    417393        kernel = self.kernel[1 if magnetic else 0] 
    418394        args = [ 
    … …  
    425401            self.result.ctypes.data,   # results 
    426402            self.real(cutoff), # cutoff 
     403            effective_radius_type, # cutoff 
    427404        ] 
    428405        #print("Calling DLL") 

  • sasmodels/kernelpy.py

    r108e70e r5399809  
    1212 
    1313import numpy as np  # type: ignore 
     14from numpy import pi 
     15try: 
     16    from numpy import cbrt 
     17except ImportError: 
     18    def cbrt(x): return x ** (1.0/3.0) 
    1419 
    1520from .generate import F64 
    … …  
    158163 
    159164        # Generate a closure which calls the form_volume if it exists. 
    160         form_volume = model_info.form_volume 
    161         self._volume = ((lambda: form_volume(*volume_args)) if form_volume else 
    162                         (lambda: 1.0)) 
    163  
    164     def __call__(self, call_details, values, cutoff, magnetic): 
     165        self._volume_args = volume_args 
     166        volume = model_info.form_volume 
     167        radius = model_info.effective_radius 
     168        self._volume = ((lambda: volume(*volume_args)) if volume 
     169                        else (lambda: 1.0)) 
     170        self._radius = ((lambda mode: radius(mode, *volume_args)) if radius 
     171                        else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume 
     172                        else (lambda mode: 1.0)) 
     173 
     174 
     175 
     176    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    165177        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    166178        if magnetic: 
    … …  
    168180        #print("Calling python kernel") 
    169181        #call_details.show(values) 
    170         res = _loops(self._parameter_vector, self._form, self._volume, 
    171                      self.q_input.nq, call_details, values, cutoff) 
    172         return res 
     182        radius = ((lambda: 0.0) if effective_radius_type == 0 
     183                  else (lambda: self._radius(effective_radius_type))) 
     184        self.result = _loops( 
     185            self._parameter_vector, self._form, self._volume, radius, 
     186            self.q_input.nq, call_details, values, cutoff) 
    173187 
    174188    def release(self): 
    … …  
    183197           form,          # type: Callable[[], np.ndarray] 
    184198           form_volume,   # type: Callable[[], float] 
     199           form_radius,   # type: Callable[[], float] 
    185200           nq,            # type: int 
    186201           call_details,  # type: details.CallDetails 
    187202           values,        # type: np.ndarray 
    188            cutoff         # type: float 
     203           cutoff,        # type: float 
    189204          ): 
    190205    # type: (...) -> None 
    … …  
    198213    #                                                              # 
    199214    ################################################################ 
     215 
     216    # WARNING: Trickery ahead 
     217    # The parameters[] vector is embedded in the closures for form(), 
     218    # form_volume() and form_radius().  We set the initial vector from 
     219    # the values for the model parameters. As we loop through the polydispesity 
     220    # mesh, we update the components with the polydispersity values before 
     221    # calling the respective functions. 
    200222    n_pars = len(parameters) 
    201223    parameters[:] = values[2:n_pars+2] 
     224 
    202225    if call_details.num_active == 0: 
    203         pd_norm = float(form_volume()) 
    204         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    205         background = values[1] 
    206         return scale*form() + background 
    207  
    208     pd_value = values[2+n_pars:2+n_pars + call_details.num_weights] 
    209     pd_weight = values[2+n_pars + call_details.num_weights:] 
    210  
    211     pd_norm = 0.0 
    212     partial_weight = np.NaN 
    213     weight = np.NaN 
    214  
    215     p0_par = call_details.pd_par[0] 
    216     p0_length = call_details.pd_length[0] 
    217     p0_index = p0_length 
    218     p0_offset = call_details.pd_offset[0] 
    219  
    220     pd_par = call_details.pd_par[:call_details.num_active] 
    221     pd_offset = call_details.pd_offset[:call_details.num_active] 
    222     pd_stride = call_details.pd_stride[:call_details.num_active] 
    223     pd_length = call_details.pd_length[:call_details.num_active] 
    224  
    225     total = np.zeros(nq, 'd') 
    226     for loop_index in range(call_details.num_eval): 
    227         # update polydispersity parameter values 
    228         if p0_index == p0_length: 
    229             pd_index = (loop_index//pd_stride)%pd_length 
    230             parameters[pd_par] = pd_value[pd_offset+pd_index] 
    231             partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
    232             p0_index = loop_index%p0_length 
    233  
    234         weight = partial_weight * pd_weight[p0_offset + p0_index] 
    235         parameters[p0_par] = pd_value[p0_offset + p0_index] 
    236         p0_index += 1 
    237         if weight > cutoff: 
    238             # Call the scattering function 
    239             # Assume that NaNs are only generated if the parameters are bad; 
    240             # exclude all q for that NaN.  Even better would be to have an 
    241             # INVALID expression like the C models, but that is too expensive. 
    242             Iq = np.asarray(form(), 'd') 
    243             if np.isnan(Iq).any(): 
    244                 continue 
    245  
    246             # update value and norm 
    247             total += weight * Iq 
    248             pd_norm += weight * form_volume() 
    249  
    250     scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    251     background = values[1] 
    252     return scale*total + background 
     226        total = form() 
     227        weight_norm = 1.0 
     228        weighted_volume = form_volume() 
     229        weighted_radius = form_radius() 
     230 
     231    else: 
     232        pd_value = values[2+n_pars:2+n_pars + call_details.num_weights] 
     233        pd_weight = values[2+n_pars + call_details.num_weights:] 
     234 
     235        weight_norm = 0.0 
     236        weighted_volume = 0.0 
     237        weighted_radius = 0.0 
     238        partial_weight = np.NaN 
     239        weight = np.NaN 
     240 
     241        p0_par = call_details.pd_par[0] 
     242        p0_length = call_details.pd_length[0] 
     243        p0_index = p0_length 
     244        p0_offset = call_details.pd_offset[0] 
     245 
     246        pd_par = call_details.pd_par[:call_details.num_active] 
     247        pd_offset = call_details.pd_offset[:call_details.num_active] 
     248        pd_stride = call_details.pd_stride[:call_details.num_active] 
     249        pd_length = call_details.pd_length[:call_details.num_active] 
     250 
     251        total = np.zeros(nq, 'd') 
     252        for loop_index in range(call_details.num_eval): 
     253            # update polydispersity parameter values 
     254            if p0_index == p0_length: 
     255                pd_index = (loop_index//pd_stride)%pd_length 
     256                parameters[pd_par] = pd_value[pd_offset+pd_index] 
     257                partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
     258                p0_index = loop_index%p0_length 
     259 
     260            weight = partial_weight * pd_weight[p0_offset + p0_index] 
     261            parameters[p0_par] = pd_value[p0_offset + p0_index] 
     262            p0_index += 1 
     263            if weight > cutoff: 
     264                # Call the scattering function 
     265                # Assume that NaNs are only generated if the parameters are bad; 
     266                # exclude all q for that NaN.  Even better would be to have an 
     267                # INVALID expression like the C models, but that is expensive. 
     268                Iq = np.asarray(form(), 'd') 
     269                if np.isnan(Iq).any(): 
     270                    continue 
     271 
     272                # update value and norm 
     273                total += weight * Iq 
     274                weight_norm += weight 
     275                weighted_volume += weight * form_volume() 
     276                weighted_radius += weight * form_radius() 
     277 
     278    result = np.hstack((total, weight_norm, weighted_volume, weighted_radius)) 
     279    return result 
    253280 
    254281 

  • sasmodels/models/barbell.c

    r71b751d rd277229  
    6262} 
    6363 
     64static double 
     65radius_from_volume(double radius_bell, double radius, double length) 
     66{ 
     67    const double vol_barbell = form_volume(radius_bell,radius,length); 
     68    return cbrt(0.75*vol_barbell/M_PI); 
     69} 
     70 
     71static double 
     72radius_from_totallength(double radius_bell, double radius, double length) 
     73{ 
     74    const double hdist = sqrt(square(radius_bell) - square(radius)); 
     75    return 0.5*length + hdist + radius_bell; 
     76} 
     77 
     78static double 
     79effective_radius(int mode, double radius_bell, double radius, double length) 
     80{ 
     81    if (mode == 1) { 
     82        return radius_from_volume(radius_bell, radius , length); 
     83    } else if (mode == 2) { 
     84        return radius; 
     85    } else if (mode == 3) { 
     86        return 0.5*length; 
     87    } else { 
     88        return radius_from_totallength(radius_bell,radius,length); 
     89    } 
     90} 
     91 
    6492static void 
    6593Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 

  • sasmodels/models/barbell.py

    r71b751d rd277229  
    117117source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] 
    118118have_Fq = True 
     119effective_radius_type = ["equivalent sphere","radius","half length","half total length"] 
    119120 
    120121def random(): 

  • sasmodels/models/capped_cylinder.c

    r71b751d rd277229  
    8484} 
    8585 
     86static double 
     87radius_from_volume(double radius, double radius_cap, double length) 
     88{ 
     89    const double vol_cappedcyl = form_volume(radius,radius_cap,length); 
     90    return cbrt(0.75*vol_cappedcyl/M_PI); 
     91} 
     92 
     93static double 
     94radius_from_totallength(double radius, double radius_cap, double length) 
     95{ 
     96    const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); 
     97    return 0.5*length + hc; 
     98} 
     99 
     100static double 
     101effective_radius(int mode, double radius, double radius_cap, double length) 
     102{ 
     103    if (mode == 1) { 
     104        return radius_from_volume(radius, radius_cap, length); 
     105    } else if (mode == 2) { 
     106        return radius; 
     107    } else if (mode == 3) { 
     108        return 0.5*length; 
     109    } else { 
     110        return radius_from_totallength(radius, radius_cap,length); 
     111    } 
     112} 
     113 
    86114static void 
    87115Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 

  • sasmodels/models/capped_cylinder.py

    r71b751d rd277229  
    137137source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"] 
    138138have_Fq = True 
     139effective_radius_type = ["equivalent sphere","radius","half length","half total length"] 
    139140 
    140141def random(): 

  • sasmodels/models/core_multi_shell.c

    r71b751d rd277229  
    1919} 
    2020 
     21static double 
     22outer_radius(double core_radius, double fp_n, double thickness[]) 
     23{ 
     24  double r = core_radius; 
     25  int n = (int)(fp_n+0.5); 
     26  for (int i=0; i < n; i++) { 
     27    r += thickness[i]; 
     28  } 
     29  return r; 
     30} 
     31 
     32static double 
     33effective_radius(int mode, double core_radius, double fp_n, double thickness[]) 
     34{ 
     35    if (mode == 1) { 
     36        double r = core_radius; 
     37        int n = (int)(fp_n+0.5); 
     38        for (int i=0; i < n; i++) { 
     39            r += thickness[i]; 
     40        } 
     41        return r; 
     42        //return outer_radius(core_radius,fp_n,thickness); 
     43    } else { 
     44        return core_radius; 
     45    } 
     46} 
     47 
    2148static void 
    2249Fq(double q, double *F1, double *F2, double core_sld, double core_radius, 

  • sasmodels/models/core_multi_shell.py

    r71b751d rd277229  
    100100source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 
    101101have_Fq = True 
     102effective_radius_type = ["outer radius", "core radius"] 
    102103 
    103104def random(): 
    … …  
    144145    return np.asarray(z), np.asarray(rho) 
    145146 
    146 def ER(radius, n, thickness): 
    147     """Effective radius""" 
    148     n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse 
    149     return np.sum(thickness[:n], axis=0) + radius 
     147#def ER(radius, n, thickness): 
     148#    """Effective radius""" 
     149#    n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse 
     150#    return np.sum(thickness[:n], axis=0) + radius 
    150151 
    151152demo = dict(sld_core=6.4, 

  • sasmodels/models/core_shell_bicelle.c

    r71b751d rd277229  
    3434 
    3535    return t; 
     36} 
     37 
     38static double 
     39radius_from_volume(double radius, double thick_rim, double thick_face, double length) 
     40{ 
     41    const double volume_bicelle = form_volume(radius,thick_rim,thick_face,length); 
     42    return cbrt(0.75*volume_bicelle/M_PI); 
     43} 
     44 
     45static double 
     46radius_from_diagonal(double radius, double thick_rim, double thick_face, double length) 
     47{ 
     48    const double radius_tot = radius + thick_rim; 
     49    const double length_tot = length + 2.0*thick_face; 
     50    return sqrt(radius_tot*radius_tot + 0.25*length_tot*length_tot); 
     51} 
     52 
     53static double 
     54effective_radius(int mode, double radius, double thick_rim, double thick_face, double length) 
     55{ 
     56    if (mode == 1) { 
     57        return radius_from_volume(radius, thick_rim, thick_face, length); 
     58    } else if (mode == 2) { 
     59        return radius + thick_rim; 
     60    } else if (mode == 3) { 
     61        return 0.5*length + thick_face; 
     62    } else { 
     63        return radius_from_diagonal(radius,thick_rim,thick_face,length); 
     64    } 
    3665} 
    3766 

  • sasmodels/models/core_shell_bicelle.py

    r71b751d rd277229  
    155155          "core_shell_bicelle.c"] 
    156156have_Fq = True 
     157effective_radius_type = ["equivalent sphere","outer rim radius", "half outer thickness","half diagonal"] 
    157158 
    158159def random(): 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r71b751d rd277229  
    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_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     14{ 
     15    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length); 
     16    return cbrt(0.75*volume_bicelle/M_PI); 
     17} 
     18 
     19static double 
     20radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     21{ 
     22    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor); 
     23    const double radius_max_tot = radius_max + thick_rim; 
     24    const double length_tot = length + 2.0*thick_face; 
     25    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot); 
     26} 
     27 
     28static double 
     29effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     30{ 
     31    if (mode == 1) { 
     32        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
     33    } else if (mode == 2) { 
     34        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
     35    } else if (mode == 3) { 
     36        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     37    } else if (mode == 4) { 
     38        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     39    } else if (mode ==5) { 
     40        return 0.5*length + thick_face; 
     41    } else { 
     42        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
     43    } 
    1044} 
    1145 

  • sasmodels/models/core_shell_bicelle_elliptical.py

    r71b751d rd277229  
    147147          "core_shell_bicelle_elliptical.c"] 
    148148have_Fq = True 
     149effective_radius_type = ["equivalent sphere","outer rim average radius","outer rim min radius", 
     150                         "outer max radius", "half outer thickness","half diagonal"] 
    149151 
    150152def random(): 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r71b751d rd277229  
    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_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     15{ 
     16    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length); 
     17    return cbrt(0.75*volume_bicelle/M_PI); 
     18} 
     19 
     20static double 
     21radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     22{ 
     23    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor); 
     24    const double radius_max_tot = radius_max + thick_rim; 
     25    const double length_tot = length + 2.0*thick_face; 
     26    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot); 
     27} 
     28 
     29static double 
     30effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     31{ 
     32    if (mode == 1) { 
     33        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
     34    } else if (mode == 2) { 
     35        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
     36    } else if (mode == 3) { 
     37        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     38    } else if (mode == 4) { 
     39        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     40    } else if (mode ==5) { 
     41        return 0.5*length + thick_face; 
     42    } else { 
     43        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
     44    } 
    1145} 
    1246 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    r71b751d rd277229  
    160160          "core_shell_bicelle_elliptical_belt_rough.c"] 
    161161have_Fq = True 
     162effective_radius_type = ["equivalent sphere","outer rim average radius","outer rim min radius", 
     163                         "outer max radius", "half outer thickness","half diagonal"] 
    162164 
    163165demo = dict(scale=1, background=0, 

  • sasmodels/models/core_shell_cylinder.c

    r71b751d rd277229  
    1111{ 
    1212    return M_PI*square(radius+thickness)*(length+2.0*thickness); 
     13} 
     14 
     15static double 
     16radius_from_volume(double radius, double thickness, double length) 
     17{ 
     18    const double volume_outer_cyl = form_volume(radius,thickness,length); 
     19    return cbrt(0.75*volume_outer_cyl/M_PI); 
     20} 
     21 
     22static double 
     23radius_from_diagonal(double radius, double thickness, double length) 
     24{ 
     25    const double radius_outer = radius + thickness; 
     26    const double length_outer = length + thickness; 
     27    return sqrt(radius_outer*radius_outer + 0.25*length_outer*length_outer); 
     28} 
     29 
     30static double 
     31effective_radius(int mode, double radius, double thickness, double length) 
     32{ 
     33    if (mode == 1) { 
     34        return radius_from_volume(radius, thickness, length); 
     35    } else if (mode == 2) { 
     36        return radius + thickness; 
     37    } else if (mode == 3) { 
     38        return 0.5*length + thickness; 
     39    } else if (mode == 4) { 
     40        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness); 
     41    } else if (mode == 5) { 
     42        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness); 
     43    } else { 
     44        return radius_from_diagonal(radius,thickness,length); 
     45    } 
    1346} 
    1447 

  • sasmodels/models/core_shell_cylinder.py

    r71b751d rd277229  
    126126source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] 
    127127have_Fq = True 
     128effective_radius_type = ["equivalent sphere","outer radius","half outer length","half min outer dimension", 
     129                         "half max outer dimension","half outer diagonal"] 
    128130 
    129 def ER(radius, thickness, length): 
    130     """ 
    131     Returns the effective radius used in the S*P calculation 
    132     """ 
    133     radius = radius + thickness 
    134     length = length + 2 * thickness 
    135     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    136     return 0.5 * (ddd) ** (1. / 3.) 
     131#def ER(radius, thickness, length): 
     132#    """ 
     133#    Returns the effective radius used in the S*P calculation 
     134#    """ 
     135#    radius = radius + thickness 
     136#    length = length + 2 * thickness 
     137#    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
     138#    return 0.5 * (ddd) ** (1. / 3.) 
    137139 
    138140def VR(radius, thickness, length): 

  • sasmodels/models/core_shell_ellipsoid.c

    r71b751d r3c60146  
    3636    double vol = M_4PI_3*equat_shell*equat_shell*polar_shell; 
    3737    return vol; 
     38} 
     39 
     40static double 
     41radius_from_volume(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     42{ 
     43    const double volume_ellipsoid = form_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     44    return cbrt(0.75*volume_ellipsoid/M_PI); 
     45} 
     46 
     47static double 
     48radius_from_curvature(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     49{ 
     50    // Trivial cases 
     51    if (1.0 == x_core && 1.0 == x_polar_shell) return radius_equat_core + thick_shell; 
     52    if ((radius_equat_core + thick_shell)*(radius_equat_core*x_core + thick_shell*x_polar_shell) == 0.)  return 0.; 
     53 
     54    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     55    const double radius_equat_tot = radius_equat_core + thick_shell; 
     56    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 
     57    const double ratio = (radius_polar_tot < radius_equat_tot 
     58                          ? radius_polar_tot / radius_equat_tot 
     59                          : radius_equat_tot / radius_polar_tot); 
     60    const double e1 = sqrt(1.0 - ratio*ratio); 
     61    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     62    const double bL = (1.0 + e1) / (1.0 - e1); 
     63    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     64    const double delta = 0.75 * b1 * b2; 
     65    const double ddd = 2.0 * (delta + 1.0) * radius_polar_tot * radius_equat_tot * radius_equat_tot; 
     66    return 0.5 * cbrt(ddd); 
     67} 
     68 
     69static double 
     70effective_radius(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     71{ 
     72    const double radius_equat_tot = radius_equat_core + thick_shell; 
     73    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 
     74 
     75    if (mode == 1) { 
     76        return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     77    } else if (mode == 2) { 
     78        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     79    } else if (mode == 3) { 
     80        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
     81    } else { 
     82        return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
     83    } 
    3884} 
    3985 

  • sasmodels/models/core_shell_ellipsoid.py

    r71b751d rd277229  
    146146source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    147147have_Fq = True 
    148  
    149 def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 
    150     """ 
    151         Returns the effective radius used in the S*P calculation 
    152     """ 
    153     from .ellipsoid import ER as ellipsoid_ER 
    154     polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell 
    155     equat_outer = radius_equat_core + thick_shell 
    156     return ellipsoid_ER(polar_outer, equat_outer) 
     148effective_radius_type = ["equivalent sphere","average outer curvature", "min outer radius", "max outer radius"] 
     149 
     150#def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 
     151#    """ 
     152#        Returns the effective radius used in the S*P calculation 
     153#    """ 
     154#    from .ellipsoid import ER as ellipsoid_ER 
     155#    polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell 
     156#    equat_outer = radius_equat_core + thick_shell 
     157#    return ellipsoid_ER(polar_outer, equat_outer) 
    157158 
    158159def random(): 

  • sasmodels/models/core_shell_parallelepiped.c

    r71b751d rd277229  
    2525        2.0 * length_a * length_b * thick_rim_c; 
    2626#endif 
     27} 
     28 
     29static double 
     30radius_from_volume(double length_a, double length_b, double length_c, 
     31                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     32{ 
     33    const double volume_paral = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     34    return cbrt(0.75*volume_paral/M_PI); 
     35} 
     36 
     37static double 
     38radius_from_crosssection(double length_a, double length_b, double thick_rim_a, double thick_rim_b) 
     39{ 
     40    const double area_xsec_paral = length_a*length_b + 2.0*thick_rim_a*length_b + 2.0*thick_rim_b*length_a; 
     41    return sqrt(area_xsec_paral/M_PI); 
     42} 
     43 
     44static double 
     45effective_radius(int mode, double length_a, double length_b, double length_c, 
     46                 double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     47{ 
     48    if (mode == 1) { 
     49        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     50    } else if (mode == 2) { 
     51        return 0.5 * (length_a + thick_rim_a); 
     52    } else if (mode == 3) { 
     53        return 0.5 * (length_b + thick_rim_b); 
     54    } else if (mode == 4) { 
     55        return 0.5 * (length_c + thick_rim_c); 
     56    } else if (mode == 5) { 
     57        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b); 
     58    } else if (mode == 6) { 
     59        return 0.5*sqrt(square(length_a+thick_rim_a) + square(length_b+thick_rim_b)); 
     60    } else { 
     61        return 0.5*sqrt(square(length_a+thick_rim_a) + square(length_b+thick_rim_b) + square(length_c+thick_rim_c)); 
     62    } 
    2763} 
    2864 

  • sasmodels/models/core_shell_parallelepiped.py

    r71b751d rd277229  
    227227source = ["lib/gauss76.c", "core_shell_parallelepiped.c"] 
    228228have_Fq = True 
    229  
    230  
    231 def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c): 
    232     """ 
    233         Return equivalent radius (ER) 
    234     """ 
    235     from .parallelepiped import ER as ER_p 
    236  
    237     a = length_a + 2*thick_rim_a 
    238     b = length_b + 2*thick_rim_b 
    239     c = length_c + 2*thick_rim_c 
    240     return ER_p(a, b, c) 
     229effective_radius_type = ["equivalent sphere","half outer length_a", "half outer length_b", "half outer length_c", 
     230                         "equivalent circular cross-section","half outer ab diagonal","half outer diagonal"] 
     231 
     232#def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c): 
     233#    """ 
     234#        Return equivalent radius (ER) 
     235#    """ 
     236#    from .parallelepiped import ER as ER_p 
     237# 
     238#    a = length_a + 2*thick_rim_a 
     239#    b = length_b + 2*thick_rim_b 
     240#    c = length_c + 2*thick_rim_c 
     241#    return ER_p(a, b, c) 
    241242 
    242243# VR defaults to 1.0 

  • sasmodels/models/core_shell_sphere.c

    r71b751d rd277229  
    33{ 
    44    return M_4PI_3 * cube(radius + thickness); 
     5} 
     6 
     7static double 
     8effective_radius(int mode, double radius, double thickness) 
     9{ 
     10    if (mode == 1) { 
     11        return radius + thickness; 
     12    } else { 
     13        return radius; 
     14    } 
    515} 
    616 

  • sasmodels/models/core_shell_sphere.py

    r71b751d rd277229  
    7878source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] 
    7979have_Fq = True 
     80effective_radius_type = ["outer radius", "core radius"] 
    8081 
    8182demo = dict(scale=1, background=0, radius=60, thickness=10, 
    8283            sld_core=1.0, sld_shell=2.0, sld_solvent=0.0) 
    8384 
    84 def ER(radius, thickness): 
    85     """ 
    86         Equivalent radius 
    87         @param radius: core radius 
    88         @param thickness: shell thickness 
    89     """ 
    90     return radius + thickness 
     85#def ER(radius, thickness): 
     86#    """ 
     87#        Equivalent radius 
     88#        @param radius: core radius 
     89#        @param thickness: shell thickness 
     90#    """ 
     91#    return radius + thickness 
    9192 
    9293def VR(radius, thickness): 

  • sasmodels/models/cylinder.c

    r71b751d rd277229  
    1111{ 
    1212    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
     13} 
     14 
     15static double 
     16radius_from_volume(double radius, double length) 
     17{ 
     18    return cbrt(0.75*radius*radius*length); 
     19} 
     20 
     21static double 
     22radius_from_diagonal(double radius, double length) 
     23{ 
     24    return sqrt(radius*radius + 0.25*length*length); 
     25} 
     26 
     27static double 
     28effective_radius(int mode, double radius, double length) 
     29{ 
     30    if (mode == 1) { 
     31        return radius_from_volume(radius, length); 
     32    } else if (mode == 2) { 
     33        return radius; 
     34    } else if (mode == 3) { 
     35        return 0.5*length; 
     36    } else if (mode == 4) { 
     37        return (radius < 0.5*length ? radius : 0.5*length); 
     38    } else if (mode == 5) { 
     39        return (radius > 0.5*length ? radius : 0.5*length); 
     40    } else { 
     41        return radius_from_diagonal(radius,length); 
     42    } 
    1343} 
    1444 

  • sasmodels/models/cylinder.py

    r71b751d rd277229  
    139139source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 
    140140have_Fq = True 
     141effective_radius_type = ["equivalent sphere","radius","half length","half min dimension","half max dimension","half diagonal"] 
    141142 
    142 def ER(radius, length): 
    143     """ 
    144         Return equivalent radius (ER) 
    145     """ 
    146     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    147     return 0.5 * (ddd) ** (1. / 3.) 
     143#def ER(radius, length): 
     144#    """ 
     145#        Return equivalent radius (ER) 
     146#    """ 
     147#    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
     148#    return 0.5 * (ddd) ** (1. / 3.) 
    148149 
    149150def random(): 

  • sasmodels/models/ellipsoid.c

    r71b751d rd277229  
    44    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    55} 
     6 
     7static double 
     8radius_from_volume(double radius_polar, double radius_equatorial) 
     9{ 
     10    return cbrt(radius_polar*radius_equatorial*radius_equatorial); 
     11} 
     12 
     13static double 
     14radius_from_curvature(double radius_polar, double radius_equatorial) 
     15{ 
     16    // Trivial cases 
     17    if (radius_polar == radius_equatorial) return radius_polar; 
     18    if (radius_polar * radius_equatorial == 0.)  return 0.; 
     19 
     20    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     21    const double ratio = (radius_polar < radius_equatorial 
     22                          ? radius_polar / radius_equatorial 
     23                          : radius_equatorial / radius_polar); 
     24    const double e1 = sqrt(1.0 - ratio*ratio); 
     25    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     26    const double bL = (1.0 + e1) / (1.0 - e1); 
     27    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     28    const double delta = 0.75 * b1 * b2; 
     29    const double ddd = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial * radius_equatorial; 
     30    return 0.5 * cbrt(ddd); 
     31} 
     32 
     33static double 
     34effective_radius(int mode, double radius_polar, double radius_equatorial) 
     35{ 
     36    if (mode == 1) { 
     37        return radius_from_volume(radius_polar, radius_equatorial); 
     38    } else if (mode == 2) { 
     39        return radius_from_curvature(radius_polar, radius_equatorial); 
     40    } else if (mode == 3) { 
     41        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial); 
     42    } else { 
     43        return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial); 
     44    } 
     45} 
     46 
    647 
    748static void 

  • sasmodels/models/elliptical_cylinder.c

    r71b751d rd277229  
    33{ 
    44    return M_PI * radius_minor * radius_minor * r_ratio * length; 
     5} 
     6 
     7static double 
     8radius_from_volume(double radius_minor, double r_ratio, double length) 
     9{ 
     10    const double volume_ellcyl = form_volume(radius_minor,r_ratio,length); 
     11    return cbrt(0.75*volume_ellcyl/M_PI); 
     12} 
     13 
     14static double 
     15radius_from_min_dimension(double radius_minor, double r_ratio, double length) 
     16{ 
     17    const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
     18    return (rad_min < length ? rad_min : length); 
     19} 
     20 
     21static double 
     22radius_from_max_dimension(double radius_minor, double r_ratio, double length) 
     23{ 
     24    const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
     25    return (rad_max > length ? rad_max : length); 
     26} 
     27 
     28static double 
     29radius_from_diagonal(double radius_minor, double r_ratio, double length) 
     30{ 
     31    const double radius_max = (r_ratio > 1.0 ? radius_minor*r_ratio : radius_minor); 
     32    return sqrt(radius_max*radius_max + 0.25*length*length); 
     33} 
     34 
     35static double 
     36effective_radius(int mode, double radius_minor, double r_ratio, double length) 
     37{ 
     38    if (mode == 1) { 
     39        return radius_from_volume(radius_minor, r_ratio, length); 
     40    } else if (mode == 2) { 
     41        return 0.5*radius_minor*(1.0 + r_ratio); 
     42    } else if (mode == 3) { 
     43        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
     44    } else if (mode == 4) { 
     45        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
     46    } else if (mode == 5) { 
     47        return sqrt(radius_minor*radius_minor*r_ratio); 
     48    } else if (mode == 6) { 
     49        return 0.5*length; 
     50    } else if (mode == 7) { 
     51        return radius_from_min_dimension(radius_minor,r_ratio,length); 
     52    } else if (mode == 8) { 
     53        return radius_from_max_dimension(radius_minor,r_ratio,length); 
     54    } else { 
     55        return radius_from_diagonal(radius_minor,r_ratio,length); 
     56    } 
    557} 
    658 

  • sasmodels/models/elliptical_cylinder.py

    r71b751d rd277229  
    123123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
    124124have_Fq = True 
     125effective_radius_type = ["equivalent sphere","average radius","min radius","max radius", 
     126                         "equivalent circular cross-section","half length","half min dimension","half max dimension","half diagonal"] 
    125127 
    126128demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 
    … …  
    128130            theta_pd=10, phi_pd=2, psi_pd=3) 
    129131 
    130 def ER(radius_minor, axis_ratio, length): 
    131     """ 
    132         Equivalent radius 
    133         @param radius_minor: Ellipse minor radius 
    134         @param axis_ratio: Ratio of major radius over minor radius 
    135         @param length: Length of the cylinder 
    136     """ 
    137     radius = sqrt(radius_minor * radius_minor * axis_ratio) 
    138     ddd = 0.75 * radius * (2 * radius * length 
    139                            + (length + radius) * (length + pi * radius)) 
    140     return 0.5 * (ddd) ** (1. / 3.) 
     132#def ER(radius_minor, axis_ratio, length): 
     133#    """ 
     134#        Equivalent radius 
     135#        @param radius_minor: Ellipse minor radius 
     136#        @param axis_ratio: Ratio of major radius over minor radius 
     137#        @param length: Length of the cylinder 
     138#    """ 
     139#    radius = sqrt(radius_minor * radius_minor * axis_ratio) 
     140#    ddd = 0.75 * radius * (2 * radius * length 
     141#                           + (length + radius) * (length + pi * radius)) 
     142#    return 0.5 * (ddd) ** (1. / 3.) 
    141143 
    142144def random(): 

  • sasmodels/models/fuzzy_sphere.py

    r71b751d rd277229  
    7878              ["sld_solvent", "1e-6/Ang^2",  3, [-inf, inf], "sld",    "Solvent scattering length density"], 
    7979              ["radius",      "Ang",        60, [0, inf],    "volume", "Sphere radius"], 
    80               ["fuzziness",   "Ang",        10, [0, inf],    "",       "std deviation of Gaussian convolution for interface (must be << radius)"], 
     80              ["fuzziness",   "Ang",        10, [0, inf],    "volume",       "std deviation of Gaussian convolution for interface (must be << radius)"], 
    8181             ] 
    8282# pylint: enable=bad-whitespace,line-too-long 
    8383 
    84 source = ["lib/sas_3j1x_x.c"] 
     84source = ["lib/sas_3j1x_x.c","fuzzy_sphere.c"] 
    8585have_Fq = True 
     86effective_radius_type = ["radius","radius + fuzziness"] 
    8687 
    87 c_code = """ 
    88 static double form_volume(double radius) 
    89 { 
    90     return M_4PI_3*cube(radius); 
    91 } 
    92  
    93 static void Fq(double q, double *F1, double *F2, double sld, double sld_solvent, 
    94                double radius, double fuzziness) 
    95 { 
    96     const double qr = q*radius; 
    97     const double bes = sas_3j1x_x(qr); 
    98     const double qf = exp(-0.5*square(q*fuzziness)); 
    99     const double contrast = (sld - sld_solvent); 
    100     const double form = contrast * form_volume(radius) * bes * qf; 
    101     *F1 = 1.0e-2*form; 
    102     *F2 = 1.0e-4*form*form; 
    103 } 
    104 """ 
    105  
    106 def ER(radius): 
    107     """ 
    108     Return radius 
    109     """ 
    110     return radius 
     88#def ER(radius): 
     89#    """ 
     90#    Return radius 
     91#    """ 
     92#    return radius 
    11193 
    11294# VR defaults to 1.0 

  • sasmodels/models/hollow_cylinder.c

    r71b751d rd277229  
    2222} 
    2323 
     24static double 
     25radius_from_volume(double radius, double thickness, double length) 
     26{ 
     27    const double volume_outer_cyl = M_PI*square(radius + thickness)*length; 
     28    return cbrt(0.75*volume_outer_cyl/M_PI); 
     29} 
     30 
     31static double 
     32radius_from_diagonal(double radius, double thickness, double length) 
     33{ 
     34    return sqrt(square(radius + thickness) + 0.25*square(length)); 
     35} 
     36 
     37static double 
     38effective_radius(int mode, double radius, double thickness, double length) 
     39{ 
     40    if (mode == 1) { 
     41        return radius_from_volume(radius, thickness, length); 
     42    } else if (mode == 2) { 
     43        return radius + thickness; 
     44    } else if (mode == 3) { 
     45        return 0.5*length; 
     46    } else if (mode == 4) { 
     47        return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length); 
     48    } else if (mode == 5) { 
     49        return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length); 
     50    } else { 
     51        return radius_from_diagonal(radius,thickness,length); 
     52    } 
     53} 
    2454 
    2555static void 

  • sasmodels/models/hollow_cylinder.py

    r71b751d rd277229  
    9090source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
    9191have_Fq = True 
     92effective_radius_type = ["equivalent sphere","outer radius","half length", 
     93                         "half outer min dimension","half outer max dimension","half outer diagonal"] 
    9294 
    9395# pylint: disable=W0613 
    94 def ER(radius, thickness, length): 
    95     """ 
    96     :param radius:      Cylinder core radius 
    97     :param thickness:   Cylinder wall thickness 
    98     :param length:      Cylinder length 
    99     :return:            Effective radius 
    100     """ 
    101     router = radius + thickness 
    102     if router == 0 or length == 0: 
    103         return 0.0 
    104     len1 = router 
    105     len2 = length/2.0 
    106     term1 = len1*len1*2.0*len2/2.0 
    107     term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0) 
    108     ddd = 3.0*term1*term2 
    109     diam = pow(ddd, (1.0/3.0)) 
    110     return diam 
     96#def ER(radius, thickness, length): 
     97#    """ 
     98#    :param radius:      Cylinder core radius 
     99#    :param thickness:   Cylinder wall thickness 
     100#    :param length:      Cylinder length 
     101#    :return:            Effective radius 
     102#    """ 
     103#    router = radius + thickness 
     104#    if router == 0 or length == 0: 
     105#        return 0.0 
     106#    len1 = router 
     107#    len2 = length/2.0 
     108#    term1 = len1*len1*2.0*len2/2.0 
     109#    term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0) 
     110#    ddd = 3.0*term1*term2 
     111#    diam = pow(ddd, (1.0/3.0)) 
     112#    return diam 
    111113 
    112114def VR(radius, thickness, length): 

  • sasmodels/models/hollow_rectangular_prism.c

    r71b751d rd277229  
    1111    double vol_shell = vol_total - vol_core; 
    1212    return vol_shell; 
     13} 
     14 
     15static double 
     16effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     17{ 
     18    if (mode == 1) { 
     19        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
     20    } else if (mode == 2) { 
     21        return 0.5 * length_a; 
     22    } else if (mode == 3) { 
     23        return 0.5 * length_a*b2a_ratio; 
     24    } else if (mode == 4) { 
     25        return 0.5 * length_a*c2a_ratio; 
     26    } else if (mode == 5) { 
     27        return length_a*sqrt(b2a_ratio/M_PI); 
     28    } else if (mode == 6) { 
     29        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     30    } else { 
     31        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     32    } 
    1333} 
    1434 

  • sasmodels/models/hollow_rectangular_prism.py

    r71b751d rd277229  
    143143source = ["lib/gauss76.c", "hollow_rectangular_prism.c"] 
    144144have_Fq = True 
    145  
    146 def ER(length_a, b2a_ratio, c2a_ratio, thickness): 
    147     """ 
    148     Return equivalent radius (ER) 
    149     thickness parameter not used 
    150     """ 
    151     b_side = length_a * b2a_ratio 
    152     c_side = length_a * c2a_ratio 
    153  
    154     # surface average radius (rough approximation) 
    155     surf_rad = sqrt(length_a * b_side / pi) 
    156  
    157     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    158     return 0.5 * (ddd) ** (1. / 3.) 
     145effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
     146                         "equivalent outer circular cross-section","half ab diagonal","half diagonal"] 
     147 
     148#def ER(length_a, b2a_ratio, c2a_ratio, thickness): 
     149#    """ 
     150#    Return equivalent radius (ER) 
     151#    thickness parameter not used 
     152#    """ 
     153#    b_side = length_a * b2a_ratio 
     154#    c_side = length_a * c2a_ratio 
     155# 
     156#    # surface average radius (rough approximation) 
     157#    surf_rad = sqrt(length_a * b_side / pi) 
     158# 
     159#    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
     160#    return 0.5 * (ddd) ** (1. / 3.) 
    159161 
    160162def VR(length_a, b2a_ratio, c2a_ratio, thickness): 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    r71b751d rd277229  
    66    double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
    77    return vol_shell; 
     8} 
     9 
     10static double 
     11effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio) 
     12{ 
     13    if (mode == 1) { 
     14        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
     15    } else if (mode == 2) { 
     16        return 0.5 * length_a; 
     17    } else if (mode == 3) { 
     18        return 0.5 * length_a*b2a_ratio; 
     19    } else if (mode == 4) { 
     20        return 0.5 * length_a*c2a_ratio; 
     21    } else if (mode == 5) { 
     22        return length_a*sqrt(b2a_ratio/M_PI); 
     23    } else if (mode == 6) { 
     24        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     25    } else { 
     26        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     27    } 
    828} 
    929 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    r71b751d rd277229  
    103103source = ["lib/gauss76.c", "hollow_rectangular_prism_thin_walls.c"] 
    104104have_Fq = True 
     105effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
     106                         "equivalent outer circular cross-section","half ab diagonal","half diagonal"] 
    105107 
    106 def ER(length_a, b2a_ratio, c2a_ratio): 
    107     """ 
    108         Return equivalent radius (ER) 
    109     """ 
    110     b_side = length_a * b2a_ratio 
    111     c_side = length_a * c2a_ratio 
    112  
    113     # surface average radius (rough approximation) 
    114     surf_rad = sqrt(length_a * b_side / pi) 
    115  
    116     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    117     return 0.5 * (ddd) ** (1. / 3.) 
     108#def ER(length_a, b2a_ratio, c2a_ratio): 
     109#    """ 
     110#        Return equivalent radius (ER) 
     111#    """ 
     112#    b_side = length_a * b2a_ratio 
     113#    c_side = length_a * c2a_ratio 
     114# 
     115#    # surface average radius (rough approximation) 
     116#    surf_rad = sqrt(length_a * b_side / pi) 
     117# 
     118#    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
     119#    return 0.5 * (ddd) ** (1. / 3.) 
    118120 
    119121def VR(length_a, b2a_ratio, c2a_ratio): 

  • sasmodels/models/mono_gauss_coil.py

    r2d81cfe rd277229  
    5454 
    5555import numpy as np 
    56 from numpy import inf, exp, errstate 
     56from numpy import inf 
    5757 
    5858name = "mono_gauss_coil" 
    … …  
    6969parameters = [ 
    7070    ["i_zero", "1/cm", 70.0, [0.0, inf], "", "Intensity at q=0"], 
    71     ["rg", "Ang", 75.0, [0.0, inf], "", "Radius of gyration"], 
     71    ["rg", "Ang", 75.0, [0.0, inf], "volume", "Radius of gyration"], 
    7272    ] 
     73 
     74source = ["mono_gauss_coil.c"] 
     75have_Fq = False 
     76effective_radius_type = ["R_g","2R_g","3R_g","(5/3)^0.5*R_g"] 
     77 
     78 
    7379# pylint: enable=bad-whitespace, line-too-long 
    7480 
    75 # NB: Scale and Background are implicit parameters on every model 
    76 def Iq(q, i_zero, rg): 
    77     # pylint: disable = missing-docstring 
    78     z = (q * rg)**2 
    79  
    80     with errstate(invalid='ignore'): 
    81         inten = (i_zero * 2.0) * (exp(-z) + z - 1.0)/z**2 
    82         inten[q == 0] = i_zero 
    83     return inten 
    84 Iq.vectorized = True # Iq accepts an array of q values 
     81## NB: Scale and Background are implicit parameters on every model 
     82#def Iq(q, i_zero, rg): 
     83#    # pylint: disable = missing-docstring 
     84#    z = (q * rg)**2 
     85# 
     86#    with errstate(invalid='ignore'): 
     87#        inten = (i_zero * 2.0) * (exp(-z) + z - 1.0)/z**2 
     88#        inten[q == 0] = i_zero 
     89#    return inten 
     90#Iq.vectorized = True # Iq accepts an array of q values 
    8591 
    8692def random(): 

  • sasmodels/models/multilayer_vesicle.c

    r71b751d rd277229  
    4747} 
    4848 
     49static double 
     50effective_radius(int mode, double radius, double thick_shell, double thick_solvent, double fp_n_shells) 
     51{ 
     52    return radius + fp_n_shells*thick_shell + (fp_n_shells - 1.0)*thick_solvent; 
     53} 
     54 
    4955static void 
    5056Fq(double q, 

  • sasmodels/models/multilayer_vesicle.py

    r71b751d rd277229  
    146146source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 
    147147have_Fq = True 
     148effective_radius_type = ["outer radius"] 
    148149 
    149 def ER(radius, thick_shell, thick_solvent, n_shells): 
    150     n_shells = int(n_shells+0.5) 
    151     return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent 
     150#def ER(radius, thick_shell, thick_solvent, n_shells): 
     151#    n_shells = int(n_shells+0.5) 
     152#    return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent 
    152153 
    153154def random(): 

  • sasmodels/models/onion.c

    r71b751d rd277229  
    4040} 
    4141 
     42static double 
     43effective_radius(int mode, double radius_core, double n_shells, double thickness[]) 
     44{ 
     45  int n = (int)(n_shells+0.5); 
     46  double r = radius_core; 
     47  for (int i=0; i < n; i++) { 
     48    r += thickness[i]; 
     49  } 
     50  return r; 
     51} 
     52 
    4253static void 
    4354Fq(double q, double *F1, double *F2, double sld_core, double radius_core, double sld_solvent, 

  • sasmodels/models/onion.py

    r71b751d rd277229  
    316316single = False 
    317317have_Fq = True 
     318effective_radius_type = ["outer radius"] 
    318319 
    319320profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
    … …  
    366367    return np.asarray(z), np.asarray(rho) 
    367368 
    368 def ER(radius_core, n_shells, thickness): 
    369     """Effective radius""" 
    370     n = int(n_shells[0]+0.5) 
    371     return np.sum(thickness[:n], axis=0) + radius_core 
     369#def ER(radius_core, n_shells, thickness): 
     370#    """Effective radius""" 
     371#    n = int(n_shells[0]+0.5) 
     372#    return np.sum(thickness[:n], axis=0) + radius_core 
    372373 
    373374demo = { 

  • sasmodels/models/parallelepiped.c

    r71b751d rd277229  
    33{ 
    44    return length_a * length_b * length_c; 
     5} 
     6 
     7static double 
     8effective_radius(int mode, double length_a, double length_b, double length_c) 
     9{ 
     10    if (mode == 1) { 
     11        return cbrt(0.75*length_a*length_b*length_c/M_PI); 
     12    } else if (mode == 2) { 
     13        return 0.5 * length_a; 
     14    } else if (mode == 3) { 
     15        return 0.5 * length_b; 
     16    } else if (mode == 4) { 
     17        return 0.5 * length_c; 
     18    } else if (mode == 5) { 
     19        return sqrt(length_a*length_b/M_PI); 
     20    } else if (mode == 6) { 
     21        return 0.5*sqrt(length_a*length_a + length_b*length_b); 
     22    } else { 
     23        return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c); 
     24    } 
    525} 
    626 

  • sasmodels/models/parallelepiped.py

    r71b751d rd277229  
    231231source = ["lib/gauss76.c", "parallelepiped.c"] 
    232232have_Fq = True 
    233  
    234 def ER(length_a, length_b, length_c): 
    235     """ 
    236     Return effective radius (ER) for P(q)*S(q) 
    237     """ 
    238     # now that axes can be in any size order, need to sort a,b,c 
    239     # where a~b and c is either much smaller or much larger 
    240     abc = np.vstack((length_a, length_b, length_c)) 
    241     abc = np.sort(abc, axis=0) 
    242     selector = (abc[1] - abc[0]) > (abc[2] - abc[1]) 
    243     length = np.where(selector, abc[0], abc[2]) 
    244     # surface average radius (rough approximation) 
    245     radius = sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi) 
    246  
    247     ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius)) 
    248     return 0.5 * (ddd) ** (1. / 3.) 
     233effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
     234                         "equivalent circular cross-section","half ab diagonal","half diagonal"] 
     235 
     236#def ER(length_a, length_b, length_c): 
     237#    """ 
     238#    Return effective radius (ER) for P(q)*S(q) 
     239#    """ 
     240#    # now that axes can be in any size order, need to sort a,b,c 
     241#    # where a~b and c is either much smaller or much larger 
     242#    abc = np.vstack((length_a, length_b, length_c)) 
     243#    abc = np.sort(abc, axis=0) 
     244#    selector = (abc[1] - abc[0]) > (abc[2] - abc[1]) 
     245#    length = np.where(selector, abc[0], abc[2]) 
     246#    # surface average radius (rough approximation) 
     247#    radius = sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi) 
     248# 
     249#    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius)) 
     250#    return 0.5 * (ddd) ** (1. / 3.) 
    249251 
    250252# VR defaults to 1.0 

  • sasmodels/models/pringle.c

    r74768cb rd277229  
    104104} 
    105105 
     106static double 
     107effective_radius(int mode, double radius, double thickness, double alpha, double beta) 
     108{ 
     109    if (mode == 1) { 
     110        return cbrt(0.75*radius*radius*thickness); 
     111    } else { 
     112        return radius; 
     113    } 
     114} 
     115 
    106116double Iq( 
    107117    double q, 

  • sasmodels/models/pringle.py

    ref07e95 rd277229  
    7474source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", \ 
    7575          "lib/sas_JN.c", "lib/gauss76.c", "pringle.c"] 
     76effective_radius_type = ["equivalent sphere","radius"] 
    7677 
    77 def ER(radius, thickness, alpha, beta): 
    78     """ 
    79     Return equivalent radius (ER) 
    80     """ 
    81     ddd = 0.75 * radius * (2 * radius * thickness + (thickness + radius) \ 
    82                            * (thickness + pi * radius)) 
    83     return 0.5 * (ddd) ** (1. / 3.) 
     78#def ER(radius, thickness, alpha, beta): 
     79#    """ 
     80#    Return equivalent radius (ER) 
     81#    """ 
     82#    ddd = 0.75 * radius * (2 * radius * thickness + (thickness + radius) \ 
     83#                           * (thickness + pi * radius)) 
     84#    return 0.5 * (ddd) ** (1. / 3.) 
    8485 
    8586def random(): 

  • sasmodels/models/raspberry.c

    r71b751d rd277229  
    1 double form_volume(double radius_lg); 
     1double form_volume(double radius_lg, double radius_sm, double penetration); 
    22 
    33double Iq(double q, 
    … …  
    66          double radius_lg, double radius_sm, double penetration); 
    77 
    8 double form_volume(double radius_lg) 
     8double form_volume(double radius_lg, double radius_sm, double penetration) 
    99{ 
    1010    //Because of the complex structure, volume normalization must 
    … …  
    1212    double volume=1.0; 
    1313    return volume; 
     14} 
     15 
     16static double 
     17effective_radius(int mode, double radius_lg, double radius_sm, double penetration) 
     18{ 
     19    if (mode == 1) { 
     20        return radius_lg; 
     21    } else { 
     22        return radius_lg + 2.0*radius_sm - penetration; 
     23    } 
    1424} 
    1525 

  • sasmodels/models/raspberry.py

    ref07e95 rd277229  
    145145              ["radius_lg", "Ang", 5000, [0, inf], "volume", 
    146146               "radius of large spheres"], 
    147               ["radius_sm", "Ang", 100, [0, inf], "", 
     147              ["radius_sm", "Ang", 100, [0, inf], "volume", 
    148148               "radius of small spheres"], 
    149               ["penetration", "Ang", 0, [-1, 1], "", 
     149              ["penetration", "Ang", 0, [-1, 1], "volume", 
    150150               "fractional penetration depth of small spheres into large sphere"], 
    151151             ] 
    152152 
    153153source = ["lib/sas_3j1x_x.c", "raspberry.c"] 
     154effective_radius_type = ["radius_large","radius_outer"] 
    154155 
    155156def random(): 

  • sasmodels/models/rectangular_prism.c

    r71b751d rd277229  
    33{ 
    44    return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio); 
     5} 
     6 
     7static double 
     8effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio) 
     9{ 
     10    if (mode == 1) { 
     11        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
     12    } else if (mode == 2) { 
     13        return 0.5 * length_a; 
     14    } else if (mode == 3) { 
     15        return 0.5 * length_a*b2a_ratio; 
     16    } else if (mode == 4) { 
     17        return 0.5 * length_a*c2a_ratio; 
     18    } else if (mode == 5) { 
     19        return length_a*sqrt(b2a_ratio/M_PI); 
     20    } else if (mode == 6) { 
     21        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     22    } else { 
     23        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     24    } 
    525} 
    626 

  • sasmodels/models/rectangular_prism.py

    r71b751d rd277229  
    136136source = ["lib/gauss76.c", "rectangular_prism.c"] 
    137137have_Fq = True 
     138effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
     139                         "equivalent circular cross-section","half ab diagonal","half diagonal"] 
    138140 
    139 def ER(length_a, b2a_ratio, c2a_ratio): 
    140     """ 
    141         Return equivalent radius (ER) 
    142     """ 
    143     b_side = length_a * b2a_ratio 
    144     c_side = length_a * c2a_ratio 
    145  
    146     # surface average radius (rough approximation) 
    147     surf_rad = sqrt(length_a * b_side / pi) 
    148  
    149     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    150     return 0.5 * (ddd) ** (1. / 3.) 
     141#def ER(length_a, b2a_ratio, c2a_ratio): 
     142#    """ 
     143#        Return equivalent radius (ER) 
     144#    """ 
     145#    b_side = length_a * b2a_ratio 
     146#    c_side = length_a * c2a_ratio 
     147# 
     148#    # surface average radius (rough approximation) 
     149#    surf_rad = sqrt(length_a * b_side / pi) 
     150# 
     151#    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
     152#    return 0.5 * (ddd) ** (1. / 3.) 
    151153 
    152154def random(): 

  • sasmodels/models/sphere.py

    r71b751d rd277229  
    6767             ] 
    6868 
    69 source = ["lib/sas_3j1x_x.c"] 
     69source = ["lib/sas_3j1x_x.c","sphere.c"] 
    7070have_Fq = True 
     71effective_radius_type = ["radius"] 
    7172 
    72 c_code = """ 
    73 static double form_volume(double radius) 
    74 { 
    75     return M_4PI_3*cube(radius); 
    76 } 
    77  
    78 static void Fq(double q, double *f1, double *f2, double sld, double sld_solvent, double radius) 
    79 { 
    80     const double bes = sas_3j1x_x(q*radius); 
    81     const double contrast = (sld - sld_solvent); 
    82     const double form = contrast * form_volume(radius) * bes; 
    83     *f1 = 1.0e-2*form; 
    84     *f2 = 1.0e-4*form*form; 
    85 } 
    86 """ 
    87  
    88 def ER(radius): 
    89     """ 
    90     Return equivalent radius (ER) 
    91     """ 
    92     return radius 
     73#def ER(radius): 
     74#    """ 
     75#    Return equivalent radius (ER) 
     76#    """ 
     77#    return radius 
    9378 
    9479# VR defaults to 1.0 

  • sasmodels/models/spherical_sld.c

    r71b751d rd277229  
    1010    } 
    1111    return M_4PI_3*cube(r); 
     12} 
     13 
     14static double 
     15effective_radius(int mode, double fp_n_shells, double thickness[], double interface[]) 
     16{ 
     17    int n_shells= (int)(fp_n_shells + 0.5); 
     18    double r = 0.0; 
     19    for (int i=0; i < n_shells; i++) { 
     20        r += thickness[i] + interface[i]; 
     21    } 
     22    return r; 
    1223} 
    1324 

  • sasmodels/models/spherical_sld.py

    r71b751d rd277229  
    210210single = False  # TODO: fix low q behaviour 
    211211have_Fq = True 
     212effective_radius_type = ["outer radius"] 
    212213 
    213214profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
    … …  
    255256 
    256257 
    257 def ER(n_shells, thickness, interface): 
    258     """Effective radius""" 
    259     n_shells = int(n_shells + 0.5) 
    260     total = (np.sum(thickness[:n_shells], axis=1) 
    261              + np.sum(interface[:n_shells], axis=1)) 
    262     return total 
     258#def ER(n_shells, thickness, interface): 
     259#    """Effective radius""" 
     260#    n_shells = int(n_shells + 0.5) 
     261#    total = (np.sum(thickness[:n_shells], axis=1) 
     262#             + np.sum(interface[:n_shells], axis=1)) 
     263#    return total 
    263264 
    264265 

  • sasmodels/models/triaxial_ellipsoid.c

    r71b751d rd277229  
    77} 
    88 
     9static double 
     10radius_from_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     11{ 
     12    return cbrt(radius_equat_minor*radius_equat_major*radius_polar); 
     13} 
     14 
     15static double 
     16radius_from_min_dimension(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     17{ 
     18    const double rad_equat_min = (radius_equat_minor < radius_equat_major ? radius_equat_minor : radius_equat_major); 
     19    return (rad_equat_min < radius_polar ? rad_equat_min : radius_polar); 
     20} 
     21 
     22static double 
     23radius_from_max_dimension(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     24{ 
     25    const double rad_equat_max = (radius_equat_minor < radius_equat_major ? radius_equat_major : radius_equat_minor); 
     26    return (rad_equat_max > radius_polar ? rad_equat_max : radius_polar); 
     27} 
     28 
     29static double 
     30effective_radius(int mode, double radius_equat_minor, double radius_equat_major, double radius_polar) 
     31{ 
     32    if (mode == 1) { 
     33        return radius_from_volume(radius_equat_minor,radius_equat_major, radius_polar); 
     34    } else if (mode == 2) { 
     35        return radius_from_min_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
     36    } else { 
     37        return radius_from_max_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
     38    } 
     39} 
    940 
    1041static void 

  • sasmodels/models/triaxial_ellipsoid.py

    r71b751d rd277229  
    158158source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"] 
    159159have_Fq = True 
     160effective_radius_type = ["equivalent sphere","min radius", "max radius"] 
    160161 
    161162def ER(radius_equat_minor, radius_equat_major, radius_polar): 

  • sasmodels/models/vesicle.c

    r71b751d rd277229  
    66} 
    77 
     8static double 
     9effective_radius(int mode, double radius, double thickness) 
     10{ 
     11    return radius + thickness; 
     12} 
    813 
    914static void 

  • sasmodels/models/vesicle.py

    r71b751d rd277229  
    9595source = ["lib/sas_3j1x_x.c", "vesicle.c"] 
    9696have_Fq = True 
     97effective_radius_type = ["outer radius"] 
    9798 
    98 def ER(radius, thickness): 
    99     ''' 
    100     returns the effective radius used in the S*P calculation 
    101  
    102     :param radius: core radius 
    103     :param thickness: shell thickness 
    104     ''' 
    105     return radius + thickness 
     99#def ER(radius, thickness): 
     100#    ''' 
     101#    returns the effective radius used in the S*P calculation 
     102# 
     103#    :param radius: core radius 
     104#    :param thickness: shell thickness 
     105#    ''' 
     106#    return radius + thickness 
    106107 
    107108def VR(radius, thickness): 
Note: See TracChangeset for help on using the changeset viewer.