From e24bc88b8d379dad40782942000f124c5e622509 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 21:14:06 +0200 Subject: [PATCH 1/9] 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/9] 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/9] 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/9] 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], From 78d13134b7ce7742bd4207e616118d41cef8e335 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 22:17:07 +0200 Subject: [PATCH 5/9] remove row-wise applies --- src/prx/precise_corrections/sp3/evaluate.py | 5 ++- src/prx/rinex_nav/evaluate.py | 37 +++++++++------------ src/prx/rinex_nav/test/test_evaluate.py | 3 +- 3 files changed, 19 insertions(+), 26 deletions(-) diff --git a/src/prx/precise_corrections/sp3/evaluate.py b/src/prx/precise_corrections/sp3/evaluate.py index bc8338b..2794de6 100644 --- a/src/prx/precise_corrections/sp3/evaluate.py +++ b/src/prx/precise_corrections/sp3/evaluate.py @@ -7,6 +7,7 @@ from scipy.interpolate import KroghInterpolator from prx import constants, util from prx.precise_corrections.antex import antex_processing as atx_processing +from prx.util import timedelta_2_seconds log = logging.getLogger(__name__) @@ -32,9 +33,7 @@ def cached_load(file_path: Path, file_hash: str): if "position" not in col and "velocity" not in col: df.rename(columns={col: col.replace("_x", "")}, inplace=True) # Convert timestamps to seconds since GPST epoch - df["gpst_s"] = (df["time"] - constants.cGpstUtcEpoch).apply( - util.timedelta_2_seconds - ) + df["gpst_s"] = timedelta_2_seconds(df["time"] - constants.cGpstUtcEpoch) df.drop(columns=["time"], inplace=True) df["sat_clock_offset_m"] = ( constants.cGpsSpeedOfLight_mps * df["clock"] diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 4034848..1487491 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -9,7 +9,7 @@ import georinex from prx import util from prx import constants -from prx.util import timeit, try_repair_with_gfzrnx +from prx.util import timeit, try_repair_with_gfzrnx, timedelta_2_seconds log = logging.getLogger(__name__) @@ -40,7 +40,9 @@ def cached_load(rinex_file_path: Path, file_hash: str): return cached_load(rinex_file, file_content_hash) -def time_scale_integer_second_offset_wrt_gpst(time_scale, utc_gpst_leap_seconds=None): +def time_scale_integer_second_offset_wrt_gpst( + time_scale: str, utc_gpst_leap_seconds=None +): if time_scale in ["GPST", "SBAST", "QZSST", "IRNSST", "GST"]: return pd.Timedelta(seconds=0) if time_scale == "BDT": @@ -552,22 +554,13 @@ def compute_gal_inav_fnav_indicators(df): return df -def to_isagpst(time, timescale, gpst_utc_leapseconds): - if (isinstance(time, pd.Timedelta) or isinstance(time, pd.Series)) and isinstance( - timescale, str - ): - return time - time_scale_integer_second_offset_wrt_gpst( - timescale, gpst_utc_leapseconds - ) - if isinstance(time, pd.Series) and isinstance(timescale, pd.Series): - return time - timescale.apply( - lambda element: time_scale_integer_second_offset_wrt_gpst( - element, gpst_utc_leapseconds - ) - ) - - assert False, ( - f"Unexpected types: time is {type(time)}, timescale is {type(timescale)}" +def to_isagpst( + time: pd.Timedelta | pd.Series, + timescale: str | pd.Series, + gpst_utc_leapseconds: int, +) -> pd.Timedelta | pd.Series: + return time - time_scale_integer_second_offset_wrt_gpst( + timescale, gpst_utc_leapseconds ) @@ -593,12 +586,12 @@ def select_ephemerides(df, query): direction="backward", ) # Compute times w.r.t. orbit and clock reference times used by downstream computations - query["query_time_wrt_ephemeris_reference_time_s"] = ( + query["query_time_wrt_ephemeris_reference_time_s"] = timedelta_2_seconds( query["query_time_isagpst"] - query["ephemeris_reference_time_isagpst"] - ).apply(util.timedelta_2_seconds) - query["query_time_wrt_clock_reference_time_s"] = ( + ) + query["query_time_wrt_clock_reference_time_s"] = timedelta_2_seconds( query["query_time_isagpst"] - query["clock_reference_time_isagpst"] - ).apply(util.timedelta_2_seconds) + ) query["ephemeris_valid"] = (query["query_time_isagpst"] < query["validity_end"]) & ( query["query_time_isagpst"] > query["validity_start"] ) diff --git a/src/prx/rinex_nav/test/test_evaluate.py b/src/prx/rinex_nav/test/test_evaluate.py index e45f7b3..9283af4 100644 --- a/src/prx/rinex_nav/test/test_evaluate.py +++ b/src/prx/rinex_nav/test/test_evaluate.py @@ -11,7 +11,7 @@ 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, disk_cache import shutil import pytest import itertools @@ -108,6 +108,7 @@ def input_for_test_2023(tmp_path_factory): def test_parse_nav_file(input_for_test): path_to_rnx3_nav_file = converters.anything_to_rinex_3(input_for_test["nav"]) + disk_cache.clear() df = parse_rinex_nav_file(path_to_rnx3_nav_file) assert not df.empty From 8956d60755061ccd52c95f4654c85339aad41b2a Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 22:24:58 +0200 Subject: [PATCH 6/9] test it --- src/prx/test/test_util.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 src/prx/test/test_util.py diff --git a/src/prx/test/test_util.py b/src/prx/test/test_util.py new file mode 100644 index 0000000..ebad492 --- /dev/null +++ b/src/prx/test/test_util.py @@ -0,0 +1,38 @@ +import logging +from datetime import timedelta + +import numpy as np + +import pandas as pd +import polars as pl + +from prx.constants import cSecondsPerDay +from prx.util import timedelta_2_seconds + +log = logging.getLogger(__name__) + + +def test_timedelta_2_seconds(): + expected_timedelta_s = cSecondsPerDay + 1.23456789 + assert np.isclose( + timedelta_2_seconds(pd.Timedelta(days=1, seconds=1.23456789)), + expected_timedelta_s, + atol=1e-9, + ) + assert np.isclose( + timedelta_2_seconds(pd.Series([pd.Timedelta(days=1, seconds=1.23456789)])).iloc[ + 0 + ], + expected_timedelta_s, + atol=1e-9, + ) + assert np.isclose( + timedelta_2_seconds( + pl.Series( + values=[timedelta(days=1, seconds=1.23456789)], + dtype=pl.Duration(time_unit="ns"), + ) + )[0], + expected_timedelta_s, + atol=1e-9, + ) From c65ef4aa73fd4150cb2f9933f7276c065162ee94 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 18 Jun 2026 23:02:20 +0200 Subject: [PATCH 7/9] oups --- src/prx/util.py | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/prx/util.py b/src/prx/util.py index 57df626..06f3b5e 100644 --- a/src/prx/util.py +++ b/src/prx/util.py @@ -9,7 +9,7 @@ from pathlib import Path import importlib.metadata as md import git - +import polars as pl import georinex import joblib import numpy as np @@ -233,21 +233,20 @@ def week_and_seconds_2_timedelta(weeks, seconds): return pd.Timedelta(weeks * constants.cSecondsPerWeek + seconds, "seconds") -def timedelta_2_seconds(time_delta: pd.Timedelta): - if pd.isnull(time_delta): - return np.nan - assert isinstance(time_delta, pd.Timedelta), ( - "time_delta must be of type pd.Timedelta" - ) - integer_seconds = np.float64(round(time_delta.total_seconds())) - fractional_seconds = ( - np.float64( - timedelta_2_nanoseconds(time_delta) - - integer_seconds * constants.cNanoSecondsPerSecond - ) - / constants.cNanoSecondsPerSecond +def timedelta_2_seconds( + time_delta: pd.Timedelta | pd.Series | pl.Series, +) -> float | pd.Series | pl.Series: + if isinstance(time_delta, pd.Timedelta): + return timedelta_2_seconds( + pl.Series([time_delta.value], dtype=pl.Duration(time_unit="ns")) + )[0] + if isinstance(time_delta, pd.Series): + return timedelta_2_seconds(pl.from_pandas(time_delta)).to_pandas() + assert isinstance(time_delta, pl.Series) + assert time_delta.dtype.time_unit == "ns" + return ( + time_delta.dt.total_nanoseconds().cast(float) / constants.cNanoSecondsPerSecond ) - return integer_seconds + fractional_seconds def timedelta_2_nanoseconds(time_delta: pd.Timedelta): From bbcd9f907eaadf7ec02d946525a9dd7696102506 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Fri, 19 Jun 2026 10:02:57 +0200 Subject: [PATCH 8/9] fix --- src/prx/rinex_nav/evaluate.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 1487491..11b4be8 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -41,7 +41,7 @@ def cached_load(rinex_file_path: Path, file_hash: str): def time_scale_integer_second_offset_wrt_gpst( - time_scale: str, utc_gpst_leap_seconds=None + time_scale: str, utc_gpst_leap_seconds: int = None ): if time_scale in ["GPST", "SBAST", "QZSST", "IRNSST", "GST"]: return pd.Timedelta(seconds=0) @@ -556,8 +556,8 @@ def compute_gal_inav_fnav_indicators(df): def to_isagpst( time: pd.Timedelta | pd.Series, - timescale: str | pd.Series, - gpst_utc_leapseconds: int, + timescale: str, + gpst_utc_leapseconds: int | None, ) -> pd.Timedelta | pd.Series: return time - time_scale_integer_second_offset_wrt_gpst( timescale, gpst_utc_leapseconds From 05606afe1c78a2c8dff6be7549fc3d1a72bac34c Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 25 Jun 2026 08:50:08 +0200 Subject: [PATCH 9/9] fix naming --- src/prx/test/test_helpers.py | 346 ---------------------------------- src/prx/test/test_util.py | 349 ++++++++++++++++++++++++++++++++++- 2 files changed, 343 insertions(+), 352 deletions(-) delete mode 100644 src/prx/test/test_helpers.py diff --git a/src/prx/test/test_helpers.py b/src/prx/test/test_helpers.py deleted file mode 100644 index 63d37bf..0000000 --- a/src/prx/test/test_helpers.py +++ /dev/null @@ -1,346 +0,0 @@ -import logging -import numpy as np -import pytest -from prx import constants, util -import pandas as pd -from pathlib import Path -import shutil -import os - -from prx.converters import anything_to_rinex_3 - -log = logging.getLogger(__name__) - - -@pytest.fixture -def input_for_test(): - test_directory = Path(f"./tmp_test_directory_{__name__}").resolve() - if test_directory.exists(): - # Make sure the expected file has not been generated before and is still on disk due to e.g. a previous - # test run having crashed: - shutil.rmtree(test_directory) - os.makedirs(test_directory) - compressed_compact_rinex_file = "TLSE00FRA_R_20230010100_10S_01S_MO.crx.gz" - test_file = {"obs": test_directory.joinpath(compressed_compact_rinex_file)} - shutil.copy( - Path(__file__).parent - / f"datasets/TLSE_2023001/{compressed_compact_rinex_file}", - test_file["obs"], - ) - assert test_file["obs"].exists() - - # Also provide ephemerides so the test does not have to download them: - ephemerides_file = "BRDC00IGS_R_20230010000_01D_MN.rnx.gz" - test_file["nav"] = test_directory.joinpath(ephemerides_file) - shutil.copy( - Path(__file__).parent / f"datasets/TLSE_2023001/{ephemerides_file}", - test_file["nav"].parent.joinpath(ephemerides_file), - ) - assert test_file["nav"].parent.joinpath(ephemerides_file).exists() - - # sp3 file - sp3_file = "GFZ0MGXRAP_20230010000_01D_05M_ORB.SP3" - test_file["sp3"] = test_directory.joinpath(sp3_file) - shutil.copy( - Path(__file__).parent / f"datasets/TLSE_2023001/{sp3_file}", - test_file["sp3"].parent.joinpath(sp3_file), - ) - assert test_file["sp3"].parent.joinpath(sp3_file).exists() - - yield test_file - shutil.rmtree(test_directory) - - -coords = { - "equator_0": {"ecef": [6378137.0, 0.0, 0.0], "geodetic": [0.0, 0.0, 0.0]}, - "equator_90": { - "ecef": [0.0, 6378137.0, 0.0], - "geodetic": [np.deg2rad(0.0), np.deg2rad(90), 0.0], - }, - "toulouse": { - "ecef": [4624518, 116590, 4376497], - "geodetic": [np.deg2rad(43.604698100243), np.deg2rad(1.444193786348), 151.9032], - }, - "Sidney": { - "ecef": [-4.646050004314417e06, 2.553206120634516e06, -3.534374202256767e06], - "geodetic": [np.deg2rad(-33.8688197), np.deg2rad(151.2092955), 0], - }, - "Ushuaia": { - "ecef": [1.362205559782862e06, -3.423584689115747e06, -5.188704112366104e06], - "geodetic": [np.deg2rad(-54.8019121), np.deg2rad(-68.3029511), 0], - }, -} - - -def test_rinex_header_time_string_2_timestamp_ns(): - assert util.timestamp_2_timedelta( - util.rinex_header_time_string_2_timestamp_ns( - " 1980 1 6 0 0 0.0000000 GPS" - ), - "GPST", - ) == pd.Timedelta(0) - assert util.timestamp_2_timedelta( - util.rinex_header_time_string_2_timestamp_ns( - " 1980 1 6 0 0 1.0000000 GPS" - ), - "GPST", - ) == pd.Timedelta(constants.cNanoSecondsPerSecond, unit="ns") - timestamp = util.rinex_header_time_string_2_timestamp_ns( - " 1980 1 6 0 0 1.0000001 GPS" - ) - timedelta = util.timestamp_2_timedelta(timestamp, "GPST") - - assert timedelta == pd.Timedelta(constants.cNanoSecondsPerSecond + 100, unit="ns") - assert util.timestamp_2_timedelta( - util.rinex_header_time_string_2_timestamp_ns( - " 1980 1 7 0 0 0.0000000 GPS" - ), - "GPST", - ) == pd.Timedelta( - constants.cSecondsPerDay * constants.cNanoSecondsPerSecond, unit="ns" - ) - - -def test_ecef_to_geodetic(): - tolerance_rad = 1e-3 / 6400e3 # equivalent to one mm at Earth surface - tolerance_alt = 1e-3 - - for coord in coords.values(): - computed_geodetic = util.ecef_2_geodetic(coord["ecef"]) - assert ( - np.abs(np.array(coord["geodetic"][:2]) - np.array(computed_geodetic[:2])) - < tolerance_rad - ).all() - assert np.abs(coord["geodetic"][2] - computed_geodetic[2]) < tolerance_alt - - -def test_geodetic_2_ecef(): - tolerance_ecef = 1e-3 - - for coord in coords.values(): - computed_ecef = util.geodetic_2_ecef( - coord["geodetic"][0], coord["geodetic"][1], coord["geodetic"][2] - ) - assert ( - np.linalg.norm(np.array(coord["ecef"]) - np.array(computed_ecef)) - < tolerance_ecef - ) - - -def test_satellite_elevation_and_azimuth(): - tolerance = np.deg2rad(1e-3) - - sat_pos_ecef = np.array([[26600e3, 0.0, 0.0]]) - rx_pos_ecef = np.array([6400e3, 0.0, 0.0]) - expected_el, expected_az = np.deg2rad(90), np.deg2rad(0) - computed_el, computed_az = util.compute_satellite_elevation_and_azimuth( - sat_pos_ecef, rx_pos_ecef - ) - assert np.abs(expected_el - computed_el) < tolerance - assert np.abs(expected_az - computed_az) < tolerance - - sat_pos_ecef = np.array([[2.066169397996826e07, 0.0, 1.428355697996826e07]]) - rx_pos_ecef = np.array([6378137.0, 0.0, 0.0]) - expected_el, expected_az = np.deg2rad(45), np.deg2rad(0) - computed_el, computed_az = util.compute_satellite_elevation_and_azimuth( - sat_pos_ecef, rx_pos_ecef - ) - assert np.abs(expected_el - computed_el) < tolerance - assert np.abs(expected_az - computed_az) < tolerance - - sat_pos_ecef = np.array( - [[2.066169397996826e07, 7.141778489984130e06, 1.236992320105505e07]] - ) - rx_pos_ecef = np.array([6378137.0, 0.0, 0.0]) - expected_el, expected_az = np.deg2rad(45), np.deg2rad(30) - computed_el, computed_az = util.compute_satellite_elevation_and_azimuth( - sat_pos_ecef, rx_pos_ecef - ) - assert np.abs(expected_el - computed_el) < tolerance - assert np.abs(expected_az - computed_az) < tolerance - - -def test_sagnac_effect(): - # load validation data - path_to_validation_file = ( - util.prx_src_directory() / "tools/validation_data/sagnac_effect.csv" - ) - - # satellite position (from reference CSV header) - sat_pos = np.array([[28400000, 0, 0]]) - - # read data - data = np.loadtxt( - path_to_validation_file, - delimiter=",", - skiprows=3, - ) - - sagnac_effect_reference = np.zeros((data.shape[0],)) - sagnac_effect_computed = np.zeros((data.shape[0],)) - for ind in range(data.shape[0]): - rx_pos = data[ind, 0:3] - sagnac_effect_reference[ind] = data[ind, 3] - sagnac_effect_computed[ind] = util.compute_sagnac_effect(sat_pos, rx_pos) - - # errors come from the approximation of cos and sin for small angles - # millimeter accuracy should be sufficient - tolerance = 1e-3 - assert np.max(np.abs(sagnac_effect_computed - sagnac_effect_reference)) < tolerance - - -def test_is_sorted(): - assert util.is_sorted([1, 2, 3, 4, 5]) - assert not util.is_sorted([1, 2, 3, 5, 4]) - assert not util.is_sorted([5, 4, 3, 2, 1]) - assert util.is_sorted([1, 1, 1, 1, 1]) - assert util.is_sorted([1]) - assert util.is_sorted([]) - - -def test_row_wise_dot_product(): - # Check whether the way we compute the row-wise dot product with numpy yields the expected result - A = np.array([[1, 2], [4, 5], [7, 8]]) - B = np.array([[10, 20], [30, 40], [50, 60]]) - row_wise_dot = np.sum(A * B, axis=1).reshape(-1, 1) - assert (row_wise_dot == np.array([[10 + 40], [120 + 200], [350 + 480]])).all() - - -def test_timedelta_2_weeks_and_seconds(): - tested_timestamps = [ - pd.Timestamp("1981-01-21T06:00:00"), - pd.Timestamp("1981-01-21T06:10:00"), - pd.Timestamp("1999-09-01T16:10:00"), # after first week number rollover - pd.Timestamp("2019-04-10T12:00:00"), # after 2nd week number rollover - pd.Timestamp(""), # yields a NaT - ] - - week_computed = [] - seconds_of_week_computed = [] - for timestamp in tested_timestamps: - tested_timedelta = ( - timestamp - constants.system_time_scale_rinex_utc_epoch["GPST"] - ) - w, s = util.timedelta_2_weeks_and_seconds(tested_timedelta) - week_computed.append(w) - seconds_of_week_computed.append(s) - - # expected values come from https://gnsscalc.com/ - week_expected = [54, 54, 1025, 2048, np.nan] - seconds_of_week_expected = [280800, 281400, 317400, 302400, np.nan] - - np.testing.assert_array_equal(week_computed, week_expected) - np.testing.assert_array_equal(seconds_of_week_computed, seconds_of_week_expected) - - # We also expect the function to work for Series of timestamps - week_series_computed, seconds_of_week_series_computed = ( - util.timedelta_2_weeks_and_seconds( - pd.Series(tested_timestamps) - - constants.system_time_scale_rinex_utc_epoch["GPST"] - ) - ) - np.testing.assert_array_equal(week_series_computed, week_expected) - np.testing.assert_array_equal( - seconds_of_week_series_computed, seconds_of_week_series_computed - ) - - -def test_compute_gps_leap_seconds(): - # expected GPS leap second come from https://gnsscalc.com/ - test_cases = [ - (1982, 181, 2), # year, doy (corresponds to 01-Jul), expected leap seconds - (1989, 181, 5), - (1992, 182, 8), # leap year, still corresponds to 01-Jul - (2008, 182, 14), - (2016, 182, 17), - (2022, 181, 18), - (1979, 1, np.nan), - ] - for year, doy, expected_leap_second in test_cases: - np.testing.assert_equal( - util.compute_gps_utc_leap_seconds(year, doy), expected_leap_second - ) - - -def test_timestamp_to_mid_day(): - assert util.timestamp_to_mid_day( - pd.Timestamp("2023-01-01T01:02:03") - ) == pd.Timestamp("2023-01-01T12:00:00") - - -def test_sun_pos(): - pos_sun = ( - util.compute_sun_ecef_position( - np.array( - [ - pd.Timestamp( - year=2020, month=6, day=20, hour=6 - ), # approx time when the sun is "close" to [0,1,xx] AU (Equinox)pd.Timestamp( - pd.Timestamp( - year=2020, month=6, day=20, hour=12 - ), # approx time when the sun is "close" to [1,0,xx] AU (Equinox) - pd.Timestamp( - year=2020, month=6, day=20, hour=18 - ), # approx time when the sun is "close" to [0,-1,xx] AU - ] - ) - ) - / np.array([149_597_870_700] * 3) - ) - assert pos_sun[0:2, 0] == pytest.approx(np.array([0, 1]), abs=0.1) - assert pos_sun[0:2, 1] == pytest.approx(np.array([1, 0]), abs=0.1) - assert pos_sun[0:2, 2] == pytest.approx(np.array([0, -1]), abs=0.1) - - -def test_sat_frame(): - """ - Some tests to validate the conversion between ECEF and satellite-fixed frame. - - In order to control the sun position, the function util.ecef_2_satellite is not used, but its content is - replicated here. - """ - pos_sat_ecef = np.array([26_600_000, 0, 0]).reshape(1, -1) - pos_ecef = np.array([6_400_000, 0, 0]).reshape(1, -1) - - k = -pos_sat_ecef / np.linalg.norm(pos_sat_ecef, axis=1).reshape(1, -1) - pos_sun_ecef = np.array([0, 149_597_870_700, 0]).reshape(1, -1) # 1 AU on y - unit_vector_sun_ecef = (pos_sun_ecef - pos_sat_ecef) / np.linalg.norm( - pos_sun_ecef - pos_sat_ecef, axis=1 - ).reshape(-1, 1) - j = np.cross(k, unit_vector_sun_ecef) - i = np.cross(j, k) - - rot_mat_ecef2sat = np.stack([i, j, k], axis=1) - pos_sat_frame = np.stack( - [ - rot_mat_ecef2sat[i, :, :] @ (pos_ecef[i, :] - pos_sat_ecef[i, :]) - for i in range(pos_sat_ecef.shape[0]) - ] - ) - assert pos_sat_frame[0, :] == pytest.approx(np.array([0, 0, 20_200_000])) - - # other simple examples - pos_sun = np.array([0, 149_597_870_700, 0]) # 1 AU on y - pos_rx = np.array([6_400_000, 0, 0]) # on +x - pos_sat = np.array([26_600_000, 0, 0]) # on +x - - k = -pos_sat / np.linalg.norm(pos_sat) - unit_vector_sun_ecef = (pos_sun - pos_sat) / np.linalg.norm(pos_sun - pos_sat) - j = np.cross(k, unit_vector_sun_ecef) - i = np.cross(j, k) - - pos_rx_sat = np.stack([i, j, k], axis=0) @ (pos_rx - pos_sat) - assert (pos_rx_sat == np.array([0, 0, 20_200_000])).all() - - pos_rx = np.array([0, 6_400_000, 0]) # on +y - pos_rx_sat = np.stack([i, j, k], axis=0) @ (pos_rx - pos_sat) - assert pos_rx_sat == pytest.approx(np.array([6_400_000, 0, 20_200_000 + 6_400_000])) - - -def test_rinex_type_based_on_first_line(input_for_test): - obs = anything_to_rinex_3(input_for_test["obs"]) - nav = anything_to_rinex_3(input_for_test["nav"]) - assert util.is_rinex_3_obs_file(obs) - assert util.is_rinex_3_nav_file(nav) - assert not util.is_rinex_3_obs_file(nav) - assert not util.is_rinex_3_nav_file(obs) diff --git a/src/prx/test/test_util.py b/src/prx/test/test_util.py index ebad492..0acddca 100644 --- a/src/prx/test/test_util.py +++ b/src/prx/test/test_util.py @@ -1,17 +1,354 @@ import logging -from datetime import timedelta - import numpy as np - +import pytest +from prx import constants, util import pandas as pd +from pathlib import Path +import shutil +import os import polars as pl - +import datetime from prx.constants import cSecondsPerDay +from prx.converters import anything_to_rinex_3 from prx.util import timedelta_2_seconds log = logging.getLogger(__name__) +@pytest.fixture +def input_for_test(): + test_directory = Path(f"./tmp_test_directory_{__name__}").resolve() + if test_directory.exists(): + # Make sure the expected file has not been generated before and is still on disk due to e.g. a previous + # test run having crashed: + shutil.rmtree(test_directory) + os.makedirs(test_directory) + compressed_compact_rinex_file = "TLSE00FRA_R_20230010100_10S_01S_MO.crx.gz" + test_file = {"obs": test_directory.joinpath(compressed_compact_rinex_file)} + shutil.copy( + Path(__file__).parent + / f"datasets/TLSE_2023001/{compressed_compact_rinex_file}", + test_file["obs"], + ) + assert test_file["obs"].exists() + + # Also provide ephemerides so the test does not have to download them: + ephemerides_file = "BRDC00IGS_R_20230010000_01D_MN.rnx.gz" + test_file["nav"] = test_directory.joinpath(ephemerides_file) + shutil.copy( + Path(__file__).parent / f"datasets/TLSE_2023001/{ephemerides_file}", + test_file["nav"].parent.joinpath(ephemerides_file), + ) + assert test_file["nav"].parent.joinpath(ephemerides_file).exists() + + # sp3 file + sp3_file = "GFZ0MGXRAP_20230010000_01D_05M_ORB.SP3" + test_file["sp3"] = test_directory.joinpath(sp3_file) + shutil.copy( + Path(__file__).parent / f"datasets/TLSE_2023001/{sp3_file}", + test_file["sp3"].parent.joinpath(sp3_file), + ) + assert test_file["sp3"].parent.joinpath(sp3_file).exists() + + yield test_file + shutil.rmtree(test_directory) + + +coords = { + "equator_0": {"ecef": [6378137.0, 0.0, 0.0], "geodetic": [0.0, 0.0, 0.0]}, + "equator_90": { + "ecef": [0.0, 6378137.0, 0.0], + "geodetic": [np.deg2rad(0.0), np.deg2rad(90), 0.0], + }, + "toulouse": { + "ecef": [4624518, 116590, 4376497], + "geodetic": [np.deg2rad(43.604698100243), np.deg2rad(1.444193786348), 151.9032], + }, + "Sidney": { + "ecef": [-4.646050004314417e06, 2.553206120634516e06, -3.534374202256767e06], + "geodetic": [np.deg2rad(-33.8688197), np.deg2rad(151.2092955), 0], + }, + "Ushuaia": { + "ecef": [1.362205559782862e06, -3.423584689115747e06, -5.188704112366104e06], + "geodetic": [np.deg2rad(-54.8019121), np.deg2rad(-68.3029511), 0], + }, +} + + +def test_rinex_header_time_string_2_timestamp_ns(): + assert util.timestamp_2_timedelta( + util.rinex_header_time_string_2_timestamp_ns( + " 1980 1 6 0 0 0.0000000 GPS" + ), + "GPST", + ) == pd.Timedelta(0) + assert util.timestamp_2_timedelta( + util.rinex_header_time_string_2_timestamp_ns( + " 1980 1 6 0 0 1.0000000 GPS" + ), + "GPST", + ) == pd.Timedelta(constants.cNanoSecondsPerSecond, unit="ns") + timestamp = util.rinex_header_time_string_2_timestamp_ns( + " 1980 1 6 0 0 1.0000001 GPS" + ) + timedelta = util.timestamp_2_timedelta(timestamp, "GPST") + + assert timedelta == pd.Timedelta(constants.cNanoSecondsPerSecond + 100, unit="ns") + assert util.timestamp_2_timedelta( + util.rinex_header_time_string_2_timestamp_ns( + " 1980 1 7 0 0 0.0000000 GPS" + ), + "GPST", + ) == pd.Timedelta( + constants.cSecondsPerDay * constants.cNanoSecondsPerSecond, unit="ns" + ) + + +def test_ecef_to_geodetic(): + tolerance_rad = 1e-3 / 6400e3 # equivalent to one mm at Earth surface + tolerance_alt = 1e-3 + + for coord in coords.values(): + computed_geodetic = util.ecef_2_geodetic(coord["ecef"]) + assert ( + np.abs(np.array(coord["geodetic"][:2]) - np.array(computed_geodetic[:2])) + < tolerance_rad + ).all() + assert np.abs(coord["geodetic"][2] - computed_geodetic[2]) < tolerance_alt + + +def test_geodetic_2_ecef(): + tolerance_ecef = 1e-3 + + for coord in coords.values(): + computed_ecef = util.geodetic_2_ecef( + coord["geodetic"][0], coord["geodetic"][1], coord["geodetic"][2] + ) + assert ( + np.linalg.norm(np.array(coord["ecef"]) - np.array(computed_ecef)) + < tolerance_ecef + ) + + +def test_satellite_elevation_and_azimuth(): + tolerance = np.deg2rad(1e-3) + + sat_pos_ecef = np.array([[26600e3, 0.0, 0.0]]) + rx_pos_ecef = np.array([6400e3, 0.0, 0.0]) + expected_el, expected_az = np.deg2rad(90), np.deg2rad(0) + computed_el, computed_az = util.compute_satellite_elevation_and_azimuth( + sat_pos_ecef, rx_pos_ecef + ) + assert np.abs(expected_el - computed_el) < tolerance + assert np.abs(expected_az - computed_az) < tolerance + + sat_pos_ecef = np.array([[2.066169397996826e07, 0.0, 1.428355697996826e07]]) + rx_pos_ecef = np.array([6378137.0, 0.0, 0.0]) + expected_el, expected_az = np.deg2rad(45), np.deg2rad(0) + computed_el, computed_az = util.compute_satellite_elevation_and_azimuth( + sat_pos_ecef, rx_pos_ecef + ) + assert np.abs(expected_el - computed_el) < tolerance + assert np.abs(expected_az - computed_az) < tolerance + + sat_pos_ecef = np.array( + [[2.066169397996826e07, 7.141778489984130e06, 1.236992320105505e07]] + ) + rx_pos_ecef = np.array([6378137.0, 0.0, 0.0]) + expected_el, expected_az = np.deg2rad(45), np.deg2rad(30) + computed_el, computed_az = util.compute_satellite_elevation_and_azimuth( + sat_pos_ecef, rx_pos_ecef + ) + assert np.abs(expected_el - computed_el) < tolerance + assert np.abs(expected_az - computed_az) < tolerance + + +def test_sagnac_effect(): + # load validation data + path_to_validation_file = ( + util.prx_src_directory() / "tools/validation_data/sagnac_effect.csv" + ) + + # satellite position (from reference CSV header) + sat_pos = np.array([[28400000, 0, 0]]) + + # read data + data = np.loadtxt( + path_to_validation_file, + delimiter=",", + skiprows=3, + ) + + sagnac_effect_reference = np.zeros((data.shape[0],)) + sagnac_effect_computed = np.zeros((data.shape[0],)) + for ind in range(data.shape[0]): + rx_pos = data[ind, 0:3] + sagnac_effect_reference[ind] = data[ind, 3] + sagnac_effect_computed[ind] = util.compute_sagnac_effect(sat_pos, rx_pos) + + # errors come from the approximation of cos and sin for small angles + # millimeter accuracy should be sufficient + tolerance = 1e-3 + assert np.max(np.abs(sagnac_effect_computed - sagnac_effect_reference)) < tolerance + + +def test_is_sorted(): + assert util.is_sorted([1, 2, 3, 4, 5]) + assert not util.is_sorted([1, 2, 3, 5, 4]) + assert not util.is_sorted([5, 4, 3, 2, 1]) + assert util.is_sorted([1, 1, 1, 1, 1]) + assert util.is_sorted([1]) + assert util.is_sorted([]) + + +def test_row_wise_dot_product(): + # Check whether the way we compute the row-wise dot product with numpy yields the expected result + A = np.array([[1, 2], [4, 5], [7, 8]]) + B = np.array([[10, 20], [30, 40], [50, 60]]) + row_wise_dot = np.sum(A * B, axis=1).reshape(-1, 1) + assert (row_wise_dot == np.array([[10 + 40], [120 + 200], [350 + 480]])).all() + + +def test_timedelta_2_weeks_and_seconds(): + tested_timestamps = [ + pd.Timestamp("1981-01-21T06:00:00"), + pd.Timestamp("1981-01-21T06:10:00"), + pd.Timestamp("1999-09-01T16:10:00"), # after first week number rollover + pd.Timestamp("2019-04-10T12:00:00"), # after 2nd week number rollover + pd.Timestamp(""), # yields a NaT + ] + + week_computed = [] + seconds_of_week_computed = [] + for timestamp in tested_timestamps: + tested_timedelta = ( + timestamp - constants.system_time_scale_rinex_utc_epoch["GPST"] + ) + w, s = util.timedelta_2_weeks_and_seconds(tested_timedelta) + week_computed.append(w) + seconds_of_week_computed.append(s) + + # expected values come from https://gnsscalc.com/ + week_expected = [54, 54, 1025, 2048, np.nan] + seconds_of_week_expected = [280800, 281400, 317400, 302400, np.nan] + + np.testing.assert_array_equal(week_computed, week_expected) + np.testing.assert_array_equal(seconds_of_week_computed, seconds_of_week_expected) + + # We also expect the function to work for Series of timestamps + week_series_computed, seconds_of_week_series_computed = ( + util.timedelta_2_weeks_and_seconds( + pd.Series(tested_timestamps) + - constants.system_time_scale_rinex_utc_epoch["GPST"] + ) + ) + np.testing.assert_array_equal(week_series_computed, week_expected) + np.testing.assert_array_equal( + seconds_of_week_series_computed, seconds_of_week_series_computed + ) + + +def test_compute_gps_leap_seconds(): + # expected GPS leap second come from https://gnsscalc.com/ + test_cases = [ + (1982, 181, 2), # year, doy (corresponds to 01-Jul), expected leap seconds + (1989, 181, 5), + (1992, 182, 8), # leap year, still corresponds to 01-Jul + (2008, 182, 14), + (2016, 182, 17), + (2022, 181, 18), + (1979, 1, np.nan), + ] + for year, doy, expected_leap_second in test_cases: + np.testing.assert_equal( + util.compute_gps_utc_leap_seconds(year, doy), expected_leap_second + ) + + +def test_timestamp_to_mid_day(): + assert util.timestamp_to_mid_day( + pd.Timestamp("2023-01-01T01:02:03") + ) == pd.Timestamp("2023-01-01T12:00:00") + + +def test_sun_pos(): + pos_sun = ( + util.compute_sun_ecef_position( + np.array( + [ + pd.Timestamp( + year=2020, month=6, day=20, hour=6 + ), # approx time when the sun is "close" to [0,1,xx] AU (Equinox)pd.Timestamp( + pd.Timestamp( + year=2020, month=6, day=20, hour=12 + ), # approx time when the sun is "close" to [1,0,xx] AU (Equinox) + pd.Timestamp( + year=2020, month=6, day=20, hour=18 + ), # approx time when the sun is "close" to [0,-1,xx] AU + ] + ) + ) + / np.array([149_597_870_700] * 3) + ) + assert pos_sun[0:2, 0] == pytest.approx(np.array([0, 1]), abs=0.1) + assert pos_sun[0:2, 1] == pytest.approx(np.array([1, 0]), abs=0.1) + assert pos_sun[0:2, 2] == pytest.approx(np.array([0, -1]), abs=0.1) + + +def test_sat_frame(): + """ + Some tests to validate the conversion between ECEF and satellite-fixed frame. + + In order to control the sun position, the function util.ecef_2_satellite is not used, but its content is + replicated here. + """ + pos_sat_ecef = np.array([26_600_000, 0, 0]).reshape(1, -1) + pos_ecef = np.array([6_400_000, 0, 0]).reshape(1, -1) + + k = -pos_sat_ecef / np.linalg.norm(pos_sat_ecef, axis=1).reshape(1, -1) + pos_sun_ecef = np.array([0, 149_597_870_700, 0]).reshape(1, -1) # 1 AU on y + unit_vector_sun_ecef = (pos_sun_ecef - pos_sat_ecef) / np.linalg.norm( + pos_sun_ecef - pos_sat_ecef, axis=1 + ).reshape(-1, 1) + j = np.cross(k, unit_vector_sun_ecef) + i = np.cross(j, k) + + rot_mat_ecef2sat = np.stack([i, j, k], axis=1) + pos_sat_frame = np.stack( + [ + rot_mat_ecef2sat[i, :, :] @ (pos_ecef[i, :] - pos_sat_ecef[i, :]) + for i in range(pos_sat_ecef.shape[0]) + ] + ) + assert pos_sat_frame[0, :] == pytest.approx(np.array([0, 0, 20_200_000])) + + # other simple examples + pos_sun = np.array([0, 149_597_870_700, 0]) # 1 AU on y + pos_rx = np.array([6_400_000, 0, 0]) # on +x + pos_sat = np.array([26_600_000, 0, 0]) # on +x + + k = -pos_sat / np.linalg.norm(pos_sat) + unit_vector_sun_ecef = (pos_sun - pos_sat) / np.linalg.norm(pos_sun - pos_sat) + j = np.cross(k, unit_vector_sun_ecef) + i = np.cross(j, k) + + pos_rx_sat = np.stack([i, j, k], axis=0) @ (pos_rx - pos_sat) + assert (pos_rx_sat == np.array([0, 0, 20_200_000])).all() + + pos_rx = np.array([0, 6_400_000, 0]) # on +y + pos_rx_sat = np.stack([i, j, k], axis=0) @ (pos_rx - pos_sat) + assert pos_rx_sat == pytest.approx(np.array([6_400_000, 0, 20_200_000 + 6_400_000])) + + +def test_rinex_type_based_on_first_line(input_for_test): + obs = anything_to_rinex_3(input_for_test["obs"]) + nav = anything_to_rinex_3(input_for_test["nav"]) + assert util.is_rinex_3_obs_file(obs) + assert util.is_rinex_3_nav_file(nav) + assert not util.is_rinex_3_obs_file(nav) + assert not util.is_rinex_3_nav_file(obs) + + def test_timedelta_2_seconds(): expected_timedelta_s = cSecondsPerDay + 1.23456789 assert np.isclose( @@ -29,10 +366,10 @@ def test_timedelta_2_seconds(): assert np.isclose( timedelta_2_seconds( pl.Series( - values=[timedelta(days=1, seconds=1.23456789)], + values=[datetime.timedelta(days=1, seconds=1.23456789)], dtype=pl.Duration(time_unit="ns"), ) )[0], expected_timedelta_s, - atol=1e-9, + atol=1e-12, )