`PolyhedralGravity`#

class PolyhedralGravity(origin: CelestialBody, base_frame: Frame, shape=None, faces: ndarray = None, density=None, gm=None, body_fixed_frame=None, scale: float = 1.0)#

Bases: DynamicModel

Models the gravitational acceleration of a constant-density polyhedron.

For a small, irregular body the shape is the gravity field: a spherical-harmonic expansion diverges inside the circumscribing sphere, which for an elongated asteroid is most of the space a spacecraft actually occupies. The polyhedral formulation of Werner & Scheeres [1] instead integrates the potential of a homogeneous polyhedron in closed form, one solid-angle term per triangular facet and one line term per edge:

\[U = \frac{1}{2} G \rho \sum_{e \in edges} \mathbf{r}_e \cdot \mathbf{E}_e \mathbf{r}_e L_e - \frac{1}{2} G \rho \sum_{f \in facets} \mathbf{r}_f \cdot \mathbf{F}_f \mathbf{r}_f \omega_f\]

\[\nabla U = -G \rho \sum_e \mathbf{E}_e \mathbf{r}_e L_e + G \rho \sum_f \mathbf{F}_f \mathbf{r}_f \omega_f\]

\[\nabla\nabla U = G \rho \sum_e \mathbf{E}_e L_e - G \rho \sum_f \mathbf{F}_f \omega_f\]

where \(L_e = \ln\frac{a+b+e}{a+b-e}\) is the edge factor, \(\omega_f\) the signed solid angle the facet subtends at the field point, and \(\mathbf{E}_e\) / \(\mathbf{F}_f\) the edge and facet dyads. The result is exact everywhere outside the body — no truncation degree, no divergence — and \(\nabla^2 U = -4\pi G \rho\) inside it, which the model also exposes as a containment test.

The summations run in Rust (PolyhedronGravityCore); Python keeps the shape bookkeeping, the body-fixed to inertial rotation and the estimated-density handling.

Parameters:

originCelestialBody

Body the propagation is centred on, and the body the shape belongs to.

base_frameFrame

Reference frame the propagator integrates in.

shapestr or Path or object or np.ndarray

The body’s shape, in whichever form is at hand:

a path to a Wavefront .obj file, or to an ASCII vertex/facet table (.tab, .txt, .dat) of the kind the PDS distributes for small bodies;
any object exposing vertices and faces (a trimesh.Trimesh and friends drop straight in);
a (vertices, faces) pair;
an (N, 3) vertex array, with faces given separately.

facesnp.ndarray, optional

(M, 3) zero-based vertex indices, one row per triangular facet, wound counter-clockwise as seen from outside the body. Only needed when shape is a bare vertex array.

densityfloat or ArrayWUnits, optional

Bulk density in kg/km^3. Supply this or gm, not both.

gmfloat or ArrayWUnits, optional

Gravitational parameter in km^3/s^2. The density follows from the mesh volume, which is the usual way a shape model is published: a measured GM plus a shape, with the density inferred.

body_fixed_framestr or Frame, optional

Frame the vertices are expressed in. When given, positions are rotated into it through SPICE at each evaluation, so a rotating body is handled properly. Defaults to None, meaning the mesh is already expressed in base_frame and no rotation is applied.

scalefloat, optional

Multiplier taking the file’s length unit to kilometers. Shape models are published in whatever unit the team preferred — pass 1e-3 for a mesh in meters. Defaults to 1.0.

Attributes:

body_fixed_frame: Frame the mesh is expressed in, or None when it is ref_frame.
density: Bulk density in kg/km^3.
faces: (M, 3) zero-based facet vertex indices.
gm: Gravitational parameter in km^3/s^2, from the density and mesh volume.
n_edges: Number of unique edges.
n_facets: Number of triangular facets.
origin_name: Name of the body the shape belongs to.
ref_frame: Reference frame the acceleration is returned in.
vertices: (N, 3) body-fixed vertex coordinates in km.

Methods

`center_of_mass`()	Centre of mass of the shape, in its own frame, km.
`compute_acceleration`(position, ...)	Computes the gravitational acceleration of the polyhedral body.
`compute_potential`(position[, epoch])	Gravitational potential at a field point.
`contains`(position[, epoch])	Whether a point lies inside the body.
`volume`()	Enclosed volume of the shape, in km^3.

Raises:

ShapeModelError: If the mesh is unusable, or if neither/both of density and gm are supplied.
RustCoreUnavailableError: If the compiled extension is missing, or autograd tracing is active.

See also

scarabaeus.SphericalHarmonicsGravity: Harmonic expansion, for bodies whose field is wanted outside the Brillouin sphere.
scarabaeus.PointMassGravity: Single-body Newtonian gravity.

Notes

The mesh must be a closed, consistently wound surface: every edge is shared by exactly two facets. The constructor checks this and reports the offending edge, since a hole or a duplicated facet otherwise produces a plausible but wrong field.

A negative volume() means the facets are wound inwards; flip the last two indices of every facet.

References

[1]

Werner, R. A., & Scheeres, D. J. (1997). Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia. Celestial Mechanics and Dynamical Astronomy, 65(3), 313-344.

Examples

# straight from a shape file
gravity = PolyhedralGravity(eros, J2000, "eros_shape.obj", gm=4.463e-4)

# from anything mesh-shaped already in memory
gravity = PolyhedralGravity(eros, J2000, trimesh.load("eros.obj"),
                            density=2670e9)

# or from explicit arrays
gravity = PolyhedralGravity(eros, J2000, vertices, faces, density=2670e9)

center_of_mass() → ndarray#

Centre of mass of the shape, in its own frame, km.

Worth checking on any shape model before trusting a comparison against a point mass: a mesh whose centroid is off the origin carries a dipole term, and its field will not approach GM/r^2 about the origin at any distance.

Returns:

np.ndarray: 3-vector centroid in km.

compute_acceleration(position: ndarray, current_state: dict, epoch: float, body: CelestialBody) → ndarray#

Computes the gravitational acceleration of the polyhedral body.

Parameters:

positionnp.ndarray: 3x1 spacecraft position relative to the central body, in km, expressed in the reference frame.
current_statedict: Dictionary containing the current state parameters. A ("rho_poly", spice_id) entry overrides the bulk density, so the density can be estimated.
epochfloat: Current simulation epoch in seconds.
bodyCelestialBody: Included for interface consistency; not used.

Returns:

np.ndarray: 3x1 acceleration in km/s^2, in the reference frame.

compute_potential(position: ndarray, epoch: float = None) → float#

Gravitational potential at a field point.

Parameters:

positionnp.ndarray: 3-vector in ref_frame, km.
epochfloat, optional: Epoch used for the body-fixed rotation. Required only when body_fixed_frame is set.

Returns:

float: Potential in km^2/s^2, positive by the convention \(\nabla U = \mathbf{a}\).

contains(position: ndarray, epoch: float = None) → bool#

Whether a point lies inside the body.

Uses the summed solid angle, which is exact for a closed mesh: zero outside, \(4\pi\) inside.

Parameters:

positionnp.ndarray: 3-vector in ref_frame, km.
epochfloat, optional: Epoch used for the body-fixed rotation.

Returns:

bool: True when the point is inside the shape.

volume() → float#

Enclosed volume of the shape, in km^3.

Signed: a negative value means the facets are wound inwards, which is the cheapest check that a mesh follows the outward convention.

Returns:

float: Enclosed volume, km^3.

property body_fixed_frame#: Frame the mesh is expressed in, or None when it is ref_frame.

property density: float#: Bulk density in kg/km^3.

property faces: ndarray#: (M, 3) zero-based facet vertex indices.

property gm: float#: Gravitational parameter in km^3/s^2, from the density and mesh volume.

property n_edges: int#: Number of unique edges.

property n_facets: int#: Number of triangular facets.

property origin_name: str#: Name of the body the shape belongs to.

property ref_frame: Frame#: Reference frame the acceleration is returned in.

property vertices: ndarray#: (N, 3) body-fixed vertex coordinates in km.

exception RustCoreUnavailableError#

Bases: Exception

Raised when the compiled polyhedral kernel is not available.

exception ShapeModelError#

Bases: ValueError

Raised when the supplied shape model cannot define a gravity field.

PolyhedralGravity#

`PolyhedralGravity`#