Skip to content

FIX: correct default Earth radius (10x too large) in geopotential height conversions #1193

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
zjk199506 wants to merge 1 commit into from
+30 −12
Closed

FIX: correct default Earth radius (10x too large) in geopotential height conversions #1193

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 11 additions & 11 deletions rocketpy/tools.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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.

Expand All @@ -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
Expand All @@ -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.

Expand All @@ -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

Expand Down
20 changes: 19 additions & 1 deletion tests/unit/test_tools.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -22,7 +23,6 @@
tuple_handler,
)


WEB_MERCATOR_EARTH_RADIUS = 6378137.0


Expand Down Expand Up @@ -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)