The analytical computation of the full gravity tensor from a polyhedral source of homogeneous density is presented, with emphasis on its algorithmic implementation. The theoretical development is based on the subsequent transition of the general expressions from volume to surface and from surface to line integrals, defined along the closed polygons building each polyhedral face. However, the accurate numerical computation of the obtained transcendental expressions is linked with the relative position of the computation point and its corresponding projections on the plane of each face and on the line of each segment with respect to the polygons defining each face. Depending on this geometric setup, the application of the divergence theorem of Gauss leads to the appearance of additional correction terms, valid only for these boundary conditions and crucial for the correct numerical evaluation of the polyhedral-related gravity quantities at those locations of the computation point. A program in Fortran is supplied and thoroughly documented; it computes the gravitational potential, its first-order derivatives, and the full gradiometric tensor at arbitrary space points due to a general polyhedral source of constant density.

The study of the gravity field of arbitrary polyhedral bodies of homogeneous density has provoked a series of publications over the last decades. Some of the researchers represented an arbitrary three dimensional body in terms of contours obtained by the intersection of horizontal planes with the body. They replaced each contour by a horizontal irregular n-sided polygonal lamina and chose the combination of the analytical solution for the attraction of an arbitrary polygonal lamina (Jung, 1961; Talwani, 1973) with a numerical integration along the z-direction for the computation of the body's gravity field (Talwani and Ewing, 1960; Hellinger, 1983). Others devoted their efforts to the derivation of purely analytical expressions, based on the technique of reducing the volume Newtonian integral to a summation of line integrals along the segments of each facet, applying Gauss' divergence theorem twice: for the conversion first to surface and then to line integrals. The pioneering work in this field came from Paul (1974), Plouff (1976), and Barnett (1976). Analytical formulas for the potential due to a homogeneous polyhedral body of arbitrary shape and its derivatives up to second order were presented in a compact manner, suitable for efficient computer programming, by Petrović (1996) and Werner and Scheeres (1997). Optimized analytical formulas for the gravity field of a homogeneous polyhedral body, appropriate for efficient calculations, have been presented by Pohanka (1988). An error analysis and further optimization of these expressions were undertaken by Holstein and Ketteridge (1996) and Holstein et al. (1999). Recently, the research in this field was extended to the nonhomogeneous case (Pohanka, 1998; Hansen, 1999). The goal of this research note is to draw attention to the singularities appearing in the numerical evaluation of the closed expressions for the gravitational field of a homogeneous polyhedron. These appear in the transition from surface to line integrals when the orthogonal projection of the computation point onto the plane defined by a facet of the polyhedron lies inside the polygon defining the facet or falls on the line segment describing the edge. Indeed, certain singularities arise when one of these cases occur. We will show, however, that the gravity field is nevertheless defined, if one treats the singularities in the usual way of potential field theory, excluding a small sphere or circle around the singular point. Thus, certain correction terms arise, which should be taken into account when the respective singularities occur.