Changes in / [44bd2be:a8bd500] in sasmodels
- Files:
-
- 8 added
- 2 deleted
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
doc/_extensions/dollarmath.py
r19dcb933 r103ea45 12 12 import re 13 13 14 _dollar = re.compile(r"(?:^|(?<=\s ))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;\\]))")14 _dollar = re.compile(r"(?:^|(?<=\s|[(]))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;)\\]))") 15 15 _notdollar = re.compile(r"\\[$]") 16 16 … … 51 51 assert replace_dollar(u"dollar\$ escape")==u"dollar$ escape" 52 52 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" 53 54 assert replace_dollar(u"emb\ $ed$\ ed")==u"emb\ :math:`ed`\ ed" 54 55 assert replace_dollar(u"$first$a")==u"$first$a" 55 56 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)" 56 61 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" 57 64 58 65 if __name__ == "__main__": -
sasmodels/generate.py
r0a4628d rfa8011eb 191 191 The function :func:`make` loads the metadata from the module and returns 192 192 the kernel source. The function :func:`doc` extracts the doc string 193 and adds the parameter table to the top. The function :func:` sources`193 and adds the parameter table to the top. The function :func:`model_sources` 194 194 returns a list of files required by the model. 195 195 """ … … 206 206 import numpy as np 207 207 208 __all__ = ["make", "doc", "sources", "convert_type"]208 #__all__ = ["make", "doc", "model_sources", "convert_type"] 209 209 210 210 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') -
sasmodels/kernel_template.c
r840b859 r960cd80 12 12 #ifndef USE_OPENCL 13 13 # 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) 18 18 #include <limits> 19 19 #include <float.h> 20 20 #define kernel extern "C" __declspec( dllexport ) 21 21 inline double trunc(double x) { return x>=0?floor(x):-floor(-x); } 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 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 50 50 #else 51 51 #define kernel extern "C" … … 260 260 const double vol_weight = VOLUME_WEIGHT_PRODUCT; 261 261 vol += vol_weight*form_volume(VOLUME_PARAMETERS); 262 #endif263 262 norm_vol += vol_weight; 263 #endif 264 264 } 265 265 //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); } -
sasmodels/models/capped_cylinder.c
rc138211 r139c528 22 22 const double upper = 1.0; 23 23 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); 24 26 // cos term in integral is: 25 27 // cos (q (R t - h + L/2) cos(alpha)) … … 36 38 // translate a point in [-1,1] to a point in [lower,upper] 37 39 //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; 39 41 const double radical = 1.0 - t*t; 40 42 const double arg = qrst*sqrt(radical); // cap bessel function arg -
sasmodels/models/core_shell_bicelle.py
r8007311 rfa8011eb 7 7 The form factor is normalized by the particle volume. 8 8 9 .. _core-shell- cylinder-geometry:9 .. _core-shell-bicelle-geometry: 10 10 11 11 .. figure:: img/core_shell_bicelle_geometry.png -
sasmodels/models/core_shell_sphere.py
r8c9dbc9 rfa8011eb 15 15 16 16 .. math:: 17 17 18 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}+ 18 19 V_s(\rho_s-\rho_{solv})\frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3}\right] … … 41 42 our model and the output of the NIST software. 42 43 43 .. image:: img/core_shell_sphere_1d.jpg44 .. figure:: img/core_shell_sphere_1d.jpg 44 45 45 Figure 1:Comparison of the SasView scattering intensity for a core-shell sphere with46 Comparison of the SasView scattering intensity for a core-shell sphere with 46 47 the output of the NIST SANS analysis software. The parameters were set to: 47 48 *scale* = 1.0, *radius* = 60 , *contrast* = 1e-6 |Ang^-2|, and -
sasmodels/models/ellipsoid.py
r9c461c7 r431caae 111 111 ---------- 112 112 113 L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum, 114 New York, 1987. 113 L A Feigin and D I Svergun. 114 *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 115 Plenum Press, New York, 1987. 115 116 """ 116 117 -
sasmodels/models/elliptical_cylinder.py
rb7c2fce rfa8011eb 10 10 to any of the orientation angles, and also for the minor radius and the ratio of the ellipse radii. 11 11 12 .. image:: img/elliptical_cylinder_geometry.gif12 .. figure:: img/elliptical_cylinder_geometry.gif 13 13 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$). 15 15 16 16 The function calculated is 17 17 18 18 .. math:: 19 19 20 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} 20 21 … … 22 23 23 24 .. math:: 25 24 26 F(\mathbf{q},\alpha,\psi)=2\frac{J_1(a)\sin(b)}{ab} 25 27 \\ … … 40 42 P(q) = scale <F^2> / V 41 43 42 The returned value is scaled to units of |cm^-1|.43 44 44 To provide easy access to the orientation of the elliptical cylinder, we define the axis of the cylinder using two 45 45 angles |theta|, |phi| and |bigpsi|. As for the case of the cylinder, the angles |theta| and |phi| are defined on … … 49 49 All angle parameters are valid and given only for 2D calculation; ie, an oriented system. 50 50 51 .. image:: img/elliptical_cylinder_geometry_2d.jpg51 .. figure:: img/elliptical_cylinder_geometry_2d.jpg 52 52 53 *Figure. Definition of angles for 2D*53 Definition of angles for 2D 54 54 55 .. image:: img/core_shell_bicelle_fig2.jpg55 .. figure:: img/core_shell_bicelle_fig2.jpg 56 56 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. 58 58 59 59 NB: The 2nd virial coefficient of the cylinder is calculated based on the averaged radius (= sqrt(*r_minor*\ :sup:`2` \* *r_ratio*)) … … 61 61 62 62 63 .. image:: img/elliptical_cylinder_comparison_1d.jpg63 .. figure:: img/elliptical_cylinder_comparison_1d.jpg 64 64 65 *Figure. 1D plot using the default values (w/1000 data point).*65 1D plot using the default values (w/1000 data point). 66 66 67 67 Validation … … 73 73 and 76 degrees are taken for the angles of |theta|, |phi|, and |bigpsi| respectively). 74 74 75 .. image:: img/elliptical_cylinder_validation_1d.gif75 .. figure:: img/elliptical_cylinder_validation_1d.gif 76 76 77 *Figure. Comparison between 1D and averaged 2D.*77 Comparison between 1D and averaged 2D. 78 78 79 79 In the 2D average, more binning in the angle |phi| is necessary to get the proper result. The following figure shows 80 80 the results of the averaging by varying the number of angular bins. 81 81 82 .. image:: img/elliptical_cylinder_averaging.gif82 .. figure:: img/elliptical_cylinder_averaging.gif 83 83 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. 85 85 86 86 Reference -
sasmodels/models/guinier_porod.py
r21d1031 rfa8011eb 51 51 q = \sqrt{q_x^2+q_y^2} 52 52 53 .. image:: img/guinier_porod_model.jpg53 .. figure:: img/guinier_porod_model.jpg 54 54 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$. 56 56 57 57 -
sasmodels/models/hardsphere.py
r093f754 r97e6d3c 69 69 return(HARDSPH); 70 70 } 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; 73 76 X=fabs(q*effect_radius*2.0); 74 77 … … 77 80 return(HARDSPH); 78 81 } 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; 82 86 G=0.5*volfraction*A; 83 87 84 88 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. 87 91 // 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 90 99 // combining the terms makes things worse at smallest Q in single precision 91 100 //FF = (8-0.8*X2)*A +(3.0-X2/3.)*2*B + (4+0.5*X2)*G +(A/35. +B/40. +G/50.)*X4; 92 101 // 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; 94 103 HARDSPH= 1./(1. + volfraction*FF ); 95 104 return(HARDSPH); 96 105 } 106 X4=X2*X2; 97 107 SINCOS(X,S,C); 98 108 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 ; 101 114 HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); 102 115 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 104 123 // 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 ); 105 124 // HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); -
sasmodels/models/line.py
re66075f rfa8011eb 13 13 .. note:: 14 14 For 2D plots intensity has different definition than other shape independent models 15 15 16 .. math:: 16 17 I(q) = I(qx) \cdot I(qy) 17 18 .. figure:: None19 18 20 19 References -
sasmodels/models/rpa.py
r8dd6914 rfa8011eb 43 43 component. 44 44 45 .. figure:: img/ image215.jpg45 .. figure:: img/rpa_1d.jpg 46 46 47 47 1D plot using the default values (w/500 data points). -
sasmodels/models/spherepy.py
reb69cce rd2950f4 118 118 g[low] = sqrt(1 - dlow2 / 4.) * (1 + dlow2 / 8.) + dlow2 / 2.*(1 - dlow2 / 16.) * log(dlow / (2. + sqrt(4. - dlow2))) 119 119 return g 120 sesans.vectorized = True # sesans accepts an darray of z values120 sesans.vectorized = True # sesans accepts an array of z values 121 121 122 122 def ER(radius): -
sasmodels/models/vesicle.py
r068cebd rfa8011eb 24 24 is a flat background level (due for example to incoherent scattering in the 25 25 case of neutrons), and $j_1$ is the spherical bessel function 26 $j_1 = ( sin(x) - xcos(x))/ x^2$.26 $j_1 = (\sin(x) - x \cos(x))/ x^2$. 27 27 28 28 The functional form is identical to a "typical" core-shell structure, except … … 35 35 thickness = $R_{\text{tot}} - R_{\text{core}}$. 36 36 37 .. figure: img/vesicle_geometry.jpg 37 .. figure:: img/vesicle_geometry.jpg 38 39 Vesicle geometry. 38 40 39 41 The 2D scattering intensity is the same as *P(q)* above, regardless of the … … 48 50 radius for *S(Q)* when *P(Q)* \* *S(Q)* is applied. 49 51 50 .. image:: img/vesicle_1d.jpg52 .. figure:: img/vesicle_1d.jpg 51 53 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. 55 57 56 58 REFERENCE
Note: See TracChangeset
for help on using the changeset viewer.