Repository: fatiando/boule

Issue:
Fix normal gravity of Sphere
Is the centrifugal acceleration for the sphere we have implemented correct? Shouldn't it be projected to the normal of the sphere? This is where it gets a bit confusing. The "gravity" we want is the gradient of the combined potential. But how does that translate to the combined acceleration vector? Do we need to do a vector sum and then take the norm? I would imagine so but that doesn't seem to be what we do. Or am I missing something?

Discussion:
santisoler: > @santisoler maybe we should merge #51 first (with the failing test) and then use it here to make sure we're doing the right thing this time. What do you think?

I see the touching point between them and I don't think we should merge both features into a single PR, so if you agree, I'm ok with merging #51 even with failing CIs and then work on this PR.
santisoler: Bad news: Somigliana test fails :disappointed: 
santisoler: @MarkWieczorek I implemented some of the things we discussed on #50.
Although I agreed with your proposal, I forgot that we want `normal_gravity` to be in agreement with Somigliana.
Therefore I created three different methods for computing normal gravity of a Sphere:
- `normal_gravity`: computes the radial component of the gradient of the combined potential (gravitational + centrifugal)
- `normal_gravity_norm`: computes the norm of the gradient of the combined potential (gravitational + centrifugal)
- `normal_gravity_no_rotation`: computes the norm of the gradient of the gravitational potential (no rotating sphere).

I don't have any use case for the last two, so I'm not entirely sure if they are needed, but I would do like `normal_gravity` to agree with Somigliana.
Would you mind taking a look at it?
MarkWieczorek: Just some quick comments without looking at the code. As a reminder, this was what my suggestion was 

1. Define the normal_gravity() as GM/r2, which is the hydrostatic value for a non-rotating sphere. This will make it consistent with Somliagno's formula and the definition used in Ellipsoid.normal_gravity(). I.e., Ellipsoid->Sphere as f and angular_velocity go to zero.
2. Define a new function surface_acceleration() which is GM/r2 + centrifugal term. The docsring would need to explicitly state that the sphere is assumed to be perfectly rigid, and that it does not deform under rotation.

If you want normal_gravity to agree with Somliagno's formula (in the limit of zero rotation and flattening), then you have to define normal gravity as in (1) above. I don't think that there is anyway to avoid this. 

I don't think that there is much use for returning only the radial component as in your first definition above. If you really need this, I would instead probably return the full acceleration vector.
santisoler: Thanks @MarkWieczorek for the quick response.

> If you want normal_gravity to agree with Somliagno's formula (in the limit of zero rotation and flattening), then you have to define normal gravity as in (1) above. I don't think that there is anyway to avoid this.

I respectfully disagree. Somigliana equation does take into account the rotation of the ellipsoid. On #50 I referenced Somigliana equation and the expression I got after replacing the parameters with the ones of a sphere: the gravitational attraction minus the radial component of the centrifugal acceleration. So if we want `normal_gravity` to agree with Somigliana, we actually need it to return the radial component of the gradient of the combined potential.

> I don't think that there is much use for returning only the radial component as in your first definition above. If you really need this, I would instead probably return the full acceleration vector.

I have never worked with spherical ellipsoids, maybe that's why I'm facing this issues now. I don't know which is the common usage of such ellipsoids, and how normal gravity is usually computed. What I do want is our `Sphere` class to match the `Ellipsoid` class as much as possible, so I would like `normal_gravity` to return a single value (or one array) rather than a tuple with vector components, just to agree with the `Ellipsoid.normal_gravity`. 

Because other users might need other gravity computations, I though that would be clever to add them, like the norm of the gradient or the normal gravity without considering rotation.
But maybe we can just include `normal_gravity` on this PR and discuss later if we need some of the other methods.
MarkWieczorek: Here are some more comments on the way I see this. Its possible that we are using different definitions, so it probably best to start from the basics and to see at which point we disagree. 

## Conceptual framework
1. The "normal gravity" is defined as the magnitude of the gravitational acceleration on the surface of a rotating hydrostatic body.
2. For Earth, the normal gravity can be approximated using Somigliana's formula (among others), which depends on the semiminor axis, the semimajor axis, and the rotation rate (alternatively, instead of the rotation rate, you can use the gravity at the pole and equator).
3. A hydrostatic body that is not rotating corresponds to a sphere.
4. For a non-rotating hydrostatic body, the semimajor and semiminor axes are equal, and the normal gravity on the surface is everywhere GM/r**2. This is the limiting value of Somigliana's formula when the rotation rate goes to zero and the semimajor and semiminor axes approach each other.
5. A rigid sphere that is rotating is not in hydrostatic equilibrium.
6. For a rigid rotating sphere, the concept of "normal gravity" does not apply, because normal gravity assumes the object is in hydrotsatic equilibrium.
7. For a rotating rigid sphere, the magnitude of the gravity on the surface of the sphere can be easily computed from |GM/r^2 r-hat + g_centrigual|


## Proposal for Sphere class:
1. Normal gravity is defined as the magnitude of the gravitational acceleration on a non-rotating hydrostatic sphere. This is equal to GM/r^2.
2. The sphere class should provide the magnitude of the centrifugal acceleration. As on option [`vector=True`?], the three vector components could be returned.
3. The sphere class should provide the magnitude of the gravitational acceleration on the surface of the rigid rotating sphere. As on option, the three vector components could be returned. The name of this function should be descriptive, and one possible name is `surface_acceleration`.
santisoler: Thanks @MarkWieczorek for the knowledge and the advices.

I completely agree with the items of the conceptual framework.
And I like the first proposal for the Sphere class. As `Sphere` class is a subclass of `Ellipsoid`, I want it to have the `normal_gravity` method which should return a single value (or a single array).
But we can surely add methods for the centrifugal acceleration and the `surface_acceleration`, which can return either the three vector components or the magnitude, as we wish.

On this PR I will fix the `normal_gravity` method then. And maybe leave the extra methods for separated PRs.
I will ask you for a code review when I finish. I will take some of your words for writing the description of the class and the `normal_gravity` method, they are very good and clear!
leouieda: Hi Mark, thanks for the input! A few thoughts:

We seem to have different definitions of normal gravity here. We've been going with the one from physical geodesy "the magnitude of the gravity potential of a reference ellipsoid". This is also what is used in Geophysics for the Earth at least. There is no restriction on the computation point (e.g. surface) and it's gravity = gravitation + centrifugal. 

I'd like to keep things consistent and avoid having the ellipsoid return gravity and the sphere return gravitation. That is bound to be the source of bugs in the future. So this rules out the first part of proposal 1. 

> A hydrostatic body that is not rotating corresponds to a sphere.

Sure, but the moon is rotating and the reference surface for it is a sphere. So in practice this doesn't seem to apply.


> For a rigid rotating sphere, the concept of "normal gravity" does not apply, because normal gravity assumes the object is in hydrotsatic equilibrium.

Again, we seem to be using conflicting definitions here. 

I like the idea of providing other methods for calculating the magnitude and vectors of gravitation and centrifugal. But we also need a normal_gravity method that is consistent with the ellipsoid definition. Otherwise, we can't calculate gravity anomalies and disturbances from observed gravity (gravitational + centrifugal). 

Unless the data for the Moon etc don't include the centrifugal components. In that case, there would be issues with the ellipsoids is an angular velocity is provided.
MarkWieczorek: Could you let me know the first numbered point that you disagree with? If it is number 1, could you let me know which part you disagree with?

> We seem to have different definitions of normal gravity here. We've been going with the one from physical geodesy "the magnitude of the gravity potential of a reference ellipsoid". This is also what is used in Geophysics for the Earth at least. There is no restriction on the computation point (e.g. surface) and it's gravity = gravitation + centrifugal.

I don't think that there is any disagreement here. The only possible difference is that my definition in (1) requires the reference ellipsoid to be hydrostatic (i.e., to correpsond to a constant gravitational potential, as in Somigliana's derivation). 

santisoler: Hi @MarkWieczorek and @leouieda .

After reading more and doing some math I have new thoughts about this.

Firstly, I think we should stick with the gravity = gravitational + centrifugal definition. Therefore the _normal gravity_ of a sphere **must** include a centrifugal component.

Secondly, we can define the normal gravity as the norm of the gradient of the gravity potential (gravitational + centrifugal) of an equipotential ellipsoid of revolution. In the case of a rotating sphere is impossible to define normal gravity because a rotating sphere cannot be its own equipotential surface. 

So, in order to define a `normal_gravity` for the `Sphere` class we must sacrifice one of these hypotheses:
1. We ignore that the sphere is not its own equipotential surface, or
2. We return the gravity of a non-rotating sphere.

Following my idea that we **must** consider the centrifugal component, I would go with the first option: ignoring the fact that the sphere is not its own equipotential surface.

This opens another question.
When computing the normal gravity for an ellipsoid, the gradient has only one non-zero component: the normal component.
Therefore, the normal gravity is just the absolute value of the normal component.
On the sphere, the gradient has two non-zero components: the radial and the latitudinal one (due to a component of the centrifugal acceleration).
So, we have two options:
1. Return the norm of the gradient:
![image](https://user-images.githubusercontent.com/11541317/93121746-bed1b680-f69b-11ea-8394-e91440cc9631.png)
2. Return only the radial component
![image](https://user-images.githubusercontent.com/11541317/93122265-94ccc400-f69c-11ea-81b5-441ed7e9dc94.png)

If we choose the first option, the `normal_gravity` won't be in agreement with Somigliana.
But it will with the second one:
![image](https://user-images.githubusercontent.com/11541317/93122334-add57500-f69c-11ea-94d8-2fb8547f36e9.png)

## My opinion

I personally prefer the second one for two reasons:
1. It will make it compatible with Somigliana, which is the way I would use to compute the normal gravity of a spherical body if Boule wouldn't exists.
2. I think we can safely neglect the latitudinal component of the centrifugal acceleration. Spherical bodies has a very small angular velocity, therefore the centrifugal acceleration is very small in comparison with the gravitational acceleration. Therefore, not considering the latitudinal component of the small centrifugal acceleration won't introduce a considerable error.


I would like to hear your opinions on these @MarkWieczorek and @leouieda . Also would be nice to hear what @birocoles thinks of this. I'm open to suggestions, and remember that if any of you need a method that computes anything else (like the gravitational acceleration, for example) we can add it on a future PR as a separated method.
birocoles: Hi everyone, sorry for the delayed response. In my opinion, the hydrostatic equilibrium is not relevant for geophysicists. Actually, the geophysical normal Earth should be defined with the purposes of (1) keeping the difference from the actual gravity field as small as possible and (2) incorporate the internal layered structure of the Earth. It means that, for geophysical purposes, it is not necessary a normal Earth having a limiting surface defined by an equipotential of its own gravity field. There is a discussion about this topic in : Roman Pašteka, Ján Mikuška, and Bruno Meurers (2017). "Understanding the Bouguer Anomaly: A Gravimetry Puzzle". Elsevier. ISBN 978-0-12-812913-5. Based on these ideas, I think we should ignore that the sphere is not its own equipotential surface.
santisoler: Thanks @birocoles for your opinion!
I agree with your point of view: in geophysics we shouldn't care too much about the sphere not being in hydrostatic equilibrium.

I ended up reverting the commits so `normal_gravity` returns the norm of the gradient of the full gravity potential (gravitational + centrifugal). I also added an if statement on the function that tests ellipsoids against Somigliana, so the test is carried out only for non-zero flattening ellipsoids.

If you all agree, I think we could merge this. As I said before, if any of us wants to add any other method that returns some other computation, we can do it on a separate PR.
leouieda: I agree with @birocoles and @santisoler on the definitions here. That makes the normal gravity of ellipsoids and spheres have the same definition and avoids confusion when switching from one to the other (users shouldn't have to care which one they use). For reference, we might want to add a note about this to the Normal Gravity documentation page and reference the Bouguer anomaly book.

@MarkWieczorek we will make another issue to add a `Sphere.normal_gravitation` method so that there is a way to get the norm of the gravitation vector (without the centrifugal component). From what I understand, this is something that is useful for you to have, is that right? We can also add methods for the centrifugal part only, as you suggested. It would be nice to have these methods for the ellipsoid as well but that's for a later discussion.

@santisoler I made a few changes to the docstrings but nothing major. Will merge as soon as CIs are happy. 👍🏽 



Patch:
diff --git a/boule/sphere.py b/boule/sphere.py
index c099e62..166c15f 100644
--- a/boule/sphere.py
+++ b/boule/sphere.py
@@ -12,11 +12,11 @@ from . import Ellipsoid
 @attr.s(frozen=True)
 class Sphere(Ellipsoid):
     """
-    Reference spherical ellipsoid
+    Reference sphere (zero flattening ellipsoids)
 
-    Represents ellipsoids with zero flattening (spheres). Inherits methods and
-    properties of the :class:`boule.Ellipsoid`, guaranteeing no singularities
-    due to zero flattening (and thus zero eccentricity).
+    Inherits methods and properties of the :class:`boule.Ellipsoid` while
+    guaranteeing no singularities due to zero flattening (and thus zero
+    eccentricity) in gravity calculations.
 
     All parameters are in SI units.
 
@@ -56,6 +56,13 @@ class Sphere(Ellipsoid):
     1.00
     >>> print("{:.2f}".format(sphere.mean_radius))
     1.00
+    >>> print(sphere.normal_gravity(latitude=0, height=1))
+    25000.0
+    >>> print(sphere.normal_gravity(latitude=90, height=1))
+    25000.0
+
+    The flattening and eccentricities will all be zero:
+
     >>> print("{:.2f}".format(sphere.flattening))
     0.00
     >>> print("{:.2f}".format(sphere.linear_eccentricity))
@@ -64,10 +71,7 @@ class Sphere(Ellipsoid):
     0.00
     >>> print("{:.2f}".format(sphere.second_eccentricity))
     0.00
-    >>> print(sphere.normal_gravity(latitude=0, height=1))
-    25000.0
-    >>> print(sphere.normal_gravity(latitude=90, height=1))
-    25000.0
+
     """
 
     name = attr.ib()
@@ -88,19 +92,39 @@ class Sphere(Ellipsoid):
         return self.radius
 
     def normal_gravity(self, latitude, height):
-        """
+        r"""
         Calculate normal gravity at any latitude and height
 
         Computes the magnitude of the gradient of the gravity potential
         (gravitational + centrifugal) generated by the sphere at the given
-        latitude and height.
+        latitude :math:`\theta` and height :math:`h`:
+
+        .. math::
+
+            \gamma(\theta, h) =
+            \sqrt{\left( \frac{GM}{(R + h)^2} \right)^2
+            + \left(\omega^2 (R + h) - 2\frac{GM}{(R + h)^2} \right)
+            \omega^2 (R + h) \cos^2 \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.
+
+        .. note::
+
+            A sphere under rotation is not in hydrostatic equilibrium.
+            Therefore, it is not it's own equipotential gravity surface (as is
+            the case for the ellipsoid) and the normal gravity vector is not
+            normal to the surface of the sphere.
 
         Parameters
         ----------
         latitude : float or array
-            The latitude where the normal gravity will be computed (in degrees).
+            The latitude where the normal gravity will be computed (in
+            degrees). For a reference sphere there is no difference between
+            geodetic and spherical latitudes.
         height : float or array
-            The height (above the sphere) of computation point (in meters).
+            The height (above the sphere) of the computation point (in meters).
 
         Returns
         -------
@@ -110,24 +134,17 @@ class Sphere(Ellipsoid):
         References
         ----------
         [Heiskanen-Moritz]_
-        """
-        return self._gravity_sphere(height) + self._centrifugal_force(latitude, height)
 
-    def _gravity_sphere(self, height):
-        """
-        Calculate the gravity generated by a solid sphere (mGal)
-        """
-        return 1e5 * self.geocentric_grav_const / (self.radius + height) ** 2
-
-    def _centrifugal_force(self, latitude, height):
-        """
-        Calculate the centrifugal force due to the rotation of the sphere (mGal)
         """
-        return 1e5 * (
-            (-1)
+        radial_distance = self.radius + height
+        gravity_acceleration = self.geocentric_grav_const / (radial_distance) ** 2
+        return 1e5 * np.sqrt(
+            gravity_acceleration ** 2
+            + (self.angular_velocity ** 2 * radial_distance - 2 * gravity_acceleration)
             * self.angular_velocity ** 2
-            * (self.radius + height)
-            * np.cos(np.radians(latitude))
+            * radial_distance
+            # replace cos^2 with (1 - sin^2) for more accurate results on the pole
+            * (1 - np.sin(np.radians(latitude)) ** 2)
         )
 
     @property


Test patch:
diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py
index bc74dac..cf833f9 100644
--- a/boule/tests/test_ellipsoid.py
+++ b/boule/tests/test_ellipsoid.py
@@ -248,7 +248,10 @@ def test_normal_gravity_against_somigliana(ellipsoid):
     Check if normal gravity on the surface satisfies Somigliana equation
     """
     latitude = np.linspace(-90, 90, 181)
-    npt.assert_allclose(
-        ellipsoid.normal_gravity(latitude, height=0),
-        normal_gravity_surface(latitude, ellipsoid),
-    )
+    # Somigliana equation applies only to ellipsoids that are their own
+    # equipotential gravity surface. Spheres (with zero flattening) aren't.
+    if ellipsoid.flattening != 0:
+        npt.assert_allclose(
+            ellipsoid.normal_gravity(latitude, height=0),
+            normal_gravity_surface(latitude, ellipsoid),
+        )
diff --git a/boule/tests/test_sphere.py b/boule/tests/test_sphere.py
index cb7676e..17f0c1a 100644
--- a/boule/tests/test_sphere.py
+++ b/boule/tests/test_sphere.py
@@ -13,7 +13,10 @@ from .. import Sphere
 def sphere():
     "A spherical ellipsoid"
     ellipsoid = Sphere(
-        name="unit_sphere", radius=1.0, geocentric_grav_const=2.0, angular_velocity=1.3,
+        name="unit_sphere",
+        radius=1.0,
+        geocentric_grav_const=2.0,
+        angular_velocity=1.3,
     )
     return ellipsoid
 
@@ -141,22 +144,29 @@ def test_normal_gravity_only_rotation():
     """
     radius = 1
     omega = 2
+    heights = [1, 100, 1000]
     sphere = Sphere(
         name="sphere", radius=radius, geocentric_grav_const=0, angular_velocity=omega
     )
     # Check normal gravity on the equator
-    for height in [1, 2, 3, 4]:
-        expected_value = -1e5 * (omega ** 2) * (radius + height)
+    # Expected value is positive because normal gravity is the norm of the
+    # vector.
+    for height in heights:
+        expected_value = 1e5 * (omega ** 2) * (radius + height)
         npt.assert_allclose(
-            expected_value, sphere.normal_gravity(latitude=0, height=height),
+            expected_value,
+            sphere.normal_gravity(latitude=0, height=height),
         )
     # Check normal gravity on the poles (must be equal to zero)
-    for height in [1, 2, 3, 4]:
+    for height in heights:
         assert sphere.normal_gravity(latitude=90, height=height) < 1e-15
         assert sphere.normal_gravity(latitude=-90, height=height) < 1e-15
     # Check normal gravity at 45 degrees latitude
-    for height in [1, 2, 3, 4]:
-        expected_value = -1e5 * (omega ** 2) * (radius + height) * np.sqrt(2) / 2
+    # Expected value is positive because normal gravity is the norm of the
+    # vector.
+    for height in heights:
+        expected_value = 1e5 * (omega ** 2) * (radius + height) * np.sqrt(2) / 2
         npt.assert_allclose(
-            expected_value, sphere.normal_gravity(latitude=45, height=height),
+            expected_value,
+            sphere.normal_gravity(latitude=45, height=height),
         )