%
% 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, \Tmp)}
\newcommand{\DirFrePr} {\ensuremath{(\PDir, \PDir^\prime, \Frq, \Tmp)}}


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

The radiative transfer equation for a medium with thermal emission
comprising sparseley and randomly distributed, spherical or
arbitrarily oriented nonspherical particles 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 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 you could use the electromagetic 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 convienient 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$ signify the state of linear
polarization and the last component $V$ signifies the state of circular
polarization. More details about the physical meaning of the Stokes
vector and about possibilities to measure the Stokes components
can be
found in Chapter \ref{sec:polarization}.\\
In Equation  (\ref{eq:scattering:RTE})
\PDir\ denotes the propagation direction and
\DiffD $s$ the pathlength element along this direction. The radiation
field depends on frequency \Frq\ and temperature \Tmp.

\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 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 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, \Tmp) = \sum_i{\PDen_i\PhaMat_i(\PDir,
  \PDir', \Frq, \Tmp)}.
\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 nonspherical particles in the scattering
medium. If there are nonspherical 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 scattering 
matrix, commonly only its first element, which is usually called the
phase function, is used although there is no physical justification for
doing this.

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, \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 vector 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. For this purpose the radiation field has to be
discretised. This is done by defining cloudbox (cp. Section
\ref{sec:fm_defs:cloudbox}). In the 3D atmosphere the cloudbox
consists of a specific number of grid points which are the corners of
the so called grid cells. Inside the grid cells we assume that the
extinction matrix, the absorption vector and the scattering matrix are
constant. 

\begin{itemize}
\item {\bf Scattering integral}

The first step  of the iterative method is solving the scattering integral
\begin{eqnarray}
  \ScaInt^{(0)} = \int_{4\pi} \DiffD \PDir^\prime
  \SEfMat \StoVec^{(0)}
\end{eqnarray}
using a first guess initial field \StoVec$^{(0)}$. To do the
integration we also need a discretised angular grid to define the
directions. The initial field depends on the
position inside the cloudbox and on the incident direction of the
radiation. For
convenience the frequency and temperature dependence is omitted in
this equation and in the following equations. 

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 $\ScaInt^{(0)}$
which has the same dimension as the intensity field \StoVec$^{(0)}$. 




\item{\bf Radiative transfer with fixed scattered field}

We can substitute  $\ScaInt^{(0)}$ as a fixed scattered field into the
radiative transfer 
Equation (\ref{eq:scattering:RTE}):   
\begin{eqnarray}
     \frac{\DiffD \StoVec^{(1)}}{\DiffD s} =
     -\ExtMat \StoVec^{(1)} + \AbsVec \Planck
     +\ScaInt^{(0)}
\label{eq:rad_fs_sol}
\end{eqnarray} 
As the coefficients are constant in each
grid cell Equation (\ref{eq:rad_fs_sol}) is formally a linear differential
equation with constant coefficients which can be solved analytically
leading to the result (cp. FIXME Appendix: exact solution method) 
\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}
Here \IdnMat\ denotes the identity matrix. 
$\StoVec^{(0)}$ is given by the initial condition, which in this
case is given by the first guess initial field. For each grid
point inside the cloudbox the intersection points with the adjacent
grid cell boundaries are calculated and then the initial field is
interpolated on these points.

$\StoVec^{(1)}$ is called the first iteration field.


\item {\bf Iteration}

The first iteration field $\StoVec^{(1)}$ 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 field $\StoVec^{(2)}$
\begin{eqnarray}
     \frac{\DiffD \StoVec^{(2)}}{\DiffD s} =
     -\ExtMat \StoVec^{(2)} + \AbsVec \Planck
     +\ScaInt^{(1)}
\end{eqnarray} 
which is solved by
\begin{eqnarray}
   \StoVec^{(2)} = e^{-\ExtMat s}\cdot\StoVec^{(1)} + (\IdnMat - e^{-\ExtMat
    s}) \ExtMat\Inv (\AbsVec \Planck + \ScaInt^{(1)})
\end{eqnarray}
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 is given by
\begin{eqnarray}
     \frac{\DiffD \StoVec^{(n)}}{\DiffD s} =
     -\ExtMat \StoVec^{(n)} + \AbsVec \Planck
     +\ScaInt^{(n-1)}.
\end{eqnarray} 
So the  $n$-th order iteration field is given by:
\begin{eqnarray}
   \StoVec^{(n)} = e^{-\ExtMat s}\cdot\StoVec^{(n-1)} + (\IdnMat - e^{-\ExtMat
    s}) \ExtMat\Inv (\AbsVec \Planck + \ScaInt^{(n-1)})
\end{eqnarray}


\item {\bf Convergence test}

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. 
\end{itemize}


\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
leading to the result
\begin{eqnarray}
   \label{eq:scattering:scalar_rte_sol}
I^{(1)} = I^{(0)} e^{-K_{11}s} + \frac{a_1
  B + S^{(0)}}{K_{11}}\left(1-e^{-K_{11}s}\right)
\end{eqnarray}
where $I^{(0)}$ again is specified by the initial field.

Applying this equation leads to first iteration scalar intensity
field  $I^{(1)}$.  

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

There are 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
[\cite{Mishchenko:98}] and [\cite{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 are 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 is one 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. 
All 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">
<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="18" nbooks="1" npages="1" 
                                 nrows="1" ncols="1">
  0.          
  10.
  20.         // incoming zenith angle grid
   .
   .
   .
 </Tensor6>
 <Tensor6 nvitrines="1" nshelves="1" nbooks="36" npages="1"
                                 nrows="1" ncols="1">
   0.
  10.
   .          // incoming azimuth angle grid   
   .
   .
 </Tensor6>
 <Tensor6 nvitrines="1" nshelves="1" nbooks="1" npages="18"
                                 nrows="1" ncols="1">
  0.
  10.
   .          // scattered zenith angle grid
   .
   .
 </Tensor6>
 <Tensor6 nvitrines="1" nshelves="1" nbooks="1" npages="18"
                                 nrows="1" 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="18" nbooks="36" 
                        npages="18" nrows="36" 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. The data files contain more 
information about the particles, for example the particle shape, the
particle shape and the aspect ratio for non-spherical particles.


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

%==================================================================
\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. This is the physical extinction matrix, that means
that it includes gaseous extinction and 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:

\vspace{1ex}
\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}). Second input is the
local particle number density field for all particle types
\artsstyle{pnd\_field}, which
 is a Tensor4 with dimension 
[$N_{pt}, N_{p}, N_{lat}, N_{lon}$]. Furthermore the  workspace
variables
\wsvindex{atmosphere\_dim}, \wsvindex{scat\_p\_index},
\wsvindex{scat\_lat\_index} and \wsvindex{scat\_lon\_index} are
required as input. \\
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_{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_{41} &=& -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_{41} &=& -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}






\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 density field \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 agenda
\funcindex{ext\_mat\_agenda} and the method
\funcindex{pha\_matCalc} 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 is 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}. 

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: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}. 

\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:solution_rte}.  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}, the phase matrix for single 
particle \artsstyle{pha\_mat\_spt}, 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} ). \\ 
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. 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 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.
\end{itemize}


\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
