diff --git a/docs/source/developer/images/Sensor.pdf b/docs/source/developer/images/Sensor.pdf index c9592e7..e4c1092 100644 Binary files a/docs/source/developer/images/Sensor.pdf and b/docs/source/developer/images/Sensor.pdf differ diff --git a/docs/source/developer/images/class_diagram.pdf b/docs/source/developer/images/class_diagram.pdf index dbfa869..080b711 100644 Binary files a/docs/source/developer/images/class_diagram.pdf and b/docs/source/developer/images/class_diagram.pdf differ diff --git a/esbo_etc/classes/sensor/Heterodyne.py b/esbo_etc/classes/sensor/Heterodyne.py index c03775e..baa4aef 100644 --- a/esbo_etc/classes/sensor/Heterodyne.py +++ b/esbo_etc/classes/sensor/Heterodyne.py @@ -1,12 +1,14 @@ -from astropy import units as u from .ASensor import ASensor from ..IRadiant import IRadiant from ..Entry import Entry -import numpy as np -from astropy.constants import k_B -from typing import Union from ...lib.logger import logger 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): @@ -75,24 +77,25 @@ class Heterodyne(ASensor): """ # Calculate the signal and background temperatures 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) # Calculate the noise bandwidth - delta_nu = self.__lambda_line.to(u.Hz, equivalencies=u.spectral()) / ( - self.__lambda_line / self.__common_conf.wl_delta() + 1) - # Calculate the RMS background temperature - if self.__n_on is None: - t_rms = 2 * t_sys * self.__kappa / np.sqrt(exp_time * delta_nu) - else: - t_rms = t_sys * self.__kappa * np.sqrt(1 + 1 / np.sqrt(self.__n_on)) / np.sqrt(exp_time * delta_nu) - # Calculate the SNR - snr = t_signal / t_rms - # Print details - if exp_time.size > 1: - for i in range(exp_time.size): - self.__printDetails(t_sys, delta_nu, t_rms[i], t_signal, "t_exp=%.2f s: " % exp_time[i].value) - else: - self.__printDetails(t_sys, delta_nu, t_rms, t_signal, "t_exp=%.2f s: " % exp_time.value) - return snr + delta_nu = t_signal.wl.to(u.Hz, equivalencies=u.spectral()) / (t_signal.wl / 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 + if self.__n_on is None: + t_rms = 2 * t_sys * self.__kappa / np.sqrt(exp_time_ * delta_nu) + else: + t_rms = t_sys * self.__kappa * np.sqrt(1 + 1 / np.sqrt(self.__n_on)) / np.sqrt(exp_time_ * delta_nu) + # Calculate the SNR + snr_ = t_signal / t_rms + snr.append(snr_.qty[line_ind]) + # Print details + self.__printDetails(t_sys.qty[line_ind], delta_nu[line_ind], t_rms.qty[line_ind], t_signal.qty[line_ind], + "t_exp=%.2f s: " % exp_time_.value) + self.__output(t_signal, t_background, t_rms, "texp_%.2f" % exp_time_.value, snr=snr_) + return u.Quantity(snr) if len(snr) > 1 else u.Quantity(snr[0]) @u.quantity_input(snr=u.dimensionless_unscaled) 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 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) # Calculate the noise bandwidth - delta_nu = self.__lambda_line.to(u.Hz, equivalencies=u.spectral()) / ( - self.__lambda_line / self.__common_conf.wl_delta() + 1) - # Calculate the RMS background temperature - t_rms = t_signal / snr - # Calculate the exposure time - if self.__n_on is None: - exp_time = ((2 * t_sys * self.__kappa / t_rms) ** 2 / delta_nu).decompose() - else: - exp_time = ((t_sys * self.__kappa / t_rms) ** 2 * (1 + 1 / np.sqrt(self.__n_on)) / delta_nu).decompose() - # Print details - if snr.size > 1: - for i in range(snr.size): - self.__printDetails(t_sys, delta_nu, t_rms[i], t_signal, "SNR=%.2f: " % snr[i].value) - else: - self.__printDetails(t_sys, delta_nu, t_rms, t_signal, "SNR=%.2f: " % snr.value) - return exp_time + delta_nu = t_signal.wl.to(u.Hz, equivalencies=u.spectral()) / (t_signal.wl / self.__common_conf.wl_delta() + 1) + exp_time = [] + for snr_ in snr if snr.size > 1 else [snr]: + # Calculate the RMS background temperature + t_rms = t_signal / snr_ + # Calculate the exposure time + if self.__n_on is None: + exp_time_ = ((2 * t_sys * self.__kappa / t_rms) ** 2 / delta_nu) + else: + 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 + self.__printDetails(t_sys.qty[line_ind], delta_nu[line_ind], t_rms.qty[line_ind], t_signal.qty[line_ind], + "SNR=%.2f: " % snr_.value) + self.__output(t_signal, t_background, t_rms, "snr_%.2f" % snr_.value, exp_time=exp_time_) + return u.Quantity(exp_time) if len(exp_time) > 1 else u.Quantity(exp_time[0]) - @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, - 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. @@ -155,7 +162,7 @@ class Heterodyne(ASensor): snr : Quantity The SNR for which the sensitivity time shall be calculated. target_brightness : Quantity - The target brightness in magnitudes. + The target brightness in mag or mag / sr. Returns ------- @@ -164,26 +171,30 @@ class Heterodyne(ASensor): """ # Calculate the signal and background temperatures 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) # Calculate the noise bandwidth - delta_nu = self.__lambda_line.to(u.Hz, equivalencies=u.spectral()) / ( - self.__lambda_line / self.__common_conf.wl_delta() + 1) - # Calculate the RMS background temperature - if self.__n_on is None: - t_rms = 2 * t_sys * self.__kappa / np.sqrt(exp_time * delta_nu) - else: - 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 - t_signal_lim = t_rms * snr - # Print details - if snr.size > 1: - for i in range(snr.size): - self.__printDetails(t_sys, delta_nu, t_rms[i], t_signal_lim[i], - "SNR=%.2f t_exp=%.2f s: " % (snr[i].value, exp_time[i].value)) - else: - 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 + delta_nu = t_signal.wl.to(u.Hz, equivalencies=u.spectral()) / (t_signal.wl / 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 + if self.__n_on is None: + t_rms = 2 * t_sys * self.__kappa / np.sqrt(exp_time_ * delta_nu) + else: + 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 + 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 + self.__printDetails(t_sys.qty[line_ind], delta_nu[line_ind], t_rms.qty[line_ind], + t_signal_lim.qty[line_ind], "SNR=%.2f t_exp=%.2f s: " % (snr_.value, exp_time_.value)) + self.__output(t_signal, t_background, t_rms, "snr_%.2f_texp_%.2f" % (snr_.value, exp_time_.value), + sensitivity=sensitivity_) + return u.Quantity(sensitivity) if len(sensitivity) > 1 else u.Quantity(sensitivity[0]) @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, @@ -214,6 +225,53 @@ class Heterodyne(ASensor): logger.info(prefix + "Antenna temperature: %1.2e K" % t_signal.value) 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): """ Calculate the noise temperatures of the signal and the background radiation. @@ -229,32 +287,42 @@ class Heterodyne(ASensor): Returns ------- - t_signal : u.Quantity - The signal temperature in Kelvins. - t_background : u.Quantity - The background temperature in Kelvins. + t_signal : SpectralQty + The spectral signal temperature in Kelvins. + t_background : SpectralQty + The spectral signal temperature in Kelvins. """ logger.info("Calculating the system temperature.") - t_background = (background.rebin(self.__lambda_line).qty.to( - u.W / (u.m ** 2 * u.Hz * u.sr), equivalencies=u.spectral_density(self.__lambda_line)) * - self.__lambda_line ** 2 / (2 * k_B) * u.sr).decompose() + # Add desired wavelength to wavelength bins + wl_bins = np.sort(np.append(self.__common_conf.wl_bins(), self.__lambda_line)).view(u.Quantity) + 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 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" if size == "point": - signal = signal.rebin(self.__lambda_line).qty.to(u.W / (u.m ** 2 * u.Hz), - equivalencies=u.spectral_density(self.__lambda_line)) - t_signal = (signal * self.__aperture_efficiency * self.__common_conf.d_aperture() ** 2 * - np.pi / 4 / (2 * k_B) * self.__eta_fss).decompose() + signal = SpectralQty(signal.wl, signal.qty.to(u.W / (u.m ** 2 * u.Hz), + equivalencies=u.spectral_density(signal.wl))) + t_signal = signal * (self.__aperture_efficiency * self.__common_conf.d_aperture() ** 2 * + np.pi / 4 / (2 * k_B) * self.__eta_fss) + t_signal = SpectralQty(t_signal.wl, t_signal.qty.decompose()) else: - signal = signal.rebin(self.__lambda_line).qty.to(u.W / (u.m ** 2 * u.Hz * u.sr), - equivalencies=u.spectral_density(self.__lambda_line)) - t_signal = (signal * u.sr * self.__main_beam_efficiency * self.__lambda_line ** 2 / ( - 2 * k_B) * self.__eta_fss).decompose() - logger.debug("Signal temperature: %1.2e K" % t_signal.value) + signal = SpectralQty(signal.wl, signal.qty.to(u.W / (u.m ** 2 * u.Hz * u.sr), + equivalencies=u.spectral_density(signal.wl))) + t_signal = signal * (u.sr * self.__main_beam_efficiency * self.__lambda_line ** 2 / ( + 2 * k_B) * self.__eta_fss) + 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("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 @staticmethod