Changeset 68dd6a9 in sasmodels for sasmodels/models/triaxial_ellipsoid.c

Timestamp:

Mar 22, 2017 9:46:16 PM (7 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:

d86967b

Parents:

01ea374

Message:

triaxial ellipsoid: adjust to new asymmetric orientation definition; reimplement 1D model from Finnigan and Jacobs, 1971. Refs #890.

File:

• sasmodels/models/triaxial_ellipsoid.c (modified) (2 diffs)

sasmodels/models/triaxial_ellipsoid.c

@@ -20 +20 @@
     double radius_polar)
 {
-    double sn, cn;
-    // translate a point in [-1,1] to a point in [0, 1]
-    const double zm = 0.5;
-    const double zb = 0.5;
+    const double pa = square(radius_equat_minor/radius_equat_major) - 1.0;
+    const double pc = square(radius_polar/radius_equat_major) - 1.0;
+    // translate a point in [-1,1] to a point in [0, pi/2]
+    const double zm = M_PI_4;
+    const double zb = M_PI_4;
     double outer = 0.0;
     for (int i=0;i<76;i++) {
-        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;
-        const double x = 0.5*(Gauss76Z[i] + 1.0);
-        SINCOS(M_PI_2*x, sn, cn);
-        const double acosx2 = radius_equat_minor*radius_equat_minor*cn*cn;
-        const double bsinx2 = radius_equat_major*radius_equat_major*sn*sn;
-        const double c2 = radius_polar*radius_polar;
+        //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2;
+        const double phi = Gauss76Z[i]*zm + zb;
+        const double pa_sinsq_phi = pa*square(sin(phi));
 
         double inner = 0.0;
+        const double um = 0.5;
+        const double ub = 0.5;
         for (int j=0;j<76;j++) {
-            const double ysq = square(Gauss76Z[j]*zm + zb);
-            const double t = q*sqrt(acosx2 + bsinx2*(1.0-ysq) + c2*ysq);
-            const double fq = sas_3j1x_x(t);
-            inner += Gauss76Wt[j] * fq * fq ;
+            // translate a point in [-1,1] to a point in [0, 1]
+            const double usq = square(Gauss76Z[j]*um + ub);
+            const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq);
+            const double fq = sas_3j1x_x(q*r);
+            inner += Gauss76Wt[j] * fq * fq;
         }
-        outer += Gauss76Wt[i] * 0.5 * inner;
+        outer += Gauss76Wt[i] * inner;  // correcting for dx later
     }
-    // translate dx in [-1,1] to dx in [lower,upper]
-    const double fqsq = outer*zm;
+    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi
+    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI);
     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);
     return 1.0e-4 * s * s * fqsq;
@@ -58 +59 @@
     double psi)
 {
-    double q, calpha, cmu, cnu;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu);
+    double q, xhat, yhat, zhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
 
-    const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu
-                          + radius_equat_major*radius_equat_major*cmu*cmu
-                          + radius_polar*radius_polar*calpha*calpha);
-    const double fq = sas_3j1x_x(t);
+    const double r = sqrt(square(radius_equat_minor*xhat)
+                          + square(radius_equat_major*yhat)
+                          + square(radius_polar*zhat));
+    const double fq = sas_3j1x_x(q*r);
     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar);