%
% To start the document, use
%  \levela{...}
% For lover level, sections use
%  \levelb{...}
%  \levelc{...}
%
\levela{Scattering}
 \label{sec:scattering}

%
% Document history, format:
%  \starthistory
%    date1 & text .... \\
%    date2 & text .... \\
%    ....
%  \stophistory
%
 \starthistory
 100502 & Created and written by Claudia Emde.\\
 120103 & Updated.
 \stophistory


%
% Symbol table, format:
%  \startsymbols
%    ... & \artsstyle{...} & text ... \\
%    ... & \artsstyle{...} & text ... \\
%    ....
%  \stopsymbols
%
 \startsymbols
 \StoVec       & \artsstyle{i\_field}       & monochromatic intensity field\/ Stokes Vector\\
 & \artsstyle{i\_field\_old}  & needed for the iteration\\
 i$_I$         & -                          & Stokes component \\
 \ExtMat       & -                        & total extinction matrix \\
 \SExMat & \artsstyle{ext\_mat} & extinction matrix for a
 single particle type\\
 \AbsVec       & -                        & total absorption vector \\
 \SAbVec & \artsstyle{abs\_vec} & absorption vector for a
 single particle type\\
 \SEfMat & - & scattering efficiency
 matrix\\
 \PhaMat & \artsstyle{pha\_mat} & phase matrix for a particle
 type\\
 \AmpMat       & \artsstyle{amp\_mat} & amplitude matrix\\
 \PDir         & -                        & propagation direction \\
 \Frq          & \artsstyle{f\_grid}       & frequency grid\\
 \Tmp          & \artsstyle{t\_field}         & temperature field\\
 \PDen         & \artsstyle{pnd\_field} & particle density field \\
 \ScaInt & \artsstyle{sca\_vec} & scattered vector,
 evaluated scattering integral\\
 \Prs & \artsstyle{p\_grid} & pressure grid inside cloud
 box\\
 \Lat & \artsstyle{lat\_grid} & latitude grid inside cloud
 box\\
 \Lon & \artsstyle{lon\_grid} & longitude grid inside cloud
 box\\
 \ScaZa        & \artsstyle{scat\_za\_grid}  & zenith angle \\
 \ScaAa        & \artsstyle{scat\_aa\_grid}  & azimuth angle  \\
 & \artsstyle{scat\_f\_index}  & frequency index\\
 \Planck & \artsstyle{a\_planck\_value} & Planck function
\label{symtable:scattering}
\stopsymbols
%=====================================================================


%=====================================================================
% Definition of new commands:
% ====================================================================
\newcommand{\DirFre} {(\PDir, \Frq, \Tmp)} \newcommand{\DirFrePr}
{\ensuremath{(\PDir, \PDir^\prime, \Frq, \Tmp)}}


\levelb{Vector radiative transfer equation (VRTE)}
%=====================================================================
\label{sec:scattering:general_rte}
 
The radiative transfer equation for a medium with thermal emission
comprising sparsely and randomly distributed, spherical or arbitrarily
oriented non-spherical particles is according to \cite{sree02}:
\begin{eqnarray}
\label{eq:scattering:RTE} 
     \frac{\DiffD \StoVec}{\DiffD s}(\PDir, \Frq) =
     -\ExtMat\DirFre\StoVec(\PDir, \Frq)+\AbsVec\DirFre
     \Planck(\Frq) \\ \nonumber
     +\int_{4\pi} \DiffD \PDir^\prime \SEfMat\DirFrePr
     \StoVec(\PDir^\prime, \Frq) 
\end{eqnarray} 
Formally it is an inhomogeneous four dimensional vector differential
equation for the Stokes vector \StoVec $(I,Q,U,V)$. To describe the
full state of an electromagnetic wave one could use the
electromagnetic field vector, which consists of two complex numbers.
This is the common notation in Maxwell's theory. An alternative
notation is the Stokes vector which is more convenient for
experimental purposes, as it consists of real numbers which can be
associated more directly with the measurements.  The first component
of the Stokes vector $I$ signifies the monochromatic flux or the
intensity of the radiation. The other components describe the
polarization state: the second and third components $Q$ and $U$
characterize the state of linear polarization and the last component
$V$
characterizes the state of circular polarization.\\
In Equation (\ref{eq:scattering:RTE}) \PDir\ denotes the propagation
direction and \DiffD $s$ the path-length element along this direction.
The radiation field depends on position, propagation direction {\bf
  n}, and frequency \Frq.  \vspace{1ex} The first term of Equation
(\ref{eq:scattering:RTE}) corresponds to the extinction of radiation
determined by the extinction coefficient matrix \ExtMat . For
microwave radiation traveling through the atmosphere there is
extinction due to gaseous absorption, particle absorption and particle
scattering. Lambert's law states that the extinction process is
proportional to the amount of matter. Therefore \ExtMat\ can be
written as a sum of two matrices, the first for particle extinction (
\aExtMat{p}) and the second for gaseous extinction (\aExtMat{g}):
\begin{eqnarray}
  \ExtMat\DirFre &=&
  \aExtMat{p}\DirFre+\aExtMat{g}\DirFre
\end{eqnarray}

The particle extinction term can be written as a sum, each term
corresponding to one particle type:
\begin{eqnarray}
  \aExtMat{p}\DirFre = \sum_i \PDen_i \aSExMat{p}_i \DirFre
\end{eqnarray}

Here $\PDen_i$ is the particle number density of the {\sl i}th
particle type and $\aSExMat{p}_i\DirFre$ the particle extinction cross
section matrix of the {\sl i}th particle type.


The gaseous extinction matrix is derived from the scalar gas
absorption. As extinction by gas absorption does not cause any
polarization effects, the off-diagonal elements of the gaseous
extinction matrix have to be zero. Gaseous extinction is determined
only by absorption, so the coefficients on the diagonal correspond to
the gaseous absorption coeffient ({\bf FIXME: Why are the coefficients
  equal, may be only K11 = a, and all others 0 ?}):

\[
\aExtMat{g}_{i,j}(\Prs, \Frq, \Tmp) =
\left\{\begin{array}{r@{\quad:\quad}l} \alpha(\Prs, \Frq, \Tmp) & i =
    j \\ 0 & i \not= j \end{array} \right.
\]

Here \Prs and \Tmp are pressure and temperature. $\alpha$ is the
scalar gaseous absorption coefficient.

The second term in Equation (\ref{eq:scattering:RTE}) describes the
thermal emission. \Planck\ is the Planck function and \AbsVec\ is the
absorption coefficient vector which can be written as
\begin{eqnarray}
  \AbsVec \DirFre  &=& \aAbsVec{p} \DirFre + \aAbsVec{g} \DirFre 
\end{eqnarray}
with \aAbsVec{p} and \aAbsVec{g} as the particle absorption
coefficient and the gaseous absorption coefficient, respectively.

The particle absorption term can also be written as a sum over all
particle types:
\begin{eqnarray}
  \aAbsVec{p}\DirFre = \sum_i \PDen_i \aSAbVec{p}_i \DirFre
\end{eqnarray}
Here $\aAbsVec{p}_i$ is the particle absorption cross section vector
of the {\sl i}th particle type.  The gaseous absorption vector can be
written as follows:

\begin{eqnarray}
  \aAbsVec{g}(\Prs, \Frq, \Tmp)  = (\alpha(\Prs, \Frq, \Tmp), 0, 0, 0) 
\end{eqnarray}


The last term in Equation (\ref{eq:scattering:RTE}) is the scattering
source term. It is the amount of radiation which is scattered from all
directions \PDir$^\prime$ into the propagation direction \PDir.  The
scattering efficiency matrix \SEfMat\ is the sum of products of the
particle densities \PDen\ and the phase matrices \PhaMat
\begin{eqnarray}
\SEfMat(\PDir, \PDir', \Frq) = \sum_i{\PDen_i\PhaMat_i(\PDir,
  \PDir', \Frq)}.
\end{eqnarray}

Molecular scattering is neglected as it is not important in the
considered frequency range from 10 to 1000 GHz.


\levelb{Scalar radiative transfer equation}
%======================================================
\label{sec:scattering:scalar_rte}

If we know that the radiation field is unpolarized we may use the
scalar radiative transfer equation, which is easier to handle than the
vector equation. It can directly be derived from the general equation
provided that there are no non-spherical particles in the scattering
medium. If there are non-spherical particles the polarization state of
the radiation will be changed, in this case the medium is called
dichroistic.  For spherical particles the extinction matrix has only
elements on the diagonal and all of them have equal values. Only the
first component of the absorption vector will be nonzero.

Instead of taking the full phase matrix, commonly only its first
element, which is usually called the phase function, is used, for
aspherical particles as well as for spherical particles. This concept
has no physical justification as there are non-diagonal terms in the
phase matrix.

Thus to write down the scalar radiative transfer equation the
extinction matrix in (\ref{eq:scattering:RTE}) is replaced by the
first element of the extinction matrix $K_{11}$, the scalar absorption
coefficient is replaced by the first element of the absorption vector
$a_1$ and the scalar scattering coefficient by the first element of
the scattering efficiency matrix $Y_{11}$:
\begin{eqnarray}
  \label{eq:scattering:scalar_rte}
\frac{\DiffD I}{\DiffD s} (\PDir, \Frq) = -K_{11}(\PDir, \Frq)I
(\PDir, \Frq) + a_1(\PDir, \Frq)\Planck(\Frq) + \\ \nonumber
 +\int_{4\pi} \DiffD \PDir^\prime Y_{11}\DirFrePr I(\PDir^\prime,
  \Frq)
\end{eqnarray}

\levelb{Iteration scheme for solving the vector RTE}
%===============================
\label{sec:scattering:solution_rte}

There exists no analytical solution for differential equations as the
general radiative transfer equation. A numerical solution can be
obtained using for example an iterative method. The method implemented
in ARTS is described in this section.

\levelc{Basic definitions}
\begin{itemize}
\item{\bf Cloudbox}
  
  It is not necessary to solve the VRTE in the whole atmosphere. In
  ARTS it is possible to define a region called cloudbox, in which
  cloud particles may exist. Only in this particular region the VRTE
  is solved, this saves computation time and memory.  The cloudbox is
  defined by its corner points which have to correspond to atmospheric
  grid points (cp. Figure \ref{fig:scattering:cloudbox}). For a 3D
  calulation the user has to define lower and upper pressure limit
  (\aPrs{1} and \aPrs{2}), two latitude limits (\aLat{1} and \aLat{2})
  and two longitude limits (\aLon{1} and \aLon{2}), the cloudbox
  concept is described mor detailed in section
  \ref{sec:fm_defs:cloudbox}.

\begin{figure}[htbp]
 \begin{center}
  \begin{minipage}[c]{0.65\textwidth}
   \begin{center}
     \includegraphics*[width=0.9\hsize]{Figs/scattering/cloudbox}
   \end{center}
  \end{minipage}%
  \begin{minipage}[c]{0.35\textwidth}
   \caption{2D cloudbox defined by lower and upper pressure limit (\aPrs{1} and \aPrs{2}) and two latitude limits (\aLat{1} and \aLat{2}).}
   \label{fig:scattering:cloudbox}
  \end{minipage}
 \end{center}
\end{figure}   

The aim is to calculate the radiation at all grid points inside the
cloudbox for all viewing angles, this quantity we call the radiation
field.

\item{\bf Radiation field}
  
  To describe the radiation field inside the cloudbox completely we
  need the Stokes vectors at all grid points. The Stokes vector
  depends not only on the position in the cloudbox but also on the
  propagation direction of the radiation (cp. equation
  \ref{eq:scattering:RTE}). Thus the radiation field is defined to be
  a set including the Stokes vectors at all grid points in the
  cloudbox, which are elements of the atmospheric grids, and for all
  propagation directions defined by special angular grids for the
  scattering calculation:

\begin{eqnarray}
{\mathcal I} = \left\{ \StoVec \left( p_i, \alpha_i, \beta_i, \theta_i, \phi_i\right)\right\}   \quad
\forall \quad p_i &\in& p\_grid,\\
 \alpha_i &\in& lat\_grid, \\
 \beta_i &\in& lon\_grid, \\
 \theta_i &\in& scat\_za\_grid, \\
 \phi_i &\in& scat\_aa\_grid 
\end{eqnarray}

Of course the Stokes vector also depends on the frequency of the
radiation. But in order to save computation memory, the scattering
calculation is only done for one frequency at a time.

\end{itemize}


\levelc{First guess field}
%===========================================================================

The starting point of the iteration method is the first guess field
\aIFld{0}. This is partly determined by the boundary condition given
by the radiation coming from the clear sky part of the atmosphere
traveling into the cloudbox. To take into account this condition, the
first guess field on the cloudbox boundary has to be set to the clear
sky field on the cloudbox boundary.

Inside the cloudbox any field can be chosen as a first guess. In order
to minimize the number of iterations it should be close to the
iterative solution. Our test calculations have shown, that taking the
clearsky field as first guess is reasonable.


\levelc{First scattering integral field}
%============================================================================

The next step is solving the scattering integral
\begin{eqnarray}
  \ScaInt^{(0)} = \int_{4\pi} \DiffD \PDir^\prime
  \SEfMat \StoVec^{(0)}
\label{eq:scat_int}
\end{eqnarray}
using the first guess initial field. For the integration we use again
the angular grid defined for scattering calculations.

The integration has to be performed over all incident directions
$\PDir^\prime$ for each propagation direction \PDir{}. The evaluation
of the scattering integral is done for all grid points inside the
cloudbox. The result is the first guess scattered field \aSFld{0}
which is a set of scatterd field vectors at all positions in the
cloudbox for all directions:

\begin{eqnarray}
{\mathcal S}^{(0)} = \left\{ \bf{S}^{(0)} \left( p_i, \alpha_i, \beta_i, \theta_i, \phi_i\right)\right\}  \qquad  
&\forall& {\rm positions } (p_i,\alpha_i, \beta_i)\\
&\forall& {\rm directions } (\theta_i,\phi_i)
\end{eqnarray}


\levelc{Averaging of the coefficients of the VRTE}
%============================================================================
\label{scattering:av_coeff}

In the following we are looking at one grid point inside the cloudbox
into one specified direction.  We can substitute $\ScaInt^{(0)}$, the
scattering integral vector calculated in the previous step at this
point for the specified direction, into Equation
(\ref{eq:scattering:RTE}), in order to get the VRTE for this point and
direction:
\begin{eqnarray}
     \frac{\DiffD \StoVec^{(1)}}{\DiffD s} =
     -\ExtMat \StoVec^{(1)} + \AbsVec \Planck
     +\ScaInt^{(0)}
\label{eq:scattering:vrte_fs}
\end{eqnarray} 
Here $\StoVec^{(1)}$ is the Stokes vector at the considered grid point
for the specified direction.

\begin{figure}[htbp]
 \begin{center}
  \begin{minipage}[c]{0.65\textwidth}
   \begin{center}
     \includegraphics*[width=0.9\hsize]{Figs/scattering/average}
   \end{center}
  \end{minipage}%
  \begin{minipage}[c]{0.35\textwidth}
   \caption{Path from one grid point (black) to the intersection point with the next grid cell boundary. Multilinear interpolation is done from the grid points on the intersection point. The coefficients are averaged between the black grid point and the intersection point.}
   \label{fig:scattering:averaging}
  \end{minipage}
 \end{center}
\end{figure}   

As mentioned before, all optical properties are assumed to vary
linearly between the grid points. Equation
(\ref{eq:scattering:vrte_fs}) is formally a linear differential
equation, which could be solved analytically if the coefficients were
constant.  To find a constant approximation for the coefficients, we
follow a propagation path from the considered point into the specified
direction until the path intersects a grid cell boundary (cp. Figure
\ref{fig:scattering:averaging}). By multi-linear interpolation we
obtain all coefficients and the temperature at the intersection point.
Also the scattering integral vector is obtained for the intersection
point by multi-linear interpolation of the first guess scattered
field. Now we approximate the coefficients (extinction matrix and
absorption vector) in Equation (\ref{eq:scattering:vrte_fs}) by taking
the average of the coefficients at the intersection point and the
coefficients at the considered point element-vise for each
matrix-/vector-element. The same averaging is done for the scattering
integral vector and the temperature. The average value of the
temperature is used to get the averaged Planck function.

Then Equation (\ref{eq:scattering:vrte_fs}) can be written as follows:
 \begin{eqnarray}
     \frac{\DiffD \StoVec^{(1)}}{\DiffD s} =
     -\bar{\ExtMat} \StoVec^{(1)} + \bar{\AbsVec} \bar{\Planck}
     +\bar{\ScaInt^{(0)}}
\label{eq:scattering:vrte_fs_av}
\end{eqnarray} 
Here $\bar{\ExtMat}$, $\bar{\AbsVec}$, $\bar{\Planck}$ and
$\bar{\ScaInt^{(0)}}$ are the averaged quantities.


\levelc{Radiative transfer with fixed scattered field and averaged
  coefficients}
%===================================================================

Equation (\ref{eq:scattering:vrte_fs_av}) can be solved analytically
using the following matrix exponential approach:
 \begin{eqnarray}
 \label{eq:scattering:ansatz}
&   \StoVec^{(1)} = e^{-\bar{\ExtMat}s}{\bf C_1} + {\bf C_2}
\end{eqnarray}
Here ${\bf C_1}$ and ${\bf C_2}$ are constants which have to be
determined. Substituting (\ref{eq:scattering:ansatz}) into
(\ref{eq:scattering:vrte_fs_av}) gives the constant ${\bf C_2}$:
\begin{eqnarray}
 \Rightarrow& -\bar{\ExtMat}e^{-\bar{\ExtMat}s}{\bf C_1} = -\bar{\ExtMat}e^{-\bar{\ExtMat}s}{\bf C_1} - \bar{\ExtMat} {\bf C_2}+ \bar{\AbsVec} \bar{\Planck} +\bar{\ScaInt}^{(0)}\\
\Rightarrow&  \bar{\ExtMat} {\bf C_2} = \bar{\AbsVec} \bar{\Planck} +\bar{\ScaInt}^{(0)}\\
\Rightarrow& {\bf C_2} =  \bar{\ExtMat}\Inv\left(\bar{\AbsVec} \bar{\Planck} +\bar{\ScaInt}^{(0)}\right)
 \label{eq:scattering:c1}
\end{eqnarray}
${\bf C_1}$ can be determined using the initial condition, which is
the radiation at the intersection point ($s=0$) travelling towards the
considered gridpoint:
\begin{eqnarray}
   \StoVec^{(1)}(s=0) &=& \StoVec^{(0)} {\rm(at\ intersection\ point)}\\
\rm{From\ Ansatz:} &&\\
\StoVec^{(0)} &=& {\bf C_1} +  \bar{\ExtMat}\Inv\left(\bar{\AbsVec} \bar{\Planck} +\bar{\ScaInt^{(0)}}\right)\\
\Rightarrow   {\bf C_1} &=& \bar{\StoVec^{(0)}} -  \bar{\ExtMat}\Inv \left(\bar{\AbsVec} \bar{\Planck} +\bar{\ScaInt^{(0)}}\right)
\label{eq:scattering:c2}
\end{eqnarray}
Substituting (\ref{eq:scattering:c1}) and (\ref{eq:scattering:c2})
into Equation (\ref{eq:scattering:ansatz}) leads to the solution:
\begin{eqnarray}
  \StoVec^{(1)} =  e^{-\bar{\ExtMat}s}\cdot\left(\StoVec^{(0)} -\bar{\ExtMat}\Inv\left(\bar{\AbsVec} \bar{\Planck} +\bar{\ScaInt}^{(0)}\right)\right) + 
 \bar{\ExtMat}\Inv\left(\bar{\AbsVec} \bar{\Planck} +\bar{\ScaInt}^{(0)}\right)
\end{eqnarray}
This can be resorted to the following form:
\begin{eqnarray}
\label{eq:scattering:vrte_sol}
 \StoVec^{(1)} =  e^{-\bar{\ExtMat}s} \StoVec^{(0)} + \left( {\bf E} -  e^{-\bar{\ExtMat}s}\right)  \bar{\ExtMat}\Inv\left(\bar{\AbsVec} \bar{\Planck} +\bar{\ScaInt}^{(0)}\right)
\end{eqnarray}

Here \IdnMat\ denotes the identity matrix.  $\StoVec^{(0)}$ is given
by the initial condition. It is the Stokes vector at the intersection
point which is obtained by interpolation of the first guess field.

This radiative transfer step calculation is done for all points inside
the cloudbox in all directions. The resulting Stokes Vectors (
$\StoVec^{(1)}$ for all points in all directions) form the first
iteration field \aIFld{1}:

\begin{eqnarray}
{\mathcal I}^{(1)} = \left\{ \bf{I}^{(1)} \left( p_i, \alpha_i, \beta_i, \theta_i, \phi_i\right)\right\}  \qquad  
&\forall& {\rm positions } (p_i,\alpha_i, \beta_i)\\
&\forall& {\rm directions } (\theta_i,\phi_i)
\end{eqnarray}


\levelc{Iterations}

The first iteration field is used to evaluate the scattered field
\ScaInt\ again at all grid points inside the cloudbox:
\begin{eqnarray}
  \ScaInt^{(1)} = \int_{4\pi} \DiffD \PDir^\prime
  \SEfMat \StoVec^{(1)} 
\end{eqnarray}
Now $\ScaInt^{(1)}$ is used as fixed scattering term in the radiative
transfer equation for the second iteration to get the Stokes vectors
$\StoVec^{(2)}$
\begin{eqnarray}
     \frac{\DiffD \StoVec^{(2)}}{\DiffD s} =
     -\bar{\ExtMat} \StoVec^{(2)} + \bar{\AbsVec} \bar{\Planck}
     +\bar{\ScaInt}^{(1)}.
\end{eqnarray} 
This equation contains already the averaged values and is valid for
specified position and direction.
 
The solution is given by:
\begin{eqnarray}
   \StoVec^{(2)} = e^{-\bar{\ExtMat} s}\cdot\StoVec^{(1)} + (\IdnMat - e^{-\bar{\ExtMat}
    s}) \bar{\ExtMat}\Inv (\bar{\AbsVec} \bar{\Planck} + \bar{\ScaInt}^{(1)})
\end{eqnarray}
Thus the Stokes vectors ( $\StoVec^{(2)}$ for all points in all
directions) form the second iteration field.

After that the scattering integral and higher order iteration fields
are calculated alternately. We can formulate a differential equation
for the $n$-th order iteration field. The fixed scattering integral
term in this equation is
 \begin{eqnarray}
  \ScaInt^{(n-1)} = \int_{4\pi} \DiffD \PDir^\prime
  \SEfMat \StoVec^{(n-1)}
\end{eqnarray}
and the differential equation for specified position and direction is
given by
\begin{eqnarray}
     \frac{\DiffD \StoVec^{(n)}}{\DiffD s} =
     -\bar{\ExtMat} \StoVec^{(n)} + \bar{\AbsVec} \bar{\Planck}
     +\bar{\ScaInt}^{(n-1)}.
\end{eqnarray} 
Thus the $n$-th order iteration is given by:
\begin{eqnarray}
   \StoVec^{(n)} = e^{-\bar{\ExtMat} s}\cdot\StoVec^{(n-1)} + (\IdnMat - e^{-\bar{\ExtMat}
    s}) \bar{\ExtMat}\Inv (\bar{\AbsVec} \bar{\Planck} + \bar{\ScaInt}^{(n-1)})
\end{eqnarray}


\levelc{Convergence test}

After each iteration the convergence is checked. If
\begin{eqnarray}
|\StoVec^{(m)} \DirFre -  \StoVec^{(m-1)} \DirFre| < {\bf \epsilon}
\end{eqnarray}
for all points inside the cloudbox and all directions a solution to
the vector radiative transfer Equation (\ref{eq:scattering:RTE}) has
been found:
\begin{eqnarray}
{\mathcal I}^{(m)} = \left\{ \bf{I}^{(m)} \left( p_i, \alpha_i, \beta_i, \theta_i, \phi_i\right)\right\}  \qquad  
&\forall& {\rm positions } (p_i,\alpha_i, \beta_i)\\
&\forall& {\rm directions } (\theta_i,\phi_i)
\end{eqnarray}



\levelb{Iteration scheme for solving the scalar RTE}
%===============================
\label{sec:scattering:solution_rte_scalar}
\begin{itemize}
\item {\bf Scattering integral}
  
  In analogy to the scattering integral vector field the scalar
  scattering integral field is obtained:
\begin{eqnarray}
  S^{(0)}  = \int_{4\pi} \DiffD \PDir^\prime Y_{11} I^{(0)} 
\end{eqnarray}
As first guess field only the first Stokes component of the initial
field is needed.

\item{\bf Radiative transfer with fixed scattered field}
  
  The scalar radiative transfer equation
  (\ref{eq:scattering:scalar_rte}) with fixed scattering integral is
\begin{eqnarray}
  \label{eq:scattering:scalar_rte_scatint}
\frac{\DiffD I^{(1)}}{\DiffD s} = -K_{11} I^{(1)}
 + a_1 \Planck + S^{(0)}.
\end{eqnarray} 
Assuming constant coefficients this equation is solved analytically
after averaging extinction coefficient, absorption coefficient,
scattered vector and the temperature. The averaging procedure is done
analogously to the procedure described for solving the VRTE (see
Section \ref{sec:scattering:solution_rte}).  The solution of the
averaged differential equation is
\begin{eqnarray}
   \label{eq:scattering:scalar_rte_sol}
I^{(1)} = I^{(0)} e^{-\bar{K_{11}}s} + \frac{\bar{a_1}
  \bar{B} + \bar{S^{(0)}}}{\bar{K_{11}}}\left(1-e^{-\bar{K_{11}}s}\right)
\end{eqnarray}
where $I^{(0)}$ again obtained by interpolating the initial field.
$\bar{K_{11}}$, $\bar{a_1}$, $\bar{B}$ and $\bar{S^{(0)}}$ are the
averaged values for extinction coefficient, absorption coefficient,
Planck function and the scattered integral respectively.

Applying this equation leads to first iteration scalar intensity
field, consisting of the intensities $I^{(1)}$ at all points in the
cloudbox for all directions.

\item{\bf Iteration and Convergence}
  
  As the solution to the vector radiative transfer equation the
  solution to the scalar radiative transfer equation is found
  numerically by iterations.
  
  The convergence test for the scalar equation compares the values of
  the calculated intensities of two successive radiation fields.

\end{itemize}


\levelb{Scattering databases}
%============================================================
\label{sec:scattering:database}

Two databases which are required for scattering calculations inside
the cloudbox. One contains the optical properties of different
particle types and the other contains the corresponding particle
number density fields.

\levelc{Database for optical properties}
%============================================================
\label{sec:scattering:amp_mat_data}

The required quantities for the radiative transfer calculations are
the extinction cross section matrix, the absorption coefficient vector
and the phase matrix. All these quantities can be calculated for
spherical, cylindrical and spheroidal particles using the T-matrix
method. Applying Mishchenko's FORTRAN code, which is explained in
\citet{Mishchenko:98} and \citet{Mishchenko:00}, the amplitude matrix
can be calculated. The amplitude matrix ${\bf S}$ is a 2x2 complex
matrix, which linearly transforms the electric field vector components
of the incident wave into the electric field vector components of the
scattered wave:
\begin{eqnarray}
  \label{eq:ampl_matrix}
  \left[{E^{\rm sca}_{\theta}\atop E^{\rm sca}_{\phi}}\right] =
  \frac{\exp({\rm i} kR)}{R}{\bf S}({\bf n};{\bf
      n'};\epsilon_1,\epsilon_2,\epsilon_3)\left[{E^{\rm inc}_{\theta
          0}\atop E^{\rm inc}_{\phi 0}}\right] 
\end{eqnarray}
It depends on the directions of incidence $(\theta, \phi)$ and
scattering $(\theta', \phi')$ as well as on size, morphology and
composition of the scattering material, furthermore on the orientation
of the particles which is specified by the Euler angles of rotation
$\epsilon _i$.  The equation above is valid for a monochromatic wave,
generally the amplitude matrix is also frequency dependent.

All information about the optical properties of the scattering medium
is contained in the amplitude matrix. Extinction cross section matrix,
absorption cross section vector and phase matrix can be obtained by
simple additions or multiplications of amplitude matrix components
(see Sections \ref{sec:scattering:ext_mat_spt},
\ref{sec:scattering:pha_mat_spt} and
\ref{sec:scattering:abs_vec_spt}).  For this reason, it is sufficient
to store only the amplitude matrices in the scattering data base
instead of storing all coefficients.

For each particle type determined by the particle size, the particle
shape and the orientation there exists the corresponding file in the
database containing the amplitude matrix for this type.  So far we
have only considered ice particles as we are mainly
interested in cirrus clouds.\\
The methods to read from the database allow to construct homogeneous
clouds as well as inhomogeneous clouds, which may consist of various
particle types.  The data files are stored in XML-format. The datatype
of each file is a 6D gridded field (cp. Section FIXME
%\ref{sec:gridded_fields}
).  The following example shows the datafile for spherical particles
having 200 $\mu$m equal volume sphere radius:

\begin{verbatim}
<?xml version="1.0"?>
<arts format="ascii" version="1">
<comment>
  ------------------------------------------
  Describtion of the particle:
  =========================================
  ice particle                                
  radius : 200 micrometer 
  shape: spherical          
  aspect ratio: 1
 ------------------------------------------- 
</comment>
<Array type="Tensor6" nelem="7">
<Tensor6 nvitrines="1" nshelves="1" nbooks="1" npages="1"
                                nrows="1" ncols="1">
  3.25E+11    // frequency grid, this file contains only one
</Tensor6>
<Tensor6 nvitrines="1" nshelves="19" nbooks="1" npages="1" 
                                 nrows="1" ncols="1">
  0.          
  10.
  20.         // incoming zenith angle grid
   .
   .
   .
 </Tensor6>
 <Tensor6 nvitrines="1" nshelves="1" nbooks="37" npages="1"
                                 nrows="1" ncols="1">
   0.
  10.
   .          // incoming azimuth angle grid   
   .
   .
 </Tensor6>
 <Tensor6 nvitrines="1" nshelves="1" nbooks="1" npages="19"
                                 nrows="1" ncols="1">
  0.
  10.
   .          // scattered zenith angle grid
   .
   .
 </Tensor6>
 <Tensor6 nvitrines="1" nshelves="1" nbooks="1" npages="1"
                                 nrows="37" ncols="1">
  0.
  10.
   .           // scattered azimuth angle grid
   .
   .
  </Tensor6>
 <Tensor6 nvitrines="1" nshelves="1" nbooks="1" npages="1"
                                 nrows="1" ncols="8"> 
  1.
  1.           // dummy grid for the amplitue matrix 
  .            // components
  .
  .
  </Tensor6>
 <Tensor6 nvitrines="1" nshelves="19" nbooks="37" 
                        npages="19" nrows="37" ncols="8">
   0.204E+03   0.102E+03  -0.382E-37   0.186E-37   0.000E+00  
                    0.000E+00   0.204E+03   0.102E+03
   0.201E+03   0.100E+03  -0.354E+02  -0.177E+02   0.354E+02 
                    0.177E+02   0.201E+03   0.100E+03
   .
   .          // the amplitude matrix data 
   .    
   .          // the order of the elements is:
   .          // S11_real, S11_imag, S12_real, S12_imag, 
   .          // S21_real, S21_imag, S22_real, S22_imag    
   .
   0.201E+03   0.100E+03  -0.348E+02  -0.174E+02   0.348E+02  
                    0.174E+02   0.201E+03   0.100E+03
   0.204E+03   0.102E+03  -0.234E-22  -0.446E-22   0.234E-22   
                    0.446E-22   0.204E+03   0.102E+03
</Tensor6>
</Array>
</arts>
\end{verbatim}

This example file shows that each file in the database contains for
each frequency, for each propagation direction (\ScaZa, \ScaAa) and
each scattered direction (\ScaZa$^\prime$, \ScaAa$^\prime$) the real
components of the amplitude matrix ($S^{real}_{11}$, $S^{real}_{12}$,
$S^{real}_{21}$, $S^{real}_{22}$) and the imaginary components
($S^{imag}_{11}$, $S^{imag}_{12}$,$S^{imag}_{21}$, $S^{imag}_{22}$).

\levelc{Database for particle number density fields}
%==========================================================
\label{sec:scattering:pnd_data}

A second database contains the particle number density fields.  The
datatype contained in each data file is a 3D gridded field. It
contains three elements corresponding to pressure, latitude and
longitude grid. The forth element contains the data for the particle
number density. Each chosen particle type has to be defined together
with a particle number density field.

\levelc{Reading the databases}
%===========================================================
\label{sec:scattering:read_data}

Two workspace methods are implemented for reading the database. These
functions are used to define the clouds, the horizontal and vertical
extension and the particle distribution inside the cloud. They are
called \wsmindex{ParticleTypeInit} and \wsmindex{ParticleTypeAdd}.

The workspace variable containing the amplitude matrix data is
\wsvindex{amp\_mat\_raw} which is an array of gridded fields. The
array contains for each particle type one gridded field having the
same dimensions as the files in the database. Therefore the first six
elements of a gridded field contain the frequency grid, the angular
grids and a dummy grid corresponding to the number of elements of the
amplitude matrix (i.e. 8) and the seventh element contains the data
itself.

The analogous workspace variable for the particle number density field
is called \wsvindex{pnd\_field\_raw} which as well is an array of
gridded fields containing one element for each particle type. The
first three elements of each gridded field contain the pressure, the
latitude and the longitude grids and the fourth element the particle
number density data.

The workspace method \wsmindex{ParticleTypeInit} is used for
initializing \wsvindex{amp\_mat\_raw} and \wsvindex{pnd\_field\_raw}.
The method \wsmindex{ParticleTypeAdd} puts the data for one particle
type into \wsvindex{amp\_mat\_raw} and \wsvindex{pnd\_field\_raw}. The
filenames are passed into the methods by keywords which are specified
in the control file.

The following example shows a part of a control file where
\wsvindex{amp\_mat\_raw} and \wsvindex{pnd\_field\_raw} are created
for two particle types:

\begin{verbatim}

# Initialize variables
# --------------------------------------------
ParticleTypeInit{}

# Add spherical particles, size 200 mirons
# --------------------------------------------
ParticleTypeAdd{
        filename_amp_mat = "sph_200_ampmat.xml" 
        filename_pnd_field = "sph_200_pnd.xml"
        }       

# Add cylindrical particles, size 100 mirons
#---------------------------------------------
ParticleTypeAdd{
        filename_amp_mat = "cyl_100_ampmat.xml" 
        filename_pnd_field = "cyl_100_pnd.xml"
        }   

\end{verbatim}
 
There is no convention for the filenames. All needed information about
the particles is included in the data files.


\levelb{Radiative transfer implementation in the cloudbox}
%=================================================================
\label{sec:scattering:rt_cloudbox}


\levelc{Scattering main function}
%=====================================
\label{sec:scattering:main_function}
For a scattering calculation only the workspace method
\wsmindex{ScatteringMain} has to be executed in the control file.  But
before executing the main fuction, agendas required for a scattering
calculation have to be defined.  The following example shows a part of
a controlfile where a scattering calculation is performed:
\begin{verbatim}

#===========================================================
# Agendas for calculating the single scattering properties
# from the amplitude matrix data:
#===========================================================

  AgendaSet(spt_calc_agenda){
      pha_mat_sptCalc{}
      ext_mat_sptCalc{}
      abs_vec_sptCalc{}
  }

  AgendaSet( opt_prop_part_agenda ) {
      ext_matAddPart{}
      abs_vecAddPart{}
  }
 
#===========================================================
# Define the method for the convergence test:
#===========================================================
  AgendaSet( convergence_test_agenda) {
        convergence_flagAbs{
        epsilon = [1e-17, 1e-18, 1e-18, 1e-18] 
        # epsilon = [1e-16]
   }
  
#==========================================================
# Methods for monochromatic scattering calculation:
#==========================================================
  AgendaSet(scat_mono_agenda){
  amp_matCalc{}
  i_fieldSetClearsky{}
  WriteXML(i_field){"initial_field.xml"}        
  i_fieldIterate{}
  scat_iPut{}
 }
                
#==========================================================
# Execute the scattering main function:
#==========================================================       
        
  ScatteringMain{}

\end{verbatim}


All listed agendas are described in detail in the next sections.\\
\vspace*{1ex}

The first step in the scattering main function is the calculation of
the clearsky field on the cloudbox boundary using the method
\wsmindex{CloudboxGetIncoming}.

Then it executes \wsaindex{scat\_mono\_agenda}(see Section
\ref{sec:scattering:scat_mono_ag}) inside a loop over all frequencies
defined in \wsvindex{f\_grid}.

\levelc{Numerical grids and fields}
%====================================================================
\label{sec:scattering:grids}
The calculations inside the cloudbox are performed on the common
atmospheric grids: the pressure grid (\wsvindex{p\_grid}), the
latitude grid (\wsvindex{lat\_grid}) and the longitude grid
(\wsvindex{lon\_grid}).  For the calculation gaseous absorption (see
Section \ref{sec:absorption}) and atmospheric fields
(\wsvindex{t\_field} and \wsvindex{vmr\_field}) are required. More
atmospheric grids and fields (cp. section \ref{sec:fm_defs:grids}) are
needed for the ray tracing.  Furthermore anglular grids are necessary
for the scattering calculations as the radiation field, the scattering
efficiency matrix etc. depend on the propagation and incident
direction of the radiation. For this purpose the workspace variables
\wsvindex{scat\_za\_grid} and \wsvindex{scat\_aa\_grid} have to be
defined in the control file. They hold a zenith angle grid and a
azimuthal angle grid. These angles are defined in the same manner as
the angles which describe the viewing direction of the sensor (cp.
section \ref{sec:fm_defs:los}).  The user has to find a compromise.
The grids should not be too coarse for accuracy reasons but they also
can not be very fine as there is only limited random acces memory
available.  The grids do not neccesarily have to be equidistant. It is
reasonble to take a finer resolution around zenith angle equal
90\degree\ because here the radiation field shows a high variability.
The zenith angles up to an angle a bit higher than 90\degree\ 
correspond to uplooking angles and the angles above correspond to
downlooking angles. Of course the intensity of the radiation coming
from below is much higher than the intensity of the radiation coming
from above as in the lower atmosphere the density is much higher,
which corresponds to higher thermal emission.\newline

\vspace{1ex} The radiation field in the cloudbox is stored in the
workspace variable \wsvindex{i\_field} which has the dimension:
\begin{center}
  \wsvindex{i\_field} = \wsvindex{i\_field} (\Prs, \Lat, \Lon, \ScaZa,
  \ScaAa, i$_I$).
\end{center}
It has to be noted that \wsvindex{i\_field} is defined only in the
cloudbox which means that its size is
\[ \begin{array}{rl}
 N(\artsstyle{i\_field}) =& [\artsstyle{cloudbox\_limits}[1] - \artsstyle{cloudbox\_limits}[0] + 1,\\
              &\artsstyle{cloudbox\_limits}[3] -\artsstyle{cloudbox\_limits}[2] + 1, \\
              &\artsstyle{cloudbox\_limits}[5] -\artsstyle{cloudbox\_limits}[4] + 1,  \\
              & N_{za}, N_{aa}, N_{I} ] 
\end{array}\]
The first three numbers give the sizes of pressure grid, latitude grid
and longitude grid respectively inside the cloud box. $N_{za}$ and
$N_{aa}$ are the size of the zenith and azimuth angle grids and
$N_{I}$ is the Stokes dimension.


\levelc{Interface between cloudbox and clearsky}
%================================================================
\label{sec:scattering:interface}
The interface between the clearsky part and the cloudbox part of the
radiative transfer calulation is the cloudbox boundary. The scattering
calculation requires the incloming clearsky radiation field on the
boundary as initial condition.For the clearsky calculation the
outcoming radiation field on the boundary, which is propagated from
there to the sensor, is needed.  The interface variables are
\wsvindex{scat\_i\_p}, \wsvindex{scat\_i\_lat} and
\wsvindex{scat\_i\_lon}. They have the dimensions
\begin{center}
  \artsstyle{scat\_i\_p} = \artsstyle{scat\_i\_p} (\Frq, 2(two surfaces), \Lat, \Lon, \ScaZa, \ScaAa, i$_I$)\\
  \artsstyle{scat\_i\_lat} = \artsstyle{scat\_i\_lat} (\Frq, \Prs,
  2(two surfaces), \Lon, \ScaZa,
  \ScaAa, i$_I$ )\\
  \artsstyle{scat\_i\_lon} = \artsstyle{scat\_i\_lon} (\Frq, \Prs,
  \Lat, 2(two surfaces), \ScaZa, \ScaAa, i$_I$).
\end{center}
where \Frq\ is the frequency, \Prs\ the pressure, \Lat\ the latitude,
\Lon\ the longitude, \ScaZa\ and \ScaAa\ the zenith of the propagation
the azimuthal angles of the propagation direction respectively and
i$_I$ is the Stokes component.  In 3D geometry the variable
\wsvindex{scat\_i\_p} for example has the size:
\begin{center}
  N(\artsstyle{scat\_i\_p}) = $[N_\Frq, 2, N_\Lat, N_\Lon, N_\ScaZa,
  N_\ScaAa, N_I]$
\end{center} 

Here $N_\Lat$ and $N_\Lon$ are the number of latitude and longitude
gridponts inside the cloudbox region. $N_\ScaZa$ and $N_\ScaAa$ are
the number of grid points in the angular grids for scattering
calculations, $N_I$ is the number of Stokes components, and $N_\Frq$
the number of points in the frequency grid.  There are no special
workspace variables for 1D calculations. The variables which are not
needed, e.g. \wsvindex{scat\_i\_lat} and \wsvindex{scat\_i\_lon}, are
still in the workspace but they are empty. For the interface only
\wsvindex{scat\_i\_p} is needed and its size is in the 1D case
\begin{center}
  N( \artsstyle{scat\_i\_p} ) = $[N_\Frq, 2, 1, 1, N_\ScaZa, 1, N_I]$
\end{center}  
Inside the scattering box we need to specify only
\wsvindex{scat\_za\_grid}.



\leveld{Clearsky field on cloudbox boundary} The method
\wsmindex{CloudboxGetIncoming} calculates for each grid point on the
cloudbox boundary in each propagation direction the radiation field
using the method \wsmindex{RteCalc}, which performs a clearsky
radiative transfer calculation.  Output of this method are the
interface variables listed above.

\leveld{Scattered field on cloudbox boundary} The scattered field on
the cloudbox boundary for given position and direction is obtained by
using the method \wsmindex{CloudboxGetOutgoing}. The method performs
an interpolation from the interface variables on the requested
position and direction.


\levelc{Solution of the monochromatic VRTE}
%========================================================================
\label{sec:scattering:scat_mono_ag}

The method \wsmindex{ScatteringMain} executes
\wsaindex{scat\_mono\_agenda} for each frequency defined in
\wsvindex{f\_grid}. The frequency index \wsvindex{scat\_f\_index} is
used in the agenda as the frequncy for which the calculations are done
has to be known.  The agenda is a list of methods which solve the
monochromatic VRTE and it must be defined in the control file before
executing \wsmindex{ScatteringMain} for examle in the following way:

\begin{verbatim}
AgendaSet(scat_mono_agenda) {
  amp_matCalc{}
  i_fieldSetConst{}
  i_fieldIterate1D{}
  scat_iPut{}
}
\end{verbatim}

\noindent
The methods used by this agenda are described below.

 
\leveld{First guess field}
%=======================================================================

The workspace method \wsmindex{i\_fieldSetClearsky} uses a linear 3D
interpolation scheme to obtain the radiation field on all grid points
inside the cloud box from the clear sky field on the cloudbox
boundary.  This can be taken as a first guess for the iterative
solution method of the RTE.  The method picks only the monochromatic
radiation field corresponding to the actual frequency out of the
variables \wsvindex{scat\_i\_p}, \wsvindex{scat\_i\_lat} and
\wsvindex{scat\_i\_lon}. Output of the method is the initial field
stored in the workspace variable \wsvindex{i\_field}.  The workspace
method \wsmindex{i\_fieldSetClearsky} is implemented in such a way
that the method should adapt itself for 1D or 3D atmosphere. In this
case the main input is the interface (\artsstyle{scat\_i\_p},
\artsstyle{scat\_i\_lat}, \artsstyle{scat\_i\_lon}) the pressure grid
(\artsstyle{p\_grid}), the frequency grid (\artsstyle{f\_grid}), the
frequency index (\artsstyle{scat\_f\_index}) and of course the
boundary of the cloudbox defined in \artsstyle{cloudbox\_limits}.  An
interpolation of \artsstyle{scat\_i\_p} which is defined only on the
cloudbox boundary to all points inside the cloud box is
done to get the initial field. \\

\vspace*{1ex} Another method to set the first guess field is
\wsmindex{i\_fieldSetConst}. Using this method the user can assign a
constant value for each Stokes component of the first guess field. On
the boundary of course the clearsky field has to be taken, as this is
the boundary condition for the scattering calculation.  Setting the
first guess field to a constant unrealistic value is useful for
testing the solution. A calculation can be done with different first
guess fields. If the radiation field converges each time towards the
same solution field, it can be assumed that the solution field is the
correct one.


\leveld{Iteration}
%=====================================================================
The method \wsmindex{i\_fieldIterate} solves the VRTE using the
iterative method.  The function has included switches to adapt
automatically to the atmospheric dimensionality specified by the
workspace variable \wsvindex{atmosphere\_dim}. Note that only 1D or 3D
scattering calculations are possible.  The following steps are
performed in each iteration until the solution of the radiation field
converges.

\begin{itemize}
\item \artsstyle{i\_field} is copied to \artsstyle{i\_field\_old}.
\item The scattered field is calculated using the method
  \wsmindex{scat\_fieldCalc} (see Sections
  \ref{sec:scattering:sca_fieldCalc} and
  \ref{sec:scattering:solution_rte}).
\item Calculate the new radiation field using\\
  \wsmindex{i\_fieldUpdate1D} or \wsmindex{i\_fieldUpdate3D}.  (see
  Sections \ref{sec:scattering:RT_methods} and
  \ref{sec:scattering:solution_rte}).
\item Do the convergence test (see Sections
  \ref{sec:scattering:conv_method} and
  \ref{sec:scattering:solution_rte}).
\end{itemize}
The solution of the radiative transfer equation is returned as output
using the variable \artsstyle{i\_field}.


\levelc{Computation of the scattered field}
%=========================================
\label{sec:scattering:sca_fieldCalc}
The method \wsmindex{scat\_fieldCalc} does exactly what is described
in Section \ref{sec:scattering:solution_rte}.  Output of this method
is the scattered field (\wsvindex{scat\_field}) which has the same
dimension as the intensity field (\artsstyle{i\_field}).  It is
generated by integrating the product of intensity field
(\artsstyle{i\_field}) and phase matrix \artsstyle{pha\_mat} over all
incident zenith and azimuth angles.  In addition to
\artsstyle{i\_field}, the phase matrix for the single particles
(\artsstyle{pha\_mat\_spt}) and the particle number density field
(\artsstyle{pnd\_field}) are inputs to this method.  They are used to
calculate the phase matrix \artsstyle{pha\_mat} by calling the method
\artsstyle{pha\_mat\_partCalc}\ (see Section
\ref{sec:scattering:pha_mat} ). How the phase matrices for the individual particles are calculated is desribed in detail in Section \ref{sec:scattering:pha_mat_spt} \\
\artsstyle{scat\_fieldCalc} uses the trapezoidal integration method
described in Section \ref{sec:integration}. More sophisticated
integration methods are not implemented yet.
%The integration method used is the trapezoidal integration implemented
%by the function \artsstyle{AngIntegrate\_trapezoid} in the file
%\artsstyle{math\_funcs.cc}.  This explains why we need this
%computation of scattered field to be
%put in as an agenda since there can be different methods by which we
%can perform the integral. 
%The output \artsstyle{scat\_field} is used in the
%methods \artsstyle{i\_fieldUpdate1D} and \artsstyle{i\_fieldIterate}.


\levelc{Solving the VRTE with fixed scattering integral}
%==================================================================
\label{sec:scattering:RT_methods}

The methods \wsmindex{i\_fieldUpdate1D} and
\wsmindex{i\_fieldUpdate3D} are used to evaluate the RT Equation
(\ref{eq:scattering:vrte_fs}) and update the radiation field
\wsvindex{i\_field} for each iteration. The theoretical description of
these functions can be found in Section
\ref{sec:scattering:solution_rte}.

Boths functions consist basically of two parts:

\leveld{Calculate coefficients of VRTE} In the first part, the
coefficients required for the radiative transfer calculation are
determined for all gridpoints inside the cloudbox and on the cloudbox
boundary. These coefficients are the total absorption vector, the
total extrinction matrix and the scattered field vector. They are
obtained in the following way:
\begin{itemize}
\item Scattered field vector: Scattered field
  (\artsstyle{scat\_field}) is input variable.
\item Total extinction matrix total and absorption vector:
  Contributions from gaseous species and from particles have to be
  taken into accout.  The first step is calculating the scatar gas
  absorption using \wsaindex{scalar\_gas\_absorption\_agenda}
  described in Section \ref{sec:absorption}.  Then
  \wsaindex{opt\_prop\_part\_agenda} and
  \wsaindex{opt\_prop\_gas\_agenda} are used to compute total
  extinction and absorption. They can be defined as follows:
  
  \vspace{1ex}
\begin{minipage}{0.9\hsize}
\begin{verbatim}
AgendaSet(opt_prop_gas_agenda)
{
        ext_matInit{} \\ Initialize ext_mat and abs_vec.
        abs_vecInit{} \\ 
      ext_matAddGas{} \\ Add gaseous extinction.              
      abs_vecAddGas{} \\ Add gaseous absorption.
}
\end{verbatim}
\end{minipage} 

\vspace{1ex}
\begin{minipage}{0.9\hsize}
\begin{verbatim}
AgendaSet(opt_prop_part_agenda)
{
      ext_matAddPart{} \\ Add particle extinction.            
      abs_vecAddPart{} \\ Add particle absorption.
}
\end{verbatim}
\end{minipage} 

Note: It is important, that the initialization is done in
\wsaindex{opt\_prop\_gas\_agenda} as this agenda is called before
\wsaindex{opt\_prop\_part\_agenda}.It is also important to put the
initializing functions inside the agendas and not somewhere else in
the controlfile as these agendas are executed inside loops over all
positions in the cloudbox. For each loop index \wsvindex{abs\_vec} and
\wsvindex{ext\_mat} have to be reinitialized, otherwise all extinction
matrices and absorption vectors would be added up.

\item The Planck function at each point can be calculated using
  \wsvindex{t\_field} which is an input variable.

\end{itemize}

All coefficients are stored in arrays during the monochromatic
radiative transfer equation. They aer replaced when the calculation
for the next frequency starts.

\leveld{Radiative transfer step with constant coefficients}

In this part several steps are executed inside a loop over all
positions inside the cloudbox and all directions to obtain the
radiation field:
\begin{itemize}
\item Determine the intersection point with the next layer or grid
  cell boundary using \wsaindex{ppath\_step\_agenda}
\item Interpolate optical properties and the Stokes vector
  (\wsvindex{i\_field}) on the intersection point.
\item Average optical properties to find the coefficients for the VRTE
  with constant coefficients (see Section \ref{scattering:av_coeff})
\item Do a radative transfer calculation with constant coefficients
  using the function \wsmindex{RteCalc}.
\end{itemize}




\levelc{Single scattering properties}
%==================================================================
\label{sec:scattering:single_scat_prop}

\leveld{Method which generates \wsvindex{amp\_mat} from
  \artsstyle{amp\_mat\_raw}}
%==================================================================
\label{sec:scattering:gen_ampmat}

As we know from the previous section \ref{sec:scattering:read_data}
the workspace variable containing the amplitude matrix data is
\artsstyle{amp\_mat\_raw} which is an array of gridded fields. The
amplitude matrix data given by the seventh element of the workspace
variable \artsstyle{amp\_mat\_raw} is calculated for frequencies given
by the first element and angles given by the second, third, fourth and
fifth elements.  The method \wsmindex{amp\_matCalc} takes
\artsstyle{amp\_mat\_raw} as input and interpolates the amplitude
matrix data onto the respective frequency and angular grids necessary
for the calculation.  The output of this method is
\artsstyle{amp\_mat}.  The size of \artsstyle{amp\_mat} is [$N_{pt},
N_{za}, N_{aa}, N_{za}, N_{aa}, 8$]. Here $N_{pt}$ gives the number of
particle types considered for the calculation which is specified from
the input \artsstyle{amp\_mat\_raw}. $N_{za}$ gives the number of
zenith angles considered for scattering calculation which is specified
by the input \wsvindex{scat\_za\_grid}, and $N_{aa}$ gives the number
of azimuth angle considered for scattering calculation which is
specified by the input \wsvindex{scat\_aa\_grid}.  The workspace
variable \artsstyle{amp\_mat} is calculated for one frequency
specified by the input variables \artsstyle{f\_grid} and
\artsstyle{scat\_za\_grid}. This is in compliance with our radiative
transfer scheme which has frequency as its outermost loop.

Chapter \ref{sec:interpolation} has a more detailed documentation of
the different interpolation schemes implemented in ARTS. A couple of
examples are also discussed there.  In this case we interpolate the
seventh element of \artsstyle{amp\_mat\_raw} which is a Tensor6 having
size [$N_{f}, N_{za}, N_{aa}, N_{za}, N_{aa}, 8$] from one gridded
field to another.  This is a green type interpolation of all columns
of this Tensor6.  The interpolation is done simultaneously for the 5
dimensions.  The first step is to set up the grid position arrays by
calling the function \artsstyle{gridpos}.  The output of this function
is the grid positions for the new grid stored in ArrayOfGridPos and
the inputs are the original grid and the new grid where we want to
have the interpolated values. Interpolation weight tensors can be
computed by the function called \artsstyle{interpweights}.  This
computes the interpolation weights simultaneously for all the five
dimensions simultaneously. For this step also we do not need the
actual fields, just the grid positions. Now, the last step is to do a
green interpolation for all columns of the Tensor6.  This is by
calling the function \artsstyle{interps} which takes in as input the
gridpositions, interpolation weights and the source field and output
is the new field.

\leveld{Generate scattering efficiency matrix}
%========================================================================
\label{sec:scattering:tot_pha_mat}

\levele{Calculate the total scattering efficiency matrix:}
%=======================================================================
\label{sec:scattering:pha_mat}

Output of the workspace method \wsmindex{pha\_matCalc} is the phase
matrix stored in the variable \wsvindex{pha\_mat} with dimension
[$N_{\ScaZa}, N_{\ScaAa}, N_{I}, N_{I}$], where $N_{I}$ denotes the
chosen Stokes dimension.  $N_{\ScaZa}$ and $N_{\ScaAa}$ are the number
of angles in the angular grids defined for
the scattering calculations.\\
Input is first the workspace variable \wsvindex{pha\_mat\_spt} (see
Section \ref{sec:scattering:pha_mat_spt}). The second input variable
is the local particle number density vector \artsstyle{pnd\_field} for
all
particle types chosen in the control file.\\
The function sums up the phase matrices for all particle types
weighted with the particle number density.
\begin{eqnarray}
\SEfMat = \sum_i{\PDen\PhaMat}.
\end{eqnarray}

\levele{Generate phase matrix for single particle types:}
%==============================================================
\label{sec:scattering:pha_mat_spt}

From the amplitude matrix the phase matrix for each particle type is
calculated using the workspace method \wsmindex{pha\_mat\_sptCalc}.
The phase matrix depends on the direction of the incoming radiation,
so it is stored for all incoming directions specified by
\wsvindex{scat\_za\_grid} and
\wsvindex{scat\_aa\_grid}.\\
Output of the method is a tensor5 including phase matrices for all
particle types. Therefore it is of dimension [$N_{pt}, N_{\ScaZa},
N_{\ScaAa}, N_{I}, N_{I}$], where $N_{pt}$ is the number of particle
types and $N_{I}$ denotes the Stokes dimension. $N_{\ScaZa}$ and
$N_{\ScaAa}$ are the number of angles in the angular grids defined for
the scattering calculations.\\
As input variables it needs the amplitude matrix (\wsvindex{amp\_mat})
and the propagation direction given by the local zenith and azimuth
angle specified by the indices \wsvindex{scat\_za\_index} and
\artsstyle{scat\_aa\_index}.  The phase matrices \PhaMat{} are
calculated for all particle types chosen in the control file using the
following equations:
\begin{eqnarray}
  Z_{11} &=& \frac{1}{2}(|S_{11}|^2+|S_{12}|^2+|S_{21}|^2+|S_{22}|^2)\\
  Z_{12} &=& \frac{1}{2}(|S_{11}|^2-|S_{12}|^2+|S_{21}|^2-|S_{22}|^2)\\
  Z_{13} &=& -\Re(S_{11}S_{12}^*+S_{22}S_{21}^*)\\
  Z_{14} &=& -\Im(S_{11}S_{12}^*-S_{22}S_{21}^*)\\
  Z_{21} &=& \frac{1}{2}(|S_{11}|^2+|S_{12}|^2-|S_{21}|^2-|S_{22}|^2)\\
  Z_{22} &=& \frac{1}{2}(|S_{11}|^2-|S_{12}|^2-|S_{21}|^2+|S_{22}|^2)\\
  Z_{23} &=& -\Re(S_{11}S_{12}^*-S_{22}S_{21}^*)\\
  Z_{24} &=& -\Im(S_{11}S_{12}^*+S_{22}S_{21}^*)\\
  Z_{31} &=& -\Re(S_{11}S_{21}^*+S_{22}S_{12}^*)\\
  Z_{32} &=& -\Re(S_{11}S_{21}^*-S_{22}S_{12}^*)\\
  Z_{33} &=& \Re(S_{11}S_{22}^*+S_{12}S_{21}^*)\\
  Z_{34} &=& \Im(S_{11}S_{22}^*+S_{21}S_{12}^*)\\
  Z_{41} &=& -\Im(S_{21}S_{11}^*+S_{22}S_{12}^*)\\
  Z_{42} &=& -\Im(S_{21}S_{11}^*-S_{22}S_{12}^*)\\
  Z_{43} &=& \Im(S_{22}S_{11}^*-S_{12}S_{21}^*)\\
  Z_{44} &=& \Re(S_{22}S_{11}^*-S_{12}S_{21}^*)
\end{eqnarray}


\leveld{Agenda for calculating the optical properties}
%==================================================================
\label{sec:scattering:opt_prop_part_agenda}

The agenda \wsaindex{opt\_prop\_part\_agenda} calculates the particle
extinction matrix \ExtMat\ and the particle absorption vector \AbsVec
. Output and input are the workspace variables \wsvindex{ext\_mat} and
\wsvindex{abs\_vec}.  The dimensions of these variables are:
\begin{center}
  \artsstyle{ext\_mat} = \artsstyle{ext\_mat}(\Frq, i$_I$, i$_I$)\\
  \artsstyle{abs\_vec} = \artsstyle{abs\_vec}(\Frq, i$_I$ )
\end{center}
where \Frq is the frequency and i$_I$ is the Stokes component.

As input for calculating the particle extinction the following
variables are required: \wsvindex{pnd\_field}, i.e. the particle
number density, and \wsvindex{ext\_mat\_spt}, which contains the
extinction cross section matrices for each chosen particle type. How
this variable is generated is explained in Section
\ref{sec:scattering:ext_mat_spt}. For calculating the gaseous
extinction the volume mixing ratio (\wsvindex{vmrs}), the pressure
(\wsvindex{p\_grid}), the temperature (\wsvindex{t\_field}) and the
propagation direction given by zenith angle
(\wsvindex{scat\_za\_grid}) and azimuth angle
(\wsvindex{scat\_aa\_grid}) are needed. An example how the agenda can
be defined is given in Section \ref{sec:scattering:RT_methods}.  The
methods \wsmindex{ext\_matAddPart} and \wsvindex{abs\_vecAddPart} are
used in that example. In that example the properties for gases are
calculated first and \wsaindex{opt\_prop\_part\_agenda} adds the
particle properties the the gaseous properties to obtain the total
extinction matrix.

The workspace methods \wsmindex{ext\_matAddPart} and
\wsmindex{abs\_vecAddPart} sum up the extinction matrices and
absorption vector respectively for all particle types weighted with
the particle number density. \begin{eqnarray}
  \aExtMat{p} &=& \sum_i \PDen_i \aSExMat{p}_i\\
  \aAbsVec{p} &=& \sum_i \PDen_i \aSAbVec{p}
\end{eqnarray}

\levele{Generate extinction matrix for single particle types:}
%==============================================================
\label{sec:scattering:ext_mat_spt}

From the database for single scattering properties (see Section
\ref{sec:scattering:database}) the extinction matrices for each
particle type can be calculated using the workspace method
\wsmindex{ext\_mat\_sptCalc}.\\
Output is a tensor including local extinction matrices for each
particle type, so it is of the dimension [$N_{pt}, N_{I}, N_{I}$],
where $N_{I}$ denotes the Stokes dimension.\\
As input variables it needs the amplitude matrix (\wsvindex{amp\_mat})
and the propagation direction specified by the angular grid indices
\wsvindex{scat\_za\_grid} and \wsvindex{scat\_aa\_grid}. Furthermore
the frequency is required. It is passed into the function by the
workspace variables \wsvindex{scat\_f\_index} and \wsvindex{f\_grid}.
\wsvindex{amp\_mat} is obtained using the reading routine
\wsmindex{get\_amp\_fromDb} which is explained in Section
\ref{sec:scattering:database}.


The extinction cross section matrices \SExMat{} are calculated for all
particle types chosen in the control file. The following formulas are
used:
\begin{eqnarray}
  \label{eq:gen_ext_matrix}
  L^p_{jj}({\bf n}) &=& \frac{c}{\nu} \Im[S_{11}({\bf n},{\bf
    n})+S_{22}({\bf n},{\bf n})], \quad j=1,...,4 \\
  L^p_{12}({\bf n}) = L^p_{21}({\bf n}) &=& \frac{c}{\nu} \Im[S_{11}({\bf n},{\bf
    n})-S_{22}({\bf n},{\bf n})]\\
  L^p_{13}({\bf n}) = L^p_{31}({\bf n}) &=& -\frac{c}{\nu} \Im[S_{12}({\bf n},{\bf
    n})+S_{21}({\bf n},{\bf n})]\\
  L^p_{14}({\bf n}) = L^p_{41}({\bf n}) &=& \frac{c}{\nu} \Re[S_{21}({\bf n},{\bf
    n})-S_{12}({\bf n},{\bf n})]\\
  L^p_{23}({\bf n}) = -L^p_{32}({\bf n}) &=& \frac{c}{\nu} \Im[S_{21}({\bf n},{\bf
    n})-S_{12}({\bf n},{\bf n})]\\
  L^p_{24}({\bf n}) = -L^p_{42}({\bf n}) &=& -\frac{c}{\nu} \Re[S_{12}({\bf n},{\bf
    n})+S_{21}({\bf n},{\bf n})]\\
  L^p_{34}({\bf n}) = -L^p_{43}({\bf n}) &=& \frac{c}{\nu} \Re[S_{22}({\bf n},{\bf
    n})-S_{11}({\bf n},{\bf n})]
\end{eqnarray}


\levele{Generate absorption vector for single particle types:}
%==============================================================
\label{sec:scattering:abs_vec_spt}

From the amplitude matrix the absorption vector for each particle type
is calculated using the workspace method
\wsmindex{abs\_vec\_sptCalc}.\\
Output is the absorption vector tensor, which is of dimension
[$N_{pt}, N_{I}$], where $N_{pt}$ is the number of particle
types  and $N_{I}$ denotes the Stokes dimension.\\
As input it uses the extinction matrix (\wsvindex{ext\_mat}) and the
phase matrix
(\wsvindex{pha\_mat}).\\
The absorption cross section vectors \SAbVec{} are calculated for all
particle types chosen in the control file. The following formulas are
used:
\begin{eqnarray}
  \label{eq:gen_abs_vector}
  b^p_1({\bf n}) &=&  L^p_{11}({\bf n}) - \int_{4\pi} d{\bf n'}
  Z_{11}({\bf n, n'})\\
  b^p_2({\bf n}) &=&  L^p_{21}({\bf n}) - \int_{4\pi} d{\bf n'}
  Z_{21}({\bf n, n'})\\
  b^p_3({\bf n}) &=&  L^p_{31}({\bf n}) - \int_{4\pi} d{\bf n'}
  Z_{31}({\bf n, n'})\\
  b^p_4({\bf n}) &=&  L^p_{41}({\bf n}) - \int_{4\pi} d{\bf n'}
  Z_{41}({\bf n, n'}) \\
\end{eqnarray}


\levelc{Radiative transfer step calculation with averaged
  coefficients}
%=================================================================
\label{sec:scattering:rte_step}
The method for doing a radiative step calculation through one grid
cell for a 3D atmosphere or one layer for a 1D atmosphere is called
\wsmindex{RteCalc}. The VRTE with constant coefficients and fixed
scattering integral vector is solved.  \wsmindex{RteCalc}
automatically chooses the most efficient way to compute the radiative
transfer. There are three possibilities:

\begin{itemize}
\item {\bf General method}\\
  If the Stokes dimension is greater than one and the extinction
  matrix non-diagonal this method is chosen.  The two terms in the
  radiative transfer equation \ref{eq:scattering:vrte_fs_av} are
  computed separately using different functions. The first term
\begin{equation}
 e^{-\ExtMat s}\cdot\StoVec_0
\end{equation}
includes the matrix exponential function. This is solved numerically
using the Pad\'e approximation as implemented in the function
\artsstyle{matrix\_exp} (see Section \ref{sec:lin_alg:mat_exp}).  The
second term
\begin{equation}
(\IdnMat - e^{-\ExtMat
    s}) \ExtMat\Inv (\AbsVec \Planck + \ScaInt_0)
\end{equation}
includes the inverse of \ExtMat\ multiplied with a vector. This is
computed numerically by a LU decomposition as described in the
Sections \ref{sec:lin_alg:backsub} and \ref{sec:lin_alg:lu_decomp}.
The functions used here are \artsstyle{ludcmp} and
\artsstyle{lubacksub}. The matrix exponential is
computed using again \artsstyle{matrix\_exp}.\\
Finally the method sums up the two terms to get the updated Stokes
Vector.
 
\item {\bf Scalar case}\\
  If the Stokes dimension is one, only the scalar RTE has to be
  solved.  In this case it is straightforward to compute a radiative
  transfer step through one grid cell/layer according to equation
  (\ref{eq:scattering:scalar_rte_sol}). We only have to compute the
  scalar exponential function and do not need the matrix inverse, such
  that the standard C math library is sufficient.

  
\item{\bf Diagonal exinction matrix and no scattering}
\end{itemize}
{\Large: FIXME: The description of RteCalc should be shifted to
  another place}


\levelc{Convergence test method}
%====================================
\label{sec:scattering:conv_method}

Different methods to do the convergence test are implemented. The user
has to chose one by setting the \wsvindex{convergence\_test\_agenda}.
The following methods can be used:

\begin{itemize}
\item {\bf Method \artsstyle{convergence\_flagAbs}:}\\
  The function calculates the absolute differences for two successive
  iteration fields. It picks out the maximum values for each Stokes
  component separately. The convergence test is fulfilled under the
  following conditions:
\begin{eqnarray}
|I_{m+1} - I_m| < \epsilon_1    \\
|Q_{m+1} - Q_m| < \epsilon_2    \\
|U_{m+1} - U_m| < \epsilon_3    \\
|V_{m+1} - V_m| < \epsilon_4     
\end{eqnarray}

The limits for convergence have to be set in the control file by
setting the vector \artsstyle{epsilon} as a keyword to appropriate
values.  The conditions have to be valid for all positions in the
cloudbox and for all directions.  Then the workspace variable
\wsvindex{convergence\_flag} is set to 1.


\item {\bf Method \artsstyle{convergence\_flagLsq}:}\\
  This method performs a least square test.
\end{itemize}



 


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "uguide"
%%% End: 
% LocalWords:  Emde ext matrix abs vec pha pnd sca lat lon za aa pt FIXME Eq
% LocalWords:  Eqs mishchenko scatt nonsp partic RTE
