diff --git a/aiapy/response/channel.py b/aiapy/response/channel.py index da07c8c..ddc9e19 100644 --- a/aiapy/response/channel.py +++ b/aiapy/response/channel.py @@ -5,15 +5,15 @@ import collections from urllib.parse import urljoin -import astropy.constants as const import astropy.units as u import numpy as np +from sunkit_instruments.response.abstractions import AbstractChannel from sunpy.data import manager from sunpy.io.special import read_genx from sunpy.util.metadata import MetaDict +import aiapy.calibrate from aiapy import _SSW_MIRRORS -from aiapy.calibrate import degradation from aiapy.calibrate.util import _select_epoch_from_correction_table, get_correction_table from aiapy.util import telescope_number from aiapy.util.decorators import validate_channel @@ -38,7 +38,7 @@ } -class Channel: +class Channel(AbstractChannel): """ Interface to AIA channel properties and response functions. @@ -53,6 +53,21 @@ class Channel: instrument_file : `str`, optional Path to AIA instrument file. If not specified, the latest version will be downloaded from SolarSoft. + include_eve_correction : `bool`, optional + If true, include correction to EVE calibration in effective area. + include_crosstalk : `bool`, optional + If true, include the effect of crosstalk between channels that share a telescope + correction_table : `~astropy.table.Table` or `str`, optional + Table of correction parameters or path to correction table file. + If not specified, it will be queried from JSOC. + If you are calling this function repeatedly, it is recommended to + read the correction table once and pass it with this argument to avoid + multiple redundant network calls. + calibration_version : `int`, optional + The version of the calibration to use when calculating the + degradation. By default, this is the most recent version available + from JSOC. If you are using a specific calibration response file, + you may need to specify this according to the version in that file. Examples -------- @@ -69,9 +84,21 @@ class Channel: @u.quantity_input @validate_channel("channel") - def __init__(self, channel: u.angstrom, *, instrument_file=None): + def __init__( + self, + channel: u.angstrom, + *, + instrument_file=None, + include_eve_correction=False, + include_crosstalk=True, + **kwargs, + ): self._channel = channel self._instrument_data = self._get_instrument_data(instrument_file) + self.include_eve_correction = include_eve_correction + self.include_crosstalk = include_crosstalk + self.correction_table = kwargs.get("correction_table", None) + self.calibration_version = kwargs.get("calibration_version", None) @property def is_fuv(self): @@ -147,148 +174,120 @@ def wavelength( @property @u.quantity_input - def primary_reflectance( + def primary_mirror_reflectance( self, ) -> u.dimensionless_unscaled: return u.Quantity(self._data["primary"]) @property @u.quantity_input - def secondary_reflectance( + def secondary_mirror_reflectance( self, ) -> u.dimensionless_unscaled: return u.Quantity(self._data["secondary"]) @property @u.quantity_input - def focal_plane_filter_efficiency( + def mirror_reflectance(self): + """ + Combined reflectance of the primary and secondary mirrors. + """ + return self.primary_reflectance * self.secondary_reflectance + + @property + @u.quantity_input + def focal_plane_filter_transmittance( self, ) -> u.dimensionless_unscaled: return u.Quantity(self._data["fp_filter"]) @property @u.quantity_input - def entrance_filter_efficiency( + def entrance_filter_transmittance( self, ) -> u.dimensionless_unscaled: return u.Quantity(self._data["ent_filter"]) @property @u.quantity_input - def geometrical_collecting_area( + def filter_transmittance(self): + "Combined transmittance of the focal plane and entrance filters" + return self.entrance_filter_transmittance * self.focal_plane_filter_transmittance + + @property + @u.quantity_input + def geometrical_area( self, ) -> u.cm**2: return u.Quantity(self._data["geoarea"], u.cm**2) @property @u.quantity_input - def quantum_efficiency( + def effective_quantum_efficiency( self, ) -> u.dimensionless_unscaled: return u.Quantity(self._data["ccd"]) @property @u.quantity_input - def contamination( - self, - ) -> u.dimensionless_unscaled: - # Contamination missing for FUV channels - if "contam" in self._data: - return u.Quantity(self._data["contam"]) - return u.Quantity([1]) + def _preflight_contamination(self): + # NOTE: Preflight contamination calibration data missing for FUV channels + return u.Quantity(self._data.get("contam"), np.ones(self.wavelength.shape)) + + @u.quantity_input + def energy_per_electron(self) -> u.eV / u.electron: + return 1 / self._data["elecperev"] * u.eV / u.electron + + @u.quantity_input + def camera_gain(self) -> u.DN / u.electron: + return 1 / self._data["elecperdn"] * u.DN / u.electron - @property @u.quantity_input - def plate_scale( + def pixel_solid_angle( self, ) -> u.steradian / u.pixel: return u.Quantity(self._data["platescale"], u.steradian / u.pixel) - @property @u.quantity_input - def effective_area( - self, - ) -> u.cm**2: - r""" - Uncorrected effective area as a function of wavelength. - - According to Section 2 of [boerner]_, the effective area is given by, - - .. math:: - - A_{eff}(\lambda) = A_{geo}R_P(\lambda)R_S(\lambda)T_E(\lambda)T_F(\lambda)D(\lambda)Q(\lambda), - - where, - - - :math:`A_{geo}`: geometrical collecting area - - :math:`R_P`, :math:`R_S`: reflectances of primary and secondary mirrors, respectively - - :math:`T_E`, :math:`T_F`: transmission efficiency of the entrance and focal-plane filters, respectively - - :math:`D`: contaminant transmittance of optics - - :math:`Q`: quantum efficiency of the CCD - - The effective area contains information about the efficiency of the telescope optics and its sensitivity as a function of wavelength. - All of the telescope properties are read from the AIA instrument files available in SolarSoft. + def _get_eve_correction(self, obstime) -> u.dimensionless_unscaled: + """ + Calculate degradation factor from EVE rocket flight cross-calibration. - References - ---------- - .. [boerner] Boerner et al., 2012, Sol. Phys., `275, 41 `__ + This function is adapted directly from the + `aia_bp_corrections.pro `__ + routine in SolarSoft. """ - return ( - self.primary_reflectance - * self.secondary_reflectance - * self.focal_plane_filter_efficiency - * self.entrance_filter_efficiency - * self.geometrical_collecting_area - * self.quantum_efficiency - * self.contamination + table = _select_epoch_from_correction_table( + self.channel, + obstime, + get_correction_table(correction_table=self.correction_table), + version=self.calibration_version, ) + effective_area_interp = np.interp(table["EFF_WVLN"][-1], self.wavelength, self.effective_area) + return table["EFF_AREA"][0] / effective_area_interp - @property @u.quantity_input - def crosstalk( - self, - ) -> u.cm**2: - """ - Contamination of effective area from crosstalk between channels. + def degradation(self, obstime=None) -> u.dimensionless_unscaled: + r""" + Wavelength- and time-dependent degradation factor. - On telescopes 1, 3, and 4, both channels are always illuminated. - This can lead to contamination in a channel from the channel with which it shares a telescope. - This impacts the 94 and 304 Å channels as well as 131 and 335 Å. - See Section 2.2.1 of [1]_ for more details. + The AIA wavelength- and time-dependent degradation is modeled as, - References - ---------- - .. [1] Boerner et al., 2012, Sol. Phys., `275, 41 `__ - """ - crosstalk_lookup = { - 94 * u.angstrom: 304 * u.angstrom, - 304 * u.angstrom: 94 * u.angstrom, - 131 * u.angstrom: 335 * u.angstrom, - 335 * u.angstrom: 131 * u.angstrom, - } - if self.channel in crosstalk_lookup: - cross = Channel(crosstalk_lookup[self.channel], instrument_file=self._instrument_data) - return ( - cross.primary_reflectance - * cross.secondary_reflectance - * self.focal_plane_filter_efficiency - * cross.entrance_filter_efficiency - * cross.geometrical_collecting_area - * cross.quantum_efficiency - * cross.contamination - ) - return u.Quantity(np.zeros(self.wavelength.shape), u.cm**2) + .. math:: - @u.quantity_input - def eve_correction(self, obstime, **kwargs) -> u.dimensionless_unscaled: - r""" - Correct effective area to give good agreement with full-disk EVE data. + D(\lambda, t) = D_c(\lambda)C_e(t)C_d(t) + + where :math:`D_c(\lambda)` is the pre-flight contamination model + as described in Section 2.1.6 of [boerner], + :math:`C_e(t)` is the EVE cross-calibration correction factor, + and :math:`C_d(t)` is the time-varying telescope degradation computed + by `aiapy.calibrate.degradation`. The EVE correction factor is given by, .. math:: - \frac{A_{eff}(\lambda_n,t_0)}{A_{eff}(\lambda_E,t_e)} + C_e(t) = \frac{A_{eff}(\lambda_n,t_0)}{A_{eff}(\lambda_E,t_e)} where :math:`A_{eff}(\lambda_n,t_0)` is the effective area at the nominal wavelength of the channel (:math:`\lambda_n`) at the first @@ -296,25 +295,10 @@ def eve_correction(self, obstime, **kwargs) -> u.dimensionless_unscaled: area at the ```obstime``` calibration epoch interpolated to the effective wavelength (:math:`\lambda_E`). - .. note:: This function is adapted directly from the - `aia_bp_corrections.pro `__ - routine in SolarSoft. - Parameters ---------- obstime : `~astropy.time.Time` The time of the observation. - correction_table : `~astropy.table.Table` or `str`, optional - Table of correction parameters or path to correction table file. - If not specified, it will be queried from JSOC. - If you are calling this function repeatedly, it is recommended to - read the correction table once and pass it with this argument to avoid - multiple redundant network calls. - calibration_version : `int`, optional - The version of the calibration to use when calculating the - degradation. By default, this is the most recent version available - from JSOC. If you are using a specific calibration response file, - you may need to specify this according to the version in that file. Returns ------- @@ -322,106 +306,96 @@ def eve_correction(self, obstime, **kwargs) -> u.dimensionless_unscaled: See Also -------- - aiapy.calibrate.util.get_correction_table - """ - table = _select_epoch_from_correction_table( - self.channel, - obstime, - get_correction_table(correction_table=kwargs.get("correction_table")), - version=kwargs.get("calibration_version"), - ) - effective_area_interp = np.interp(table["EFF_WVLN"][-1], self.wavelength, self.effective_area) - return table["EFF_AREA"][0] / effective_area_interp - - @property - @u.quantity_input - def gain( - self, - ) -> u.DN / u.ph: - r""" - Gain of the CCD camera system. - - According to Section 2 of [boerner1]_, the gain of the CCD-camera system, in - DN per photon, is given by, - - .. math:: - - G(\lambda) = \frac{hc}{\lambda}\frac{g}{a} - - where :math:`g` is the camera gain in DN per electron - and :math:`\approx 3.65` eV per electron is a conversion factor. + aiapy.calibrate.degradation References ---------- - .. [boerner1] Boerner et al., 2012, Sol. Phys., `275, 41 `__ + .. [boerner] Boerner et al., 2012, Sol. Phys., `275, 41 `__ """ - _e = u.electron # Avoid rewriting u.electron a lot - electron_per_ev = self._data["elecperev"] * _e / u.eV - energy_per_photon = const.h * const.c / self.wavelength / u.ph - electron_per_photon = (electron_per_ev * energy_per_photon).to(_e / u.ph) - return electron_per_photon / (self._data["elecperdn"] * _e / u.DN) + contamination = self._preflight_contamination + if obstime is not None: + contamination *= aiapy.calibrate.degradation( + self.channel, + obstime, + correction_table=self.correction_table, + calibration_version=self.calibration_version, + ) + if self.include_eve_correction: + contamination *= self.eve_correction(obstime) + return contamination @u.quantity_input - def wavelength_response( + def _get_crosstalk(self, obstime) -> u.cm**2: + """ + Compute cross-talk component of effective area. + """ + crosstalk_lookup = { + 94 * u.angstrom: 304 * u.angstrom, + 304 * u.angstrom: 94 * u.angstrom, + 131 * u.angstrom: 335 * u.angstrom, + 335 * u.angstrom: 131 * u.angstrom, + } + if self.channel not in crosstalk_lookup: + return u.Quantity(np.zeros(self.wavelength.shape), u.cm**2) + cross = Channel(crosstalk_lookup[self.channel], instrument_file=self._instrument_data) + effective_area = ( + cross.geometrical_area + * cross.mirror_reflectance + * self.focal_plane_filter_transmittance + * cross.entrance_filter_transmittance + * cross.effective_quantum_efficiency + * cross._preflight_contamination # noqa: SLF001 + ) + # NOTE: This applies only the time-dependent component of the degradation for + # the nominal channel. The logic here is that the nominal time-dependent correction + # should be applied to the total (nominal + crosstalk) effective area since the + # cross-calibration is done on the whole channel postflight which includes the crosstalk + effective_area *= self.degradation(obstime=obstime) / self._preflight_contamination + return effective_area + + @u.quantity_input + def effective_area( self, *, obstime=None, - include_eve_correction=False, - include_crosstalk=True, - **kwargs, - ) -> u.DN / u.ph * u.cm**2: + ) -> u.cm**2: r""" - The wavelength response function is the product of the gain and the - effective area. + Effective area as a function of wavelength. - The wavelength response as a function of time and wavelength is given - by, + According to Section 2 of [boerner]_, the effective area is given by, .. math:: - R(\lambda,t) = (A_{eff}(\lambda) + A_{cross}(\lambda))G(\lambda)C_T(t)C_E(t) + A_{eff}(\lambda) = A_{geo}R_P(\lambda)R_S(\lambda)T_E(\lambda)T_F(\lambda)Q(\lambda)D(\lambda,t), where, - - :math:`A_{eff}(\lambda)` is the effective area as a function of - wavelength - - :math:`A_{cross}(\lambda)` is the effective area of the crosstalk - channel - - :math:`G(\lambda)` is the gain of the telescope - - :math:`C_T(t)` is the time-dependent correction factor for the - instrument degradation - - :math:`C_E(t)` is the time-dependent EVE correction factor + - :math:`A_{geo}`: geometrical collecting area + - :math:`R_P`, :math:`R_S`: reflectance of primary and secondary mirrors, respectively + - :math:`T_E`, :math:`T_F`: transmissitance of the entrance and focal-plane filters, respectively + - :math:`Q`: effective quantum efficiency of the CCD + - :math:`D`: degradation of optics, including preflight contaminant transmittance and time-dependent corrections - Parameters - ---------- - obstime : `~astropy.time.Time`, optional - If specified, a time-dependent correction is applied to account for degradation. - include_eve_correction : `bool`, optional - If true and ``obstime`` is not `None`, include correction to EVE calibration. - The time-dependent correction is also included. - include_crosstalk : `bool`, optional - If true, include the effect of crosstalk between channels that share a telescope - correction_table : `~astropy.table.Table` or `str`, optional - Table of correction parameters or path to correction table file. - If not specified, it will be queried from JSOC. - See `aiapy.calibrate.util.get_correction_table` for more information. + The effective area contains information about the efficiency of the telescope optics and its sensitivity + as a function of wavelength. + All of the telescope properties are read from the AIA instrument files available in SolarSoft. - Returns - ------- - `~astropy.units.Quantity` + On telescopes 1, 3, and 4, both channels are always illuminated. + This can lead to contamination in a channel from the channel with which it shares a telescope. + This impacts the 94 and 304 Å channels as well as 131 and 335 Å. + For these channels, additional component is added to the effective area to model the cross-talk + between these channels. + See Section 2.2.1 of [boerner]_ for more details. See Also -------- - effective_area - gain - crosstalk - eve_correction - aiapy.calibrate.degradation + degradation + + References + ---------- + .. [boerner] Boerner et al., 2012, Sol. Phys., `275, 41 `__ """ - eve_correction, time_correction = 1, 1 - if obstime is not None: - time_correction = degradation(self.channel, obstime, **kwargs) - if include_eve_correction: - eve_correction = self.eve_correction(obstime, **kwargs) - crosstalk = self.crosstalk if include_crosstalk else 0 * u.cm**2 - return (self.effective_area + crosstalk) * self.gain * time_correction * eve_correction + effective_area = super().effective_area(obstime=obstime) + if self.include_crosstalk: + effective_area += self._get_crosstalk(obstime) + return effective_area diff --git a/pyproject.toml b/pyproject.toml index a008b9f..31e1b34 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -42,6 +42,8 @@ classifiers = [ ] dependencies = [ 'sunpy[net,image,map]>=6.0', + # FIXME: Change to just sunkit-instruments once released or to main when merged + 'sunkit-instruments@git+https://github.com/namurphy/sunkit-instruments#egg=add-abstract-channel', ] [project.urls]