# Source code for boule._triaxialellipsoid

# Copyright (c) 2019 The Boule Developers.
#
# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
#
"""
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
get in touch <https://www.fatiando.org/contact>__.

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:

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
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
)