Adding torus elliptical shell model

sasmodels/models/torus_elliptical_shell.c

@@ -0,0 +1,66 @@
+static double form_volume(double radius, double core_radius, double thickness,
+                          double nu) {
+  double ao = core_radius + thickness;
+  double bo = nu * ao;
+  double area = M_PI * ao * bo;
+  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 = 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 double Iq(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 = 0.0, theta = 0.0, 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);
+  }
+
+  // multiply by pi/4 to get the integral over [0, pi/2]
+  double result = I_total / form_volume(radius, core_radius, thickness, nu) *
+                  M_PI_4;  // A^-1
+  return 1e8 * result;     // convert from A^-1 to cm^-1
+}

sasmodels/models/torus_elliptical_shell.py

@@ -0,0 +1,115 @@
+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.
+
+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"],
+    [
+        "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"],
+]
+
+valid = "radius >= core_radius + 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
+
+
+tests = [
+    [{}, 1.002266990452620e-03, 2.390925896904364e07],
+    [{}, 2.009078240301904e-03, 2.366687319304617e07],
+    [{}, 1.017875911130553e-02, 1.684835184636866e07],
+]