Fix normal gravity equation of Sphere

[santisoler]

Change how normal gravity of a sphere is computed.
It now computes it as the norm of the gradient of the combined potential (gravitational potential + centrifugal potential).

Fixes #50

Reminders:

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst and the base __init__.py file for the package.
• Write detailed docstrings for all functions/classes/methods. It often helps to design better code if you write the docstrings first.
• If adding new functionality, add an example to the docstring, gallery, and/or tutorials.
• Add your full name, affiliation, and ORCID (optional) to the AUTHORS.md file (if you haven't already) in case you'd like to be listed as an author on the Zenodo archive of the next release.

[leouieda]

@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?

[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?

Sure, I want to have the Somigliana test before merging this, so we can be sure that our normal gravity fix works fine.

[santisoler]

Bad news: Somigliana test fails 😞

[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 1 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.

Proposal 2 for Sphere class:

1. Do not include a method called normal_gravity as this will just confuse people.
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:
   
2. Return only the radial component
   

If we choose the first option, the normal_gravity won't be in agreement with Somigliana.
But it will with the second one:

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 think we should go with the first option and return the norm of the gradient because that's in agreement with the definition of normal gravity.
Besides, we can safely keep it not in agreement with Somigliana equation, which is only valid for ellipsoids with non-zero flattening. We are tempted to use it as a replacement because it doesn't have any singularity when eccentricity is zero, but it doesn't apply to spherical bodies.

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 @santisoler , @leouieda and @MarkWieczorek , sorry for my delayed response. There is a lot of discussion here and I would like to give you my humble opinion.

For me, the hydrostatic equilibrium is not relevant for geophysicists in this case. 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. 👍🏽