Changeset 513efc5 in sasmodels


Ignore:
Timestamp:
Jan 20, 2016 4:50:50 AM (9 years ago)
Author:
piotr
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
7ed702f
Parents:
30b4ddf
Message:

Code review issues from PK addressed.

Added Taylor expansion utility function.

Location:
sasmodels/models
Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/ellipsoid.c

    r3f832f9 r513efc5  
    77double _ellipsoid_kernel(double q, double rpolar, double requatorial, double sin_alpha) 
    88{ 
    9     double sn, cn; 
    109    double ratio = rpolar/requatorial; 
    1110    const double u = q*requatorial*sqrt(1.0 
    1211                   + sin_alpha*sin_alpha*(ratio*ratio - 1.0)); 
    13     SINCOS(u, sn, cn); 
    14     //const double f = ( u==0.0 ? 1.0 : 3.0*(sn-u*cn)/(u*u*u) ); 
    15     const double usq = u*u; 
    16     const double f = (u < 1.e-1) 
    17         ? 1.0 + usq*(-3./30. + usq*(3./840. + usq*(-3./45360.)))// + qrsq*(3./3991680.)))) 
    18         : 3.0*(sn/u - cn)/usq; 
     12    const double f = J1c(u); 
     13 
    1914    return f*f; 
    2015} 

  • sasmodels/models/ellipsoid.py

    rcade620 r513efc5  
    152152             ] 
    153153 
    154 source = ["lib/J1.c", "lib/gauss76.c", "ellipsoid.c"] 
     154source = ["lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "ellipsoid.c"] 
    155155 
    156156def ER(rpolar, requatorial): 

  • sasmodels/models/gel_fit.py

    r30b4ddf r513efc5  
    1818.. math:: 
    1919 
    20     I(Q) = I(0)_L \frac{1}{\left( 1+\left[ ((D+1/3)Q^2a_{1}^2 \right]\right)^{D/2}} + I(0)_G exp\left( -Q^2a_{2}^2\right) + B 
     20    I(Q) = I(0)_L \frac{1}{\left( 1+\left[ ((D+1/3)Q^2a_{1}^2 
     21    \right]\right)^{D/2}} + I(0)_G exp\left( -Q^2a_{2}^2\right) + B 
    2122 
    2223where 
    … …  
    2627    a_{2}^2 \approx \frac{R_{g}^2}{3} 
    2728 
    28 Note that the first term reduces to the Ornstein-Zernicke equation when $D = 2$; 
    29 ie, when the Flory exponent is 0.5 (theta conditions). 
    30 In gels with significant hydrogen bonding $D$ has been reported to be ~2.6 to 2.8. 
     29Note that the first term reduces to the Ornstein-Zernicke equation 
     30when $D = 2$; ie, when the Flory exponent is 0.5 (theta conditions). 
     31In gels with significant hydrogen bonding $D$ has been reported to be 
     32~2.6 to 2.8. 
    3133 
    3234 
    … …  
    3840--------- 
    3941 
    40 Mitsuhiro Shibayama, Toyoichi Tanaka, Charles C Han, J. Chem. Phys. 1992, 97 (9), 6829-6841 
     42Mitsuhiro Shibayama, Toyoichi Tanaka, Charles C Han, J. Chem. Phys. 1992, 97 (9), 
     436829-6841 
    4144 
    42 Simon Mallam, Ferenc Horkay, Anne-Marie Hecht, Adrian R Rennie, Erik Geissler, Macromolecules 1991, 24, 543-548 
     45Simon Mallam, Ferenc Horkay, Anne-Marie Hecht, Adrian R Rennie, Erik Geissler, 
     46Macromolecules 1991, 24, 543-548 
    4347 
    4448""" 

  • sasmodels/models/polymer_excl_volume.py

    r23d3b41 r513efc5  
    11r""" 
    2 This model describes the scattering from polymer chains subject to excluded volume effects 
    3 and has been used as a template for describing mass fractals. 
     2This model describes the scattering from polymer chains subject to excluded 
     3volume effects and has been used as a template for describing mass fractals. 
    44 
    55Definition 
    66---------- 
    77 
    8 The form factor was originally presented in the following integral form (Benoit, 1957) 
     8The form factor was originally presented in the following integral form 
     9(Benoit, 1957) 
    910 
    1011.. math:: 
    … …  
    1617$a$ is the statistical segment length of the polymer chain, 
    1718and $n$ is the degree of polymerization. 
    18 This integral was later put into an almost analytical form as follows (Hammouda, 1993) 
     19This integral was later put into an almost analytical form as follows 
     20(Hammouda, 1993) 
    1921 
    2022.. math:: 
    … …  
    4345Note that this model applies only in the mass fractal range (ie, $5/3<=m<=3$ ) 
    4446and **does not apply** to surface fractals ( $3<m<=4$ ). 
    45 It also does not reproduce the rigid rod limit (m=1) because it assumes chain flexibility 
    46 from the outset. It may cover a portion of the semi-flexible chain range ( $1<m<5/3$ ). 
     47It also does not reproduce the rigid rod limit (m=1) because it assumes chain 
     48flexibility from the outset. It may cover a portion of the semi-flexible chain 
     49range ( $1<m<5/3$ ). 
    4750 
    48 A low-Q expansion yields the Guinier form and a high-Q expansion yields the Porod form 
    49 which is given by 
     51A low-Q expansion yields the Guinier form and a high-Q expansion yields the 
     52Porod form which is given by 
    5053 
    5154.. math:: 
    5255 
    53     P(Q\rightarrow \infty) = \frac{1}{\nu U^{1/2\nu}}\Gamma\left(\frac{1}{2\nu}\right) - 
    54     \frac{1}{\nu U^{1/\nu}}\Gamma\left(\frac{1}{\nu}\right) 
     56    P(Q\rightarrow \infty) = \frac{1}{\nu U^{1/2\nu}}\Gamma\left( 
     57    \frac{1}{2\nu}\right) - \frac{1}{\nu U^{1/\nu}}\Gamma\left( 
     58    \frac{1}{\nu}\right) 
    5559 
    5660Here $\Gamma(x) = \gamma(x,\infty)$ is the gamma function. 
    … …  
    102106name = "polymer_excl_volume" 
    103107title = "Polymer Excluded Volume model" 
    104 description = """Compute the scattering intensity from polymers with excluded volume effects. 
     108description = """Compute the scattering intensity from polymers with excluded 
     109                volume effects. 
    105110                rg:         radius of gyration 
    106111                porod_exp:  Porod exponent 
    … …  
    127132 
    128133    intensity = ((1.0/(nu*power(u, o2nu))) * (gamma(o2nu)*gammainc(o2nu, u) - 
    129                  1.0/power(u, o2nu) * gamma(porod_exp)*gammainc(porod_exp, u)))*(q > 0) + \ 
    130                  1.0*(q <= 0) 
     134                  1.0/power(u, o2nu) * gamma(porod_exp) * 
     135                  gammainc(porod_exp, u))) * (q > 0) + 1.0*(q <= 0) 
    131136 
    132137    return intensity 
    … …  
    156161         [{'rg': 2.2, 'porod_exp': 22.0, 'background': 100.0}, 5.0, 100.0], 
    157162 
    158          [{'rg': 1.1, 'porod_exp': 1, 'background': 10.0, 'scale': 1.25}, 20000., 10.0000712097] 
     163         [{'rg': 1.1, 'porod_exp': 1, 'background': 10.0, 'scale': 1.25}, 
     164         20000., 10.0000712097] 
    159165         ] 

  • sasmodels/models/sphere.py

    rd53d3cd r513efc5  
    7878             ] 
    7979 
     80source = ["lib/J1c.c"] 
    8081 
    8182# No volume normalization despite having a volume parameter 
    … …  
    8788Iq = """ 
    8889    const double qr = q*radius; 
    89     const double qrsq = qr*qr; 
    90     double sn, cn; 
    91     SINCOS(qr, sn, cn); 
    92     // Use taylor series for low q to avoid cancellation error.  Tested against 
    93     // the following expression in quad precision: 
    94     //     3.0*(sn-qr*cn)/(qr*qr*qr); 
    95     // Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but 
    96     // in this case it is likely cancellation errors in the original expression 
    97     // using double precision that are the source.  Single precision only 
    98     // requires the first 3 terms.  Double precision requires the 4th term. 
    99     // The fifth term is not needed, and is commented out below. 
    100     const double bes = (qr < 1.e-1) 
    101         ? 1.0 + qrsq*(-3./30. + qrsq*(3./840. + qrsq*(-3./45360.)))// + qrsq*(3./3991680.)))) 
    102         : 3.0*(sn/qr - cn)/qrsq; 
     90    const double bes = J1c(qr); 
    10391    const double fq = bes * (sld - solvent_sld) * form_volume(radius); 
    10492    return 1.0e-4*fq*fq; 

  • sasmodels/models/surface_fractal.c

    rda84551 r513efc5  
    11double form_volume(double radius); 
    22 
    3 double Iq(double q, double radius, double surface_dim, double cutoff_length); 
    4 double Iqxy(double qx, double qy, double radius, double surface_dim, double cutoff_length); 
     3double Iq(double q, 
     4          double radius, 
     5          double surface_dim, 
     6          double cutoff_length); 
     7 
     8double Iqxy(double qx, double qy, 
     9            double radius, 
     10            double surface_dim, 
     11            double cutoff_length); 
    512 
    613static double _gamln(double q) 
    … …  
    2633} 
    2734 
    28 static double surface_fractal_kernel(double q, 
     35static double _surface_fractal_kernel(double q, 
    2936    double radius, 
    3037    double surface_dim, 
    … …  
    3340    double pq, sq, mmo, result; 
    3441 
    35     //calculate P(q) for the spherical subunits; not normalized 
    36         pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2); 
     42    //Replaced the original formula with Taylor expansion near zero. 
     43    //pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2); 
     44 
     45    pq = J1c(q*radius); 
     46    pq = pq*pq; 
    3747 
    3848    //calculate S(q) 
    … …  
    5969    ) 
    6070{ 
    61     return surface_fractal_kernel(q, radius, surface_dim, cutoff_length); 
     71    return _surface_fractal_kernel(q, radius, surface_dim, cutoff_length); 
    6272} 
    6373 
    … …  
    6979{ 
    7080    double q = sqrt(qx*qx + qy*qy); 
    71     return surface_fractal_kernel(q, radius, surface_dim, cutoff_length); 
     81    return _surface_fractal_kernel(q, radius, surface_dim, cutoff_length); 
    7282} 
    7383 

  • sasmodels/models/surface_fractal.py

    r30b4ddf r513efc5  
    2222.. math:: 
    2323 
    24     S(q) = \frac{\Gamma(5-D_S)\zeta^{5-D_S}}{\left[1+(q\zeta)^2\right]^{(5-D_S)/2}} 
     24    S(q) = \frac{\Gamma(5-D_S)\zeta^{5-D_S}}{\left[1+(q\zeta)^2 
     25    \right]^{(5-D_S)/2}} 
    2526    \frac{sin\left[(D_S - 5) tan^{-1}(q\zeta) \right]}{q} 
    2627 
    … …  
    3334    V = \frac{4}{3}\pi R^3 
    3435 
    35 where $R$ is the radius of the building block, $D_S$ is the **surface** fractal dimension, 
    36 $\zeta$ is the cut-off length, $\rho_{solvent}$ is the scattering length density of the solvent, 
     36where $R$ is the radius of the building block, $D_S$ is the **surface** fractal 
     37dimension,$\zeta$ is the cut-off length, $\rho_{solvent}$ is the scattering 
     38length density of the solvent, 
    3739and $\rho_{particle}$ is the scattering length density of particles. 
    3840 
    3941.. note:: 
    4042    The surface fractal dimension $D_s$ is only valid if $1<surface\_dim<3$. 
    41     It is also only valid over a limited $q$ range (see the reference for details) 
     43    It is also only valid over a limited $q$ range (see the reference for 
     44    details) 
    4245 
    4346 
    … …  
    7881 
    7982#             ["name", "units", default, [lower, upper], "type","description"], 
    80 parameters = [["radius",        "Ang", 10.0, [0, inf],   "", "Particle radius"], 
    81               ["surface_dim",   "",    2.0,  [0, inf],   "", "Surface fractal dimension"], 
    82               ["cutoff_length", "Ang", 500., [0.0, inf], "",  "Cut-off Length"], 
     83parameters = [["radius",        "Ang", 10.0, [0, inf],   "", 
     84               "Particle radius"], 
     85              ["surface_dim",   "",    2.0,  [0, inf],   "", 
     86               "Surface fractal dimension"], 
     87              ["cutoff_length", "Ang", 500., [0.0, inf], "", 
     88               "Cut-off Length"], 
    8389              ] 
    8490 
    8591 
    86 source = ["surface_fractal.c"] 
     92source = ["lib/J1c.c", "surface_fractal.c"] 
    8793 
    8894demo = dict(scale=1, background=0, 

  • sasmodels/models/triaxial_ellipsoid.c

    rd53d3cd r513efc5  
    3636            const double y = 0.5*(Gauss76Z[j] + 1.0); 
    3737            const double t = q*sqrt(acosx2 + bsinx2*(1.0-y*y) + c2*y*y); 
    38             SINCOS(t, st, ct); 
    39             //const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) ); 
    40             const double tsq = t*t; 
    41             const double fq = (t < 1.e-1) 
    42                 ? 1.0 + tsq*(-3./30. + tsq*(3./840. + tsq*(-3./45360.)))// + tsq*(3./3991680.)))) 
    43                 : 3.0*(st/t - ct)/tsq; 
     38            const double fq = J1c(t); 
    4439            inner += Gauss76Wt[j] * fq * fq ; 
    4540        } 

  • sasmodels/models/triaxial_ellipsoid.py

    reb69cce r513efc5  
    117117             ] 
    118118 
    119 source = ["lib/J1.c", "lib/gauss76.c", "triaxial_ellipsoid.c"] 
     119source = ["lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "triaxial_ellipsoid.c"] 
    120120 
    121121def ER(req_minor, req_major, rpolar): 

  • sasmodels/models/two_lorentzian.py

    r421e55c r513efc5  
    11r""" 
    2 This model calculates an empirical functional form for SAS data characterized by two Lorentzian-type functions. 
     2This model calculates an empirical functional form for SAS data characterized 
     3by two Lorentzian-type functions. 
    34 
    45Definition 
    … …  
    1213 
    1314where $A$ = Lorentzian scale factor #1, $C$ = Lorentzian scale #2, 
    14 $\xi_1$ and $\xi_2$ are the corresponding correlation lengths, and $n$ and $m$ are the respective 
    15 power law exponents (set $n = m = 2$ for Ornstein-Zernicke behaviour). 
     15$\xi_1$ and $\xi_2$ are the corresponding correlation lengths, and $n$ and 
     16$m$ are the respective power law exponents (set $n = m = 2$ for 
     17Ornstein-Zernicke behaviour). 
    1618 
    1719For 2D data the scattering intensity is calculated in the same way as 1D, 
    … …  
    5254category = "shape-independent" 
    5355 
    54 #             ["name", "units", default, [lower, upper], "type", "description"], 
    55 parameters = [["lorentz_scale_1",  "",     10.0, [-inf, inf], "", "First power law scale factor"], 
    56               ["lorentz_length_1", "Ang", 100.0, [-inf, inf], "", "First Lorentzian screening length"], 
    57               ["lorentz_exp_1",    "",      3.0, [-inf, inf], "", "First exponent of power law"], 
    58               ["lorentz_scale_2",  "",      1.0, [-inf, inf], "", "Second scale factor for broad Lorentzian peak"], 
    59               ["lorentz_length_2", "Ang",  10.0, [-inf, inf], "", "Second Lorentzian screening length"], 
    60               ["lorentz_exp_2",    "",      2.0, [-inf, inf], "", "Second exponent of power law"], 
     56#            ["name", "units", default, [lower, upper], "type", "description"], 
     57parameters = [["lorentz_scale_1",  "",     10.0, [-inf, inf], "", 
     58               "First power law scale factor"], 
     59              ["lorentz_length_1", "Ang", 100.0, [-inf, inf], "", 
     60               "First Lorentzian screening length"], 
     61              ["lorentz_exp_1",    "",      3.0, [-inf, inf], "", 
     62               "First exponent of power law"], 
     63              ["lorentz_scale_2",  "",      1.0, [-inf, inf], "", 
     64               "Second scale factor for broad Lorentzian peak"], 
     65              ["lorentz_length_2", "Ang",  10.0, [-inf, inf], "", 
     66               "Second Lorentzian screening length"], 
     67              ["lorentz_exp_2",    "",      2.0, [-inf, inf], "", 
     68               "Second exponent of power law"], 
    6169              ] 
    6270 
    6371 
    64 def Iq(q, lorentz_scale_1, lorentz_length_1, lorentz_exp_1, lorentz_scale_2, lorentz_length_2, lorentz_exp_2): 
     72def Iq(q, 
     73       lorentz_scale_1, 
     74       lorentz_length_1, 
     75       lorentz_exp_1, 
     76       lorentz_scale_2, 
     77       lorentz_length_2, 
     78       lorentz_exp_2): 
    6579 
    6680    """ 
    … …  
    7589    """ 
    7690 
    77     intensity  = lorentz_scale_1/(1.0 + power(q*lorentz_length_1, lorentz_exp_1)) 
    78     intensity += lorentz_scale_2/(1.0 + power(q*lorentz_length_2, lorentz_exp_2)) 
     91    intensity  = lorentz_scale_1/(1.0 + 
     92                                  power(q*lorentz_length_1, lorentz_exp_1)) 
     93    intensity += lorentz_scale_2/(1.0 + 
     94                                  power(q*lorentz_length_2, lorentz_exp_2)) 
    7995 
    8096    return intensity 
    … …  
    101117               lorentz_exp_1='exponent_1',  lorentz_exp_2='exponent_2') 
    102118 
    103 tests = [[{'lorentz_scale_1': 10, 'lorentz_length_1': 100.0, 'lorentz_exp_1' : 3.0, 
    104           'lorentz_scale_2': 1, 'lorentz_length_2': 10.0, 'lorentz_exp_2' : 2.0, 
     119tests = [[{'lorentz_scale_1':   10.0, 
     120           'lorentz_length_1': 100.0, 
     121           'lorentz_exp_1':      3.0, 
     122           'lorentz_scale_2':    1.0, 
     123           'lorentz_length_2':  10.0, 
     124           'lorentz_exp_2':      2.0, 
    105125           }, 0.000332070182643, 10.9996228107], 
    106126 
    107          [{'lorentz_scale_1': 0, 'lorentz_length_1': 0.0, 'lorentz_exp_1' : 0.0, 
    108           'lorentz_scale_2': 0, 'lorentz_length_2': 0.0, 'lorentz_exp_2' : 0.0, 
    109            'background':100. 
     127         [{'lorentz_scale_1':  0.0, 
     128           'lorentz_length_1': 0.0, 
     129           'lorentz_exp_1':    0.0, 
     130           'lorentz_scale_2':  0.0, 
     131           'lorentz_length_2': 0.0, 
     132           'lorentz_exp_2':    0.0, 
     133           'background':     100.0 
    110134           }, 5.0, 100.0], 
    111135 
    112          [{'lorentz_scale_1': 200, 'lorentz_length_1': 10.0, 'lorentz_exp_1' : 0.1, 
    113           'lorentz_scale_2': 0.1, 'lorentz_length_2': 5.0, 'lorentz_exp_2' : 2.0 
     136         [{'lorentz_scale_1': 200.0, 
     137           'lorentz_length_1': 10.0, 
     138           'lorentz_exp_1':     0.1, 
     139           'lorentz_scale_2':   0.1, 
     140           'lorentz_length_2':  5.0, 
     141           'lorentz_exp_2':     2.0 
    114142           }, 20000., 45.5659201896], 
    115143         ] 
Note: See TracChangeset for help on using the changeset viewer.