# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: healpyvenv
#     language: python
#     name: healpyvenv
# ---

# # Unit conversion with broadband detectors looking at the CMB
# - categories: [cmbs4, pysm, simonsobs]

import healpy as hp
import numpy as np
import matplotlib.pyplot as plt
from pysm3 import units as u
import pysm3 as pysm
# %matplotlib inline

dip = hp.synfast([0,1], lmax=1, nside=128) * u.V

hp.mollview(dip, unit=dip.unit)

# We measure the sky with out broadband instrument, we assume we only measure the CMB solar dipole,
# initially the units are arbitrary, for example Volts of our instrument.
#
# Next we calibrate on the solar dipole, which is known to be 3.3 mK.

calibration_factor = 2 * 3.3 * u.mK_CMB / (dip.max() - dip.min())

calibration_factor

calibrated_dip = calibration_factor * dip

calibrated_dip

hp.mollview(calibrated_dip, unit=calibrated_dip.unit)

# First we simplify and consider a delta-frequency instrument at 300 GHz

center_frequency = 300 * u.GHz

dip_peak = calibrated_dip.max()

calibrated_dip.max()

calibrated_dip.max().to(u.mK_RJ, equivalencies=u.cmb_equivalencies(center_frequency))

calibrated_dip.max().to(u.MJy/u.sr, equivalencies=u.cmb_equivalencies(center_frequency))

# Next we assume instead that we have a broadband instrument, of 20% bandwidth, with uniform response in that range.
# For simplicity, we only take 4 points.

freq = [270, 290, 310, 330] * u.GHz

weights = [1, 1, 1, 1]

weights /= np.trapz(weights, freq)

weights

# The instrument bandpass is defined in power so we can transform our signal in MJy/sr at the 4 reference frequencies,
# then integrate.

dip_peak_MJysr = dip_peak.to(u.MJy/u.sr, equivalencies=u.cmb_equivalencies(freq))

dip_peak_MJysr

integrated_SR = np.trapz(dip_peak_MJysr * weights, freq)

integrated_SR

# This is different than assuming uniform bandpass in $K_{CMB}$, where instead we would recover the same result of the delta-bandpass:

np.trapz(dip_peak * weights, freq)

SR = u.MJy/u.sr

# We use the PySM 3 function to compute unit conversion given a bandpass
#
# $ \tilde{I}[unit_{out}] = \tilde{I}[unit_{in}] \dfrac{ \int C_{unit_{in}}^{Jy~sr^{-1}}(\nu) g(\nu) d\nu}  { \int C_{unit_{out}}^{Jy~sr^{-1}}(\nu) g(\nu) d\nu} $
#
# which comes from equating in power:
#
# $ \tilde{I}[unit_{out}]{ \int C_{unit_{out}}^{Jy~sr^{-1}}(\nu) g(\nu) d\nu} = \tilde{I}[unit_{in}]\int C_{unit_{in}}^{Jy~sr^{-1}}(\nu) g(\nu) d\nu   $

SR

pysm.utils.bandpass_unit_conversion(freq, weights=weights, output_unit=u.mK_CMB, input_unit=SR)

integrated_SR * _

1 * u.mK_CMB / (1 * SR)

# We can doublecheck the implementation of the PySM function by executing it here:

K_CMB_to_MJysr = ((1*SR).to(u.mK_CMB, equivalencies=u.cmb_equivalencies(freq)))/(1*SR)

K_CMB_to_MJysr

# Integrating the `K_CMB_to_MJysr` conversion factor is wrong, we always need to do the integral in power,
# therefore we integrate the inverse and then take its inverse.

np.trapz(K_CMB_to_MJysr * weights, freq)

1/np.trapz(1/K_CMB_to_MJysr * weights, freq)



StarcoderdataJupyterScriptsDedupFiltered