fix the units on FHD reference frequencies

CHANGELOG.md

@@ -8,6 +8,9 @@ information with a value per component.
 - Updated the pyuvdata requirement to >= 2.4.1
 - Updated the scipy requirement to >= 1.5
 
+### Fixed
+- A bug where FHD frequencies were interpreted as being in Hz rather than MHz.
+
 ## [0.3.0] - 2023-04-10
 
 ### Added

docs/tutorial.rst

@@ -148,10 +148,8 @@ a) using extended_model_group attribute
   point
   >>> print(sm.spectral_type)
   spectral_index
-  >>> # correction done since catalog reference frequencies had wrong power
-  >>> sm.reference_frequency = sm.reference_frequency*10**6
   >>> print(np.unique(sm.reference_frequency))
-  [1.82435013e+08 2.15675003e+08] Hz
+  [1.82435e+08 2.15675e+08] Hz
   >>> print(np.unique(sm.spectral_index))
   [-0.8]
   >>> print(np.unique(sm.extended_model_group))
@@ -184,7 +182,7 @@ a) using extended_model_group attribute
 
   >>> # confirming that there is one reference frequency for this extended model group
   >>> print(np.unique(sm.reference_frequency[index_32768]))
-  [2.15675003e+08] Hz
+  [2.15675e+08] Hz
 
   >>> # plots of fluxes are sensible at one frequency since fluxes can change with frequency, plots below provide fluxes
   >>> # when frequency = reference frequency (more on this in at_frequencies section)
@@ -870,15 +868,13 @@ b) spectral index spectral type
   >>> filename = os.path.join(DATA_PATH, "fhd_catalog.sav")
   >>> sm.read_fhd_catalog(filename)
 
-  >>> # correction done since catalog reference frequencies had wrong power
-  >>> sm.reference_frequency = sm.reference_frequency*10**6
-  >>> print(np.unique(sm.reference_frequency))
-  [7.40000000e+07 1.80000000e+08 1.81000000e+08 2.15675003e+08] Hz
+  >>> print(np.unique(sm.reference_frequency.to("MHz")))
+  [ 74.    180.    181.    215.675] MHz
 
   >>> print(sm.stokes.value[0,0,8235])
   0.5017849802970886
   >>> print(sm.reference_frequency[8235])
-  215675003.0517578 Hz
+  215675008.0 Hz
   >>> # last component (at index 8325) was chosen due to nonzero spectral index
   >>> print(sm.spectral_index[8235])
   -0.8

src/pyradiosky/skymodel.py

@@ -4337,7 +4337,8 @@ def read_fhd_catalog(
         ids = catalog["id"].astype(str)
         ra = catalog["ra"]
         dec = catalog["dec"]
-        source_freqs = catalog["freq"]
+        # FHD catalogs frequencies are in MHz
+        source_freqs = catalog["freq"] * 1e6 * units.Hz
         spectral_index = catalog["alpha"]
         Nsrcs = len(catalog)
         extended_model_group = np.full(Nsrcs, "", dtype="<U10")
@@ -4428,7 +4429,9 @@ def read_fhd_catalog(
                         beam_amp_new[:, use_index : use_index + Ncomps] = beam_amp_ext
                         beam_amp_new[:, use_index + Ncomps :] = beam_amp[:, use_index:]
                         beam_amp = beam_amp_new
-                    source_freqs = np.insert(source_freqs, use_index, src["freq"])
+                    source_freqs = np.insert(
+                        source_freqs, use_index, src["freq"] * 1e6 * units.Hz
+                    )
                     spectral_index = np.insert(spectral_index, use_index, src["alpha"])
 
         ra = Longitude(ra, units.deg)
@@ -4443,7 +4446,7 @@ def read_fhd_catalog(
             frame="icrs",
             stokes=stokes,
             spectral_type="spectral_index",
-            reference_frequency=Quantity(source_freqs, "hertz"),
+            reference_frequency=source_freqs,
             spectral_index=spectral_index,
             beam_amp=beam_amp,
             extended_model_group=extended_model_group,

tests/test_skymodel.py

@@ -2956,6 +2956,8 @@ def test_fhd_catalog_reader():
     catalog = scipy.io.readsav(catfile)["catalog"]
     assert skyobj.Ncomponents == len(catalog)
 
+    assert np.all(skyobj.reference_frequency > 50 * units.MHz)
+
 
 @pytest.mark.parametrize("extended", [True, False])
 def test_fhd_catalog_reader_extended_sources(extended):