Changeset d277229 in sasmodels for sasmodels/models/core_shell_ellipsoid.c

Timestamp:

Sep 7, 2018 5:29:38 AM (6 years ago)

Author:

grethevj

Branches:

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

Children:

3c60146

Parents:

2a12351b

Message:

Models updated to include choices for effective interaction radii

File:

• sasmodels/models/core_shell_ellipsoid.c (modified) (1 diff)

sasmodels/models/core_shell_ellipsoid.c

@@ -36 +36 @@
     double vol = M_4PI_3*equat_shell*equat_shell*polar_shell;
     return vol;
+}
+
+static double
+radius_from_volume(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
+{
+    const double volume_ellipsoid = form_volume(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    return cbrt(0.75*volume_ellipsoid/M_PI);
+}
+
+static double
+radius_from_curvature(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
+{
+    // Trivial cases
+    if (x_core == x_polar_shell == 1.0) return radius_equat_core + thick_shell;
+    if ((radius_equat_core + thick_shell)*(radius_equat_core*x_core + thick_shell*x_polar_shell) == 0.)  return 0.;
+
+    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    const double radius_equat_tot = radius_equat_core + thick_shell;
+    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell;
+    const double ratio = (radius_polar_tot < radius_equat_tot
+                          ? radius_polar_tot / radius_equat_tot
+                          : radius_equat_tot / radius_polar_tot);
+    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_tot * radius_equat_tot * radius_equat_tot;
+    return 0.5 * cbrt(ddd);
+}
+
+static double
+effective_radius(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
+{
+    if (mode == 1) {
+        return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    } else if (mode == 2) {
+        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    } else if (mode == 3) {
+        const double radius_equat_tot = radius_equat_core + thick_shell;
+        const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell;
+        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot);
+    } else {
+        const double radius_equat_tot = radius_equat_core + thick_shell;
+        const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell;
+        return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot);
+    }
 }