fix bug with RM synthesis failing to handle cal scans, fixes #38

Author: joshoewahp

dstools/dynamic_spectrum.py

@@ -627,7 +627,8 @@ def rm_synthesis(self, I, L):
         # to pull Q / U from the cleaned RM synthesis FDF spectrum.
         fdf_peak_timeseries = np.argmax(np.abs(fdf_spectrum), axis=0)
         peak_rm_spectrum = phis[fdf_peak_timeseries]
-        tslice = np.argmax(np.nanmean(I.real, axis=1))
+        tslice = np.nanargmax(np.nanmean(I.real, axis=1))
+
         RM = peak_rm_spectrum[tslice]
         logger.debug(f"Peak RM of {RM:.1f} rad/m2")
 

tests/test_dynamic_spectrum.py

@@ -155,6 +155,22 @@ def test_rm_synthesis(ds_paths):
     assert round(ds.RM, 1) == -830.5
 
 
+def test_rm_synthesis_handles_calibrator_scan(ds_paths):
+    ds_path = ds_paths.get("atca_pulse")
+    ds = DynamicSpectrum(ds_path)
+
+    # Simulate early cal scan
+    ds.data["I"][5, :] = np.nan
+    ds.data["L"][5, :] = np.nan
+
+    # Test that RM extraction runs on the pulse at later time index
+    I = ds.data["I"]
+    L = ds.data["L"].T
+    RM = ds.rm_synthesis(I, L)
+
+    assert round(RM, 1) == -830.5
+
+
 def test_ds_acf(ds_paths):
     ds_path = ds_paths.get("atca_pulse")
     ds = DynamicSpectrum(ds_path)