Changeset 3b571ae in sasmodels
 Timestamp:
 Mar 22, 2017 6:32:33 PM (9 months ago)
 Branches:
 master, boltzmann, costrafo411, doc_update, generic_integration_loop, ticket1043, ticket786
 Children:
 61104c8
 Parents:
 b00a646
 Location:
 sasmodels/models
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

sasmodels/models/ellipsoid.c
r130d4c7 r3b571ae 3 3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 4 4 double radius_polar, double radius_equatorial, double theta, double phi); 5 6 static double7 _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha)8 {9 double ratio = radius_polar/radius_equatorial;10 // Using ratio v = Rp/Re, we can expand the following to match the11 // form given in Guinier (1955)12 // r = Re * sqrt(1 + cos^2(T) (v^2  1))13 // = Re * sqrt( (1  cos^2(T)) + v^2 cos^2(T) )14 // = Re * sqrt( sin^2(T) + v^2 cos^2(T) )15 // = sqrt( Re^2 sin^2(T) + Rp^2 cos^2(T) )16 //17 // Instead of using pythagoras we could pass in sin and cos; this may be18 // slightly better for 2D which has already computed it, but it introduces19 // an extra sqrt and square for 1D not required by the current form, so20 // leave it as is.21 const double r = radius_equatorial22 * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio  1.0));23 const double f = sas_3j1x_x(q*r);24 25 return f*f;26 }27 5 28 6 double form_volume(double radius_polar, double radius_equatorial) … … 37 15 double radius_equatorial) 38 16 { 17 // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955) 18 // i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT 19 // = int_0^pi/2 Phi^2(h a sqrt((1sin^2) + v^2 sin^2) cos dT 20 // = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^21)) cos dT 21 // usubstitution of 22 // u = sin, du = cos dT 23 // i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^21)) du 24 const double v_square_minus_one = square(radius_polar/radius_equatorial)  1.0; 25 39 26 // translate a point in [1,1] to a point in [0, 1] 27 // const double u = Gauss76Z[i]*(upperlower)/2 + (upper+lower)/2; 40 28 const double zm = 0.5; 41 29 const double zb = 0.5; 42 30 double total = 0.0; 43 31 for (int i=0;i<76;i++) { 44 //const double cos_alpha = (Gauss76Z[i]*(upperlower) + upper + lower)/2; 45 const double cos_alpha = Gauss76Z[i]*zm + zb; 46 total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 32 const double u = Gauss76Z[i]*zm + zb; 33 const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 34 const double f = sas_3j1x_x(q*r); 35 total += Gauss76Wt[i] * f * f; 47 36 } 48 37 // translate dx in [1,1] to dx in [lower,upper] … … 62 51 double q, sin_alpha, cos_alpha; 63 52 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 64 const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 53 const double r = sqrt(square(radius_equatorial*sin_alpha) 54 + square(radius_polar*cos_alpha)); 55 const double f = sas_3j1x_x(q*r); 65 56 const double s = (sld  sld_solvent) * form_volume(radius_polar, radius_equatorial); 66 57 67 return 1.0e4 * form * s * s;58 return 1.0e4 * square(f * s); 68 59 } 69 60 
sasmodels/models/ellipsoid.py
r925ad6e r3b571ae 18 18 .. math:: 19 19 20 F(q,\alpha) = \frac{3 \Delta \rho V (\sin[qr(R_p,R_e,\alpha)] 21  \cos[qr(R_p,R_e,\alpha)])} 22 {[qr(R_p,R_e,\alpha)]^3} 20 F(q,\alpha) = \Delta \rho V \frac{3(\sin qr  qr \cos qr)}{(qr)^3} 23 21 24 and 22 for 25 23 26 24 .. math:: 27 25 28 r(R_p,R_e,\alpha) = \left[ R_e^2 \sin^2 \alpha 29 + R_p^2 \cos^2 \alpha \right]^{1/2} 26 r = \left[ R_e^2 \sin^2 \alpha + R_p^2 \cos^2 \alpha \right]^{1/2} 30 27 31 28 32 29 $\alpha$ is the angle between the axis of the ellipsoid and $\vec q$, 33 $V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid , $R_p$ is the polar radius along the 34 rotational axis of the ellipsoid, $R_e$ is the equatorial radius perpendicular 35 to the rotational axis of the ellipsoid and $\Delta \rho$ (contrast) is the 36 scattering length density difference between the scatterer and the solvent. 30 $V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid, $R_p$ is the polar 31 radius along the rotational axis of the ellipsoid, $R_e$ is the equatorial 32 radius perpendicular to the rotational axis of the ellipsoid and 33 $\Delta \rho$ (contrast) is the scattering length density difference between 34 the scatterer and the solvent. 37 35 38 For randomly oriented particles :36 For randomly oriented particles use the orientational average, 39 37 40 38 .. math:: 41 39 42 F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}40 \langle F^2(q) \rangle = \int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)\,d\alpha} 43 41 42 43 computed via substitution of $u=\sin(\alpha)$, $du=\cos(\alpha)\,d\alpha$ as 44 45 .. math:: 46 47 \langle F^2(q) \rangle = \int_0^1{F^2(q, u)\,du} 48 49 with 50 51 .. math:: 52 53 r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2  1\right)\right]^{1/2} 44 54 45 55 To provide easy access to the orientation of the ellipsoid, we define … … 48 58 :ref:`cylinder orientation figure <cylinderangledefinition>`. 49 59 For the ellipsoid, $\theta$ is the angle between the rotational axis 50 and the $z$ axis. 60 and the $z$ axis in the $xz$ plane followed by a rotation by $\phi$ 61 in the $xy$ plane. 51 62 52 63 NB: The 2nd virial coefficient of the solid ellipsoid is calculated based … … 90 101 than 500. 91 102 103 Model was also tested against the triaxial ellipsoid model with equal major 104 and minor equatorial radii. It is also consistent with the cyclinder model 105 with polar radius equal to length and equatorial radius equal to radius. 106 92 107 References 93 108  … … 96 111 *Structure Analysis by SmallAngle XRay and Neutron Scattering*, 97 112 Plenum Press, New York, 1987. 113 114 Authorship and Verification 115  116 117 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 118 * **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014 119 * **Last Modified by:** Paul Kienzle **Date:** March 22, 2017 98 120 """ 99 121
Note: See TracChangeset
for help on using the changeset viewer.