nirwals.physics.exposure

Functions for signal-to-noise ratio calculations.

detection_rates(area, grating_angle, grating_constant, observation)

Return the count rates of an observation, binned for the wavelength resolution.

The function bins the observation's wavelength bins together, so that the new bins cover the wavelength resolution element as tightly as possible. For example, if each observation wavelength bin covers 5 A and the wavelength resolution element is 23 A, then the new bins will consist of 5 of the original bins each. The first five original bins become the first new bin, the sixth to tenth original bin become the second new bin, and so on.

The wavelengths corresponding to the new bins are taken to be the midpoints of the new bins.

The flux in each new bin is integrated and multiplied by the given area. These rates are returned together with the new wavelengths.

This function assumes that the bin set of the observation is equidistant, and that each wavelength bin corresponds to the wavelength range covered by a pixel.

The rates of the first and last bin should be considered to have arbitrary values, as they are affected by boundary effects.

Parameters:

Name	Type	Description	Default
area	Quantity	Effective mirror area.	required
grating_angle	deg	Grating angle.	required
grating_constant	Quantity	The grating constant, i.e. the spacing between grooves.	required
observation	Observation	Observation.	required

Returns:

Type	Description
tuple	A tuple of wavelengths and corresponding rates.

Source code in nirwals/physics/exposure.py

def detection_rates(
    area: Quantity,
    grating_angle: u.deg,
    grating_constant: Quantity,
    observation: Observation,
) -> tuple[Quantity, Quantity]:
    """
    Return the count rates of an observation, binned for the wavelength resolution.

    The function bins the observation's wavelength bins together, so that the new bins
    cover the wavelength resolution element as tightly as possible. For example, if each
    observation wavelength bin covers 5 A and the wavelength resolution element is
    23 A, then the new bins will consist of 5 of the original bins each. The first five
    original bins become the first new bin, the sixth to tenth original bin become the
    second new bin, and so on.

    The wavelengths corresponding to the new bins are taken to be the midpoints of the
    new bins.

    The flux in each new bin is integrated and multiplied by the given area. These rates
    are returned together with the new wavelengths.

    This function assumes that the bin set of the observation is equidistant, and that
    each wavelength bin corresponds to the wavelength range covered by a pixel.

    The rates of the first and last bin should be considered to have arbitrary values,
    as they are affected by boundary effects.

    Parameters
    ----------
    area: Quantity
        Effective mirror area.
    grating_angle: Angle
        Grating angle.
    grating_constant: u.micron
        The grating constant, i.e. the spacing between grooves.
    observation: Observation
        Observation.

    Returns
    -------
    tuple
        A tuple of wavelengths and corresponding rates.
    """

    wavelengths = observation.binset
    wavelength_values = wavelengths.to(u.AA).value
    flux_values = observation(wavelengths).to(units.PHOTLAM).value
    delta_lambda = wavelengths[1] - wavelengths[0]
    delta_lambda_value = delta_lambda.to(u.AA).value
    pixel_flux_values = _bin_integrals(wavelength_values, flux_values)

    # Find the binning so that each bin covers the wavelength resolution element as
    # tightly as possible.
    binning = 1
    wre = wavelength_resolution_element(
        grating_angle=grating_angle, grating_constant=grating_constant
    )
    while binning * delta_lambda < wre:
        binning += 1

    # Get the wavelengths for the bin wavelengths and the fluxes in the bins.
    bin_wavelength_values: list[float] = []
    bin_flux_values: list[float] = []
    pixel = 0
    while pixel < len(wavelengths):
        # Get the wavelength in the middle of the bin. The current pixel is the first in
        # the bin and the distance between the midpoints of the bin's outer pixels is
        # the binsize less one pixel.
        bin_wavelength_value = (
            wavelength_values[0] + (pixel + 0.5 * (binning - 1)) * delta_lambda_value
        )
        bin_wavelength_values.append(bin_wavelength_value)

        # Add up the fluxes of the pixels in the bin. If we are running out of pixels
        # for the last bin, we assume a flux of 0 for the "missing" pixels.
        bin_flux_value = 0
        for _ in range(binning):
            pixel_flux = (
                pixel_flux_values[pixel] if pixel < len(pixel_flux_values) else 0
            )
            bin_flux_value += pixel_flux
            pixel += 1
        bin_flux_values.append(bin_flux_value)

    # Add the units, and convert fluxes to rates.
    bin_wavelengths = np.array(bin_wavelength_values) * u.AA
    rates = np.array(bin_flux_values) * area * u.AA * units.PHOTLAM

    # Returns wavelengths and rates.
    return bin_wavelengths, rates

electrons(area, exposures, exposure_time, grating_angle, grating_constant, observation)

Return the number of electrons accumulated for an observation.

The electron counts are calculated for the same bins which are used by the detection_rates function. They should not be confused with the detector counts, for which you still would have to divide by the gain.

Parameters:

Name	Type	Description	Default
area	Quantity	Effective mirror area.	required
exposures	int	Number of exposures.	required
exposure_time	s	Exposure time per exposure.	required
grating_angle	deg	Grating angle.	required
grating_constant	Quantity	The grating constant, i.e. the spacing between grooves.	required
observation	Observation	Observation.	required

Returns:

Type	Description
tuple[Quantity, Quantity]	A tuple with an array of wavelengths and an array of the corresponding electron counts.

Source code in nirwals/physics/exposure.py

def electrons(
    area: Quantity,
    exposures: int,
    exposure_time: u.s,
    grating_angle: u.deg,
    grating_constant: Quantity,
    observation: Observation,
) -> tuple[Quantity, Quantity]:
    """
    Return the number of electrons accumulated for an observation.

    The electron counts are calculated for the same bins which are used by the
    detection_rates function. They should not be confused with the detector counts,
    for which you still would have to divide by the gain.

    Parameters
    ----------
    area: Quantity
        Effective mirror area.
    exposures: int
        Number of exposures.
    exposure_time: Quantity
        Exposure time per exposure.
    grating_angle: Angle
        Grating angle.
    grating_constant: u.micron
        The grating constant, i.e. the spacing between grooves.
    observation: Observation
        Observation.

    Returns
    -------
    tuple[Quantity, Quantity]
        A tuple with an array of wavelengths and an array of the corresponding electron
        counts.
    """
    wavelengths, rates = detection_rates(
        area=area,
        grating_angle=grating_angle,
        grating_constant=grating_constant,
        observation=observation,
    )
    return wavelengths, exposures * exposure_time * rates

exposure_time(configuration)

Calculate the exposure time as a function of the signal-to-noise ratio.

The given configuration must include a wavelength lambda and a signal-to-noise ratio snr, and the exposure time are calculated at lambda for 101 equidistant SNR values in the interval [0, 2 * snr]. The choice of 101 (rather than 100) values is deliberate as it means that the given snr is one of the values. It has the index 50.

The function returns a tuple of the signal-to-noise ratios and corresponding exposure times.

Parameters:

Name	Type	Description	Default
configuration	Configuration	The configuration.	required

Returns:

Type	Description
tuple	A tuple of 101 signal-to-noise ratios and corresponding exposure times.

Source code in nirwals/physics/exposure.py

def exposure_time(configuration: Configuration) -> tuple[Quantity, Quantity]:
    """
    Calculate the exposure time as a function of the signal-to-noise ratio.

    The given configuration must include a wavelength lambda and a signal-to-noise ratio
    snr, and the exposure time are calculated at lambda for 101 equidistant SNR values
    in the interval [0, 2 * snr]. The choice of 101 (rather than 100) values is
    deliberate as it means that the given snr is one of the values. It has the index 50.

    The function returns a tuple of the signal-to-noise ratios and corresponding
    exposure times.

    Parameters
    ----------
    configuration: Configuration
        The configuration.

    Returns
    -------
    tuple
        A tuple of 101 signal-to-noise ratios and corresponding exposure times.
    """
    # Get the source and sky observation.
    source = source_observation(configuration)
    sky = sky_observation(configuration)

    # Get the source and sky detection rates,
    grating = cast(Grating, configuration.telescope.grating)
    wavelengths_source, rates_source = detection_rates(
        area=configuration.telescope.effective_mirror_area,
        grating_angle=grating.grating_angle,
        grating_constant=grating.grating_constant,
        observation=source,
    )
    wavelengths_sky, rates_sky = detection_rates(
        area=configuration.telescope.effective_mirror_area,
        grating_angle=grating.grating_angle,
        grating_constant=grating.grating_constant,
        observation=sky,
    )

    # Collect the remaining relevant configuration parameters.
    exposure = cast(Exposure, configuration.exposure)
    e = exposure.exposures
    snr_ = cast(SNR, exposure.snr)
    requested_wavelength = snr_.wavelength
    detector = cast(Detector, configuration.detector)
    read_noise = detector.read_noise
    samplings = detector.samplings
    sampling_mode = detector.sampling_mode
    r = readout_noise(
        read_noise=read_noise, samplings=samplings, sampling_mode=sampling_mode
    )

    # Define the SNR values for which to calculate the exposure time.
    sigma = np.linspace(0, 2 * float(snr_.snr), 101)

    # Find the source and sky rate for the requested wavelength.
    rates_source_spectrum = SpectralElement(
        Empirical1D,
        points=wavelengths_source,
        lookup_table=rates_source.to(units.PHOTLAM * u.cm**2 * u.AA).value,
    )
    rate_source = rates_source_spectrum(requested_wavelength)
    rates_sky_spectrum = SpectralElement(
        Empirical1D,
        points=wavelengths_sky,
        lookup_table=rates_sky.to(units.PHOTLAM * u.cm**2 * u.AA).value,
    )
    rate_sky = rates_sky_spectrum(requested_wavelength)

    # The exposure time t is defined by a quadratic equation r^2 + p t + q = 0.
    p = -(sigma**2) * (rate_source + rate_sky) / (e * rate_source**2)
    q = -(sigma**2) * r / (e * rate_source**2)
    t = -(p / 2) + np.sqrt((p / 2) ** 2 - q)

    # Add units and return the result.
    return sigma * u.dimensionless_unscaled, t * u.s

pixel_wavelength_range(grating_angle, grating_constant)

Return the wavelength range covered by a single CCD pixel.

Parameters:

Name	Type	Description	Default
grating_angle	deg	The grating angle, i.e. the angle of the incoming rays to the grating normal.	required
grating_constant	micron	The grating constant, i.e. the groove spacing.	required

Returns:

Type	Description
Quantity	The wavelength range covered by a single pixel.

Source code in nirwals/physics/exposure.py

def pixel_wavelength_range(grating_angle: u.deg, grating_constant: u.micron) -> u.AA:
    """
    Return the wavelength range covered by a single CCD pixel.

    Parameters
    ----------
    grating_angle: Angle
        The grating angle, i.e. the angle of the incoming rays to the grating normal.
    grating_constant: Quantity
        The grating constant, i.e. the groove spacing.

    Returns
    -------
    Quantity
        The wavelength range covered by a single pixel.
    """
    return (
        grating_constant
        * CCD_PIXEL_SIZE
        * math.cos(grating_angle.to(u.rad).value)
        / CAMERA_FOCAL_LENGTH
    )

readout_noise(read_noise, samplings, sampling_mode)

Return the readout noise for a single exposure.

Parameters:

Name	Type	Description	Default
read_noise	float	Read noise.	required
samplings	int	Number of samplings (per exposure).	required
sampling_mode	SamplingMode	Sampling mode, such as "Fowler".	required

Returns:

Type	Description
float	The readout noise.

Source code in nirwals/physics/exposure.py

def readout_noise(
    read_noise: float, samplings: int, sampling_mode: SamplingMode
) -> float:
    """
    Return the readout noise for a single exposure.

    Parameters
    ----------
    read_noise: float
        Read noise.
    samplings: int
        Number of samplings (per exposure).
    sampling_mode: SamplingMode
        Sampling mode, such as "Fowler".

    Returns
    -------
    float
        The readout noise.
    """
    match sampling_mode:
        case "Fowler":
            return read_noise**2 / (samplings / 2)
        case "Up-the-Ramp":
            return read_noise**2 / (samplings / 12)

sky_observation(configuration)

Returns the sky background observation for a given configuration.

The observation, represented by a synphot Observation object, is the photon rate detected, which means it is the sky background flux with the following losses applied:

• Telescope throughput
• Fibre throughput
• Filter transmission
• Grating efficiency
• Detector quantum efficiency

No atmospheric extinction is applied, as it is assumed that it is included in the background spectrum already.

The bin set of the Observation object is chosen so that the bin just size is equal to the wavelength range covered by a single pixel.

The bin set covers the wavelength range from lambda_min - 100 A to at least lambda_max + 100 A, where lambda_min and lambda_max are the minimum and maximum requested wavelengths. The extra 100 A are added to avoid artifacts at the boundaries of the requested wavelength range.

Parameters:

Name	Type	Description	Default
configuration	Configuration	Simulator configuration.	required

Returns:

Type	Description
Observation	The observation for the sky background.

Source code in nirwals/physics/exposure.py

def sky_observation(configuration: Configuration) -> Observation:
    """
    Returns the sky background observation for a given configuration.

    The observation, represented by a synphot Observation object, is the photon rate
    detected, which means it is the sky background flux with the following losses
    applied:

    - Telescope throughput
    - Fibre throughput
    - Filter transmission
    - Grating efficiency
    - Detector quantum efficiency

    No atmospheric extinction is applied, as it is assumed that it is included in the
    background spectrum already.

    The bin set of the Observation object is chosen so that the bin just size is equal
    to the wavelength range covered by a single pixel.

    The bin set covers the wavelength range from lambda_min - 100 A to at least
    lambda_max + 100 A, where lambda_min and lambda_max are the minimum and maximum
    requested wavelengths. The extra 100 A are added to avoid artifacts at the
    boundaries of the requested wavelength range.

    Parameters
    ----------
    configuration: Configuration
        Simulator configuration.

    Returns
    -------
    Observation
        The observation for the sky background.
    """
    sky = sky_spectrum()
    grating = cast(Grating, configuration.telescope.grating)
    bandpass = (
        telescope_throughput()
        * fibre_throughput(
            seeing=configuration.seeing,
            source_extension="Diffuse",
            zenith_distance=configuration.zenith_distance,
        )
        * filter_transmission(filter_name=cast(Filter, configuration.telescope.filter))
        * grating_efficiency(
            grating_angle=grating.grating_angle,
            grating_name=grating.name,
        )
        * detector_quantum_efficiency()
    )
    return Observation(
        sky,
        bandpass,
        binset=_binset(grating.grating_angle, grating.grating_constant),
    )

snr(configuration)

Source code in nirwals/physics/exposure.py

def snr(configuration: Configuration) -> tuple[Quantity, Quantity]:
    # Get the source and sky observation.
    source = source_observation(configuration)
    sky = sky_observation(configuration)

    # Get the source and sky detection rates,
    grating = cast(Grating, configuration.telescope.grating)
    wavelengths_source, rates_source = detection_rates(
        area=configuration.telescope.effective_mirror_area,
        grating_angle=grating.grating_angle,
        grating_constant=grating.grating_constant,
        observation=source,
    )
    wavelengths_sky, rates_sky = detection_rates(
        area=configuration.telescope.effective_mirror_area,
        grating_angle=grating.grating_angle,
        grating_constant=grating.grating_constant,
        observation=sky,
    )

    # Collect the remaining relevant configuration parameters.
    exposure = cast(Exposure, configuration.exposure)
    e = exposure.exposures
    t = exposure.exposure_time
    detector = cast(Detector, configuration.detector)
    read_noise = detector.read_noise
    samplings = detector.samplings
    sampling_mode = detector.sampling_mode
    r = readout_noise(
        read_noise=read_noise, samplings=samplings, sampling_mode=sampling_mode
    )

    # Calculate the SNR.
    source_counts = (
        (rates_source * e * t).to(units.PHOTLAM * u.AA * u.cm**2 * u.s).value
    )
    sky_counts = (rates_sky * e * t).to(units.PHOTLAM * u.AA * u.cm**2 * u.s).value
    snr_values = source_counts / np.sqrt(source_counts + sky_counts + r * e)

    # Return the wavelengths and SNR values.
    return wavelengths_source, snr_values * u.dimensionless_unscaled

source_electrons(configuration)

Return the number of electrons accumulated due to source photons.

The electron counts are calculated for the same bins which are used by the detection_rates function. They should not be confused with the detector counts, for which you still would have to divide by the gain.

Parameters:

Name	Type	Description	Default
configuration	Configuration	Simulator configuration.	required

Returns:

Type	Description
tuple[Quantity, Quantity]	A tuple with an array of wavelengths and an array of the corresponding electron counts.

Source code in nirwals/physics/exposure.py

def source_electrons(configuration: Configuration) -> tuple[Quantity, Quantity]:
    """
    Return the number of electrons accumulated due to source photons.

    The electron counts are calculated for the same bins which are used by the
    detection_rates function. They should not be confused with the detector counts,
    for which you still would have to divide by the gain.

    Parameters
    ----------
    configuration: Configuration
        Simulator configuration.

    Returns
    -------
    tuple[Quantity, Quantity]
        A tuple with an array of wavelengths and an array of the corresponding electron
        counts.
    """
    # Collect the required parameters.
    area = configuration.telescope.effective_mirror_area
    exposure = cast(Exposure, configuration.exposure)
    grating = cast(Grating, configuration.telescope.grating)
    observation = source_observation(configuration)

    return electrons(
        area=area,
        exposures=exposure.exposures,
        exposure_time=exposure.exposure_time,
        grating_angle=grating.grating_angle,
        grating_constant=grating.grating_constant,
        observation=observation,
    )

source_observation(configuration)

Returns the source observation for a given configuration.

The observation, represented by a synphot Observation object, is the photon flux detected, which means it is the source spectrum with the following losses applied:

• Atmospheric extinction
• Telescope throughput
• Fibre throughput
• Filter transmission
• Grating efficiency
• Detector quantum efficiency

The bin set of the Observation object is chosen so that the bin just size is equal to the wavelength range covered by a single pixel.

The bin set covers the wavelength range from lambda_min - 100 A to at least lambda_max + 100 A, where lambda_min and lambda_max are the minimum and maximum requested wavelengths. The extra 100 A are added to avoid artifacts at the boundaries of the requested wavelength range.

Parameters:

Name	Type	Description	Default
configuration	Configuration	Simulator configuration.	required

Returns:

Type	Description
Observation	The observation for the source spectrum.

Source code in nirwals/physics/exposure.py

def source_observation(configuration: Configuration) -> Observation:
    """
    Returns the source observation for a given configuration.

    The observation, represented by a synphot Observation object, is the photon flux
    detected, which means it is the source spectrum with the following losses applied:

    - Atmospheric extinction
    - Telescope throughput
    - Fibre throughput
    - Filter transmission
    - Grating efficiency
    - Detector quantum efficiency

    The bin set of the Observation object is chosen so that the bin just size is equal
    to the wavelength range covered by a single pixel.

    The bin set covers the wavelength range from lambda_min - 100 A to at least
    lambda_max + 100 A, where lambda_min and lambda_max are the minimum and maximum
    requested wavelengths. The extra 100 A are added to avoid artifacts at the
    boundaries of the requested wavelength range.

    Parameters
    ----------
    configuration: Configuration
        Simulator configuration.

    Returns
    -------
    Observation
        The observation for the source spectrum.
    """
    source = source_spectrum(configuration=configuration)
    grating = cast(Grating, configuration.telescope.grating)
    configuration_source = cast(Source, configuration.source)
    bandpass = (
        atmospheric_transmission(zenith_distance=configuration.zenith_distance)
        * telescope_throughput()
        * fibre_throughput(
            seeing=configuration.seeing,
            source_extension=configuration_source.extension,
            zenith_distance=configuration.zenith_distance,
        )
        * filter_transmission(filter_name=cast(Filter, configuration.telescope.filter))
        * grating_efficiency(
            grating_angle=grating.grating_angle,
            grating_name=grating.name,
        )
        * detector_quantum_efficiency()
    )
    return Observation(
        source,
        bandpass,
        binset=_binset(
            grating_angle=grating.grating_angle,
            grating_constant=grating.grating_constant,
        ),
        force="taper"
    )

wavelength_resolution_element(grating_angle, grating_constant)

Return the wavelength resolution element for a grating setup.

Parameters:

Name	Type	Description	Default
grating_constant	micron	The grating constant, i.e. the groove spacing.	required
grating_angle	deg	The grating angle, i.e. the angle of the incoming rays to the grating normal.	required

Returns:

Type	Description
Quantity	The wavelength resolution element.

Source code in nirwals/physics/exposure.py

def wavelength_resolution_element(
    grating_angle: u.deg, grating_constant: u.micron
) -> u.AA:
    """
    Return the wavelength resolution element for a grating setup.

    Parameters
    ----------
    grating_constant: Quantity
        The grating constant, i.e. the groove spacing.
    grating_angle: Angle
        The grating angle, i.e. the angle of the incoming rays to the grating normal.

    Returns
    -------
    Quantity
        The wavelength resolution element.
    """
    # Angle of a fibre on the sky, in radians
    phi_fibre = 2 * FIBRE_RADIUS.to(u.rad).value

    return (
        phi_fibre
        * (TELESCOPE_FOCAL_LENGTH / COLLIMATOR_FOCAL_LENGTH)
        * grating_constant
        * math.cos(grating_angle.to(u.rad).value)
    )