`ShapeHarmonics`#

class ShapeHarmonics#

Bases: object

Convert a polyhedral shape model into spherical harmonic coefficients.

A shape model and a harmonic expansion are two descriptions of the same field, each convenient somewhere different: the polyhedron is exact right down to the surface but costs one term per facet, while the expansion is cheap and truncatable but only converges outside the sphere that circumscribes the body. Turning the first into the second is what lets a shape model be flown with SphericalHarmonicsGravity.

The coefficients follow Werner’s exact recursion over the tetrahedra formed by each facet with the origin, so for a constant-density polyhedron they are exact at every degree — the only error is the truncation you choose and the fidelity of the mesh itself.

For a body of constant density and volume \(V\), the fully normalized coefficients of reference radius \(R\) are [1]

\[\bar{C}_{nm} + i\,\bar{S}_{nm} = \frac{1}{(2n+1)\,V} \int_V \left(\frac{r}{R}\right)^{n} \bar{P}_{nm}(\sin\phi)\, e^{i m \lambda}\, dV.\]

The integrand is a homogeneous polynomial of degree \(n\) in \((x, y, z)/R\), built by a recursion on \((n, m)\) (see scarabaeus.rust.PolyhedronHarmonics). The polyhedron is split into the tetrahedra \((\mathbf{0}, \mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2)\) that each facet forms with the origin. Inside one, \(\mathbf{r} = \alpha\mathbf{v}_0 + \beta\mathbf{v}_1 + \gamma\mathbf{v}_2\), so the integral becomes a sum of monomials in \((\alpha, \beta, \gamma)\). Each of these integrates exactly over the unit simplex:

\[\int_V f\, dV = \sum_{\text{facets}} \det[\mathbf{v}_0\ \mathbf{v}_1\ \mathbf{v}_2] \int_{\alpha+\beta+\gamma \le 1} f\, d\alpha\, d\beta\, d\gamma, \qquad \int_{\alpha+\beta+\gamma \le 1} \alpha^i \beta^j \gamma^k\, d\alpha\, d\beta\, d\gamma = \frac{i!\, j!\, k!}{(i+j+k+3)!}.\]

The output is written in the JSON layout Scarabaeus reads through sph_harm_cs_file, so a shape file goes to a flyable gravity model in one call.

Methods

`center_of_mass`(vertices, faces)	Centre of mass of a constant-density mesh, km.
`degree_power`(Cnm, Snm)	RMS power per degree, the usual way to pick a truncation.
`from_shape`(shape, degree[, faces, ...])	Compute the harmonic coefficients of a constant-density shape model.
`to_coefficient_file`(shape, path, degree[, ...])	Write a shape model's harmonic coefficients as a Scarabaeus file.
`volume`(vertices, faces)	Signed volume enclosed by a mesh, km^3.

See also

scarabaeus.PolyhedralGravity: The exact shape-model field.
scarabaeus.SphericalHarmonicsGravity: Consumer of the file written here.

Notes

Coefficients are normalised so that \(C_{00} = 1\), which is what the loader expects, and the mesh is centred on its centre of mass by default so the degree-1 terms vanish as they should.

A harmonic expansion diverges inside the Brillouin sphere — the sphere of radius ref_radius about the origin. For an elongated body that is a large region, and it is exactly where the polyhedral model should be used instead. The reference radius defaults to the circumscribing radius so the boundary is explicit in the file.

References

[1]

Werner, R. A. (1997). Spherical harmonic coefficients for the potential of a constant-density polyhedron. Computers & Geosciences, 23(10), 1071-1077.

Examples

# shape file in, coefficient file out
ShapeHarmonics.to_coefficient_file(
    "bennu.obj", "bennu_deg16.json", degree=16,
    gm=4.891e-9, body_name="Bennu",
)

# then fly it
force_model = ForceModelTranslation(
    primary_body=sc, sph_harm=True, sph_harm_body=bennu,
    sph_harm_order=16, sph_harm_cs_file="bennu_deg16.json",
    sph_harm_norm_flag=True,
)

classmethod center_of_mass(vertices, faces) → ndarray#

Centre of mass of a constant-density mesh, km.

Parameters:

verticesnp.ndarray: (N, 3) coordinates.
facesnp.ndarray: (M, 3) facet indices.

Returns:

np.ndarray: 3-vector centroid.

Raises:

ShapeHarmonicsError: If the mesh encloses no volume.

classmethod from_shape(shape, degree: int, faces=None, ref_radius=None, normalized: bool = True, center: bool = True, scale: float = 1.0)#

Compute the harmonic coefficients of a constant-density shape model.

Parameters:

shapestr or Path or object or np.ndarray: The body’s shape, in any form PolyhedralGravity accepts: a Wavefront .obj or ASCII vertex/facet table path, a mesh object exposing vertices/faces, a (vertices, faces) pair, or a vertex array with faces supplied separately.
degreeint: Maximum degree and order to compute.
facesnp.ndarray, optional: Facet indices, when shape is a bare vertex array.
ref_radiusfloat or ArrayWUnits, optional: Normalising radius in km. Defaults to the circumscribing radius of the mesh, which is the Brillouin sphere the expansion is valid outside of.
normalizedbool, optional: Fully normalised coefficients when True (the default and what Scarabaeus prefers), raw when False.
centerbool, optional: Shift the mesh onto its centre of mass first. Defaults to True; without it the degree-1 terms carry the offset instead of vanishing.
scalefloat, optional: Multiplier taking the shape’s length unit to km.

Returns:

Cnm, Snmnp.ndarray: (degree + 1, degree + 1) coefficient arrays, normalised to C[0, 0] = 1; only entries with n >= m are populated.
ref_radiusfloat: The reference radius the coefficients are expressed against, km.

Raises:

ShapeHarmonicsError: If the degree is negative or the mesh encloses no volume.
RustCoreUnavailableError: If the compiled extension is missing.

classmethod to_coefficient_file(shape, path: str, degree: int, gm=None, density=None, faces=None, ref_radius=None, normalized: bool = True, center: bool = True, scale: float = 1.0, body_name: str = None, source: str = None) → dict#

Write a shape model’s harmonic coefficients as a Scarabaeus file.

The one call that takes a published shape model to something SphericalHarmonicsGravity can load.

Parameters:

shapestr or Path or object or np.ndarray: The body’s shape, in any form from_shape() accepts.
pathstr: Destination .json file.
degreeint: Maximum degree and order to write.
gmfloat or ArrayWUnits, optional: Gravitational parameter in km^3/s^2. Supply this or density: the coefficients are scale-free, but the file records the field strength the model will be flown with.
densityfloat or ArrayWUnits, optional: Bulk density in kg/km^3, from which GM follows via the mesh volume.
facesnp.ndarray, optional: Facet indices, when shape is a bare vertex array.
ref_radiusfloat or ArrayWUnits, optional: Reference radius in km; defaults to the circumscribing radius.
normalizedbool, optional: Write fully normalised coefficients. Defaults to True.
centerbool, optional: Centre the mesh on its centre of mass first. Defaults to True.
scalefloat, optional: Multiplier taking the shape’s length unit to km.
body_namestr, optional: Name recorded in the file. Defaults to the shape file’s stem.
sourcestr, optional: Provenance note recorded in the file.

Returns:

dict: The record written, so the caller can inspect it without re-reading the file.

Raises:

ShapeHarmonicsError: If neither or both of gm and density are supplied, or the shape cannot be reduced to coefficients.

static degree_power(Cnm, Snm) → ndarray#

RMS power per degree, the usual way to pick a truncation.

\[P_n = \sqrt{\sum_{m=0}^{n} \left(\bar C_{nm}^2 + \bar S_{nm}^2\right)}\]

Plotting this against degree shows where the field stops carrying information, which is the degree worth keeping.

Parameters:

Cnm, Snmnp.ndarray: Coefficient arrays.

Returns:

np.ndarray: P_n for each degree.

static volume(vertices, faces) → float#

Signed volume enclosed by a mesh, km^3.

Parameters:

verticesnp.ndarray: (N, 3) coordinates.
facesnp.ndarray: (M, 3) facet indices.

Returns:

float: Enclosed volume; negative when the facets are wound inwards.

exception RustCoreUnavailableError#

Bases: Exception

Raised when the compiled harmonics kernel is not available.

exception ShapeHarmonicsError#

Bases: ValueError

Raised when a shape cannot be reduced to harmonic coefficients.