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: [email protected]

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 Mathematical equation is

Mathematical equation

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

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 Mathematical equation; 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

Mathematical equation

and mid-to-end vector

Mathematical equation

The normal vector Mathematical equation of the plane spanned by Vj shall be oriented such that ∂Γ has the winding number +1 (fulfills the right-hand rule with respect to Mathematical equation). The oriented plane characterized by Mathematical equation induces a decomposition of any vector Mathematical equation into a component perpendicular to the plane,

Mathematical equation

and an in-plane component,

Mathematical equation

This decomposition will be applied to position vectors r and to wavevectors q. The oriented plane is fully specified by its normal vector Mathematical equation and its distance from the origin, r.

Complex conjugation is denoted by a superscript asterisk. The absolute value of a complex vector is written Mathematical equation.

Note that the in-plane unit vector

Mathematical equation

differs from the in-plane component Mathematical equation of the unit vector Mathematical equation. In this work, we shall only use Mathematical equation and Mathematical equation, not Mathematical equation.

The triple product is denoted

Mathematical equation

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 Mathematical equation has the analytic continuation Mathematical equation. 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

Mathematical equation

Proposition

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

Mathematical equation

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 Γ,

Mathematical equation

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

Mathematical equation

Proof

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

Mathematical equation

To prove (9)[link], we choose Mathematical equation. The left-hand side of (12)[link] is

Mathematical equation

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

Mathematical equation

Each edge can be written as a parametric curve Mathematical equation so that

Mathematical equation

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 Mathematical equation. □

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 Mathematical equation should be swapped). In our notation, it reads

Mathematical equation

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

[Figure 1]
Figure 1
Modulus of the form factor of an equilateral triangle as function of wavenumber q for three different wavevector directions Mathematical equation. 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 Mathematical equation. 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

Mathematical equation

with the function

Mathematical equation

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

Mathematical equation

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 Mathematical equation. 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 Mathematical equation, 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

Mathematical equation

As a consistency check, we note that the limit Mathematical equation, at constant in-plane direction Mathematical equation, is the n = 1 term Mathematical equation. 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 Mathematical equation. 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 Mathematical equation and Mathematical equation. For q ≠ 0, the form factor is

Mathematical equation

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 Mathematical equation shall point to the outside of Π; this then determines the order of the vertices in the sequence Mathematical equation.

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

Mathematical equation

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

Mathematical equation

Proof

For a polyhedron, the divergence theorem takes the form

Mathematical equation

With the choice Mathematical equation, this yields

Mathematical equation

With the notation (8)[link], this proves (21)[link]. With the choice Mathematical equation, 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 Mathematical equation.

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 Mathematical equation 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

Mathematical equation

The leading, apparently singular term is identically zero because Mathematical equation. 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 Mathematical equation versus q for the off-symmetric direction Mathematical equation. 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

Mathematical equation

with the function

Mathematical equation

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

Mathematical equation

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 Mathematical equation to write the form factor as

Mathematical equation

where

Mathematical equation

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

Mathematical equation

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 Mathematical equation, we obtain the form factor

Mathematical equation

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

Mathematical equation

For large q it has the asymptotic envelope |F| ≤ 8/|qxqyqz|, which goes as q−3 for fixed direction Mathematical equation, provided none of the three components Mathematical equation, Mathematical equation, Mathematical equation is zero. If Mathematical equation 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 Mathematical equation is perpendicular to one of the faces of the cube, then (33)[link] has two constant factors Mathematical equation. 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

Mathematical equation

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

Mathematical equation

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

Mathematical equation

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 Mathematical equation 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 Mathematical equation in (34)[link] takes the q-independent value Mathematical equation; 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 Mathematical equation, 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 Mathematical equation 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 Mathematical equation that determines which algorithm to use. For large q, it still can happen that q is small. Therefore, a second metaparameter Mathematical equation 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 Mathematical equation and Mathematical equation 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 Mathematical equation.

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 Mathematical equation:

Mathematical equation

Expanding the sinc function, and using (2)[link] and (3)[link] to simplify Mathematical equation, we obtain

Mathematical equation

Using Mathematical equation, and hence qE = qE,

Mathematical equation

Using Mathematical equation, 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 Mathematical equation limit of the generic form factor (9)[link], and the surveyor's formula (10)[link]. For any constant direction Mathematical equation with Mathematical equation, we have from (16)[link] and (19)[link]

Mathematical equation

Splitting q = q + q, drawing the constant qr in front of the sum and using Mathematical equation we obtain

Mathematical equation

Inserting the definitions of Ej and Rj,

Mathematical equation

Multiplying out and shuffling indices Mathematical equation for some terms under the sum gives

Mathematical equation

Using Mathematical equation, 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 (Mathematical equation), 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 Mathematical equation 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,

Mathematical equation

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

Mathematical equation

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