Load and scale the maps of the noise from the atmosphere
import toml
config = toml.load("s4_design.toml")
log.basicConfig(level=log.INFO)
def get_telecope_years(config, site, channel):
"""Compute the number of telescope/years in the CMB-S4 configuration
config_telescopes : dict
CMB-S4 telescopes configuration,
generally loaded from a TOML file
site : str
'Pole' or 'Chile', case doesn't matter
channel : str
Channel tag, e.g. 'MFHS1'
"""
telescope_years = 0
for telescope_name, telescope_config in config["telescopes"][
get_telescope(channel)
].items():
if telescope_config["site"].lower() == site.lower():
has_band = telescope_config.get(channel[:-1], 0) > 0
telescope_years += has_band * telescope_config.get(
"years", config["experiment"]["total_experiment_length_years"]
)
return telescope_years
s4 = read_instrument_model()
!tail -30 s4_design.toml
for site in ["Pole", "Chile"]:
for channel in s4:
telescope_years = get_telecope_years(config, site, channel["band"])
print(site, channel["band"], telescope_years)
telescope = get_telescope(channel["band"])
if site == "Chile":
if telescope == "SAT":
assert telescope_years == 0, "All SAT at Pole"
elif channel["band"].startswith("ULFL"):
assert telescope_years == 0, "No ULFL in Chile"
else:
if "P" in channel["band"]:
assert telescope_years == 0, "Pole only LAT"
else:
assert telescope_years == 14, "2 LAT each band"
if site == "Pole":
if telescope == "SAT":
assert telescope_years == 14, "2 SAT telescopes each band"
else:
if "P" in channel["band"]:
assert telescope_years == 7, "1 LAT with all 4 bands"
else:
assert telescope_years == 0, "Chile only LAT"
print(30 * "=")
def load_atmosphere(config, site, channel, realization=0, raw=False):
"""Load foreground maps for a channel
Parameters
----------
config : dict
CMB-S4 configuration,
generally loaded from a TOML file
site : str
'Pole' or 'Chile', case doesn't matter
channel : str
Channel tag, e.g. 'MFHS1'
realization : int
Choose one of the available 8 realizations
Returns
-------
output_map : numpy array
Output map with all emissions combined, uses nan for missing pixels
"""
telescope = get_telescope(channel)
channel_noP = channel.replace("P", "")
map_filename = (
Path(f"{realization:08d}")
/ f"{site.lower()}_atmosphere_{telescope}_{channel_noP}_{mapmaking_naming[telescope]}"
)
log.info(f"Reading {map_filename}")
atmosphere_map = hp.read_map(
Path(base_folder) / map_filename, (0, 1, 2), dtype=None, verbose=False
)
if raw:
atmosphere_map[atmosphere_map == 0] = hp.UNSEEN
return atmosphere_map
atmosphere_map[atmosphere_map == hp.UNSEEN] = np.nan
atmosphere_map[atmosphere_map == 0] = np.nan
# input map is 10 days at 100% efficiency
atmosphere_map *= np.sqrt(
10
* simulations_observing_efficiency[site.lower()][channel]
/ (
365.25
* get_observing_efficiency(
config["experiment"]["observing_efficiency"], site, telescope, channel
)
)
)
atmosphere_map /= np.sqrt(get_telecope_years(config, site, channel))
atmosphere_map[0] *= config["experiment"].get("atmosphere_scaling_T", 1)
atmosphere_map[1:] *= config["experiment"].get("atmosphere_scaling_P", 1)
atmosphere_map[1:] /= np.sqrt(get_thinfp(channel))
return atmosphere_map
filenames = !ls $base_folder/00000000/*atmo*
import os.path
for f in map(os.path.basename, filenames):
print(f)
channel = "HFS1"
site = "Pole"
telescope = get_telescope(channel)
input_map = hp.read_map(
base_folder + "/00000000/" + \
f"{site.lower()}_atmosphere_{telescope}_{channel}_telescope_all_time_all_filtered.fits.gz"
, (0,1,2))
input_map[input_map == 0] = hp.UNSEEN
output_map = load_atmosphere(config, site, channel, realization=0)
np.nanmin(output_map), np.nanmax(output_map)
output_map[np.isnan(output_map)] = hp.UNSEEN
np.testing.assert_allclose(input_map == hp.UNSEEN, output_map == hp.UNSEEN)
channel
observing_efficiency = get_observing_efficiency(
config["experiment"]["observing_efficiency"], site, telescope, channel
)
input_map[input_map == hp.UNSEEN] = np.nan
expected_map = input_map * np.sqrt(10 * simulations_observing_efficiency[site.lower()][channel] / (7 * 365.25 * observing_efficiency) / 2)
expected_map[1:] /= np.sqrt(get_thinfp(channel))
np.testing.assert_allclose(expected_map[np.logical_not(np.isnan(expected_map))],
output_map[output_map != hp.UNSEEN], rtol=1e-6)
%matplotlib inline
hp.mollview(output_map[0], min=-1e-5, max=1e-5, unit="K", title="Atmosphere I")
There are no systematics in the half-wave plate or bandpass mismatch, so almost all the atmosphere signal is rejected in polarization.
hp.mollview(output_map[1], min=-1e-9, max=1e-9, unit="K", title="Atmosphere Q")
hp.mollview(output_map[2], min=-1e-9, max=1e-9, unit="K", title="Atmosphere U")
config["experiment"]["atmosphere_scaling_T"] = 3
config["experiment"]["atmosphere_scaling_P"] = 0.1
output_map = load_atmosphere(config, site, channel, realization=0)
output_map[np.isnan(output_map)] = hp.UNSEEN
np.testing.assert_allclose(config["experiment"]["atmosphere_scaling_T"] * expected_map[0, np.logical_not(np.isnan(expected_map[0]))],
output_map[0][output_map[0] != hp.UNSEEN], rtol=1e-6)
np.testing.assert_allclose(config["experiment"]["atmosphere_scaling_P"] * expected_map[1:][np.logical_not(np.isnan(expected_map[1:]))],
output_map[1:][output_map[1:] != hp.UNSEEN], rtol=1e-6)