Changeset 63d4dd1 in sasmodels for sasmodels/models/triaxial_ellipsoid.c

Timestamp:

Nov 9, 2018 3:25:10 PM (5 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

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

Children:

7126c04

Parents:

599993b9 (diff), cf3d0ce (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'beta_approx' into ticket-1015-gpu-mem-error

File:

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

sasmodels/models/triaxial_ellipsoid.c

@@ -8 +8 @@
 
 static double
-Iq(double q,
+radius_from_curvature(double radius_equat_minor, double radius_equat_major, double radius_polar)
+{
+    // Trivial cases
+    if (radius_equat_minor == radius_equat_major == radius_polar) return radius_polar;
+    if (radius_equat_minor * radius_equat_major * radius_polar == 0.)  return 0.;
+
+
+    double r_equat_equiv, r_polar_equiv;
+    double radii[3] = {radius_equat_minor, radius_equat_major, radius_polar};
+    double radmax = fmax(radii[0],fmax(radii[1],radii[2]));
+
+    double radius_1 = radmax;
+    double radius_2 = radmax;
+    double radius_3 = radmax;
+
+    for(int irad=0; irad<3; irad++) {
+        if (radii[irad] < radius_1) {
+            radius_3 = radius_2;
+            radius_2 = radius_1;
+            radius_1 = radii[irad];
+            } else {
+                if (radii[irad] < radius_2) {
+                        radius_2 = radii[irad];
+                }
+            }
+    }
+    if(radius_2-radius_1 > radius_3-radius_2) {
+        r_equat_equiv = sqrt(radius_2*radius_3);
+        r_polar_equiv = radius_1;
+    } else  {
+        r_equat_equiv = sqrt(radius_1*radius_2);
+        r_polar_equiv = radius_3;
+    }
+
+    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    const double ratio = (r_polar_equiv < r_equat_equiv
+                          ? r_polar_equiv / r_equat_equiv
+                          : r_equat_equiv / r_polar_equiv);
+    const double e1 = sqrt(1.0 - ratio*ratio);
+    const double b1 = 1.0 + asin(e1) / (e1 * ratio);
+    const double bL = (1.0 + e1) / (1.0 - e1);
+    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL);
+    const double delta = 0.75 * b1 * b2;
+    const double ddd = 2.0 * (delta + 1.0) * r_polar_equiv * r_equat_equiv * r_equat_equiv;
+    return 0.5 * cbrt(ddd);
+}
+
+static double
+radius_from_volume(double radius_equat_minor, double radius_equat_major, double radius_polar)
+{
+    return cbrt(radius_equat_minor*radius_equat_major*radius_polar);
+}
+
+static double
+radius_from_min_dimension(double radius_equat_minor, double radius_equat_major, double radius_polar)
+{
+    const double rad_equat_min = (radius_equat_minor < radius_equat_major ? radius_equat_minor : radius_equat_major);
+    return (rad_equat_min < radius_polar ? rad_equat_min : radius_polar);
+}
+
+static double
+radius_from_max_dimension(double radius_equat_minor, double radius_equat_major, double radius_polar)
+{
+    const double rad_equat_max = (radius_equat_minor < radius_equat_major ? radius_equat_major : radius_equat_minor);
+    return (rad_equat_max > radius_polar ? rad_equat_max : radius_polar);
+}
+
+static double
+effective_radius(int mode, double radius_equat_minor, double radius_equat_major, double radius_polar)
+{
+    switch (mode) {
+    default:
+    case 1: // equivalent biaxial ellipsoid average curvature
+        return radius_from_curvature(radius_equat_minor,radius_equat_major, radius_polar);
+    case 2: // equivalent volume sphere
+        return radius_from_volume(radius_equat_minor,radius_equat_major, radius_polar);
+    case 3: // min radius
+        return radius_from_min_dimension(radius_equat_minor,radius_equat_major, radius_polar);
+    case 4: // max radius
+        return radius_from_max_dimension(radius_equat_minor,radius_equat_major, radius_polar);
+    }
+}
+
+static void
+Fq(double q,
+    double *F1,
+    double *F2,
     double sld,
     double sld_solvent,
@@ -20 +106 @@
     const double zm = M_PI_4;
     const double zb = M_PI_4;
-    double outer = 0.0;
+    double outer_sum_F1 = 0.0;
+    double outer_sum_F2 = 0.0;
     for (int i=0;i<GAUSS_N;i++) {
         //const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2;
@@ -26 +113 @@
         const double pa_sinsq_phi = pa*square(sin(phi));
 
-        double inner = 0.0;
+        double inner_sum_F1 = 0.0;
+        double inner_sum_F2 = 0.0;
         const double um = 0.5;
         const double ub = 0.5;
@@ -34 +122 @@
             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 += GAUSS_W[j] * fq * fq;
+            inner_sum_F1 += GAUSS_W[j] * fq;
+            inner_sum_F2 += GAUSS_W[j] * fq * fq;
         }
-        outer += GAUSS_W[i] * inner;  // correcting for dx later
+        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1;  // correcting for dx later
+        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2;  // correcting for dx later
     }
     // 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 vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar);
-    const double drho = (sld - sld_solvent);
-    return 1.0e-4 * square(vol*drho) * fqsq;
+    outer_sum_F1 *= 0.25;  // = outer*um*zm*8.0/(4.0*M_PI);
+    outer_sum_F2 *= 0.25;  // = outer*um*zm*8.0/(4.0*M_PI);
+
+    const double volume = form_volume(radius_equat_minor, radius_equat_major, radius_polar);
+    const double contrast = (sld - sld_solvent);
+    *F1 = 1.0e-2 * contrast * volume * outer_sum_F1;
+    *F2 = 1.0e-4 * square(contrast * volume) * outer_sum_F2;
 }
+
 
 static double