diff --git a/CHANGELOG.md b/CHANGELOG.md index 6a173965..13d31be5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,7 @@ # HEAD +- Add conversion from Temperature to Spectral Radiance units for the non linearity factor [#524](https://github.com/litebird/litebird_sim/pull/524) + - Add complete HWP Jones formalism, including band integration [#499](https://github.com/litebird/litebird_sim/pull/499) - Fix missing parallelizations in multiple places [#521](https://github.com/litebird/litebird_sim/pull/521) diff --git a/litebird_sim/hwp_harmonics/hwp_harmonics.py b/litebird_sim/hwp_harmonics/hwp_harmonics.py index 98af4d35..edbc6222 100644 --- a/litebird_sim/hwp_harmonics/hwp_harmonics.py +++ b/litebird_sim/hwp_harmonics/hwp_harmonics.py @@ -1,21 +1,17 @@ +import logging +import os + import numpy as np import numpy.typing as npt from astropy import constants as const from astropy.cosmology import Planck18 as cosmo from ducc0.healpix import Healpix_Base from numba import njit -import os -import logging from litebird_sim.hwp_jones_parameters import HWPJonesParams -from .jones_methods import ( - compute_signal_for_one_detector as compute_signal_for_one_detector_jones, - integrate_inband_signal_for_one_detector as integrate_inband_signal_for_one_detector_jones, -) -from .mueller_methods import ( - compute_signal_for_one_detector as compute_signal_for_one_detector_mueller, -) + from ..bandpass_template_module import bandpass_profile +from ..constants import NUM_THREADS_ENVVAR from ..coordinates import CoordinateSystem from ..hwp_non_ideal import HWPFormalism, NonIdealHWP from ..input_sky import SkyGenerationParams @@ -24,7 +20,15 @@ from ..pointings_in_obs import ( _get_pointings_array, ) -from ..constants import NUM_THREADS_ENVVAR +from .jones_methods import ( + compute_signal_for_one_detector as compute_signal_for_one_detector_jones, +) +from .jones_methods import ( + integrate_inband_signal_for_one_detector as integrate_inband_signal_for_one_detector_jones, +) +from .mueller_methods import ( + compute_signal_for_one_detector as compute_signal_for_one_detector_mueller, +) COND_THRESHOLD = 1e10 @@ -34,11 +38,11 @@ Tcmb0 = getattr(cosmo, "Tcmb0").value # CMB temperature today in K -def _dBodTrj(nu: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: +def _dBodTrj(nu: npt.NDArray[np.float64] | np.float64): return 2 * k_B * nu * nu * 1e18 / c / c -def _dBodTth(nu: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: +def _dBodTth(nu: npt.NDArray[np.float64] | np.float64): x = h * nu * 1e9 / k_B / Tcmb0 ex = np.exp(x) exm1 = ex - 1.0e0 diff --git a/litebird_sim/non_linearity.py b/litebird_sim/non_linearity.py index dbea2fd1..7f1c04a4 100644 --- a/litebird_sim/non_linearity.py +++ b/litebird_sim/non_linearity.py @@ -1,6 +1,8 @@ -import numpy as np from dataclasses import dataclass +import numpy as np + +from .hwp_harmonics.hwp_harmonics import _dBodTth from .observations import Observation from .seeding import regenerate_or_check_detector_generators @@ -57,8 +59,11 @@ def apply_quadratic_nonlin_for_one_sample( def apply_quadratic_nonlin_for_one_detector( tod_det, + det_bandcenter_ghz: np.float64 | None = None, + det_bandwidth_ghz: np.float64 | None = None, nl_params: NonLinParams | None = None, random: np.random.Generator | None = None, + conv_K_to_SR: bool = False, ): """This function applies the quadratic non-linearity on the TOD corresponding to only one detector. @@ -81,6 +86,16 @@ def apply_quadratic_nonlin_for_one_detector( random (np.random.Generator, optional): A random number generator. Defaults to None. + + conv_K_to_SR (bool, optional): Flag for temperature to spectral radiance + units conversion. Defaults to False. + + det_freq_ghz (np.float64, optional): Detector central frequency in GHz. + Used in the K to SR conversion. Defaults to None. + + det_bandwidth_ghz (np.float64, optional): Detector bandwidth in GHz. + Used in the K to SR conversion. Defaults to None. + """ assert random is not None, ( "You should pass a random number generator which implements the `normal` method." @@ -88,19 +103,30 @@ def apply_quadratic_nonlin_for_one_detector( if nl_params is None: nl_params = NonLinParams() - g_one_over_k = random.normal( + g_nonlin = random.normal( loc=nl_params.sampling_gaussian_loc, scale=nl_params.sampling_gaussian_scale, ) + if conv_K_to_SR: + assert det_bandcenter_ghz is not None and det_bandwidth_ghz is not None, ( + "You should pass det_bandcenter_ghz and det_bandwidth_ghz when conv_K_to_SR is set to True." + ) + conv_factor = _dBodTth(det_bandcenter_ghz) + # convert sampled g (1/K units) to (1/spectral radiance) units + g_nonlin = 1 / (conv_factor * (1 / g_nonlin) * det_bandwidth_ghz * 1e9) + for i in range(len(tod_det)): - tod_det[i] += g_one_over_k * tod_det[i] ** 2 + tod_det[i] += g_nonlin * tod_det[i] ** 2 def apply_quadratic_nonlin( tod: np.ndarray, + bandcenter_ghz: np.ndarray, + bandwidth_ghz: np.ndarray, nl_params: NonLinParams | None = None, dets_random: list[np.random.Generator] | None = None, + conv_K_to_SR: bool = False, ): """Apply a quadratic nonlinearity to some time-ordered data @@ -115,6 +141,9 @@ def apply_quadratic_nonlin( tod_det=tod[detector_idx], nl_params=nl_params, random=dets_random[detector_idx], + conv_K_to_SR=conv_K_to_SR, + det_bandcenter_ghz=bandcenter_ghz[detector_idx], + det_bandwidth_ghz=bandwidth_ghz[detector_idx], ) @@ -124,6 +153,7 @@ def apply_quadratic_nonlin_to_observations( component: str = "tod", user_seed: int | None = None, dets_random: list[np.random.Generator] | None = None, + conv_K_to_SR: bool = False, ): """ Apply a quadratic nonlinearity to some time-ordered data @@ -162,6 +192,8 @@ def apply_quadratic_nonlin_to_observations( List of per-detector random number generators. If not provided, and `user_seed` is given, generators are created internally. One of `user_seed` or `dets_random` must be provided. + conv_K_to_SR (bool, optional): Flag for temperature to spectral radiance + units conversion. Defaults to False. Raises ------ @@ -192,9 +224,14 @@ def apply_quadratic_nonlin_to_observations( # iterate through each observation for cur_obs in obs_list: tod = getattr(cur_obs, component) + bandcenter_ghz = getattr(cur_obs, "bandcenter_ghz") + bandwidth_ghz = getattr(cur_obs, "bandwidth_ghz") apply_quadratic_nonlin( tod=tod, nl_params=nl_params, dets_random=dets_random, + conv_K_to_SR=conv_K_to_SR, + bandcenter_ghz=bandcenter_ghz, + bandwidth_ghz=bandwidth_ghz, ) diff --git a/litebird_sim/simulations.py b/litebird_sim/simulations.py index 98f1e260..6c252140 100644 --- a/litebird_sim/simulations.py +++ b/litebird_sim/simulations.py @@ -2094,6 +2094,7 @@ def apply_quadratic_nonlin( component: str = "tod", append_to_report: bool = False, rng_hierarchy: RNGHierarchy | None = None, + conv_K_to_SR: bool = False, ): """A method to apply non-linearity to the observation. @@ -2101,6 +2102,9 @@ def apply_quadratic_nonlin( :func:`.apply_quadratic_nonlin_to_observations()` that applies non-linearity to a list of :class:`.Observation` instance. Random number generators are obtained from the detector-level layer. As default it uses the `dets_random` field of a :class:`.Simulation` object for this. + + conv_K_to_SR (bool, optional) is a flag for temperature to spectral radiance + units conversion. Defaults to False. """ if rng_hierarchy is None: rng_hierarchy = self.rng_hierarchy @@ -2116,6 +2120,7 @@ def apply_quadratic_nonlin( user_seed=user_seed, component=component, dets_random=dets_random, + conv_K_to_SR=conv_K_to_SR, ) if append_to_report and MPI_COMM_WORLD.rank == 0: @@ -2132,6 +2137,7 @@ def apply_quadratic_nonlin( self.append_to_report( markdown_template, g=g, + conv_K_to_SR=conv_K_to_SR, ) @_profile diff --git a/litebird_sim/templates/report_quad_nonlin.md b/litebird_sim/templates/report_quad_nonlin.md index f1969628..cc9cf03d 100644 --- a/litebird_sim/templates/report_quad_nonlin.md +++ b/litebird_sim/templates/report_quad_nonlin.md @@ -1,5 +1,14 @@ -## Addition of quadratic non-lineariy +# Detector Non Linearity -Quadratic non-linearity has been added to the TOD +{% if g -%} +- Quadratic Non Linearity was applied to the detectors. +- Detector Non-Linearity params: {{ g }} +{% else %} +- Quadratic Non Linearity was NOT applied to the detectors. +{% endif %} -- detector non-linearity factor : {{g}} +{% if conv_K_to_SR -%} +- Non-Linearity g factor was converted from temperature to spectral radiance units. +{% else %} +- Non-Linearity g factor was used in temperature units. +{% endif %}