From 475613e8b4b39c6f4bb1e1ac5284a374e9e274fa Mon Sep 17 00:00:00 2001
From: "Christopher H. Jordan" <christopherjordan87@gmail.com>
Date: Mon, 6 Feb 2023 16:08:43 +0800
Subject: [PATCH] Fix raw data timestamp bug.

For some obsids (large ones?), treating an mwalib GPS time as a float, adding
half an integration time and dividing by 1000 to get units of seconds was enough
to introduce a float precision issue. hifitime allows us to use its types from
nanoseconds, so we keep things as integers and fix these issues for good. This
commit also changes a couple of other potentially-problematic bits of code.
---
 CHANGELOG.md                      |  2 ++
 src/cli/vis_utils/simulate/mod.rs |  6 +++++-
 src/io/read/raw/mod.rs            | 35 ++++++++++++++++++++++---------
 3 files changed, 32 insertions(+), 11 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index bacbbf56..2c18ef24 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -24,6 +24,8 @@ Versioning](https://semver.org/spec/v2.0.0.html).
   now read properly.
 - RTS source lists with multiple SHAPELET components are now removed properly
   (not to be confused with the favoured SHAPELET2 component type).
+- Some MWA raw data observations were not handled correctly due to a
+  floating-point error.
 
 ### Changed
 - The performance of CPU visibility modelling has been dramatically improved.
diff --git a/src/cli/vis_utils/simulate/mod.rs b/src/cli/vis_utils/simulate/mod.rs
index a962e392..b7cbb3ad 100644
--- a/src/cli/vis_utils/simulate/mod.rs
+++ b/src/cli/vis_utils/simulate/mod.rs
@@ -416,7 +416,11 @@ impl VisSimParams {
         let time_res = Duration::from_seconds(*time_res);
         let timestamps = {
             let mut timestamps = Vec::with_capacity(*num_timesteps);
-            let start = Epoch::from_gpst_seconds(metafits.sched_start_gps_time_ms as f64 / 1e3)
+            let start_ns = metafits
+                .sched_start_gps_time_ms
+                .checked_mul(1_000_000)
+                .expect("does not overflow u64");
+            let start = Epoch::from_gpst_nanoseconds(start_ns)
                 + time_res / 2
                 + Duration::from_seconds(*time_offset);
             for i in 0..*num_timesteps {
diff --git a/src/io/read/raw/mod.rs b/src/io/read/raw/mod.rs
index 7178695f..c270217c 100644
--- a/src/io/read/raw/mod.rs
+++ b/src/io/read/raw/mod.rs
@@ -215,7 +215,12 @@ impl RawDataReader {
             }
         }
 
-        let time_res = Duration::from_milliseconds(metafits_context.corr_int_time_ms as f64);
+        let time_res = {
+            let int_time_ns = i128::from(metafits_context.corr_int_time_ms)
+                .checked_mul(1_000_000)
+                .expect("does not overflow i128");
+            Duration::from_total_nanoseconds(int_time_ns)
+        };
 
         let all_timesteps = Vec1::try_from_vec(mwalib_context.provided_timestep_indices.clone())
             .map_err(|_| RawReadError::NoTimesteps)?;
@@ -223,9 +228,11 @@ impl RawDataReader {
             .timesteps
             .iter()
             .map(|t| {
-                Epoch::from_gpst_seconds(
-                    (t.gps_time_ms + metafits_context.corr_int_time_ms / 2) as f64 / 1e3,
-                )
+                let gps_nanoseconds = t
+                    .gps_time_ms
+                    .checked_mul(1_000_000)
+                    .expect("does not overflow u64");
+                Epoch::from_gpst_nanoseconds(gps_nanoseconds) + time_res / 2
             })
             .collect();
         let timestamps = Vec1::try_from_vec(timestamps).map_err(|_| RawReadError::NoTimesteps)?;
@@ -359,12 +366,20 @@ impl RawDataReader {
             // cotter has a nasty bug that can cause the start time listed in
             // mwaf files to be offset from data HDUs. Warn the user if this is
             // noticed.
-            let data_start =
-                Epoch::from_gpst_seconds(mwalib_context.common_start_gps_time_ms as f64 / 1e3)
-                    + time_res / 2;
-            let data_end =
-                Epoch::from_gpst_seconds(mwalib_context.common_end_gps_time_ms as f64 / 1e3)
-                    + time_res / 2;
+            let data_start = {
+                let gps_nanoseconds = mwalib_context
+                    .common_start_gps_time_ms
+                    .checked_mul(1_000_000)
+                    .expect("does not overflow u64");
+                Epoch::from_gpst_nanoseconds(gps_nanoseconds) + time_res / 2
+            };
+            let data_end = {
+                let gps_nanoseconds = mwalib_context
+                    .common_end_gps_time_ms
+                    .checked_mul(1_000_000)
+                    .expect("does not overflow u64");
+                Epoch::from_gpst_nanoseconds(gps_nanoseconds) + time_res / 2
+            };
             let flags_start = f.start_time;
             let flags_end = flags_start + f.num_time_steps as f64 * time_res;
             let diff = (flags_start - data_start).to_seconds() / time_res.to_seconds();