\chapter{Clear-sky Jacobians}
 \label{sec:wfuns}

 \starthistory
 110826 & First complete version by Patrick Eriksson.\\
 \stophistory

\graphicspath{{Figs/wfuns/}}


Inversions of both OEM and Tikhonov type require that the Jacobian can be
provided by the forward model \citep[see e.g.][]{eriksson:analy:00}. A
retrieval characterisation following \citet{rodgers:90,rodgers:00} raises the
same demand. A column of the Jacobian, \aWfnMtr{\SttVct}, is defined as
\begin{equation}
  \frac{\partial\MsrVct}{\partial\aSttElm{p}},
  \label{eq:wfuns:ki}
\end{equation}
where \MsrVct\ is the vector of measurement data and \aSttElm{p}\ is one forward
model (scalar) variable. See further Section~\ref{T-sec:formalism:wfuns} of
\theory. The nomenclature of that section is also used here.

The quantity in Eq.~\ref{eq:wfuns:ki} is in the atmospheric sounding
community frequently denoted as a ``weighting function'', and accordingly
\aWfnMtr{\SttVct} is called the weighting function matrix. In the documentation
of ARTS both terms (Jacobian and weighting functions) are used. These names
refer normally to \aWfnMtr{\SttVct}, the partial derivatives with respect to
the variables to be retrieved, forming the state vector \SttVct. However, in
the context of retrieval characterisation, the same matrix for the remaining
model parameters is of equally high interest, denoted as \aWfnMtr{\FrwMdlVct}.
In the same manner, the terms inversion and retrieval are used interchangeably.

The main task of the user is to select which quantities that shall be
retrieved, and to define the associated retrieval grids. These aspects must be
considered for successful inversions, but are out of scope for this document.
Beside for the most simple retrievals, it is further important to understand
how the Jacobian is calculated. A practical point is the calculation speed,
primarily determined if perturbations or analytical expressions are used
(Sec.~\ref{sec:wfuns:intro}). The derivation of the different Jacobians
involves some approximations due to theoretical and practical considerations.
Such approximations can be accepted, if of low or moderate size, but will
result in a slower convergence (the inversion will require more iterations).
Due to these later aspects, and to meet the needs of more experienced users,
this section is relatively detailed and contains a (high?) number of equations.

This section is restricted to Jacobians for clear-sky conditions, i.e.\ to be
applied outside the cloudbox. So far none of the scattering methods provide
Jacobians. Sections \ref{sec:wfuns:intro}\,-\,\ref{sec:wfuns:basis} contain
information of general character, while the available quantities are discussed
in the remaining sections (Section~\ref{sec:wfuns:absspecies} and forward).



\section{Introduction}
%==============================================================================
\label{sec:wfuns:intro}
%
There are two main approaches for calculating Jacobians, by analytical
expressions and by perturbations. We start with the conceptually simplest one,
but also the more inefficient approach.



\subsection{Perturbations}
%==============================================================================
\label{sec:wfuns:pert}
%
The most straightforward method to determine the Jacobian is by perturbing the
model parameter of concern. For example, the Jacobian corresponding to state
variable $p$ can always be calculated as
\begin{equation}
  \aWfnMtr{\SttVct}^p = \frac{\FrwMdl(\SttVct+\Delta\SttElm\VctStl{e}^p,
                      \FrwMdlVct)-\FrwMdl(\SttVct,\FrwMdlVct)} {\Delta\SttElm}
 \label{eq:wfuns:perturb}
\end{equation}
where $(\SttVct,\FrwMdlVct)$ is the linearization state, $\VctStl{e}^p$ is a
vector of zeros except for the $p$:th component that is unity, and
$\Delta\SttElm$ is a small disturbance (but sufficiently large to avoid
numerical instabilities).

However, it is normally not needed to make a recalculation using the total
forward model as the variables are either part of the atmospheric or the sensor
state, but not both. In addition, in many cases it is possibly to find
short-cuts. For example, the perturbed state can be approximated by an
interpolation of existing data (such as for a perturbed zenith angle). Such
short-cuts are discussed separately for each retrieval quantity.


\subsection{Analytical expressions}
%==============================================================================
\label{sec:wfuns:anal}
%
For most atmospheric variables, such as species abundance, it is possible to
derive an analytical expression for the Jacobians. This is advantageous because
they result in faster and more accurate calculations. Such expressions are
derived below. Some of the terms involved are calculated as a perturbation.
This is because some underlying functions of ARTS are best-fit methods that
are internally interpolated to the correct grid anyways. 
The analytical calculations are introduced in Sec.~\ref{sec:wfuns:atmvars}. 


\subsection{Workspace variables and methods}
%==============================================================================
\label{sec:wfuns:wsm}
%
As a workspace variable, the complete Jacobian is denoted as
\wsvindex{jacobian}. Auxiliary information is provided by
\wsvindex{jacobian\_quantities}. The actual calculations are made as part of
\builtindoc{yCalc}.

The retrieval quantities are defined separately, before calling
\builtindoc{yCalc}. This process is started by calling \wsmindex{jacobianInit}.
The retrieval quantities are then introduced through workspace methods named as
jacobianAdd{\it Something}. For example, for atmospheric temperature the method
is \builtindoc{jacobianAddTemperature}. It does not matter in which order these
``add methods'' are called.

The definition of retrieval quantities is finalised by calling
\wsmindex{jacobianClose}. To disable the calculation of the Jacobian, skip all
above, and just use \wsmindex{jacobianOff}. The methods named jacobianCalc{\it
  Something} shall never be used directly. Neither needs the user to consider
\wsaindex{jacobian\_agenda}.

The input to the ``add methods'' differs. In some cases you can select between
the analytical and perturbation options. For all perturbation calculations you
must specify the size of the perturbation. For atmospheric gases you can use
different units. For atmospheric fields, and some other quantities, you must
define the retrieval grid(s) to use.



 
\section{Basis functions}
%==============================================================================
\label{sec:wfuns:basis}

A forward model must use a discrete representation: it describes each quantity
with one or several variables. This is unproblematic for quantities that are of
discrete nature (including scalar variables). However, for atmospheric fields
and other continuous model quantities, the discrete representation inside the
forward model requires consideration. To avoid inconsistencies between model
input and output it is important that the mapping from the discrete variables
to the ``continuous view'' of the quantity is well defined, and applied
consistently through the forward model . This mapping is given by the basis
functions\index{basis function}. Similar arguments and nomenclature are found
in \citet{read:thecl:06}.

The basis functions are discussed explicitly in few places in this user guide,
but it shall be noted that all interpolations imply an underlying set of basis
functions. On the other hand, an understanding of both the derivation and the
obtained Jacobians require direct consideration of the basis functions. ARTS
operates with two types of basis functions.




\subsection{Basis functions for piece-wise linear quantities}
\label{sec:wfuns:basis1}
%
To treat an one-dimensional quantity to be piece-wise linear, or to say that a
linear interpolation shall be applied, are identical definitions. The basis
functions matching this definition have triangular shape, sometimes denoted as
``tent functions''. Such functions are exemplified in
Fig.~\ref{fig:wfuns:zbasis}, see also \citet{buehler:artst:05}.

\begin{figure}[t]
 \begin{center}
  \includegraphics*[width=0.7\hsize]{fig_absbasis_z}
  \caption{Examples on 1d basis functions for a vertical grid with a 1 km
           spacing: \lsolid~30~km, \ldashed~31~km and \ldashdot~32~km.}
  \label{fig:wfuns:zbasis}  
 \end{center}
\end{figure}

In most cases, the quantity is considered to be undefined outside the end
points of the grid. Hence, the basis function for a grid end point is then just
``half a tent''. The exception to this rule is retrieval grids of piece-wise
linear variables. To avoid that retrieval grids must cover the complete
atmosphere, end point values are assumed to be valid to the end of the
atmosphere (or data range of concern). That is, the basis functions for end
points of retrieval grids follow the tent shape inside the grid range, and
have a constant value of 1 outside. In terms of interpolation, this matches to
allow extrapolation, the applying a ``nearest'' interpolation for positions
outside the covered range (the end values are valid all the way to
$\pm$infinity).
The basis functions are defined likewise for higher dimensions, but the
tent functions are then 2D or 3D ``tents''.




\subsection{Polynomial basis functions}
\label{sec:wfuns:basis2}
%
Some retrieval quantities are expressed using a polynomial basis. Sensor zenith
angle pointing off-set is one such quantity. The off-set is then treated to
have a polynomial variation as a function of time. If the offset is assumed to
be constant in time, a zero order polynomial shall be selected. If there is
also a linear drift with time, use a first order polynomial, etc.

For these basis functions, the explanatory variable (time in the example above)
is normalised to cover the range [-1,1], here denoted as $z$, and the
continuous representation ($f$) of the variable of concern can be written as
\begin{equation}
  f(z) = \aSttElm{0} + \aSttElm{1}(z-b_1) + \aSttElm{2}(z^2-b_2) + 
                     \aSttElm{3}(z^3-b_3) + \aSttElm{4}(z^3-b_4) + \dots  
\end{equation}
where $\aSttElm{0}, \aSttElm{1}, \dots$ are the coefficients to be retrieved
(elements of \SttVct). The interpretation of a retrieval is simplified if the
average of $f$ equals \aSttElm{p0}, and the scalars $b_1, b_2, \dots$ are
selected, schematically, as
\begin{equation}
  0 = \int_{-1}^1 \!\left(z^n-b_n\right) \, \DiffD z, \quad n>0.
\end{equation}
According to this expression, $b_n$ is zero for odd $n$. However, $z$ is in
practise a discrete variable ($z_i$, not necessarily symmetric around 0), and
$b_n$ is taken as the average of $z_i^n$: all $b_n$ can be non-zero. The
normalisation of $z$ is not only made for interpretation reasons, it can be
required for pure numerical reasons, such as when $z$ represent frequency (in
Hz).

In practise, the basis functions are vectors, denoted below as $\VctStl{z}_i$.
Element $j$ of $\VctStl{z}_i$ is
\begin{equation}
  \VctStl{z}_i(j) = z(j)^i - b_i.
\end{equation}




\section{Atmospheric variables, common expressions}
%==============================================================================
\label{sec:wfuns:atmvars}
%
The analytically-oriented calculation procedure to obtain the Jacobian for
atmospheric quantities is here outlined. The expressions are based on the
chain rule, and can be applied for absorption constituents, atmospheric
temperatures and winds. It is important to notice that only local effects are
considered, and the expressions have limitations, as discussed below in
Sec.~\ref{sec:wfuns:atmvars:limit}.


\subsection{Matrix derivatives}
\label{sec:wfuns:atmvars:matder}
%
Matrix exponents are calculated using the Pad\'{e} approximation with scaling \& squaring.
Likewise, matrix exponent derivatives are calculated using the same method.
The implementation in ARTS follows \citet{brancik:06}.  As a brief reminder, following
the same notation as \citet{brancik:06}, but with
altered input/output to match the rest of this text, the transmission
matrix is
\begin{equation}
  \MtrStl{T}_i(x_i) = \exp\left(-\MtrStl{K}_i(x_i)r\right) \approx \MtrStl{D}_{pq}^{-1}(x_i)\MtrStl{N}_{pq}(x_i),\label{eq: pade}
\end{equation}
with
\begin{equation}
  \MtrStl{N}_{pq}(x_i)=\sum_{j=0}^p\frac{\left(p+q-j\right)!}{\left(p+q\right)!j!\left(p-j\right)}\left(-\MtrStl{K}_i(x_i)r\right)^j,
\end{equation}
and
\begin{equation}
  \MtrStl{D}_{pq}(x_i)=\sum_{j=0}^q\frac{\left(p+q-j\right)!}{\left(p+q\right)!j!\left(p-j\right)}\left(\MtrStl{K}_i(x_i)r\right)^j.
\end{equation}
This means that the matrix derivation of the transmission of the Pad\'{e} approximation method is
\begin{equation}
  \frac{\PartD \MtrStl{T}_i}{\PartD x_i} \approx \frac{\PartD  \MtrStl{D}_{pq}^{-1}}{\PartD x_i}\MtrStl{N}_{pq} + \MtrStl{D}_{pq}^{-1}\frac{\PartD  \MtrStl{N}_{pq}}{\PartD x_i},
\end{equation}
where the input arguments have been dropped (and will not appear again in this subsection) to make the equation(s) easier to read.
To be explicit, the above can be rewritten as
\begin{equation}
  \frac{\PartD \MtrStl{T}_i}{\PartD x_i} \approx \MtrStl{D}_{pq}^{-1}\left(\frac{\PartD  \MtrStl{N}_{pq}}{\PartD x_i} - \frac{\PartD \MtrStl{D}_{pq}}{\PartD x_i} \MtrStl{T}_i\right) \label{eq: pade derivative}
\end{equation}
in calculations.
Finally, the internal partial derivatives are just
\begin{equation}
  \frac{\PartD \MtrStl{N}_{pq}}{\PartD x_i}=\sum_{j=0}^p\frac{\left(p+q-j\right)!}{\left(p+q\right)!j!\left(p-j\right)}\left(-\frac{\PartD \MtrStl{K}_i}{\PartD x_i}r\right)^j,\label{eq: pade dn/dx}
\end{equation}
\begin{equation}
  \frac{\PartD \MtrStl{D}_{pq}}{\PartD x_i}=\sum_{j=0}^q\frac{\left(p+q-j\right)!}{\left(p+q\right)!j!\left(p-j\right)}\left(\frac{\PartD \MtrStl{K}_i}{\PartD x_i}r\right)^j.\label{eq: pade dd/dx}
\end{equation}

Important to note in ARTS is that we average layer properties from properties of the surrounding path points simply by
averaging the property at the two path points.
That is
\begin{equation}
   \MtrStl{K}_i = \frac{ \MtrStl{K}_{i-1/2}+\MtrStl{K}_{i+1/2}}{2},
\end{equation}
where the half denotes the propagation matrix at the two levels.  
So for every layer, we need the partial derivation of the layer for both the influence
by the level above and the level below. Thus, 
\begin{equation}
  \frac{\PartD \MtrStl{K}_i}{\PartD x_i} = \frac{1}{2}\left(\frac{\PartD \MtrStl{K}_{i-1/2}}{\PartD x_i}+\frac{\PartD \MtrStl{K}_{i+1/2}}{\PartD x_i}\right)
\end{equation}
These calculations are done in parallel to save memory and/or time.  In fact, only
Equations~\ref{eq: pade dn/dx}~and~\ref{eq: pade dd/dx} are recalculated for each analytical Jacobian entity that has been requested
by the user.

\subsection{Analytical expression for partial derivation of the propagation matrix}

All lines and contributions to the propagation matrix are summed for a total
propagation.  This means that
\begin{equation}
  \MtrStl{K}_i = \sum_j \MtrStl{K}_j;\;\;\;\frac{\PartD \MtrStl{K}_i}{\PartD x_i} = \sum_j \frac{\PartD \MtrStl{K}_j}{\PartD x_i},
  \label{eq: general dk/dx}
\end{equation}
for all $j$ lines, continua contributions, collision induced absorption contributions, etc..
This means that to calculate each individual ${\PartD \MtrStl{K}_j}/{\PartD x_i}$ is enough to get 
${\PartD \MtrStl{K}}/{\PartD x_i}$.  The method of calculating ${\PartD \MtrStl{K}_j}/{\PartD x_i}$
depends on how the absorption is calculated.  


\subsubsection{Parameterized absorption models}

Several methods in ARTS are best-fit or parameterized models.  If $\MtrStl{K}_j$ is 
from one of these functions, then a perturbation approach is used in a
low-level function to get ${\PartD \MtrStl{K}_j}/{\PartD x_i}$.  The user sets a 
small perturbation that is then used internally.  
This way is used when continua, collision induced absorption and lookup table 
is used since these are only available in parameterized form in ARTS.

\subsubsection{Line-by-line absorption}

Expanding $\MtrStl{K}_j$ for line-by-line calculations,
it can roughly be written
\begin{equation}
  \MtrStl{K}_j = S_jF_j\MtrStl{\Phi}_j, \label{eq: simple propagation}
\end{equation}
where $S_j$ is the scaled line strength, $F_j$ is the shape of the line, 
and $\MtrStl{\Phi}_j$ is the polarization matrix.
This rough expression is useful for defining the partial derivatives since
we can write
\begin{equation}
  \frac{\PartD \MtrStl{K}_j}{\PartD x_i} = 
  \frac{\PartD S_j}{\PartD x_i}F_j\MtrStl{\Phi}_j+
  S_j\frac{\PartD F_j}{\PartD x_i}\MtrStl{\Phi}_j+
  S_jF_j\frac{\PartD \MtrStl{\Phi}_j}{\PartD x_i}.
  \label{eq: line dk/dx}
\end{equation}
An important thing to remember about the equation above is that
everything in Equation~\ref{eq: simple propagation} is  already 
known (to the level of certainty in the model).  Thus, 
Equation~\ref{eq: line dk/dx} can, when it is convenient, be written as
\begin{equation}
  \frac{\PartD \MtrStl{K}_j}{\PartD x_i} = \MtrStl{K}_j
  \left(
  S_j^{-1}\frac{\PartD S_j}{\PartD x_i}+
  F_j^{-1}\frac{\PartD F_j}{\PartD x_i}+
  \MtrStl{\Phi}_j^{-1}\frac{\PartD \MtrStl{\Phi}_j}{\PartD x_i}
  \right). \label{eq: line dk/dx special}
\end{equation}
Both Equation~\ref{eq: line dk/dx} or \ref{eq: line dk/dx special} are
used internally in the code base.

To extend the rough expression a little bit, the line strength 
variable can be thought of as
\begin{equation}
  S_j = n_j(T,p) S_{j,0} f(T,\cdots),\label{eq: rough line strength}
\end{equation}
where $n_j$ is the number density, $S_{j,0}$ is the line strength 
at a reference temperature in LTE, and
$f(T,\cdots)$ is a function that mostly depends on the temperature
(but also other line parameters and non-LTE effects are included here).
The line shape variable can be thought of as
\begin{equation}
  F_j = F_n F_m F_u F_s,\label{eq: rough line shape}
\end{equation}
where $F_n$ is a renormalization term far from the line center (e.g., van Vleck and Huber renormalization),
$F_m$ is the line mixing,
$F_u$ is a unit conversion term, and
$F_s$ is the classical definition of line shape (e.g., Lorentz or Voigt distributions).
The polarization matrix is the unit matrix except for the Zeeman effect and for Faraday rotation.
In scalar simulations, it is simply the $1\times1$ unit matrix, or a scalar 1.

\subsection{Separation of terms}
\label{sec:wfuns:atmvars:sepe}
%
The overall task is to calculate (Eq.~\ref{eq:wfuns:ki} with subscript $p$ left
out)
\begin{displaymath}
  \frac{\PartD\MsrVct}{\PartD\SttElm}  
\end{displaymath}
where \SttElm\ is the element of $\SttVct$ for which we want to obtain the
weighting function. This is a column of the complete weighting function matrix.

A first step is to identify how sensor characteristics can be incorporated? To
make the nomenclature simpler, we assume here that the simulations cover only a
single measurement block and we have (cf.\ Eq.~\ref{eq:fm_defs:measseq})
\begin{equation}
  \MsrVct = \SnsMtr\aMpiVct{b}.
  \label{eq:wfuns:Hi}
\end{equation}
We can then apply the chain-rule for a first time to obtain 
\begin{equation}
  \frac{\PartD\MsrVct}{\PartD\SttElm} = 
  \frac{\PartD\MsrVct}{\PartD\aMpiVct{b}} 
  \frac{\PartD\aMpiVct{b}}{\PartD\SttElm} = 
  \SnsMtr \frac{\PartD\aMpiVct{b}}{\PartD\SttElm} = \SnsMtr\VctStl{k}_\MpiVct.
  \label{eq:wfuns:Hk}
\end{equation}
Hence, sensor characteristics can be handled by calculating the weighting
function column matching all monochromatic pencil beam calculation of the
measurement block ($\VctStl{k}_\MpiVct$) and perform a multiplication with
\SnsMtr, a parallel procedure compared to how \aMpiVct{b}\ is compiled to obtain
\MsrVct. The calculation procedure expands to allow that the complete
weighting function matrix (for quantities covered by the analytical calculation
procedure) is calculated as
\begin{equation}
  \frac{\PartD\MsrVct}{\PartD\SttVct} = \SnsMtr\aWfnMtr{\MpiVct},
  \label{eq:wfuns:HK}
\end{equation}
where $\aWfnMtr{\MpiVct} = \PartD\aMpiVct{b}/\PartD\SttVct$. One column of this
matrix is $\VctStl{k}_\MpiVct$.

The vector \aMpiVct{b}\ consists of a number of Stokes vectors appended,
$\aMpiVct{b}=[\aStoVec{1}^T, \aStoVec{2}^T,\dots,\aStoVec{n}^T]^T$, and the
calculation of $\VctStl{k}_\MpiVct$ can schematically be written as
\begin{equation}
  \VctStl{k}_\MpiVct = \sum_{j=1}^n \frac{\PartD\aMpiVct{b}}{\PartD\aStoVec{j}}
    \frac{\PartD\aStoVec{j}}{\PartD\SttElm}
  \label{eq:wfuns:kisum}
\end{equation}
The terms $\PartD\MpiVct/\PartD\aStoVec{j}$ are formally matrices. However,
these matrices are not calculated explicitly as they only contain information
on where \aStoVec{i}\ is stored inside \MpiVct. That is, these matrices are of
bookkeeping character, consisting only of zeros and ones. In practice, results
matching $\PartD\aStoVec{j}/\PartD\SttElm$ are simply inserted in correct place
of $\VctStl{k}_\MpiVct$, mimicking how is \aStoVec{i}\ put into \MpiVct.

Accordingly, the core task is to calculate $\PartD\aStoVec{j}/\PartD\SttElm$,
where for simplicity the subscript $j$ is dropped below. This calculation is
expanded as
\begin{equation}
  \frac{\PartD\StoVec}{\PartD\SttElm} = \sum_{i=1}^n
    \frac{\PartD\StoVec}{\PartD\aSttElm{i}}
    \frac{\PartD\aSttElm{i}}{\PartD\SttElm},
\end{equation}
where \aSttElm{i}\ is the value of the quantity of concern, at each point of
the propagation path. That is, $i$ indexes the path points. The actual
radiative transfer enters by the terms $\PartD\StoVec/\PartD\aSttElm{i}$, and
they are discussed separately below.

The term $\PartD\aSttElm{i}/\PartD\SttElm$ appears due to the fact that
\aSttElm{i}\ and \SttElm\ are placed at different positions, and the
representation of the atmospheric fields must be considered here. In practise,
the term is calculated as the value of the basis function for \SttElm\ at the
location of \aSttElm{i}\ (further discussed in \citet{buehler:artst:05}). This
is a slight approximation with respect to the goal of fully incorporating the
piece-wise linear representation in the weighting functions
\citep{buehler:artst:05}. A low value of $\Delta\aPpathLng{i}$ decreases the
degree of approximation.

It can be noted that $\PartD\aSttElm{i} / \PartD\SttElm$ is normally non-zero
for more than one element of \SttVct. The exception is if the positions of
\aSttElm{i}\ and \SttElm\ are identical. Reversely, the weighting function for
element \SttElm\ can have contributions from several propagation path points
(\aSttElm{i}), as well as from several radiance spectra. Accordingly, the
practical calculations are done by first determine all
$\PartD\StoVec/\PartD\aSttElm{i}$, of the given propagation path. These data
are then used to determine $\PartD\StoVec/\PartD\SttVct$, for the retrieval
quantity of concern. That is, each $\PartD\StoVec/\PartD\aSttElm{i}$ is
combined with $\PartD\aSttElm{i}/\PartD\SttElm$ for all elements of \SttVct.
Most of these combinations yield a zero result. The terms
$\PartD\aSttElm{i}/\PartD\SttElm$ are determined with help of ARTS' internal
interpolation grid position routines.



\subsection{$\PartD\StoVec/\PartD\aSttElm{i}$, general 
case}
\label{sec:wfuns:atmvars:gene}
%
The term $\PartD\StoVec/\PartD\aSttElm{i}$ is here outlined for the general case
of vector radiative transfer. In this case, the Stokes elements can not be
treated separately, and matrix-vector notation is required. The final Stokes
vector obtained through Eq.~\ref{eq:fm_defs:vrte_step} can be expressed as
\begin{equation}
  \label{eq:s4J}
  \StoVec = \StoVec' + \TrnMat' \left[ \aTrnMat{i} \aStoVec{i} +
         (\IdnMtr-\aTrnMat{i}) \left(\bar{\EmsVec}_i + \bar{\ExtMat}_i^{-1}\bar{\aSrcVec{n,i}}\right) \right],
\end{equation}
where $\StoVec'$ is the Stokes vector for all emission generated between the
sensor and point $i+1$, $\TrnMat'$ is the transmission Mueller matrix for the
same part of the propagation path and \aTrnMat{i} is defined in
Eq.~\ref{eq:rte:transmat}. The quantities $\StoVec'$, $\TrnMat'$ and
\aStoVec{i}\ are not function of \aSttElm{i}. This gives
\begin{eqnarray}
  \frac{\PartD\StoVec}{\PartD\aSttElm{i}} &=& \TrnMat'
  \left[
    \frac{\PartD\aTrnMat{i}}{\PartD\aSttElm{i}}
      \left( \aStoVec{i}-\bar{\EmsVec}_i -\bar{\ExtMat}_i^{-1}\bar{\aSrcVec{n,i}}\right)+\right. \nonumber\\ &&\left.
    \left( \IdnMtr-\aTrnMat{i} \right) \left(
      \frac{\PartD\bar{\EmsVec}_i}{\PartD\aSttElm{i}}+\bar{\ExtMat}_i^{-1}\frac{\PartD\bar{\ExtMat}_i}{\PartD\aSttElm{i}}\bar{\ExtMat}_i^{-1}\bar{\aSrcVec{n,i}} + \bar{\ExtMat}_i^{-1}\frac{\PartD\bar{\aSrcVec{n,i}}}{\PartD\aSttElm{i}} \right)
  \right].
  \label{eq:wfuns:dsdxi1}
\end{eqnarray}
The terms $\PartD\aTrnMat{i}/\PartD\aSttElm{i}$ and
$\PartD\bar{\EmsVec}_i/\PartD\aSttElm{i}$ are both calculated in a pure
numerical manner, for different reasons discussed below.


\subsection{Including the surface}
%
For scattered down-welling radiation the effective transmission matrix is
\begin{equation}
  \TrnMat' = \aTrnMat{2}\MtrStl{R}_i\aTrnMat{1},
  \label{eq:wfuns:surface}
\end{equation}
where \aTrnMat{2}\ is the transmission between the surface and the sensor,
$\MtrStl{R}_i$ is defined in Eq.~\ref{eq:fm_defs:surfacerefl} and
\aTrnMat{1}\ is the transmission between the point $i+1$ and the surface.



\subsection{$\PartD\StoVec/\PartD\aSttElm{i}$, locally unpolarized absorption}
\label{sec:wfuns:atmvars:unpol}
%
Eq.~\ref{eq:wfuns:dsdxi1} can be simplified for some conditions. A first case
is scalar radiative transfer, i.e. only the first element of the Stokes vector
is considered and all terms of Eq.~\ref{eq:wfuns:dsdxi1} are scalar quantities.
For all other cases, in principle vector radiative transfer must be performed,
as $\TrnMat'$ can always have off-diagonal elements. Even if atmospheric
absorption is totally unpolarized, $\TrnMat'$ can be non-diagonal due to the
surface (Eq.~\ref{eq:wfuns:surface}). However, if the absorption locally lacks
polarization, the calculations can be handled by analytical expressions in a
higher degree.

To focus on the analytical expressions dealing with local radiative transfer
effects, let us write Eq.~\ref{eq:s4J} as
\begin{equation}
  \StoVec = \StoVec' + \TrnMat'\aStoVec{s}.
\end{equation}
In the nomenclature of Eq.~\ref{eq:fm_defs:vrte_step},
$\aStoVec{s}=\aStoVec{i+1}$ (but would be confusing to here use \aStoVec{i+1}).
Hence
\begin{displaymath}
  \frac{\PartD\StoVec}{\PartD\aSttElm{i}} = \TrnMat'
    \frac{\PartD\aStoVec{s}}{\PartD\aSttElm{i}}.
\end{displaymath}
The first element of \aStoVec{s}, \StoI, can be written as (cf.\
Eq.~\ref{eq:fm_defs:rte_step})
\begin{equation}
  \label{eq:I4J}
  \StoI = e^{-\aOth{i}}\aStoI{i} + (1-e^{-\aOth{i}})\left(\bar{\Planck}_i+{\bar{j_{n,i}}}/{\bar{\AbsCoef}}\right).
\end{equation}
By again using the chain rule, the derivative of \StoI\ with respect to
\aSttElm{i}\ can be written as
\begin{equation}
  \frac{\PartD\StoI}{\PartD\aSttElm{i}} =
  \frac{\PartD\StoI}{\PartD\aOth{i}}\frac{\PartD\aOth{i}}{\PartD\aSttElm{i}} + 
  \frac{\PartD\StoI}{\PartD j_i}
  \frac{\PartD j_i}{\PartD\aSttElm{i}} = 
   e^{-\aOth{i}}\frac{\PartD\aOth{i}}{\PartD\aSttElm{i}}
      (j_i-\aStoI{i}) +
  (1-e^{-\aOth{i}})\frac{\PartD j_i}{\PartD\aSttElm{i}},
  \label{eq:wfuns:dsdxi2}
\end{equation}
where $j_i$ is short for $\bar{\Planck}_i+{\bar{j_{n,i}}}/{\bar{\AbsCoef}}$.
The two remaining derivatives $\PartD\aOth{i}/\PartD\aSttElm{i}$ and
$\PartD j_i/\PartD\aSttElm{i}$
depend on the quantity considered, and are discussed further below.

For higher Stokes components of \aStoVec{s}, or in general if emission is
totally ignored, Eq.~\ref{eq:I4J} is simplified to (here exemplified for the
second Stokes element, \StoQ)
\begin{equation}
  \label{eq:Q4J}
  \StoQ = e^{-\aOth{i}}\aStoQ{i},
\end{equation}
and the chain rule expression is correspondingly shorter:
\begin{equation}
  \frac{\PartD\StoQ}{\PartD\aSttElm{i}} =
   -e^{-\aOth{i}}
   \frac{\PartD\aOth{i}}{\PartD\aSttElm{i}}\aStoQ{i}.
  \label{eq:wfuns:dsdxi3}
\end{equation}


\subsection{Limitations}
\label{sec:wfuns:atmvars:limit}
%
A constraint for the analytical expressions above is that the effect of the
variable must only be local. Main examples on non-local effects should occur
through hydrostatic equilibrium and refraction. Significant impact of a gas
through these mechanisms should only be found for water vapour in the lower
troposphere, but is a general concern for temperature as discussed in
Sec.~\ref{sec:wfuns:atmtemp}.





\section{Absorption species}
%==============================================================================
\label{sec:wfuns:absspecies}

\subsection{Common practicalities}
%
To obtain the Jacobian for absorption species, use
\wsmindex{jacobianAddAbsSpecies}. The method handles one species at the time.
The calculations can either be done in an ``analytical'' or ``perturbation''
manner. For gases, weighting functions (WFs) can be provided for several units
of the gas abundance:
\begin{description}
\item[vmr] Volume mixing ratio (a value between 0 and 1). The WF divided by
  $10^6$ corresponds to that 1 ppm of the gas is added to the atmospheric
  volume of concern.
\item[nd] Number density. The WF corresponds here to that one molecule is added.
\item[rel] Relative/fractional change. In a perfectly linear case, the WF
  corresponds here to that the gas amount is doubled.
%\item[logrel] This option returns the same WFs as ``rel'', but is included to
%  flag that the natural logarithm of the ``rel'' unit is retrieved.
\end{description}
For the ``rel'' option it is important to note that ARTS
calculate the WFs with respect to the given state, ARTS does not know anything
about the actual reference state for which the ``rel'' unit is valid (where
normally the a priori state is selected). For iterative inversions, a rescaling
of the WFs provided by ARTS is likely needed, to make the WFs valid with
respect to the (original) reference state. 
%For the assumption made inside ARTS,
%the WFs for ``rel'' and ``logrel'' are identical.

A second main consideration is to select the retrieval grids. For analytical
calculations there are no other selections to be made. 


\subsection{Perturbation calculations}
%
For pure numerical calculations, also the size of the perturbation must be
specified ($\Delta\SttElm$ in Eq.~\ref{eq:wfuns:perturb}). The perturbation
shall given following the unit selected. The same value is applied for all WFs
(which can cause practical problems \dots).


\subsection{Analytical expressions}
%
If not made clear above, the only term that differs between the Jacobian
quantities is $\PartD\StoVec/\PartD\aSttElm{i}$.

All parameterized models in ARTS are parameterized on basis of cross-section.
These models' contribution to Equation\ref{eq: general dk/dx} are simply going to be
\begin{equation}
 \frac{\PartD \MtrStl{K}_j}{\PartD x_i}=\frac{1}{x_i} \MtrStl{K}_j,
\end{equation}
for the parameterized model contributions $\MtrStl{K}_j$ that are functions of $x_i$.
One thing of importance to note here is that in collision induced absorption, the 
absorption caused by the secondary species is also considered in the above expression.
 
For line-by-line contributions in Equation~\ref{eq: line dk/dx},
we only consider $\PartD S_j/\PartD x_i$ as important (when the line is from the absorption species).
The same expression as for parameterized models still apply for the partial derivative.

Note that we ignore that $\PartD F_j/\PartD x_i$ is mathematically non-trivial since pressure broadening 
and line mixing is changed if the atmosphere is changed.  In fact, $\PartD F_j/\PartD x_i$ is 
mathematically non-trivial even for lines other than those from the absorption species itself
since it changes the total pressure.  We justify ignoring this effect because the term is very 
small in most reasonable atmospheric setups.

\section{Winds}
%==============================================================================
\label{sec:wfuns:winds}
%
Calculation of wind weighting functions are triggered by
\wsmindex{jacobianAddWind}. Each wind component (see Sec.~\ref{sec:winds:defs})
is treated as an individual retrieval variable. That is, if you want to
retrieve all three wind components, \builtindoc{jacobianAddWind} must be called
three times. If you want the total wind Jacobian, this can also be attained.
Only the analytically inclined calculation approach is at hand for winds.

Theoretically, the Doppler shift induced by winds affects the emission source
term, but this impact is extremely small
and the related terms are ignored. This gives a case basically identical to the
one above for absorption species.

By design, the absorption calculations of ARTS are unaware of wind velocities.  
Therefore, only the frequency derivative is calculated at the lowest level.  That is,
for each component, this equation is evaluated
\begin{equation}
   \frac{\PartD \MtrStl{K}_j}{\PartD x_i}=\frac{\PartD \MtrStl{K}_j}{\PartD \nu} \frac{\PartD \nu}{\PartD x_i}.
\end{equation}

The frequency partial derivative term is
\begin{equation}
  \frac{\PartD \nu}{\PartD x_i} = \frac{1}{c}\frac{\PartD v}{\PartD x_i},
\end{equation}
where $c$ the speed of light and $v$ is the wind velocity relative to the observational path of the component of interest.
The last partial derivative is there due to the influence of the angle of observation on the shift in frequency.
Look along the direction of the wind component and it is the only component that matters so 
the partial derivative above is unity.  The other two components are therefore perpendicular
to the observational path and have partial derivatives that are nil.  We will not write the full expression since it
depends on atmospheric dimensionality and is different for each of the three components so there are a total of 
ten different expressions depending on user input.

For Equation~\ref{eq: line dk/dx}, only ${\PartD F_j}/{\PartD \nu}$ is non-zero.
The contribution is
\begin{equation}
   \frac{\PartD \MtrStl{K}_j}{\PartD \nu} = \MtrStl{K}_j \left(\frac{1}{F_n}\frac{\PartD F_n}{\PartD \nu}+\frac{1}{F_s}\frac{\PartD F_s}{\PartD \nu}\right).
   \label{eq: wind propagation derivative}
\end{equation}
For efficiency, the second term is calculated slightly differently due to efficiency.

There is a way to parameterize classical line shapes by two terms and write
\begin{equation}
  F_s = w(a+ib) = F_A + iF_B,
\end{equation}
where $a$ and $b$ are parameters used to determine how much pressure and temperature influence the line shape 
(see details in the theory guide),
$F_A$ is the Voigt line shape and $F_B$ is the Faraday-Voigt line shape.
This is usually called the Faddeeva line shape.
Using the Faddeeva line shape,
\begin{equation}
  \frac{\PartD w}{\PartD x_i} = \frac{\PartD w}{\PartD a}\frac{\PartD a}{\PartD x_i}+i\frac{\PartD w}{\PartD b}\frac{\PartD b}{\PartD x_i},
  \label{eq: line shape derivative}
\end{equation}
and a great convenience is that 
\begin{equation}
  \frac{\PartD w}{\PartD b} = - \frac{\PartD w}{\PartD a},
\end{equation}
with
\begin{equation}
  \frac{\PartD w}{\PartD a} = 2b \cdot F_B-2a \cdot F_A - i \left[ 2b \cdot F_A+2a \cdot F_B-\frac{2}{\Delta\nu_D\cdot\pi} \right],
\end{equation}
where the last term is from $F_u$ and $\Delta\nu_D$ is the Doppler half width in frequency units.
Since the partial derivation within the line shape is often repeated, using just one calculation of $\PartD w/\PartD a$ 
and then recalculating ${\PartD a}/{\PartD x_i}$ and ${\PartD b}/{\PartD x_i}$ for different request Jacobian quantities
is efficient.

\section{Atmospheric temperatures}
%==============================================================================
\label{sec:wfuns:atmtemp}

\subsection{Common practicalities}
%
To obtain WFs for temperatures, use \wsmindex{jacobianAddTemperature}.
The calculations can either be done in ``analytical'' or ``perturbation''
manner. Retrieval grids must be specified.

A special consideration for temperature is hydrostatic equilibrium. If effects
originating in hydrostatic equilibrium shall be included in the WFs, or
not, is selected by an argument denoted as \verb|hse|. A full account of
hydrostatic equilibrium is possible for perturbation calculations, while the
analytical approach only treats the local effect (see further below).


\subsection{Perturbation calculations}
%
The size of the perturbation must be selected (in K).
Complete radiative transfer calculations are done after perturbing the
temperature field. Hence, all possible effects are included, such as changed
propagation paths through the impact of temperature on the refractive index.
Please, note that hydrostatic equilibrium comes in during the perturbation. If
\verb|hse| is set to ``on'', also \builtindoc{z\_field} is recalculated as part
of the temperature perturbation (Section~\ref{sec:fm_defs:hse}). If set to
``off'', there is no change of \builtindoc{z\_field}. That is, you must make an
active choice regarding hydrostatic equilibrium, while others effects are
included automatically.


\subsection{Analytical expressions}
%
\subsubsection{Unpolarized absorption}
%
Compared to atmospheric species, the expressions become here more complex as
temperature also affects the propagation path length ($\Delta \aPpathLng{i}$)
and the emission source term ($j_i$). Accordingly, all terms of
Eq.~\ref{eq:wfuns:dsdxi2} are relevant, and the expansion of
$\PartD\aOth{i}/\PartD\aSttElm{i}$ generates additional terms
\begin{equation}
  \frac{\PartD\StoI}{\PartD\aSttElm{i}} =
   e^{-\aOth{i}}\left[
      \frac{\PartD\aOth{i}}{\PartD\aAbsCoef{i}}
      \frac{\PartD\aAbsCoef{i}}{\PartD\aSttElm{i}} + 
      \frac{\PartD\aOth{i}}{\PartD\Delta\aPpathLng{i}}
     \frac{\PartD\Delta\aPpathLng{i}}{\PartD\aSttElm{i}}
  \right](j_i-\aStoI{i})+
  (1-e^{-\aOth{i}})\frac{\PartD j_i}{\PartD\aSttElm{i}}.  
\end{equation}
Terms part of expressions found above are not discussed separately here. The
term, $\PartD\Delta\aPpathLng{i}/\PartD\aSttElm{i}$, originates in the
constrain of hydrostatic equilibrium, and is set to zero when \verb|hse| 
is set to ``off''. Otherwise it is set as derived below.
The term $\PartD\aAbsCoef{i}/\PartD\aSttElm{i}$ is calculated in a pure
numerical manner, by perturbing the temperature. Eq.~\ref{eq:taustep} gives
\begin{equation}
  \frac{\PartD \aOth{i}}{\PartD \Delta \aPpathLng{i}} = 
  \frac{\aAbsCoef{i}+\aAbsCoef{i+1}}{2}. 
\end{equation}
The path length ($\Delta \aPpathLng{i}$), for a given pressure, is linearly
proportional to the temperature, and if $\bar{\Tmp}$ is the average temperature
along the path step:
\begin{equation}
  \frac{\PartD \Delta \aPpathLng{i}}{\PartD \bar{\Tmp}} =   
                                    \frac{\Delta \aPpathLng{i}}{\bar{\Tmp}}.
  \label{eq:wfuns:ddldtbar}
\end{equation}
Following the other variables, we set $\bar{\Tmp}=(\aTmp{i}+\aTmp{i+1})/2$, and
\begin{equation}
  \frac{\PartD \Delta \aPpathLng{i}}{\PartD\aSttElm{i}} = 
                                    \frac{\Delta \aPpathLng{i}}{2\aTmp{i}}.
  \label{eq:wfuns:ddldti}
\end{equation}
In summary (assuming \verb|hse| set to ``on''):
\begin{equation}
  \frac{\PartD\StoI}{\PartD\aSttElm{i}} =
   e^{-\aOth{i}}\frac{\Delta \aPpathLng{i}}{2}\left[
      \frac{\PartD\aAbsCoef{i}}{\PartD\aSttElm{i}} + 
  \frac{\aAbsCoef{i}+\aAbsCoef{i+1}}{2\aTmp{i}}
  \right](j_i-\aStoI{i})+
  (1-e^{-\aOth{i}})\frac{\PartD j_i}{\PartD\aSttElm{i}}.  
\end{equation}
The derivative of the Planck function ($\PartD\Planck/\PartD\aSttElm{i}$) can
be expressed analytically \citep{eriksson:studi:02}. 

The expression for higher Stokes elements, or if emission is totally ignored,
is obtained by setting $j_i$ (and $\PartD j_i/\PartD\aSttElm{i}$) to zero.


\subsubsection{General case}
%
Both line strength and line shape change with temperature.
The rough expression is
\begin{equation}
  \frac{\PartD \MtrStl{K}_i}{\PartD x_i} = \MtrStl{K}_i
  \left(
  S_j^{-1}\frac{\PartD S_j}{\PartD x_i}+
  F_j^{-1}\frac{\PartD F_j}{\PartD x_i}
  \right),
\end{equation}
where individual terms of both line shape and line strength are
treated separately.

\subsubsection{Hydrostatic equilibrium and limitations}
%
A changed temperature has non-local effects, originating from refraction and
hydrostatic equilibrium. The expressions above ignore totally refraction
effects. 

As  mentioned,  if \verb|hse|  is  set  to  ``off'',  the term  $\PartD  \Delta
\aPpathLng{i}/\PartD\aSttElm{i}$  is set  to  zero. That  is,  the path  length
through the layer is not affected  by a temperature change. With \verb|hse| set
to ``on'', the complete expressions above are used.

If this treatment of hydrostatic equilibrium is sufficient or not depends on
the observation geometry. It should be insufficient for limb sounding,
where changes even at altitudes below the tangent point can have an influence
as the geometrical altitudes of all higher layers is changed through
hydrostatic equilibrium.  However, this effect vanishes for
ground-based observations at zenith and satellite measurements at nadir, giving
a full account of hydrostatic equilibrium even with the analytical expressions.
In practise it should be possible to apply the expressions outside zenith and
nadir, as long as the observations are of ``up'' or ``down-ward'' type. The
same applies to measurements from inside the atmosphere, (e.g.\ aircraft ones),
if the reference pressure for hydrostatic equilibrium (\builtindoc{p\_hse}) is
matched to the pressure of the observation point.

Non-LTE is presently not supported for HSE.

\section{Magnetic field}
%==============================================================================
\label{sec:wfuns:magneticfield}

\subsection{Common practicalities}
%
To obtain WFs for magnetic field components, use \wsmindex{jacobianAddMagField}.
 Retrieval grids must be specified.  There are six input options for component.
 These are ``u'', ``v'', ``w'', ``theta'', ``eta'', and ``strength''.  The 
 last of these three are derived analytically, the other three small
 perturbations.  Note that the Faraday effect is not supporting magnetic field
 weighting function calculations at the time of writing this.

\subsection{The analytical solutions}

\subsubsection{Strength component}
Magnetic field strength affects the center of the line.  This means it affects the line 
strength through changing the state distribution by changing the energy differences 
between states.  It also means it affects the line shape by moving the line away from its 
original line center.  The former effect is ignored because it is very small 
compared to the latter effect for reasonable atmospheric magnetic fields.

In rough estimation, the magnetic field change the energy state of the molecule by
\begin{equation}
  \Delta E = -g H,
\end{equation}
where $g$ is a constant (see the theory guide) and $H$ is the field strength.  The frequency change
is then
\begin{equation}
  \Delta \nu_0 = \frac{1}{h}\left(g''-g'\right)H,
\end{equation}
where $h$ is Planck's constant and the primes stand for upper (single prime)
and lower (double prime) of the state.

Since we ignore line strength partial derivations,
both Equations~\ref{eq: wind propagation derivative}~and~\ref{eq: line shape derivative}
applies for the magnetic field strength. The remaining partial derivative is then just
\begin{equation}
  \frac{\PartD a}{\PartD x_i} = \frac{\Delta \nu_0}{\Delta\nu_D};\;\;\;\frac{\PartD b}{\PartD x_i}=0.
\end{equation}

\subsubsection{Angular components}
The angle along the path of observation ($\theta$) only changes the polarization state. 
The angle of linear polarization rotation ($\eta$) also only changes the polarization state.
So
\begin{equation}
   \frac{\PartD \MtrStl{K}_i}{\PartD x_i} = \MtrStl{K}_i \MtrStl{\Phi}^{-1}\frac{\PartD \MtrStl{\Phi}}{\PartD x_i}.
\end{equation}
The derivation of a matrix with respect to a variable is just the derivation of individual parameters, so see
the theory guide for the original shape and content of $\MtrStl{\Phi}$.

\section{Non-LTE effects}
%==============================================================================
\label{sec:wfuns:nlte}

This is handled passively whenever it is encountered.  If non-LTE is off, then the source function is
the Planck function, so $\aSrcVec{n,i}=0$ and $\PartD \aSrcVec{n,i} /\PartD x_i=0$ in Equation~\ref{eq:wfuns:dsdxi1}.
Otherwise, we must note two things.  First is that the non-LTE effect is included in $f(T,\cdots)$ as described above,
so we know $\MtrStl{K}_i$ and ${\PartD \MtrStl{K}_i}/{\PartD x_i}$ equally well as before.  The second part is that 
we have defined $\aSrcVec{n,i}$ from Equation~\ref{eq:nlte_source}.  By this definition, for all effects the NLTE
contribution is 
\begin{equation}
  \frac{\PartD \aSrcVec{n,i}}{\PartD x_i} = \Planck
   \left[\frac{\PartD \aAbsVec{s}}{\PartD x_i}\oslash\AbsVec - \frac{\PartD \AbsVec}{\PartD x_i}\odot\aAbsVec{s}\oslash\left(\AbsVec\odot\AbsVec\right)\right]
  +\frac{\PartD\Planck}{\PartD x_i}\left(\aAbsVec{s}\oslash\AbsVec-1\right).
\end{equation}

In practice, the above is only an overview and the calculations are done at a much lower level.  In fact
\begin{equation}
  \aSrcVec{n,i} = \left(\frac{f_s(T,\cdots)}{f(T,\cdots)}-1\right)\Planck\MtrStl{\Phi}\cdot[1,\,0,\,0,\,0],
\end{equation}
where the ratio of $f_s/f$ gives the ratio of the source emission to the absorption due to non-LTE effects.
These $f$s are otherwise as in Equation~\ref{eq: rough line strength} for calculations of $\MtrStl{K}_j$.

\section{Sensor pointing}
%==============================================================================
\label{sec:wfuns:sensorpointing}

The term ``sensor pointing'' refers to deviations between nominal and
actual viewing direction of the sensor. So far, only deviations in zenith angle
can be considered. The workspace method to initiate such Jacobians is
\wsmindex{jacobianAddPointingZa}.

The pointing deviation is treated as a time varying variable, then having a
polynomial variation. hence, the basis functions described in
Section~\ref{sec:wfuns:basis2} are applied. The time is taken from
\wsvindex{sensor\_time}. If the pointing error is assumed to be constant with
time, the polynomial order to select is 0, and so on. As a special case, the
polynomial order -1 signifies here that the pointing off-set is so highly
varying that an off-set must be assigned to each spectrum
(sometime called ``pointing jitter'').

The Jacobian can be calculated in two manners:
\begin{description}
\item[recalc] If this option is selected, radiative transfer calculations are
  performed for a shift of \builtindoc{sensor\_los} (perturbation size selected
  by \verb|dza|). Only a ``one-sided'' perturbation is applied.
\item[interp] The Jacobian is derived from existing data, by an interpolation
  of existing data. This achieved by interpolating pencil beam data to a
  shifted zenith angle grid. This will involve some extrapolation of the data,
  and this aspect should be considered when selecting the zenith angles in
  \builtindoc{mblock\_dlos\_grid}. The average of a positive and negative shift
  is determined. The shift to apply (\verb|dza|) should be smaller than the
  minimum spacing of the zenith angles in\builtindoc{mblock\_dlos\_grid} for
  accurate results. As interpolation is a relative fast operation and a
  ``double-sided'' disturbance is used, this option should in general be
  preferred.
\end{description}




\section{Sensor frequencies}
%==============================================================================
\label{sec:wfuns:sensorfreq}

This class of Jacobians treats deviations between nominal and actually recorded
frequencies. Such differences can originate in several ways, but the exact
origin can normally be ignored and the effect can be modelled as the backend
(spectrometer, filter bank \dots) channels are shifted from their nominal
position. The workspace methods to define such Jacobians are
\wsmindex{jacobianAddFreqShift} and \wsmindex{jacobianAddFreqStretch}. These
Jacobians can so far only be determined by applying an interpolation of existing
monochromatic data, then shifted \verb|df| from the nominal
values.

The methods treat either the ``shift'' or ``stretch'' effects. This follows
standard nomenclature. A ``shift'' is an off-set that is of the same size for
all backend channels. That is, if only a shift is assumed, the nominal
distances between the channels are assumed to be valid. The ``stretch'' term
considers the distance between the channels. For a backend with all channels
equally spaced, a stretch signifies that the spacing deviates from the nominal
value (but all channels still equally spaced). More generally, a stretch means
that the deviation from the nominal channel position increases linearly from
the middle point of the backend. In terms of the basis functions
Section~\ref{sec:wfuns:basis2}, shift and stretch correspond to polynomial
order 0 and 1.

Both frequency shift and stretch can be assumed to be time varying, where
exactly the same polynomial approach as for pointing is applied. This
including the case of setting the order to -1.




\section{Polynomial baseline fit}
%==============================================================================
\label{sec:wfuns:polyfit}

A ``baseline'' is microwave jargon for a disturbance of the spectrum that is
not covered by the common sensor characteristics. The most common case is that
the local oscillator signal leaks into the measurement by reflections occurring
inside the sensor, causing a pattern in the spectrum of standing-wave type.
Such effects are difficult to model in a physical manner, and a more general
fitting procedure must be applied. A common option is then to model the
baseline as a polynomial, of a specified order. That is, assuming a measurement
giving a single spectrum, the measured spectrum \MsrVct\ is modelled as
\begin{equation}
  \label{eq:polyfit}
  \MsrVct = \MsrVct' + \sum_{i=0}^n \aSttElm{i}\VctStl{z}_i,
\end{equation}
where $\MsrVct'$ is the ``baseline free'' spectrum, and \aSttElm{i}\ and
$\VctStl{z}_i$ are introduced in Section~\ref{sec:wfuns:basis2}.

The Jacobian for such baseline models are obtained by
\wsmindex{jacobianAddPolyfit}. For single spectra measurements, the only
consideration is the polynomial order to use. For measurements where several
spectra are appended to form the measurement vector, the default option is that
baseline can vary between all spectra. In some cases, it could be assumed that
the baseline is common between data for different polarizations, viewing
directions or measurement blocks, and flags can be set to mimic such
assumptions.

For a given set of retrieved \aSttElm{i}, the simplest way to determine the
estimated baseline is to perform a multiplication between the relevant parts of
the Jacobian and the state vector:
\begin{equation}
  \aWfnMtr{p}\aSttVct{p},
  \label{eq:bline:recreate}
\end{equation}
where $p$ indicates the $n+1$ columns/elements corresponding to \aSttElm{i}.


\section{Sinusoidal baseline fit}
%==============================================================================
\label{sec:wfuns:sinefit}

If the baseline has components of sinusoidal character there is also a second
option, provided the method \wsmindex{jacobianAddSinefit}. The baseline is
in this method modelled as \citep{kuntz:97}
\begin{equation}
  \label{eq:sinefit}
  \MsrVct = \MsrVct' + \sum_{i=1}^n 
       a_i\sin\left(\frac{2\pi(\Frq-\aFrq{1})}{l_i}\right)+
       b_i\cos\left(\frac{2\pi(\Frq-\aFrq{1})}{l_i}\right)
\end{equation}
where $a_i$ and $b_i$ are the coefficients to be retrieved, \Frq\ is frequency,
\aFrq{i} is a reference frequency and $l_i$ is period length. The reference
frequency (\aFrq{i}) is in practice taken as the first frequency of the
spectrometer. 

The period lengths ($l_i$) are user input. For each given period length, the
corresponding sine and cosine terms are included in the Jacobian. As for the
polynomial fit, there exist options to set the baseline to be common between
different polarizations, viewing directions or measurement blocks.
Equation~\ref{eq:bline:recreate} is also applicable for sinusoidal baseline
fits.

An alternative way to express the expression above is
\begin{equation}
  \label{eq:sinefit2}
  \MsrVct = \MsrVct' + \sum_{i=1}^n A_i
             \sin\left(\frac{2\pi(\Frq-\aFrq{1})}{l_i} + \phi_i\right)
\end{equation}
where
\begin{equation}
  A_i = \sqrt{a_i^2+b_i^2}
\end{equation}
and
\begin{equation}
  \tan{\phi_i} = \frac{a_i}{b_i}.
\end{equation}
Equation~\ref{eq:sinefit} is used to define the weighting functions as this
gives a linear retrieval problem, in contrast to if Equation~\ref{eq:sinefit2}
would be used, that would require an iterative process to determine $A_i$
and $\phi_i$.



