Adding torus elliptical shell model

sasmodels/models/torus_elliptical_shell.c

@@ -0,0 +1,137 @@
+static double form_volume(double radius, double radius_core, double thickness,
+                          double nu_core, double nu_shell) {
+  double ao = radius_core + thickness;
+  double bo = nu_core * radius_core + nu_shell * thickness;
+  double area = M_PI * ao * bo;
+  return 2.0 * M_PI * radius * area;
+}
+
+static double radius_from_volume(double radius_core, double thickness,
+                                 double nu_core, double nu_shell) {
+  return cbrt(form_volume(1.0, radius_core, thickness, nu_core, nu_shell) /
+              M_4PI_3);
+}
+
+static double radius_from_diagonal(double radius, double radius_core,
+                                   double thickness, double nu_core,
+                                   double nu_shell) {
+  return radius + radius_core + thickness;
+}
+
+static double max_radius(double radius_core, double thickness, double nu_core,
+                         double nu_shell, double torus_radius) {
+  // equatorial radius of the outer ellipse
+  double r_e = radius_core + thickness;
+
+  // polar radius of the outer ellipse
+  double r_p = nu_core * radius_core + nu_shell * thickness;
+
+  double denom = square(r_p) - square(r_e);
+
+  if (denom > 0.0) {
+    double cos_theta = (torus_radius * r_e) / denom;
+
+    if (fabs(cos_theta) <= 1.0) {
+      double nu = r_p / r_e;
+
+      double r_max_sq = square(torus_radius) + square(r_e) +
+                        square(r_p) * ((square(nu) - 1.0) / square(nu));
+
+      return sqrt(r_max_sq);
+    }
+  }
+
+  return torus_radius + r_e;
+}
+
+static double radius_effective(int mode, double radius, double radius_core,
+                               double thickness, double nu_core,
+                               double nu_shell) {
+  switch (mode) {
+    case 1:
+      return radius_from_diagonal(radius, radius_core, thickness, nu_core,
+                                  nu_shell);
+    case 2:
+      return radius_from_volume(radius_core, thickness, nu_core, nu_shell);
+    case 3:
+      return max_radius(radius_core, thickness, nu_core, nu_shell, radius);
+    case 4:
+    default:
+      return radius;
+  }
+}
+
+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 = 0, int_r_delta = 0, r = 0, f_total = 0;
+  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;
+
+  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(square_x - square(r - R));
+
+    int_r_delta =
+        r * sas_J0(r * Q_sin_theta) * sas_sinx_x(Q_cos_theta * gamma) * gamma;
+    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 void Fq(double q, double* F1, double* F2, double radius,
+               double radius_core, double thickness, double nu_core,
+               double nu_shell, double sld_core, double sld_shell,
+               double sld_solvent) {
+  // F2 = integral_{0}^{pi / 2}
+  //  | F_torus(Q, theta, R, radius_core + thickness, nu, sld_shell -
+  //  sld_solvent)   // shell
+  //    - F_torus(Q, theta, R, radius_core, nu, sld_core - sld_solvent) //
+  //    core
+  //  |^2 dtheta
+  double F_diff = 0.0, theta = 0.0, F1_total = 0.0, F2_total = 0.0;
+
+  double nu_outer =
+      (nu_shell * thickness + nu_core * radius_core) /
+      (thickness + radius_core);  // aspect ratio of the outer ellipse
+
+  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, radius_core + thickness, nu_outer,
+                sld_shell - sld_solvent) -
+        F_torus(q, theta, radius, radius_core, nu_core, sld_shell - sld_core);
+
+    F1_total += GAUSS_W[i] * F_diff * sin(theta);
+    F2_total += GAUSS_W[i] * square(F_diff) * sin(theta);
+  }
+
+  // multiply by pi/4 to get the integral over [0, pi/2]
+
+  *F1 = 1e-2 * F1_total * M_PI_4;
+  *F2 = 1e-4 * F2_total * M_PI_4;
+}
+
+static double Iq(double q, double radius, double radius_core, double thickness,
+                 double nu_core, double nu_shell, double sld_core,
+                 double sld_shell, double sld_solvent) {
+  double F1 = 0.0, F2 = 0.0;
+  Fq(q, &F1, &F2, radius, radius_core, thickness, nu_core, nu_shell, sld_core,
+     sld_shell, sld_solvent);
+  return F2;
+}

sasmodels/models/torus_elliptical_shell.py

@@ -0,0 +1,117 @@
+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_{core}$: aspect ratio of the elliptical core cross-section
+* $\nu_{shell}$: aspect ratio of the elliptical shell cross-section
+* $t$: shell thickness
+
+.. figure:: img/torus_elliptical_shell_geometry.png
+
+    Schematic geometry of the torus with elliptical tube cross-section.
+
+The core and shell have independent aspect ratios, so the core semi-axes are
+$a$ and $\nu_{core} \cdot a$, and the outer semi-axes are $a+t$ and
+$\nu_{core} \cdot a + \nu_{shell} \cdot 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, \nu)
+    = 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},\nu_{outer}\right)
+    -F\!\left(q,\theta;a,\rho_{shell}-\rho_{core},\nu_{core}\right)
+
+
+where $\nu_{outer} = \frac{\nu_{shell}t + \nu_{core}a}{a+t}$ is the aspect ratio of the outer ellipse.
+
+and the orientationally averaged intensity is
+
+.. math::
+
+    I(q) = \int_0^{\pi/2}\!\left|F_{cs}(q,\theta)\right|^2\sin\theta\,d\theta
+
+
+References
+----------
+#.  T. Kawaguchi, *J. Appl. Crystallogr*, 34(2001) 580-584
+#.  S. Förster, *J. Phys. Chem.*, 103(1999) 6657-6668
+
+Authorship and Verification
+---------------------------
+
+* **Author:** Itsuki Tajima and Yuhei Yamada (Github user name: 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"],
+    [
+        "radius_core",
+        "Ang",
+        5.0,
+        [0, inf],
+        "volume",
+        "Elliptical core minor radius a",
+    ],
+    ["thickness", "Ang", 2.0, [0, inf], "volume", "Shell thickness t"],
+    ["core_nu", "", 1.0, [0.1, 10.0], "volume", "Aspect ratio of core cross-section"],
+    ["shell_nu", "", 1.0, [0.1, 10.0], "volume", "Aspect ratio of shell cross-section"],
+    ["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"],
+]
+
+valid = "radius >= radius_core + thickness"
+
+# -- 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
+
+radius_effective_models = ["outer_radius", "equivalent_volume_sphere", "max_distance_from_center", "maijor_radius"]
+have_Fq = True
+
+tests = [
+    [{}, 1.047452236000633e-03, 2.312834927839701e00],
+    [{}, 2.099653600438444e-03, 2.287247365381240e00],
+    [{}, 4.208826990208894e-03, 2.186990768206136e00],
+]