Changeset a8bd500 in sasmodels


Ignore:
Timestamp:
Feb 26, 2016 3:05:27 PM (9 years ago)
Author:
gonzalezm
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:
deac08c
Parents:
44bd2be (diff), 139c528 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' of ​https://github.com/SasView/sasmodels

Files:
17 added
1 deleted
17 edited
3 moved

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 

  • sasmodels/convert.py

    r34d6cab r44bd2be  
    132132        _remove_pd(oldpars, 'num_pearls', name) 
    133133        _remove_pd(oldpars, 'thick_string', name) 
     134    elif name == 'core_shell_parallelepiped': 
     135        _remove_pd(oldpars, 'rimA', name) 
     136        _remove_pd(oldpars, 'rimB', name) 
     137        _remove_pd(oldpars, 'rimC', name) 
    134138    elif name == 'rpa': 
    135139        # convert scattering lengths from femtometers to centimeters 

  • sasmodels/models/parallelepiped.c

    rc5b7d07 rdeb7ee0  
    99    double argA,argB,argC,tmp1,tmp2,tmp3; 
    1010    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    11     argA = a*ala/2.0; 
    12     argB = b*alb/2.0; 
    13     argC = c*alc/2.0; 
     11    argA = 0.5*a*ala; 
     12    argB = 0.5*b*alb; 
     13    argC = 0.5*c*alc; 
    1414    if(argA==0.0) { 
    1515        tmp1 = 1.0; 
    … …  
    3131} 
    3232 
     33 
    3334double form_volume(double a_side, double b_side, double c_side) 
    3435{ 
    3536    return a_side * b_side * c_side; 
    3637} 
     38 
    3739 
    3840double Iq(double q, 
    … …  
    4345    double c_side) 
    4446{ 
    45     double tmp1, tmp2, yyy; 
     47    double tmp1, tmp2; 
    4648     
    4749    double mu = q * b_side; 
    … …  
    5052    double a_scaled = a_side / b_side; 
    5153    double c_scaled = c_side / b_side; 
    52      
    53     // outer integral (with 76 gauss points), integration limits = 0, 1 
     54         
     55    //Order of integration 
     56    int nordi=76;                                
     57    int nordj=76; 
     58 
     59    // outer integral (with nordi gauss points), integration limits = 0, 1 
    5460    double summ = 0; //initialize integral 
    5561 
    56     for( int i=0; i<76; i++) { 
     62    for( int i=0; i<nordi; i++) { 
    5763                 
    58         // inner integral (with 76 gauss points), integration limits = 0, 1 
     64        // inner integral (with nordj gauss points), integration limits = 0, 1 
    5965         
    6066        double summj = 0.0; 
    61         double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );              
     67            double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );          
    6268                 
    63         for(int j=0; j<76; j++) { 
     69            for(int j=0; j<nordj; j++) { 
    6470 
    6571            double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    6672            double mudum = mu * sqrt(1.0-sigma*sigma); 
    67             double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu); 
    68             double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu); 
     73                double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu); 
     74                double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu); 
    6975            if(arg1==0.0) { 
    7076                tmp1 = 1.0; 
    … …  
    7783                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; 
    7884            } 
    79             yyy = Gauss76Wt[j] * tmp1 * tmp2; 
    8085 
    81             summj += yyy; 
     86            summj += Gauss76Wt[j] * tmp1 * tmp2; 
    8287        } 
    8388                 
    … …  
    9297        } 
    9398                 
    94         // sum of outer integral 
    95         yyy = Gauss76Wt[i] * answer; 
    96         summ += yyy; 
     99            // sum of outer integral 
     100        summ += Gauss76Wt[i] * answer; 
    97101         
    98102    }    
    99103    
    100104    const double vd = (sld-solvent_sld) * form_volume(a_side, b_side, c_side); 
     105     
     106    // convert from [1e-12 A-1] to [cm-1] and 0.5 factor for outer integral 
    101107    return 1.0e-4 * 0.5 * vd * vd * summ; 
    102108     

  • sasmodels/models/parallelepiped.py

    reb69cce rdeb7ee0  
    77---------- 
    88 
    9 This model provides the form factor, $P(q)$, for a rectangular parallelepiped 
    10 (below) where the form factor is normalized by the volume of the 
    11 parallelepiped. If you need to apply polydispersity, see also 
    12 rectangular_prism_. 
    13  
    14 The calculated form factor is: 
    15  
    16 .. math:: 
    17  
    18     P(q) = \frac{\text{scale}}{V} F^2(q) + \text{background} 
    19  
    20 where the volume $V = A B C$ and the averaging $\left<\ldots\right>$ is 
    21 applied over all orientations for 1D. 
     9| This model calculates the scattering from a rectangular parallelepiped (:num:`Figure #parallelepiped-image`). 
     10| If you need to apply polydispersity, see also :ref:`rectangular-prism`. 
     11 
     12.. _parallelepiped-image: 
    2213 
    2314.. figure:: img/parallelepiped.jpg 
    … …  
    2516   Parallelepiped with the corresponding definition of sides. 
    2617 
    27 The edge of the solid must satisfy the condition that $A < B < C$. 
    28 Then, assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
    29 form factor is 
    30  
    31 .. math:: 
    32  
    33     P(q) = \frac{\text{scale}}{V}\int_0^1 
    34         \phi\left(\mu \sqrt{1-\sigma^2},a\right) 
     18.. note:: 
     19 
     20   The edge of the solid must satisfy the condition that $A < B < C$. 
     21   This requirement is not enforced in the model, so it is up to the 
     22   user to check this during the analysis. 
     23 
     24The 1D scattering intensity $I(q)$ is calculated as: 
     25 
     26.. Comment by Miguel Gonzalez: 
     27   I am modifying the original text because I find the notation a little bit 
     28   confusing. I think that in most textbooks/papers, the notation P(Q) is 
     29   used for the form factor (adim, P(Q=0)=1), although F(q) seems also to 
     30   be used. But here (as for many other models), P(q) is used to represent 
     31   the scattering intensity (in cm-1 normally). It would be good to agree on 
     32   a common notation. 
     33 
     34.. math:: 
     35 
     36    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 \left< P(q, \alpha) \right> 
     37    + \text{background} 
     38 
     39where the volume $V = A B C$, the contrast is defined as 
     40:math:`\Delta\rho = \rho_{\textstyle p} - \rho_{\textstyle solvent}`, 
     41$P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented 
     42at an angle $\alpha$ (angle between the long axis C and :math:`\vec q`), 
     43and the averaging $\left<\ldots\right>$ is applied over all orientations. 
     44 
     45Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
     46form factor is given by (Mittelbach and Porod, 1961) 
     47 
     48.. math:: 
     49 
     50    P(q, \alpha) = \int_0^1 \phi_Q\left(\mu \sqrt{1-\sigma^2},a\right) 
    3551        \left[S(\mu c \sigma/2)\right]^2 d\sigma 
    3652 
    … …  
    3955.. math:: 
    4056 
    41     \phi(\mu,a) = \int_0^1 
    42         \left\{S\left[\frac{\mu}{2}\cos(\frac{\pi}{2}u)\right] 
    43                S\left[\frac{\mu a}{2}\sin(\frac{\pi}{2}u)\right] 
     57    \phi_Q(\mu,a) = \int_0^1 
     58        \left\{S\left[\frac{\mu}{2}\cos\left(\frac{\pi}{2}u\right)\right] 
     59               S\left[\frac{\mu a}{2}\sin\left(\frac{\pi}{2}u\right)\right] 
    4460               \right\}^2 du 
    4561 
    … …  
    4864    \mu = qB 
    4965 
    50 and the contrast is defined as 
    51  
    52 .. math:: 
    53  
    54     \Delta\rho = \rho_{\textstyle p} - \rho_{\textstyle solvent} 
    55  
    56 The scattering intensity per unit volume is returned in units of |cm^-1|; 
    57 i.e., $I(q) = \phi P(q)$. 
    58  
    59 NB: The 2nd virial coefficient of the parallelpiped is calculated based on 
     66 
     67The scattering intensity per unit volume is returned in units of |cm^-1|. 
     68 
     69NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    6070the averaged effective radius $(=\sqrt{A B / \pi})$ and 
    6171length $(= C)$ values, and used as the effective radius for 
    … …  
    6575three angles $\theta$, $\phi$ and $\Psi$. The definition of $\theta$ and 
    6676$\phi$ is the same as for the cylinder model (see also figures below). 
    67 The angle $\Psi$ is the rotational angle around the $C$ axis against 
    68 the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel 
    69 to the $x$-axis of the detector. 
    70  
     77 
     78.. Comment by Miguel Gonzalez: 
     79   The following text has been commented because I think there are two 
     80   mistakes. Psi is the rotational angle around C (but I cannot understand 
     81   what it means against the q plane) and psi=0 corresponds to a||x and b||y. 
     82 
     83   The angle $\Psi$ is the rotational angle around the $C$ axis against 
     84   the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel 
     85   to the $x$-axis of the detector. 
     86 
     87The angle $\Psi$ is the rotational angle around the $C$ axis. 
     88For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis 
     89oriented parallel to the y-axis of the detector with $A$ along the z-axis. 
     90For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated 
     91$\theta$ degrees around $z$ and $\phi$ degrees around $y$, 
     92before doing a final rotation of $\Psi$ degrees around the resulting $C$ to 
     93obtain the final orientation of the parallelepiped. 
     94For example, for $\theta = 0$ and $\phi = 90$, we have that $\Psi = 0$ 
     95corresponds to $A$ along $x$ and $B$ along $y$, 
     96while for $\theta = 90$ and $\phi = 0$, $\Psi = 0$ corresponds to 
     97$A$ along $z$ and $B$ along $x$. 
    7198 
    7299.. _parallelepiped-orientation: 
    73100 
    74 .. figure:: img/orientation.jpg 
     101.. figure:: img/parallelepiped_angles_definition.jpg 
    75102 
    76103    Definition of the angles for oriented parallelepipeds. 
    77104 
    78 .. figure:: img/orientation2.jpg 
     105.. figure:: img/parallelepiped_angles_examples.jpg 
    79106 
    80107    Examples of the angles for oriented parallelepipeds against the detector plane. 
     108 
     109For a given orientation of the parallelepiped, the 2D form factor is calculated as 
     110 
     111.. math:: 
     112 
     113    P(q_x, q_y) = \left[\frac{sin(qA\cos\alpha/2)}{(qA\cos\alpha/2)}\right]^2 
     114                  \left[\frac{sin(qB\cos\beta/2)}{(qB\cos\beta/2)}\right]^2 
     115                  \left[\frac{sin(qC\cos\gamma/2)}{(qC\cos\gamma/2)}\right]^2 
     116 
     117with 
     118 
     119.. math:: 
     120 
     121    \cos\alpha = \hat A \cdot \hat q, \; 
     122    \cos\beta  = \hat B \cdot \hat q, \; 
     123    \cos\gamma = \hat C \cdot \hat q 
     124 
     125and the scattering intensity as: 
     126 
     127.. math:: 
     128 
     129    I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y) + \text{background} 
     130 
     131.. Comment by Miguel Gonzalez: 
     132   This reflects the logic of the code, as in parallelepiped.c the call 
     133   to _pkernel returns P(q_x, q_y) and then this is multiplied by V^2 * (delta rho)^2. 
     134   And finally outside parallelepiped.c it will be multiplied by scale, normalized by 
     135   V and the background added. But mathematically it makes more sense to write 
     136   I(q_x, q_y) = \text{scale} V \Delta\rho^2 P(q_x, q_y) + \text{background}, 
     137   with scale being the volume fraction. 
    81138 
    82139 
    … …  
    97154   Comparison between 1D and averaged 2D. 
    98155 
    99 This model reimplements the form factor calculations implemented in a c-library 
     156This model is based on form factor calculations implemented in a c-library 
    100157provided by the NIST Center for Neutron Research (Kline, 2006). 
    101158 
    102159""" 
    103160 
     161import numpy as np 
    104162from numpy import pi, inf, sqrt 
    105163 
    … …  
    107165title = "Rectangular parallelepiped with uniform scattering length density." 
    108166description = """ 
    109      P(q)= scale/V*integral from 0 to 1 of ... 
     167    I(q)= scale*V*(sld - solvent_sld)^2*P(q,alpha)+background 
     168        P(q,alpha) = integral from 0 to 1 of ... 
    110169           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma 
    111  
     170        with 
    112171            phi(mu,a) = integral from 0 to 1 of .. 
    113         (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du 
     172            (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du 
    114173            S(x) = sin(x)/x 
    115         mu = q*B 
     174            mu = q*B 
     175        V: Volume of the rectangular parallelepiped 
     176        alpha: angle between the long axis of the  
     177            parallelepiped and the q-vector for 1D 
    116178""" 
    117 category = "shape:parallelpiped" 
     179category = "shape:parallelepiped" 
    118180 
    119181#             ["name", "units", default, [lower, upper], "type","description"], 
    … …  
    139201 
    140202def ER(a_side, b_side, c_side): 
     203    """ 
     204        Return effective radius (ER) for P(q)*S(q) 
     205    """ 
    141206 
    142207    # surface average radius (rough approximation) 
    143208    surf_rad = sqrt(a_side * b_side / pi) 
    144209 
    145     # DiamCyl recoded here (to check and possibly put in a library?) 
    146     a = surf_rad 
    147     b = 0.5 * c_side 
    148     t1 = a * a * b 
    149     t2 = 1.0 + (b / a) * (1.0 + a / b / 2.0) * (1.0 + pi * a / b / 2.0) 
    150     ddd = 3.0 * t1 * t2 
    151  
     210    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    152211    return 0.5 * (ddd) ** (1. / 3.) 
     212 
     213# VR defaults to 1.0 
    153214 
    154215# parameters for demo 
    … …  
    171232               sld='sldPipe', solvent_sld='sldSolv') 
    172233 
     234 
     235qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
     236tests = [[{}, 0.2, 0.17658004974], 
     237         [{}, [0.2], [0.17658004974]], 
     238         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.00460296014], 
     239         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.00460296014]], 
     240        ] 
     241del qx, qy  # not necessary to delete, but cleaner 
Note: See TracChangeset for help on using the changeset viewer.