Changeset 6e7ba14 in sasmodels for sasmodels/models

Timestamp:

Aug 20, 2018 8:52:21 AM (6 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:

aa44a6a

Parents:

bad3093

Message:

allow for different forms of effective_radius

Location:

sasmodels/models

Files:

• ellipsoid.c (modified) (1 diff)
• ellipsoid.py (modified) (1 diff)

sasmodels/models/ellipsoid.c

@@ -4 +4 @@
     return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial;
 }
+
+
+static double
+radius_from_volume(double radius_polar, double radius_equatorial)
+{
+    return cbrt(radius_polar*radius_equatorial*radius_equatorial);
+}
+
+static double
+radius_from_curvature(double radius_polar, double radius_equatorial)
+{
+    // Trivial cases
+    if (radius_polar == radius_equatorial) return radius_polar;
+    if (radius_polar * radius_equatorial == 0.)  return 0.;
+
+    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    const double ratio = (radius_polar < radius_equatorial 
+                          ? radius_polar / radius_equatorial
+                          : radius_equatorial / radius_polar);
+    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) * radius_polar * radius_equatorial * radius_equatorial;
+    return 0.5 * cbrt(ddd);
+}
+
+static double
+effective_radius(int mode, double radius_polar, double radius_equatorial)
+{
+    if (mode == 1) {
+        return radius_from_curvature(radius_polar, radius_equatorial);
+    } else if (mode == 2) {
+        return radius_from_volume(radius_polar, radius_equatorial);
+    } else if (mode == 3) {
+        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial);
+    } else {
+        return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial);
+    }
+}
+
 
 static void

sasmodels/models/ellipsoid.py

@@ -164 +164 @@
 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"]
 have_Fq = True
-
-def ER(radius_polar, radius_equatorial):
-    # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
-    ee = np.empty_like(radius_polar)
-    idx = radius_polar > radius_equatorial
-    ee[idx] = (radius_polar[idx] ** 2 - radius_equatorial[idx] ** 2) / radius_polar[idx] ** 2
-    idx = radius_polar < radius_equatorial
-    ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2
-    idx = radius_polar == radius_equatorial
-    ee[idx] = 2 * radius_polar[idx]
-    valid = (radius_polar * radius_equatorial != 0)
-    bd = 1.0 - ee[valid]
-    e1 = np.sqrt(ee[valid])
-    b1 = 1.0 + np.arcsin(e1) / (e1 * np.sqrt(bd))
-    bL = (1.0 + e1) / (1.0 - e1)
-    b2 = 1.0 + bd / 2 / e1 * np.log(bL)
-    delta = 0.75 * b1 * b2
-
-    ddd = np.zeros_like(radius_polar)
-    ddd[valid] = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial ** 2
-    return 0.5 * ddd ** (1.0 / 3.0)
+effective_radius_type = ["average curvature", "equivalent sphere", "min radius", "max radius"]
 
 def random():