%
% 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
%  \Ind           & -                 & vector/matrix/tensor index           \\
%  \aInd{\Lat}    & -                 & the \Ind:th latitude                 \\
%  \VctLng        & -                 & vector length or size of matrix/tensor for a dimension \\
%  \aVctLng{\Lat} & -                 & length of the latitude grid \\
%  \Prs           & \artsstyle{p}     & pressure                             \\
%  \PrsAlt        & \artsstyle{pz}    & pressure altitude                    \\
%  \Rds           & \artsstyle{r}     & radius from the geoid centre         \\
%  \Alt           & \artsstyle{z}     & geometrical altitude above the geoid \\
%  \Lat           & \artsstyle{alpha} & latitude                             \\
%  \Lon           & \artsstyle{beta}  & longitude                            \\
%  \ZntAng        & \artsstyle{psi}   & zenith angle                         \\
%  \AzmAng        & \artsstyle{omega} & azimuthal angle                      \\
% \label{symtable:ppath_defs}     
%\stopsymbols

\startsymbols
\StoVec       & \artsstyle{i\_field}       & monochromatic intensity field\/ Stokes Vector\\
              & \artsstyle{i\_field\_old}  & needed only internally
              for the iteration\\
\ExtMat       & -                        & total extinction matrix \\
\SExMat        & \artsstyle{ext\_mat} & extinction martrix 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

\newcommand{\DirFre} {(\PDir, \Frq)}
\newcommand{\DirFrePr} {\ensuremath{(\PDir^\prime, \Frq)}}


\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+\AbsVec\DirFre \Planck(\Tmp, \Frq)
     +\int_{4\pi} \DiffD \PDir^\prime \SEfMat\DirFrePr \StoVec\DirFrePr 
\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
polarisation and the last component $V$ denotes the state of circular
polarisation 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 travelling 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 K can be written as
a sum of two matrices, the first for particle extinction ( \aExtMat{p})
and the second for gaseous extinction (\aExtMat{g}). Then
\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 T 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{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\DirFrePr 
\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. 

\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(\PDir, \Frq, \Tmp)}{\DiffD s} =
     -\ExtMat\DirFre\StoVec + \AbsVec\DirFre \Planck(\Tmp, \Frq)
     +\ScaInt_0\DirFre
\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 aone
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 field \StoVec$_1$.

\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 (\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.

\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 the 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:sca_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. 


%=============================this part has to be updated !!!!!!!!
For each considered frequency there is one file for each particle
type in the database. The particle
type depends on the particle size, the particle shape and the
orientation. We consider only ice crystals, so 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 propagation direction (\ScaZa,
\ScaAa) and each scattered direction (\ScaZa', \ScaAa') 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
\wsfindex{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{Extract 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 \wsfindex{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}, 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{vmr}), the
pressure  (\wsvindex{p}), the temperature  (\wsvindex{T}) and the
propagation direction given by zenith angle (\wsvindex{za}) and
azimuth angle (\wsvindex{aa}) are needed.\\
The agenda can for example be defined in the following manner:

\verb|AgendaDefine(ext_mat_agenda)|\\
\begin{tabular}[h]{l l}
\verb|{ |\\
\verb|ext_mat_partCalc{}| & \verb| //function to calculate particle |\\
& \verb|extinction | (see section \ref{sec:scattering:ext_mat_part}).\\
\verb|ext_mat_gasCalc{}| & \verb| //function to calculate gaseous | \\
& \verb|extinction | (see section \ref{sec:scattering:ext_mat_gas}).\\ 
\verb|ext_matCalc{}| &\verb| // function which sums up | \\
& \verb| particle extinction and gaseous extinction.|\\
\verb| }|
\end{tabular}

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

The workspace method \wsfindex{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.\\
Input is first the workspace variable \wsvindex{ext\_mat\_spt} (see section
\ref{sec:scattering:ext_mat_spt}). The second input variable are the
local particle number densities \artsstyle{pnd} for all particle types,
that means \artsstyle{pnd} is of dimension $N_{pt}$ (number of
particle types). \\
The function 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 derived using the workspace method
\wsfindex{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.\\

% =================update !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
As input variables it needs the amplitude matrix (\wsvindex{amp\_mat})
and the propagation direction given by the local zenith and azimuth
angle (\artsstyle{za} and \artsstyle{aa}). \wsvindex{amp\_mat} is
obtained using the reading routine \wsfindex{get\_amp\_fromDb} which is
explained in section \ref{sec:scattering:database}.\\
%============================================================

The extinction cross section matrices \SExMat{} are
calculated for all patricle 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_sca_mat}

\levelc{Calculate the total scattering efficiency matrix}
%=======================================================================
\label{sec:scattering:sca_mat}

Output of the workspace method \wsfindex{sca\_matCalc} is the phase
matrix stored in the variable \wsvindex{sca\_mat}  with dimension
according to the chosen Stokes
dimension, so in general it is a 4x4 matrix.\\
Input is first the workspace variable \wsvindex{sca\_mat\_spt} (see section
\ref{sec:scattering:sca_mat_spt}). The second input variable is the
local particle number density vector \artsstyle{pnd} 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:sca_mat_spt}

From the amplitude matrix the phase matrix for each
particle type is calculated using the workspace method
\wsfindex{pha\_mat\_sptCalc}.\\
Output is a tensor including phase matrices for all particle
types. Therefore it is  of  dimension
[$N_{pt}, N_{I}, N_{I}$], where $N_{pt}$ is the number of particle
types  and $N_{I}$ denotes the Stokes dimension.\\
As input variables it needs the amplitude matrix (\wsvindex{amp\_mat})
and the propagation direction given by the local zenith and azimuth
angle (\artsstyle{za} and \artsstyle{aa}). Furthermore it needs the
angles describing the scattering direction  (\artsstyle{za\_inc} and
\artsstyle{aa\_inc}). 
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 \wsfindex{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} 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}), the temperature  (\wsvindex{T}) and the
propagation direction given by zenith angle (\wsvindex{za}) and
azimuth angle (\wsvindex{aa}) are needed.\\
The definition of the agenda is listed below:

\verb|AgendaDefine(abs_vec_agenda)|\\
\begin{tabular}[h]{l l}
\verb|{ |\\
\verb|abs_vec_partCalc{}| & \verb| //function to calculate particle |\\
& \verb|absorption | (see section \ref{sec:scattering:abs_vec_part}).\\
\verb|abs_vec_gasCalc{}| & \verb| //function to calculate gaseous | \\
& \verb|absorption | (see section \ref{sec:scattering:abs_vec_gas}).\\ 
\verb|abs_vecCalc{}| &\verb| // function which sums up | \\
& \verb| particle absorption and gaseous absorption.|\\
\verb| }|
\end{tabular}



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

The workspace method \wsfindex{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.\\
Input is first the workspace variable \wsvindex{abs\_vec\_spt} (see section
\ref{sec:scattering:abs_vec_spt}). The second input variable are the
local particle number densities \artsstyle{pnd} for all particle types
chosen in the control file.
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
\wsfindex{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 variables it needs the amplitude matrix (\wsvindex{amp\_mat})
and the propagation direction given by the local zenith and azimuth
angle (\artsstyle{za} and \artsstyle{aa}). Furthermore 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}
This 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{sca\_mat\_agenda} before executing
\funcindex{abs\_vec\_agenda}. 



\levelb{Methods and Agendas for the RT inside the Cloud Box}
%==============================================================
\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 cloud box}
%==========================================
\label{sec:scattering:cloudbox}

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

The concept of a cloud box is explained in
\ref{sec:fm_defs:cloudbox}. It is activated by setting the variable
\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 upper and lower latitude index and the lower and
upper longitude index. 

The calculations inside the cloud box are performed an 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 direction. 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 scattering
box. The radiation field \StoVec\ on the boundary of the scattering box 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}

\leveld{1D Geometry}
%=========================================================
For 1D scattering calculations the same workspace variables are 
used as for the 3D calculations. The dimensions which are not needed are empty. 
As special grid inside the scattering box we need to specify only  
\wsvindex{scat\_za\_grid}. 

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 \wsfindex{scat\_iterateCalc} 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 \wsfindex{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{Agenda \wsfindex{scat\_mono\_agenda}}
%===================================
\label{sec:scattering:scat_mono_ag}

The iteration scheme outlined in \ref{sec:scattering:solution_rte} is
performed by \wsfindex{scat\_mono\_agenda} for one frequency specified
by the frequency index \artsstyle{scat\_f\_index}. 
The agenda consists of a number of methods which have to be executed in
the correct order.


\leveld{Initial field}
Input to the agenda is the radiation field on the cloud
box boundary.
Using a linear 3D interpolation scheme, we obtain a
radiation field on all grid points inside the cloud box.
This can be taken as a first guess for the iteration. The method
\funcindex{i\_fieldInterpolateFromBoundary} performs this
interpolation. It 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 internal variable
\begin{center}
  \artsstyle{i\_field} $\in$ (\Prs, \Lat, \Lon, \ScaZa,
\ScaAa). 
\end{center}

\leveld{Iteration}

The method \funcindex{i\_fieldCalc} performs the following steps 
until convergence is obtained. As input the initial field stored in
\artsstyle{i\_field} is required. Furthermore the method needs two
agendas, the \wsfindex{scattered\_field\_agenda} and the
\wsfindex{i\_field\_update\_agenda}.

\begin{itemize}
\item \artsstyle{i\_field} is copied to \artsstyle{i\_field\_old}.
\item The scattered field is calculated using
  \wsfindex{scattered\_field\_agenda}
    (see sections \ref{sec:scattering:scat_field_ag} and
    \ref{sec:scattering:scat_int}).
\item Calculate the new radiation field using
  \wsfindex{i\_field\_update\_agenda}
 (see sections \ref{sec:scattering:RT_ag} 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\_iPutIn} is doing this task.
Of course it has to use the frequency index to know where
to put \artsstyle{i\_field}.
 

\levelc{Agenda to compute the Scattered Field}
%=========================================
\label{sec:scattering:scat_field_ag}

\levelc{Agenda for RT with fixed Scattering Integral}
%=====================================
\label{sec:scattering:RT_ag}
The \wsfindex{i\_field\_update\_agenda} needs the variables
\artsstyle{i\_field\_old}, \artsstyle{scattered\_field}, \wsvindex{ext\_mat}
and  \wsvindex{abs\_vec} as input.\\ 
According to equation(\ref{eq:scattering:RTE_sol}) it computes the new
radiation field. The agenda starts with a loop over all grid points
inside the scattering box. Then for each grid point the new intensity
vector
 \StoVec\  is obtained for all directions by calculating the
radiative transfer through the adjacent grid cells. We consider paths
starting at the grid point into all directions and ending at the next
intersection points with a grid cell boundary. The paths are calculated
using the method \ref{sec:ppath}.\\
 The intensity vector at
the intersection points \StoVec{0} is obtained by linear interpolation
between the grid points using the interpolation routine (*** give 
reference to section here ***)
%\ref{sec:interpolation}. \\
Computing the matrix exponential is done using the method
\wsfindex{matrix\_exp} explained in section
\ref{sec:lin_alg:mat_exp}.\\
 The term
\ExtMat$^{-1}$[\AbsVec\Planck+\ScaInt] can be computed numerically  
performing a LU decomposition. This is implemented in the workspace
methods \wsfindex{ludcmp} and \wsfindex{lubacksub} which are described
in the sections \ref{sec:lin_alg:backsub} and
\ref{sec:lin_alg:lu_decomp}. \\
Output of the \wsfindex{i\_field\_update\_agenda} is the variable 
\artsstyle{i\_field} which contains the new radiation field at all
grid points in the cloudbox for all directions defined by 
\wsvindex{scat\_za\_grid} and \wsvindex{scat\_aa\_grid}.

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





 


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "uguide"
%%% End: 
