# Source code for boule._sphere

# Copyright (c) 2019 The Boule Developers.
#
# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
#
"""
Define the reference sphere (ellipsoid with 0 flattening).
"""
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 Sphere:
r"""
A rotating sphere (zero-flattening ellipsoid).

The ellipsoid is defined by three parameters: radius, geocentric
gravitational constant, and angular velocity. The internal density
structure can be either homogeneous or vary radially (e.g. in homogeneous
concentric spherical shells). The gravity potential of the sphere is not
constant on its surface because of the latitude-dependent centrifugal
potential.

**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.

Parameters
----------
name : str
A short name for the sphere, for example "Moon".
Definition: :math:R.
Units: :math:m.
geocentric_grav_const : float
The geocentric gravitational constant. The product of the mass of the
sphere :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 sphere.
Definition: :math:\omega.
Units: :math:\\rad.s^{-1}.
long_name : str or None
A long name for the sphere, for example "Moon Reference System"
(optional).
reference : str or None
Citation for the sphere parameter values (optional).

.. caution::

Must be used instead of :class:boule.Ellipsoid with zero flattening
for gravity calculations because it is impossible for a rotating sphere
to have constant gravity (gravitational + centrifugal) potential on its
surface. So the underlying ellipsoid gravity calculations don't apply
and are in fact singular when the flattening is zero.

Examples
--------

We can define a sphere by specifying the 3 key numerical parameters:

>>> sphere = Sphere(
...     name="Moon",
...     long_name="That's no moon",
...     geocentric_grav_const=2,
...     angular_velocity=0.5,
... )
>>> print(sphere) # doctest: +ELLIPSIS
Sphere(name='Moon', ...)
>>> print(sphere.long_name)
That's no moon

The sphere defines semi-axis, flattening, and some eccentricities similar
to :class:~bould.Ellipsoid for compatibility with the coordinate
conversion functions of pymap3d:

>>> print(sphere.semiminor_axis)
1
>>> print(sphere.semimajor_axis)
1
>>> print(sphere.first_eccentricity)
0
>>> print(sphere.eccentricity)
0
>>> print(sphere.flattening)
0
>>> print(sphere.thirdflattening)
0

"""

name = attr.ib()
geocentric_grav_const = attr.ib()
angular_velocity = attr.ib()
long_name = attr.ib(default=None)
reference = attr.ib(default=None)

"Check if the radius is positive."
if not value > 0:
raise ValueError(f"Invalid radius '{value}'. Should be greater than zero.")

@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 semiminor_axis(self):
"""
The semiminor axis of the sphere is equal to its radius. Added for
compatibility with pymap3d.
Definition: :math:b = R.
Units: :math:m.
"""

@property
def semimajor_axis(self):
"""
The semimajor axis of the sphere is equal to its radius. Added for
compatibility with pymap3d.
Definition: :math:a = R.
Units: :math:m.
"""

@property
def flattening(self):
r"""
The flattening of the sphere is equal to zero. Added for compatibility
with pymap3d.
Definition: :math:f = \dfrac{a - b}{a}.
"""
return 0

@property
def thirdflattening(self):
r"""
The third flattening of the sphere is equal to zero. Added for
compatibility with pymap3d
Definition: :math:f^{\prime\prime}= \dfrac{a -b}{a + b}.
"""
return 0

@property
def eccentricity(self):
"Alias for the first eccentricity."
return self.first_eccentricity

@property
def first_eccentricity(self):
r"""
The (first) eccentricity of the sphere is equal to zero. Added for
compatibility with pymap3d.
Definition: :math:e = \dfrac{\sqrt{a^2 - b^2}}{a} = \sqrt{2f - f^2}.
"""
return 0

[docs]    def normal_gravity(self, latitude, height, si_units=False):
r"""
Normal gravity of the sphere at the given latitude and height.

Computes the magnitude of the gradient of the gravity potential
(gravitational + centrifugal; see [HofmannWellenhofMoritz2006]_)
generated by the sphere at the given spherical latitude :math:\theta
and height above the surface of the sphere :math:h:

.. math::

\gamma(\theta, h) = \|\vec{\nabla}U(\theta, h)\|

in which :math:U = V + \Phi is the gravity potential of the sphere,
:math:V is the gravitational potential of the sphere, and
:math:\Phi is the centrifugal potential.

.. caution::

The current implementation is only valid for heights on or above
the surface of the sphere.

Parameters
----------
latitude : float or array
The spherical latitude where the normal gravity will be computed
(in degrees).
height : float or array
The height above the surface of the sphere of the computation point
(in meters).
si_units : bool
Return the value in mGal (False, default) or m/s² (True)

Returns
-------
gamma : float or array
The normal gravity in mGal or m/s².

Examples
--------

Normal gravity can be calculated at any spherical latitude and height
above the sphere:

>>> sphere = Sphere(
...     name="Moon",
...     long_name="That's no moon",
...     geocentric_grav_const=2,
...     angular_velocity=0.5,
... )
>>> gamma_equator = sphere.normal_gravity(latitude=0, height=0)
>>> print(f"{gamma_equator:.2f} mGal")
175000.00 mGal
>>> gamma_pole = sphere.normal_gravity(latitude=90, height=0)
>>> print(f"{gamma_pole:.2f} mGal")
200000.00 mGal

Notes
-----

The gradient of the gravity potential is the sum of the gravitational
:math:\vec{g} and centrifugal :math:\vec{f} accelerations for a
rotating sphere:

.. math::

\vec{\nabla}U(\theta, h) = \vec{g}(\theta, h)
+ \vec{f}(\theta, h)

The radial and latitudinal components of the two acceleration vectors
are:

.. math::

g_r = -\dfrac{GM}{(R + h)^2}

.. math::

g_\theta = 0

and

.. math::

f_r = \omega^2 (R + h) \cos^2 \theta

.. math::

f_\theta = \omega^2 (R + h) \cos\theta\sin\theta

in which :math:R is the sphere radius, :math:G is the gravitational
constant, :math:M is the mass of the sphere, and :math:\omega is
the angular velocity.

The norm of the combined gravitational and centrifugal accelerations
is:

.. math::

\gamma(\theta, h) = \sqrt{
\left( \dfrac{GM}{(R + h)^2} \right)^2
+ \left( \omega^2 (R + h) - 2\dfrac{GM}{(R + h)^2} \right)
\omega^2 (R + h) \cos^2 \theta
}

It's worth noting that a sphere under rotation is not in hydrostatic
equilibrium. Therefore unlike the oblate ellipsoid, it is not it's own
equipotential gravity surface (as is the case for the ellipsoid), the
gravity potential is not constant at the surface, and  the normal
gravity vector is not normal to the surface of the sphere.

"""
if np.any(height < 0):
warn(
"Formulas used are valid for points outside the sphere. "
"Height must be greater than or equal to zero."
)

gravity_acceleration = self.geocentric_grav_const / (radial_distance) ** 2
gamma = np.sqrt(
gravity_acceleration**2
+ (self.angular_velocity**2 * radial_distance - 2 * gravity_acceleration)
* self.angular_velocity**2
# Use cos^2 = (1 - sin^2) for more accurate results on the pole
* (1 - np.sin(np.radians(latitude)) ** 2)
)

# Convert gamma from SI to mGal
if not si_units:
gamma *= 1e5

return gamma

[docs]    def normal_gravitation(self, height, si_units=False):
r"""
Calculate normal gravitation at any height.

Computes the magnitude of the gradient of the gravitational potential
generated by the sphere at the given height :math:h:

.. math::

\gamma(h) = \|\vec{\nabla}V(h)\| = \dfrac{GM}{(R + h)^2}

in which :math:R is the sphere radius, :math:G is the gravitational
constant, and :math:M is the mass of the sphere.

.. caution::

These expressions are only valid for heights on or above the
surface of the sphere.

Parameters
----------
height : float or array
The height above the surface of the sphere of the computation point
(in meters).
si_units : bool
Return the value in mGal (False, default) or m/s² (True)

Returns
-------
gamma : float or array
The normal gravitation in mGal.

Examples
--------

Normal gravitation can be calculated at any point. However as this is a
sphere, only the height is used in the calculation.

>>> sphere = Sphere(
...     name="Moon",
...     long_name="That's no moon",
...     geocentric_grav_const=2,
...     angular_velocity=0.5,
... )
>>> g = sphere.normal_gravitation(height=1)
>>> print(f"{g:.2f} mGal")
50000.00 mGal

"""