From 6c716ac984749659254a36b02c9aac1b72d46685 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Mon, 8 Apr 2024 14:29:49 +0200 Subject: [PATCH 01/29] Add centrifugal potential to Ellipsoid and Sphere --- boule/_ellipsoid.py | 44 ++++++++++++++++++++++++++++++++++++++++++++ boule/_sphere.py | 36 ++++++++++++++++++++++++++++++++++++ 2 files changed, 80 insertions(+) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index c2f9a820..faaf1f1f 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -619,3 +619,47 @@ def normal_gravity(self, latitude, height, si_units=False): gamma *= 1e5 return gamma + + def centrifugal_potential(self, latitude, height): + r""" + Centrifugal potential at the given geodetic latitude and height above + the ellipsoid. + + Parameters + ---------- + latitude : float or array + The geodetic latitude where the centrifugal potential will be + computed (in degrees). + height : float or array + The ellipsoidal (geometric) height of the computation point (in + meters). + + Returns + ------- + Phi : float or array + The centrifugal potential in m²/s². + + Notes + ----- + + The centrifugal potential :math:`\Phi` at geodetic latitude + :math:`\phi` and height above the ellipsoid :math:`h` (geometric + height) is + + .. math:: + + \Phi(\phi, h) = \dfrac{1}{2} + \omega^2 \left(N(\phi) + h\right)^2 \cos^2(\phi) + + in which :math:`N(\phi)` is the prime vertical radius of curvature of + the ellipsoid. + """ + # Pre-compute to avoid repeated calculations + sinlat = np.sin(np.radians(latitude)) + coslat = np.sqrt(1 - sinlat**2) + + return (1 / 2) * ( + self.angular_velocity + * (self.prime_vertical_radius(sinlat) + height) + * coslat + ) ** 2 diff --git a/boule/_sphere.py b/boule/_sphere.py index 462a1b04..30c28cbb 100644 --- a/boule/_sphere.py +++ b/boule/_sphere.py @@ -379,3 +379,39 @@ def normal_gravitation(self, height, si_units=False): gamma *= 1e5 return gamma + + def centrifugal_potential(self, latitude, height): + r""" + Centrifugal potential at the given latitude and height. + + Parameters + ---------- + latitude : float or array + The latitude where the centrifugal potential will be computed + (in degrees). + height : float or array + The height above the sphere of the computation point (in meters). + + Returns + ------- + Phi : float or array + The centrifugal potential in m²/s². + + Notes + ----- + + The centrifugal potential :math:`\Phi` at latitude :math:`\phi` and + height above the sphere :math:`h` is + + .. math:: + + \Phi(\phi, h) = \dfrac{1}{2} + \omega^2 \left(R + h\right)^2 \cos^2(\phi) + + in which :math:`R` is the sphere radius. + """ + return (1 / 2) * ( + self.angular_velocity + * (self.radius + height) + * np.cos(np.radians(latitude)) + ) ** 2 From 6409055b8ee5405d867a1e6deb823db2fec917a8 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Wed, 10 Apr 2024 14:11:16 +0200 Subject: [PATCH 02/29] Add geodetic_to_ellipsoidal_harmonic to Ellipsoid class --- boule/_ellipsoid.py | 77 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 75 insertions(+), 2 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index c2f9a820..a6d77adc 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -435,8 +435,8 @@ def geodetic_to_spherical(self, longitude, latitude, height): ------- longitude : array Longitude coordinates on geocentric spherical coordinate system in - degrees. - The longitude coordinates are not modified during this conversion. + degrees. The longitude coordinates are not modified during this + conversion. spherical_latitude : array Converted latitude coordinates on geocentric spherical coordinate system in degrees. @@ -501,6 +501,79 @@ def spherical_to_geodetic(self, longitude, spherical_latitude, radius): height = (k + self.eccentricity**2 - 1) / k * hypot_dz return longitude, latitude, height + def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): + """ + Convert from geodetic to ellipsoidal-harmonic coordinates. + + The geodetic datum is defined by this ellipsoid, and the coordinates + are converted following [Lakshmanan1991]_ and [LiGotze2001]. + + Parameters + ---------- + longitude : array + Longitude coordinates on geodetic coordinate system in degrees. + latitude : array + Latitude coordinates on geodetic coordinate system in degrees. + height : array + Ellipsoidal heights in meters. + + Returns + ------- + longitude : array + Longitude coordinates on ellipsoidal-harmonic coordinate system in + degrees. The longitude coordinates are not modified during this + conversion. + reduced_latitude : array + The reduced (or parametric) latitude in degrees. + u : array + The coordinate u, which is the semiminor axis of the ellipsoid that + passes through the input coordinates. + """ + if height == 0: + # Use simple formula that relates geodetic and reduced latitude + beta = np.degrees( + np.arctan( + self.semiminor_axis + / self.semimajor_axis + * np.tan(np.radians(latitude)) + ) + ) + u = self.semiminor_axis + + else: + # The variable names follow Li and Goetze (2001). The prime terms + # (*_p) refer to quantities on an ellipsoid passing through the + # computation point. + sinlat = np.sin(np.radians(latitude)) + coslat = np.sqrt(1 - sinlat**2) + + # Reduced latitude of the projection of the point on the + # reference ellipsoid + beta = np.arctan2( + self.semiminor_axis * sinlat, self.semimajor_axis * coslat + ) + sinbeta = np.sin(beta) + cosbeta = np.sqrt(1 - sinbeta**2) + + # Distance squared between computation point and equatorial plane + 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 + + # Auxiliary variables + big_d = (r_p2 - z_p2) / self.linear_eccentricity**2 + big_r = (r_p2 + z_p2) / self.linear_eccentricity**2 + + # 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) + + # Semiminor axis of ellipsoid passing through the computation point + u = np.sqrt(r_p2 + z_p2 - self.linear_eccentricity**2 * cosbeta_p2) + + beta = np.degrees(np.arccos(np.sqrt(cosbeta_p2))) + + return longitude, beta, u + def normal_gravity(self, latitude, height, si_units=False): r""" Normal gravity of the ellipsoid at the given latitude and height. From 13d6a6fc21eeb6e903a7684aa90a8680d0439760 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Wed, 10 Apr 2024 15:49:43 +0200 Subject: [PATCH 03/29] Add inverse --- boule/_ellipsoid.py | 50 ++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 47 insertions(+), 3 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index a6d77adc..9e966b78 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -540,6 +540,8 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): ) u = self.semiminor_axis + return longitude, beta, u + else: # The variable names follow Li and Goetze (2001). The prime terms # (*_p) refer to quantities on an ellipsoid passing through the @@ -568,11 +570,53 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): cosbeta_p2 = 0.5 + big_r / 2 - np.sqrt(0.25 + big_r**2 / 4 - big_d / 2) # Semiminor axis of ellipsoid passing through the computation point - u = np.sqrt(r_p2 + z_p2 - self.linear_eccentricity**2 * cosbeta_p2) + b_p = np.sqrt(r_p2 + z_p2 - self.linear_eccentricity**2 * cosbeta_p2) + + beta_p = np.degrees(np.arccos(np.sqrt(cosbeta_p2))) + + return longitude, beta_p, b_p + + def ellipsoidal_harmonic_to_geodetic(self, longitude, reduced_latitude, u): + """ + Convert from ellipsoidal-harmonic coordinates to geodetic coordinates. + + The geodetic datum is defined by this ellipsoid. + + Parameters + ---------- + longitude : array + Longitude coordinates on ellipsoidal-harmonic coordinate system in + degrees. + latitude : array + Reduced (parametric) latitude coordinates on ellipsoidal-harmonic + coordinate system in degrees. + u : array + The coordinate u, which is the semiminor axes of the ellipsoid that + passes through the input coordinates. - beta = np.degrees(np.arccos(np.sqrt(cosbeta_p2))) + Returns + ------- + longitude : array + Longitude coordinates on geodetic coordinate system in degrees. The + longitude coordinates are not modified during this conversion. + latitude : array + Latitude coordinates on geodetic coordinate system in degrees. + height : array + Ellipsoidal heights in meters. + """ + # semimajor axis of the ellipsoid that passes through the input + # coordinates + a_p = np.sqrt(u**2 + self.linear_eccentricity**2) + # latitude = np.arctan(a_p / u * np.tan(np.radians(reduced_latitude))) + latitude = np.arctan(a_p / u * np.tan(np.radians(reduced_latitude))) + + # Compute height as the difference of the prime_vertical_radius of the + # input ellipsoid and reference ellipsoid + height = self.prime_vertical_radius(np.sin(latitude)) * ( + a_p / self.semimajor_axis - 1 + ) - return longitude, beta, u + return longitude, np.degrees(latitude), height def normal_gravity(self, latitude, height, si_units=False): r""" From 113bcc0d782da947a70ea1b059fe6638c6d0ebc5 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Wed, 10 Apr 2024 16:22:34 +0200 Subject: [PATCH 04/29] minor cleanup --- boule/_ellipsoid.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 9e966b78..aea2c28d 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -607,7 +607,8 @@ def ellipsoidal_harmonic_to_geodetic(self, longitude, reduced_latitude, u): # semimajor axis of the ellipsoid that passes through the input # coordinates a_p = np.sqrt(u**2 + self.linear_eccentricity**2) - # latitude = np.arctan(a_p / u * np.tan(np.radians(reduced_latitude))) + + # geodetic latitude latitude = np.arctan(a_p / u * np.tan(np.radians(reduced_latitude))) # Compute height as the difference of the prime_vertical_radius of the From 78bd56bf7cfe00b87d54fda8c72e4625192f5e69 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Wed, 10 Apr 2024 22:34:05 +0200 Subject: [PATCH 05/29] improve checking of height==0 case --- boule/_ellipsoid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index aea2c28d..2bee5d75 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -529,7 +529,7 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): The coordinate u, which is the semiminor axis of the ellipsoid that passes through the input coordinates. """ - if height == 0: + if np.array(height).all() == 0: # Use simple formula that relates geodetic and reduced latitude beta = np.degrees( np.arctan( From d2fe1cbf209fdf656f17ba22a5fb531b459a5c26 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 11 Apr 2024 14:22:23 +0200 Subject: [PATCH 06/29] Add normal_gravitational_potential --- boule/_ellipsoid.py | 79 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 2bee5d75..326e2adb 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -737,3 +737,82 @@ def normal_gravity(self, latitude, height, si_units=False): gamma *= 1e5 return gamma + + def normal_gravitational_potential(self, latitude, height): + r""" + Normal gravitational potential of the ellipsoid at the given latitude + and height. + + Computes the gravitational potential generated by the ellipsoid at the + given geodetic latitude :math:`\phi` and height above the ellipsoid + :math:`h` (geometric height). + + .. math:: + + \V(\beta, u) = \dfrac{GM}{E} \arctan{\dfrac{E}{u}} + \dfrac{1}{3} + \omega^2 a^2 \dfrac{q}{q_0} P_2 (\sin \beta) + + in which :math:`V` is the gravitational potential of the + ellipsoid (no centrifugal term), :math:`GM` is the geocentric + gravitational constant, :math:`E` is the linear eccentricty, + :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 + 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. + + .. caution:: + + These expressions are only valid for heights on or above the + surface of the ellipsoid. + + Parameters + ---------- + latitude : float or array + The geodetic latitude where the normal gravitational potential will + be computed (in degrees). + height : float or array + The ellipsoidal (geometric) height of the computation the point (in + meters). + + Returns + ------- + V : float or array + The normal gravitational potential in m²/s². + + """ + # Warn if height is negative + if np.any(height < 0): + warn( + "Formulas used are valid for points outside the ellipsoid." + "Height must be greater than or equal to zero." + ) + + # Convert geodetic latitude and height to ellipsoidal-harmonic coords + longitude, beta, u = self.geodetic_to_ellipsoidal_harmonic( + None, latitude, height + ) + + # compute the auxiliary functions q and q_0 (eq 2-113 of + # HofmannWellenhofMoritz2006) + q_0 = 0.5 * ( + (1 + 3 * (self.semiminor_axis / self.linear_eccentricity) ** 2) + * np.arctan2(self.linear_eccentricity, self.semiminor_axis) + - 3 * self.semiminor_axis / self.linear_eccentricity + ) + q = 0.5 * ( + (1 + 3 * (u / self.linear_eccentricity) ** 2) + * np.arctan2(self.linear_eccentricity, u) + - 3 * u / self.linear_eccentricity + ) + + big_V = self.geocentric_grav_const / self.linear_eccentricity * np.arctan( + self.linear_eccentricity / u + ) + (1 / 3) * (self.angular_velocity * self.semimajor_axis) ** 2 * q / q_0 * ( + 1.5 * np.sin(np.radians(beta)) ** 2 - 0.5 + ) + + return big_V From 58bf1079398a79a30c99d8932cf257246b3f375c Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 11 Apr 2024 14:43:51 +0200 Subject: [PATCH 07/29] Add normal_gravity_potential --- boule/_ellipsoid.py | 89 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 89 insertions(+) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 326e2adb..08bf073e 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -816,3 +816,92 @@ def normal_gravitational_potential(self, latitude, height): ) return big_V + + def normal_gravity_potential(self, latitude, height): + r""" + Normal gravity potential of the ellipsoid at the given latitude and + height. + + Computes the gravity potential generated by the ellipsoid at the + given geodetic latitude :math:`\phi` and height above the ellipsoid + :math:`h` (geometric height). + + .. math:: + + \U(\beta, u) = \dfrac{GM}{E} \arctan{\dfrac{E}{u}} + \dfrac{1}{2} + \omega^2 a^2 \dfrac{q}{q_0} \left(\sin^2 \beta + - \dfrac{1}{3}\right) + \dfrac{1}{2} \omega^2 \left(u^2 + E^2) + \cos^2 \beta + + in which :math:`U` is the gravity potential of the ellipsoid, + :math:`GM` is the geocentric gravitational constant, :math:`E` is the + linear eccentricty, :math:`\omega` is the angular rotation rate, + :math:`q` and :math:`q_0` are auxiliary functions, and :math:`u` and + :math:`beta` are ellipsoidal-harmonic coordinates corresponding to the + input geodetic latitude and and ellipsoidal height. See eq. 2-126 of + [HofmannWellenhofMoritz2006]_. + + Assumes that the internal density distribution of the ellipsoid is such + that the gravity potential is constant at its surface. + + .. caution:: + + These expressions are only valid for heights on or above the + surface of the ellipsoid. + + Parameters + ---------- + latitude : float or array + The geodetic latitude where the normal gravity potential will be + computed (in degrees). + height : float or array + The ellipsoidal (geometric) height of the computation the point (in + meters). + + Returns + ------- + V : float or array + The normal gravitational potential in m²/s². + + """ + # Warn if height is negative + if np.any(height < 0): + warn( + "Formulas used are valid for points outside the ellipsoid." + "Height must be greater than or equal to zero." + ) + + # Convert geodetic latitude and height to ellipsoidal-harmonic coords + longitude, beta, u = self.geodetic_to_ellipsoidal_harmonic( + None, latitude, height + ) + + # compute the auxiliary functions q and q_0 (eq 2-113 of + # HofmannWellenhofMoritz2006) + q_0 = 0.5 * ( + (1 + 3 * (self.semiminor_axis / self.linear_eccentricity) ** 2) + * np.arctan2(self.linear_eccentricity, self.semiminor_axis) + - 3 * self.semiminor_axis / self.linear_eccentricity + ) + q = 0.5 * ( + (1 + 3 * (u / self.linear_eccentricity) ** 2) + * np.arctan2(self.linear_eccentricity, u) + - 3 * u / self.linear_eccentricity + ) + + big_U = ( + self.geocentric_grav_const + / self.linear_eccentricity + * np.arctan(self.linear_eccentricity / u) + + 0.5 + * (self.angular_velocity * self.semimajor_axis) ** 2 + * q + / q_0 + * (np.sin(np.radians(beta)) ** 2 - 1 / 3) + + 0.5 + * self.angular_velocity**2 + * (u**2 + self.linear_eccentricity**2) + * np.cos(np.radians(beta)) ** 2 + ) + + return big_U From 5ac3321b63abafe451f4528179166e2b3e065c48 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 11 Apr 2024 15:39:14 +0200 Subject: [PATCH 08/29] Add tests to check gravity, gravitationa, centrifugal potentials --- boule/_ellipsoid.py | 72 +++++++++++++++++------------------ boule/tests/test_ellipsoid.py | 18 +++++++++ 2 files changed, 54 insertions(+), 36 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index ae6588f4..d5c76940 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -906,46 +906,46 @@ def normal_gravity_potential(self, latitude, height): return big_U -def centrifugal_potential(self, latitude, height): - r""" - Centrifugal potential at the given geodetic latitude and height above - the ellipsoid. + def centrifugal_potential(self, latitude, height): + r""" + Centrifugal potential at the given geodetic latitude and height above + the ellipsoid. - Parameters - ---------- - latitude : float or array - The geodetic latitude where the centrifugal potential will be - computed (in degrees). - height : float or array - The ellipsoidal (geometric) height of the computation point (in - meters). + Parameters + ---------- + latitude : float or array + The geodetic latitude where the centrifugal potential will be + computed (in degrees). + height : float or array + The ellipsoidal (geometric) height of the computation point (in + meters). - Returns - ------- - Phi : float or array - The centrifugal potential in m²/s². + Returns + ------- + Phi : float or array + The centrifugal potential in m²/s². - Notes - ----- + Notes + ----- - The centrifugal potential :math:`\Phi` at geodetic latitude - :math:`\phi` and height above the ellipsoid :math:`h` (geometric - height) is + The centrifugal potential :math:`\Phi` at geodetic latitude + :math:`\phi` and height above the ellipsoid :math:`h` (geometric + height) is - .. math:: + .. math:: - \Phi(\phi, h) = \dfrac{1}{2} - \omega^2 \left(N(\phi) + h\right)^2 \cos^2(\phi) + \Phi(\phi, h) = \dfrac{1}{2} + \omega^2 \left(N(\phi) + h\right)^2 \cos^2(\phi) - in which :math:`N(\phi)` is the prime vertical radius of curvature of - the ellipsoid. - """ - # Pre-compute to avoid repeated calculations - sinlat = np.sin(np.radians(latitude)) - coslat = np.sqrt(1 - sinlat**2) - - return (1 / 2) * ( - self.angular_velocity - * (self.prime_vertical_radius(sinlat) + height) - * coslat - ) ** 2 + in which :math:`N(\phi)` is the prime vertical radius of curvature of + the ellipsoid. + """ + # Pre-compute to avoid repeated calculations + sinlat = np.sin(np.radians(latitude)) + coslat = np.sqrt(1 - sinlat**2) + + return (1 / 2) * ( + self.angular_velocity + * (self.prime_vertical_radius(sinlat) + height) + * coslat + ) ** 2 diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py index 546e577b..fca5d34f 100644 --- a/boule/tests/test_ellipsoid.py +++ b/boule/tests/test_ellipsoid.py @@ -337,6 +337,24 @@ def test_normal_gravity_computed_on_internal_point(ellipsoid): assert len(warn) >= 1 +@pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES) +def test_normal_gravity_gravitational_centrifugal_potential(ellipsoid): + """ + Test that the normal gravity potential is equal to the sum of the normal + gravitational potential and centrifugal potential. + """ + latitude = np.array([-90, -45, 0, 45, 90]) + big_u0 = ellipsoid.normal_gravity_potential(latitude, height=0) + big_v0 = ellipsoid.normal_gravitational_potential(latitude, height=0) + big_phi0 = ellipsoid.centrifugal_potential(latitude, height=0) + height = 1000 + big_u = ellipsoid.normal_gravity_potential(latitude, height) + big_v = ellipsoid.normal_gravitational_potential(latitude, height) + big_phi = ellipsoid.centrifugal_potential(latitude, height) + npt.assert_allclose(big_u0, big_v0 + big_phi0) + npt.assert_allclose(big_u, big_v + big_phi) + + def test_emm_wgs84(): "The _emm property should be equal this value for WGS84" npt.assert_allclose(WGS84._emm, 0.00344978650684) From 62cb3fc745deb198f928acda1dc092f734d02ffa Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 11 Apr 2024 15:50:48 +0200 Subject: [PATCH 09/29] update for code-style --- boule/_ellipsoid.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index d5c76940..b757ebba 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -809,13 +809,13 @@ def normal_gravitational_potential(self, latitude, height): - 3 * u / self.linear_eccentricity ) - big_V = self.geocentric_grav_const / self.linear_eccentricity * np.arctan( + big_v = self.geocentric_grav_const / self.linear_eccentricity * np.arctan( self.linear_eccentricity / u ) + (1 / 3) * (self.angular_velocity * self.semimajor_axis) ** 2 * q / q_0 * ( 1.5 * np.sin(np.radians(beta)) ** 2 - 0.5 ) - return big_V + return big_v def normal_gravity_potential(self, latitude, height): r""" @@ -889,7 +889,7 @@ def normal_gravity_potential(self, latitude, height): - 3 * u / self.linear_eccentricity ) - big_U = ( + big_u = ( self.geocentric_grav_const / self.linear_eccentricity * np.arctan(self.linear_eccentricity / u) @@ -904,7 +904,7 @@ def normal_gravity_potential(self, latitude, height): * np.cos(np.radians(beta)) ** 2 ) - return big_U + return big_u def centrifugal_potential(self, latitude, height): r""" From 3bc08380b6aad807b92803d93036403a97c30631 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 11 Apr 2024 16:39:58 +0200 Subject: [PATCH 10/29] Fix sign of output lat in geodetic_to_ellipsoidal_harmonic --- boule/_ellipsoid.py | 4 +++- boule/tests/test_ellipsoid.py | 38 ++++++++++++++++++++++++++++++++++- 2 files changed, 40 insertions(+), 2 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index b757ebba..3850e501 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -572,7 +572,9 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): # Semiminor axis of ellipsoid passing through the computation point b_p = np.sqrt(r_p2 + z_p2 - self.linear_eccentricity**2 * cosbeta_p2) - beta_p = np.degrees(np.arccos(np.sqrt(cosbeta_p2))) + # Note that the sign of beta_p needs to be fixed as it is defined + # from -90 to 90 degrees. + beta_p = np.sign(latitude) * np.degrees(np.arccos(np.sqrt(cosbeta_p2))) return longitude, beta_p, b_p diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py index fca5d34f..f977e0b6 100644 --- a/boule/tests/test_ellipsoid.py +++ b/boule/tests/test_ellipsoid.py @@ -168,6 +168,35 @@ def test_spherical_to_geodetic_on_poles(ellipsoid): npt.assert_allclose(radius, height + ellipsoid.semiminor_axis, rtol=rtol) +@pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES) +def test_geodetic_to_ellipsoidal_conversions(ellipsoid): + "Test the geodetic to ellipsoidal-harmonic coordinate conversions by + going from geodetic to ellipsoidal to geodetic." + rtol = 1e-6 # The conversion is not too accurate for large height + size = 5 + geodetic_latitude_in = np.linspace(-90, 90, size) + height_in = np.zeros(size) + + longitude, reduced_latitude, u = ellipsoid.geodetic_to_ellipsoidal_harmonic( + None, geodetic_latitude_in, height_in + ) + longitude, geodetic_latitude_out, height_out = ellipsoid.ellipsoidal_harmonic_to_geodetic( + None, reduced_latitude, u + ) + npt.assert_allclose(geodetic_latitude_in, geodetic_latitude_out) + npt.assert_allclose(height_in, height_out) + + height_in = np.array(size * [1000]) + longitude, reduced_latitude, u = ellipsoid.geodetic_to_ellipsoidal_harmonic( + None, geodetic_latitude_in, height_in + ) + longitude, geodetic_latitude_out, height_out = ellipsoid.ellipsoidal_harmonic_to_geodetic( + None, reduced_latitude, u + ) + npt.assert_allclose(geodetic_latitude_in, geodetic_latitude_out, rtol=rtol) + npt.assert_allclose(height_in, height_out, rtol=rtol) + + @pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES) def test_spherical_to_geodetic_identity(ellipsoid): "Test if applying both conversions in series is the identity operator" @@ -329,12 +358,19 @@ def test_normal_gravity_against_somigliana(ellipsoid): @pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES) def test_normal_gravity_computed_on_internal_point(ellipsoid): """ - Check if warn is raised if height is negative + Check if warn is raised if height is negative for normal_gravity, + normal_gravity_potential, and normal_gravitational_potential. """ latitude = np.linspace(-90, 90, 100) with warnings.catch_warnings(record=True) as warn: ellipsoid.normal_gravity(latitude, height=-10) assert len(warn) >= 1 + with warnings.catch_warnings(record=True) as warn: + ellipsoid.normal_gravity_potential(latitude, height=-10) + assert len(warn) >= 1 + with warnings.catch_warnings(record=True) as warn: + ellipsoid.normal_gravitational_potential(latitude, height=-10) + assert len(warn) >= 1 @pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES) From 6cd8b3981302680e8db9773bffd00ba22fb50794 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 11 Apr 2024 23:32:06 +0200 Subject: [PATCH 11/29] Add normal gravity and gravitational potential to sphere class. --- boule/_ellipsoid.py | 13 ++-- boule/_sphere.py | 140 ++++++++++++++++++++++++++++++---- boule/tests/test_ellipsoid.py | 26 +++---- boule/tests/test_sphere.py | 20 +++++ 4 files changed, 163 insertions(+), 36 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 3850e501..af44037c 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -529,7 +529,7 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): The coordinate u, which is the semiminor axis of the ellipsoid that passes through the input coordinates. """ - if np.array(height).all() == 0: + if (np.array(height) == 0).all(): # Use simple formula that relates geodetic and reduced latitude beta = np.degrees( np.arctan( @@ -538,7 +538,7 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): * np.tan(np.radians(latitude)) ) ) - u = self.semiminor_axis + u = np.full_like(height, fill_value=self.semiminor_axis) return longitude, beta, u @@ -569,12 +569,13 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): # 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) - # Semiminor axis of ellipsoid passing through the computation point + # Semiminor axis of the ellipsoid passing through the computation + # point. This is the coordinate u b_p = np.sqrt(r_p2 + z_p2 - self.linear_eccentricity**2 * cosbeta_p2) # Note that the sign of beta_p needs to be fixed as it is defined - # from -90 to 90 degrees. - beta_p = np.sign(latitude) * np.degrees(np.arccos(np.sqrt(cosbeta_p2))) + # from -90 to 90 degrees, but arccos is from 0 to 180. + beta_p = np.copysign(np.degrees(np.arccos(np.sqrt(cosbeta_p2))), latitude) return longitude, beta_p, b_p @@ -940,7 +941,7 @@ def centrifugal_potential(self, latitude, height): \omega^2 \left(N(\phi) + h\right)^2 \cos^2(\phi) in which :math:`N(\phi)` is the prime vertical radius of curvature of - the ellipsoid. + the ellipsoid and :math:`\omega` is the angular velocity. """ # Pre-compute to avoid repeated calculations sinlat = np.sin(np.radians(latitude)) diff --git a/boule/_sphere.py b/boule/_sphere.py index 30c28cbb..398d65de 100644 --- a/boule/_sphere.py +++ b/boule/_sphere.py @@ -204,8 +204,8 @@ def normal_gravity(self, latitude, height, si_units=False): .. caution:: - The current implementation is only valid for heights on or above - the surface of the sphere. + These expressions are only valid for heights on or above the + surface of the sphere. Parameters ---------- @@ -292,10 +292,9 @@ def normal_gravity(self, latitude, height, si_units=False): } 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. + equilibrium. Therefore, unlike the oblate 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): @@ -380,9 +379,113 @@ def normal_gravitation(self, height, si_units=False): return gamma + def normal_gravity_potential(self, latitude, height): + r""" + Normal gravity potential of the sphere at the given latitude and + height. + + Computes the normal gravity potential (gravitational + centrifugal) + generated by the sphere at the given spherical latitude :math:`\theta` + and height above the surface of the sphere :math:`h`: + + .. math:: + + U(\theta, h) = V(h) + \Phi(\theta, h) = \dfrac{GM}{(R + h)} + + \dfrac{1}{2} \omega^2 \left(R + h\right)^2 \cos^2(\theta) + + 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:: + + These expressions are 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). + + Returns + ------- + U : float or array + The normal gravity potential in m²/s². + + Notes + ----- + A sphere under rotation is not in hydrostatic equilibrium. Therefore, + unlike the oblate ellipsoid, it is not its 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." + ) + + radial_distance = self.radius + height + big_u = self.geocentric_grav_const / radial_distance + big_phi = ( + 0.5 + * ( + self.angular_velocity + * (self.radius + height) + * np.cos(np.radians(latitude)) + ) + ** 2 + ) + + return big_u + big_phi + + def normal_gravitational_potential(self, height): + r""" + Normal gravitational potential at a given height above a sphere. + + .. math:: + + V(h) = \dfrac{GM}{(R + h)} + + in which :math:`R` is the sphere radius, :math:`h` is the height above + the sphere, :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). + + Returns + ------- + V : float or array + The normal gravitational potential in m²/s². + + """ + if np.any(height < 0): + warn( + "Formulas used are valid for points outside the sphere. " + "Height must be greater than or equal to zero." + ) + radial_distance = self.radius + height + return self.geocentric_grav_const / radial_distance + def centrifugal_potential(self, latitude, height): r""" - Centrifugal potential at the given latitude and height. + Centrifugal potential at the given latitude and height above the + sphere. Parameters ---------- @@ -400,18 +503,23 @@ def centrifugal_potential(self, latitude, height): Notes ----- - The centrifugal potential :math:`\Phi` at latitude :math:`\phi` and + The centrifugal potential :math:`\Phi` at latitude :math:`\theta` and height above the sphere :math:`h` is .. math:: - \Phi(\phi, h) = \dfrac{1}{2} - \omega^2 \left(R + h\right)^2 \cos^2(\phi) + \Phi(\theta, h) = \dfrac{1}{2} + \omega^2 \left(R + h\right)^2 \cos^2(\theta) - in which :math:`R` is the sphere radius. + in which :math:`R` is the sphere radius and :math:`\omega` is the + angular velocity. """ - return (1 / 2) * ( - self.angular_velocity - * (self.radius + height) - * np.cos(np.radians(latitude)) - ) ** 2 + return ( + 0.5 + * ( + self.angular_velocity + * (self.radius + height) + * np.cos(np.radians(latitude)) + ) + ** 2 + ) diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py index f977e0b6..6465c0fe 100644 --- a/boule/tests/test_ellipsoid.py +++ b/boule/tests/test_ellipsoid.py @@ -170,28 +170,29 @@ def test_spherical_to_geodetic_on_poles(ellipsoid): @pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES) def test_geodetic_to_ellipsoidal_conversions(ellipsoid): - "Test the geodetic to ellipsoidal-harmonic coordinate conversions by - going from geodetic to ellipsoidal to geodetic." - rtol = 1e-6 # The conversion is not too accurate for large height + """ + Test the geodetic to ellipsoidal-harmonic coordinate conversions by + going from geodetic to ellipsoidal and back. + """ size = 5 geodetic_latitude_in = np.linspace(-90, 90, size) height_in = np.zeros(size) - longitude, reduced_latitude, u = ellipsoid.geodetic_to_ellipsoidal_harmonic( None, geodetic_latitude_in, height_in ) - longitude, geodetic_latitude_out, height_out = ellipsoid.ellipsoidal_harmonic_to_geodetic( - None, reduced_latitude, u + longitude, geodetic_latitude_out, height_out = ( + ellipsoid.ellipsoidal_harmonic_to_geodetic(None, reduced_latitude, u) ) npt.assert_allclose(geodetic_latitude_in, geodetic_latitude_out) npt.assert_allclose(height_in, height_out) + rtol = 1e-5 # The conversion is not too accurate for large heights height_in = np.array(size * [1000]) longitude, reduced_latitude, u = ellipsoid.geodetic_to_ellipsoidal_harmonic( None, geodetic_latitude_in, height_in ) - longitude, geodetic_latitude_out, height_out = ellipsoid.ellipsoidal_harmonic_to_geodetic( - None, reduced_latitude, u + longitude, geodetic_latitude_out, height_out = ( + ellipsoid.ellipsoidal_harmonic_to_geodetic(None, reduced_latitude, u) ) npt.assert_allclose(geodetic_latitude_in, geodetic_latitude_out, rtol=rtol) npt.assert_allclose(height_in, height_out, rtol=rtol) @@ -379,15 +380,12 @@ def test_normal_gravity_gravitational_centrifugal_potential(ellipsoid): Test that the normal gravity potential is equal to the sum of the normal gravitational potential and centrifugal potential. """ - latitude = np.array([-90, -45, 0, 45, 90]) - big_u0 = ellipsoid.normal_gravity_potential(latitude, height=0) - big_v0 = ellipsoid.normal_gravitational_potential(latitude, height=0) - big_phi0 = ellipsoid.centrifugal_potential(latitude, height=0) - height = 1000 + size = 5 + latitude = np.array([np.linspace(-90, 90, size)] * 2) + height = np.array([[0] * size, [1000] * size]) big_u = ellipsoid.normal_gravity_potential(latitude, height) big_v = ellipsoid.normal_gravitational_potential(latitude, height) big_phi = ellipsoid.centrifugal_potential(latitude, height) - npt.assert_allclose(big_u0, big_v0 + big_phi0) npt.assert_allclose(big_u, big_v + big_phi) diff --git a/boule/tests/test_sphere.py b/boule/tests/test_sphere.py index bf1c9cfd..9b18b00a 100644 --- a/boule/tests/test_sphere.py +++ b/boule/tests/test_sphere.py @@ -56,6 +56,12 @@ def test_normal_gravity_computed_on_internal_point(sphere): with pytest.warns(UserWarning) as warn: sphere.normal_gravity(latitude, height=-10) assert len(warn) >= 1 + with warnings.catch_warnings(record=True) as warn: + sphere.normal_gravity_potential(latitude, height=-10) + assert len(warn) >= 1 + with warnings.catch_warnings(record=True) as warn: + sphere.normal_gravitational_potential(height=-10) + assert len(warn) >= 1 def test_check_geocentric_grav_const(): @@ -146,3 +152,17 @@ def test_normal_gravity_only_rotation(): expected_value, sphere.normal_gravity(latitude=45, height=height), ) + + +def test_normal_gravity_gravitational_centrifugal_potential(sphere): + """ + Test that the normal gravity potential is equal to the sum of the normal + gravitational potential and centrifugal potential. + """ + size = 5 + latitude = np.array([np.linspace(-90, 90, size)] * 2) + height = np.array([[0] * size, [1000] * size]) + big_u = sphere.normal_gravity_potential(latitude, height) + big_v = sphere.normal_gravitational_potential(height) + big_phi = sphere.centrifugal_potential(latitude, height) + npt.assert_allclose(big_u, big_v + big_phi) From 749bd1f8346cbf380be74bc63b58bebf909fbf24 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Fri, 12 Apr 2024 11:29:24 +0200 Subject: [PATCH 12/29] Fix dtype of output for geodetic_to_ellipsoidal_harmonic when height is an integer --- boule/_ellipsoid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index af44037c..c4e84df7 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -538,7 +538,7 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): * np.tan(np.radians(latitude)) ) ) - u = np.full_like(height, fill_value=self.semiminor_axis) + u = np.full_like(height, fill_value=self.semiminor_axis, dtype=np.float_) return longitude, beta, u From d3b34e921b77438e2d75adbb0085e14c1873381a Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Fri, 12 Apr 2024 12:00:26 +0200 Subject: [PATCH 13/29] Minor documation cleanup --- boule/_ellipsoid.py | 40 +++++++++++++++++++--------------------- boule/_sphere.py | 41 ++++++++++++++++++++--------------------- 2 files changed, 39 insertions(+), 42 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index c4e84df7..6a345ce7 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -607,7 +607,7 @@ def ellipsoidal_harmonic_to_geodetic(self, longitude, reduced_latitude, u): height : array Ellipsoidal heights in meters. """ - # semimajor axis of the ellipsoid that passes through the input + # Semimajor axis of the ellipsoid that passes through the input # coordinates a_p = np.sqrt(u**2 + self.linear_eccentricity**2) @@ -752,7 +752,7 @@ def normal_gravitational_potential(self, latitude, height): .. math:: - \V(\beta, u) = \dfrac{GM}{E} \arctan{\dfrac{E}{u}} + \dfrac{1}{3} + V(\beta, u) = \dfrac{GM}{E} \arctan{\dfrac{E}{u}} + \dfrac{1}{3} \omega^2 a^2 \dfrac{q}{q_0} P_2 (\sin \beta) in which :math:`V` is the gravitational potential of the @@ -762,7 +762,7 @@ def normal_gravitational_potential(self, latitude, height): 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 - and ellipsoidal height. See eq. 2-124 of [HofmannWellenhofMoritz2006]_. + 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. @@ -799,7 +799,7 @@ def normal_gravitational_potential(self, latitude, height): None, latitude, height ) - # compute the auxiliary functions q and q_0 (eq 2-113 of + # Compute the auxiliary functions q and q_0 (eq 2-113 of # HofmannWellenhofMoritz2006) q_0 = 0.5 * ( (1 + 3 * (self.semiminor_axis / self.linear_eccentricity) ** 2) @@ -831,7 +831,7 @@ def normal_gravity_potential(self, latitude, height): .. math:: - \U(\beta, u) = \dfrac{GM}{E} \arctan{\dfrac{E}{u}} + \dfrac{1}{2} + U(\beta, u) = \dfrac{GM}{E} \arctan{\dfrac{E}{u}} + \dfrac{1}{2} \omega^2 a^2 \dfrac{q}{q_0} \left(\sin^2 \beta - \dfrac{1}{3}\right) + \dfrac{1}{2} \omega^2 \left(u^2 + E^2) \cos^2 \beta @@ -864,7 +864,7 @@ def normal_gravity_potential(self, latitude, height): Returns ------- V : float or array - The normal gravitational potential in m²/s². + The normal gravity potential in m²/s². """ # Warn if height is negative @@ -879,7 +879,7 @@ def normal_gravity_potential(self, latitude, height): None, latitude, height ) - # compute the auxiliary functions q and q_0 (eq 2-113 of + # Compute the auxiliary functions q and q_0 (eq 2-113 of # HofmannWellenhofMoritz2006) q_0 = 0.5 * ( (1 + 3 * (self.semiminor_axis / self.linear_eccentricity) ** 2) @@ -914,6 +914,18 @@ def centrifugal_potential(self, latitude, height): Centrifugal potential at the given geodetic latitude and height above the ellipsoid. + The centrifugal potential :math:`\Phi` at geodetic latitude + :math:`\phi` and height above the ellipsoid :math:`h` (geometric + height) is + + .. math:: + + \Phi(\phi, h) = \dfrac{1}{2} + \omega^2 \left(N(\phi) + h\right)^2 \cos^2(\phi) + + in which :math:`N(\phi)` is the prime vertical radius of curvature of + the ellipsoid and :math:`\omega` is the angular velocity. + Parameters ---------- latitude : float or array @@ -928,20 +940,6 @@ def centrifugal_potential(self, latitude, height): Phi : float or array The centrifugal potential in m²/s². - Notes - ----- - - The centrifugal potential :math:`\Phi` at geodetic latitude - :math:`\phi` and height above the ellipsoid :math:`h` (geometric - height) is - - .. math:: - - \Phi(\phi, h) = \dfrac{1}{2} - \omega^2 \left(N(\phi) + h\right)^2 \cos^2(\phi) - - in which :math:`N(\phi)` is the prime vertical radius of curvature of - the ellipsoid and :math:`\omega` is the angular velocity. """ # Pre-compute to avoid repeated calculations sinlat = np.sin(np.radians(latitude)) diff --git a/boule/_sphere.py b/boule/_sphere.py index 398d65de..a99b4658 100644 --- a/boule/_sphere.py +++ b/boule/_sphere.py @@ -419,10 +419,9 @@ def normal_gravity_potential(self, latitude, height): Notes ----- A sphere under rotation is not in hydrostatic equilibrium. Therefore, - unlike the oblate ellipsoid, it is not its 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. + unlike the oblate 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): @@ -447,15 +446,17 @@ def normal_gravity_potential(self, latitude, height): def normal_gravitational_potential(self, height): r""" - Normal gravitational potential at a given height above a sphere. + Normal gravitational potential at the given height above a sphere. + + Computes the normal gravitational potential generated by the sphere + at the given height above the surface of the sphere :math:`h`: .. math:: V(h) = \dfrac{GM}{(R + h)} - in which :math:`R` is the sphere radius, :math:`h` is the height above - the sphere, :math:`G` is the gravitational constant, and :math:`M` is - the mass of the sphere. + in which :math:`R` is the sphere radius and :math:`GM` is the + geocentric gravitational constant of the sphere. .. caution:: @@ -487,6 +488,17 @@ def centrifugal_potential(self, latitude, height): Centrifugal potential at the given latitude and height above the sphere. + The centrifugal potential :math:`\Phi` at latitude :math:`\theta` and + height above the sphere :math:`h` is + + .. math:: + + \Phi(\theta, h) = \dfrac{1}{2} + \omega^2 \left(R + h\right)^2 \cos^2(\theta) + + in which :math:`R` is the sphere radius and :math:`\omega` is the + angular velocity. + Parameters ---------- latitude : float or array @@ -500,19 +512,6 @@ def centrifugal_potential(self, latitude, height): Phi : float or array The centrifugal potential in m²/s². - Notes - ----- - - The centrifugal potential :math:`\Phi` at latitude :math:`\theta` and - height above the sphere :math:`h` is - - .. math:: - - \Phi(\theta, h) = \dfrac{1}{2} - \omega^2 \left(R + h\right)^2 \cos^2(\theta) - - in which :math:`R` is the sphere radius and :math:`\omega` is the - angular velocity. """ return ( 0.5 From 6efc961a3fe6fd92ee58fc642bfee1fcb0e8d602 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 18 Apr 2024 11:08:42 +0200 Subject: [PATCH 14/29] Add centrifugal potential to triaxial ellipsoid --- boule/_triaxialellipsoid.py | 40 +++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/boule/_triaxialellipsoid.py b/boule/_triaxialellipsoid.py index d3fe2d6e..5cf0e882 100644 --- a/boule/_triaxialellipsoid.py +++ b/boule/_triaxialellipsoid.py @@ -295,3 +295,43 @@ def geocentric_radius(self, longitude, latitude, longitude_semimajor_axis=0.0): * np.cos(longitude_rad - longitude_semimajor_axis_rad) ** 2 ) return radius + + def centrifugal_potential(self, longitude, latitude, height): + r""" + Centrifugal potential at the given latitude and height above the + ellipsoid. + + The centrifugal potential :math:`\Phi` at spherical latitude + :math:`\phi`, spherical longitude :math:`\lambda` and spherical height + above the ellipsoid :math:`h` is + + .. math:: + + \Phi(\phi, \lambda, h) = \dfrac{1}{2} + \omega^2 \left(R(\phi, lambda) + h\right)^2 \cos^2(\phi) + + in which :math:`R(\phi, \lambda)` is the radius of the ellipsoid and + :math:`\omega` is the angular velocity. + + Parameters + ---------- + longitude : float or array + The spherical longitude where the centrifugal potential will be + computed (in degrees). + latitude : float or array + The spherical latitude where the centrifugal potential will be + computed (in degrees). + height : float or array + The spherical height of the computation point (in meters). + + Returns + ------- + Phi : float or array + The centrifugal potential in m²/s². + + """ + return (1 / 2) * ( + self.angular_velocity + * (self.geocentric_radius(longitude, latitude) + height) + * np.cos(np.radians(latitude)) + ) ** 2 From 53eb08602ccb552140cd1849e73e74901c56fcfa Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 18 Apr 2024 11:10:49 +0200 Subject: [PATCH 15/29] fix docstring --- boule/_triaxialellipsoid.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/boule/_triaxialellipsoid.py b/boule/_triaxialellipsoid.py index c61fed65..c1ab8a4e 100644 --- a/boule/_triaxialellipsoid.py +++ b/boule/_triaxialellipsoid.py @@ -415,8 +415,8 @@ def geocentric_radius(self, longitude, latitude): def centrifugal_potential(self, longitude, latitude, height): r""" - Centrifugal potential at the given latitude and height above the - ellipsoid. + Centrifugal potential at the given latitude, longitude and height above + the ellipsoid. The centrifugal potential :math:`\Phi` at spherical latitude :math:`\phi`, spherical longitude :math:`\lambda` and spherical height From 894d0ae30caa8311b578a919b5abc665ca2e96fd Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 18 Apr 2024 11:30:59 +0200 Subject: [PATCH 16/29] Add test for triaxial centrifugal potential --- boule/tests/test_triaxialellipsoid.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/boule/tests/test_triaxialellipsoid.py b/boule/tests/test_triaxialellipsoid.py index f97c21b6..d1575603 100644 --- a/boule/tests/test_triaxialellipsoid.py +++ b/boule/tests/test_triaxialellipsoid.py @@ -295,3 +295,26 @@ def test_geocentric_radius_biaxialellipsoid(triaxialellipsoid): npt.assert_allclose( radius_true, biaxialellipsoid.geocentric_radius(longitude, latitude) ) + + +def test_centrifugal_potential(triaxialellipsoid): + """ + Test that the centrifugal potential on the surface is greater along the + semimajor axis than the semimedium axis, and that the centrifugal potential + is zero at the poles. + """ + height = 2 * [0, 1, 1000] + latitude = 3 * [-90] + 3 * [90] + longitude = np.linspace(0, 360, 6) + assert ( + triaxialellipsoid.centrifugal_potential( + longitude=longitude, latitude=latitude, height=height + ) + < 1e-15 + ).all() + assert ( + triaxialellipsoid.centrifugal_potential(longitude=0, latitude=0, height=height) + >= triaxialellipsoid.centrifugal_potential( + longitude=90, latitude=0, height=height + ) + ).all() From 48c07635d0b6fc3cd7adf93416a946403ea72807 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Tue, 30 Apr 2024 12:41:28 +0200 Subject: [PATCH 17/29] Add small flattening approximation for Ellipsoid.reference_normal_gravity_potential --- boule/_ellipsoid.py | 45 +++++++++++++++++++++++---------------------- 1 file changed, 23 insertions(+), 22 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 4aa6c2d1..11917475 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -24,9 +24,9 @@ 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. @@ -65,9 +65,13 @@ class Ellipsoid: .. caution:: - Use :class:`boule.Sphere` if you desire zero flattening because there - are singularities for this particular case in the normal gravity - calculations. + 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. 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): From 2122f533235520fc14492b6d7b4ff0424ef94dc7 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Tue, 30 Apr 2024 13:06:11 +0200 Subject: [PATCH 18/29] Add small flattening approximation for Ellipsoid.gravity_equator --- boule/_ellipsoid.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 11917475..1d6bac08 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -363,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 From 4da541f763e2f754c3e87b151a100f8e4b21aae7 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Tue, 30 Apr 2024 13:10:49 +0200 Subject: [PATCH 19/29] Add small flattening approximation for Ellipsoid.gravity_pole --- boule/_ellipsoid.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 1d6bac08..bdd1ff04 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -391,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 From f6c09c1130664e55f1b021c0ff2646e7f01c307d Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Tue, 30 Apr 2024 13:39:41 +0200 Subject: [PATCH 20/29] Add small flattening approximation for Ellipsoid.normal_gravitational_potential --- boule/_ellipsoid.py | 53 ++++++++++++++++++++++++++++++--------------- 1 file changed, 36 insertions(+), 17 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index bdd1ff04..32a87547 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -927,29 +927,48 @@ 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 ) - # 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 - self.linear_eccentricity**2 / (3 * u**3) + ) + else: + var = ( + self.geocentric_grav_const + / self.linear_eccentricity + * np.arctan(self.linear_eccentricity / 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 / self.linear_eccentricity) ** 2) - * np.arctan2(self.linear_eccentricity, self.semiminor_axis) - - 3 * self.semiminor_axis / self.linear_eccentricity - ) - q = 0.5 * ( - (1 + 3 * (u / self.linear_eccentricity) ** 2) - * np.arctan2(self.linear_eccentricity, u) - - 3 * u / self.linear_eccentricity - ) + 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 / self.linear_eccentricity) ** 2) + * np.arctan2(self.linear_eccentricity, self.semiminor_axis) + - 3 * self.semiminor_axis / self.linear_eccentricity + ) + q = 0.5 * ( + (1 + 3 * (u / self.linear_eccentricity) ** 2) + * np.arctan2(self.linear_eccentricity, u) + - 3 * u / self.linear_eccentricity + ) + ratio = q / q_0 - big_v = self.geocentric_grav_const / self.linear_eccentricity * np.arctan( - self.linear_eccentricity / u - ) + (1 / 3) * (self.angular_velocity * self.semimajor_axis) ** 2 * q / q_0 * ( - 1.5 * np.sin(np.radians(beta)) ** 2 - 0.5 - ) + big_v = var + (1 / 3) * ( + self.angular_velocity * self.semimajor_axis + ) ** 2 * ratio * (1.5 * np.sin(np.radians(beta)) ** 2 - 0.5) return big_v From c9744d60dd8c2ab07c98dd731003a10b842f01d8 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Tue, 30 Apr 2024 13:43:21 +0200 Subject: [PATCH 21/29] Add small flattening approximation for Ellipsoid.normal_gravity_potential --- boule/_ellipsoid.py | 52 ++++++++++++++++++++++++++++++--------------- 1 file changed, 35 insertions(+), 17 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 32a87547..1d52e786 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -1026,32 +1026,50 @@ 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 ) - # 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 - self.linear_eccentricity**2 / (3 * u**3) + ) + else: + var = ( + self.geocentric_grav_const + / self.linear_eccentricity + * np.arctan(self.linear_eccentricity / 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 / self.linear_eccentricity) ** 2) - * np.arctan2(self.linear_eccentricity, self.semiminor_axis) - - 3 * self.semiminor_axis / self.linear_eccentricity - ) - q = 0.5 * ( - (1 + 3 * (u / self.linear_eccentricity) ** 2) - * np.arctan2(self.linear_eccentricity, u) - - 3 * u / self.linear_eccentricity - ) + 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 / self.linear_eccentricity) ** 2) + * np.arctan2(self.linear_eccentricity, self.semiminor_axis) + - 3 * self.semiminor_axis / self.linear_eccentricity + ) + q = 0.5 * ( + (1 + 3 * (u / self.linear_eccentricity) ** 2) + * np.arctan2(self.linear_eccentricity, u) + - 3 * u / self.linear_eccentricity + ) + ratio = q / q_0 big_u = ( - self.geocentric_grav_const - / self.linear_eccentricity - * np.arctan(self.linear_eccentricity / 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 From 26342067224a0092381bdac0d3aedf260f34be7d Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 2 May 2024 18:37:59 +0200 Subject: [PATCH 22/29] Add small flattening approximation to Ellipsoid.geodetic_to_ellipsoidal_harmonic --- boule/_ellipsoid.py | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 1d52e786..591a0e27 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -694,13 +694,31 @@ 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 - - # Auxiliary variables - big_d = (r_p2 - z_p2) / self.linear_eccentricity**2 - big_r = (r_p2 + z_p2) / self.linear_eccentricity**2 - - # 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) + # 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 + + self.linear_eccentricity**2 + * 0.25 + * (z_pp2**2 / r_pp2**3 - 1 / r_pp2) + ) + else: + # Auxiliary variables + big_d = z_pp2 / self.linear_eccentricity**2 + big_r = r_pp2 / self.linear_eccentricity**2 + # 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 sqrt 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 From e4c95f474e8119d08791d6c9db3ccdd672a93a81 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 2 May 2024 22:47:31 +0200 Subject: [PATCH 23/29] Add small flattening approximation to Ellipsoid.normal_gravity --- boule/_ellipsoid.py | 63 +++++++++++++++++++++++++++++++-------------- 1 file changed, 43 insertions(+), 20 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 591a0e27..bb809341 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -716,8 +716,8 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): 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 sqrt of a negative number, - # use the absolute value of this quantity. + # 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 @@ -825,34 +825,52 @@ def normal_gravity(self, latitude, height, si_units=False): "Height must be greater than or equal to zero." ) - # Pre-compute to avoid repeated calculations + # The variable names follow Li and Goetze (2001). The prime terms + # (*_p) refer to quantities on an ellipsoid passing through the + # computation point. sinlat = np.sin(np.radians(latitude)) coslat = np.sqrt(1 - sinlat**2) - # The terms below follow the variable names from Li and Goetze (2001). - # The prime terms (*_p) refer to quantities on an ellipsoid passing - # through the computation point. - - # The reduced latitude of the projection of the point on the ellipsoid + # Reduced latitude of the projection of the point on the + # reference ellipsoid beta = np.arctan2(self.semiminor_axis * sinlat, self.semimajor_axis * coslat) sinbeta = np.sin(beta) cosbeta = np.sqrt(1 - sinbeta**2) - # Distance between the computation point and the equatorial plane + # Distance squared between computation point and equatorial plane z_p2 = (self.semiminor_axis * sinbeta + height * sinlat) ** 2 - # Distance between the computation point and the spin axis + # 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 - # Auxiliary variables - big_d = (r_p2 - z_p2) / self.linear_eccentricity**2 - big_r = (r_p2 + z_p2) / self.linear_eccentricity**2 - - # Reduced latitude of the computation point - cosbeta_p2 = 0.5 + big_r / 2 - np.sqrt(0.25 + big_r**2 / 4 - big_d / 2) + if self.flattening < 1.25e-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 + + self.linear_eccentricity**2 * 0.25 * (z_pp2**2 / r_pp2**3 - 1 / r_pp2) + ) + else: + # Auxiliary variables + big_d = z_pp2 / self.linear_eccentricity**2 + big_r = r_pp2 / self.linear_eccentricity**2 + # 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) sinbeta_p2 = 1 - cosbeta_p2 - # Auxiliary variables + # Semiminor axis of the ellipsoid passing through the computation + # point. This is the coordinate u b_p = np.sqrt(r_p2 + z_p2 - self.linear_eccentricity**2 * cosbeta_p2) + + # Auxiliary variables q_0 = 0.5 * ( (1 + 3 * (self.semiminor_axis / self.linear_eccentricity) ** 2) * np.arctan2(self.linear_eccentricity, self.semiminor_axis) @@ -874,14 +892,19 @@ def normal_gravity(self, latitude, height, si_units=False): / (b_p**2 + self.linear_eccentricity**2) ) + # Compute the ratio E q_p / q_0 + if self.flattening <= 1.25e-5: + aux = 3 * self.semiminor_axis**3 / b_p**2 + else: + aux = self.linear_eccentricity * q_p / q_0 + # Put together gamma using 3 separate terms term1 = self.geocentric_grav_const / (b_p**2 + self.linear_eccentricity**2) term2 = (0.5 * sinbeta_p2 - 1 / 6) * ( self.semimajor_axis**2 - * self.linear_eccentricity - * q_p + * aux * self.angular_velocity**2 - / ((b_p**2 + self.linear_eccentricity**2) * q_0) + / (b_p**2 + self.linear_eccentricity**2) ) term3 = -cosbeta_p2 * b_p * self.angular_velocity**2 gamma = (term1 + term2 + term3) / big_w From e1b650afaae24bb3e81f1c7cf08233538de532ae Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Thu, 2 May 2024 23:24:26 +0200 Subject: [PATCH 24/29] Add tests for small flattenings --- boule/_ellipsoid.py | 45 +++++++++++++++++------------------ boule/tests/test_ellipsoid.py | 23 ++++++++++++++++++ 2 files changed, 45 insertions(+), 23 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index bb809341..9aaaefd3 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -845,7 +845,7 @@ def normal_gravity(self, latitude, height, si_units=False): z_pp2 = r_p2 - z_p2 r_pp2 = r_p2 + z_p2 - if self.flattening < 1.25e-5: + if self.flattening < 5.e-5: # Use the Taylor series approximation for flattenings close to # zero to avoid numerical issues. cosbeta_p2 = ( @@ -870,34 +870,33 @@ def normal_gravity(self, latitude, height, si_units=False): # point. This is the coordinate u b_p = np.sqrt(r_p2 + z_p2 - self.linear_eccentricity**2 * cosbeta_p2) - # Auxiliary variables - q_0 = 0.5 * ( - (1 + 3 * (self.semiminor_axis / self.linear_eccentricity) ** 2) - * np.arctan2(self.linear_eccentricity, self.semiminor_axis) - - 3 * self.semiminor_axis / self.linear_eccentricity - ) - q_p = ( - 3 - * (1 + (b_p / self.linear_eccentricity) ** 2) - * ( - 1 - - b_p - / self.linear_eccentricity - * np.arctan2(self.linear_eccentricity, b_p) + # Compute auxiliary variables and the ratio E q_p / q_0 + if self.flattening <= 1.25e-5: + aux = 3 * self.semiminor_axis**3 / b_p**2 + else: + q_0 = 0.5 * ( + (1 + 3 * (self.semiminor_axis / self.linear_eccentricity) ** 2) + * np.arctan2(self.linear_eccentricity, self.semiminor_axis) + - 3 * self.semiminor_axis / self.linear_eccentricity ) - - 1 - ) + q_p = ( + 3 + * (1 + (b_p / self.linear_eccentricity) ** 2) + * ( + 1 + - b_p + / self.linear_eccentricity + * np.arctan2(self.linear_eccentricity, b_p) + ) + - 1 + ) + aux = self.linear_eccentricity * q_p / q_0 + big_w = np.sqrt( (b_p**2 + self.linear_eccentricity**2 * sinbeta_p2) / (b_p**2 + self.linear_eccentricity**2) ) - # Compute the ratio E q_p / q_0 - if self.flattening <= 1.25e-5: - aux = 3 * self.semiminor_axis**3 / b_p**2 - else: - aux = self.linear_eccentricity * q_p / q_0 - # Put together gamma using 3 separate terms term1 = self.geocentric_grav_const / (b_p**2 + self.linear_eccentricity**2) term2 = (0.5 * sinbeta_p2 - 1 / 6) * ( diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py index ec9be376..ef802429 100644 --- a/boule/tests/test_ellipsoid.py +++ b/boule/tests/test_ellipsoid.py @@ -392,3 +392,26 @@ 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) + + +def test_small_flattenings(): + "Check that no errors arise when using small values for the flattening" + flattening = 10.**(-np.arange(1, 20, .5)) + latitude = np.linspace(-90, 90, 13 ) + height = np.array([[1.e3] * len(latitude)]) + with warnings.catch_warnings(record=True) as warn: + for f in flattening: + ellipsoid = Ellipsoid( + name="WGS84 with small flattening", + semimajor_axis=WGS84.semimajor_axis, + flattening=f, + geocentric_grav_const=WGS84.geocentric_grav_const, + angular_velocity=WGS84.angular_velocity, + ) + big_u_0 = ellipsoid.reference_normal_gravity_potential + g_b = ellipsoid.gravity_pole + g_a = ellipsoid.gravity_equator + big_v = ellipsoid.normal_gravitational_potential(latitude, height) + big_u = ellipsoid.normal_gravity_potential(latitude, height) + gamma = ellipsoid.normal_gravity(latitude, height) + assert warn == None From bdf750d972b6c787f7a7cacd897e9e9c8b41c2df Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Fri, 3 May 2024 14:23:23 +0200 Subject: [PATCH 25/29] Update normal gravity web documentation --- boule/_ellipsoid.py | 2 +- boule/tests/test_ellipsoid.py | 37 ++++++--------------- doc/user_guide/normal_gravity.rst | 55 ++++++++++++++++++++++--------- 3 files changed, 51 insertions(+), 43 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 9aaaefd3..4746edd7 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -845,7 +845,7 @@ def normal_gravity(self, latitude, height, si_units=False): z_pp2 = r_p2 - z_p2 r_pp2 = r_p2 + z_p2 - if self.flattening < 5.e-5: + if self.flattening < 5.0e-5: # Use the Taylor series approximation for flattenings close to # zero to avoid numerical issues. cosbeta_p2 = ( diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py index ef802429..9260d2d7 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(): @@ -396,9 +379,9 @@ def test_emm_wgs84(): def test_small_flattenings(): "Check that no errors arise when using small values for the flattening" - flattening = 10.**(-np.arange(1, 20, .5)) - latitude = np.linspace(-90, 90, 13 ) - height = np.array([[1.e3] * len(latitude)]) + flattening = 10.0 ** (-np.arange(1, 20, 0.5)) + latitude = np.linspace(-90, 90, 13) + height = np.array([[1.0e3] * len(latitude)]) with warnings.catch_warnings(record=True) as warn: for f in flattening: ellipsoid = Ellipsoid( @@ -408,10 +391,10 @@ def test_small_flattenings(): geocentric_grav_const=WGS84.geocentric_grav_const, angular_velocity=WGS84.angular_velocity, ) - big_u_0 = ellipsoid.reference_normal_gravity_potential - g_b = ellipsoid.gravity_pole - g_a = ellipsoid.gravity_equator - big_v = ellipsoid.normal_gravitational_potential(latitude, height) - big_u = ellipsoid.normal_gravity_potential(latitude, height) - gamma = ellipsoid.normal_gravity(latitude, height) - assert warn == None + _ = 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) + assert len(warn) == 0 diff --git a/doc/user_guide/normal_gravity.rst b/doc/user_guide/normal_gravity.rst index 25559bb2..78a9a557 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:`1-^{-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`: From 965e9036b1f9a9e60ebd6fbf5733fa6a58ba04dc Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Wed, 9 Oct 2024 16:06:43 +0200 Subject: [PATCH 26/29] Fix typo in docstring --- boule/_ellipsoid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index a5123392..e0909c8a 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -876,7 +876,7 @@ def normal_gravitational_potential(self, latitude, height): in which :math:`V` is the gravitational potential of the ellipsoid (no centrifugal term), :math:`GM` is the geocentric - gravitational constant, :math:`E` is the linear eccentricty, + 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- From af53cf95c4530ab5819241aa79377b20c2253d58 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 10 Oct 2024 10:21:09 -0700 Subject: [PATCH 27/29] Move caution above parameters in Ellipsoid docstring --- boule/_ellipsoid.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index e0909c8a..33ef076c 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -33,6 +33,17 @@ class Ellipsoid: **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,17 +73,6 @@ class Ellipsoid: comments : str or None Additional comments regarding the ellipsoid (optional). - - .. 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. - Examples -------- From 40513e39226bf0fc4f32fb05102bb0c178956d19 Mon Sep 17 00:00:00 2001 From: markwieczorek Date: Thu, 10 Oct 2024 21:38:19 +0200 Subject: [PATCH 28/29] Update doc/user_guide/normal_gravity.rst Fix typo in docstring Co-authored-by: Santiago Soler --- doc/user_guide/normal_gravity.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/user_guide/normal_gravity.rst b/doc/user_guide/normal_gravity.rst index 78a9a557..ddcac00a 100644 --- a/doc/user_guide/normal_gravity.rst +++ b/doc/user_guide/normal_gravity.rst @@ -119,7 +119,7 @@ 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:`1-^{-7}`, boule makes use of approximate equations in the low flattening +: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 From 4c17377f769373f31f10f64481f43ca951235267 Mon Sep 17 00:00:00 2001 From: markwieczorek Date: Thu, 10 Oct 2024 21:42:02 +0200 Subject: [PATCH 29/29] Update boule/tests/test_ellipsoid.py Improve test function using parameterize and simplefilter Co-authored-by: Santiago Soler --- boule/tests/test_ellipsoid.py | 41 +++++++++++++++++++---------------- 1 file changed, 22 insertions(+), 19 deletions(-) diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py index 9260d2d7..27ba909e 100644 --- a/boule/tests/test_ellipsoid.py +++ b/boule/tests/test_ellipsoid.py @@ -377,24 +377,27 @@ def test_emm_wgs84(): npt.assert_allclose(WGS84._emm, 0.00344978650684) -def test_small_flattenings(): - "Check that no errors arise when using small values for the flattening" - flattening = 10.0 ** (-np.arange(1, 20, 0.5)) +@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(record=True) as warn: - for f in flattening: - ellipsoid = Ellipsoid( - name="WGS84 with small flattening", - semimajor_axis=WGS84.semimajor_axis, - flattening=f, - 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) - assert len(warn) == 0 + 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)