esbo_ds/ESBO-ETC
esbo_ds/ESBO-ETC
11
0
Fork 1
Code Issues 1 Pull Requests Releases 1 Wiki Activity

Heterodyne output added.

Browse Source
This commit is contained in:
Lukas Klass 2020-07-14 12:01:28 +02:00
parent 367568e0d3
commit fa911d2b19
3 changed files with 142 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
docs/source/developer/images/Sensor.pdf
View File

Binary file not shown.

BIN
docs/source/developer/images/class_diagram.pdf
View File

Binary file not shown.

186
esbo_etc/classes/sensor/Heterodyne.py
View File

@ -1,12 +1,14 @@
from astropy import units as u
from .ASensor import ASensor from .ASensor import ASensor
from ..IRadiant import IRadiant from ..IRadiant import IRadiant
from ..Entry import Entry from ..Entry import Entry
import numpy as np
from astropy.constants import k_B
from typing import Union
from ...lib.logger import logger from ...lib.logger import logger
from ..SpectralQty import SpectralQty from ..SpectralQty import SpectralQty
import numpy as np
from astropy import units as u
from astropy.constants import k_B
from astropy.table import QTable
from typing import Union
import os
class Heterodyne(ASensor): class Heterodyne(ASensor):
@ -75,24 +77,25 @@ class Heterodyne(ASensor):
""" """
# Calculate the signal and background temperatures # Calculate the signal and background temperatures
t_signal, t_background = self.calcTemperatures(background, signal, obstruction) t_signal, t_background = self.calcTemperatures(background, signal, obstruction)
line_ind = np.where(t_signal.wl == self.__lambda_line)[0][0]
t_sys = 2 * (t_background + self.__receiver_temp) t_sys = 2 * (t_background + self.__receiver_temp)
# Calculate the noise bandwidth # Calculate the noise bandwidth
delta_nu = self.__lambda_line.to(u.Hz, equivalencies=u.spectral()) / ( delta_nu = t_signal.wl.to(u.Hz, equivalencies=u.spectral()) / (t_signal.wl / self.__common_conf.wl_delta() + 1)
self.__lambda_line / self.__common_conf.wl_delta() + 1) snr = []
for exp_time_ in exp_time if exp_time.size > 1 else [exp_time]:
# Calculate the RMS background temperature # Calculate the RMS background temperature
if self.__n_on is None: if self.__n_on is None:
t_rms = 2 * t_sys * self.__kappa / np.sqrt(exp_time * delta_nu) t_rms = 2 * t_sys * self.__kappa / np.sqrt(exp_time_ * delta_nu)
else: else:
t_rms = t_sys * self.__kappa * np.sqrt(1 + 1 / np.sqrt(self.__n_on)) / np.sqrt(exp_time * delta_nu) t_rms = t_sys * self.__kappa * np.sqrt(1 + 1 / np.sqrt(self.__n_on)) / np.sqrt(exp_time_ * delta_nu)
# Calculate the SNR # Calculate the SNR
snr = t_signal / t_rms snr_ = t_signal / t_rms
snr.append(snr_.qty[line_ind])
# Print details # Print details
if exp_time.size > 1: self.__printDetails(t_sys.qty[line_ind], delta_nu[line_ind], t_rms.qty[line_ind], t_signal.qty[line_ind],
for i in range(exp_time.size): "t_exp=%.2f s: " % exp_time_.value)
self.__printDetails(t_sys, delta_nu, t_rms[i], t_signal, "t_exp=%.2f s: " % exp_time[i].value) self.__output(t_signal, t_background, t_rms, "texp_%.2f" % exp_time_.value, snr=snr_)
else: return u.Quantity(snr) if len(snr) > 1 else u.Quantity(snr[0])
self.__printDetails(t_sys, delta_nu, t_rms, t_signal, "t_exp=%.2f s: " % exp_time.value)
return snr
@u.quantity_input(snr=u.dimensionless_unscaled) @u.quantity_input(snr=u.dimensionless_unscaled)
def calcExpTime(self, background: SpectralQty, signal: SpectralQty, obstruction: float, snr: u.Quantity) -> u.s: def calcExpTime(self, background: SpectralQty, signal: SpectralQty, obstruction: float, snr: u.Quantity) -> u.s:
@ -117,28 +120,32 @@ class Heterodyne(ASensor):
""" """
# Calculate the signal and background temperatures # Calculate the signal and background temperatures
t_signal, t_background = self.calcTemperatures(background, signal, obstruction) t_signal, t_background = self.calcTemperatures(background, signal, obstruction)
line_ind = np.where(t_signal.wl == self.__lambda_line)[0][0]
t_sys = 2 * (t_background + self.__receiver_temp) t_sys = 2 * (t_background + self.__receiver_temp)
# Calculate the noise bandwidth # Calculate the noise bandwidth
delta_nu = self.__lambda_line.to(u.Hz, equivalencies=u.spectral()) / ( delta_nu = t_signal.wl.to(u.Hz, equivalencies=u.spectral()) / (t_signal.wl / self.__common_conf.wl_delta() + 1)
self.__lambda_line / self.__common_conf.wl_delta() + 1) exp_time = []
for snr_ in snr if snr.size > 1 else [snr]:
# Calculate the RMS background temperature # Calculate the RMS background temperature
t_rms = t_signal / snr t_rms = t_signal / snr_
# Calculate the exposure time # Calculate the exposure time
if self.__n_on is None: if self.__n_on is None:
exp_time = ((2 * t_sys * self.__kappa / t_rms) ** 2 / delta_nu).decompose() exp_time_ = ((2 * t_sys * self.__kappa / t_rms) ** 2 / delta_nu)
else: else:
exp_time = ((t_sys * self.__kappa / t_rms) ** 2 * (1 + 1 / np.sqrt(self.__n_on)) / delta_nu).decompose() exp_time_ = ((t_sys * self.__kappa / t_rms) ** 2 *
(1 + 1 / np.sqrt(self.__n_on)) / delta_nu)
exp_time_ = SpectralQty(exp_time_.wl, exp_time_.qty.decompose())
exp_time.append(exp_time_.qty[line_ind])
# Print details # Print details
if snr.size > 1: self.__printDetails(t_sys.qty[line_ind], delta_nu[line_ind], t_rms.qty[line_ind], t_signal.qty[line_ind],
for i in range(snr.size): "SNR=%.2f: " % snr_.value)
self.__printDetails(t_sys, delta_nu, t_rms[i], t_signal, "SNR=%.2f: " % snr[i].value) self.__output(t_signal, t_background, t_rms, "snr_%.2f" % snr_.value, exp_time=exp_time_)
else: return u.Quantity(exp_time) if len(exp_time) > 1 else u.Quantity(exp_time[0])
self.__printDetails(t_sys, delta_nu, t_rms, t_signal, "SNR=%.2f: " % snr.value)
return exp_time
@u.quantity_input(exp_time="time", snr=u.dimensionless_unscaled, target_brightness=u.mag) # @u.quantity_input(exp_time="time", snr=u.dimensionless_unscaled,
# target_brightness=[u.mag, u.mag / u.sr])
def calcSensitivity(self, background: SpectralQty, signal: SpectralQty, obstruction: float, exp_time: u.Quantity, def calcSensitivity(self, background: SpectralQty, signal: SpectralQty, obstruction: float, exp_time: u.Quantity,
snr: u.Quantity, target_brightness: u.Quantity) -> u.mag: snr: u.Quantity, target_brightness: u.Quantity) -> [u.mag, u.mag / u.sr]:
""" """
Calculate the sensitivity of the telescope detector combination. Calculate the sensitivity of the telescope detector combination.
@ -155,7 +162,7 @@ class Heterodyne(ASensor):
snr : Quantity snr : Quantity
The SNR for which the sensitivity time shall be calculated. The SNR for which the sensitivity time shall be calculated.
target_brightness : Quantity target_brightness : Quantity
The target brightness in magnitudes. The target brightness in mag or mag / sr.
Returns Returns
------- -------
@ -164,26 +171,30 @@ class Heterodyne(ASensor):
""" """
# Calculate the signal and background temperatures # Calculate the signal and background temperatures
t_signal, t_background = self.calcTemperatures(background, signal, obstruction) t_signal, t_background = self.calcTemperatures(background, signal, obstruction)
line_ind = np.where(t_signal.wl == self.__lambda_line)[0][0]
t_sys = 2 * (t_background + self.__receiver_temp) t_sys = 2 * (t_background + self.__receiver_temp)
# Calculate the noise bandwidth # Calculate the noise bandwidth
delta_nu = self.__lambda_line.to(u.Hz, equivalencies=u.spectral()) / ( delta_nu = t_signal.wl.to(u.Hz, equivalencies=u.spectral()) / (t_signal.wl / self.__common_conf.wl_delta() + 1)
self.__lambda_line / self.__common_conf.wl_delta() + 1) sensitivity = []
for snr_, exp_time_ in zip(snr, exp_time) if snr.size > 1 else zip([snr], [exp_time]):
# Calculate the RMS background temperature # Calculate the RMS background temperature
if self.__n_on is None: if self.__n_on is None:
t_rms = 2 * t_sys * self.__kappa / np.sqrt(exp_time * delta_nu) t_rms = 2 * t_sys * self.__kappa / np.sqrt(exp_time_ * delta_nu)
else: else:
t_rms = t_sys * self.__kappa * np.sqrt(1 + 1 / np.sqrt(self.__n_on)) / np.sqrt(exp_time * delta_nu) t_rms = t_sys * self.__kappa * np.sqrt(1 + 1 / np.sqrt(self.__n_on)) / np.sqrt(exp_time_ * delta_nu)
# Calculate the limiting signal temperature # Calculate the limiting signal temperature
t_signal_lim = t_rms * snr t_signal_lim = t_rms * snr_
# Calculate the sensitivity
signal_ratio = t_signal_lim / t_signal
sensitivity_ = SpectralQty(signal_ratio.wl,
target_brightness - 2.5 * np.log10(signal_ratio.qty) * target_brightness.unit)
sensitivity.append(sensitivity_.qty[line_ind])
# Print details # Print details
if snr.size > 1: self.__printDetails(t_sys.qty[line_ind], delta_nu[line_ind], t_rms.qty[line_ind],
for i in range(snr.size): t_signal_lim.qty[line_ind], "SNR=%.2f t_exp=%.2f s: " % (snr_.value, exp_time_.value))
self.__printDetails(t_sys, delta_nu, t_rms[i], t_signal_lim[i], self.__output(t_signal, t_background, t_rms, "snr_%.2f_texp_%.2f" % (snr_.value, exp_time_.value),
"SNR=%.2f t_exp=%.2f s: " % (snr[i].value, exp_time[i].value)) sensitivity=sensitivity_)
else: return u.Quantity(sensitivity) if len(sensitivity) > 1 else u.Quantity(sensitivity[0])
self.__printDetails(t_sys, delta_nu, t_rms, t_signal_lim,
"SNR=%.2f t_exp=%.2f s: " % (snr.value, exp_time.value))
return target_brightness - 2.5 * np.log10(t_signal_lim / t_signal) * u.mag
@u.quantity_input(t_sys=u.K, delta_nu=u.Hz, t_rms=u.K, t_signal=u.K) @u.quantity_input(t_sys=u.K, delta_nu=u.Hz, t_rms=u.K, t_signal=u.K)
def __printDetails(self, t_sys: u.Quantity, delta_nu: u.Quantity, t_rms: u.Quantity, def __printDetails(self, t_sys: u.Quantity, delta_nu: u.Quantity, t_rms: u.Quantity,
@ -214,6 +225,53 @@ class Heterodyne(ASensor):
logger.info(prefix + "Antenna temperature: %1.2e K" % t_signal.value) logger.info(prefix + "Antenna temperature: %1.2e K" % t_signal.value)
logger.info("-------------------------------------------------------------------------------------------------") logger.info("-------------------------------------------------------------------------------------------------")
@u.quantity_input(signal=u.electron, background=u.electron, read_noise=u.electron ** 0.5, dark=u.electron)
def __output(self, t_signal: SpectralQty, t_background: SpectralQty, t_rms: SpectralQty,
name: str, snr: SpectralQty = None, exp_time: SpectralQty = None, sensitivity: SpectralQty = None):
"""
Write the signal and the noise in electrons to files.
Parameters
----------
t_signal : SpectralQty
The signal temperature in Kelvins.
t_background : SpectralQty
The background temperature in Kelvins.
t_rms : SpectralQty
The RMS noise temperature in Kelvins.
name : str
The name of the configuration.
snr : SpectralQty
The calculated signal-to-noise ratio per wavelength.
exp_time : SpectralQty
The calculated exposure time per wavelength.
sensitivity : SpectralQty
The calculated sensitivity per wavelength.
Returns
-------
"""
# Concatenate the paths
path = os.path.join(self.__common_conf.output.path, name)
try:
os.makedirs(path, exist_ok=True)
except FileExistsError:
logger.warning("Output directory '" + path + "' already exists.")
res = QTable([t_signal.wl, t_signal.qty, t_background.qty, t_rms.qty],
names=('Wavelength [' + t_signal.wl.unit.to_string() + ']',
'Signal Temperature [' + t_signal.qty.unit.to_string() + ']',
'Background Temperature [' + t_background.qty.unit.to_string() + ']',
'RMS Noise Temperature [' + t_rms.qty.unit.to_string() + ']'),
meta={'name': 'first table'})
if snr is not None:
res['SNR [-]'] = snr.qty
if exp_time is not None:
res['Exposure Time [' + exp_time.qty.unit.to_string() + ']'] = exp_time.qty
if sensitivity is not None:
res['Sensitivity [' + sensitivity.qty.unit.to_string() + ']'] = sensitivity.qty
res.write(os.path.join(path, "result.csv"), format='ascii.csv', overwrite=True)
def calcTemperatures(self, background: SpectralQty, signal: SpectralQty, obstruction: float): def calcTemperatures(self, background: SpectralQty, signal: SpectralQty, obstruction: float):
""" """
Calculate the noise temperatures of the signal and the background radiation. Calculate the noise temperatures of the signal and the background radiation.
@ -229,32 +287,42 @@ class Heterodyne(ASensor):
Returns Returns
------- -------
t_signal : u.Quantity t_signal : SpectralQty
The signal temperature in Kelvins. The spectral signal temperature in Kelvins.
t_background : u.Quantity t_background : SpectralQty
The background temperature in Kelvins. The spectral signal temperature in Kelvins.
""" """
logger.info("Calculating the system temperature.") logger.info("Calculating the system temperature.")
t_background = (background.rebin(self.__lambda_line).qty.to( # Add desired wavelength to wavelength bins
u.W / (u.m ** 2 * u.Hz * u.sr), equivalencies=u.spectral_density(self.__lambda_line)) * wl_bins = np.sort(np.append(self.__common_conf.wl_bins(), self.__lambda_line)).view(u.Quantity)
self.__lambda_line ** 2 / (2 * k_B) * u.sr).decompose() signal = signal.rebin(wl_bins)
background = background.rebin(wl_bins)
t_background = SpectralQty(background.wl, background.qty.to(u.W / (u.m ** 2 * u.Hz * u.sr),
equivalencies=u.spectral_density(
background.wl)))
t_background = t_background * (t_background.wl ** 2 / (2 * k_B) * u.sr)
t_background = SpectralQty(t_background.wl, t_background.qty.decompose())
# Calculate the incoming photon current of the target # Calculate the incoming photon current of the target
logger.info("Calculating the signal temperature.") logger.info("Calculating the signal temperature.")
size = "extended" if signal.qty.unit.is_equivalent(u.W / (u.m ** 2 * u.nm * u.sr)) else "point" size = "extended" if signal.qty.unit.is_equivalent(u.W / (u.m ** 2 * u.nm * u.sr)) else "point"
if size == "point": if size == "point":
signal = signal.rebin(self.__lambda_line).qty.to(u.W / (u.m ** 2 * u.Hz), signal = SpectralQty(signal.wl, signal.qty.to(u.W / (u.m ** 2 * u.Hz),
equivalencies=u.spectral_density(self.__lambda_line)) equivalencies=u.spectral_density(signal.wl)))
t_signal = (signal * self.__aperture_efficiency * self.__common_conf.d_aperture() ** 2 * t_signal = signal * (self.__aperture_efficiency * self.__common_conf.d_aperture() ** 2 *
np.pi / 4 / (2 * k_B) * self.__eta_fss).decompose() np.pi / 4 / (2 * k_B) * self.__eta_fss)
t_signal = SpectralQty(t_signal.wl, t_signal.qty.decompose())
else: else:
signal = signal.rebin(self.__lambda_line).qty.to(u.W / (u.m ** 2 * u.Hz * u.sr), signal = SpectralQty(signal.wl, signal.qty.to(u.W / (u.m ** 2 * u.Hz * u.sr),
equivalencies=u.spectral_density(self.__lambda_line)) equivalencies=u.spectral_density(signal.wl)))
t_signal = (signal * u.sr * self.__main_beam_efficiency * self.__lambda_line ** 2 / ( t_signal = signal * (u.sr * self.__main_beam_efficiency * self.__lambda_line ** 2 / (
2 * k_B) * self.__eta_fss).decompose() 2 * k_B) * self.__eta_fss)
logger.debug("Signal temperature: %1.2e K" % t_signal.value) t_signal = SpectralQty(t_signal.wl, t_signal.qty.decompose())
logger.debug("Spectral signal temperature")
logger.debug(t_signal)
logger.debug("Target size: " + size) logger.debug("Target size: " + size)
logger.debug("Obstruction: %.2f" % obstruction) logger.debug("Obstruction: %.2f" % obstruction)
logger.debug("Background temperature: %1.2e K" % t_background.value) logger.debug("Spectral background temperature")
logger.debug(t_background)
return t_signal, t_background return t_signal, t_background
@staticmethod @staticmethod