FIX: correct default Earth radius in geopotential height conversions

The default radius of geopotential_height_to_geometric_height,
geopotential_to_height_asl and geopotential_to_height_agl was
63781370 m - 10x the WGS-84 semi-major axis - contradicting their
own docstrings ("defaulting to 6378137.0"). The docstring examples
were anchored to the 10x-radius values, masking the bug.

The main Environment processing paths already pass earth_radius
explicitly and are unaffected; the default only bites when the
helpers are reused directly (e.g. EnvironmentAnalysis surface
elevation from NOAA reanalysis 'z'), biasing heights low by
~14 m at 10 km and ~57 m at 20 km.

Defaults now use 6378137.0 m; docstring examples and a regression
test pin the correct values.

Author: zhengjinkun

rocketpy/tools.py

@@ -969,7 +969,7 @@ def time_num_to_date_string(time_num, units, timezone, calendar="gregorian"):
     return date_string, hour_string, date_time
 
 
-def geopotential_height_to_geometric_height(geopotential_height, radius=63781370.0):
+def geopotential_height_to_geometric_height(geopotential_height, radius=6378137.0):
     """Converts geopotential height to geometric height.
 
     Parameters
@@ -992,14 +992,14 @@ def geopotential_height_to_geometric_height(geopotential_height, radius=63781370
     >>> geopotential_height_to_geometric_height(0)
     0.0
     >>> geopotential_height_to_geometric_height(10000)
-    10001.568101798659
+    10015.70317975257
     >>> geopotential_height_to_geometric_height(20000)
-    20006.2733909262
+    20062.91151008542
     """
     return radius * geopotential_height / (radius - geopotential_height)
 
 
-def geopotential_to_height_asl(geopotential, radius=63781370, g=9.80665):
+def geopotential_to_height_asl(geopotential, radius=6378137.0, g=9.80665):
     """Compute height above sea level from geopotential.
 
     Source: https://en.wikipedia.org/wiki/Geopotential
@@ -1010,7 +1010,7 @@ def geopotential_to_height_asl(geopotential, radius=63781370, g=9.80665):
         Geopotential in m^2/s^2. It is the geopotential value at a given
         pressure level, to be converted to height above sea level.
     radius : float, optional
-        Earth radius in m. Default is 63781370 m.
+        Earth radius in m. Default is 6378137.0 m (WGS-84 semi-major axis).
     g : float, optional
         Gravity acceleration in m/s^2. Default is 9.80665 m/s^2.
 
@@ -1025,15 +1025,15 @@ def geopotential_to_height_asl(geopotential, radius=63781370, g=9.80665):
     >>> geopotential_to_height_asl(0)
     0.0
     >>> geopotential_to_height_asl(100000)
-    10198.792680243916
+    10213.491133844715
     >>> geopotential_to_height_asl(200000)
-    20400.84750449947
+    20459.74503595946
     """
     geopotential_height = geopotential / g
     return geopotential_height_to_geometric_height(geopotential_height, radius)
 
 
-def geopotential_to_height_agl(geopotential, elevation, radius=63781370, g=9.80665):
+def geopotential_to_height_agl(geopotential, elevation, radius=6378137.0, g=9.80665):
     """Compute height above ground level from geopotential and elevation.
 
     Parameters
@@ -1044,7 +1044,7 @@ def geopotential_to_height_agl(geopotential, elevation, radius=63781370, g=9.806
     elevation : float
         Surface elevation in m
     radius : float, optional
-        Earth radius in m. Default is 63781370 m.
+        Earth radius in m. Default is 6378137.0 m (WGS-84 semi-major axis).
     g : float, optional
         Gravity acceleration in m/s^2. Default is 9.80665 m/s^2.
 
@@ -1059,9 +1059,9 @@ def geopotential_to_height_agl(geopotential, elevation, radius=63781370, g=9.806
     >>> geopotential_to_height_agl(0, 0)
     0.0
     >>> geopotential_to_height_agl(100000, 0)
-    10198.792680243916
+    10213.491133844715
     >>> geopotential_to_height_agl(100000, 1000)
-    9198.792680243916
+    9213.491133844715
     """
     return geopotential_to_height_asl(geopotential, radius, g) - elevation
 

tests/unit/test_tools.py

@@ -12,6 +12,7 @@
     euler313_to_quaternions,
     find_roots_cubic_function,
     generate_monte_carlo_ellipses,
+    geopotential_height_to_geometric_height,
     haversine,
     inverted_haversine,
     mercator_to_wgs84,
@@ -22,7 +23,6 @@
     tuple_handler,
 )
 
-
 WEB_MERCATOR_EARTH_RADIUS = 6378137.0
 
 
@@ -347,3 +347,21 @@ def test_mercator_extent_to_local_preserves_offset_sign(
     assert local_extent[0] < local_extent[1]
     assert local_extent[2] < local_extent[3]
     assert all(expected_sign * value > 0 for value in local_extent)
+
+
+def test_geopotential_height_to_geometric_height_default_radius():
+    """The default radius must be the WGS-84 semi-major axis (6378137.0 m).
+    A previous default of 63781370.0 (10x the Earth radius) biased the
+    conversion low by ~14 m at 10 km and ~57 m at 20 km.
+    """
+    assert geopotential_height_to_geometric_height(0) == 0.0
+    assert geopotential_height_to_geometric_height(10000) == pytest.approx(
+        10015.70317975257, rel=1e-10
+    )
+    assert geopotential_height_to_geometric_height(20000) == pytest.approx(
+        20062.91151008542, rel=1e-10
+    )
+    # Explicit radius still overrides the default
+    assert geopotential_height_to_geometric_height(
+        10000, radius=6378137.0
+    ) == pytest.approx(10015.70317975257, rel=1e-10)