Add tetrapod model implementation

sasmodels/models/tetrapod.c

@@ -0,0 +1,77 @@
+// Half of the tetrahedral angle acos(-1/3)/2 ~54.7356 deg.
+// Precompute its cosine and sine: cos = sqrt(1/3), sin = sqrt(2/3).
+#define COS_HALF_ANGLE 0.57735026918962573  // sqrt(1/3)
+#define SIN_HALF_ANGLE 0.81649658092772603  // sqrt(2/3)
+
+static double u_n(int n, double theta, double alpha) {
+  const double phi[4] = {0.0, M_PI_2, M_PI, 3.0 * M_PI_2};
+  const double sign[4] = {1.0, -1.0, 1.0, -1.0};
+  return sign[n] * COS_HALF_ANGLE * cos(theta) +
+         SIN_HALF_ANGLE * sin(theta) * cos(alpha - phi[n]);
+}
+
+// L: arm length, R: core radius, t: shell thickness, R_o = R + t: outer radius
+static double Fq_n(double q, double u, double L, double R, double t,
+                   double contrast_core, double contrast_shell) {
+  double R_o = R + t;
+  double quL2 = q * u * L * 0.5;
+  double mu = sqrt(fmax(0.0, 1.0 - u * u));
+  double V_c = M_PI * R * R * L;
+  double V_o = M_PI * R_o * R_o * L;
+  return sas_sinx_x(quL2) * (contrast_core * V_c * sas_2J1x_x(q * mu * R) +
+                             contrast_shell * V_o * sas_2J1x_x(q * mu * R_o));
+}
+
+static double form_volume(double L, double R, double t) {
+  // V = 4 * pi * (R + t)^2 * L  (outer volume of 4 arms)
+  return 4.0 * M_PI * (R + t) * (R + t) * L;
+}
+
+static double radius_effective(int mode, double L, double R, double t) {
+  switch (mode) {
+    default:
+    case 1:  // equivalent volume sphere
+      return cbrt(form_volume(L, R, t) / M_4PI_3);
+    case 2:  // length of tetrapod arms (L)
+      return L;
+  }
+}
+
+static double Iq(double q, double L, double R, double t, double sld_core,
+                 double sld_shell, double sld_solvent) {
+  double contrast_core = sld_core - sld_shell;
+  double contrast_shell = sld_shell - sld_solvent;
+  double total = 0.0;
+
+  for (int dtheta = 0; dtheta < GAUSS_N; dtheta++) {
+    double theta =
+        M_PI_2 * (GAUSS_Z[dtheta] + 1.0);  // map from [-1, 1] to [0, pi]
+    double w_theta =
+        GAUSS_W[dtheta] * M_PI_2;  // adjust weight for the new range
+
+    double integral_alpha = 0.0;
+    for (int dalpha = 0; dalpha < GAUSS_N; dalpha++) {
+      double alpha =
+          M_PI * (GAUSS_Z[dalpha] + 1.0);  // map from [-1, 1] to [0, 2*pi]
+      double w_alpha =
+          GAUSS_W[dalpha] * M_PI;  // adjust weight for the new range
+
+      double u[4], F[4];
+      for (int n = 0; n < 4; n++) {
+        u[n] = u_n(n, theta, alpha);
+        F[n] = Fq_n(q, u[n], L, R, t, contrast_core, contrast_shell);
+      }
+      double sum_arms = 0.0;
+      for (int n = 0; n < 4; n++) {
+        sum_arms += F[n] * F[n];
+        for (int m = n + 1; m < 4; m++) {
+          sum_arms += 2.0 * F[n] * F[m] * cos(q * (u[n] - u[m]) * L / 2.0);
+        }
+      }
+      sum_arms *= sin(theta);
+      integral_alpha += sum_arms * w_alpha;
+    }
+    total += integral_alpha * w_theta;
+  }
+  return 1e-4 * total / (4.0 * M_PI);
+}

sasmodels/models/tetrapod.py

@@ -0,0 +1,138 @@
+r"""
+Definition
+----------
+
+Calculates the scattering from a tetrapod-shaped structure. A tetrapod consists
+of four cylindrical arms radiating from a central point, oriented along the
+(1,1,1), (-1,-1,1), (-1,1,-1), and (1,-1,-1) directions.
+
+.. figure:: img/tetrapod.png
+
+    Tetrapod schematic.  Each of the four arms is a core--shell cylinder of
+    length $L$, core radius $R$ and shell thickness $t$;
+    the arm cross-section is shown in the inset.  The arms radiate from a
+    central junction along tetrahedral directions.
+
+The scattering intensity is calculated as an average over all
+orientations:
+
+.. math::
+
+    I(q) = \frac{(\Delta \rho)^2}{4\pi}
+        \int_0^{\pi} \int_0^{2\pi}
+        \left|\sum_{n=1}^{4} F_n(q, \theta, \varphi)\right|^2
+        \sin\theta \, d\theta \, d\varphi
+
+where $F_n$ is the core-shell form factor amplitude of the $n$-th arm:
+
+.. math::
+
+    F_n(q, \theta, \varphi) = \text{sinc}\!\left(\frac{q u_n L}{2}\right)
+        \left[
+            (\rho_\text{core} - \rho_\text{shell})\, V_\text{core}\,
+            \frac{2 J_1(q \mu_n R)}{q \mu_n R}
+          + (\rho_\text{shell} - \rho_\text{solvent})\, V_\text{outer}\,
+            \frac{2 J_1(q \mu_n (R+t))}{q \mu_n (R+t)}
+        \right]
+
+with $u_n = \hat{q} \cdot \hat{a}_n$, $\mu_n = \sqrt{1 - u_n^2}$, $L$ the arm
+length, $R$ the arm core radius, $R + t$ the outer radius ($t$ = shell
+thickness), and $V_\text{core} = \pi R^2 L$, $V_\text{outer} = \pi (R+t)^2 L$.
+Expanding the squared modulus into a double sum gives:
+
+.. math::
+
+    I(q) = \frac{1}{4\pi}
+        \int \sum_{n=1}^{4}\sum_{m=1}^{4}
+        F_n F_m \cos\!\left(\frac{q(u_n-u_m)L}{2}\right)
+        \sin\theta \, d\theta \, d\varphi
+
+The cosine factor is the interference term between the centres of arms $n$
+and $m$, which are displaced by $\tfrac{L}{2}\hat{a}_n$ from the junction.
+
+Geometry
+--------
+
+The four arms are oriented along tetrahedral directions.  With
+$A = 109.5 /2$ (the half-angle between arms), the arm unit vectors and the
+corresponding projections $u_n$ are
+
+.. math::
+
+    u_n = s_n \cos A \cos\theta + \sin A \sin\theta \cos(\varphi - \varphi_n)
+
+where $(s_n, \varphi_n) = (+1,\ 0),\ (-1,\ \pi/2),\ (+1,\ \pi),\ (-1,\ 3\pi/2)$
+for $n = 1, 2, 3, 4$ respectively.
+
+Each arm has length $L$, core radius $R$, shell thickness $t$, and hence outer
+radius $R + t$.
+
+.. note::
+
+    Each of the four arms is modelled as a complete cylinder of length $L$
+    extending from the central junction, so the arms overlap near the origin.
+    This overlap is neglected in two ways, following the treatment used in the
+    reference. First, the particle volume is overestimated (the effect being
+    largest when the arms are short and wide), which affects the calculated
+    intensity through the volume normalisation. Second, the scattering amplitude
+    is the plain sum of the four cylinder amplitudes, so the overlapping region
+    near the origin is counted more than once; because this region is small
+    compared with the arms, the resulting artefacts are expected to appear
+    mainly in the high-$q$ range. Consequently the model is only valid for
+    long-arm tetrapods, i.e. for arm lengths much larger than the arm width,
+    $L \gg R + t$.
+
+References
+----------
+
+#. Seoki Kyoo Seo *Korean J. Chem. Eng.* 34(2017) 1192-1198 DOI:10.1007/s11814-016-0341-x
+
+Authorship and Verification
+----------------------------
+
+* **Author:** Yuhei Yamada (Github user name: Indigo Carmine, https://orcid.org/0009-0003-9780-4135)
+* **Last Modified by:**
+"""
+
+from numpy import inf
+
+name = "tetrapod"
+title = "Core-shell tetrapod with four cylindrical arms"
+description = """
+    Calculates the scattering from a core-shell tetrapod structure with four
+    cylindrical arms radiating from a central point. Each arm has a cylindrical
+    core and a coaxial shell.
+"""
+category = "shape:cylinder"
+
+#             [ "name", "units", default, [lower, upper], "type", "description"],
+parameters = [
+    ["length", "Ang", 1000, [0, inf], "volume", "Arm length"],
+    ["radius", "Ang", 50, [0, inf], "volume", "Arm core radius"],
+    ["thickness", "Ang", 10, [0, inf], "volume", "Arm shell thickness"],
+    ["sld_core", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Arm core scattering length density"],
+    ["sld_shell", "1e-6/Ang^2", 0.5, [-inf, inf], "sld", "Arm shell scattering length density"],
+    ["sld_solvent", "1e-6/Ang^2", 0, [-inf, inf], "sld", "Solvent scattering length density"],
+]
+
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "tetrapod.c"]
+have_Fq = False
+opencl = True
+
+
+radius_effective_modes = [
+    "equivalent volume sphere",
+    "length of tetrapod arms (L)",
+]
+
+
+test = [
+    # thickness=0 reduces to uniform cylinder: contrast = sld_core - sld_solvent = 1
+    [{"thickness": 0}, 1.099643429303014e-03, 2.746377685546875e03],
+    [{"thickness": 0}, 1.033727919136565e-01, 4.374285638332367e-01],
+    [{"thickness": 0}, 3.145051218117226e-01, 3.527995198965073e-03],
+    # default core-shell parameters
+    [{}, 1.099643429303014e-03, 1.846798618097820e03],
+    [{}, 1.033727919136565e-01, 1.811197691988146e-01],
+    [{}, 3.145051218117226e-01, 1.080396557066592e-03],
+]