Merge pull request #1547 from astrojuanlu/fix-spheroid-location-take-2

Fix spheroid location, take 2

docs/source/examples/Detecting Events.myst.md

@@ -163,7 +163,7 @@ from poliastro.core.spheroid_location import cartesian_to_ellipsoidal
 latitudes = []
 for r in rr:
     position_on_body = (r / norm(r)) * Earth.R
-    lat, _, _ = cartesian_to_ellipsoidal(
+    _, lat, _ = cartesian_to_ellipsoidal(
         Earth.R, Earth.R_polar, *position_on_body
     )
     latitudes.append(np.rad2deg(lat))

src/poliastro/core/spheroid_location.py

@@ -7,45 +7,46 @@
 
 
 @jit
-def cartesian_cords(_a, _c, _lon, _lat, _h):
+def cartesian_cords(a, c, lon, lat, h):
     """Calculates cartesian coordinates.
 
     Parameters
     ----------
-    _a : float
+    a : float
         Semi-major axis
-    _c : float
+    c : float
         Semi-minor axis
-    _lon : float
+    lon : float
         Geodetic longitude
-    _lat : float
+    lat : float
         Geodetic latitude
-    _h : float
+    h : float
         Geodetic height
 
     """
-    e2 = 1 - (_c / _a) ** 2
-    N = _a / np.sqrt(1 - e2 * np.sin(_lon) ** 2)
 
-    x = (N + _h) * np.cos(_lon) * np.cos(_lat)
-    y = (N + _h) * np.cos(_lon) * np.sin(_lat)
-    z = ((1 - e2) * N + _h) * np.sin(_lon)
+    e2 = 1 - (c / a) ** 2
+    N = a / np.sqrt(1 - e2 * np.sin(lat) ** 2)
+
+    x = (N + h) * np.cos(lat) * np.cos(lon)
+    y = (N + h) * np.cos(lat) * np.sin(lon)
+    z = ((1 - e2) * N + h) * np.sin(lat)
     return x, y, z
 
 
 @jit
-def f(_a, _c):
+def f(a, c):
     """Get first flattening.
 
     Parameters
     ----------
-    _a : float
+    a : float
         Semi-major axis
-    _c : float
+    c : float
         Semi-minor axis
 
     """
-    return 1 - _c / _a
+    return 1 - c / a
 
 
 @jit
@@ -89,21 +90,21 @@ def tangential_vecs(N):
 
 
 @jit
-def radius_of_curvature(_a, _c, _lat):
+def radius_of_curvature(a, c, lat):
     """Radius of curvature of the meridian at the latitude of the given location.
 
     Parameters
     ----------
-    _a : float
+    a : float
         Semi-major axis
-    _c : float
+    c : float
         Semi-minor axis
-    _lat : float
+    lat : float
         Geodetic latitude
 
     """
-    e2 = 1 - (_c / _a) ** 2
-    rc = _a * (1 - e2) / (1 - e2 * np.sin(_lat) ** 2) ** 1.5
+    e2 = 1 - (c / a) ** 2
+    rc = a * (1 - e2) / (1 - e2 * np.sin(lat) ** 2) ** 1.5
     return rc
 
 
@@ -156,17 +157,17 @@ def is_visible(cartesian_cords, px, py, pz, N):
 
 
 @jit
-def cartesian_to_ellipsoidal(_a, _c, x, y, z):
+def cartesian_to_ellipsoidal(a, c, x, y, z):
     """
     Converts cartesian coordinates to ellipsoidal coordinates for the given ellipsoid.
     Instead of the iterative formula, the function uses the approximation introduced in
-    Bowring, B. R. (1976). TRANSFORMATION FROM SPATIAL TO GEOGRAPHICAL COORDINATES
+    Bowring, B. R. (1976). TRANSFORMATION FROM SPATIAL TO GEOGRAPHICAL COORDINATES.
 
     Parameters
     ----------
-    _a : float
+    a : float
         Semi-major axis
-    _c : float
+    c : float
         Semi-minor axis
     x : float
         x coordinate
@@ -176,18 +177,22 @@ def cartesian_to_ellipsoidal(_a, _c, x, y, z):
         z coordinate
 
     """
-    e2 = 1 - (_c / _a) ** 2
+    e2 = 1 - (c / a) ** 2
     e2_ = e2 / (1 - e2)
     p = np.sqrt(x**2 + y**2)
-    th = np.arctan(z * _a / (p * _c))
+    th = np.arctan(z * a / (p * c))
     lon = np.arctan2(
         y, x
     )  # Use `arctan2` so that lon lies in the range: [-pi, +pi]
     lat = np.arctan(
-        (z + e2_ * _c * np.sin(th) ** 3) / (p - e2 * _a * np.cos(th) ** 3)
+        (z + e2_ * c * np.sin(th) ** 3) / (p - e2 * a * np.cos(th) ** 3)
     )
 
-    v = _a / np.sqrt(1 - e2 * np.sin(lat) ** 2)
-    h = x / np.cos(lat) - v if lat == 0.0 else z / np.sin(lat) - (1 - e2) * v
+    v = a / np.sqrt(1 - e2 * np.sin(lat) ** 2)
+    h = (
+        np.sqrt(x**2 + y**2) / np.cos(lat) - v
+        if lat < abs(1e-18)  # to avoid errors very close and at zero
+        else z / np.sin(lat) - (1 - e2) * v
+    )
 
-    return lat, lon, h
+    return lon, lat, h

src/poliastro/spheroid_location.py

@@ -141,19 +141,19 @@ def is_visible(self, px, py, pz):
 
     def cartesian_to_ellipsoidal(self, x, y, z):
         """
-        Converts ellipsoidal coordinates to the Cartesian coordinate system for the given ellipsoid.
+        Converts cartesian coordinates to ellipsoidal coordinates for this ellipsoid.
 
         Parameters
         ----------
         x : ~astropy.units.quantity.Quantity
-            x coordinate
+            x-coordinate
         y : ~astropy.units.quantity.Quantity
-            y coordinate
+            y-coordinate
         z : ~astropy.units.quantity.Quantity
-            z coordinate
+            z-coordinate
 
         """
         _a, _c = self._a.to_value(u.m), self._c.to_value(u.m)
         x, y, z = x.to_value(u.m), y.to_value(u.m), z.to_value(u.m)
-        lat, lon, h = cartesian_to_ellipsoidal_fast(_a, _c, x, y, z)
-        return lat * u.rad, lon * u.rad, h * u.m
+        lon, lat, h = cartesian_to_ellipsoidal_fast(_a, _c, x, y, z)
+        return lon * u.rad, lat * u.rad, h * u.m

src/poliastro/twobody/events.py

@@ -122,7 +122,7 @@ def __init__(self, orbit, lat, terminal=False, direction=0):
     def __call__(self, t, u_, k):
         self._last_t = t
         pos_on_body = (u_[:3] / norm(u_[:3])) * self._R
-        lat_, _, _ = cartesian_to_ellipsoidal_fast(
+        _, lat_, _ = cartesian_to_ellipsoidal_fast(
             self._R, self._R_polar, *pos_on_body
         )
 

tests/test_spheroid_location.py

@@ -14,9 +14,9 @@ def with_units(draw, elements, unit):
 
 def test_cartesian_coordinates():
     expected_cords = [
-        3764258.64785411 * u.m,
-        3295359.33856106 * u.m,
-        3942945.28570563 * u.m,
+        3764859.30127275 * u.m,
+        2987201.67496698 * u.m,
+        4179160.71540021 * u.m,
     ]
 
     el_cords = (38.43 * u.deg, 41.2 * u.deg, 0 * u.m)
@@ -78,7 +78,7 @@ def test_radius_of_curvature():
 
 
 def test_distance():
-    expected_distance = 6369864.745418392 * u.m
+    expected_distance = 6368850.150294118 * u.m
     el_cords = (38.43 * u.deg, 41.2 * u.deg, 0 * u.m)
     point_cords = (10.5 * u.m, 35.5 * u.m, 45.5 * u.m)
 
@@ -90,7 +90,7 @@ def test_distance():
 
 
 def test_cartesian_conversion_approximate():
-    el_cords = (0.670680 * u.rad, 0.7190227 * u.rad, 0 * u.m)
+    el_cords = (0.7190227 * u.rad, 0.670680 * u.rad, 0 * u.m)
 
     c_cords = [
         3764258.64785411 * u.m,
@@ -120,6 +120,6 @@ def test_h_calculation_near_lat_singularity(lat):
     h = 5 * u.m
     p = SpheroidLocation(lon, lat, h, body)
     cartesian_coords = p.cartesian_cords
-    lat_, lon_, h_ = p.cartesian_to_ellipsoidal(*cartesian_coords)
+    lon_, lat_, h_ = p.cartesian_to_ellipsoidal(*cartesian_coords)
 
     assert_quantity_allclose(h_, h)