Changes in sasmodels/models/ellipsoid.c [130d4c7:3b571ae] in sasmodels

File:

• sasmodels/models/ellipsoid.c (modified) (3 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);
 }