jtec · plutonheaven · Jun 25, 2025 · May 6, 2025 · May 7, 2025 · May 7, 2025
diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py
@@ -392,7 +392,7 @@ def kepler_orbit_position_and_velocity(eph):
 
 def set_time_of_validity(df):
     def set_for_one_constellation(group):
-        group_constellation = group["constellation"].iloc[0]
+        group_constellation = group["constellation"].iat[0]
         group["validity_start"] = (
             group["ephemeris_reference_time_isagpst"]
             + constants.constellation_2_ephemeris_validity_interval[
@@ -407,7 +407,11 @@ def set_for_one_constellation(group):
         )
         return group
 
-    df = df.groupby("constellation").apply(set_for_one_constellation)
+    df = (
+        df.groupby("constellation")[df.columns]
+        .apply(set_for_one_constellation)
+        .reset_index(drop=True)
+    )
     return df
 
 
@@ -597,20 +601,27 @@ def compute_clock_offsets(df):
     return df
 
 
-def compute_parallel(rinex_nav_file_path, per_signal_query):
+def compute_parallel(
+    rinex_nav_file_path, per_signal_query, is_query_corrected_by_sat_clock_offset=False
+):
     # Warm up nav file parser cache so that we don't parse the file multiple times
     _ = parse_rinex_nav_file(rinex_nav_file_path)
     parallel = Parallel(n_jobs=round(multiprocessing.cpu_count() / 2), return_as="list")
     # split dataframe into `n_chunks` smaller dataframes
     n_chunks = min(len(per_signal_query.index), 4)
     chunks = np.array_split(per_signal_query, n_chunks)
     processed_chunks = parallel(
-        delayed(compute)(rinex_nav_file_path, chunk) for chunk in chunks
+        delayed(compute)(
+            rinex_nav_file_path, chunk, is_query_corrected_by_sat_clock_offset
+        )
+        for chunk in chunks
     )
     return pd.concat(processed_chunks)
 
 
-def compute(rinex_nav_file_path, per_signal_query):
+def compute(
+    rinex_nav_file_path, per_signal_query, is_query_corrected_by_sat_clock_offset=False
+):
     query_columns = per_signal_query.columns
     # per_signal_query is a pd.DataFrame with the following columns
     #   - time_of_reception_in_receiver_time
@@ -625,7 +636,22 @@ def compute(rinex_nav_file_path, per_signal_query):
     # Example: Galileo transmits E5a clock and group delay parameters in the F/NAV message, but parameters for other
     # signals in the I/NAV message
     per_signal_query = select_ephemerides(ephemerides, per_signal_query)
-    per_signal_query = compute_clock_offsets(per_signal_query)
+
+    # compute satellite clock bias
+    if is_query_corrected_by_sat_clock_offset:
+        per_signal_query = compute_clock_offsets(per_signal_query)
+    else:  # compute satellite clock offset iteratively
+        t = per_signal_query.query_time_wrt_clock_reference_time_s
+        for _ in range(2):
+            per_signal_query = compute_clock_offsets(per_signal_query)
+            per_signal_query.query_time_wrt_clock_reference_time_s = (
+                t - per_signal_query.sat_clock_offset_m / constants.cGpsSpeedOfLight_mps
+            )
+        # Apply sat clock correction to the query time for satellite position computation
+        per_signal_query.query_time_wrt_ephemeris_reference_time_s -= (
+            per_signal_query.sat_clock_offset_m / constants.cGpsSpeedOfLight_mps
+        )
+
     # Compute orbital states for each (satellite,ephemeris) pair only once:
     per_sat_eph_query = (
         per_signal_query.groupby(["sv", "query_time_isagpst", "ephemeris_hash"])
@@ -651,7 +677,9 @@ def evaluate_orbit(sub_df):
             sub_df[["x_m", "y_m", "z_m", "dx_mps", "dy_mps", "dz_mps"]] = np.nan
         return sub_df
 
-    per_sat_eph_query = per_sat_eph_query.groupby("orbit_type").apply(evaluate_orbit)
+    per_sat_eph_query = per_sat_eph_query.groupby("orbit_type")[
+        per_sat_eph_query.columns
+    ].apply(evaluate_orbit)
     per_sat_eph_query = per_sat_eph_query.reset_index(drop=True)
     columns_to_keep = [
         "sv",
@@ -765,7 +793,9 @@ def compute_tgds(df):
         df["sat_code_bias_m"] = df.tgd * df.gamma * df.speedOfLightIcd_mps
         return df
 
-    query = query.groupby(["signal", "constellation"]).apply(compute_tgds)
+    query = query.groupby(["signal", "constellation"])[query.columns].apply(
+        compute_tgds
+    )
     return query