Octahedron_truncated_txtytz - octahedron_truncated_txtytz.py

    # 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 the possibility to be truncated differently along 3 directions (x,y,z).
This model is constructed in a similar way as the octahedron_truncated model.
For example, it is suited to model square bipyramids or square bifrustum.

Definition
----------

The general octahedron is defined by its dimensions along its three perpendicular two-fold axis 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)

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)

tx, ty and tz are the three truncation parameters. 
No truncation corresponds to tx=ty=tz=1.
Truncation adds a square facet for each vertice that is perpendicular to a 2-fold axis.
The resulting shape is made of 6 squares and 8 hexagons, that may be non regular depending on the dimensions.
A square facet crosses the x,y,z directions at distances equal to tx*length_a, ty*length_b and tz*length_c.


The model contains 6 parameters: length_a, the two ratios b2a_ratio and c2a_ratio and tx, ty and tz:
    length_a
    b2a_ratio = length_b/length_a
    c2a_ratio = length_c/length_a
    0 < tx < 1
    0 < ty < 1
    0 < tz < 1

For a regular shape:
    b2a_ratio = c2a_ratio = 1

For example, a regular square bifrustum shape with the 4-fold axis along $z$ corresponds to :
    b2a_ratio = c2a_ratio = 1
    tx=ty=1
    0<tz<1

An elongated square bipyramid shape without truncation with the 4-fold axis along $z$ corresponds to :
    b2a_ratio = 1
    c2a_ratio > 1
    tx=ty=tz=1

.. math::

Volume of the shape is:
V = (4/3) * length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio) * (1.-(1.-tx)^3-(1.-ty)^3-(1.-(1.-tz)^3)

Lengths of edges are equal to :
    A_edge^2 = length_a^2+length_b^2
    B_edge^2 = length_a^2+length_c^2
    C_edge^2 = length_b^2+length_c^2


For a regular octahedron :
    b2a_ratio = c2a_ratio = 1
    A_edge = B_edge = C_edge = length_a*sqrt(2)
    length_a = length_b = length_c = A_edge/sqrt(2)
    V = (4/3) * length_a^3 = (sqrt(2)/3) * A_edge^3 * (1.-(1.-tx)^3-(1.-ty)^3-(1.-(1.-tz)^3)

Amplitude of the form factor AP is calculated with a scaled scattering vector (qx,qy,qz) :
    AP = 6./(1.-(1.-tx)*(1.-tx)*(1.-tx)-(1.-ty)*(1.-ty)*(1.-ty)-(1.-tz)*(1.-tz)*(1.-tz))*(AA+BB+CC);

    with:

    AA = 1./((qy*qy-qz*qz)*(qy*qy-qx*qx))*1./2.*((qy-qx)*sin(qy*(1.-tx)-qx*tx)+(qy+qx)*sin(qy*(1.-tx)+qx*tx))+
                                1./((qz*qz-qx*qx)*(qz*qz-qy*qy))*1./2.*((qz-qx)*sin(qz*(1.-tx)-qx*tx)+(qz+qx)*sin(qz*(1.-tx)+qx*tx));

    BB = 1./((qz*qz-qx*qx)*(qz*qz-qy*qy))*1./2.*((qz-qy)*sin(qz*(1.-ty)-qy*ty)+(qz+qy)*sin(qz*(1.-ty)+qy*ty))+
                                1./((qx*qx-qy*qy)*(qx*qx-qz*qz))*1./2.*((qx-qy)*sin(qx*(1.-ty)-qy*ty)+(qx+qy)*sin(qx*(1.-ty)+qy*ty));

    CC = 1./((qx*qx-qy*qy)*(qx*qx-qz*qz))*1./2.*((qx-qz)*sin(qx*(1.-tz)-qz*tz)+(qx+qz)*sin(qx*(1.-tz)+qz*tz))+
                                1./((qy*qy-qz*qz)*(qy*qy-qx*qx))*1./2.*((qy-qz)*sin(qy*(1.-tz)-qz*tz)+(qy+qz)*sin(qy*(1.-tz)+qz*tz));

    and :
            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

$\theta$ is the angle between the scattering vector and the $z$ axis of the octahedron ($length_c$).
$\phi$ is the angle between the scattering vector (lying in the $xy$ plane) and the $x$ axis ($length_a$).

The normalized form factor in 1D is obtained averaging over all possible
orientations. Same code as for a rectangular prism.

.. 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

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 octahedron, $\rho_\text{p}$
is the scattering length inside the volume, $\rho_\text{solvent}$
is the scattering length of the solvent, and (if the data are in absolute
units) *scale* represents the volume fraction (which is unitless).

For 2d data the orientation of the particle is required, described using
angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
of the calculation and angular dispersions see :ref:`orientation` .
The angle $\Psi$ is the rotational angle around the long *C* axis. For example,
$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector.

For 2d, constraints must be applied during fitting to ensure that the inequality
$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error,
but the results may be not correct.

.. figure:: img/parallelepiped_angle_definition.png

    Definition of the angles for oriented core-shell parallelepipeds.
    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then
    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder.
    The neutron or X-ray beam is along the $z$ axis.

.. figure:: img/parallelepiped_angle_projection.png

    Examples of the angles for oriented rectangular prisms against the
    detector plane.



Validation
----------

Validation of the code is made using numerical checks.


References
----------
1. Wei-Ren Chen et al. “Scattering functions of Platonic solids”.
In:Journal of AppliedCrystallography - J APPL CRYST44 (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


Authorship and Verification
----------------------------

* **Authors: Marianne Imperor-Clerc (marianne.imperor@universite-paris-saclay.fr)
             Helen Ibrahim (helenibrahim1@outlook.com)**
* **Last Modified by MI: 9 October 2020**
* **Last Reviewed by MI: March 2023:**
"""

import numpy as np
from numpy import inf

name = "octahedron_truncated_txtytz"
title = "Truncated Octahedron txtytz."
description = """
            AP = 6./(1.-(1.-tx)*(1.-tx)*(1.-tx)-(1.-ty)*(1.-ty)*(1.-ty)-(1.-tz)*(1.-tz)*(1.-tz))*(AA+BB+CC);
           const double AA = 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx))*((qy-qx)*sin(qy*(1.-tx)-qx*tx)+(qy+qx)*sin(qy*(1.-tx)+qx*tx))+
                                1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy))*((qz-qx)*sin(qz*(1.-tx)-qx*tx)+(qz+qx)*sin(qz*(1.-tx)+qx*tx));

            const double BB = 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy))*((qz-qy)*sin(qz*(1.-ty)-qy*ty)+(qz+qy)*sin(qz*(1.-ty)+qy*ty))+
                                1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz))*((qx-qy)*sin(qx*(1.-ty)-qy*ty)+(qx+qy)*sin(qx*(1.-ty)+qy*ty));

            const double CC = 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz))*((qx-qz)*sin(qx*(1.-tz)-qz*tz)+(qx+qz)*sin(qx*(1.-tz)+qz*tz))+
                                1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx))*((qy-qz)*sin(qy*(1.-tz)-qz*tz)+(qy+qz)*sin(qy*(1.-tz)+qz*tz));

normalisation to 1. of AP at q = 0. Division by a Factor 4/3.
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;
    	    0 < tx < 1
    	    0 < ty < 1
    	    0 < tz < 1


"""
category = "shape:parallelepiped"

#             ["name", "units", default, [lower, upper], "type","description"],
parameters = [["sld", "1e-6/Ang^2", 128, [-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", 35, [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"],
              ["tx", "", 0.99, [0.0, 1.0], "volume",
               "truncation along x axis"],
              ["ty", "", 0.99, [0.0, 1.0], "volume",
               "truncation along y axis"],
              ["tz", "", 0.99, [0.0, 1.0], "volume",
               "truncation along z axis"],
              ["theta", "degrees", 0, [-360, 360], "orientation",
               "c axis to beam angle"],
              ["phi", "degrees", 0, [-360, 360], "orientation",
               "rotation about beam"],
              ["psi", "degrees", 0, [-360, 360], "orientation",
               "rotation about c axis"],
             ]

source = ["lib/gauss20.c", "octahedron_truncated_txtytz.c"]
# change to "lib/gauss76.c" or "lib/gauss150.c" to increase the number of integration points
# for the orientational average. Note that it will increase calculation times.
have_Fq = True
Back to Model Download