%
% 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.\\
\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\\
\Tmp          & \artsstyle{t\_field}         & temperature\\
\PDen         & \artsstyle{pnd\_field} & particle density field \\
\ScaInt       & \artsstyle{sca\_vec} & scattered field,
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\\
\IPart        & \artsstyle{scat\_pt\_index} & particle type index\\
\Planck       &                             & Planck function
\label{symtable:scattering}
\stopsymbols
%=====================================================================


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


\levelb{Polarization and Stokes Parameters}
%=====================================================================
\label{sec:scattering:pol_stokes}

FIXME: some nice introductory sentence (motivation) here.

Plane monochromatic electromagnetic waves are commonly written in the form
\begin{equation} 
  \label{eq:scattering:e_field}
  \VctStl{E}(\VctStl{x},t) = \left[E_v \atop E_h\right] 
  e^{i(\VctStl{k}\VctStl{x} - \omega t)}
\end{equation}
where $\VctStl{E}$ is the electric field vector, the subscripts $v$
and $h$ denote the components with vertical and horizontal
polarization, respectively. $E_v$ and $E_h$, the amplitudes, are
complex numbers, $\VctStl{k}$ and $\omega$ are the wavenumber vector
and the angular frequency, respectively, of the plane wave.  It is
always implicitly understood that the actual, physical, electric field
is the real part of the above expression. Rewriting the complex
amplitudes $E_v$ and $E_h$ using real, positive amplitudes $a_v$ and
$a_h$, and phases $\delta_v$ and $\delta_v$, 
\begin{equation}
  \label{eq:scattering:compl_ampl}
  E_v=a_v e^{i\delta_v}\mbox{ , }
  E_h=a_h e^{i\delta_h}
\end{equation}
the actual electric field vector $\VctStl{\tilde{E}}$ is
\begin{equation}
  \label{eq:scattering:actual_efield}
  \VctStl{\tilde{E}}(\VctStl{x},t) = \Re[\VctStl{E}(\VctStl{x},t)] 
    = \left[a_v\cdot \cos(\VctStl{k}\VctStl{x} - \omega t + \delta_v) 
                       \atop 
            a_h\cdot \cos(\VctStl{k}\VctStl{x} - \omega t + \delta_h) 
       \right] 
\end{equation}

In general, instruments do not measure the electric or magnetic field
vectors of an electromagnetic wave, but rather the time-averaged
intensity, i.e., the energy flux, $F$. This is the time-averaged Poynting
vector (which, in turn, is proportional to the square of the electric
field), thus:
\begin{eqnarray}
  F%(\VctStl{x}) 
  &=& 
  \sqrt{\frac{\epsilon}{\mu}}\overline{(\VctStl{E}(\VctStl{x},t))^2}\\
   \nonumber
  &=&
  \sqrt{\frac{\epsilon}{\mu}}\left(
    a_v^2\overline{\cos^2(\VctStl{k}\VctStl{x} - \omega t + \delta_v)}
    + a_h^2\overline{\cos^2(\VctStl{k}\VctStl{x} - \omega t + \delta_h)}
   \right)
  \label{eq:scattering:intensity}
\end{eqnarray}
The overline denotes the time average
%, which is taken over one oscillation period $T$ ($T=2\pi\omega$). 
%Since the time average of
which for cosine squares is 1/2, thus:
\begin{equation}
  \label{eq:scattering:intensity_}
 F(\VctStl{x}) = 
  \frac{1}{2}\sqrt{\frac{\epsilon}{\mu}}(
    a_v^2 + a_h^2)  
\end{equation}
Taking into account that for plane, monochromatic waves 
the time average always results in a factor
$\frac{1}{2}$, we can also directly write the intensity using the
electric field vector in complex notation
(Eq.~\ref{eq:scattering:e_field}).
\begin{eqnarray}
  \label{eq:scattering:intensity_from_complex}
  F(\VctStl{x}) &=&  \frac{1}{2}\sqrt{\frac{\epsilon}{\mu}}
                      \VctStl{E} (\VctStl{x},t) 
                    \cdot \VctStl{E}^\ast (\VctStl{x},t)\\ \nonumber
               &=&     \frac{1}{2}\sqrt{\frac{\epsilon}{\mu}}
                     (E_v E_v^\ast + E_h E_h^\ast)
\end{eqnarray}
where the asterisk denotes complex conjugation%
\footnote{Note that the same argumentation works for ???} .

Commonly the factor  $\frac{1}{2}\sqrt{\frac{\epsilon}{\mu}}$ is
omitted, and in addition to the flux, three more intensity quantities
are defined as in the following equations. They are called 
\emph{Stokes parameters}:
\begin{eqnarray}
  \label{eq:scattering:stokesparam_I}
  I &=&   E_v E_v^\ast + E_h E_h^\ast \\
  \label{eq:scattering:stokesparam_Q}
  Q &=&   E_v E_v^\ast - E_h E_h^\ast \\
  \label{eq:scattering:stokesparam_U}
  U &=& - E_v E_h^\ast - E_h E_v^\ast \\
  \label{eq:scattering:stokesparam_V}
  V &=& i(E_h E_v^\ast - E_v E_h^\ast)
\end{eqnarray}
Using the amplitude/phase notation from
Eq.~(\ref{eq:scattering:compl_ampl}),
we can rewrite the Stokes parameters as 
\begin{eqnarray}
  \label{eq:scattering:stokesparam_alt_I}
 I &=& a_v^2 + a_h^2\\
  \label{eq:scattering:stokesparam_alt_Q}
 Q &=& a_v^2 - a_h^2\\
  \label{eq:scattering:stokesparam_alt_U}
 U &=&  - a_v a_h \cos(\delta_v-\delta_h)\\
  \label{eq:scattering:stokesparam_alt_V}
 V &=&   a_v a_h \sin(\delta_v-\delta_h)
\end{eqnarray}
The Stokes parameters fully characterize the electromagnetic wave and
therefore contain the same information as the electric field vector
(except for one absolute phase).  Since instruments generally measure
intensities (fluxes), describing electromagnetic radiation by the
Stokes parameters is more practical than describing it by the electric
(or magnetic) field vector. Furthermore, the Stokes parameters are
always real numbers.

 
In order understand what the Stokes parameters mean, we have to go
back to the electric field vector and see what polarization state it
describes.  To do so, we look at the curve that the tip of the
physical electric field vector $\VctStl{\tilde{E}}$ describes with
time at a fixed position $\VctStl{x_0}$:
\begin{eqnarray}
  \tilde{E}_v (t) &=& a_v \cos(\Delta_v - \omega t)\\
  \tilde{E}_h (t) &=& a_h \cos(\Delta_h - \omega t)
\end{eqnarray}
where $\Delta_{v,h} = \VctStl{k}\VctStl{x_0} + \delta_{v,h}$. 
To see that this is an ellipse, we first split the cosines using
the addition theorem:
\begin{eqnarray}
  \label{eq:scattering:tip_of_fieldvec1}
  \tilde{E}_v (t) &=&   a_v \cos\Delta_v \cos(\omega t)
                      + a_v \sin\Delta_v \sin(\omega t)\\
  \label{eq:scattering:tip_of_fieldvec2}
  \tilde{E}_h (t) &=&   a_h \cos\Delta_h \cos(\omega t)
                      + a_h \sin\Delta_h \sin(\omega t)
\end{eqnarray}

In order to have the tip of $\VctStl{\tilde{E}}$ describe an ellipse with semi-major axis $a_0 \cos\beta$
and
semi-minor axis $a_0 \sin\beta$, where $a_0^2 = a_v^2 + a_h^2$, 
it should have the following form
\begin{eqnarray}
  \label{eq:scattering:ellipse_parallel}
 \tilde{E}_v (t) &=&   a_0 \sin\beta \cos(\omega t)\\
 \tilde{E}_h (t) &=&   a_0 \cos\beta \sin(\omega t)
\end{eqnarray}
Here $\beta$ must be between $-45\degree$ and $45\degree$: the tip of
the vector $\VctStl{\tilde{E}}$ describes a circle for $\beta = \pm
45\degree$ (circular polarization), oscillates along the $h$-axis for
$\beta = 0$ (linear polarization) and else describes an ellipse. The
sense of rotation is counterclockwise for positive and clockwise for
negative $\beta$. Since $|\tan\beta|$ is the ratio of the semi-minor
and semi-major axes of the ellipse, $\beta$ is called the ellipticity
angle.

Note that the semi-major axis is oriented along the positive $h$-axis.  To
have the major axis of the ellipse enclose an arbitrary angle $\zeta$
($0 \leq \zeta < 180\degree$) with the $h$-axis, we apply a rotation
matrix and get the equation for an ellipse with arbitrary shape
(ellipticity) and orientation:
\begin{eqnarray}
  \label{eq:scattering:ellipse_rotated1}
 \tilde{E}_v (t) &=&  a_0(  \sin\beta \cos(\omega t) \cos\zeta
                           +\cos\beta \sin(\omega t) \sin\zeta )\\
  \label{eq:scattering:ellipse_rotated2}
 \tilde{E}_h (t) &=&  a_0( -\sin\beta \cos(\omega t) \sin\zeta
                           +\cos\beta \sin(\omega t) \cos\zeta )
\end{eqnarray}

Now we want to establish a direct connection between the parameters
$\beta$ and $\zeta$ describing the shape (ellipticity) and orientation
of the polarization ellipse on the one hand, and the amplitudes $a_v$
and $a_h$ and phases $\delta_v$ and $\delta_v$ of the components of
the electric field vector on the other hand.  Comparing the
$\sin(\omega t)$ and $\cos(\omega t)$ terms in
Eqs.~(\ref{eq:scattering:ellipse_rotated1})%
-(\ref{eq:scattering:ellipse_rotated2}) with the corresponding terms
in Eqs.~(\ref{eq:scattering:tip_of_fieldvec1})%
-(\ref{eq:scattering:tip_of_fieldvec2}), we get:
\begin{eqnarray}
  \label{eq:scattering:corresp1a}
 a_v \cos\Delta_v &=& a_0 \sin\beta \cos\zeta\\
  \label{eq:scattering:corresp1b}
 a_v \sin\Delta_v &=& a_0 \cos\beta \sin\zeta
\end{eqnarray}
and 
\begin{eqnarray}
  \label{eq:scattering:corresp2a}
 a_h \cos\Delta_h &=& -a_0 \sin\beta \sin\zeta\\
  \label{eq:scattering:corresp2b}
 a_h \sin\Delta_h &=&  a_0 \cos\beta \cos\zeta
\end{eqnarray}
Multiplying Eq.~(\ref{eq:scattering:corresp1a}) with
Eq.~(\ref{eq:scattering:corresp2a}), and
Eq.~(\ref{eq:scattering:corresp1b}) with
Eq.~(\ref{eq:scattering:corresp2b}) and adding up the results, we get
\begin{equation}
  a_v a_h (\cos\Delta_v\cos\Delta_h + \sin\Delta_v\sin\Delta_h)
  = a_0^2 \sin\zeta\cos\zeta (\cos^2\beta - \sin^2\beta) 
\end{equation}
Using the addition theorems for sinusoidals and taking into account
that
$\Delta_v-\Delta_h = \delta_v-\delta_h$:
\begin{equation}
  \label{eq:scattering:sinzetacosbeta}
  \frac{a_v a_h}{a_0^2} \cos(\delta_v-\delta_h)
  = \frac{1}{2}\sin(2\zeta)\cos(2\beta)
\end{equation}
In a similar way, subtracting the product of
Eq.~(\ref{eq:scattering:corresp1b}) with
Eq.~(\ref{eq:scattering:corresp2a}) from the product of
Eq.~(\ref{eq:scattering:corresp1a}) with
Eq.~(\ref{eq:scattering:corresp2b}) and adding up the results, we get
\begin{equation}
  \label{eq:scattering:sinbeta}
  -\frac{a_v a_h}{a_0^2} \sin(\delta_v-\delta_h)
  = \frac{1}{2}\sin(2\beta)
\end{equation}
The above two equations tell us how to translate the amplitudes
($a_v$, $a_h$) and phases ($\delta_v$, $\delta_h$) of the vertical and
horizontal component of the electric field into the orientation and
shape of the ellipse that the tip of the electric field vector describes with
time.  We can obtain one further relation by subtracting the sum of
the squares of Eq.~(\ref{eq:scattering:corresp2a}) and
Eq.~(\ref{eq:scattering:corresp2b}) from the sum of the squares of
Eq.~(\ref{eq:scattering:corresp1a}) and
Eq.~(\ref{eq:scattering:corresp1b}):
\begin{equation}
  \label{eq:scattering:cosbetazeta}
 a_v^2 - a_h^2 =  -a_0^2 \cos(2\zeta)\cos(2\beta)
\end{equation}
Finally, we use the above 3 equations
(\ref{eq:scattering:sinzetacosbeta}), 
(\ref{eq:scattering:sinbeta}) and 
(\ref{eq:scattering:cosbetazeta}) to rewrite the Stokes parameters
(Eqs.~\ref{eq:scattering:stokesparam_alt_I}%
-\ref{eq:scattering:stokesparam_alt_V}) 
as
\begin{eqnarray}
  \label{eq:scattering:stokesparam_alt2_I}
 I &=& a_0^2\\
  \label{eq:scattering:stokesparam_alt2_Q}
 Q &=&  -a_0^2 \cos(2\zeta)\cos(2\beta)\\ 
  \label{eq:scattering:stokesparam_alt2_U}
 U &=& -a_o^2 \sin(2\zeta)\cos(2\beta)\\
  \label{eq:scattering:stokesparam_alt2_V}
 V &=& -a_0^2 \sin(2\beta)
\end{eqnarray}

Thus, we can get the orientation angle $\zeta$ of the ellipse from
\begin{equation}
  \label{eq:scattering:tan2zeta}
 \tan(2\zeta) = \frac{U}{V}
\end{equation}
Since $0 \leq 2\zeta < 360\degree$, there are 2 solutions for $\zeta$ for a
given pair $U,V$. This ambiguity is resolved by looking at
Eq.~\ref{eq:scattering:stokesparam_alt2_Q}, taking into account that
$|\beta| \leq 45\degree$ and thus $\cos(2\beta) \geq 0$:
The sign of $\cos(2\zeta)$ must be the same as the sign of $-Q$.

We get the ellipticity angle $\beta$  from
\begin{equation}
  \label{eq:scattering:tan2beta}
 \tan(2\beta) = - \frac{V}{(Q^2 + U^2)^{1/2}}  
\end{equation}

\texttt{FIXME: Measuring the Stokes parameters}

\texttt{FIXME: Partial Polarization}



\levelb{General radiative transfer equation}
%=====================================================================
\label{sec:scattering:general_rte}
 

The radiative transfer equation is \citep{mishchenko00:_light_scatt_nonsp_partic}: 
\begin{eqnarray}
     \frac{\DiffD \StoVec}{\DiffD s}(\PDir, \Frq, \Tmp) =
     -\ExtMat\DirFre\StoVec(\PDir, \Frq, \Tmp)+\AbsVec\DirFre
     \Planck(\Tmp, \Frq) \\ \nonumber
     +\int_{4\pi} \DiffD \PDir^\prime \SEfMat\DirFrePr
     \StoVec(\PDir^\prime, \Frq, \Tmp) 
\label{eq:scattering:RTE} 
\end{eqnarray} 
Formally it is a four dimensional vector equation where \StoVec
$(I,Q,U,V)$ is the Stokes vector.
Its first component $I$ is the monochromatic flux, the
second and third components $Q$ and $U$ describe the state of linear
polarization and the last component $V$ denotes the state of circular
polarization of the radiation propagating in the direction \PDir.
\DiffD $s$ denotes the pathlength element along direction \PDir. \Frq\ 
is the frequency of the radiation.

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 we
have 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. Accordingly \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 second term in Equation (\ref{eq:scattering:RTE})  describes the thermal
emission. Therefore \Planck\  is the Planck
function at temperature \Tmp 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 $\aSAbVec{p}_i$ is the single particle absorption cross section
vector. 

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 other directions \PDir$^\prime$   
into 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 are not interested in polarization effects we may use the scalar
radiative transfer equation, which is easier to handle than the vector
equation. It follows directly from the general equation. The scalar
extinction coefficient is the first element of the extinction matrix
$K_{11}$, the scalar absorption coefficient is given by the first
element of the absorption vector $a_1$ and the scalar scattering
coefficient is the first element of the scattering efficiency matrix
$Y_{11}$.\\
Substituting these definitions the scalar radiative transfer equation
into (\ref{eq:scattering:RTE}) is obtained:
\begin{eqnarray}
  \label{eq:scattering:scalar_rte}
\frac{\DiffD I}{\DiffD s} (\PDir, \Frq, \Tmp) = -K_{11}(\PDir, \Frq)I
(\PDir, \Frq, \Tmp) + a_1(\PDir, \Frq)\Planck(\Tmp,\Frq) + \\ \nonumber
 +\int_{4\pi} \DiffD \PDir^\prime Y_{11}\DirFrePr I(\PDir^\prime, \Frq, \Tmp)
\end{eqnarray}

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

It is not possible to find an analytical solution to the general
radiative transfer equation. We use an iterative method to find a
numerical solution.

\levelc{Scattering integral}
%===================================
\label{sec:scattering:scat_int}

The first step is solving the scattering integral
\begin{eqnarray}
  \ScaInt_0 \DirFrePr  = \int_{4\pi} \DiffD \PDir^\prime
  \SEfMat\DirFre \StoVec_0(\PDir^\prime, \Frq, T) 
\end{eqnarray}
using a first guess initial field for \StoVec$_0$\DirFrePr . The integration
has to be performed over all incident directions $\PDir^\prime$ for each
propagation direction \PDir{}. The result is the first guess scattered field $\ScaInt_0\DirFre$
which has the same dimension as the intensity field \StoVec$_0$\DirFrePr. 

In analogy the scalar scattering integral is obtained:
\begin{eqnarray}
  S_0^{int} \DirFrePr  = \int_{4\pi} \DiffD \PDir^\prime Y_{11}\DirFre I_0(\PDir^\prime, \Frq, T) 
\end{eqnarray}
As first guess field only the first Stokes component of the clear sky
field is needed.


\levelc{Radiative transfer with fixed scattered field}
%======================================
\label{sec:scattering:RTE}
With a fixed scattered field \ScaInt\DirFre{} the radiative transfer
Equation (\ref{eq:scattering:RTE}) can be written as
\begin{eqnarray}
     \frac{\DiffD \StoVec}{\DiffD s}(\PDir, \Frq, \Tmp) =
     -\ExtMat\DirFre\StoVec + \AbsVec\DirFre \Planck(\Tmp, \Frq)
     +\ScaInt_0(\PDir^\prime, \Frq, \Tmp)
\label{eq:rad_fs_sol}
\end{eqnarray} 

Inside a  grid box (see Section
\ref{sec:fm_defs:grids}) we assume that all
coefficients and the Planck function are constant. Then inside one
grid box Equation (\ref{eq:rad_fs_sol}) is a linear differential
equation with constant coefficients which can be solved analytically
leading to the result  
\begin{eqnarray}
  \label{eq:scattering:RTE_sol}
  \StoVec_1 = e^{-\ExtMat s}\cdot\StoVec_0 + (\IdnMat - e^{-\ExtMat
    s}) \ExtMat\Inv (\AbsVec \Planck + \ScaInt_0)
\end{eqnarray}
where \StoVec$_0$ is given by the initial condition
\begin{eqnarray}
  \StoVec_0 =  \StoVec(s = 0).
\end{eqnarray}
Here we take for \StoVec$_0$ the initial field and we find the first
iteration radiation field \StoVec$_1$.

The scalar radiative transfer equation
(\ref{eq:scattering:scalar_rte})  with fixed scattering integral
is accordingly
\begin{eqnarray}
  \label{eq:scattering:scalar_rte_scatint}
\frac{\DiffD I}{\DiffD s} (\PDir, \Frq, \Tmp) = -K_{11}(\PDir, \Frq)I
(\PDir, \Frq, \Tmp) + a_1(\PDir, \Frq)\Planck(\Tmp,\Frq) + S_0^{int}
(\PDir, \Frq, \Tmp).
\end{eqnarray} 
Assuming constant coefficients the analytical solution is
\begin{eqnarray}
   \label{eq:scattering:scalar_rte_sol}
I_1 = I_0 e^{-K_{11}s} + \frac{a_1
  B + S_0^{int}}{K_11}\left(1-e^{-K_{11}s}\right)
\end{eqnarray}
where $I_0$ again is specified by the initial condition:
\begin{eqnarray}
  I_0 = I(s=0)
\end{eqnarray}
Applying this equation leads to first iteration scalar intensity
field.  



\levelc{Iteration and Convergence}
%===============================
\label{sec:scattering:conv}

The procedure above is repeated, i.e. \StoVec$_1$ is taken to evaluate
the scattering integral leading to the first iterated scattered field
$\ScaInt_1$. Then Equation (\ref{eq:scattering:RTE_sol}) is again used to
compute the
second iteration field \StoVec$_2$. \\
After each iteration the convergence is checked. If 
\begin{eqnarray}
|\StoVec_m \DirFre -  \StoVec_{m-1} \DirFre| < {\bf \epsilon}
\end{eqnarray}
a solution to the vector radiative transfer Equation (\ref{eq:scattering:RTE})
has been
found. An analogous convergence test is done for the scalar case.

\levelb{Database for the single scattering properties}
%============================================================
\label{sec:scattering:database}

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% [\cite{Mishchenko:98}]
%[\cite{Mishchenko:00}] //check in references!!!!  
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}
The amplitude matrix depends on the directions of incidence  $(\theta,
\phi)$  and
scattering $(\theta', \phi')$ as  well as on size, morphology and composition of the
scattering material, also 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. In general, the
amplitude matrix is also frequency dependent.

The amplitude matrix contains all information about the optical
properties of the scattering medium. Extinction cross section
matrix, absorption cross section vector and phase matrix can be
obtained from the amplitude matrix (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. 

For each particle type there is one file in the database.
The particle type depends on the particle size, the particle shape and the
orientation. We consider only ice crystals, which means the material of the
particles is always the same. If there is a cloud consisting of
different particles, the scattering properties for that particles
distribution can be derived from the single scattering properties.

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}$).

You can read the database using the function
\wsmindex{read\_amp\_fromDb}.\\
Output of the function is a tensor containing amplitude matrices for
a fixed angle combination stored in the variable
\wsvindex{amp\_mat}. Its dimension is [$N_{sp}, N_{I}, N_{I}$]
where $N_{pt}$ is the number of particle
types  and $N_{I}$ denotes the Stokes dimension.\\
Input to this method is a fixed angle
combination  (\artsstyle{za},\artsstyle{aa},\artsstyle{za}' and
\artsstyle{aa}') for propagation and scattered direction and the vector
\artsstyle{scat\_pt} which contains indices for the  particle types
which are considered in the calculation. This vector has to be defined
in the control file.


\levelb{Get local atmospheric properties from atmospheric fields}
%==================================================================
\label{sec:scattering:gen_atmprop}

\levelb{Generate extinction matrix}
%==================================================================
\label{sec:scattering:gen_ext}

\levelc{Agenda for calculating the extinction matrix}
%==================================================================
\label{sec:scattering:ext_mat_agenda}

The agenda \wsvindex{ext\_mat\_agenda} calculates the total extinction
matrix \ExtMat{}
for each grid point. It is the physical extinction matrix, that means
that it includes the gaseous extinction and the particle
extinction.\\
Output is the variable \wsvindex{ext\_mat} which is given
as a $N_I$x $N_I$ Matrix, i.e. one dimension for each Stokes component. If the
Stokes dimension is only one, only the element $\ExtMat_{11}$ will be
nonzero, it corresponds to the scalar extinction cross section. \\
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.\\
The agenda can for example be defined in the following manner:

\begin{minipage}{0.9\hsize}
\begin{verbatim}
AgendaDefine(ext_mat_agenda)
{
ext_mat_partCalc{} // Function to calculate particle
                   // extinction.
ext_mat_gasCalc{}  // Function to calculate gaseous
                   // extinction. 
ext_matCalc{}      // Function which sums up
                   // particle extinction and 
                   // gaseous extinction.
}
\end{verbatim}
\end{minipage}

\vspace{2ex}
\noindent
The workspace methods \wsmindex{ext\_mat\_partCalc} and
\wsmindex{ext\_mat\_gasCalc} are explained in Sections
\ref{sec:scattering:ext_mat_part} and \ref{sec:scattering:ext_mat_gas}
respectively.

\levelc{Generate particle extinction matrix}
%==================================================================0
\label{sec:scattering:ext_mat_part}

The workspace method \wsmindex{ext\_mat\_partCalc} gives as output the
particle extinction matrix \aExtMat{p} stored in the variable
\wsvindex{ext\_mat\_part} which has the chosen Stokes dimension. So in
the most general case it is a 4x4 matrix.\\
First input is the workspace variable \wsvindex{ext\_mat\_spt} (see Section
\ref{sec:scattering:ext_mat_spt}). The second input variable is the
local particle number densities \artsstyle{pnd\_field} for all particle types,
that means \artsstyle{pnd\_field} is a Tensor4 with dimension 
[$N_{pt}, N_{p}, N_{lat}, N_{lon}$]. The other workspace variable
inputs required are \wsvindex{atmosphere\_dim}, \wsvindex{scat\_p\_index},
\wsvindex{scat\_lat\_index}, and \wsvindex{scat\_lon\_index}. \\
The method sums up the extinction matrices for all particle types
weighted with the particle number density:
\begin{eqnarray}
  \aExtMat{p} = \sum_i \PDen_i \aSExMat{p}_i 
\end{eqnarray}


\levelc{Generate gaseous extinction matrix}
%===============================================================
\label{sec:scattering:ext_mat_gas}

\levelc{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}


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

\levelc{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}

\levelc{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_{21} &=& \frac{1}{2}(|S_{11}|^2+|S_{12}|^2-|S_{21}|^2-|S_{22}|^2)\\
   Z_{31} &=& -Re(S_{11}S_{21}^*+S_{22}S_{12}^*)\\
   Z_{41} &=& -Im(S_{21}S_{11}^*+S_{22}S_{12}^*)\\
\cdots ... {\rm to be completed}
\end{eqnarray}






\levelb{Generate absorption coefficient vector}
%==========================================================================
\label{sec:scattering:abs_vec}

Agendas and methods to generate the absorption coefficient vector are
defined in analogy to those generating the extinction coefficient matrix.

\levelc{Agenda for calculating the absorption vector}
%==================================================================
\label{sec:scattering:abs_vec_agenda}

The agenda \wsvindex{abs\_vec\_agenda} calculates the total absorption
matrix \AbsVec{}
for each grid point. It is the physical absorption vector including
particle absorption and gaseous absorption.\\
Output is the variable \wsvindex{abs\_vec} which is given
as a $N_I$ component vector, i.e. one dimension for each Stokes component. If the
Stokes dimension is only one, it will only have one component
corresponding to the scalar absorption coefficient. \\
Input for this agenda are the following variables:
the particle number density \wsvindex{pnd\_field} and 
\wsvindex{abs\_vec\_spt}, which contains the absorption vectors  for
each chosen particle type. How this variable is obtained
is described in Section \ref{sec:scattering:abs_vec_spt}. For calculating 
gaseous absorption the volume mixing ratio (\wsvindex{vmr}), 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.\\
The definition of the agenda is listed below:

\vspace{2ex} 
\begin{minipage}[h]{.9\hsize}
\begin{verbatim}
AgendaDefine(abs_vec_agenda)
{
abs_vec_partCalc{}  // Function to calculate particle
                    // absorption.
abs_vec_gasCalc{}   // Function to calculate gaseous
                    // absorption.
abs_vecCalc{}       // Function which sums up 
                    // particle absorption and
                    // gaseous absorption.
}
\end{verbatim}
\end{minipage}
\vspace{2ex} 

\noindent
The workspace methods \wsmindex{abs\_vec\_partCalc} and
\wsmindex{abs\_vec\_gasCalc} are described in the sections
\ref{sec:scattering:abs_vec_part} and
\ref{sec:scattering:abs_vec_gas}. 


\levelc{Calculate particle absorption coefficient vector}
%=======================================================================
\label{sec:scattering:abs_vec_part}

The workspace method \wsmindex{abs\_vec\_partCalc} gives as output the
particle absorption vector \aAbsVec{p} stored in the variable
\wsvindex{abs\_vec\_part} which has the chosen Stokes dimension, so in
general a four component vector.\\
First input is the workspace variable \wsvindex{abs\_vec\_spt} (see Section
\ref{sec:scattering:abs_vec_spt}). The second input variable is the
local particle number densities \artsstyle{pnd\_field} for all particle types,
that means \artsstyle{pnd\_field} is a Tensor4 with dimension 
[$N_{pt}, N_{p}, N_{lat}, N_{lon}$]. The other workspace variable
inputs required are \wsvindex{atmosphere\_dim}, \wsvindex{scat\_p\_index},
\wsvindex{scat\_lat\_index}, and \wsvindex{scat\_lon\_index}. \\  
The function sums up the absorption vectors for all particle types
weighted with the particle number density.
\begin{eqnarray}
  \aAbsVec{p} = \sum_i \PDen_i \aSAbVec{p}
\end{eqnarray}


\levelc{Calculate gaseous absorption coefficient vector}
%=======================================================================
\label{sec:scattering:abs_vec_gas}


\levelc{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}
These equations show, that the first column of extinction and phase
matrix are needed. So it is obligatory to execute the agendas
\funcindex{ext\_mat\_agenda} and 
\funcindex{pha\_mat\_agenda} before executing
\funcindex{abs\_vec\_agenda}. 



\levelb{Radiative transfer inside the cloudbox}
%==============================================================
\label{sec:scattering:scat_meth_rt}

The scattering radiative transfer calculation is performed
independently from the clear sky radiative transfer calculation. It is
limited only to a small part of the atmosphere defined by the cloud box.

\levelc{Define the cloudbox}
%==========================================
\label{sec:scattering:cloudbox}

\leveld{3D geometry}
%==========================================

The concept of a cloud box is explained in Section
\ref{sec:fm_defs:cloudbox}. It is activated by setting the flag
\wsvindex{cloudbox\_on} to 1. The
limits of the cloud box are stored in \wsvindex{cloudbox\_limits}
which in an array of indices, containing the lower and the upper
pressure index, the lower and upper latitude index and the lower and
upper longitude index. 

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}. 
Furthermore angle grids have to be defined for the scattering calculations  
as the radiation field, the scattering efficiency matrix etc. depend
on the propagation and incident directions. For this purpose  the workspace  
variables \wsvindex{scat\_za\_grid} and \wsvindex{scat\_aa\_grid}
have to be defined.

The interface between the clear sky RT calculation and the cloudbox RT
calculation is the radiation field on the boundary of the cloudbox
which is stored
in the variables \wsvindex{scat\_i\_p}, \wsvindex{scat\_i\_lat} and
\wsvindex{scat\_i\_lon}. The dimensions are 
\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}  


\leveld{1D geometry}
%=========================================================
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} (???? really ????). 

As interface between the clear sky calculation and the cloud box calculation only 
  \begin{center}
  \artsstyle{scat\_i\_p} = \artsstyle{scat\_i\_p} (\Frq, 2(two surfaces), 1 , 1, \ScaZa,
1, i$_I$)
\end{center}
is required. The other interface variables are empty. 

\levelc{Scattering main function}
%=====================================
\label{sec:scattering:main_function}

The scattering main function \wsmindex{scatCalc} requires the clear sky radiation 
field on the cloud box boundary as input.
It loops over the frequencies and in each loop it executes the agenda 
\wsvindex{scat\_mono\_agenda}
(see Section \ref{sec:scattering:scat_mono_ag}).\\
Output of \wsmindex{ScatCalc} is the scattered radiation field on the
boundary also stored in the variables \wsvindex{scat\_i\_p},
\wsvindex{scat\_i\_lat}
and \wsvindex{scat\_i\_lon}.


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

The agenda \wsvindex{scat\_mono\_agenda} solves the radiative transfer
equation inside the cloudbox for one frequency specified
by the frequency index \wsvindex{scat\_f\_index}. 
The agenda consists of a number of methods which have to be executed in
the correct order.
The agenda is set in the control file. To solve the RTE iteratively it
is set in the following manner:

\begin{verbatim}
AgendaSet(scat_mono_agenda) {
  i_fieldSet{}
  i_fieldIterate{}
  scat_iPut{}
}
\end{verbatim}

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

 
\leveld{Initial field}
%=======================================================================
The workspace method \wsmindex{i\_fieldSet} 
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 
which is stored in the workspace variable
\begin{center}
 \wsvindex{i\_field} = \wsvindex{i\_field} (\Prs, \Lat, \Lon, \ScaZa,
\ScaAa, i$_I$). 
\end{center}


\leveld{Iteration}
%=====================================================================
The method \wsmindex{i\_fieldIterate} solves the RTE using the
iterative method.  The function has
included switches to adapt automatically to the atmospheric
dimensionality specified
in 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{sca\_fieldCalc}
    (see Sections \ref{sec:scattering:sca_fieldCalc} and
    \ref{sec:scattering:scat_int}).
\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:RTE}).
\item Do the convergence test (see Sections
  \ref{sec:scattering:conv_method} and \ref{sec:scattering:conv}).
\end{itemize}
The solution of the radiative transfer equation is returned as output
using the variable \artsstyle{i\_field}. 

\leveld{Generate full radiation field inside cloud box}
The variable \artsstyle{i\_field} contains the radiation field only
for one frequency. So there has to be a method which puts
\artsstyle{i\_field} into the variables \wsvindex{scat\_i\_p},
\wsvindex{scat\_i\_lat} and \wsvindex{scat\_i\_lon}. The method
\funcindex{scat\_iPut} is doing this task.
It has to use the frequency index to know where
to put \artsstyle{i\_field}.
 

\levelc{Method to compute the scattered field}
%=========================================
\label{sec:scattering:sca_fieldCalc}
The method \wsmindex{scat\_integralCalc} does exactly what is
described in Section 
\ref{sec:scattering:scat_int}.  The output of this method is
\wsvindex{scat\_field} which has the same dimension as
\artsstyle{i\_field}.  This 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}, phase matrix for single 
particle \artsstyle{pha\_mat\_spt}, the particle number density field
\artsstyle{pnd\_field} are also 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} ). \\ 
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{Methods for solving the RT with fixed scattering integral}
%==================================================================
\label{sec:scattering:RT_methods}

\leveld{Workspace method \artsstyle{i\_fieldUpdate1D}}
%===================================================================
For a 1D atmosphere the method  \wsmindex{i\_fieldUpdate1D} is used to
evaluate the RT Equation (\ref{eq:scattering:RTE_sol}) and update the
radiation field \wsvindex{i\_field} for each iteration. 
The function loops over pressures given in \wsvindex{p\_grid} and
propagation directions given in \wsvindex{scat\_za\_grid}.  
The lower and upper indices for the loop over pressures are passed
into the function by the workspace variable
\wsvindex{cloudbox\_limits}. Inside the loops the following steps are
performed:
\begin{itemize}
\item Calculate the coefficients of the RT Equation using appropriate agendas 
  (\ref{eq:scattering:RTE_sol}):
  \begin{itemize}
  \item Extinction coefficient matrix: Execute \wsvindex{ext\_mat\_agenda} (see
    Section \ref{sec:scattering:ext_mat_agenda}).
  \item Scattering efficiency matrix: Execute \wsvindex{sca\_matCalc} (see
    Section \ref{sec:scattering:pha_mat}).
  \item Absorption coefficient vector:  Execute \wsvindex{abs\_vec\_agenda} (see
    Section \ref{sec:scattering:abs_vec_agenda}).
  \end{itemize}
\item A  propagation path starting at the current point specified by
  the internal variable \artsstyle{p\_index} is initialized. Its
  direction is specified by \artsstyle{scat\_za\_index}. Using the
  \wsvindex{ppath\_step\_agenda} (cf. Section
  \ref{sec:ppath:stepcalc}) 
  the intersection point with the next
  layer is determined as well as the pathlength d{\bf s} from the starting
  point to the intersection point.  
\item Solve Equation (\ref{eq:scattering:RTE_sol}) using  the agenda
  \wsvindex{scat\_rte\_agenda} (see description below).
\end{itemize}

\leveld{Agenda \wsvindex{scat\_rte\_agenda}}
%=================================================================
In the agenda \wsvindex{scat\_rte\_agenda} the method for doing a
radiative step calculation through one grid cell for a 3D atmosphere
or one layer for a 1D atmosphere is defined. The following methods can
be chosen:

\begin{itemize}
\item {\bf Method \artsstyle{stokes\_vecGeneral}:}\\
This method is the most general method. It can be applied for all
Stokes dimensions but it should not be used for the scalar equation as
for this case the method is numerically inefficient. 

There are two terms in the radiative transfer equation
\ref{eq:scattering:RTE_sol} which 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 Method \artsstyle{stokes\_vecScalar}:}\\
For the scalar case it is straightforward to compute a radiative
transfer step through one grid cell/layer. Only equation
(\ref{eq:scattering:scalar_rte_sol}) has to be computed. 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.
\end{itemize}



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





 


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