\documentclass[11pt]{article}       
\usepackage{a4wide,amsmath}
\usepackage{graphicx}
\usepackage{chicago}
\bibliographystyle{aa}
\def\aap{A\&A}
\def\apj{ApJ}
\def\pasp{PASP}
\def\mnras{MNRAS}
\def\icarus{ICARUS}
\def\araa{ARA\&A}
\def\apjs{ApJS}

\setlength{\textheight}{25cm}
\setlength{\textwidth}{15cm}
\oddsidemargin=0mm
\evensidemargin=0mm
%\setlength{\parindent}{0mm}
\setlength{\itemsep}{-5pt plus-5pt minus-5pt}
\setlength{\parskip}{4pt plus1pt minus1pt}
\setlength{\parsep}{2pt plus1pt minus1pt}
\voffset -1.5cm

    \def\abl#1#2{\frac {d #1}{d #2}}
    \def\pabl#1#2{\frac {\partial #1}{\partial #2}}
    \def\eg{e.\,g.\ }
    \def\ie{i.\,e.\ }

\begin{document}

\sloppy

\centerline{\huge\bf For Rosina}
{\ }\\[-1.7ex]
\centerline{\large Peter Woitke, June 2013}

\section{A flux-conserving scheme to convert images in polar coordinates
  to regular carthesian grids}


Coninuum images and channel maps in ProDiMo are initially calculated
on a polar grid (see Thi et al.~2011\nocite{Thi2011} for details) with
roughly\footnote{Some adjustments inside the inner rim and increased
  resolution towards the outer disk radius and beyond.}  logarithmic
equidistant radial grid points $\{r_i\,|\,i\!\in\![0...N_r]\}$ and linear
equidistant angular grid points $\{\theta_j\,|\,j\!\in\![0...N_\theta]\}$.
This is necessary to resolve the tiny inner disk rim (important for
the short wavelengths) as well as resolving the outer regions
(important for the long wavelengths) at the same time. The polar intensities
$I_\nu(i,j)$ have associated areas in the image plane as
\begin{equation}
  A_{ij} = \pi(r_i^2-r_{i-1}^2)/N_\theta \qquad\mbox{for\ }
                                       i\!\in\![1...N_r] ,
                                       j\!\in\![1...N_\theta]
\end{equation}
The intensities are assumed to be constant on the polar
``pixels'', \ie the area $A_{ij}$ bracketed by $r_i$, $r_{i-1}$,
$\theta_j$ and $\theta_{j-1}$, see Fig.~\ref{GridConvert}.

\begin{figure}
\centering
\includegraphics[width=85mm]{GridConvert2.eps}
\caption{Conversion from polar intensities $I_\nu(r,\theta)$ to a
  regular carthesian grid $I'_\nu(x,y)$. The $r_i$ and $\theta_j$
  gridpoints are drawn as black lines, the orange box represents a
  regular image pixel between carthesian grid points $x_n$ and
  $y_k$. The dashed box shows the construction of a rectangle that
  exactly contains the polar cell $(i,j)$.}
\label{GridConvert}
\end{figure}

We seek a fast numerical method to convert these polar images 
onto a regular grid with equidistant carthesian image coordinates 
$\{x_n\,|\,n\!\in\![0...N_x]\}$ and $\{y_k\,|\,k\!\in\![0...N_y]\}$
with associated pixel areas 
\begin{equation}
  A_{nk} = (x_n-x_{n-1})(y_k-y_{k-1}) \qquad\mbox{for\ } n\!\in\![1...N_x], 
                                               k\!\in\![1...N_y]
\end{equation}
The method is flux-conservative if
\begin{equation}
   \sum_{i=1}^{N_r}\sum_{j=1}^{N_\theta} A_{ij}\,I_\nu(i,j)
 = \sum_{n=1}^{N_x}\sum_{k=1}^{N_y} A_{nk}\,I'_\nu(n,k)
\end{equation}
where $I'_\nu(n,k)$ are the desired intensities on the regular
carthesian pixels.  The exact solution of this problem would be to
calculate the overlap areas, $O(A_{nk},A_{ij})$, between the regular
pixels $A_{nk}$ and any polar pixel $A_{ij}$, then sum up the fluxes 
and devide by the pixel area as
\begin{equation}
  I'_\nu(n,k) = \sum_{i=1}^{N_r}\sum_{j=1}^{N_\theta} 
                \frac{O(A_{nk},A_{ij})}{A_{nk}} \,I_\nu(i,j) \ ,
\end{equation}
but to calculate those overlap areas is painful, see
Fig.~\ref{GridConvert}. A more practical idea is to create a
carthesian rectangular area that exactly contains the polar pixel, see
Fig.~\ref{GridConvert}, by taking the minimum and maximum of
the four corner points of $A_{ij}$
\begin{eqnarray}
  \begin{array}{rp{1mm}lp{1mm}rp{1mm}l}
  x_l &=& \min\{x_1,x_2,x_3,x_4\} && y_l &=& \min\{y_1,y_2,y_3,y_4\} \\
  x_r &=& \max\{x_1,x_2,x_3,x_4\} && y_r &=& \max\{y_1,y_2,y_3,y_4\} \\
  x_1 &=& r_{i-1}\sin(\theta_{j-1}) && y_1 &=& r_{i-1}\cos(\theta_{j-1}) \\
  x_2 &=& r_{i-1}\sin(\theta_j)    && y_2 &=& r_{i-1}\cos(\theta_j) \\
  x_3 &=& r_i\sin(\theta_{j-1})    && y_3 &=& r_i\cos(\theta_{j-1}) \\
  x_4 &=& r_i\sin(\theta_j)        && y_4 &=& r_i\cos(\theta_j) 
  \end{array}
  \label{ConstructPixel}
\end{eqnarray}
The area of this rectangular pixel, 
\begin{equation}
  A'_{ij} = (x_r-x_l)(y_r-y_l) \ ,
\end{equation}
is always larger than the area of the original polar pixel $A_{ij}$,
resulting in a correction factor in Eq.~(\ref{PixelConvert})
below. The area overlaps between $A_{nk}$ and $A'_{ij}$ are now easy
to calculate, and the resulting conversion formula is
\begin{eqnarray}
  C_{nk} &=& \sum_{i=1}^{N_r}\sum_{j=1}^{N_\theta} 
             O(A_{nk},A'_{ij})\frac{A_{ij}}{A'_{ij}}\nonumber\\
  I'_\nu(n,k) &=& \frac{1}{C_{nk}}\sum_{i=1}^{N_r}\sum_{j=1}^{N_\theta} 
             O(A_{nk},A'_{ij})\frac{A_{ij}}{A'_{ij}}\,I_\nu(i,j) \ ,
  \label{PixelConvert}
\end{eqnarray}
The area coverage $C_{nk}$ is close, but not exactly equal to
$A_{nk}$, because some regular pixels will be slightly oversampled 
by the polar pixels, and others slightly undersampled, and this 
effect is automatically corrected for in Eq.~(\ref{PixelConvert}).

This conversion method is fast, reliable, and exactly flux
conservative.  However, the resulting spatial resolution in the
regular image will suffer somewhat in areas where the rectangular
pixels $A'_{ij}$, as constructed from the polar pixels $A_{ij}$, are
larger than the regular pixels $A_{nk}$.  Here, the simulation can be
improved by introducing polar sub-pixels in the first place, i.e.\ by
subdividing large polar pixels into a suitable number of smaller polar
sub-pixels, equally spaced in $r$ and $\theta$, assuming the
intensities to be constant over all sub-pixels, before applying
Eqs.~(\ref{ConstructPixel}) to (\ref{PixelConvert}).

\bibliography{references}

\end{document}