import warnings
import numpy as np
from astropy import units as u
from scipy.signal import fftconvolve
from .io import (
_SENSORFILTER_FOCUS,
_SENSORFILTER_IMPLEMENTED,
get_sensor_config,
resolve_bandpass,
)
from .meta import _MetaHolder_
from .scene import Scene
from .sensor import Sensor
from .telescope import Telescope
from .utils import list_of_quantity_to_array
# Default aperture radius for the diffraction-limited Airy core. Defocused PSFs
# get an SNR-optimized radius instead; see Simulation.default_aperture.
DEFAULT_R_APER_MAS = 70
# import logging
# logging.basicConfig(level=logging.INFO)
# logger = logging.getLogger(__name__)
[docs]
def calculate_bg_normalization_magnitude(bg_surface_brightness, psf_area):
"""
Convert the Background Surface Brightness into the total magnitude given the PSF area.
The area needs to be in square arcseconds since the typical definition of
Surface Brightness is in units of magnitudes per arcseconds^2.
Parameters
----------
bg_surface_brightness : float
The background surface brightness in mag/arcsec^2.
psf_area : float
The area of the PSF in arcsec^2.
Returns
-------
float
The total background magnitude.
"""
bg_magnitude = bg_surface_brightness - 2.5 * np.log10(psf_area)
return bg_magnitude
def _psf_from_focus_level(focus_level):
"""
Create a PSF object from a focus level string.
Parameters
----------
focus_level : str
The focus level, one of "0wave", "1wave", or "2wave".
Returns
-------
AiryPSF or DefocusPSF
The PSF object corresponding to the focus level.
Raises
------
ValueError
If focus_level is not one of the expected values.
"""
from .psfsim import DEFOCUS_1WAVE_PATH, DEFOCUS_2WAVE_PATH, AiryPSF, DefocusPSF
if focus_level == "0wave":
return AiryPSF()
if focus_level == "1wave":
return DefocusPSF(DEFOCUS_1WAVE_PATH)
if focus_level == "2wave":
return DefocusPSF(DEFOCUS_2WAVE_PATH)
raise ValueError(
f"Unknown focus_level {focus_level!r}. Expected '0wave', '1wave', or '2wave'."
)
[docs]
class Simulation(_MetaHolder_):
"""
A class to manage and run image exposure simulations.
Attributes
----------
telescope : Telescope
The telescope used for the simulation.
sensor : Sensor
The sensor used for the simulation.
scene : Scene
The scene being observed.
psf_profile : dict
Calculated PSF profile parameters.
"""
_mutable_parameters = ["time", "r_aper_mas", "n_reads"]
def __init__(
self,
telescope,
sensor,
scene=None,
time=90,
r_aper_mas=None,
n_reads=1,
meta={},
):
"""
Initialize a Simulation object.
Parameters
----------
telescope : Telescope
Telescope to be used for this simulation.
sensor : Sensor
Sensor to be used for this simulation.
scene : Scene, optional
Scene to be used for this simulation. Default is None.
time : float or array_like, optional
Exposure time(s) in seconds. Default is 90.
r_aper_mas : float, optional
Radius of the aperture in milliarcseconds. None (default) picks it
per PSF: 70 mas for the Airy PSF, the SNR-optimized radius for a
defocused or custom PSF (see default_aperture).
n_reads : int, optional
Number of coadded frames (read noise is incurred per frame; saturation
is evaluated per-frame). Default is 1.
meta : dict, optional
Additional parameters to store in the meta dictionary. Default is {}.
"""
self._telescope = telescope
self._sensor = sensor
self.set_scene(scene)
self._default_psf = None # set by from_sensorfilter only
self._psf_profile = {}
self._image_render_bundle_cache = {}
self._effective_wavelength = None
input_parameters = {
key: value
for key, value in locals().items()
if key not in ["self", "telescope", "sensor", "scene", "meta"]
and value is not None
and not key.startswith("__")
}
non_attr_meta = {
key: value
for key, value in meta.items()
if key not in ["telescope", "sensor", "scene"]
}
super().__init__(non_attr_meta | input_parameters)
# ============= #
# properties #
# ============= #
[docs]
@classmethod
def from_config(cls, config):
"""
Initialize a Simulation from a configuration dictionary.
Parameters
----------
config : dict
Configuration dictionary containing 'telescope', 'sensor', and 'scene' keys.
Returns
-------
Simulation
The initialized Simulation object.
"""
# Telescope
config_telescope = config.get("telescope", None)
if config_telescope is not None:
telescope = Telescope.from_config(config_telescope)
else:
telescope = None
# Sensor
config_sensor = config.get("sensor", None)
if config_sensor is not None:
sensor = Sensor.from_config(config_sensor)
else:
sensor = None
# Scene
config_scene = config.get("scene", None)
if config_scene is not None:
scene = Scene.from_config(config_scene)
else:
scene = None
return cls(telescope=telescope, sensor=sensor, scene=scene)
[docs]
@classmethod
def from_sensor(cls, sensor, scene=None):
"""
Naming simplification of from_sensor_and_scene().
Parameters
----------
sensor : str, list, tuple, or Sensor
Sensor specification.
scene : Scene, optional
Scene object.
Returns
-------
Simulation
"""
return cls.from_sensor_and_scene(sensor, scene=scene)
[docs]
@classmethod
def from_sensor_and_scene(cls, sensor, scene, gain_mode=None):
"""
Create a Simulation from a sensor specification and a scene.
Parameters
----------
sensor : str, list, tuple, or Sensor
Sensor specification. If str, format 'kind:band'.
scene : Scene
The scene to observe.
gain_mode : str, optional
Gain mode for detectors that have several (qCMOS: 'high' or 'low').
Ignored when `sensor` is already a Sensor.
Returns
-------
Simulation
"""
if type(sensor) is str or type(sensor) in [list, tuple]:
if type(sensor) is str:
kind, band = sensor.split(":", 1)
else:
kind, band = sensor
config = get_sensor_config(kind, band, gain_mode=gain_mode)
this = cls.from_config(config) # this has no scene
elif isinstance(sensor, Sensor):
this = cls(sensor=sensor, telescope=sensor.telescope)
else:
raise ValueError(f"I cannot parse the input {sensor=}")
this.set_scene(scene)
return this
[docs]
@classmethod
def from_sensorfilter(cls, sensorfilter, scene, gain_mode=None):
"""
Create a Simulation from a sensorfilter label, auto-selecting the PSF.
The label must be a key in sensor_info (e.g. 'zwo:r', 'zwo:r+1',
'qcmos:bb'). The PSF matching the sensor's focus_level is stored as
`default_psf`; get_image_snr uses it when no psf= argument is given.
Parameters
----------
sensorfilter : str
Label from sensor_info, format 'kind:band' (e.g. 'zwo:r',
'zwo:r+1', 'qcmos:bb'). Only canonical labels are accepted;
nicknames like 'sony:r' are not resolved here (use
from_sensor_and_scene for those).
scene : Scene
gain_mode : str, optional
Gain mode for detectors that have several (qCMOS: 'high' or 'low').
Returns
-------
Simulation
"""
if sensorfilter not in _SENSORFILTER_FOCUS:
known = sorted(_SENSORFILTER_FOCUS)
raise ValueError(
f"Unknown sensorfilter {sensorfilter!r}. Known labels: {known}"
)
if not _SENSORFILTER_IMPLEMENTED[sensorfilter]:
raise NotImplementedError(
f"Sensorfilter {sensorfilter!r} is not yet implemented "
f"(no throughput curve available)."
)
focus_level = _SENSORFILTER_FOCUS[sensorfilter]
kind, band = sensorfilter.split(":", 1)
config = get_sensor_config(kind, band, gain_mode=gain_mode)
this = cls.from_config(config)
this.set_scene(scene)
this._default_psf = _psf_from_focus_level(focus_level)
return this
# ================ #
# methods #
# ================ #
[docs]
def set_scene(self, scene_or_config):
"""
Set the scene for the simulation.
Parameters
----------
scene_or_config : Scene or dict
The Scene object or its configuration.
"""
if isinstance(scene_or_config, dict):
scene = Scene.from_config(scene_or_config)
else:
scene = scene_or_config
self._scene = scene
# this should move inside scene eventually
self._h_spec_observation = None
self._h_bkgd_observation = None
# scene drives the count rates / PSF profile: drop derived caches
self._psf_profile = {}
self._image_render_bundle_cache = {}
self._effective_wavelength = None
[docs]
def set_sensor(self, sensor_or_config):
"""
Set the sensor for the simulation.
Parameters
----------
sensor_or_config : Sensor or dict
The Sensor object or its configuration.
"""
if isinstance(sensor_or_config, dict):
sensor = Sensor.from_config(sensor_or_config)
else:
sensor = sensor_or_config
self._sensor = sensor
self._psf_profile = {} # reset the psf profile
self._image_render_bundle_cache = {}
self._effective_wavelength = None
# these following entries might depend on sensor for the bandpass
self._h_spec_observation = None
self._h_bkgd_observation = None
[docs]
def set_telescope(self, telescope_or_config):
"""
Set the telescope for the simulation.
Parameters
----------
telescope_or_config : Telescope or dict
The Telescope object or its configuration.
"""
if isinstance(telescope_or_config, dict):
telescope = Telescope.from_config(telescope_or_config)
else:
telescope = telescope_or_config
self._telescope = telescope
self._psf_profile = {} # reset the psf profile
self._image_render_bundle_cache = {}
self._effective_wavelength = None
# ------- #
# update #
# ------- #
@staticmethod
def _fullkey_to_element_and_key(fullkey):
"""
Split a double-underscore key into element and attribute key.
Parameters
----------
fullkey : str
Key like 'sensor__gain'.
Returns
-------
tuple
(element, key)
"""
element, *keys = fullkey.split("__")
if len(keys) == 0:
key = element
element = None
else:
key = "__".join(
keys
) # trick to allow scene__background__mag => scene, background__mag
return element, key
def _fullkey_to_value(self, fullkey):
"""
Get value from a double-underscore key.
Parameters
----------
fullkey : str
Key like 'sensor__gain'.
Returns
-------
object
The value of the parameter.
"""
*origin, baseparam = fullkey.split("__")
origin = "__".join(origin)
if len(origin) == 0:
return self.meta.get(baseparam)
else:
return eval(f"self.{origin.replace('__', '.')}").meta.get(baseparam)
[docs]
def reset(self):
"""
Reset the simulation and its components.
"""
super().reset() # this resets the meta
for element in [self.telescope, self.sensor, self.scene]:
if element is not None:
element.reset()
self._psf_profile = {}
self._image_render_bundle_cache = {}
self._effective_wavelength = None
[docs]
def update(self, **kwargs):
"""
Update simulation parameters.
The keyword should in principle have the following structure
`scene__source__mag = 21` to change the `self.scene.source.mag`
parameter to 21.
To simplify the use, if the name has no ambiguity, you
can skip the first structure elements.
For instance, only `scene__` has a `source__mag`. So `source__mag=21`
will automatically be associated with `scene__source__mag=21`.
Same would work, for instance for `dark_current`: only `detector__`
has a `dark_current` mutable_parameter.
However `mag=21` is not clear enough, as `scene__source__mag`,
`scene__background__mag` or `scene__host__mag` exist.
Parameters
----------
**kwargs
Parameters to update. Can use double-underscore for sub-elements.
"""
update_this = {}
to_update = {"telescope": {}, "sensor": {}, "scene": {}}
for key, value in kwargs.items():
key = key.replace(".", "__") # generi trick, sensor.gain == sensor__gain.
# is that a fully defined name like sensor__dark_current ?
if np.any(
[
key.startswith(f"{element}__")
for element in ["telescope", "sensor", "scene"]
]
):
# yes ? easy then
element, down_key = self._fullkey_to_element_and_key(key)
to_update[element][down_key] = value
continue
# see if it is missing a key
fetch_key = [
fullkey for fullkey in self.mutable_parameters if fullkey.endswith(key)
]
if len(fetch_key) == 1:
# ok, well defined, easy:
element, down_key = self._fullkey_to_element_and_key(fetch_key[0])
if element is None:
# means not a mutable parameter of a sub-element:
update_this[down_key] = value
else:
# it is a sub-element
to_update[element][down_key] = value
continue
if len(fetch_key) > 1:
# not well defined, more than one entry exist for they key
raise KeyError(
f"several entries found matching {key=} : {fetch_key}. Please clarify."
)
# if we are here, it means fetch_key didn't match.
raise KeyError(f"no entries found matching {key=}")
self.update_telescope(**to_update["telescope"])
self.update_sensor(**to_update["sensor"])
self.update_scene(**to_update["scene"])
self._meta |= update_this
# any parameter change can affect the PSF/count rates: drop caches
self._psf_profile = {}
self._image_render_bundle_cache = {}
self._effective_wavelength = None
[docs]
def update_telescope(self, **kwargs):
"""
Update telescope parameters.
Parameters
----------
**kwargs
Telescope parameters.
"""
shortcuts = {"jitter": "jitter_sigma", "diameter": "diameter_primary"}
to_update = {shortcuts.get(key, key): value for key, value in kwargs.items()}
self.telescope.update(**to_update)
[docs]
def update_scene(self, **kwargs):
"""
Update scene parameters.
Parameters
----------
**kwargs
Scene parameters.
"""
self.scene.update(**kwargs)
[docs]
def update_sensor(self, **kwargs):
"""
Update sensor parameters.
Parameters
----------
**kwargs
Sensor parameters.
"""
self.sensor.update(**kwargs)
[docs]
def get_parameter(self, name, as_dict=False):
"""
Get simulation or sub-element parameters.
Parameters
----------
name : str or list of str
Parameter name(s).
as_dict : bool, optional
Whether to return as a dictionary. Default is False.
Returns
-------
list or dict
The requested parameters.
"""
values = []
names = np.atleast_1d(name)
for name_ in names:
name_ = name_.replace(".", "__") # accept this generic way
fetch_key = [
fullkey
for fullkey in self.mutable_parameters
if fullkey.endswith(name_)
]
if len(fetch_key) == 1:
values.append(self._fullkey_to_value(fetch_key[0]))
else:
if len(fetch_key) == 0:
warnings.warn(f"not matching found for {name_=}")
else:
warnings.warn(f"several matching found for {name_=} ; {fetch_key}")
values.append(None)
if as_dict:
return dict(zip(names, values))
return values
# ------- #
# GETTER #
# ------- #
[docs]
def get_countrates(self, scene=None, band=None, units="adu/s", as_dict=True):
"""DEPRECATED (Airy, in-aperture). Use _count_rate_components for the 2D
path or get_image_snr. Retained for the analytic Airy methods."""
import warnings
warnings.warn(
"get_countrates is deprecated (Airy in-aperture model); "
"use get_image_snr / _count_rate_components.",
DeprecationWarning,
stacklevel=2,
)
return self._countrates_in_aperture(
scene=scene, band=band, units=units, as_dict=as_dict
)
def _countrates_in_aperture(
self, scene=None, band=None, units="adu/s", as_dict=True
):
"""
Get the countrates for each element in the scene (Airy in-aperture).
Parameters
----------
scene : Scene, optional
The scene to use. Defaults to self.scene.
band : SpectralElement, optional
The bandpass to use. Defaults to sensor bandpass.
units : str, optional
Units of the countrate ('adu/s', 'e/s', 'e-/s'). Default is "adu/s".
as_dict : bool, optional
Whether to return as a dictionary. Default is True.
Returns
-------
dict or ndarray
The countrates.
"""
if units not in ["adu/s", "e/s", "e-/s"]:
raise ValueError(
f"unknown countrate units. Should be 'adu/s' or 'e/s'. {units=} given"
)
# INFO: self.psf_profile is computed automatically if needed.
if band is None:
band = self.sensor.bandpass
if scene is None:
scene = self.scene
if any(p is not None for p in scene.call_down("profile")):
warnings.warn(
"spatial profiles are ignored on the analytic Airy path; "
"use get_snr / get_image_snr for a Sersic host.",
stacklevel=2,
)
# these are the countrate in e/s
if np.any(scene.call_down("mag_is_surface_brightness")):
area = self.psf_profile["psf_area"].value # area in arcsec**2
else:
area = None
scene_observations = scene.get_observation(band=band, area=area, as_dict=True)
# loop over the scene elements and get the countrate for each.
countrates = {}
for element, observation in scene_observations.items():
# - scene
count_rate_total = (
observation.countrate(area=self.telescope.surface) * u.electron / u.ct
)
count_rate = count_rate_total * self.psf_profile["ee_at_aper"]
if units in ["adu/s"]: # in [] enables short cut.
count_rate /= self.sensor.gain # ADU/s
countrates[element] = count_rate
if as_dict:
return countrates
return list_of_quantity_to_array(countrates.values())
def _count_rate_components(self, band=None):
"""PSF-independent count-rate primitives for the 2D path (electron/s).
Elements are classified by their own ``surface_brightness`` flag, never
by their name (issue #63). synphot applies ``area`` only to
surface-brightness elements, so those come back per pixel; every other
element comes back as a TOTAL rate that must be distributed spatially.
Returns
-------
dict with keys
``source_rate_total`` : the source, total electron/s.
``background_rate_per_pix`` : the sky element, electron/s/pixel.
``diffuse_rate_per_pix`` : other surface-brightness elements (e.g. a
resolved host given in mag/arcsec^2), electron/s/pixel.
``contaminant_rate_total`` : non-surface-brightness elements other
than the source (e.g. an unresolved host), total electron/s.
Co-located with the source and rendered through the same PSF.
``extended`` : list of (rate, profile, is_surface_brightness) for
elements with a spatial profile (e.g. a Sersic host). Rendered
by extended_rate_image; rate is total electron/s, or per pixel
at mu_e when is_surface_brightness.
Unlike get_countrates, this does NOT apply the Airy ee_at_aper, so it is
independent of compute_psf_profile and correct for any PSF. The rendered
PSF (in _image_render_bundle) handles spatial distribution.
"""
if band is None:
band = self.sensor.bandpass
plate_scale_arcsec = (
self.sensor.get_plate_scale(self.telescope).to(u.arcsec / u.pix).value
)
pixel_area_arcsec2 = plate_scale_arcsec**2 # one detector pixel, arcsec^2
surf = self.telescope.surface
# One call at one-pixel area: non-surface-brightness elements (the point
# source) ignore `area` -> total rate; surface-brightness elements (sky)
# use it -> per-pixel rate.
obs = self.scene.get_observation(
band=band, area=pixel_area_arcsec2, as_dict=True
)
sb_flags = self.scene.call_down("mag_is_surface_brightness", as_dict=True)
profiles = self.scene.call_down("profile", as_dict=True)
extended = []
source_rate_total = 0.0
background_rate_per_pix = 0.0
diffuse_rate_per_pix = 0.0
contaminant_rate_total = 0.0
for name, o in obs.items():
if o is None:
continue
rate = (
(o.countrate(area=surf) * u.electron / u.ct).to(u.electron / u.s).value
)
if profiles.get(name) is not None:
# Spatially resolved element: rendered by extended_rate_image.
extended.append((rate, profiles[name], bool(sb_flags.get(name))))
continue
if sb_flags.get(name, False):
# mag/arcsec^2: synphot already evaluated this at one pixel.
if name == "background":
background_rate_per_pix = rate
else:
diffuse_rate_per_pix += rate
elif name == "source":
source_rate_total = rate
else:
# A total rate for an element with no spatial extent: treat it
# as unresolved and co-located with the source, so it is spread
# by the PSF. Replicating this total into every pixel is the
# #63 bug. An extended host must be given as a surface
# brightness (mag/arcsec^2) to be treated per-pixel.
contaminant_rate_total += rate
return {
"source_rate_total": source_rate_total,
"background_rate_per_pix": background_rate_per_pix,
"diffuse_rate_per_pix": diffuse_rate_per_pix,
"contaminant_rate_total": contaminant_rate_total,
"extended": extended,
}
[docs]
def get_signal_and_variance(self, time=None, units="e-", n_reads=None):
"""
Get the signal and total variance for a given exposure time.
Parameters
----------
time : float or Quantity, optional
Exposure time. Defaults to self.meta['time'].
units : str, optional
Units of the signal ('e-', 'adu'). Default is "e-".
n_reads : int, optional
Number of coadded frames; the read-noise variance is incurred
n_reads times. Defaults to self.meta['n_reads'] (or 1).
Returns
-------
tuple
(source_signal, total_variance)
"""
import warnings
warnings.warn(
"get_signal_and_variance is the analytic Airy approximation and is "
"deprecated; use the PSF-aware 2D path (get_image_snr / "
"get_image_exptime_for_snr / get_peak_pixel).",
DeprecationWarning,
stacklevel=2,
)
if time is None:
time = self._meta.get("time", None)
if time is None:
raise ValueError("no time given, none set to meta")
# make sure the time is in the current units.
elif not isinstance(time, u.Quantity):
time = time * u.second
n_reads = self._resolve_n_reads(n_reads)
# get the count rates in {adu,e-}/s
count_rates = self._countrates_in_aperture(
units="e/s", as_dict=True
) # this is an array
# logging.info(f"Scene count rate: {count_rates:.2f} e/s")
# logging.info(f"Background count rate: {sky_count_rate:.2f} e/s")
source_signal = count_rates["source"] * time # e-
# poisson noise is at the electron level, not adu.
# poisson noise comes from the full scene source.
all_countrates = list_of_quantity_to_array(count_rates.values())
scene_signal = np.nansum(all_countrates) * time # e-
dark_signal = self.sensor.dark_current * time # e-/pix
# u.electron/u.pixel as variance, so unit square
detector_variance = (
dark_signal * u.electron / u.pixel + n_reads * self.sensor.read_noise**2
) * self.psf_profile["num_psf_pixels"] # e-**2
# logging.info(f"signal: {scene_signal:.2f} e-")
# logging.info(f"sky signal: {sky_signal:.2f} e-") part of the scene signal?
# logging.info(f"dark signal: {dark_signal:.2f} e-")
# logging.info(f"detector variance: {detector_variance:.2f} e-^2.")
# total noise
# * u.electron as photon noise ; already there in etector_variance
total_variance = (
scene_signal * u.electron + detector_variance
) # variance in e-**2
# units
if units.lower() == "adu":
scene_signal /= self.sensor.gain
total_variance /= self.sensor.gain**2
elif units not in ["e", "e-", "electron"]:
raise ValueError(f"unknown units {units=}. adu or electron/e- expected.")
return source_signal, total_variance
[docs]
def get_peak_pixel(
self,
time=None,
units="adu",
n_reads=None,
*,
psf=None,
jitter_sigma_mas=None,
npix=128,
oversample=11,
):
"""Brightest-pixel value for the actual (possibly defocused) PSF.
Peak = max over the per-frame clean image (_per_frame_clean_image_e):
(source + unresolved contaminants) through the rendered PSF, plus
extended profiles, sky, diffuse and dark, so defocused configs and an
offset host are handled correctly. `psf`
defaults to `default_psf` (diffraction-limited unless the simulation
was built by from_sensorfilter).
Saturation is per-frame (per-frame time = time / n_reads). Linear in
time, so scalar or array `time` both work.
"""
if time is None:
time = self._meta.get("time", None)
if time is None:
raise ValueError("no time given, none set to meta")
if not isinstance(time, u.Quantity):
time = time * u.second
n_reads = self._resolve_n_reads(n_reads)
tf = (time / n_reads).to(u.second).value # per-frame seconds (scalar or array)
if psf is None:
psf = self.default_psf
b = self._image_render_bundle(psf, jitter_sigma_mas, npix, oversample)
# Brightest pixel of the shared per-frame budget at unit time -- one
# image, so an offset host and the source peak are not naively summed.
# Linear in time.
peak_rate = float(self._per_frame_clean_image_e(b, 1.0).max())
peak_e = (peak_rate * tf) * u.electron
if units in ["e", "e-", "electron"]:
return peak_e
if units.lower() == "adu":
return (peak_e / self.sensor.gain).to(u.ct) + self.sensor.bias_level
raise ValueError(f"unknown units {units=}. 'adu' or electron/'e-' expected.")
[docs]
def is_saturated(
self,
time=None,
n_reads=None,
*,
psf=None,
jitter_sigma_mas=None,
npix=128,
oversample=11,
):
"""Whether the brightest pixel of the clean per-frame image reaches
sensor.saturation_e (min of full well and ADC clip after bias), for the
actual (possibly defocused) PSF. See get_peak_pixel and
Sensor.saturation_limit for which limit binds. Same threshold as
psfsim.saturation_mask_from_image_e and ImageSimulator.simulate (#64).
"""
peak_e = self.get_peak_pixel(
time,
units="e-",
n_reads=n_reads,
psf=psf,
jitter_sigma_mas=jitter_sigma_mas,
npix=npix,
oversample=oversample,
)
return peak_e >= self.sensor.saturation_e
[docs]
def peak_pixel_fraction(
self, psf=None, jitter_sigma_mas=None, npix=128, oversample=11
):
"""Fraction of total source flux in the brightest detector pixel for the
actual (possibly defocused) PSF.
PSF-aware replacement for the retired psf_profile['peak_pixel_fraction']
(which was in-focus-Airy only). `psf` defaults to `default_psf`.
"""
if psf is None:
psf = self.default_psf
b = self._image_render_bundle(psf, jitter_sigma_mas, npix, oversample)
return float(b["psf_norm"].max())
[docs]
def extended_rate_image(self, psf_norm, ctx, comps=None):
"""PSF-convolved e-/s/pix image of the profiled (extended) elements.
Each profile is rendered at the integer centre of a grid padded by
``npix`` (rounded to even) and convolved with ``psf_norm`` in
``mode="same"``, whose kernel origin is that same integer centre
(``psfsim.grid_center``). The host therefore lands wherever the PSF put
the source -- ``ctx.center`` or the PSF class's natural centre, odd or
even grid -- and ``dx``/``dy`` are offsets from the source. The padding
lets light from just outside the detector grid convolve inward -- no
edge loss. Zeros when the scene has no extended element.
"""
from .extended import render_sersic
from .psfsim import center_crop_or_pad
if comps is None:
comps = self._count_rate_components()
npix = ctx.npix
image = np.zeros((npix, npix))
if not comps["extended"]:
return image
big = npix + 2 * (npix // 2)
for rate, profile, is_sb in comps["extended"]:
image_big = rate * render_sersic(
profile,
ctx.plate_scale_mas / 1000.0,
big,
ctx.oversample,
total=not is_sb,
)
image += center_crop_or_pad(
fftconvolve(image_big, psf_norm, mode="same"), npix
)
return image
def _image_render_bundle(self, psf, jitter_sigma_mas, npix, oversample):
"""
Cached, time-independent inputs for the PSF-aware SNR/exptime path.
Returns a dict with the rendered normalized PSF, the plate scale (mas),
and the source/diffuse count *rates* (electrons / s). Memoized on
render-affecting state; cleared by update()/set_sensor()/set_telescope().
"""
from .psfsim import ImageSimulator
if jitter_sigma_mas is None:
jitter_sigma_mas = self.telescope.jitter_sigma.to("mas").value
cache = getattr(self, "_image_render_bundle_cache", None)
if cache is None:
cache = self._image_render_bundle_cache = {}
key = (psf.cache_key(), jitter_sigma_mas, int(npix), int(oversample))
if key in cache:
return cache[key]
imsim = ImageSimulator(self, npix=npix, oversample=oversample)
ctx = imsim._context(jitter_sigma_mas=jitter_sigma_mas)
psf_norm = psf.render(ctx)
comps = self._count_rate_components()
source_rate_total = comps["source_rate_total"]
diffuse_rate_per_pix = comps["diffuse_rate_per_pix"]
background_rate_per_pix = comps["background_rate_per_pix"]
contaminant_rate_total = comps["contaminant_rate_total"]
bundle = {
"psf_norm": psf_norm,
"plate_scale_mas": ctx.plate_scale_mas,
"source_rate_total": source_rate_total,
"diffuse_rate_per_pix": diffuse_rate_per_pix,
"background_rate_per_pix": background_rate_per_pix,
"contaminant_rate_total": contaminant_rate_total,
"extended_rate_image": self.extended_rate_image(psf_norm, ctx, comps),
}
cache[key] = bundle
return bundle
def _per_frame_clean_image_e(self, b, tf, source_scale=1.0):
"""Per-frame clean electron image for the saturation test.
Budget is source + unresolved contaminants (both through the PSF) +
extended profiles + background + diffuse + dark — the same charge every other API sees
(ImageSimulator.simulate, get_peak_pixel, the get_image_snr saturation
check). Issue #63: these used to disagree about the host. `b` is an
_image_render_bundle dict; `tf` is the per-frame integration time
(total time / n_reads); `source_scale` scales the source flux.
"""
dark_rate_per_pix = self.sensor.dark_current.to(
u.electron / (u.s * u.pix)
).value
psf_rate = b["source_rate_total"] * source_scale + b["contaminant_rate_total"]
return (
(psf_rate * tf) * b["psf_norm"]
+ b["background_rate_per_pix"] * tf
+ b["diffuse_rate_per_pix"] * tf
+ b["extended_rate_image"] * tf
+ dark_rate_per_pix * tf
)
[docs]
def get_image_snr(
self,
time=None,
mags=None,
psf=None,
r_aper_mas=None,
ee_frac=None,
optimize=False,
jitter_sigma_mas=None,
n_reads=None,
npix=128,
oversample=11,
warn=True,
):
"""
PSF-aware aperture signal-to-noise ratio.
Unlike get_snr_airy (the analytic Airy-disk approximation), this renders
the given PSF (default AiryPSF) on the detector grid and computes the SNR
for a circular aperture. get_snr now delegates to this method; the
in-focus default-aperture case reproduces get_snr_airy to within ~1%.
Aperture precedence: optimize > r_aper_mas > ee_frac; if none is
given, default_aperture decides (70 mas for the Airy PSF, optimized
otherwise).
Parameters
----------
time : float, array_like, or Quantity, optional
Exposure time(s) in seconds (bare floats are interpreted as
seconds). Defaults to meta['time']. Must be finite and >= 0; zero
time gives snr = 0. A scalar yields a dict of Python float/int; an
array yields a dict of equal-length ndarrays.
mags : float, array_like, or None, optional
Source magnitude(s) to evaluate at, sweeping the source brightness
relative to its set magnitude (host/background/dark/read held fixed).
None (default) leaves the source unchanged. A scalar yields a dict of
floats; a 1-D array yields a dict of arrays over magnitude. `time` and
`mags` may not both be arrays. Requires a scene with a source.
psf : PSFSource, optional
PSF model. If None, uses `default_psf` (the focus-level PSF when
built by from_sensorfilter, otherwise a diffraction-limited AiryPSF).
r_aper_mas : float, optional
Fixed aperture radius (mas).
ee_frac : float, optional
Aperture enclosing this fraction of the PSF.
optimize : bool, optional
If True, use the radius that maximizes SNR.
jitter_sigma_mas : float, optional
Override the telescope jitter (mas).
npix, oversample : int, optional
Render grid size and oversampling. `npix` must be large enough to
contain the PSF; the default (128) contains the bundled defocus PSFs
for the current sensors, but very small pixels or stronger defocus
may need a larger value.
warn : bool, optional
If True (default), emit a warning when any pixel in the rendered
per-frame image saturates. The 'n_saturated'/'saturated' values are
returned regardless of this flag.
Returns
-------
dict
'snr', 'signal_e', 'noise_e', 'enclosed_fraction', 'r_aper_mas',
'n_pix', plus saturation info 'n_saturated' (count of pixels in the
rendered per-frame image at or above full well / ADC clip) and
'saturated' (bool). Values are Python float/int/bool for scalar
`time` and `mags`, or ndarrays (n_pix/n_saturated as int, saturated
as bool) when `time` or `mags` is an array.
"""
from .psfsim import (
aperture_snr_radial,
saturation_mask_from_image_e,
select_aperture,
)
if time is None:
time = self._meta.get("time", None)
if time is None:
raise ValueError("no time given, none set to meta")
if not isinstance(time, u.Quantity):
time = time * u.second
t_val = np.asarray(time.to(u.second).value, dtype=float)
if not np.all(np.isfinite(t_val)) or np.any(t_val < 0):
raise ValueError(f"time must be finite and >= 0, got {time}")
n_reads = self._resolve_n_reads(n_reads)
if psf is None:
psf = self.default_psf
# validate mags inputs before the (expensive) render bundle
m0 = None
if mags is not None:
if not self.scene.has_source():
raise ValueError("mags sweep requires a scene with a source")
m0 = self.scene.source.mag
if m0 is None:
raise ValueError("mags sweep requires a source with a set magnitude")
m0 = m0.value
if np.ndim(mags) > 0 and not time.isscalar:
raise ValueError(
"time and mags cannot both be arrays; sweep one axis at a time"
)
b = self._image_render_bundle(psf, jitter_sigma_mas, npix, oversample)
# default to the ETC aperture if no mode was requested
if not optimize and r_aper_mas is None and ee_frac is None:
r_aper_mas, optimize = self.default_aperture(psf)
read_noise = self.sensor.read_noise.to(u.electron / u.pix).value * np.sqrt(
n_reads
)
dark_rate_per_pix = self.sensor.dark_current.to(
u.electron / (u.s * u.pix)
).value
def _snr_at(t_sec, source_scale=1.0):
source_e_total = b["source_rate_total"] * source_scale * t_sec
# Sky and any *resolved* (surface-brightness) element contribute
# per-pixel shot noise. The sky term lives in background_rate_per_pix
# and MUST be included here (it is what makes faint sources
# sky-limited); omitting it overstates SNR. The PSF-convolved
# extended-host image makes this a 2D per-pixel array.
diffuse_per_pix = (
b["diffuse_rate_per_pix"]
+ b["background_rate_per_pix"]
+ b["extended_rate_image"]
) * t_sec
# An unresolved host rides the source PSF instead (#63).
contaminant_e_total = b["contaminant_rate_total"] * t_sec
dark_per_pix = dark_rate_per_pix * t_sec
prof = aperture_snr_radial(
b["psf_norm"],
b["plate_scale_mas"],
source_e_total,
diffuse_per_pix,
dark_per_pix,
read_noise,
contaminant_e_total=contaminant_e_total,
)
idx = select_aperture(
prof, r_aper_mas=r_aper_mas, ee_frac=ee_frac, optimize=optimize
)
# per-frame clean electron image -> saturated-pixel count
# (matches get_peak_pixel/is_saturated: saturation is per-frame).
image_e = self._per_frame_clean_image_e(b, t_sec / n_reads, source_scale)
mask = saturation_mask_from_image_e(self.sensor, image_e)
return {
"snr": float(prof["snr"][idx]),
"signal_e": float(prof["signal_e"][idx]),
"noise_e": float(prof["noise_e"][idx]),
"enclosed_fraction": float(prof["enclosed_fraction"][idx]),
"r_aper_mas": float(prof["r_mas"][idx]),
"n_pix": int(prof["n_pix"][idx]),
"n_saturated": int(mask.sum()),
"saturated": bool(mask.any()),
"read_noise_e": read_noise,
}
def _assemble_array(results):
out = {
k: np.array([r[k] for r in results])
for k in (
"snr",
"signal_e",
"noise_e",
"enclosed_fraction",
"r_aper_mas",
"read_noise_e",
)
}
out["n_pix"] = np.array([r["n_pix"] for r in results], dtype=int)
out["n_saturated"] = np.array(
[r["n_saturated"] for r in results], dtype=int
)
out["saturated"] = np.array([r["saturated"] for r in results], dtype=bool)
return out
def _maybe_warn(result):
if warn and bool(np.any(result["saturated"])):
n = result["n_saturated"]
n_max = int(np.max(n)) if np.ndim(n) else int(n)
warnings.warn(
f"detector saturates — up to {n_max} pixel(s) "
f"at or above full well / ADC clip (per-frame, n_reads={n_reads}); "
f"SNR is unreliable. See is_saturated/get_peak_pixel.",
stacklevel=2,
)
return result
# resolve the source-flux scale from mags (validated above)
if mags is None:
source_scale = 1.0
elif np.ndim(mags) > 0:
scales = 10 ** (-0.4 * (np.asarray(mags, dtype=float) - m0))
t_sec = time.to(u.second).value
return _maybe_warn(_assemble_array([_snr_at(t_sec, s) for s in scales]))
else:
source_scale = 10 ** (-0.4 * (float(mags) - m0))
if time.isscalar:
return _maybe_warn(_snr_at(time.to(u.second).value, source_scale))
results = [_snr_at(t, source_scale) for t in time.to(u.second).value]
return _maybe_warn(_assemble_array(results))
[docs]
def get_image_exptime_for_snr(
self,
snr,
psf=None,
r_aper_mas=None,
ee_frac=None,
optimize=False,
jitter_sigma_mas=None,
n_reads=None,
npix=128,
oversample=11,
warn=True,
):
"""
Exposure time (s) to reach a target SNR on the PSF-aware path.
Inverse of get_image_snr. Aperture precedence: optimize > r_aper_mas >
ee_frac; if none is given default_aperture decides. Note
optimize here picks the radius that reaches the target SNR *fastest*
(minimum time), the inverse of get_image_snr's max-SNR optimize. Returns
a dict {'time_s', 'snr', 'r_aper_mas', 'enclosed_fraction', 'n_pix',
'n_saturated', 'saturated'}.
Parameters
----------
psf : PSFSource, optional
PSF model. If None, uses `default_psf` (the focus-level PSF when
built by from_sensorfilter, otherwise a diffraction-limited AiryPSF).
warn : bool, optional
If True (default), warn when the rendered image saturates at the
solved time; 'n_saturated'/'saturated' are returned regardless.
"""
from .psfsim import aperture_time_for_snr, saturation_mask_from_image_e
n_reads = self._resolve_n_reads(n_reads)
if psf is None:
psf = self.default_psf
b = self._image_render_bundle(psf, jitter_sigma_mas, npix, oversample)
dark_rate_per_pix = self.sensor.dark_current.to(
u.electron / (u.s * u.pix)
).value
read_noise = self.sensor.read_noise.to(u.electron / u.pix).value
if not optimize and r_aper_mas is None and ee_frac is None:
r_aper_mas, optimize = self.default_aperture(psf)
# Sky + resolved (surface-brightness) elements contribute per-pixel shot
# noise; the sky term must be included or the solved time is too short.
# An unresolved host rides the source PSF instead (#63).
diffuse_rate_per_pix = (
b["diffuse_rate_per_pix"]
+ b["background_rate_per_pix"]
+ b["extended_rate_image"]
)
result = aperture_time_for_snr(
b["psf_norm"],
b["plate_scale_mas"],
b["source_rate_total"],
diffuse_rate_per_pix,
dark_rate_per_pix,
read_noise,
n_reads=n_reads,
snr=snr,
r_aper_mas=r_aper_mas,
ee_frac=ee_frac,
optimize=optimize,
contaminant_rate_total=b["contaminant_rate_total"],
)
# saturated-pixel count at the solved time, on the per-frame clean image
# (one shared budget, matching get_image_snr / simulate / get_peak_pixel).
image_e = self._per_frame_clean_image_e(b, result["time_s"] / n_reads)
mask = saturation_mask_from_image_e(self.sensor, image_e)
result["n_saturated"] = int(mask.sum())
result["saturated"] = bool(mask.any())
if warn and result["saturated"]:
warnings.warn(
f"detector saturates at the solved time ({result['time_s']:.3g} s) "
f"— {result['n_saturated']} pixel(s) at or above full well / ADC clip "
f"(per-frame, n_reads={n_reads}).",
stacklevel=2,
)
return result
[docs]
def get_snr(
self,
time=None,
psf=None,
r_aper_mas=None,
ee_frac=None,
optimize=False,
jitter_sigma_mas=None,
n_reads=None,
npix=128,
oversample=11,
warn=True,
):
"""
Signal-to-noise ratio via the 2D image simulation (PSF-aware default).
Delegates to get_image_snr; see it for parameter details and the
returned dict. `time` may be a scalar or an array (returns a dict of
arrays). For the legacy analytic Airy approximation use get_snr_airy
(deprecated).
Returns
-------
dict
Same as get_image_snr: 'snr', 'signal_e', 'noise_e',
'enclosed_fraction', 'r_aper_mas', 'n_pix', 'n_saturated',
'saturated' (scalar values for scalar `time`, ndarrays for array
`time`). Pass warn=False to silence the saturation warning.
"""
return self.get_image_snr(
time=time,
psf=psf,
r_aper_mas=r_aper_mas,
ee_frac=ee_frac,
optimize=optimize,
jitter_sigma_mas=jitter_sigma_mas,
n_reads=n_reads,
npix=npix,
oversample=oversample,
warn=warn,
)
# ------- #
# PLOTS #
# ------- #
[docs]
def plot_snr(self, time, **kwargs):
"""SNR vs exposure time; see plotting.plot_snr."""
from . import plotting
return plotting.plot_snr(self, time, **kwargs)
plot = plot_snr
[docs]
def plot_snr_vs_mag(self, mags, **kwargs):
"""SNR vs source magnitude; see plotting.plot_snr_vs_mag."""
from . import plotting
return plotting.plot_snr_vs_mag(self, mags, **kwargs)
[docs]
def plot_sweep(self, key, values, **kwargs):
"""SNR vs any update() key; see plotting.plot_sweep."""
from . import plotting
return plotting.plot_sweep(self, key, values, **kwargs)
[docs]
def plot_peak_pixel(self, time, **kwargs):
"""Peak pixel vs time with the saturation ceiling; see plotting.plot_peak_pixel."""
from . import plotting
return plotting.plot_peak_pixel(self, time, **kwargs)
[docs]
def plot_exptime_for_snr(self, snr, **kwargs):
"""Exposure time vs target SNR; see plotting.plot_exptime_for_snr."""
from . import plotting
return plotting.plot_exptime_for_snr(self, snr, **kwargs)
[docs]
def plot_snr_vs_aperture(self, radii, **kwargs):
"""SNR vs aperture radius; see plotting.plot_snr_vs_aperture."""
from . import plotting
return plotting.plot_snr_vs_aperture(self, radii, **kwargs)
[docs]
def get_snr_airy(self, time=None, n_reads=None, warn=True):
"""
DEPRECATED analytic Airy-disk SNR approximation.
Use get_snr, which computes the SNR via the 2D image simulation.
Parameters
----------
time : float or Quantity, optional
Exposure time.
n_reads : int, optional
Number of reads. Defaults to meta['n_reads'] if set, else 1.
warn : bool, optional
If True (default), warn when the detector saturates at this time.
Returns
-------
Quantity
The SNR.
"""
warnings.warn(
"get_snr_airy (analytic Airy approximation) is deprecated; "
"use get_snr, which now uses the 2D image simulation.",
DeprecationWarning,
stacklevel=2,
)
signal, variance = self.get_signal_and_variance(time, n_reads=n_reads)
if warn:
try:
sat = self.is_saturated(time, n_reads=n_reads)
except ValueError:
sat = False # no time resolvable; skip the saturation check
if bool(np.any(sat)):
warnings.warn(
"detector saturates (per-frame); the analytic Airy SNR is "
"unreliable. Use get_snr for the saturated-pixel count.",
stacklevel=2,
)
return signal / np.sqrt(variance)
[docs]
def get_exptime_for_snr(self, snr, n_reads=None, warn=True):
"""
Exposure time (seconds, Quantity) to reach a target SNR on the analytic
(Airy) path. Inverse of get_snr_airy. Returns inf*u.s if the source rate is 0.
Parameters
----------
snr : float
Target SNR.
n_reads : int, optional
Number of reads. Defaults to meta['n_reads'] if set, else 1.
warn : bool, optional
If True (default), warn when the detector saturates at the solved
exposure time.
"""
import warnings
warnings.warn(
"get_exptime_for_snr is the analytic Airy approximation and is "
"deprecated; use the PSF-aware 2D path (get_image_snr / "
"get_image_exptime_for_snr / get_peak_pixel).",
DeprecationWarning,
stacklevel=2,
)
from .psfsim import solve_time_for_snr
A, B, C = self._snr_coefficients(n_reads=n_reads)
t = solve_time_for_snr(snr, A, B, C) * u.second
if (
warn
and np.isfinite(t.value)
and bool(np.any(self.is_saturated(t, n_reads=n_reads)))
):
warnings.warn(
f"detector saturates at the solved exposure time ({t:.3g}); "
"use get_image_exptime_for_snr for the saturated-pixel count.",
stacklevel=2,
)
return t
# -------------- #
# Internal #
# -------------- #
def _resolve_n_reads(self, n_reads):
"""n_reads from the argument, else meta['n_reads'], else 1. Integer >= 1."""
if n_reads is None:
n_reads = self._meta.get("n_reads", 1)
if not float(n_reads).is_integer() or n_reads < 1:
raise ValueError(f"n_reads must be an integer >= 1, got {n_reads!r}")
return int(n_reads)
def _snr_coefficients(self, n_reads=None):
"""
(A, B, C) floats for SNR(t) = A*t / sqrt(B*t + C), electrons & seconds.
A = source rate; B = total scene rate + n_pix*dark; C = n_pix*N*RN**2.
"""
n_reads = self._resolve_n_reads(n_reads)
count_rates = self._countrates_in_aperture(units="e/s", as_dict=True)
A = count_rates["source"].to(u.electron / u.s).value
all_rates = float(
np.nansum([r.to(u.electron / u.s).value for r in count_rates.values()])
)
n_pix = self.psf_profile["num_psf_pixels"]
n_pix = n_pix.value if isinstance(n_pix, u.Quantity) else float(n_pix)
dark = self.sensor.dark_current.to(u.electron / (u.s * u.pix)).value
rn = self.sensor.read_noise.to(u.electron / u.pix).value
B = all_rates + n_pix * dark
C = n_pix * n_reads * rn**2
return A, B, C
def _parse_bandpass(self, bandpass):
"""
Parse a bandpass name or object.
Parameters
----------
bandpass : str or SpectralElement
The bandpass. If 'sensor', returns the sensor's bandpass.
Returns
-------
SpectralElement
"""
if bandpass == "sensor":
return self.sensor.bandpass
return resolve_bandpass(bandpass)
def _compute_psf_profile_impl(self):
"""Internal implementation of the Airy PSF profile computation (no warning)."""
from .airy import get_airy_and_ee_curve
wavelength = self.effective_wavelength.to("m")
r_psf_mas, psf1d, ee, ee_at_aper = get_airy_and_ee_curve(
wavelength,
r_aper_mas=self._meta.get("r_aper_mas", DEFAULT_R_APER_MAS),
jitter_sigma_mas=self.telescope.jitter_sigma.to("mas"),
fnum=self.telescope.f_num,
D=self.telescope.diameter_primary.value,
pixel_size=self.sensor.pixel_size.value,
verbose=False,
)
# compute the number of pixels associated to the PSF
plate_scale = self.sensor.get_plate_scale(self.telescope) # in arcsec/pix
num_pixels_at_r = (
self._meta.get("r_aper_mas", DEFAULT_R_APER_MAS)
* u.arcsec
/ (plate_scale * 1000)
) # pix
num_psf_pixels = np.pi * num_pixels_at_r**2 # in pixels**2
# area of the psf in angular units
psf_area = num_psf_pixels * plate_scale**2 # in arcsec**2
# logging.info(f"PSF profile computed:")
# logging.info(f"EE={ee_at_aper:.2f} at {self.meta['r_aper_mas']} mas aperture")
# logging.info(f"num_psf_pixels={num_psf_pixels:.1f} pixels")
# logging.info(f"PSF area={psf_area:.2f} arcsec^2")
return {
"wavelength": wavelength,
"r_psf_mas": r_psf_mas,
"psf1d": psf1d,
"ee": ee,
"ee_at_aper": ee_at_aper,
"num_psf_pixels": num_psf_pixels,
"psf_area": psf_area,
}
[docs]
def compute_psf_profile(self):
"""
Compute the PSF profile and associated metrics.
DEPRECATED: use psf_profile (the cached property) or the PSF-aware 2D
path (get_image_snr / get_image_exptime_for_snr / get_peak_pixel).
Returns
-------
dict
Dictionary containing 'wavelength', 'r_psf_mas', 'psf1d', 'ee',
'ee_at_aper', 'num_psf_pixels', and 'psf_area'.
"""
import warnings
warnings.warn(
"compute_psf_profile is the analytic Airy approximation and is "
"deprecated; use the PSF-aware 2D path (get_image_snr / "
"get_image_exptime_for_snr / get_peak_pixel).",
DeprecationWarning,
stacklevel=2,
)
return self._compute_psf_profile_impl()
[docs]
def has_element(self, which):
"""
Check if the simulation has a specific component.
Parameters
----------
which : str
Component name ('telescope', 'sensor', 'scene').
Returns
-------
bool
"""
return getattr(self, which) is not None
# ================= #
# Properties #
# ================= #
@property
def scene(self):
"""
The Scene object.
"""
return self._scene
@property
def telescope(self):
"""
The Telescope object.
"""
return self._telescope
@property
def sensor(self):
"""
The Sensor object.
"""
return self._sensor
@property
def meta(self):
"""
Combined metadata from simulation and components.
"""
return self._meta | {
element_name: element.meta
for element_name in ["telescope", "sensor", "scene"]
if (element := getattr(self, element_name)) is not None
}
[docs]
def default_aperture(self, psf):
"""Aperture mode used when a call gives no r_aper_mas/ee_frac/optimize.
Returns ``(r_aper_mas, optimize)``: the Simulation's own ``r_aper_mas``
if one was set; else 70 mas for an ``AiryPSF``; else ``optimize=True``,
since the fixed Airy-core aperture holds only a few percent of a
defocused PSF (#95).
"""
from .psfsim import AiryPSF
r_aper_mas = self._meta.get("r_aper_mas")
if r_aper_mas is not None:
return r_aper_mas, False
if isinstance(psf, AiryPSF):
return DEFAULT_R_APER_MAS, False
return None, True
@property
def default_psf(self):
"""The PSF this simulation uses when no ``psf=`` argument is given.
``from_sensorfilter`` sets this from the filter's focus level (an
``AiryPSF`` for ``0wave``, a ``DefocusPSF`` for ``1wave`` / ``2wave``);
every other construction path leaves it diffraction limited.
"""
from .psfsim import AiryPSF
return self._default_psf if self._default_psf is not None else AiryPSF()
@property
def mutable_parameters(self):
"""
List of all mutable parameters (including sub-elements).
"""
return self._mutable_parameters + [
f"{element_name}__{k}"
for element_name in ["telescope", "sensor", "scene"]
if (element := getattr(self, element_name)) is not None
for k in element.mutable_parameters
]
@property
def effective_wavelength(self):
"""Photon-weighted effective wavelength of the bandpass times the source SED.
This is the wavelength the PSF is rendered at (see
:class:`wcc_etc.psfsim.AiryPSF`). Unlike ``sensor.wavelength`` -- the
pivot wavelength, a property of the throughput curve alone -- it moves
with the colour of the source, which matters because the diffraction
scale is linear in wavelength. Through the WCC broad band the pivot is
582 nm for every source, while the photons actually delivered average
541 nm for an O5V and 716 nm for an M5V; rendering the M5V PSF at the
pivot overpredicts its central-pixel fraction by ~26% and so predicts
premature saturation for red sources.
Falls back to ``sensor.wavelength`` when the simulation carries no
source spectrum. Cached, and invalidated by ``update`` / ``set_scene`` /
``set_sensor`` / ``set_telescope`` / ``reset`` like the render caches.
Returns
-------
Quantity
The effective wavelength, in nm.
"""
if getattr(self, "_effective_wavelength", None) is None:
self._effective_wavelength = self._compute_effective_wavelength()
return self._effective_wavelength
def _compute_effective_wavelength(self):
"""Photon-weight the bandpass by the source SED; fall back to the pivot."""
from .spectral import effective_wavelength
source = None if self.scene is None else self.scene.source
if source is None:
return self.sensor.wavelength
# apply_mag=False: the weighting is scale-free, so skip the normalization.
spectrum = source.get_spectrum(apply_mag=False, as_array=False)
if spectrum is None:
return self.sensor.wavelength
return effective_wavelength(self.sensor.bandpass, spectrum)
[docs]
def polychromatic_psf(self, n_sub=7, base=None):
"""A photon-weighted polychromatic PSF for this simulation's band and source.
The monochromatic default (rendered at :attr:`effective_wavelength`)
already removes the colour bias in the PSF *width*; this goes further and
coadds ``n_sub`` renders across the band, so the wings and the
diffraction ring structure are broadened the way a real broadband PSF is.
Parameters
----------
n_sub : int, optional
Number of equal-photon-weight sub-bands. Default 7. Cost is linear
in this; 5-9 is plenty for a broad optical band.
base : PSFSource, optional
PSF rendered at each sub-band wavelength. Defaults to the
simulation's own default PSF (``AiryPSF`` unless
``from_sensorfilter`` selected a defocused one).
Returns
-------
PolychromaticPSF
Pass it as ``psf=`` to ``get_image_snr`` / ``get_peak_pixel`` /
``peak_pixel_fraction`` / ``ImageSimulator.simulate``.
Examples
--------
>>> sim.get_peak_pixel(10, psf=sim.polychromatic_psf()) # doctest: +SKIP
"""
from .psfsim import PolychromaticPSF
if base is None:
base = self.default_psf
source = None if self.scene is None else self.scene.source
spectrum = (
source.get_spectrum(apply_mag=False, as_array=False)
if source is not None
else None
)
return PolychromaticPSF.from_bandpass(
self.sensor.bandpass, spectrum, n_sub=n_sub, base=base
)
# ---------- #
# cashed #
# ---------- #
@property
def psf_profile(self):
"""
The calculated PSF profile.
"""
if (
not hasattr(self, "_psf_profile")
or self._psf_profile is None
or len(self._psf_profile) == 0
): # like {}
self._psf_profile = self._compute_psf_profile_impl()
return self._psf_profile