Changeset ee60aa7 in sasmodels


Ignore:
Timestamp:
Sep 10, 2018 2:16:46 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
d299327
Parents:
3f818b2
Message:

clean up effective radius functions; improve mono_gauss_coil accuracy; start moving VR into C

Files:
59 edited

Legend:

Unmodified
Added
Removed

  • explore/precision.py

    r2a7e20e ree60aa7  
    357357) 
    358358add_function( 
    359     name="debye", 
     359    name="gauss_coil", 
    360360    mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4, 
    361361    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4, 
    362362    ocl_function=make_ocl(""" 
    363363    const double qsq = q*q; 
    364     if (qsq < 1.0) { // Pade approximation 
     364    // For double: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide) 
     365    // For single: use O(7) Taylor with 0.8 cutoff (7 mad) 
     366    if (qsq < 0.0) { 
    365367        const double x = qsq; 
    366368        if (0) { // 0.36 single 
    … …  
    372374            const double B1=3./8., B2=3./56., B3=1./336.; 
    373375            return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.); 
    374         } else if (1) { // 1.0 for single, 0.25 for double 
     376        } else if (0) { // 1.0 for single, 0.25 for double 
    375377            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
    376378            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 
    … …  
    385387                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 
    386388        } 
    387     } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double 
     389    } else if (qsq < 0.8) { 
    388390        const double x = qsq; 
    389391        const double C0 = +1.; 

  • sasmodels/core.py

    r5399809 ree60aa7  
    364364    actual = [p.name for p in model.info.parameters.kernel_parameters] 
    365365    target = ("sld sld_solvent radius length theta phi" 
    366               " radius_effective volfraction structure_factor_mode" 
     366              " radius_effective volfraction " 
     367              " structure_factor_mode radius_effective_mode" 
    367368              " A_sld A_sld_solvent A_radius").split() 
    368369    assert target == actual, "%s != %s"%(target, actual) 

  • sasmodels/kernel.py

    r5399809 ree60aa7  
    9292    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    9393        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    94         raise NotImpmentedError() 
     94        raise NotImplementedError() 

  • sasmodels/models/barbell.c

    rd277229 ree60aa7  
    7979effective_radius(int mode, double radius_bell, double radius, double length) 
    8080{ 
    81     if (mode == 1) { 
     81    switch (mode) { 
     82    case 1: // equivalent sphere 
    8283        return radius_from_volume(radius_bell, radius , length); 
    83     } else if (mode == 2) { 
     84    case 2: // radius 
    8485        return radius; 
    85     } else if (mode == 3) { 
     86    case 3: // half length 
    8687        return 0.5*length; 
    87     } else { 
     88    case 4: // half total length 
    8889        return radius_from_totallength(radius_bell,radius,length); 
    8990    } 

  • sasmodels/models/barbell.py

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

  • sasmodels/models/capped_cylinder.c

    rd277229 ree60aa7  
    101101effective_radius(int mode, double radius, double radius_cap, double length) 
    102102{ 
    103     if (mode == 1) { 
     103    switch (mode) { 
     104    case 1: // equivalent sphere 
    104105        return radius_from_volume(radius, radius_cap, length); 
    105     } else if (mode == 2) { 
     106    case 2: // radius 
    106107        return radius; 
    107     } else if (mode == 3) { 
     108    case 3: // half length 
    108109        return 0.5*length; 
    109     } else { 
     110    case 4: // half total length 
    110111        return radius_from_totallength(radius, radius_cap,length); 
    111112    } 

  • sasmodels/models/capped_cylinder.py

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

  • sasmodels/models/core_multi_shell.c

    ra94046f ree60aa7  
    66  const double vol = M_4PI_3 * cube(r); 
    77  return sld * vol * bes; 
    8 } 
    9  
    10 static double 
    11 form_volume(double core_radius, double fp_n, double thickness[]) 
    12 { 
    13   double r = core_radius; 
    14   int n = (int)(fp_n+0.5); 
    15   for (int i=0; i < n; i++) { 
    16     r += thickness[i]; 
    17   } 
    18   return M_4PI_3 * cube(r); 
    198} 
    209 
    … …  
    3120 
    3221static double 
     22form_volume(double core_radius, double fp_n, double thickness[]) 
     23{ 
     24  return M_4PI_3 * cube(outer_radius(core_radius, fp_n, thickness)); 
     25} 
     26 
     27static double 
    3328effective_radius(int mode, double core_radius, double fp_n, double thickness[]) 
    34 // this seems regardless to always give the result for outer radius for n=1 shells; why?? 
    35 // printf shows fp_n is always 1, not 0,1,2 
    3629{ 
    37 //        printf("fp_n =%g \n",fp_n); 
    38         if (mode == 1) { 
    39         double r = core_radius; 
    40         int n = (int)(fp_n+0.5); 
    41         if ( n > 0) { 
    42             for (int i=0; i < n; i++) { 
    43                 r += thickness[i]; 
    44             } 
    45         } 
    46         return r; 
    47         //return outer_radius(core_radius,fp_n,thickness); 
    48     } else { 
    49         return core_radius; 
    50     } 
     30  switch (mode) { 
     31  case 1: // outer radius 
     32    return outer_radius(core_radius, fp_n, thickness); 
     33  case 2: // core radius 
     34    return core_radius; 
     35  } 
    5136} 
    5237 

  • sasmodels/models/core_multi_shell.py

    rd277229 ree60aa7  
    145145    return np.asarray(z), np.asarray(rho) 
    146146 
    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 
    151  
    152147demo = dict(sld_core=6.4, 
    153148            radius=60, 

  • sasmodels/models/core_shell_bicelle.c

    rd277229 ree60aa7  
    5454effective_radius(int mode, double radius, double thick_rim, double thick_face, double length) 
    5555{ 
    56     if (mode == 1) { 
     56    switch (mode) { 
     57    case 1: // equivalent sphere 
    5758        return radius_from_volume(radius, thick_rim, thick_face, length); 
    58     } else if (mode == 2) { 
     59    case 2: // outer rim radius 
    5960        return radius + thick_rim; 
    60     } else if (mode == 3) { 
     61    case 3: // half outer thickness 
    6162        return 0.5*length + thick_face; 
    62     } else { 
     63    case 4: // half diagonal 
    6364        return radius_from_diagonal(radius,thick_rim,thick_face,length); 
    6465    } 

  • sasmodels/models/core_shell_bicelle.py

    rd277229 ree60aa7  
    155155          "core_shell_bicelle.c"] 
    156156have_Fq = True 
    157 effective_radius_type = ["equivalent sphere","outer rim radius", "half outer thickness","half diagonal"] 
     157effective_radius_type = [ 
     158    "equivalent sphere", "outer rim radius", 
     159    "half outer thickness", "half diagonal", 
     160    ] 
    158161 
    159162def random(): 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    rd277229 ree60aa7  
    2929effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
    3030{ 
    31     if (mode == 1) { 
     31    switch (mode) { 
     32    case 1: // equivalent sphere 
    3233        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
    33     } else if (mode == 2) { 
     34    case 2: // outer rim average radius 
    3435        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
    35     } else if (mode == 3) { 
     36    case 3: // outer rim min radius 
    3637        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
    37     } else if (mode == 4) { 
     38    case 4: // outer max radius 
    3839        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
    39     } else if (mode ==5) { 
     40    case 5: // half outer thickness 
    4041        return 0.5*length + thick_face; 
    41     } else { 
     42    case 6: // half diagonal 
    4243        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
    4344    } 

  • sasmodels/models/core_shell_bicelle_elliptical.py

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

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    rd277229 ree60aa7  
    3030effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
    3131{ 
    32     if (mode == 1) { 
     32    switch (mode) { 
     33    case 1: // equivalent sphere 
    3334        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
    34     } else if (mode == 2) { 
     35    case 2: // outer rim average radius 
    3536        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
    36     } else if (mode == 3) { 
     37    case 3: // outer rim min radius 
    3738        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
    38     } else if (mode == 4) { 
     39    case 4: // outer max radius 
    3940        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
    40     } else if (mode ==5) { 
     41    case 5: // half outer thickness 
    4142        return 0.5*length + thick_face; 
    42     } else { 
     43    case 6: // half diagonal 
    4344        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
    4445    } 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

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

  • sasmodels/models/core_shell_cylinder.c

    ra94046f ree60aa7  
    3030static double 
    3131effective_radius(int mode, double radius, double thickness, double length) 
    32 //effective_radius_type = ["equivalent sphere","outer radius","half outer length","half min outer dimension", 
    33 //                         "half max outer dimension","half outer diagonal"] 
    3432{ 
    35     if (mode == 1) { 
     33    switch (mode) { 
     34    case 1: // equivalent sphere 
    3635        return radius_from_volume(radius, thickness, length); 
    37     } else if (mode == 2) { 
     36    case 2: // outer radius 
    3837        return radius + thickness; 
    39     } else if (mode == 3) { 
     38    case 3: // half outer length 
    4039        return 0.5*length + thickness; 
    41     } else if (mode == 4) { 
     40    case 4: // half min outer length 
    4241        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness); 
    43     } else if (mode == 5) { 
     42    case 5: // half max outer length 
    4443        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness); 
    45     } else { 
     44    case 6: // half outer diagonal 
    4645        return radius_from_diagonal(radius,thickness,length); 
    4746    } 

  • sasmodels/models/core_shell_cylinder.py

    rc44b611 ree60aa7  
    132132source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] 
    133133have_Fq = True 
    134 effective_radius_type = ["equivalent sphere","outer radius","half outer length","half min outer dimension", 
    135                          "half max outer dimension","half outer diagonal"] 
    136  
    137 #def ER(radius, thickness, length): 
    138 #    """ 
    139 #    Returns the effective radius used in the S*P calculation 
    140 #    """ 
    141 #    radius = radius + thickness 
    142 #    length = length + 2 * thickness 
    143 #    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    144 #    return 0.5 * (ddd) ** (1. / 3.) 
     134effective_radius_type = [ 
     135    "equivalent sphere", "outer radius", "half outer length", 
     136    "half min outer dimension", "half max outer dimension", 
     137    "half outer diagonal", 
     138    ] 
    145139 
    146140def random(): 

  • sasmodels/models/core_shell_ellipsoid.c

    r3c60146 ree60aa7  
    7373    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    7474 
    75     if (mode == 1) { 
     75    switch (mode) { 
     76    case 1: // equivalent sphere 
    7677        return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 
    77     } else if (mode == 2) { 
     78    case 2: // average outer curvature 
    7879        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell); 
    79     } else if (mode == 3) { 
     80    case 3: // min outer radius 
    8081        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
    81     } else { 
     82    case 4: // max outer radius 
    8283        return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
    8384    } 

  • sasmodels/models/core_shell_ellipsoid.py

    rd277229 ree60aa7  
    146146source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    147147have_Fq = True 
    148 effective_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) 
     148effective_radius_type = [ 
     149    "equivalent sphere", "average outer curvature", 
     150    "min outer radius", "max outer radius", 
     151    ] 
    158152 
    159153def random(): 

  • sasmodels/models/core_shell_parallelepiped.c

    ra94046f ree60aa7  
    3131                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    3232{ 
    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); 
     33    const double volume = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     34    return cbrt(volume/M_4PI_3); 
    3535} 
    3636 
    … …  
    4545effective_radius(int mode, double length_a, double length_b, double length_c, 
    4646                 double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    47 //effective_radius_type = ["equivalent sphere","half outer length_a", "half outer length_b", "half outer length_c", 
    48 //                         "equivalent circular cross-section","half outer ab diagonal","half outer diagonal"] 
    49 // note the core box is A*B*C with slabs ta, tb & tc on each face but missing the corners, though that fact is ignored here 
    50 // in the equvalent sphere option 
    5147{ 
    52     if (mode == 1) { 
     48    switch (mode) { 
     49    case 1: // equivalent sphere 
    5350        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
    54     } else if (mode == 2) { 
     51    case 2: // half outer length a 
    5552        return 0.5 * length_a + thick_rim_a; 
    56     } else if (mode == 3) { 
     53    case 3: // half outer length b 
    5754        return 0.5 * length_b + thick_rim_b; 
    58     } else if (mode == 4) { 
     55    case 4: // half outer length c 
    5956        return 0.5 * length_c + thick_rim_c; 
    60     } else if (mode == 5) { 
     57    case 5: // equivalent circular cross-section 
    6158        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b); 
    62     } else if (mode == 6) { 
     59    case 6: // half outer ab diagonal 
    6360        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b)); 
    64     } else { 
     61    case 7: // half outer diagonal 
    6562        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b) + square(length_c+ 2.0*thick_rim_c)); 
    6663    } 

  • sasmodels/models/core_shell_parallelepiped.py

    rd277229 ree60aa7  
    9898is the scattering length of the solvent. 
    9999 
    100 .. note::  
     100.. note:: 
    101101 
    102102   the code actually implements two substitutions: $d(cos\alpha)$ is 
    … …  
    227227source = ["lib/gauss76.c", "core_shell_parallelepiped.c"] 
    228228have_Fq = True 
    229 effective_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) 
    242  
    243 # VR defaults to 1.0 
     229effective_radius_type = [ 
     230    "equivalent sphere", 
     231    "half outer length_a", "half outer length_b", "half outer length_c", 
     232    "equivalent circular cross-section", 
     233    "half outer ab diagonal", "half outer diagonal", 
     234    ] 
    244235 
    245236def random(): 

  • sasmodels/models/core_shell_sphere.c

    rd277229 ree60aa7  
    88effective_radius(int mode, double radius, double thickness) 
    99{ 
    10     if (mode == 1) { 
     10    switch (mode) { 
     11    case 1: // outer radius 
    1112        return radius + thickness; 
    12     } else { 
     13    case 2: // core radius 
    1314        return radius; 
    1415    } 

  • sasmodels/models/core_shell_sphere.py

    r2cc8aa2 ree60aa7  
    8383            sld_core=1.0, sld_shell=2.0, sld_solvent=0.0) 
    8484 
    85 #def ER(radius, thickness): 
    86 #    """ 
    87 #        Equivalent radius 
    88 #        @param radius: core radius 
    89 #        @param thickness: shell thickness 
    90 #    """ 
    91 #    return radius + thickness 
    92  
    9385def random(): 
    9486    outer_radius = 10**np.random.uniform(1.3, 4.3) 

  • sasmodels/models/cylinder.c

    rd277229 ree60aa7  
    2828effective_radius(int mode, double radius, double length) 
    2929{ 
    30     if (mode == 1) { 
     30    switch (mode) { 
     31    case 1: 
    3132        return radius_from_volume(radius, length); 
    32     } else if (mode == 2) { 
     33    case 2: 
    3334        return radius; 
    34     } else if (mode == 3) { 
     35    case 3: 
    3536        return 0.5*length; 
    36     } else if (mode == 4) { 
     37    case 4: 
    3738        return (radius < 0.5*length ? radius : 0.5*length); 
    38     } else if (mode == 5) { 
     39    case 5: 
    3940        return (radius > 0.5*length ? radius : 0.5*length); 
    40     } else { 
     41    case 6: 
    4142        return radius_from_diagonal(radius,length); 
    4243    } 

  • sasmodels/models/cylinder.py

    rd277229 ree60aa7  
    139139source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 
    140140have_Fq = True 
    141 effective_radius_type = ["equivalent sphere","radius","half length","half min dimension","half max dimension","half diagonal"] 
    142  
    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.) 
     141effective_radius_type = [ 
     142    "equivalent sphere", "radius", 
     143    "half length", "half min dimension", "half max dimension", "half diagonal", 
     144    ] 
    149145 
    150146def random(): 

  • sasmodels/models/ellipsoid.c

    rd277229 ree60aa7  
    3434effective_radius(int mode, double radius_polar, double radius_equatorial) 
    3535{ 
    36     if (mode == 1) { 
     36    switch (mode) { 
     37    case 1: // equivalent sphere 
    3738        return radius_from_volume(radius_polar, radius_equatorial); 
    38     } else if (mode == 2) { 
     39    case 2: // average curvature 
    3940        return radius_from_curvature(radius_polar, radius_equatorial); 
    40     } else if (mode == 3) { 
     41    case 3: // min radius 
    4142        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial); 
    42     } else { 
     43    case 4: // max radius 
    4344        return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial); 
    4445    } 

  • sasmodels/models/ellipsoid.py

    r2773c66 ree60aa7  
    169169source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
    170170have_Fq = True 
    171 effective_radius_type = ["equivalent sphere","average curvature", "min radius", "max radius"] 
     171effective_radius_type = [ 
     172    "equivalent sphere", "average curvature", "min radius", "max radius", 
     173    ] 
    172174 
    173175def random(): 

  • sasmodels/models/elliptical_cylinder.c

    rfbaef04 ree60aa7  
    3535static double 
    3636effective_radius(int mode, double radius_minor, double r_ratio, double length) 
    37 //effective_radius_type = ["equivalent sphere","average radius","min radius","max radius", 
    38 //                         "equivalent circular cross-section","half length","half min dimension","half max dimension","half diagonal"] 
    3937{ 
    40     if (mode == 1) { 
     38    switch (mode) { 
     39    case 1: // equivalent sphere 
    4140        return radius_from_volume(radius_minor, r_ratio, length); 
    42     } else if (mode == 2) { 
     41    case 2: // average radius 
    4342        return 0.5*radius_minor*(1.0 + r_ratio); 
    44     } else if (mode == 3) { 
     43    case 3: // min radius 
    4544        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
    46     } else if (mode == 4) { 
     45    case 4: // max radius 
    4746        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
    48     } else if (mode == 5) { 
     47    case 5: // equivalent circular cross-section 
    4948        return sqrt(radius_minor*radius_minor*r_ratio); 
    50     } else if (mode == 6) { 
     49    case 6: // half length 
    5150        return 0.5*length; 
    52     } else if (mode == 7) { 
     51    case 7: // half min dimension 
    5352        return radius_from_min_dimension(radius_minor,r_ratio,0.5*length); 
    54     } else if (mode == 8) { 
     53    case 8: // half max dimension 
    5554        return radius_from_max_dimension(radius_minor,r_ratio,0.5*length); 
    56     } else { 
     55    case 9: // half diagonal 
    5756        return radius_from_diagonal(radius_minor,r_ratio,length); 
    5857    } 

  • sasmodels/models/elliptical_cylinder.py

    r2cc8aa2 ree60aa7  
    123123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
    124124have_Fq = True 
    125 effective_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"] 
     125effective_radius_type = [ 
     126    "equivalent sphere", "average radius", "min radius", "max radius", 
     127    "equivalent circular cross-section", 
     128    "half length", "half min dimension", "half max dimension", "half diagonal", 
     129    ] 
    127130 
    128131demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 
    129132            sld=4.0, sld_solvent=1.0, theta=10.0, phi=20, psi=30, 
    130133            theta_pd=10, phi_pd=2, psi_pd=3) 
    131  
    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.) 
    143134 
    144135def random(): 

  • sasmodels/models/fuzzy_sphere.c

    rd277229 ree60aa7  
    77effective_radius(int mode, double radius, double fuzziness) 
    88{ 
    9     if (mode == 1) { 
     9    switch (mode) { 
     10    case 1: // radius 
    1011        return radius; 
    11     } else { 
     12    case 2: // radius + fuzziness 
    1213        return radius + fuzziness; 
    1314    } 

  • sasmodels/models/fuzzy_sphere.py

    rd277229 ree60aa7  
    8484source = ["lib/sas_3j1x_x.c","fuzzy_sphere.c"] 
    8585have_Fq = True 
    86 effective_radius_type = ["radius","radius + fuzziness"] 
    87  
    88 #def ER(radius): 
    89 #    """ 
    90 #    Return radius 
    91 #    """ 
    92 #    return radius 
    93  
    94 # VR defaults to 1.0 
     86effective_radius_type = ["radius", "radius + fuzziness"] 
    9587 
    9688def random(): 

  • sasmodels/models/hollow_cylinder.c

    rd277229 ree60aa7  
    1515} 
    1616 
     17// TODO: interface to form_volume/shell_volume not yet settled 
     18static double 
     19shell_volume(double *total, double radius, double thickness, double length) 
     20{ 
     21    *total = M_PI*length*square(radius+thickness); 
     22    return *total - M_PI*length*radius*radius; 
     23} 
     24 
    1725static double 
    1826form_volume(double radius, double thickness, double length) 
    1927{ 
    20     double v_shell = M_PI*length*(square(radius+thickness) - radius*radius); 
    21     return v_shell; 
     28    double total; 
     29    return shell_volume(&total, radius, thickness, length); 
    2230} 
    2331 
    … …  
    3846effective_radius(int mode, double radius, double thickness, double length) 
    3947{ 
    40     if (mode == 1) { 
     48    switch (mode) { 
     49    case 1: // equivalent sphere 
    4150        return radius_from_volume(radius, thickness, length); 
    42     } else if (mode == 2) { 
     51    case 2: // outer radius 
    4352        return radius + thickness; 
    44     } else if (mode == 3) { 
     53    case 3: // half length 
    4554        return 0.5*length; 
    46     } else if (mode == 4) { 
     55    case 4: // half outer min dimension 
    4756        return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length); 
    48     } else if (mode == 5) { 
     57    case 5: // half outer max dimension 
    4958        return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length); 
    50     } else { 
     59    case 6: // half outer diagonal 
    5160        return radius_from_diagonal(radius,thickness,length); 
    5261    } 

  • sasmodels/models/hollow_cylinder.py

    r2cc8aa2 ree60aa7  
    100100source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
    101101have_Fq = True 
    102 effective_radius_type = ["equivalent sphere","outer radius","half length", 
    103                          "half outer min dimension","half outer max dimension","half outer diagonal"] 
    104  
    105 # pylint: disable=W0613 
    106 #def ER(radius, thickness, length): 
    107 #    """ 
    108 #    :param radius:      Cylinder core radius 
    109 #    :param thickness:   Cylinder wall thickness 
    110 #    :param length:      Cylinder length 
    111 #    :return:            Effective radius 
    112 #    """ 
    113 #    router = radius + thickness 
    114 #    if router == 0 or length == 0: 
    115 #        return 0.0 
    116 #    len1 = router 
    117 #    len2 = length/2.0 
    118 #    term1 = len1*len1*2.0*len2/2.0 
    119 #    term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0) 
    120 #    ddd = 3.0*term1*term2 
    121 #    diam = pow(ddd, (1.0/3.0)) 
    122 #    return diam 
    123  
    124 def VR(radius, thickness, length): 
    125     """ 
    126     :param radius:      Cylinder radius 
    127     :param thickness:   Cylinder wall thickness 
    128     :param length:      Cylinder length 
    129     :return:            Volf ratio for P(q)*S(q) 
    130     """ 
    131     router = radius + thickness 
    132     vol_core = pi*radius*radius*length 
    133     vol_total = pi*router*router*length 
    134     vol_shell = vol_total - vol_core 
    135     return vol_total, vol_shell 
     102effective_radius_type = [ 
     103    "equivalent sphere", "outer radius", "half length", 
     104    "half outer min dimension", "half outer max dimension", 
     105    "half outer diagonal", 
     106    ] 
    136107 
    137108def random(): 

  • sasmodels/models/hollow_rectangular_prism.c

    ra94046f ree60aa7  
     1// TODO: interface to form_volume/shell_volume not yet settled 
    12static double 
    2 form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     3shell_volume(double *total, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    34{ 
    45    double length_b = length_a * b2a_ratio; 
    … …  
    89    double c_core = length_c - 2.0*thickness; 
    910    double vol_core = a_core * b_core * c_core; 
    10     double vol_total = length_a * length_b * length_c; 
    11     double vol_shell = vol_total - vol_core; 
    12     return vol_shell; 
     11    *total = length_a * length_b * length_c; 
     12    return *total - vol_core; 
     13} 
     14 
     15static double 
     16form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     17{ 
     18    double total; 
     19    return shell_volume(&total, length_a, b2a_ratio, c2a_ratio, thickness); 
    1320} 
    1421 
    … …  
    1926// NOTE length_a is external dimension! 
    2027{ 
    21     if (mode == 1) { 
     28    switch (mode) { 
     29    case 1: // equivalent sphere 
    2230        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
    23     } else if (mode == 2) { 
     31    case 2: // half length_a 
    2432        return 0.5 * length_a; 
    25     } else if (mode == 3) { 
     33    case 3: // half length_b 
    2634        return 0.5 * length_a*b2a_ratio; 
    27     } else if (mode == 4) { 
     35    case 4: // half length_c 
    2836        return 0.5 * length_a*c2a_ratio; 
    29     } else if (mode == 5) { 
     37    case 5: // equivalent outer circular cross-section 
    3038        return length_a*sqrt(b2a_ratio/M_PI); 
    31     } else if (mode == 6) { 
     39    case 6: // half ab diagonal 
    3240        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
    33     } else { 
     41    case 7: // half diagonal 
    3442        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
    3543    } 

  • sasmodels/models/hollow_rectangular_prism.py

    rc44b611 ree60aa7  
    149149source = ["lib/gauss76.c", "hollow_rectangular_prism.c"] 
    150150have_Fq = True 
    151 effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
    152                          "equivalent outer circular cross-section","half ab diagonal","half diagonal"] 
    153  
    154 #def ER(length_a, b2a_ratio, c2a_ratio, thickness): 
    155 #    """ 
    156 #    Return equivalent radius (ER) 
    157 #    thickness parameter not used 
    158 #    """ 
    159 #    b_side = length_a * b2a_ratio 
    160 #    c_side = length_a * c2a_ratio 
    161 # 
    162 #    # surface average radius (rough approximation) 
    163 #    surf_rad = sqrt(length_a * b_side / pi) 
    164 # 
    165 #    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    166 #    return 0.5 * (ddd) ** (1. / 3.) 
    167  
    168 def VR(length_a, b2a_ratio, c2a_ratio, thickness): 
    169     """ 
    170     Return shell volume and total volume 
    171     """ 
    172     b_side = length_a * b2a_ratio 
    173     c_side = length_a * c2a_ratio 
    174     a_core = length_a - 2.0*thickness 
    175     b_core = b_side - 2.0*thickness 
    176     c_core = c_side - 2.0*thickness 
    177     vol_core = a_core * b_core * c_core 
    178     vol_total = length_a * b_side * c_side 
    179     vol_shell = vol_total - vol_core 
    180     return vol_total, vol_shell 
    181  
     151effective_radius_type = [ 
     152    "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     153    "equivalent outer circular cross-section", 
     154    "half ab diagonal", "half diagonal", 
     155    ] 
    182156 
    183157def random(): 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    rd277229 ree60aa7  
     1// TODO: interface to form_volume/shell_volume not yet settled 
    12static double 
    2 form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     3shell_volume(double *total, double length_a, double b2a_ratio, double c2a_ratio) 
    34{ 
    45    double length_b = length_a * b2a_ratio; 
    56    double length_c = length_a * c2a_ratio; 
    67    double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
     8    *total = length_a * length_b * length_c; 
    79    return vol_shell; 
    810} 
    911 
    1012static double 
     13form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     14{ 
     15    double total; 
     16    return shell_volume(&total, length_a, b2a_ratio, c2a_ratio); 
     17} 
     18 
     19 
     20static double 
    1121effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio) 
    1222{ 
    13     if (mode == 1) { 
     23    switch (mode) { 
     24    case 1: // equivalent sphere 
    1425        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
    15     } else if (mode == 2) { 
     26    case 2: // half length_a 
    1627        return 0.5 * length_a; 
    17     } else if (mode == 3) { 
     28    case 3: // half length_b 
    1829        return 0.5 * length_a*b2a_ratio; 
    19     } else if (mode == 4) { 
     30    case 4: // half length_c 
    2031        return 0.5 * length_a*c2a_ratio; 
    21     } else if (mode == 5) { 
     32    case 5: // equivalent outer circular cross-section 
    2233        return length_a*sqrt(b2a_ratio/M_PI); 
    23     } else if (mode == 6) { 
     34    case 6: // half ab diagonal 
    2435        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
    25     } else { 
     36    case 7: // half diagonal 
    2637        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
    2738    } 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

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

  • sasmodels/models/mono_gauss_coil.c

    rd277229 ree60aa7  
    77effective_radius(int mode, double rg) 
    88{ 
    9     if (mode == 1) { 
     9    switch (mode) { 
     10    case 1: // R_g 
    1011        return rg; 
    11     } else if (mode == 2) { 
     12    case 2: // 2R_g 
    1213        return 2.0*rg; 
    13     } else if (mode == 3) { 
     14    case 3: // 3R_g 
    1415        return 3.0*rg; 
    15     } else { 
     16    case 4: // (5/3)^0.5*R_g 
    1617        return sqrt(5.0/3.0)*rg; 
    1718    } 
    1819} 
    1920 
     21static double 
     22gauss_coil(double qr) 
     23{ 
     24    const double x = qr*qr; 
     25 
     26    // Use series expansion at low q for higher accuracy. We could use 
     27    // smaller polynomials if we sacrifice some digits of precision or 
     28    // introduce an additional series expansion around x == 1. 
     29    // See explore/precision.py, gauss_coil function. 
     30#if FLOAT_SIZE>4 // DOUBLE_PRECISION 
     31    // For double precision: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide) 
     32    if (x < 0.5) { 
     33        // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
     34        const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.; 
     35        const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.; 
     36        return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.) 
     37                / (((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 
     38    } 
     39#else 
     40    // For single precision: use O(7) Taylor with 0.8 cutoff (7 mad) 
     41    if (x < 0.8) { 
     42        const double C0 = +1.; 
     43        const double C1 = -1./3.; 
     44        const double C2 = +1./12.; 
     45        const double C3 = -1./60.; 
     46        const double C4 = +1./360.; 
     47        const double C5 = -1./2520.; 
     48        const double C6 = +1./20160.; 
     49        const double C7 = -1./181440.; 
     50        return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
     51    } 
     52#endif 
     53 
     54    return 2.0 * (expm1(-x) + x)/(x*x); 
     55} 
     56 
    2057double Iq(double q, double i_zero, double rg) 
    2158{ 
    22     const double uarg = square(q*rg); 
    23     const double inten; 
    24     if (q == 0) { 
    25         inten = i_zero; 
    26     } else { 
    27         inten = 2.0*i_zero * (exp(-uarg) + uarg - 1.0)/square(uarg); 
    28     } 
    29  
    30     return inten; 
     59    return i_zero * gauss_coil(q*rg); 
    3160} 

  • sasmodels/models/mono_gauss_coil.py

    rd277229 ree60aa7  
    7171    ["rg", "Ang", 75.0, [0.0, inf], "volume", "Radius of gyration"], 
    7272    ] 
     73# pylint: enable=bad-whitespace, line-too-long 
    7374 
    7475source = ["mono_gauss_coil.c"] 
    7576have_Fq = False 
    76 effective_radius_type = ["R_g","2R_g","3R_g","(5/3)^0.5*R_g"] 
     77effective_radius_type = ["R_g", "2R_g", "3R_g", "(5/3)^0.5*R_g"] 
    7778 
    78  
    79 # pylint: enable=bad-whitespace, line-too-long 
    80  
    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 
    9179 
    9280def random(): 

  • sasmodels/models/multilayer_vesicle.c

    rd277229 ree60aa7  
    5050effective_radius(int mode, double radius, double thick_shell, double thick_solvent, double fp_n_shells) 
    5151{ 
     52    // case 1: outer radius 
    5253    return radius + fp_n_shells*thick_shell + (fp_n_shells - 1.0)*thick_solvent; 
    5354} 

  • sasmodels/models/multilayer_vesicle.py

    rd277229 ree60aa7  
    148148effective_radius_type = ["outer radius"] 
    149149 
    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 
    153  
    154150def random(): 
    155151    volfraction = 10**np.random.uniform(-3, -0.5)  # scale from 0.1% to 30% 

  • sasmodels/models/onion.c

    rd277229 ree60aa7  
    3030 
    3131static double 
    32 form_volume(double radius_core, double n_shells, double thickness[]) 
    33 { 
    34   int n = (int)(n_shells+0.5); 
    35   double r = radius_core; 
    36   for (int i=0; i < n; i++) { 
    37     r += thickness[i]; 
    38   } 
    39   return M_4PI_3*cube(r); 
    40 } 
    41  
    42 static double 
    43 effective_radius(int mode, double radius_core, double n_shells, double thickness[]) 
     32outer_radius(double radius_core, double n_shells, double thickness[]) 
    4433{ 
    4534  int n = (int)(n_shells+0.5); 
    … …  
    4938  } 
    5039  return r; 
     40} 
     41 
     42static double 
     43form_volume(double radius_core, double n_shells, double thickness[]) 
     44{ 
     45  return M_4PI_3*cube(outer_radius(radius_core, n_shells, thickness)); 
     46} 
     47 
     48static double 
     49effective_radius(int mode, double radius_core, double n_shells, double thickness[]) 
     50{ 
     51  // case 1: outer radius 
     52  return outer_radius(radius_core, n_shells, thickness); 
    5153} 
    5254 

  • sasmodels/models/onion.py

    rd277229 ree60aa7  
    367367    return np.asarray(z), np.asarray(rho) 
    368368 
    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 
    373  
    374369demo = { 
    375370    "sld_solvent": 2.2, 

  • sasmodels/models/parallelepiped.c

    rd277229 ree60aa7  
    88effective_radius(int mode, double length_a, double length_b, double length_c) 
    99{ 
    10     if (mode == 1) { 
     10    switch (mode) { 
     11    case 1: // equivalent sphere 
    1112        return cbrt(0.75*length_a*length_b*length_c/M_PI); 
    12     } else if (mode == 2) { 
     13    case 2: // half length_a 
    1314        return 0.5 * length_a; 
    14     } else if (mode == 3) { 
     15    case 3: // half length_b 
    1516        return 0.5 * length_b; 
    16     } else if (mode == 4) { 
     17    case 4: // half length_c 
    1718        return 0.5 * length_c; 
    18     } else if (mode == 5) { 
     19    case 5: // equivalent circular cross-section 
    1920        return sqrt(length_a*length_b/M_PI); 
    20     } else if (mode == 6) { 
     21    case 6: // half ab diagonal 
    2122        return 0.5*sqrt(length_a*length_a + length_b*length_b); 
    22     } else { 
     23    case 7: // half diagonal 
    2324        return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c); 
    2425    } 
    2526} 
    26  
    2727 
    2828static void 

  • sasmodels/models/parallelepiped.py

    rd277229 ree60aa7  
    140140   ensuring that the inequality $A < B < C$ is not violated,  The calculation 
    141141   will not report an error, but the results may be not correct. 
    142     
     142 
    143143.. _parallelepiped-orientation: 
    144144 
    … …  
    231231source = ["lib/gauss76.c", "parallelepiped.c"] 
    232232have_Fq = True 
    233 effective_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.) 
    251  
    252 # VR defaults to 1.0 
    253  
     233effective_radius_type = [ 
     234    "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     235    "equivalent circular cross-section", "half ab diagonal", "half diagonal", 
     236    ] 
    254237 
    255238def random(): 

  • sasmodels/models/pringle.c

    rd277229 ree60aa7  
    107107effective_radius(int mode, double radius, double thickness, double alpha, double beta) 
    108108{ 
    109     if (mode == 1) { 
     109    switch (mode) { 
     110    case 1: // equivalent sphere 
    110111        return cbrt(0.75*radius*radius*thickness); 
    111     } else { 
     112    case 2: // radius 
    112113        return radius; 
    113114    } 

  • sasmodels/models/pringle.py

    rd277229 ree60aa7  
    7272 
    7373 
    74 source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", \ 
     74source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", 
    7575          "lib/sas_JN.c", "lib/gauss76.c", "pringle.c"] 
    76 effective_radius_type = ["equivalent sphere","radius"] 
    77  
    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.) 
     76effective_radius_type = ["equivalent sphere", "radius"] 
    8577 
    8678def random(): 

  • sasmodels/models/raspberry.c

    rd277229 ree60aa7  
    1717effective_radius(int mode, double radius_lg, double radius_sm, double penetration) 
    1818{ 
    19     if (mode == 1) { 
     19    switch (mode) { 
     20    case 1: // radius_large 
    2021        return radius_lg; 
    21     } else { 
     22    case 2: // radius_outer 
    2223        return radius_lg + 2.0*radius_sm - penetration; 
    2324    } 

  • sasmodels/models/raspberry.py

    rd277229 ree60aa7  
    152152 
    153153source = ["lib/sas_3j1x_x.c", "raspberry.c"] 
    154 effective_radius_type = ["radius_large","radius_outer"] 
     154effective_radius_type = ["radius_large", "radius_outer"] 
    155155 
    156156def random(): 

  • sasmodels/models/rectangular_prism.c

    rd277229 ree60aa7  
    88effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio) 
    99{ 
    10     if (mode == 1) { 
     10    switch (mode) { 
     11    case 1: // equivalent sphere 
    1112        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
    12     } else if (mode == 2) { 
     13    case 2: // half length_a 
    1314        return 0.5 * length_a; 
    14     } else if (mode == 3) { 
     15    case 3: // half length_b 
    1516        return 0.5 * length_a*b2a_ratio; 
    16     } else if (mode == 4) { 
     17    case 4: // half length_c 
    1718        return 0.5 * length_a*c2a_ratio; 
    18     } else if (mode == 5) { 
     19    case 5: // equivalent circular cross-section 
    1920        return length_a*sqrt(b2a_ratio/M_PI); 
    20     } else if (mode == 6) { 
     21    case 6: // half ab diagonal 
    2122        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
    22     } else { 
     23    case 7: // half diagonal 
    2324        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
    2425    } 

  • sasmodels/models/rectangular_prism.py

    rd277229 ree60aa7  
    136136source = ["lib/gauss76.c", "rectangular_prism.c"] 
    137137have_Fq = True 
    138 effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
    139                          "equivalent circular cross-section","half ab diagonal","half diagonal"] 
    140  
    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.) 
     138effective_radius_type = [ 
     139    "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     140    "equivalent circular cross-section", "half ab diagonal", "half diagonal", 
     141    ] 
    153142 
    154143def random(): 

  • sasmodels/models/sphere.c

    rd277229 ree60aa7  
    77effective_radius(int mode, double radius) 
    88{ 
     9    // case 1: radius 
    910    return radius; 
    1011} 

  • sasmodels/models/sphere.py

    r2cc8aa2 ree60aa7  
    7171effective_radius_type = ["radius"] 
    7272 
    73 #def ER(radius): 
    74 #    """ 
    75 #    Return equivalent radius (ER) 
    76 #    """ 
    77 #    return radius 
    78  
    79 # VR defaults to 1.0 
    80  
    8173def random(): 
    8274    radius = 10**np.random.uniform(1.3, 4) 

  • sasmodels/models/spherical_sld.c

    rd277229 ree60aa7  
    1 static double form_volume( 
    2     double fp_n_shells, 
    3     double thickness[], 
    4     double interface[]) 
    5 { 
    6     int n_shells= (int)(fp_n_shells + 0.5); 
    7     double r = 0.0; 
    8     for (int i=0; i < n_shells; i++) { 
    9         r += thickness[i] + interface[i]; 
    10     } 
    11     return M_4PI_3*cube(r); 
    12 } 
    13  
    141static double 
    15 effective_radius(int mode, double fp_n_shells, double thickness[], double interface[]) 
     2outer_radius(double fp_n_shells, double thickness[], double interface[]) 
    163{ 
    174    int n_shells= (int)(fp_n_shells + 0.5); 
    … …  
    218    } 
    229    return r; 
     10} 
     11 
     12static double form_volume( 
     13    double fp_n_shells, 
     14    double thickness[], 
     15    double interface[]) 
     16{ 
     17    return M_4PI_3*cube(outer_radius(fp_n_shells, thickness, interface)); 
     18} 
     19 
     20static double 
     21effective_radius(int mode, double fp_n_shells, double thickness[], double interface[]) 
     22{ 
     23    // case 1: outer radius 
     24    return outer_radius(fp_n_shells, thickness, interface); 
    2325} 
    2426 

  • sasmodels/models/spherical_sld.py

    re8eff7b ree60aa7  
    2222 
    2323    0: erf($\nu z$) 
    24      
     24 
    2525    1: Rpow($z^\nu$) 
    26      
     26 
    2727    2: Lpow($z^\nu$) 
    28      
     28 
    2929    3: Rexp($-\nu z$) 
    30      
     30 
    3131    4: Lexp($-\nu z$) 
    3232 
    … …  
    270270 
    271271 
    272 #def ER(n_shells, thickness, interface): 
    273 #    """Effective radius""" 
    274 #    n_shells = int(n_shells + 0.5) 
    275 #    total = (np.sum(thickness[:n_shells], axis=1) 
    276 #             + np.sum(interface[:n_shells], axis=1)) 
    277 #    return total 
    278  
    279  
    280272demo = { 
    281273    "n_shells": 5, 

  • sasmodels/models/triaxial_ellipsoid.c

    rd277229 ree60aa7  
    3030effective_radius(int mode, double radius_equat_minor, double radius_equat_major, double radius_polar) 
    3131{ 
    32     if (mode == 1) { 
     32    switch (mode) { 
     33    case 1: // equivalent sphere 
    3334        return radius_from_volume(radius_equat_minor,radius_equat_major, radius_polar); 
    34     } else if (mode == 2) { 
     35    case 2: // min radius 
    3536        return radius_from_min_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
    36     } else { 
     37    case 3: // max radius 
    3738        return radius_from_max_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
    3839    } 

  • sasmodels/models/triaxial_ellipsoid.py

    rd277229 ree60aa7  
    158158source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"] 
    159159have_Fq = True 
    160 effective_radius_type = ["equivalent sphere","min radius", "max radius"] 
    161  
    162 def ER(radius_equat_minor, radius_equat_major, radius_polar): 
    163     """ 
    164     Returns the effective radius used in the S*P calculation 
    165     """ 
    166     from .ellipsoid import ER as ellipsoid_ER 
    167  
    168     # now that radii can be in any size order, radii need sorting a,b,c 
    169     # where a~b and c is either much smaller or much larger 
    170     radii = np.vstack((radius_equat_major, radius_equat_minor, radius_polar)) 
    171     radii = np.sort(radii, axis=0) 
    172     selector = (radii[1] - radii[0]) > (radii[2] - radii[1]) 
    173     polar = np.where(selector, radii[0], radii[2]) 
    174     equatorial = np.sqrt(np.where(~selector, radii[0]*radii[1], radii[1]*radii[2])) 
    175     return ellipsoid_ER(polar, equatorial) 
     160effective_radius_type = ["equivalent sphere", "min radius", "max radius"] 
    176161 
    177162def random(): 

  • sasmodels/models/vesicle.c

    rd277229 ree60aa7  
     1// TODO: interface to form_volume/shell_volume not yet settled 
     2static double 
     3shell_volume(double *total, double radius, double thickness) 
     4{ 
     5    //note that for the vesicle model, the volume is ONLY the shell volume 
     6    *total = M_4PI_3 * cube(radius+thickness); 
     7    return *total - M_4PI_3*cube(radius); 
     8} 
     9 
    110static double 
    211form_volume(double radius, double thickness) 
    312{ 
    413    //note that for the vesicle model, the volume is ONLY the shell volume 
    5     return M_4PI_3*(cube(radius+thickness) - cube(radius)); 
     14    double total; 
     15    return shell_volume(&total, radius, thickness); 
    616} 
    717 
    … …  
    919effective_radius(int mode, double radius, double thickness) 
    1020{ 
     21    // case 1: outer radius 
    1122    return radius + thickness; 
    1223} 

  • sasmodels/models/vesicle.py

    r2cc8aa2 ree60aa7  
    102102effective_radius_type = ["outer radius"] 
    103103 
    104 #def ER(radius, thickness): 
    105 #    ''' 
    106 #    returns the effective radius used in the S*P calculation 
    107 # 
    108 #    :param radius: core radius 
    109 #    :param thickness: shell thickness 
    110 #    ''' 
    111 #    return radius + thickness 
    112  
    113 def VR(radius, thickness): 
    114     ''' 
    115     returns the volumes of the shell and of the whole sphere including the 
    116     core plus shell - is used to normalize when including polydispersity. 
    117  
    118     :param radius: core radius 
    119     :param thickness: shell thickness 
    120     :return whole: volume of core and shell 
    121     :return whole-core: volume of the shell 
    122     ''' 
    123  
    124     whole = 4./3. * pi * (radius + thickness)**3 
    125     core = 4./3. * pi * radius**3 
    126     return whole, whole - core 
    127  
    128104def random(): 
    129105    total_radius = 10**np.random.uniform(1.3, 5) 
Note: See TracChangeset for help on using the changeset viewer.