From e24bc88b8d379dad40782942000f124c5e622509 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 21:14:06 +0200 Subject: [PATCH 1/4] compare old and new PASSED [100%]2026-06-17 21:13:31,815 [DEBUG] [prx.util]: Function handle_bds_geos took 0 days 00:00:24.540819 to run. 2026-06-17 21:13:34,538 [DEBUG] [prx.util]: Function handle_bds_geos_faster took 0 days 00:00:02.708934 to run. --- src/prx/rinex_nav/evaluate.py | 60 +++++++++++++++++++++++-- src/prx/rinex_nav/test/test_evaluate.py | 20 ++++++++- src/prx/test/benchmark.py | 2 +- 3 files changed, 75 insertions(+), 7 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index b221cc2..2c15fa1 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -314,7 +314,7 @@ def orbital_plane_to_earth_centered_cartesian(eph): eph["dZ_k"] = eph.y_k * eph.di_k * np.cos(eph.i_k) + eph.dy_k * np.sin(eph.i_k) pass - +@timeit def handle_bds_geos(eph): # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see # Beidou_ICD_B3I_v1.0, Table 5-11 @@ -323,10 +323,10 @@ def handle_bds_geos(eph): return P_GK = np.reshape(geos[["X_k", "Y_k", "Z_k"]].to_numpy(), (-1, 1)) V_GK = np.reshape(geos[["dX_k", "dY_k", "dZ_k"]].to_numpy(), (-1, 1)) - z_angles = geos.OmegaEarthIcd_rps * geos.t_k + z_angles = geos["OmegaEarthIcd_rps"] * geos["t_k"] rotation_matrices = [] + x_angle = util.deg_2_rad(-5.0) for i, z_angle in enumerate(z_angles): - x_angle = util.deg_2_rad(-5.0) Rx = np.array( [ [1, 0, 0], @@ -366,6 +366,58 @@ def frozen_to_rotating_bdcs(row): geos = geos.apply(frozen_to_rotating_bdcs, axis=1) eph[eph.is_bds_geo] = geos + return eph + + +@timeit +def handle_bds_geos_faster(eph): + # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see + # Beidou_ICD_B3I_v1.0, Table 5-11 + geos = eph[eph.is_bds_geo] + if geos.empty: + return + P_GK = np.reshape(geos[["X_k", "Y_k", "Z_k"]].to_numpy(), (-1, 1)) + V_GK = np.reshape(geos[["dX_k", "dY_k", "dZ_k"]].to_numpy(), (-1, 1)) + z_angles = geos["OmegaEarthIcd_rps"] * geos["t_k"] + rotation_matrices = [] + x_angle = util.deg_2_rad(-5.0) + for i, z_angle in enumerate(z_angles): + Rx = np.array( + [ + [1, 0, 0], + [0, np.cos(x_angle), np.sin(x_angle)], + [0, -np.sin(x_angle), np.cos(x_angle)], + ] + ) + Rz = np.array( + [ + [np.cos(z_angle), np.sin(z_angle), 0], + [-np.sin(z_angle), np.cos(z_angle), 0], + [0, 0, 1], + ] + ) + rotation_matrices.append(np.matmul(Rz, Rx)) + R = scipy.sparse.block_diag(rotation_matrices) + P_K = R @ P_GK + P_K = np.reshape(P_K, (-1, 3)) + geos["X_k"] = P_K[:, 0] + geos["Y_k"] = P_K[:, 1] + geos["Z_k"] = P_K[:, 2] + # Velocity in inertial frame that coincides with BDCS at this time, ie a "frozen" ECEF frame + V_K_frozen = R @ V_GK + V_K_frozen = np.reshape(V_K_frozen, (-1, 3)) + geos["dX_k"] = V_K_frozen[:, 0] + geos["dY_k"] = V_K_frozen[:, 1] + geos["dZ_k"] = V_K_frozen[:, 2] + + # Add term due to ECEFs angular velocity w.r.t. the frozen frame + # Leverage the fact that there are only BDS GEOs + assert geos["OmegaEarthIcd_rps"].nunique() == 1 + OmegaEarthIcd_rps = geos["OmegaEarthIcd_rps"].iloc[0] + geos[["dX_k", "dY_k", "dZ_k"]] += np.cross(np.array([0, 0, - OmegaEarthIcd_rps]), geos[["X_k", "Y_k", "Z_k"]].to_numpy()) + + eph[eph.is_bds_geo] = geos + return eph # Adapted from gnss_lib_py's find_sat() @@ -392,7 +444,7 @@ def kepler_orbit_position_and_velocity(eph): ) position_in_orbital_plane(eph) orbital_plane_to_earth_centered_cartesian(eph) - handle_bds_geos(eph) + eph = handle_bds_geos(eph) eph = eph.rename( columns={ "X_k": "sat_pos_x_m", diff --git a/src/prx/rinex_nav/test/test_evaluate.py b/src/prx/rinex_nav/test/test_evaluate.py index e45f7b3..8babc17 100644 --- a/src/prx/rinex_nav/test/test_evaluate.py +++ b/src/prx/rinex_nav/test/test_evaluate.py @@ -5,13 +5,13 @@ from prx.rinex_nav.evaluate import ( select_ephemerides, set_time_of_validity, - parse_rinex_nav_file, + parse_rinex_nav_file, handle_bds_geos, handle_bds_geos_faster, ) from prx.rinex_obs.parser import parse_rinex_obs_file from prx.precise_corrections.sp3 import evaluate as sp3_evaluate from prx.rinex_nav import evaluate as rinex_nav_evaluate from prx import constants, converters -from prx.util import week_and_seconds_2_timedelta +from prx.util import week_and_seconds_2_timedelta, configure_logging import shutil import pytest import itertools @@ -833,3 +833,19 @@ def test_compute_health_flag(input_for_test_2): assert (values == test[2]).all() print("done") + + +def test_benchmark_geo_rotation(): + # GIVEN the following - not physically meaningful - satellite positions and velocities + df = pd.DataFrame(np.random.rand(int(1e6), 6)) + df.columns = ["X_k", "Y_k", "Z_k", "dX_k", "dY_k", "dZ_k"] + # With roughly half of them belonging to BDS GEO satellites + df["is_bds_geo"] = df["X_k"] > 0.5 + df["t_k"] = 1.23 + df["OmegaEarthIcd_rps"] = constants.cBdsOmegaDotEarth_rps + configure_logging("DEBUG") + reference = handle_bds_geos(df.copy()) + candidate = handle_bds_geos_faster(df.copy()) + assert reference.equals(candidate) + #assert len(result_df) == len(df) + pass \ No newline at end of file diff --git a/src/prx/test/benchmark.py b/src/prx/test/benchmark.py index f5483e4..a50b7e0 100644 --- a/src/prx/test/benchmark.py +++ b/src/prx/test/benchmark.py @@ -86,7 +86,7 @@ def main(ram: bool, obs_file: Path, warm_parser_cache: bool): configure_logging("DEBUG") cases = generate_inputs( - n_steps=10, + n_steps=6, obs_file=obs_file, root=obs_file.parent / "benchmark_datasets" if obs_file is not None else None, ) From b94bc408d7989cddb5de4d932bc31413862372b8 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 21:19:47 +0200 Subject: [PATCH 2/4] cleanup --- src/prx/rinex_nav/evaluate.py | 5 ++++- src/prx/rinex_nav/test/test_evaluate.py | 20 ++------------------ 2 files changed, 6 insertions(+), 19 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 2c15fa1..9ad3a10 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -314,6 +314,7 @@ def orbital_plane_to_earth_centered_cartesian(eph): eph["dZ_k"] = eph.y_k * eph.di_k * np.cos(eph.i_k) + eph.dy_k * np.sin(eph.i_k) pass + @timeit def handle_bds_geos(eph): # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see @@ -414,7 +415,9 @@ def handle_bds_geos_faster(eph): # Leverage the fact that there are only BDS GEOs assert geos["OmegaEarthIcd_rps"].nunique() == 1 OmegaEarthIcd_rps = geos["OmegaEarthIcd_rps"].iloc[0] - geos[["dX_k", "dY_k", "dZ_k"]] += np.cross(np.array([0, 0, - OmegaEarthIcd_rps]), geos[["X_k", "Y_k", "Z_k"]].to_numpy()) + geos[["dX_k", "dY_k", "dZ_k"]] += np.cross( + np.array([0, 0, -OmegaEarthIcd_rps]), geos[["X_k", "Y_k", "Z_k"]].to_numpy() + ) eph[eph.is_bds_geo] = geos return eph diff --git a/src/prx/rinex_nav/test/test_evaluate.py b/src/prx/rinex_nav/test/test_evaluate.py index 8babc17..e45f7b3 100644 --- a/src/prx/rinex_nav/test/test_evaluate.py +++ b/src/prx/rinex_nav/test/test_evaluate.py @@ -5,13 +5,13 @@ from prx.rinex_nav.evaluate import ( select_ephemerides, set_time_of_validity, - parse_rinex_nav_file, handle_bds_geos, handle_bds_geos_faster, + parse_rinex_nav_file, ) from prx.rinex_obs.parser import parse_rinex_obs_file from prx.precise_corrections.sp3 import evaluate as sp3_evaluate from prx.rinex_nav import evaluate as rinex_nav_evaluate from prx import constants, converters -from prx.util import week_and_seconds_2_timedelta, configure_logging +from prx.util import week_and_seconds_2_timedelta import shutil import pytest import itertools @@ -833,19 +833,3 @@ def test_compute_health_flag(input_for_test_2): assert (values == test[2]).all() print("done") - - -def test_benchmark_geo_rotation(): - # GIVEN the following - not physically meaningful - satellite positions and velocities - df = pd.DataFrame(np.random.rand(int(1e6), 6)) - df.columns = ["X_k", "Y_k", "Z_k", "dX_k", "dY_k", "dZ_k"] - # With roughly half of them belonging to BDS GEO satellites - df["is_bds_geo"] = df["X_k"] > 0.5 - df["t_k"] = 1.23 - df["OmegaEarthIcd_rps"] = constants.cBdsOmegaDotEarth_rps - configure_logging("DEBUG") - reference = handle_bds_geos(df.copy()) - candidate = handle_bds_geos_faster(df.copy()) - assert reference.equals(candidate) - #assert len(result_df) == len(df) - pass \ No newline at end of file From d57a8c8f6aaad7a5aa0b2deb9af3d3bbe1a6b04c Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 21:24:33 +0200 Subject: [PATCH 3/4] fix --- src/prx/rinex_nav/evaluate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 9ad3a10..4a0f248 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -371,12 +371,12 @@ def frozen_to_rotating_bdcs(row): @timeit -def handle_bds_geos_faster(eph): +def handle_bds_geos(eph): # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see # Beidou_ICD_B3I_v1.0, Table 5-11 geos = eph[eph.is_bds_geo] if geos.empty: - return + return eph P_GK = np.reshape(geos[["X_k", "Y_k", "Z_k"]].to_numpy(), (-1, 1)) V_GK = np.reshape(geos[["dX_k", "dY_k", "dZ_k"]].to_numpy(), (-1, 1)) z_angles = geos["OmegaEarthIcd_rps"] * geos["t_k"] From 2c2f36c7fb912703cf10f87d410e33c1af73b9af Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 21:26:44 +0200 Subject: [PATCH 4/4] oups --- src/prx/rinex_nav/evaluate.py | 69 ++++------------------------------- 1 file changed, 7 insertions(+), 62 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 4a0f248..4034848 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -315,61 +315,6 @@ def orbital_plane_to_earth_centered_cartesian(eph): pass -@timeit -def handle_bds_geos(eph): - # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see - # Beidou_ICD_B3I_v1.0, Table 5-11 - geos = eph[eph.is_bds_geo] - if geos.empty: - return - P_GK = np.reshape(geos[["X_k", "Y_k", "Z_k"]].to_numpy(), (-1, 1)) - V_GK = np.reshape(geos[["dX_k", "dY_k", "dZ_k"]].to_numpy(), (-1, 1)) - z_angles = geos["OmegaEarthIcd_rps"] * geos["t_k"] - rotation_matrices = [] - x_angle = util.deg_2_rad(-5.0) - for i, z_angle in enumerate(z_angles): - Rx = np.array( - [ - [1, 0, 0], - [0, np.cos(x_angle), np.sin(x_angle)], - [0, -np.sin(x_angle), np.cos(x_angle)], - ] - ) - Rz = np.array( - [ - [np.cos(z_angle), np.sin(z_angle), 0], - [-np.sin(z_angle), np.cos(z_angle), 0], - [0, 0, 1], - ] - ) - rotation_matrices.append(np.matmul(Rz, Rx)) - R = scipy.sparse.block_diag(rotation_matrices) - P_K = R @ P_GK - P_K = np.reshape(P_K, (-1, 3)) - geos["X_k"] = P_K[:, 0] - geos["Y_k"] = P_K[:, 1] - geos["Z_k"] = P_K[:, 2] - # Velocity in inertial frame that coincides with BDCS at this time, ie a "frozen" ECEF frame - V_K_frozen = R @ V_GK - V_K_frozen = np.reshape(V_K_frozen, (-1, 3)) - geos["dX_k"] = V_K_frozen[:, 0] - geos["dY_k"] = V_K_frozen[:, 1] - geos["dZ_k"] = V_K_frozen[:, 2] - - # Add term due to ECEFs angular velocity w.r.t. the frozen frame - - def frozen_to_rotating_bdcs(row): - p = np.array([row["X_k"], row["Y_k"], row["Z_k"]]) - v_frozen = np.array([row["dX_k"], row["dY_k"], row["dZ_k"]]) - v_rotating = v_frozen + np.cross(np.array([0, 0, -row.OmegaEarthIcd_rps]), p) - row[["dX_k", "dY_k", "dZ_k"]] = v_rotating - return row - - geos = geos.apply(frozen_to_rotating_bdcs, axis=1) - eph[eph.is_bds_geo] = geos - return eph - - @timeit def handle_bds_geos(eph): # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see @@ -382,14 +327,14 @@ def handle_bds_geos(eph): z_angles = geos["OmegaEarthIcd_rps"] * geos["t_k"] rotation_matrices = [] x_angle = util.deg_2_rad(-5.0) + Rx = np.array( + [ + [1, 0, 0], + [0, np.cos(x_angle), np.sin(x_angle)], + [0, -np.sin(x_angle), np.cos(x_angle)], + ] + ) for i, z_angle in enumerate(z_angles): - Rx = np.array( - [ - [1, 0, 0], - [0, np.cos(x_angle), np.sin(x_angle)], - [0, -np.sin(x_angle), np.cos(x_angle)], - ] - ) Rz = np.array( [ [np.cos(z_angle), np.sin(z_angle), 0],