From 96b1ed62ae6ea2354b86e149de22b1bb3cad45bb Mon Sep 17 00:00:00 2001
From: Paul Kienzle <pkienzle@nist.gov>
Date: Tue, 2 Jun 2026 12:27:31 -0400
Subject: [PATCH 1/5] truncated octahedron model copied from
 adding_trOh_model_purePython branch

---
 .../models/octahedron_truncated_fibonacci.py  | 492 ++++++++++++++++++
 1 file changed, 492 insertions(+)
 create mode 100644 explore/models/octahedron_truncated_fibonacci.py

diff --git a/explore/models/octahedron_truncated_fibonacci.py b/explore/models/octahedron_truncated_fibonacci.py
new file mode 100644
index 00000000..6647bcba
--- /dev/null
+++ b/explore/models/octahedron_truncated_fibonacci.py
@@ -0,0 +1,492 @@
+# octahedron_truncated model
+# Note: model title and parameter table are inserted automatically
+r"""
+This model provides the form factor P(q) for a general octahedron with orientational averaging (Fibonacci quadrature).
+It can be a regular octahedron shape with all edges of the same length.
+Or a general shape with different elongations along the three perpendicular two-fold axes.
+It includes the possibility to add an adjustable square truncation at each of the six vertices.
+This model includes the general cuboctahedron shape for the maximum value of truncation.
+The form factor expression is obtained by analytical integration over the volume of the shape.
+This model is constructed in a similar way as the rectangular prism model.
+It contains both the form factor for a reference orientation and the 1D form factor after orientation average
+using the Fibonacci quadrature. This quadrature provides a quasi-uniform distribution of points on the unit sphere
+using the golden ratio. The only new parameter is the number of points to generate on the unit sphere.
+The default value (around 400 points) provides a good balance between accuracy and computational efficiency.
+
+Definition
+----------
+
+The general octahedron is defined by its dimensions along its three perpendicular two-fold axes along x, y and z directions.
+length_a, length_b and length_c are the distances from the center of the general octahedron to its 6 vertices.
+
+Coordinates of the six vertices are:
+    (length_a, 0, 0),
+    (-length_a, 0, 0),
+    (0, length_b, 0),
+    (0, -length_b, 0),
+    (0, 0, length_c),
+    (0, 0, -length_c)
+
+t is the truncation parameter.
+Truncation adds a square facet for each vertex that is perpendicular to a 2-fold axis.
+The resulting shape consists of six squares and eight hexagons, which may be irregular depending on the three dimensions
+A square facet crosses the x, y, z directions at distances equal to t length_a, t length_b and t length_c.
+
+A regular octahedron corresponds to:
+
+.. math::
+
+    length_a = length_b = length_c, \quad t = 1
+
+A regular cuboctahedron shape with 6 squares and 8 triangles corresponds to:
+
+.. math::
+
+    length_a = length_b = length_c, \quad t = \frac{1}{2}
+
+The model contains 4 parameters: length_a, the two ratios b2a_ratio and c2a_ratio and t:
+
+.. math::
+
+    b2a_{\text{ratio}} = \frac{length_b}{length_a}, \quad
+    c2a_{\text{ratio}} = \frac{length_c}{length_a}, \quad
+    \frac{1}{2} < t < 1
+
+For a regular shape:
+
+.. math::
+
+    b2a_{\text{ratio}} = c2a_{\text{ratio}} = 1
+
+Volume of the general shape including truncation is given by:
+
+.. math::
+
+    V = \frac{4}{3}\, length_{\text{a}}^{3}\, b2a_{\text{ratio}}\, c2a_{\text{ratio}}\,\bigl(1 - 3(1 - t)^{3}\bigr)
+
+The general octahedron is made of eight triangular faces. The three edge lengths
+are:
+
+.. math::
+
+    A_{\text{edge}}^{2} = length_{\text{a}}^{2} + length_{\text{b}}^{2},\qquad
+    B_{\text{edge}}^{2} = length_{\text{a}}^{2} + length_{\text{c}}^{2},\qquad
+    C_{\text{edge}}^{2} = length_{\text{b}}^{2} + length_{\text{c}}^{2}
+
+For a regular shape (no elongation):
+
+.. math::
+
+    b2a_{\text{ratio}} = c2a_{\text{ratio}} = 1,\qquad
+    A_{\text{edge}} = B_{\text{edge}} = C_{\text{edge}} = length_{\text{a}} \sqrt{2},\qquad
+    length_{\text{a}} = length_{\text{b}} = length_{\text{c}} = A_{\text{edge}} / \sqrt{2}
+
+.. math::
+
+    V = \frac{4}{3} \, length_{\text{a}}^{3} \, \bigl(1 - 3 (1 - t)^3 \bigr)
+
+The reference orientation of the shape is: a along x, b along y and c along z.
+Amplitude of the form factor AP for the reference orientation of the shape reads
+
+.. math::
+
+    AP(q,\theta,\phi) = \frac{3}{1 - 3(1 - t)^3}\,(AA + BB + CC)
+
+.. math::
+
+    AA = \frac{1}{2\,(q_y^2 - q_z^2)\,(q_y^2 - q_x^2)}\Big[(q_y - q_x)\sin\big(q_y(1 - t) - q_x t\big)
+    + (q_y + q_x)\sin\big(q_y(1 - t) + q_x t\big)\Big]
+    + \frac{1}{2\,(q_z^2 - q_x^2)\,(q_z^2 - q_y^2)}\Big[(q_z - q_x)\sin\big(q_z(1 - t) - q_x t\big)
+    + (q_z + q_x)\sin\big(q_z(1 - t) + q_x t\big)\Big]
+
+.. math::
+
+    BB = \frac{1}{2\,(q_z^2 - q_x^2)\,(q_z^2 - q_y^2)}\Big[(q_z - q_y)\sin\big(q_z(1 - t) - q_y t\big)
+    + (q_z + q_y)\sin\big(q_z(1 - t) + q_y t\big)\Big]
+    + \frac{1}{2\,(q_x^2 - q_y^2)\,(q_x^2 - q_z^2)}\Big[(q_x - q_y)\sin\big(q_x(1 - t) - q_y t\big)
+    + (q_x + q_y)\sin\big(q_x(1 - t) + q_y t\big)\Big]
+
+.. math::
+
+    CC = \frac{1}{2\,(q_x^2 - q_y^2)\,(q_x^2 - q_z^2)}\Big[(q_x - q_z)\sin\big(q_x(1 - t) - q_z t\big)
+    + (q_x + q_z)\sin\big(q_x(1 - t) + q_z t\big)\Big]
+    + \frac{1}{2\,(q_y^2 - q_z^2)\,(q_y^2 - q_x^2)}\Big[(q_y - q_z)\sin\big(q_y(1 - t) - q_z t\big)
+    + (q_y + q_z)\sin\big(q_y(1 - t) + q_z t\big)\Big]
+
+Capital Qx Qy Qz are the three components in [A-1] of the scattering vector.
+qx qy qz are rescaled components (no unit) for computing AA, BB and CC terms.
+
+.. math::
+
+    Q_x = q\,\sin\theta\,\cos\phi, \qquad
+    Q_y = q\,\sin\theta\,\sin\phi, \qquad
+    Q_z = q\,\cos\theta
+
+.. math::
+
+    q_x = Q_x \, length_{\text{a}},\qquad
+    q_y = Q_y \, length_{\text{b}},\qquad
+    q_z = Q_z \, length_{\text{c}}
+
+
+θ is the angle between the scattering vector and the z axis.
+ϕ is the rotation angle in the xy plane.
+
+The octahedron is in its reference orientation, with the c-axis aligned along z and the a-axis aligned along x.
+
+The 1D form factor P(q) corresponds to the orientation average with all the possible orientations having the same probability.
+Instead of rotating the shape through all the possible orientations in the integral,
+it is equivalent to integrate the 3D scattering vector over a sphere of radius q with the shape in its reference orientation.
+
+.. math::
+
+    P(q) =  \frac{2}{\pi} \int_0^{\frac{\pi}{2}} \,
+    \int_0^{\frac{\pi}{2}} A_P^2(q) \, \sin\theta \, d\theta \, d\phi
+
+The sphere is sampled using Fibonacci quadrature to provide a quasi-uniform distribution of points on the unit sphere.
+The repartition of the points is computed using the golden ratio (see fibonacci.py).
+
+
+.. figure:: img/fibonacci_sphere.png
+
+    Fibonacci sphere using 5810 points.
+
+And the 1D scattering intensity is calculated as
+
+.. math::
+
+    I(q) = \text{scale} \times V \times (\rho_\text{p} -
+    \rho_\text{solvent})^2 \times P(q)
+
+where V is the volume of the truncated octahedron, ρ
+is the scattering length inside the volume, ρ *solvent*
+is the scattering length of the solvent, and (if the data are in absolute
+units) *scale* represents the volume fraction (which is unitless).
+
+
+.. figure:: img/octa-truncated.png
+
+    Truncated octahedron shape for different truncation.
+
+.. figure:: img/octahedrons_intensity_plot.png
+
+    Scattering intensity of a cuboctahedron (t=0.5) and a regular octahedron (t=1) of a = 300 Angstroms.
+
+Validation
+----------
+
+Validation of the code is made using numerical checks.
+Comparisons with Debye formula calculations were made using DebyeCalculator library (https://github.com/FrederikLizakJohansen/DebyeCalculator).
+Good agreement was found at q < 0.1 1/Angstrom.
+
+References
+----------
+
+1. Wei-Ren Chen et al. "Scattering functions of Platonic solids".
+   In: Journal of Applied Crystallography - J APPL CRYST 44 (June 2011).
+   DOI: 10.1107/S0021889811011691
+
+2. Croset, Bernard, "Form factor of any polyhedron: a general compact
+   formula and its singularities" In: J. Appl. Cryst. (2017). 50, 1245–1255
+   https://doi.org/10.1107/S1600576717010147
+
+3. Wuttke, J. Numerically stable form factor of any polygon and polyhedron
+   J Appl Cryst 54, 580-587 (2021)
+   https://doi.org/10.1107/S160057672100171
+
+4. Álvaro González. Measurement of areas on a sphere using Fibonacci and latitude–longitude lattices.
+Mathematical geosciences 42 (2010), pp. 49–64
+
+5. MB Kozin, VV Volkov, and DI Svergun. ASSA, a program for three-dimensional rendering in solution
+scattering from biopolymers. Applied Crystallography 30.5 (1997), pp. 811–815
+
+Authorship and Verification
+----------------------------
+
+* **Authors:** Marianne Imperor-Clerc (marianne.imperor@cnrs.fr)
+             Sara Mokhtari (smokhtari@insa-toulouse.fr)
+
+* **Last Modified by:** SM **Date:** 13th January 2026
+
+* **Last Reviewed by:** MIC **Date:** 13th January 2026
+
+"""
+
+import numpy as np
+from numpy import inf
+
+from sasmodels.quadratures.fibonacci import fibonacci_sphere
+
+name = "octahedron_truncated_fibonacci"
+title = "General octahedron pure python model using Fibonacci quadrature"
+description = """
+The amplitude AP is defined as follows.
+
+AP = 6./(1.-3*(1.-t)*(1.-t)*(1.-t))*(AA+BB+CC)
+
+AA = 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx)) * ((qy-qx)*sin(qy*(1.-t)-qx*t) + (qy+qx)*sin(qy*(1.-t)+qx*t)) + 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy)) * ((qz-qx)*sin(qz*(1.-t)-qx*t) + (qz+qx)*sin(qz*(1.-t)+qx*t))
+
+BB = 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy)) * ((qz-qy)*sin(qz*(1.-t)-qy*t) + (qz+qy)*sin(qz*(1.-t)+qy*t)) + 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz)) * ((qx-qy)*sin(qx*(1.-t)-qy*t) + (qx+qy)*sin(qx*(1.-t)+qy*t))
+
+CC = 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz)) * ((qx-qz)*sin(qx*(1.-t)-qz*t) + (qx+qz)*sin(qx*(1.-t)+qz*t)) + 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx)) * ((qy-qz)*sin(qy*(1.-t)-qz*t) + (qy+qz)*sin(qy*(1.-t)+qz*t))
+
+With capital QX QY QZ are the three components in [A-1] of the scattering vector
+and qx qy qz are the rescaled components (no unit) for computing AP term.
+
+Qx = q * sin_theta * cos_phi
+Qy = q * sin_theta * sin_phi
+Qz = q * cos_theta
+qx = Qx * length_a
+qy = Qy * length_b
+qz = Qz * length_c
+
+Reference orientation is with a along x axis, b along y axis and c along z axis
+
+Valid truncation parameter range: 0.5 < t < 1.
+
+t=1 is for octahedron
+t=0.5 is for cuboctahedron
+"""
+category = "shape:polyhedron"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["sld", "1e-6/Ang^2", 130, [-inf, inf], "sld",
+               "Octahedron scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2", 9.4, [-inf, inf], "sld",
+               "Solvent scattering length density"],
+              ["length_a", "Ang", 500, [0, inf], "volume",
+               "half height along a axis"],
+              ["b2a_ratio", "", 1, [0, inf], "volume",
+               "Ratio b/a"],
+              ["c2a_ratio", "", 1, [0, inf], "volume",
+               "Ratio c/a"],
+              ["t", "", 0.99, [0.0, 1.0], "volume",
+               "truncation along axis"],
+              ["npoints_fibonacci",     "",           400, [1, 1e6],   "",
+                       "Number of points on the sphere for the Fibonacci integration"]
+              ]
+
+
+def form_volume(length_a, b2a_ratio, c2a_ratio, t):
+    """
+    Calculate the volume of the truncated octahedron.
+
+    Parameters
+    ----------
+    length_a: Half height along the a axis of the octahedron without truncature
+    b2a_ratio: Half height along the b axis of the octahedron without truncature
+    c2a_ratio: Half height along the c axis of the octahedron without truncature
+    t: Truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
+
+    Returns
+    -------
+    volume: Volume of the truncated octahedron
+    """
+    return (
+        (4 / 3)
+        * length_a**3
+        * b2a_ratio
+        * c2a_ratio
+        * (1.0 - 3 * (1.0 - t) ** 3)
+    )
+
+
+
+def A(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
+    """
+    Computes the AA term of the amplitude of the form factor AP at (qa,qb,qc).
+
+    qnx, qny, qnz are the rescaled components (no unit) for computing AA term.
+
+    Parameters
+    ----------
+    qa: the component of the scattering vector along a axis, units 1/Ang
+    qb: the component of the scattering vector along b axis, units 1/Ang
+    qc: the component of the scattering vector along c axis, units 1/Ang
+    length_a: half height along the a axis of the octahedron without truncature
+    b2a_ratio: half height along the b axis of the octahedron without truncature
+    c2a_ratio: half height along the c axis of the octahedron without truncature
+    t: truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
+
+    Returns
+    -------
+    aa: AA term of the amplitude of the form factor at (qa,qb,qc)
+    """
+    length_b = length_a * b2a_ratio
+    length_c = length_a * c2a_ratio
+    qnx = qa * length_a  # conversion to dimensionless coordinate
+    qny = qb * length_b  # conversion to dimensionless coordinate
+    qnz = qc * length_c  # conversion to dimensionless coordinate
+
+    # Protect denominators against exact zeros by adding tiny epsilon
+    eps_den = 1e-18
+    denom1 = (qny * qny - qnz * qnz) * (qny * qny - qnx * qnx)
+    denom2 = (qnz * qnz - qnx * qnx) * (qnz * qnz - qny * qny)
+    denom1 += eps_den * (np.abs(denom1) < eps_den)  # add epsilon where denom is zero
+    denom2 += eps_den * (np.abs(denom2) < eps_den)  # add epsilon where denom is zero
+
+    term1 = (
+        (qny - qnx) * np.sin(qny * (1.0 - t) - qnx * t)
+        + (qny + qnx) * np.sin(qny * (1.0 - t) + qnx * t)
+    ) / denom1
+    term2 = (
+        (qnz - qnx) * np.sin(qnz * (1.0 - t) - qnx * t)
+        + (qnz + qnx) * np.sin(qnz * (1.0 - t) + qnx * t)
+    ) / denom2
+
+    return term1 + term2
+
+def Fqabc(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
+    """
+    Computes the amplitude of the form factor AP at (qa,qb,qc).
+
+    Parameters
+    ----------
+    qa: component of the scattering vector along a axis, units 1/Ang
+    qb: component of the scattering vector along b axis, units 1/Ang
+    qc: component of the scattering vector along c axis, units 1/Ang
+    length_a: half height along the a axis of the octahedron without truncature
+    b2a_ratio: half height along the b axis of the octahedron without truncature
+    c2a_ratio: half height along the c axis of the octahedron without truncature
+    t: truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
+
+    Returns
+    -------
+    ap: Amplitude of the form factor at (qa,qb,qc)
+    """
+    # Code taking into account the circular permutation between AA, BB and CC terms in AP.
+    # amp3D is scaled to 1 near the origin for a regular shape.
+    AA = A(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t)
+    BB = A(qb, qc, qa, length_a, b2a_ratio, c2a_ratio, t)
+    CC = A(qc, qa, qb, length_a, b2a_ratio, c2a_ratio, t)
+
+    # Normalisation to 1 of AP at q = 0. Division by a factor 4/3 and global 1/2 coefficient.
+    ap = 3.0 / (1.0 - 3.0 * (1.0 - t) ** 3) * (AA + BB + CC)
+
+    # Guard against numerical issues: mirror the C fallback which sets
+    # invalid (NaN/inf) values to zero so tests and downstream code
+    # do not receive NaNs. Use np.nan_to_num for vector/scalar safety.
+    return np.nan_to_num(ap, nan=0.0, posinf=0.0, neginf=0.0)
+
+def Iqabc(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
+    """
+    Computes the 3D scattering intensity.
+
+    Parameters
+    ----------
+    qa: component of the scattering vector along a axis, units 1/Ang
+    qb: component of the scattering vector along b axis, units 1/Ang
+    qc: component of the scattering vector along c axis, units 1/Ang
+    length_a: half height along the a axis of the octahedron without truncature
+    b2a_ratio: half height along the b axis of the octahedron without truncature
+    c2a_ratio: half height along the c axis of the octahedron without truncature
+    t: truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
+
+    Returns
+    -------
+    intensity: 3D scattering intensity at (qa,qb,qc), units 1/cm
+    """
+    # int3D is scaled to 1 near the origin for a regular shape
+    amp = Fqabc(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t)
+    intensity = amp**2
+
+    # If numerical issues produced invalid numbers, return 0.0 (like the C
+    # implementation does) rather than NaN/inf which break tests and usage.
+    if np.isnan(intensity) or np.isinf(intensity):
+        return 0.0
+
+    return intensity
+
+
+# Cache for Fibonacci sphere points/weights to avoid recomputing across calls
+_fibonacci_sphere_cache = {}
+
+def _get_fibonacci_sphere(npoints):
+    """
+    Return (q_unit, w) for given npoints using module cache.
+    """
+    key = int(npoints)
+    if key in _fibonacci_sphere_cache:
+        return _fibonacci_sphere_cache[key]
+
+    q_unit, w = fibonacci_sphere(int(npoints))
+    _fibonacci_sphere_cache[key] = (q_unit, w)
+    return q_unit, w
+
+
+def Iq(
+    q,
+    sld,
+    sld_solvent,
+    length_a=500,
+    b2a_ratio=1,
+    c2a_ratio=1,
+    t=0.99,
+    npoints_fibonacci: int = 400,
+):
+    """
+    Computes the 1D scattering intensity I(q) using Fibonacci quadrature.
+
+    Parameters
+    ----------
+    q : float or ndarray
+        Magnitude of the scattering vector, units 1/Ang
+    sld : float
+        Octahedron scattering length density, units 1e-6/Ang^2
+    sld_solvent : float
+        Solvent scattering length density, units 1e-6/Ang^2
+    length_a : float
+        Half height along the a axis of the octahedron without truncature, units Angstrom
+    b2a_ratio : float
+        Ratio b/a
+    c2a_ratio : float
+        Ratio c/a
+    t : float
+        Truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
+    npoints_fibonacci : int
+        Number of points on the sphere for the Fibonacci integration
+
+    Returns
+    -------
+    1D scattering intensity at q, units 1/cm
+    """
+    q = np.atleast_1d(q)
+    q_unit, w = _get_fibonacci_sphere(npoints_fibonacci)  # shape (npoints, 3)
+
+    # build qx,qy,qz arrays with correct broadcasting -> shape (nq, npoints)
+    qx = q[:, None] * q_unit[None, :, 0]
+    qy = q[:, None] * q_unit[None, :, 1]
+    qz = q[:, None] * q_unit[None, :, 2]
+
+    # compute intensity grid using existing amp3D (vectorized)
+    amp_grid = Fqabc(
+        qx, qy, qz, length_a, b2a_ratio, c2a_ratio, t
+    )  # shape (nq, npoints)
+
+    # Replace any NaN/inf in amplitude grid (temporary fallback mirroring C
+    # implementation). Use 0.0 to avoid producing huge/NaN intensities.
+    amp_grid = np.nan_to_num(amp_grid, nan=0.0, posinf=0.0, neginf=0.0)
+    # Also clip extremely large values to avoid overflow in square
+    amp_grid = np.clip(amp_grid, -1e12, 1e12)
+
+    intensity_grid = np.abs(amp_grid) ** 2
+
+    integral = np.sum(intensity_grid * w[None, :], axis=1)  # summation over all points
+
+    # Convert from [1e-12 A-1] to [cm-1]
+    return (
+        integral
+        * 0.0001
+        * (sld - sld_solvent) ** 2
+        * form_volume(length_a, b2a_ratio, c2a_ratio, t) ** 2
+    )
+
+
+Iq.vectorized = True
+
+tests = [
+    [{"background": 0, "scale": 1, "length_a": 100, "t": 1, "sld": 1., "sld_solvent": 0.},
+     0.01, 120.57219],
+    [{"background": 0, "scale": 1, "length_a": 100, "t": 1, "sld": 1., "sld_solvent": 0.},
+     [0.01, 0.1], [120.57219, 0.37414]],
+]

From 4bcbbb4e2b6f7431640140679f624284c8199981 Mon Sep 17 00:00:00 2001
From: Paul Kienzle <pkienzle@nist.gov>
Date: Tue, 2 Jun 2026 12:34:07 -0400
Subject: [PATCH 2/5] rename octahedron_truncated to truncated_octahedron

---
 ...cated_fibonacci.py => truncated_octahedron_fibonacci.py} | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)
 rename explore/models/{octahedron_truncated_fibonacci.py => truncated_octahedron_fibonacci.py} (99%)

diff --git a/explore/models/octahedron_truncated_fibonacci.py b/explore/models/truncated_octahedron_fibonacci.py
similarity index 99%
rename from explore/models/octahedron_truncated_fibonacci.py
rename to explore/models/truncated_octahedron_fibonacci.py
index 6647bcba..127efe71 100644
--- a/explore/models/octahedron_truncated_fibonacci.py
+++ b/explore/models/truncated_octahedron_fibonacci.py
@@ -1,4 +1,4 @@
-# octahedron_truncated model
+# truncated octahedron model
 # Note: model title and parameter table are inserted automatically
 r"""
 This model provides the form factor P(q) for a general octahedron with orientational averaging (Fibonacci quadrature).
@@ -215,9 +215,9 @@
 import numpy as np
 from numpy import inf
 
-from sasmodels.quadratures.fibonacci import fibonacci_sphere
+from sasmodels.special.fibonacci import fibonacci_sphere
 
-name = "octahedron_truncated_fibonacci"
+name = "truncated_octahedron_fibonacci"
 title = "General octahedron pure python model using Fibonacci quadrature"
 description = """
 The amplitude AP is defined as follows.

From 168b1dfc2b3ee832fe6f9687fad7907a72f18297 Mon Sep 17 00:00:00 2001
From: Paul Kienzle <pkienzle@nist.gov>
Date: Tue, 2 Jun 2026 12:41:52 -0400
Subject: [PATCH 3/5] rename length_a to radius_a and t to truncation

---
 .../models/truncated_octahedron_fibonacci.py  | 108 +++++++++---------
 1 file changed, 54 insertions(+), 54 deletions(-)

diff --git a/explore/models/truncated_octahedron_fibonacci.py b/explore/models/truncated_octahedron_fibonacci.py
index 127efe71..8e20714f 100644
--- a/explore/models/truncated_octahedron_fibonacci.py
+++ b/explore/models/truncated_octahedron_fibonacci.py
@@ -17,39 +17,39 @@
 ----------
 
 The general octahedron is defined by its dimensions along its three perpendicular two-fold axes along x, y and z directions.
-length_a, length_b and length_c are the distances from the center of the general octahedron to its 6 vertices.
+radius_a, radius_b and radius_c are the distances from the center of the general octahedron to its 6 vertices.
 
 Coordinates of the six vertices are:
-    (length_a, 0, 0),
-    (-length_a, 0, 0),
-    (0, length_b, 0),
-    (0, -length_b, 0),
-    (0, 0, length_c),
-    (0, 0, -length_c)
+    (radius_a, 0, 0),
+    (-radius_a, 0, 0),
+    (0, radius_b, 0),
+    (0, -radius_b, 0),
+    (0, 0, radius_c),
+    (0, 0, -radius_c)
 
 t is the truncation parameter.
 Truncation adds a square facet for each vertex that is perpendicular to a 2-fold axis.
 The resulting shape consists of six squares and eight hexagons, which may be irregular depending on the three dimensions
-A square facet crosses the x, y, z directions at distances equal to t length_a, t length_b and t length_c.
+A square facet crosses the x, y, z directions at distances equal to truncation radius_a, truncation radius_b and truncation radius_c.
 
 A regular octahedron corresponds to:
 
 .. math::
 
-    length_a = length_b = length_c, \quad t = 1
+    radius_a = radius_b = radius_c, \quad t = 1
 
 A regular cuboctahedron shape with 6 squares and 8 triangles corresponds to:
 
 .. math::
 
-    length_a = length_b = length_c, \quad t = \frac{1}{2}
+    radius_a = radius_b = radius_c, \quad truncation = \frac{1}{2}
 
-The model contains 4 parameters: length_a, the two ratios b2a_ratio and c2a_ratio and t:
+The model contains 4 parameters: radius_a, the two ratios b2a_ratio and c2a_ratio and truncation:
 
 .. math::
 
-    b2a_{\text{ratio}} = \frac{length_b}{length_a}, \quad
-    c2a_{\text{ratio}} = \frac{length_c}{length_a}, \quad
+    b2a_{\text{ratio}} = \frac{radius_b}{radius_a}, \quad
+    c2a_{\text{ratio}} = \frac{radius_c}{radius_a}, \quad
     \frac{1}{2} < t < 1
 
 For a regular shape:
@@ -236,9 +236,9 @@
 Qx = q * sin_theta * cos_phi
 Qy = q * sin_theta * sin_phi
 Qz = q * cos_theta
-qx = Qx * length_a
-qy = Qy * length_b
-qz = Qz * length_c
+qx = Qx * radius_a
+qy = Qy * radius_b
+qz = Qz * radius_c
 
 Reference orientation is with a along x axis, b along y axis and c along z axis
 
@@ -254,29 +254,29 @@
                "Octahedron scattering length density"],
               ["sld_solvent", "1e-6/Ang^2", 9.4, [-inf, inf], "sld",
                "Solvent scattering length density"],
-              ["length_a", "Ang", 500, [0, inf], "volume",
+              ["radius_a", "Ang", 500, [0, inf], "volume",
                "half height along a axis"],
               ["b2a_ratio", "", 1, [0, inf], "volume",
                "Ratio b/a"],
               ["c2a_ratio", "", 1, [0, inf], "volume",
                "Ratio c/a"],
-              ["t", "", 0.99, [0.0, 1.0], "volume",
+              ["truncation", "", 0.99, [0.0, 1.0], "volume",
                "truncation along axis"],
               ["npoints_fibonacci",     "",           400, [1, 1e6],   "",
                        "Number of points on the sphere for the Fibonacci integration"]
               ]
 
 
-def form_volume(length_a, b2a_ratio, c2a_ratio, t):
+def form_volume(radius_a, b2a_ratio, c2a_ratio, truncation):
     """
     Calculate the volume of the truncated octahedron.
 
     Parameters
     ----------
-    length_a: Half height along the a axis of the octahedron without truncature
+    radius_a: Half height along the a axis of the octahedron without truncature
     b2a_ratio: Half height along the b axis of the octahedron without truncature
     c2a_ratio: Half height along the c axis of the octahedron without truncature
-    t: Truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
+    truncation: Truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
 
     Returns
     -------
@@ -284,15 +284,15 @@ def form_volume(length_a, b2a_ratio, c2a_ratio, t):
     """
     return (
         (4 / 3)
-        * length_a**3
+        * radius_a**3
         * b2a_ratio
         * c2a_ratio
-        * (1.0 - 3 * (1.0 - t) ** 3)
+        * (1.0 - 3 * (1.0 - truncation) ** 3)
     )
 
 
 
-def A(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
+def A(qa, qb, qc, radius_a, b2a_ratio, c2a_ratio, truncation):
     """
     Computes the AA term of the amplitude of the form factor AP at (qa,qb,qc).
 
@@ -303,20 +303,20 @@ def A(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
     qa: the component of the scattering vector along a axis, units 1/Ang
     qb: the component of the scattering vector along b axis, units 1/Ang
     qc: the component of the scattering vector along c axis, units 1/Ang
-    length_a: half height along the a axis of the octahedron without truncature
+    radius_a: half height along the a axis of the octahedron without truncature
     b2a_ratio: half height along the b axis of the octahedron without truncature
     c2a_ratio: half height along the c axis of the octahedron without truncature
-    t: truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
+    truncation: truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
 
     Returns
     -------
     aa: AA term of the amplitude of the form factor at (qa,qb,qc)
     """
-    length_b = length_a * b2a_ratio
-    length_c = length_a * c2a_ratio
-    qnx = qa * length_a  # conversion to dimensionless coordinate
-    qny = qb * length_b  # conversion to dimensionless coordinate
-    qnz = qc * length_c  # conversion to dimensionless coordinate
+    radius_b = radius_a * b2a_ratio
+    radius_c = radius_a * c2a_ratio
+    qnx = qa * radius_a  # conversion to dimensionless coordinate
+    qny = qb * radius_b  # conversion to dimensionless coordinate
+    qnz = qc * radius_c  # conversion to dimensionless coordinate
 
     # Protect denominators against exact zeros by adding tiny epsilon
     eps_den = 1e-18
@@ -326,17 +326,17 @@ def A(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
     denom2 += eps_den * (np.abs(denom2) < eps_den)  # add epsilon where denom is zero
 
     term1 = (
-        (qny - qnx) * np.sin(qny * (1.0 - t) - qnx * t)
-        + (qny + qnx) * np.sin(qny * (1.0 - t) + qnx * t)
+        (qny - qnx) * np.sin(qny * (1.0 - truncation) - qnx * truncation)
+        + (qny + qnx) * np.sin(qny * (1.0 - truncation) + qnx * truncation)
     ) / denom1
     term2 = (
-        (qnz - qnx) * np.sin(qnz * (1.0 - t) - qnx * t)
-        + (qnz + qnx) * np.sin(qnz * (1.0 - t) + qnx * t)
+        (qnz - qnx) * np.sin(qnz * (1.0 - truncation) - qnx * truncation)
+        + (qnz + qnx) * np.sin(qnz * (1.0 - truncation) + qnx * truncation)
     ) / denom2
 
     return term1 + term2
 
-def Fqabc(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
+def Fqabc(qa, qb, qc, radius_a, b2a_ratio, c2a_ratio, truncation):
     """
     Computes the amplitude of the form factor AP at (qa,qb,qc).
 
@@ -345,10 +345,10 @@ def Fqabc(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
     qa: component of the scattering vector along a axis, units 1/Ang
     qb: component of the scattering vector along b axis, units 1/Ang
     qc: component of the scattering vector along c axis, units 1/Ang
-    length_a: half height along the a axis of the octahedron without truncature
+    radius_a: half height along the a axis of the octahedron without truncature
     b2a_ratio: half height along the b axis of the octahedron without truncature
     c2a_ratio: half height along the c axis of the octahedron without truncature
-    t: truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
+    truncation: truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
 
     Returns
     -------
@@ -356,19 +356,19 @@ def Fqabc(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
     """
     # Code taking into account the circular permutation between AA, BB and CC terms in AP.
     # amp3D is scaled to 1 near the origin for a regular shape.
-    AA = A(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t)
-    BB = A(qb, qc, qa, length_a, b2a_ratio, c2a_ratio, t)
-    CC = A(qc, qa, qb, length_a, b2a_ratio, c2a_ratio, t)
+    AA = A(qa, qb, qc, radius_a, b2a_ratio, c2a_ratio, truncation)
+    BB = A(qb, qc, qa, radius_a, b2a_ratio, c2a_ratio, truncation)
+    CC = A(qc, qa, qb, radius_a, b2a_ratio, c2a_ratio, truncation)
 
     # Normalisation to 1 of AP at q = 0. Division by a factor 4/3 and global 1/2 coefficient.
-    ap = 3.0 / (1.0 - 3.0 * (1.0 - t) ** 3) * (AA + BB + CC)
+    ap = 3.0 / (1.0 - 3.0 * (1.0 - truncation) ** 3) * (AA + BB + CC)
 
     # Guard against numerical issues: mirror the C fallback which sets
     # invalid (NaN/inf) values to zero so tests and downstream code
     # do not receive NaNs. Use np.nan_to_num for vector/scalar safety.
     return np.nan_to_num(ap, nan=0.0, posinf=0.0, neginf=0.0)
 
-def Iqabc(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
+def Iqabc(qa, qb, qc, radius_a, b2a_ratio, c2a_ratio, truncation):
     """
     Computes the 3D scattering intensity.
 
@@ -377,17 +377,17 @@ def Iqabc(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t):
     qa: component of the scattering vector along a axis, units 1/Ang
     qb: component of the scattering vector along b axis, units 1/Ang
     qc: component of the scattering vector along c axis, units 1/Ang
-    length_a: half height along the a axis of the octahedron without truncature
+    radius_a: half height along the a axis of the octahedron without truncature
     b2a_ratio: half height along the b axis of the octahedron without truncature
     c2a_ratio: half height along the c axis of the octahedron without truncature
-    t: truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
+    truncation: truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
 
     Returns
     -------
     intensity: 3D scattering intensity at (qa,qb,qc), units 1/cm
     """
     # int3D is scaled to 1 near the origin for a regular shape
-    amp = Fqabc(qa, qb, qc, length_a, b2a_ratio, c2a_ratio, t)
+    amp = Fqabc(qa, qb, qc, radius_a, b2a_ratio, c2a_ratio, truncation)
     intensity = amp**2
 
     # If numerical issues produced invalid numbers, return 0.0 (like the C
@@ -418,10 +418,10 @@ def Iq(
     q,
     sld,
     sld_solvent,
-    length_a=500,
+    radius_a=500,
     b2a_ratio=1,
     c2a_ratio=1,
-    t=0.99,
+    truncation=0.99,
     npoints_fibonacci: int = 400,
 ):
     """
@@ -435,13 +435,13 @@ def Iq(
         Octahedron scattering length density, units 1e-6/Ang^2
     sld_solvent : float
         Solvent scattering length density, units 1e-6/Ang^2
-    length_a : float
+    radius_a : float
         Half height along the a axis of the octahedron without truncature, units Angstrom
     b2a_ratio : float
         Ratio b/a
     c2a_ratio : float
         Ratio c/a
-    t : float
+    truncation : float
         Truncation parameter, varies from 0.5 (cuboctahedron) to 1 (octahedron)
     npoints_fibonacci : int
         Number of points on the sphere for the Fibonacci integration
@@ -460,7 +460,7 @@ def Iq(
 
     # compute intensity grid using existing amp3D (vectorized)
     amp_grid = Fqabc(
-        qx, qy, qz, length_a, b2a_ratio, c2a_ratio, t
+        qx, qy, qz, radius_a, b2a_ratio, c2a_ratio, truncation
     )  # shape (nq, npoints)
 
     # Replace any NaN/inf in amplitude grid (temporary fallback mirroring C
@@ -478,15 +478,15 @@ def Iq(
         integral
         * 0.0001
         * (sld - sld_solvent) ** 2
-        * form_volume(length_a, b2a_ratio, c2a_ratio, t) ** 2
+        * form_volume(radius_a, b2a_ratio, c2a_ratio, truncation) ** 2
     )
 
 
 Iq.vectorized = True
 
 tests = [
-    [{"background": 0, "scale": 1, "length_a": 100, "t": 1, "sld": 1., "sld_solvent": 0.},
+    [{"background": 0, "scale": 1, "radius_a": 100, "truncation": 1, "sld": 1., "sld_solvent": 0.},
      0.01, 120.57219],
-    [{"background": 0, "scale": 1, "length_a": 100, "t": 1, "sld": 1., "sld_solvent": 0.},
+    [{"background": 0, "scale": 1, "radius_a": 100, "truncation": 1, "sld": 1., "sld_solvent": 0.},
      [0.01, 0.1], [120.57219, 0.37414]],
 ]

From 1b1ea726f91c1aaf8eb2c3e1676dc4301a9cfe75 Mon Sep 17 00:00:00 2001
From: Paul Kienzle <pkienzle@nist.gov>
Date: Tue, 2 Jun 2026 12:54:19 -0400
Subject: [PATCH 4/5] change truncation to amount of truncation instead of
 amount not truncated

---
 .../models/truncated_octahedron_fibonacci.py  | 192 +++++++++---------
 1 file changed, 92 insertions(+), 100 deletions(-)

diff --git a/explore/models/truncated_octahedron_fibonacci.py b/explore/models/truncated_octahedron_fibonacci.py
index 8e20714f..7b7c31b6 100644
--- a/explore/models/truncated_octahedron_fibonacci.py
+++ b/explore/models/truncated_octahedron_fibonacci.py
@@ -1,117 +1,119 @@
 # truncated octahedron model
 # Note: model title and parameter table are inserted automatically
 r"""
-This model provides the form factor P(q) for a general octahedron with orientational averaging (Fibonacci quadrature).
-It can be a regular octahedron shape with all edges of the same length.
+This model provides the form factor P(q) for a general octahedron.
+It can be a regular octahedron, with all edges of equal length.
 Or a general shape with different elongations along the three perpendicular two-fold axes.
 It includes the possibility to add an adjustable square truncation at each of the six vertices.
 This model includes the general cuboctahedron shape for the maximum value of truncation.
 The form factor expression is obtained by analytical integration over the volume of the shape.
 This model is constructed in a similar way as the rectangular prism model.
-It contains both the form factor for a reference orientation and the 1D form factor after orientation average
-using the Fibonacci quadrature. This quadrature provides a quasi-uniform distribution of points on the unit sphere
-using the golden ratio. The only new parameter is the number of points to generate on the unit sphere.
-The default value (around 400 points) provides a good balance between accuracy and computational efficiency.
+It contains both the form factor for a reference orientation and the 1D form factor after orientation average (Gauss-Legendre).
 
 Definition
 ----------
 
-The general octahedron is defined by its dimensions along its three perpendicular two-fold axes along x, y and z directions.
-radius_a, radius_b and radius_c are the distances from the center of the general octahedron to its 6 vertices.
+This model computes the form factor of a general octahedron by defining its size through
+its circumradius :math:`radius_a` (parameter called *radius_a* in the model),
+the elongations through the ratios :math:`\frac{b}{a}` and :math:`\frac{c}{a}` (parameters called *b2a_ratio* and *c2a_ratio* in the model) and the truncation level through the
+truncation ratio *t* (parameter called *truncation* in the model).
+
+Indeed, the general octahedron is defined by its dimensions along its three perpendicular two-fold axes along x, y and z directions.
+:math:`radius_a` (parameter called *radius_a* in the model), :math:`radius_b` and :math:`radius_c` are the distances from the center of the general octahedron to its 6 vertices,
+which are equivalent to the circumradiuses of the general octahedron along the three directions.
 
 Coordinates of the six vertices are:
-    (radius_a, 0, 0),
-    (-radius_a, 0, 0),
-    (0, radius_b, 0),
-    (0, -radius_b, 0),
-    (0, 0, radius_c),
-    (0, 0, -radius_c)
-
-t is the truncation parameter.
+
+.. math::
+
+    (radius_a,\ 0,\ 0) \\
+    (-radius_a,\ 0,\ 0) \\
+    (0,\ radius_b,\ 0) \\
+    (0,\ -radius_b,\ 0) \\
+    (0,\ 0,\ radius_c) \\
+    (0,\ 0,\ -radius_c)
+
 Truncation adds a square facet for each vertex that is perpendicular to a 2-fold axis.
-The resulting shape consists of six squares and eight hexagons, which may be irregular depending on the three dimensions
-A square facet crosses the x, y, z directions at distances equal to truncation radius_a, truncation radius_b and truncation radius_c.
+The resulting shape consists of six squares and eight hexagons, which may be irregular depending on the three dimensions.
+The truncation ratio *t* (parameter called `truncation` in the model) is defined as: 0 ≤ t ≤ 0.5, 0 corresponding to no truncation
+(full octahedron) and 0.5 corresponding to the maximum truncation (cuboctahedron).
+For the following formulas, we will use the notation :math:`t_inv = 1 - t`.
+Indeed, a square facet crosses the x, y, z directions at distances equal to
+:math:`t_{\mathrm{inv}} \, radius_a`, :math:`t_{\mathrm{inv}} \, radius_b` and :math:`t_{\mathrm{inv}} \, radius_c`.
 
 A regular octahedron corresponds to:
 
 .. math::
 
-    radius_a = radius_b = radius_c, \quad t = 1
+    radius_a = radius_b = radius_c, \quad t = 0
 
 A regular cuboctahedron shape with 6 squares and 8 triangles corresponds to:
 
 .. math::
 
-    radius_a = radius_b = radius_c, \quad truncation = \frac{1}{2}
+    radius_a = radius_b = radius_c, \quad t  = \frac{1}{2}
+
 
-The model contains 4 parameters: radius_a, the two ratios b2a_ratio and c2a_ratio and truncation:
+The volume of the general shape including truncation is given by:
 
 .. math::
 
-    b2a_{\text{ratio}} = \frac{radius_b}{radius_a}, \quad
-    c2a_{\text{ratio}} = \frac{radius_c}{radius_a}, \quad
-    \frac{1}{2} < t < 1
+    V = \frac{4}{3}\, radius_{\text{a}}^{3}\, \frac{b}{a}\, \frac{c}{a}\,\bigl(1 - 3t^{3}\bigr)
 
-For a regular shape:
+The general octahedron is made of eight triangular faces. The three edge lengths
+are:
 
 .. math::
 
-    b2a_{\text{ratio}} = c2a_{\text{ratio}} = 1
+    A_{\text{edge}}^{2} = radius_{\text{a}}^{2} + radius_{\text{b}}^{2},\qquad
+    B_{\text{edge}}^{2} = radius_{\text{a}}^{2} + radius_{\text{c}}^{2},\qquad
+    C_{\text{edge}}^{2} = radius_{\text{b}}^{2} + radius_{\text{c}}^{2}
 
-Volume of the general shape including truncation is given by:
+For a regular shape (no elongation):
 
 .. math::
 
-    V = \frac{4}{3}\, length_{\text{a}}^{3}\, b2a_{\text{ratio}}\, c2a_{\text{ratio}}\,\bigl(1 - 3(1 - t)^{3}\bigr)
-
-The general octahedron is made of eight triangular faces. The three edge lengths
-are:
+    \frac{b}{a} = \frac{c}{a} = 1
 
 .. math::
 
-    A_{\text{edge}}^{2} = length_{\text{a}}^{2} + length_{\text{b}}^{2},\qquad
-    B_{\text{edge}}^{2} = length_{\text{a}}^{2} + length_{\text{c}}^{2},\qquad
-    C_{\text{edge}}^{2} = length_{\text{b}}^{2} + length_{\text{c}}^{2}
-
-For a regular shape (no elongation):
+    A_{\text{edge}} = B_{\text{edge}} = C_{\text{edge}} = radius_{\text{a}} \sqrt{2}
 
 .. math::
 
-    b2a_{\text{ratio}} = c2a_{\text{ratio}} = 1,\qquad
-    A_{\text{edge}} = B_{\text{edge}} = C_{\text{edge}} = length_{\text{a}} \sqrt{2},\qquad
-    length_{\text{a}} = length_{\text{b}} = length_{\text{c}} = A_{\text{edge}} / \sqrt{2}
+    radius_{\text{a}} = radius_{\text{b}} = radius_{\text{c}} = A_{\text{edge}} / \sqrt{2}
 
 .. math::
 
-    V = \frac{4}{3} \, length_{\text{a}}^{3} \, \bigl(1 - 3 (1 - t)^3 \bigr)
+    V = \frac{4}{3} \, radius_{\text{a}}^{3} \, \bigl(1 - 3 t^3 \bigr)
 
 The reference orientation of the shape is: a along x, b along y and c along z.
 Amplitude of the form factor AP for the reference orientation of the shape reads
 
 .. math::
 
-    AP(q,\theta,\phi) = \frac{3}{1 - 3(1 - t)^3}\,(AA + BB + CC)
+    AP(q,\theta,\phi) = \frac{6}{1 - 3t^3}\,(AA + BB + CC)
 
 .. math::
 
-    AA = \frac{1}{2\,(q_y^2 - q_z^2)\,(q_y^2 - q_x^2)}\Big[(q_y - q_x)\sin\big(q_y(1 - t) - q_x t\big)
-    + (q_y + q_x)\sin\big(q_y(1 - t) + q_x t\big)\Big]
-    + \frac{1}{2\,(q_z^2 - q_x^2)\,(q_z^2 - q_y^2)}\Big[(q_z - q_x)\sin\big(q_z(1 - t) - q_x t\big)
-    + (q_z + q_x)\sin\big(q_z(1 - t) + q_x t\big)\Big]
+    AA = \frac{1}{2\,(q_y^2 - q_z^2)\,(q_y^2 - q_x^2)}\Big[(q_y - q_x)\sin\big(q_y t - q_x t_{\text{inv}}\big)
+    + (q_y + q_x)\sin\big(q_y t + q_x t_{\text{inv}}\big)\Big] \\
+    + \frac{1}{2\,(q_z^2 - q_x^2)\,(q_z^2 - q_y^2)}\Big[(q_z - q_x)\sin\big(q_z t - q_x t_{\text{inv}}\big)
+    + (q_z + q_x)\sin\big(q_z t + q_x t_{\text{inv}}\big)\Big]
 
 .. math::
 
-    BB = \frac{1}{2\,(q_z^2 - q_x^2)\,(q_z^2 - q_y^2)}\Big[(q_z - q_y)\sin\big(q_z(1 - t) - q_y t\big)
-    + (q_z + q_y)\sin\big(q_z(1 - t) + q_y t\big)\Big]
-    + \frac{1}{2\,(q_x^2 - q_y^2)\,(q_x^2 - q_z^2)}\Big[(q_x - q_y)\sin\big(q_x(1 - t) - q_y t\big)
-    + (q_x + q_y)\sin\big(q_x(1 - t) + q_y t\big)\Big]
+    BB = \frac{1}{2\,(q_z^2 - q_x^2)\,(q_z^2 - q_y^2)}\Big[(q_z - q_y)\sin\big(q_z t - q_y t_{\text{inv}}\big)
+    + (q_z + q_y)\sin\big(q_z t + q_y t_{\text{inv}}\big)\Big] \\
+    + \frac{1}{2\,(q_x^2 - q_y^2)\,(q_x^2 - q_z^2)}\Big[(q_x - q_y)\sin\big(q_x t - q_y t_{\text{inv}}\big)
+    + (q_x + q_y)\sin\big(q_x t + q_y t_{\text{inv}}\big)\Big]
 
 .. math::
 
-    CC = \frac{1}{2\,(q_x^2 - q_y^2)\,(q_x^2 - q_z^2)}\Big[(q_x - q_z)\sin\big(q_x(1 - t) - q_z t\big)
-    + (q_x + q_z)\sin\big(q_x(1 - t) + q_z t\big)\Big]
-    + \frac{1}{2\,(q_y^2 - q_z^2)\,(q_y^2 - q_x^2)}\Big[(q_y - q_z)\sin\big(q_y(1 - t) - q_z t\big)
-    + (q_y + q_z)\sin\big(q_y(1 - t) + q_z t\big)\Big]
+    CC = \frac{1}{2\,(q_x^2 - q_y^2)\,(q_x^2 - q_z^2)}\Big[(q_x - q_z)\sin\big(q_x t - q_z t_{\text{inv}}\big)
+    + (q_x + q_z)\sin\big(q_x t + q_z t_{\text{inv}}\big)\Big] \\
+    + \frac{1}{2\,(q_y^2 - q_z^2)\,(q_y^2 - q_x^2)}\Big[(q_y - q_z)\sin\big(q_y t - q_z t_{\text{inv}}\big)
+    + (q_y + q_z)\sin\big(q_y t + q_z t_{\text{inv}}\big)\Big]
 
 Capital Qx Qy Qz are the three components in [A-1] of the scattering vector.
 qx qy qz are rescaled components (no unit) for computing AA, BB and CC terms.
@@ -124,9 +126,9 @@
 
 .. math::
 
-    q_x = Q_x \, length_{\text{a}},\qquad
-    q_y = Q_y \, length_{\text{b}},\qquad
-    q_z = Q_z \, length_{\text{c}}
+    q_x = Q_x \, radius_{\text{a}},\qquad
+    q_y = Q_y \, radius_{\text{b}},\qquad
+    q_z = Q_z \, radius_{\text{c}}
 
 
 θ is the angle between the scattering vector and the z axis.
@@ -143,14 +145,6 @@
     P(q) =  \frac{2}{\pi} \int_0^{\frac{\pi}{2}} \,
     \int_0^{\frac{\pi}{2}} A_P^2(q) \, \sin\theta \, d\theta \, d\phi
 
-The sphere is sampled using Fibonacci quadrature to provide a quasi-uniform distribution of points on the unit sphere.
-The repartition of the points is computed using the golden ratio (see fibonacci.py).
-
-
-.. figure:: img/fibonacci_sphere.png
-
-    Fibonacci sphere using 5810 points.
-
 And the 1D scattering intensity is calculated as
 
 .. math::
@@ -163,6 +157,8 @@
 is the scattering length of the solvent, and (if the data are in absolute
 units) *scale* represents the volume fraction (which is unitless).
 
+The sphere is sampled using Fibonacci quadrature to provide a quasi-uniform distribution of points on the unit sphere.
+The repartition of the points is computed using the golden ratio (see fibonacci.py).
 
 .. figure:: img/octa-truncated.png
 
@@ -170,40 +166,34 @@
 
 .. figure:: img/octahedrons_intensity_plot.png
 
-    Scattering intensity of a cuboctahedron (t=0.5) and a regular octahedron (t=1) of a = 300 Angstroms.
+    Scattering intensity of a cuboctahedron (t=0.5) and a regular octahedron (t=0) of radius_a = 400 Å.
 
 Validation
 ----------
 
 Validation of the code is made using numerical checks.
 Comparisons with Debye formula calculations were made using DebyeCalculator library (https://github.com/FrederikLizakJohansen/DebyeCalculator).
-Good agreement was found at q < 0.1 1/Angstrom.
+Good agreement was found at q < 0.1 1/Å.
 
 References
 ----------
 
-1. Wei-Ren Chen et al. "Scattering functions of Platonic solids".
-   In: Journal of Applied Crystallography - J APPL CRYST 44 (June 2011).
-   DOI: 10.1107/S0021889811011691
-
-2. Croset, Bernard, "Form factor of any polyhedron: a general compact
-   formula and its singularities" In: J. Appl. Cryst. (2017). 50, 1245–1255
-   https://doi.org/10.1107/S1600576717010147
+1. Li, X., Shew, C., He, L., Meilleur, F., Myles, D. A. A., Liu, E., Zhang, Y., Smith, G. S.,
+   Herwig, K. W., Pynn, R., & Chen, W. (2011). Scattering functions of Platonic solids.
+   *Journal Of Applied Crystallography*, 44(3), 545‑557. https://doi.org/10.1107/s0021889811011691
 
-3. Wuttke, J. Numerically stable form factor of any polygon and polyhedron
-   J Appl Cryst 54, 580-587 (2021)
-   https://doi.org/10.1107/S160057672100171
+2. Croset, B. (2017). Form factor of any polyhedron : a general compact formula and its singularities.
+   *Journal Of Applied Crystallography*, 50(5), 1245‑1255. https://doi.org/10.1107/s1600576717010147
 
-4. Álvaro González. Measurement of areas on a sphere using Fibonacci and latitude–longitude lattices.
-Mathematical geosciences 42 (2010), pp. 49–64
+3. Wuttke, J. (2021). Numerically stable form factor of any polygon and polyhedron.
+   *Journal Of Applied Crystallography*, 54(2), 580‑587. https://doi.org/10.1107/s1600576721001710
 
-5. MB Kozin, VV Volkov, and DI Svergun. ASSA, a program for three-dimensional rendering in solution
-scattering from biopolymers. Applied Crystallography 30.5 (1997), pp. 811–815
 
 Authorship and Verification
 ----------------------------
 
 * **Authors:** Marianne Imperor-Clerc (marianne.imperor@cnrs.fr)
+             Helen Ibrahim (helenibrahim1@outlook.com)
              Sara Mokhtari (smokhtari@insa-toulouse.fr)
 
 * **Last Modified by:** SM **Date:** 13th January 2026
@@ -221,14 +211,15 @@
 title = "General octahedron pure python model using Fibonacci quadrature"
 description = """
 The amplitude AP is defined as follows.
+t is the truncation parameter, tinv is 1 - t and AA, BB and CC are the three terms of the form factor for the reference orientation of the shape.
 
-AP = 6./(1.-3*(1.-t)*(1.-t)*(1.-t))*(AA+BB+CC)
+AP = 6./(1.-3*(t)**3)*(AA+BB+CC)
 
-AA = 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx)) * ((qy-qx)*sin(qy*(1.-t)-qx*t) + (qy+qx)*sin(qy*(1.-t)+qx*t)) + 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy)) * ((qz-qx)*sin(qz*(1.-t)-qx*t) + (qz+qx)*sin(qz*(1.-t)+qx*t))
+AA = 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx)) * ((qy-qx)*sin(qy*(t)-qx*tinv) + (qy+qx)*sin(qy*(t)+qx*tinv)) + 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy)) * ((qz-qx)*sin(qz*(t)-qx*tinv) + (qz+qx)*sin(qz*(t)+qx*tinv))
 
-BB = 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy)) * ((qz-qy)*sin(qz*(1.-t)-qy*t) + (qz+qy)*sin(qz*(1.-t)+qy*t)) + 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz)) * ((qx-qy)*sin(qx*(1.-t)-qy*t) + (qx+qy)*sin(qx*(1.-t)+qy*t))
+BB = 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy)) * ((qz-qy)*sin(qz*(t)-qy*tinv) + (qz+qy)*sin(qz*(t)+qy*tinv)) + 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz)) * ((qx-qy)*sin(qx*(t)-qy*tinv) + (qx+qy)*sin(qx*(t)+qy*tinv))
 
-CC = 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz)) * ((qx-qz)*sin(qx*(1.-t)-qz*t) + (qx+qz)*sin(qx*(1.-t)+qz*t)) + 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx)) * ((qy-qz)*sin(qy*(1.-t)-qz*t) + (qy+qz)*sin(qy*(1.-t)+qz*t))
+CC = 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz)) * ((qx-qz)*sin(qx*(t)-qz*tinv) + (qx+qz)*sin(qx*(t)+qz*tinv)) + 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx)) * ((qy-qz)*sin(qy*(t)-qz*tinv) + (qy+qz)*sin(qy*(t)+qz*tinv))
 
 With capital QX QY QZ are the three components in [A-1] of the scattering vector
 and qx qy qz are the rescaled components (no unit) for computing AP term.
@@ -242,26 +233,26 @@
 
 Reference orientation is with a along x axis, b along y axis and c along z axis
 
-Valid truncation parameter range: 0.5 < t < 1.
+Valid truncation parameter range: 0 ≤ t ≤ 0.5.
 
-t=1 is for octahedron
+t=0 is for octahedron
 t=0.5 is for cuboctahedron
 """
 category = "shape:polyhedron"
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["sld", "1e-6/Ang^2", 130, [-inf, inf], "sld",
+parameters = [["sld", "1e-6/Ang^2", 126., [-inf, inf], "sld",
                "Octahedron scattering length density"],
               ["sld_solvent", "1e-6/Ang^2", 9.4, [-inf, inf], "sld",
                "Solvent scattering length density"],
-              ["radius_a", "Ang", 500, [0, inf], "volume",
+              ["radius_a", "Ang", 400, [0, inf], "volume",
                "half height along a axis"],
               ["b2a_ratio", "", 1, [0, inf], "volume",
                "Ratio b/a"],
               ["c2a_ratio", "", 1, [0, inf], "volume",
                "Ratio c/a"],
-              ["truncation", "", 0.99, [0.0, 1.0], "volume",
-               "truncation along axis"],
+              ["truncation", "", 0, [0, 0.5], "volume",
+               "truncation ratio, 0 for octahedron and 0.5 for cuboctahedron"],
               ["npoints_fibonacci",     "",           400, [1, 1e6],   "",
                        "Number of points on the sphere for the Fibonacci integration"]
               ]
@@ -287,7 +278,7 @@ def form_volume(radius_a, b2a_ratio, c2a_ratio, truncation):
         * radius_a**3
         * b2a_ratio
         * c2a_ratio
-        * (1.0 - 3 * (1.0 - truncation) ** 3)
+        * (1.0 - 3 * truncation ** 3)
     )
 
 
@@ -314,6 +305,7 @@ def A(qa, qb, qc, radius_a, b2a_ratio, c2a_ratio, truncation):
     """
     radius_b = radius_a * b2a_ratio
     radius_c = radius_a * c2a_ratio
+    tinv = 1.0 - truncation
     qnx = qa * radius_a  # conversion to dimensionless coordinate
     qny = qb * radius_b  # conversion to dimensionless coordinate
     qnz = qc * radius_c  # conversion to dimensionless coordinate
@@ -326,12 +318,12 @@ def A(qa, qb, qc, radius_a, b2a_ratio, c2a_ratio, truncation):
     denom2 += eps_den * (np.abs(denom2) < eps_den)  # add epsilon where denom is zero
 
     term1 = (
-        (qny - qnx) * np.sin(qny * (1.0 - truncation) - qnx * truncation)
-        + (qny + qnx) * np.sin(qny * (1.0 - truncation) + qnx * truncation)
+        (qny - qnx) * np.sin(qny * truncation - qnx * tinv)
+        + (qny + qnx) * np.sin(qny * truncation + qnx * tinv)
     ) / denom1
     term2 = (
-        (qnz - qnx) * np.sin(qnz * (1.0 - truncation) - qnx * truncation)
-        + (qnz + qnx) * np.sin(qnz * (1.0 - truncation) + qnx * truncation)
+        (qnz - qnx) * np.sin(qnz * truncation - qnx * tinv)
+        + (qnz + qnx) * np.sin(qnz * truncation + qnx * tinv)
     ) / denom2
 
     return term1 + term2
@@ -361,7 +353,7 @@ def Fqabc(qa, qb, qc, radius_a, b2a_ratio, c2a_ratio, truncation):
     CC = A(qc, qa, qb, radius_a, b2a_ratio, c2a_ratio, truncation)
 
     # Normalisation to 1 of AP at q = 0. Division by a factor 4/3 and global 1/2 coefficient.
-    ap = 3.0 / (1.0 - 3.0 * (1.0 - truncation) ** 3) * (AA + BB + CC)
+    ap = 3.0 / (1.0 - 3.0 * truncation ** 3) * (AA + BB + CC)
 
     # Guard against numerical issues: mirror the C fallback which sets
     # invalid (NaN/inf) values to zero so tests and downstream code
@@ -418,10 +410,10 @@ def Iq(
     q,
     sld,
     sld_solvent,
-    radius_a=500,
-    b2a_ratio=1,
-    c2a_ratio=1,
-    truncation=0.99,
+    radius_a=500.,
+    b2a_ratio=1.,
+    c2a_ratio=1.,
+    truncation=0.,
     npoints_fibonacci: int = 400,
 ):
     """

From c8821e006a0b25cd18a9e2b3275b7a00412d8bcd Mon Sep 17 00:00:00 2001
From: marimperorclerc <95975874+marimperorclerc@users.noreply.github.com>
Date: Wed, 3 Jun 2026 12:13:28 +0200
Subject: [PATCH 5/5] Update truncated_octahedron_fibonacci.py

states at the beginning that this is a pure python version using Fibonacci integration.
---
 explore/models/truncated_octahedron_fibonacci.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/explore/models/truncated_octahedron_fibonacci.py b/explore/models/truncated_octahedron_fibonacci.py
index 7b7c31b6..fb292dea 100644
--- a/explore/models/truncated_octahedron_fibonacci.py
+++ b/explore/models/truncated_octahedron_fibonacci.py
@@ -1,4 +1,5 @@
 # truncated octahedron model
+# pure python version using Fibonacci integration.
 # Note: model title and parameter table are inserted automatically
 r"""
 This model provides the form factor P(q) for a general octahedron.
@@ -8,7 +9,7 @@
 This model includes the general cuboctahedron shape for the maximum value of truncation.
 The form factor expression is obtained by analytical integration over the volume of the shape.
 This model is constructed in a similar way as the rectangular prism model.
-It contains both the form factor for a reference orientation and the 1D form factor after orientation average (Gauss-Legendre).
+It contains both the form factor for a reference orientation and the 1D form factor after orientation average (Fibonacci integration).
 
 Definition
 ----------