Bands and dipole in MBS

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:

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

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

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"])
 

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,
     ):

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(