diff --git a/esbo_etc/classes/sensor/Heterodyne.py b/esbo_etc/classes/sensor/Heterodyne.py new file mode 100644 index 0000000..01c28dd --- /dev/null +++ b/esbo_etc/classes/sensor/Heterodyne.py @@ -0,0 +1,281 @@ +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 + + +class Heterodyne(ASensor): + """ + A class for modelling the behaviour of a superheterodyne spectrometer. + """ + + def __init__(self, parent: IRadiant, aperture_efficiency: float, main_beam_efficiency: float, + receiver_temp: u.Quantity, eta_fss: float, lambda_line: u.Quantity, kappa: float, common_conf: Entry, + n_on: float = None): + """ + Initialize a new heterodyne detector + + Parameters + ---------- + parent : IRadiant + The parent element of the optical component from which the electromagnetic radiation is received. + aperture_efficiency : float + The aperture efficiency of the antenna. + main_beam_efficiency : float + The main beam efficiency of the telescope. + receiver_temp : u.Quantity in Kelvins + The intrinsic noise temperature of all receiver components. + eta_fss : float + The forward scattering efficiency of the antenna. + lambda_line : u.Quantity + The wavelength to be used for calculating the SNR. + kappa : float + The backend degradation factor. + common_conf : Entry + The common-Entry of the configuration. + n_on : float + The number of on source observations. + """ + self.aperture_efficiency = aperture_efficiency + self.main_beam_efficiency = main_beam_efficiency + self.receiver_temp = receiver_temp + self.eta_fss = eta_fss + self.lambda_line = lambda_line + self.kappa = kappa + self.common_conf = common_conf + self.n_on = n_on + super().__init__(parent) + + @u.quantity_input(exp_time="time") + def getSNR(self, exp_time: u.Quantity) -> u.dimensionless_unscaled: + """ + Calculate the signal to background ratio (SNR) for the given exposure time using the CCD-equation. + + Parameters + ---------- + exp_time : time-Quantity + The exposure time to calculate the SNR for. + + Returns + ------- + snr : Quantity + The calculated SNR as dimensionless quantity + """ + # Calculate the signal and background temperatures + t_signal, t_background = self.calcTemperatures() + 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 + + @u.quantity_input(snr=u.dimensionless_unscaled) + def getExpTime(self, snr: u.Quantity) -> u.s: + """ + Calculate the necessary exposure time in order to achieve the given SNR. + + Parameters + ---------- + snr : Quantity + The SNR for which the necessary exposure time shall be calculated as dimensionless quantity. + + Returns + ------- + exp_time : Quantity + The necessary exposure time in seconds. + """ + # Calculate the signal and background temperatures + t_signal, t_background = self.calcTemperatures() + 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 + + @u.quantity_input(exp_time="time", snr=u.dimensionless_unscaled, target_brightness=u.mag) + def getSensitivity(self, exp_time: u.Quantity, snr: u.Quantity, target_brightness: u.Quantity) -> u.mag: + """ + Calculate the sensitivity of the telescope detector combination. + + Parameters + ---------- + exp_time : Quantity + The exposure time in seconds. + snr : Quantity + The SNR for which the sensitivity time shall be calculated. + target_brightness : Quantity + The target brightness in magnitudes. + + Returns + ------- + sensitivity: Quantity + The sensitivity as limiting apparent star magnitude in mag. + """ + # Calculate the signal and background temperatures + t_signal, t_background = self.calcTemperatures() + 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 + + @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, + t_signal: u.Quantity, prefix: str = ""): + """ + Print details on the signal and noise composition. + + Parameters + ---------- + t_sys : Quantity + The system temperature. + delta_nu : Quantity + The noise bandwidth. + t_rms : Quantity + The RMS antenna temperature. + t_signal : Quantity + The antenna temperature. + prefix : str + The prefix to be used for printing. + + Returns + ------- + """ + logger.info("-------------------------------------------------------------------------------------------------") + logger.info(prefix + "System temperature: %1.2e K" % t_sys.value) + logger.info(prefix + "Noise bandwidth: %1.2e K" % delta_nu.value) + logger.info(prefix + "RMS antenna temperature: %1.2e K" % t_rms.value) + logger.info(prefix + "Antenna temperature: %1.2e K" % t_signal.value) + logger.info("-------------------------------------------------------------------------------------------------") + + def calcTemperatures(self): + """ + Calculate the noise temperatures of the signal and the background radiation. + + Returns + ------- + t_signal : u.Quantity + The signal temperature in Kelvins. + t_background : u.Quantity + The background temperature in Kelvins. + """ + logger.info("Calculating the system temperature.") + t_background = (self._parent.calcBackground().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() + # Calculate the incoming photon current of the target + logger.info("Calculating the signal temperature.") + signal, obstruction = self._parent.calcSignal() + 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() + 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) + logger.debug("Target size: " + size) + logger.debug("Obstruction: %.2f" % obstruction) + logger.debug("Background temperature: %1.2e K" % t_background.value) + return t_signal, t_background + + @staticmethod + def check_config(sensor: Entry, conf: Entry) -> Union[None, str]: + """ + Check the configuration for this class + + Parameters + ---------- + sensor : Entry + The configuration entry to be checked. + conf: Entry + The complete configuration. + + Returns + ------- + mes : Union[None, str] + The error message of the check. This will be None if the check was successful. + """ + if not hasattr(sensor, "aperture_efficiency"): + return "Missing container 'aperture_efficiency'." + mes = sensor.aperture_efficiency.check_float("val") + if mes is not None: + return "aperture_efficiency: " + mes + if not hasattr(sensor, "main_beam_efficiency"): + return "Missing container 'main_beam_efficiency'." + mes = sensor.main_beam_efficiency.check_float("val") + if mes is not None: + return "main_beam_efficiency: " + mes + if not hasattr(sensor, "receiver_temp"): + return "Missing container 'receiver_temp'." + mes = sensor.receiver_temp.check_quantity("val", u.K) + if mes is not None: + return "receiver_temp: " + mes + if not hasattr(sensor, "eta_fss"): + return "Missing container 'eta_fss'." + mes = sensor.eta_fss.check_float("val") + if mes is not None: + return "eta_fss: " + mes + if not hasattr(sensor, "lambda_line"): + return "Missing container 'lambda_line'." + mes = sensor.lambda_line.check_quantity("val", u.nm) + if mes is not None: + return "lambda_line: " + mes + if not hasattr(sensor, "kappa"): + return "Missing container 'kappa'." + mes = sensor.kappa.check_float("val") + if mes is not None: + return "kappa: " + mes + if hasattr(sensor, "n_on") and isinstance(sensor.n_on, Entry): + mes = sensor.n_on.check_float("val") + if mes is not None: + return "n_on: " + mes diff --git a/esbo_etc/classes/sensor/SensorFactory.py b/esbo_etc/classes/sensor/SensorFactory.py index fa66564..7084dee 100644 --- a/esbo_etc/classes/sensor/SensorFactory.py +++ b/esbo_etc/classes/sensor/SensorFactory.py @@ -2,6 +2,7 @@ from ..IRadiant import IRadiant from ..Entry import Entry from .ASensor import ASensor from .Imager import Imager +from .Heterodyne import Heterodyne from ...lib.logger import logger @@ -49,5 +50,14 @@ class SensorFactory: options.photometric_aperture.contained_pixels, Entry): args["contained_pixels"] = options.photometric_aperture.contained_pixels() return Imager(**args) + elif options.type == "Heterodyne": + args = dict(parent=self.__parent, aperture_efficiency=options.aperture_efficiency(), + main_beam_efficiency=options.main_beam_efficiency(), receiver_temp=options.receiver_temp(), + eta_fss=options.eta_fss(), lambda_line=options.lambda_line(), kappa=options.kappa(), + common_conf=self.__common_conf) + if hasattr(options, "n_on"): + # noinspection PyCallingNonCallable + args["n_on"] = options.n_on() + return Heterodyne(**args) else: logger.error("Wrong sensor type: " + options.type) diff --git a/esbo_etc/classes/sensor/__init__.py b/esbo_etc/classes/sensor/__init__.py index f050240..3772c5f 100644 --- a/esbo_etc/classes/sensor/__init__.py +++ b/esbo_etc/classes/sensor/__init__.py @@ -1,3 +1,4 @@ from .ASensor import * from .Imager import * +from .Heterodyne import * from .SensorFactory import * diff --git a/tests/sensor/test_Heterodyne.py b/tests/sensor/test_Heterodyne.py new file mode 100644 index 0000000..217dbda --- /dev/null +++ b/tests/sensor/test_Heterodyne.py @@ -0,0 +1,50 @@ +from unittest import TestCase +import astropy.units as u +import numpy as np +from esbo_etc.classes.Config import Configuration +from esbo_etc.classes.target.FileTarget import FileTarget +from esbo_etc.classes.target.BlackBodyTarget import BlackBodyTarget +from esbo_etc.classes.optical_component.Atmosphere import Atmosphere +from esbo_etc.classes.optical_component.CosmicBackground import CosmicBackground +from esbo_etc.classes.optical_component.Mirror import Mirror +from esbo_etc.classes.sensor.Heterodyne import Heterodyne + + +class TestHeterodyne(TestCase): + def setUp(self): + self.config = Configuration("data/esbo-etc_defaults_heterodyne.xml").conf + self.heterodyne_args = dict(aperture_efficiency=0.55, main_beam_efficiency=0.67, + receiver_temp=1050 * u.K, eta_fss=0.97, lambda_line=157.774 * u.um, kappa=1.0, + common_conf=self.config.common) + self.target = FileTarget("data/target/line.csv", self.config.common.wl_bins()) + self.atmosphere = Atmosphere(self.target, "data/atmosphere/transmittance_great.csv") + self.cosmic = CosmicBackground(self.atmosphere, temp=220 * u.K, emissivity=0.14) + self.mirror = Mirror(self.cosmic, reflectance="data/mirror/reflectance_great.csv", emissivity=0.08, + temp=230 * u.K) + self.heterodyne = Heterodyne(self.mirror, **self.heterodyne_args) + + def test_getSNR(self): + snr = self.heterodyne.getSNR(1900 * u.s) + self.assertAlmostEqual(snr.value, 10.059625717085497) + + def test_getExpTime(self): + exp_time = 1900 * u.s + snr = self.heterodyne.getSNR(exp_time) + exp_time_ = self.heterodyne.getExpTime(snr) + self.assertAlmostEqual(exp_time.value, exp_time_.value) + + def test_getSensitivity(self): + exp_time = 1900 * u.s + target = BlackBodyTarget(self.config.common.wl_bins(), mag=20 * u.mag) + atmosphere = Atmosphere(target, "data/atmosphere/transmittance_great.csv") + cosmic = CosmicBackground(atmosphere, temp=220 * u.K, emissivity=0.14) + mirror = Mirror(cosmic, reflectance="data/mirror/reflectance_great.csv", emissivity=0.08, temp=230 * u.K) + heterodyne = Heterodyne(mirror, **self.heterodyne_args) + snr = heterodyne.getSNR(exp_time) + target = BlackBodyTarget(self.config.common.wl_bins(), mag=10 * u.mag) + atmosphere = Atmosphere(target, "data/atmosphere/transmittance_great.csv") + cosmic = CosmicBackground(atmosphere, temp=220 * u.K, emissivity=0.14) + mirror = Mirror(cosmic, reflectance="data/mirror/reflectance_great.csv", emissivity=0.08, temp=230 * u.K) + heterodyne = Heterodyne(mirror, **self.heterodyne_args) + sensitivity = heterodyne.getSensitivity(exp_time, snr, 10 * u.mag) + self.assertAlmostEqual(sensitivity.value, 20)