\chapter{Propagation paths and the geoid}
 \label{sec:ppath}


\starthistory
  030310 & First complete version written by Patrick Eriksson.\\
\stophistory



A propagation path is the way the radiation travels to reach the
sensor for a specified line-of-sight. A general description of
propagation paths is given in Section~\ref{sec:fm_defs:ppaths} and it
can be a good idea to read that section before continuing here. This
section describes how propagation paths are described and calculated.
In addition, at the end of the section some geodetic issues are
discussed, such as the choice of reference ellipsoid for the geoid.


\section{Implementation files}
%===
Variables and functions related to propagation paths are defined in the files:
\begin{itemize}
\item \fileindex{ppath.h}
\item \fileindex{ppath.cc}
\item \fileindex{m\_ppath.cc}
\item \fileindex{m\_atmosphere.cc}
\end{itemize}
The first file, \artsstyle{ppath.h}, contains the definition of the
structure to describe propagation paths, \typeindex{Ppath}. The
second file, \artsstyle{ppath.cc}, contains functions to perform
calculations to determine propagation paths. The third file,
\artsstyle{m\_ppath.cc}, contains the workspace methods related to
propagation paths, but these methods mainly check the input and the
actual calculations are performed by sub-functions in
\artsstyle{ppath.cc}. The fourth file, \artsstyle{m\_atmosphere.cc},
contains methods to set the geoid radius.




\section{Calculation approach}
%===================
\label{sec:ppath:approach}

The propagation paths are calculated in steps, as outlined below in
this section. The path steps are normally from one crossing of the
atmospheric grids to next. This solution is necessary to allow that
the same code is used throughout the program. To introduce propagation
paths steps was necessary to handle the iterative solution for
scattering inside the cloud box, as made clear from
Figure~\ref{fig:scattering:averaging}.

A full propagation path is stored in the workspace variable
\wsvindex{ppath}, that is of the type \artsstyle{Ppath} (see
Section~\ref{sec:ppath:Ppath}). The paths are determined by
calculating a number of path steps. A path step is the path from a
point to the next crossing of either the pressure, latitude or
longitude grid (Figure~\ref{fig:ppath:ex1}). There is one exception to
this definition of a path step, and that is when there is an
intersection with the surface, which ends the propagation path at that
point. The starting point for the calculation of a path step is
normally a grid crossing point, but can also be an arbitrary point
inside the atmosphere, such as the sensor position. Only points inside
the model atmosphere are handled. The path steps are stored in the
workspace variable \wsvindex{ppath\_step}, that is of the same type as
\artsstyle{ppath}. The path steps are calculated by an agenda called
\wsaindex{ppath\_step\_agenda}. Example on methods that can be used in
\artsstyle{ppath\_step\_agenda} are \artsstyle{ppath\_stepGeometric}
and \artsstyle{ppath\_stepRefractionEuler}.

\begin{figure}[!t]
 \begin{center}
  \includegraphics*[width=0.80\hsize]{Figs/ppath/ppath_ex1}
  \caption{Tracking of propagation paths. For legend, see 
    Figure~\ref{fig:ppath:ex2}. The figure tries to visualize how the
    calculations of propagation paths are performed from one grid cell
    to next. In this example, the calculations start directly at the
    sensor position $(\ast)$ as it placed inside the model
    atmosphere. The circles give the points defining the propagation
    path. Path points are always included at the crossings of the grid
    cell boundaries. Such a point is then used as the starting point
    for the calculations inside the next grid cell. }
  \label{fig:ppath:ex1}  
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_ppath

\begin{figure}[!t]
 \begin{center}
   \includegraphics*[width=0.98\hsize]{Figs/ppath/ppath_ex2}
  \caption{As Figure~\ref{fig:ppath:ex1}, but with a length criterion 
    for the distance between the points defining the path.
    The inclusion of the tangent point is not a result of this length
    criterion, it is always included among the path points.}
  \label{fig:ppath:ex2}  
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_ppath


Propagation paths are calculated with the workspace method
\wsmindex{ppathCalc}. The communication between this method and
\artsstyle{ppath\_step\_agenda} is handled by \artsstyle{ppath\_step}.
That variable is used both as input and output to
\artsstyle{ppath\_step\_agenda}.  The agenda gets back
\artsstyle{ppath\_step} as returned to \wsmindex{ppathCalc} and the
last path point hold by the structure is accordingly the starting
point for the new calculations. If a total propagation path shall be
determined, the agenda is called repeatedly until the starting point
of the propagation path is found and \artsstyle{ppath\_step} will hold
all path steps that together make up \artsstyle{ppath}. The starting
point is included in the returned structure.

The path is determined by starting at the end point and moving
backwards to the starting point. The calculations are initiated by
filling \artsstyle{ppath\_step} with the practical end point of the
path. This is either the position of the sensor (true or
hypothetical), or some point at the top of the atmosphere (determined
by geometrical calculations starting at the sensor). This
initialization is not handled by \artsstyle{ppath\_step\_agenda}. 
The field \artsstyle{constant} is set by \artsstyle{ppathCalc}
to the correct value if the sensor is above the model atmosphere.
Otherwise, the field is set to be negative and is corrected by
\artsstyle{ppath\_step\_agenda} at the first call. This procedure is
needed as the propagation path constant changes if refraction is
considered, or not, when the sensor is placed inside the atmosphere.

The agenda performs only calculations to next crossing of a grid, all
other tasks are performed by \wsmindex{ppathCalc}, with one exception.
If there is an intersection with the surface, the calculations stop at
this point. This is flagged by setting the background field of
\artsstyle{ppath\_step}. Beside this, \wsmindex{ppathCalc} checks if
the starting point of the calculations is inside the scattering box or
below the surface level, and check if the last point of the path has
been reached. 

The \artsstyle{ppath\_step\_agenda} put in points along the
propagation path at all crossings with the grids, tangent points and
points of surface reflection. Additional points can be included in the
propagation paths. For example, an upper distance between the points
defining the path can be set for \artsstyle{ppath\_stepGeometric} by
the keyword \artsstyle{lmax} (see Figure~\ref{fig:ppath:ex2}).



\section{The propagation path data structure}
%===================
\label{sec:ppath:Ppath}

A propagation path is represented by a structure of type
\typeindex{Ppath}. This structure holds also auxiliary variables to
facilitate the radiative transfer calculations and to speed up the
interpolation. The fields of \artsstyle{Ppath}, where
the data type is given inside square brackets, are:
\begin{description}

  \item[dim] [Index] The atmospheric dimensionality. This field shall always 
     be equal to the workspace variable \artsstyle{atmosphere\_dim}.
     
   \item[np] [Index] Number of positions to define the propagation
     path. Allowed values are $\geq 0$. The number of rows of
     \artsstyle{pos} and \artsstyle{los}, and the length of
     \artsstyle{z}, \artsstyle{gp\_p}, \artsstyle{gp\_lat} and
     \artsstyle{gp\_lon}, shall be equal to \artsstyle{np}. The length
     of \artsstyle{l\_step} is \artsstyle{np} - 1. If \artsstyle{np}
     $\leq$ 1, the observed spectrum is identical to the radiative
     background. For cases where the sensor is placed inside the model
     atmosphere and \artsstyle{np} = 1, the stored position is
     identical to the sensor position and that position can be used to
     determinate the radiative background (see below).

   \item[refraction] [Index] A flag (0 or 1) to indicate if refraction
     has been considered when determining the path. A value of 1 means
     that refraction has been considered.

   \item[method] [String] A string describing the calculation approach.
     For example, '1D basic geometrical'.
     
   \item[constant] [Numeric] The propagation path constant. Such a
     constant can be assigned to all geometrical paths and for 1D
     cases (with or without refraction). See
     Sections~\ref{sec:ppath:basicgeom} and [**]. This field can be
     initiated to a negative value to indicate that the constant is
     undefined or not yet set. For cases where the constant applies,
     \artsstyle{ppath\_step\_agenda} sets this constant at the first
     call of the agenda if the given value is negative.

   \item[pos] [Matrix] The position of the propagation path points.
     This matrix has \artsstyle{np} rows and up to 3 columns. Each row
     holds a position where column 1 is the radius, column 2 the
     latitude and column 3 the longitude (cf.
     Section~\ref{sec:fm_defs:sensorpos}). The number of columns for
     1D and 2D is 2, while for 3D it is 3. The latitudes are stored
     for 1D cases as these can be of interest for some applications
     and are useful if the propagation path shall be plotted. The
     latitudes for 1D give the angular distance to the sensor (see
     further Section~\ref{sec:fm_defs:atmdim}).
     
     The propagation path is stored in reversed order, that is, the
     position with index 0 is the path point closest to the sensor
     (and equals the sensor position if it is inside the atmosphere).
     The full path is stored also for 1D cases with symmetry around a
     tangent point (in contrast to ARTS-1). 
     
  \item[z] [Vector] The geometrical altitude for each path position. The
     length of this vector is accordingly \artsstyle{np}. This is a help
     variable for plotting and similar purposes. It shall not be used to
     interpolate the atmospheric fields, as pressure is the main altitude
     coordinate.
     
   \item[l\_step] [Vector] The length along the propagation path
     between the positions in \artsstyle{pos}. The first value is the
     length between the first and second point etc. For \artsstyle{np}
     $\geq 2$, the length of the vector is \artsstyle{np} - 1.
     Otherwise it is 0.

   \item[gp\_p] [ArrayOfGridPos] Index position with respect to the
     pressure grid. The structure for grid positions is described in
     Section~\ref{sec:interpolation:gridpos}. 
     
   \item[gp\_lat] [ArrayOfGridPos] As \artsstyle{gp\_p} but with
     respect to the latitude grid.

   \item[gp\_lon] [ArrayOfGridPos] As \artsstyle{gp\_p} but with
     respect to the longitude grid.
     
   \item[los] [Matrix] The line-of-sight of the propagation path at
     each point. The number of rows of the matrix is \artsstyle{np}.
     For 1D and 2D, the matrix has a single column holding the zenith
     angle. For 3D there is an additional column giving the azimuth
     angle. The zenith and azimuth angles are defined in
     Section~\ref{sec:fm_defs:los}. If the radiative background is the
     cloud box, the last position (in \artsstyle{pos}) and
     line-of-sight give the relevant information needed when
     extracting the radiative background from the cloud box intensity
     field.
     
   \item[background] [String] The radiative background for the
     propagation path. The possible
     options for this field are 'space', 'blackbody surface', 'cloud
     box interior' and 'cloud box surface', where the source of
     radiation should be clear the content of the strings.
     
   \item[tan\_pos] [Vector] The position of the tangent point. This
     vector is only set if there exists a tangent point (above the
     surface level), the length of the vector is otherwise 0. The
     tangent point is defined as the point with the lowest radius
     along the path. This means that (the absolute value of) the
     zenith angle at the tangent point is always 90\degree. For 2D
     and 3D this point can deviate from the point with lowest
     geometrical altitude.
     
   \item[geom\_tan\_pos] [Vector] The position of the geometrical
     tangent point. This vector is set for all downward observations.
     Refraction and surface reflections are neglected when calculating
     this tangent point position. This field is not handled by
     \artsstyle{ppath\_step\_agenda}. Definition of the tangent point
     as for \artsstyle{tan\_pos}.
\end{description}



\section{Structure of implementation}
%===================
\label{sec:ppath:structure}

The workspace method for calculating propagation paths is
\wsmindex{ppathCalc}, but this is just a getaway function for
\artsstyle{ppath\_calc}. The main use of \artsstyle{ppathCalc} is to
debug and test the path calculations, and that WSM should normally not
be part of the control file. Propagation paths, or steps, are
generated from inside other functions.


\subsection{Main functions for clear sky paths}
%===================

The master function to calculate full clear sky propagation paths is
\funcindex{ppath\_calc}. This function is outlined in
Algorithm~\ref{alg:ppath:ppath_calc}. The function can be divided into
three main parts, initialization (handled by
\artsstyle{ppath\_start\_stepping}), a repeated call of
\artsstyle{ppath\_step\_agenda} and putting data into the return
structure (\artsstyle{ppath}).

\begin{algorithm}[!t]
 \begin{algorithmic}
  \STATE{check consistency of function input}
  \STATE{call \artsstyle{ppath\_start\_stepping} to set 
    \artsstyle{ppath\_step}}
  \STATE{create an array of \artsstyle{Ppath} structures, 
         \artsstyle{ppath\_array}}
  \STATE{add \artsstyle{ppath\_step} to \artsstyle{ppath\_array}}
  \WHILE{radiative background not reached}
   \STATE{call \wsaindex{ppath\_step\_agenda}}
   \IF{path is at the highest pressure surface}
    \STATE{radiative background is space}
   \ELSIF{path is at either end point of latitude or longitude grid}
    \STATE{this is not allowed, issue an runtime error}
   \ENDIF
   \IF{cloud box is active}
    \IF{path is at the surface of the cloud box}
     \STATE{radiative background is the cloud box surface}
    \ENDIF
   \ENDIF
   \STATE{add \artsstyle{ppath\_step} to \artsstyle{ppath\_array}}
  \ENDWHILE
  \STATE{initialize the WSV \artsstyle{ppath} to hold found number of 
         path points} 
  \STATE{copy data from \artsstyle{ppath\_array} to \artsstyle{ppath}}
 \end{algorithmic}
 \caption{Outline of the function \artsstyle{ppath\_calc}.}
 \label{alg:ppath:ppath_calc}
\end{algorithm}

The main task of the function \funcindex{ppath\_start\_stepping} is to
set up \wsvindex{ppath\_step} for the first call of
\artsstyle{ppath\_step\_agenda}, which means that the practical
starting point for the path calculations must be determined. If the
sensor is placed inside the model atmosphere, the sensor position
gives directly the starting point. For cases when the sensor is found
outside the atmosphere, the point where the path exits the atmosphere
must be determined. The exit point can be determined by pure
geometrical calculations (see Sections \ref{sec:ppath:basicgeom} and
\ref{sec:ppath:stepcalc}) as the refractive index is assumed to have the
constant value of 1 outside the atmosphere. The problem is accordingly
to find the geometrical crossing between the limit of the atmosphere
and the sensor line-of-sight (LOS). The function performs further some
other tasks, which include:
\begin{itemize}
\item For all LOS with a zenith angle $\geq 90\degree$ the position of
  the geometrical tangent point is calculated.
\item If the sensor is placed inside the model atmosphere
  \begin{itemize}
  \item Checks that the sensor is placed above the surface level. If
    not, an error is issued.
  \item Checks for 2D and 3D and when the sensor position as at an end
    point of the latitude or longitude grid, that the LOS is inwards
    with respect to the atmospheric limit.
  \item If the sensor and surface altitudes are equal, and the sensor
    LOS is downward, the radiative background is set to be the
    surface. For 2D and 3D, the tilt of the surface radius is considered
    when determining if the LOS is downward.
  \item If the cloud box is active and the sensor position is inside
    the cloud box, the radiative backsurface is set to be ``cloud box
    interior''. All sensor positions on the cloud box surface are for
    2D and 3D treated as points inside the box (for simplicity
    reasons), while for 1D the behavior is as expected.
  \end{itemize}
\item If the sensor is placed outside the model atmosphere
  \begin{itemize}
  \item Checks that the zenith angle is $\geq  90\degree$.  Upward
    observations are here not allowed.
  \item If it is found for 2D and 3D that the exit point of the path
    not is at the top of the atmosphere, but is either at a latitude
    or longitude end face of the atmosphere, an error is issued. This
    problem can not appear for 1D.
  \end{itemize}
\end{itemize}
For further details, see the code.


\subsection{Main functions for propagation path steps}
%===================

Example on workspace methods to calculate propagation path steps are
\wsmindex{ppath\_stepGeometric} and
\wsmindex{ppath\_stepRefractionEuler}. All such methods adapt
automatically to the atmospheric dimensionality, but the different
dimensionalities are handled by separate internal functions. For
example, the sub-functions to \artsstyle{ppath\_stepGeometric} are
\funcindex{ppath\_step\_geom\_1d}, \funcindex{ppath\_step\_geom\_2d}
and \funcindex{ppath\_step\_geom\_3d}. See \artsstyle{m\_ppath.cc} to
get the names of the sub-functions for other propagation path step
workspace methods. The variables to describe the atmosphere are
compacted for 1D and 2D when handed over to atmospheric dimensionality
specific sub-functions. For example, the variable in
\artsstyle{ppath\_step\_geom\_2d} for the geoid radius is a vector,
while the workspace variable is a matrix.

\begin{algorithm}[!t]
 \begin{algorithmic}
  \STATE{call \funcindex{ppath\_start\_2d}}
  \IF{\artsstyle{ppath\_step.ppc} $<1$}
   \STATE{calculate the path constant}
   \COMMENT{this is then first path step}
  \ENDIF
  \STATE{call \funcindex{do\_gridcell\_2d}}
  \STATE{call \funcindex{ppath\_end\_2d}}
  \IF{calculated step ends with tangent point}
   \STATE{call \artsstyle{ppath\_step\_geom\_2d} with temporary 
     \artsstyle{Ppath} structure}
   \STATE{append temporary \artsstyle{Ppath} structure to 
     \artsstyle{ppath\_step}}
  \ENDIF
 \end{algorithmic}
 \caption{Outline of the function \funcindex{ppath\_step\_geom\_2d}.}
 \label{alg:ppath:ppath_step_geom_2d}
\end{algorithm}

Many tasks are independent of the algorithm for refraction that is
used, or if refraction is considered at all. These tasks are solved by
two functions for each atmospheric dimensionality. For 1D the
functions are \funcindex{ppath\_start\_1d} and
\funcindex{ppath\_end\_1d}, and the corresponding functions for 2D and
3D are named in the same way. The functions to calculate geometrical
path steps are denoted as \funcindex{do\_gridrange\_1d},
\funcindex{do\_gridcell\_2d} and \funcindex{do\_gridcell\_3d}. Paths
steps passing a tangent point are handled by a recursive call of the
step function. Algorithm~\ref{alg:ppath:ppath_step_geom_2d} summerizes
this for geometrical 2D steps.


\section{General comments}
%===================
\label{sec:ppath:comments}

The calculation of propagation paths involves a number of mathematical
expressions and they are presented in
Sections~\ref{sec:ppath:basicgeom}~-~\ref{sec:ppath:refreuler}. In
addition, the path calculations present a number of practical
problems. These practical problems are discussed briefly in this
section. For further details, see the code.


\subsection{Numerical precision}
%===================

The aim here is not to make a complete discussion around the limited
numerical accuracy, but just to point out some of the problems caused.
We can start by noticing that the precision with which atmospheric
positions can be given is about 0.5~m when the numeric type is
\textindex{float} and 2\topowerten{-8}~m for \textindex{double}
(assuming that the mantissa has 24 and 48 bits, respectively). The
numbers given correspond to the change of the position for a change of
1 bit, in either radius, latitude and longitude. Already these numbers
cause problems for the approach taken to calculate propagation paths.
For any path along the border of a grid cell, any rounding error in
the wrong direction will move the position outside the grid cell,
which would lead to a crash of the code without countermeasures.

The values above give the representation precision. The precision will
be even poorer if a position is obtained by calculations as numerical
problems tend to accumulate. The calculation precision depends on what
mathematical expressions that are involved.  For example, a radius or
length obtained by the Pythagorean relation will have a relatively
high uncertainty as the calculations involve taking the square of a
radius in the order of 6400~km. It was found that for calculations
performed using only float as numeric type, could lead to
displacements from the true position up to 10~m. It was first tried to
hard-code double as the numerical type for the most critical passages
of the calculations, but a total success was not achieved and some
code had to be duplicated (to be used with either the float or double
option by if-statements for the pre-compiler) to avoid compiler
warnings. A step further was then taken, and double is now hard-coded
for all internal variables of \fileindex{ppath.cc}. This deviation
from the rule to have an uniform numeric type inside ARTS was
introduced to avoid more complicated coding and it has a very small
impact on the overall calculation speed. However, this measure will
not lead to that the precision of the path calculations will be the
same for float and double, as the results will be converted to float
between each propagation path step when copied to
\artsstyle{ppath\_step}.

As pointed out above, the most critical cases are when the path goes
along the boundary of a grid cell. This situation is not common for
arbitrary observation positions, but it is a standard case for 3D
scattering calculations as the starting point for the calculations
there is always a crossing point of the atmospheric grids. The
solution to this problem is to introduce special treatment for such
geometrical paths. For strictly vertical 2D and 3D paths, the
latitude, and also longitude for 3D, of the start and end points shall
be identical. Paths in 3D with an azimuth angle of 0\degree or
180\degree\ have a constant longitude; the paths are in the north-south
plane, and this should also then be valid for the longitude value of
the start and end positions of the path step.

The variables connected to different problems associated with the
numerical inaccuracy and singularity of mathematical expressions are
defined at the top of the file \artsstyle{ppath.cc}. The variables
include the accepted tolerance when making asserts in internal
functions that the given point is inside the specified grid cell.
Another example is the latitude limit to use the special mathematical
expressions needed for positions on the poles.



\subsection{Propagation paths and grid positions}
%===================

The grid positions are calculated on the same time as the path is
determined. The main reason to this is that the grid positions make it
possible to quickly determine inside which grid box the path step is
found. Without the grid positions, each call of the functions would
need a costly search to locate the starting position with respect to
the grids. If you are not familiar with grid positions, it is
recommended to read Section~\ref{sec:interpolation} before you
continue here.

The limited numerical accuracy requires some care when setting the
grid positions. First of all, rounding errors can give a fractional
distance $< 0$ or $> 1$ and this must be avoided. The function
\funcindex{gridpos\_check\_fd} was created for this purpose, and
should be called for each grid position. This function just
sets all values below 0 to 0 and all value above 1 to 1. In addition,
the grid position for the end point of a path step (beside when there
is an intersection with the ground) must have one fractional
distance of exactly 0 or 1, but this is not ensured by
\artsstyle{gridpos\_check\_fd} and for end points the function
\funcindex{gridpos\_force\_end\_fd} shall also be called.

Some care is needed to determine in which grid range a path step is
found. First of all, there exists an ambiguity for the fractional
distance at the grid points. It can either be 0 or 1. In addition, if
a position is exactly on top of a grid point, the observation
direction determines the interesting grid range. As an help to resolve
these question there is the function \funcindex{gridpos2gridrange}.
This function takes an argument describing the direction of the
line-of-sight with respect to the grids. This argument shall be set to
1 if the viewing direction is towards higher indexes. The direction
argument can be set with the following logical expressions, for the
different combinations of atmospheric dimensionality and grid of
interest:

 {\bf 1D-3D, pressure}: $\quad |\ZntAng| \leq 90\degree$

 {\bf 2D, latitude}: $\quad \ZntAng \geq 0\degree$

 {\bf 3D, latitude}: $\quad \AzmAng \leq 90\degree$

 {\bf 3D, longitude}: $\quad \AzmAng \geq 0\degree$





\section{Some basic geometrical relationships for 1D and 2D}
%===================
\label{sec:ppath:basicgeom}

This section gives some expressions to determine positions along a
propagation path when refraction is neglected. The expressions deal
only with propagation path inside a plane, where the latitude angle is
the angular distance from an arbitrary point. This means that the
expressions given here can be directly applied for 1D and 2D. Some of
the expression are also of interest for 3D. The ARTS method for making
the calculation of concern is given inside parenthesis above each
equation, if not stated explicitly. A part of a geometrical
propagation path is shown in Figure~\ref{fig:ppath:1d2dgeom}.

The law of sines gives that the product must $\Rds\sin(\ZntAng)$ be
constant along the propagation path:
\begin{equation}
  p_c = \Rds\sin(\ZntAng)
  \label{eq:ppath:geomconst}
\end{equation}
where the absolute value is taken for 2D zenith angles as they can for
such cases be negative. The propagation path constant, $p_c$, is
determined by the position and line-of-sight of the sensor, a
calculation done by the function \funcindex{geometrical\_ppc}. The
constant equals also the radius of the tangent point of the path (that
is found along an imaginary prolongation of the path behind the sensor
if the viewing direction is upwards). The expressions below are based
on $p_c$ as the usage of a global constant for the path should
decrease the sensitivity to numerical inaccuracies. If the
calculations are based solely on the values for the neighboring
point, a numerical inaccuracy can accumulate when going from one point
to next. The propagation path constant is stored in the field
\artsstyle{constant} of \wsvindex{ppath} and \wsvindex{ppath\_step}.

\begin{figure}[!t]
 \begin{center}
  \begin{minipage}[c]{0.65\textwidth}
   \begin{center}
    \includegraphics*[width=0.9\hsize]{ppath/geom1d}
   \end{center}
  \end{minipage}%
  \begin{minipage}[c]{0.35\textwidth}
   \caption{The radius (\Rds) and zenith angle (\ZntAng) for two points along
     the propagation path, and the distance along the path ($\Delta\PpathLng$)
     and the latitude difference ($\Delta\Lat$) between these points.}
   \label{fig:ppath:1d2dgeom}
  \end{minipage}
 \end{center}
\end{figure}   

The relationship between the distance along the path for an
infinitesimal change in radius is here denoted as the
\textindex{geometrical factor}, $g$. If refraction is neglected, valid
expressions for the geometrical factor are
\begin{equation}
  g = \frac{\DiffD l}{\DiffD r} 
           = \frac{1}{\cos(\ZntAng)} = \frac{1}{\sqrt{1-\sin^2(\ZntAng)}}
                                            = \frac{\Rds}{\sqrt{\Rds^2-p_c^2}}
  \label{eq:ppath:g_geom}
\end{equation}
For the radiative transfer calculations, only the distance between the
points, $\Delta \PpathLng$, is of interest, but for the internal
propagation path calculations the length from the tangent point (real
or imaginary), \PpathLng, is used. By integrating
Equation~\ref{eq:ppath:g_geom}, we get that
(\funcindex{geomppath\_l\_at\_r})
\begin{equation}
  \PpathLng(\Rds) = \sqrt{\Rds^2-p_c^2} 
  \label{eq:ppath:r2l}
\end{equation}
As refraction is here neglected, the tangent point, the point of
concern and the center of the coordinate system make up a right
triangle and Equation~\ref{eq:ppath:r2l} corresponds to the
Pythagorean relation where $p_c$ is the radius of the tangent point.
The distance between two points ($\Delta \PpathLng$) is obtained by
taking the difference of Equation~\ref{eq:ppath:r2l} for the two
radii.

The radius for a given \PpathLng\ is simply (\funcindex{geomppath\_r\_at\_l})
\begin{equation}
  \Rds(\PpathLng) = \sqrt{\PpathLng^2+p_c^2} 
  \label{eq:ppath:l2r}
\end{equation}
The radius for a given zenith angle is simply obtained by rearranging 
Equation~\ref{eq:ppath:geomconst} (\funcindex{geomppath\_r\_at\_za})
\begin{equation}
  \Rds(\ZntAng) = \frac{p_c}{sin(\ZntAng)}
  \label{eq:ppath:za2r}
\end{equation}
The zenith angle for a given radius is (\funcindex{geomppath\_za\_at\_r})
\begin{equation}
  \ZntAng(\Rds) = \left\{
   \begin{array}{ll}
    180 - \sin^{-1}(p_c/\Rds) & 
                   \textrm{for}\quad 90\degree < \aZntAng{a} \leq 180\degree\\
    \sin^{-1}(p_c/\Rds) & 
                   \textrm{for}\quad 0\degree \leq \aZntAng{a} \leq 90\degree\\
    -\sin^{-1}(p_c/\Rds) & 
                   \textrm{for}\quad -90\degree \leq \aZntAng{a} < 0\degree\\
    \sin^{-1}(p_c/\Rds) - 180 & 
                  \textrm{for}\quad -180\degree \leq \aZntAng{a} < -90\degree\\
   \end{array}   \right.
  \label{eq:ppath:r2psi}
\end{equation}
where \aZntAng{a} is any zenith angle valid for the path on the same
side of the tangent point. For example, for a 1D case, the part of the
path between the tangent point and the sensor has zenith angles
$90\degree < \aZntAng{a} \leq 180\degree$.

The latitude for a point (\funcindex{geomppath\_lat\_at\_za}) is most 
easily determined by its zenith angle \\

\begin{equation}
  \Lat(\ZntAng) = \aLat{0} + \aZntAng{0} - \ZntAng
  \label{eq:ppath:za2lat}
\end{equation}
where \aZntAng{0} and \aLat{0} are the zenith angle and latitude of some 
other point of the path. Equation~\ref{eq:ppath:za2lat} is based on the 
fact that the quantities \aZntAng{1}, \aZntAng{2} and $\Delta\Lat$
fulfill the relationship
\begin{equation}
  \Delta\Lat = \aZntAng{1} - \aZntAng{2},
  \label{eq:ppath:dlat}
\end{equation}
this independently of the sign of the zenith angles. The definitions
used here result in that the absolute value of the zenith angle always
decreases towards zero when following the path in the line-of-sight
direction, that is, when going away from the sensor. It should then be
remembered that the latitudes for 1D measures the angular distance to
the sensor, and for 2D a positive zenith angle means observation
towards higher latitudes.

The radius for a given latitude (\funcindex{geomppath\_r\_at\_lat})
is obtained by combining Equations \ref{eq:ppath:za2lat} and
\ref{eq:ppath:za2r}.



\section{Calculation of geometrical propagations paths}
%===================
\label{sec:ppath:stepcalc}

This section describes the calculation of geometrical propagation
paths for different atmospheric dimensionalities. That is, the effect
of refraction is neglected. These calculations are performed by the
workspace method \artsstyle{ppath\_stepGeometric}. This method, as all
methods for propagation path steps, adjust automatically to the
atmospheric dimensionality, but the actual calculations are performed
a sub-function for each dimensionality.


\subsection{1D}
%===================
\label{sec:ppath:1Dgeom}

The core function for this case is \funcindex{do\_gridrange\_1d}. The
lowest and highest radius value along the path step is first
determined. If the line-of-sight is upwards ($\|\ZntAng| \leq
90\degree$), then the start point of the step gives the lowest radius,
and the radius of the pressure surface above gives the highest value.
In the case of a downwards line-of-sight, the lowest radius is either
the tangent point, the pressure surface below or the surface. The
needed quantities to describe the propagation path between the two
found radii are calculated by the function
\funcindex{geompath\_from\_r1\_to\_r2}, that has the option to
introduce more points to fulfill a length criterion between the path
points. The mathematics of \artsstyle{geompath\_from\_r1\_to\_r2} are
given by Equations~\ref{eq:ppath:geomconst}~-~\ref{eq:ppath:za2lat}.


\subsection{2D}
%===================
\label{sec:ppath:2Dgeom}

The definitions given in Sections~\ref{sec:fm_defs:atmdim} results in
that for a 2D case the radius of a pressure surface varies linearly
from one point of the latitude grid to next. This is the main
additional problem to solve, compared to the 1D case.
Figure~\ref{fig:ppath:psurf_crossing} gives a schematic description of
the problem at hand, which is handled by the internal function
\funcindex{psurface\_crossing\_2d}.

The law of sine gives the following relationship for the crossing
point:
\begin{equation}
  \frac{\sin\Theta_p}{\aRds{0}+c\Lat} = 
                                \frac{\sin(\pi-\Lat-\Theta_p)}{\aRds{p}}
\end{equation}
which can be re-written to

\begin{figure}
 \begin{minipage}[c]{0.45\textwidth}
 \includegraphics*[width=0.92\textwidth]{ppath/psurf_crossing}
 \end{minipage}%
 \begin{minipage}[c]{0.55\textwidth}
  \caption{Quantities used to describe how to find the crossing between a 
    geometrical propagation path and a tilted pressure surface. The
    angle \Lat\ is the angular distance from a reference point on the
    path. The problem at hand is to find \Lat\ for the crossing
    point. The radius of the pressure surface at $\Lat = 0$ is
    denoted as $r_0$. The tilt of the pressure surface is $c$.}
  \label{fig:ppath:psurf_crossing}
 \end{minipage}%
\end{figure}   

\begin{equation}
   \aRds{p} \sin(\Theta_p) = (\aRds{0}+c\Lat) 
           (\sin\Theta_p \cos\Lat + \cos\Theta_p \sin\Lat)
 \label{eq:ppath:psurf1}
\end{equation}
This equation has no analytical solution. A first step to find an
approximative solution is to note that \Lat\ will be limited to
relatively small values. For example, if it shall be possible for the
angular distance \Lat\ to reach the value of 3\degree, the vertical
spacing between the pressure surfaces must be about 8~km, while it
normally is below 2~km. For angles $\Lat \leq 3\degree$, the sine and
cosine terms can be replaced with the two first terms of their Taylor
expansions with a relative accuracy of $< 4\topowerten{-7}$. That is,
\begin{eqnarray}
  \cos \Lat & \approx & 1 - \Lat^2 / 2    \nonumber \\
  \sin \Lat & \approx & \Lat - \Lat^3 / 6 \nonumber
\end{eqnarray}
Equation~\ref{eq:ppath:psurf1} becomes with these replacements a
polynomial equation of order 4:
\begin{eqnarray}
  \label{eq:ppath:psurf2}
  0 & = & p_0 + p_1 \Lat + p_2 \Lat^2 + p_3 \Lat^3 + p_4 \Lat^4 \\
  p_0 & = & (\aRds{0}-\aRds{p}) \sin\Theta_p \nonumber \\ 
  p_1 & = & \aRds{0}\cos\Theta_p + c \sin\Theta_p \nonumber \\ 
  p_2 & = & -(\aRds{0}\sin\Theta_p)/2 + c\cos\Theta_p  \nonumber \\ 
  p_3 & = & -(\aRds{0}\cos\Theta_p)/6 - (c\sin\Theta_p)/2 \nonumber \\ 
  p_4 & = & -(c \cos\Theta_p)/6  \nonumber 
\end{eqnarray}
This equation is solved numerically with the root finding algorithm
implemented in the function \funcindex{poly\_root\_solve}. Solutions of
interest shall not be imaginary.

\begin{figure}[!t]
 \begin{center}
  \includegraphics*[width=0.80\hsize]{Figs/ppath/ppath_ex3}
  \caption{Example on propagation path steps starting from a latitude end face 
    (solid lines), or the lower pressure surface (dashed lines), to
    all other grid cell faces. The distortion of the grid cell from
    cylinder segment is highly exaggerated compared to a real case.
    The relationship between vertical and horisontal size deviates
    also from normal real cases.  Typical values for the vertical
    extension is around 500~m, while the horisontal length is
    normally $>$~10~km.}
  \label{fig:ppath:ex3}  
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_ppath

Geometrical 2D propagation path steps are determined by
\funcindex{do\_gridcell\_2d}. This function uses
\artsstyle{psurface\_crossing\_2d} to calculate the latitude distance
to a crossing of the pressure surface below and above the present path
point. If the closest crossing point with the pressure surfaces is
outside the latitude range of the grid cell, it is the crossing of the
path with the end latitude (in the viewing direction) that is of
interest (Figure~\ref{fig:ppath:ex3}).



\subsection{3D}
%===================
\label{sec:ppath:3Dgeom}

Geometrical 3D propagation path steps are determined by the function
\funcindex{do\_gridcell\_3d}. It was first tested to use different
analytical expressions to calculate the length between a point and
the crossing of some radius, latitude or longitude. However, the
expressions found include the trigonometric functions and the squaring
of radii, which resulted in a high sensitivity to the numerical
inaccuracy. It was found that the numerical problems made the created
algorithm impossible to use in practice.
Equation~\ref{eq:ppath:rcoss3Db} below is a reminiscence of that work.
In addition, no simple solution to the problem of finding the crossing
with a tilted 3D pressure surface using the analytical expressions was
found.

\begin{algorithm}[!t]
 \begin{algorithmic}
  \STATE{calculate the spherical position $(x_0,y_0,z_0)$ and LOS vector 
         $(\DiffD x,\DiffD y,\DiffD z)$}
  \STATE{calculate $(\Rds_c,\Lat_c,\Lon_c)=S(x_0,y_0,z_0)-(\Rds_0,\Lat_0,\Lon_0)$, the position correction term}
  \STATE{set $l_{in} = 0$} 
  \STATE{set $l_{out} = 1$}
  \IF{LOS is downwards}
   \STATE{calculate length to the tangent point, $l_{tan}$}
  \ELSE
   \STATE{set $l_{tan} = 99\topowerten{6}$ m}
  \ENDIF
  \WHILE{$S(x_0+l_{out}\DiffD x,y_0+l_{out}\DiffD y,z_0+l_{out}\DiffD z)-(\Rds_c,\Lat_c,\Lon_c)$ is inside grid cell}
   \IF{$l_{out}<l_{tan}$ and $10l_{out}>l_{tan}$}
    \STATE{$l_{out} = l_{tan}$}
    \COMMENT{to assure that tangent point is included in search}
   \ELSE
    \STATE{$l_{out} \gets 10 * l_{out}$ }
   \ENDIF
  \ENDWHILE
  \STATE{set $l_{end} = (l_{in}+l_{out})/2$}
  \STATE{set accuracy flag to false}
  \WHILE{accuracy flag is false}
   \STATE{calculate $(\Rds,\Lat,\Lon)=S(x_0+l_{end}\DiffD x,y_0+l_{end}\DiffD y,z_0+l_{end}\DiffD z)-(\Rds_c,\Lat_c,\Lon_c)$}
   \IF{$(\Rds,\Lat,\Lon)$ is inside grid cell}
    \STATE{$l_{in} = l_{end}$}
   \ELSE
    \STATE{$l_{out} = l_{end}$}
   \ENDIF
   \IF{$(l_{out}-l_{in})$ smaller than specified accuracy}
    \STATE{set accuracy flag to true}
   \ELSE
    \STATE{$l_{end} = (l_{in}+l_{out})/2$}
   \ENDIF
  \ENDWHILE
  \STATE{$(\Rds,\Lat,\Lon)\gets(\Rds,\Lat,\Lon)+(\Rds_c,\Lat_c,\Lon_c)$}
 \end{algorithmic}
 \caption{The method applied in \artsstyle{do\_gridcell\_3d} to find the 
          total length of the path step to be calculated.
          The symbol $S$ signifies here conversion from cartesian to
          spherical coordinates (Equation \ref{eq:ppath:cart2sph}).}
 \label{alg:ppath:dogridcell3d}
\end{algorithm}

A straightforward trail-and-error algorithm was then tested
(Algorithm~\ref{alg:ppath:dogridcell3d} and
Figure~\ref{fig:ppath:3Dsearch}). The main advantage of the algorithm
is that a correction for the shift in position caused by the
transformations back and fourth to a cartesian coordinate system can
be applied. The correction term assures that the position is not
changed for a step of zero length, and is not moved outside the grid
cell due to the numerical problems. The algorithm was further found to
be sufficiently fast to be accepted. A simple bisection search to find
the length of the propagation path step is used. Both the position and
the line-of-sight for the other end point of the path step are
calculated using a transformation to cartesian coordinates.  The
cartesian coordinate system used here is defined as:
\begin{description}
\item[x-axis] is along latitude 0\degree and longitude 0\degree
\item[y-axis] is along latitude +90\degree
\item[z-axis] is along latitude 0\degree and longitude +90\degree
\end{description}
This definition results in the following relationships between the
spherical $(\Rds,\Lat,\Lon)$ and cartesian $(x,y,z)$ coordinates
\begin{eqnarray}
  x &=& \Rds \cos(\Lat) \cos(\Lon) \nonumber \\
  y &=& \Rds \sin(\Lat)            \\
  z &=& \Rds \cos(\Lat) \sin(\Lon) \nonumber
 \label{eq:ppath:sph2cart}
\end{eqnarray}
and
\begin{eqnarray}
 \label{eq:ppath:cart2sph}
  \Rds &=& \sqrt{x^2+y^2+z^2}  \nonumber  \\
  \Lat &=& \arcsin(y/\Rds)                \\
  \Lon &=& \arctan(z/x) \qquad \mathrm{(implemented\ by\ the\ atan2\ function)}
                               \nonumber
\end{eqnarray}
The functions performing these transformations are \funcindex{sph2cart} and
\funcindex{cart2sph}.

\begin{figure}[!t]
 \begin{center}
  \includegraphics*[width=0.80\hsize]{Figs/ppath/ppath_3Dsearch}
  \caption{Schematic of Algorithm~\ref{alg:ppath:dogridcell3d}. The
    figure shows two iterations of the algorith to search for the
    total length of the path step. The asterisk $(\ast)$ gives the
    start point for the calculations and the circles $(\circ)$ are the
    final end points of the path step. The plus signs $(+)$ shows the
    position of the different lengths tested during the iterations.}
  \label{fig:ppath:3Dsearch}  
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_ppath


The first step to transform a line-of-sight, given by the zenith
($\ZntAng$) and the azimuth ($\AzmAng$) angle, to cartesian
coordinates is to determine the corresponding vector with unit length
in the spherical coordinate system:
\begin{equation}
 \left[ \begin{array}{c}
  \DiffD \Rds \\
  \DiffD \Lat \\
  \DiffD \Lon
 \end{array} \right] =
 \left[ \begin{array}{c}
   \cos(\ZntAng) \\
   \sin(\ZntAng) \cos(\AzmAng) / \Rds \\
   \sin(\ZntAng) \sin(\AzmAng) / ( \Rds \cos(\Lat) )
 \end{array} \right]
 \label{eq:ppath:los2sphvec}
\end{equation}
This vector is then translated to the cartesian coordinate system as
\begin{equation}
 \left[ \begin{array}{c}
  \DiffD x \\
  \DiffD y \\
  \DiffD z
 \end{array} \right] =
 \left[ \begin{array}{ccc}
  \cos(\Lat)\cos(\Lon) & -\Rds\sin(\Lat)\cos(\Lon) & 
                                                   -\Rds\cos(\Lat)\sin(\Lon) \\
  \sin(\Lat)           & \Rds\cos(\Lat)            & 0                     \\ 
  \cos(\Lat)\sin(\Lon) & -\Rds\sin(\Lat)\sin(\Lon) & \Rds\cos(\Lat)\cos(\Lon) 
 \end{array} \right] 
 \left[ \begin{array}{c}
  \DiffD \Rds \\
  \DiffD \Lat \\
  \DiffD \Lon
 \end{array} \right]
 \label{eq:ppath:los2cart}
\end{equation}
Note that the radial terms (\Rds) in Equations
\ref{eq:ppath:los2sphvec} and \ref{eq:ppath:los2cart} cancel each
other.  These calculations are performed in \funcindex{poslos2cart}.
Special expressions must be used for positions at the north and south
pole (see the code) as the azimuth angle has there a special
definition (Section~\ref{sec:fm_defs:los}).

The cartesian position of a point along the geometrical path at a
distance $l$ is then simply
\begin{equation}
 \left[ \begin{array}{c}
  x_2 \\
  y_2 \\
  z_2
 \end{array} \right] =
 \left[ \begin{array}{c}
  x_1 + l\DiffD x \\
  y_1 + l\DiffD y \\
  z_1 + l\DiffD z
 \end{array} \right]
  \label{eq:ppath:xdl}
\end{equation}
The cartesian viewing vector $[\DiffD x, \DiffD y, \DiffD z]^T$ is
constant along a geometrical path. The new position is converted to
spherical coordinates by Equation~\ref{eq:ppath:cart2sph} and the new
spherical viewing vector is calculated as
\begin{equation}
 \left[ \begin{array}{c}
  \DiffD \Rds \\
  \DiffD \Lat \\
  \DiffD \Lon
 \end{array} \right] =
 \left[ \begin{array}{ccc}
  \cos(\Lat)\cos(\Lon) & \sin(\Lat)            & \cos(\Lat)\sin(\Lon) \\ 
  -\sin(\Lat)\cos(\Lon)/\Rds & \cos(\Lat)/\Rds & -\sin(\Lat)\sin(\Lon)/\Rds \\ 
  -\sin(\Lon)/(\Rds\cos(\Lat)) & 0             & \cos(\Lon)/(\Rds\cos(\Lat)) 
 \end{array} \right] 
 \left[ \begin{array}{c}
  \DiffD x \\
  \DiffD y \\
  \DiffD z
 \end{array} \right]
 \label{eq:ppath:los2sph}
\end{equation}
which is converted to a zenith and azimuth angle as
\begin{eqnarray}
  \ZntAng &=& \arccos(\DiffD \Rds) \nonumber  \\
  \AzmAng &=& \arccos(\Rds\DiffD\Lat/\sin(\ZntAng)), 
                      \qquad \textrm{for} \quad \DiffD\Lon >= 0 \\
  \AzmAng &=& -\arccos(\Rds\DiffD\Lat/\sin(\ZntAng)), 
                      \qquad \textrm{for} \quad \DiffD\Lon < 0  \nonumber
 \label{eq:ppath:conv2zaaa}
\end{eqnarray}
Special expressions must be used for positions at the north and south
pole (see the code) as the azimuth angle has there a special
definition (Section~\ref{sec:fm_defs:los}). These calculations are performed in
\funcindex{cart2poslos}.

For sensor positions outside the atmosphere, the calculations made in
\funcindex{ppath\_start\_stepping} involve the problem of finding the
position where the path leaves the atmosphere. This position is found
by an iterative search. The maximum radius of the uppermost pressure
surface is taken as first guess for the radius of the exit point.  The
exit latitude and longitude for this radius is determined (as
discussed below), and the radius for the top of the atmosphere for the
found position is used as radius for next iteration. This procedure is
repeated until the change from one iteration to next for both latitude
and longitude is smaller than 1\topowerten{-6}. The exit position for
a given radius, \Rds, is found by solving the following equation
system:
\begin{eqnarray}
  \Rds \cos(\Lat) \cos(\Lon)  & = & x+l\DiffD x \nonumber \\
  \Rds \sin(\Lat)             & = & y+l\DiffD y \\
  \Rds \cos(\Lat) \sin(\Lon)  & = & z+l\DiffD z \nonumber
  \label{eq:ppath:rcoss3Da}
\end{eqnarray}
where $(x,y,z)$ is the position of the sensor, $(\DiffD x,\DiffD
y,\DiffD z)$ the sensor LOS, and $l$, \Lat\ and \Lon\ are the variables
to be determined. The first step is to determine the distance $l$ to
the exit point, which is found by adding the square of all three
equations:
\begin{equation}
  \Rds^2 = (x+l\DiffD x)^2 + (y+l\DiffD y)^2 + (z+l\DiffD z)^2
  \label{eq:ppath:rcoss3Db}
\end{equation}
Once $l$ is determined, the latitude and longitude are easily
calculated by Equations~\ref{eq:ppath:xdl} and
\ref{eq:ppath:cart2sph}. These calculations are implemented in the
function \funcindex{psurface\_crossing\_3d}.  Similar expressions were
derived to find the position for the crossing of a given latitude or
longitude but those expressions were removed from the code as they are
not used with present algorithms.\footnote{The expressions mentioned
  can be extracted from the function
  \artsstyle{gridcell\_crossing\_3d} in ARTS version 1-1-440.}



\section{Refraction with simple Euler scheme}
%===================
\label{sec:ppath:refreuler}

\begin{figure}[!t]
 \begin{center}
  \includegraphics*[width=0.80\hsize]{Figs/ppath/ppath_refr1}
  \caption{Comparison of propagation paths calculated geometrically and 
    with refraction considered, for the same zenith angle of the
    sensor line-of-sight. The figure include two pair of paths, with
    refracted tangent altitude of about 0 and 10~km, respectively.
    The horisontal coordinate is the latitude distance from the point
    where the path exits the model atmosphere (at 80~km). The model
    atmosphere used had a spherical symmetry (that is, 1 D case, but
    the calculations were performed in 2D mode).}
  \label{fig:ppath:ppath_refr1}  
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_refraction

Refraction affects the radiative transfer in several ways. The
distance through a layer of a fixed vertical thickness will be
changed, and for a limb sounding observation the tangent point is
moved both vertically and horizontally. If the atmosphere is assumed
to be horizontally stratified (1D), a horizontal displacement is of no
importance but for 2D and 3D calculations this effect must be
considered. For limb sounding and a fixed zenith angle, the tangent
point is moved downwards compared to the pure geometrical case
(Figure~\ref{fig:ppath:ppath_refr1}), resulting in that inclusion of
refraction in general gives higher intensities. However, the
propagation path is still symmetric around tangent and surface points.

The refraction causes a bending of the path, which gives a deviation
from the geometrical approximation of propagation along a straight
line. The bending of the path is obtained by the relationship
\begin{equation}
  \frac{\DiffD x}{\DiffD l} = \frac{1}{n} \left( \frac{\PartD n}{\PartD y} \right)_x
  \label{eq:ppath:refrbend}
\end{equation}
where $x$ is the direction of propagation, $l$ the distance along the
path, $n$ the refractive index\footnote{The refractive index is here
  assumed to have no imaginary part}, and $y$ is the coordinate
perpendicular to the path. See further Section~9.4 in
\citet{rodgers:00}.

\begin{figure}[!t]
 \begin{center}
  \includegraphics*{Figs/ppath/euler}
  \caption{Schematic of the Euler ray tracing scheme. The ray tracing step 
    length is $l_r$. }
  \label{fig:ppath:euler}  
 \end{center}
\end{figure}

The workspace method \wsmindex{ppath\_stepRefractionEuler} takes
refraction into consideration by probably the most simple (from the
viewpoint of implementation) algorithm possible. This does not mean
that it is the best way to consider refraction, it is rather
inefficient regarding computational burden, and if the step length for
the ray tracing (see below) is made very small, the result can be
completely wrong due to numerical problems.

The approach taken in \artsstyle{ppath\_stepRefractionEuler} is to
take a geometrical ray tracing step from the present point of the path
(and in the direction of present line-of-sight). Refraction is
considered only when the line-of-sight at the new point is determined
(Figure~\ref{fig:ppath:euler}). The found line-of-sight is used to
calculate the next ray tracing step etc. This can be seen as an Euler
solution to the differential problem given by
Equation~\ref{eq:ppath:refrbend}. The main difference between handling
1D, 2D or 3D cases is how the line-of-sight for the new point is
corrected to compensate for the bending due to refraction. The
calculation of propagation paths including the effect of refraction is
often denoted as \textindex{ray tracing}.

The length of the calculation steps is set by the keyword
\artsstyle{lraytrace}. This length shall not be confused with the
final distance between the points that define the path, which is
controlled by the keyword \artsstyle{lmax}. The path is first
determined in steps of \artsstyle{lraytrace}. The found ray tracing
points are then used for an interpolation to create a path step
defined exactly as for geometrical calculations. The normal situation
is that the ray tracing step length is considerably shorter than the
final spacing between the path points. Suitable values for
\artsstyle{lraytrace} have not yet been investigated in detail, but
for limb sounding values in around 1~-~10~km should be appropriate.
Shorter ray tracing steps (down to a level where rounding errors will
start to have an impact) will of course give a propagation path more
accurately determined, but on the cost of more time consuming
calculations.


\subsection{1D}
\label{sec:ppath:refr1D}
%-----------------------
\begin{figure}[!p]
  \begin{center}
    \includegraphics*{Figs/ppath/snell}
    \caption{Geometry to derive Snell's law for a spherical atmosphere. }
    \label{fig:ppath:snell} 
  \end{center} 
\end{figure}
\begin{figure}[!p]
 \begin{center}
  \includegraphics*[width=0.80\hsize]{Figs/ppath/ppath_N}
  \caption{Vertical variation of refractivity $(n-1)\topowerten{6}$.
     Calculated for a mid-latitude summer climatology (FASCODE), where
     the dashed line is for a completely dry atmosphere, and the solid line
     includes also contribution from water vapour.}
  \label{fig:ppath:N}  
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_refraction

When determining the propagation path through the atmosphere
geometrical optics can be applied because the change of the refractive
index over a wavelength can be neglected. Applying Snell's law
to the geometry shown in Figure
\ref{fig:ppath:snell} gives
\begin{equation}
  n_i \sin (\aZntAng{i}) = n_{i+1} \sin (\aZntAng{i'})
\end{equation}
Using the same figure, the law of sines gives the relationship
\begin{equation}
  \frac{\sin(\aZntAng{i+1})}{\aRds{i}} = 
  \frac{\sin(180^\circ-\aZntAng{i+1}')}{\aRds{i+1}} =
  \frac{\sin(\aZntAng{i'})}{\aRds{i+1}} 
\end{equation}
By combining the two equations above, the Snell's law for a spherical
atmosphere (that is, 1D cases) is derived
\citep[e.g.][]{kyle:91,balluch:97}:
\begin{equation}
  p_c = \aRds{i} n_i \sin(\aZntAng{i}) = \aRds{i+1} n_{i+1}\sin(\aZntAng{i+1}) 
 \label{eq:ppath:snellspherical}
\end{equation}
where $c$ is a constant. With other words, the Snell's law for
spherical atmospheres states that the product of $n$, \Rds\ and
$\sin(\ZntAng)$ is constant along the propagation path. It is
noteworthy that with $n=1$, Equations~\ref{eq:ppath:geomconst} and
\ref{eq:ppath:snellspherical} are identical.

The Snell's law for a spherical atmosphere makes it very easy to
determine the zenith angle of the path for a given radius. A
rearrangement of Equation~\ref{eq:ppath:snellspherical} gives
\begin{equation}
  \ZntAng = \arcsin( \Rds n / p_c )
 \label{eq:ppath:za1D}
\end{equation}
This relationship makes it possible to handle refraction for 1D
without calculating any gradients of the refractive index, which is
needed for 2D and 3D. These calculations are implemented in the
function \funcindex{raytrace\_1d\_linear\_euler}.
Figure~\ref{fig:ppath:N} shows the vertical variation of the
refractive index.





\subsection{2D}
\label{sec:ppath:refr2D}
%-----------------------
\begin{figure}[!t]
 \begin{center}
  \includegraphics*[width=0.85\hsize]{Figs/ppath/ppath_dndr}
  \caption{Vertical gradient of the refractive index.
     Calculated for a mid-latitude summer climatology (FASCODE), where
     the dashed line is for a completely dry atmosphere, and the solid line
     includes also contribution from water vapour.}
  \label{fig:ppath:dndr}  
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_refraction

\begin{figure}[!t]
 \begin{center}
  \includegraphics*[width=0.85\hsize]{Figs/ppath/ppath_dndlat}
  \caption{Latitude gradient of the refractive index due to varying radius 
    of the geoid. The gradient is given as the change in refractive
    index over 1~m, which allows direct comparison with the values in
    Figure~\ref{fig:ppath:dndr}e. The wet atmosphere from
    Figure~\ref{fig:ppath:dndr} was used for all latitudes, and the
    the plotted gradient is only caused by the fact that the radius of
    the geoid is not constant.  The gradient is positive on the
    southern hemisphere (shown), and negative on the northern
    hemisphere.}
  \label{fig:ppath:dndlat}  
 \end{center}
\end{figure}
% This figure was produced by the Matlab function mkfigs_refraction

Equation~\ref{eq:ppath:refrbend} expressed in polar coordinates
is \citep[Eq. 9.30]{rodgers:00}
\begin{equation}
  \frac{\DiffD(\Lat+\ZntAng)}{\DiffD l} = 
    -\frac{\sin\ZntAng}{n} \left( \frac{\PartD n}{\PartD \Rds} \right)_\Lat
    +\frac{\cos\ZntAng}{n\Rds} \left( \frac{\PartD n}{\PartD \Lat} \right)_\Rds
  \label{eq:ppath:refrbend2d}
\end{equation}
If the gradients are zero (corresponding to the geometrical case) we
find that the sum of the zenith angle and the latitude is constant
along a 2D geometrical path, which is also made clear by
Equation~\ref{eq:ppath:za2lat}. The geometrical zenith angle at ray
tracing point $i+1$ is accordingly $\aZntAng{i+1} = \aZntAng{i} -
(\aLat{i+1}-\aLat{i})$. If then also the refraction is considered, we
get the following expression:
\begin{equation}
  \aZntAng{i+1} = \aZntAng{i} - (\aLat{i+1}-\aLat{i}) + \frac{l_g}{n_i}
   \left[
    -\sin\aZntAng{i} \left( \frac{\PartD n}{\PartD \Rds} \right)_{\aLat{i}}
    +\frac{\cos\aZntAng{i}}{\aRds{i}} 
                \left( \frac{\PartD n}{\PartD \Lat} \right)_{\aRds{i}}
  \right]  
  \label{eq:ppath:za2d}
\end{equation}
The gradients of the refractive index for 2D are calculated by the
function \funcindex{refr\_gradients\_2d}. This function returns the
gradients as the change of the refractive index over 1~m. The
conversion for the latitude gradient corresponds to the $1/\Rds$ term
found in Equation~\ref{eq:ppath:za2d}, and this term is accordingly
left out in \funcindex{raytrace\_2d\_linear\_euler}, which is the
function of this section. 

The radial and latitudinal gradients of the refractive index are
calculated in pure numerical way, by shifting the position slightly
from the position of concern. Figures~\ref{fig:ppath:dndr} and
\ref{fig:ppath:dndlat} show example on gradients of the refractive
index.




\subsection{3D}
\label{sec:ppath:refr3D}
%----------------------
For 3D, the geomtrical expressions are used to
calculate the geometrical zenith and azimuth angles at the end of the
ray tracing step. Following the methodology for 2D, the geometrical
zenith and azimuth angles are then corrected to incorporate the
influence of refraction. The zenith angle is calculated as
\begin{eqnarray}
  \aZntAng{i+1} &=& \aZntAng{g} - \frac{l_g\sin\aZntAng{i}}{n_i} 
     \left( \frac{\PartD n}{\PartD \Rds} \right)_{(\aLat{i},\aLon{i})} + \\
     & &  + \frac{l_g\cos\aZntAng{i}}{\aRds{i}n_i}
   \left[
     \cos\aAzmAng{i} 
          \left( \frac{\PartD n}{\PartD \Lat} \right)_{(\aRds{i},\aLon{i})}
    +\frac{\sin\aAzmAng{i}}{\cos\aLat{i}} 
          \left( \frac{\PartD n}{\PartD \Lon} \right)_{(\aRds{i},\aLat{i})}
  \right]  \nonumber
  \label{eq:ppath:za3d} 
\end{eqnarray}
where \aZntAng{g} is the zenith angle obtained from the geometrical
expressions. In similar manner, the geometrical azimuth angle,
\aAzmAng{g}, is corrected as
\begin{equation}
  \aAzmAng{i+1} = \aAzmAng{g} + \frac{l_g\sin\aZntAng{i}}{r_in_i} 
   \left[
    -\sin\aAzmAng{i} 
          \left( \frac{\PartD n}{\PartD \Lat} \right)_{(\aRds{i},\aLon{i})}
    +\frac{\cos\aAzmAng{i}}{\cos\aLat{i}}
          \left( \frac{\PartD n}{\PartD \Lon} \right)_{(\aRds{i},\aLat{i})}
  \right]  
  \label{eq:ppath:aa3d} 
\end{equation}
This expression, slightly modified, is found in
\funcindex{raytrace\_3d\_linear\_euler}. The terms of
Equation~\ref{eq:ppath:aa3d} missing in that function, are part of
\funcindex{refr\_gradients\_3d} to convert the gradients to the same
unit. The longitude gradient is converted to the unit [1/m] by
multiplication with the term $1/(\Rds\cos\Lat)$.




\section{Geoid ellipsoids and geodetic datums}
%===================
\label{sec:ppath:geoids}

This section defines the geoid ellipsoid and discusses related
issues. The geoid is introduced in Section~\ref{sec:fm_defs:geoid}.
The workspace variable representing the geoid is \wsvindex{r\_geoid}.


\subsection{Geoid ellipsoids}
%===================
\label{sec:ppath:geoid}

All geodetic datums are based on a reference ellipsoid\index{geoid
  ellipsoid}. The ellipsoid is rotationally symmetric around the
north-south axis. That is, the ellipsoid radius has no longitude
variation, it is only a function of latitude. The ellipsoid is
described by an equatorial radius, \aRds{e}, and a polar radius,
\aRds{p}. These radii are indicated in Figure~\ref{fig:ppath:lats}.
The radius of the ellipsoid for a given latitude is
\begin{equation}
 \aRds{\odot}(\Lat) = \sqrt{\frac{\aRds{e}^2\aRds{p}^2}
                    {\aRds{e}^2\sin^2\Lat+\aRds{p}^2\cos^2\Lat}}
 \label{eq:ppath:ellipsradius} 
\end{equation}
The radius given by Equation~\ref{eq:ppath:ellipsradius} can be
directly applied for 2D and 3D cases. On the other hand, for 1D cases
the reference geoid is by definition a sphere and the radius of this
sphere shall be selected in such way that it represents the local
shape of a reference ellipsoid. This is achieved by setting
\aRds{\odot} to the radius of curvature of the ellipsoid. The
curvature radius differs from the local radius except at the equator
and an east-west direction. For example, at the equator and a
north-south direction, the curvature radius is smaller then the local
radius, while at the poles (for all directions) it is greater (see
further Figure~\ref{fig:ppath:wgs84radii}). 

%The implications of
%approximating the geoid by the local curvature radius, instead of
%using a complete representation, is exemplified in
%Figure~\ref{fig:ppath:wgs84_dz}.

The \textindex{curvature radius}, \aRds{c}, of an ellipsoid is 
\citep{rodgers:00}
\begin{equation}
 \aRds{c} = \frac{1}{\aRds{ns}^{-1}\cos^2 \Lat + \aRds{ew}^{-1}\sin^2 \Lat}
 \label{eq:ppath:curvradius} 
\end{equation}
where \aRds{ns} and \aRds{ew} are the north-south and east-west curvature radius, respectively,
\begin{eqnarray}
 \aRds{ns} &=& \aRds{e}^2\aRds{p}^2 (
           \aRds{e}^2\cos^2\AzmAng+\aRds{p}^2\sin^2\AzmAng )^{-\frac{3}{2}} \\
 \aRds{ew} &=& \aRds{e}^2 (
           \aRds{e}^2\cos^2\AzmAng+\aRds{p}^2\sin^2\AzmAng )^{-\frac{1}{2}} 
 \label{eq:ppath:rew} 
\end{eqnarray}
The azimuth angle, \AzmAng, is defined in
Section~\ref{sec:fm_defs:los}. The latitude and azimuth angle to
apply in Equations \ref{eq:ppath:curvradius}~-~\ref{eq:ppath:rew}
shall rather be valid for a middle point of the propagation paths
(such as some tangent point), instead of the sensor position. 

\begin{figure}[!p]
 \begin{center}
  \begin{minipage}[c]{0.65\textwidth}
   \begin{center}
    \includegraphics*[width=0.9\hsize]{ppath/latitudes}
   \end{center}
  \end{minipage}%
  \begin{minipage}[c]{0.35\textwidth}
   \caption{Definition of the ellipsoid radii, \aRds{e} and \aRds{p}, 
     geocentric latitude, \Lat, and geodetic latitude, \Lat$^*$. The
     dotted line is the normal to the local tangent of the geoid
     ellipsoid. The zenith and nadir directions, and geometrical
     altitudes, are here defined to follow the solid line.}
   \label{fig:ppath:lats}
  \end{minipage}
 \end{center}
\end{figure}   

\begin{figure}[!p]
 \begin{minipage}[c]{0.65\textwidth}
 \includegraphics*[width=0.96\textwidth]{ppath/wgs84_radii}
 \end{minipage}%
 \begin{minipage}[c]{0.35\textwidth}
  \caption{The ellipsoid radius (\aRds{\odot}) and curvature radius (\aRds{c})
    for the
    WGS-84 reference ellipsoid. The curvature radii are valid for the
    north-south direction.}
  \label{fig:ppath:wgs84radii}
 \end{minipage}%
\end{figure}   
        
\begin{figure}[!p]
 \begin{minipage}[c]{0.65\textwidth}
 \includegraphics*[width=0.96\textwidth]{ppath/wgs84_latdiff}
 \end{minipage}%
 \begin{minipage}[c]{0.35\textwidth}
  \caption{The change of the WGS-84 ellipsoid radius for  1\degree\ 
            latitude differences.}
  \label{fig:ppath:latdiff}
 \end{minipage}%
\end{figure}   

%\begin{figure}[!t]
% \begin{center}
%  \includegraphics*[width=0.80\hsize]{Figs/ppath/wgs84_dz}
%  \caption{The altitude above the geoid as a function of latitude
%    when using the WGS84 reference ellipsiod (solid line) and when
%    using a spherical geoid (dashed line). The radius in the latter
%    case is set to the curvature radius of WGS84 at 45\degree, in the
%    direction (north-south) of the simulated measurement. The two
%    propagation paths share the same tangent point (defined by a
%    zenith angle of 90\degree) at latitude 45\degree, but the lowest
%    geometrical altitude is slightly shifted from that position with
%    an ellipsiodal geoid.}
%  \label{fig:ppath:wgs84_dz}  
% \end{center}
%\end{figure}
% This figure was produced by the Matlab function mkfigs_geoid


\subsection{Geocentric and geodetic latitudes}
%===================
\label{sec:ppath:geolat}

The fact that the geoid is an ellipsoid, instead of a sphere, opens up
for the two different definitions of the latitude. The
\textindex{geocentric latitude}, which is the the one used here, is the
angle between the equatorial plane and the vector from the coordinate
system center to the position of concern. The \textindex{geodetic
  latitude} is also defined with respect to the equatorial plane, but
the angle to the normal to the reference ellipsoid is considered here, as
shown in Figure~\ref{fig:ppath:lats}. It could be mentioned that a
geocentric latitude does not depend on the geoid ellipsoid used, while
the geodetic latitudes change if another reference ellipsoid is
selected. An approximative relationship between the geodetic
($\Lat^*$) and geocentric (\Lat) latitudes is \citep{montenbruck:00}
\begin{equation}
 \Lat^* = \Lat + f\,\sin(2\Lat)  
 \label{eq:ppath:lats}
\end{equation}
where $f$ is the flattening of the ellipse:
\begin{equation}
 f = \frac{\aRds{e}-\aRds{p}}{\aRds{e}}
 \label{eq:ppath:flattening}
\end{equation}
The value of $f$ for the Earth is about 1/298.26. This means that the
largest differences between \Lat\ and $\Lat^*$ are found at
mid-latitudes and the maximum value is about 12 arc-minutes.

The \textindex{zenith} and \textindex{nadir} directions shall normally be
defined to follow the normal to the reference ellipsoid, but, if
nothing else is mentioned, these directions are here treated to go
along the vector the center of the coordinate system, as indicated in
Figure~\ref{fig:ppath:lats}. This latter definition is preferred
as it results in that a propagation path in the zenith/nadir direction
can be described by a single latitude and longitude value. The
difference in geometrical altitude when using these two possible
definitions on the zenith direction is proportional to the deviation
between geocentric and geodetic latitude (Equation~\ref{eq:ppath:lats}).
For an altitude of 100~km around $\Lat=45\degree$, the difference is
about 350~m.


\subsection{Geodetic datums}
%===================
\label{sec:ppath:geodatums}

Table~\ref{tab:ppath:geodatums} gives the equatorial and polar radii
of the reference ellipsoid for the geodetic datums handled by ARTS.

\begin{table}[!h]
  \begin{center}
    \begin{tabular}{c c c c l}
     Datum & \aRds{e} & \aRds{p} & $1/f$ & Reference \vspace*{1mm} \\ 
     \hline 
     WGS-84 & 6378.137 km & \emph{6356.752 km} & 298.2572235 & {\small \citet{montenbruck:00}}  \rule{0mm}{5mm} \vspace*{1mm} \\
     \hline
    \end{tabular}
    \caption{Equatorial and polar radius of reference ellipsoids. Values 
      given as \emph{italic} are 
      derived by the other two values and Equation~\ref{eq:ppath:flattening}.}
    \label{tab:ppath:geodatums}
  \end{center}
\end{table}



\section{Control file examples}
%===================
\label{sec:ppath:cfile}


Some examples on how the geoid radius and the surface altitude can be set:
\begin{verbatim}

# Set the geoid to model WGS84 for 2D and 3D.
r_geoidWGS84 {
}

# For 1D, *lat_1d* and *meridian_angle_1d* must be specified 
# to define the position and direction for which the curvature 
# radius shall be extracted.
NumericSet( lat_1d ) {
  value = 45
}
NumericSet( meridian_angle_1d ) {
  value = 0
}
r_geoidWGS84 {
}

# Set the geoid to be spherical, with a radius of 6370 km
r_geoidSpherical {
  r = 6370e3
}

# Set the geoid to be spherical, with standard radius (defined 
# in arts.h)
r_geoidSpherical {
  r = -1
}

# Set a constant surface altitude of 1 km
MatrixSetTakingSizeFromMatrix( z_surface, r_geoid ) {
  value = 1e3
}

\end{verbatim}

\noindent
Different possibilities for the \artsstyle{ppath\_step\_agenda}:
\begin{verbatim}

# Select geometric calculations, with no length criterion 
# for path points.
AgendaSet( ppath_step_agenda ) {
  ppath_stepGeometric{
    lmax = -1
  }
}

# Consider refraction. The ray tracing step length is 2 km 
# and the length criterion for path points is 10 km.
AgendaSet( refr_index_agenda ) {
  refr_indexThayer{
  }
}
AgendaSet( ppath_step_agenda ) {
  ppath_stepRefractionEuler{
    lraytrace = 2e3
    lmax      = 10e3
  }
}
\end{verbatim}





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

% LocalWords:  ppath cc stepGeometric stepGeometricWithLmax ppathCalc pos los
% LocalWords:  ArrayOfGridPos geom ppc geomppath gridpos fd gridrange Eq Eqs
% LocalWords:  rodgers WGS montenbruck
