Merge pull request #114 from zonca/fix_units_co

Fix units in CO
galsci · Mar 12, 2022 · 9f51b6c · 9f51b6c
2 parents 7c55a86 + 3f0cfd1
commit 9f51b6c
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 13 deletions.
diff --git a/pysm3/models/co_lines.py b/pysm3/models/co_lines.py
@@ -139,10 +139,10 @@ def get_emission(self, freqs: u.GHz, weights=None) -> u.uK_RJ:
             if line_freq >= freqs[0] and line_freq <= freqs[-1]:
                 weight = np.interp(line_freq, freqs, weights)
                 convert_to_uK_RJ = (1 * u.K_CMB).to_value(
-                    u.K_RJ,
+                    u.uK_RJ,
                     equivalencies=u.cmb_equivalencies(line_freq * u.GHz),
                 )
-                I_map = self.planck_templatemap[line]
+                I_map = self.planck_templatemap[line].copy()
                 if self.include_high_galactic_latitude_clouds:
                     I_map += self.simulate_high_galactic_latitude_CO(line)
 

diff --git a/pysm3/tests/test_co.py b/pysm3/tests/test_co.py
@@ -36,29 +36,27 @@ def test_co(include_high_galactic_latitude_clouds):
     expected_map_filename = remote_data.get(
         "co/testing/CO10_TQUmaps_{}_nside16_ring.fits.zip".format(tag)
     )
+
     expected_co_map = (
-        hp.read_map(expected_map_filename, field=(0, 1, 2), dtype=np.float64) * u.uK_CMB
+        hp.read_map(expected_map_filename, field=(0, 1, 2), dtype=np.float64) * u.K_CMB
     ).to(u.uK_RJ, equivalencies=u.cmb_equivalencies(line_freq))
 
     assert_quantity_allclose(co_map, expected_co_map, rtol=1e-5)
 
-    # to avoid numerical errors, only compare part of the map where signal is more than 1e-3 uK
-    nonzero_values = expected_co_map > 1e-3 * u.uK_RJ
-
     co_map = co.get_emission([114.271, 115.271, 116.271, 117.271] * u.GHz)
     # weight is about 1/3 as bandwidth is 3 GHz
     assert_quantity_allclose(
-        co_map[nonzero_values],
-        expected_co_map[nonzero_values] * 0.3304377039951613,
-        rtol=1e-4,
+        co_map,
+        expected_co_map * 0.3304377039951613,
+        rtol=1e-3
     )
 
     co_map = co.get_emission([100, 120] * u.GHz)
     # weight is about 1/20, consider normalization also includes conversion to power units
     assert_quantity_allclose(
-        co_map[nonzero_values],
-        expected_co_map[nonzero_values] * 0.05475254098360655,
-        rtol=1e-4,
+        co_map,
+        expected_co_map * 0.05475254098360655,
+        rtol=1e-4
     )
 
 
@@ -76,7 +74,7 @@ def test_co_model(model_tag):
         "co/testing/CO10_TQUmaps_{}_nside16_ring.fits.zip".format(tag)
     )
     expected_co_map = (
-        hp.read_map(expected_map_filename, field=(0, 1, 2), dtype=np.float64) * u.uK_CMB
+        hp.read_map(expected_map_filename, field=(0, 1, 2), dtype=np.float64) * u.K_CMB
     )
 
     assert_quantity_allclose(co_map, expected_co_map, rtol=1e-5)