Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Use power units under the hood #57

Merged
merged 1 commit into from
Mar 5, 2024
Merged

Use power units under the hood #57

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions docs/source/supported-sites.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
###############

===============
Supported sites
###############
===============

+---------------------+--------------------------------------------------------------+------------------+-------------------+--------------------+---------------------+------------------------------------------------------------------------+
| site | description | region | altitude (meters) | latitude (degrees) | longitude (degrees) | documentation |
Expand Down
48 changes: 19 additions & 29 deletions maria/atmosphere/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
from tqdm import tqdm

from .. import utils
from ..constants import k_B
from ..sim.base import BaseSimulation
from ..site import InvalidRegionError, all_regions
from .turbulent_layer import TurbulentLayer
Expand Down Expand Up @@ -84,12 +85,6 @@ def __init__(
from_cache=weather_from_cache,
)

def emission(self, nu):
return

def transmission(self, nu):
return


class AtmosphereMixin:
def _initialize_atmosphere(self):
Expand Down Expand Up @@ -196,32 +191,18 @@ def _simulate_atmospheric_emission(self, units="K_RJ"):
)

for band in bands:
band_index = self.instrument.dets.subset(band=band.name).index
band_index = self.instrument.dets.subset(band_name=band.name).index

if self.verbose:
bands.set_description(f"Computing atm. emission ({band.name})")

# multiply by 1e9 to go from GHz to Hz
# det_power_grid = (
# 1e9
# * 1.380649e-23
# * np.trapz(
# self.atmosphere.emission_spectrum
# * band.passband(self.atmosphere.spectrum_side_nu),
# self.atmosphere.spectrum_side_nu,
# axis=-1,
# )
# )

# this is NOT power
det_power_grid = np.sum(
# in watts. the 1e9 is for GHz -> Hz
det_power_grid = k_B * np.trapz(
self.atmosphere.emission_spectrum
* band.passband(self.atmosphere.spectrum_side_nu),
1e9 * self.atmosphere.spectrum_side_nu,
axis=-1,
)
det_power_grid /= np.sum(
band.passband(self.atmosphere.spectrum_side_nu), axis=-1
)

band_power_interpolator = sp.interpolate.RegularGridInterpolator(
(
Expand All @@ -240,14 +221,19 @@ def _simulate_atmospheric_emission(self, units="K_RJ"):
)
)

self.atmospheric_transmission = np.empty(
(self.instrument.n_dets, self.plan.n_time), dtype=np.float32
)

# to make a new progress bar
bands = (
tqdm(self.instrument.dets.bands)
if self.verbose
else self.instrument.dets.bands
)

for band in bands:
band_index = self.instrument.dets.subset(band=band.name).index
band_index = self.instrument.dets.subset(band_name=band.name).index

if self.verbose:
bands.set_description(f"Computing atm. transmission ({band.name})")
Expand All @@ -257,14 +243,13 @@ def _simulate_atmospheric_emission(self, units="K_RJ"):
* self.atmosphere.spectrum_side_nu**2
)

# multiply by 1e9 to go from GHz to Hz
self.det_transmission_grid = np.trapz(
rel_T_RJ_spectrum * self.atmosphere.transmission_spectrum,
self.atmosphere.spectrum_side_nu,
1e9 * self.atmosphere.spectrum_side_nu,
axis=-1,
) / np.trapz(
rel_T_RJ_spectrum,
self.atmosphere.spectrum_side_nu,
1e9 * self.atmosphere.spectrum_side_nu,
axis=-1,
)

Expand All @@ -277,7 +262,12 @@ def _simulate_atmospheric_emission(self, units="K_RJ"):
self.det_transmission_grid,
)

self.atmospheric_transmission = band_transmission_interpolator(
# what's happening here? the atmosphere blocks some of the light from space.
# we want to calibrate to the stuff in space, so we make the atmosphere *hotter*

self.atmospheric_transmission[
band_index
] = band_transmission_interpolator(
(
self.zenith_scaled_pwv[band_index],
self.atmosphere.weather.temperature[0],
Expand Down
2 changes: 1 addition & 1 deletion maria/atmosphere/turbulent_layer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -312,7 +312,7 @@ def sample(self):
for band in self.instrument.bands:
# we assume the atmosphere looks the same for every nu in the band

band_index = self.instrument.dets.subset(band=band.name).index
band_index = self.instrument.dets.subset(band_name=band.name).index
band_angular_fwhm = self.instrument.angular_fwhm(z=self.depth)[
band_index
].mean()
Expand Down
2 changes: 1 addition & 1 deletion maria/atmosphere/weather.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@
import pytz
import scipy as sp

from ..constants import g
from ..site import InvalidRegionError, all_regions, supported_regions_table
from ..utils import get_utc_day_hour, get_utc_year_day
from ..utils.constants import g
from . import utils

here, this_filename = os.path.split(__file__)
Expand Down
2 changes: 2 additions & 0 deletions maria/utils/constants.py → maria/constants.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,3 +8,5 @@
k_B = 1.380649e-23

T_CMB = 2.7255

c = 299_792_458.0
0 maria/utils/errors.py → maria/errors.py
View file
Open in desktop
File renamed without changes.
4 changes: 2 additions & 2 deletions maria/instrument/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@ def __repr__(self):

@property
def ubands(self):
return self.dets.ubands
return self.dets.bands.names

@property
def offset_x(self):
Expand Down Expand Up @@ -178,7 +178,7 @@ def plot_dets(self, units="deg"):

legend_handles = []
for iub, uband in enumerate(self.ubands):
band_mask = self.dets.band == uband
band_mask = self.dets.band_name == uband

fwhm = np.degrees(self.fwhm[band_mask])
offsets = np.degrees(self.offsets[band_mask])
Expand Down
38 changes: 23 additions & 15 deletions maria/instrument/bands.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,14 @@
import pandas as pd

from .. import utils
from ..constants import k_B
from ..utils.io import flatten_config, read_yaml

BAND_FIELD_TYPES = {
"center": "float",
"width": "float",
"shape": "str",
"tau": "float",
"time_constant": "float",
"white_noise": "float",
"pink_noise": "float",
"efficiency": "float",
Expand Down Expand Up @@ -102,18 +103,16 @@ class Band:
name: str
center: float
width: float
tau: float = 0.0
time_constant: float = 0.0
white_noise: float = 0.0
pink_noise: float = 0.0
shape: str = "gaussian"
shape: str = "top_hat"
efficiency: float = 1.0

@classmethod
def from_passband(cls, name, nu, pb, pb_err=None):
center = np.round(np.sum(pb * nu), 3)
width = np.round(
nu[pb > pb.max() / np.e**2].ptp(), 3
) # width is the two-sigma interval
width = np.round(nu[pb > 1e-2 * pb.max()].ptp(), 3)

band = cls(name=name, center=center, width=width, shape="custom")

Expand All @@ -126,34 +125,43 @@ def from_passband(cls, name, nu, pb, pb_err=None):
def nu_min(self) -> float:
if self.shape == "flat":
return self.center - 0.5 * self.width
if self.shape == "gaussian":
return self.center - self.width
if self.shape == "custom":
return self._nu[self._pb > 1e-2 * self._pb.max()].min()

return self.center - self.width

@property
def nu_max(self) -> float:
if self.shape == "flat":
return self.center + 0.5 * self.width
if self.shape == "gaussian":
return self.center + self.width
if self.shape == "custom":
return self._nu[self._pb > 1e-2 * self._pb.max()].max()

return self.center + self.width

def passband(self, nu):
"""
Passband response as a function of nu (in GHz). These integrate to one.
Passband response as a function of nu (in GHz).
"""

_nu = np.atleast_1d(nu)

if self.shape == "gaussian":
band_sigma = self.width / 4
if self.shape == "top_hat":
return np.exp(-np.log(2) * (2 * (_nu - self.center) / self.width) ** 8)

return np.exp(-0.5 * np.square((_nu - self.center) / band_sigma))
if self.shape == "gaussian":
return np.exp(-np.log(2) * (2 * (_nu - self.center) / self.width) ** 2)

if self.shape == "flat":
return np.where((_nu > self.nu_min) & (_nu < self.nu_max), 1.0, 0.0)

elif self.shape == "custom":
return np.interp(_nu, self._nu, self._pb)

@property
def abs_cal_rj(self):
"""
Absolute calibration (i.e. temperature per power) in Kelvins per Watt.
"""
nu = np.linspace(self.nu_min, self.nu_max, 256)
dP_dT = self.efficiency * np.trapz(self.passband(nu) * k_B, 1e9 * nu)
return 1 / dP_dT
2 changes: 1 addition & 1 deletion maria/instrument/bands/abs.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
f150:
center: 150
width: 30
tau: 0
time_constant: 0
white_noise: 266.e-6 # in K_RJ s^-1/2
pink_noise: 1.e-1 # in K_RJ s^-1/2
efficiency: 0.5
12 changes: 6 additions & 6 deletions maria/instrument/bands/act.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,15 +3,15 @@ pa4:
center: 150
width: 30
shape: gaussian
tau: 0
time_constant: 0
white_noise: 266.e-6
pink_noise: 1.e-1
efficiency: 0.5
f220:
center: 220
width: 30
shape: gaussian
tau: 0
time_constant: 0
white_noise: 266.e-6
pink_noise: 1.e-1
efficiency: 0.5
Expand All @@ -20,15 +20,15 @@ pa5:
center: 90
width: 30
shape: gaussian
tau: 0
time_constant: 0
white_noise: 266.e-6
pink_noise: 1.e-1
efficiency: 0.5
f150:
center: 150
width: 30
shape: gaussian
tau: 0
time_constant: 0
white_noise: 266.e-6
pink_noise: 1.e-1
efficiency: 0.5
Expand All @@ -37,15 +37,15 @@ pa6:
center: 90
width: 30
shape: gaussian
tau: 0
time_constant: 0
white_noise: 266.e-6
pink_noise: 1.e-1
efficiency: 0.5
f150:
center: 150
width: 30
shape: gaussian
tau: 0
time_constant: 0
white_noise: 266.e-6
pink_noise: 1.e-1
efficiency: 0.5
Expand Down
Loading
Loading