From 6c716ac984749659254a36b02c9aac1b72d46685 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Mon, 8 Apr 2024 14:29:49 +0200 Subject: [PATCH 01/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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/22] 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 aa5ca212793fe3732c478a9f0f7df696671b15e9 Mon Sep 17 00:00:00 2001 From: Leonardo Uieda Date: Fri, 7 Jun 2024 15:41:08 -0300 Subject: [PATCH 17/22] Fix equations in docstrings of potential methods --- boule/_ellipsoid.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 4aa6c2d1..5302cc07 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -874,12 +874,13 @@ 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- - 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. @@ -950,15 +951,15 @@ def normal_gravity_potential(self, latitude, height): 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) + - \dfrac{1}{3}\right) + \dfrac{1}{2} \omega^2 \left(u^2 + E^2\right) \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, + linear eccentricity, :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 + :math:`\beta` are ellipsoidal-harmonic coordinates corresponding to the + input geodetic latitude and ellipsoidal height. See eq. 2-126 of [HofmannWellenhofMoritz2006]_. Assumes that the internal density distribution of the ellipsoid is such From 7fc90719a0c927ecb18d878cd5b8d68cd1ea8d3b Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Sat, 8 Jun 2024 17:23:48 +0200 Subject: [PATCH 18/22] precompute variables and remove redundant code in normal grav computations --- boule/_ellipsoid.py | 120 +++++++++++++++----------------------------- 1 file changed, 40 insertions(+), 80 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 5302cc07..c8de1463 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -665,6 +665,7 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): # computation point. sinlat = np.sin(np.radians(latitude)) coslat = np.sqrt(1 - sinlat**2) + big_e2 = self.linear_eccentricity**2 # Reduced latitude of the projection of the point on the # reference ellipsoid @@ -680,15 +681,15 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): 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 + 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) # 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) + b_p = np.sqrt(r_p2 + z_p2 - big_e2 * cosbeta_p2) # Note that the sign of beta_p needs to be fixed as it is defined # from -90 to 90 degrees, but arccos is from 0 to 180. @@ -791,65 +792,34 @@ def normal_gravity(self, latitude, height, si_units=False): "Height must be greater than or equal to zero." ) - # Pre-compute to avoid repeated calculations - 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 - 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 - z_p2 = (self.semiminor_axis * sinbeta + height * sinlat) ** 2 - # Distance between the computation point and the 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 - - # Reduced latitude of the computation point - cosbeta_p2 = 0.5 + big_r / 2 - np.sqrt(0.25 + big_r**2 / 4 - big_d / 2) - sinbeta_p2 = 1 - cosbeta_p2 + # Convert geodetic latitude and height to ellipsoidal-harmonic coords + longitude, beta, u = self.geodetic_to_ellipsoidal_harmonic( + None, latitude, height + ) + sinbeta2 = np.sin(np.radians(beta)) ** 2 + cosbeta2 = 1 - sinbeta2 + big_e = self.linear_eccentricity - # Auxiliary variables - b_p = np.sqrt(r_p2 + z_p2 - self.linear_eccentricity**2 * cosbeta_p2) + # 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_p = ( - 3 - * (1 + (b_p / self.linear_eccentricity) ** 2) - * ( - 1 - - b_p - / self.linear_eccentricity - * np.arctan2(self.linear_eccentricity, b_p) - ) - - 1 - ) - big_w = np.sqrt( - (b_p**2 + self.linear_eccentricity**2 * sinbeta_p2) - / (b_p**2 + self.linear_eccentricity**2) + (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)) # 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) * ( + term1 = self.geocentric_grav_const / (u**2 + big_e**2) + term2 = (0.5 * sinbeta2 - 1 / 6) * ( self.semimajor_axis**2 - * self.linear_eccentricity + * big_e * q_p * self.angular_velocity**2 - / ((b_p**2 + self.linear_eccentricity**2) * q_0) + / ((u**2 + big_e**2) * q_0) ) - term3 = -cosbeta_p2 * b_p * self.angular_velocity**2 + term3 = -cosbeta2 * u * self.angular_velocity**2 gamma = (term1 + term2 + term3) / big_w # Convert gamma from SI to mGal @@ -916,25 +886,20 @@ def normal_gravitational_potential(self, latitude, height): longitude, beta, u = self.geodetic_to_ellipsoidal_harmonic( None, latitude, height ) + big_e = self.linear_eccentricity # 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 + (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) - 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 = self.geocentric_grav_const / big_e * np.arctan(big_e / 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 @@ -951,8 +916,8 @@ def normal_gravity_potential(self, latitude, height): 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\right) - \cos^2 \beta + - \dfrac{1}{3}\right) + \dfrac{1}{2} \omega^2 \left(u^2 + + E^2\right) \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 @@ -996,24 +961,19 @@ def normal_gravity_potential(self, latitude, height): longitude, beta, u = self.geodetic_to_ellipsoidal_harmonic( None, latitude, height ) + big_e = self.linear_eccentricity # 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 + (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) big_u = ( - self.geocentric_grav_const - / self.linear_eccentricity - * np.arctan(self.linear_eccentricity / u) + self.geocentric_grav_const / big_e * np.arctan(big_e / u) + 0.5 * (self.angular_velocity * self.semimajor_axis) ** 2 * q @@ -1021,7 +981,7 @@ def normal_gravity_potential(self, latitude, height): * (np.sin(np.radians(beta)) ** 2 - 1 / 3) + 0.5 * self.angular_velocity**2 - * (u**2 + self.linear_eccentricity**2) + * (u**2 + big_e**2) * np.cos(np.radians(beta)) ** 2 ) From 452318d5967dc58322643f825ee9151ed9c3467c Mon Sep 17 00:00:00 2001 From: markwieczorek Date: Tue, 2 Jul 2024 16:41:45 +0200 Subject: [PATCH 19/22] Update boule/_ellipsoid.py Fix typo in documentation Co-authored-by: Santiago Soler --- boule/_ellipsoid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index c8de1463..c07a24ce 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -623,7 +623,7 @@ 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]. + are converted following [Lakshmanan1991]_ and [LiGotze2001]_. Parameters ---------- From 1df518d4c5963ba52501ab4e7388544de399e470 Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Tue, 2 Jul 2024 16:49:52 +0200 Subject: [PATCH 20/22] Update np.float_ to np.float64 form numpy 2.0 --- boule/_ellipsoid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index c07a24ce..a746538f 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -655,7 +655,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, dtype=np.float_) + u = np.full_like(height, fill_value=self.semiminor_axis, dtype=np.float64) return longitude, beta, u From dd465a748d34cff9c658c5df7a092cfb61ed5ea9 Mon Sep 17 00:00:00 2001 From: markwieczorek Date: Tue, 27 Aug 2024 21:54:20 +0200 Subject: [PATCH 21/22] Update boule/_ellipsoid.py Co-authored-by: Santiago Soler --- boule/_ellipsoid.py | 61 ++++++++++++++++++++++----------------------- 1 file changed, 30 insertions(+), 31 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index a746538f..1b4912bb 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -659,43 +659,42 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): 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 - # computation point. - sinlat = np.sin(np.radians(latitude)) - coslat = np.sqrt(1 - sinlat**2) - big_e2 = self.linear_eccentricity**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) + # 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) + big_e2 = self.linear_eccentricity**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 + # 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) / big_e2 - big_r = (r_p2 + z_p2) / big_e2 + # 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) + # 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 the ellipsoid passing through the computation - # point. This is the coordinate u - b_p = np.sqrt(r_p2 + z_p2 - big_e2 * cosbeta_p2) + # Semiminor axis of the ellipsoid passing through the computation + # point. This is the coordinate u + b_p = np.sqrt(r_p2 + z_p2 - big_e2 * cosbeta_p2) - # Note that the sign of beta_p needs to be fixed as it is defined - # 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) + # Note that the sign of beta_p needs to be fixed as it is defined + # 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 + return longitude, beta_p, b_p def ellipsoidal_harmonic_to_geodetic(self, longitude, reduced_latitude, u): """ From 5744b5a50d9f0be76a328a4cfe7e5662fcc3dfcf Mon Sep 17 00:00:00 2001 From: Mark Wieczorek Date: Tue, 27 Aug 2024 22:08:28 +0200 Subject: [PATCH 22/22] format code --- boule/_ellipsoid.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/boule/_ellipsoid.py b/boule/_ellipsoid.py index 1b4912bb..6568d029 100644 --- a/boule/_ellipsoid.py +++ b/boule/_ellipsoid.py @@ -668,9 +668,7 @@ def geodetic_to_ellipsoidal_harmonic(self, longitude, latitude, height): # Reduced latitude of the projection of the point on the # reference ellipsoid - beta = np.arctan2( - self.semiminor_axis * sinlat, self.semimajor_axis * coslat - ) + beta = np.arctan2(self.semiminor_axis * sinlat, self.semimajor_axis * coslat) sinbeta = np.sin(beta) cosbeta = np.sqrt(1 - sinbeta**2)