\chapter{Propagation paths}
 \label{sec:ppaththeory}


 \starthistory
 120227 & Created by splitting and revising the corresponding chapter 
          in \user\ (Patrick Eriksson).\\
 \stophistory


\graphicspath{{Figs/ppath/}}



\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
\funcindex{ppath\_calc}. The main use of \builtindoc{ppathCalc} is to
debug and test the path calculations, and that WSM should normally not
be part of the control file. 


\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, initialisation (handled by
\funcindex{ppath\_start\_stepping}), a repeated call of
\builtindoc{ppath\_step\_agenda} and putting data into the return
structure (\builtindoc{ppath}).

\begin{algorithm}
 \begin{algorithmic}
  \STATE{check consistency of function input}
  \STATE{call \funcindex{ppath\_start\_stepping} to set 
    \builtindoc{ppath\_step}}
  \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
  \ENDWHILE
  \STATE{initialise the WSV \builtindoc{ppath} to hold found number of 
         path points} 
 \end{algorithmic}
 \caption{Outline of the function \funcindex{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 \builtindoc{ppath\_step\_agenda},
which means that the practical starting point for the path calculations must be
determined. As the propagation path is followed in the backward direction, the
calculation starting points equals the end point of the path. 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 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 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 background is set to be ``cloud box
    interior''. 
  \end{itemize}
\item If the sensor is placed outside the model atmosphere
  \begin{itemize}
  \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\_stepRefractionBasic}. 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 \builtindoc{ppath\_stepGeometric} are
\funcindex{ppath\_step\_geom\_1d}, \funcindex{ppath\_step\_geom\_2d}
and \funcindex{ppath\_step\_geom\_3d}. See \shortcode{m\_ppath.cc} to
get the names of the sub-functions for other propagation path step
workspace methods. 

\begin{algorithm}
 \begin{algorithmic}
  \STATE{call \funcindex{ppath\_start\_2d}}
  \IF{\shortcode{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 \shortcode{ppath\_step\_geom\_2d} with temporary 
     \shortcode{Ppath} structure}
   \STATE{append temporary \shortcode{Ppath} structure to 
     \shortcode{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\_byltest}. Paths
steps passing a tangent point are handled by a recursive call of the
step function. Algorithm \ref{alg:ppath:ppath_step_geom_2d} summarises
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
% \builtindoc{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 \shortcode{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 \developer, Section \ref{D-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
% \shortcode{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 $\Rds\sin(\ZntAng)$ must 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 neighbouring
point, a numerical inaccuracy can accumulate when going from one point
to next. The propagation path constant is stored in the field
\shortcode{constant} of \wsvindex{ppath} and \wsvindex{ppath\_step}.

\begin{figure}
 \begin{center}
  \begin{minipage}[c]{0.65\textwidth}
   \begin{center}
    \includegraphics*[width=0.9\hsize]{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 centre 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$
fulfil 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 propagation 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 \builtindoc{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 fulfil a length criterion between the path
points. The mathematics of \shortcode{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 Chapter~\ref{U-sec:atmosphere} of \user\ results in
that for a 2D case the radius of a pressure surface varies linearly from one
point of the latitude grid to next. Compared to the 1D case, this is the main
additional problem to solve, handled by \funcindex{plevel\_crossing\_2d}. A two
step procedure is applied. In the first step the propagation path is moved
towards the pressure level as far as exact expressions can be used. For
example, if the level is approached from above the path is moved down to the
maximum radius of the level inside the gridbox. An approximative solution is
needed for the second step. Figure \ref{fig:ppath:psurf_crossing} gives a
schematic description of the problem at hand, which is handled by the internal
function \funcindex{rslope\_crossing}.
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]{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 approximate
solution is to note that \Lat\ is 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 distance between $\aRds{p}$ and $\Rds$ must be
about 8\,km. For angles $\Lat \leq 3\degree$, the sine and cosine terms can be
replaced with the three first (non-constant) terms of their Taylor expansions
maintaining a high accuracy. That is,
\begin{eqnarray}
  \cos \Lat & \approx & 1 - \Lat^2 / 2 + \Lat^4 / 24 + \Lat^6 / 720 \nonumber\\
  \sin \Lat & \approx & \Lat - \Lat^3 / 6 + \Lat^5 / 120 \nonumber
\end{eqnarray}
Equation \ref{eq:ppath:psurf1} becomes with these replacements a polynomial
equation of order 6:
\begin{eqnarray}
    0 & = & p_0 + p_1 \Lat + p_2 \Lat^2 + p_3 \Lat^3 + p_4 \Lat^4 + 
                  p_5 \Lat^5 + p_6 \Lat^6, \\
  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 & = &  \aRds{0}\sin(\Theta_p)/24  - c \cos(\Theta_p)/6,  \nonumber \\
  p_5 & = &  \aRds{0}\cos(\Theta_p)/120 + c \sin(\Theta_p)/24,  \nonumber \\
  p_6 & = & -\aRds{0}\sin(\Theta_p)/720 + c \cos(\Theta_p)/120.  \nonumber 
  \label{eq:ppath:psurf2}
\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. Several issues associated with numerical accuracy must be
considered, see the code (\funcindex{rslope\_crossing2d}) for details.

\begin{figure}
 \begin{center}
  \includegraphics*[width=0.80\hsize]{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 horizontal size deviates
    also from normal real cases.  Typical values for the vertical
    extension is around 500\,m, while the horizontal 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
\shortcode{plevel\_crossing\_2d} to calculate the latitude distance to a
crossing of the pressure surface below and above the present path point, as
well as the planets surface if it is found inside the grid box. 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}
\subsubsection{Conversion between polar and Cartesian coordinates}
%-
\label{sec:ppath:cart2sph}

The Cartesian coordinate system used follows the (standard?) Earth-centred
earth-fixed (ECEF) system (\url{http://en.wikipedia.org/wiki/ECEF}), with the
axes defined as:
\begin{description}
\item[x-axis] is along latitude 0\degree and longitude 0\degree
\item[y-axis] is along latitude 0\degree and longitude +90\degree
\item[z-axis] is along latitude +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 \cos(\Lat) \sin(\Lon) \\
  z &=& \Rds \sin(\Lat)            \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(z/\Rds)                \\
  \Lon &=& \arctan(y/x) \qquad \mathrm{(implemented\ by\ the\ atan2\ function)}
                               \nonumber
\end{eqnarray}
The functions performing these transformations are \funcindex{sph2cart} and
\funcindex{cart2sph}.

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) \\
  \cos(\Lat)\sin(\Lon) & -\Rds\sin(\Lat)\sin(\Lon) & \Rds\cos(\Lat)\cos(\Lon)\\ 
  \sin(\Lat)           & \Rds\cos(\Lat)            & 0                     
 \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 (see Section \ref{U-sec:fm_defs:los} of \user).

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) & \cos(\Lat)\sin(\Lon) & \sin(\Lat) \\ 
  -\sin(\Lat)\cos(\Lon)/\Rds & -\sin(\Lat)\sin(\Lon)/\Rds & \cos(\Lat)/\Rds \\ 
  -\sin(\Lon)/(\Rds\cos(\Lat)) & \cos(\Lon)/(\Rds\cos(\Lat)) & 0 
 \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}
These calculations are performed in \funcindex{cart2poslos}. Again special
expressions must be used for positions at the north and south pole (see the
code).


\subsubsection{Finding the crossing of a specified \Rds, \Lat\ or \Lon}
%-
\label{sec:ppath:3dcross}

The starting point in for all three cases is the following
equation system:
\begin{eqnarray}
  \Rds \cos(\Lat) \cos(\Lon)  & = & x+l\DiffD x, \nonumber \\
  \Rds \cos(\Lat) \sin(\Lon)  & = & y+l\DiffD y,           \\
  \Rds \sin(\Lat)             & = & 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 either \Rds, \Lat\ or \Lon\ is given.

The distance $l$ to a given \Rds\ 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 can easily be
calculated by Equations \ref{eq:ppath:xdl} and
\ref{eq:ppath:cart2sph}. These calculations are implemented in the
function \funcindex{r\_crossing\_3d}.  

If instead \Lat\ is given, the length to the point of interest can found by
again squaring the three equations, but now summing the x- and y-terms
and diving with the z-term:
\begin{equation}
  \tan^2(\Lat) = \frac{(z+l\DiffD z)^2}{(x+l\DiffD x)^2+(y+l\DiffD y)^2}.
\end{equation}
The solution of this quadratic equation is implemented in the function
\funcindex{lat\_crossing\_3d}\footnote{This function is presently not used. The
  source code can be found in \texttt{ppath\_NotUsed.cc}}. The solution for
$\Lat=0\degree$ is particularly simple ($l=-z/\DiffD z$). The case of
$\Lat=90\degree$ is set to have no solution ($\tan(90\degree)=\infty$), and is
instead assumed to be picked up as a crossing with one of the two longitudes
defining the grid box. Another complication is that, as the $\tan$-term is
squared, both $\pm\Lat$ can show up as possible solutions, and it must be
tested that the found length gives a \Lat\ with the correct sign.

For a given longitude, the x- and y-equations can be combined to give:
\begin{equation}
  l = \frac{y-x\tan(\Lon)}{\DiffD x\tan(\Lon)-\DiffD y}.
\end{equation}
This case is handled by \funcindex{lon\_crossing\_3d}\footnote{This function is
  presently not used. The source code can be found in
  \texttt{ppath\_NotUsed.cc}}. If the zenith or azimuth angle equals 0\degree\
or 180\degree, or if the start and target longitudes are equal, there is no
valid solution.



\subsubsection{Finding the crossing with a pressure level}
%-
\label{sec:ppath:3dplevel}

The same approach as for 2D is applied. The difference is that for 3D the
additional dimension gives a more complex variation of the radius of the
pressure level. For 2D, the variation can be expressed as a first order
polynomial ($r=\aRds{0}+c\Lat$), while for 3D a second order polynomial must be
used
\begin{equation}
  r = \aRds{0}+c_1\Lat+c_2\Lat^2.
\end{equation}
The coefficients $c_1$ and $c_2$ are detwermined in a purely numerical way, by
\funcindex{plevel\_slope\_3d}. The change in radius, $\Delta r_1$ and
$\Delta r_2$, at a distance of $\Delta\Lat$ and $2\Delta\Lat$, respectively, are
determined. These values give
\begin{equation}
  c_1 = \frac{4\Delta r_1 - \Delta r_2}{2\Delta\Lat}
\end{equation}
and
\begin{equation}
  c_2 = \frac{4\Delta r_1 - c_1\Delta\Lat}{(\Delta\Lat)^2}.
\end{equation}
The polynomial to solve becomes (cf.\ Eq.~\ref{eq:ppath:psurf2})
\begin{eqnarray}
    0 & = & p_0 + p_1 \Lat + p_2 \Lat^2 + p_3 \Lat^3 + p_4 \Lat^4 + 
                  p_5 \Lat^5 + p_6 \Lat^6, \\
  p_0 & = & (\aRds{0}-\aRds{p}) \sin(\Theta_p),                    \nonumber \\ 
  p_1 & = &   \aRds{0}\cos(\Theta_p)    + c_1 \sin(\Theta_p),      \nonumber \\ 
  p_2 & = & -\aRds{0}\sin(\Theta_p)/2   + c_1 \cos(\Theta_p)  
                                        + c_2 \sin(\Theta_p),      \nonumber \\ 
  p_3 & = & -\aRds{0}\cos(\Theta_p)/6   - c_1 \sin(\Theta_p)/2    
                                        + c_2 \cos(\Theta_p),      \nonumber \\ 
  p_4 & = &  \aRds{0}\sin(\Theta_p)/24  - c_1 \cos(\Theta_p)/6
                                        - c_2 \sin(\Theta_p)/2,    \nonumber \\
  p_5 & = &  \aRds{0}\cos(\Theta_p)/120 + c_1 \sin(\Theta_p)/24   
                                        - c_2 \cos(\Theta_p)/6,    \nonumber \\
  p_6 & = & -\aRds{0}\sin(\Theta_p)/720 + c_1 \cos(\Theta_p)/120  
                                        + c_2 \sin(\Theta_p)/24.   \nonumber
  \label{eq:ppath:psurf2}
\end{eqnarray}
The solution of this polynomial is handled by
\funcindex{rslope\_crossing3d}\footnote{This function is presently not used.
  The source code can be found in \texttt{ppath\_NotUsed.cc}}.



\subsubsection{A robust 3D algorithm}
%-
\label{sec:ppath:3dgeom}
%
\begin{algorithm}
 \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$} 
   \IF{$l_s>0$}
    \STATE{$l_{out} = l_s$}
    \COMMENT{$l_s$ is a function input}
   \ELSE
    \STATE{set $l_s$ to 3*vertical thickness of gid cell }
   \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}
  \STATE{$l_{out} \gets 5 * l_{out}$ }
  \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)$}
  \STATE{A recursive call can be neeeded, see the text.}
 \end{algorithmic}
 \caption{The method applied in \shortcode{do\_gridcell\_3d\_byltest} 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}
%
\begin{figure}[!b]
 \begin{center}
  \includegraphics*[width=0.80\hsize]{ppath_3Dsearch}
  \caption{Schematic of Algorithm \ref{alg:ppath:dogridcell3d}. The
    figure shows two iterations of the algorithm 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

\noindent
Some of the expressions presented above, for finding the crossing of a
specified \Rds, \Lat\ or \Lon, wehre found to be sensitive to numerical
inaccuracy and an algorith that avoids those expressions have been devised. It
applies a straightforward ``length-search'' algorithm (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 the 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. This algorithm is implemented by the function
\funcindex{do\_gridcell\_3d\_byltest}. 

The core task is to find the length of the path step. The search algorithm is
safe with respect to all grid cell boundaries, except the lower pressure level
where it can fail for zenith angles around 90\degree. In this case, the path
can pass the lower pressure level and re-enter the grid cell after a short
distance. For 1D cases, the part inside the lower cell would hold the tangent
point, but for a non-spherical reference ellipsoid and ``titled'' pressure
levels the tangent point can be found elsewhere. 

The bisection algorithm can miss such excursions to the lower grid cell.
Analytic approaches to handle this was rejected due to numerical problems.
Instead, all final points of the path step are checked with respect to this
issue and if any point is found to be below the lower pressure level, the
function is called recursively with the distance to the problematic point as
maximum search length ($l_s$ in Algorithm~\ref{alg:ppath:dogridcell3d}).


\section{Basic treatment of refraction}
%===================
\label{sec:ppath:refreuler}

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. 

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}[!p]
 \begin{center}
  \includegraphics*[width=0.70\hsize]{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 horizontal 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


\begin{figure}[!p]
 \begin{center}
  \includegraphics*{euler}
  \caption{Schematic of the ``basic'' ray tracing scheme. The ray tracing step 
    length is $l_r$. }
  \label{fig:ppath:euler}  
 \end{center}
\end{figure}

The workspace method \wsmindex{ppath\_stepRefractionBasic} takes
refraction into consideration by probably the most simple (from the
viewpoint of implementation) algorithm possible. 

The approach taken in \builtindoc{ppath\_stepRefractionBasic} 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. 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 generic input
\shortcode{lraytrace}. This length shall not be confused with the
final distance between the points that define the path, which is
controlled by \shortcode{lmax}. The path is first
determined in steps of \shortcode{lraytrace}. The normal situation
is that the ray tracing step length is considerably shorter than the
final spacing between the path points. Suitable values for
\shortcode{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}[tb!]
  \begin{center}
    \includegraphics*{snell}
    \caption{Geometry to derive Snell's law for a spherical atmosphere. }
    \label{fig:ppath:snell} 
  \end{center} 
\end{figure}

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 $p_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.

\begin{figure}[tb!]
 \begin{center}
  \includegraphics*[width=0.7\hsize]{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

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}[!p]
 \begin{center}
  \includegraphics*[width=0.65\hsize]{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}[!p]
 \begin{center}
  \includegraphics*[width=0.65\hsize]{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}
These calculations are handled by \funcindex{raytrace\_2d\_linear\_euler}.

The gradients of the refractive index for 2D are calculated by the function
\funcindex{refr\_gradients\_2d}. 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. This
function returns both gradients as the change of the refractive index over
1\,m. The conversion for the latitude gradient, from rad$^{-1}$ to m$^{-1}$,
corresponds to the $1/\Rds$ term found in Equation \ref{eq:ppath:za2d}, and
this term is accordingly left out in \shortcode{raytrace\_2d\_linear\_euler}.





\subsection{3D}
\label{sec:ppath:refr3D}
%----------------------
For 3D, the geometrical 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)$.










%%% 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
