diff --git a/sasmodels/models/tetrapod.c b/sasmodels/models/tetrapod.c
new file mode 100644
index 00000000..b0189fd9
--- /dev/null
+++ b/sasmodels/models/tetrapod.c
@@ -0,0 +1,57 @@
+// half of the angle between arms in radians
+const double A = acos(-1 / 3) / 2.0 * M_PI / 180.0;
+
+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(A) * cos(theta) +
+         sin(A) * sin(theta) * cos(alpha - phi[n]);
+}
+
+static double Fq_n(double q, double u, double L, double R) {
+  double quL2 = q * u * L * 0.5;
+
+  double mu = sqrt(fmax(0.0, 1.0 - u * u));
+  double qmuR = q * mu * R;
+
+  return sas_sinx_x(quL2) * sas_2J1x_x(qmuR);
+}
+
+static double form_volume(double L, double R) {
+  // V = 4 * \pi * R^2 * L
+  return 4.0 * M_PI * R * R * L;
+}
+
+static double Iq(double q, double L, double R, double sld_particle,
+                 double sld_solvent) {
+  double contrast = sld_particle - 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 sum_arms = 0.0;
+      for (int n = 0; n < 4; n++) {
+        for (int m = 0; m < 4; m++) {
+          double u = u_n(n, theta, alpha);
+          double v = u_n(m, theta, alpha);
+          sum_arms += Fq_n(q, u, L, R) * Fq_n(q, v, L, R) *
+                      cos(q * (u - v) * L / 2.0) * sin(theta);
+        }
+      }
+      integral_alpha += sum_arms * w_alpha;
+    }
+    total += integral_alpha * w_theta;
+  }
+  return 1e-4 * contrast * contrast * total * 2 / M_PI;
+}
\ No newline at end of file
diff --git a/sasmodels/models/tetrapod.py b/sasmodels/models/tetrapod.py
new file mode 100644
index 00000000..eae3521f
--- /dev/null
+++ b/sasmodels/models/tetrapod.py
@@ -0,0 +1,90 @@
+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.
+
+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 $\Delta\rho$ is the SLD contrast and $F_n$ is the form factor amplitude
+of the $n$-th cylindrical arm:
+
+.. math::
+
+    F_n(q, \theta, \varphi) = \text{sinc}\!\left(\frac{q u_n L}{2}\right)
+        \cdot \frac{2 J_1(q \mu_n R)}{q \mu_n R}
+
+with $u_n = \hat{q} \cdot \hat{a}_n$ the projection of $\hat{q}$ onto the
+arm axis, $\mu_n = \sqrt{1 - u_n^2}$, $L$ the arm length, and $R$ the arm
+radius.  Expanding the squared modulus into a double sum and exploiting the
+reality of $F_n$ gives:
+
+.. math::
+
+    I(q) = \frac{(\Delta \rho)^2}{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$ and radius $R$.
+
+References
+----------
+
+#. Seoki Kyoo Seo *Korean J. Chem. Eng.* 34(2017) 1192-1198
+
+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 = "Tetrapod with four cylindrical arms"
+description = """
+    Calculates the scattering from a tetrapod structure with four cylindrical
+    arms radiating from a central point.
+"""
+category = "shape:cylinder"
+
+#             [ "name", "units", default, [lower, upper], "type", "description"],
+parameters = [
+    ["length", "Ang", 400, [0, inf], "volume", "Cylindrical arm length"],
+    ["radius", "Ang", 30, [0, inf], "volume", "Cylindrical arm radius"],
+    ["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Tetrapod scattering length density"],
+    ["sld_solvent", "1e-6/Ang^2", 1, [-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