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)