Changeset 130d4c7 in sasmodels

Timestamp:

Feb 14, 2017 12:43:14 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:

9ed53e7, a38b065

Parents:

1105286

Message:

restore broken ellipsoid 2-D calculation to 4.0.1 equivalent

File:

• sasmodels/models/ellipsoid.c (modified) (4 diffs)

sasmodels/models/ellipsoid.c

@@ -4 +4 @@
     double radius_polar, double radius_equatorial, double theta, double phi);
 
-double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha);
-double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha)
+static double
+_ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha)
 {
     double ratio = radius_polar/radius_equatorial;
-    // Given the following under the radical:
-    //     1 + sin^2(T) (v^2 - 1)
-    // we can expand to match the form given in Guinier (1955)
-    //     = (1 - sin^2(T)) + v^2 sin^2(T) = cos^2(T) + sin^2(T)
+    // 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
@@ -17 +20 @@
     // leave it as is.
     const double r = radius_equatorial
-                     * sqrt(1.0 + sin_alpha*sin_alpha*(ratio*ratio - 1.0));
+                     * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0));
     const double f = sas_3j1x_x(q*r);
 
@@ -39 +42 @@
     double total = 0.0;
     for (int i=0;i<76;i++) {
-        //const double sin_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;
-        const double sin_alpha = Gauss76Z[i]*zm + zb;
-        total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha);
+        //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);
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -59 +62 @@
     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, sin_alpha);
+    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha);
     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);