Sky emission

load_sky_emission[source]

load_sky_emission(config_sky_emission, site, channel)

Load foreground maps for a channel

Parameters
----------
config_sky_emission : dict
    CMB-S4 sky emission configuration,
    generally config['sky_emission'] for a configuration
    loaded from a TOML file
    dictionary with standard emission names and their weights
site : str
    'Pole' or 'Chile', case doesn't matter
channel : str
    Channel tag, e.g. 'MFHS1'

Returns
-------
output_map : numpy array
    Output map with all emissions combined, uses nan for missing pixels

def load_sky_emission(config_sky_emission, site, channel):
    """Load foreground maps for a channel

    Parameters
    ----------
    config_sky_emission : dict
        CMB-S4 sky emission configuration,
        generally config['sky_emission'] for a configuration
        loaded from a TOML file
        dictionary with standard emission names and their weights
    site : str
        'Pole' or 'Chile', case doesn't matter
    channel : str
        Channel tag, e.g. 'MFHS1'

    Returns
    -------
    output_map : numpy array
        Output map with all emissions combined, uses nan for missing pixels
    """
    log.info("Configuration %s", str(config_sky_emission))
    telescope = get_telescope(channel)
    nside = 512 if telescope == "SAT" else 4096
    npix = hp.nside2npix(nside)
    output_map = np.zeros((3, npix), dtype=np.float32)
    realization = 0  # foregrounds are deterministic
    channel_noP = channel.replace("P", "")
    for emission, weight in config_sky_emission.items():
        if weight == 0:
            log.info("Skip %s", emission)
            continue
        log.info("Processing %s", emission)

        emission_map = hp.read_map(
            Path(base_folder)
            / f"{realization:08d}"
            / f"{site.lower()}_{emission_naming[emission]}_{telescope}_{channel_noP}_{mapmaking_naming[telescope]}",
            (0, 1, 2),
            dtype=None,
            verbose=False,
        )
        emission_map[emission_map == 0] = np.nan
        emission_map[emission_map == hp.UNSEEN] = np.nan

        if emission == "CMB_tensor_to_scalar_ratio":
            # input maps are simulated with r=3e-3
            # tensor-to-scalar ratio is defined in the power spectrum, we need to weight the maps with the square root
            weight = np.sqrt(weight / 3e-3)
        output_map += emission_map * weight

    output_map /= 1e6  # uK -> K
    return output_map

import toml

config = toml.load("s4_design.toml")

s4 = read_instrument_model()

s4

QTable length=21

band	telescope	center_frequency	fwhm	detectors_per_tube	NET	chile_temp_sensitivity_DSR	chile_pol_sensitivity_DSR	pole_pol_sensitivity_DSR
		GHz	arcmin		s(1/2) uK	arcmin uK	arcmin uK	arcmin uK
str6	str3	float64	float64	int64	float64	float64	float64	float64
ULFPL1	LAT	20.0	11.4	54	332.2	0.0	0.0	8.4
LFL1	LAT	25.75	7.4	96	275.0	21.8	30.8	5.0
LFPL1	LAT	25.75	8.4	96	278.6	21.8	30.8	5.0
LFL2	LAT	38.75	5.1	96	215.0	12.4	17.6	4.5
LFPL2	LAT	38.75	5.8	96	268.6	12.4	17.6	4.5
MFPL1	LAT	91.5	2.5	864	285.2	2.0	2.9	0.68
MFL1	LAT	91.5	2.2	864	251.0	2.0	2.9	0.68
MFL2	LAT	148.5	1.4	864	280.0	2.0	2.8	0.96
MFPL2	LAT	148.5	1.6	864	264.6	2.0	2.8	0.96
HFL1	LAT	227.0	1.0	864	625.0	6.9	9.8	5.7
HFPL1	LAT	227.0	1.1	864	534.9	6.9	9.8	5.7
HFL2	LAT	285.5	0.9	864	1528.0	16.7	23.6	9.8
HFPL2	LAT	285.5	1.0	864	1228.3	16.7	23.6	9.8
LFS1	SAT	25.75	72.8	288	169.2	0.0	0.0	3.5
LFS2	SAT	38.75	72.8	288	204.3	0.0	0.0	4.5
MFLS1	SAT	85.0	25.5	3528	290.0	0.0	0.0	0.88
MFHS1	SAT	95.0	22.7	3528	248.2	0.0	0.0	0.78
MFLS2	SAT	145.0	25.5	3528	280.4	0.0	0.0	1.2
MFHS2	SAT	155.1	22.7	3528	297.4	0.0	0.0	1.3
HFS1	SAT	227.0	13.0	8592	609.1	0.0	0.0	3.5
HFS2	SAT	285.5	13.0	8592	1434.9	0.0	0.0	6.0

log.basicConfig(level=log.INFO)

ch = "LFS1"
site = "Pole"

output_map = load_sky_emission(config["sky_emission"], site, ch)

INFO:root:Configuration {'foreground_emission': 1, 'CMB_unlensed': 1, 'CMB_lensing_signal': 1, 'CMB_tensor_to_scalar_ratio': 0.005}
INFO:root:Processing foreground_emission
INFO:root:Processing CMB_unlensed
INFO:root:Processing CMB_lensing_signal
INFO:root:Processing CMB_tensor_to_scalar_ratio

assert len(output_map[output_map == hp.UNSEEN]) == 0, \
    "We should be using nan not UNSEEN to simplify summing"

np.nanmin(output_map), np.nanmax(output_map)

(-0.0015811161, 0.0018167086)

assert np.nanmax(output_map) > 1e-4 and np.nanmax(output_map) < 1e-2, "Rough unit check"

%matplotlib inline

hp.mollview(output_map[0], min=-1e-4, max=1e-4, unit="K", title="Foregrounds I")

/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/projaxes.py:920: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
  newcm.set_over(newcm(1.0))
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/projaxes.py:921: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
  newcm.set_under(bgcolor)
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/projaxes.py:922: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
  newcm.set_bad(badcolor)
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/projaxes.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
  **kwds

hp.mollview(output_map[1], min=-1e-5, max=1e-5, unit="K", title="Foregrounds Q")

hp.mollview(output_map[2], min=-1e-5, max=1e-5, unit="K", title="Foregrounds U")