Import COLines into PySM 3

docs/co_lines.rst

@@ -0,0 +1,43 @@
+.. _colines:
+
+COLines
+=======
+
+This class implements simulations for Galactic CO emission involving the first 3 CO rotational lines, i.e. :math:`J=1-0,2-1,3-2` whose center frequency is respectively at :math:`\nu_0 = 115.3, 230.5,345.8` GHz. The CO emission map templates are the CO Planck maps obtained with ``MILCA`` component separation algorithm (See `Planck paper <https://www.aanda.org/articles/aa/abs/2014/11/aa21553-13/aa21553-13.html>`). The CO maps have been released at the nominal resolution (10 and 5 arcminutes). However, to reduce  noise contamination from template maps (especially at intermediate and high Galactic latitudes), we  convolved them with a 1 deg gaussian beam.
+
+The Stokes I map is computed from the template one as it follows:
+
+if target :math:`N_{side}` <= 512:
+
+    #. The template map at a ``nside=512``  is downgraded at the target :math:`N_{side}`
+
+if target :math:`N_{side}` > 512 :
+
+    #. The template map at a ``nside=2048``  is downgraded(eventually upgraded) at the target :math:`N_{side}`
+
+Q and U maps can be computed from the template CO emission  map, :math:`I_{CO}`,  assuming a constant  fractional polarization, as:
+
+.. math::
+
+    Q = f_{pol} I_{CO}  g_d \cos( 2 \psi)
+
+    U  = f_{pol} I_{CO}  g_d \sin( 2 \psi)
+
+with :math:`g_d` and :math:`\psi` being respectively the depolarization and polarization angle maps estimated from a dust map as :
+
+.. math::
+
+    g_d = \frac{ \sqrt{Q^2_{d,353}    + U^2_{d,353}   } }{f_{pol} I_{d,353} }
+
+    \psi = \frac{1}{2} \arctan {\frac{U_{d,353}}{Q_{d,353}}}
+
+
+Most of the CO emission is expected to be confined in the  Galactic midplane. However, there are still regions at high Galactic latitudes  where the CO emission has been purely assessed (by current surveys) and where the Planck signal-to-noise was not enough to detect any emission.
+
+The PySM user can include the eventuality of molecular emission (both unpolarized and polarized) at High Gal. Latitudes by co-adding to the emission maps one realization of CO emission simulated with MCMole3D together with  the Planck CO map. The polarization is simulated similarly as above.
+
+The ``MCMole3D`` input parameters  are are obtained from best fit with the Planck CO 1-0 map (see Puglisi et al. 2017 and the `documentation <http://giuspugl.github.io/mcmole/index.html>`). If ``include_high_galactic_latitude_clouds=True``, a mock CO cloud map is simulated with ``MCMole3D``, encoding high Galactic latitudes clouds at latitudes above and below  than 20 degrees. The mock emission map is then co-added to the Planck CO emission map. The polarization is simulated similarly as above.
+
+The installation of ``mcmole3d`` is not required, HGL clouds can be input to the CO emission by setting ``run_mcmole3d=False``  (which is the default). However, if one wants to run several mock CO  realizations observing high Galactic latitude patches we encourage to run ``mcmole3d`` by changing ``random_seed`` in the CO class constructor. The parameter ``theta_high_galactic_latitude_deg`` set the latitude above which CO emission from high Galactic latitudes can be included and it has an impact **only when** ``run_mcmole3d=True``.
+
+The level of polarization in **co2** and **co3** is 0.1%, which on average is the expected level on 10% of the sky. However, polarization from CO emission have been detected at larger fluxes in  Orion and Taurus complexes (Greaves et al.1999 )

docs/models.rst

@@ -80,3 +80,10 @@ CMB
 
 - **c1**: A lensed CMB realisation is computed using Taylens, a code to compute a lensed CMB realisation using nearest-neighbour Taylor interpolation (`taylens <https://github.com/amaurea/taylens>`_; Naess, S. K. and Louis, T. JCAP 09 001, 2013, astro-ph/1307.0719). This code takes, as an input, a set of unlensed Cl's generated using `CAMB <http://www.camb.info/>`_. The params.ini is in the Ancillary directory. There is a pre-computed CMB map provided at Nside 512.
 
+CO line emission
+================
+For more details see :ref:`colines`.
+
+- **co1**: Galactic CO emission involving the first 3 CO rotational lines, i.e. :math:`J=1-0,2-1,3-2` whose center frequency is respectively at :math:`\nu_0 = 115.3, 230.5,345.8` GHz. The CO emission map templates are the CO Planck maps obtained with ``MILCA`` component separation algorithm (See `Planck paper <https://www.aanda.org/articles/aa/abs/2014/11/aa21553-13/aa21553-13.html>`). The CO maps have been released at the nominal resolution (10 and 5 arcminutes). However, to reduce  noise contamination from template maps (especially at intermediate and high Galactic latitudes), we  convolved them with a 1 deg gaussian beam.
+- **co2**: like **co1** with polarized emission at the level of 0.1%.
+- **co3**: like **co2** with a mock CO clouds map 20 degrees off the Galactic plane simulated with ``MCMole3D``.

pysm3/__init__.py

@@ -1,6 +1,5 @@
 # flake8: noqa
 from ._astropy_init import *  # noqa
-
 import warnings
 from astropy.utils.exceptions import AstropyDeprecationWarning
 

pysm3/data/presets.cfg

@@ -247,3 +247,15 @@ delensing_ells = "pysm_2/delens_ells.txt"
 [c2]
 class = "CMBMap"
 map_IQU = "pysm_2/lensed_cmb.fits"
+[co1]
+class = "COLines"
+has_polarization = false
+[co2]
+class = "COLines"
+has_polarization = true
+polarization_fraction = 0.001
+[co3]
+class = "COLines"
+has_polarization = true
+polarization_fraction = 0.001
+include_high_galactic_latitude_clouds = true

pysm3/models/__init__.py

@@ -8,3 +8,4 @@
 from .interpolating import InterpolatingComponent
 from .cmb import CMBMap, CMBLensed
 from .dust_layers import ModifiedBlackBodyLayers
+from .co_lines import COLines

pysm3/models/co_lines.py

@@ -0,0 +1,244 @@
+import numpy as np
+import healpy as hp
+
+from .. import units as u
+from .. import utils
+from .template import Model
+
+
+def build_lines_dict(lines, maps):
+    """Build a dictionary for lines and maps
+
+    Takes a list of tags (strings) and a map or set of maps
+    and returns a dictionary where each tag is associated with
+    a map
+    """
+    return dict(zip(lines, np.atleast_2d(maps)))
+
+
+class COLines(Model):
+    def __init__(
+        self,
+        nside,
+        has_polarization=True,
+        lines=["10", "21", "32"],
+        include_high_galactic_latitude_clouds=False,
+        polarization_fraction=0.001,
+        theta_high_galactic_latitude_deg=20.0,
+        random_seed=1234567,
+        verbose=False,
+        run_mcmole3d=False,
+        map_dist=None,
+    ):
+
+        """Class defining attributes for CO line emission.
+        CO templates are extracted from Type 1 CO Planck maps.
+        See further details in:
+        https://www.aanda.org/articles/aa/abs/2014/11/aa21553-13/aa21553-13.html
+
+        Parameters
+        ----------
+        nside : int
+            HEALPix NSIDE of the output maps
+        has_polarization : bool
+            whether or not to simulate also polarization maps
+        lines : list of strings
+            CO rotational transitions to consider.
+            Accepted values : 10, 21, 32
+        polarization_fraction: float
+            polarisation fraction for polarised CO emission.
+        include_high_galactic_latitude_clouds: bool
+            If True it includes a simulation from MCMole3D to include
+            high Galactic Latitude clouds.
+            (See more details at http://giuspugl.github.io/mcmole/index.html)
+        run_mcmole3d: bool
+            If True it simulates  HGL cluds by running MCMole3D, otherwise it coadds
+            a map of HGL emission.
+        random_seed: int
+            set random seed for mcmole3d simulations.
+        theta_high_galactic_latitude_deg : float
+            Angle in degree  to identify High Galactic Latitude clouds
+            (i.e. clouds whose latitude b is `|b|> theta_high_galactic_latitude_deg`).
+        map_dist : mpi4py communicator
+            Read inputs across a MPI communicator, see pysm.read_map
+        """
+
+        self.lines = lines
+        self.line_index = {"10": 0, "21": 1, "32": 2}
+        self.line_frequency = {
+            "10": 115.271 * u.GHz,
+            "21": 230.538 * u.GHz,
+            "32": 345.796 * u.GHz,
+        }
+        self.nside = nside
+
+        self.template_nside = 512 if self.nside <= 512 else 2048
+
+        super().__init__(nside=nside, map_dist=map_dist)
+
+        self.remote_data = utils.RemoteData()
+
+        self.planck_templatemap_filename = (
+            "co/HFI_CompMap_CO-Type1_{}_R2.00_ring.fits".format(self.template_nside)
+        )
+        self.planck_templatemap = build_lines_dict(
+            self.lines,
+            hp.ud_grade(
+                map_in=self.read_map(
+                    self.remote_data.get(self.planck_templatemap_filename),
+                    field=[self.line_index[line] for line in self.lines],
+                    unit=u.K_CMB,
+                ),
+                nside_out=self.nside,
+            )
+            << u.K_CMB,
+        )
+
+        self.include_high_galactic_latitude_clouds = (
+            include_high_galactic_latitude_clouds
+        )
+        self.has_polarization = has_polarization
+        if self.has_polarization:
+            self.polangle = self.read_map(
+                self.remote_data.get(
+                    "co/psimap_dust90_{}.fits".format(self.template_nside)
+                )
+            ).value
+            self.depolmap = self.read_map(
+                self.remote_data.get(
+                    "co/gmap_dust90_{}.fits".format(self.template_nside)
+                )
+            ).value
+        self.polarization_fraction = polarization_fraction
+        self.theta_high_galactic_latitude_deg = theta_high_galactic_latitude_deg
+        self.random_seed = random_seed
+        self.run_mcmole3d = run_mcmole3d
+
+        if include_high_galactic_latitude_clouds and not run_mcmole3d:
+            # Dictionary where keys are "10", "21" and values
+            self.mapclouds = build_lines_dict(
+                self.lines,
+                self.read_map(
+                    self.remote_data.get(
+                        "co/mcmoleCO_HGL_{}.fits".format(self.template_nside)
+                    ),
+                    field=[self.line_index[line] for line in self.lines],
+                    unit=u.K_CMB,
+                ),
+            )
+
+        self.verbose = verbose
+
+    @u.quantity_input
+    def get_emission(self, freqs: u.GHz, weights=None) -> u.uK_RJ:
+        freqs = utils.check_freq_input(freqs)
+        weights = utils.normalize_weights(freqs, weights)
+        out = np.zeros((3, hp.nside2npix(self.nside)), dtype=np.double)
+        for line in self.lines:
+            line_freq = self.line_frequency[line].to_value(u.GHz)
+            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,
+                    equivalencies=u.cmb_equivalencies(line_freq * u.GHz),
+                )
+                I_map = self.planck_templatemap[line]
+                if self.include_high_galactic_latitude_clouds:
+                    I_map += self.simulate_high_galactic_latitude_CO(line)
+
+                if self.has_polarization:
+                    out[1:] += (
+                        self.simulate_polarized_emission(I_map).value
+                        * convert_to_uK_RJ
+                        * weight
+                    )
+                out[0] += I_map.value * convert_to_uK_RJ * weight
+
+        return out << u.uK_RJ
+
+    def simulate_polarized_emission(self, I_map):
+        """
+        Add polarized emission by means of:
+        * an overall constant polarization fraction,
+        * a depolarization map to mimick the line of sight depolarization effect
+          at low Galactic latitudes
+        * a polarization angle map coming from a dust template (we exploit the observed correlation
+        between polarized dust and molecular emission in star forming regions).
+        """
+
+        cospolangle = np.cos(2.0 * self.polangle)
+        sinpolangle = np.sin(2.0 * self.polangle)
+
+        P_map = self.polarization_fraction * self.depolmap * I_map
+        return P_map * np.array([cospolangle, sinpolangle])
+
+    def simulate_high_galactic_latitude_CO(self, line):
+        """
+        Coadd High Galactic Latitude CO emission, simulated with  MCMole3D.
+        """
+        if self.run_mcmole3d:
+            import mcmole3d as cl
+
+            # params to MCMole
+            N = 40000
+            L_0 = 20.4  # pc
+            L_min = 0.3
+            L_max = 60.0
+            R_ring = 5.8
+            sigma_ring = 2.7  # kpc
+            R_bulge = 3.0
+            R_z = 10  # kpc
+            z_0 = 0.1
+            Em_0 = 240.0
+            R_em = 6.6
+            model = "LogSpiral"
+
+            nside = self.nside
+            Itot_o, _ = cl.integrate_intensity_map(
+                self.planck_templatemap[line],
+                hp.get_nside(self.planck_templatemap[line]),
+                planck_map=True,
+            )
+            Pop = cl.Cloud_Population(N, model, randseed=self.random_seed)
+
+            Pop.set_parameters(
+                radial_distr=[R_ring, sigma_ring, R_bulge],
+                typical_size=L_0,
+                size_range=[L_min, L_max],
+                thickness_distr=[z_0, R_z],
+                emissivity=[Em_0, R_em],
+            )
+            Pop()
+
+            if self.verbose:
+                Pop.print_parameters()
+            # project into  Healpix maps
+            mapclouds = cl.do_healpy_map(
+                Pop,
+                nside,
+                highgalcut=np.deg2rad(90.0 - self.theta_high_galactic_latitude_deg),
+                apodization="gaussian",
+                verbose=self.verbose,
+            )
+            Itot_m, _ = cl.integrate_intensity_map(mapclouds, nside)
+            # convert simulated map into the units of the Planck one
+            rescaling_factor = Itot_m / Itot_o
+            mapclouds /= rescaling_factor
+            hglmask = np.zeros_like(mapclouds)
+            # Apply mask to low galactic latitudes
+            listhgl = hp.query_strip(
+                nside,
+                np.deg2rad(90.0 + self.theta_high_galactic_latitude_deg),
+                np.deg2rad(90 - self.theta_high_galactic_latitude_deg),
+            )
+            hglmask[listhgl] = 1.0
+            rmsplanck = self.planck_templatemap[line][listhgl].std()
+            rmssim = mapclouds[listhgl].std()
+            if rmssim == 0.0:
+                belowplanck = 1.0
+            else:
+                belowplanck = rmssim / rmsplanck
+
+            return mapclouds * hglmask / belowplanck
+        else:
+            return self.mapclouds[line]