2020-05-14 17:33:02 +02:00
|
|
|
from astropy import units as u
|
|
|
|
from .ASensor import ASensor
|
|
|
|
from ..IRadiant import IRadiant
|
|
|
|
from ..Entry import Entry
|
|
|
|
import numpy as np
|
2020-05-15 14:34:55 +02:00
|
|
|
from typing import Union, Tuple
|
2020-05-14 17:33:02 +02:00
|
|
|
from ..psf.Airy import Airy
|
|
|
|
from ..psf.Zemax import Zemax
|
|
|
|
from ..SpectralQty import SpectralQty
|
|
|
|
from .PixelMask import PixelMask
|
2020-05-29 09:37:05 +02:00
|
|
|
from ...lib.logger import logger
|
2020-05-14 17:33:02 +02:00
|
|
|
import astropy.constants as const
|
2020-05-23 16:25:18 +02:00
|
|
|
import os
|
|
|
|
import astropy.io.fits as fits
|
2020-05-14 17:33:02 +02:00
|
|
|
|
|
|
|
|
|
|
|
class Imager(ASensor):
|
|
|
|
"""
|
|
|
|
A class for modelling a Image-sensor
|
|
|
|
"""
|
|
|
|
__encircled_energy: Union[str, float, u.Quantity]
|
|
|
|
|
2020-09-10 15:02:52 +02:00
|
|
|
@u.quantity_input(pixel_geometry=u.pixel, pixel_size="length", sigma_read_out=u.electron ** 0.5 / u.pix,
|
2020-05-19 10:32:03 +02:00
|
|
|
center_offset=u.pix, dark_current=u.electron / u.pix / u.second, well_capacity=u.electron)
|
2020-05-14 17:33:02 +02:00
|
|
|
def __init__(self, parent: IRadiant, quantum_efficiency: Union[str, u.Quantity],
|
2020-09-10 15:02:52 +02:00
|
|
|
pixel_geometry: u.Quantity, pixel_size: u.Quantity, sigma_read_out: u.Quantity, dark_current: u.Quantity,
|
2020-05-15 11:16:42 +02:00
|
|
|
well_capacity: u.Quantity, f_number: Union[int, float], common_conf: Entry,
|
2020-05-19 10:32:03 +02:00
|
|
|
center_offset: u.Quantity = np.array([0, 0]) << u.pix, shape: str = "circle",
|
2020-09-09 13:14:01 +02:00
|
|
|
contained_energy: Union[str, int, float] = "FWHM", aperture_size: u.Quantity = None):
|
2020-05-14 17:33:02 +02:00
|
|
|
"""
|
|
|
|
Initialize a new Image-sensor model.
|
|
|
|
|
|
|
|
Parameters
|
|
|
|
----------
|
|
|
|
parent : IRadiant
|
|
|
|
The parent element of the optical component from which the electromagnetic radiation is received.
|
|
|
|
quantum_efficiency : Union[str, u.Quantity]
|
|
|
|
The quantum efficiency of the detector. This can be either the path to the file containing the values of
|
2020-05-19 10:32:03 +02:00
|
|
|
the spectral quantum efficiency or the overall quantum efficiency as astropy quantity.
|
2020-05-14 17:33:02 +02:00
|
|
|
pixel_geometry : u.Quantity
|
|
|
|
The geometry of the pixel array as Quantity in pixels with two entries:
|
|
|
|
[number of pixels in x-direction, number of pixels in y-direction]
|
|
|
|
pixel_size : length-Quantity
|
|
|
|
The edge length of a pixel (assumed to be square).
|
2020-09-10 15:02:52 +02:00
|
|
|
sigma_read_out : Quantity
|
2020-05-14 17:33:02 +02:00
|
|
|
The RMS-read noise per detector pixel in electrons^0.5 / pixel.
|
|
|
|
dark_current : Quantity
|
|
|
|
The dark current of a detector pixel in electrons / (pixels * s).
|
2020-05-15 11:16:42 +02:00
|
|
|
well_capacity : Quantity
|
|
|
|
The pixel's well capacity in electrons.
|
2020-05-14 17:33:02 +02:00
|
|
|
f_number : Union[int, float]
|
|
|
|
The f-number of the optical system.
|
|
|
|
common_conf : Entry
|
|
|
|
The common-Entry of the configuration.
|
|
|
|
center_offset : u.Quantity
|
|
|
|
The offset of the PSF-center relative to the center of the detector array as length-quantity with two
|
|
|
|
entries: [offset in x-direction, offset in y-direction]
|
|
|
|
shape : str
|
|
|
|
The shape of the photometric aperture. Can be either square or circle
|
|
|
|
contained_energy : Union[str, int, float]
|
|
|
|
The energy contained within the photometric aperture.
|
2020-09-09 13:14:01 +02:00
|
|
|
aperture_size : u.Quantity
|
|
|
|
The radius respectively the edge length of the photometric aperture.
|
2020-05-14 17:33:02 +02:00
|
|
|
"""
|
|
|
|
super().__init__(parent)
|
2020-05-19 10:32:03 +02:00
|
|
|
if type(quantum_efficiency) == str:
|
2020-05-14 17:33:02 +02:00
|
|
|
self.__quantum_efficiency = SpectralQty.fromFile(quantum_efficiency, u.nm, u.electron / u.photon)
|
2020-05-19 15:59:07 +02:00
|
|
|
else:
|
2020-05-14 17:33:02 +02:00
|
|
|
self.__quantum_efficiency = quantum_efficiency
|
|
|
|
self.__pixel_geometry = pixel_geometry
|
|
|
|
self.__pixel_size = pixel_size
|
2020-09-10 15:02:52 +02:00
|
|
|
self.__sigma_read_out = sigma_read_out
|
2020-05-14 17:33:02 +02:00
|
|
|
self.__dark_current = dark_current
|
2020-05-15 11:16:42 +02:00
|
|
|
self.__well_capacity = well_capacity
|
2020-05-14 17:33:02 +02:00
|
|
|
self.__f_number = f_number
|
|
|
|
self.__center_offset = center_offset
|
|
|
|
self.__shape = shape
|
|
|
|
self.__contained_energy = contained_energy
|
2020-09-09 13:14:01 +02:00
|
|
|
self.__aperture_size = aperture_size
|
2020-05-14 17:33:02 +02:00
|
|
|
self.__common_conf = common_conf
|
|
|
|
# Calculate central wavelength
|
|
|
|
self.__central_wl = self.__common_conf.wl_min() + (
|
|
|
|
self.__common_conf.wl_max() - self.__common_conf.wl_min()) / 2
|
|
|
|
# Parse PSF
|
|
|
|
if hasattr(common_conf, "psf") and common_conf.psf().lower() == "airy":
|
|
|
|
# Use an airy disk as PSF
|
|
|
|
self.__psf = Airy(self.__f_number, self.__central_wl, common_conf.d_aperture(), common_conf.psf.osf,
|
|
|
|
pixel_size)
|
|
|
|
else:
|
|
|
|
# Read PSF from Zemax file
|
|
|
|
self.__psf = Zemax(common_conf.psf(), self.__f_number, self.__central_wl, common_conf.d_aperture(),
|
|
|
|
common_conf.psf.osf, pixel_size)
|
|
|
|
|
|
|
|
@u.quantity_input(exp_time="time")
|
2020-07-07 09:11:21 +02:00
|
|
|
def calcSNR(self, background: SpectralQty, signal: SpectralQty, obstruction: float,
|
|
|
|
exp_time: u.Quantity) -> u.dimensionless_unscaled:
|
2020-05-14 17:33:02 +02:00
|
|
|
"""
|
|
|
|
Calculate the signal to background ratio (SNR) for the given exposure time using the CCD-equation.
|
|
|
|
|
|
|
|
Parameters
|
|
|
|
----------
|
2020-07-07 09:11:21 +02:00
|
|
|
background : SpectralQty
|
|
|
|
The received background radiation
|
|
|
|
signal : SpectralQty
|
|
|
|
The received signal radiation
|
|
|
|
obstruction : float
|
|
|
|
The obstruction factor of the aperture as ratio A_ob / A_ap
|
2020-05-14 17:33:02 +02:00
|
|
|
exp_time : time-Quantity
|
|
|
|
The exposure time to calculate the SNR for.
|
|
|
|
|
|
|
|
Returns
|
|
|
|
-------
|
2020-05-15 14:34:55 +02:00
|
|
|
snr : Quantity
|
|
|
|
The calculated SNR as dimensionless quantity
|
2020-05-14 17:33:02 +02:00
|
|
|
"""
|
2020-05-15 14:34:55 +02:00
|
|
|
# Calculate the electron currents
|
2020-07-07 09:11:21 +02:00
|
|
|
signal_current, background_current, read_noise, dark_current = self.__exposePixels(background, signal,
|
|
|
|
obstruction)
|
2020-05-18 14:44:08 +02:00
|
|
|
# Calculate the SNR using the CCD-equation
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("Calculating the SNR...", extra={"spinning": True})
|
2020-05-18 14:44:08 +02:00
|
|
|
snr = signal_current.sum() * exp_time / np.sqrt(
|
|
|
|
(signal_current + background_current + dark_current).sum() * exp_time + (read_noise ** 2).sum())
|
2020-05-19 10:34:10 +02:00
|
|
|
# Print information
|
2020-05-29 09:53:10 +02:00
|
|
|
for exp_time_ in exp_time if exp_time.size > 1 else [exp_time]:
|
|
|
|
self.__printDetails(signal_current * exp_time_, background_current * exp_time_, read_noise,
|
|
|
|
dark_current * exp_time_, "t_exp=%.2f s: " % exp_time_.value)
|
|
|
|
self.__output(signal_current * exp_time_, background_current * exp_time_, read_noise,
|
|
|
|
dark_current * exp_time_, "texp_%.2f" % exp_time_.value)
|
2020-05-18 14:44:08 +02:00
|
|
|
# Return the value of the SNR, ignoring the physical units (electrons^0.5)
|
|
|
|
return snr.value * u.dimensionless_unscaled
|
2020-05-14 17:33:02 +02:00
|
|
|
|
2020-05-15 14:58:47 +02:00
|
|
|
@u.quantity_input(snr=u.dimensionless_unscaled)
|
2020-07-07 09:11:21 +02:00
|
|
|
def calcExpTime(self, background: SpectralQty, signal: SpectralQty, obstruction: float, snr: u.Quantity) -> u.s:
|
2020-05-14 17:33:02 +02:00
|
|
|
"""
|
|
|
|
Calculate the necessary exposure time in order to achieve the given SNR.
|
|
|
|
|
|
|
|
Parameters
|
|
|
|
----------
|
2020-07-07 09:11:21 +02:00
|
|
|
background : SpectralQty
|
|
|
|
The received background radiation
|
|
|
|
signal : SpectralQty
|
|
|
|
The received signal radiation
|
|
|
|
obstruction : float
|
|
|
|
The obstruction factor of the aperture as ratio A_ob / A_ap
|
2020-05-15 14:34:55 +02:00
|
|
|
snr : Quantity
|
|
|
|
The SNR for which the necessary exposure time shall be calculated as dimensionless quantity.
|
2020-05-14 17:33:02 +02:00
|
|
|
|
|
|
|
Returns
|
|
|
|
-------
|
|
|
|
exp_time : Quantity
|
|
|
|
The necessary exposure time in seconds.
|
|
|
|
"""
|
2020-05-15 14:34:55 +02:00
|
|
|
# Calculate the electron currents
|
2020-07-07 09:11:21 +02:00
|
|
|
signal_current, background_current, read_noise, dark_current = self.__exposePixels(background, signal,
|
|
|
|
obstruction)
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("Calculating the exposure time...", extra={"spinning": True})
|
2020-05-15 14:34:55 +02:00
|
|
|
# Calculate the electron currents for all pixels
|
|
|
|
signal_current_tot = signal_current.sum()
|
|
|
|
# Fix the physical units of the SNR
|
2020-05-18 14:44:08 +02:00
|
|
|
snr = snr * u.electron ** 0.5
|
2020-05-15 14:34:55 +02:00
|
|
|
|
|
|
|
# Calculate the ratio of the background- and dark-current to the signal current as auxiliary variable
|
2020-05-19 10:45:42 +02:00
|
|
|
current_ratio = (background_current.sum() + dark_current.sum()) / signal_current_tot
|
2020-05-15 14:34:55 +02:00
|
|
|
# Calculate the necessary exposure time as inverse of the CCD-equation
|
2020-05-15 14:58:47 +02:00
|
|
|
exp_time = snr ** 2 * (
|
2020-05-19 10:45:42 +02:00
|
|
|
1 + current_ratio + np.sqrt((1 + current_ratio) ** 2 + 4 * (read_noise ** 2).sum() / snr ** 2)) / (
|
2020-05-18 14:44:08 +02:00
|
|
|
2 * signal_current_tot)
|
2020-05-19 10:34:10 +02:00
|
|
|
# Print information
|
2020-05-29 09:53:10 +02:00
|
|
|
for snr_, exp_time_ in zip(snr, exp_time) if snr.size > 1 else zip([snr], [exp_time]):
|
|
|
|
self.__printDetails(signal_current * exp_time_, background_current * exp_time_, read_noise,
|
|
|
|
dark_current * exp_time_, "SNR=%.2f: " % snr_.value)
|
|
|
|
self.__output(signal_current * exp_time_, background_current * exp_time_, read_noise,
|
|
|
|
dark_current * exp_time_, "snr_%.2f" % snr_.value)
|
2020-05-15 14:34:55 +02:00
|
|
|
return exp_time
|
|
|
|
|
2020-07-14 12:01:08 +02:00
|
|
|
# @u.quantity_input(exp_time="time", snr=u.dimensionless_unscaled, target_brightness=[u.mag, u.mag / u.sr])
|
2020-07-07 09:11:21 +02:00
|
|
|
def calcSensitivity(self, background: SpectralQty, signal: SpectralQty, obstruction: float, exp_time: u.Quantity,
|
2020-07-14 12:01:08 +02:00
|
|
|
snr: u.Quantity, target_brightness: u.Quantity) -> [u.mag, u.mag / u.sr]:
|
2020-05-16 15:52:27 +02:00
|
|
|
"""
|
|
|
|
Calculate the sensitivity of the telescope detector combination.
|
|
|
|
|
|
|
|
Parameters
|
|
|
|
----------
|
2020-07-07 09:11:21 +02:00
|
|
|
background : SpectralQty
|
|
|
|
The received background radiation
|
|
|
|
signal : SpectralQty
|
|
|
|
The received signal radiation
|
|
|
|
obstruction : float
|
|
|
|
The obstruction factor of the aperture as ratio A_ob / A_ap
|
2020-05-16 15:52:27 +02:00
|
|
|
exp_time : Quantity
|
|
|
|
The exposure time in seconds.
|
|
|
|
snr : Quantity
|
|
|
|
The SNR for which the sensitivity time shall be calculated.
|
|
|
|
target_brightness : Quantity
|
2020-07-14 12:01:08 +02:00
|
|
|
The target brightness in mag or mag / sr.
|
2020-05-16 15:52:27 +02:00
|
|
|
|
|
|
|
Returns
|
|
|
|
-------
|
|
|
|
sensitivity: Quantity
|
|
|
|
The sensitivity as limiting apparent star magnitude in mag.
|
|
|
|
"""
|
|
|
|
# Calculate the electron currents
|
2020-07-07 09:11:21 +02:00
|
|
|
signal_current, background_current, read_noise, dark_current = self.__exposePixels(background, signal,
|
|
|
|
obstruction)
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("Calculating the sensitivity...", extra={"spinning": True})
|
2020-05-16 15:52:27 +02:00
|
|
|
# Fix the physical units of the SNR
|
|
|
|
snr = snr * u.electron ** 0.5
|
|
|
|
signal_current_lim = snr * (snr + np.sqrt(
|
2020-05-18 14:44:08 +02:00
|
|
|
snr ** 2 + 4 * (exp_time * (background_current.sum() + dark_current.sum()) +
|
2020-05-19 09:40:47 +02:00
|
|
|
(read_noise ** 2).sum()))) / (2 * exp_time)
|
2020-05-19 10:34:10 +02:00
|
|
|
# Print information
|
2020-05-29 09:53:10 +02:00
|
|
|
for snr_, exp_time_, signal_current_lim_ in zip(snr, exp_time, signal_current_lim) if snr.size > 1 else zip(
|
|
|
|
[snr], [exp_time], [signal_current_lim]):
|
|
|
|
self.__printDetails(signal_current_lim_ * exp_time_, background_current * exp_time_, read_noise,
|
|
|
|
dark_current * exp_time_, "SNR=%.2f t_exp=%.2f s: " % (snr_.value, exp_time_.value))
|
|
|
|
self.__output(signal_current * signal_current_lim_ / signal_current.sum() * exp_time_,
|
|
|
|
background_current * exp_time_, read_noise, dark_current * exp_time_,
|
|
|
|
"snr_%.2f_texp_%.2f" % (snr_.value, exp_time_.value))
|
2020-07-14 12:01:08 +02:00
|
|
|
return target_brightness - 2.5 * np.log10(signal_current_lim / signal_current.sum()) * target_brightness.unit
|
2020-05-16 15:52:27 +02:00
|
|
|
|
2020-05-15 14:34:55 +02:00
|
|
|
@u.quantity_input(signal=u.electron, background=u.electron, read_noise=u.electron ** 0.5, dark=u.electron)
|
2020-05-18 14:44:08 +02:00
|
|
|
def __printDetails(self, signal: u.Quantity, background: u.Quantity, read_noise: u.Quantity,
|
|
|
|
dark: u.Quantity, prefix: str = ""):
|
2020-05-15 14:34:55 +02:00
|
|
|
"""
|
2020-05-19 09:40:47 +02:00
|
|
|
Print details on the signal and noise composition.
|
2020-05-15 14:34:55 +02:00
|
|
|
|
|
|
|
Parameters
|
|
|
|
----------
|
|
|
|
signal : Quantity
|
|
|
|
The collected electrons from the target in electrons.
|
|
|
|
background : Quantity
|
|
|
|
The collected electrons from the background in electrons.
|
|
|
|
read_noise : Quantity
|
|
|
|
The read noise in electrons.
|
|
|
|
dark : Quantity
|
|
|
|
The electrons from the dark current in electrons.
|
2020-05-23 16:25:18 +02:00
|
|
|
prefix : str
|
|
|
|
The prefix to be used for printing
|
2020-05-15 14:34:55 +02:00
|
|
|
|
|
|
|
Returns
|
|
|
|
-------
|
|
|
|
"""
|
|
|
|
# Calculate the total collected electrons per pixel
|
2020-05-15 11:18:03 +02:00
|
|
|
total = signal + background + dark
|
2020-05-15 14:34:55 +02:00
|
|
|
# Check for overexposed pixels
|
2020-05-15 11:18:03 +02:00
|
|
|
overexposed = total > self.__well_capacity
|
|
|
|
if np.any(overexposed):
|
2020-10-13 11:47:27 +02:00
|
|
|
# Calculate number of overexposed pixels
|
|
|
|
n_overexposed = np.count_nonzero(overexposed)
|
2020-05-15 14:34:55 +02:00
|
|
|
# Show a warning for the overexposed pixels
|
2020-10-13 11:47:27 +02:00
|
|
|
if n_overexposed == 1:
|
|
|
|
logger.warning(prefix + str(n_overexposed) + " pixel is overexposed.")
|
|
|
|
else:
|
|
|
|
logger.warning(prefix + str(n_overexposed) + " pixels are overexposed.")
|
2020-10-16 09:45:48 +02:00
|
|
|
logger.info("--------------------------------------------------------------------------")
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info(prefix + "Collected electrons from target: %1.2e electrons" % signal.sum().value)
|
|
|
|
logger.info(prefix + "Collected electrons from background: %1.2e electrons" % background.sum().value)
|
|
|
|
logger.info(prefix + "Electrons from dark current: %1.2e electrons" % dark.sum().value)
|
|
|
|
logger.info(prefix + "Read noise: %1.2e electrons" % (read_noise ** 2).sum().value)
|
|
|
|
logger.info(prefix + "Total collected electrons: %1.2e electrons" % total.sum().value)
|
2020-10-16 09:45:48 +02:00
|
|
|
logger.info("--------------------------------------------------------------------------")
|
2020-05-20 09:13:59 +02:00
|
|
|
|
|
|
|
@u.quantity_input(signal=u.electron, background=u.electron, read_noise=u.electron ** 0.5, dark=u.electron)
|
|
|
|
def __output(self, signal: u.Quantity, background: u.Quantity, read_noise: u.Quantity,
|
2020-05-23 16:25:18 +02:00
|
|
|
dark: u.Quantity, name: str):
|
2020-05-20 09:13:59 +02:00
|
|
|
"""
|
|
|
|
Write the signal and the noise in electrons to files.
|
|
|
|
|
|
|
|
Parameters
|
|
|
|
----------
|
|
|
|
signal : Quantity
|
|
|
|
The collected electrons from the target in electrons.
|
|
|
|
background : Quantity
|
|
|
|
The collected electrons from the background in electrons.
|
|
|
|
read_noise : Quantity
|
|
|
|
The read noise in electrons.
|
|
|
|
dark : Quantity
|
|
|
|
The electrons from the dark current in electrons.
|
2020-05-23 16:25:18 +02:00
|
|
|
name : str
|
|
|
|
The name of the configuration.
|
2020-05-20 09:13:59 +02:00
|
|
|
|
|
|
|
Returns
|
|
|
|
-------
|
|
|
|
"""
|
2020-05-23 16:25:18 +02:00
|
|
|
# Concatenate the paths
|
|
|
|
path = os.path.join(self.__common_conf.output.path, name)
|
|
|
|
try:
|
2020-07-13 15:23:49 +02:00
|
|
|
os.makedirs(path, exist_ok=True)
|
2020-05-23 16:25:18 +02:00
|
|
|
except FileExistsError:
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.warning("Output directory '" + path + "' already exists.")
|
2020-05-23 16:25:18 +02:00
|
|
|
# Calculate the indices of nonzero values and create a bounding rectangle
|
|
|
|
y, x = np.nonzero(signal)
|
|
|
|
y_min = min(y)
|
|
|
|
y_max = max(y)
|
|
|
|
x_min = min(x)
|
|
|
|
x_max = max(x)
|
|
|
|
# Write arrays to file
|
|
|
|
if self.__common_conf.output.format.lower() == "csv":
|
|
|
|
mes = "Range reduced to nonzero values.\nThe origin is in the top left corner, starting with 0.\n" + \
|
|
|
|
"Column index range: %d - %d\nRow index range: %d - %d\n" % (y_min, y_max, x_min, x_max)
|
|
|
|
np.savetxt(os.path.join(path, "signal.csv"), signal[y_min:(y_max + 1), x_min:(x_max + 1)].value,
|
|
|
|
delimiter=",", header="Signal in electrons\n" + mes)
|
|
|
|
np.savetxt(os.path.join(path, "background.csv"), background[y_min:(y_max + 1), x_min:(x_max + 1)].value,
|
|
|
|
delimiter=",", header="Background in electrons\n" + mes)
|
|
|
|
np.savetxt(os.path.join(path, "read_noise.csv"), read_noise[y_min:(y_max + 1), x_min:(x_max + 1)].value,
|
|
|
|
delimiter=",", header="Read noise in electrons\n" + mes)
|
|
|
|
np.savetxt(os.path.join(path, "dark_noise.csv"), dark[y_min:(y_max + 1), x_min:(x_max + 1)].value,
|
|
|
|
delimiter=",", header="Dark noise in electrons\n" + mes)
|
|
|
|
elif self.__common_conf.output.format.lower() == "fits":
|
|
|
|
mes = "Range reduced to nonzero values. The origin is in the top left corner, starting with 0. " + \
|
|
|
|
"Column index range: %d - %d Row index range: %d - %d " % (y_min, y_max, x_min, x_max)
|
|
|
|
hdu = fits.PrimaryHDU(header=fits.Header(dict(COMMENT="Simulation data created by ESBO-ETC.",
|
|
|
|
TELESCOP="ESBO-ETC")))
|
|
|
|
signal_hdu = fits.ImageHDU(signal[y_min:(y_max + 1), x_min:(x_max + 1)].value, name="signal",
|
|
|
|
header=fits.Header(dict(COMMENT="Signal in electrons. " + mes,
|
|
|
|
TELESCOP="ESBO-ETC")))
|
|
|
|
background_hdu = fits.ImageHDU(background[y_min:(y_max + 1), x_min:(x_max + 1)].value, name="background",
|
|
|
|
header=fits.Header(dict(COMMENT="Background in electrons. " + mes,
|
|
|
|
TELESCOP="ESBO-ETC")))
|
|
|
|
read_noise_hdu = fits.ImageHDU(read_noise[y_min:(y_max + 1), x_min:(x_max + 1)].value, name="read noise",
|
|
|
|
header=fits.Header(dict(COMMENT="Read noise in electrons. " + mes,
|
|
|
|
TELESCOP="ESBO-ETC")))
|
|
|
|
dark_hdu = fits.ImageHDU(dark[y_min:(y_max + 1), x_min:(x_max + 1)].value, name="dark noise",
|
|
|
|
header=fits.Header(dict(COMMENT="Dark noise in electrons. " + mes,
|
|
|
|
TELESCOP="ESBO-ETC")))
|
|
|
|
hdul = fits.HDUList([hdu, signal_hdu, background_hdu, read_noise_hdu, dark_hdu])
|
|
|
|
hdul.writeto(os.path.join(path, "results.fits"), overwrite=True)
|
2020-05-14 17:33:02 +02:00
|
|
|
|
2020-07-07 09:11:21 +02:00
|
|
|
def __exposePixels(self, background: SpectralQty, signal: SpectralQty,
|
|
|
|
obstruction: float) -> Tuple[u.Quantity, u.Quantity, u.Quantity, u.Quantity]:
|
2020-05-15 14:34:55 +02:00
|
|
|
"""
|
|
|
|
Expose the pixels and calculate the signal and noise electron currents per pixel.
|
|
|
|
|
2020-07-07 09:11:21 +02:00
|
|
|
Parameters
|
|
|
|
----------
|
|
|
|
background : SpectralQty
|
|
|
|
The received background radiation
|
|
|
|
signal : SpectralQty
|
|
|
|
The received signal radiation
|
|
|
|
obstruction : float
|
|
|
|
The obstruction factor of the aperture as ratio A_ob / A_ap
|
|
|
|
|
2020-05-15 14:34:55 +02:00
|
|
|
Returns
|
|
|
|
-------
|
|
|
|
signal_current : Quantity
|
|
|
|
The electron current from the target as PixelMask in electrons / s
|
|
|
|
background_current : Quantity
|
|
|
|
The electron current from the background as PixelMask in electrons / s
|
|
|
|
read_noise : Quantity
|
|
|
|
The read noise per pixel in electrons
|
|
|
|
dark_current : Quantity
|
|
|
|
The electron current from the dark noise as PixelMask in electrons / s
|
|
|
|
"""
|
|
|
|
# Calculate the total incoming electron current
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("Calculating incoming electron current...", extra={"spinning": True})
|
2020-07-07 09:11:21 +02:00
|
|
|
signal_current, size, obstruction, background_current = self.__calcIncomingElectronCurrent(background, signal,
|
|
|
|
obstruction)
|
2020-05-29 09:37:05 +02:00
|
|
|
# info("Finished calculating incoming electron current", extra={"spinning": False})
|
2020-05-15 14:34:55 +02:00
|
|
|
# Initialize a new PixelMask
|
2020-05-14 17:33:02 +02:00
|
|
|
mask = PixelMask(self.__pixel_geometry, self.__pixel_size, self.__center_offset)
|
|
|
|
if size.lower() == "extended":
|
2020-05-15 14:34:55 +02:00
|
|
|
# Target is extended, a diameter of 0 pixels results in a mask with one pixel marked
|
2020-05-14 17:33:02 +02:00
|
|
|
d_photometric_ap = 0 * u.pix
|
2020-05-15 14:34:55 +02:00
|
|
|
# Mask the pixels to be exposed
|
2020-05-14 17:33:02 +02:00
|
|
|
mask.createPhotometricAperture("circle", d_photometric_ap / 2, np.array([0, 0]) << u.pix)
|
|
|
|
else:
|
2020-05-15 14:34:55 +02:00
|
|
|
# Target is a point source
|
2020-09-09 13:14:01 +02:00
|
|
|
if self.__aperture_size is not None:
|
2020-10-05 20:45:05 +02:00
|
|
|
# Us ethe aperture size as diameter
|
|
|
|
d_photometric_ap = self.__aperture_size
|
2020-05-14 17:33:02 +02:00
|
|
|
else:
|
2020-05-15 14:34:55 +02:00
|
|
|
# Calculate the diameter of the photometric aperture from the given contained energy
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("Calculating the diameter of the photometric aperture...",
|
2020-07-07 09:11:21 +02:00
|
|
|
extra={"spinning": True})
|
2020-05-14 17:33:02 +02:00
|
|
|
d_photometric_ap = self.__calcPhotometricAperture(obstruction)
|
2020-10-05 20:45:13 +02:00
|
|
|
# Mask the pixels to be exposed
|
|
|
|
mask.createPhotometricAperture(self.__shape, d_photometric_ap / 2)
|
2020-05-15 14:34:55 +02:00
|
|
|
# Calculate the background current PixelMask
|
|
|
|
background_current = mask * background_current * u.pix
|
|
|
|
# Calculate the read noise PixelMask
|
2020-09-10 15:02:52 +02:00
|
|
|
read_noise = mask * self.__sigma_read_out * u.pix
|
2020-05-15 14:34:55 +02:00
|
|
|
# Calculate the dark current PixelMask
|
2020-05-14 17:33:02 +02:00
|
|
|
dark_current = mask * self.__dark_current * u.pix
|
2020-09-09 13:14:01 +02:00
|
|
|
if self.__aperture_size is None and size.lower() != "extended":
|
2020-05-15 16:22:18 +02:00
|
|
|
if type(self.__contained_energy) == str:
|
|
|
|
if self.__contained_energy.lower() == "peak":
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("The radius of the photometric aperture is %.2f pixels. This equals the peak value" % (
|
2020-05-15 16:22:18 +02:00
|
|
|
d_photometric_ap.value / 2))
|
|
|
|
elif self.__contained_energy.lower() == "fwhm":
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("The radius of the photometric aperture is %.2f pixels. This equals the FWHM" % (
|
2020-05-15 16:22:18 +02:00
|
|
|
d_photometric_ap.value / 2))
|
|
|
|
elif self.__contained_energy.lower() == "min":
|
2020-07-07 09:11:21 +02:00
|
|
|
logger.info(
|
|
|
|
"The radius of the photometric aperture is %.2f pixels. This equals the first minimum" % (
|
|
|
|
d_photometric_ap.value / 2))
|
2020-05-15 16:22:18 +02:00
|
|
|
else:
|
2020-07-07 09:11:21 +02:00
|
|
|
logger.info(
|
|
|
|
"The radius of the photometric aperture is %.2f pixels. This equals %.0f%% encircled energy" % (
|
|
|
|
d_photometric_ap.value / 2, self.__contained_energy))
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("The photometric aperture contains " + str(np.count_nonzero(mask)) + " pixels.")
|
2020-05-14 17:33:02 +02:00
|
|
|
if size.lower() != "extended":
|
2020-05-15 14:34:55 +02:00
|
|
|
# Map the PSF onto the pixel mask in order to get the relative irradiance of each pixel
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("Mapping the PSF onto the pixel grid...", extra={"spinning": True})
|
2020-05-14 17:33:02 +02:00
|
|
|
mask = self.__psf.mapToPixelMask(mask,
|
|
|
|
getattr(getattr(self.__common_conf, "jitter_sigma", None), "val", None),
|
|
|
|
obstruction)
|
2020-05-15 14:34:55 +02:00
|
|
|
# Calculate the signal current PixelMask
|
|
|
|
signal_current = mask * signal_current
|
|
|
|
return signal_current, background_current, read_noise, dark_current
|
2020-05-14 17:33:02 +02:00
|
|
|
|
|
|
|
def __calcPhotometricAperture(self, obstruction: float) -> u.Quantity:
|
2020-05-15 14:34:55 +02:00
|
|
|
"""
|
|
|
|
Calculate the diameter of the photometric aperture
|
|
|
|
|
|
|
|
Parameters
|
|
|
|
----------
|
|
|
|
obstruction : float
|
|
|
|
The obstruction factor as A_ob / A_ap.
|
|
|
|
|
|
|
|
Returns
|
|
|
|
-------
|
|
|
|
d_photometric_ap : Quantity
|
|
|
|
The diameter of the photometric aperture in pixels.
|
|
|
|
"""
|
2020-05-14 17:33:02 +02:00
|
|
|
# Calculate the reduced observation angle
|
|
|
|
jitter_sigma = getattr(getattr(self.__common_conf, "jitter_sigma", None), "val", None)
|
|
|
|
reduced_observation_angle = self.__psf.calcReducedObservationAngle(self.__contained_energy, jitter_sigma,
|
|
|
|
obstruction)
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.debug("Reduced observation angle: %.2f" % reduced_observation_angle.value)
|
2020-05-14 17:33:02 +02:00
|
|
|
# Calculate angular width of PSF
|
|
|
|
observation_angle = (reduced_observation_angle * self.__central_wl / self.__common_conf.d_aperture() *
|
|
|
|
180.0 / np.pi * 3600).decompose() * u.arcsec
|
|
|
|
# Calculate FOV of a single pixel
|
|
|
|
pixel_fov = (self.__pixel_size / (self.__f_number * self.__common_conf.d_aperture()) * 180.0 /
|
|
|
|
np.pi * 3600).decompose() * u.arcsec
|
|
|
|
# Calculate the radius of the photometric aperture in pixels
|
|
|
|
d_photometric_ap = observation_angle / pixel_fov
|
|
|
|
return d_photometric_ap * u.pix
|
|
|
|
|
2020-07-07 09:11:21 +02:00
|
|
|
def __calcIncomingElectronCurrent(self, background: SpectralQty, signal: SpectralQty,
|
|
|
|
obstruction: float) -> Tuple[u.Quantity, str, float, u.Quantity]:
|
2020-05-14 17:33:02 +02:00
|
|
|
"""
|
|
|
|
Calculate the detected electron current of the signal and the background.
|
|
|
|
|
2020-07-07 09:11:21 +02:00
|
|
|
Parameters
|
|
|
|
----------
|
|
|
|
background : SpectralQty
|
|
|
|
The received background radiation
|
|
|
|
signal : SpectralQty
|
|
|
|
The received signal radiation
|
|
|
|
obstruction : float
|
|
|
|
The obstruction factor of the aperture as ratio A_ob / A_ap
|
|
|
|
|
2020-05-14 17:33:02 +02:00
|
|
|
Returns
|
|
|
|
-------
|
|
|
|
signal_current : Quantity
|
|
|
|
The electron current on the detector caused by the target in electrons / s.
|
|
|
|
size : str
|
|
|
|
The size of the target.
|
|
|
|
obstruction : float
|
2020-05-15 14:34:55 +02:00
|
|
|
The obstruction factor as A_ob / A_ap.
|
2020-05-14 17:33:02 +02:00
|
|
|
background_current : Quantity
|
|
|
|
The electron current on the detector caused by the background in electrons / (s * pix).
|
|
|
|
"""
|
|
|
|
# Calculate the photon current of the background
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("Calculating the background photon current.")
|
2020-07-07 09:11:21 +02:00
|
|
|
background_photon_current = background * np.pi * (
|
2020-05-14 17:33:02 +02:00
|
|
|
self.__pixel_size.to(u.m) ** 2 / u.pix) / (4 * self.__f_number ** 2 + 1) * (1 * u.sr)
|
|
|
|
# Calculate the incoming photon current of the target
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.info("Calculating the signal photon current.")
|
2020-06-26 19:15:43 +02:00
|
|
|
size = "extended" if signal.qty.unit.is_equivalent(u.W / (u.m ** 2 * u.nm * u.sr)) else "point"
|
|
|
|
if size == "point":
|
|
|
|
signal_photon_current = signal * np.pi * (self.__common_conf.d_aperture() / 2) ** 2
|
|
|
|
else:
|
|
|
|
signal_photon_current = signal * np.pi * self.__pixel_size.to(u.m) ** 2 / (
|
|
|
|
4 * self.__f_number ** 2 + 1) * (1 * u.sr)
|
2020-05-14 17:33:02 +02:00
|
|
|
# Calculate the electron current of the background and thereby handling the photon energy as lambda-function
|
|
|
|
background_current = (
|
|
|
|
background_photon_current / (lambda wl: (const.h * const.c / wl).to(u.W * u.s) / u.photon) *
|
2020-06-08 11:36:04 +02:00
|
|
|
self.__quantum_efficiency).integrate().decompose()
|
2020-05-14 17:33:02 +02:00
|
|
|
# Calculate the electron current of the signal and thereby handling the photon energy as lambda-function
|
|
|
|
signal_current = (signal_photon_current / (lambda wl: (const.h * const.c / wl).to(u.W * u.s) / u.photon) *
|
2020-06-08 11:36:04 +02:00
|
|
|
self.__quantum_efficiency).integrate().decompose()
|
2020-05-29 09:37:05 +02:00
|
|
|
logger.debug("Signal current: %1.2e e-/s" % signal_current.value)
|
|
|
|
logger.debug("Target size: " + size)
|
|
|
|
logger.debug("Obstruction: %.2f" % obstruction)
|
|
|
|
logger.debug("Background current: %1.2e e-/s" % background_current.value)
|
2020-05-14 17:33:02 +02:00
|
|
|
return signal_current, size, obstruction, background_current
|
|
|
|
|
|
|
|
@staticmethod
|
2020-05-15 11:15:18 +02:00
|
|
|
def check_config(sensor: Entry, conf: Entry) -> Union[None, str]:
|
2020-05-14 17:33:02 +02:00
|
|
|
"""
|
|
|
|
Check the configuration for this class
|
|
|
|
|
|
|
|
Parameters
|
|
|
|
----------
|
2020-05-15 11:15:18 +02:00
|
|
|
sensor : Entry
|
2020-05-14 17:33:02 +02:00
|
|
|
The configuration entry to be checked.
|
2020-05-15 11:15:18 +02:00
|
|
|
conf: Entry
|
|
|
|
The complete configuration.
|
2020-05-14 17:33:02 +02:00
|
|
|
|
|
|
|
Returns
|
|
|
|
-------
|
|
|
|
mes : Union[None, str]
|
|
|
|
The error message of the check. This will be None if the check was successful.
|
|
|
|
"""
|
2020-05-15 11:15:18 +02:00
|
|
|
if not hasattr(sensor, "f_number"):
|
2020-05-14 17:33:02 +02:00
|
|
|
return "Missing container 'f_number'."
|
2020-05-15 11:15:18 +02:00
|
|
|
mes = sensor.f_number.check_float("val")
|
2020-05-14 17:33:02 +02:00
|
|
|
if mes is not None:
|
|
|
|
return "f_number: " + mes
|
2020-05-15 11:15:18 +02:00
|
|
|
if not hasattr(sensor, "pixel_geometry"):
|
2020-05-14 17:33:02 +02:00
|
|
|
return "Missing container 'pixel_geometry'."
|
2020-05-15 11:15:18 +02:00
|
|
|
mes = sensor.pixel_geometry.check_quantity("val", u.pix)
|
2020-05-14 17:33:02 +02:00
|
|
|
if mes is not None:
|
|
|
|
return "pixel_geometry: " + mes
|
2020-05-15 11:15:18 +02:00
|
|
|
if hasattr(sensor, "center_offset") and isinstance(sensor.center_offset, Entry):
|
|
|
|
mes = sensor.center_offset.check_quantity("val", u.pix)
|
2020-05-14 17:33:02 +02:00
|
|
|
if mes is not None:
|
|
|
|
return "center_offset: " + mes
|
|
|
|
|
|
|
|
# Check pixel
|
2020-05-15 11:15:18 +02:00
|
|
|
if not hasattr(sensor, "pixel"):
|
2020-05-14 17:33:02 +02:00
|
|
|
return "Missing container 'pixel'."
|
2020-05-15 11:15:18 +02:00
|
|
|
if not hasattr(sensor.pixel, "quantum_efficiency"):
|
2020-05-14 17:33:02 +02:00
|
|
|
return "Missing container 'quantum_efficiency'."
|
2020-05-19 15:59:07 +02:00
|
|
|
mes = sensor.pixel.quantum_efficiency.check_quantity("val", u.electron / u.photon)
|
2020-05-14 17:33:02 +02:00
|
|
|
if mes is not None:
|
2020-05-15 11:15:18 +02:00
|
|
|
mes = sensor.pixel.quantum_efficiency.check_file("val")
|
2020-05-14 17:33:02 +02:00
|
|
|
if mes is not None:
|
|
|
|
return "pixel -> quantum_efficiency: " + mes
|
2020-05-15 11:16:58 +02:00
|
|
|
if not hasattr(sensor.pixel, "pixel_size"):
|
2020-05-14 17:33:02 +02:00
|
|
|
return "Missing container 'pixel_size'."
|
2020-05-15 11:15:18 +02:00
|
|
|
mes = sensor.pixel.pixel_size.check_quantity("val", u.m)
|
2020-05-14 17:33:02 +02:00
|
|
|
if mes is not None:
|
|
|
|
return "pixel -> pixel_size: " + mes
|
2020-05-15 11:15:18 +02:00
|
|
|
if not hasattr(sensor.pixel, "dark_current"):
|
2020-05-14 17:33:02 +02:00
|
|
|
return "Missing container 'dark_current'."
|
2020-05-15 11:15:18 +02:00
|
|
|
mes = sensor.pixel.dark_current.check_quantity("val", u.electron / (u.pix * u.s))
|
2020-05-14 17:33:02 +02:00
|
|
|
if mes is not None:
|
|
|
|
return "pixel -> dark_current: " + mes
|
2020-05-15 11:15:18 +02:00
|
|
|
if not hasattr(sensor.pixel, "sigma_read_out"):
|
2020-05-14 17:33:02 +02:00
|
|
|
return "Missing container 'sigma_read_out'."
|
2020-05-15 11:15:18 +02:00
|
|
|
mes = sensor.pixel.sigma_read_out.check_quantity("val", u.electron ** 0.5 / u.pix)
|
2020-05-14 17:33:02 +02:00
|
|
|
if mes is not None:
|
|
|
|
return "pixel -> sigma_read_out: " + mes
|
2020-05-15 11:15:18 +02:00
|
|
|
if not hasattr(sensor.pixel, "well_capacity"):
|
|
|
|
return "Missing container 'well_capacity'."
|
|
|
|
mes = sensor.pixel.well_capacity.check_quantity("val", u.electron)
|
|
|
|
if mes is not None:
|
|
|
|
return "pixel -> well_capacity: " + mes
|
2020-05-14 17:33:02 +02:00
|
|
|
|
|
|
|
# Check photometric aperture
|
2020-06-26 19:15:43 +02:00
|
|
|
if not hasattr(sensor, "photometric_aperture"):
|
|
|
|
setattr(sensor, "photometric_aperture", Entry(shape=Entry(val="circle"),
|
|
|
|
contained_energy=Entry(val="FWHM")))
|
2020-10-15 08:38:51 +02:00
|
|
|
if not hasattr(sensor.photometric_aperture, "shape"):
|
|
|
|
return "Missing container 'shape'."
|
|
|
|
mes = sensor.photometric_aperture.shape.check_selection("val", ["square", "circle"])
|
|
|
|
if mes is not None:
|
|
|
|
return "photometric_aperture -> shape: " + mes
|
2020-09-09 13:14:01 +02:00
|
|
|
if hasattr(sensor.photometric_aperture, "aperture_size"):
|
|
|
|
mes = sensor.photometric_aperture.aperture_size.check_quantity("val", u.pix)
|
2020-06-26 19:15:43 +02:00
|
|
|
if mes is not None:
|
2020-09-09 13:14:01 +02:00
|
|
|
return "photometric_aperture -> aperture_size: " + mes
|
2020-06-26 19:15:43 +02:00
|
|
|
else:
|
|
|
|
if not hasattr(sensor.photometric_aperture, "contained_energy"):
|
|
|
|
return "Missing container 'contained_energy'."
|
|
|
|
mes = sensor.photometric_aperture.contained_energy.check_float("val")
|
|
|
|
if mes is not None:
|
|
|
|
if conf.common.psf().lower() == "airy":
|
|
|
|
mes = sensor.photometric_aperture.contained_energy.check_selection("val",
|
|
|
|
["peak", "FWHM", "fwhm",
|
|
|
|
"min"])
|
|
|
|
if mes is not None:
|
|
|
|
return "photometric_aperture -> contained_energy: " + mes
|
|
|
|
else:
|
|
|
|
mes = sensor.photometric_aperture.contained_energy.check_selection("val", ["FWHM", "fwhm"])
|
|
|
|
if mes is not None:
|
|
|
|
return "photometric_aperture -> contained_energy: " + mes
|