\chapter{Reference ellipsoid and surface properties}
 \label{sec:surface}


\starthistory
  110614 & Extended and revised (Patrick Eriksson). \\
  050613 & First version finished by Patrick Eriksson. \\
\stophistory


\graphicspath{{Figs/ppath/}}


\section{The reference ellipsoid}
%===================
\label{sec:fm_defs:geoid}

The \textindex{ellipsoid} is an imaginary surface used as a reference when
specifying the surface altitude and the altitude of pressure levels. As the
name indicates, this reference surface is defined as an ellipsoid. The
reference ellipsoid should normally be set to some global geodetic datum, such
as WGS-84.
Inside ARTS, the ellipsoid is represented as a vector denoted as
\wsvindex{refellipsoid}. This vector must have length two, where the two
elements are equatorial radius (\aRds{e}) and eccentricity ($e$), respectively.




\subsection{Ellipsoid models}
%===================
\label{sec:ppath:geoid}

A geodetic datums is based on a reference 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
defintion of the ellispoid used in ARTS is based on \aRds{e} and the
eccentricity, $e$:
\begin{equation}
 e = \sqrt{1-\aRds{e}^2/\aRds{p}^2}.
 \label{eq:ppath:eccentricity} 
\end{equation}
The radius of the ellipsoid for a given (geocentric) latitude, \Lat, 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}} =
  \frac{\aRds{p}}{\sqrt{\sin^2\Lat+(1-e^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
ellipsoid 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{e} to the radius of curvature of
the ellipsoid and $e$ to zero. 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 \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}
 \begin{center}
  \begin{minipage}[c]{0.65\textwidth}
   \begin{center}
    \includegraphics*[width=0.9\hsize]{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
     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}
 \begin{minipage}[c]{0.65\textwidth}
 \includegraphics*[width=0.96\textwidth]{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}
 \begin{minipage}[c]{0.65\textwidth}
 \includegraphics*[width=0.96\textwidth]{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}
% \begin{center}
%  \includegraphics*[width=0.80\hsize]{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

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

% \begin{table}
%   \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}


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

The fact that Earth is rather an ellipsoid than 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 centre 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 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.






\section{Surface altitude}
%===================

The surface altitude, \aAlt{g}, is given as the geometrical altitude above the
reference ellispoid. The radius for the surface is accordingly
\begin{equation}
  \aRds{s} = \aRds{\odot} + \aAlt{s}
 \label{eq:fm_defs:zsurface}
\end{equation}
As also mentioned in Section~\ref{sec:fm_defs:surf}, a gap between the surface
and the lowermost pressure level is not allowed.

The ARTS variable for the \textindex{surface altitude} (\wsvindex{z\_surface})
is a matrix. For 1D, the surface constitutes a sphere by definition (as the
ellipsoid), while for 2D and 3D any shape is allowed and a rough model of the
surface topography can be made.




\section{Surface emission and reflection}
\label{sec:surf:eandr}
%===================

An introduction to the treatment of surface emission and reflections is given
in Section \ref{sec:fm_defs:surface}. As described in that section, the surface
properties are described by three workspace variables and the methods developed
to handle different surface properties set these variables:
\wsvindex{surface\_emission}, \wsvindex{surface\_los} and
\wsvindex{surface\_rmatrix}. The sections below outlines the methods available,
how these set the output variables.
Section~\ref{T-sec:surface} of \theory\ provides the theoretical
background.



\subsection{Blackbody surface}
%
If the surface can be assumed to act as a blackbody, the workspace method
\wsmindex{surfaceBlackbody} can be used. This method sets
\builtindoc{surface\_emission} to $[B,0,0,0]^T$, and \builtindoc{surface\_los}
and \builtindoc{surface\_rmatrix} to be empty.




\subsection{Specular reflections}
%
Several methods to incorporate a flat surface exist, including
\wsmindex{surfaceFlatRefractiveIndex} and
\wsmindex{surfaceFlatSingleEmissivity}. The methods differ in how the dielectric
properties of the surface are given, and if these are constant or not with
frequency.

In the case of specular reflections, \builtindoc{surface\_los} has the length 1.
The specular direction is calculated by the internal function
\funcindex{surface\_specular\_los}\footnote{Any tilt of the surface is
  neglected when determining the specular direction. If there would be any need
  to consider surface tilt, almost complete code for this task existed in
  \builtindoc{surface\_specular\_los} but was removed in version 1-1-876. The
  code can be obtained by e.g.\ checking out version 1-1-875.}. Equations
\ref{T-eq:surface:specular_matrix}-\ref{T-eq:surface:specular_emission}
in \theory give the values of \builtindoc{surface\_rmatrix} and
\builtindoc{surface\_emission}.



\subsection{Lambertian surface}
%---
A basic treatment of Lambertian surfaces is provided by the method
\wsmindex{surfaceLambertianSimple}. This method assumes that the down-welling
radiation has no azimuthal dependency, which fits the assumptions for 1D
atmospheres. The number of angles to apply in \builtindoc{surface\_los} is
selected by the user. 

For a Lambertian surface the reflected radiation is unpolarised (thus
independent of the polarisation of the down-welling radiation). That is,
each \builtindoc{surface\_rmatrix} has the structure:
\begin{equation}
  \MtrStl{R} =
     \left[\begin{array}{cccc}
       w&0&0&0\\
       0&0&0&0\\
       0&0&0&0\\
       0&0&0&0\\
     \end{array}
     \right].
\end{equation}
When determining the ``weight'' $w$ above, the method assumes that the
down-welling radiance ($I$) is constant inside each zenith angle range:
$[\theta_a,\theta_b]$. Hence, $w$ equals (cf.\
Equation~\ref{T-eq:surface:brdf1} of \theory)
\begin{equation}
  \label{eq:surface:brdf2}
  w = \int_{\theta_a}^{\theta_b} \! \int_{\phi_a}^{\phi_b} 
  \cos(\theta) f(\theta,\phi,\theta_1,\phi_1)
  \sin(\theta) \, \DiffD\phi \, \DiffD\theta.
\end{equation}
that gives
\begin{equation}
  w = \frac{r_d}{2}\left[\cos(2\theta_a)-\cos(2\theta_b)\right].
\end{equation}
Thus, this value is a combination of the surface reflectivity and an solid
angle weight.


The emission (\builtindoc{surface\_emission}) becomes:
\begin{equation}
  \label{eq:surface:blambertian} 
  \VctStl{b} =  \left[\begin{array}{c} r_dB \\ 0 \\0\\0 \end{array}\right].
\end{equation}



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