Source code for boule._triaxialellipsoid

# Copyright (c) 2019 The Boule Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
# This code is part of the Fatiando a Terra project (
Define a reference triaxial ellipsoid.
from warnings import warn

import attr
import numpy as np

# Don't let ellipsoid parameters be changed to avoid messing up calculations
# accidentally.
[docs]@attr.s(frozen=True) class TriaxialEllipsoid: r""" A rotating triaxial ellipsoid. The ellipsoid is defined by five parameters: semimajor axis, semimedium axis, semiminor axis, geocentric gravitational constant, and angular velocity The thee semi-axis are different and the ellipsoid spins around it's largest moment of inertia. **This class is read-only:** Input parameters and attributes cannot be changed after instantiation. **Units:** All input parameters and derived attributes are in SI units. .. attention:: Gravity calculations have not been implemented yet for triaxial ellipsoids. If you're interested in this feature or would like to help implement it, please `get in touch <>`__. Parameters ---------- name : str A short name for the ellipsoid, for example ``"WGS84"``. semimajor_axis : float The semimajor (largest) axis of the ellipsoid. Definition: :math:`a`. Units: :math:`m`. semimedium_axis : float The semimedium (middle) axis of the ellipsoid. Definition: :math:`b`. Units: :math:`m`. semiminor_axis : float The semiminor (smallest) axis of the ellipsoid. Definition: :math:`c`. Units: :math:`m`. geocentric_grav_const : float The geocentric gravitational constant. The product of the mass of the ellipsoid :math:`M` and the gravitational constant :math:`G`. Definition: :math:`GM`. Units: :math:`m^3.s^{-2}`. angular_velocity : float The angular velocity of the rotating ellipsoid. Definition: :math:`\omega`. Units: :math:`\\rad.s^{-1}`. long_name : str or None A long name for the ellipsoid, for example ``"World Geodetic System 1984"`` (optional). reference : str or None Citation for the ellipsoid parameter values (optional). Examples -------- We can define an ellipsoid by setting the 5 key numerical parameters: >>> ellipsoid = TriaxialEllipsoid( ... name="VESTA", ... long_name="Vesta Triaxial Ellipsoid", ... semimajor_axis=286_300, ... semimedium_axis=278_600, ... semiminor_axis=223_200, ... geocentric_grav_const=1.729094e10, ... angular_velocity=326.71050958367e-6, ... reference=( ... "Russell, C. T., Raymond, C. A., Coradini, A., McSween, " ... "H. Y., Zuber, M. T., Nathues, A., et al. (2012). Dawn at " ... "Vesta: Testing the Protoplanetary Paradigm. Science. " ... "doi:10.1126/science.1219381" ... ), ... ) >>> print(ellipsoid) # doctest: +ELLIPSIS TriaxialEllipsoid(name='VESTA', ...) >>> print(ellipsoid.long_name) Vesta Triaxial Ellipsoid The class then defines several derived attributes based on the input parameters: >>> print(f"{ellipsoid.mean_radius:.0f} m") 262700 m >>> print(f"{ellipsoid.volume * 1e-9:.0f} km³") 74573626 km³ """ name = attr.ib() semimajor_axis = attr.ib() semimedium_axis = attr.ib() semiminor_axis = attr.ib() geocentric_grav_const = attr.ib() angular_velocity = attr.ib() long_name = attr.ib(default=None) reference = attr.ib(default=None) def _raise_invalid_axis(self): "Raise a ValueError informing that the axis are invalid." raise ValueError( "Invalid triaxial ellipsoid axis: " f"major={self.semimajor_axis} " f"medium={self.semimedium_axis} " f"minor={self.semiminor_axis}. " "Must be major > medium > minor." ) @semimajor_axis.validator def _check_semimajor_axis(self, semimajor_axis, value): """ Check if semimajor_axis is positive and is the largest of the axis. """ if not value > 0: raise ValueError( f"Invalid semi-major axis '{value}'. Should be greater than zero." ) if self.semimedium_axis > value: self._raise_invalid_axis() @semimedium_axis.validator def _check_semimedium_axis(self, semimedium_axis, value): """ Check if semimedium_axis is positive and larger than semi-minor axis. """ if not value > 0: raise ValueError( f"Invalid semi-medium axis '{value}'. Should be greater than zero." ) if self.semiminor_axis > value: self._raise_invalid_axis() @semiminor_axis.validator def _check_semiminor_axis(self, semiminor_axis, value): "Check if semiminor_axis is positive." if not value > 0: raise ValueError( f"Invalid semi-minor axis '{value}'. Should be greater than zero." ) # Don't need to check here because if the two checks for major and # medium pass it means that this is the smallest. @geocentric_grav_const.validator def _check_geocentric_grav_const(self, geocentric_grav_const, value): "Warn if geocentric_grav_const is negative." if value < 0: warn(f"The geocentric gravitational constant is negative: '{value}'") @property def mean_radius(self): r""" The arithmetic mean radius of the ellipsoid. Definition: :math:`R = \dfrac{a + b + c}{3}`. Units: :math:`m`. """ return (self.semimajor_axis + self.semimedium_axis + self.semiminor_axis) / 3 @property def volume(self): r""" The volume bounded by the ellipsoid. Definition: :math:`V = \dfrac{4}{3} \pi a b c`. Units: :math:`m^3`. """ return ( (4 / 3 * np.pi) * self.semimajor_axis * self.semimedium_axis * self.semiminor_axis )