\chapter{Clear-sky radiative transfer}
 \label{sec:rte}


 \starthistory
 130220 & Revised after parts moved to a new chapter (Patrick Eriksson).\\
 120831 & Added flowchart and sections on polarised absorption, 
       \builtindoc{iyCalc}, auxiliary data and dispersion (Patrick Eriksson).\\
 110611 & Extended and general revision (Patrick Eriksson).\\
 050613 & First complete version by Patrick Eriksson.\\
 \stophistory

\graphicspath{{Figs/rte/}}

This section discusses variables and the approach used to handle the actual
radiative transfer calculations. This includes how effects caused by the sensor
and surface are incorporated. Measurements of thermal emission in absence of
particle scattering are used as example, and the basic theory for such
simulations is also covered. The first ARTS version was developed for emission
measurements, and such observations remain the standard case in ARTS.

A basic assumption for this chapter is thus that there is no particle
scattering. This is denoted as clear-sky calculations. 
Scattering is restricted to 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 that domain is mainly treated by dedicated workspace methods. For
pure transmission measurements (where scattering into the line-of-sight is
neglected), see Chapter~\ref{sec:trans}. This chapter discusses only the direct
radiative transfer, partial derivatives (i.e.\ the
Jacobian or weighting functions) are discussed in Section~\ref{sec:wfuns}. 

Absorption by atmospheric gases does normally not depend on polarisation but
exceptions exist, where Zeeman splitting is one example. Both polarised and
unpolarised absorption is handled. Even if the gaseous absorption in itself is
unpolarised, the expressions to apply must allow that polarisation signals
from the surface and the cloud box are correctly propagated to the sensor.

For an introduction to a complete radiative transfer calculations, see
Chapter~\ref{sec:complcalcs}. For example, the content of this chapter
corresponds roughly to the flowchart displayed in Figure~\ref{fig:ycalc_flow},
outlining a standard radiative transfer emission calculation. In fact, this
chapter can be seen as a direct continuation of Chapter~\ref{sec:complcalcs}.



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

The structure handling complete radiative transfer calculations is fixed, where
the main workspace method is denoted as \wsmindex{yCalc}
(Fig.~\ref{fig:ycalc_flow}). 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. As workspace variable denoted as
  \wsvindex{iyb}, but can be considered as a pure internal variable and should
  not be of concern for the user.
\item[\aMpiVct{y}]: Monochromatic data for one line-of-sight, i.e.\ a single
  pencil beam calculation. The corresponding workspace variable is
  \wsvindex{iy}. (\aMpiVct{b} consists of one or several \aMpiVct{y} appended.)
\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\_main\_agenda}, giving \aMpiVct{y}}
    \COMMENT{Algorithm \ref{alg:fm_defs:iyCSagenda}}
    \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\_agenda}}
   \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}}
   \STATE{unit conversion of \aMpiVct{y} following \wsvindex{iy\_unit}}
   \COMMENT{Section \ref{sec:fm_defs:unit}}
 \end{algorithmic}
 \caption{The main operations for methods to be part of
   \wsaindex{iy\_main\_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 \wsaindex{iy\_main\_agenda}
(Algorithm~\ref{alg:fm_defs:iyCSagenda}) to perform a single clear sky
radiative transfer calculation. This agenda, in its turn, makes us of other
agendas, such as \builtindoc{ppath\_agenda}. All methods developed for
\builtindoc{iy\_main\_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\_main\_agenda}. 
The three following sections describes the main calculation steps of 
\builtindoc{iy\_main\_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
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}[p]
 \begin{center}
  \includegraphics*[width=0.99\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)$; but note that these no longer are
    explicitly included as path point (in contrast to ARTS-2.0 and earlier).
    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}[p]
 \begin{center}
  \includegraphics*[width=0.99\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 or to 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. In short, 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.

Controlled by \builtindoc{ppath\_step\_agenda}, propagation paths can be
calculated purely geometrically or considering refraction. When considering
refraction, the refractive index is determined at each point along the path
according to \builtindoc{refr\_index\_air\_agenda}. Details about different
methods applicable within \builtindoc{refr\_index\_air\_agenda} are given in
Chapter~\ref{sec:rindex}.

If nothing else is stated, it assumed that all frequency components share a
single propagation path. Another way to express this assumption is that
\textindex{dispersion} is neglected. See Section~\ref{sec:fm_defs:dispersion}
for how to consider dispersion. In the non-dispersive case, the propagation
path is valid for average of the first and last element in \wsvindex{f\_grid},
as this is the frequency given to \wsaindex{refr\_index\_air\_agenda}.

Propagation paths can be calculated separately by the method
\wsmindex{ppathCalc}, but for standard calculations the propagation paths are
calculated internally by \builtindoc{yCalc}. Methods and variables to control
the path calculations are discussed in Section~\ref{sec:ppath:usage}. 



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

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. It is the task of \wsaindex{iy\_surface\_agenda} to return this
  up-welling radiation from the surface, see further
  Chapter~\ref{sec: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 could be
determined by a new call of \builtindoc{iy\_main\_agenda} and the radiative
background for that calculation is then space or the cloud box. The intensity
field entering the cloud box is in some cases calculated by calls of
\builtindoc{iy\_main\_agenda} (with cloud box deactivated) and the radiative
background for these calculations 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.




\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. As mentioned, in this chapter focus is put on emission
measurements. 
The equation to solve is Equation~\ref{eq:VRTE1}:
\begin{equation}
  \frac{\DiffD\StoVec}{\DiffD \PpathLng} = \aAbsMat{a}\left[ \EmsVec- \StoVec
  \right] = -\aAbsMat{a}\StoVec + \Planck\aAbsVec{a},\nonumber
\end{equation}
where the involved quantities are defined and discussed in
Section~\ref{sec:rteq}.


\subsection{Unpolarised absorption}

Let's start with the simpler case of non-polarised absorption (that is, the
absorption is independent of polarisation state). For unpolarised absorption the
matrix \aAbsMat{a}\ is diagonal, with all diagonal elements equal, and only the
first of the elements of \aAbsVec{a})\ 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}
  \aStoI{i+1} = \aStoI{i}e^{-\aOth{i}} + \left(\bar{\Planck}_i+{\bar{j_{n,i}}}/{\bar{\AbsCoef}}\right)(1-e^{-\aOth{i}}),
\end{equation}
with
\begin{eqnarray}
  \bar{\Planck}_i &=& (\Planck(\aTmp{i})+\Planck(\aTmp{i+1}))/2, \\
  \label{eq:bbarstep}
  \bar{\AbsCoef}&=&(\aAbsCoef{i}+\aAbsCoef{i+1})/2\\
  \label{eq:absbar}
  \bar{j_{n,i}}&=&(j_{n,i}+j_{n,i+1})/2\\
  \label{eq:jnbar}
  \aOth{i} &=& \Delta\aPpathLng{i}\bar{\AbsCoef},  
  \label{eq:taustep}
\end{eqnarray}
where \aStoI{i}, \aTmp{i}\ and \aAbsCoef{i}\ are the radiance, temperature and
absorption coefficient, respectively, at point $i$ of the propagation path, and
$\Delta\aPpathLng{i}$ is the distance along the path between point $i$ and
$i+1$. That is, $\bar{\Planck}_i$ is an average of the Planck function at the
path step end points, and the absorption is assumed to vary linearly between
the two points. The start value of \StoI\ is governed by the radiative
background (Section \ref{sec:fm_defs:rad_bkgr}).

A consequence of unpolarised absorption is that also the emission is
unpolarised, and the emission term vanishes for higher Stokes elements.
Accordingly, the expression for the second Stokes component is
\begin{equation}
  \aStoQ{i+1}(\Frq) = \aStoQ{i}(\Frq)e^{-\aOth{i}}.
  \label{eq:fm_defs:rte_step2}
\end{equation}
The third and forth Stokes component are handled likewise. The expressions
above are implemented in the workspace method \wsmindex{iyEmissionStandard},
intended to be part of \builtindoc{iy\_main\_agenda}.  The non-LTE term is
as in Equation~\ref{eq:nlte_source}.

An alternative way to perform the calculations for the first Stokes element
would be
\begin{equation}
  \label{eq:fm_defs:rte_alt}
  \StoI = \sum_i \aTrn{i+1}  \left(\bar{\Planck}_i+{\bar{j_{n,i}}}/{\bar{\AbsCoef}}\right)(1-e^{-\aOth{i}}),
\end{equation}
where \StoI\ is the final intensity and \aTrn{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 (\aStoI{i} must be
known, Section~\ref{sec:wfuns}). However, the calculation of weighting
functions is simplified if \aTrnMat{i}\ is at hand, and this quantity is
also tracked by \builtindoc{iyEmissionStandard}.


\subsection{Polarised absorption}

The overall calculation procedure is the same with polarised absorption, the
only difference is the radiative transfer expression applied. The calculations
for the different Stokes components can here not be separated, and
matrix-vector notation is required: 
\begin{equation}
  \label{eq:fm_defs:vrte_step}
  \aStoVec{i+1} = e^{-\Delta\aPpathLng{i}\bar{\ExtMat}_i} \aStoVec{i} + 
                  (\IdnMtr-e^{-\Delta\aPpathLng{i}\bar{\ExtMat}_i})\left(\bar{\EmsVec}_i+\bar{\ExtMat}_i^{-1}\bar{\aSrcVec{n,i}}\right),
\end{equation}
where \IdnMtr\ is the identity matrix. The \ExtMat\ and \EmsVec\
at point $i$ and $i+1$ are averaged (element-wise) to give $\bar{\ExtMat}$ and
$\bar{\EmsVec}$, respectively, in line with Equation~\ref{eq:bbarstep}. The
calculation of the transmission matrix,
\begin{equation}
  \aTrnMat{i} = e^{-\Delta\aPpathLng{i}\bar{\ExtMat}_i},
  \label{eq:rte:transmat}
\end{equation}
involves a matrix exponential. This calculation step is handled for simpler
cases with analytical expressions, while for more complex cases the Pad\'e
approximation (\developer, Section \ref{D-sec:lin_alg:mat_exp}) is applied.
Only the first element of $\bar{\EmsVec}$ is non-zero, and only the first
column of the matrix corresponding to the term
$(\IdnMtr-e^{-\Delta\aPpathLng{i}\bar{\ExtMat}_i})$ is of interest.



\subsection{Blackbody and cosmic background radiation}

As mentioned, the term \Planck\ is the Planck function.
In ARTS the following version of the Planck function is used:
\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.) The expression in
Equation~\ref{eq:planck} deviates from the exact definition (see
Eq.~\ref{T-eq:rtetheory:Planck} in \theory) as it includes \speedoflight\
instead of the local propagation speed ($v$). The reason for this is the
n$^2$-law of radiance, discussed in the section below.

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
\wsaindex{iy\_space\_agenda} is \wsmindex{MatrixCBR} (together with some calls
of \builtindoc{Ignore}).



\section{Output unit and the n$^2$-law}
%==============================================================================
\label{sec:fm_defs:unit}

First of all, it should be noticed that ARTS does not enforce any 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}). Conversion to other units is selected by the
\wsvindex{iy\_unit} workspace variable. The standard manner is to apply the
unit conversion as part of the calculations performed inside
\builtindoc{yCalc}. See the built-in documentation of the workspace method you
have selected for \builtindoc{iy\_main\_agenda} for comments on practical
aspects and available output units. The most extensive support for conversion
to other units is provided by \builtindoc{iyEmissionStandard}, while other
methods have no support at all (ie.\ they ignore \builtindoc{iy\_unit}). It is
also possible to change the unit as a post-processing step by
\wsmindex{yApplyUnit} (or \wsmindex{iyApplyUnit}), but some restrictions apply
and there are no automatic checks if the input data have correct unit. Further
considerations and expressions for the unit conversion are discussed in the
ARTS-2 journal paper \citep[][Sec.~5.7]{eriksson:arts2:11}.

The n$^2$-law of radiance\index{n2-law of radiance} is introduced in
Section~\ref{T-sec:n2law} of \theory. As shown in that section, the main impact
of the law is handled by consistently using the vacuum speed in the definition
of the Planck radiation law, as done inside ARTS (Eq.~\ref{eq:planck}). This
suffices if the sensor is placed in space (where the refractive index is 1), or
if you use brightness temperatures. Remaining cases are also handled exactly if
\builtindoc{iyEmissionStandard} is used. For those remaining cases, radiance
data shall be scaled with the refractive index squared at the observation
position. For Earth, the maximum value of this factor is about 0.1\,\%, and can
anyhow normally be neglected.

In summary, there is normally no need for you as an user to consider the
n$^2$-law. The exception is if you extract radiance data for a point inside an
atmosphere, and the refractive index deviates significantly from 1 at this
point.




\section{Single pencil beam calculations}
%==============================================================================
\label{sec:fm_defs:single_beams}

The text above assumes that \wsmindex{yCalc} is used. This method can always be
used, but \builtindoc{yCalc} is not mandatory if the simulations only deal with
monochromatic data for a single line-of-sight. In this case, it could be more
handy to use \wsmindex{iyCalc}, which basically is a direct call of
\wsaindex{iy\_main\_agenda}. A reason for selecting \builtindoc{iyCalc} is that
a larger set of auxiliary quantities can be extracted
(Sec.~\ref{sec:fm_defs:aux}). 

On the input side, the main difference when using \builtindoc{iyCalc} is that
the observation position and line-of-sight are specified by \wsvindex{rte\_pos}
and \wsvindex{rte\_los} (instead of \builtindoc{sensor\_pos} and
\builtindoc{sensor\_los}). The calculated radiances are returned as the matrix
\wsvindex{iy} (instead of the vector \builtindoc{y}). 
No automatic unit conversion is made inside \builtindoc{iyCalc}. This is
instead handled separately by \wsmindex{iyApplyUnit}.




\section{Dispersion}
%==============================================================================
\label{sec:fm_defs:dispersion}

The clear-sky radiative transfer methods handle all frequencies in
\builtindoc{f\_grid} in parallel, for efficiency reasons. One consequence of
this feature is that only a single propagation path is calculated, that is
assumed to be common for all frequencies. With other words,
\textindex{dispersion} is not considered. This is in general an acceptable
simplification, but exceptions exist where one example is radiative
transfer through the ionosphere at frequencies approaching the ``plasma
frequency''.

When dispersion is expected to give a significant impact on the results, ARTS
offers a general solution. Dispersion can be handled by setting
\wsaindex{iy\_main\_agenda} as:
\begin{code}
AgendaSet( iy_main_agenda ){
  iyLoopFrequencies
}
\end{code}
The radiative transfer method you put in \builtindoc{iy\_main\_agenda} for
non-dispersive calculations are now moved to \wsaindex{iy\_sub\_agenda}. For
example, if \builtindoc{iyEmissionStandard} is the method of your choice:
\begin{code}
AgendaSet( iy_sub_agenda ){
  iyEmissionStandard
}
\end{code}
The approach is simple, \wsmindex{iyLoopFrequencies} calls
\builtindoc{iy\_sub\_agenda} for each single frequency in \builtindoc{f\_grid}
and appends the output. With some details, \builtindoc{iyLoopFrequencies}
performs a loop over the \builtindoc{f\_grid}, creates an internal
\builtindoc{f\_grid} of length 1 holding the frequency of concern and calls
\builtindoc{iy\_sub\_agenda} with this length-1 frequency grid. This has the
result that a propagation path is calculated for each frequency component.

Some more steps are required to correctly include dispersion. A basic demand is
that \wsaindex{ppath\_agenda} considers refraction. Further,
\wsaindex{refr\_index\_air\_agenda} must provide a dispersive refractive
index. Most methods aimed for \builtindoc{refr\_index\_air\_agenda} give a
refractive index that does not varies with frequency. An example on the
opposite is \wsmindex{refr\_index\_airFreeElectrons}. If a method with
dispersive refractive index is used for non-dispersive calculations, it
receives the mean of the first and last element in \builtindoc{f\_grid} (as
already commented above).

A limitation of \builtindoc{iyLoopFrequencies} is that it can not be combined
with auxiliary data of along-the-path character (Sec.~\ref{sec:fm_defs:aux})






\section{Auxiliary data}
%==============================================================================
\label{sec:fm_defs:aux}

The core output of the radiative calculations is \builtindoc{y}
(\builtindoc{iy} if \builtindoc{iyCalc} is used, jacobians discussed in
Sec.~\ref{sec:wfuns}), but different auxiliary data can be extracted. First of
all, \wsmindex{yCalc} outputs automatically \wsvindex{y\_f}, \wsvindex{y\_pol},
\wsvindex{y\_pos} and \wsvindex{y\_los}. These data give information about the
frequency, polarisation, sensor position and sensor bore-sight, respectively,
corresponding to each value in \builtindoc{y}. The content of the variables are
governed by the sensor settings and the order calculated radiances are stored
(discussed in Sec.~\ref{sec:fm_defs:seqsandblocks}).

A more general mechanism for extracting auxiliary data is controlled by the
\wsvindex{iy\_aux\_vars} workspace variable. This mechanism is most useful
together with \builtindoc{iyCalc}, and for the moment we assume that this
method is used (limitations for \builtindoc{yCalc} are discussed below). The
quantities that can be extracted differ, see the built-in documentation for
the options for each workspace method of concern, e.g.:
\begin{code}
arts -d iyEmissionStandard
\end{code}
The options for this particular method  (\wsmindex{iyEmissionStandard}) can be
divided into different groups (more variables will/can be added):
\begin{description}
\item[Atmosphere, along-the-path] The pressure, temperature and volume mixing
  rations along the propagation path.
\item[Attenuation, along-the-path] Total and species specific absorption
  coefficients along the propagation path.
\item[Radiative properties, along-the-path] The radiance at each propagation
  path point.
\item[Overall radiative properties] The total (clear-sky) optical depth along
  the path and flag giving the radiative background.
\end{description}
``Along-the-path'' means that data are provided for each point of the
propagation path. The path is described by \wsvindex{ppath}, that is also
returned by \wsaindex{iy\_main\_agenda}. The \builtindoc{ppath} variable
contains the information needed to geo-position, for example, ``along-the path
temperatures''. 

Example on setting of \wsvindex{iy\_aux\_vars} (again valid for
\builtindoc{iyEmissionStandard}):
\begin{code}
ArrayOfStringSet( iy_aux_vars,  
    [ "Temperature", 
      "VMR, species 0",
      "Absorption, summed", 
      "Absorption, species 0",
      "Absorption, species 2",
      "iy", 
      "Optical depth" ] )
\end{code}
The data are outputted in a single variable, \wsvindex{iy\_aux}. This variable
is an array of Tensor4. All dimensions are used when storing e.g.\ the
propagation matrix along the path (for all frequencies of
\builtindoc{f\_grid}). For other types of quantities, one or several dimensions
are set to have length 1. See the built-in documentation for further details,
such as the order of the data dimensions.

Storage of quantities of ``along-the-path'' type assumes that there exists a
common propagation path. This is necessarily not the case for calculations by
\builtindoc{yCalc}. This is the case as a calculation considering an antenna
response includes radiative transfer along several propagation paths. The
points of these paths do not end up on common altitude grid, neither are at a
fixed distance from the sensor. In fact, the number of points of the paths will
likely differ. For this reason, \wsmindex{yCalc} will issue an error if you
in \wsvindex{iy\_aux\_vars} include a quantity of ``along-the-path'' character.

The same applies to dispersion calculations (here the propagation path differs
already between the frequencies), and also \wsmindex{iyLoopFrequencies} gives
also an error if ``along-the-path'' auxiliary data are selected.




\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 require 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 piece-wise (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
underlying 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. See also Section~\ref{sec:ppath:lmax}.


