"""
Data models
"""
# Copyright CFHT/CNRS/CEA/UParisSaclay
# Licensed under the MIT licence
from typing import Annotated, Optional, Tuple
from astropy import units as u #type: ignore[import-untyped]
import numpy as np
from pydantic import BaseModel, ConfigDict, Field, model_validator
from synphot import ( #type: ignore[import-untyped]
BlackBody1D,
ConstFlux1D,
Observation,
SourceSpectrum,
SpectralElement
)
from synphot.spectrum import BaseSpectrum #type: ignore[import-untyped]
from ... import package
from ..photsys import PhotSys
from ..types.quantity import AnnotatedQuantity
# Setup AB photometric system
abphotsys = PhotSys('abmag')
[docs]
def spectral_to_arrays(spectral: BaseSpectrum) -> Tuple[np.ndarray, np.ndarray]:
"""Convert a spectrum to wavelength and value arrays.
One zero-valued sample is retained on either side of the non-zero region
when available. An entirely zero-valued spectrum is returned unchanged.
Parameters
----------
spectral: synphot.spectrum.BaseSpectrum
Input spectrum with a defined wavelength set.
Returns
-------
wave: ~numpy.ndarray
Wavelength samples, normally an AstroPy quantity array.
values: ~numpy.ndarray
Spectrum values at ``wave``, normally an AstroPy quantity array.
"""
w = spectral.waveset
x = spectral(w)
# Trim extra 0 values at beginning and at the end
idx = np.where(x.value != 0.)[0]
if len(idx):
start = idx[0] - 1 if idx[0] > 0 else 0
end = idx[-1] + 2 if idx[-1] < w.size - 1 else w.size
w = w[start:end]
x = x[start:end]
return w, x
[docs]
class CacheModel(BaseModel):
'''
Pydantic model for cached data.
'''
tpeaks: dict[str, float]
zp_abmags: dict[str, AnnotatedQuantity( #type: ignore[valid-type]
unit = "s / ct",
decimals = 4,
description = "Instrumental AB magnitude zero-point"
)]
transmissions: dict[str, 'TransmissionModel']
emissions_ct: dict[str, float]
[docs]
class DetectorModel(BaseModel):
"""Physical model for an instrument detector, such as a CMOS or CCD.
Attributes
----------
gain: ~astropy.units.Quantity
Positive conversion factor in electrons per ADU.
ron: ~astropy.units.Quantity
Nonnegative RMS readout noise in electrons.
scale: ~astropy.units.Quantity
Two positive angular pixel scales in arcseconds per pixel.
full_well: ~astropy.units.Quantity
Positive full-well capacity in electrons; defaults to 1,000,000 electrons.
range: ~astropy.units.Quantity
Digital range in ADUs, greater than one ADU; defaults to 65,536 ADUs.
bias: ~astropy.units.Quantity
Nonnegative bias level in ADUs; defaults to zero ADUs.
Saturation calculations use the smaller
of ``range - 1 adu - bias`` and ``full_well / gain`` as signal headroom.
transmissions: dict[str, TransmissionModel]
Loaded quantum-efficiency curves.
emissions: dict[str, SBSEDModel]
Loaded detector emission spectra.
Examples
--------
>>> detector = DetectorModel(
... gain="2 electron/adu", ron="4 electron",
... scale="[0.3, 0.3] arcsec/pix", full_well="100000 electron",
... range="65536 adu", bias="1200 adu",
... transmissions={}, emissions={}
... )
>>> detector.range - 1 * u.adu - detector.bias
<Quantity 64335. adu>
"""
gain: AnnotatedQuantity( #type: ignore[valid-type]
unit = "electron / adu",
gt = 0. * u.electron / u.adu,
decimals = 4,
description = "Detector conversion factor."
)
ron: AnnotatedQuantity( #type: ignore[valid-type]
unit = "electron",
ge = 0. * u.electron,
decimals = 4,
description = "Read-out noise RMS amplitude."
)
scale: AnnotatedQuantity( #type: ignore[valid-type]
unit = "arcsec/pix",
gt = [0., 0.] * u.arcsec / u.pix,
min_shape = (2),
max_shape = (2),
decimals = 4,
description = "Angular pixel scale along each axis."
)
full_well: AnnotatedQuantity( #type: ignore[valid-type]
default = 1e6 * u.electron,
unit = "electron",
gt = 0. * u.electron,
decimals = 4,
description = "Detector full-well capacity in electrons."
)
range: AnnotatedQuantity( #type: ignore[valid-type]
default = 65536 * u.adu,
unit = "adu",
gt = 1. * u.adu,
decimals = 4,
description = "Digital range in ADUs; the highest level is range minus one ADU."
)
bias: AnnotatedQuantity( #type: ignore[valid-type]
default = 0. * u.adu,
unit = "adu",
ge = 0. * u.adu,
decimals = 4,
description = "Detector bias level in ADUs."
)
transmissions: dict[str, 'TransmissionModel']
emissions: dict[str, 'SBSEDModel']
[docs]
class InstrumentModel(BaseModel):
'''
Pydantic model for a PyDIET instrument.
'''
id: str
name: str
description: str
wavelength_range: AnnotatedQuantity( #type: ignore[valid-type]
unit = "nm",
gt = [0., 0.] * u.nm,
min_shape = (2),
max_shape = (2),
decimals = 4,
description = "Instrument minimum and maximum wavelengths."
)
obstruction_area: AnnotatedQuantity( #type: ignore[valid-type]
unit = "m2",
gt = 0. * u.m**2,
decimals = 3,
description = "Effective instrument obstruction area, resolved from the telescope when omitted in configuration."
)
overhead: AnnotatedQuantity( #type: ignore[valid-type]
unit = "s",
ge = 0. * u.second,
decimals = 3,
description = "Total instrument time overhead between exposures."
)
psf_beta: float = Field(
default = 3.2,
gt = 1.,
le = 10.,
description="Moffat beta parameter of the instrument PSF model."
)
filters: 'FiltersModel'
optics: 'OpticsModel'
detector: 'DetectorModel'
telescope: 'TelescopeModel'
site: 'SiteModel'
default: bool = False
cache: Optional['CacheModel'] = Field(default=None)
@model_validator(mode="after")
def _update_cache(self) -> 'InstrumentModel':
# Compute extra parameters during initialization
area = self.telescope.collecting_area - self.obstruction_area
# Transmission peaks (essentially to check that some flux goes through)
tpeaks: dict[str, float] = {}
# Instrumental magnitude zero-points in the AB system
zp_abmags: dict[str, float] = {}
# Filter emissions and transmissions
transmissions : dict[str, TransmissionModel] = {} #type: ignore[annotation-unchecked]
emissions_ct : dict[str, u.Quantity[u.ct/u.s]] = {} #type: ignore[annotation-unchecked]
for mirror_status in self.telescope.transmissions:
upstream_transmission: float | SpectralElement = 1.
mirror_transmission = self.telescope.transmissions[mirror_status]
mirror_emission = self.telescope.emissions[mirror_status]
# Pre-filter list of transmissions
transmission_list = [
mirror_transmission,
*self.optics.transmissions.values()
]
upstream_emission = SourceSpectrum(ConstFlux1D, amplitude=0.)
emission_list = [
mirror_emission,
*self.optics.emissions.values()
]
for i, v in enumerate(transmission_list):
emission = emission_list[i].spectral
transmission: SpectralElement = transmission_list[i].spectral
assert transmission is not None
upstream_transmission *= transmission
upstream_emission = upstream_emission * transmission + emission
for f in self.filters.transmissions:
filter = self.filters.transmissions[f]
assert filter.vars is not None
scale = filter.vars['scale'] if 'scale' in filter.vars else 1.
assert isinstance(scale, (int, float))
assert filter.spectral is not None
filter_transmission = scale * filter.spectral
filter_emission = self.filters.emissions[f].spectral
assert filter_transmission is not None
transmission = upstream_transmission * filter_transmission
emission = upstream_emission * transmission + filter_emission
detector_transmission = self.detector.transmissions["0"].spectral
assert detector_transmission is not None
transmission *= detector_transmission
emission *= detector_transmission
wave, response = spectral_to_arrays(transmission)
# Configuration ID includes mirror status ID and filter ID
config_id = f"{mirror_transmission.id}+{filter.id}" \
if mirror_transmission.id != "" \
else filter.id
# Compute transmission maximum to verify that light goes through
tpeaks[config_id] = transmission.tpeak()
# Compute instrumental magnitude zero-point as mag(s/ADU)
zp_abmags[config_id] = u.Magnitude(
self.detector.gain.value / \
Observation(
abphotsys.spectrum,
transmission
).countrate(area=area, binned=False)
) if tpeaks[config_id] > 0. else -100. * u.Magnitude(u.s / u.ct)
transmissions[config_id] = TransmissionModel(
id = config_id,
name = filter.name,
description = filter.description,
vars = filter.vars,
wave = wave.astype(np.float32),
response = response.astype(np.float32),
spectral = SpectralElement.from_spectrum1d(transmission.to_spectrum1d())
)
# Compute emission countrate
emissions_ct[config_id] = Observation(
emission,
transmission,
force='taper'
).countrate(area=area).astype(np.float32).value if tpeaks[config_id] > 0. else 0.
self.cache = CacheModel(
tpeaks=tpeaks,
zp_abmags=zp_abmags,
transmissions=transmissions,
emissions_ct=emissions_ct
)
del transmissions, emissions_ct
return self
model_config = ConfigDict(arbitrary_types_allowed=True)
[docs]
class OpticsModel(BaseModel):
'''
Pydantic model for optics.
'''
transmissions: dict[str, 'TransmissionModel']
emissions: dict[str, 'SBSEDModel']
[docs]
class FiltersModel(OpticsModel):
'''
Pydantic model for a filter set.
'''
pass
[docs]
class SBSEDModel(BaseModel):
'''
Pydantic model for a Surface Brightness Spectral Energy Distribution (SBSED).
'''
id: str
name: str
description: str
vars: dict[str, float|str] = {}
wave: AnnotatedQuantity( #type: ignore[valid-type]
unit = "nm",
ge = 100. * u.nm,
le = 100. * u.micron,
min_shape = (2),
max_shape = (100000),
decimals = 4
) | None = None
sbsed: AnnotatedQuantity( #type: ignore[valid-type]
unit = "Jy / arcsec2",
ge = 0. * u.Jy / u.arcsec**2,
min_shape = (2),
max_shape = (100000),
decimals = 6
) | None = None
spectral: SourceSpectrum = Field(exclude=True)
default: bool = False
model_config = ConfigDict(arbitrary_types_allowed=True)
[docs]
class SEDModel(BaseModel):
'''
Pydantic model for a Spectral Energy Distribution (SED).
'''
id: str
name: str
description: str
vars: dict[str, float] = {}
wave: AnnotatedQuantity( #type: ignore[valid-type]
unit = "nm",
ge = 100. * u.nm,
le = 100. * u.micron,
min_shape = (2),
max_shape = (100000),
decimals = 4
) | None = None
sed: AnnotatedQuantity( #type: ignore[valid-type]
unit = "Jy",
ge = 0. * u.Jy,
min_shape = (2),
max_shape = (100000),
decimals = 6
) | None = None
spectral: SourceSpectrum = Field(exclude=True)
default: bool = False
model_config = ConfigDict(arbitrary_types_allowed=True)
[docs]
class SiteModel(BaseModel):
'''
Pydantic model for an observing site.
'''
id: str
name: str
description: str
sky_transmissions: dict[str, 'TransmissionModel']
sky_emissions: dict[str, 'SBSEDModel']
default: bool = False
[docs]
class TelescopeModel(BaseModel):
'''
Pydantic model for a telescope.
'''
id: str
name: str
description: str
collecting_area: AnnotatedQuantity( #type: ignore[valid-type]
unit = "m**2",
gt = 0. * u.m**2,
decimals = 4,
description = "Full collecting area, ignoring obstructions."
)
obstruction_area: AnnotatedQuantity( #type: ignore[valid-type]
unit = "m**2",
gt = 0. * u.m**2,
decimals = 4,
description = "Default obstruction area for instruments without an explicit override."
)
transmissions: dict[str, 'TransmissionModel']
emissions: dict[str, 'SBSEDModel']
default: bool = False
[docs]
class TransmissionModel(BaseModel):
'''
Pydantic model for a transmission curve (with wavelength).
'''
id: str
name: str
description: str = ""
vars: Optional[dict[str, float | str]] = None
wave_range: AnnotatedQuantity( #type: ignore[valid-type]
unit = "nm",
gt = [0., 0.] * u.nm,
min_shape = (2),
max_shape = (2),
decimals = 4,
description = "Instrument minimum and maximum wavelengths."
) | None = None
wave: AnnotatedQuantity( #type: ignore[valid-type]
unit = "nm",
ge = 100. * u.nm,
le = 100. * u.micron,
min_shape = (2),
max_shape = (100000),
decimals = 3
) | None = None
response: AnnotatedQuantity( #type: ignore[valid-type]
unit = "",
ge = -100.,
le = 100.,
min_shape = (2),
max_shape = (100000),
decimals = 4
) | None = None
spectral: Optional[SpectralElement] = Field(default=None, exclude=True)
default: bool = False
model_config = ConfigDict(arbitrary_types_allowed=True)