\graphicspath{{Figs/scattering/}}

\chapter{Scattering calculations -- The DOIT module}
 \label{sec:scattering:doit}

\starthistory
 020601 & Created and written by Claudia Emde\\ 
 050223 & Rewritten by Claudia Emde, mostly taken from
 Chapter 4 of Claudia's PhD thesis \\
 050929 & Included technical part, example contol file \\
 141009 & Moved general model parts to theory guide
\stophistory

The \textindex{Discrete Ordinate ITerative (DOIT) method} is one of the
scattering algorithms in ARTS. The DOIT method is unique because a discrete
ordinate iterative method is used to solve the scattering problem in a
spherical atmosphere. Although the DOIT module is implemented for 1D and 3D
atmospheres, it is strongly recommended to use it only for 1D, because the
Monte Carlo module (Chapter \ref{sec:scattering:mc}) is much more
appropriate for 3D calulations. More appropriate in the sense that it is much
more efficient. A literature review about scattering models for the microwave
region, which is presented in \citet{sreerekha04:_devel_rt_ghz_wp1}, shows that
former implementations of discrete ordinate schemes are only applicable for
(1D-)plane-parallel or 3D-cartesian atmospheres. All of these algorithms can
not be used for the simulation of limb radiances. A description of the DOIT
method is given in \theory, Chapter \ref{T-sec:doit} and has been published in
\citet{emde04:_doit_jgr} and in \citet{emde05:_phdthesis}.

The workspace methods required for DOIT calculations are implemented
in the files \fileindex{m\_scatrte.cc}, \fileindex{m\_cloudbox.cc} and
\fileindex{m\_optproperties.cc}. Here we introduce the steps to be performed in
a DOIT calculation along with the relevant workspace methods by means of a
controlfile example.

\section{The 1D control file example}
This example demonstrates how to set up a 1D DOIT calculation. A full
running controlfile example for a DOIT calculation can be found in the
ARTS package in the \fileindex{controlfiles/artscomponents/doit} directory. The file
is called \shortcode{TestDOIT.arts}. For detailed descriptions of the
workspace methods and variables please refer also to the online documentation
(start doc-server on your own computer with \shortcode{arts -s} or use
\artsurl{\docserver}).
 
\section{DOIT frame}
\label{sec:scattering:frame}

The first step for a DOIT calculation is the initialization of variables
required for a DOIT calculation using \wsmindex{DoitInit}:
\begin{code}
DoitInit
\end{code}
As the next step we have to calculate the
incoming field on the boundary of the cloudbox. This is done using the
workspace method \wsmindex{DoitGetIncoming}:
\begin{code}
DoitGetIncoming
\end{code}
The method \wsmindex{doit\_i\_fieldSetClearsky} interpolates the
incoming radiation field on all points inside the cloudbox to obtain
the initial field (\wsvindex{doit\_i\_field}) for the DOIT calculation. 
As a test one can alternatively start with a constant radiation field
using the method \wsmindex{doit\_i\_fieldSetConst}. 
\begin{code}
doit_i_fieldSetClearsky
\end{code}

The grid discretization plays a very significant role in discrete
ordinate methods. In spherical geometry the zenith angular grid is of
particular importance (cf. \theory, Section \ref{T-sec:doit:grid_opt_interp}). 
The angular discretization is defined in the workspace method
\wsmindex{DOAngularGridsSet}:
\begin{code} 
DOAngularGridsSet( doit_za_grid_size,
                   scat_aa_grid, scat_za_grid,
                   19, 10, "doit_za_grid_opt.xml" )
\end{code}
For down-looking geometries it is sufficient to define the generic inputs:
\begin{tabular}{ll}
  \shortcode{N\_za\_grid}& Number of grid points in zenith angle grid,
  recommended value: 19\\
  \shortcode{N\_aa\_grid}& Number of grid points in azimuth angle grid,
  recommended value: 37\\ 
\end{tabular}
From these numbers equally spaced grids are created and stored in the
work space variables \wsvindex{scat\_za\_grid} and
\wsvindex{scat\_aa\_grid}. 

For limb simulations it is important to use an optimized zenith angle
grid with a very fine resolution about 90\degree\ for the RT calculations.
Such a grid can be generated using the workspace method
\wsmindex{doit\_za\_grid\_optCalc}. Please refer to the online
documentation of this method. 
The filename of the optimized zenith angle grid can be given
as a generic input. If a filename is given, the equidistant grid is
taken for the calculation of the scattering integrals and the
optimized grid is taken for the radiative transfer part.
Otherwise, if no filename is specified
(\shortcode{za\_grid\_opt\_file = ""}) the equidistant grid is
taken for the calculation of the scattering integrals and for
the radiative transfer calculations. This option makes sense for
down-looking cases to speed up the calculation.

The main agenda for a DOIT calculation is
\wsaindex{doit\_mono\_agenda}. 
The agenda is executed by the workspace method
\wsmindex{DoitCalc}:
\begin{code}
DoitCalc
\end{code}

\subsection{The DOIT main agenda}
\label{sec:scattering:doit_main_agenda}

Although there are alternatives, the most elegant usage of DOIT involves
specifying it within the agenda \wsaindex{doit\_mono\_agenda}, which calculates
the radiation field inside the cloudbox.
The agenda requires the incoming
clearsky field on the cloudbox boundary as an input and gives as
output the scattered field on the cloudbox boundary if the sensor is
placed outside the cloudbox or the full scattered field in the
cloudbox if the sensor is placed inside the cloudbox.
\begin{code}
AgendaSet( doit_mono_agenda ){
 # Prepare scattering data for DOIT calculation (Optimized method):
   DoitScatteringDataPrepare
 # Alternative method (needs less memory):
 # scat_data_monoCalc
 # Perform iterations: 1. scattering integral. 2. RT calculations with 
 #   fixed scattering integral field, 3. convergence test 
   doit_i_field_monoIterate
 }		
\end{code}
The first method \wsmindex{DoitScatteringDataPrepare} prepares the single
scattering data for use in a DOIT scattering calculation. Namely, it
interpolates the data on the requested
frequency and performs the transformation from the scattering
frame into the laboratory frame. Alternatively the method
\wsmindex{scat\_data\_monoCalc} can be used. In this case only the
frequency interpolation is done and the transformations are done
later. The advantage is that this method needs less memory. For 1D
calculation it is recommended to use
\wsmindex{DoitScatteringDataPrepare} because it is much more
efficient. 

The iteration is performed in the method
\wsmindex{doit\_i\_field\_monoIterate}, which includes
\begin{itemize}
\item the calculation of the scattering integral field
\wsvindex{doit\_scat\_field} (requires agendas \wsaindex{pha\_mat\_spt\_agenda}
and \wsaindex{doit\_scat\_field\_agenda}),
\item the radiative transfer calculations in the cloudbox with fixed scattering
integral (requires agendas \wsaindex{spt\_calc\_agenda},
\wsaindex{ppath\_step\_agenda}, and \wsaindex{doit\_rte\_agenda}), and
\item the convergence test (requires agenda
\wsaindex{doit\_conv\_test\_agenda}).
\end{itemize}
For details on the agendas involved see Section
\ref{sec:scattering:doit_iterate_agendas}.

\subsection{Agendas used in \wsmindex{doit\_i\_field\_monoIterate}}
\label{sec:scattering:doit_iterate_agendas}
There are several methods which can be used in
\wsmindex{doit\_i\_field\_monoIterate}, for instance for the calculation of
the scattering integral. The methods are selected in the control-file by
defining several agendas.

\subsubsection{Calculation of the scattering integral:}

To calculate the scattering integral (\theory, Equation \ref{T-eq:scattering:scat-int})
the phase matrix (\wsvindex{pha\_mat}) is required. How 
the phase matrix is calculated is defined in the agenda
\wsaindex{pha\_mat\_spt\_agenda}:
\begin{code}
# Calculation of the phase matrix
AgendaSet( pha_mat_spt_agenda ){
 # Optimized option:
   pha_mat_sptFromDataDOITOpt
 # Alternative option:
 #  pha_mat_sptFromMonoData
}
\end{code}
If in \wsaindex{doit\_mono\_agenda}
the optimized method \wsmindex{DoitScatteringDataPrepare} is used we
have to use here the corresponding method
\wsmindex{pha\_mat\_sptFromDataDOITOpt}. Otherwise we have to use
\wsmindex{pha\_mat\_sptFromMonoData}.

To do the integration itself, we have to define
\wsaindex{doit\_scat\_field\_agenda}:
\begin{code}
AgendaSet( doit_scat_field_agenda ){
 doit_scat_fieldCalcLimb
 # Alternative: 
 # doit_scat_fieldCalc
}
\end{code}
Here we have two options. One is \wsmindex{doit\_scat\_fieldCalcLimb},
which should be used for limb simulations, for which we need a fine
zenith angle grid resolution to represent the radiation field. This
method has to be used if a zenith angle grid file is given in
\wsmindex{DOAngularGridsSet}. The
scattering integral can be calculated on a coarser grid resolution, hence
in \wsmindex{doit\_scat\_fieldCalcLimb}, the radiation field is
interpolated on the equidistant angular grids specified in
\wsmindex{DOAngularGridsSet} by the generic inputs \shortcode{Nza} and
\shortcode{Naa}. Alternatively, one can use 
\wsmindex{doit\_scat\_fieldCalc}, where this interpolation is not 
performed. This function is efficient for simulations in up- or
down-looking geometry, where a fine zenith angle grid resolution around
90\degree is not needed.

\subsubsection{Radiative transfer with fixed scattering integral term:}

With a fixed scattering integral field the radiative transfer equation
can be solved (\theory, Equation \ref{T-eq:scattering:vrte-fs-av}). The workspace
method to be used for this calculation is defined in
\wsaindex{doit\_rte\_agenda}. The most efficient and recommended
workspace method is \wsmindex{doit\_i\_fieldUpdateSeq1D} where the
sequential update which is described in \theory,
Section \ref{T-sec:doit:sequential_update} is applied. The workspace
method \wsmindex{doit\_i\_fieldUpdate1D} does the same calculation
without 
sequential update and is therefore much less efficient because the
number of iterations depends in this case on the number of pressure
levels in the cloudbox.  Other options
are to use a plane-parallel approximation implemented in the workspace
method \wsmindex{doit\_i\_fieldUpdateSeq1DPP}. This method is not
much more efficient than \wsmindex{doit\_i\_fieldUpdateSeq1D},
therefore it is usually better to use
\wsmindex{doit\_i\_fieldUpdateSeq1D} since it is more accurate. 
\begin{code}
AgendaSet( doit_rte_agenda ){
     doit_i_fieldUpdateSeq1D
   # Alternatives:
   # doit_i_fieldUpdateSeq1DPP
   # i_fieldUpdate1D
   }
\end{code}
The optical properties of the particles, i.e., extinction matrix and
absorption vector (for all scattering elements) are required for solving the
radiative transfer equation. How they are calculated is specified in
\wsaindex{spt\_calc\_agenda}. The workspace method
\wsmindex{opt\_prop\_sptFromMonoData} requires that the raw data is
already interpolated on the frequency of the monochromatic
calculation. This requirement is fulfilled when
\wsmindex{DoitScatteringDataPrepare} or \wsmindex{scat\_data\_monoCalc}
is executed before \wsmindex{doit\_i\_field\_monoIterate} (see
Section \ref{sec:scattering:doit_main_agenda}).
\begin{code}
AgendaSet( spt_calc_agenda ){
   opt_prop_sptFromMonoData
}
\end{code}
The work space method
\wsmindex{opt\_prop\_bulkCalc} is then used internally to derived the
bulk absorption vector \wsvindex{abs\_vec} and extinction
matrix \wsvindex{ext\_mat} from the workspace variable
\wsvindex{ext\_mat\_spt}. The gas absorption is added internally.

\subsubsection{Convergence test:}
After the radiative transfer calculations with a fixed scattering
integral field are complete the newly obtained radiation field is
compared to the old radiation field by a convergence test. The
functions and parameters for the convergence test are defined in the
agenda \wsaindex{doit\_conv\_test\_agenda}. There are several
options. The workspace methods \wsmindex{doit\_conv\_flagAbsBT} and
\wsmindex{doit\_conv\_flagAbs} compare the absolute differences of the
radiation field element-wise as described in
\theory, Equation \ref{T-eq:scattering:conv_test}. The convergence limits are specified by the
generic input \shortcode{epsilon} which specifies the convergence limit. A limit
must be given for each Stokes component. In
\wsmindex{doit\_conv\_flagAbsBT} the limits must be specified in
Rayleigh Jeans brightness temperatures whereas in
\wsmindex{doit\_conv\_flagAbs} they must be defined in the basic radiance unit
([W/(m$^2$Hz sr)]. Another option is to perform a least square
convergence test using the workspace method
\wsmindex{doit\_conv\_flagLsq}. Test calculations have shown that this
test is not safe, therefore the least square convergence test should
only be used for test purposes.
\begin{code}
AgendaSet( doit_conv_test_agenda) {
    doit_conv_flagAbsBT( doit_conv_flag, doit_iteration_counter,
                         doit_i_field, doit_i_field_old,
                         f_grid, f_index,
                         [0.1, 0.01, 0.01, 0.01] )
    # Alternative: Give limits in radiances
    #   doit_conv_flagAbs( doit_conv_flag, doit_iteration_counter,
    #                      doit_i_field, doit_i_field_old,
    #                      [0.1e-15, 0.1e-18, 0.1e-18, 0.1e-18] )
    #
    # If you want to look at several iteration fields, for example 
    # to investigate the convergence behavior, you can use
    # the following workspace method:
    # DoitWriteIterationFields( doit_iteration_counter, doit_i_field,
    #                           [2, 4] )
}
\end{code}
\subsection{Propagation of the DOIT result towards the sensor}
In order to propagate the result of the scattering calculation towards
the sensor, the fields needs to be interpolated on the direction of
the sensor's line of sight. This is done in the workspace method
\wsmindex{iyInterpCloudboxField}, which has to be put into the agenda
\wsaindex{iy\_cloudbox\_agenda}:

\begin{code}
AgendaSet( iy_cloudbox_agenda ){
    iyInterpCloudboxField
}
\end{code}

\section{3D DOIT calculations}

The DOIT method is implemented for 1D and 3D spherical atmospheres, but it is
strongly recommended to use it only for 1D calculations, because there
are several numerical difficulties related to the grid
discretizations. It is difficult to find appropriate discretizations
to get sufficiently accurate results in reasonable computation
time. Therefore only experienced ARTS users should use DOIT for 3D
calculations only for smaller cloud scenarios. Please refer to the online
documentation for the workspace method for  3D scattering calculations
(\wsmindex{doit\_i\_fieldUpdateSeq3D}). All other workspace
methods adapt automatically to the atmospheric dimensionality.





\chapter{Scattering calculations -- The Monte Carlo scattering module}
 \label{sec:scattering:mc}

\starthistory
 140924 & Started by Patrick Eriksson\\ 
\stophistory

(This is just a stub for this chapter. So far just rescuing some text from an
old (and now deleted) Wiki page.)\\

\noindent
The ARTS \textindex{Monte Carlo scattering module} offers an efficient method
for single viewing direction radiative transfer calculations in arbitrarily
complex 3D cloudy cases. When simulating the observation of detailed 3D cloudy
scenarios, reversed Monte Carlo algorithms have several advantages over other
methods, such as discrete ordinate and Forward Monte Carlo methods. These
features include:

\begin{itemize}
\item All computational effort is dedicated to calculating the Stokes vector at
  the location of interest and in the direction of interest. This is
  particularly relevant for space-borne remote sensing, where we are only
  interested in a narrow field of view. This is contrast to DOM methods where
  the whole radiation field is calculated.
\item CPU and memory cost scales more slowly than other methods with grid size.
  so that large or detailed 3D scenarios are not a problem. This stems from the
  suitability of Monte Carlo Integration (MCI) for evaluating integrals over
  highly dimensioned spaces. As well as CPU cost increasing dramatically in 3D
  DOM applications with the number of grid-points in each dimension, the memory
  requirements becomes prohibitive at moderate grid sizes due to the
  requirement that the radiance in every direction must be stored at each grid
  point.
\item Optically thick media are no problem. A feature of reversed Monte Carlo
  algorithms is that only parts of the atmosphere that actually contribute to
  the observed radiance are considered in the computation. So where the medium
  is optically thick due to absorption or scattering, only the parts of the
  atmosphere closest to the sensor are visited by the algorithm. This contrasts
  with DOM methods, where, as mentioned above, the whole radiation field is
  computed. Also, a requirement of DOM methods is that the optical thickness
  between adjacent grid points must be <= 1, which increases the grid-size, and
  hence cost, of the method.
\end{itemize}

For the theory behind the ARTS MC module see Section \ref{T-sec:montecarlo} of
\theory.





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