Changeset 3b571ae in sasmodels

Timestamp:

Mar 22, 2017 6:32:33 PM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

61104c8

Parents:

b00a646

Message:

ellipsoid: fix docs so equations are consistent with code; rearrange code for speed and readability

Location:

sasmodels/models

Files:

• ellipsoid.c (modified) (3 diffs)
• ellipsoid.py (modified) (4 diffs)

sasmodels/models/ellipsoid.c

@@ -3 +3 @@
 double Iqxy(double qx, double qy, double sld, double sld_solvent,
     double radius_polar, double radius_equatorial, double theta, double phi);
-
-static double
-_ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha)
-{
-    double ratio = radius_polar/radius_equatorial;
-    // Using ratio v = Rp/Re, we can expand the following to match the
-    // form given in Guinier (1955)
-    //     r = Re * sqrt(1 + cos^2(T) (v^2 - 1))
-    //       = Re * sqrt( (1 - cos^2(T)) + v^2 cos^2(T) )
-    //       = Re * sqrt( sin^2(T) + v^2 cos^2(T) )
-    //       = sqrt( Re^2 sin^2(T) + Rp^2 cos^2(T) )
-    //
-    // Instead of using pythagoras we could pass in sin and cos; this may be
-    // slightly better for 2D which has already computed it, but it introduces
-    // an extra sqrt and square for 1-D not required by the current form, so
-    // leave it as is.
-    const double r = radius_equatorial
-                     * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0));
-    const double f = sas_3j1x_x(q*r);
-
-    return f*f;
-}
 
 double form_volume(double radius_polar, double radius_equatorial)
@@ -37 +15 @@
     double radius_equatorial)
 {
+    // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955)
+    //     i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT
+    //          = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT
+    //          = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT
+    // u-substitution of
+    //     u = sin, du = cos dT
+    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du
+    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0;
+
     // translate a point in [-1,1] to a point in [0, 1]
+    // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2;
     const double zm = 0.5;
     const double zb = 0.5;
     double total = 0.0;
     for (int i=0;i<76;i++) {
-        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;
-        const double cos_alpha = Gauss76Z[i]*zm + zb;
-        total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha);
+        const double u = Gauss76Z[i]*zm + zb;
+        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one);
+        const double f = sas_3j1x_x(q*r);
+        total += Gauss76Wt[i] * f * f;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -62 +51 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha);
+    const double r = sqrt(square(radius_equatorial*sin_alpha)
+                          + square(radius_polar*cos_alpha));
+    const double f = sas_3j1x_x(q*r);
     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
 
-    return 1.0e-4 * form * s * s;
+    return 1.0e-4 * square(f * s);
 }
 

sasmodels/models/ellipsoid.py

@@ -18 +18 @@
 .. math::
 
-    F(q,\alpha) = \frac{3 \Delta \rho V (\sin[qr(R_p,R_e,\alpha)]
-                - \cos[qr(R_p,R_e,\alpha)])}
-                {[qr(R_p,R_e,\alpha)]^3}
+    F(q,\alpha) = \Delta \rho V \frac{3(\sin qr  - qr \cos qr)}{(qr)^3}
 
-and
+for
 
 .. math::
 
-    r(R_p,R_e,\alpha) = \left[ R_e^2 \sin^2 \alpha
-        + R_p^2 \cos^2 \alpha \right]^{1/2}
+    r = \left[ R_e^2 \sin^2 \alpha + R_p^2 \cos^2 \alpha \right]^{1/2}
 
 
 $\alpha$ is the angle between the axis of the ellipsoid and $\vec q$,
-$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid , $R_p$ is the polar radius along the
-rotational axis of the ellipsoid, $R_e$ is the equatorial radius perpendicular
-to the rotational axis of the ellipsoid and $\Delta \rho$ (contrast) is the
-scattering length density difference between the scatterer and the solvent.
+$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid, $R_p$ is the polar
+radius along the rotational axis of the ellipsoid, $R_e$ is the equatorial
+radius perpendicular to the rotational axis of the ellipsoid and
+$\Delta \rho$ (contrast) is the scattering length density difference between
+the scatterer and the solvent.
 
-For randomly oriented particles:
+For randomly oriented particles use the orientational average,
 
 .. math::
 
-   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}
+   \langle F^2(q) \rangle = \int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)\,d\alpha}
 
+
+computed via substitution of $u=\sin(\alpha)$, $du=\cos(\alpha)\,d\alpha$ as
+
+.. math::
+
+    \langle F^2(q) \rangle = \int_0^1{F^2(q, u)\,du}
+
+with
+
+.. math::
+
+    r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2}
 
 To provide easy access to the orientation of the ellipsoid, we define
@@ -48 +58 @@
 :ref:`cylinder orientation figure <cylinder-angle-definition>`.
 For the ellipsoid, $\theta$ is the angle between the rotational axis
-and the $z$ -axis.
+and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$
+in the $xy$ plane.
 
 NB: The 2nd virial coefficient of the solid ellipsoid is calculated based
@@ -90 +101 @@
 than 500.
 
+Model was also tested against the triaxial ellipsoid model with equal major
+and minor equatorial radii.  It is also consistent with the cyclinder model
+with polar radius equal to length and equatorial radius equal to radius.
+
 References
 ----------
@@ -96 +111 @@
 *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
 Plenum Press, New York, 1987.
+
+Authorship and Verification
+----------------------------
+
+* **Author:** NIST IGOR/DANSE **Date:** pre 2010
+* **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014
+* **Last Modified by:** Paul Kienzle **Date:** March 22, 2017
 """