Changes in / [44bd2be:a8bd500] in sasmodels


Ignore:
Files:
8 added
2 deleted
14 edited

Legend:

Unmodified
Added
Removed
  • doc/_extensions/dollarmath.py

    r19dcb933 r103ea45  
    1212import re 
    1313 
    14 _dollar = re.compile(r"(?:^|(?<=\s))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;\\]))") 
     14_dollar = re.compile(r"(?:^|(?<=\s|[(]))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;)\\]))") 
    1515_notdollar = re.compile(r"\\[$]") 
    1616 
     
    5151    assert replace_dollar(u"dollar\$ escape")==u"dollar$ escape" 
    5252    assert replace_dollar(u"dollar \$escape\$ too")==u"dollar $escape$ too" 
     53    assert replace_dollar(u"spaces $in the$ math")==u"spaces :math:`in the` math" 
    5354    assert replace_dollar(u"emb\ $ed$\ ed")==u"emb\ :math:`ed`\ ed" 
    5455    assert replace_dollar(u"$first$a")==u"$first$a" 
    5556    assert replace_dollar(u"a$last$")==u"a$last$" 
     57    assert replace_dollar(u"$37")==u"$37" 
     58    assert replace_dollar(u"($37)")==u"($37)" 
     59    assert replace_dollar(u"$37 - $43")==u"$37 - $43" 
     60    assert replace_dollar(u"($37, $38)")==u"($37, $38)" 
    5661    assert replace_dollar(u"a $mid$dle a")==u"a $mid$dle a" 
     62    assert replace_dollar(u"a ($in parens$) a")==u"a (:math:`in parens`) a" 
     63    assert replace_dollar(u"a (again $in parens$) a")==u"a (again :math:`in parens`) a" 
    5764 
    5865if __name__ == "__main__": 
  • sasmodels/generate.py

    r0a4628d rfa8011eb  
    191191The function :func:`make` loads the metadata from the module and returns 
    192192the kernel source.  The function :func:`doc` extracts the doc string 
    193 and adds the parameter table to the top.  The function :func:`sources` 
     193and adds the parameter table to the top.  The function :func:`model_sources` 
    194194returns a list of files required by the model. 
    195195""" 
     
    206206import numpy as np 
    207207 
    208 __all__ = ["make", "doc", "sources", "convert_type"] 
     208#__all__ = ["make", "doc", "model_sources", "convert_type"] 
    209209 
    210210C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 
  • sasmodels/kernel_template.c

    r840b859 r960cd80  
    1212#ifndef USE_OPENCL 
    1313#  ifdef __cplusplus 
    14      #include <cstdio> 
    15      #include <cmath> 
    16      using namespace std; 
    17      #if defined(_MSC_VER) 
     14      #include <cstdio> 
     15      #include <cmath> 
     16      using namespace std; 
     17      #if defined(_MSC_VER) 
    1818         #include <limits> 
    1919         #include <float.h> 
    2020         #define kernel extern "C" __declspec( dllexport ) 
    2121         inline double trunc(double x) { return x>=0?floor(x):-floor(-x); } 
    22              inline double fmin(double x, double y) { return x>y ? y : x; } 
    23              inline double fmax(double x, double y) { return x<y ? y : x; } 
    24              inline double isnan(double x) { return _isnan(x); } 
    25              #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN 
    26              static double cephes_expm1(double x) { 
    27                  // Adapted from the cephes math library. 
    28                  // Copyright 1984 - 1992 by Stephen L. Moshier 
    29                  if (x != x || x == 0.0) { 
    30                      return x; // NaN and +/- 0 
    31                  } else if (x < -0.5 || x > 0.5) { 
    32                      return exp(x) - 1.0; 
    33                  } else { 
    34                      const double xsq = x*x; 
    35                      const double p = ((( 
    36                           +1.2617719307481059087798E-4)*xsq 
    37                       +3.0299440770744196129956E-2)*xsq 
    38                       +9.9999999999999999991025E-1); 
    39                  const double q = (((( 
    40                       +3.0019850513866445504159E-6)*xsq 
    41                       +2.5244834034968410419224E-3)*xsq 
    42                       +2.2726554820815502876593E-1)*xsq 
    43                       +2.0000000000000000000897E0); 
    44                  double r = x * p; 
    45                      r =  r / (q - r); 
    46                      return r+r; 
    47                  } 
    48              } 
    49              #define expm1 cephes_expm1 
     22         inline double fmin(double x, double y) { return x>y ? y : x; } 
     23         inline double fmax(double x, double y) { return x<y ? y : x; } 
     24         inline double isnan(double x) { return _isnan(x); } 
     25         #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN 
     26         static double cephes_expm1(double x) { 
     27            // Adapted from the cephes math library. 
     28            // Copyright 1984 - 1992 by Stephen L. Moshier 
     29            if (x != x || x == 0.0) { 
     30               return x; // NaN and +/- 0 
     31            } else if (x < -0.5 || x > 0.5) { 
     32               return exp(x) - 1.0; 
     33            } else { 
     34               const double xsq = x*x; 
     35               const double p = ((( 
     36                  +1.2617719307481059087798E-4)*xsq 
     37                  +3.0299440770744196129956E-2)*xsq 
     38                  +9.9999999999999999991025E-1); 
     39               const double q = (((( 
     40                  +3.0019850513866445504159E-6)*xsq 
     41                  +2.5244834034968410419224E-3)*xsq 
     42                  +2.2726554820815502876593E-1)*xsq 
     43                  +2.0000000000000000000897E0); 
     44               double r = x * p; 
     45               r =  r / (q - r); 
     46               return r+r; 
     47             } 
     48         } 
     49         #define expm1 cephes_expm1 
    5050     #else 
    5151         #define kernel extern "C" 
     
    260260        const double vol_weight = VOLUME_WEIGHT_PRODUCT; 
    261261        vol += vol_weight*form_volume(VOLUME_PARAMETERS); 
    262       #endif 
    263262        norm_vol += vol_weight; 
     263      #endif 
    264264      } 
    265265      //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); } 
  • sasmodels/models/capped_cylinder.c

    rc138211 r139c528  
    2222    const double upper = 1.0; 
    2323    const double lower = h/cap_radius; // integral lower bound 
     24    const double zm = 0.5*(upper-lower); 
     25    const double zb = 0.5*(upper+lower); 
    2426    // cos term in integral is: 
    2527    //    cos (q (R t - h + L/2) cos(alpha)) 
     
    3638        // translate a point in [-1,1] to a point in [lower,upper] 
    3739        //const double t = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    38         const double t = 0.5*(Gauss76Z[i]*(upper-lower)+upper+lower); 
     40        const double t = Gauss76Z[i]*zm + zb; 
    3941        const double radical = 1.0 - t*t; 
    4042        const double arg = qrst*sqrt(radical); // cap bessel function arg 
  • sasmodels/models/core_shell_bicelle.py

    r8007311 rfa8011eb  
    77The form factor is normalized by the particle volume. 
    88 
    9 .. _core-shell-cylinder-geometry: 
     9.. _core-shell-bicelle-geometry: 
    1010 
    1111.. figure:: img/core_shell_bicelle_geometry.png 
  • sasmodels/models/core_shell_sphere.py

    r8c9dbc9 rfa8011eb  
    1515 
    1616.. math:: 
     17 
    1718    F^2(q)=\frac{3}{V_s}\left[V_c(\rho_c-\rho_s)\frac{\sin(qr_c)-qr_c\cos(qr_c)}{(qr_c)^3}+ 
    1819    V_s(\rho_s-\rho_{solv})\frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3}\right] 
     
    4142our model and the output of the NIST software. 
    4243 
    43 .. image:: img/core_shell_sphere_1d.jpg 
     44.. figure:: img/core_shell_sphere_1d.jpg 
    4445 
    45     Figure 1: Comparison of the SasView scattering intensity for a core-shell sphere with 
     46    Comparison of the SasView scattering intensity for a core-shell sphere with 
    4647    the output of the NIST SANS analysis software. The parameters were set to: 
    4748    *scale* = 1.0, *radius* = 60 , *contrast* = 1e-6 |Ang^-2|, and 
  • sasmodels/models/ellipsoid.py

    r9c461c7 r431caae  
    111111---------- 
    112112 
    113 L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum, 
    114 New York, 1987. 
     113L A Feigin and D I Svergun. 
     114*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 
     115Plenum Press, New York, 1987. 
    115116""" 
    116117 
  • sasmodels/models/elliptical_cylinder.py

    rb7c2fce rfa8011eb  
    1010to any of the orientation angles, and also for the minor radius and the ratio of the ellipse radii. 
    1111 
    12 .. image:: img/elliptical_cylinder_geometry.gif 
     12.. figure:: img/elliptical_cylinder_geometry.gif 
    1313 
    14     *Figure.* *a* = *r_minor* and |nu|\ :sub:`n` = $r_ratio$ (i.e., $r_major / r_minor$). 
     14    *a* = *r_minor* and |nu|\ :sub:`n` = $r_ratio$ (i.e., $r_major / r_minor$). 
    1515 
    1616The function calculated is 
    1717 
    1818.. math:: 
     19 
    1920    I(\mathbf{q})=\frac{1}{V_{cyl}}\int{d\psi}\int{d\phi}\int{p(\theta,\phi,\psi)F^2(\mathbf{q},\alpha,\psi)\sin(\theta)d\theta} 
    2021 
     
    2223 
    2324.. math:: 
     25 
    2426    F(\mathbf{q},\alpha,\psi)=2\frac{J_1(a)\sin(b)}{ab} 
    2527    \\ 
     
    4042    P(q) = scale  <F^2> / V 
    4143 
    42 The returned value is scaled to units of |cm^-1|. 
    43  
    4444To provide easy access to the orientation of the elliptical cylinder, we define the axis of the cylinder using two 
    4545angles |theta|, |phi| and |bigpsi|. As for the case of the cylinder, the angles |theta| and |phi| are defined on 
     
    4949All angle parameters are valid and given only for 2D calculation; ie, an oriented system. 
    5050 
    51 .. image:: img/elliptical_cylinder_geometry_2d.jpg 
     51.. figure:: img/elliptical_cylinder_geometry_2d.jpg 
    5252 
    53     *Figure. Definition of angles for 2D* 
     53    Definition of angles for 2D 
    5454 
    55 .. image:: img/core_shell_bicelle_fig2.jpg 
     55.. figure:: img/core_shell_bicelle_fig2.jpg 
    5656 
    57     *Figure. Examples of the angles for oriented elliptical cylinders against the detector plane.* 
     57    Examples of the angles for oriented elliptical cylinders against the detector plane. 
    5858 
    5959NB: The 2nd virial coefficient of the cylinder is calculated based on the averaged radius (= sqrt(*r_minor*\ :sup:`2` \* *r_ratio*)) 
     
    6161 
    6262 
    63 .. image:: img/elliptical_cylinder_comparison_1d.jpg 
     63.. figure:: img/elliptical_cylinder_comparison_1d.jpg 
    6464 
    65     *Figure. 1D plot using the default values (w/1000 data point).* 
     65    1D plot using the default values (w/1000 data point). 
    6666 
    6767Validation 
     
    7373and 76 degrees are taken for the angles of |theta|, |phi|, and |bigpsi| respectively). 
    7474 
    75 .. image:: img/elliptical_cylinder_validation_1d.gif 
     75.. figure:: img/elliptical_cylinder_validation_1d.gif 
    7676 
    77     *Figure. Comparison between 1D and averaged 2D.* 
     77    Comparison between 1D and averaged 2D. 
    7878 
    7979In the 2D average, more binning in the angle |phi| is necessary to get the proper result. The following figure shows 
    8080the results of the averaging by varying the number of angular bins. 
    8181 
    82 .. image:: img/elliptical_cylinder_averaging.gif 
     82.. figure:: img/elliptical_cylinder_averaging.gif 
    8383 
    84     *Figure. The intensities averaged from 2D over different numbers of bins and angles.* 
     84    The intensities averaged from 2D over different numbers of bins and angles. 
    8585 
    8686Reference 
  • sasmodels/models/guinier_porod.py

    r21d1031 rfa8011eb  
    5151    q = \sqrt{q_x^2+q_y^2} 
    5252 
    53 .. image:: img/guinier_porod_model.jpg 
     53.. figure:: img/guinier_porod_model.jpg 
    5454 
    55     Figure 1: Guinier-Porod model for $R_g=100$ |Ang|, $s=1$, $m=3$, and $background=0.1$. 
     55    Guinier-Porod model for $R_g=100$ |Ang|, $s=1$, $m=3$, and $background=0.1$. 
    5656 
    5757 
  • sasmodels/models/hardsphere.py

    r093f754 r97e6d3c  
    6969               return(HARDSPH); 
    7070      } 
    71       D=pow((1.-volfraction),2); 
    72       A=pow((1.+2.*volfraction)/D, 2); 
     71      // removing use of pow(xxx,2) and rearranging the calcs of A, B & G cut ~40% off execution time ( 0.5 to 0.3 msec) 
     72      X = 1.0/( 1.0 -volfraction); 
     73      D= X*X; 
     74      A= (1.+2.*volfraction)*D; 
     75      A *=A; 
    7376      X=fabs(q*effect_radius*2.0); 
    7477 
     
    7780                 return(HARDSPH); 
    7881      } 
    79       X2=pow(X,2); 
    80       X4=pow(X2,2); 
    81       B=-6.*volfraction* pow((1.+0.5*volfraction)/D ,2); 
     82      X2 =X*X; 
     83      B = (1.0 +0.5*volfraction)*D; 
     84      B *= B; 
     85      B *= -6.*volfraction; 
    8286      G=0.5*volfraction*A; 
    8387 
    8488      if(X < 0.2) { 
    85       // use Taylor series expansion for small X, IT IS VERY PICKY ABOUT THE X CUT OFF VALUE, ought to be lower in double.  
    86       // No obvious way to rearrange the equations to avoid needing a very high number of significant figures.  
     89      // RKH Feb 2016, use Taylor series expansion for small X, IT IS VERY PICKY ABOUT THE X CUT OFF VALUE, ought to be lower in double.  
     90      // else no obvious way to rearrange the equations to avoid needing a very high number of significant figures.  
    8791      // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient  
    88       // for 5 or 6 significant figures but I put the X^4 one in anyway  
    89             FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4; 
     92      // for 5 or 6 significant figures, but I put the X^4 one in anyway  
     93            //FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4; 
     94            // refactoring the polynomial makes it very slightly faster (0.5 not 0.6 msec) 
     95            //FF = 8*A +6*B + 4*G + ( -0.8*A -2.0*B/3.0 -0.5*G +(A/35. +B/40. +G/50.)*X2)*X2; 
     96 
     97            FF = 8.0*A +6.0*B + 4.0*G + ( -0.8*A -B/1.5 -0.5*G +(A/35. +0.0125*B +0.02*G)*X2)*X2; 
     98 
    9099            // combining the terms makes things worse at smallest Q in single precision 
    91100            //FF = (8-0.8*X2)*A +(3.0-X2/3.)*2*B + (4+0.5*X2)*G +(A/35. +B/40. +G/50.)*X4; 
    92101            // note that G = -volfraction*A/2, combining this makes no further difference at smallest Q 
    93             //FF = (8 +2.*volfraction + ( volfraction/4. -0.8 +(volfraction/100. -1./35.)*X2 )*X2 )*A  + (3.0 -X2/3. +X4/40)*2*B; 
     102            //FF = (8 +2.*volfraction + ( volfraction/4. -0.8 +(volfraction/100. -1./35.)*X2 )*X2 )*A  + (3.0 -X2/3. +X4/40.)*2.*B; 
    94103            HARDSPH= 1./(1. + volfraction*FF ); 
    95104            return(HARDSPH); 
    96105      } 
     106      X4=X2*X2; 
    97107      SINCOS(X,S,C); 
    98108 
    99 // RKH Feb 2016, use version from FISH code as it is better than original sasview one at small Q in single precision 
    100       FF=A*(S-X*C)/X + B*(2.*X*S -(X2-2.)*C -2.)/X2 + G*( (4.*X2*X -24.*X)*S -(X4 -12.*X2 +24.)*C +24. )/X4; 
     109// RKH Feb 2016, use version FISH code as is better than original sasview one at small Q in single precision, and more than twice as fast in double. 
     110      //FF=A*(S-X*C)/X + B*(2.*X*S -(X2-2.)*C -2.)/X2 + G*( (4.*X2*X -24.*X)*S -(X4 -12.*X2 +24.)*C +24. )/X4; 
     111      // refactoring the polynomial here & above makes it slightly faster 
     112 
     113      FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S -(X2-2.)*C -2.) )/X + A*(S-X*C))/X ; 
    101114      HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); 
    102115 
    103 // rearrange the terms, is now about same as sasmodels 
     116      // changing /X and /X2 to *MX1 and *MX2, no significantg difference? 
     117      //MX=1.0/X; 
     118      //MX2=MX*MX; 
     119      //FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )*MX2 + B*(2.*X*S -(X2-2.)*C -2.) )*MX + A*(S-X*C)) ; 
     120      //HARDSPH= 1./(1. + 24.*volfraction*FF*MX2*MX ); 
     121 
     122// grouping the terms, was about same as sasmodels for single precision issues 
    104123//     FF=A*(S/X-C) + B*(2.*S/X - C +2.0*(C-1.0)/X2) + G*( (4./X -24./X3)*S -(1.0 -12./X2 +24./X4)*C +24./X4 ); 
    105124//     HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); 
  • sasmodels/models/line.py

    re66075f rfa8011eb  
    1313.. note:: 
    1414    For 2D plots intensity has different definition than other shape independent models 
     15 
    1516.. math:: 
    1617    I(q) = I(qx) \cdot I(qy) 
    17  
    18 .. figure:: None 
    1918 
    2019References 
  • sasmodels/models/rpa.py

    r8dd6914 rfa8011eb  
    4343component. 
    4444 
    45 .. figure:: img/image215.jpg 
     45.. figure:: img/rpa_1d.jpg 
    4646 
    4747    1D plot using the default values (w/500 data points). 
  • sasmodels/models/spherepy.py

    reb69cce rd2950f4  
    118118    g[low] = sqrt(1 - dlow2 / 4.) * (1 + dlow2 / 8.) + dlow2 / 2.*(1 - dlow2 / 16.) * log(dlow / (2. + sqrt(4. - dlow2))) 
    119119    return g 
    120 sesans.vectorized = True  # sesans accepts and array of z values 
     120sesans.vectorized = True  # sesans accepts an array of z values 
    121121 
    122122def ER(radius): 
  • sasmodels/models/vesicle.py

    r068cebd rfa8011eb  
    2424is a flat background level (due for example to incoherent scattering in the 
    2525case of neutrons), and $j_1$ is the spherical bessel function 
    26 $j_1 = (sin(x) - x cos(x))/ x^2$. 
     26$j_1 = (\sin(x) - x \cos(x))/ x^2$. 
    2727 
    2828The functional form is identical to a "typical" core-shell structure, except 
     
    3535thickness = $R_{\text{tot}} - R_{\text{core}}$. 
    3636 
    37 .. figure: img/vesicle_geometry.jpg 
     37.. figure:: img/vesicle_geometry.jpg 
     38 
     39    Vesicle geometry. 
    3840 
    3941The 2D scattering intensity is the same as *P(q)* above, regardless of the 
     
    4850radius for *S(Q)* when *P(Q)* \* *S(Q)* is applied. 
    4951 
    50 .. image:: img/vesicle_1d.jpg 
     52.. figure:: img/vesicle_1d.jpg 
    5153 
    52 *Figure. 1D plot using the default values given in the table 
    53 (w/200 data point). Polydispersity and instrumental resolution normally 
    54 will smear out most of the rapidly oscillating features.* 
     54    1D plot using the default values given in the table (w/200 data point). 
     55    Polydispersity and instrumental resolution normally will smear out most 
     56    of the rapidly oscillating features. 
    5557 
    5658REFERENCE 
Note: See TracChangeset for help on using the changeset viewer.