%
%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}


