% %Paolo Grigis 2008/01/07 imported template %Modified by Yingna Su on 2011/08/25 %Modified by Mark Weber on 2012/01/06 %Modified by Mark Cheung on 2012/01/13 %Modified by Mark Weber on 2012/01/17 \documentclass[12pt,twoside]{article} \usepackage{fancyheadings,graphicx,setspace} \usepackage[pdftex]{color} \makeatletter \renewcommand\normalsize{% \@setfontsize\normalsize{13pt}{14.5pt}% \abovedisplayskip 12\p@ \@plus3\p@ \@minus7\p@ \abovedisplayshortskip \z@ \@plus3\p@ \belowdisplayshortskip 6.5\p@ \@plus3.5\p@ \@minus3\p@ \belowdisplayskip \abovedisplayskip \let\@listi\@listI} \normalsize \makeatother % Page format \oddsidemargin -0.00 in \evensidemargin -0.00 in \textwidth 6.50 in % US Letter is 8.5 X 11 inches, so 6.5 X 9 max for text \columnsep 0.25 in \topmargin 0.00 in \textheight 8.6 in % Fancy headings (put AFTER page format) \pagestyle{fancyplain} \lhead[\fancyplain{}{\thepage}]% {\fancyplain{}{\it AIA PSF Characterization}} \rhead[\fancyplain{}{\it AIA PSF Characterization}]% {\fancyplain{}{\thepage}} %\cfoot{} \headrulewidth 2pt % Handy figure magic. \renewcommand{\bottomfraction}{1.0} \renewcommand{\topfraction}{1.0} \renewcommand{\textfraction}{0.} % Mathematical symbols \newcommand{\vm}[1]{\ensuremath{\mathbf{#1}}} % vector bold, e.g. \vm{B} => \mathbf{B} \newcommand{\Cdot}{\boldsymbol{\cdot}} \newcommand{\dotprod}[2]{\ensuremath{#1\Cdot#2}} % dot product \newcommand{\crossprod}[2]{\ensuremath{#1\times#2}} % cross product \newcommand{\curl}[1]{\ensuremath{ % \crossprod{\boldsymbol{\nabla}}{#1} % }} % curl e.g. \curl{a} => nabla x a \newcommand{\del}{\ensuremath{\partial}} % partial derivative sign \newcommand{\diverg}[1]{\ensuremath{ % \dotprod{\boldsymbol{\nabla}}{#1} % }} % divergence e.g. \diverg{a} => nabla . a \newcommand{\grad}{\ensuremath{\boldsymbol{\nabla}}} % gradient \newcommand{\gradsq}{\ensuremath{\nabla^{2}}} % gradient squared \newcommand{\gradient}[1]{ % \ensuremath{\grad #1} % } % gradient e.g. \gradient{u} => \grad u \newcommand{\laplacian}[1]{\ensuremath{\gradsq {#1}}} % laplacian \newcommand{\dderiv}[2]{\ensuremath{(#1 \Cdot \boldsymbol{\nabla})#2}} % directional derivative e.g. \dderiv{a}{b} => (a . \nabla) b \newcommand{\sinc}{\ensuremath{\mathrm{sinc}}} % sinc function \newcommand{\Bf}{{\bf B}} \newcommand{\ep}{{\bf e}_\phi} \newcommand{\vf}{{\bf v}} \newcommand{\pa}{\partial} \newcommand{\Rs}{R_{\odot}} \newcommand{\er}{\mbox{erf}} \newcommand{\Bc}{B_c} \newcommand{\miss}[1]{{\it #1}} %%% JOURNALS THAT HAVE A COMMAND IN AASTEX \newcommand{\apj}{ApJ} \newcommand{\apjl}{ApJ Letters} \newcommand{\apjs}{ApJ Supplement} \newcommand{\nat}{Nature} \newcommand{\aap}{A\&A} \newcommand{\araa}{ARA\&A} \newcommand{\solphys}{Sol.~Phys.} \newcommand{\mnras}{MNRAS} \newcommand{\jgr}{Journ.~Geophys.~Res.} \begin{document} %\relax \mbox{} \vskip 2cm \centerline{\Huge \sffamily AIA PSF Characterization and Image Deconvolution} \vskip 1cm \centerline{\Large \sffamily Paolo Grigis, Yingna Su, and Mark Weber for the AIA team} \vskip 1cm \centerline{\Large \sffamily Version 2012-Jan-17} \vskip 3cm \centerline{\includegraphics[width=60mm]{./AIA_logo_medium.jpg}} \newpage %\tableofcontents % \section{Introduction} This document describes the calibration work performed to build a model for the Point Spread Function (PSF) for the telescopes on the Atmospheric Imaging Assembly (AIA) on SDO in the Extreme UltraViolet (EUV) channels (94 \AA, 131 \AA, 171 \AA, 193 \AA, 211 \AA, 304 \AA, 335 \AA). \section{Components of the PSF} EUV radiation incident on the AIA telescope encounters the following optical components. \begin{enumerate} \item \emph{Entrance Filter} (EF), a thin layer of Aluminum or Zirconium supported by a 70 lines/inch mesh. The main function of the EF is to absorb visible light and off-band EUV radiation. The EUV light is diffracted by the mesh. \item \emph{Primary and secondary mirrors}. Their multi-layer coating enables the reflection of the EUV light. The mirrors can scatter some of the incoming EUV radiation. \item \emph{Focal Plane Filter} (FPF), a thin layer of Aluminum or Zirconium supported by a 70 lines/inch mesh. The main function of the FPF is to absorb visible light and off-band EUV radiation. The EUV light is diffracted by the mesh. \item \emph{CCD}, where the EUV radiation is converted into electrons. The electrons can diffuse before reaching the potential well that traps them in each pixel. \end{enumerate} \begin{figure}[h!] \centering\includegraphics[width=10cm]{listcomp.png} \caption{PSF components} \label{listcomp} \end{figure} A schematic representation of these components is shown in Figure~\ref{listcomp}. The core of the model PSF is determined by the RMS (Root Mean Square) spot diameter ($R$) which includes: jitter, CCD pixelization, charge spreading, and scatter by the mirror surfaces and the filter-support meshes (for details see Table~\ref{tbl_errbudg}). The wings of the PSF are dominated by the filter diffraction pattern. The following sections examine those components in more detail. \begin{table}[t!] \begin{center} \caption{A sample image-resolution error budgets for the 171~\AA~channel. The effects are added in quadrature to give an estimate of the system-level optical performance. \label{tbl_errbudg}} \vspace{0.5cm} %\begin{minipage}[b]{\hsize} %\parbox[t]{12cm}{ \mbox{}\\ \begin{tabular}{lc}\hline\hline & Contribution to \\ & RMS Spot diameter \\ Item & (arcsec) \\ \hline Optical Prescription & 0.60 \\ Fabrication, alignment and assembly effects & 1.21\\ Launch shift effects & 0.10 \\ On-orbit thermal effects & 0.21 \\ Focus error (with on-orbit correction) & 0.10 \\ Jitter residual & 0.48 \\ Detector Pixelization & 0.48 \\ CCD Charge spreading & 0.80 \\ \hline On-orbit performance prediction & 1.73 \\ \hline \end{tabular} %\end{minipage} %\vspace{0.5cm} \end{center} \end{table} \section{Diffraction from the mesh} The mesh used in the EF and FPF has a density $\rho=70 \pm 1$ lines/inch and the fractional geometrical area not covered by the mesh (i.e. free for the passage of the EUV radiation) is $A=0.82\pm0.02$. One inch is $n=25400\,\mu\mathrm{m}$. That means that the spacing $d$ between two consecutive mesh wires is $$d=\frac{n}{\rho} \pm \frac{n}{\rho}\frac{\delta \rho}{\rho} =362.9\pm 5.2\,\mu\mathrm{m}\, .$$ The width $w$ of the mesh wires can be found by using $$ A=\left(1-\frac{w}{d}\right)^2 $$ and therefore $$ w=d\left(1-\sqrt{A}\right) \pm \sqrt{\left(1-\sqrt{A}\right)^2{\delta d}^2+\frac{d^2}{4A}{\delta A}^2}=34.3 \pm 4.0 \,\mu\mathrm{m}\, . $$ A graphical representation of those numbers is shown in Figure~\ref{meshpars}. \begin{figure}[h!] \centering\includegraphics[width=10cm]{meshpar.png} \caption{Mesh parameters used in both the entrance and focal plane filters} \label{meshpars} \end{figure} \subsection{Diffraction from a mesh} The wings of the PSF are dominated by the filter diffraction pattern. The diffraction pattern from a 1-dimensional grid consists of a set of \emph{diffraction spikes} with a regular spacing as a function of the angle $\theta$. Using the approximation $\sin\theta\simeq\theta$, the spikes are at positions $$\theta=m\lambda/d \label{pitch1}$$ for integer $m \ne 0$, where $\lambda$ is the wavelength. The intensity of the spikes is governed by the function: \begin{equation} I(\theta)=\sinc^2(\theta \frac{w}{\lambda})\, , \label{eqone} \end{equation} where \begin{equation} \sinc(x)=\frac{\sin(\pi x)}{\pi x}\, . \end{equation} This is shown in Figure~\ref{onedimpattern}. \begin{figure}[h!] \centering\includegraphics[width=12cm]{diffraction1diminfo.png} \caption{Diffraction pattern from a grid for a one dimensional case} \label{onedimpattern} \end{figure} The 2-dimensional case is analogous: there are 2 main arms, spaced at 90 degrees from each other, whose intensity profile is the same as in the 1-dimensional case. In addition there are weaker spikes between the arms. This can be neglected in first approximation. For AIA, the front filters consist of two distinct pieces. They are mounted at different angles of (nominally) 40$^\circ$ and 50$^\circ$, as shown in Figure~\ref{fig_meshcfg}. Therefore the diffraction will consist of 4 arms located (nominally) at 50$^\circ\!\!$, 40$^\circ\!\!$, $-40^\circ\!\!$, $-50^\circ\!\!$. The \emph{actual} values of the angles can be measured from AIA images directly. Table~\ref{tableangle} shows the results of these measurements. The spikes belonging to the diffraction patter from the EF are further diffracted by the FPF. Therefore we can describe the total diffraction effect as a convolution of the diffraction of the EF and the FPF. The FPF filter is located at a distance of 9.6 cm from the CCD. \begin{figure}[h!] \centering\includegraphics[width=10cm]{meshmount.png} \caption{Mesh configuration (Entrance Filter) \label{fig_meshcfg}} \label{meshmounting} \end{figure} \begin{table}[t!] % \centering % \caption{Measured angles} % \vspace{0.5cm} \begin{minipage}[b]{\hsize} \parbox[t]{12cm}{ \mbox{}\\ \begin{tabular}{rrrrrrrr}\hline\hline \rule[1.4ex]{0cm}{1.4ex}Channel & Angle 1$^*$ & Angle 2$^*$ & Angle 3$^*$ & Angle 4$^*$ & Spacing$^{*a}$ & Mesh pitch$^{*b}$ & Spacing$^\star$ \\ \hline \rule[1.4ex]{0cm}{1.4ex} 94 \AA & 49.81 & 40.16 & $-$40.28 & $-$49.92 & 8.99 & 70.66 & 0.207 \\%ok double checked numebrs ofr FPF 2011/06/29 \rule[1.4ex]{0cm}{1.4ex} 131 \AA & 50.27 & 40.17 & $-$39.70 & $-$49.95 & 12.37 & 69.77 & 0.289 \\ \rule[1.4ex]{0cm}{1.4ex} 171 \AA & 49.81 & 39.57 & $-$40.13 & $-$50.38 & 16.26 & 70.26 & 0.377 \\ \rule[1.4ex]{0cm}{1.4ex} 193 \AA & 49.82 & 39.57 & $-$40.12 & $-$50.37 & 18.39 & 70.40 & 0.425 \\ \rule[1.4ex]{0cm}{1.4ex} 211 \AA & 49.78 & 40.08 & $-$40.34 & $-$49.95 & 19.97 & 69.93 & 0.465 \\ \rule[1.4ex]{0cm}{1.4ex} 304 \AA & 49.76 & 40.18 & $-$40.14 & $-$49.90 & 28.87 & 70.17 & 0.670 \\ \rule[1.4ex]{0cm}{1.4ex} 335 \AA & 50.40 & 39.80 & $-$39.64 & $-$50.25 & 31.83 & 70.20 & 0.738 \\ \hline \end{tabular} $*$: Entrance Filter\\ $a$: in pixels\\ $b$: lines/inch, computed from the measured spacing in pixels\\ $\star$: Focal Plane Filter, assuming 70 lines/inch mesh pitch\\ \label{tableangle} } \end{minipage} \vspace{0.5cm} \end{table} \section{PSF Core} The PSF core distribution is modelled as a 2-dimensional Gaussian function of the radial distance ($r$) from the center, as given by \begin{equation} I(r,\theta)=I_0\exp\left(\frac{-r^2}{2\sigma^2}\right)\, . \end{equation} As mentioned in Section~2, the width of the PSF core is determined by several instrument characteristics which contribute to the size of the RMS spot diameter: jitter, CCD pixelization, charge spreading, and scatter by the mirror surfaces and the filter-support meshes. A complete list of the components contributing to the RMS spot diameter is shown in Table~\ref{tbl_errbudg}, which is taken from Table~7 in Boerner et al.\ (2011). The AIA PSF software ({\tt aia\_calc\_psf.pro} in the SolarSoft AIA tree), provides a keyword option \linebreak ({\tt /use\_preflightcore}) to directly use the RMS spot diameters to fit the Gaussian core. However, the values represented in the table are conservative estimates of the imaging performance based on pre-flight measurements. They were originally derived in order to verify that the instrument met its performance requirements, and thus should be viewed as upper limits on the spot diameter rather than as most-accurate predictions. For example, it is known from the GT-to-science telescope alignment shifts that the telescope mirrors did not move as much as was allocated, and thus the contribution to the PSF from optical misalignment is overstated. Therefore, the default behavior of {\tt aia\_calc\_psf.pro} (i.e., {\it not} using the {\tt /use\_preflightcore} keyword) is to use a somewhat narrower core that has been evaluated using deconvolutions of on-orbit images. Now we discuss how to use the RMS spot diameters to fit the Gaussian model of the PSF cores, as it is implemented for the {\tt /use\_preflightcore} keyword of {\tt aia\_calc\_psf.pro}. The standard deviation ($\mu$) of the Gaussian model ($I(r,\theta)$) is set equal to one-half of the RMS spot diameter ($\frac{1}{2}R$). For a 2D Gaussian, \begin{equation} \mu = \sqrt 2 \sigma\, . \end{equation} The pixel size for each EUV channel is 0.6 arcsec, so we may write \begin{equation} \sigma = \frac{5\sqrt 2R}{12}\, , \label{eq_sigma} \end{equation} in which $R$ is in units of arcseconds, and $\sigma$ is in units of pixels. The $\sigma$ for each EUV channel is listed in Table~\ref{tbl_rmsspot}. In the PSF software, we compute the PSF by doing a convolution of the two diffraction patterns from the focal plane and entrance filter meshes. Each of these two patterns is itself a convolution of a diffraction sinc function with its own Gaussian core ($\sigma_{1}$ and $\sigma_{2}$). When the two patterns are convolved together to generate the full PSF, the sigmas of the cores obey the relation $\sigma^2 = \sigma_1^2 + \sigma_2^2$. In particular, the convolution of two identical Gaussians with $\sigma_{1} = \sigma_{2}$ gives $\sigma=\sqrt{2}\sigma_{1}$. So the sigmas used in the program for each of the two patterns ($\sigma_{1}$) should be ${\sigma}/\sqrt{2}$ such that {\it after convolution of the two patterns together} we get the proper value of the full $\sigma$ given by Eq.~\ref{eq_sigma}. In the program calculating the PSF, the variable {\tt wx} is defined as $1/2\sigma^{2}_{1}$. Therefore, \begin{equation} {\tt wx} = \frac{72}{25R^{2}}\, . \end{equation} The {\tt wx} for each EUV channel is listed in Table~\ref{tbl_rmsspot}, under ``Gaussian width''. \begin{table}[t!] \begin{center} \caption{Predicted image resolution (RMS spot diameter, $R$), $\sigma$, and Gaussian width ($wx$) used in the program modeling PSF for seven AIA EUV channels. \label{tbl_rmsspot}} \vspace{0.5cm} %\begin{minipage}[b]{\hsize} %\parbox[t]{12cm}{ \mbox{}\\ \begin{tabular}{cccc}\hline\hline & Predicted resolution & $\sigma$ & Gaussian width \\ Channel & (arcsec) & (pixels) & (pixels) \\ \hline 94 \AA & 1.74 & 1.025 & 0.951 \\ 131 \AA & 1.67 & 0.984 & 1.033 \\ 171 \AA & 1.73 & 1.019 & 0.962 \\ 193 \AA & 1.38 & 0.813 & 1.512\\ 211 \AA & 1.55 & 0.913 & 1.199\\ 304 \AA & 1.52 & 0.896 & 1.247\\ 335 \AA & 1.73 & 1.019 & 0.962\\ \hline \end{tabular} %\end{minipage} %\vspace{0.5cm} \end{center} \end{table} \section{Model PSF} We build a model PSF including the following effects: 1) diffraction from the EF, 2) diffraction from the FPF, 3) RMS spot diameter. Other effects are not taken into account. The PSF is built in 3 steps: first, an array is created with a series of Gaussians arranged into the pattern expected from the diffraction by the EF. We use the measured angles from Table~\ref{tableangle} and the intensity profile of Eq.~\ref{eqone}, using the plate scale value of 0.6~arcsec~pix$^{-1}$ to convert angles to pixels. These Gaussians have a narrower width than the $\sigma$ defined in Eq.~\ref{eq_sigma}. A second array is created with a series of Gaussians arranged into the pattern expected from the diffraction by the FPF. The conversion between angles and pixels is done using a distance between the filter and the CCD of 9.6~cm. These Gaussians also have a narrower width than the $\sigma$ defined above. The two patterns are then convolved with each other, such that the convolution corresponds to the cumulative effect of the two diffractions, and the width of its Gaussians is the proper value of $\sigma$ as defined in Eq.~\ref{eq_sigma}. Figures~\ref{psflog} and \ref{psflogwide} show the PSF for the 211 \AA\ channel displayed on a logarithmic color scaling. \begin{figure}[h!] \centering\includegraphics[width=10cm]{psfsmall.png} \caption{Central part of the model PSF, displayed in a logarithmic color scale.} \label{psflog} \end{figure} \begin{figure}[h!] \centering\includegraphics[width=10cm]{psflarge.png} \caption{Model PSF displayed in a logarithmic color scale. Compared with Figure~\ref{psflog}, this is a wider view of the strongest parts only.} \label{psflogwide} \end{figure} \section{IMAGE DECONVOLUTION ALGORITHMS} There are different options available to deconvolve images. A popular algorithm is the Richardson-Lucy (RL) algorithm, which has been very successful in deconvolving Hubble Space Telescope images. The RL algorithm is an iterative procedure for recovering a latent image that has been blurred by a known point spread function. %(Richardson 1972) Pixels in the observed image can be represented in terms of the point spread function and the latent image as : \begin{equation} c_{i} = \sum_{j} p_{ij} u_{j}\, \end{equation} where $p_{ij}$ is the point spread function (the fraction of light coming from true location $j$ that is observed at position $i$), $u_{j}$ is the pixel value at location $j$ in the latent image, and $c_{i}$ is the observed value at pixel location $i$. The statistics are performed under the assumption that the $u_j$ are Poisson distributed, which is appropriate for photon noise in the data. The basic idea is to calculate the most likely $u_j$ given the observed $c_i$ and known $p_{ij}$. This leads to an equation for $u_j$ which can be solved iteratively according to : \begin{equation} {u}_{j}^{(t+1)} ={u}_j^{(t)} \sum_{i}\frac{c_{i}}{\hat{c}_{i}}p_{ij} \end{equation} where: \begin{equation} \hat{c}_{i} = \sum_{j}{u}_{j}^{(t)} p_{ij} \end{equation} It has been shown empirically that if this iteration converges, it converges to the maximum likelihood solution for $u_j$. %(Shepp 1982) \section{PSF software in SSW} %We provide 3 main routines in SSWIDL to yield the model PSF and perform the %deconvolution: For the deconvolution of level 1 AIA images, we provide two main routines in the SSW AIA package: \begin{enumerate} \item {\ttfamily aia\_calc\_psf} returns PSFs for the different EUV channels, and\\ \item {\ttfamily aia\_deconvolve\_richardsonlucy} performs iterative RL deconvolution %\item {\ttfamily pg\_aia\_getmeshinfo} returns the mesh parameters\\ %\item {\ttfamily pg\_aia\_diffractionpattern } builds the PSF in different channels\\ %\item {\ttfamily pg\_richardsonlucy\_deconv} performs RL deconvolution\\ \end{enumerate} \noindent Here's an example code: %psf_parameters=pg_aia_getmeshinfo(wavelength) %;this returns a structure with the information %;needed to build a PSF model \begin{verbatim} wavelength=`171' ;this computes the PSF model from the information ;structure - be patient ;This version uses a PSF core within the error budgets and optimized ;for on-orbit performance. psf=aia_calc_psf(wavelength) ;This version uses a PSF core that exactly uses the channel error budgets ;like Table 1. psf=aia_calc_psf(wavelength, /use_preflightcore) ;perform deconvolution - assume `image' is a level 1 ;4096x4096 AIA image of the same wavelength ;as the psf ; cast image to a floating point array image = float(image) ;Now do RL deconvolution - be even more patient ;niter is an optional keyword. The default is 25 iterations. deconvolved_image=aia_richardsonlucy_deconvolve(image,psf,niter=25) \end{verbatim} If you want to create level 1.6 images, you should perform deconvolution from level 1 data and then use aia\_prep. \begin{thebibliography}{} \bibitem{Boe11} Boerner, P., Edwards, C., Lemen, J. et al.\ 2011, ``Initial Calibration of the Atmospheric Imaging Assembly (AIA) on the Solar Dynamics Observatory (SDO)", Solar Physics, submitted. \end{thebibliography}{} \end{document}