diff --git a/docs/source/detectors.rst b/docs/source/detectors.rst index f5221b34..e356dd0e 100644 --- a/docs/source/detectors.rst +++ b/docs/source/detectors.rst @@ -168,7 +168,6 @@ with four detectors. Here is a summary of its contents: foo3: band center at 67.0 GHz foo4: band center at 68.0 GHz - Now, let's turn back to the problem of specifying a set of detectors in a parameter file. The following TOML file shows some of the possibilities granted by the framework: @@ -215,7 +214,6 @@ The following code will read the TOML file above and produce a list of 5. planck30GHz: band center at 28.4 GHz 6. foo1: band center at 65.0 GHz - You are not forced to use ``detectors`` as the name of the parameter in the TOML file, as :func:`.detector_list_from_parameters` accepts a generic list. As an example, consider the following TOML file: diff --git a/litebird_sim/constants.py b/litebird_sim/constants.py new file mode 100644 index 00000000..1d268669 --- /dev/null +++ b/litebird_sim/constants.py @@ -0,0 +1,14 @@ +# -*- encoding: utf-8 -*- +from astropy.constants import c as c_light +from astropy.constants import h, k_B + +C_LIGHT_KM_S = c_light.value / 1e3 +H_OVER_K_B = h.value / k_B.value + +T_CMB_K = 2.72548 # Fixsen 2009 http://arxiv.org/abs/0911.1955 + +SOLAR_VELOCITY_KM_S = 369.8160 +SOLAR_VELOCITY_GAL_LAT_RAD = 0.842_173_724 +SOLAR_VELOCITY_GAL_LON_RAD = 4.608_035_744_4 + +EARTH_L2_DISTANCE_KM = 1_496_509.30522 diff --git a/litebird_sim/dipole.py b/litebird_sim/dipole.py index 30807c00..438bdb49 100644 --- a/litebird_sim/dipole.py +++ b/litebird_sim/dipole.py @@ -7,14 +7,10 @@ from typing import Union, List -from astropy.constants import c as c_light -from astropy.constants import h, k_B - from .observations import Observation from .spacecraft import SpacecraftPositionAndVelocity -C_LIGHT_KM_S = c_light.value / 1e3 -H_OVER_K_B = h.value / k_B.value +from litebird_sim import constants as c # We use a IntEnum class so that comparisons are much faster than with strings @@ -65,7 +61,7 @@ class DipoleType(IntEnum): @njit def planck(nu_hz, t_k): """Return occupation number at frequency nu_hz and temperature t_k""" - return 1 / (np.exp(H_OVER_K_B * nu_hz / t_k) - 1) + return 1 / (np.exp(c.H_OVER_K_B * nu_hz / t_k) - 1) @njit @@ -79,15 +75,15 @@ def compute_scalar_product(theta, phi, v): @njit def calculate_beta(theta, phi, v_km_s): """Return a 2-tuple containing β·n and β""" - beta_dot_n = compute_scalar_product(theta, phi, v_km_s) / C_LIGHT_KM_S - beta = np.sqrt(v_km_s[0] ** 2 + v_km_s[1] ** 2 + v_km_s[2] ** 2) / C_LIGHT_KM_S + beta_dot_n = compute_scalar_product(theta, phi, v_km_s) / c.C_LIGHT_KM_S + beta = np.sqrt(v_km_s[0] ** 2 + v_km_s[1] ** 2 + v_km_s[2] ** 2) / c.C_LIGHT_KM_S return beta_dot_n, beta @njit def compute_dipole_for_one_sample_linear(theta, phi, v_km_s, t_cmb_k): - beta_dot_n = compute_scalar_product(theta, phi, v_km_s) / C_LIGHT_KM_S + beta_dot_n = compute_scalar_product(theta, phi, v_km_s) / c.C_LIGHT_KM_S return t_cmb_k * beta_dot_n @@ -212,7 +208,7 @@ def add_dipole( nu_hz = frequency_ghz[detector_idx] * 1e9 # freq in GHz # Note that x is a dimensionless parameter - x = h.value * nu_hz / (k_B.value * t_cmb_k) + x = c.H_OVER_K_B * nu_hz / t_cmb_k f_x = x * np.exp(x) / (np.exp(x) - 1) @@ -235,7 +231,7 @@ def add_dipole_to_observations( obs: Union[Observation, List[Observation]], pos_and_vel: SpacecraftPositionAndVelocity, pointings: Union[np.ndarray, List[np.ndarray], None] = None, - t_cmb_k: float = 2.72548, # Fixsen 2009 http://arxiv.org/abs/0911.1955 + t_cmb_k: float = c.T_CMB_K, dipole_type: DipoleType = DipoleType.TOTAL_FROM_LIN_T, frequency_ghz: Union[ np.ndarray, None diff --git a/litebird_sim/mbs/mbs.py b/litebird_sim/mbs/mbs.py index 73ab2e36..db408bc7 100644 --- a/litebird_sim/mbs/mbs.py +++ b/litebird_sim/mbs/mbs.py @@ -15,6 +15,7 @@ import pysm3.units as u import litebird_sim as lbs +from litebird_sim import constants as c @dataclass @@ -93,6 +94,8 @@ class generates sky maps. You can choose to include the following - Foreground signal (use ``make_fg=True``); + - Solar dipole signal (use ``make_dipole=True``) + - Noise (use ``make_noise=True``); this should be used only when running map-based simulations, because if you are creating time-ordered data, chances are that you want to inject noise @@ -159,6 +162,16 @@ class generates sky maps. You can choose to include the following which associates a name to a component, e.g., ``{"ame": "pysm_ame_1"}``. + - ``make_dipole` (default: ``False``): when ``True``, a solar dipole + map is added to the maps + + - ``sun_velocity`` (default: ``None``): this list specifies the + direction, latitude (rad) and longitude (rad) in galactic coordinates, + and the amplitude in Km/s of the dipole.If ``None`` these values + are taken `.constants`. The dipole implemented is the "Linear + approximation in β using thermodynamic units" as in + :class:``.dipole.DipoleType``. + - ``output_string`` (default: ``""``): a string used to build the file names of the Healpix FITS files saved by the :class:`.Mbs` class. @@ -191,6 +204,8 @@ class generates sky maps. You can choose to include the following seed_cmb: Union[int, None] = None make_fg: bool = False fg_models: Union[Dict[str, Any], None] = None + make_dipole: bool = False + sun_velocity: Union[list, None] = None output_string: Union[str, None] = None units: str = "K_CMB" maps_in_ecliptic: bool = False @@ -324,13 +339,13 @@ def _parse_instrument_from_ch_list(self): len(self.ch_list) except TypeError: self.ch_list = [self.ch_list] - for c in self.ch_list: - name = c.channel.replace(" ", "_") + for ch in self.ch_list: + name = ch.channel.replace(" ", "_") self.instrument[name] = _InstrumentFreq( - bandcenter_ghz=c.bandcenter_ghz, - bandwidth_ghz=c.bandwidth_ghz, - fwhm_arcmin=c.fwhm_arcmin, - p_sens_ukarcmin=c.pol_sensitivity_channel_ukarcmin, + bandcenter_ghz=ch.bandcenter_ghz, + bandwidth_ghz=ch.bandwidth_ghz, + fwhm_arcmin=ch.fwhm_arcmin, + p_sens_ukarcmin=ch.pol_sensitivity_channel_ukarcmin, ) def _parse_instrument(self): @@ -589,14 +604,23 @@ def generate_cmb(self): ) for Nchnl, chnl in enumerate(channels): - freq = instr[chnl].bandcenter_ghz + if hasattr(instr[chnl], "band"): + freq = instr[chnl].band.bandcenter_ghz + else: + freq = instr[chnl].bandcenter_ghz + if self.params.bandpass_int: - band = instr[chnl].bandwidth_ghz - fmin = freq - band / 2.0 - fmax = freq + band / 2.0 - fsteps = int(np.ceil(fmax - fmin) + 1) - bandpass_frequencies = np.linspace(fmin, fmax, fsteps) * u.GHz - weights = np.ones(len(bandpass_frequencies)) + if hasattr(instr[chnl], "band"): + bandpass_frequencies = instr[chnl].band.freqs_ghz + weights = instr[chnl].band.weights + else: + band = instr[chnl].bandwidth_ghz + fmin = freq - band / 2.0 + fmax = freq + band / 2.0 + fsteps = int(np.ceil(fmax - fmin) + 1) + bandpass_frequencies = np.linspace(fmin, fmax, fsteps) * u.GHz + weights = np.ones(len(bandpass_frequencies)) + cmb_map = sky.get_emission(bandpass_frequencies, weights) cmb_map = cmb_map * pysm3.bandpass_unit_conversion( bandpass_frequencies, weights, self.pysm_units @@ -680,15 +704,23 @@ def generate_fg(self): fg_map_matrix = np.zeros((n_channels, 3, npix)) for Nchnl, chnl in enumerate(channels): - freq = instr[chnl].bandcenter_ghz + if hasattr(instr[chnl], "band"): + freq = instr[chnl].band.bandcenter_ghz + else: + freq = instr[chnl].bandcenter_ghz fwhm_arcmin = instr[chnl].fwhm_arcmin if self.params.bandpass_int: - band = instr[chnl].bandwidth_ghz - fmin = freq - band / 2.0 - fmax = freq + band / 2.0 - fsteps = int(np.ceil(fmax - fmin) + 1) - bandpass_frequencies = np.linspace(fmin, fmax, fsteps) * u.GHz - weights = np.ones(len(bandpass_frequencies)) + if hasattr(instr[chnl], "band"): + bandpass_frequencies = instr[chnl].band.freqs_ghz + weights = instr[chnl].band.weights + else: + band = instr[chnl].bandwidth_ghz + fmin = freq - band / 2.0 + fmax = freq + band / 2.0 + fsteps = int(np.ceil(fmax - fmin) + 1) + bandpass_frequencies = np.linspace(fmin, fmax, fsteps) * u.GHz + weights = np.ones(len(bandpass_frequencies)) + sky_extrap = sky.get_emission(bandpass_frequencies, weights) sky_extrap = sky_extrap * pysm3.bandpass_unit_conversion( bandpass_frequencies, weights, self.pysm_units @@ -698,6 +730,7 @@ def generate_fg(self): sky_extrap = sky_extrap.to( self.pysm_units, equivalencies=u.cmb_equivalencies(freq * u.GHz) ) + if smooth: sky_extrap_smt = hp.smoothing( sky_extrap, @@ -724,6 +757,76 @@ def generate_fg(self): else: return (None, saved_maps) + def generate_dipole(self): + parallel = self.params.parallel_mc + instr = self.instrument + nside = self.params.nside + npix = hp.nside2npix(nside) + root_dir = self.sim.base_path + output_directory = root_dir / "dipole" + sun_velocity = self.params.sun_velocity + file_str = self.params.output_string + channels = instr.keys() + n_channels = len(channels) + col_units = [self.params.units] + saved_maps = [] + + if parallel: + comm = lbs.MPI_COMM_WORLD + rank = comm.Get_rank() + else: + comm = None + rank = 0 + + if rank == 0 and self.params.save: + output_directory.mkdir(parents=True, exist_ok=True) + + if sun_velocity is None: + sun_velocity = ( + c.SOLAR_VELOCITY_GAL_LAT_RAD, + c.SOLAR_VELOCITY_GAL_LON_RAD, + c.SOLAR_VELOCITY_KM_S, + ) + + dipvec = ( + hp.ang2vec(sun_velocity[0], sun_velocity[1]) + * sun_velocity[2] + / c.C_LIGHT_KM_S + * c.T_CMB_K + ) + + dipole = np.dot(dipvec, hp.pix2vec(nside, np.arange(npix))) * u.K_CMB + + if not self.params.save: + dipole_map_matrix = np.zeros((n_channels, npix)) + + for Nchnl, chnl in enumerate(channels): + if hasattr(instr[chnl], "band"): + freq = instr[chnl].band.bandcenter_ghz + else: + freq = instr[chnl].bandcenter_ghz + + sky_dipole = dipole.to( + self.pysm_units, equivalencies=u.cmb_equivalencies(freq * u.GHz) + ) + + if self.params.save: + file_name = f"{chnl}_dipole_{file_str}.fits" + cur_map_path = output_directory / file_name + lbs.write_healpix_map_to_file( + cur_map_path, sky_dipole, column_units=col_units + ) + saved_maps.append( + MbsSavedMapInfo(path=cur_map_path, component="dipole", channel=chnl) + ) + else: + dipole_map_matrix[Nchnl] = sky_dipole + + if not self.params.save: + return (dipole_map_matrix, saved_maps) + else: + return (None, saved_maps) + def write_coadded_maps(self, saved_maps): root_dir = self.sim.base_path fg_dir = root_dir / "foregrounds" @@ -847,6 +950,14 @@ def run_all(self): for cmp in fg.keys(): tot += fg[cmp] + if self.params.make_dipole: + log.info("generating and saving dipole") + dipole, dipole_map = self.generate_dipole() + saved_maps += dipole_map + + if not self.params.save: + tot[:, 0, :] += dipole + if self.params.maps_in_ecliptic: r = hp.Rotator(coord=["G", "E"]) diff --git a/litebird_sim/simulations.py b/litebird_sim/simulations.py index 90f7248b..ff9b8027 100644 --- a/litebird_sim/simulations.py +++ b/litebird_sim/simulations.py @@ -13,6 +13,7 @@ from shutil import copyfile, copytree, SameFileError import litebird_sim +from litebird_sim import constants as c from . import HWP from .detectors import DetectorInfo, InstrumentInfo from .distribute import distribute_evenly, distribute_optimally @@ -1129,9 +1130,9 @@ def compute_pointings( def compute_pos_and_vel( self, delta_time_s=86400.0, - solar_velocity_km_s: float = 369.8160, - solar_velocity_gal_lat_rad: float = 0.842_173_724, - solar_velocity_gal_lon_rad: float = 4.608_035_744_4, + solar_velocity_km_s: float = c.SOLAR_VELOCITY_KM_S, + solar_velocity_gal_lat_rad: float = c.SOLAR_VELOCITY_GAL_LAT_RAD, + solar_velocity_gal_lon_rad: float = c.SOLAR_VELOCITY_GAL_LON_RAD, ): """Computes the position and the velocity of the spacescraft for computing the dipole. @@ -1202,7 +1203,7 @@ def fill_tods( def add_dipole( self, - t_cmb_k: float = 2.72548, # Fixsen 2009 http://arxiv.org/abs/0911.1955 + t_cmb_k: float = c.T_CMB_K, dipole_type: DipoleType = DipoleType.TOTAL_FROM_LIN_T, append_to_report: bool = True, ): diff --git a/litebird_sim/spacecraft.py b/litebird_sim/spacecraft.py index 6bf1c95f..bbfef8c0 100644 --- a/litebird_sim/spacecraft.py +++ b/litebird_sim/spacecraft.py @@ -10,11 +10,10 @@ from numba import njit import numpy as np +from litebird_sim import constants as c from .coordinates import DEFAULT_COORDINATE_SYSTEM, DEFAULT_TIME_SCALE from .observations import Observation -EARTH_L2_DISTANCE_KM = 1_496_509.30522 - def cycles_per_year_to_rad_per_s(x: float) -> float: """Convert an angular speed from cycles/yr to rad/s""" @@ -27,7 +26,7 @@ def get_ecliptic_vec(vec): def compute_l2_pos_and_vel( - time0: astropy.time.Time, earth_l2_distance_km: float = EARTH_L2_DISTANCE_KM + time0: astropy.time.Time, earth_l2_distance_km: float = c.EARTH_L2_DISTANCE_KM ): """ Compute the position and velocity of the L2 Sun-Earth point at a given time. @@ -236,15 +235,15 @@ class SpacecraftOrbit: """ start_time: astropy.time.Time - earth_l2_distance_km: float = EARTH_L2_DISTANCE_KM + earth_l2_distance_km: float = c.EARTH_L2_DISTANCE_KM radius1_km: float = 244_450.0 radius2_km: float = 137_388.0 ang_speed1_rad_s: float = cycles_per_year_to_rad_per_s(2.02104) ang_speed2_rad_s: float = cycles_per_year_to_rad_per_s(1.98507) phase_rad: float = np.deg2rad(-47.944) - solar_velocity_km_s: float = 369.8160 - solar_velocity_gal_lat_rad: float = 0.842_173_724 - solar_velocity_gal_lon_rad: float = 4.608_035_744_4 + solar_velocity_km_s: float = c.SOLAR_VELOCITY_KM_S + solar_velocity_gal_lat_rad: float = c.SOLAR_VELOCITY_GAL_LAT_RAD + solar_velocity_gal_lon_rad: float = c.SOLAR_VELOCITY_GAL_LON_RAD solar_velocity_ecl_xyz_km_s = ( SkyCoord(