Adding torus elliptical shell model

sasmodels/models/torus_elliptical_shell.c

@@ -0,0 +1,89 @@
+static double form_volume(double radius, double core_radius, double thickness,
+                          double nu) {
+  double a = core_radius;
+  double b = nu * a;
+  double ao = a + thickness;
+  double bo = nu * ao;
+  double area = M_PI * ao * bo - M_PI * a * b;
+  return 2.0 * M_PI * radius * area;
+}
+
+static double F_torus(double Q, double theta, double R, double x, double nu,
+                      double delta_eta) {
+  // integral_{R - x}^{R + x}
+  // 4 * pi * r * J0(Q * r * sin(theta)) * ganmma
+  // sinx_x(Q * ganmma * cos(theta)) dr
+
+  // Set lower integration bound to R-x (not min(R-x,0)), since a torus has a
+  // hole and R-x > 0 is always valid.
+
+  double gamma, int_r_delta, r, f_total, gamma_arg;
+  const double square_x = square(x);
+  const double Q_sin_theta = Q * sin(theta);
+  const double Q_cos_theta = Q * cos(theta);
+
+  f_total = 0.0;
+
+  if (Q_sin_theta == 0.0) {
+    // If Q * sin(theta) is zero, the integral simplifies to a constant
+    // value since J0(0) = 1 and sin(0) = 0.
+    return 4.0 * M_PI * R * x *
+           delta_eta;  // integral over [R - x, R + x] is just a constant
+  }
+
+  for (int i = 0; i < GAUSS_N; i++) {
+    // translate a point in[-1, 1] to a point in[R - x, R + x]
+    r = GAUSS_Z[i] * x + R;
+
+    // gamma = nu * sqrt(fmax(square_x - square(r - R), 0.0));
+#ifdef USE_OPENCL  // if you can use OpenCL, use the step function.
+    gamma_arg = square_x - square(r - R);
+    gamma = nu * sqrt(fabs(gamma_arg)) * step(0.0, gamma_arg);
+#else
+    gamma = nu * sqrt(fmax(square_x - square(r - R), 0.0));
+#endif
+
+    int_r_delta =
+        r * sas_J0(r * Q_sin_theta) * sin(gamma * Q_cos_theta) / Q_cos_theta;
+    f_total += GAUSS_W[i] * int_r_delta;
+  }
+
+  return 4.0 * M_PI * f_total * x *
+         delta_eta;  // multiply by x to get the integral over [R - x, R + x]
+}
+
+static double Iq_non_correlation(double q, double radius, double core_radius,
+                                 double thickness, double nu, double sld_core,
+                                 double sld_shell, double sld_solvent) {
+  // integral_{0}^{pi / 2}
+  //  | F_torus(Q, theta, R, core_radius + thickness, nu, sld_shell -
+  //  sld_solvent)   // shell
+  //    - F_torus(Q, theta, R, core_radius, nu, sld_core - sld_solvent) //
+  //    core
+  //  |^2 dtheta
+  double F_diff, theta;
+  double I_total = 0.0;
+
+  for (int i = 0; i < GAUSS_N; i++) {
+    // translate a point in[-1, 1] to a point in[0, pi / 2]
+    theta = GAUSS_Z[i] * M_PI_4 + M_PI_4;
+
+    F_diff = F_torus(q, theta, radius, core_radius + thickness, nu,
+                     sld_shell - sld_solvent) -
+             F_torus(q, theta, radius, core_radius, nu, sld_shell - sld_core);
+
+    I_total += GAUSS_W[i] * square(F_diff) * sin(theta);
+  }
+
+  return I_total *
+         M_PI_4;  // multiply by pi/4 to get the integral over [0, pi/2]
+}
+
+static double Iq(double q, double radius, double core_radius, double thickness,
+                 double nu, double sld_core, double sld_shell,
+                 double sld_solvent, double zeta, double kappa) {
+  double Sq = 1 + (kappa / (1 + square(q * zeta)));
+
+  return Sq * Iq_non_correlation(q, radius, core_radius, thickness, nu,
+                                 sld_core, sld_shell, sld_solvent);
+}

sasmodels/models/torus_elliptical_shell.py

@@ -0,0 +1,122 @@
+r"""
+Definition
+----------
+
+This model computes the scattering from a torus with an elliptical tube
+cross-section and a concentric shell.
+
+The geometric parameters are:
+
+* $R$: major (ring) radius of the torus centerline
+* $a$: core minor radius in the radial direction
+* $\nu$: aspect ratio of the elliptical cross-section
+* $t$: shell thickness
+
+.. figure:: img/torus_elliptical_shell_geometry.png
+
+    Schematic geometry of the torus with elliptical tube cross-section.
+
+The same aspect ratio $\nu$ is used for both the core and outer cross-sections,
+so the outer semi-axes are $a+t$ and $\nu(a+t)$.
+
+For a given orientation angle $\theta$ between the torus symmetry axis and
+$\vec q$, the kernel evaluates the amplitude using a numerical integral over
+the tube section:
+
+.. math::
+
+    F(q,\theta; x, \Delta\rho)
+    = 4\pi\,\Delta\rho\int_{R-x}^{R+x}
+      r\,J_0\!\left(qr\sin\theta\right)
+      \frac{\sin\!\left(q\,\gamma\cos\theta\right)}{q\cos\theta}\,dr
+
+with
+
+.. math::
+
+    \gamma = \nu\sqrt{x^2-(r-R)^2}
+
+The core-shell amplitude is formed from outer and inner contributions:
+
+.. math::
+
+    F_{cs}(q,\theta) =
+    F\!\left(q,\theta;a + t,\rho_{shell}-\rho_{solvent}\right)
+    -F\!\left(q,\theta;a,\rho_{shell}-\rho_{core}\right)
+
+and the orientationally averaged intensity is
+
+.. math::
+
+    I(q) = S(q)\int_0^{\pi/2}\!\left|F_{cs}(q,\theta)\right|^2\sin\theta\,d\theta
+
+where
+
+.. math::
+
+    S(q)=1+\frac{\kappa}{1+(q\zeta)^2}
+
+Here $S(q)$ represents the structure factor accounting for
+interparticle correlations.
+
+.. note::
+
+    The C kernel integrates over $r \in [R-x, R+x]$, so the model assumes an
+    open torus hole and is only geometrically meaningful when $R > a+t$.
+    (Probably, you don't call it a donut if the hole is closed!)
+
+References
+----------
+#.  T. Kawaguchi, *J. Appl. Crystallogr*, 34(2001) 580-584
+#.  S. Förster, *J. Phys. Chem.*, 103(1999) 6657-6668
+#.  M. J. Hollamby, *Angew. Chem. Int. Ed.* 55(2016) 9890
+
+Authorship and Verification
+---------------------------
+
+* **Author:** Itsuki Tajima and Yuhei Yamada (Github username: Indigo Carmine, https://orcid.org/0009-0003-9780-4135)
+* **Last Modified by:**
+* **Last Reviewed by:**
+"""
+
+from numpy import inf
+
+# ---- SasView model metadata -------------------------------------------------
+name = "torus_elliptical_shell"
+title = "core-shell torus with elliptical cross-section"
+description = "Core-shell torus with elliptical tube cross-section"
+category = "shape:cylinder"
+parameters = [
+    # name          units        default  [min,   max]  type     description
+    ["radius", "Ang", 100.0, [0, inf], "volume", "Torus major radius R"],
+    [
+        "core_radius",
+        "Ang",
+        5.0,
+        [0, inf],
+        "volume",
+        "Elliptical core minor radius a",
+    ],
+    ["thickness", "Ang", 2.0, [0, inf], "volume", "Shell thickness"],
+    ["nu", "", 1.0, [0.1, 10.0], "volume", "Aspect ratio b/a"],
+    ["sld_core", "1e-6/Ang^2", 0.0, [-inf, inf], "sld", "Core SLD"],
+    ["sld_shell", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", "Shell SLD"],
+    ["sld_solvent", "1e-6/Ang^2", 0.0, [-inf, inf], "sld", "Solvent SLD"],
+    ["zeta", "Ang", 0.0, [-inf, inf], "", "a correlation length"],
+    [
+        "kappa",
+        "",
+        0.0,
+        [-inf, inf],
+        "",
+        "the strength of interactions via the isothermal compressibility",
+    ],
+]
+
+# -- tell sasmodels that a C kernel is provided -------------------------------
+source = [
+    "lib/polevl.c",
+    "lib/sas_j0.c",
+    "lib/gauss76.c",
+    "torus_elliptical_shell.c",
+]  # compiled together with default libs