GravDyn: A Python Framework for Gravitational Modeling of Irregular Celestial Bodies Using Polyhedral, Mascon, and Series Expansion Methods

Safwan Aljbaae; Valerio Carruba; Allan K. de Almeida Jr; Gabriel Antonio Caritá; Antonio F. B. A. Prado; Marcelo L. Mota; Carlos E. Ferreira Lopes

doi:10.5334/jors.743

1 Overview

1.1 Introduction

Irregularly shaped small bodies do not settle into hydrostatic equilibrium (1), so their gravity fields differ markedly from the simple point-mass picture. Those differences matter when we design spacecraft trajectories, assess orbital stability, or plan missions that will operate close to an asteroid (2, 3, 4, 5). Several approaches have been developed to address this problem. The classical polyhedral method provides an accurate analytical solution for the gravitational potential of a homogeneous body by transforming volume integrals into surface and line integrals, making it suitable for modeling the gravitational field near the surface of irregular bodies (6, 7). Despite its accuracy, this method is computationally expensive, especially for high-resolution shape models, as the number of faces and edges directly impacts the evaluation cost. Alternatively, a different philosophy replaces the continuous mass with a collection of point masses distributed throughout the volume (8). Building on that, later research formalized mascons using tetrahedral decompositions of a polyhedron. By placing a point mass at each tetrahedron’s centroid, this approach preserves much of the original density field while reducing the computational load (9, 10, 11). In addition, the method naturally allows the incorporation of layered internal structures with varying densities. More recently, the Potential Series Expansion Method (PSEM) has been proposed (12). Here, the potential is expressed as a convergent polynomial series derived analytically from the tetrahedral elements. PSEM derives a polynomial expansion that provides a fast-to-evaluate approximation outside the Brillouin sphere.

Although these formulations are well established, their implementations are often tied to individual studies. For instance, several tools for the polyhedral gravity model are available in different languages such as C++17 with Python interface¹ or Fortran.² This isolation between different methods makes it difficult to reproduce previous calculations, compare methods under the same assumptions, or reuse common shape-processing steps across projects. GravDyn addresses this practical gap by placing polyhedral, mascon, and PSEM gravity models in a single Python package. Each method presents different advantages and limitations depending on the application. GravDyn helps users determine which method is most appropriate for a given application.

The package follows a common workflow. A user provides a triangular shape model and the physical parameters of the body. GravDyn then verifies and rescales the geometry, builds the selected gravity representation, and evaluates the potential or acceleration at prescribed field points. This paper describes the structure of the package, the implemented models, and the validation tests used to compare them.

1.2 Implementation and architecture

GravDyn is implemented in Python and is organized into modular components. The package provides tools for efficient numerical evaluation of gravitational potentials and accelerations in small-body environments. As shown in Figure 1, the system consists of two pillars: a core computational library and an optional graphical user interface (GUI). The core lives in the gravdyn package and contains all numerical routines, data structures, and processing tools. The GUI offers a user-friendly interface for visualization, parameter configuration, and interaction with the underlying models. The GUI is intentionally decoupled from the core computational layer. This separation ensures that the scientific functionality remains independent of the interface, allowing the software to be used both programmatically and interactively.

Structure of the **GravDyn** package. The software is organized into two main modules: shape verification and potential calculation. The potential calculation module includes polyhedral, expansion, and mascon-based approaches.

1.2.1 Shape-model preprocessing

The implementation is essentially based on a functional decomposition which constitutes the standard workflow for gravitational modeling. The first step is geometry preprocessing, which is handled by the shape_verification module. This makes sure that the polyhedral shape model is geometrically and physically consistent before any gravitational analysis is performed. First, the mesh is loaded and the body is aligned to its principal axes of inertia, which is helpful to ensure that the following calculations are consistent with the physical structure of the object. The shape is then rescaled to match some reference volume derived from the specified mass and density and translated so that its center of mass coincides with the origin. Although these steps are standard from a scientific point of view, they are essential, as small inconsistencies at this stage can propagate into the dynamical results. The processed geometry is stored for later use, and a set of diagnostic plots is generated to verify that the preprocessing stage has been performed correctly.

The second stage of the workflow is devoted to building the gravitational model and evaluating the resulting field. In GravDyn, the polyhedral method, the mascon approach, and the PSEM are implemented within a single computational framework with a common interface. This design allows users to switch between gravity models, compare their performance and accuracy under the same conditions, and integrate different modeling strategies into the same dynamical workflow, rather than treating each method as a separate standalone implementation.

1.2.2 Polyhedral gravity model

The classical polyhedral method gives an exact analytical expression for the gravitational potential of a homogeneous body. In GravDyn, this is done by combining the preprocessing routine masconprepare_werner_model with the evaluation function batched_werner_potential. The preprocessing step follows the work of Werner and Scheeres (6), Scheeres (2), to perform calculations characterizing the body, such as face and edge dyads, normal vectors, and other geometric coefficients. The polyhedral formulation generally provides the most accurate results, especially near the surface where the approximations are often not valid.

The package implements the constant-density polyhedron formulation of Werner and Scheeres (6), Scheeres (2). In this model, the body is described by a closed triangular mesh, and the gravitational potential is evaluated as a sum of edge and face contributions. For a field point $r$ , we define $r_{e}$ as the vector from an arbitrary point on edge e to $r$ , and $r_{f}$ as the vector from an arbitrary point on face f to $r$ . Each triangular facet is specified by three ordered vertices, for example $v_{1}$ , $v_{2}$ , and $v_{3}$ , with associated edges E₁₂, E₂₃, and E₃₁, and an outward unit normal ${\hat{n}}_{123}$ . A schematic representation of these geometric quantities is shown in Figure 2.

Redrawn schematic of the geometric entities associated with a triangular facet in the polyhedral gravity model of Werner and Scheeres (6), Scheeres (2). The facet is defined by the ordered vertices $v_{1}$ , $v_{2}$ , and $v_{3}$ , with edges E₁₂, E₂₃, and E₃₁, and outward unit normal ${\hat{n}}_{123}$ . The field point is represented in the body-fixed frame *(x,y,z)*.

The potential can then be written as

1

U (r) = \frac{G σ}{2} [\underset{e \in edges}{Σ} r_{e} E_{e} r_{e} L_{e} - \underset{f \in faces}{Σ} r_{f} F_{f} r_{f} ω_{f}]

The corresponding first and second derivatives are

2

\nabla U (r) = - G σ [\underset{e \in edges}{Σ} E_{e} r_{e} L_{e} - \underset{f \in faces}{Σ} F_{f} r_{f} ω_{f}],

3

\nabla_{r}^{2} U (r) = G σ [\underset{e \in edges}{Σ} E_{e} L_{e} - \underset{f \in faces}{Σ} F_{f} ω_{f}] .

where G is the gravitational constant and $σ$ is the bulk density. The matrices $E_{e}$ and $F_{f}$ are the edge and face dyads, respectively, while L_e is the logarithmic edge factor and $ω_{f}$ is the signed solid angle subtended by face f at the evaluation point. These quantities are defined as

4

E_{e} = {\hat{n}}_{f} {\hat{n}}_{e}^{f} + {\hat{n}}_{f^{’}} {\hat{n}}_{e}^{f^{’}},

5

F_{f} = {\hat{n}}_{f} {\hat{n}}_{f},

6

L_{e} = \ln (\frac{r_{1}^{e} + r_{2}^{e} + e_{e}}{r_{1}^{e} + r_{2}^{e} - e_{e}}),

7

ω_{f} = 2 \arctan (\frac{r_{1}^{f} \cdot r_{2}^{f} \cdot r_{3}^{f}}{r_{1}^{f} r_{2}^{f} r_{3}^{f} + r_{1}^{f} r_{2}^{f} \cdot r_{3}^{f} + r_{2}^{f} r_{3}^{f} \cdot r_{1}^{f} + r_{3}^{f} r_{1}^{f} \cdot r_{2}^{f}}) .

Here f and f’ are the two faces adjacent to edge e, ${\hat{n}}_{f}$ and ${\hat{n}}_{f^{’}}$ are their outward unit normals, and ${\hat{n}}_{e}^{f}$ and ${\hat{n}}_{e}^{f^{’}}$ are the corresponding edge-normal unit vectors lying in each face and pointing outward from the face. The scalar e_e is the length of edge e, while $r_{1}^{e}$ and $r_{2}^{e}$ are the distances from the evaluation point to the two vertices of that edge. Similarly, $r_{1}^{f}$ , $r_{2}^{f}$ , and $r_{3}^{f}$ are the vectors from the evaluation point to the three vertices of face f, with $r_{i}^{f} = ‖ r_{i}^{f} ‖$ .

In the preprocessing stage, GravDyn therefore computes and stores the geometric quantities required by these sums: face normals, edge normals, face dyads, edge dyads, edge lengths, and the connectivity between adjacent faces and edges. During field evaluation, these precomputed quantities are reused so that only the point-dependent terms L_e, $ω_{f}$ , $r_{e}$ , and $r_{f}$ must be evaluated. This is also why the cost of the polyhedral method scales with the number of faces and edges in the mesh. The formulation is especially valuable near the surface because it remains valid down to the boundary of the body, unlike spherical-harmonic expansions, which may diverge close to or inside the circumscribing sphere.

The numerical evaluation of closed-form polyhedral gravity models near singular geometrical configurations was analyzed by Tsoulis and Petrovic (7). Their work complements the formulation of Werner and Scheeres (6) by treating the cases in which the orthogonal projection of the computation point onto a face lies inside the polygon, on an edge, or at a vertex. In such situations, apparent singularities arise in the line-integral expressions, although the potential and first-order derivatives remain well defined when appropriate limiting correction terms are included. Incorporating this treatment would require additional geometrical tests and correction terms during the field evaluation, increasing the computational cost of the method. The present version of GravDyn therefore keeps the Werner and Scheeres (6) edge-face formulation as the baseline polyhedral model, while the explicit Tsoulis and Petrovic (7) singularity treatment is left for a future implementation stage.

1.2.3 Layered mascon model

The mascon gravity model with a polyhedral source is implemented by the functions generate_layered_mascons and batched_pot_mascon. As in Chanut et al. (9), Aljbaae et al. (10, 11), we connect each face of the surface mesh to the body-fixed origin to form a tetrahedron. If the three vertices of a face are denoted by $v_{1}$ , $v_{2}$ , and $v_{3}$ , then the signed volume of the associated tetrahedron is obtained as

8

V_{T} = \frac{1}{6} v_{1} \cdot (v_{2} \times v_{3}) .

The tetrahedron is then subdivided radially according to the number of density layers provided by the user. For a model with $N_{ℓ}$ layers, the scaled vertices of the k-th radial boundary are written as

9

v_{j, k} = λ_{k} v_{j}, λ_{k} = \frac{k}{N_{ℓ}}, j = 1, 2, 3 .

The cumulative volume enclosed by this scaled tetrahedron is

10

V_{T, k} = \frac{1}{6} v_{1, k} \cdot (v_{2, k} \times v_{3, k}),

so that the volume assigned to layer k is

11

Δ V_{T, k} = {\begin{matrix} V_{T, 1}, & if k = 1, \\ V_{T, k} - V_{T, k - 1}, & if k > 1 . \end{matrix}

The preliminary mass and gravitational parameter of each mascon are then computed from the layer density as

12

m_{T, k}^{0} = ρ_{k} Δ V_{T, k}, μ_{T, k}^{0} = G m_{T, k}^{0} .

The mascon position is assigned using the same approximation adopted in the original layered implementation. For the innermost layer, the point mass is placed at

13

c_{T, 1} = \frac{v_{1, 1} + v_{2, 1} + v_{3, 1}}{4},

while for the outer layers it is placed at the average position of the three vertices of the outer and inner layer boundaries,

14

c_{T, k} = \frac{v_{1, k} + v_{2, k} + v_{3, k} + v_{1, k - 1} + v_{2, k - 1} + v_{3, k - 1}}{6}, k > 1 .

After all mascons have been generated, the total preliminary mass is compared with the prescribed asteroid mass M. If necessary, all masses are rescaled by a common factor,

15

α = \frac{M}{Σ_{i} m_{i}^{0}}, m_{i} = α m_{i}^{0}, μ_{i} = G m_{i} .

This step conserves the total mass of the body,

16

\underset{i}{Σ} m_{i} = M,

while preserving the relative mass distribution imposed by the selected density layers. The potential at a field point $r$ is then approximated as

17

U (r) = Σ_{i = 1}^{N_{m}} \frac{μ_{i}}{r_{i}}, r_{i} = ‖ r - c_{i} ‖,

where N_m is the total number of mascons. The function batched_pot_mascon evaluates these point-mass contributions for many field points simultaneously, making the method suitable for large grids and long orbital integrations. The main advantage of this implementation is that it keeps the geometry tied to the original polyhedral shape, supports layered internal density models, and preserves the prescribed total mass.

1.2.4 Potential Series Expansion Method

The PSEM is based on the analytical expansion of the Newtonian potential kernel $1 / r$ into a convergent power series of Legendre polynomials. Considering a field point $r$ and a source point inside the body, the inverse distance can be written as

18

\frac{1}{r} = \frac{1}{ρ} Σ_{i = 0}^{\infty} P_{i} (u) {(\frac{ρ^{’}}{ρ})}^{i},

where $ρ = ‖ r ‖$ , $ρ^{’}$ is the distance from the origin to the mass element, and $u = \cos γ$ , with $γ$ being the angle between the vectors $r$ (from the origin to the field point) and $r^{’}$ (from the origin to the mass element). This angle satisfies

19

\cos γ = \frac{r \cdot r^{’}}{‖ r ‖ ‖ r^{’} ‖} .

This expansion is uniformly convergent for points located outside the Brillouin sphere of the body.

In PSEM, the irregular body is first represented as a polyhedron and decomposed into tetrahedral elements, similarly to the classical polyhedral method. The total gravitational potential is then obtained by integrating the series term-by-term over each tetrahedron and summing all contributions, yielding

20

U (r) = G σ Σ_{k = 1}^{N_{T}} Σ_{i = 0}^{m} ∭_{Q_{k}} P_{i} (u) \frac{ρ^{’, i}}{ρ^{i + 1}} 𝑑 V,

where Q_k denotes each tetrahedral element, N_T is the total number of tetrahedra, and m is the truncation order of the expansion. As shown in Mota et al. (12), each term of the integrand becomes a polynomial function of the coordinates, allowing the volume integrals to be evaluated analytically using a coordinate transformation to a reference tetrahedron. The resulting potential can be expressed as a finite polynomial expansion in Cartesian coordinates of the form

21

U (r) = Σ_{l = 0}^{m} \underset{i + j + k = l}{Σ} U_{i, j, k}^{l} x^{i} y^{j} z^{k} r^{- (2 l + 1) / 2} \forall i, j, k, l \in N : i + j + k = l,

which enables fast evaluation of the potential and its derivatives.

In GravDyn, the function build_potential_derivatives is used to obtain the PSEM coefficients. Deriving these coefficients symbolically becomes computationally expensive for high truncation orders and high-resolution shape models. The computational cost depends strongly on the number of faces in the polyhedral shape model, the corresponding number of tetrahedral elements, and the selected truncation order. For some high-resolution asteroid models, the coefficient-generation stage may require several days of computation.

For this reason, the current release of GravDyn provides PSEM as a runtime evaluation tool based on precomputed and validated coefficient files. These polynomial expansions are generated offline on dedicated servers and stored in the GravDyn repository for selected bodies. At runtime, build_potential_derivatives loads the precomputed expansion for the selected asteroid and returns both the symbolic expressions and optimized numerical callables, enabling very fast evaluation of the potential and its derivatives from Eq. (21).

If a precomputed model is not yet available, the function notifies the user and provides instructions to request its generation. The generation of new coefficient sets is currently maintained as an offline preprocessing workflow by the authors, rather than as a standard user-facing routine, because each new model requires substantial symbolic computation and validation against reference solutions. Once generated and validated, the corresponding model can be added to the repository for future reuse by the community.

1.3 Quality control

1.3.1 Automated tests

GravDyn includes automated tests written with pytest. These cover shape preprocessing, polyhedral precomputation, mascon generation, PSEM model loading, plotting utilities, and point-mass evaluations against analytical solutions. The aim of the test suite is to verify that the implemented routines return consistent arrays, reuse cached data correctly, and reproduce reference values for representative cases.

1.3.2 Input checks

Input validation is applied before field evaluation. The package checks array dimensions, evaluation-point shapes, mass and density inputs, and the availability of external data files, and it raises informative errors when a configuration is inconsistent.

1.3.3 Validation and performance

GravDyn is validated through unit-level tests and also by comparing the gravitational modeling approaches against the classical polyhedral solution, which is taken as the reference. The mascon and PSEM implementations are evaluated at varying distances from the body.

Figure 3 compares the mascon and PSEM models against the polyhedral reference as a function of radial distance. Both approximations converge toward the reference in the external region, while Table 1 shows the expected cost trade-off: PSEM is fastest once its coefficients have been precomputed, whereas the mascon model remains useful when layered interiors must be represented explicitly. The results shown in Figure 3 and Table 1 refer specifically to asteroid (99942) Apophis, which is used here as an illustrative benchmark. Similar comparisons for other bodies, including (21) Lutetia and (87) Sylvia, have been performed in previous work (12). The row labeled “Tsoulis” is included only as an external reference benchmark. It does not correspond to a released gravdyn module, but to a standalone Fortran implementation of the polyhedral gravity formulation following Tsoulis and Petrović (7). This implementation was not written using the vectorized, parallel, or JAX-based execution paths used in the Python routines. It is therefore not intended as a direct software-performance comparison with the released package, but as a reference point for a well-established and robust polyhedral formulation.

Example comparison of relative error in the gravitational potential as a function of radial distance r. The polyhedral method is used as a reference solution. The expansion method (green) and mascon approach (blue) show convergence toward the reference as the distance increases. The vertical red line indicates the approximate radius of the body.

Table 1

Performance comparison of gravitational modeling methods for asteroid (99942) Apophis. The test consists of evaluating the gravitational potential at 40,000 field points using a shape model with 3,996 faces. Timings were obtained on a system equipped with an Intel(R) Core(TM) i7-14700 processor (20 physical cores, 28 logical CPUs), with a maximum clock frequency of 5.4 GHz and 33 MiB of L3 cache.

METHOD	EXECUTION TIME (s)	RELATIVE SPEED	MEAN RELATIVE ERROR (%)
Tsoulis	155.21	$\times 1$	Reference
Werner	9.96	$\times 15.6$	Reference
Mascon (30 layers)	20.14	$\times 7.7$	$6.86 \times 10^{- 5}$
PSEM	1.50	$\times 103.5$	$1.84 \times 10^{- 4}$

The “Mean Relative Error” reported in Table 1 is the average of the pointwise relative potential errors shown as a function of radial distance in Figure 3, computed with respect to the polyhedral reference solution over the 40,000 evaluation points.

The PSEM evaluation is essentially the direct computation of an analytical polynomial, which makes it very fast and inexpensive to evaluate, even when higher orders are required. The mascon approach, by contrast, involves summing the contributions of a large number of discrete masses; this becomes increasingly demanding as the number of layers grows, particularly when a refined internal structure is modeled. Very close to the surface, the mascon formulation can behave more steadily provided the discretization is sufficiently fine, whereas PSEM becomes less reliable due to truncation effects and its formal validity being restricted to the region outside the Brillouin sphere.³ Even so, both methods reach a level of accuracy that is more than adequate from a dynamical perspective: the residual errors are small enough that they do not meaningfully affect the qualitative or quantitative evolution of orbital trajectories in practical short-term trajectory integrations.

1.3.4 Reproducibility and execution environment

For general use, no special hardware is required. For computationally intensive applications, in particular for generating mascon models with multiple layered internal structure, it is advantageous to have a system with sufficient memory (suggested $\geq$ 16 GB RAM) and multi-core CPU support. The generation of layered mascon distributions involves processing large polyhedral meshes and assigning mass elements for multiple density layers, which can greatly increase memory usage and computational time. The required disk space depends on the size of the polyhedral shape model and associated data sets. Typically, less than 1 GB of space is required, but larger shape models and pre-computed data sets may require more.

For typical use, GravDyn does not require specialized hardware. Larger layered-mascon calculations benefit from additional memory and multi-core CPUs. Detailed system requirements are summarized in the Availability section.

2 Availability

2.1 Operating system

GravDyn is platform-independent and has been tested on Linux (Ubuntu 22.04 and later). The software is expected to run on macOS and Windows systems with a compatible Python environment.

2.2 Programming language

Python (version 3.9 or higher)

2.3 System requirements

For general use, no special hardware is required. For computationally intensive applications, in particular for generating mascon models with multiple layered internal structure, it is advantageous to have a system with sufficient memory (suggested $\geq$ 16 GB RAM) and multi-core CPU support. The generation of layered mascon distributions involves processing large polyhedral meshes and assigning mass elements for multiple density layers, which can greatly increase memory usage and computational time. Disk space requirements depend on the size of the polyhedral shape model and the associated data sets. Typical use requires less than 1 GB of space; however, larger shape models and precomputed data sets may require more.

2.4 Dependencies

GravDyn relies on the following core Python libraries:

NumPy (numerical computations)
SymPy (symbolic computation for PSEM)
Matplotlib (visualization)
Trimesh (geometry handling and mesh processing)

Optional dependencies:

JAX (accelerated computation and automatic differentiation)
Pandas (data handling for mascon datasets)

All dependencies can be installed via standard package managers such as pip.

2.5 List of contributors

Safwan Aljbaae
Conceived and led the project, developed the GravDyn software, and wrote the manuscript.Affiliation: Universidad de Atacama, Instituto de Astronomía y Ciencias Planetarias, Copayapu 485, Copiapó, Chile
Valerio Carruba
Contributed to testing and validation of the software, assisted in refining the implementation, and contributed to the revision and improvement of the manuscript.
Affiliation: School of Natural Sciences and Engineering, São Paulo State University (UNESP), Guaratinguetá, SP 12516-410, BrazilLaboratório Interinstitucional de e-Astronomia, Rio de Janeiro, RJ 20765-000, Brazil
Allan K. de Almeida Jr
Contributed to testing and validation of the software, assisted in refining the implementation, and contributed to the revision and improvement of the manuscript.
Affiliation: CFisUC, Departamento de Física, Universidade de Coimbra, Coimbra, 3004-516, Portugal
Gabriel Caritá
Contributed to testing and validation of the software, assisted in refining the implementation, and contributed to the revision and improvement of the manuscript.
Affiliation: São Paulo State University (UNESP), Instituto de Geociências e Ciências Exatas, Rio Claro, SP 13506-900, Brazil
Antonio F. B. A. Prado
Contributed to testing and validation of the software, assisted in refining the implementation, and contributed to the revision and improvement of the manuscript.
Affiliation: Divisão de Pós-Graduação, Instituto Nacional de Pesquisas Espaciais (INPE), Avenida dos Astronautas, 515, São José dos Campos, 12227-310, São Paulo, Brazil
Marcelo L. Mota
Contributed to the theoretical development of the Potential Series Expansion Method (PSEM).
Affiliation: Federal Institute of São Paulo, IFSP, Avenida Thereza Ana Cecon Breda, s/n - Vila São Pedro, Hortolândia, 13183-250, SP, Brazil
Carlos E. Ferreira Lopes
Contributed to testing and validation of the software, assisted in refining the implementation, and contributed to the revision and improvement of the manuscript.
Affiliation: Universidad de Atacama, Instituto de Astronomía y Ciencias Planetarias, Copayapu 485, Copiapó, Chile

2.6 Software location

Archive

Name: Zenodo.
Persistent identifier: https://doi.org/10.5281/zenodo.20210316.
Licence: MIT License.
Publisher: Zenodo.
Version published: v0.1.2.
Date published: 15-07-2026.

Code repository

Name: GravDyn.
Identifier: https://github.com/safwanaljbaae/GravDyn.git
Licence: MIT License.
Date published: 2026-04-20.

Project website

The official GravDyn website is available at https://gravdyn.linea.org.br. The website is hosted by the Laboratório Interinstitucional de e-Astronomia (LIneA).

Package distribution

GravDyn is also distributed via the Python Package Index (PyPI) and can be installed using:

pip install gravdyn

2.7 Language

All documentation is provided in English.

3 Reuse Potential

GravDyn can be reused in studies that require gravitational fields for irregular small bodies. The most direct applications are spacecraft trajectory design, orbital stability analysis, equilibrium-point calculations, and close-proximity motion near asteroids. Because the same input shape can be used to build polyhedral, mascon, and PSEM models, the package is also useful for method comparison under common geometric and physical assumptions.

The layered mascon implementation is particularly useful when internal structure matters. By assigning densities to radial layers, users can test how simplified heterogeneous interiors modify the potential and acceleration field.

The PSEM implementation serves a different role. Once its coefficients have been generated and stored, the potential and its derivatives can be evaluated rapidly outside the Brillouin sphere. This makes the method useful for large grids of field points and repeated evaluations in long integrations.

GravDyn can also serve as a validation environment for new approximations, since the polyhedral model is available as a reference solution within the same workflow. Additional precomputed PSEM models may be added to the repository as they are generated and validated.

Users can report bugs, request support, or suggest improvements through the GitHub issue tracker associated with the GravDyn repository. For scientific questions or collaboration requests, users may also contact the corresponding author.

Notes

[1] https://github.com/esa/polyhedral-gravity-model.

[2] https://github.com/a-amarante/minor-gravity/tree/v2.1.

[3] The Brillouin sphere is defined as the smallest sphere centered at the center of mass that completely encloses the body. Series expansions of the gravitational potential, such as spherical harmonics or PSEM, are guaranteed to converge only outside this sphere; inside it, the series may diverge or lose accuracy (see, e.g., (13)).

Acknowledgements

The authors acknowledge the use of AI-based writing tools to assist with language editing and readability. All scientific content, results, and interpretations remain the responsibility of the authors. The official GravDyn project website (https://gravdyn.linea.org.br) is hosted by the Laboratório Interinstitucional de e-Astronomia (LIneA), whose infrastructure is supported by FINEP, CNPq, FAPERJ, and the INCT do e-Universo program. GAC is grateful to the São Paulo Research Foundation (FAPESP), grant 2025/07596-3. AKAJ acknowledges support from project SPACE, ref. COMPETE2030-FEDER-00860300, funded by COMPETE 2030 and FCT, Portugal.

Author Contributions

Safwan Aljbaae: Conceived and led the project, developed the GravDyn software, and wrote the manuscript.

Valerio Carruba: Contributed to testing and validation of the software, assisted in refining the implementation, and contributed to the revision and improvement of the manuscript.

Allan K. de Almeida Jr: Contributed to testing and validation of the software, assisted in refining the implementation, and contributed to the revision and improvement of the manuscript.

Gabriel Antonio Caritá: Contributed to testing and validation of the software, assisted in refining the implementation, and contributed to the revision and improvement of the manuscript.

Antonio F. B. A. Prado: Contributed to testing and validation of the software, assisted in refining the implementation, and contributed to the revision and improvement of the manuscript.

Marcelo L. Mota: Contributed to the theoretical development of the Potential Series Expansion Method (PSEM).

Carlos E. Ferreira Lopes: Contributed to testing and validation of the software, assisted in refining the implementation, and contributed to the revision and improvement of the manuscript.