research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoJOURNAL OF
APPLIED
CRYSTALLOGRAPHY
ISSN: 1600-5767

Numerically stable form factor of any polygon and polyhedron

CROSSMARK_Color_square_no_text.svg

aForschungszentrum Jülich GmbH, Jülich Centre for Neutron Science (JCNS) at Heinz Maier-Leibnitz Zentrum (MLZ), Lichtenbergstrasse 1, 85748 Garching, Germany
*Correspondence e-mail: j.wuttke@fz-juelich.de

Edited by D. I. Svergun, European Molecular Biology Laboratory, Hamburg, Germany (Received 15 January 2021; accepted 11 February 2021; online 25 March 2021)

Coordinate-free expressions for the form factors of arbitrary polygons and polyhedra are derived using the divergence theorem and Stokes's theorem. Apparent singularities, all removable, are discussed in detail. Cancellation near the singularities causes a loss of precision that can be avoided by using series expansions. An important application domain is small-angle scattering by nanocrystals.

1. Introduction

1.1. Overview

The term `form factor' has different meanings in science and in engineering. Here, we are concerned with the form factor of a geometric figure as defined in the physical sciences, namely the Fourier transform of the figure's indicator function, also called the `shape transform' of the figure.

This form factor has important applications in the emission, detection and scattering of radiation. Two-dimensional shape transforms are used in the theory of reflector antennas (Lee & Mittra, 1983[Lee, S. W. & Mittra, R. (1983). IEEE Trans. Antennas Propagat. 31, 99-103.]). The three-dimensional form factors of the sphere and the cylinder go back to Lord Rayleigh (1881[Rayleigh, Lord (1881). Philos. Mag. 12, 81-101.]). Shapes of three-dimensional nanoparticles are investigated by neutron and X-ray small-angle scattering (Hammouda, 2010[Hammouda, B. (2010). Probing Nanoscale Structures - The SANS Toolbox, https://www.ncnr.nist.gov/staff/hammouda/the_SANS_toolbox.pdf.]). Particles grown on a substrate (Henry, 2005[Henry, C. R. (2005). Prog. Surf. Sci. 80, 92-116.]) develop many different shapes, especially polyhedral ones, as observed by grazing-incidence neutron and X-ray small-angle scattering (GISAS, GISANS, GISAXS) (Renaud et al., 2009[Renaud, G., Lazzari, R. & Leroy, F. (2009). Surf. Sci. Rep. 64, 255-380.]). Large collections of particle shape transforms have therefore been derived for and implemented in GISAS software (Lazzari, 2006[Lazzari, R. (2006). IsGISAXS, Version 2.6 of 4 May 2006, https://www.insp.jussieu.fr/oxydes/IsGISAXS/isgisaxs.htm.]; Pospelov et al., 2020[Pospelov, G., Van Herck, W., Burle, J., Carmona Loaiza, J. M., Durniak, C., Fisher, J. M., Ganeva, M., Yurov, D. & Wuttke, J. (2020). J. Appl. Cryst. 53, 262-276.]); another GISAS software package uses surface triangulation for computing approximative form factors (Chourou et al., 2013[Chourou, S. T., Sarje, A., Li, X. S., Chan, E. R. & Hexemer, A. (2013). J. Appl. Cryst. 46, 1781-1795.]). For attempts at direct reconstruction of polyhedral shapes from scattering patterns see the article by Engel & Laasch (2020[Engel, K. & Laasch, B. (2020). arxiv:2011.06971v1 [Math-ph].]) and literature cited therein.

In this paper, we derive a numerically stable algorithm for computing the form factor of any polygon or polyhedron, as implemented in the GISAS software BornAgain (Pospelov et al., 2020[Pospelov, G., Van Herck, W., Burle, J., Carmona Loaiza, J. M., Durniak, C., Fisher, J. M., Ganeva, M., Yurov, D. & Wuttke, J. (2020). J. Appl. Cryst. 53, 262-276.]). Originally, this algorithm was documented in a terse mathematical note (Wuttke, 2017[Wuttke, J. (2017). arxiv:1703.00255.]). In the present paper, derivations and results have been simplified, the material has been completely reorganized for better readability, and additional literature is taken into account.

1.2. Different ways to compute form factors

The form factor of a three-dimensional solid body [\Pi\subset{\bb R}^{3}] is

[F({{\bf q}},\Pi):= \textstyle\int\!\!\int\limits_{\Pi}\!\!\int {\rm d}^{3}r\,\exp({i{\bf q}{\bf r}}).\eqno(1)]

In most applications, the wavevectors q are real. In GISAS, however, the incident and scattered radiation may undergo substantial absorption, which can be modeled by an imaginary part of q. Therefore, we admit complex wavevectors [{{\bf q}}\in{\bb C}^{3}].

For any polyhedron, (1)[link] can be evaluated analytically by successive integration in the three coordinates. This is straightforward for a cuboid with edges along the coordinate axes. In most other cases, the algebra is cumbersome, and the resulting expressions are complicated and unattractive in that they do not reflect symmetries of the underlying shape. Striking examples are the form factors of the Platonic solids worked out in a tour de force by Li et al. (2011[Li, X., Shew, C.-Y., He, L., Meilleur, F., Myles, D. A. A., Liu, E., Zhang, Y., Smith, G. S., Herwig, K. W., Pynn, R. & Chen, W.-R. (2011). J. Appl. Cryst. 44, 545-557.]).

It is therefore preferable to derive a coordinate-free solution of (1)[link] that expresses the form factor of a generic polyhedron in terms of its topology and vertex coordinates. This has been undertaken in different ways by Senesi & Lee (2015[Senesi, A. & Lee, B. (2015). J. Appl. Cryst. 48, 565-577.]), Croset (2017[Croset, B. (2017). J. Appl. Cryst. 50, 1245-1255.]) and Wuttke (2017[Wuttke, J. (2017). arxiv:1703.00255.]). Senesi & Lee (2015[Senesi, A. & Lee, B. (2015). J. Appl. Cryst. 48, 565-577.]) decomposed the polyhedron into pyramids and wrote the polyhedral form factor as the sum of the pyramidal form factor evaluated at different rotated q. Croset (2017[Croset, B. (2017). J. Appl. Cryst. 50, 1245-1255.]) decomposed the polyhedron into simplexes, as explored previously by Lien & Kajiya (1984[Lien, S. & Kajiya, J. T. (1984). IEEE Comput. Graph. 4, 35-42.]) and most recently by Li & Xie (2020[Li, S. Z. & Xie, Z. Q. (2020). Appl. Math. Comput. 376, 121540.]), for the integration of multinomials. Following Wuttke (2017[Wuttke, J. (2017). arxiv:1703.00255.]), we here present a different derivation that is based on use of the divergence theorem and Stokes's theorem to reduce the volume integral to integrals over polygonal faces, and further reduce these surface integrals to line integrals over straight edges. This approach has been previously used for the integration of polynomials (Cattani & Paoluzzi, 1990[Cattani, C. & Paoluzzi, A. (1990). Comput. Aided Des. 22, 130-135.]; Bernadini, 1991[Bernadini, F. (1991). Comput. Aided Des. 23, 51-58.]) and for the computation of inertia moments (Mirtich, 1996[Mirtich, B. (1996). J. Graphics Tools 1, 31-50.]).

Applications to nanoparticle assemblies typically require some averaging over particle sizes or/and orientations. How to compute these averages efficiently and with sufficient accuracy is an interesting and important question, which however is beyond the scope of the present work.

1.3. Singularities and asymptotes

All analytical expressions for polyhedral form factors, derived by whatever method, contain denominators that vanish at q = 0. Croset (2017[Croset, B. (2017). J. Appl. Cryst. 50, 1245-1255.]) suggested, and we will confirm, that the degree of this singularity is closely related to the asymptotic envelope of F(q, Π) for large q, which goes as q−1, q−2 or q−3 depending on whether q is perpendicular to a face or an edge or points in an off-symmetric direction.

However, there is nothing fundamental about the singularities at q = 0: From the definition (1)[link] in conjunction with Leibniz's integral rule of differentiation we see that F(q, Π) is infinitely many times differentiable for all [{{\bf q}}\in{\bb C}^{3}]; therefore F is a holomorphic function of each of the Cartesian components of q; therefore any apparent singularity is removable. Croset (2018[Croset, B. (2018). J. Appl. Cryst. 51, 1005-1012.]) rederived the asymptotic envelopes by classifying the endpoint singularities of the section normal to q as function of height. In Section 3.6[link], we obtain them directly from our form factors.

The main purpose of this paper is to overcome numeric instabilities for small q and q. The latter is the wavevector component in the plane of a polygonal face. We will explain how rounding errors can grossly distort form factors when q or q is of the order of ε/a, where ε is the machine precision and a is a typical edge length.

At this point the reader may wonder whether wavevectors with extremely small, but nonzero, q or q have any practical importance. If wavevectors were drawn at random from an entire Brillouin zone, then the chance of ever hitting numerically problematic values would indeed be negligible. Often, however, q is chosen along a face normal. Roundoff errors then easily yield a tiny nonzero q, which causes huge, and symmetry breaking, errors in the form factor. Actually, this entire study started from the unexpected discovery of such artifacts in conventionally computed form factors.

2. Polygon form factor

2.1. Notation

A flat polygon Γ, embedded in three-dimensional space, shall be specified by its J vertices Vj (j = 1, …, J). Vertex indices shall be understood modulo J so that V0VJ. With this convention, the vertex sequence forms the closed loop ∂Γ. Edge j of the polygon is a straight line from Vj−1 to Vj. In most of this work it is more advantageously specified through its position

[{{\bf R}}_{j}:=({{\bf V}}_{j}+{{\bf V}}_{j-1})/2\eqno(2)]

and mid-to-end vector

[{{\bf E}}_{j}:=({{\bf V}}_{j}-{{\bf V}}_{j-1})/2.\eqno(3)]

The normal vector [{{\hat{\bf n}}}] of the plane spanned by Vj shall be oriented such that ∂Γ has the winding number +1 (fulfills the right-hand rule with respect to [{{ \hat{\bf n}}}]). The oriented plane characterized by [{{\hat{\bf n}}}] induces a decomposition of any vector [{{\bf v}}\in{\bb C}^{3}] into a component perpendicular to the plane,

[{{\bf v}}_{\perp}:=({{\bf v}}{{ \hat{\bf n}}}){{\hat{\bf n}}},\eqno(4)]

and an in-plane component,

[{{\bf v}}_{\parallel}:={{\bf v}}-{{\bf v}}_{\perp}.\eqno(5)]

This decomposition will be applied to position vectors r and to wavevectors q. The oriented plane is fully specified by its normal vector [{{\hat{\bf n}}}] and its distance from the origin, r.

Complex conjugation is denoted by a superscript asterisk. The absolute value of a complex vector is written [q:=|{{\bf q}}| = ({{\bf q}}{{\bf q}}^{*})^{1/2}].

Note that the in-plane unit vector

[\widehat{{{\bf q}}_{\parallel}}:={{\bf q}}_{\parallel}/q_{\parallel}\eqno(6)]

differs from the in-plane component [{{\hat{\bf q}}}_{\parallel}] of the unit vector [{{\hat{\bf q}}}]. In this work, we shall only use [\widehat{{{\bf q}}_{\parallel}}] and [{{\hat{\bf q}}}], not [{{\hat{\bf q}}}_{\parallel}].

The triple product is denoted

[[{{\bf a}},{{\bf b}},{{\bf c}}]:={{\bf a}}\cdot({{\bf b} }\times{{\bf c}}),\eqno(7)]

with the standard operators dot (·) for the scalar product and cross (×) for the vector product. Between adjacent vector symbols, as in the parentheses in (4)[link], we omit the dot.

The cardinal sine function [{\rm sinc}(z):=\sin(z)/z] has the analytic continuation [{\rm sinc}(0) = 1]. The numeric implementation for |z| → 0 is unproblematic: as sin(z) has full floating-point accuracy, so has sin(z)/z.

2.2. Form factor

We define the form factor of a flat figure Γ, embedded in three-dimensional space, as

[f({{\bf q}},\Gamma):= \mathop{\textstyle\int\limits_{\,\,\,\Gamma}\!\!\!\int} {\rm d}^{2}r\,\exp({i{{\bf q}}{{\bf r}}}).\eqno(8)]

Proposition

The form factor (8)[link] of a flat J-gon Γ is

[f({{\bf q}},\Gamma) = {{2} \over {iq_{\parallel}}}\sum_{j = 1}^{J}\left[{{\hat{\bf n}}},\widehat{{{\bf q}}_{\parallel}}^{*},{{\bf E}}_{j}\right]\left [{\rm sinc}({{\bf q}}{{\bf E}}_{j})\exp({i{{\bf q}}{ {\bf R}}_{j}})-c\right]\eqno(9)]

for q ≠ 0, with notations from Section 2.1[link], and with an arbitrary constant c that can be chosen for computational convenience. The value at q = 0 is the area of Γ,

[f(0,\Gamma)\equiv{\rm Ar}(\Gamma) = \textstyle{{1} \over {2}}\sum\limits_{j = 1}^{J}\left[{ {\hat{\bf n}}},{{\bf V}}_{j-1},{{\bf V}}_{j}\right].\eqno(10)]

Values at q ≠ 0, q = 0 can be obtained from

[f({{\bf q}},\Gamma) = \exp({i{{\bf q}}_{\perp}{{\bf r}}_{\perp} })\,f({{\bf q}}_{\parallel},\Gamma).\eqno(11)]

Proof

For any vector field G, we have Stokes's theorem:

[\textstyle\int\limits_{\,\,\Gamma}\!\!\!\int{\rm d}r^{2}\,{{\hat{\bf n}}}\cdot({\nabla}\times{ {\bf G}}) = \oint\limits_{\partial\Gamma} {\rm d}{{\bf r}}\cdot{{\bf G}}.\eqno(12)]

To prove (9)[link], we choose [{{\bf G}}({{\bf r}}):={{\hat{\bf n}}}\times{{\bf q}}^{*} [\exp(i{{\bf q}}{{\bf r}})-c ]]. The left-hand side of (12)[link] is

[\eqalignno{\textstyle\int\limits_{\,\,\,\Gamma}\!\!\!\int {\rm d}r^{2}\,{{\hat{\bf n}}}\cdot({ \nabla}\times{{\bf G}}) & = {{\hat{\bf n}}}\cdot[i{{\bf q}}\times({ {\hat{\bf n}}}\times{{\bf q}}^{*})]\textstyle\int\limits_{\,\,\,\,\Gamma}\!\!\!\int {\rm d}r^{2}\, \exp({i{{\bf q}}{{\bf r}}}) \cr & = i{|{{\hat{\bf n}}}\times{{\bf q}}|}^{2}\,f({{\bf q}},\Gamma) \cr & = iq_{\parallel}^{2}\,f({{\bf q}},\Gamma).&(13a)}]

The right-hand side of (12)[link] is

[\textstyle\oint\limits_{\partial\Gamma}{\rm d}{{\bf r}}\cdot{{\bf G}} = \sum\limits_{j = 1}^{J} \int\limits_{{{\bf V}}_{j-1}}^{{{\bf V}}_{j}}\!\rm {\rm d}{{\bf r}}\cdot{ {\bf G}}.\eqno(13b)]

Each edge can be written as a parametric curve [{{\bf r}}_{j}(\lambda):={{\bf R}}_{j}+{{\hat{\bf E}}}_{j}\lambda] so that

[\eqalignno{\textstyle\oint\limits_{\partial\Gamma}{\rm d}{{\bf r}}\cdot{{\bf G}}& = \textstyle\sum\limits_{j = 1}^{J}{{\hat{\bf E}}}_{j}\cdot\int\limits_{-E}^{+E}\, {\rm d}\lambda\, {{\bf G}}[{{\bf r}}_{j}(\lambda)] \cr & = \textstyle\sum\limits_{j = 1}^{J}\left[{{\hat{\bf E}}}_{j},{{\hat{\bf n}}},{ {\bf q}}^{*}\right]\int\limits_{-E}^{+E}\, {\rm d}\lambda\,\left\{\exp[ i{\bf q}{\bf r}_{j}(\lambda)]-c\right\} \cr & \textstyle = \textstyle\sum\limits_{j = 1}^{J}\left[\hat{\bf n},{\bf q}_{\parallel} ^{*},{ \hat{\bf E}}_{j}\right]2E\left[{\rm sinc}({{\bf q}}{ {\bf E}}_{j})\exp(i{\bf q}{\bf R}_{j})-c\right] \cr & \textstyle = 2q_{\parallel}\sum\limits_{j = 1}^{J}\left[{\hat{\bf n}},\widehat{ {\bf q}_{\parallel}}^{*},{{\bf E}}_{j}\right]\left[{\rm sinc}({\bf q}{\bf E}_{j})\exp(i{\bf q}{\bf R}_{j})-c \right].&(14)}]

With q ≠ 0, we obtain (9)[link]. Equation (11)[link] follows directly from (8)[link] and the fact that Γ is a flat figure with constant r. To prove (10)[link], we use Stokes's theorem (12)[link] with [{{\bf G}} = {{\hat{\bf n}}}\times{{\bf r}}/2]. □

2.3. Remarks and example

A closed expression for the form factor of the polygon has long since been known (Lee & Mittra, 1983[Lee, S. W. & Mittra, R. (1983). IEEE Trans. Antennas Propagat. 31, 99-103.], equation 6). A more symmetric expression was obtained by Croset (2017[Croset, B. (2017). J. Appl. Cryst. 50, 1245-1255.]) (equation 4, where i+1 and [i-1] should be swapped). In our notation, it reads

[f({{\bf q}},\Gamma) = {{\hat{\bf n}}}\cdot\sum_{j = 1}^{J}{{{{\bf E}} _{j-1}\times{{\bf E}}_{j}} \over {({{\bf q}}{{\bf E}}_{j-1})({{\bf q}}{ {\bf E}}_{j})}}\exp(i{\bf q}{\bf V}_{j-1}).\eqno(15)]

The equivalence with our equation (9)[link] is proven in Appendix A[link]. Equation (15)[link] is esthetically more pleasing than (9)[link], but (15)[link] is problematic for computer implementation and ill suited for the theoretical study of singularities, because for each j there are two q planes for which the denominator vanishes.

Equation (10)[link] is well known as the `surveyor's formula'. The standard proof uses triangular tessellation (Braden, 1986[Braden, B. (1986). Collect. Math. J. 17, 326-337.]). See Appendix B[link] for a derivation of (10)[link] from the q → 0 limit of (9)[link].

Fig. 1[link] shows for an equilateral triangle how strongly the form factor, plotted as |f(q)| versus q, varies with the wavevector direction [{{\hat{\bf q}}}].

[Figure 1]
Figure 1
Modulus of the form factor of an equilateral triangle as function of wavenumber q for three different wavevector directions [{{\hat{\bf q}}}]. The triangle lies in the xy plane and has a symmetry axis along x. The center of gravity is at the origin. The edge length is L = 1.

2.4. Removable singularity

The closed expression (9)[link] for the polygonal form factor f(q, Γ) has a singular factor [q_{\parallel}^{-1}]. As discussed in Section 1.3[link], the definition of f guarantees its analyticity. So we know that the apparent singularity at q → 0 is removable. We also know the value at q = 0, given by (10)[link] and (11)[link]. Nonetheless, the presence of a divergent factor may cause numeric instabilities for small values of q. To investigate this more closely, let us write (9)[link] as

[f({{\bf q}},\Gamma) = 2\textstyle\sum\limits_{j = 1}^{J}\left[{{\hat{\bf n}}},\widehat{{ {\bf q}}_{\parallel}}^{*},{{\bf E}}_{j}\right]\tau_{j}({{\bf q}},c)\eqno(16)]

with the function

[\tau_{j}({{\bf q}},c):={{{\rm sinc}({{\bf q}}{ {\bf E}}_{j})\exp({i{{\bf q}}{{\bf R}}_{j}})-c} \over {iq_{\parallel}}}.\eqno(17)]

The constant c can be chosen differently for different q. At large q, c = 0 prevents roundoff errors in the numerator of (17)[link]. For small q, c = 1 is the better choice. To see this, we expand τj as

[\tau_{j}({{\bf q}},c) = {{1-c} \over {iq_{\parallel}}}+{{\hat{\bf q}}}{\bf {R}}_{j}+{\cal O}(q).\eqno(18)]

The leading, j-independent term in the expansion, with the apparent singularity (iq)−1, contributes nothing to the sum in (16)[link], whatever the value of c, because [\textstyle\sum_{j}{{\bf E}}_{j} = 0]. This, however, holds only in exact arithmetics; in a floating-point implementation, roundoff errors can make the sum nonzero. For q → 0, any such error outgrows all other terms in the expansion (Fig. 2[link]). Therefore, in the small-q regime, the only sensible choice is c = 1, which lets the leading term vanish.

[Figure 2]
Figure 2
Same form factor modulus as in Fig. 1[link], for [{{\hat{\bf q}}} = (0.6,0.8,0)], now in a double logarithmic plot, computed in double-precision arithmetic according to the closed expression (9)[link] with c = 0. For some small q, the results are totally wrong owing to imperfect cancellation in the leading-order term that ought to vanish.

Unfortunately though, the subtraction of c = 1 can cause roundoff errors in the numerator of (17)[link]. As a straightforward remedy, we compute τj for small q from its series expansion

[\tau_{j}({{\bf q}},1) = {(iq_{\parallel})}^{-1}\sum_{n = 1}^{\infty}{(iq)}^{n} \sum_{l = 0}^{2l\leq n}{{{\left({{\hat{\bf q}}}{{\bf R}}_{j}\right)}^{ n-2l}} \over {(n-2l)!}}{{{\left({{\hat{\bf q}}}{{\bf E}}_{j}\right)}^{2l}} \over {(2l+1)!}}.\eqno(19)]

As a consistency check, we note that the limit [\tau_{j}(q\to0)], at constant in-plane direction [{{\hat{\bf q}}}], is the n = 1 term [{{\hat{\bf q}}}{{\bf R}}_{j}]. Plugging this value into (16)[link], we recover the surveyor's formula (10)[link]. The algebra is quite lengthy and therefore relegated to Appendix B.

Fig. 3[link] compares the series expansion with the closed expression (9)[link]. The series expansion works well even beyond the first minima in f(|q|). In practice, the series expansion is only needed for qL << 1, and therefore only a few expansion orders are needed to keep errors close to machine precision.

[Figure 3]
Figure 3
Modulus of the form factor of the equilateral triangle of Fig. 1[link], as function of wavenumber q for wavevector direction [{{\hat{\bf q}}} = (0,1,0)]. The black chain is computed using the analytical expression (9)[link]. The colored curves are computed using the series expansion (19)[link] up to the indicated order n.

2.5. Polygon with inversion center

Computations can be simplified, and the numeric instability at q → 0 avoided, if a polygon has a perpendicular twofold symmetry axis (Schoenflies group S2) or, equivalently, an inversion center at ρ. The form factor has the symmetry f(q + q, Γ) = f(qq, 2ρΓ). As the number of vertices is even, we can write J = 2J′ and use [{{\bf R}}_{j+J\,^{\prime}}-{{ \boldrho}} = -({{\bf R}}_{j}-{{ \boldrho}})] and [{{\bf E}}_{j+J^{\prime}} = -{{\bf E}}_{j}]. For q ≠ 0, the form factor is

[f({{\bf q}},\Gamma) = \exp({i{{\bf q}}\,{{\boldrho}}}){{4} \over {q_{ \parallel}}}\sum_{j = 1}^{J^{\prime}}\left[{{\hat{\bf n}}},\widehat{{{\bf q }}_{\parallel}}^{*},{{\bf E}}_{j}\right]\,{\rm sinc}({{\bf q}}{ {\bf E}}_{j})\sin[{{\bf q}}({{\bf R}}_{j}-{{ \boldrho}})].\eqno(20)]

In contrast to (9)[link], the summand in (20)[link] has no constant contribution but is linear in q. There is no cancellation for q → 0 and no need to use a series expansion for the accurate computation of f.

3. Polyhedron form factors

3.1. Notation and parameterization

An orientable polyhedron Π shall be specified by its K polygonal faces Γk (k = 1, …, K). For each face Γk, the normal [{{\hat{\bf n}}}_{k}:={{\hat{\bf n}}}(\Gamma_{k})] shall point to the outside of Π; this then determines the order of the vertices in the sequence [{{\bf V}}_{k1},\ldots,{{\bf V}}_{kJ_{k}}].

In a computer implementation, the topology and geometry of a polyhedron can be specified through two arrays: Array C holds one coordinate triple Vα for each of the polyhedron's vertices. Array T holds one array γk for each of the polyhedron's faces Γk; γk holds the global indices αkj of the vertices that belong to face Γk, such that [{{\bf V}}_{\alpha_{kj}} = {{\bf V}}_{kj}]. In short, array C holds the coordinates and array T holds the topology of the polyhedron. For the latter, Schlegel diagrams (Fig. 4[link]) provide a helpful visualization. In physical simulations, C is typically generated by a parametric function, whereas T is static. An assertion in the computer code should ensure that all faces are planar for any geometry parameters.

[Figure 4]
Figure 4
Schlegel diagram of a facetted cube. This representation of polyhedral topology can help to assign vertex indices in a systematic way, which then facilitates the coding of C and T. The topology array T has elements (0, 1, 2, 3, 4, 5, 6, 7), (0, 8, 1) etc. The coordinate array C, parameterized on lengths a and b, has elements (a, b, a), (b, a, a), (−b, a, a) etc.

Additionally, it is advantageous to foresee boolean parameters to indicate the presence or absence of inversion centers. One needs one such parameter for the entire polyhedron and one for each of its polygonal faces.

3.2. Form factor

Proposition

If q ≠ 0, then the form factor (1)[link] of a K-hedron Π is

[F({{\bf q}},\Pi) = {{1} \over {iq}}\,\sum\limits_{k = 1}^{K}{{\hat{\bf q}}}^{*}{{\hat{\bf n}}}_{k}\left[f({{\bf q}},\Gamma_{k})-Cf(0,\Gamma_{k})\right]\eqno(21)]

with an arbitrary constant C that can be chosen for computational convenience. Otherwise, for q = 0, the form factor is just the volume of Π,

[F(0,\Pi)\equiv{\rm Vol}(\Pi) = \textstyle{{1} \over {3}}\sum\limits_{k}{\rm Ar}(\Gamma_{k})\,r_{\perp k}.\eqno(22)]

Proof

For a polyhedron, the divergence theorem takes the form

[\textstyle\int\!\!\int\limits_{\Pi}\!\!\int {\rm d}^{3}r\,\,{\nabla}{{\bf H}} = \int\limits_{\,\,\partial\Pi}\!\!\!\!\int {\rm d}^{2}r\,\,{{\hat{\bf n}}}\,{{\bf H}} = \sum\limits_{k}\int\!\!\int\limits_{\!\!\Gamma\limits_{k}} {\rm d}^{2}r\,\,{{\hat{\bf n}}}_{k}{{\bf H}}.\eqno(23)]

With the choice [{{\bf H}}:={{\bf q}}^{*} [\exp(i{{\bf q}}{{\bf r}})-C ]], this yields

[iq^{2}F({{\bf q}},\Pi) = \textstyle\sum\limits_{k = 1}^{K}{{\bf q}}^{*}{{\hat{\bf n}}}_{k} \int \!\!\int\limits_{\!\!\Gamma_{k}}{\rm d}^{2}r\,\left[\exp({i{{\bf q}}{ {\bf r}}})-C\right].\eqno(24)]

With the notation (8)[link], this proves (21)[link]. With the choice [{{\bf H}}:={{\bf r}}/3], we obtain the volume formula (22)[link]. □

Similar to c in Section 2[link], the constant C can and should be chosen differently for different q domains. At large q, the best choice is C = 0. The small-q case is discussed in Section 3.3[link].

The volume formula (22)[link] has previously been derived by tetrahedral tessellation (Comessatti, 1930[Comessatti, A. A. (1930). Lezioni di Geometria Analitica e Proiettiva I. Padova: Cedam.], Cap. II, §3, III 171).

To see the equivalence of (12)[link] with equation 15 of Croset (2017[Croset, B. (2017). J. Appl. Cryst. 50, 1245-1255.]), we let C = 0, take f(q, Γ) from (15)[link] and use the fact that Ej−1 × Ej is colinear with [{{\hat{\bf n}}}].

3.3. Removable singularity

The closed expression (21)[link] for the polyhedral form factor F(q, Π) contains two removable singularities: the explicit factor q−1, and the factor [q_{\parallel}^{-1}] contained in the polygonal form factors f(q, Γk). For the case that only q, but not q, is close to 0, we rely on the numerically stable computation of f derived in Section 2.4[link]. Here we address the case q → 0.

In analogy to Section 2.4[link], it is sufficient to invoke analyticity to convince ourselves that the singularity of F is removable. The value of F in the limit q → 0 is just the volume of Π. The expansion of (21)[link] starts with

[F({{\bf q}},\Pi) = \sum_{k = 1}^{K}{{\hat{\bf q}}}^{*}{{\hat{\bf n}}}_{k} \left[{{1-C} \over {iq}}{\rm Ar}(\Gamma_{k})+{\cal O}(q^{0})\right].\eqno(25)]

The leading, apparently singular term is identically zero because [\textstyle\sum_{k}{{\hat{\bf n}}}_{k}{\rm Ar}(\Gamma_{k}) = 0]. This, however, holds only in exact arithmetics; in a floating-point implementation, roundoff errors can make the sum nonzero. For q → 0, any such error outgrows all other terms (Fig. 5[link]). Therefore, in the small-q regime, the only sensible choice is C = 1, which lets the leading term vanish.

[Figure 5]
Figure 5
Modulus of the form factor of a truncated tetrahedron (trigonal pyramidal frustum). The base is an equilateral triangle in the xy plane, oriented so that an edge points in the y direction, with edge length L = 1; the dihedral angle is 72°; the height H = L/2. The plot shows [{\rm Re}F] versus q for the off-symmetric direction [{{\hat{\bf q}}} = (1,2,3)/ {14} ^{1/2}]. Blue spots are computed using the analytical expressions (9)[link] and (21)[link]. For qL < 10−6, roundoff errors dominate. The orange line is computed according to (26)[link] with summation (28)[link] up to n = 19.

Unfortunately though, this can lead to roundoff errors in the difference f(q) − f(0) in (21)[link]. As a remedy, we compute the form factor from a series expansion as follows: Combine (21)[link] and (16)[link] to write the form factor as

[F({{\bf q}},\Pi) = 2\textstyle\sum\limits_{k = 1}^{K}{{\hat{\bf q}}}^{*}{{\hat{\bf n}}}_{k }\sum\limits_{j = 1}^{J_{k}}\left[{{\hat{\bf n}}},\widehat{{{\bf q}}_{\parallel}} ^{*},{{\bf E}}_{kj}\right]\phi_{kj}({{\bf q}})\eqno(26)]

with the function

[\phi_{\alpha}({{\bf q}}):={{\tau_{\alpha}({{\bf q}},1)-\tau_{ \alpha}(0,1)} \over {iq}}.\eqno(27)]

The constant τα(0, 1) neutralizes the first term in the expansion (19)[link] of τα(q, 1) so that

[\phi_{\alpha}({{\bf q}}) = \sum_{n = 2}^{\infty}{(iq)}^{n-2}\sum_{l = 0}^{2l\leq n }{{{\left({{\hat{\bf q}}}{{\bf R}}_{\alpha}\right)}^{n-2l}} \over {(n-2l)!}} {{{\left({{\hat{\bf q}}}{{\bf E}}_{\alpha}\right)}^{2l}} \over {(2l+1)!}}.\eqno(28)]

Fig. 5[link] shows that there is good overlap between the domains of the closed expression and the series expansion.

3.4. Polyhedron with inversion center

If a polyhedron has an inversion center at ρ (Schoenflies group Ci), then the form factor has the symmetry F(q, Π) = F(−q, 2ρΠ). As the number of faces is even, we can write K = 2K′. We require that faces numbered k and k + K′ be opposite to each other. We use [{{\hat{\bf n}}}_{k+K^{\prime}} = -{{\hat{\bf n}}}_{k}] to write the form factor as

[F({{\bf q}},\Pi) = \exp(i{\bf q}\,{\boldrho}){{1} \over {iq}} \sum_{k = 1}^{K^{\prime}}{{\hat{\bf q}}}{{\hat{\bf n}}}_{k}\,\tilde{f}({ {\bf q}},\Gamma_{k}-{{\bf \rho}}),\eqno(29)]

where

[\tilde{f}({{\bf q}},\Gamma):= f({{\bf q}},\Gamma)-f({{\bf q}}, -\Gamma)\eqno(30)]

is the form factor of a pair of opposite faces. The symmetry f(q, −Γ) = f(−q, Γ) allows some economy in computing F from the generic closed expression.

In the small-q case, the expansion (26)[link] is symmetrized as

[F({{\bf q}},\Pi) = 2\textstyle\sum\limits_{k = 1}^{K^{\prime}}{{\hat{\bf q}}}^{*}{{\hat{\bf n}}}_{k}\sum\limits_{j = 1}^{J_{k}}\left[{{\hat{\bf n}}},\widehat{{{\bf q}}_ {\parallel}}^{*},{{\bf E}}_{kj}\right]\left[\phi_{kj}({{\bf q}})+\phi_{ kj}(-{{\bf q}})\right],\eqno(31)]

and in consequence in (28)[link] the terms with odd n cancel.

3.5. Prism

For a prism Π = {r | rΓ, |r| < h/2} a much simpler solution is available. We return to the definition (1)[link]. Applying Fubini's theorem to factorize the triple integral from the onset into an integral (8)[link] over the base Γ of the prism and an integral along the normal direction [{{\hat{\bf n}}}], we obtain the form factor

[F({{\bf q}},\Pi) = h\,{\rm sinc}({{\bf q}}_{\perp}{{\hat{\bf n}} }\,h/2)\,f({{\bf q}}_{\parallel},\Gamma_{\parallel})\eqno(32)]

for all q. Thanks to the sinc function in (32)[link], there is no singularity in q and therefore no series expansion is needed for q → 0.

3.6. Asymptotic envelopes

We now come back to the asymptotic power-law envelopes for large q discussed in Section 1.3[link]. A cube Π with side lengths L, centered at the origin and oriented along the coordinate axes, has the form factor

[F({{\bf q}},\Pi) = L^{3}\,{\rm sinc}(q_{x}L/2)\,{\rm sinc}(q_{y }L/2)\,{\rm sinc}(q_{z}L/2).\eqno(33)]

For large q it has the asymptotic envelope |F| ≤ 8/|qxqyqz|, which goes as q−3 for fixed direction [{{\hat{\bf q}}}], provided none of the three components [\hat{q}_{x}], [\hat{q}_{y}], [\hat{q}_{z}] is zero. If [{{\hat{\bf q}}}] is perpendicular to one of the edges of the cube, then one of the three sinc functions has the fixed argument 0 and value 1. And if [{{\hat{\bf q}}}] is perpendicular to one of the faces of the cube, then (33)[link] has two constant factors [{\rm sinc}(0) = 1]. As Croset (2017[Croset, B. (2017). J. Appl. Cryst. 50, 1245-1255.]) has worked out, these observations can be generalized to any polygon. Within our present formalism, this can be confirmed as follows.

For q ≠ 0, the form factor (21)[link] of any K-hedron Π is limited by

[\left|{F({{\bf q}},\Pi)}\right|\leq({{1} / {q}})K\max_{k}|{f({{\bf q}}, \Gamma_{k})}|.\eqno(34)]

For q ≠ 0, the form factor (9)[link] of any J-gon Γ is limited by

[\left|{f({{\bf q}},\Gamma)}\right|\leq({{2} / {q_{\parallel}}})\,J\max_{j}| {\rm sinc}({{\bf q}}_{\parallel}{{\bf E}}_{j})|.\eqno(35)]

For qEj ≠ 0, the sinc function in (35)[link] is limited by

[\left|{\rm sinc}({{\bf q}}{{\bf E}}_{j})\right|\leq{{1} \over {{ {\bf q}}{{\bf E}}_{j}}}.\eqno(36)]

So if all the above conditions are fulfilled, then the polyhedron form factor F has an asymptotic envelope ∼q−3. If there is any edge perpendicular to q, then (36)[link] is not applicable, and [\max|{\rm sinc}({{\bf q}}{{\bf E}}_{j})|] takes the q-independent value 1, so that the envelope of F goes with q−2. If there is any face perpendicular to q, then (35)[link] is not applicable, and [\max|\,f\,|] in (34)[link] takes the q-independent value [{\rm Ar}(\Gamma)]; so the envelope of F goes with q−1.

4. Concluding remarks

4.1. Implementation

Code for computing the form factor of any polygon or polyhedron, based on all the above, has been implemented as part of the open-source GISAS simulation package Born­Again (Pospelov et al., 2020[Pospelov, G., Van Herck, W., Burle, J., Carmona Loaiza, J. M., Durniak, C., Fisher, J. M., Ganeva, M., Yurov, D. & Wuttke, J. (2020). J. Appl. Cryst. 53, 262-276.]). All floating-point numbers, internal and external, have double precision. A summary of the algorithm is given in Appendix C[link].

The code underwent extensive tests for internal consistency and for compatibility with conventional form factor formulae. Checks of BornAgain against the reference code IsGISAXS (Renaud et al., 2009[Renaud, G., Lazzari, R. & Leroy, F. (2009). Surf. Sci. Rep. 64, 255-380.]; Lazzari, 2006[Lazzari, R. (2006). IsGISAXS, Version 2.6 of 4 May 2006, https://www.insp.jussieu.fr/oxydes/IsGISAXS/isgisaxs.htm.]) have been documented by Pospelov et al. (2020[Pospelov, G., Van Herck, W., Burle, J., Carmona Loaiza, J. M., Durniak, C., Fisher, J. M., Ganeva, M., Yurov, D. & Wuttke, J. (2020). J. Appl. Cryst. 53, 262-276.]). In the following, we describe form factor consistency checks that have been permanently added to the BornAgain unit tests.

4.2. Tests

The internal consistency tests address symmetry, specialization and continuity. Symmetry tests are performed for particle shapes that are invariant under some rotation or reflection R. For a suite of wavevectors q, it is checked that the relative deviation of form factors F(q) and F(Rq) stays below a given bound.

The specialization tests address pairs of figures Π1, Π2 with different topologies that coincide for certain geometry parameters. For instance, if Π1 is a box with side lengths a, b, c, and Π2 a truncated cube with side length a and truncation length t, then the choices a = b = c and t = 0 reduce Π1 and Π2 to the same cube. For a suite of wavevectors q, it is checked that the relative deviation of form factors F(q, Π1) and F(q, Π2) stays below a given bound.

The continuity tests search for possible discontinuities due to a change in the computational method. They need special instrumentation of the code, activated through a CMake option and a precompiler macro. Under this option, additional variables tell us whether the analytical expression or the series expansion has been used in the latest form factor computation, and, if applicable, at which expansion order the summation was terminated. For a given direction [{{\hat{\bf q}}}], bisection is used to determine wavevectors where one of these variables changes. Then, the form factor F is computed for wavevectors slightly before and slightly after the transition, and it is checked that the relative step in F remains below a given bound.

All these tests are performed for a suite of particle shapes, for different wavevector directions [{{\hat{\bf q}}}] with different degrees of symmetry, for a logarithmically wide range of magnitudes q and for a range of complex phases.

4.3. Crossover metaparameters

For large q, the polyhedral form factor is computed from (21)[link] with C = 0. For small q, we use (26)[link] with the expansion (28)[link]. Therefore, we need a heuristic metaparameter [\tilde{q}] that determines which algorithm to use. For large q, it still can happen that q is small. Therefore, a second metaparameter [\tilde{q}_{\parallel}] is needed to determine whether face form factors are computed from the closed expression (9)[link] or from (16)[link] with the expansion (19)[link]. As [\tilde{q}] and [\tilde{q}_{\parallel}] are dimensionless, the choice of algorithm is based on qr and qr, where r is the radius of the circumscribing sphere of figure Π. Under a multitude of tests, we obtained the best results with [\tilde{q} = \tilde{q}_{\parallel}\simeq 10^{-2}].

4.4. Accuracy

Currently, the bounds for maximum relative form factor discrepancies are 10−11 in symmetry tests, 6 × 10−12 in specialization tests and 6 × 10−9 in continuity tests. Discrepancies reaching the order of magnitude of these bounds are only observed for a few out of hundreds of thousands of test cases. Most often, errors are smaller than 10−15, i.e. a small multiple of the machine precision. Some of the larger discrepancies are compiler or processor dependent.

The cases of relatively large discrepancy that we have investigated so far all involve special wavevectors that make the integral (1)[link] more symmetric than the underlying figure Π. Appendix D[link] presents one such case: a pyramid that acquires the inversion symmetry of a bipyramid if q lies in the base plane.

It remains to be seen whether such cases warrant closer attention and improved code. So far, we have not encountered a single q, Π combination where symmetry, specialization or continuity tests revealed numeric errors larger than single-precision machine error.

APPENDIX A

Polygon form factor in the literature

As a complement to Section 2.3[link], we demonstrate the equivalence of our form factor (9)[link] with the result (15)[link] of Lee & Mittra (1983[Lee, S. W. & Mittra, R. (1983). IEEE Trans. Antennas Propagat. 31, 99-103.]) and Croset (2017[Croset, B. (2017). J. Appl. Cryst. 50, 1245-1255.]). We start from (9)[link], choose c = 0 and expand the unit in-plane vector [\widehat{{{\bf q}}_{\parallel}}]:

[f({{\bf q}},\Gamma) = {{2} \over {iq_{\parallel}^{2}}}\sum_{j = 1}^{J}[{{\hat{\bf n}}},{{\bf q}}_{\parallel}^{*},{{\bf E}}_{j}]\,{\rm sinc}({ {\bf q}}{{\bf E}}_{j})\exp(i{\bf q}{\bf R}_{j}).\eqno(37)]

Expanding the sinc function, and using (2)[link] and (3)[link] to simplify [{{\bf R}}_{j}\pm{{\bf E}}_{j}], we obtain

[\eqalignno{ q_{\parallel}^{2}\,f({{\bf q}},\Gamma) & = \sum_{j = 1}^{J}{{[{ {\hat{\bf n}}},{{\bf q}}_{\parallel}^{*},{{\bf E}}_{j}]} \over {{{\bf q}} {{\bf E}}_{j}}}\exp(i{{\bf q}}{{\bf V}}_{j-1})\cr &\quad-\sum_{j = 1}^{J} {{[{{\hat{\bf n}}},{{\bf q}}_{\parallel}^{*},{{\bf E}}_{j}]} \over {{ {\bf q}}{{\bf E}}_{j}}}\exp(i{{\bf q}}{{\bf V})_{j}}\cr & = \sum_{j = 1}^{J}\left({{[{{\hat{\bf n}}},{{\bf q}}_{ \parallel}^{*},{{\bf E}}_{j}]} \over {{{\bf q}}{{\bf E}}_{j}}} -{{[{ {\hat{\bf n}}},{{\bf q}}_{\parallel}^{*},{{\bf E}}_{j-1}]} \over {{{\bf q }}{{\bf E}}_{j-1}}}\right) \exp({i{{\bf q}}{{\bf V}}_{j-1}}) \cr & = \left({{\hat{\bf n}}}\times{{\bf q}}_{\parallel}^{*} \right)\cdot\sum_{j = 1}^{J}\left({{{{\bf E}}_{j}} \over {{{\bf q}}{{\bf E}}_{j}}}-{{{{\bf E}}_{j-1}} \over {{{\bf q}}{{\bf E}}_{j-1}}}\right)\exp({i{{\bf q}}{{\bf V}}_{j-1}}).&(38)}]

Using [{{\hat{\bf n}}}{{\bf E}} = 0], and hence qE = qE,

[\eqalignno{ q_{\parallel}^{2}\,f({{\bf q}},\Gamma) & = \left({{\hat{\bf n}} }\times{{\bf q}}_{\parallel}^{*}\right)\cdot\sum_{j = 1}^{J}{{{{\bf E }}_{j}({{\bf q}}_{\parallel}{{\bf E}}_{j-1})-{{\bf E}}_{j-1}({\bf {q}}_{\parallel}{{\bf E}}_{j})} \over {({{\bf q}}{{\bf E}}_{j-1})({{\bf q }}{{\bf E}}_{j})}}\exp({i{{\bf q}}{{\bf V}}_{j-1}})\cr & = \left({{\hat{\bf n}}}\times{{\bf q}}_{\parallel}^{*} \right)\cdot\sum_{j = 1}^{J}{{({{\bf E}}_{j-1}\times{{\bf E}}_{j}) \times{{\bf q}}_{\parallel}} \over {({{\bf q}}{{\bf E}}_{j-1})({{\bf q}}{ {\bf E}}_{j})}}\exp({i{{\bf q}}{{\bf V}}_{j-1}})\cr & = \left[{{\bf q}}_{\parallel}\times\left({{\hat{\bf n}}} \times{{\bf q}}_{\parallel}^{*}\right)\right]\cdot\sum_{j = 1}^{J}{{{ {\bf E}}_{j-1}\times{{\bf E}}_{j}} \over {({{\bf q}}{{\bf E}}_{j-1})({ {\bf q}}{{\bf E}}_{j})}}\exp({i{{\bf q}}{{\bf V}}_{j-1}}).\cr & &(39)}]

Using [{{\hat{\bf n}}}{{\bf q}}_{\parallel} = 0], we obtain (15)[link].

APPENDIX B

Polygon area as lowest expansion coefficient

As a complement to Section 2.4[link], we demonstrate the equivalence of two different expressions for the area of a polygon: the [q_{\parallel}\!\to\!0] limit of the generic form factor (9)[link], and the surveyor's formula (10)[link]. For any constant direction [{{\hat{\bf q}}}] with [q_{\parallel}\!\neq\!0], we have from (16)[link] and (19)[link]

[f(0,\Gamma) = {{2} \over {q_{\parallel}}}\left({{\hat{\bf n}}}\times\widehat{{ {\bf q}}_{\parallel}}^{*}\right)\cdot\sum_{j = 1}^{J}{{\bf E}}_{j}({ {\bf q}}{{\bf R}}_{j}).\eqno(40)]

Splitting q = q + q, drawing the constant qr in front of the sum and using [\textstyle\sum{{\bf E}}_{j} = 0] we obtain

[f(0,\Gamma) = 2\left({{\hat{\bf n}}}\times\widehat{{{\bf q}}_{\parallel}}^ {*}\right)\cdot\textstyle\sum\limits_{j = 1}^{J}{{\bf E}}_{j}(\widehat{{{\bf q}}_{ \parallel}}\,{{\bf R}}_{j}).\eqno(41)]

Inserting the definitions of Ej and Rj,

[\phi_{0}(\widehat{{{\bf q}}_{\parallel}},\Gamma) = {{{{\hat{\bf n}}} \times\widehat{{{\bf q}}_{\parallel}}^{*}} \over {2}}\cdot\sum_{j = 1}^{J}\left({ {\bf V}}_{j}-{{\bf V}}_{j-1}\right)\left[\widehat{{{\bf q}}_{ \parallel}}\cdot({{\bf V}}_{j}+{{\bf V}}_{j-1})\right].\eqno(42)]

Multiplying out and shuffling indices [(j\!-\!1)\to j] for some terms under the sum gives

[\eqalignno{\phi_{0}(\widehat{{{\bf q}}_{\parallel}},\Gamma) & = {{{ {\hat{\bf n}}}\times\widehat{{{\bf q}}_{\parallel}}^{*}} \over {2}}\cdot\sum_{j = 1}^{J}\left[{{\bf V}}_{j}(\widehat{{{\bf q}}_{\parallel}}{{\bf V}}_{j -1})-{{\bf V}}_{j-1}(\widehat{{{\bf q}}_{\parallel}}{{\bf V}}_{j})\right] \cr & = {{{{\hat{\bf n}}}\times\widehat{{{\bf q}}_{\parallel}} ^{*}} \over {2}}\cdot\sum_{j = 1}^{J}\left({{\bf V}}_{j-1}\times{{\bf V}}_{j} \right)\times\widehat{{{\bf q}}_{\parallel}} \cr & = {{\widehat{{{\bf q}}_{\parallel}}\times({{\hat{\bf n}} }\times\widehat{{{\bf q}}_{\parallel}}^{*})} \over {2}}\cdot\sum_{j = 1}^{J}\left({ {\bf V}}_{j-1}\times{{\bf V}}_{j}\right).&(43)}]

Using [{{\hat{\bf n}}}{{\bf q}}_{\parallel} = 0], we recover (10)[link].

APPENDIX C

Algorithm summary

In this appendix, we summarize the algorithm for the computation of the form factor F(q; Π) as derived in this work. For clarity and brevity, we only consider a polyhedron Π without any inversion symmetry. For polyhedra with inversion symmetry, and for other details omitted here, see the actual implementation in the open-source code BornAgain.

Read the wavevector q, the topology T, the vertex coordinates C and the symmetry flags (Section 3.1[link], ignored here). Discard faces with zero or negligible area. For each face, merge adjacent vertices with zero or negligible distance. Assert that all remaining faces are planar. Compute the circumscribing radius r of Π (Section 4.3[link]). For each face Γk, compute the circumscribing rk.

If qr < 10−2 (Section 4.3[link]), then compute F according to (26)[link] and (28)[link], with changed summation order so that the outermost summation runs over the expansion power n. Terminate if the relative contribution of two subsequent expansion terms is less than ε = 2 × 10−16.

Otherwise ([qr\geq\tilde{q}]), compute F according to (21)[link] with C = 0. This is a weighted sum of polygonal face form factors f(q, Γk). They are computed as follows:

If qrk < 10−2, then compute f according to (16)[link] and (19)[link] with changed summation order so that the outermost summation runs over the expansion power n. Terminate if the relative contribution of two subsequent expansion terms is less than ε = 2 × 10−16.

Otherwise (qrk ≥ 10−2), compute f according to (9)[link] with c = 0.

APPENDIX D

Additional symmetry at special wavevectors

As a complement to Section 4.4[link], we give one example of how special wavevectors cause extra symmetries in the integral (1)[link] that lead to cancellation and roundoff errors. Consider a pyramid Π with a base Γ0 that has at least a twofold rotation axis, as considered in Section 2.5[link]. Let [{{\hat{\bf n}}}] be the normal of Γ0, pointing towards the outside of Π. With (4)[link] and (5)[link] it induces the decomposition r = r + r.

Now consider a wavevector q in the plane of Γ0 so that q = 0. To see how this causes an extra symmetry in the integral (1)[link], consider the bipyramid Π2 that is the union of Π and the mirror image of Π under reflection about the Γ0 plane,

[\Pi_{2} = \{{{\bf r}}\,|\,{{\bf r}}_{\parallel}+{{\bf r}}_{\perp}\in\Pi \lor{{\bf r}}_{\parallel}-{{\bf r}}_{\perp}\in\Pi\}.\eqno(44)]

Obviously, for the special in-plane q under consideration,

[F({{\bf q}},\Pi) = \textstyle{{1} \over {2}}F({{\bf q}},\Pi_{2}).\eqno(45)]

Thereby, the form factor integral (1)[link] for F(q, Π) has the extra symmetries of the bipyramid Π2. The combination of the S2 symmetry assumed for Γ0 and the mirror symmetry that constitutes the bipyramid implies an inversion symmetry around the center of Γ0. As discussed in Section 3.4[link] this leads to cancellation of some terms. Unless these terms are removed from the implemented formulae, severe roundoff errors must be expected.

Acknowledgements

I thank Céline Durniak, Walter Van Herck and Gennady Pospelov for reference code and for help with the test framework. I thank Bernard Croset for helpful suggestions and for confirming the sign of (15)[link]. Open access funding enabled and organized by Projekt DEAL.

References

First citationBernadini, F. (1991). Comput. Aided Des. 23, 51–58.  Google Scholar
First citationBraden, B. (1986). Collect. Math. J. 17, 326–337.  Google Scholar
First citationCattani, C. & Paoluzzi, A. (1990). Comput. Aided Des. 22, 130–135.  Google Scholar
First citationChourou, S. T., Sarje, A., Li, X. S., Chan, E. R. & Hexemer, A. (2013). J. Appl. Cryst. 46, 1781–1795.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationComessatti, A. A. (1930). Lezioni di Geometria Analitica e Proiettiva I. Padova: Cedam.  Google Scholar
First citationCroset, B. (2017). J. Appl. Cryst. 50, 1245–1255.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationCroset, B. (2018). J. Appl. Cryst. 51, 1005–1012.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationEngel, K. & Laasch, B. (2020). arxiv:2011.06971v1 [Math-ph].  Google Scholar
First citationHammouda, B. (2010). Probing Nanoscale Structures – The SANS Toolbox, https://www.ncnr.nist.gov/staff/hammouda/the_SANS_toolbox.pdfGoogle Scholar
First citationHenry, C. R. (2005). Prog. Surf. Sci. 80, 92–116.  CAS Google Scholar
First citationLazzari, R. (2006). IsGISAXS, Version 2.6 of 4 May 2006, https://www.insp.jussieu.fr/oxydes/IsGISAXS/isgisaxs.htmGoogle Scholar
First citationLee, S. W. & Mittra, R. (1983). IEEE Trans. Antennas Propagat. 31, 99–103.  Google Scholar
First citationLi, S. Z. & Xie, Z. Q. (2020). Appl. Math. Comput. 376, 121540.  Google Scholar
First citationLi, X., Shew, C.-Y., He, L., Meilleur, F., Myles, D. A. A., Liu, E., Zhang, Y., Smith, G. S., Herwig, K. W., Pynn, R. & Chen, W.-R. (2011). J. Appl. Cryst. 44, 545–557.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationLien, S. & Kajiya, J. T. (1984). IEEE Comput. Graph. 4, 35–42.  Google Scholar
First citationMirtich, B. (1996). J. Graphics Tools 1, 31–50.  Google Scholar
First citationPospelov, G., Van Herck, W., Burle, J., Carmona Loaiza, J. M., Durniak, C., Fisher, J. M., Ganeva, M., Yurov, D. & Wuttke, J. (2020). J. Appl. Cryst. 53, 262–276.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationRayleigh, Lord (1881). Philos. Mag. 12, 81–101.  Google Scholar
First citationRenaud, G., Lazzari, R. & Leroy, F. (2009). Surf. Sci. Rep. 64, 255–380.  Web of Science CrossRef CAS Google Scholar
First citationSenesi, A. & Lee, B. (2015). J. Appl. Cryst. 48, 565–577.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationWuttke, J. (2017). arxiv:1703.00255.  Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

Journal logoJOURNAL OF
APPLIED
CRYSTALLOGRAPHY
ISSN: 1600-5767
Follow J. Appl. Cryst.
Sign up for e-alerts
Follow J. Appl. Cryst. on Twitter
Follow us on facebook
Sign up for RSS feeds