\chapter{Overview of clear-sky radiative transfer 
calculations}
 \label{sec:rte}


 \starthistory
 110611 & Extended and general revision (Patrick Eriksson).\\
 050613 & First complete version by Patrick Eriksson.\\
 \stophistory

\graphicspath{{Figs/rte/}}

This section gives an overview of the variables and the approach used to handle
the actual radiative transfer calculations. This includes an overview of how
effects caused by the sensor and surface are incorporated. The default
assumption is that there is no particle scattering and the (atmospheric)
absorption/emission is unpolarised. This is for simplicity denoted as clear-sky
calculations. Cases involving scattering or polarised absorption are handled by
the ``cloud box'' (Sec.~\ref{sec:fm_defs:cloudbox}). In short, the more
demanding calculations are restricted to a smaller domain of the model
atmosphere, and the radiative transfer in this domain is treated by dedicated
workspace methods.

The section deals with general aspects of radiative transfer and the
algorithms applied outside the cloud box. Even though the atmospheric absorption
and emission outside the cloud box are unpolarised, the expressions to apply
must allow that polarisation signals from the surface and the cloud box are
correctly propagated to the sensor.




\section{Stokes dimensionality}
%==============================================================================
\label{sec:fm_defs:polarisation}

To full polarisation state of radiation can be described by the Stokes
vector. The vector can be defined in different ways, but it has always
four elements. The Stokes vector, \StoVec, is here written as
\begin{equation}
  \label{eq:fm_defs:stokevec}
  \StoVec = \left[
  \begin{array}{c}
   \StoI\\ \StoQ \\ \StoU\\ \StoV
  \end{array}
  \right],
\end{equation}
where the first component (\StoI) is the full intensity of the
radiation, the second component (\StoQ) is the difference between
vertical and horizontal polarisation, the third component (\StoU) is the
difference for $\pm$45$^\circ$ polarisation and the last component
(\StoV) is the difference between left and right circular polarisation.
That is:
\begin{eqnarray}
  \StoI &=&   \Iv + \Ih = \Ipff + \Imff = \Ilhc + \Irhc, \\
  \StoQ &=&   \Iv - \Ih,                                 \\
  \StoU &=&   \Ipff - \Imff,                             \\
  \StoV &=&   \Ilhc - \Irhc,                             
\end{eqnarray}
where \Iv, \Ih, \Ipff, and \Imff\ are the intensity of the component linearly
polarised at the vertical, horizontal, +45\degree\ and -45\degree\ direction,
respectively, and \Irhc, and \Ilhc\ are the intensity for the right- and
left-hand circular components. Further details on polarisation and the
definition of the Stokes vector are found in \theory,
Section~\ref{T-sec:polarization}.

ARTS is a fully polarised forward model, but can be run with a smaller
number of Stokes components. The selection is made with the workspace
variable \wsvindex{stokes\_dim}. For example, gaseous absorption and
emission are in general unpolarised, and if not scattering has to be
considered it is sufficient to only include the first Stokes
components in the simulations. To include higher order Stokes
components results only, in this case, in slower calculations. The
general case is here denoted as \textindex{vector radiative transfer},
while \textindex{scalar radiative transfer} refers to the case when
only the first Stokes component is considered.
 



\section{Overall calculation procedure}
%===================
\label{sec:fm_defs:calcproc}

The structure of the part handling complete radiative transfer calculations is
fixed, where the main workspace method is denoted as \wsmindex{yCalc}. That is,
most ARTS control files include a call of \builtindoc{yCalc} and this section
outlines this method and the associated main variables.

The calculation approach fits with the formalism presented in Sections
\ref{T-sec:formalism:fm}-\ref{T-sec:formalism:sensor} of \theory, where the
separation between atmospheric radiative transfer and inclusion of sensor
effects shall be noted especially, and a similar nomenclature is used here:
\begin{description}
\item[\MsrVct]: Complete measurement vector. In addition to atmospheric
  radiative transfer, the vector can include effects by sensor characteristics
  and data reduction operations. The corresponding workspace variable is 
  \wsvindex{y}.
\item[\aMpiVct{b}]: Monochromatic pencil beam data for a measurement block. The
  definition of a measurement block is found in
  Section~\ref{sec:fm_defs:seqsandblocks}. This vector is only affected by
  atmospheric radiative transfer. Only used internally and there is no
  corresponding workspace variable.
\item[\aMpiVct{y}]: Monochromatic data for a single (pencil beam)
  line-of-sight. \aMpiVct{b} consists of one or several \aMpiVct{y} appended.
  The corresponding workspace variable is \wsvindex{iy}.
\item[\aSnsMtr{b}]: The complete sensor response matrix, for a measurement
  block. Can include data reduction. The corresponding workspace variable is
  \wsvindex{sensor\_response}.
\end{description}
\begin{algorithm}[t]
 \begin{algorithmic}
  \STATE{allocate memory for the matrix $\MsrVct$}
  \COMMENT{Equation \ref{eq:fm_defs:measseq}}
  \STATE{allocate memory for the matrix \aMpiVct{b}}
  \COMMENT{Equation \ref{eq:fm_defs:freqs_of_ib}}
  \FORALL{sensor positions}
   \FORALL[Section \ref{sec:fm_defs:seqsandblocks}]
                                    {pencil beam directions of the block}
    \STATE{call \builtindoc{iy\_clearsky\_agenda}, giving \aMpiVct{y}}
    \COMMENT{Algorithm \ref{alg:fm_defs:iyCSagenda}}
    \STATE{unit conversion of \aMpiVct{y} following \wsvindex{y\_unit}}
    \COMMENT{Section \ref{sec:fm_defs:unit}}
    \STATE{copy \aMpiVct{y} to correct part of \aMpiVct{b}}
   \ENDFOR
   \STATE{put the product \aSnsMtr{b}\aMpiVct{b} in correct part of 
          $\MsrVct$}
  \ENDFOR
 \end{algorithmic}
 \caption{Outline of the overall clear sky radiative transfer calculations
   (\builtindoc{yCalc}).}
 \label{alg:fm_defs:yCalc}
\end{algorithm}
\begin{algorithm}[t]
 \begin{algorithmic}
   \STATE{determine the propagation path by \builtindoc{ppath\_calc}}
   \COMMENT{Section \ref{sec:fm_defs:ppaths}}
   \STATE{determine the radiation at the start of the propagation path}
   \COMMENT{Section \ref{sec:fm_defs:rad_bkgr}}
   \STATE{perform radiative transfer along the propagation path}
   \COMMENT{Section \ref{sec:fm_defs:rte}}
 \end{algorithmic}
 \caption{The main operations for methods to be part of
   \wsvindex{iy\_clearsky\_agenda}. The same applies to methods for
   \wsvindex{iy\_clearsky\_basic\_agenda}.}
 \label{alg:fm_defs:iyCSagenda}
\end{algorithm}
The \builtindoc{yCalc} method is outlined in Algorithm~\ref{alg:fm_defs:yCalc}.
For further details of each calculation step, see the indicated equation or
section. In summary, \builtindoc{yCalc} appends data from different pencil beam
calculations and applies the sensor response matrix (\aSnsMtr{b}). The actual
radiative transfer calculations are not part of \builtindoc{yCalc}.

Atmospheric radiative transfer is solved for each pencil beam direction
(line-of-sight) separately. It is the task of \wsvindex{iy\_clearsky\_agenda}
(Algorithm~\ref{alg:fm_defs:iyCSagenda}) to perform a single such clear sky
radiative transfer calculation. This agenda, in its turn, makes us of other
agendas, such as \builtindoc{ppath\_step\_agenda}. All methods developed for
\builtindoc{iy\_clearsky\_agenda} adapt automatically to the value of
\wsvindex{stokes\_dim}.

That is, \builtindoc{yCalc} is a common method, independent of the details of
the radiative transfer. For example, \builtindoc{yCalc} is used both if emission
measurements or pure transmission data are simulated, that choice is made
inside \builtindoc{iy\_clearsky\_agenda} (see further Section
\ref{sec:fm_defs:rte}). 

The three following sections describes the main calculation steps of 
\builtindoc{iy\_clearsky\_agenda}, in the order they are executed.


\section{Propagation paths}
%===================
\label{sec:fm_defs:ppaths}

A pencil beam path through the atmosphere to reach a position along a specific
line-of-sight is denoted as the \textindex{propagation path}. Propagation paths
are described by a set of points on the path, and the distance along the path
between the points. These quantities, and a number of auxiliary variables, are
stored together in a structure described in Section~\ref{sec:ppath:Ppath}. The
path points are primarily placed at the crossings of the path with the
atmospheric grids (\builtindoc{p\_grid}, \builtindoc{lat\_grid} and
\builtindoc{lon\_grid}). A path point is also placed at the sensor if it is
placed inside the atmosphere. Points of surface reflections and tangent points
are also included if such exist. More points can also be added to the
propagation path, for example, by setting an upper limit for the distance along
the path between the points. This is achieved by the variable
\builtindoc{ppath\_lmax}, see further Sections \ref{sec:fm_defs:accuracy} and
\ref{sec:ppath:usage}.

\begin{figure}
 \begin{center}
  \includegraphics*[width=0.95\hsize]{ppath_cases2}
  \caption{Examples on allowed propagation paths for a 2D atmosphere. 
    The atmosphere is plotted as in Figure~\ref{fig:fm_defs:2d} beside
    that the points for the atmospheric fields are not emphasised.
    The position of the sensor is indicated by an asterisk $(\ast)$,
    the points defining the paths are plotted as circles $(\circ)$,
    joined by a solid line. The part of the path outside the
    atmosphere, not included in the path structure, is shown by a
    dashed line. Path points corresponding to a tangent point are
    marked by an extra plus sign $(\oplus)$. The shown paths include
    the minimum set of definition points. There exists also the
    possibility to add points inside the grid cells, for example, to
    ensure that the distance between the path points does not exceed
    a specified limit.}
  \label{fig:fm_defs:ppath_cases2}
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_ppath_cases.

\begin{figure}
 \begin{center}
  \includegraphics*[width=0.95\hsize]{ppath_cases1}
  \caption{Examples on allowed propagation paths for a 1D atmosphere
    with an activated cloud box. Plotting symbols as in
    Figure~\ref{fig:fm_defs:ppath_cases2}. When the sensor is placed 
    inside the cloud box, the path is defined with a single point, 
    to know for which position and line-of-sight the intensity field of
    the cloud box shall be interpolated. }
  \label{fig:fm_defs:ppath_cases1}
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_ppath_cases.

\begin{figure}
 \begin{center}
  \includegraphics*[width=0.95\hsize]{ppath_badcases}
  \caption{Examples on \emph{not} allowed propagation paths for a 2D 
    atmosphere. The constraints for allowed paths are discussed in the
    text.}
  \label{fig:fm_defs:ppath_badcases}
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_ppath_cases.


The propagation paths are determined basically by starting at the
sensor and following the path backwards by some \textindex{ray
  tracing} technique. If the sensor is placed above the model
atmosphere, geometrical calculations are used (as there is no
refraction in space) to find the crossing between the path and the top
of the atmosphere where the ray tracing then starts. Paths are tracked
backwards until the top of the atmosphere is reached, or there is an
intersection with the cloud box or the surface. The propagation path
(or paths) before a surface reflection is calculated when determining
the up-welling radiation from the surface
(Section~\ref{sec:fm_defs:surface}). Example on propagation
paths are shown in Figures~\ref{fig:fm_defs:ppath_cases2} and 
\ref{fig:fm_defs:ppath_cases1}.
 
Not all propagation paths are allowed for 2D and 3D. The paths can
only enter and leave the model atmosphere at the top of the
atmosphere, as the atmospheric fields are treated to be undefined
outside the covered latitude and longitude ranges
(Figure~\ref{fig:fm_defs:ppath_badcases}). In addition, if the sensor
is placed outside the model atmosphere, the line-of-sight zenith angle
must be $\geq90\degree$, and the tangent point position of the
propagation paths must be inside the end points of the latitude and
longitude grids, but can be above the top of the atmosphere. Hence, it
is allowed that the propagation path is totally outside the
atmosphere, as long as the viewing direction is downward and the
lowest point of the path, the tangent point, is inside the latitude
and longitude limits of the model atmosphere.

Propagation paths can be calculated separately by the method
\wsmindex{ppathCalc}. However, for standard calculations the propagation paths
are calculated internally by \builtindoc{yCalc}. The calculation of the path
from one crossing of the grids to next crossing is defined by
\wsaindex{ppath\_step\_agenda}. Depending on which function that is selected
for \builtindoc{ppath\_step\_agenda}, refraction will be considered etc.
Available workspace methods are presented in Section~\ref{sec:ppath:usage}.




\section{The radiative background}
%===================
\label{sec:fm_defs:rad_bkgr}

\subsection{Overview}
%===================
%
The radiative intensity at the starting point of the path, and in the
direction of the line-of-sight at that point, is denoted as the
\textindex{radiative background}. Four possible radiative backgrounds
exist:
\begin{description}
\item[Space] When the propagation path starts at the top of the
  atmosphere, space is the radiative background. The normal case
  should be to set the radiation at the top of the atmosphere to be
  cosmic background radiation. An exception is when the sensor is
  directed towards the sun. The radiative background at the top of the
  atmosphere is determined by \wsaindex{iy\_space\_agenda}. If a
  propagation path is totally outside the model atmosphere, the
  observed monochromatic pencil beam intensity (\aMpiVct{y}\ in
  Algorithm~\ref{alg:fm_defs:yCalc}) equals the output of
  \builtindoc{iy\_space\_agenda}.
\item[The surface] The sum of surface emission and radiation reflected by the
  surface is the radiative background when the propagation path intersects with
  the surface. The calculation of the up-welling radiation from the surface is
  treated by a dedicated section below (Section~\ref{sec:fm_defs:surface}.)
\item[Surface of cloud box] For cases when the propagation path enters
  the cloud box the radiative background is the intensities leaving
  the cloud box. This radiation is obtained by
  \wsaindex{iy\_cloudbox\_agenda}. 
\item[Interior of cloud box] If the sensor is situated inside the
  cloud box, there is basically no propagation path. The radiative
  background, and also the final spectrum, equals the internal
  intensity field of the cloud box at the position of the sensor, in
  the direction of the sensor line-of-sight. This case is also handled
  by \builtindoc{iy\_cloudbox\_agenda}.
\end{description}
It should be noted that except for the first case above, the determination of
the radiative background involves further radiative transfer calculations. For
example, in the case of surface reflection, the down-welling radiation is
determined by a new call of \builtindoc{iy\_clearsky\_agenda} and the radiative
background for that calculation is then space or the cloud box. The intensity
field entering the cloud box is calculated by calls of
\builtindoc{iy\_clearsky\_basic\_agenda} (with cloud box deactivated) and the
radiative background is then space or the surface. This results in that space
is normally the ultimate radiative background for the calculations. The
exception is for propagation paths that intersects with the surface, and the
surface is treated to act as a blackbody. For such cases, the propagation path
effectively starts at the surface.


\subsection{Surface scattering and emission}
%===================
\label{sec:fm_defs:surface}

If there is an interception of the propagation path by the surface,
emission and scattering by the surface must be considered. The overall
treatment of these effects is fixed. The upwelling radiation from the surface
can be written as (Figure~\ref{fig:fm_defs:surface_refl})
\begin{equation}
  \MpiVct_s^u = \MpiVct_e + \sum_l^{} \mathbf{R}_l \MpiVct_l^d
  \label{eq:fm_defs:surfacerefl}
\end{equation}
where \MpiVct\ is the Stokes vector for one frequency, $\MpiVct_s^u$
is the total upward travelling intensity from the surface along the
propagation path, $\MpiVct_e$ is the emission from the surface,
$\MpiVct_l^d$ is the downward travelling intensity reaching the
surface along direction $l$, and $\mathbf{R}_l$ is the reflection
coefficient matrix from direction $l$ to the present propagation path.
The emission from the surface $(\MpiVct_e)$ is stored in
\wsvindex{surface\_emission}, the directions $l$ for which downward
travelling intensities are given by \wsvindex{surface\_los}, and the
reflection coefficients $(\mathbf{R})$ are stored in
\wsvindex{surface\_rmatrix}. These workspace variables are handled by
\wsaindex{surface\_prop\_agenda}. Surface reflections and emission are
discussed further in Section~\ref{sec:surf:eandr}.
\begin{figure}
 \begin{center}
  \includegraphics*[width=0.95\hsize]{ground_refl}
  \caption{Schematic of Equation \ref{eq:fm_defs:surfacerefl}.}
  \label{fig:fm_defs:surface_refl}
 \end{center}
\end{figure}



\section{Basic radiative transfer variables and expressions}
%---
\label{sec:fm_defs:rte}

This section describes how the core radiative transfer equation is solved
practically in ARTS. Focus is put on emission measurements as ARTS is
intended primarily for such observations. However, simulations of transmission
measurements are also possible.

Beside the actual radiances, \builtindoc{iy\_clearsky\_agenda} can provide
weighting functions and auxiliary data. These later variables are not supported
by all parts of ARTS, and for efficiency reasons there exists a simpler version
of the agenda. This version is denoted as
\wsvindex{iy\_clearsky\_basic\_agenda} and returns only radiances.


\subsection{Cases with emission}
%===================
\label{sec:fm_defs:emission}

The complete vector radiative transfer equation, including scattering, is given
by Equation \ref{T-eq:rtetheory:VRTE} in \theory. If scattering can be
neglected, the equation can be written as
\begin{equation}
  \label{eq:rte:vrte}
  \frac {\DiffD\StoVec}{\DiffD s} = -\ExtMat\StoVec + \AbsVec B,
\end{equation}
where \StoVec\ is the intensity vector (the Stokes vector), $s$ is the distance
along the propagation path, \ExtMat\ is the extinction matrix, \AbsVec\ is the
absorption vector and $B$ is the source function (a scalar). If local
thermodynamic equilibrium applies, $B$ equals the Planck function describing
blackbody radiation. See further \theory, Section~\ref{T-sec:rte_theory}. For
``clear-sky conditions'' the matrix \ExtMat\ is diagonal, with all diagonal
elements equal, and only the first of the elements of \AbsVec\ is non-zero.

The radiative transfer equation above can be solved in many ways, and with
different level of refinement. The standard approach in ARTS is to solve the
radiative transfer from one point of the propagation path to next. For the
first Stokes element the following expression is applied (compare \theory,
Equation~\ref{T-eq:rtetheory:layer})
\begin{equation}
  \label{eq:fm_defs:rte_step}
  I_{i+1} = I_ie^{-\aOth{i}} + \bar{B}_i(1-e^{-\aOth{i}}),
\end{equation}
with
\begin{eqnarray}
  \bar{B}_i &=& (B(T_i)+B(T_{i+1}))/2, \\
  \aOth{i} &=& \Delta s_i(k_i+k_{i+1})/2, 
  \label{eq:taustep}
\end{eqnarray}
where $I_i$, $T_i$ and $k_i$ are the radiance, temperature and absorption
coefficient, respectively, at point $i$ of the propagation path, and $\Delta
s_i$ is the distance along the path between point $i$ and $i+1$. That is,
$\bar{B}_i$ is a Planck function average of the path step, and the absorption
is assumed to vary linearly between the two points. The start value of \Mpi\ is
governed by the radiative background (Section \ref{sec:fm_defs:rad_bkgr}).

As mentioned above, the emission is unpolarised for the conditions assumed
here, and the emission term vanishes for higher Stokes elements. Accordingly,
the expression for the second Stokes component is
\begin{equation}
  Q_{i+1}(\Frq) = Q_i(\Frq)e^{-\aOth{i}}.
  \label{eq:fm_defs:rte_step2}
\end{equation}
The third and forth Stokes component are handled likewise.

An alternative way to perform the calculations for the first Stokes element
would be
\begin{equation}
  \label{eq:fm_defs:rte_alt}
  I = \sum_i \mathcal{T}_{i+1} \bar{B}_i(1-e^{-\aOth{i}}),
\end{equation}
where $I$ is the final intensity and $\mathcal{T}_i$ is the transmission
between the sensor and point $i$. This calculation approach is not used as it
fits poorer with the calculation of weighting functions
(Section~\ref{sec:wfuns}). However, the calculation of weighting functions is
simplified if $\mathcal{T}_i$ is at hand, and this quantity is also tracked by
\builtindoc{iyEmissionStandardClearsky}. The value of $\mathcal{T}$ at the
radiative background can also be outputted as auxiliary data
(\wsvindex{iy\_aux}).

The expressions above are implemented in the workspace method
\wsmindex{iyEmissionStandardClearsky}, intended to be part of
\builtindoc{iy\_clearsky\_agenda}. For inclusion in
\builtindoc{iy\_clearsky\_basic\_agenda}, there is a version denoted as
\wsmindex{iyEmissionStandardClearskyBasic}. This version is not calculating
$\mathcal{T}_i$, or provides any auxiliary data. 

The term $B$ is set by \wsvindex{emission\_agenda}. The setting of this agenda
basically determines the unit of the final outcome of
Eq.~\ref{eq:fm_defs:rte_step}, see further Sec.~\ref{sec:fm_defs:unit}.
For radiance calculations, the standard workspace method to use inside 
\builtindoc{emission\_agenda} is \wsmindex{emissionPlanck}. The Planck function
is in this method, and in ARTS generally, defined as
\begin{equation}
  \label{eq:planck}
  \Planck(\Tmp) = \frac{2\planckCns\Frq^3}
                  {\speedoflight^2(exp(\planckCns\Frq/\boltzmannCns\Tmp)-1)}
\end{equation}
where \planckCns\ is the Planck constant, \speedoflight\ the speed of light and
\boltzmannCns\ the Boltzmann constant. This expression gives the total power,
per unit frequency per unit area per solid angle. (The Planck function can also
be defined as a function of wavelength.)

As long as cosmic background radiation is the only type of non-telluric
radiation that has to be considered, the standard method for inclusion in
\wsvindex{iy\_space\_agenda} is \wsmindex{MatrixCBR} (together with some calls
of \builtindoc{Ignore}). Please, note that \builtindoc{emission\_agenda} and 
\builtindoc{iy\_space\_agenda} must be defined in a consistent manner, that
they use the same unit for $B$.


\subsection{Pure transmission calculations}
%===================
\label{sec:fm_defs:transmission}

If only the attenuation of a signal is of concern (i.e.\ atmospheric and
surface emissions are neglected), Equation~\ref{eq:fm_defs:rte_step2} can be
applied for all Stokes elements. The initiation of the Stokes vector by the
radiative background must also be adopted. The standard choice for
\wsvindex{iy\_space\_agenda} should here be \wsmindex{MatrixUnitIntensity}
(instead of \builtindoc{MatrixCBR}). Otherwise the calculations are performed
exactly as for cases with emission (\builtindoc{emission\_agenda} can be left
undefined).

Calculations of pure transmission type are performed by the method denoted as
\wsmindex{iyBeerLambertStandardClearsky}. This method is complemented by
\wsmindex{iyBeerLambertStandardCloudbox}, to be applied inside the cloud box.


\section{Output unit}
%==============================================================================
\label{sec:fm_defs:unit}

There is no fixed unit for calculated spectra (\wsvindex{y}), it depends on the
calculation set-up. For example, if emission is considered, or if just
transmissions are calculated. 

The primary unit for emission data (radiances) is [W/(Hz$\cdot$m$^2\cdot$sr)],
The emission intensity corresponds directly with the definition of the Planck
function (Eq.~\ref{eq:planck}). This radiance can be converted to other
emission units by the workspace variable \wsvindex{y\_unit}. For comments on
the practical aspects, see the built-in documentation of \builtindoc{y\_unit}.
Otherwise, considerations and expressions for the unit conversion are mainly
discussed in ARTS-2 journal paper \citep[][Sec.~5.7]{eriksson:arts2:11}.



\section{Compulsory sensor and data reduction variables}
%==============================================================================
\label{sec:fm_defs:sensor1}

The instrument that detects the simulated radiation is denoted as the
sensor\index{sensor, the}. The forward model is constructed in such
way that a sensor must exist. For cases when only monochromatic pencil
beam radiation is of interest, the positions and directions for which
the radiation shall be calculated are given by specifying an imaginary
sensor with infinite frequency and angular resolution. The workspace
variables for the sensor that always must be specified are
\builtindoc{sensor\_response}, \builtindoc{sensor\_pos},
\builtindoc{sensor\_los}, \builtindoc{antenna\_dim},
\builtindoc{mblock\_za\_grid} and \builtindoc{mblock\_aa\_grid}. These
variables are presented separately
below. 


\subsection{Sensor position\index{sensor position}}
%===================
\label{sec:fm_defs:sensorpos}

The observation positions of the sensor are stored in
\wsvindex{sensor\_pos}. This is a matrix where each row corresponds to
a sensor position. The number of columns in the matrix equals the
atmospheric dimensionality (1 column for 1D etc.). The columns of the
matrix (from first to last) are radius, latitude and longitude.
Accordingly, row $i$ of \builtindoc{sensor\_pos} for a 3D case is
$(\aRds{i},\aLat{i},\aLon{i})$. The sensor position can be set to any
value, but the resulting propagation paths (also dependent on
\builtindoc{sensor\_los}) must be valid with respect to the model
atmosphere (see Section~\ref{sec:fm_defs:ppaths}). An obviously
incorrect choice is to place the senor below the surface altitude. If
the sensor is placed inside the model atmosphere, any sensor
line-of-sight is allowed, this including the cases that the sensor is
placed on the surface looking down, and that the sensor is placed
inside the cloud box.

The fact that the sensor position can be given any value implies that
the radius must be used in \builtindoc{sensor\_pos}, in contrast to
\builtindoc{z\_surface} and \builtindoc{z\_field} where the altitude
above the geoid is applied. This is the case as, for 2D and 3D, the
sensor can be placed outside the covered latitude and longitude
ranges, thus outside the defined geoid, and the geometrical altitude is
undefined. 

The sensor is treated to be motionless when calculating the spectrum,
or spectra, for each given observation position. One or several
spectra can be calculated for each position as described in
Section~\ref{sec:fm_defs:seqsandblocks}.


\subsection{Line-of-sight\index{line-of-sight}}
%===================
\label{sec:fm_defs:los}

The viewing direction of the sensor, the line-of-sight, is described
by two angles, the zenith angle (\ZntAng) and the azimuth angle
(\AzmAng). The zenith angle exists for all atmospheric
dimensionalities, while the azimuth angle is defined only for 3D.
The term line-of-sight is not only used in connection with the sensor,
it is also used to describe the local propagation direction along the
path taken by the observed radiation
(Section~\ref{sec:fm_defs:ppaths}).  The zenith and azimuthal angles
are defined in an identical way in both of these contexts (sensor
pointing direction; local propagation direction). This is expected as
the position of the sensor is the end point of the propagation path.
The sensor line-of-sight is the direction the antenna is pointed to
receive the radiation. The line-of-sight for propagation paths is
defined likewise, it is the direction in which a hypothetical sensor
must be placed to receive the radiation along the propagation path at
the point of interest. This means that the line-of-sight and the
photons are going in opposite directions. As a true sensor has a
finite spatial resolution (described by the antenna pattern),
theoretically there is an infinite number of line-of-sights associated
with the sensor, but in the forward model, spectra are only calculated
for a discrete set of directions. If a sensor line-of-sight is
mentioned without any comments, it refers to the direction in which
the centre of the antenna pattern is directed.

\begin{figure}
 \begin{center}
  \begin{minipage}[c]{0.6\textwidth}
   \includegraphics*[width=0.99\hsize]{za_and_aa_angles}
  \end{minipage}%
  \begin{minipage}[c]{0.4\textwidth}
   \caption{Definition of zenith angle, \ZntAng, and azimuth angle, 
       \AzmAng, for a line-of-sight. The figure shows a line-of-sight
       with a negative azimuth angle.}
   \label{fig:fm_defs:los}
  \end{minipage}
 \end{center}
\end{figure}           
 
The \textindex{zenith angle}, \ZntAng, is simply the angle between the
line-of-sight and the zenith direction (Figure~\ref{fig:fm_defs:los}).
The valid range for 1D and 3D cases is $[0,180\degree]$. In the case
of 2D, zenith angles down to -180\degree\ are also allowed, where the
distinction is that positive angles mean a viewing direction towards
higher latitudes, and negative angles mean a viewing direction towards
lower latitudes. It should be mentioned that the zenith and nadir
directions are here defined to be along the line passing the centre of
the coordinate system and the point of concern
(Section~\ref{sec:ppath:geoid}). A nadir observation,
$\ZntAng=180\degree$, is thus a measurement towards the centre of the
coordinate system.

The \textindex{azimuth angle}, \AzmAng, is given with respect to the
\textindex{meridian plane}.  That is, the plane going through the
north and south poles of the coordinate system $(\Lat=\pm90\degree)$
and the sensor. The valid range is $[-180\degree,180\degree]$ where
angles are counted clockwise; 0\degree means that the viewing or
propagation direction is north-wise and +90\degree\ means that the
direction of concern goes eastward. This definition does not work for
position on the poles. To cover these special cases, the definition is
extended to say that for positions on the poles the azimuth angle
equals the longitude along the viewing direction. For example, if
standing on any of the poles and the viewing direction is towards
Greenwich, the azimuth angle is 0\degree.

The sensor line-of-sights are stored in \wsvindex{sensor\_los}. This
workspace variable is a matrix, where the first column holds zenith
angles and the second column is azimuth angles. For 1D and 2D there is
only one column in the matrix, while for 3D a row $i$ of the matrix is
$(\aZntAng{i},\aAzmAng{i})$. The number of rows for
\builtindoc{sensor\_los} must be the same as for
\builtindoc{sensor\_pos}.


\subsection{Sensor characteristics and data reduction}
%===================
\label{sec:fm_defs:sensorchar}

The term ``sensor characteristics''\index{sensor characteristics} is
used here as a comprehensive term for the response of all sensor parts
that affect how the field of monochromatic pencil beam intensities are
translated to the recorded spectrum. For example, the antenna pattern,
the side-band filtering and response of the spectrometer channels are
normally the most important characteristics for a microwave heterodyne
radiometer. Any processing of the spectral data that takes place
before the retrieval is denoted as \textindex{data reduction}. The
most common processing is to represent the original spectra with a
smaller set of values, that is, a reduction of the data size. The most
common data reduction techniques is binning and Hotelling
transformation by an eigenvector expansion.

In ARTS, the influence of sensor characteristics and data reduction is
incorporated by transfer matrices\index{sensor transfer matrix}. The
application of these transfer matrices assumes that each step is a
linear operation, which should be the case for the response of the
parts of a well designed instrument. Non-linear data reduction could
be handled by special workspace methods.

The sensor and data reduction are described as a series of units, each
having its own transfer matrix.  There is only one compulsory transfer
matrix and it is \wsvindex{sensor\_response}. There are several workspace
variables associated with this transfer matrix where
\wsvindex{antenna\_dim}, \wsvindex{mblock\_za\_grid} and
\wsvindex{mblock\_aa\_grid} are the compulsory ones.

The variable \builtindoc{antenna\_dim} gives the dimensionality of the
antenna pattern\index{antenna pattern dimensionality}, where the
options are 1 and 2, standing for 1D and 2D, respectively. A 1D
antenna dimensionality means that the azimuth extension of the
antenna pattern is neglected, there is only a zenith angle variation
of the response. A 2D antenna pattern is converted to a 1D pattern by
integrating the azimuth response for each zenith angle. For cases
with 1D antenna patterns, \builtindoc{mblock\_aa\_grid} must be set to
be an empty vector.

For each sensor position, a number of monochromatic pencil beam
spectra are calculated. The monochromatic frequencies are given by
\builtindoc{f\_grid}. The pencil
beam directions are obtained by summing the sensor line-of-sight
angles (\builtindoc{sensor\_los}) for the position and the values of
\builtindoc{mblock\_za\_grid} and \builtindoc{mblock\_aa\_grid}. For
example, pencil beam zenith angle $i$ is calculated as
\begin{equation}
  \aZntAng{i} = \aZntAng{0} + \Delta\aZntAng{i}
  \label{eq:fm_defs:psi_grid}
\end{equation}
where \aZntAng{0} is the sensor line-of-sight for the position of
concern and $\Delta\aZntAng{i}$ is value $i$ of
\builtindoc{mblock\_za\_grid}.  With other words,
\builtindoc{mblock\_za\_grid} and \builtindoc{mblock\_aa\_grid} give
the grid (relative to the sensor line-of-sight) for the calculation of
the intensity field that will be weighted with the antenna response.


\subsection{Measurement sequences and blocks}
%===================
\label{sec:fm_defs:seqsandblocks}

The series of observations modelled by the simulations is denoted as
the \textindex{measurement sequence}. That is, a measurement sequence
covers all spectra recorded at all considered sensor positions. A
measurement sequence consists of one or several \textindex{measurement
  block}s. The observations inside the various blocks differ only with
an off-set of the line-of-sight, all other factors should be common
for all blocks. A block can be treated as a measurement cycle that is
repeated, an integer number of times, to form the measurement
sequence.  The measurement blocks correspond normally to each unique
sensor position of the sequence.

A measurement block covers one or several recorded spectra, depending
on the measurement conditions and the atmospheric dimensionality. A
block can consist of several spectra when there is no effective motion
of the sensor with respect to the atmospheric fields. It should be
noted that for 1D cases, a motion along a constant radius has no
influence on the simulated spectra as the same atmospheric fields are
seen for a given viewing direction. It is favourable, if possible, to
handle all spectra as a single block, instead of using a block for
each sensor position. This is the case as the antenna patterns for the
different line-of-sights are normally overlapping and a pencil beam
spectrum can be used in connection with several measurement spectra to
estimate the intensity field. If a measurement sequence is divided
into several blocks even if a single block would be sufficient, pencil
beam spectra for basically identical propagation paths can be
calculated several times, which of course will increase the
computational time. To summarise, for cases when the sensor is not in
motion, or with a 1D atmosphere and a sensor not moving vertically,
the aim should be to use a single block for the measurement sequence.

If not a single block is used, the standard option should be that the
blocks cover one spectrum each. There could exist reasons to select an
intermediate solution, to let the extent of the blocks be several
spectra (but not the full measurement sequence). This could be the
case when the atmospheric dimensionality is 2D or 3D, and the sensor
is moving but the movement during some subsequent spectra can be
neglected. If this can be done must be judged by comparing the
movement of the sensor during the extent of the considered block size
and the spatial resolution, in the direction of the movement, that is
hoped to be achieved. If this intermediate solution shall be an
option, the difference in zenith and azimuth angles between the
spectra must be the same for all blocks, otherwise
\builtindoc{sensor\_response} cannot be applied for all blocks as done
below in Equation~\ref{eq:fm_defs:measseq}.

For each block, pencil beam spectra are calculated for the
line-of-sights obtained when summing \wsvindex{sensor\_los} and
\wsvindex{mblock\_za\_grid} (and possibly
\wsvindex{mblock\_aa\_grid}), as described in
Section~\ref{sec:fm_defs:sensorchar}. The pencil beam spectra for each
line-of-sight are appended vertically to form a common vector,
\aMpiVct{b}. Values are put in following the order in
\builtindoc{f\_grid}. Hence, the frequencies for this vector are
\begin{equation}
  \aMpiVct{b} = 
  \left[ \begin{array}{c} 
     \left[
          \begin{array}{c} \aFrq{1}\\\vdots\\\aFrq{n} \end{array} 
     \right] \\
     \vdots \\
     \left[
          \begin{array}{c} \aFrq{1}\\\vdots\\\aFrq{n} \end{array} 
     \right]
     \end{array} \right]
  \label{eq:fm_defs:freqs_of_ib}
\end{equation}
where \aFrq{i} is element $i$ of \builtindoc{f\_grid} and $n$ the length of
the same vector. The order of the angles inside \builtindoc{mblock\_za\_grid}
and \builtindoc{mblock\_aa\_grid} is followed when looping the pencil beam
directions, where the azimuth angle direction is the innermost loop.
That is, for 2D antenna patterns all azimuth angles are looped for the
first zenith angle etc. 

The workspace variable \builtindoc{sensor\_response} is here denoted as
\aSnsMtr{b}. It is applied on each \aMpiVct{b} and the results are
appended vertically, following the order of the positions in
\builtindoc{sensor\_pos}
\begin{equation}
  \MsrVct = \left[ \begin{array}{c} \aSnsMtr{b}\aMpiVct{{b,1}} \\ 
                                    \aSnsMtr{b}\aMpiVct{{b,2}} \\
                                    \vdots                     \\
                                    \aSnsMtr{b}\aMpiVct{{b,n}} 
            \end{array} \right]
  \label{eq:fm_defs:measseq}
\end{equation}
where $1$ indicates the first sensor position etc. This equation shows
that \wsvindex{sensor\_response} shall contain at least a description
of the antenna response. The matrix \shortcode{Hb} can also cover
other sensor characteristics and data reduction if the features of
concern are common for all measurement blocks. 

As the sensor line-of-sight and block grid values are just added,
there is an ambiguity of the line-of-sight. It is possible to apply a
constant off-set to the line-of-sights, if the block grids are
corrected accordingly. For example, if the simulations deal with limb
sounding and a 1D atmosphere, where normally a single block should be
used despite a number of spectra are recorded, it could be practical
to set the line-of-sight to the viewing direction of the uppermost or
lowermost spectrum, and the zenith angles in \builtindoc{mblock\_za\_grid}
will not be centred around zero which is the case when the ``true''
line-of-sight is used.

It should be noted that the compulsory sensor variables give no
information about the content of the obtained \MsrVct, as it is not
clear which parts and features the block transfer matrix covers. If
\shortcode{Hb} only incorporates the antenna pattern, the result is a set
of hypothetical spectra corresponding to a point inside the sensor. On
the other hand, if \shortcode{Hb} includes the whole of the sensor and an
eigenvector data reduction, the result is not even a spectrum in
traditional way, it is just a column of coefficients with a vague
physical meaning.



\section{Calculation accuracy}
%===================
\label{sec:fm_defs:accuracy}

The accuracy of the calculations depends on many factors. For many
factors, such as spectroscopic parameters, there is nothing else to do
than using best available data. On the other hand, for other factors
there is a trade-off between accuracy and speed. More accurate
calculations requires normally also more computer memory. All
different grids and the propagation path step length fall into this
category of accuracy factors. It could be worth discussing the
selection of atmospheric grids and the path step length as there can
be some confusion about how that affects the accuracy.

The main purpose of the atmospheric grids (\builtindoc{p\_grid},
\builtindoc{lat\_grid} and \builtindoc{lon\_grid}) is to build up the
mesh on which the atmospheric fields are defined. This means that the
spacing of these grids shall be selected having the representation of
the atmospheric fields in mind. That is, the spacing shall be fine
enough that the atmospheric field is sufficiently well approximated by
the piecewise (multi-)linear representation between the grid
crossings. The result is that a finer spacing must be used to
represent correctly atmospheric fields with a lot of structure, while
the grids can have fewer points when the atmospheric fields are
smooth. 

The accuracy when performing the actual radiative transfer calculations depends
on the refinement of the expressions used and the discretisation of the
propagation path. If Equation \ref{eq:fm_defs:rte_step} is used, the
underlaying assumption is that the Planck function and the absorption vary
linearly along the propagation path step. These assumptions are of course less
violated if the path step length is made small. An upper limit of the path step
length is set by \wsvindex{ppath\_lmax}. In many cases it should suffice to
just include path points at the crossings of the atmospheric grids
(\builtindoc{ppath\_lmax}$\leq0$). An exception can be limb sounding where the
path step length can be very long around the tangent point, but a limit of
about 25~km should suffice normally.

As points are always included in the propagation paths at the
crossings of the atmospheric grids, finer grids will give shorter path
steps. However, it is neither good practise or efficient to use the
atmospheric grids to control the accuracy of the radiative transfer
calculations. An upper limit on the path step length shall be applied
for this purpose.\footnote{Further discussion can be found in message
  399 and 410 of the ARTS developers mailing list.}


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