Skip to content

Faster GEO rotation #233

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

Merged
jtec merged 4 commits into from
Jun 19, 2026
Merged

Faster GEO rotation #233

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e24bc88
compare old and new
jtec Jun 17, 2026
b94bc40
cleanup
jtec Jun 17, 2026
d57a8c8
fix
jtec Jun 17, 2026
2c2f36c
oups
jtec Jun 17, 2026
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
38 changes: 19 additions & 19 deletions src/prx/rinex_nav/evaluate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -315,25 +315,26 @@ 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
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
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):
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)],
]
)
Rz = np.array(
[
[np.cos(z_angle), np.sin(z_angle), 0],
Expand All @@ -356,16 +357,15 @@ def handle_bds_geos(eph):
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()
)

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


# Adapted from gnss_lib_py's find_sat()
Expand All @@ -392,7 +392,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)

@jtec jtec Jun 17, 2026

Copy link
Copy Markdown
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In most parts of the code we have functions return the modified dataframe.

eph = eph.rename(
columns={
"X_k": "sat_pos_x_m",
Expand Down
2 changes: 1 addition & 1 deletion src/prx/test/benchmark.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,

@jtec jtec Jun 17, 2026

Copy link
Copy Markdown
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Delivers results faster with enough resolution.

obs_file=obs_file,
root=obs_file.parent / "benchmark_datasets" if obs_file is not None else None,
)
Expand Down