Fix normal gravity of Sphere about boule HOT 8 CLOSED

As we discussed today, we haven't implemented it correctly. The normal gravity by definition is the norm of the gradient of the potential, which includes both the gravitational potential and the centrifugal one. Therefore we cannot compute it as the sum of the norms of the acceleration vectors, like we are currently doing.

So, the gravitational acceleration vector of a sphere is:

While, the centrifugal acceleration vector is:

which always lives in the equatorial plane and has only radial and latitudinal components.

Adding these two vectors and computing its norm, we obtain the normal gravity of a sphere:

Therefore, the question turned out to be a bug and we need to fix it.

from boule.

Hurray for basic maths 🙂 🎉

from boule.

After seeing that the Somigliana tests fails when comparing against the normal gravity equation on the previous comment.
That triggered me some reading on the level ellipsoids and the definition of the normal gravity.

So, for a reference ellipsoid (flattening != 0) is defined as

[...] an ellipsoid of revolution which is an equipotential surface of a normal gravity field.
Moritz (2005). Physical Geodesy

Where the potential is defined as the sum of the gravitational potential and the centrifugal potential.

Therefore the normal gravity potential is constant on the surface of the ellipsoid.
That's why, on the calculation of the normal gravity, the author obtains it as the module of the "vertical" component of the acceleration vector (the latitudinal and longitudinal components are zero because we are on the equipotential surface).

So, for a reference ellipsoid (with flattening != 0):

• The ellipsoid is its own equipotential surface
• Normal gravity, defined as the norm of the acceleration vector, is equal to the absolute value of the "vertical" (non-latitudinal and not-longitudinal) component.

This raise me a question: if the ellipsoid is its own equipotential surface, wouldn't the flattening and the angular_velocity be dependent parameters? I mean, for a higher angular_velocity, the equipotential surface will be more oblate, and thus the ellipsoid flattening should also increase.

Leaving this question aside, these definitions introduce some problems into the definition of a spherical ellipsoid (and perhaps these are the sources of the singularities).

If the spherical ellipsoids rotates (angular_velocity != 0), its equipotential surface cannot be a sphere. This is really problematic, because with the ellipsoid = equipotential surface definition, we cannot define a spherical ellipsoid that rotates.

This problem raises a difference on the definition of the normal gravity.
Following the definition of the normal gravity as the module of the gradient of the combined potential (gravitational + centrifugal), we obtain the normal gravity expression that I wrote on the previous comment. But, obviously, this expression doesn't agree with the Somigliana equation (normal gravity at the surface of the ellipsoid), which can be expressed as:

For a sphere, the semimajor axis and the semiminor axis are equal, therefore:

Which is equal to the radial component of the gravity acceleration on the surface of the sphere.

So, what do you think @leouieda @MarkWieczorek ?
Do you agree with this reasoning? Can you find any mistake? Do you know something I might be missing?

from boule.

It is debatable whether it makes sense to define a "normal gravity" for a rotating rigid sphere in the same way that geodecists do for the Earth. An object that is rotating and that remains a sphere can not be in hydrostatic equilibirum.

I think that there are two ways to approach this problem:

1. We can simply define the normal gravity on a rotating rigid sphere as the norm of the GM/r2 + centrifugual term, or
2. We can call this something else, in order not to confuse people.

I am ok with option 1. The only use cases for are for Mercury, Venus, and the Moon. For these bodies the rotation rate is so small and the flattening is so small that people just assume that the reference ellipsoid is a sphere. When you calculate the total gravitational acceleration from a spacecraft model, you simply want to subtract out the constant GM/r2 term. The centrifugal contribution for these cases is in the noise, and it doesn't matter if you take this into account or not.

So: in the docstring/documentation, I think that we just need to add a Note at the bottom stating the normal gravity in this case is defined as the norm of GM/r2 + the centrifugual term, and that a constant potential surface does not correspond to a sphere.

from boule.

Thanks @MarkWieczorek for the swift response.

It is debatable whether it makes sense to define a "normal gravity" for a rotating rigid sphere in the same way that geodecists do for the Earth. An object that is rotating and that remains a sphere can not be in hydrostatic equilibirum.

So, with the current definition of a reference ellipsoid, shouldn't we remove the angular_velocity of the Sphere?
I mean, if it's rotating (and the rotation is not negligible), then you shouldn't use a spherical ellipsoid, right?
If we do remove the angular_velocity from the Sphere, then all the normal gravity problems are solved: we won't have centrifugal acceleration and therefore the norm of the acceleration will be equal to the radial component.

Another option could be what @MarkWieczorek proposes: leave the angular_velocity and decide whether we want to keep the normal gravity definition as the norm of the gradient of the combined potential or we want to define it as the radial component of the acceleration.
If we choose the former, spherical ellipsoids won't satisfy Somigliana equation, but with the latter, the Sphere.normal_gravity will be in agreement with the Ellipsoid.normal_gravity: the Sphere.normal_gravity will represent the same as Ellipsoid.normal_gravity in the limit of flattening = 0. We can also define a new method that computes the norm of the acceleration vector, just in case a user might need it instead.

In summary:

1. Sphere has not rotation: it will be its own equipotential surface and the normal gravity will be in agreement with the Ellipsoid one.
2. Sphere has rotation and:
   a. Normal gravity is defined as the norm of the combined acceleration (not in agreement with Somigliana)
   b. Normal gravity is defined as the radial component (in agreement with Somigliana, and it will be equivalent of the Ellipsoid.normal_gravity in the limit of flattening=0). Optional: add a new method for the norm of the acceleration vector.

I think I feel more comfortable with choosing the first option, but don't know if it satisfies the needs of the potential users.

from boule.

I think that we need to include the rotation rate in this class, even if it is never used. The logic is that there is no such thing as a non-rotating sphere in planetary science. We use spherical reference frames only when it is impossible to compute the flattening from observations. As an example: The natural variability in the degree-2 shape of Venus from geological processes is larger than the degree-2 shape you would expect based on its small rotation rate. Nevertheless, we often need to know what the rotation rate is for other calculations, so it would be useful to have access to it.

My revised proposal would be the following:

1. Include angular_velocity as a parameter,
2. Define the normal_gravity() as GM/r2, which is the hydrostatic value for a non-rotating sphere. This will make it consistent with the definition used in Ellipsoid.normal_gravity(): Ellipsoid->Sphere as f and angular_velocity go to zero.
3. 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.

How does that sound?

from boule.

I really like this approach 🚀
I think it's the best compromise between a robust API with solid foundation that also satisfies any other possible use cases.

What do you think @leouieda ?

from boule.

For the record, after following how Moritz (2005) arrives to the potential of a reference ellipsoid, now I can answer the following question:

if the ellipsoid is its own equipotential surface, wouldn't the flattening and the angular_velocity be dependent parameters? I mean, for a higher angular_velocity, the equipotential surface will be more oblate, and thus the ellipsoid flattening should also increase.

If we evaluate Equation (2-126) on the surface of the ellipsoid (u = b), the gravity potential doesn't depend on the latitude, therefore is constant for every point in the surface of the ellipsoid. So the flattening and the angular_velocity can be independent of each other.

from boule.