diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 6568d029..33ef076c 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -24,15 +24,26 @@ class Ellipsoid: The ellipsoid is defined by four parameters: semimajor axis, flattening, geocentric gravitational constant, and angular velocity. It spins around - its semiminor axis and has constant gravity potential at its surface. The + its semiminor axis and has a constant gravity potential at its surface. The internal density structure of the ellipsoid is unspecified but must be such - that the constant potential condition is satisfied. + that the constant gravity potential condition is satisfied. **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. + .. caution:: + + The :class:`boule.Ellipsoid` class with a flattening of zero is not + equivalent to the :class:`boule.Sphere` class. The internal density + structure of an Ellipsoid is unspecified, but must give rise to a + constant gravity potential on the Ellipsoid surface. For a Sphere, the + internal density depends only on radius, and when the angular velocity + is greater than zero, the gravity potential on the surface varies with + latitude. + + Parameters ---------- name : str @@ -62,13 +73,6 @@ class Ellipsoid: comments : str or None Additional comments regarding the ellipsoid (optional). - - .. caution:: - - Use :class:`boule.Sphere` if you desire zero flattening because there - are singularities for this particular case in the normal gravity - calculations. - Examples -------- @@ -147,17 +151,6 @@ def _check_flattening(self, flattening, value): f"Invalid flattening '{value}'. " "Should be greater than zero and lower than 1." ) - if value == 0: - raise ValueError( - "Flattening equal to zero will lead to errors in normal gravity. " - "Use boule.Sphere for representing ellipsoids with zero flattening." - ) - if value < 1e-7: - warn( - f"Flattening is too close to zero ('{value}'). " - "This may lead to inaccurate results and division by zero errors. " - "Use boule.Sphere for representing ellipsoids with zero flattening." - ) @semimajor_axis.validator def _check_semimajor_axis(self, semimajor_axis, value): @@ -302,12 +295,20 @@ def reference_normal_gravity_potential(self): + \dfrac{1}{3} \omega^2 a^2`. Units: :math:`m^2 / s^2`. """ - return ( - self.geocentric_grav_const - / self.linear_eccentricity - * np.arctan(self.linear_eccentricity / self.semiminor_axis) - + (1 / 3) * self.angular_velocity**2 * self.semimajor_axis**2 - ) + if self.flattening < 4.0e-16: + # Use the 3rd order arctan small-angle approximation to avoid + # numerical instabilities when the flattening is close to zero. + var = self.geocentric_grav_const * ( + 1 / self.semiminor_axis + - self.linear_eccentricity**2 / (3 * self.semiminor_axis**3) + ) + else: + var = ( + self.geocentric_grav_const + / self.linear_eccentricity + * np.arctan(self.linear_eccentricity / self.semiminor_axis) + ) + return var + (1 / 3) * (self.angular_velocity * self.semimajor_axis) ** 2 @property def area_equivalent_radius(self): @@ -362,16 +363,24 @@ def gravity_equator(self): centrifugal accelerations) at the equator on the surface of the ellipsoid. Units: :math:`m/s^2`. """ - ratio = self.semiminor_axis / self.linear_eccentricity - arctan = np.arctan2(self.linear_eccentricity, self.semiminor_axis) - aux = ( - self.second_eccentricity - * (3 * (1 + ratio**2) * (1 - ratio * arctan) - 1) - / (3 * ((1 + 3 * ratio**2) * arctan - 3 * ratio)) - ) + if self.flattening < 1.25e-5: + # Use the 5th order arctan small-angle approximation to avoid + # numerical instabilities when the flattening is close to zero. + aux = 3 + else: + ratio = self.semiminor_axis / self.linear_eccentricity + arctan = np.arctan2(self.linear_eccentricity, self.semiminor_axis) + aux = ( + self.second_eccentricity + * (3 * (1 + ratio**2) * (1 - ratio * arctan) - 1) + / (0.5 * ((1 + 3 * ratio**2) * arctan - 3 * ratio)) + ) + axis_mul = self.semimajor_axis * self.semiminor_axis result = ( - self.geocentric_grav_const * (1 - self._emm - self._emm * aux) / axis_mul + self.geocentric_grav_const + * (1 - self._emm - self._emm * aux / 6) + / axis_mul ) return result @@ -382,15 +391,22 @@ def gravity_pole(self): centrifugal accelerations) at the poles on the surface of the ellipsoid. Units: :math:`m/s^2`. """ - ratio = self.semiminor_axis / self.linear_eccentricity - arctan = np.arctan2(self.linear_eccentricity, self.semiminor_axis) - aux = ( - self.second_eccentricity - * (3 * (1 + ratio**2) * (1 - ratio * arctan) - 1) - / (1.5 * ((1 + 3 * ratio**2) * arctan - 3 * ratio)) - ) + if self.flattening < 1.25e-5: + # Use the 5th order arctan small-angle approximation to avoid + # numerical instabilities when the flattening is close to zero. + aux = 3 + else: + ratio = self.semiminor_axis / self.linear_eccentricity + arctan = np.arctan2(self.linear_eccentricity, self.semiminor_axis) + aux = ( + self.second_eccentricity + * (3 * (1 + ratio**2) * (1 - ratio * arctan) - 1) + / (0.5 * ((1 + 3 * ratio**2) * arctan - 3 * ratio)) + ) result = ( - self.geocentric_grav_const * (1 + self._emm * aux) / self.semimajor_axis**2 + self.geocentric_grav_const + * (1 + self._emm * aux / 3) + / self.semimajor_axis**2 ) return result @@ -676,13 +692,29 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): z_p2 = (self.semiminor_axis * sinbeta + height * sinlat) ** 2 # Distance squared between computation point and spin axis r_p2 = (self.semimajor_axis * cosbeta + height * coslat) ** 2 + # Auxialiary variables + z_pp2 = r_p2 - z_p2 + r_pp2 = r_p2 + z_p2 + + if self.flattening < 5.0e-5: + # Use the Taylor series approximation for flattenings close to + # zero to avoid numerical issues. + cosbeta_p2 = ( + 0.5 + + 0.5 * z_pp2 / r_pp2 + + big_e2 * 0.25 * (z_pp2**2 / r_pp2**3 - 1 / r_pp2) + ) + else: + # Auxiliary variables + big_d = z_pp2 / big_e2 + big_r = r_pp2 / big_e2 + # cos(reduced latitude) squared of the computation point + cosbeta_p2 = 0.5 + big_r / 2 - np.sqrt(0.25 + big_r**2 / 4 - big_d / 2) - # Auxiliary variables - big_d = (r_p2 - z_p2) / big_e2 - big_r = (r_p2 + z_p2) / big_e2 - - # cos(reduced latitude) squared of the computation point - cosbeta_p2 = 0.5 + big_r / 2 - np.sqrt(0.25 + big_r**2 / 4 - big_d / 2) + # Note that cosbeta_p2 can sometimes be less than 0 to within + # machine precision. To avoid taking the square root of a negative + # number, use the absolute value of this quantity. + cosbeta_p2 = np.abs(cosbeta_p2) # Semiminor axis of the ellipsoid passing through the computation # point. This is the coordinate u @@ -796,25 +828,28 @@ def normal_gravity(self, latitude, height, si_units=False): sinbeta2 = np.sin(np.radians(beta)) ** 2 cosbeta2 = 1 - sinbeta2 big_e = self.linear_eccentricity + big_e2 = big_e**2 - # Compute the auxiliary functions q and q_0 (eq 2-113 of - # HofmannWellenhofMoritz2006) - q_0 = 0.5 * ( - (1 + 3 * (self.semiminor_axis / big_e) ** 2) - * np.arctan2(big_e, self.semiminor_axis) - - 3 * self.semiminor_axis / big_e - ) - q_p = 3 * (1 + (u / big_e) ** 2) * (1 - u / big_e * np.arctan2(big_e, u)) - 1 - big_w = np.sqrt((u**2 + big_e**2 * sinbeta2) / (u**2 + big_e**2)) + # Compute auxiliary variables and the ratio E q_p / q_0 + if self.flattening <= 1.25e-5: + aux = 3 * self.semiminor_axis**3 / u**2 + else: + q_0 = 0.5 * ( + (1 + 3 * (self.semiminor_axis / big_e) ** 2) + * np.arctan2(big_e, self.semiminor_axis) + - 3 * self.semiminor_axis / big_e + ) + q_p = ( + 3 * (1 + (u / big_e) ** 2) * (1 - u / big_e * np.arctan2(big_e, u)) - 1 + ) + aux = big_e * q_p / q_0 + + big_w = np.sqrt((u**2 + big_e2 * sinbeta2) / (u**2 + big_e2)) # Put together gamma using 3 separate terms - term1 = self.geocentric_grav_const / (u**2 + big_e**2) + term1 = self.geocentric_grav_const / (u**2 + big_e2) term2 = (0.5 * sinbeta2 - 1 / 6) * ( - self.semimajor_axis**2 - * big_e - * q_p - * self.angular_velocity**2 - / ((u**2 + big_e**2) * q_0) + self.semimajor_axis**2 * aux * self.angular_velocity**2 / (u**2 + big_e2) ) term3 = -cosbeta2 * u * self.angular_velocity**2 gamma = (term1 + term2 + term3) / big_w @@ -844,10 +879,9 @@ def normal_gravitational_potential(self, latitude, height): gravitational constant, :math:`E` is the linear eccentricity, :math:`\omega` is the angular rotation rate, :math:`q` and :math:`q_0` are auxiliary functions, :math:`P_2` is the degree 2 unnormalized - Legendre Polynomial, and :math:`u` and :math:`\beta` are - ellipsoidal-harmonic coordinates corresponding to the input geodetic - latitude and ellipsoidal height. See eq. 2-124 of - [HofmannWellenhofMoritz2006]_. + Legendre Polynomial, and :math:`u` and :math:`beta` are ellipsoidal- + harmonic coordinates corresponding to the input geodetic latitude and + ellipsoidal height. See eq. 2-124 of [HofmannWellenhofMoritz2006]_. Assumes that the internal density distribution of the ellipsoid is such that the gravity potential is constant at its surface. @@ -879,24 +913,42 @@ def normal_gravitational_potential(self, latitude, height): "Height must be greater than or equal to zero." ) - # Convert geodetic latitude and height to ellipsoidal-harmonic coords + # Convert geodetic latitude and height to ellipsoidal-harmonic + # coordinates. Note that beta is in degrees. longitude, beta, u = self.geodetic_to_ellipsoidal_harmonic( None, latitude, height ) big_e = self.linear_eccentricity + big_e2 = big_e**2 - # Compute the auxiliary functions q and q_0 (eq 2-113 of + # Compute the term GM * arctan(E/u) / E + if self.flattening < 4.0e-16: + # Use the 3rd order arctan small-angle approximation to avoid + # numerical instabilities when the flattening is close to zero. + var = self.geocentric_grav_const * (1 / u - big_e2 / (3 * u**3)) + else: + var = self.geocentric_grav_const / big_e * np.arctan(big_e / u) + + # Compute the ratio of the auxiliary functions q and q_0 (eq 2-113 of # HofmannWellenhofMoritz2006) - q_0 = 0.5 * ( - (1 + 3 * (self.semiminor_axis / big_e) ** 2) - * np.arctan2(big_e, self.semiminor_axis) - - 3 * self.semiminor_axis / big_e - ) - q = 0.5 * ((1 + 3 * (u / big_e) ** 2) * np.arctan2(big_e, u) - 3 * u / big_e) + if self.flattening < 1.5e-5: + # Use the 5th order arctan small-angle approximation to avoid + # numerical instabilities when the flattening is close to zero. + ratio = (self.semiminor_axis / u) ** 3 + else: + q_0 = 0.5 * ( + (1 + 3 * (self.semiminor_axis / big_e) ** 2) + * np.arctan2(big_e, self.semiminor_axis) + - 3 * self.semiminor_axis / big_e + ) + q = 0.5 * ( + (1 + 3 * (u / big_e) ** 2) * np.arctan2(big_e, u) - 3 * u / big_e + ) + ratio = q / q_0 - big_v = self.geocentric_grav_const / big_e * np.arctan(big_e / u) + (1 / 3) * ( + big_v = var + (1 / 3) * ( self.angular_velocity * self.semimajor_axis - ) ** 2 * q / q_0 * (1.5 * np.sin(np.radians(beta)) ** 2 - 0.5) + ) ** 2 * ratio * (1.5 * np.sin(np.radians(beta)) ** 2 - 0.5) return big_v @@ -954,31 +1006,48 @@ def normal_gravity_potential(self, latitude, height): "Height must be greater than or equal to zero." ) - # Convert geodetic latitude and height to ellipsoidal-harmonic coords + # Convert geodetic latitude and height to ellipsoidal-harmonic + # coordinates. Note that beta is in degrees. longitude, beta, u = self.geodetic_to_ellipsoidal_harmonic( None, latitude, height ) big_e = self.linear_eccentricity + big_e2 = big_e**2 - # Compute the auxiliary functions q and q_0 (eq 2-113 of + # Compute the term GM * arctan(E/u) / E + if self.flattening < 4.0e-16: + # Use the 3rd order arctan small-angle approximation to avoid + # numerical instabilities when the flattening is close to zero. + var = self.geocentric_grav_const * (1 / u - big_e2 / (3 * u**3)) + else: + var = self.geocentric_grav_const / big_e * np.arctan(big_e / u) + + # Compute the ratio of the auxiliary functions q and q_0 (eq 2-113 of # HofmannWellenhofMoritz2006) - q_0 = 0.5 * ( - (1 + 3 * (self.semiminor_axis / big_e) ** 2) - * np.arctan2(big_e, self.semiminor_axis) - - 3 * self.semiminor_axis / big_e - ) - q = 0.5 * ((1 + 3 * (u / big_e) ** 2) * np.arctan2(big_e, u) - 3 * u / big_e) + if self.flattening < 1.5e-5: + # Use the 5th order arctan small-angle approximation to avoid + # numerical instabilities when the flattening is close to zero. + ratio = (self.semiminor_axis / u) ** 3 + else: + q_0 = 0.5 * ( + (1 + 3 * (self.semiminor_axis / big_e) ** 2) + * np.arctan2(big_e, self.semiminor_axis) + - 3 * self.semiminor_axis / big_e + ) + q = 0.5 * ( + (1 + 3 * (u / big_e) ** 2) * np.arctan2(big_e, u) - 3 * u / big_e + ) + ratio = q / q_0 big_u = ( - self.geocentric_grav_const / big_e * np.arctan(big_e / u) + var + 0.5 * (self.angular_velocity * self.semimajor_axis) ** 2 - * q - / q_0 + * ratio * (np.sin(np.radians(beta)) ** 2 - 1 / 3) + 0.5 * self.angular_velocity**2 - * (u**2 + big_e**2) + * (u**2 + big_e2) * np.cos(np.radians(beta)) ** 2 ) diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py index ec9be376..27ba909e 100644 --- a/boule/tests/test_ellipsoid.py +++ b/boule/tests/test_ellipsoid.py @@ -40,14 +40,6 @@ def test_check_flattening(): geocentric_grav_const=0, angular_velocity=0, ) - with pytest.raises(ValueError): - Ellipsoid( - name="zero_flattening", - semimajor_axis=1.0, - flattening=0, - geocentric_grav_const=0, - angular_velocity=0, - ) with pytest.raises(ValueError): Ellipsoid( name="almost_zero_negative_flattening", @@ -56,15 +48,6 @@ def test_check_flattening(): geocentric_grav_const=0, angular_velocity=0, ) - with warnings.catch_warnings(record=True) as warn: - Ellipsoid( - name="almost-zero-flattening", - semimajor_axis=1, - flattening=1e-8, - geocentric_grav_const=1, - angular_velocity=0, - ) - assert len(warn) >= 1 def test_check_semimajor_axis(): @@ -392,3 +375,29 @@ def test_normal_gravity_gravitational_centrifugal_potential(ellipsoid): def test_emm_wgs84(): "The _emm property should be equal this value for WGS84" npt.assert_allclose(WGS84._emm, 0.00344978650684) + + +@pytest.mark.parametrize( + "flattening", 10.0 ** (-np.arange(1, 20, 0.5)), ids=lambda f: f"{f:.1e}" +) +def test_no_warning_small_flattenings(flattening): + """ + Check that no warnings arise when using small values for the flattening + """ + latitude = np.linspace(-90, 90, 13) + height = np.array([[1.0e3] * len(latitude)]) + with warnings.catch_warnings(): + warnings.simplefilter("error") # raise error if warning is raised + ellipsoid = Ellipsoid( + name="WGS84 with small flattening", + semimajor_axis=WGS84.semimajor_axis, + flattening=flattening, + geocentric_grav_const=WGS84.geocentric_grav_const, + angular_velocity=WGS84.angular_velocity, + ) + _ = ellipsoid.reference_normal_gravity_potential + _ = ellipsoid.gravity_pole + _ = ellipsoid.gravity_equator + _ = ellipsoid.normal_gravitational_potential(latitude, height) + _ = ellipsoid.normal_gravity_potential(latitude, height) + _ = ellipsoid.normal_gravity(latitude, height) diff --git a/doc/user_guide/normal_gravity.rst b/doc/user_guide/normal_gravity.rst index 25559bb2..ddcac00a 100644 --- a/doc/user_guide/normal_gravity.rst +++ b/doc/user_guide/normal_gravity.rst @@ -91,7 +91,8 @@ And plotted with :mod:`pygmt`: particularly useful for geophysics. These calculations can be performed for any oblate ellipsoid (see -:ref:`ellipsoids`). Here is the normal gravity of the Martian ellipsoid: +:ref:`ellipsoids`). Here is an example for the normal gravity of the Martian +ellipsoid: .. jupyter-execute:: @@ -109,10 +110,18 @@ These calculations can be performed for any oblate ellipsoid (see Notice that the overall trend is the same as for the Earth (the Martian -ellipsoid is also oblate) but the range of values is different. The mean +ellipsoid is slightly more oblate than Earth) but the range of values +is different. The mean gravity on Mars is much weaker than on the Earth: around 370,000 mGal or 3.7 m/s² when compared to 970,000 mGal or 9.7 m/s² for the Earth. +The computations of the gravimetric quantities in boule are accurate for oblate +ellipsoids with flattenings that are arbitrarily small. In fact, even a +flattening of zero is permissible. Whereas the standard textbook equations +become numerically unstable when the flattening is less than about +:math:`10-^{-7}`, boule makes use of approximate equations in the low flattening +limit that do not suffer any numerical limitations. + .. admonition:: Assumptions for oblate ellipsoids :class: important @@ -132,16 +141,28 @@ Spheres ------- Method :meth:`boule.Sphere.normal_gravity` performs the normal gravity -calculations for spheres. It behaves mostly the same as the oblate ellipsoid -version except that the latitude is a *geocentric spherical latitude* instead -of a geodetic latitude (for spheres they are actually the same thing). +calculations for spheres. This method behaves mostly the same as the oblate +ellipsoid version, with two exceptions. First, spherical coordinates are +used in the case of a sphere, and the latitude coordinate corresponds to +*geocentric spherical latitude*. Geodetic and spherical latitude are, in fact, +the same for an ellipsoid with zero flattening. + +Second, boule makes the assumption that the interior density distribution of +the planet varies only as a function of radius. Because of this, the normal +gravitation potential is constant on the sphere surface, but the normal gravity +potential (which includes the centrifugal potential) is not. + +One planetary object that makes use of the Sphere class is the Moon. This +example computes the normal gravity of the Moon at 45 degrees latitude +and for heights between 0 and 1 km above the reference radius. .. jupyter-execute:: gamma = bl.Moon2015.normal_gravity(latitude=45, height=height) print(gamma) -This is what the normal gravity of Moon looks like on a map: +This is what the normal gravity of Moon looks like in map form, 10 km above +the surface: .. jupyter-execute:: @@ -162,22 +183,26 @@ This is what the normal gravity of Moon looks like on a map: Normal gravity of spheres is calculated under the following assumptions: + * The :term:`gravitational potential` is constant on the surface of the + ellipsoid. + * The internal density structure is unspecified but must be either + homogeneous or vary only as a function of radius (e.g., in concentric + layers). * The normal gravity is the magnitude of the gradient of the :term:`gravity potential` of the sphere. - * The internal density structure is unspecified but must be either - homogeneous or vary radially (e.g., in concentric layers). - A constant gravity potential on the surface of a rotating sphere is not - possible. Therefore, the normal gravity calculated for a sphere is - different than that of an oblate ellipsoid (hence why we need a separate - method of calculation). + **Important:** Unlike an ellipsoid, the normal gravity potential of a + sphere is not constant on its surface, and the normal gravity vector is + not perpendicular to the surface. + Gravity versus gravitation ++++++++++++++++++++++++++ -Notice that the variation between poles and equator is much smaller than for -the Earth or Mars. -That's because the **variation is due solely to the centrifugal acceleration**. +Notice that the variation of the normal gravity between the poles and equator +for the Moon is much smaller than for the Earth or Mars. +That's because the **variation is due solely to the centrifugal acceleration**, +and the angular rotation rate of the Moon is small. We can see this clearly when we calculate the **normal gravitation** (without the centrifugal component) using :meth:`boule.Sphere.normal_gravitation`: