"""
Calculate the point spread function (PSF) for the AIA telescopes.
"""
import astropy.units as u
import numpy as np
from sunpy import log
from aiapy.util.decorators import validate_channel
try:
import cupy
HAS_CUPY = True
except ImportError:
HAS_CUPY = False
__all__ = ["psf", "filter_mesh_parameters", "_psf"]
[docs]
def filter_mesh_parameters(*, use_preflightcore=False):
"""
Geometric parameters for meshes in AIA filters used to calculate the point
spread function.
Parameters
----------
use_preflightcore : `bool`, optional
If True, use the pre-flight values for the filter mesh parameters
Returns
-------
meshinfo : `dict`
Dictionary with filter mesh information for each channel. Each channel
entry then contains another dictionary with the following keys
describing filter mesh properties of that channel
(see Table 2 of [1]_):
* ``angle_arm``: Angles of the four entrance filter arms
* ``error_angle_arm``: Error in angle of the four entrance filter arms
* ``spacing_e``: Distance between diffraction spikes from entrance filter
* ``spacing_fp``: Distance between diffraction spikes from focal plane filter
* ``mesh_pitch``: Pitch of the mesh
* ``mesh_width``: Width of the mesh
* ``width``: Width applied to the Gaussian such that *after*
convolution we have the proper width
(:math:`4/3` at :math:`1/e` of max)
References
----------
.. [1] `Grigis, P., Su, Y., Weber M., et al., 2012,
AIA PSF Characterization and Deconvolution
<https://sohoftp.nascom.nasa.gov/solarsoft/sdo/aia/idl/psf/DOC/psfreport.pdf>`_
See Also
--------
psf : Calculate the composite point spread function
Notes
-----
These parameters were calculated from the following images and
reference background images.
94:
* image: 'AIA20101016_191039_0094.fits'
* reference: 'AIA20101016_190903_0094.fits'
131:
* image: 'AIA20101016_191035_0131.fits'
* reference: 'AIA20101016_190911_0131.fits'
171:
* image: 'AIA20101016_191037_0171.fits'
* reference: 'AIA20101016_190901_0171.fits'
193:
* image: 'AIA20101016_191056_0193.fits'
* reference: 'AIA20101016_190844_0193.fits'
211:
* image: 'AIA20101016_191038_0211.fits'
* reference: 'AIA20101016_190902_0211.fits'
304:
* image: 'AIA20101016_191021_0304.fits'
* reference: 'AIA20101016_190845_0304.fits'
335:
* image: 'AIA20101016_191041_0335.fits'
* reference: 'AIA20101016_190905_0335.fits'
"""
return {
94 * u.angstrom: {
"angle_arm": [49.81, 40.16, -40.28, -49.92] * u.deg,
"error_angle_arm": [0.02, 0.02, 0.02, 0.02] * u.deg,
"spacing_e": 8.99 * u.pixel,
"mesh_pitch": 363.0 * u.um,
"mesh_width": 34.0 * u.um,
"spacing_fp": 0.207 * u.pixel,
"width": (0.951 if use_preflightcore else 4.5) * u.pixel,
"CDELT": [0.600109, 0.600109] * u.arcsec,
},
131 * u.angstrom: {
"angle_arm": [50.27, 40.17, -39.70, -49.95] * u.deg,
"error_angle_arm": [0.02, 0.02, 0.02, 0.02] * u.deg,
"spacing_e": 12.37 * u.pixel,
"mesh_pitch": 363.0 * u.um,
"mesh_width": 34.0 * u.um,
"spacing_fp": 0.289 * u.pixel,
"width": (1.033 if use_preflightcore else 4.5) * u.pixel,
"CDELT": [0.600698, 0.600698] * u.arcsec,
},
171 * u.angstrom: {
"angle_arm": [49.81, 39.57, -40.13, -50.38] * u.deg,
"error_angle_arm": [0.02, 0.02, 0.02, 0.02] * u.deg,
"spacing_e": 16.26 * u.pixel,
"mesh_pitch": 363.0 * u.um,
"mesh_width": 34.0 * u.um,
"spacing_fp": 0.377 * u.pixel,
"width": (0.962 if use_preflightcore else 4.5) * u.pixel,
"CDELT": [0.599489, 0.599489] * u.arcsec,
},
193 * u.angstrom: {
"angle_arm": [49.82, 39.57, -40.12, -50.37] * u.deg,
"error_angle_arm": [0.02, 0.02, 0.03, 0.04] * u.deg,
"spacing_e": 18.39 * u.pixel,
"mesh_pitch": 363.0 * u.um,
"mesh_width": 34.0 * u.um,
"spacing_fp": 0.425 * u.pixel,
"width": (1.512 if use_preflightcore else 4.5) * u.pixel,
"CDELT": [0.600758, 0.600758] * u.arcsec,
},
211 * u.angstrom: {
"angle_arm": [49.78, 40.08, -40.34, -49.95] * u.deg,
"error_angle_arm": [0.02, 0.02, 0.02, 0.02] * u.deg,
"spacing_e": 19.97 * u.pixel,
"mesh_pitch": 363.0 * u.um,
"mesh_width": 34.0 * u.um,
"spacing_fp": 0.465 * u.pixel,
"width": (1.199 if use_preflightcore else 4.5) * u.pixel,
"CDELT": [0.600758, 0.600758] * u.arcsec,
},
304 * u.angstrom: {
"angle_arm": [49.76, 40.18, -40.14, -49.90] * u.degree,
"error_angle_arm": [0.02, 0.02, 0.02, 0.02] * u.deg,
"spacing_e": 28.87 * u.pixel,
"mesh_pitch": 363.0 * u.um,
"mesh_width": 34.0 * u.um,
"spacing_fp": 0.670 * u.pixel,
"width": (1.247 if use_preflightcore else 4.5) * u.pixel,
"CDELT": [0.600165, 0.600165] * u.arcsec,
},
335 * u.angstrom: {
"angle_arm": [50.40, 39.80, -39.64, -50.25] * u.degree,
"error_angle_arm": [0.02, 0.02, 0.02, 0.02] * u.deg,
"spacing_e": 31.83 * u.pixel,
"mesh_pitch": 363.0 * u.um,
"mesh_width": 34.0 * u.um,
"spacing_fp": 0.738 * u.pixel,
"width": (0.962 if use_preflightcore else 4.5) * u.pixel,
"CDELT": [0.600737, 0.600737] * u.arcsec,
},
}
[docs]
@u.quantity_input
@validate_channel("channel", valid_channels=[94, 131, 171, 193, 211, 304, 335] * u.angstrom)
def psf(channel: u.angstrom, *, use_preflightcore=False, diffraction_orders=None, use_gpu=True):
r"""
Calculate the composite PSF for a given channel, including diffraction and
core effects.
.. note:: This function has been adapted from
`aia_calc_psf.pro <https://sohoftp.nascom.nasa.gov/solarsoft/sdo/aia/idl/psf/PRO/aia_calc_psf.pro>`_.
.. note:: If the `~cupy` package is installed
and your machine has an NVIDIA GPU, the PSF calculation will
automatically be accelerated with CUDA. This can lead to
several orders of magnitude in performance increase compared to
pure `numpy` on a CPU.
The point spread function (PSF) can be modeled as a 2D Gaussian function
of the radial distance :math:`r` from the center,
.. math::
I(r, \theta) = I_0 \exp\left(\frac{-r^2}{2\sigma^2}\right)
where,
- :math:`I_0` : the intensity of a diffraction spike
- :math:`r` : the radial distance from the center
- :math:`\theta = m\lambda/d`
- :math:`m` : diffraction order
- :math:`\lambda` : the wavelength of light
- :math:`\sigma` : width of Gaussian
The intensity of a particular diffraction spike, :math:`I_0`, is given by,
.. math::
I_0 = \mathrm{sinc}^2\left(\frac{\theta w}{\lambda}\right)
where,
- :math:`w` : the width of the mesh wires
- :math:`d` : spacing between two consecutive mesh wires
The PSF for a given filter can then be calculated as,
.. math::
\mathrm{PSF} = \sum_{m=-\infty}^{+\infty}I_m(r,\theta)
where, in practice, one can approximate the summation by simply summing
over a sufficiently large number of diffraction orders. In this case, we
sum from :math:`m=--100` to :math:`m=100`.
Finally, the composite PSF of the entrance and focal plane filters is
given by,
.. math::
\mathrm{PSF}_c = \left|\mathcal{F}\left\{
\mathcal{F}\{\mathrm{PSF}_f\}
\mathcal{F}\{\mathrm{PSF}_e\}
\right\}\right|
where :math:`\mathcal{F}` denotes the Fourier transform,
:math:`\mathrm{PSF}_f` is the PSF of the focal plane filter, and
:math:`\mathrm{PSF}_e` is the PSF of the entrance filter. For a more
detailed explanation of the PSF and the above calculation, see [1]_.
Parameters
----------
channel : `~astropy.units.Quantity`
Wavelength of channel
use_preflightcore : `bool`, optional
If True, use the pre-flight values of the mesh width
diffraction_orders : array-like, optional
The diffraction orders to sum over. If None, the full
range from -100 to +100 in steps of 1 will be used.
use_gpu : `bool`, optional
If True and `~cupy` is installed, do PSF deconvolution on the GPU
with `~cupy`.
Returns
-------
`~numpy.ndarray`
The composite PSF of the entrance and focal plane filters.
See Also
--------
filter_mesh_parameters
deconvolve
References
----------
.. [1] `Grigis, P., Su, Y., Weber M., et al., 2012,
AIA PSF Characterization and Deconvolution
<https://sohoftp.nascom.nasa.gov/solarsoft/sdo/aia/idl/psf/DOC/psfreport.pdf>`__
"""
meshinfo = filter_mesh_parameters(use_preflightcore=use_preflightcore)
meshinfo = meshinfo[channel]
angles_entrance = meshinfo["angle_arm"]
angles_focal_plane = u.Quantity([45.0, -45.0], "deg")
if diffraction_orders is None:
diffraction_orders = np.arange(-100, 101, 1)
psf_entrance = _psf(meshinfo, angles_entrance, diffraction_orders, use_gpu=use_gpu)
psf_focal_plane = _psf(
meshinfo,
angles_focal_plane,
diffraction_orders,
focal_plane=True,
use_gpu=use_gpu,
)
# Composite PSF
psf = abs(np.fft.fft2(np.fft.fft2(psf_focal_plane) * np.fft.fft2(psf_entrance)))
# Center PSF in the middle of the image
psf = np.roll(np.roll(psf, psf.shape[1] // 2, axis=1), psf.shape[0] // 2, axis=0)
# Normalize by total number of pixels
psf = psf / (psf.shape[0] * psf.shape[1])
# If using cupy, cast back to a normal numpy array
if HAS_CUPY and use_gpu:
psf = cupy.asnumpy(psf)
return psf
def _psf(meshinfo, angles, diffraction_orders, *, focal_plane=False, use_gpu=True):
psf = np.zeros((4096, 4096), dtype=float)
if use_gpu and not HAS_CUPY:
log.info("cupy not installed or working, falling back to CPU")
# If cupy is available, cast to a cupy array
if HAS_CUPY and use_gpu:
psf = cupy.array(psf)
Nx, Ny = psf.shape
width_x = meshinfo["width"].value
width_y = meshinfo["width"].value
# x and y position grids
x = np.outer(np.ones(Ny), np.arange(Nx) + 0.5)
y = np.outer(np.arange(Ny) + 0.5, np.ones(Nx))
if HAS_CUPY and use_gpu:
x = cupy.array(x)
y = cupy.array(y)
area_not_mesh = 0.82 # fractional area not covered by the mesh
spacing = meshinfo["spacing_fp"] if focal_plane else meshinfo["spacing_e"]
mesh_ratio = (meshinfo["mesh_pitch"] / meshinfo["mesh_width"]).decompose().value
spacing_x = spacing * np.cos(angles)
spacing_y = spacing * np.sin(angles)
for order in diffraction_orders:
if order == 0:
continue
intensity = np.sinc(order / mesh_ratio) ** 2 # I_0
for dx, dy in zip(spacing_x.value, spacing_y.value, strict=True):
x_centered = x - (0.5 * Nx + dx * order + 0.5)
y_centered = y - (0.5 * Ny + dy * order + 0.5)
# NOTE: this step is the bottleneck and is VERY slow on a CPU
psf += np.exp(-width_x * x_centered * x_centered - width_y * y_centered * y_centered) * intensity
# Contribution from core
psf_core = np.exp(-width_x * (x - 0.5 * Nx - 0.5) ** 2 - width_y * (y - 0.5 * Ny - 0.5) ** 2)
return (1 - area_not_mesh) * psf / psf.sum() + area_not_mesh * psf_core / psf_core.sum()
