planetmodel package¶

Subpackages¶

Submodules¶

planetmodel.behaviours module¶

The free functions as methods: an adaptor, and stateless mixins.

A model type is a class derived from Model alone. What its constructor is given is its base information, PREM’s polynomial coefficients say; the model is that information and everything derived from it, fields and behaviour alike, and the mixins here are the common derivations written once: the elastic description completed in both directions (moduli from velocities and velocities from moduli), the elastic tensor and its averages, gravity, and the linear rheologies built from their static parts. A model type gets them by mixing in, never by reimplementing them and never through a hierarchy of model types. Two forms of wrapping are here.

layer_method(fn) turns a function of a layer, such as moduli(layer), into a method model.moduli(which) that resolves which (an index or a name) through model.layer. A function of a model, such as gravity(model, radii) or frozen(model, omega), needs no adaptor: assigned in a class body it is already a method.

The mixins bundle the wrapped methods a kind of model exposes. Each is a class body of such assignments and nothing else, with one exception: Elastic has a constructor hook that attaches the five Love moduli to every layer as first-class fields, since a spherically symmetric elastic medium is those five and there is no reason to hide them behind a call. No mixin holds state of its own, so class PREM(Elastic, ConstantQ, SelfGravitating, Viscoelastic, Model) is still one model type derived from Model alone. The transformations isotropic, elastic and with_gravity return copies of the same class through Model.replaced.

class planetmodel.behaviours.ConstantQ[source]¶

Bases: object

The Love moduli at a frequency under the constant-Q absorption band: kappa dispersed by qkappa and mu by qmu about a reference frequency, following the logarithmic dispersion relation, so that a layer’s five at angular frequency omega are complex fields. The model is not changed; Viscoelastic.frozen is the model that is. The reference frequency is the model’s constant omega_ref where it declares one, else 2 pi rad/s, a period of one second, in the model’s units (reference_omega).

elastic_moduli_at(which: int | str, omega: float, *, reference_omega: float | None = None) → materials.ElasticField[source]¶: The transversely isotropic tensor of a layer at omega.

moduli_at(which: int | str, omega: float, *, reference_omega: float | None = None) → dict[str, Field][source]¶: The five moduli of a layer at omega, complex where the layer holds a Q, else its static five.

reference_omega() → float¶: The reference angular frequency of the constant-Q band in the model’s units: the constant omega_ref where the model declares one (a deck’s header period, say), else 2 pi rad/s, a period of one second, which is PREM’s.

planetmodel.behaviours.ELASTIC_FAMILIES = (('vp', 'vs'), ('vpv', 'vph', 'vsv', 'vsh', 'eta'), ('kappa', 'mu'), ('A', 'C', 'F', 'L', 'N'), ('elastic_moduli',))¶: The ways a layer describes its elastic medium, each a family of names that is complete on its own: the isotropic and the transversely isotropic velocities beside rho, the isotropic moduli, the five Love moduli, and the tensor stored under elastic_moduli.

class planetmodel.behaviours.Elastic(*args: Any, **kwargs: Any)[source]¶

Bases: object

The elastic behaviour of a spherically symmetric model.

On construction the elastic description of every layer is completed in both directions: the five Love moduli A, C, F, L, N are attached where a layer was given rho with velocities or kappa and mu, exact where those are polynomial, and the velocities where a layer was given rho with the five (with_moduli, with_velocities); a layer already holding both, or its tensor under elastic_moduli, is left as it is. The methods then read the fields: moduli(which) the five, elastic_moduli the tensor (the layer’s own of any symmetry where it holds one), kappa_mu the Voigt average, is_fluid the vanishing of shear, and isotropic() the model re-described by its Voigt average, kappa and mu with the five and the velocities recomputed from them. A layer holding a general anisotropic tensor is refused by the three that read the five.

The derived fields are kept in step by with_field: replacing an elastic field on a layer makes the family it belongs to (rho with the isotropic or the transversely isotropic velocities, kappa and mu, or the five) the description, drops the other elastic names from that layer, and completes the description again, so the five follow a new vs and the velocities a new L; a field added beside the existing ones is attached as it is, and a velocity of the other symmetry than the layer’s is refused, as is any change to what a frozen model’s moduli were derived from. A change that goes through replaced directly, as the library’s own transformations do, is the caller’s to keep consistent.

coarsened(*, keep: Iterable[int | str] | None = None, drop: Iterable[int | str] | None = None, degree: int = 3, n: int = 64) → tuple[Model, CoarseningFit][source]¶

Model.coarsened with the elastic description kept consistent.

Each merged group is refit through its base description alone, rho with the velocities where every merged layer holds them (the five then follow as exact polynomial products of the fits), else the five Love moduli, else kappa and mu; the other elastic names are dropped before the fit and derived again from the fitted description, so the coarse model satisfies the elastic relations exactly. A group whose layers share no such description (one speaking only through elastic_moduli, say) is refused, as is a frozen model. Layers left whole are untouched.

elastic_moduli() → Field¶

The elastic tensor of a layer.

The layer’s own elastic_moduli field where it holds one, whatever its symmetry; otherwise an ElasticField built from the fields it holds: ISOTROPIC when the medium is, whether it holds kappa and mu, rho with vp and vs, or five moduli that satisfy the isotropic relations throughout; VTI otherwise; refused with KeyError where it holds no elastic description.

is_fluid() → bool¶: Whether every shear-bearing field the layer holds, among SHEAR_NAMES, vanishes throughout it; a layer holding none of them is refused with KeyError.

isotropic() → Model[source]¶: The model with every elastic description replaced by its Voigt average: each layer keeps rho and everything that is not an elastic name, and holds kappa and mu, and the five recomputed from them, in place of whatever it described its medium by. Exact on polynomial layers. A layer holding no elastic description is left alone.

kappa_mu() → tuple[Field, Field]¶: The Voigt average (kappa, mu) of moduli(layer), as fields.

moduli() → dict[str, Field]¶: The transversely isotropic moduli A, C, F, L, N of a layer as rank-0 fields: the layer’s own where it holds the five, else built from whatever it holds (see elastic_moduli), exact on polynomial layers. A layer holding a general tensor is refused.

with_field(which: int | str, name: str, field: Field, *, replace: bool = False) → Model[source]¶

A copy with field attached to one layer under name.

Replacing an elastic field makes its family the layer’s description and derives the other elastic fields from it again; replacing rho on a layer holding velocities makes the velocities the description, so the moduli follow the new density. A field added beside the existing ones, or a family too incomplete to describe the medium, is attached as it is. A velocity of one symmetry on a layer holding velocities of the other is refused, added or replaced: no single velocity completes the other set, and isotropic() re-describes a layer by vp and vs. On a frozen model (rheology.is_frozen) rho, the elastic names and the rheology names are refused, the frozen moduli having been derived from them.

class planetmodel.behaviours.SelfGravitating[source]¶

Bases: object

The gravity and mass of a spherically symmetric model with a radial density on every layer: the reference body’s, computed on each call and never stored, as numbers by gravity and mass, as one radial field per layer by gravity_fields, and as a copy of the model holding that field under the vocabulary name g by with_gravity; with_stratification attaches d rho / d Phi_0 = rho’(r) / g(r) to the fluid layers. A density that depends on direction is refused, and the geometry’s mapping does not enter (see mesh1d.gravity).

gravity(radii: ArrayLike) → np.ndarray¶

g(r) = G M(r) / r^2 at radii, in the model’s units.

radii broadcast to any shape and must lie in the model, to the geometry’s tolerance; g is zero at the centre and on the inner boundary of a hollow model. M(r) is the antiderivative of each layer’s rho s^2, as gravity_fields has it.

gravity_fields() → tuple[RadialField, ...]¶

g(r) = G M(r) / r^2 as one radial field per layer, named g.

Each field is exact in value and in its first derivative where the layer’s density is polynomial: M(r) is the polynomial antiderivative of 4 pi rho r^2 plus the mass below the layer, and dg/dr = G (4 pi rho - 2 M / r^3). The fields are continuous across every boundary and vanish at the centre; a model attaches them by with_field(i, "g", field).

mass(*, radius: float | None = None) → float¶: The mass inside radius, the outer boundary by default.

with_gravity(*, name: str = 'g', replace: bool = False) → Model[source]¶: The model with its gravity attached to every layer as a radial field under name, exact where the density is polynomial; a layer already holding name is refused unless replace.

with_stratification(*, name: str = 'drho_dphi0', replace: bool = False) → Model[source]¶: The model with the stratification of its fluid layers, d rho / d Phi_0 = rho’(r) / g(r), attached to each of them under name; solid layers, and layers whose fluidity cannot be read, are left alone, and a model with no fluid layer is refused. rho’ is exact where the density is polynomial; the field is zero where g vanishes (the centre). A fluid layer already holding name is refused unless replace.

class planetmodel.behaviours.Viscoelastic[source]¶

Bases: object

The rheology read from the fields: a layer’s viscoelasticity, the model frozen at a frequency, and the elastic model that drops the rheology fields.

elastic() → Model[source]¶: The model without its rheology fields (RHEOLOGY_NAMES), so that every layer is elastic; a frozen model keeps its complex moduli, which are what its rheology became.

frozen(omega: float, *, reference_omega: float | None = None) → Model¶: The model frozen at angular frequency omega: every viscoelastic layer carries its complex moduli under A, C, F, L, N, and the constant omega records the frequency. A model with no viscoelastic layer comes back as a copy with the constant alone; a model that already carries omega is refused, since its moduli are no longer the elastic ones the dispersion starts from. reference_omega defaults to the model’s, see reference_omega.

is_viscoelastic() → bool¶: Whether a layer holds any of RHEOLOGY_NAMES.

planetmodel.behaviours.elastic_family(name: str) → tuple[str, ...][source]¶: The family of ELASTIC_FAMILIES an elastic name belongs to.

planetmodel.behaviours.layer_method(fn: Callable[[...], T]) → Callable[[...], T][source]¶

fn(layer, *args, **kwargs) as a method model.fn(which, *args, **kwargs).

which is a layer index (negatives counting back) or a layer name, resolved by model.layer; the docstring and name are the function’s.

planetmodel.behaviours.with_moduli(fields: Mapping[str, Field]) → dict[str, Field][source]¶: A layer’s fields with the five Love moduli A, C, F, L, N added, read from whatever elastic description it holds; unchanged where it already holds the five, holds its tensor directly, or holds no elastic description at all.

planetmodel.behaviours.with_velocities(fields: Mapping[str, Field]) → dict[str, Field][source]¶: A layer’s fields with the velocities added from rho and the five Love moduli: vp and vs where the five are isotropic, else vpv, vph, vsv, vsh and eta. Square roots are taken pointwise, so these are composed fields, not polynomials. Unchanged where the layer already holds any velocity, lacks rho or the five, or holds its tensor directly.

planetmodel.catalogue module¶

Named model types: PREM from its polynomials, and a simple layered model for use and for testing.

Each is a class derived from Model alone, with the elastic, gravity and rheology behaviours of planetmodel.behaviours mixed in: every layer holds the five Love moduli A, C, F, L, N as exact fields beside the velocities it was built from, PREM().moduli_at("lower_mantle", omega) gives them dispersed by the constant-Q band, PREM().gravity(r) and PREM().frozen(omega) are the free functions of the library as methods. Their constructors build the fields; every copy, surgery, conversion and freezing keeps the class.

PREM is the Preliminary Reference Earth Model of Dziewonski and Anderson (1981, Phys. Earth Planet. Inter. 25, 297-356), built from the low-order polynomials in x = r / a of their Table I with a = 6371 km, so every field is an exact piecewise polynomial and evaluation, derivatives, moduli and integrals carry no interpolation error. The table below holds the coefficients in the paper’s units (g/cm^3, km/s, dimensionless) and the model is in SI. Regions with one P and one S velocity are isotropic and get vph = vpv, vsh = vsv and eta = 1 as exact constants; the fluid outer core and ocean hold no qmu; the elastic values are those at a reference period of 1 s.

LayeredIsotropicElastic builds an isotropic model of exact layers from a few numbers per layer, a constant or the coefficients of a polynomial in r / scale, a fluid layer being one with vs = 0; its homogeneous classmethod is the one-layer case.

MineosModel is the model of a mineos deck, PREM’s own tabulation (examples/data/prem.200) or any other: the table’s columns become the base fields by interpolation (planetmodel.deck), the header names the core layers and sets the reference period of the constant-Q band, and the mixins add the rest.

class planetmodel.catalogue.LayeredIsotropicElastic(boundaries: Sequence[float], *, rho: Sequence[float | Sequence[float]], vp: Sequence[float | Sequence[float]], vs: Sequence[float | Sequence[float]], scale: float = 1.0, layer_names: Sequence[str | None] | None = None, interface_names: Sequence[str | None] | None = None, scales: Scales = Scales.SI, name: str | None = None)[source]¶

Bases: Elastic, SelfGravitating, Model

An isotropic model of exact polynomial layers between boundaries.

boundaries are the skeleton’s, centre outward (an inner radius above zero gives a hollow model); rho, vp and vs give one entry per layer, a number for a constant or a sequence of ascending polynomial coefficients in r / scale, and a layer with vs = 0 throughout is fluid. name is the model’s.

classmethod homogeneous(radius: float, *, rho: float, vp: float, vs: float, name: str | None = None, scales: Scales = Scales.SI) → LayeredIsotropicElastic[source]¶: A uniform isotropic sphere of radius: constant rho, vp and vs; name names the model and its one layer.

class planetmodel.catalogue.MineosModel(source: str | PathLike[str] | Deck, *, kind: str = 'cubic')[source]¶

Bases: Elastic, ConstantQ, SelfGravitating, Viscoelastic, Tabulated, Model

The model of a mineos deck: r rho vpv vsv qkappa qmu vph vsh eta in SI, or the six-column isotropic form, under three header lines.

The columns are interpolated layer by layer by kind (see planetmodel.deck.KINDS); a qmu of zero on a layer marks no shear loss and is kept as read. The header’s nic and noc name the inner and outer core and their boundaries, its tref, in seconds, becomes the constant omega_ref that moduli_at and frozen disperse about, and its title is the model’s name. The deck is kept as deck, its knots per layer as knots and its header as header, so to_deck writes the model back out. The model is in SI, as the deck is; converted gives it in other scales.

class planetmodel.catalogue.PREM(*, ocean: bool = True)[source]¶

Bases: Elastic, ConstantQ, SelfGravitating, Viscoelastic, Model

PREM as an exact polynomial model in SI; see the module docstring.

ocean=False drops the 3 km water layer, so the model ends at the top of the upper crust and that interface is the surface; the polynomials keep their normalisation radius.

planetmodel.catalogue.PREM_RADIUS = 6371000.0¶: The normalisation radius a of PREM’s polynomials, metres.

planetmodel.character module¶

How a field transforms under a mapping, and how elastic symmetry is named.

A field is numbers plus a Character = (rank, weight), and push-forward is one generic operation driven by that character alone. For a mapping m taking the reference body to the physical one, with F[i, j] = d m_i / d X_j and J = det F, a contravariant field of rank n and weight 1 transforms as

T_phys[i1..in](x) = (1/J) F[i1, m1] … F[in, mn] T_ref[m1..mn](X)

with x = m(X); a weight-0 field carries the factors of F but no 1/J, and a rank-0 weight-0 field simply composes with the mapping. In particular rho_phys = rho_ref / J.

field rank weight rho 0 1 second elasticity tensor 4 1 equilibrium stress, stress glut 2 1 Q_kappa, Q_mu, eta, relaxation times, viscosity 0 0

The rheological entries are weight 0 on the understanding that the tensor character of the constitutive law is carried entirely by the elastic tensor: a scalar relaxation time modulating an elastic tensor is invariant, whereas a genuine viscosity tensor would be rank 4 weight 1.

Velocities have no transformation law. v_p and v_s are built from rho and the moduli, neither of which is invariant, and the combination is a tensor of no weight; the moduli are the canonical fields and the velocities are derived from them.

Products of rank-0 fields multiply their characters: weights add, and a product of two weight-1 fields (a density times a modulus, say) has no character here and is refused. Dividing by a weight-1 field is a product with weight -1 and is likewise refused; a modulus over a density is a squared velocity, which is not a field.

Symmetry classifies second elasticity tensors only: those with the full minor and major symmetries, for which the Voigt 6x6 form is faithful. The first elasticity tensor keeps only the major symmetry and has no Voigt form, which is what a Character’s voigt flag records; it is Character(4, 1, voigt=False), a different character from ELASTIC.

class planetmodel.character.Character(rank: int, weight: int, *, voigt: bool = True)[source]¶

Bases: object

The transformation law of a field: (rank, weight), and Voigt or not.

rank is the number of factors of F, weight the power of 1/J, 0 or 1. voigt says whether a rank-2 or rank-4 value carries the symmetries that make the Voigt reduction faithful; it is True for the stress and the second elasticity tensor and False for the first elasticity tensor. It takes part in equality, so a field of each kind has a different character and the algebra keeps them apart.

property component_shape: tuple[int, ...]¶

(3,) * rank.

Type:: Trailing shape of a value array

property is_invariant: bool¶: Whether the field merely composes with the mapping.

rank: int¶

voigt: bool = True¶

property voigt_shape: tuple[int, ...] | None¶

Trailing shape under Voigt reduction, or None if not reducible.

Only the symmetric ranks 2 and 4 reduce, and only where voigt says the symmetries are there; ranks 0 and 1, and a rank-4 value without the minor symmetries, keep their component shape.

weight: int¶

planetmodel.character.DENSITY = Character(rank=0, weight=1, voigt=True)¶: Density, and every modulus, which shares its factor of 1/J.

planetmodel.character.ELASTIC = Character(rank=4, weight=1, voigt=True)¶: The second elasticity tensor.

planetmodel.character.SCALAR = Character(rank=0, weight=0, voigt=True)¶: Q_kappa, Q_mu, eta, relaxation times, viscosity, velocities.

planetmodel.character.STRESS = Character(rank=2, weight=1, voigt=True)¶: Equilibrium stress, stress glut.

class planetmodel.character.Symmetry(*values)[source]¶

Bases: Enum

Material symmetry of a second elasticity tensor, by moduli count.

ISOTROPIC = 2¶

VTI = 5¶

property n_independent: int¶: The number of independent moduli.

planetmodel.deck module¶

Tabulated radial models: decks, their formats, and fields by interpolation.

A deck is a table of knots, radius first, one column per named field, in which a repeated radius marks a boundary: the layering comes from that structure, so a boundary is recorded even where the tabulated values happen to agree across it. What the columns are called and what the header lines mean is a DeckFormat: the column names (by row width, where a format admits more than one), the number of header lines, and how to read and write them. MINEOS is the format of the mineos and PREM decks, three header lines and nine or six columns; GIA that of a viscoelastic model for glacial isostatic adjustment, r rho vpv vsv vph vsh eta viscosity after two header lines (a title and the column names), SI, NaN where a layer lacks a field (the viscosity of an elastic or fluid layer); any other table is a format of its own, given its column names.

read_deck turns a file into a Deck, numbers and a header, nothing else. deck_layers turns a deck into the base fields of a model: the knots of every layer are interpolated by a piecewise polynomial (a cubic spline by default), which is a PolynomialLayer, so the field algebra downstream is exact on it and the moduli of a deck model are products of piecewise polynomials, never refits. A column that is NaN throughout a layer is absent from that layer, the way a fluid layer holds no shear quality factor; NaN on part of a layer is refused. From those base fields the rest of a model is built the normal way, by the mixins of planetmodel.behaviours.

Tabulated is the mixin of a model type built from a deck: it keeps the knots per layer and the header as knots and header, and gives the model back as a deck by to_deck, sampled on those knots, so that a model that was read can be written.

class planetmodel.deck.Deck(radius: ~numpy.ndarray, columns: ~collections.abc.Mapping[str, ~numpy.ndarray], *, header: ~collections.abc.Mapping[str, ~typing.Any] = <factory>, source: str | None = None)[source]¶

Bases: object

A table of knots: radii, one array per named column, and a header.

radius is non-decreasing with every repeat a boundary; columns maps each name to its values at the knots, NaN where a value is absent; header is what the format read from the header lines. Arrays are read-only.

property boundaries: ndarray¶

the ends and every repeated radius.

Type:: The skeleton’s boundary radii

columns: Mapping[str, ndarray]¶

header: Mapping[str, Any]¶

layers() → list[slice][source]¶: The rows of every layer, centre outward: one slice per span between repeated radii, each of at least two strictly increasing knots.

property names: tuple[str, ...]¶: The column names, in the deck’s order.

property nknots: int¶: The number of knots, repeated radii counted twice.

property nlayers: int¶: The number of layers the repeated radii cut the deck into.

radius: ndarray¶

source: str | None = None¶

class planetmodel.deck.DeckFormat(columns: tuple[str, ...] | Mapping[int, tuple[str, ...]], *, name: str = 'custom', header_lines: int = 0, parse_header: Callable[[Sequence[str]], Mapping[str, Any]] | None = None, write_header: Callable[[Deck], Sequence[str]] | None = None)[source]¶

Bases: object

What a family of deck files looks like.

columns names the columns after the radius, either one tuple or a mapping from row width to tuple where the format admits several; header_lines is how many lines precede the table, parse_header turns them into the header mapping a Deck carries and write_header turns a deck back into them. Without a parser the header lines are kept verbatim under “lines”.

columns: tuple[str, ...] | Mapping[int, tuple[str, ...]]¶

header(lines: Sequence[str]) → Mapping[str, Any][source]¶: The header mapping of the given header lines.

header_lines: int = 0¶

name: str = 'custom'¶

names(width: int) → tuple[str, ...][source]¶: The column names of a row with width values after the radius.

parse_header: Callable[[Sequence[str]], Mapping[str, Any]] | None = None¶

write_header: Callable[[Deck], Sequence[str]] | None = None¶

planetmodel.deck.GIA = DeckFormat(columns=('rho', 'vpv', 'vsv', 'vph', 'vsh', 'eta', 'viscosity'), name='gia', header_lines=2, parse_header=<function _gia_header>, write_header=<function _gia_header_lines>)¶

a title line, a line naming the columns, then r rho vpv vsv vph vsh eta viscosity in SI, NaN where a layer lacks a field.

Type:: The GIA deck

planetmodel.deck.GIA_COLUMNS = ('rho', 'vpv', 'vsv', 'vph', 'vsh', 'eta', 'viscosity')¶: The columns of a GIA deck, after the radius.

planetmodel.deck.KINDS = ('cubic', 'natural', 'pchip', 'akima', 'linear')¶: The interpolants a deck’s knots may be joined by, each a piecewise polynomial: a not-a-knot cubic spline, the parabola through three knots; a natural cubic spline, as the old Fortran codes take; a shape-preserving cubic (pchip), Akima’s cubic, or straight lines. Two knots give a line under every kind.

planetmodel.deck.MINEOS = DeckFormat(columns={8: ('rho', 'vpv', 'vsv', 'qkappa', 'qmu', 'vph', 'vsh', 'eta'), 5: ('rho', 'vp', 'vs', 'qkappa', 'qmu')}, name='mineos', header_lines=3, parse_header=<function _mineos_header>, write_header=<function _mineos_header_lines>)¶

three header lines, then r rho vpv vsv qkappa qmu vph vsh eta in SI (nine columns) or r rho vp vs qkappa qmu (six).

Type:: The mineos deck

planetmodel.deck.MINEOS_COLUMNS = ('rho', 'vpv', 'vsv', 'qkappa', 'qmu', 'vph', 'vsh', 'eta')¶: The columns of a transversely isotropic mineos deck, after the radius.

planetmodel.deck.MINEOS_ISOTROPIC_COLUMNS = ('rho', 'vp', 'vs', 'qkappa', 'qmu')¶: The columns of an isotropic mineos deck, after the radius.

class planetmodel.deck.Tabulated[source]¶

Bases: object

The behaviour of a model type built from a deck.

The type’s constructor sets knots, one array of knot radii per layer, and header, the deck’s header; both survive every copy. to_deck samples the model on those knots, every rank-0 radial field by name, NaN where a layer lacks one, so that a model that was read can be written by write_deck, and a model whose fields were changed can be written as the deck it now is. The knots are in the model’s units and follow it through converted; a deck is in SI, so to_deck converts a model in other scales first. A model whose layering no longer matches its knots (after surgery) is refused.

converted(scales: Scales) → Model[source]¶: Model.converted, the knots carried into the new length unit.

header: Mapping[str, Any] = mappingproxy({})¶

knots: tuple[ndarray, ...] = ()¶

to_deck(*, columns: Sequence[str] | None = None) → Deck[source]¶: The model as a deck on its knots, in SI: columns are the names to tabulate, every rank-0 radial real field by default.

planetmodel.deck.deck_knots(deck: Deck) → tuple[ndarray, ...][source]¶: The knot radii of every layer of a deck, read-only copies.

planetmodel.deck.deck_layers(deck: Deck, *, kind: str = 'cubic', specs: Mapping[str, FieldSpec] | None = None) → tuple[Skeleton, list[dict[str, Field]]][source]¶

The skeleton of a deck and the base fields of every layer.

Every column becomes one RadialField per layer on the piecewise polynomial of kind through that layer’s knots, with the character the vocabulary or specs give the name (SCALAR for a name neither knows). A column that is NaN throughout a layer is absent from that layer; NaN on part of one is refused.

planetmodel.deck.mineos_names(deck: Deck) → tuple[list[str | None], list[str | None]][source]¶: Layer and interface names a mineos header implies: nic and noc count the knots of the inner core and to the top of the outer core, so where each sits on a repeated radius the first two layers are inner_core and outer_core under icb and cmb; the outermost interface is surface. Every other name is None. A count of zero says there is no such core; a header whose other counts do not sit on repeated radii is warned about and ignored.

planetmodel.deck.read_deck(source: str | PathLike[str] | Iterable[str], format: DeckFormat | Sequence[str], *, header_lines: int | None = None) → Deck[source]¶

A deck from a file, or from lines: the header lines by the format, then whitespace-separated rows, radius first.

format is a DeckFormat, or the column names after the radius for a table whose header is header_lines verbatim lines (none by default). Blank lines in the table are skipped.

planetmodel.deck.write_deck(path: str | PathLike[str], deck: Deck, format: DeckFormat) → Path[source]¶: Write a deck in a format: its header lines, then one row per knot, radius first and the format’s columns in order. A column the deck lacks is refused, a column the format does not name is not written, and an absent value is written as nan.

planetmodel.displacement module¶

The radial displacement h(r, theta, phi) that drives a radial stretch.

A radial mapping is driven by one scalar function: how far each point moves along its own radius. The callable is the interface,

class RadialDisplacement(Protocol):
def __call__(self, r, theta, phi) -> array

and everything else is an optional enrichment discovered by attribute:

radial_derivative(r, theta, phi) -> dh/dr angular_gradient(r, theta, phi) -> (dh/dtheta, dh/dphi) knots -> radii where dh/dr may jump bounds() -> (min h, max h) scaled(k) -> k h(r / k, theta, phi)

the last the same displacement on a body k times the size, in its own form, so that a change of length scale keeps a stack of reliefs a stack.

A bare function is accepted wherever a displacement is wanted, with the derivatives taken by central differences; supplying them exactly gives exact deformation gradients. as_displacement(h) adapts any callable to an object carrying every method, exact where h supplied it. Two displacements are shipped: flattening, the degree-2 shape of a hydrostatic ellipsoid, and layer_linear, which gives every boundary of a skeleton an analytic relief and interpolates linearly in r between them within each layer.

Knots make smoothness checkable. A displacement that is C0 but whose radial derivative jumps declares the radii where it does, which a geometry takes anywhere within its skeleton and a radial mesh takes as element edges. A displacement declaring no knots asserts that it is smooth.

class planetmodel.displacement.CallableDisplacement(fn: SphericalFunction, *, knots: Iterable[float] = (), radial_derivative: SphericalFunction | None = None, angular_gradient: AngularGradientFunction | None = None, name: str | None = None, dstep: float = 1e-06)[source]¶

Bases: object

Any callable h(r, theta, phi), with its derivatives exact or differenced.

Derivatives passed in, or carried by fn as attributes named radial_derivative and angular_gradient, are used as given; otherwise central differences are taken, with step dstep in the angles (angular_difference) and max(dstep, dstep * |r|) in the radius (radial_difference), one-sided where they would leave the sphere or step below the centre. knots declares the radii where dh/dr jumps, by default those fn carries; leaving them empty asserts smoothness. bounds is fn’s where it has one.

angular_gradient(r: ArrayLike, theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray][source]¶: (dh/dtheta, dh/dphi), exact if supplied, else differenced.

radial_derivative(r: ArrayLike, theta: ArrayLike, phi: ArrayLike) → ndarray[source]¶: dh/dr, exact if supplied, else a central difference.

scaled(k: float) → CallableDisplacement[source]¶: k h(r / k, theta, phi): dh/dr read at r / k and the angular gradient multiplied by k, each exact where this displacement’s is and differenced as it is otherwise; the knots and the bounds multiplied by k.

class planetmodel.displacement.LayerLinear(skeleton: Skeleton, reliefs: Iterable[Relief | None], *, dstep: float = 1e-06)[source]¶

Bases: object

Reliefs on a skeleton’s boundaries, interpolated linearly in r.

reliefs gives one callable relief(theta, phi) per boundary of skeleton, or None for a boundary that stays spherical; within layer i the displacement is the linear interpolation in r between the reliefs of its two boundaries, so h is continuous, its radial derivative jumps only at the boundaries (the knots), and every boundary takes exactly its relief. A relief carrying an angular_gradient(theta, phi) attribute is used for the angular gradient; any other is differenced by angular_difference with step dstep.

angular_gradient(r: ArrayLike, theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray][source]¶: (dh/dtheta, dh/dphi), exact where the reliefs supply theirs.

property boundaries: ndarray¶: The skeleton boundaries the reliefs sit on, innermost first.

radial_derivative(r: ArrayLike, theta: ArrayLike, phi: ArrayLike) → ndarray[source]¶: dh/dr: the slope of the layer’s interpolation, exact.

property reliefs: tuple[Relief | None, ...]¶: One relief per boundary, None where a boundary carries none.

scaled(k: float) → LayerLinear[source]¶: The same reliefs on a body k times the size: the boundaries and every relief multiplied by k, a relief by its own scaled where it has one (a gridded relief stays gridded).

class planetmodel.displacement.RadialDisplacement(*args, **kwargs)[source]¶

Bases: Protocol

h(r, theta, phi): how far a point moves along its radius.

Required: __call__(r, theta, phi) -> array, broadcasting. Optional, used when present: radial_derivative, angular_gradient, knots, bounds, scaled.

class planetmodel.displacement.ZeroDisplacement[source]¶

Bases: object

No displacement: the identity mapping’s h.

angular_gradient(r: ArrayLike, theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray][source]¶: (dh/dtheta, dh/dphi): zeros of the broadcast shape.

bounds() → tuple[float, float][source]¶: (min h, max h): both zero.

knots: tuple[float, ...] = ()¶

radial_derivative(r: ArrayLike, theta: ArrayLike, phi: ArrayLike) → ndarray[source]¶: dh/dr: zeros of the broadcast shape.

scaled(k: float) → ZeroDisplacement[source]¶: No displacement on a body of any size.

planetmodel.displacement.angular_difference(f: Callable[[ndarray, ndarray], ArrayLike], theta: ndarray, phi: ndarray, *, step: float) → tuple[ndarray, ndarray][source]¶: (df/dtheta, df/dphi) of a function on the sphere by central differences of step radians; a colatitude stepped past a pole is read across it, at the opposite longitude, so that every evaluation lies on the sphere and the difference through a pole is that of the function’s smooth continuation.

planetmodel.displacement.as_displacement(fn: RadialDisplacement | SphericalFunction, *, knots: Iterable[float] = (), **kw: object) → RadialDisplacement[source]¶

Adapt any callable of (r, theta, phi) to the protocol.

An object that already carries radial_derivative, angular_gradient and knots is returned unchanged; extra arguments are then refused rather than ignored. Anything else is wrapped in a CallableDisplacement, which keeps what it carries of the three and its bounds.

planetmodel.displacement.flattening(f: float) → CallableDisplacement[source]¶

h = -(2f/3) r P2(cos theta): the degree-2 shape of flattening f.

Every sphere of the reference body becomes a spheroid of the same mean radius with polar radius r (1 - 2f/3) and equatorial radius r (1 + f/3), so that its flattening, the difference of the two over the equatorial radius, is f / (1 + f/3), f to first order. h is proportional to r, so the shape is the same at every radius and on a body of any size.

planetmodel.displacement.layer_linear(skeleton: Skeleton, reliefs: Iterable[Relief | None], *, dstep: float = 1e-06) → LayerLinear[source]¶: The displacement interpolating boundary reliefs linearly in r.

planetmodel.displacement.radial_difference(f: Callable[[ndarray], ArrayLike], r: ndarray, *, step: float) → ndarray[source]¶: df/dr by central differences of max(step, step |r|), one-sided (second order) where the lower point would fall below r = 0.

planetmodel.expansions module¶

Fields and surface functions given by spherical-harmonic expansions.

A HarmonicField is a scalar field on one interval written as

f(r, theta, phi) = sum_lm c_lm(r) Y_lm(theta, phi),

with Y_lm the real orthonormal harmonics of planetmodel.harmonics and the coefficients c_lm(r) one scipy piecewise polynomial of radius whose values have shape (2, lmax + 1, lmax + 1), in that module’s layout. It is exact in r on polynomial pieces and exact in angle to rounding: the form of a model built on spherical harmonics and radial splines, the oracle the gridded fields are held to, and what a spherical random field is. A HarmonicFunction is the same for a function of (theta, phi) alone, a fixed coefficient array, with its two angular derivatives, so that it serves as a relief of layer_linear and as a surface quantity.

Both evaluate through scipy’s harmonics (harmonics.synthesise), a block of points at a time, at a cost of order (lmax + 1)^2 a point and with the coefficients held once for each distinct radius; on the nodes of a grid, on_grid takes ducc0’s grid transforms (a shipped kind of band at least lmax) or its synthesis at points (any other grid), the fast routes of the harmonics extra. Entries of a coefficient array that belong to no harmonic (the sine slot of order zero, orders above the degree) must be zero. Everything is numbers.

class planetmodel.expansions.HarmonicField(interval: tuple[float, float], coefficients: PPoly | BSpline | ArrayLike, *, character: Character = Character(rank=0, weight=0, voigt=True), name: str | None = None, rtol: float = 1e-09)[source]¶

Bases: FieldBase

sum_lm c_lm(r) Y_lm(theta, phi) on one interval, a scalar field.

coefficients is a scipy PPoly or BSpline of radius whose values have shape (2, lmax + 1, lmax + 1), or such an array for coefficients constant in radius; beyond its breakpoints a piecewise polynomial continues its end pieces. character must be of rank 0. Values are real.

property coefficients: PPoly¶: The coefficients as a piecewise polynomial of radius, values of shape (2, lmax + 1, lmax + 1).

coefficients_at(r: ArrayLike) → ndarray[source]¶: c_lm at the radii r: shape (2, lmax + 1, lmax + 1) + r’s shape.

lateral = 'harmonic'¶

property lmax: int¶: The highest degree of the expansion.

on_grid(r: ArrayLike, grid: AngularGrid) → np.ndarray[source]¶: The values at radii r, within the interval, on the nodes of grid, (r.size, ntheta, nphi): through ducc0 where it is installed (its grid transforms on a shipped kind of band at least lmax, its synthesis at points otherwise), through scipy’s harmonics point by point without it.

on_interval(lo: float, hi: float) → HarmonicField[source]¶: The same expansion on [lo, hi], the end pieces continuing.

property ppoly: PPoly¶: The coefficients’ piecewise polynomial, as radial_pieces_of reads a layer function.

radial = 'polynomial'¶

exact pieces in radius and the expansion between nodes.

Type:: The rules a harmonic field states for itself

radial_pieces() → tuple[ndarray, int][source]¶: The coefficients’ breakpoints in the interval, its ends included, and their degree.

renamed(name: str | None) → HarmonicField[source]¶: The same field under name.

rescaled(*, k: float, v: float) → HarmonicField[source]¶: v f(r / k, theta, phi): the breakpoints and interval scaled by k, each piece’s coefficients by v / k^p for its power p.

spherical_mean(*, grid: AngularGrid | None = None) → RadialField[source]¶: c_00(r) / sqrt(4 pi), exact polynomial pieces; grid is not needed.

class planetmodel.expansions.HarmonicFunction(coefficients: ArrayLike)[source]¶

Bases: object

sum_lm c_lm Y_lm(theta, phi): a function of the angles alone.

coefficients has shape (2, lmax + 1, lmax + 1), in the layout of planetmodel.harmonics. Calling it gives the values at points; angular_gradient gives (d/dtheta, d/dphi), so that it is a relief for layer_linear. Values are real.

angular_gradient(theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray][source]¶: (d/dtheta, d/dphi) at (theta, phi), broadcast together.

property coefficients: ndarray¶: The coefficient array, (2, lmax + 1, lmax + 1).

property lmax: int¶: The highest degree of the expansion.

planetmodel.fields module¶

Static fields on one interval, and the small algebra between them.

A field is data on one layer: an interval, a character (the rank and weight that say how it transforms), a name, and evaluate(r, theta, phi, *, frame). It knows no skeleton, no sides and no units: a discontinuity is two layers asked separately, and what the numbers mean is the model’s business.

evaluate(r, theta, phi) broadcasts its coordinates and returns float64 — complex128 for a complex-valued field, such as a model frozen at a frequency holds — of the broadcast shape followed by the character’s Voigt shape for ranks 2 and 4, or its component shape otherwise. Components are given in the local spherical frame (e_r, e_theta, e_phi) at the point unless frame="cartesian" asks for Cartesian ones; a rank-1 field rotates as R v, a Voigt rank-2 as M v and a Voigt rank-4 as M C M^T with M the Bond matrix of R, anything else slot by slot. Calling a field, f(r, theta, phi), is evaluate in the spherical frame, and a radial field of rank 0 may be called with the radius alone, f(r). A radius outside the interval by more than rtol of its width is refused: stepping beyond a layer on purpose is on_interval.

Three kinds are shipped. RadialField holds one layer function per stored component and is exact when they are polynomial. AnalyticField is a formula of (r, theta, phi). ComposedField is a pointwise function of other fields, the one escape from the algebra below, which is: + and - between fields of one character on one interval, * and / between fields of rank 0 (weights add and subtract), ** of a rank-0 weight-0 field, scaling by a number, and + and - with a number, which stands for the constant field of rank 0 and weight 0 on the interval, so that 1 + delta is a field and rho + 1 is refused like any sum across characters. When every operand is radial the result is a RadialField on the combined layer functions, exact when they are polynomial; otherwise it is a ComposedField.

class planetmodel.fields.AnalyticField(interval: tuple[float, float], fn: SphericalFunction, *, character: Character = Character(rank=0, weight=0, voigt=True), name: str | None = None, frame: str = 'spherical', rtol: float = 1e-09)[source]¶

Bases: FieldBase

A formula fn(r, theta, phi) returning components in frame.

The formula receives broadcast coordinate arrays and returns the components at each point, full or Voigt for ranks 2 and 4 (a constant tensor is broadcast); the result is presented Voigt. A formula may return complex values.

property fn: SphericalFunction¶: The formula of (r, theta, phi).

property frame: str¶: The frame the formula returns its components in.

on_interval(lo: float, hi: float) → AnalyticField[source]¶: The formula on [lo, hi]: it extends wherever it is asked.

renamed(name: str | None) → AnalyticField[source]¶: The same field under name.

rescaled(*, k: float, v: float) → AnalyticField[source]¶: v fn(r / k, theta, phi) on the interval scaled by k.

class planetmodel.fields.ComposedField(fn: Callable[[...], ArrayLike], sources: Iterable[Field], *, character: Character, name: str | None = None)[source]¶

Bases: FieldBase

fn(*values) pointwise on fields of one interval.

The sources’ stored components in the spherical frame are passed to fn, which returns the stored components of character; it is never sampled or refitted. Radial when every source is.

property fn: Callable[[...], ArrayLike]¶: The pointwise function applied to the sources’ values.

on_interval(lo: float, hi: float) → ComposedField[source]¶: The same composition with every source restated on (lo, hi).

renamed(name: str | None) → ComposedField[source]¶: The same field under name.

rescaled(*, k: float, v: float) → ComposedField[source]¶: The sources on the scaled coordinate, the result times v.

property sources: tuple[Field, ...]¶: The fields whose values the function is applied to.

planetmodel.fields.FRAMES = ('spherical', 'cartesian')¶: The frames a field can present its components in.

class planetmodel.fields.Field(*args, **kwargs)[source]¶

Bases: Protocol

What every field has: three attributes and evaluate.

character: Character¶

evaluate(r: ArrayLike, theta: ArrayLike, phi: ArrayLike, *, frame: str = 'spherical') → ndarray[source]¶: The components at the broadcast coordinates, in frame: the broadcast shape followed by the character’s stored shape.

interval: tuple[float, float]¶

name: str | None¶

class planetmodel.fields.FieldBase[source]¶

Bases: object

The shared behaviour of the shipped fields.

A subclass sets _interval, _character, _name and _rtol, and implements _values(r, theta, phi) returning the stored components in the spherical frame, plus on_interval, rescaled and renamed. This class checks the coordinates, presents the values in the frame asked for, and supplies the algebra.

property character: Character¶

rank, weight, Voigt or not.

Type:: The field’s law under a mapping

property dtype: dtype¶: float64, or complex128 for a complex-valued field, from a probe at the middle of the interval.

evaluate(r: ArrayLike, theta: ArrayLike | None, phi: ArrayLike | None, *, frame: str = 'spherical') → ndarray[source]¶: The components at (r, theta, phi) in frame, float64 or complex128.

evaluate_at(X: ArrayLike, *, frame: str = 'cartesian') → ndarray[source]¶: The components at Cartesian points X of shape (…, 3), in frame.

property interval: tuple[float, float]¶: The (lo, hi) span of radius the field lives on.

is_radial = False¶

property name: str | None¶: The field’s name, or None for an unnamed one.

property rtol: float¶: The fraction of the interval’s width a radius may overshoot.

property stored_shape: tuple[int, ...]¶

stored_shape of the character.

Type:: The trailing shape of the values

class planetmodel.fields.RadialField(interval: tuple[float, float], function: LayerFunctionLike | Sequence | ndarray, *, character: Character = Character(rank=0, weight=0, voigt=True), name: str | None = None, rtol: float = 1e-09)[source]¶

Bases: FieldBase

A field of the radius alone: one layer function per stored component.

function is a layer function, a PPoly, a number or a callable of r for a rank-0 field; for a higher rank it is a nested sequence of those in the stored shape (Voigt for ranks 2 and 4), giving the components in the spherical frame. functions is that array of layer functions; function the single one of a rank-0 field.

derivative(*, nu: int = 1) → RadialField[source]¶: The nu-th radial derivative of every component; the character is kept.

property function: LayerFunction¶: The one layer function of a rank-0 field.

property functions: ndarray¶: The layer functions, an object array of the stored shape.

integrate(a: float, b: float) → float | complex | ndarray[source]¶: The signed integral from a to b: a float for rank 0, else an array.

is_radial = True¶

on_interval(lo: float, hi: float) → RadialField[source]¶: The same layer functions re-stated on [lo, hi], by their own rule.

property radial_degree: int | None¶: The least polynomial degree in r reproducing every component throughout the interval, or None where a component is not a single polynomial piece (interior breakpoints, or not a polynomial at all).

radial_pieces() → tuple[ndarray, int] | None[source]¶: The breakpoints in the interval, its ends included, and the highest degree of an exact piecewise-polynomial description of every component, or None when a component is not polynomial.

renamed(name: str | None) → RadialField[source]¶: The same field under name.

rescaled(*, k: float, v: float) → RadialField[source]¶: v f(r / k) on the interval scaled by k.

spherical_mean(*, grid: AngularGrid | None = None) → RadialField[source]¶: A radial field of rank 0 is its own spherical mean; grid is not needed.

planetmodel.fields.check_frame(frame: str) → None[source]¶: Refuse a frame outside FRAMES, naming it.

planetmodel.fields.constant_field(interval: tuple[float, float], value: ArrayLike, *, character: Character = Character(rank=0, weight=0, voigt=True), name: str | None = None) → RadialField[source]¶: The constant value (a number, or an array of the stored shape) as an exact radial field.

planetmodel.fields.radial_pieces_of(interval: tuple[float, float], functions: Iterable[object]) → tuple[ndarray, int] | None[source]¶: The breakpoints inside interval, its ends included, and the highest degree of a set of layer functions that are all piecewise polynomials (each with a ppoly), or None when one is not.

planetmodel.fields.rotate_stored(values: ndarray, R: ndarray, character: Character) → ndarray[source]¶: Stored components rotated by R, R of shape (…, 3, 3) broadcasting against the points: the Bond matrix on a Voigt form, a factor of R on every slot otherwise, and nothing for rank 0. With R the spherical frame this takes spherical components to Cartesian ones, and with its transpose back.

planetmodel.fields.stored_shape(character: Character) → tuple[int, ...][source]¶: The trailing shape a field of this character presents: Voigt for ranks 2 and 4 when the character says so, full components otherwise.

planetmodel.fields.to_stored(raw: ndarray, point_shape: tuple[int, ...], character: Character, owner: object) → ndarray[source]¶: A formula’s return value as stored components at every point.

planetmodel.frames module¶

The local spherical frame, and tensors moved between frames.

Fields speak (r, theta, phi) and give their components in the local orthonormal frame (e_r, e_theta, e_phi); mappings speak Cartesian points and give Cartesian components. This module is where the two meet. With theta the colatitude and phi the longitude, both in radians,

e_r = (sin th cos ph, sin th sin ph, cos th) e_theta = (cos th cos ph, cos th sin ph, -sin th) e_phi = (-sin ph, cos ph, 0) R = [e_r, e_theta, e_phi] as columns, so R is orthogonal

and R takes spherical-frame components to Cartesian ones: a vector as R v, a rank-2 tensor as R V R^T, and in general one factor of R on every slot, which is what rotate_slots does. The inverse direction is R^T.

Symmetric rank-2 and rank-4 tensors are also carried in Voigt form, the six index pairs in the order (11, 22, 33, 23, 13, 12) with no engineering-strain factors; bond_matrix is the (6, 6) matrix that moves such components between frames, voigt_to_tensor and tensor_to_voigt go between the two forms.

Geographic coordinates are degrees: latitude 90 - theta, north positive, and longitude phi east in [0, 360). geographic and from_geographic are the one definition of the conversion, the poles landing exactly on theta = 0 and pi and any longitude range reduced to one turn.

Nothing here depends on the rest of the package.

planetmodel.frames.MAX_RANK = 6¶

above it the contraction has 3**(2 rank) intermediate entries and is refused rather than allocated.

Type:: The rank ceiling of rotate_slots

planetmodel.frames.VOIGT_PAIRS = ((0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1))¶: The index pairs behind the six Voigt slots, in Voigt order.

planetmodel.frames.bond_matrix(R: ArrayLike) → ndarray[source]¶

The (6, 6) Bond matrix of an orthogonal R, shape R.shape[:-2] + (6, 6).

A frame change with orthogonal R takes tensor components as c’_ijkl = R_ia R_jb R_kc R_ld c_abcd, and the Voigt matrix holding those same components, with no engineering-strain factors, therefore transforms as V’ = M V M^T with

M[a, b] = R[i, k] R[j, l] + R[i, l] R[j, k] (k != l)
= R[i, k] R[j, k] (k == l)

for a <-> (i, j) and b <-> (k, l) in Voigt order. A rank-2 Voigt vector transforms by the same M applied once, v’ = M v. With R = spherical_frame(theta, phi) this takes spherical Voigt components to Cartesian ones; the inverse is bond_matrix of R^T.

planetmodel.frames.cartesian_points(r: ArrayLike, theta: ArrayLike, phi: ArrayLike) → ndarray[source]¶: The Cartesian points X = r e_r(theta, phi), shape broadcast + (3,).

planetmodel.frames.from_geographic(latitude: ArrayLike, longitude: ArrayLike) → tuple[ndarray, ndarray][source]¶

(latitude, longitude) in degrees -> (theta, phi) in radians.

theta = 90 - latitude in radians, the poles at exactly 0 and pi; phi is the longitude east in radians, from any range (-180 to 180, 0 to 360, or beyond) reduced into [0, 2 pi). Latitude must lie in [-90, 90].

planetmodel.frames.geographic(theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray][source]¶

(theta, phi) in radians -> (latitude, longitude) in degrees.

Latitude is 90 - theta in degrees, north positive; longitude is phi in degrees east, reduced into [0, 360). theta must lie in [0, pi]. With from_geographic this is the one definition of the conversion.

planetmodel.frames.rotate_slots(values: ArrayLike, M: ArrayLike, rank: int) → ndarray[source]¶

One factor of M on every slot: out_{i..} = M_{iA} … T_{A..}.

values has trailing shape (3,) * rank and M is (…, 3, 3), both broadcasting over leading axes. With M = R this takes spherical components to Cartesian ones; with M = R^T it takes them back; with M = F it is the tensor part of a push-forward. Complex values stay complex.

planetmodel.frames.rotation_subscripts(rank: int) → str[source]¶

The einsum string that puts one matrix factor on every slot.

Each factor is written [out, contracted], so the same string serves a rotation R (new components from old), a deformation gradient F (physical from reference) and their inverses without transposing:

rank 1 …iA,…A->…i rank 2 …iA,…jB,…AB->…ij rank 4 …iA,…jB,…kC,…lD,…ABCD->…ijkl

planetmodel.frames.spherical_coordinates(X: ArrayLike) → tuple[ndarray, ndarray, ndarray, ndarray][source]¶

Cartesian points -> (r, theta, phi, R), with R the frame there.

X has shape (…, 3). At the origin the direction is undefined and (theta, phi) = (pi/2, 0) is returned with the frame that implies; a caller for whom the origin is special tests r itself. No length floor is involved, so the function means the same thing in metres and in non-dimensional units.

planetmodel.frames.spherical_frame(theta: ArrayLike, phi: ArrayLike) → ndarray[source]¶

The local orthonormal frame at (theta, phi), basis vectors as columns.

R[…, :, 0] = e_r, R[…, :, 1] = e_theta, R[…, :, 2] = e_phi, in Cartesian components. theta and phi broadcast; the result has the broadcast shape followed by (3, 3).

planetmodel.frames.tensor_to_voigt(t: ArrayLike, *, rank: int = 4) → ndarray[source]¶

Reduce full components to Voigt by reading the six index pairs.

The inverse of voigt_to_tensor on a tensor with the symmetries; on any other it is a projection that keeps one slot of each symmetry class, so the caller is responsible for knowing the reduction is faithful.

planetmodel.frames.voigt_to_tensor(v: ArrayLike, *, rank: int = 4) → ndarray[source]¶

Expand Voigt components to full ones: (…, 6, 6) -> (…, 3, 3, 3, 3) for rank 4, (…, 6) -> (…, 3, 3) for rank 2.

Every symmetry-related slot receives the same number, so the result has the minor and major symmetries by construction. rank is explicit because a (6, 6) array is one Voigt matrix or six Voigt vectors, and only the caller’s character says which.

planetmodel.frames.wrapped(x: ArrayLike, *, start: float = 0.0, period: float = 6.283185307179586) → ndarray[source]¶

Angles reduced into [start, start + period), by whole periods.

np.mod alone can round a small negative offset up to the period itself; that case is brought back to start.

planetmodel.geometry module¶

A skeleton placed in the physical world by one continuous mapping.

A Geometry is a Skeleton, a Mapping from the reference ball to the physical body, and the names of the layers and interfaces. It is the object a mesher takes and the object on which fields are later hung.

The mapping is one continuous map of the whole reference domain, orientation-preserving, whose gradient may jump only at the radii the mapping declares, its knots: skeleton boundaries, or radii inside the layers where a relief sits on a surface that is no boundary of the model or a mapping read from a mesh or a grid has an element edge. Those are the invariants a Geometry checks on construction: the mapping satisfies the protocol, its Jacobian is positive on a lattice covering every layer and the whole sphere, and it is continuous across every interior boundary. The geometry’s knots are those of the mapping’s within the skeleton; a mapping may be defined on a larger body, and a truncated or hollowed geometry keeps its mapping. A piecewise construction is therefore legitimate as long as the pieces agree where they meet.

Interfaces are the boundaries that separate two layers plus the outer boundary. For a full geometry (innermost radius zero) interface k is skeleton boundary k + 1, and between = (k, k + 1) with -1 standing for the outside. For a hollow geometry the inner boundary is interface 0 with between = (-1, 0), so interface k is skeleton boundary k.

Nothing here knows about units: radii are numbers and every tolerance is rtol times the skeleton’s span.

class planetmodel.geometry.Geometry(skeleton: Skeleton, *, mapping: Mapping | None = None, layer_names: Names = None, interface_names: Names = None, rtol: float = 1e-09, check: bool = True)[source]¶

Bases: object

A skeleton, a mapping and names; see the module docstring.

check=False skips the invariant checks, for a mapping already known to satisfy them.

coarsened(*, keep: Iterable[int] | None = None, drop: Iterable[int] | None = None) → tuple[Geometry, CoarseningMap][source]¶

Interior boundaries removed; identity mapping only.

A merged layer’s name is None; kept interfaces keep theirs.

extended(radii: ArrayLike, *, names: Names = None, interface_names: Names = None) → Geometry[source]¶

Layers appended beyond the outer boundary; identity mapping only.

names name the new layers, interface_names their outer boundaries. A geometry whose mapping is not the identity is refused: extend the skeleton first and build the mapping after.

hollowed(radius: float, *, name: str | None = None) → Geometry[source]¶

The geometry cut at radius from below; the mapping is kept.

The result is hollow, with the cut as its inner interface. A cut on an existing boundary keeps that interface’s name unless name is given; the shortened layer keeps its name.

interface(which: int | str) → InterfaceInfo[source]¶: An interface by index (negatives count back) or by name.

property interfaces: tuple[InterfaceInfo, ...]¶: Index, radius, between and name of every interface, centre outward; see the module docstring for the numbering.

property is_hollow: bool¶: Whether the innermost boundary is above zero.

property is_identity: bool¶: Whether the mapping is known to move nothing.

knots() → tuple[float, ...][source]¶: The radii within the skeleton, to rtol of its span, where the mapping declares its gradient may jump, in increasing order.

layer(which: int | str) → LayerInfo[source]¶: A layer by index (negatives count back) or by name.

property layers: tuple[LayerInfo, ...]¶: Index, interval and name of every layer, centre outward.

property mapping: Mapping¶: The mapping taking the reference body to the physical one.

property nlayers: int¶: The number of layers.

refined(radii: ArrayLike, *, names: Names = None) → Geometry[source]¶

Interior boundaries inserted; the mapping is kept.

A split layer loses its name; names name the new interfaces.

renamed(*, layers: Renames = None, interfaces: Renames = None) → Geometry[source]¶

A copy with layer and interface names replaced.

Each argument is a full sequence of names, or a mapping from index or current name to the new name.

property rtol: float¶: The relative tolerance every check of this geometry uses.

scaled(k: float) → Geometry[source]¶: The same geometry with every length multiplied by k.

property skeleton: Skeleton¶: The boundary radii the geometry is built over.

stretched(h: RadialDisplacement | SphericalFunction, *, name: str | None = None, check: bool = True) → Geometry[source]¶

A copy whose mapping is the radial stretch driven by h.

h is any radial displacement; the stretch’s rmax is this geometry’s outer radius.

truncated(radius: float, *, name: str | None = None) → Geometry[source]¶

The geometry cut at radius; the mapping is kept.

A cut on an existing boundary keeps that interface’s name unless name is given; the shortened layer keeps its name.

validity(*, sample: SphericalSample | None = None) → ValidityReport[source]¶: The mapping’s validity report on sample, or on the lattice.

with_mapping(mapping: Mapping, *, check: bool = True) → Geometry[source]¶: A copy with another mapping, checked unless told otherwise.

class planetmodel.geometry.InterfaceInfo(index: int, radius: float, between: tuple[int, int], *, name: str | None = None)[source]¶

Bases: object

One interface of a geometry.

between is (layer below, layer above); -1 stands for the outside of the geometry, so the outer interface has between = (L - 1, -1) and a hollow geometry’s inner interface has between = (-1, 0).

between: tuple[int, int]¶

index: int¶

name: str | None = None¶

radius: float¶

class planetmodel.geometry.LayerInfo(index: int, interval: tuple[float, float], *, name: str | None = None)[source]¶

Bases: object

One layer of a geometry: its index, interval and optional name.

index: int¶

interval: tuple[float, float]¶

name: str | None = None¶

planetmodel.gridded module¶

Stored fields that vary with direction: values on a grid at GLL nodes.

A GriddedField is a field on one interval held as numbers: the interval is cut into radial elements at breaks, each element carries ngll Gauss-Lobatto-Legendre nodes, and at every node the field is given on one AngularGrid, so that values has shape

(nelem, ngll, ntheta, nphi) + stored shape,

the stored components in the spherical frame at each node (Voigt for the ranks that have it), which is a Sample’s layout cut to one layer. Between the nodes it is read radially by one of mesh1d.RADIAL_METHODS within the element holding the radius (the Lagrange polynomial through its GLL nodes, or straight lines between consecutive nodes), and laterally by one of the rules of planetmodel.lateral at every node of that element. A table of depth slices joined by straight lines is ngll = 2 with the breaks at the slices. The end elements continue beyond the interval, as a polynomial layer’s end pieces do, which is what on_interval asks of them.

Components of rank 1 and above are read in the Cartesian frame, where a field is single-valued at a pole, and turned back into the frame asked for at the point. The values at the poles follow the rule of lateral.polar_values and have one source each: a pole that is a node of the grid is its own row, and a pole the grid lacks is given in poles, of shape (nelem, ngll, 2) + stored shape with the components Cartesian, or found from the grid’s expansion when the grid has a band; poles holds both poles once resolved. GriddedField.sampled builds one from any field, evaluating its values at the poles exactly.

A GriddedFunction is the same thing for a function of (theta, phi) alone: values on a grid, its values at the poles, a lateral rule, and the two angular derivatives, so that it serves as a relief of layer_linear and as a surface quantity such as the thickness of a lithosphere. A GriddedDisplacement is a mapping held the same way: its displacement m(X) - X a gridded vector field on every layer of a skeleton, its deformation gradient the derivatives of the rules that read it, and every element boundary declared as a knot, since the radial rule’s derivative may jump there.

Everything is numbers and every stored array is read-only.

class planetmodel.gridded.GriddedDisplacement(skeleton: Skeleton, fields: Sequence[GriddedField])[source]¶

Bases: MappingBase

The mapping m(X) = X + u(X), with u held on grids.

fields holds one GriddedField of character VECTOR per layer of skeleton, in order, each on its layer’s interval: the displacement m(X) - X in the spherical frame at X, as a sample holds it. The deformation gradient is I + grad u from the radial and lateral rules’ derivatives, the e_phi entry taken on the polar axis as its limit along the meridian; every element boundary of every field is a knot. GriddedDisplacement.sampled builds one from any mapping.

deformation_gradient(X: ArrayLike) → ndarray[source]¶: F = I + du/dr e_r + (1/r) du/dtheta e_theta + (1/(r sin theta)) du/dphi e_phi, each a Cartesian column against a frame vector.

displacement(X: ArrayLike) → ndarray[source]¶: u(X) = m(X) - X, Cartesian, of the shape of X.

property fields: tuple[GriddedField, ...]¶: The displacement’s gridded field on each layer.

jacobian(X: ArrayLike) → ndarray[source]¶: J = det F.

classmethod sampled(mapping: Mapping, mesh: RadialMesh, grid: AngularGrid, *, lateral: str = 'cubic') → GriddedDisplacement[source]¶: mapping’s displacement sampled at the nodes of mesh’s elements on grid, layer by layer, and exactly at both poles.

scaled(k: float) → GriddedDisplacement[source]¶: k m(X / k): the skeleton and the fields’ radii multiplied by k, and the displacement with them.

property skeleton: Skeleton¶: The skeleton the displacement is held over.

class planetmodel.gridded.GriddedField(interval: tuple[float, float], breaks: ArrayLike, grid: AngularGrid, values: ArrayLike, *, character: Character = Character(rank=0, weight=0, voigt=True), name: str | None = None, radial: str = 'polynomial', lateral: str = 'linear', poles: ArrayLike | None = None, rtol: float = 1e-09)[source]¶

Bases: FieldBase

A field held on an angular grid at the GLL nodes of radial elements.

breaks are the element boundaries, strictly increasing and spanning interval (within rtol of its width at each end), and values has shape (nelem, ngll, ntheta, nphi) + stored shape, the components in the spherical frame at the nodes of the elements, with ngll at least two. radial is one of mesh1d.RADIAL_METHODS and lateral one of lateral.METHODS. poles, of shape (nelem, ngll, 2) + stored shape with components in the Cartesian frame, gives the values at the poles the grid lacks, by the rule of lateral.polar_values with rtol of the values’ largest; when it is None they are found from the grid’s expansion, which needs a band and the harmonics extra, and refused otherwise. Values are real.

property breaks: ndarray¶: The radii of the element boundaries.

cartesian_derivatives(r: ArrayLike, theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray, ndarray, ndarray][source]¶: The Cartesian components at flat points and their derivatives in r, theta and phi, each (P,) + stored shape: the radial rule’s weights and their derivatives against the lateral rule’s values and gradient at the nodes of the point’s element. Radii beyond the breaks continue the end elements; nothing is checked against the interval.

property grid: AngularGrid¶: The angular grid at every node.

property lateral: str¶: The lateral rule, one of lateral.METHODS.

property ngll: int¶: The GLL nodes per element.

on_interval(lo: float, hi: float) → GriddedField[source]¶: The same elements on [lo, hi]: beyond the breaks the end elements continue.

property poles: ndarray¶: The values at the north and south poles at every node, (nelem, ngll, 2) + stored shape, components in the Cartesian frame.

property radial: str¶: The radial rule, one of mesh1d.RADIAL_METHODS.

radial_pieces() → tuple[ndarray, int][source]¶: The breakpoints and degree that hold the field exactly on its interval: the breaks inside it and ngll - 1 under the polynomial rule, the nodes inside it and one under the linear rule, with the interval’s ends.

property radii: ndarray¶: The radius of every node, (nelem, ngll).

renamed(name: str | None) → GriddedField[source]¶: The same field under name.

rescaled(*, k: float, v: float) → GriddedField[source]¶: v f(r / k, theta, phi): the breaks and interval scaled by k, the values and the polar values by v.

classmethod sampled(field: Field, grid: AngularGrid, *, breaks: ArrayLike | None = None, ngll: int = 5, radial: str = 'polynomial', lateral: str = 'linear', name: str | None = None) → GriddedField[source]¶: field evaluated on grid at the GLL nodes of elements with breaks (the ends of its interval by default) and ngll nodes each, with its values at both poles evaluated exactly: a gridded field that equals field at every node. The name is the field’s unless name is given.

spherical_mean(*, grid: AngularGrid | None = None) → RadialField[source]¶: The spherical mean as exact polynomial pieces: at every node the exact mean of the lateral rule’s interpolant, joined by the radial rule, so the mean of the field as its rules define it; grid is not needed.

property values: ndarray¶: The stored components at the nodes, (nelem, ngll, ntheta, nphi) + stored shape, in the spherical frame.

class planetmodel.gridded.GriddedFunction(grid: AngularGrid, values: ArrayLike, *, lateral: str = 'cubic', poles: ArrayLike | None = None, rtol: float = 1e-09)[source]¶

Bases: object

A function of (theta, phi) held on an angular grid.

values has shape (ntheta, nphi); lateral is one of lateral.METHODS, the bicubic one by default; poles, shape (2,), gives the values at the poles the grid lacks, by the rule of lateral.polar_values with rtol, and when it is None they are found from the grid’s expansion or refused. Calling it gives the values at points; angular_gradient gives (d/dtheta, d/dphi), so that it is a relief for layer_linear.

angular_gradient(theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray][source]¶: (d/dtheta, d/dphi) at (theta, phi), broadcast together.

property grid: AngularGrid¶: The angular grid the values are on.

property lateral: str¶: The lateral rule, one of lateral.METHODS.

mean() → float[source]¶: The mean over the sphere of the function as its rule defines it.

property poles: ndarray¶: The values at the north and south poles.

classmethod sampled(fn: Relief, grid: AngularGrid, *, lateral: str = 'cubic') → GriddedFunction[source]¶: fn(theta, phi) evaluated on grid and, exactly, at both poles.

scaled(k: float) → GriddedFunction[source]¶: k times the function, on the same grid by the same rule: the relief of a body k times the size.

property values: ndarray¶: The values at the nodes, (ntheta, nphi).

planetmodel.harmonics module¶

Real orthonormal spherical harmonics, and synthesis from coefficients.

The basis is real and orthonormal on the unit sphere: with theta the colatitude, phi the longitude and P_l^m the associated Legendre functions without the Condon-Shortley phase,

Y_l0 = N_l0 P_l(cos theta), Y_lm^cos = sqrt 2 N_lm P_l^m(cos theta) cos(m phi), m > 0, Y_lm^sin = sqrt 2 N_lm P_l^m(cos theta) sin(m phi), m > 0, N_lm = sqrt((2l + 1) / 4 pi (l - m)! / (l + m)!),

so that a field u = sum_lm c_lm Y_lm has sum_lm c_lm^2 = int u^2 dS. Coefficient arrays are laid out as coeffs[0, l, m] for the cosine harmonic of order m and coeffs[1, l, m] for the sine one, zero for m = 0 in the sine slot and for m > l. Scipy’s complex harmonics carry the Condon-Shortley phase; it is removed here. The layout and the convention are those of pyshtools under normalization="ortho" and csphase=1, so a coefficient array here is a pyshtools cilm array as it stands. Of its 2 (lmax + 1)^2 entries, (lmax + 1)^2 belong to a harmonic; packing lists them in the one order in which they are laid out as a vector: by degree, then the cosine harmonics of orders 0 … l, then the sine harmonics of orders 1 … l.

Two routes. real_harmonics and synthesise evaluate the basis at arbitrary points through scipy, vectorised over the points: the definition, and the oracle everything else is held to. synthesise_grid and analyse_grid go between coefficients and values on a grid of a shipped kind – Gauss-Legendre, Driscoll-Healy, or equiangular with or without the poles, from planetmodel.sampling – through ducc0’s transforms: fast, exact for a field band-limited to the grid’s band, and needing the harmonics extra. The longitudes may start anywhere and number anything from 2 lmax + 1 up, provided they are equally spaced. synthesise_points is the fast counterpart of synthesise at arbitrary points, through ducc0’s non-uniform synthesis, accurate to about POINT_EPSILON of the field’s largest value rather than to rounding, and synthesise_points_gradient gives the two angular derivatives there the same way. ducc0 holds complex coefficients a_lm for m >= 0 in healpy’s order, with the Condon-Shortley phase; the real ones here are

a_l0 = c[0, l, 0], a_lm = (-1)^m (c[0, l, m] - i c[1, l, m]) / sqrt 2,

and that is the whole of the conversion.

planetmodel.harmonics.POINT_EPSILON = 1e-13¶: The relative accuracy asked of ducc0’s synthesis at arbitrary points.

planetmodel.harmonics.analyse_grid(values: ArrayLike, grid: AngularGrid, *, lmax: int | None = None) → np.ndarray[source]¶

Coefficients from values on a grid of a shipped kind, through ducc0.

values has shape (…, ntheta, nphi), any extra axes first; the result is (2, lmax + 1, lmax + 1) followed by the extra axes, with lmax the grid’s band by default or a smaller degree to stop at. The inverse of synthesise_grid for a field band-limited to the grid; a projection otherwise.

planetmodel.harmonics.packing(lmax: int) → tuple[ndarray, ndarray, ndarray][source]¶: Index arrays (s, l, m) of the (lmax + 1)^2 entries of a coefficient array that belong to a harmonic, in the order of a vector of them: by degree l, then s = 0 with m = 0 … l, then s = 1 with m = 1 … l. So coeffs[s, l, m] is the vector of an array, and out[s, l, m] = vector puts one back.

planetmodel.harmonics.real_harmonics(lmax: int, theta: ArrayLike, phi: ArrayLike) → ndarray[source]¶: The basis up to degree lmax at colatitudes theta and longitudes phi (broadcast together): shape (2, lmax + 1, lmax + 1) + points.

planetmodel.harmonics.real_harmonics_gradient(lmax: int, theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray][source]¶: (d/dtheta, d/dphi) of the basis up to degree lmax, each of the shape real_harmonics gives, through scipy’s derivatives of the complex harmonics.

planetmodel.harmonics.synthesise(coeffs: ArrayLike, theta: ArrayLike, phi: ArrayLike, *, index: ArrayLike | None = None) → ndarray[source]¶: sum_lm c_lm Y_lm(theta, phi) for coefficients of shape (2, lmax + 1, lmax + 1), or that shape followed by the points’ shape when the coefficients vary from point to point, or followed by one axis of n columns with index, of the points’ shape, naming each point’s column (so that coefficients shared by many points are held once). Through scipy, a block of points at a time.

planetmodel.harmonics.synthesise_gradient(coeffs: ArrayLike, theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray][source]¶: (d/dtheta, d/dphi) of sum_lm c_lm Y_lm at points, for coefficients of shape (2, lmax + 1, lmax + 1): the derivatives of synthesise, through scipy, a block of points at a time.

planetmodel.harmonics.synthesise_grid(coeffs: ArrayLike, grid: AngularGrid) → np.ndarray[source]¶

Values on a grid of a shipped kind from coefficients, through ducc0.

coeffs has shape (2, l + 1, l + 1) or that followed by any extra axes, a radial axis say, with l at most the grid’s band; the result has the extra axes first and then (ntheta, nphi), the layout of a sample. Exact for a field band-limited to the grid’s degree.

planetmodel.harmonics.synthesise_points(coeffs: ArrayLike, theta: ArrayLike, phi: ArrayLike) → ndarray[source]¶

sum_lm c_lm Y_lm at arbitrary points, fast, through ducc0.

coeffs has shape (2, l + 1, l + 1) or that followed by any extra axes; theta in [0, pi] and phi, in any range, broadcast together. The result has the extra axes first and then the points’ shape. It agrees with synthesise to about POINT_EPSILON of the field’s largest value, the accuracy asked of ducc0’s non-uniform synthesis, at a cost nearly independent of the number of points.

planetmodel.harmonics.synthesise_points_gradient(coeffs: ArrayLike, theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray][source]¶

(d/dtheta, d/dphi) of sum_lm c_lm Y_lm at arbitrary points, through ducc0, in the layout and to the accuracy of synthesise_points.

ducc0 gives the gradient’s components d/dtheta and (1 / sin theta) d/dphi; the second is multiplied by sin theta here, so that at a pole d/dphi is zero, as it is for any field smooth there.

planetmodel.lateral module¶

Values between the nodes of an angular grid: the lateral rules.

A field that varies with direction is stored at the nodes of an AngularGrid and read between them by one of three rules, the same here and in every reader of the library’s files:

linear bilinear in (theta, phi) on the padded grid cubic the tensor-product cubic spline on the padded grid,

periodic in longitude and clamped at the poles

harmonic analysis on the grid and synthesis at the point, exact
for a field band-limited to the grid’s band (to the accuracy of harmonics.synthesise_points); only on a grid of a shipped kind with a band, and through ducc0 (the harmonics extra)

The padded grid. The local rules are rules on a rectangle of (theta, phi), so the sphere is first laid out as one. The longitudes phi_0 < … < phi_{n-1} gain the wrap column phi_0 + 2 pi, a copy of the first, and the colatitudes gain a row at each pole the grid lacks as a node, holding the field’s value at that pole across every column; a pole that is a node keeps its own row. A point is read at its colatitude, which must lie in [0, pi], and at its longitude reduced into [phi_0, phi_0 + 2 pi). The bilinear rule is the tensor-product spline of degree one on that grid. The bicubic rule is the tensor-product cubic spline on it with these end conditions:

in longitude, periodic, so that the interpolant and its first two derivatives are continuous across the seam at phi_0;
in colatitude, clamped at each pole to a slope s(phi) = a cos phi + b sin phi, with (a, b) fitted by least squares, over the grid’s longitudes, to the slopes that the not-a-knot spline through each meridian’s values has at that pole. A function with a gradient at a pole has exactly such slopes leaving it, the gradient’s two components; with them the interpolant has a gradient there too, to the accuracy of the spline in longitude between the grid’s columns.

Both local rules interpolate the nodes and the polar values, and both are linear in the values, the bicubic one through the fitted slopes as well.

The polar values have one source each. A pole that is a node of the grid is its own row, which must hold one value, constant along the longitudes to rtol of the values’ largest: a row that varies there describes no field, and making it single-valued is for whoever read the data. A pole the grid lacks takes its value from poles (north then south), as a model file stores it, evaluated exactly by the writer; poles always has both entries, the one for a pole that is a node either agreeing with its row or NaN, meaning not given. Without poles, a grid of a shipped kind takes the values from its own expansion, those of the field band-limited to the grid, and any other grid lacking a pole is refused, so that no polar value is ever guessed. polar_values is that rule. A scalar has one value at a pole; a field with components is interpolated component by component in the Cartesian frame, where it too is single-valued at a pole, and the rotation into that frame belongs to the field, not to this module. The harmonic rule’s expansion defines its own values at the poles, so it checks poles like any rule but does not use them.

Accuracy. A local rule wants an oversampled grid. For a field of band 16 with coefficients falling as 1/(l + 1), on Gauss-Legendre grids, the largest error relative to the field’s largest value is about 0.09 for the bicubic rule at the field’s own band, 2e-4 at four times it and 1e-5 at eight, and 0.3, 2e-2 and 5e-3 for the bilinear one: fourth order against second (the cases of tests/test_lateral.py). Analysis on an oversampled grid is still exact to the field’s band, so a grid of about four times the band serves a local reader and a spectral one alike. The bicubic rule is continuously differentiable everywhere, so a relief read by it gives a mapping whose deformation gradient is single-valued on the polar axis and on the seam; its second derivatives may jump at a pole.

LateralInterpolant builds a rule once over any array whose two trailing axes are a grid’s, and gives its values and its two angular derivatives at any points, and its mean over the sphere: exactly, the mean of the interpolant itself, so that a mean is consistent with the rule that reads the values (for the local rules the coefficients against the integrals of the B-splines, sin(theta) dtheta and dphi; for the harmonic rule the degree-zero term); interpolate builds one and evaluates it. Everything is numbers.

class planetmodel.lateral.LateralInterpolant(values: ArrayLike, grid: AngularGrid, *, method: str, poles: ArrayLike | None = None, rtol: float = 1e-09)[source]¶

Bases: object

One lateral rule, built once over an array of values on a grid.

values has shape (…, ntheta, nphi), the grid’s two axes last, and must be finite and real; the leading shape is shape, and calling the interpolant at colatitudes and longitudes that broadcast to a shape P gives an array of shape shape + P. method is one of METHODS. poles and rtol are those of polar_values, whose rule fixes the values at the poles: poles has shape shape + (2,), north then south, and the harmonic rule checks it without using it.

gradient(theta: ArrayLike, phi: ArrayLike) → tuple[ndarray, ndarray][source]¶: (d/dtheta, d/dphi) of the interpolant at the points, each of the shape __call__ gives: the derivatives of the spline, piecewise constant for the bilinear rule, or of the expansion.

property grid: AngularGrid¶: The angular grid the values are on.

mean() → ndarray[source]¶: The mean over the sphere of the interpolant itself, of shape shape: exact for the local rules, whose B-spline coefficients are contracted with the integrals of the B-splines against sin(theta) dtheta and dphi over one turn, and the degree-zero term of the expansion for the harmonic rule. No weights of the grid are used, so any grid has a mean.

property method: str¶: The lateral rule, one of METHODS.

property shape: tuple[int, ...]¶: The leading shape of the values, carried through to every result.

planetmodel.lateral.METHODS = ('linear', 'cubic', 'harmonic')¶: The lateral rules, in the order of the module docstring.

planetmodel.lateral.check_rule(grid: AngularGrid, method: str) → None[source]¶: Refuse, by name, a rule the grid cannot carry: an unknown method; the harmonic rule on a grid that is not of a shipped kind with a band; the cubic rule on fewer than four rows once the poles the grid lacks are added, or fewer than three longitudes.

planetmodel.lateral.interpolate(values: ArrayLike, grid: AngularGrid, theta: ArrayLike, phi: ArrayLike, *, method: str, poles: ArrayLike | None = None, rtol: float = 1e-09) → np.ndarray[source]¶

Values on a grid read at points by a lateral rule.

LateralInterpolant(values, grid, method=method, poles=poles, rtol=rtol) evaluated at (theta, phi): the leading shape of values followed by the broadcast shape of the points. Build the interpolant itself to read the same values more than once.

planetmodel.lateral.polar_values(values: ArrayLike, grid: AngularGrid, *, poles: ArrayLike | None = None, rtol: float = 1e-09) → np.ndarray[source]¶

The values at the north and south poles, by the rule of the module docstring: shape (…, 2) for values of shape (…, ntheta, nphi).

A pole that is a node gives its row, checked to hold one value to rtol of the values’ largest. A pole the grid lacks gives its entry of poles, which must be finite; an entry of poles for a pole that is a node must agree with its row or be NaN. With poles None, the poles the grid lacks come from its expansion, which needs a grid of a shipped kind with a band and the harmonics extra; any other grid lacking a pole is refused.

planetmodel.layerfunction module¶

Functions of one radius on one interval, exact where they can be.

A layer function is what a field stores per component: a callable of the radius on one interval that also differentiates and integrates itself, re-states itself on another interval, and rescales its coordinate and value. Two kinds are shipped. PolynomialLayer wraps a scipy PPoly and does everything on the coefficients: sums, products and powers of polynomials are polynomials, a derivative or an antiderivative is exact, and the rescaling f(r) -> v f(r / k) is one multiply per coefficient (breakpoints k b, coefficients v c / k^p). NumericLayer wraps any callable and supplies the rest honestly: a central difference with a step relative to the interval’s width, and adaptive quadrature.

Arithmetic between layer functions on the same interval follows the operands: polynomial with polynomial stays polynomial, anything with a numeric layer becomes a numeric closure. Beyond its interval a polynomial continues as the same polynomial and a callable is called as given; on_interval is how that continuation is asked for on purpose, while the fields built on layer functions refuse radii outside their own interval.

class planetmodel.layerfunction.LayerFunction(*args, **kwargs)[source]¶

Bases: Protocol

A function of one radius on one interval.

__call__(r) broadcasts over any shape and returns float64, or complex128 for a complex-valued function; derivative(*, nu) and integrate(a, b) (signed) are the calculus; on_interval(lo, hi) re-states the same function on another interval by its own rule; rescaled(*, k, v) is v f(r / k) on the interval scaled by k. as_layer_function adapts a bare callable, a PPoly or a BSpline.

derivative(*, nu: int = 1) → LayerFunction[source]¶: The nu-th derivative as a layer function; nu=0 is the function.

integrate(a: float, b: float) → float | complex[source]¶: The signed integral from a to b.

interval: tuple[float, float]¶

on_interval(lo: float, hi: float) → LayerFunction[source]¶: The same function restated on [lo, hi].

rescaled(*, k: float, v: float) → LayerFunction[source]¶: v f(r / k) on the interval scaled by k.

class planetmodel.layerfunction.NumericLayer(interval: tuple[float, float], fn: Callable[[ndarray], ArrayLike], *, derivative: Callable[[ndarray], ArrayLike] | None = None, dstep: float = 1e-06)[source]¶

Bases: object

Any callable of the radius on one interval, with honest calculus.

The derivative is a central difference with step dstep times the interval’s width, or the callable given as derivative for the first one; the integral is adaptive quadrature. Arithmetic with anything gives another NumericLayer.

derivative(*, nu: int = 1) → NumericLayer[source]¶: The nu-th derivative: the supplied first derivative and then differences of it, else a finite-difference stencil for the nu-th derivative, centred inside the interval and slid inward near its ends, so that no evaluation leaves it. The step is dstep, or eps^(1 / (nu + 2)) where that is larger, times the width: the balance of truncation and rounding for a nu-th difference.

property fn: Callable[[ndarray], ArrayLike]¶: The callable underneath.

integrate(a: float, b: float) → float | complex[source]¶: The signed integral from a to b by adaptive quadrature, to a relative tolerance of 1e-10 of the integral of |f| as a few samples estimate it; the real and imaginary parts of a complex function separately.

property interval: tuple[float, float]¶: The (lo, hi) span of radius the function is stated on.

on_interval(lo: float, hi: float) → NumericLayer[source]¶: The same callable on [lo, hi].

rescaled(*, k: float, v: float) → NumericLayer[source]¶: v f(r / k) as a closure, on the interval scaled by k.

class planetmodel.layerfunction.PolynomialLayer(interval: tuple[float, float], ppoly: PPoly | BSpline)[source]¶

Bases: object

A piecewise polynomial on one interval, exact in every operation.

Wraps a scipy PPoly whose pieces are local polynomials in (r - x[i]) stored highest degree first. The interval is given separately from the breakpoints, the end pieces continuing beyond them where it is wider; as_layer_function takes the breakpoint range where none is meant.

property degree: int¶: The polynomial degree of the pieces.

derivative(*, nu: int = 1) → PolynomialLayer[source]¶: The nu-th derivative, exact on the coefficients.

integrate(a: float, b: float) → float | complex[source]¶: The signed integral from a to b, exact.

property interval: tuple[float, float]¶: The (lo, hi) span of radius the function is stated on.

is_zero() → bool[source]¶: Whether every coefficient is exactly zero.

on_interval(lo: float, hi: float) → PolynomialLayer[source]¶: The same polynomial on [lo, hi]: the end pieces continue.

property ppoly: PPoly¶: The scipy piecewise polynomial underneath.

rescaled(*, k: float, v: float) → PolynomialLayer[source]¶: v f(r / k) on the interval scaled by k, one multiply per coefficient.

planetmodel.layerfunction.as_interval(interval: tuple[float, float]) → tuple[float, float][source]¶: interval as a (lo, hi) pair of floats, refused unless it increases.

planetmodel.layerfunction.as_layer_function(interval: tuple[float, float], fn: LayerFunctionLike) → LayerFunction[source]¶

Adapt fn to the protocol on interval.

A layer function is returned as it is when its interval is the one given, and re-stated on it otherwise; a PPoly or BSpline becomes a PolynomialLayer, a number a constant, any other callable a NumericLayer. Every constructor here takes the interval first.

planetmodel.layerfunction.as_scalar(x: ArrayLike) → float | complex[source]¶: x as a float, or a complex when it is complex.

planetmodel.layerfunction.as_values(x: ArrayLike) → ndarray[source]¶: x as a float64 array, or complex128 when it is complex.

planetmodel.layerfunction.constant_layer(interval: tuple[float, float], value: float | complex) → PolynomialLayer[source]¶: The constant value on interval, as an exact polynomial.

planetmodel.layerfunction.polynomial_fit(interval: tuple[float, float], fn: Callable[[ndarray], ArrayLike], *, degree: int, n: int | None = None) → PolynomialLayer[source]¶

The least-squares polynomial of degree through fn on interval.

fn is sampled at n Chebyshev points (2 degree + 1 by default); the fit is the way a numeric layer is brought into the exact algebra on purpose, and the caller judges its residual.

planetmodel.layerfunction.polynomial_layer(interval: tuple[float, float], coeffs: ArrayLike, *, scale: float = 1.0) → PolynomialLayer[source]¶: sum_k c_k (r / scale)^k on interval, as one exact piece in r.

planetmodel.layerfunction.same_interval(a: tuple[float, float], b: tuple[float, float], *, rtol: float = 1e-09) → bool[source]¶: Whether two intervals coincide to rtol of the wider width.

planetmodel.mapping module¶

Mappings from the reference body to the physical one.

A mapping m takes points of the spherical reference body to points of the physical body. With

F = (grad m)^T, F[i, j] = d m_i / d X_j, J = det F,

the committed interface is three methods on Cartesian points of shape (…, 3):

class Mapping(Protocol):
def __call__(self, X) -> m(X), (…, 3) def deformation_gradient(self, X) -> F, (…, 3, 3) def jacobian(self, X) -> J, (…)

Anything exposing these is a Mapping. An optional tier, discovered by attribute, adds displacement(X) = m(X) - X, is_valid(...), inverse(x), linearise(delta, X=...), knots (radii where F may jump), is_identity, and scaled(k), the mapping k m(X / k) of a body k times the size in the mapping’s own form; MappingBase provides generic versions of the first two, and scaled_mapping(m, k) asks for scaled and falls back on a ScaledMapping. A radial mapping also carries h, its radial displacement, and deformation_gradient_spherical(X).

Mappings are Cartesian-first because meshes and solvers are; the conversion to and from (r, theta, phi) happens once, through planetmodel.frames. Nothing here knows about units: every tolerance is relative to a length the mapping is given.

class planetmodel.mapping.IdentityMapping[source]¶

Bases: MappingBase

The mapping that moves nothing: F = I, J = 1.

deformation_gradient(X: ArrayLike) → ndarray[source]¶: The identity on every point.

inverse(x: ArrayLike) → ndarray[source]¶: A float copy of x.

property is_identity: bool¶

this is the identity.

Type:: True

is_valid(*, X: ArrayLike | None = None, sample: SphericalSample | None = None) → ValidityReport[source]¶: Valid everywhere, whatever the sample.

jacobian(X: ArrayLike) → ndarray[source]¶: Ones of the points’ shape.

knots: tuple[float, ...] = ()¶

scaled(k: float) → IdentityMapping[source]¶: The identity, on a body of any size.

class planetmodel.mapping.Mapping(*args, **kwargs)[source]¶

Bases: Protocol

m: reference body -> physical body, with its gradient and Jacobian.

Three methods on Cartesian points of shape (…, 3): __call__(X), deformation_gradient(X) returning (…, 3, 3) with F[i, j] = d m_i / d X_j, and jacobian(X) returning det F.

deformation_gradient(X: ArrayLike) → ndarray[source]¶: F at X, (…, 3, 3), with F[i, j] = d m_i / d X_j.

jacobian(X: ArrayLike) → ndarray[source]¶: J = det F at X, of the points’ shape.

class planetmodel.mapping.MappingBase[source]¶

Bases: object

Generic implementations of the optional tier, built on m, F and J.

Mix it in for displacement, right_cauchy_green, is_valid, and refusals of inverse and linearise; nothing requires it.

deformation_gradient(X: ArrayLike) → ndarray[source]¶: F at X; the concrete mapping supplies it.

displacement(X: ArrayLike) → ndarray[source]¶: u(X) = m(X) - X, Cartesian in and out.

inverse(x: ArrayLike) → ndarray[source]¶: X such that m(X) = x, where the mapping can invert itself.

property is_identity: bool¶: Whether this mapping is known to move nothing.

is_valid(*, X: ArrayLike | None = None, sample: SphericalSample | None = None) → ValidityReport[source]¶

Whether J > 0 on the given points or on a (r, theta, phi) sample.

The verdict is only as good as the sample; validity_lattice builds one that covers every layer and the whole sphere.

jacobian(X: ArrayLike) → ndarray[source]¶: J = det F at X; the concrete mapping supplies it.

linearise(delta: RadialDisplacement | SphericalFunction, *, X: ArrayLike) → MappingPerturbation[source]¶: dF and dJ at points X for a perturbation of the mapping.

right_cauchy_green(X: ArrayLike) → ndarray[source]¶: C = F^T F.

class planetmodel.mapping.MappingPerturbation(dF: ndarray, dJ: ndarray)[source]¶

Bases: object

dF (…, 3, 3) and dJ (…) induced by a perturbation of the mapping.

dF: ndarray¶

dJ: ndarray¶

planetmodel.mapping.R_FLOOR_FRACTION = 1e-09¶: A point closer to the origin than this fraction of rmax is treated as the origin, where the radial frame is undefined.

class planetmodel.mapping.RadialStretch(h: RadialDisplacement | SphericalFunction, *, rmax: float | Skeleton | Geometry, name: str | None = None)[source]¶

Bases: MappingBase

m(X) = (r + h) e_r: every point moves along its own radius.

In the local frame (e_r, e_theta, e_phi), which is the same at X and at m(X) because direction is preserved,

F_rr = 1 + dh/dr F_r,th = (1/r) dh/dtheta F_r,ph = (1/(r sin theta)) dh/dphi F_th,th = F_ph,ph = 1 + h/r J = (1 + dh/dr) (1 + h/r)^2

and orientation is preserved exactly when both factors are positive. Cartesian components are R F R^T with R the frame matrix. rmax is the outer radius of the domain the mapping is meant for, given as a number or as a skeleton or geometry whose outer boundary it is: it sets the scale of “close to the origin” and brackets the inverse.

deformation_gradient(X: ArrayLike) → ndarray[source]¶: F in Cartesian components, R F_sph R^T.

deformation_gradient_spherical(X: ArrayLike) → ndarray[source]¶: F in the local (e_r, e_theta, e_phi) frame, shape (…, 3, 3).

inverse(x: ArrayLike, *, tol: float = 1e-12) → ndarray[source]¶

X with m(X) = x, by a safeguarded Newton iteration on every ray.

Direction is preserved, so each point is one scalar problem: find s with s + h(s, theta, phi) equal to the physical radius. Where the mapping is valid the residual increases with s and the root is unique. Each point keeps a bracket, [0, max(rmax, 1.5 times its radius)] widened until it holds the root, and takes the Newton step when it stays inside the bracket and the bisection when it does not; a point whose residual is already positive at the centre goes to the centre. tol is relative to rmax.

property is_identity: bool¶: Whether the driving displacement is the zero one.

is_valid(*, X: ArrayLike | None = None, sample: SphericalSample | None = None) → ValidityReport[source]¶

The analytic conditions 1 + dh/dr > 0 and 1 + h/r > 0.

Conditions on the mapping itself rather than on a discrete Jacobian, checked at the given points or sample. A displacement that is not zero at the origin is reported as its own failure, with margin h / rmax.

jacobian(X: ArrayLike) → ndarray[source]¶: J = (1 + dh/dr)(1 + h/r)^2.

property knots: tuple[float, ...]¶: The radii where dh/dr, and so F, may jump.

linearise(delta: RadialDisplacement | SphericalFunction, *, X: ArrayLike) → MappingPerturbation[source]¶

dF and dJ at X for a perturbation delta of the displacement.

dF_rr = d(dh)/dr dF_r,th = (1/r) d(dh)/dtheta dF_r,ph = (1/(r sin theta)) d(dh)/dphi dF_th,th = dF_ph,ph = dh / r dJ = d(dh)/dr (1 + h/r)^2 + 2 (1 + dh/dr)(1 + h/r) dh / r

dF is returned in Cartesian components.

scaled(k: float) → RadialStretch | ScaledMapping[source]¶: k m(X / k): the stretch of the displacement’s own scaled(k), rmax multiplied by k, or a ScaledMapping of this stretch where the displacement has none.

planetmodel.mapping.SIN_FLOOR = 1e-12¶: sin(theta) below this is the polar axis, where the 1/(r sin theta) entry of F is taken as its limit along the meridian.

class planetmodel.mapping.ScaledMapping(mapping: Mapping, k: float)[source]¶

Bases: MappingBase

A mapping conjugated by a change of length scale: k m(X / k).

Points of the scaled body are divided by k, mapped, and multiplied back. F and J pass through unchanged, since both are dimensionless; knots are multiplied by k. is_valid is the underlying mapping’s validity in its own coordinates; inverse and linearise are its own where it has them (a perturbation delta of the scaled body is delta(k r) / k to it), refused otherwise.

deformation_gradient(X: ArrayLike) → ndarray[source]¶: F of the underlying mapping at X / k, dimensionless as it is.

inverse(x: ArrayLike) → ndarray[source]¶: k m^-1(x / k), where the underlying mapping inverts itself.

property is_identity: bool¶: Whether the underlying mapping is the identity.

is_valid(*, X: ArrayLike | None = None, sample: SphericalSample | None = None) → ValidityReport[source]¶: The underlying mapping’s verdict, in its own coordinates.

jacobian(X: ArrayLike) → ndarray[source]¶: J of the underlying mapping at X / k.

property knots: tuple[float, ...]¶: The underlying mapping’s knots, multiplied by k.

linearise(delta: RadialDisplacement | SphericalFunction, *, X: ArrayLike) → MappingPerturbation[source]¶: The underlying mapping’s dF and dJ at X / k, for delta(k r) / k.

scaled(k: float) → Mapping[source]¶: The underlying mapping scaled by this factor times k, through scaled_mapping.

class planetmodel.mapping.ValidityReport(valid: bool, margin: float, *, worst_point: tuple[float, float, float] | None = None, reason: str = '')[source]¶

Bases: object

Whether a mapping preserves orientation, and where it is worst.

margin is dimensionless: the smallest factor of J found. For an invalid mapping worst_point is (r, theta, phi) of the failure and reason names the failing factor.

margin: float¶

reason: str = ''¶

valid: bool¶

worst_point: tuple[float, float, float] | None = None¶

planetmodel.mapping.outer_radius_of(domain: float | Skeleton | Geometry) → float[source]¶: The outer radius a mapping is meant for: a number, or the outer boundary of a skeleton or of a geometry.

planetmodel.mapping.scaled_mapping(m: Mapping, k: float) → Mapping[source]¶: m conjugated by a change of length scale, k m(X / k), for k > 0: the mapping’s own scaled(k) where it has one, which keeps its form (a stack of boundary reliefs stays a stack), else a ScaledMapping.

planetmodel.mapping.validity(m: Mapping, *, X: ArrayLike | None = None, sample: SphericalSample | None = None) → ValidityReport[source]¶: A mapping’s verdict on the points or sample: its own is_valid where it has one, else J > 0 there (MappingBase.is_valid).

planetmodel.mapping.validity_lattice(skeleton: Skeleton, *, n_r: int = 8, n_theta: int = 25, n_phi: int = 16) → tuple[np.ndarray, np.ndarray, np.ndarray][source]¶

A (r, theta, phi) sample covering a skeleton, for is_valid(sample=...).

n_r radii are laid strictly inside every layer of the skeleton, so a thin layer is sampled as densely as a thick one, and the centre is sampled where the skeleton holds it; theta covers [0, pi] including both poles and phi covers [-pi, pi). Returns broadcastable arrays of shapes (n, 1, 1), (1, n_theta, 1) and (1, 1, n_phi), n being nlayers*n_r, one more for the centre.

planetmodel.materials module¶

Elastic moduli: their conversions, the elastic tensor, and what a layer’s fields imply.

The moduli are the canonical elastic description, because velocities have no transformation law (see character.py). The general elastic medium is a second elasticity tensor, character ELASTIC, which a layer may hold directly under the name elastic_moduli; what this module builds and reads is the spherically symmetric case, where the tensor is at most transversely isotropic about e_r. A transversely isotropic medium with its symmetry axis along e_r is described by the five moduli

A = rho vph^2 C = rho vpv^2 L = rho vsv^2 N = rho vsh^2 F = eta (A - 2 L)

which collapse to the isotropic case at vph = vpv, vsh = vsv, eta = 1, where A = C = kappa + 4 mu / 3, L = N = mu and F = kappa - 2 mu / 3. The Voigt average of a transversely isotropic medium is

kappa = (4 A + C + 4 F - 4 N) / 9 mu = (A + C - 2 F + 5 N + 6 L) / 15

and reduces to (kappa, mu) exactly where the medium is isotropic.

The conversions are written on numbers and arrays with the arithmetic operators alone, so the same code runs on rank-0 fields through the field algebra: a product of polynomial layers is a polynomial layer, and the moduli of a polynomial model are exact. velocities_from_moduli takes square roots and is numbers only.

voigt_matrix builds the Voigt form (the order and the absence of engineering factors are stated in frames.py) in the spherical frame, where (1, 2, 3) = (r, theta, phi) and the symmetry axis is index 1:

[ C F F . . . ] [ F A A-2N . . . ] rows and columns in the order [ F A-2N A . . . ] (rr, tt, pp, tp, rp, rt) [ . . . N . . ] [ . . . . L . ] [ . . . . . L ]

ElasticField holds the independent moduli of its symmetry as rank-0 fields on one interval and builds the Voigt matrix when asked; its Cartesian form is the Bond rotation the base class applies to every Voigt rank-4 field.

The functions of a layer read whatever fields it holds, by name, and refuse by name when the fields they need are absent: is_fluid asks the shear-bearing fields whether they vanish, moduli builds the five moduli, elastic_moduli the tensor, and kappa_mu the Voigt average. A layer holding elastic_moduli itself is answered with that field by elastic_moduli, whatever its symmetry; moduli, the five, and the Voigt average through them are the transversely isotropic reading and refuse such a layer by name. A fully anisotropic medium is therefore stored as its tensor.

planetmodel.materials.ELASTIC_NAMES = ('vp', 'vs', 'vpv', 'vph', 'vsv', 'vsh', 'eta', 'kappa', 'mu', 'A', 'C', 'F', 'L', 'N', 'elastic_moduli')¶

what moduli reads from, and what an isotropic re-description replaces.

Type:: Every name that describes the elastic medium beside rho

class planetmodel.materials.ElasticField(symmetry: Symmetry, moduli: Mapping[str, Field], *, name: str | None = None)[source]¶

Bases: FieldBase

A second elasticity tensor, character ELASTIC, stored by its moduli.

moduli maps the independent moduli of symmetry to rank-0 fields on one interval: kappa and mu for ISOTROPIC, A, C, F, L and N for VTI. The Voigt matrix in the spherical frame is built at each evaluation from the moduli there; evaluate(voigt=False) expands it to the full (3, 3, 3, 3) components. The angles are always required: the Cartesian components depend on direction even where the moduli do not, because the frame carrying the symmetry axis does.

as_symmetry(symmetry: Symmetry) → ElasticField[source]¶

The same medium described in symmetry.

ISOTROPIC promotes to VTI through the field algebra, exactly on polynomial layers; a VTI medium is not isotropic and narrowing is refused. The same symmetry returns the field itself.

evaluate(r: ArrayLike, theta: ArrayLike | None, phi: ArrayLike | None, *, frame: str = 'spherical', voigt: bool = True) → ndarray[source]¶: The tensor at (r, theta, phi) in frame: Voigt (…, 6, 6), or the full (…, 3, 3, 3, 3) components with voigt=False.

property is_radial: bool¶: Whether every modulus is a function of the radius alone.

property moduli: dict[str, Field]¶: The independent moduli, in canonical order.

on_interval(lo: float, hi: float) → ElasticField[source]¶: The same moduli re-stated on [lo, hi], each by its own rule.

renamed(name: str | None) → ElasticField[source]¶: The same field under name.

rescaled(*, k: float, v: float) → ElasticField[source]¶: v C(r / k) on the interval scaled by k, through the moduli.

property symmetry: Symmetry¶: The symmetry whose independent moduli are stored.

planetmodel.materials.MODULI_NAMES = {Symmetry.ISOTROPIC: ('kappa', 'mu'), Symmetry.VTI: ('A', 'C', 'F', 'L', 'N')}¶: The independent moduli of each symmetry, in canonical order.

planetmodel.materials.SHEAR_NAMES = ('vs', 'vsv', 'vsh', 'mu', 'L', 'N')¶: The names whose vanishing makes a layer fluid.

planetmodel.materials.VELOCITY_NAMES = ('vp', 'vs', 'vpv', 'vph', 'vsv', 'vsh', 'eta')¶: The velocity names an elastic description may be given in, beside rho: both sets of VELOCITY_SETS.

planetmodel.materials.VELOCITY_SETS = {Symmetry.ISOTROPIC: ('vp', 'vs'), Symmetry.VTI: ('vpv', 'vph', 'vsv', 'vsh', 'eta')}¶: The velocities of each symmetry beside rho, each set complete on its own: the isotropic two and the transversely isotropic five. A layer is described by one set, never by both.

planetmodel.materials.elastic_moduli(layer: LayerLike) → Field[source]¶

The elastic tensor of a layer.

planetmodel.materials.independent_moduli(layer: LayerLike) → tuple[Symmetry, dict[str, Field]][source]¶

The symmetry a layer states and its independent moduli as fields.

The five moduli when held, read as ISOTROPIC where A = C, L = N and F = A - 2L throughout the layer to a relative 1e-10 at nine radii, and as VTI otherwise; else kappa and mu; else rho with the five transversely isotropic velocities; else rho with vp and vs, which are isotropic; a layer holding none of these, or holding its tensor directly under elastic_moduli, is refused with KeyError.

planetmodel.materials.is_fluid(layer: LayerLike) → bool[source]¶: Whether every shear-bearing field the layer holds, among SHEAR_NAMES, vanishes throughout it; a layer holding none of them is refused with KeyError.

planetmodel.materials.kappa_mu(layer: LayerLike) → tuple[Field, Field][source]¶: The Voigt average (kappa, mu) of moduli(layer), as fields.

planetmodel.materials.kappa_mu_from_moduli(A: Operand, C: Operand, F: Operand, L: Operand, N: Operand) → tuple[Field | ndarray, Field | ndarray][source]¶: The Voigt average (kappa, mu) of the moduli, on numbers or fields.

planetmodel.materials.moduli(layer: LayerLike) → dict[str, Field][source]¶: The transversely isotropic moduli A, C, F, L, N of a layer as rank-0 fields: the layer’s own where it holds the five, else built from whatever it holds (see elastic_moduli), exact on polynomial layers. A layer holding a general tensor is refused.

planetmodel.materials.moduli_from_velocities(rho: Operand, vpv: Operand, vsv: Operand, *, vph: Operand | None = None, vsh: Operand | None = None, eta: Operand | None = None) → dict[str, Field | ndarray][source]¶

The moduli A, C, F, L, N from density and velocities.

With vph, vsh and eta omitted the medium is isotropic: vph = vpv, vsh = vsv, eta = 1, giving A = C and F = A - 2 L. Numbers and arrays broadcast; rank-0 fields on one interval combine through the field algebra.

planetmodel.materials.named(field: Field, name: str) → Field[source]¶: The field under name where it can be renamed, else as it is.

planetmodel.materials.velocities_from_moduli(rho: ArrayLike, A: ArrayLike, C: ArrayLike, F: ArrayLike, L: ArrayLike, N: ArrayLike) → dict[str, ndarray][source]¶

The velocities vpv, vph, vsv, vsh and eta from density and moduli.

The inverse of moduli_from_velocities on numbers and arrays. A velocity is zero where rho is not positive or the modulus is negative, and eta is one where A = 2 L leaves it undefined.

planetmodel.materials.voigt_matrix(symmetry: Symmetry, components: Mapping[str, ArrayLike]) → ndarray[source]¶

The Voigt matrix of a second elasticity tensor in the spherical frame, shape broadcast + (6, 6), with the symmetry axis along e_r.

components maps the independent moduli of symmetry to numbers or arrays that broadcast: kappa and mu for ISOTROPIC, A, C, F, L and N for VTI. The dtype follows the inputs, so complex moduli give a complex matrix.

planetmodel.means module¶

Spherical means: a field reduced to a function of radius.

The spherical mean of a scalar field f on its layer is

fbar(r) = (1 / 4 pi) integral over S^2 of f(r, theta, phi) dOmega,

a RadialField of the field’s character. A field that can give it exactly is asked (spherical_mean, a method): a radial field is its own mean; a harmonic field’s is c_00(r) / sqrt(4 pi), exact polynomial pieces; a gridded field gives at each of its nodes the exact mean of its lateral rule’s interpolant, joined by its radial rule, the mean of the field as its rules define it; a perturbed field combines its parts, exactly from the perturbation’s mean when the perturbation is absolute or relative, and from the mean of its exponential, by quadrature, when it is logarithmic. A radial field without a mean of its own is wrapped as one. Any other field needs grid=, and its mean is that grid’s quadrature at each radius, a NumericLayer. Fields of rank 1 and above have no mean here.

spherically_averaged(model) puts the mean of every field that is not radial in its place, through with_field(replace=True) in each layer’s order, so that a model type’s own rule runs: on an Elastic model the mean of a velocity re-derives the Love moduli, rather than each modulus being averaged apart. The mean and the reference of a perturbed field differ in general; unperturbed gives the second.

planetmodel.means.quadrature_mean(field: Field, grid: AngularGrid) → RadialField[source]¶: The mean of a scalar field by grid’s quadrature at every radius it is asked at: a RadialField over a NumericLayer.

planetmodel.means.spherical_mean(field: Field, *, grid: AngularGrid | None = None) → RadialField[source]¶: The spherical mean of a scalar field, by the rule of the module docstring: the field’s own when it has one, else quadrature on grid.

planetmodel.means.spherically_averaged(model: Model, *, grid: AngularGrid | None = None) → Model[source]¶: model with the spherical mean of every field that is not radial in its place, through with_field(replace=True), layer by layer in each layer’s order of names; a field made radial by an earlier replacement is left as it now is.

planetmodel.model module¶

A model: a geometry with a bag of fields on every layer.

A Model holds a Geometry, one mapping of name to field per layer, the Scales its numbers are in, and the specs and constants that say what the names mean. Everything a field is asked goes through its layer: model.layer(i)["rho"], or model.layer("mantle")["rho"]. A discontinuity is two layers asked separately.

What a name means comes from the vocabulary, or from specs= for a name outside it; a field attached under a name with a spec must have the spec’s character, and a name with no spec is accepted, carries no dimensions, and is refused by name on conversion. Nothing is required: the function that needs a field asks the layer for it. Behaviour shared by groups of models is a free function of a layer or a model.

Units live here and nowhere else. A field is numbers; the model’s scales say what one stored unit of each base dimension is in SI, and converted re-expresses the whole model, geometry included, in other scales, by name and exactly for polynomial layers. A constant is read in the model’s units, G included.

Surgery goes through the geometry’s own and carries the fields: a split layer’s fields are re-stated on each part, a cut layer’s on its remainder, and appended shells hold the fields they are given. Every change is a copy, and every change that can break the model’s invariants (its fields, its layering) validates the copy again.

A model type is a class derived from Model alone, plus the stateless behaviour mixins of planetmodel.behaviours that wrap the free functions as methods; there is no hierarchy of model types. Every copy goes through replaced, a shallow copy with the changed parts swapped in, never through __init__, so a model type may give itself whatever constructor its construction needs and every surgery, conversion and freezing returns an instance of the same class carrying the same instance state.

class planetmodel.model.CoarseningFit(map: CoarseningMap, residuals: Mapping[tuple[int, str], float])[source]¶

Bases: object

What a model coarsening did: the map, and each fit’s misfit.

residuals maps (coarse layer index, field name) to the fit’s relative root-mean-square misfit against the fine model over the merged interval, the absolute RMS where the data vanish, measured on a mesh four times finer than the fit’s nodes. A layer left whole keeps its fields and appears in no entry.

describe() → str[source]¶: One line per fit, the worst first.

map: CoarseningMap¶

residuals: Mapping[tuple[int, str], float]¶

property worst: float¶: The largest residual; zero where nothing was fitted.

class planetmodel.model.Layer(info: LayerInfo, fields: Mapping[str, Field])[source]¶

Bases: object

One layer of a model: its place in the geometry and its fields.

index, interval and name are the geometry’s; fields is a read-only mapping of name to field, also reached by layer[name], name in layer, and names.

property fields: Mapping[str, Field]¶: The layer’s fields by name, read-only.

property index: int¶: The layer’s position, 0 at the centre.

property info: LayerInfo¶

index, interval, name.

Type:: The geometry’s record of this layer

property interval: tuple[float, float]¶: The (lo, hi) span of radius the layer occupies.

property name: str | None¶: The layer’s name, or None.

property names: tuple[str, ...]¶: The names of the fields this layer holds.

class planetmodel.model.Model(geometry: Geometry, layers: Iterable[Mapping[str, Field]], *, scales: Scales = Scales.SI, specs: Mapping[str, FieldSpec] | None = None, constants: Mapping[str, Constant] | None = None, name: str | None = None, check: bool = True)[source]¶

Bases: object

A geometry with fields on every layer; see the module docstring.

layers is one mapping of name to field per geometry layer, empty where a layer has no material. specs gives a FieldSpec for names outside the vocabulary (or overriding it) and constants a Constant beyond the shipped ones. name is what the model is called, in its repr and in the files that write it; it survives every copy. check=False skips validation for fields already known to fit.

property G: float¶: The gravitational constant in the model’s units.

coarsened(*, keep: Iterable[int | str] | None = None, drop: Iterable[int | str] | None = None, degree: int = 3, n: int = 64) → tuple[Model, CoarseningFit][source]¶

Interior interfaces removed, each merged layer’s fields refit.

Exactly one of keep and drop selects interior interfaces by name or index; the outermost interface, and the inner boundary of a hollow model, are not selectable. A layer left whole keeps its field objects. A merged layer holds, for every name all its fine layers carry, the least-squares polynomial of degree through n Chebyshev samples of the fine model, and the returned CoarseningFit reports each fit’s relative misfit for the caller to judge. Merging a fluid layer with a solid one, a name only some of the merged layers hold, and a field that is not a scalar of the radius alone are refused by name. The identity mapping is required, as for every surgery; the fits are independent, so derived fields (an attached gravity, a modulus beside its velocities) are more consistently derived again from the coarse model.

common_names() → tuple[str, ...][source]¶: The names every layer holds.

constant(name: str) → float[source]¶: A declared constant in the model’s units.

property constants: Mapping[str, Constant]¶: The declared constants by name; constant reads their values.

converted(scales: Scales) → Model[source]¶: The model in other scales: every length and every field by name.

extended(radii: ArrayLike, *, fields: str | Sequence[Mapping[str, Field]] | None = None, names: Names = None, interface_names: Names = None) → Model[source]¶

Shells appended outside, holding the fields given.

fields is None for empty shells, “extrapolate” to re-state the outermost layer’s fields on each shell by their own rule, or one mapping of name to field per shell.

field_names() → tuple[str, ...][source]¶: Every name any layer holds, in first-appearance order.

property geometry: Geometry¶: The geometry the fields hang on.

hollowed(radius: float, *, name: str | None = None) → Model[source]¶: The model cut from below at radius.

in_si() → Model[source]¶: The model converted to SI scales.

layer(which: int | str) → Layer[source]¶: A layer by index (negatives allowed) or by name.

property layers: tuple[Layer, ...]¶: The layers, centre outward, each a bag of fields.

layers_with(name: str) → tuple[int, ...][source]¶: The indices of the layers holding name.

property name: str | None¶: What the model is called, or None.

named(name: str | None) → Model[source]¶: A copy called name; None unnames it.

property nlayers: int¶: The number of layers.

nondimensionalised(*, density: float = 5515.0, length: float | None = None) → Model[source]¶

The model in geophysical scales, with G equal to one.

length is the outer radius by default; the model must be in SI.

refined(radii: ArrayLike, *, names: Names = None) → Model[source]¶: Interior boundaries inserted; a split layer’s fields on each part.

renamed(*, layers: Renames = None, interfaces: Renames = None) → Model[source]¶: The model with layers or interfaces renamed; see Geometry.renamed.

A shallow copy with the given parts replaced, validated unless check=False: the one path every copy of a model takes.

layers is one mapping of name to field per layer of the (new) geometry; specs and constants replace the model’s own tables, the vocabulary and the shipped constants being merged in again. The class and every other instance attribute are kept, the name included unless name gives another.

property scales: Scales¶: What one stored unit of each base dimension is in SI.

property skeleton: Skeleton¶: The geometry’s boundary radii.

spec(name: str) → FieldSpec | None[source]¶: The spec of a name, or None when it has none.

property specs: Mapping[str, FieldSpec]¶

the vocabulary, then specs=.

Type:: Every name with a meaning

stretched(h: RadialDisplacement | SphericalFunction, *, name: str | None = None, check: bool = True) → Model[source]¶: The model under the radial stretch driven by h; see Geometry.stretched.

truncated(radius: float, *, name: str | None = None) → Model[source]¶: The model cut at radius; the cut layer’s fields on its remainder.

validate() → None[source]¶: Every field a Field on its layer’s interval, with its spec’s character.

with_field(which: int | str, name: str, field: Field, *, replace: bool = False) → Model[source]¶: A copy with field attached to one layer under name.

with_geometry(geometry: Geometry) → Model[source]¶: A copy on another geometry over the same skeleton.

with_mapping(mapping: Mapping, *, check: bool = True) → Model[source]¶: The model placed by mapping, over the same skeleton.

without_field(name: str, *, layers: Iterable[int | str] | None = None) → Model[source]¶: A copy without name on every layer, or on the layers given.

planetmodel.modelfile module¶

The model file: a model in one netCDF file, on its reference skeleton.

One netCDF-4 file holds a model as numbers: its skeleton, its fields at the nodes of a radial layout and an angular grid, and its mapping, in the model’s own units with its scales recorded, CF-styled so that xarray, Panoply and ncview open it as it stands. The Python here is the definition of the format and the reference the C++ readers are held to; docs/source/formats/model_file.rst describes it for consumers.

Global attributes: schema (SCHEMA, checked by every reader), Conventions, the model’s name, scale_length, scale_mass, scale_time and scale_temperature (what one stored unit is in SI), rtol (the geometry’s tolerance, the default where a file has none), constant_<name> for every constant in the model’s units, and meta_<key> for the writer’s meta. The skeleton, at the root: dimensions layer, interface and bound (two), layer_bounds(layer, bound), layer_name(layer), interface_radius(interface), interface_name(interface) and layer_part(layer).

Parts. The layers are cut into parts, contiguous runs of whole layers, one group each (part_0, part_1, … from the centre), each on an angular grid of its own. A part holds the Sample layout of its layers: dimension node, the per-element GLL nodes flattened element by element, so that an interface is two nodes at one radius; radius(node) (an auxiliary coordinate, since it repeats), element_layer(element), element_bounds(element, bound) and the attribute ngll; the grid as the CF coordinate variables latitude (degrees north, decreasing) and longitude (degrees east, in [0, 360)), with the attributes grid_kind and grid_lmax and the variable colatitude_weights(latitude) where it has them. Then one variable per field, name(node) when every layer holding it is radial and name(node, latitude, longitude) otherwise, the stored components of rank 1 and above on trailing axes (component, voigt, or row and column pairs of them) in the spherical frame at the node, NaN on the nodes of the layers that lack it; attributes dimensions (the exponents of units.BASE), units (for CF tools, never read back), long_name, rank, weight, voigt, layers, radial_interpolation and, laterally, lateral_interpolation. A field that is a radial reference with a perturbation on every layer of the part that holds it perturbed is two variables: name(node), the reference, with attributes perturbation and perturbation_kind, and name_perturbation. Where the grid lacks either pole, every variable that varies with direction has name_poles(node, pole, ...), north then south, its components Cartesian, evaluated exactly by the writer.

The mapping is held in one of two forms. A mapping that offers the reliefs of its surfaces (a radial stretch by layer_linear) is the group geometry: boundary_radius(boundary), the radii of its surfaces, and boundary_displacement(boundary, latitude, longitude) on a grid of its own, with boundary_displacement_poles where that grid lacks a pole and the attributes lateral_interpolation and between (the rule between surfaces, linear). Any other mapping is the field displacement(node, latitude, longitude, component) of every part, m(X) - X in the spherical frame at X, as a sample holds it, read back as a GriddedDisplacement. The identity is neither.

The radial layout is exact where it can be: the writer builds, per part, a radial mesh whose element edges are the layer boundaries, every breakpoint a field declares (radial_pieces) and every knot of the mapping, with ngll one more than the highest declared degree; a field that declares no polynomial description subdivides the elements of its layer to the spacing a sample takes from the grid’s band (every layer, for a mapping written as a displacement field), and raises ngll to at least 5. So a deck model written and read back is the same polynomials, and a band-limited field on a grid of its band is held exactly. radial= gives the layout instead. Every variable states its intended interpolation: a gridded field its own rules, a harmonic field polynomial and harmonic, any other field the defaults (polynomial, and linear laterally) unless interpolation= names others; a reader uses the stated rules unless its caller overrides them.

read_model_file gives the numbers and the metadata as a ModelFile, the way read_deck gives a deck; file_layers and file_geometry turn them into fields and a geometry; read_model is the plain model. netCDF4 (the netcdf extra) is imported where a file is opened. Real fields only; a model is never converted on writing, and is read back in the file’s scales unless scales= asks for others.

planetmodel.modelfile.CONVENTIONS = 'CF-1.8'¶: The CF conventions the file follows; groups are admitted from 1.8.

planetmodel.modelfile.DEFAULT_INTERPOLATION = ('polynomial', 'linear')¶

(radial, lateral).

Type:: The rules a variable states when its field states none

class planetmodel.modelfile.FileGeometry(radii: np.ndarray, grid: AngularGrid, values: np.ndarray, poles: np.ndarray | None, lateral: str, between: str)[source]¶

Bases: object

The mapping as boundary displacements: the surfaces’ radii, their grid, the displacements (boundary, ntheta, nphi) and their polar values (boundary, 2) where the grid lacks a pole, and the rules.

between: str¶

grid: AngularGrid¶

lateral: str¶

poles: np.ndarray | None¶

radii: np.ndarray¶

values: np.ndarray¶

class planetmodel.modelfile.FilePart(layers: tuple[int, ...], mesh: RadialMesh, grid: AngularGrid, variables: Mapping[str, FileVariable])[source]¶

Bases: object

One part: its run of layers, their radial layout as a mesh over the skeleton, the angular grid, and the variables, by name.

grid: AngularGrid¶

layers: tuple[int, ...]¶

mesh: RadialMesh¶

variables: Mapping[str, FileVariable]¶

class planetmodel.modelfile.FileVariable(values: ndarray, attributes: Mapping[str, Any])[source]¶

Bases: object

One variable of a part: its values, as stored, and its attributes.

attributes: Mapping[str, Any]¶

values: ndarray¶

class planetmodel.modelfile.ModelFile(schema: str, name: str | None, scales: Scales, constants: Mapping[str, float], meta: Mapping[str, Any], skeleton: Skeleton, layer_names: tuple[str | None, ...], interface_names: tuple[str | None, ...], parts: tuple[FilePart, ...], geometry: FileGeometry | None, *, declared: Mapping[str, Constant] = mappingproxy({}), rtol: float = 1e-09)[source]¶

Bases: object

A model file’s numbers and metadata, as read_model_file gives them: nothing is interpolated or converted.

constants: Mapping[str, float]¶

declared: Mapping[str, Constant] = mappingproxy({})¶: The constants whose dimensions the file states, as a model holds them.

geometry: FileGeometry | None¶

interface_names: tuple[str | None, ...]¶

layer_names: tuple[str | None, ...]¶

meta: Mapping[str, Any]¶

name: str | None¶

part_of(layer: int) → FilePart[source]¶: The part holding layer.

parts: tuple[FilePart, ...]¶

rtol: float = 1e-09¶: The tolerance of the geometry the file was written from, relative to the skeleton’s span; the default where the file states none.

scales: Scales¶

schema: str¶

skeleton: Skeleton¶

planetmodel.modelfile.SCHEMA = 'planetmodel.model/1'¶: The schema string every reader checks; bumped only for an incompatible change.

planetmodel.modelfile.file_geometry(data: ModelFile, *, interpolation: Interpolation | None = None) → Geometry[source]¶: The geometry a file holds: its skeleton and names, and its mapping from the geometry group or the displacement field, read by the stated rules unless interpolation names others.

planetmodel.modelfile.file_layers(data: ModelFile, *, interpolation: Interpolation | None = None) → list[dict[str, Field]][source]¶: One mapping of name to field per layer: radial variables as exact polynomial layers through the nodes, the others as gridded fields, and references with perturbations as perturbed fields, read by the stated rules unless interpolation names others.

planetmodel.modelfile.read_model(path: str | os.PathLike[str], *, interpolation: Interpolation | None = None, scales: Scales | None = None) → Model[source]¶: The model a file holds, as a plain Model with the file’s specs, in the file’s scales or, with scales=, converted to those.

planetmodel.modelfile.read_model_file(path: str | PathLike[str]) → ModelFile[source]¶: A model file’s numbers and metadata, its structure validated first.

planetmodel.modelfile.validate_structure(source: str | PathLike[str] | Any) → None[source]¶: Check that a file, by path or as an open netCDF4 Dataset, has the structure the module docstring defines, raising ValueError naming the first thing wrong.

planetmodel.modelfile.write_model(path: str | os.PathLike[str], model: Model, *, grid: Grids, radial: RadialMesh | None = None, fields: Iterable[str] | None = None, interpolation: Interpolation | None = None, meta: Mapping[str, Any] | None = None) → Path[source]¶

Write model to a model file at path, which is returned.

grid is one AngularGrid for every layer, or a mapping from layers (by index or name) to grids, every layer given one: contiguous layers on one grid form one part, and a grid’s layers that are not contiguous are refused. radial is a RadialMesh over the model’s skeleton to lay the nodes out by, in place of the exact layout of the module docstring. fields are the names to write, every name the model holds by default; interpolation maps names (and displacement, boundary_displacement) to the (radial, lateral) rules to state in place of the fields’ own; meta becomes meta_<key> attributes. The model is written in its own units.

planetmodel.perturbation module¶

A radial reference with a perturbation: a kind of field, not two names.

A laterally varying model is often stated as a one-dimensional reference and a perturbation of it, a tomographic model as relative perturbations of a reference’s velocities, a viscosity as the logarithm of its departure from a radial profile. PerturbedField(reference, perturbation, *, kind) is such a field, and the kind is how the two combine:

absolute reference + perturbation relative reference (1 + perturbation) log reference exp(perturbation), the natural logarithm

The reference is a radial field of rank 0; the perturbation is a field of rank 0 on the same interval, of the reference’s character when it is absolute and dimensionless (SCALAR) otherwise. The field has the reference’s character and is radial exactly when the perturbation is. reference, perturbation and kind are attributes, so that a writer of files and a spherical mean ask the field for them rather than a central module inspecting types. Under a change of units the reference rescales by its dimensions and the perturbation by them only when it is absolute: a relative or logarithmic perturbation is a pure number.

perturbed(model, name, perturbation, *, kind) wraps the field of that name on every layer the perturbation’s interval covers, restating the perturbation on each by on_interval, so that one perturbation over the whole mantle perturbs each of its layers; a layer covered only in part is refused by name. It goes through with_field(replace=True), so a model type’s own rule runs: on an Elastic model a perturbed velocity re-derives the Love moduli from it. unperturbed(model) puts every perturbed field’s reference back the same way. A fully 3D field is simply a field that is not perturbed; the reference and the spherical mean of a perturbed field differ in general and are named apart.

planetmodel.perturbation.KINDS = ('absolute', 'relative', 'log')¶: How a reference and its perturbation combine.

class planetmodel.perturbation.PerturbedField(reference: Field, perturbation: Field, *, kind: str, name: str | None = None)[source]¶

Bases: FieldBase

A radial reference of rank 0 combined with a perturbation by kind.

kind is one of KINDS; name is the reference’s unless given. Both fields must lie on the same interval, and the perturbation must have the reference’s character when absolute and be SCALAR otherwise.

property is_radial: bool¶

when the perturbation does.

Type:: Whether the field depends on the radius alone

property kind: str¶: How the two combine, one of KINDS.

on_interval(lo: float, hi: float) → PerturbedField[source]¶: Both parts on [lo, hi].

property perturbation: Field¶: The perturbation of the reference.

radial_pieces() → tuple[ndarray, int] | None[source]¶: The breakpoints and degree that hold both parts exactly, each being held apart: the union of their breakpoints and the higher of their degrees, or None when either part has none.

property reference: Field¶: The radial reference.

renamed(name: str | None) → PerturbedField[source]¶: The same field under name.

rescaled(*, k: float, v: float) → PerturbedField[source]¶: v f(r / k, theta, phi): the reference by (k, v), the perturbation by (k, v) when absolute and by (k, 1) otherwise.

spherical_mean(*, grid: AngularGrid | None = None) → RadialField[source]¶: The mean from the parts: reference + mean(perturbation) when absolute, reference (1 + mean(perturbation)) when relative, and reference mean(exp(perturbation)) when logarithmic, the last by quadrature on grid, or on the perturbation’s own grid where it has weights; a radial perturbation’s mean is itself.

planetmodel.perturbation.perturbed(model: Model, name: str, perturbation: Field, *, kind: str, layers: Iterable[int | str] | None = None) → Model[source]¶

model with the field name perturbed on the layers the perturbation’s interval covers, or on layers when given.

Each layer’s field becomes PerturbedField(field, perturbation restated on the layer, kind=kind), attached through with_field(replace=True). A layer the interval covers only in part, a chosen layer it does not cover, a layer lacking name, and an interval covering no layer are refused by name.

planetmodel.perturbation.unperturbed(model: Model) → Model[source]¶: model with every perturbed field replaced by its reference, through with_field(replace=True).

planetmodel.pushforward module¶

Fields carried across a mapping, forwards and back.

For a mapping m taking the reference body to the physical one, with F[i, j] = d m_i / d X_j and J = det F, a field of rank n and weight w transforms as

T_phys[i1..in](x) = J^-w F[i1, A1] … F[in, An] T_ref[A1..An](X)

with x = m(X): one factor of F on every slot, and a factor of 1/J when the weight is 1. That is the whole rule, driven by the Character and nothing else. Rank 0 weight 1 is rho / J; rank 0 weight 0 composes with the mapping; rank 2 weight 1 carries the second Piola-Kirchhoff stress to the Cauchy stress F S F^T / J; rank 4 weight 1 carries the second elasticity tensor. The pull-back is the same rule with F^-1 in place of F and J^+w in place of J^-w, so the two are inverses at one (F, J). Identical wrapping on every slot preserves the minor and major symmetries, so a Voigt reduction after the rule is faithful.

Two levels. push_forward and pull_back act on arrays: values, F and J already evaluated at the same points, in Cartesian components, at full tensor rank. A Voigt array is a bookkeeping device with no slots to wrap and is refused; full_components is the one place a field’s Voigt presentation is expanded. Complex values pass through the array level, since the rule is linear.

PushedForwardField and PulledBackField act on fields and are lazy: they hold a source and a mapping and do the work in evaluate. Both are asked at reference coordinates and neither inverts the mapping. A pushed-forward field returns the physical tensor at m(X), in Cartesian components or in the spherical frame at m(X), the frame a physical quantity is named in. A pulled-back field takes a callable of the physical point’s spherical coordinates returning components in the frame there, evaluates it at m(X), and gives the referential tensor at X in the frame at X, so it is a field like any other.

class planetmodel.pushforward.PulledBackField(interval: tuple[float, float], physical: SphericalFunction, mapping: Mapping, *, character: Character, name: str | None = None, rtol: float = 1e-09)[source]¶

Bases: FieldBase

A field known on the physical body, stated referentially.

physical(r, theta, phi) is a callable of the physical point’s spherical coordinates returning components in the spherical frame there, full or Voigt for ranks 2 and 4 (a constant tensor is broadcast). At a reference point X the field evaluates physical at m(X), rotates to Cartesian components with the frame at m(X), pulls back with F(X) and J(X), and presents the referential tensor in the frame at X. No inverse mapping is involved. All three coordinates are required and the field is never radial. A formula returning complex values is refused. rescaled conjugates the mapping by k and reads physical at r / k, times v.

is_radial = False¶

property mapping: Mapping¶: The mapping the physical values are pulled back through.

on_interval(lo: float, hi: float) → PulledBackField[source]¶: The same physical formula and mapping on [lo, hi].

property physical: SphericalFunction¶: The formula of the physical point’s spherical coordinates.

renamed(name: str | None) → PulledBackField[source]¶: The same field under name.

rescaled(*, k: float, v: float) → PulledBackField[source]¶: v physical(r / k, theta, phi) under the mapping conjugated by k.

class planetmodel.pushforward.PushedForwardField(field: Field, mapping: Mapping, *, name: str | None = None)[source]¶

Bases: FieldBase

The image of a referential field under a mapping, asked at reference points.

evaluate(r, theta, phi) takes reference coordinates and returns the tensor the field has at the physical point m(X): the source’s Cartesian components are pushed forward with F(X) and J(X), then given in Cartesian components or in the spherical frame at m(X), Voigt-reduced when the character has a Voigt shape. Interval and character are the source’s; the name defaults to <source>_phys where the source is named. All three coordinates are required, F depending on direction even where the source does not, so the field is never radial. rescaled rescales the source and conjugates the mapping by the same length scale.

evaluate(r: ArrayLike, theta: ArrayLike | None, phi: ArrayLike | None, *, frame: str = 'spherical') → ndarray[source]¶

The physical tensor at m(X) for the reference point X = (r, theta, phi).

In Cartesian components for frame="cartesian", in the local spherical frame at m(X) for frame="spherical".

is_radial = False¶

property mapping: Mapping¶: The mapping it is carried across.

on_interval(lo: float, hi: float) → PushedForwardField[source]¶: The source re-stated on [lo, hi], under the same mapping.

renamed(name: str | None) → PushedForwardField[source]¶: The same field under name; None unnames it rather than restoring the <source>_phys default.

rescaled(*, k: float, v: float) → PushedForwardField[source]¶: The source rescaled by (k, v) under the mapping conjugated by k.

property source: Field¶: The referential field carried across.

planetmodel.pushforward.full_components(field: Field, r: ArrayLike, theta: ArrayLike, phi: ArrayLike, *, frame: str = 'cartesian') → ndarray[source]¶

A field’s values at full tensor rank, trailing shape (3,) * rank.

The one place a Voigt-presenting field is expanded: a character with a Voigt shape has its (6,) or (6, 6) values expanded to (3, 3) or (3, 3, 3, 3), anything else is asked plainly.

planetmodel.pushforward.pull_back(values: ArrayLike, F: ArrayLike, J: ArrayLike, character: Character) → ndarray[source]¶

The referential values from the physical ones at the same (F, J).

The inverse of push_forward: F^-1 on every slot and J^+w,

T_AB.. = J^w (F^-1)_Ai (F^-1)_Bj … T_ij..

values are the physical components at x = m(X); the result is referential, at X. Complex values stay complex; a Voigt array is refused.

planetmodel.pushforward.push_forward(values: ArrayLike, F: ArrayLike, J: ArrayLike, character: Character) → ndarray[source]¶

The physical values of a field of character from its referential ones.

values has shape broadcast + (3,) * rank in full Cartesian components at reference points X; F (…, 3, 3) and J (…) are the mapping’s deformation gradient and Jacobian there. Every slot receives one factor of F and a weight-1 field is divided by J:

T_ij.. = J^-w F_iA F_jB … T_AB..

A rank-0 weight-0 field is returned unchanged. Complex values stay complex; a Voigt array is refused.

planetmodel.rheology module¶

Linear viscoelastic rheologies, and a model frozen at one frequency.

A linear viscoelastic body at angular frequency omega is an elastic body with complex moduli, so a solver written for the moduli A, C, F, L, N serves any linear rheology once the model is frozen at a frequency. What a layer’s rheology is, is read from the fields it holds, in the manner of is_fluid and moduli:

viscosity: a Maxwell body in shear, mu(omega) = mu i omega tau / (1 + i omega tau) with tau = viscosity / mu, mu being the Voigt-averaged shear modulus;
qmu, qkappa: a constant-Q absorption band about the reference frequency omega_0, mu(omega) = mu [1 + (2 / pi Q_mu) ln(omega / omega_0) + i / Q_mu], and likewise kappa with Q_kappa; a non-positive Q means no loss.

A transversely isotropic layer is perturbed through its Voigt averages: with delta kappa = (f_kappa - 1) kappa and delta mu = (f_mu - 1) mu, A and C gain delta kappa + 4 delta mu / 3, F gains delta kappa - 2 delta mu / 3, and L and N gain delta mu, which is the isotropic rule when the layer is isotropic and leaves the anisotropy real. A layer holding both a viscosity and Q_mu is a Maxwell body in shear (the shear band is not applied on top) with the bulk band; a layer holding none of the three is elastic and left alone.

frozen(model, omega, *, reference_omega=None) returns the model with the complex moduli of every viscoelastic layer stored under A, C, F, L, N, which moduli reads before anything else, and the frequency recorded as the constant omega. Frequencies are in the model’s units; the reference frequency of the band is the model’s constant omega_ref where it declares one, else 2 pi rad/s, a period of one second, which is PREM’s (reference_omega). is_frozen(model) says whether a model carries that constant. The velocities of a frozen layer, and kappa and mu where it holds them, stay as the elastic description at the reference frequency; the five are what a solver reads.

planetmodel.rheology.RHEOLOGY_NAMES = ('viscosity', 'qmu', 'qkappa')¶: The fields a layer’s rheology is read from.

planetmodel.rheology.dispersive_moduli(layer: LayerLike, omega: float, *, reference_omega: float) → dict[str, Field][source]¶: The complex moduli A, C, F, L, N of a layer at omega under the constant-Q absorption band alone: kappa dispersed by qkappa and mu by qmu about reference_omega, whichever of the two the layer holds, a viscosity ignored. A layer holding neither Q comes back with its static moduli.

planetmodel.rheology.frozen(model: Model, omega: float, *, reference_omega: float | None = None) → Model[source]¶: The model frozen at angular frequency omega: every viscoelastic layer carries its complex moduli under A, C, F, L, N, and the constant omega records the frequency. A model with no viscoelastic layer comes back as a copy with the constant alone; a model that already carries omega is refused, since its moduli are no longer the elastic ones the dispersion starts from. reference_omega defaults to the model’s, see reference_omega.

planetmodel.rheology.frozen_moduli(layer: LayerLike, omega: float, *, reference_omega: float) → dict[str, Field][source]¶: The complex moduli A, C, F, L, N of a viscoelastic layer at omega as fields, from its rheology fields and its elastic moduli; an elastic layer’s moduli come back unchanged. A viscosity makes the shear a Maxwell body; otherwise the Q fields make the constant-Q band of dispersive_moduli.

planetmodel.rheology.is_frozen(model: Model) → bool[source]¶: Whether the model is frozen at a frequency: it carries the constant omega that frozen records.

planetmodel.rheology.is_viscoelastic(layer: LayerLike) → bool[source]¶: Whether a layer holds any of RHEOLOGY_NAMES.

planetmodel.rheology.reference_omega(model: Model) → float[source]¶: The reference angular frequency of the constant-Q band in the model’s units: the constant omega_ref where the model declares one (a deck’s header period, say), else 2 pi rad/s, a period of one second, which is PREM’s.

planetmodel.sampling module¶

A model evaluated on radial times angular nodes.

A consumer’s unknowns live on a product of a radial node set and an angular node set, and a model is evaluated there once, in one vectorised call per layer and field; the result is a Sample.

Layout. Every sampled array is (node, colatitude, longitude, components...), node outermost and longitude fastest. The radial nodes are the per-element GLL nodes of a RadialMesh, flattened element by element, so the flat radius array repeats every element boundary: an interface is two nodes at one radius, one on each side, and each carries its own layer’s one-sided value. That is how a discontinuous field survives sampling without a second array saying where the jumps are. A field that does not depend on direction, one whose layers all declare is_radial, is sampled on (node,) alone and the shape is the mark; a consumer broadcasts. A radial field of any rank is direction-free in this sense, since its components in the local spherical frame are functions of the radius alone.

Components are in the spherical frame (e_r, e_theta, e_phi) at the node, with ranks 2 and 4 in Voigt form where the character has one. The displacement of the geometry’s mapping, m(X) - X, is sampled the same way: a vector on the reference body in the spherical frame at X, which for a radial stretch is (h, 0, 0). It is None when the geometry is the identity.

Missing fields. A name a layer does not hold is refused, naming the layer and the name; missing="nan" fills NaN on that layer’s nodes instead, and fields=None samples the names every layer holds.

Angular grids are node arrays first. kind and lmax say what family the nodes came from and what band the grid resolves; the arrays are the truth, and a consumer with its own grid hands its nodes over as a custom grid. Colatitudes lie in [0, pi], so either pole may be a node, and longitudes in [0, 2 pi). Three families are shipped, each with its colatitude weights:

gauss_legendre(lmax) lmax + 1 Gauss-Legendre colatitudes, no pole driscoll_healy(lmax) 2 lmax + 2 equally spaced from the north

pole, which is a node, to one step short of the south pole (Driscoll and Healy’s grid)

equiangular(ntheta, nphi) equally spaced midpoints, no pole (Fejer’s
first rule), or both poles with poles=True (Clenshaw-Curtis)

The weight convention is

integral over S^2 of f dOmega ~= sum_i w_i sum_j (2 pi / nphi) f_ij,

so weights are colatitude weights against sin(theta) dtheta, the Legendre weights on x = cos(theta) for a Gauss grid, and the longitude rule is the consumer’s own trapezoid, exact for exp(i (m - m') phi) only while |m - m'| < nphi. The band is the highest degree the grid determines, so that a field of that band is recovered exactly by analysis on it: lmax for the Gauss-Legendre and Driscoll-Healy grids, ntheta - 1 for the equiangular midpoints and ntheta - 2 with the poles, and never more than (nphi - 1) // 2, since resolving orders |m| <= lmax needs nphi >= 2 lmax + 1. Every shipped weight rule integrates the polynomials in cos(theta) of degree below ntheta exactly, so the weights integrate every field of the grid’s band exactly. On the Gauss-Legendre and Driscoll-Healy grids they also integrate the product of two fields of the band; the equiangular grids reach their band only through the upsampling analysis of planetmodel.harmonics.

At a pole node the frame (e_theta, e_phi) depends on the longitude, so a component of rank 1 or above is sampled there as the limit along its own longitude’s meridian: the columns differ while the vector does not.

class planetmodel.sampling.AngularGrid(colatitudes: ndarray, longitudes: ndarray, *, kind: str = 'custom', lmax: int | None = None, weights: ndarray | None = None)[source]¶

Bases: object

A product angular node set: colatitudes times longitudes.

colatitudes are strictly increasing in [0, pi], so that either pole may be a node, and longitudes strictly increasing in [0, 2 pi), both in radians and stored read-only. kind names the family the nodes came from and lmax the band the grid resolves, when it is that kind of grid; weights, when present, are colatitude weights in the convention of the module docstring. gauss_legendre, driscoll_healy and equiangular build the shipped families; anything else is a custom grid, which the spherical-harmonic transforms (and so the harmonic lateral rule) refuse. Two grids are equal when their nodes, kind, band and weights are.

colatitudes: ndarray¶

kind: str = 'custom'¶

lmax: int | None = None¶

longitudes: ndarray¶

mean(values: ArrayLike) → ndarray[source]¶

The mean over the sphere of values at the nodes, by the grid’s quadrature: shape (…, ntheta, nphi) to (…).

The colatitude weights against sin(theta) dtheta and, in longitude, the periodic trapezoid, 2 pi / nphi for equally spaced longitudes and half the spacing to either neighbour otherwise; divided by 4 pi. A grid without weights has no quadrature and is refused.

property nphi: int¶: The number of longitudes.

property ntheta: int¶: The number of colatitudes.

property pole_nodes: tuple[bool, bool]¶: Whether the north pole and the south pole are nodes, in that order.

weights: ndarray | None = None¶

planetmodel.sampling.KINDS = ('gauss_legendre', 'driscoll_healy', 'equiangular', 'custom')¶: The families an angular grid can come from.

planetmodel.sampling.MISSING = ('refuse', 'nan')¶: What sample does on a layer lacking a field.

class planetmodel.sampling.Sample(radial: RadialMesh, angular: AngularGrid, fields: Mapping[str, ndarray], displacement: ndarray | None, characters: Mapping[str, Character], dimensions: Mapping[str, Dimensions | None], scales: Scales, layer_names: tuple[str | None, ...])[source]¶

Bases: object

A model’s fields and mapping evaluated on radial times angular nodes.

fields[name] has shape (nnode, ntheta, nphi) + c for a field that depends on direction and (nnode,) + c for one that does not, with c the stored shape of its character: Voigt where the character has one, the component shape otherwise; the components are in the spherical frame at the node. A field a layer lacks is NaN on the nodes of that layer’s elements. displacement is m(X) - X in the spherical frame at X, shape (nnode, ntheta, nphi, 3), or None for an identity geometry. characters and dimensions say what each array is, scales what its numbers are in, and layer_names names the layers the nodes refer to. Every array and mapping the sample holds or hands out is read-only; the radial mesh is the mesh’s own.

angular: AngularGrid¶

characters: Mapping[str, Character]¶

dimensions: Mapping[str, Dimensions | None]¶

displacement: ndarray | None¶

property element_layer: ndarray¶: The skeleton layer of every element, shape (nspec,), read-only.

fields: Mapping[str, ndarray]¶

is_radial(name: str) → bool[source]¶: Whether name was sampled on (node,) alone.

layer_names: tuple[str | None, ...]¶

property nnode: int¶

element boundaries counted twice.

Type:: The flat node count, nspec * ngll

property node_layer: ndarray¶: The skeleton layer of every flat node, shape (nnode,), read-only.

radial: RadialMesh¶

property radius: ndarray¶

the per-element GLL nodes flattened.

Type:: The radius of every flat node

scales: Scales¶

stored_shape(name: str) → tuple[int, ...][source]¶: The trailing axes a field’s array carries beyond the nodes.

planetmodel.sampling.driscoll_healy(lmax: int, *, sampling: int = 1) → AngularGrid[source]¶

The Driscoll-Healy grid resolving degrees and orders up to lmax.

n = 2 lmax + 2 equally spaced colatitudes theta_i = pi i / n, from the north pole, which is a node, to one step short of the south pole, which is not; n equally spaced longitudes from zero for sampling=1, or 2 n for sampling=2, the grid equally spaced in both angles. The weights are Driscoll and Healy’s, exact for the polynomials in cos(theta) of degree below n; the north pole’s is zero.

planetmodel.sampling.equiangular(ntheta: int, nphi: int, *, poles: bool = False) → AngularGrid[source]¶

Equally spaced colatitudes and longitudes, with weights and a band.

By default the colatitudes are the midpoints theta_i = pi (i + 1/2) / ntheta, both poles left out, with the weights of Fejer’s first rule; poles=True makes both poles nodes, theta_i = pi i / (ntheta - 1) for ntheta of at least two, with the Clenshaw-Curtis weights. Either rule is exact for the polynomials in cos(theta) of degree below ntheta. The longitudes are phi_j = 2 pi j / nphi. The band is ntheta - 1 for the midpoints and ntheta - 2 with the poles, and at most (nphi - 1) // 2.

planetmodel.sampling.gauss_legendre(lmax: int, *, nphi: int | None = None) → AngularGrid[source]¶

The Gauss-Legendre grid resolving degrees and orders up to lmax.

lmax + 1 Gauss-Legendre colatitudes in x = cos(theta), sorted so that theta increases, with the Legendre weights alongside, and nphi equally spaced longitudes from zero, 2 lmax + 1 by default, the fewest that resolve the orders. An explicit nphi below that bound is refused: at 2 lmax the orders +lmax and -lmax are one discrete mode and cannot be told apart.

planetmodel.sampling.sample(model: Model, grid: AngularGrid, *, fields: Iterable[str] | None = None, radial: RadialMesh | None = None, ngll: int = 5, drmax: float | None = None, missing: str = 'refuse') → Sample[source]¶

Evaluate a model’s fields and its mapping’s displacement on a grid.

grid is the angular node set. The radial one is radial when the caller has a RadialMesh over the model’s skeleton, else RadialMesh(model.geometry, ngll=ngll, drmax=drmax) when drmax is given, else RadialMesh(model.geometry, ngll=ngll, lmax=grid.lmax) when the grid carries a band; a custom grid without a band needs one of the two said explicitly.

fields is None for the names every layer holds, or a sequence of names. Each is evaluated once per layer on the nodes of that layer’s elements, so the two one-sided values at an interface come from the two layers that own them. A layer lacking a name is refused with KeyError naming the layer and the name, unless missing="nan" asks for NaN on that layer’s nodes. The displacement is sampled unless the geometry is the identity.

planetmodel.skeleton module¶

The reference geometry: ordered shell boundaries.

A Skeleton holds boundary radii b0 < b1 < … < bL defining L concentric layers and answers questions about intervals and membership. It carries no fields, no names and no physics, and knows nothing about units: radii are numbers. Every tolerance is relative to the span bL - b0 and is passed by the caller as rtol.

locate(r) reports boundary ambiguity through Location rather than resolving it: at an interior boundary both adjacent layers are candidates and the choice belongs to the caller.

class planetmodel.skeleton.CoarseningMap(fine: Skeleton, coarse: Skeleton, layers: tuple[tuple[int, ...], ...], boundaries: tuple[int, ...])[source]¶

Bases: object

What a coarsening did: which fine layers each coarse layer merges.

boundaries: tuple[int, ...]¶

coarse: Skeleton¶

property dropped_interfaces: tuple[int, ...]¶: Interior boundaries removed, indexed on the fine skeleton.

fine: Skeleton¶

fine_layer(r: float) → int[source]¶: The fine layer containing r; on a fine boundary, the layer above.

property kept_interfaces: tuple[int, ...]¶: Interior boundaries retained, indexed on the fine skeleton.

layers: tuple[tuple[int, ...], ...]¶

planetmodel.skeleton.DEFAULT_RTOL = 1e-09¶

a fraction of the skeleton’s span (of an interval’s width, for a field).

Type:: The relative tolerance to which radii are the same by default

class planetmodel.skeleton.Location(layers: tuple[int, ...], *, boundary: int | None = None)[source]¶

Bases: object

Where a radius sits: candidate layers plus a boundary flag.

Interior point -> layers=(i,), boundary=None On interior boundary j -> layers=(j-1,j), boundary=j At the innermost or outermost boundary -> layers=(i,), boundary=0 or L

boundary: int | None = None¶

property layer: int¶: The unique layer index; ValueError when two layers are candidates.

layers: tuple[int, ...]¶

class planetmodel.skeleton.Skeleton(boundaries: ArrayLike)[source]¶

Bases: object

Ordered shell boundaries; geometry only.

property boundaries: ndarray¶: The read-only boundary radii b0 < b1 < … < bL.

coarsened(*, keep: Iterable[int] | None = None, drop: Iterable[int] | None = None) → tuple[Skeleton, CoarseningMap][source]¶

A skeleton retaining a subset of the interior boundaries.

Exactly one of keep or drop is given, indexing the interior boundaries as they appear in inner_boundaries (0-based, negatives counting back). The innermost and outermost boundaries are not removable. Returns the coarse skeleton and a CoarseningMap recording which fine layers each coarse layer merges.

contains(lo: float, hi: float, *, rtol: float = 1e-09) → bool[source]¶: Whether [lo, hi] lies inside the skeleton’s span, to rtol*span.

extended(radii: ArrayLike) → Skeleton[source]¶: A skeleton with extra layers appended beyond the outer boundary.

hollowed(radius: float) → Skeleton[source]¶

The skeleton cut at radius from below, which becomes the innermost boundary.

Boundaries at or below the cut are dropped; a cut on an existing boundary, to DEFAULT_RTOL of the span, removes what lies below it, otherwise the innermost surviving layer is shortened. The result is hollow, so its inner boundary is positive.

property inner_boundaries: ndarray¶: The interior boundary radii b1 … b_{L-1}.

interval(i: int) → tuple[float, float][source]¶: The radial interval (b_i, b_{i+1}) of layer i.

property is_hollow: bool¶: Whether the innermost boundary is above the centre.

layer_index(i: int) → int[source]¶

A layer index normalised to 0 <= i < nlayers, negatives counting back.

IndexError outside that range. Every method taking a layer index accepts what this accepts, so -1 is the outermost layer.

locate(r: float, *, rtol: float = 1e-09) → Location[source]¶

Locate radius r, flagging boundary ambiguity through Location.

A radius within rtol*span of a boundary is reported on that boundary. Radii outside the skeleton by more than that raise ValueError.

property nlayers: int¶: The number of layers, L.

refined(radii: ArrayLike) → Skeleton[source]¶: A skeleton with extra interior boundaries inserted.

property span: float¶: bL - b0, the length every relative tolerance refers to.

spans(lo: float, hi: float, *, layer: int | None = None, rtol: float = 1e-09) → bool[source]¶: Whether [lo, hi] is this skeleton’s span, or one layer’s, to rtol*span.

truncated(radius: float) → Skeleton[source]¶

The skeleton cut at radius, which becomes the outer boundary.

Boundaries at or above the cut are dropped; a cut on an existing boundary, to DEFAULT_RTOL of the span, removes what lies above it, otherwise the outermost surviving layer is shortened.

planetmodel.slices module¶

Slices: a model’s field on the points one looks at it by.

A slice is a field evaluated on a set of points chosen for looking at, as plain numbers, for any model: a sphere at one radius (sphere), a great-circle section, distance along the arc by radius (section, to an end= or at an azimuth= for a length=; a whole circle is the plane through the centre), a profile at a point, radius alone (profile), and the spherical mean as a profile (mean_profile). Drawing them is planetmodel.plotting’s business; nothing here draws or needs more than numpy and scipy.

What is sliced. A field by name, from each layer holding it, NaN on a layer that lacks it, or a mapping from layer indices to fields, for a quantity the model does not hold by name. A scalar is sliced as it is; a vector or tensor by component= (an index into its stored components, flattened in C order: the Voigt slot for a Voigt form) and frame= together, one without the other refused, or component="norm" for a vector’s length. relative_to= makes the slice a departure: from the field’s spherical mean at each point’s reference radius ("mean", by means.spherical_mean, consistent with the field’s own rules), from a perturbed field’s own reference ("reference"), or from the same field of another model (a Model), as departure="percent" (100 (v / ref - 1), the default), "difference" (v - ref) or "log10" (log10(v / ref)). Every slice records what it shows (Shown) and labels itself.

Reference or physical. A slice is taken at reference points, where the fields live (at="reference"), and gives besides where the mapping takes each of them, so that it can be drawn as the body is. At physical points (at="physical": a true sphere |x| = r, a true plane, a true radius) each point is pulled back through the mapping’s inverse and the field evaluated there, which is offered where the mapping inverts itself and refused by name where it does not; a physical section or profile also needs the mapping to keep each ray, as a radial stretch does, so that each layer’s part of a ray is one run of radii. Points are given as the library gives them, colatitude and longitude in radians; frames.from_geographic converts latitudes and longitudes.

Spectra. spectrum analyses a sphere slice into the library’s real orthonormal harmonics (planetmodel.harmonics: the integral of Y_lm^2 over the unit sphere is one, with no Condon-Shortley phase), through the harmonics extra, and gives the power by degree

P_l = sum_m c_lm^2 / (4 pi),

so that the sum over degrees is the field’s mean square over the sphere. The slice is taken on gauss_legendre(2 lmax) unless a grid is given, so that a field not limited to the band does not alias into it. correlation gives, at each of a list of radii, two fields’ correlation by degree, r_l = sum_m a_lm b_lm / sqrt(sum_m a_lm^2 sum_m b_lm^2), and over the degrees from 1 to lmax, NaN where a degree has no power in either (its amplitude below NO_POWER of the field’s); the fields may be two models’ or two of one model’s.

Interfaces. A section or a profile keeps both sides of an interface, each layer’s numbers including its own ends. A sphere slice is asked for one side, side="below" or "above", required where its radius lies on an interface to the geometry’s tolerance and of no effect elsewhere, so that a loop over radii names it once; a physical sphere’s points that fall on an interface take the side the same way, and its points outside the body are NaN.

planetmodel.slices.AT = ('reference', 'physical')¶

in the reference body, or the physical one.

Type:: Where a slice’s points are

class planetmodel.slices.Correlation(radii: ndarray, degrees: ndarray, by_degree: ndarray, total: ndarray, first: tuple[Spectrum, ...], second: tuple[Spectrum, ...])[source]¶

Bases: object

Two fields’ correlation by degree and radius.

radii (k,) and degrees (lmax + 1,); by_degree (k, lmax + 1), r_l at each radius; total (k,), over the degrees from 1 to lmax; first and second the two fields’ spectra at each radius.

by_degree: ndarray¶

degrees: ndarray¶

first: tuple[Spectrum, ...]¶

property label: str¶

“correlation of vs”, or “correlation of vs with rho”.

Type:: What is correlated with what, by name

radii: ndarray¶

second: tuple[Spectrum, ...]¶

total: ndarray¶

planetmodel.slices.DEPARTURES = ('percent', 'difference', 'log10')¶: How a departure is expressed.

planetmodel.slices.NO_POWER = 1e-10¶: The amplitude of a degree, relative to the whole field’s, below which a correlation takes the degree to have no power: what the analysis gives a field with nothing at that degree.

class planetmodel.slices.ProfileSlice(theta: float | None, phi: float | None, layers: tuple[SliceLayer, ...], shown: Shown, *, at: str = 'reference')[source]¶

Bases: object

A field along a radius at (theta, phi), layer by layer; for a mean profile theta and phi are None and the values are the spherical means.

at: str = 'reference'¶

property label: str¶: The label of what the slice shows.

layers: tuple[SliceLayer, ...]¶

phi: float | None¶

shown: Shown¶

theta: float | None¶

planetmodel.slices.SIDES = ('below', 'above')¶: The sides of an interface a sphere slice may be taken on.

class planetmodel.slices.SectionSlice(distance: ndarray, path: ndarray, plane: ndarray, layers: tuple[SliceLayer, ...], shown: Shown, *, at: str = 'reference')[source]¶

Bases: object

A field on a great-circle section, layer by layer.

distance (m,) is the angle along the arc from its start, in radians; path (m, 3) the unit vectors of the arc’s points; plane (2, 3) the start and the direction of travel there, an orthonormal basis of the section’s plane, for drawing positions in it; layers each layer’s points and values.

at: str = 'reference'¶

distance: ndarray¶

property label: str¶: The label of what the slice shows.

layers: tuple[SliceLayer, ...]¶

property length: float¶: The arc’s length, in radians.

path: ndarray¶

plane: ndarray¶

shown: Shown¶

Bases: object

What a slice’s values are: the field, the component, the departure.

name is the field’s name (None for fields given by layer); component an index into the stored components or "norm", and frame the frame of a component; relative_to what a departure is from ("mean", "reference", or the other model’s name) and departure how it is expressed; statistic "spherical mean" for a mean profile; dimensions and scales the field’s and the model’s, so that the label states the unit.

component: int | str | None = None¶

departure: str | None = None¶

dimensions: Dimensions | None = None¶

frame: str | None = None¶

property is_departure: bool¶: Whether the values are a departure, centred on zero.

property label: str¶

“vs, percent from the mean [%]”, “rho [kg m-3]”, “|u| [m]”.

Type:: A label for an axis or a colour bar

name: str | None¶

relative_to: str | None = None¶

scales: Scales = Scales.SI¶

statistic: str | None = None¶

class planetmodel.slices.SliceLayer(index: int, radius: ndarray, values: ndarray, reference: ndarray, physical: ndarray)[source]¶

Bases: object

One layer’s part of a section or a profile.

radius holds the points’ radii, reference or physical as the slice’s at, and values the field there, both of the points’ shape: (n, m), n radii by m points along the arc, for a section, and (n,) for a profile; reference and physical are the points’ Cartesian positions, that shape followed by 3. Each layer’s points include its own ends, so that an interface carries both sides.

index: int¶

physical: ndarray¶

radius: ndarray¶

reference: ndarray¶

values: ndarray¶

class planetmodel.slices.Spectrum(radius: float, degrees: ndarray, power: ndarray, coefficients: ndarray, shown: Shown, *, at: str = 'reference', side: str | None = None)[source]¶

Bases: object

A field’s power by degree on a sphere.

degrees 0 .. lmax; power P_l, summing to the mean square over the sphere; coefficients the real orthonormal ones, (2, lmax + 1, lmax + 1) as planetmodel.harmonics lays them out; radius, side and at as the sphere slice they were analysed from.

property amplitude: ndarray¶: The square root of the power, by degree.

at: str = 'reference'¶

coefficients: ndarray¶

degrees: ndarray¶

property label: str¶: The label of what the spectrum is of.

power: ndarray¶

radius: float¶

shown: Shown¶

side: str | None = None¶

class planetmodel.slices.SphereSlice(radius: float, grid: AngularGrid, values: np.ndarray, reference: np.ndarray, physical: np.ndarray, shown: Shown, *, at: str = 'reference', side: str | None = None)[source]¶

Bases: object

A field on a sphere: values (ntheta, nphi) on grid at radius.

radius is a reference radius, or a physical one when at is "physical"; reference and physical are the points’ Cartesian positions (ntheta, nphi, 3) in the two bodies; side the side of an interface the slice was taken on, where it lies on one.

at: str = 'reference'¶

grid: AngularGrid¶

property label: str¶: The label of what the slice shows.

physical: np.ndarray¶

radius: float¶

reference: np.ndarray¶

shown: Shown¶

side: str | None = None¶

values: np.ndarray¶

planetmodel.slices.correlation(model: Model, fields: Fields, other: Model, *, other_fields: Fields | None = None, lmax: int, radii: ArrayLike | None = None, depths: ArrayLike | None = None, grid: AngularGrid | None = None, side: str | None = None, at: str = 'reference', component: int | str | None = None, frame: str | None = None, relative_to: RelativeTo = None, departure: str | None = None, mean_grid: AngularGrid | None = None) → Correlation[source]¶: Two fields’ correlation by degree and radius, by the rules of the module docstring: fields of model against other_fields (the same by default) of other, which may be model itself, at each of radii or depths; both sliced with the same choices.

planetmodel.slices.evaluate(model: Model, fields: Fields, layer: int, r: ArrayLike, theta: ArrayLike, phi: ArrayLike, *, component: int | str | None = None, frame: str | None = None, relative_to: RelativeTo = None, departure: str | None = None, mean_grid: AngularGrid | None = None) → np.ndarray[source]¶: What a slice of the same choices shows, at reference points (r, theta, phi) of a layer, broadcast: for a set of points of the caller’s own, as the VTK writer’s departures are; NaN where the layer lacks the field.

planetmodel.slices.mean_profile(model: Model, fields: Fields, *, n: int = 100, mean_grid: AngularGrid | None = None) → ProfileSlice[source]¶: The spherical mean of a scalar field along the radius, n points per layer, its ends included; a field without a mean of its own takes it on mean_grid.

planetmodel.slices.profile(model: Model, fields: Fields, theta: float, phi: float, *, n: int = 100, radii: tuple[float, float] | None = None, at: str = 'reference', component: int | str | None = None, frame: str | None = None, relative_to: RelativeTo = None, departure: str | None = None, mean_grid: AngularGrid | None = None) → ProfileSlice[source]¶: A field along the radius at (theta, phi), by the rules of the module docstring: n points per layer, its ends included, within radii (lo, hi), the whole body by default.

planetmodel.slices.section(model: Model, fields: Fields, start: tuple[float, float], *, end: tuple[float, float] | None = None, azimuth: float | None = None, length: float | None = None, radii: tuple[float, float] | None = None, n: int = 50, m: int = 181, at: str = 'reference', component: int | str | None = None, frame: str | None = None, relative_to: RelativeTo = None, departure: str | None = None, mean_grid: AngularGrid | None = None) → SectionSlice[source]¶

A field on a great-circle section, by the rules of the module docstring.

The arc runs from start (theta, phi) along the minor arc to end, or from start at azimuth (from north toward east, in radians) for length radians, up to a whole circle, the one or the other given; radii (lo, hi) bounds it, the whole body by default; n points per layer, its ends included, and m along the arc; at, the component, frame and departure as the module docstring says.

planetmodel.slices.spectrum(model: Model, fields: Fields, *, lmax: int, radius: float | None = None, depth: float | None = None, grid: AngularGrid | None = None, side: str | None = None, at: str = 'reference', component: int | str | None = None, frame: str | None = None, relative_to: RelativeTo = None, departure: str | None = None, mean_grid: AngularGrid | None = None) → Spectrum[source]¶: A field’s power by degree to lmax on a sphere, by the rules of the module docstring; the sphere and the field’s choices as sphere takes them. A slice with NaN in it (a layer lacking the field, a physical sphere leaving the body) is refused.

A field on a sphere, by the rules of the module docstring.

Exactly one of radius and depth, the depth below the reference body’s outer boundary; grid the angular grid, equiangular at one degree by default (180 x 360 midpoints); side the side of an interface the radius lies on; at reference or physical points; the component, frame and departure as the module docstring says, the mean of a field without one of its own taken on mean_grid or, where it has weights, on grid.

planetmodel.testing module¶

Executable contracts for the protocols of planetmodel.

Each check_* takes a candidate object, exercises what its protocol promises, and raises AssertionError naming the first violation. The shipped implementations are held to these functions in the test suite; an implementation written elsewhere is checked by the same call. The imports of the rest of the package happen inside the functions.

planetmodel.testing.check_displacement(h: RadialDisplacement, skeleton: Skeleton, *, n: int = 33) → None[source]¶

Hold a RadialDisplacement to its protocol over a skeleton’s span.

Checks finite values of the broadcast shape; single-valuedness on the polar axis; continuity across every declared knot inside the skeleton; radial_derivative, when present, against a central difference away from the knots and free of undeclared kinks; the shape of angular_gradient; and, when present, scaled(k) against k h(r / k, theta, phi), with its knots and derivatives. Every tolerance is relative to the skeleton’s span or to the value.

planetmodel.testing.check_field(field: Field, *, rng: np.random.Generator | None = None, n: int = 64) → None[source]¶

Hold a Field to its protocol.

Checks the three attributes; that a call with the radius alone is accepted exactly when the field is radial and of rank 0; float64 (or complex128) values of the broadcast shape plus the stored shape, agreeing with dtype; both ends of the interval; refusal outside it and of an unknown frame; the two frames related by the rotation the character implies (a factor of R on every slot, expanded from Voigt where the field is Voigt); evaluate_at against evaluate; on_interval on the overlap; and the rescaled round trip.

planetmodel.testing.check_geometry(geometry: Geometry) → None[source]¶

Hold a Geometry to its invariants and its numbering.

Checks that the invariants hold (by reconstructing the geometry with checks on), that names are unique, that the interfaces sit on the skeleton’s boundaries with consistent between, that scaled(k) round-trips, and that the mapping passes check_mapping on the validity lattice’s angles, at radii inside every span between the boundaries and the declared knots.

planetmodel.testing.check_layer_function(fn: LayerFunction, *, n: int = 17, rtol: float = 1e-06) → None[source]¶

Hold a LayerFunction to its protocol on its own interval.

Checks finite float64 (or complex128) values of the argument’s shape; derivative against a central difference; integrate against a fine trapezoid rule and its sign convention; on_interval agreeing on the overlap; and rescaled against v f(r / k) and its round trip. Every tolerance is relative to the value or the interval’s width.

planetmodel.testing.check_mapping(m: Mapping, points: ArrayLike, *, rtol: float = 1e-06, step: float | None = None) → None[source]¶

Hold a Mapping to its protocol on the given Cartesian points.

Checks the shapes; deformation_gradient against a central difference of m with F[i, j] = d m_i / d X_j; jacobian against det F; and, when present, displacement against m(X) - X, inverse as a round trip, and scaled(k) against k m(X / k) with the same F and J. The difference step is 1e-6 of each point’s radius unless step gives one step for every point.

planetmodel.testing.check_model(model: Model, *, rng: np.random.Generator | None = None) → None[source]¶

Hold a Model to its invariants.

Checks that replaced with checks on accepts it and keeps its class and fields; that every layer’s fields pass check_field, sit on the layer’s interval, and have their spec’s character; that layer(name) and layer(i) agree; that every constant is its SI value over the scales’ factor; that conversion to other scales and back restores every field to 1e-12; and that every surgery keeps the class and validates again.

planetmodel.testing.check_model_file(path: str | os.PathLike[str], model: Model, *, rng: np.random.Generator | None = None, n: int = 48, rtol: float = 1e-12) → None[source]¶

Hold a model file to the model it was written from.

Checks the structure (validate_structure); the skeleton, the layer and interface names and the scales; that every variable is held on the layers its layers attribute names, NaN elsewhere, by fields of the model of the stated character; that a lateral variable names in radial_layers exactly the layers where the model’s field is radial, and reads back radial there; that its values at n random entries of the file’s own nodes are the model’s there, the reference and the perturbation of a perturbed field apart, and so are its polar values; that the geometry group holds the mapping’s reliefs at its nodes, and a displacement variable the mapping’s displacement; and that the model read back from the file takes the stored values at those nodes, for every variable read by a local rule. Tolerances are rtol of a variable’s largest value. netCDF4 is imported here.

planetmodel.testing.check_refinement(refinement: Refinement, points: ArrayLike) → None[source]¶

Hold a refinement to its protocol on the given points.

Checks size and far_size as a positive floor and ceiling; that the (n, d) points give float64 sizes of shape (n,), finite and within [size, far_size]; and that a batch agrees with its points taken one at a time, so vectorisation is a promise rather than an optimisation. The points are the caller’s probes and should span the domain the refinement will be asked about, the refined region included.

planetmodel.testing.check_sample(sample: Sample, model: Model, *, rng: np.random.Generator | None = None, n: int = 64, rtol: float = 1e-12) → None[source]¶

Hold a Sample to its layout and to the model it was taken from.

Layout: the radial mesh lies over the model’s skeleton and the nodes are its per-element GLL nodes flattened; the angular nodes increase, colatitudes in [0, pi] and longitudes in [0, 2 pi); every field is float64 and read-only, shaped (nnode,) + c exactly when every layer holding it declares itself radial and (nnode, ntheta, nphi) + c otherwise, with c the stored shape of its character. Metadata: the characters and dimensions are those of the model’s fields and specs, the scales and layer names the model’s. Values: NaN on every node of a layer lacking the field and finite everywhere else; n random (node, colatitude, longitude) entries re-evaluated from the node’s own layer, so the two nodes at an interface are two questions. The displacement is None for an identity geometry and otherwise (nnode, ntheta, nphi, 3), read-only, finite, and R^T (m(X) - X) at n random entries. Tolerances are rtol of the field’s largest value, or of the outer radius for the displacement.

planetmodel.testing.check_spectral_basis(basis: SpectralBasis, *, rng: np.random.Generator | None = None, rtol: float = 1e-10) → None[source]¶

Hold a SpectralBasis to what it promises of one degree’s eigenbasis.

The spectrum: theta read-only, ascending, at least 1, nmodes of them. The modes: read-only, of shape (nglob, nmodes), M-orthonormal on the full mesh, each an eigenfunction of the family’s operator with its theta, and zero at the centre of a ball when l >= 1. The transforms: the physical synthesis is the full one read on the restriction’s nodes; analyse inverts the full synthesis and is idempotent on arbitrary nodal values, for one vector and for a stack. Evaluation: evaluate at the physical nodes reproduces the modes there, keeps the shape of its argument and refuses a radius outside the interval. Integration: weights are non-negative, the family’s mass inside the physical interval and no more than it at the two ends (where the mass counts the pad elements too), and, on a mesh of at least three GLL nodes, sum to the measure of the physical interval, (r2^3 - r1^3) / 3 or r2 - r1. pointwise_variance is the diagonal of the synthesised covariance. white_noise has the shapes it promises and unit covariance on the leading modes to six standard errors. rtol is relative to the largest kept eigenvalue where an operator is applied and to one elsewhere.

planetmodel.testing.check_spherical_basis(basis: SphericalBasis, *, rng: np.random.Generator | None = None, rtol: float = 1e-10) → None[source]¶

Hold a SphericalBasis to the order it fixes and to its degrees.

The bases are those of degrees 0 … lmax on one family and one restriction, with nmodes and size counting them. The blocks of the harmonics, taken in the order of harmonics.packing, tile range(size) without gap or overlap, and theta on each block is its degree’s. synthesise fills each harmonic of the coefficient array with its degree’s synthesis of that block and leaves the entries of no harmonic zero, physical values being the full ones on the restriction’s nodes; analyse inverts the full synthesis, and ignores the entries of no harmonic.

planetmodel.units module¶

Physical dimensions, and the scales that remove them.

Data arrives with known dimensions, is non-dimensionalised once into O(1) numbers, is worked with in that form, and is re-dimensionalised on output. The two pieces here make each step checkable rather than remembered. Nothing inside a field or a mesh consults them: a model holds a Scales and a table of dimensions by name, and conversion is the model’s operation.

Dimensions are physical: integer exponents of (mass, length, time, temperature). They are independent of a field’s character, which is its law under a mapping: Q_kappa and a relaxation time are both invariant, yet one is dimensionless and the other is seconds.

Scales say what one unit of each base dimension is, in SI. A quantity of dimensions d stored under scales s relates to its SI value by

stored = SI / s.factor(d), s.factor(d) = mass^d.mass length^d.length
time^d.time temperature^d.temperature

so Scales.SI is the identity. Scales.geophysical(a, density=rho) takes the length a and the mass rho a^3 and fixes the time by G rho T^2 = 1, so that G is one in those units: factor(G’s dimensions) is G_SI to machine precision, the time scale passing through a square root; the temperature scale is a kelvin unless it is given. A solver that takes G as a parameter is handed G_SI / scales.factor(GRAVITATIONAL_CONSTANT).

The dimension constants are module names. DENSITY here is the dimensions kg m^-3 and DENSITY in planetmodel.character is a transformation law; a module using both refers to one of them through its module.

planetmodel.units.BASE = ('mass', 'length', 'time', 'temperature')¶: The base dimensions, in the order of every exponent tuple and scale.

class planetmodel.units.Dimensions(*, mass: int = 0, length: int = 0, time: int = 0, temperature: int = 0)[source]¶

Bases: object

Physical dimensions as integer exponents of (mass, length, time, temperature).

Multiplying quantities adds exponents, so MODULUS == DENSITY * VELOCITY**2 holds here as it does on paper.

property is_dimensionless: bool¶: Whether every exponent is zero.

length: int = 0¶

mass: int = 0¶

temperature: int = 0¶

time: int = 0¶

unit_string(*, si: bool = True) → str[source]¶

A udunits-style unit string: “kg m-3” in SI, “1” otherwise.

A non-dimensional model’s numbers are pure whatever their physical dimensions, so si=False gives “1” for every quantity.

planetmodel.units.EARTH_MEAN_DENSITY = 5515.0¶: The conventional mean density of the Earth, kg m^-3. Prescribed, not computed: a calculated mean would shift with every crustal edit.

planetmodel.units.G_SI = 6.6743e-11¶: The gravitational constant, m^3 kg^-1 s^-2 (CODATA 2018). The one definition in the package; everything that needs G reads it here or from a model’s constants.

class planetmodel.units.Scales(*, length: float = 1.0, mass: float = 1.0, time: float = 1.0, temperature: float = 1.0)[source]¶

Bases: object

What one unit of each base dimension is, in SI.

A quantity of dimensions d stored under these scales relates to its SI value by stored = SI / factor(d). Scales.SI is the identity and the default everywhere: a model is SI until it says otherwise.

SI: ClassVar[Scales] = Scales.SI¶

factor(dims: Dimensions) → float[source]¶: The SI size of one stored unit of a quantity with dims.

classmethod geophysical(length: float, *, density: float = 5515.0, temperature: float = 1.0) → Scales[source]¶

Scales for self-gravitating problems: G becomes one.

length is the body’s radius in metres and density a prescribed reference density in kg m^-3, conventionally the Earth’s mean. The mass scale is density length^3 and the time scale is fixed by G density T^2 = 1; temperature is the temperature scale in kelvin, one by default.

property is_si: bool¶: Whether one stored unit of each base dimension is one SI unit.

length: float = 1.0¶

mass: float = 1.0¶

temperature: float = 1.0¶

time: float = 1.0¶

planetmodel.units.scaled_unit(dims: Dimensions | None, scales: Scales) → str[source]¶: The stored unit as a UDUNITS scaled unit: the SI unit, preceded by the size of one stored unit when the scales make it other than one (“m s-1” in SI, “3865.300165 m s-1” under the geophysical scales of the Earth); “1” for a pure number and “unknown” where no dimensions are declared.

planetmodel.units.unit_string(dims: Dimensions | None, *, si: bool = True) → str[source]¶: dims.unit_string(si=si), or “unknown” where no dimensions are declared.

planetmodel.vocabulary module¶

The field names the library itself speaks, and the constants it knows.

A model may hold a field under any name. What this table fixes is what a name means when the library reads, writes or guarantees a field by it: a model looks its names up here, refuses a field attached under a known name with the wrong character, and converts a field’s units by the dimensions the name carries. Each entry gives the name’s character (its law under a mapping), its physical dimensions, and what it is; a model supplies a FieldSpec of its own for a name outside the table.

Character and dimensions are independent. The rule for the character: a density is rank 0 weight 1; a modulus is a coefficient of the second elasticity tensor, whose density-like factor of 1/J it shares, so every modulus (kappa, mu, and the Love moduli A, C, F, L, N) is rank 0 weight 1 too, and the tensor itself is rank 4 weight 1; a velocity, a quality factor, an anisotropy ratio and a viscosity are rank 0 weight 0. A velocity is the square root of a modulus over a density, the weights cancelling, which is why it has no law of its own.

CONSTANTS holds the named constants a model reads in its own units, each declared with its SI value and dimensions.

planetmodel.vocabulary.CONSTANTS: dict[str, Constant] = {'G': Constant(value_si=6.6743e-11, dimensions=Dimensions(mass=-1, length=3, time=-2, temperature=0), meaning='the gravitational constant')}¶: The named constants, each with its SI value and dimensions.

class planetmodel.vocabulary.Constant(value_si: float, dimensions: Dimensions, *, meaning: str = '')[source]¶

Bases: object

A named constant: its SI value, its dimensions and a gloss.

dimensions: Dimensions¶

meaning: str = ''¶

value_si: float¶

class planetmodel.vocabulary.FieldSpec(character: Character, dimensions: Dimensions | None, *, meaning: str = '')[source]¶

Bases: object

What a field name means: its character, dimensions and a gloss.

dimensions is None for a name whose units are not declared; a model refuses to convert such a field by name.

character: Character¶

dimensions: Dimensions | None¶

meaning: str = ''¶

planetmodel.vocabulary.VOCABULARY: dict[str, FieldSpec] = {'A': FieldSpec(character=Character(rank=0, weight=1, voigt=True), dimensions=Dimensions(mass=1, length=-1, time=-2, temperature=0), meaning='Love modulus A = rho vph^2'), 'C': FieldSpec(character=Character(rank=0, weight=1, voigt=True), dimensions=Dimensions(mass=1, length=-1, time=-2, temperature=0), meaning='Love modulus C = rho vpv^2'), 'F': FieldSpec(character=Character(rank=0, weight=1, voigt=True), dimensions=Dimensions(mass=1, length=-1, time=-2, temperature=0), meaning='Love modulus F = eta (A - 2L)'), 'L': FieldSpec(character=Character(rank=0, weight=1, voigt=True), dimensions=Dimensions(mass=1, length=-1, time=-2, temperature=0), meaning='Love modulus L = rho vsv^2'), 'N': FieldSpec(character=Character(rank=0, weight=1, voigt=True), dimensions=Dimensions(mass=1, length=-1, time=-2, temperature=0), meaning='Love modulus N = rho vsh^2'), 'drho_dphi0': FieldSpec(character=Character(rank=0, weight=1, voigt=True), dimensions=Dimensions(mass=1, length=-5, time=2, temperature=0), meaning="stratification of a fluid layer: d rho / d Phi_0 = rho'(r) / g(r)"), 'elastic_moduli': FieldSpec(character=Character(rank=4, weight=1, voigt=True), dimensions=Dimensions(mass=1, length=-1, time=-2, temperature=0), meaning='the second elasticity tensor, static, Voigt-reduced in the spherical frame'), 'eta': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=0, time=0, temperature=0), meaning='the anisotropy parameter F / (A - 2L)'), 'g': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=1, time=-2, temperature=0), meaning='gravitational acceleration of the spherically symmetric reference, radially inward, by magnitude'), 'kappa': FieldSpec(character=Character(rank=0, weight=1, voigt=True), dimensions=Dimensions(mass=1, length=-1, time=-2, temperature=0), meaning='bulk modulus of an isotropic medium'), 'mu': FieldSpec(character=Character(rank=0, weight=1, voigt=True), dimensions=Dimensions(mass=1, length=-1, time=-2, temperature=0), meaning='shear modulus of an isotropic medium'), 'qkappa': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=0, time=0, temperature=0), meaning='bulk quality factor; zero attenuates nothing'), 'qmu': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=0, time=0, temperature=0), meaning='shear quality factor; zero attenuates nothing (a fluid)'), 'rho': FieldSpec(character=Character(rank=0, weight=1, voigt=True), dimensions=Dimensions(mass=1, length=-3, time=0, temperature=0), meaning='mass density'), 'temperature': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=0, time=0, temperature=1), meaning='absolute temperature'), 'viscosity': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=1, length=-1, time=-1, temperature=0), meaning='dynamic viscosity'), 'vp': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=1, time=-1, temperature=0), meaning='P-wave speed of an isotropic medium'), 'vph': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=1, time=-1, temperature=0), meaning='P-wave speed normal to the symmetry axis (horizontal)'), 'vpv': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=1, time=-1, temperature=0), meaning='P-wave speed along the symmetry axis (vertical)'), 'vs': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=1, time=-1, temperature=0), meaning='S-wave speed of an isotropic medium'), 'vsh': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=1, time=-1, temperature=0), meaning='S-wave speed polarised normal to the symmetry axis (horizontal)'), 'vsv': FieldSpec(character=Character(rank=0, weight=0, voigt=True), dimensions=Dimensions(mass=0, length=1, time=-1, temperature=0), meaning='S-wave speed polarised along the symmetry axis (vertical)')}¶: The shipped names, each mapped to its spec.

planetmodel.vocabulary.spec(name: str) → FieldSpec[source]¶: The spec of a shipped name; KeyError, naming it, for any other.

Module contents¶

planetmodel: spherically layered planetary models, and their meshes.

A model is a skeleton of boundary radii, a geometry that places the skeleton in the physical world through one continuous mapping, and a bag of fields on each layer. This package holds the skeleton and the geometry, the mappings that place them, the fields and the model, the catalogue of named models, and the meshers that hand a model to a solver. Every internal is numbers: the model alone carries units.

Modules reached by name, one level down:

frames the local spherical frame, tensors moved between frames, Voigt materials velocity and moduli conversions, the elastic field, what a

layer’s fields imply (fluidity, moduli)

rheology what a layer’s rheology fields imply, and a model frozen at
a frequency with complex moduli

units Dimensions, Scales and the dimension constants vocabulary the shipped field names, their characters and dimensions behaviours the free functions as methods: layer_method and the mixins catalogue PREM, LayeredIsotropicElastic and MineosModel, model types deck tabulated models: decks, their formats, fields by interpolation sampling a model on a radial mesh times an angular grid lateral values between the nodes of an angular grid: the lateral rules modelfile the model file: a model in one netCDF file (netCDF4, optional) harmonics real spherical harmonics; grid transforms via ducc0 (optional) slices a field on a sphere, a section or a profile; spectra and

correlations by degree

plotting profiles, sections, maps and spectra drawn (matplotlib and
cartopy, optional), and VTK files for ParaView

mesh1d radial spectral-element meshes (GLL), nodal values, gravity mesh3d 2D and 3D meshes via gmsh, and MFEM export randomfield Matern random fields on balls, annuli and layers testing the executable contracts, check_field and its kin