m_wfs.cc

Go to the documentation of this file.

/* Copyright (C) 2000, 2001 Patrick Eriksson <patrick@rss.chalmers.se>

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */


//   File description


//   External declarations

#ifdef SUNOS
# include <ieeefp.h>
#endif

#include <math.h>
#include "arts.h"
#include "auto_md.h"
#include "matpackI.h"
#include "messages.h"          
#include "auto_wsv.h"          
#include "math_funcs.h"          
#include "atm_funcs.h"          
#include "los.h"
extern const Numeric PLANCK_CONST;
extern const Numeric BOLTZMAN_CONST;



//   Function(s) to join two WF matrices and related variables


void k_append (
                    Matrix&          kx,
                    ArrayOfString&   kx_names,
                    ArrayOfIndex&    kx_lengths,
                    Matrix&          kx_aux,
              const Matrix&          k,
              const ArrayOfString&   k_names,
              const Matrix&          k_aux )
{
  // Size of Kx and K
  const Index  ny1  = kx.nrows();         // length of measurement vector (y)
  const Index  nx1  = kx.ncols();         // length of state vector (x)
  const Index  nri1 = kx_names.nelem();    // number of retrieval identities
  const Index  ny2  = k.nrows();  
  const Index  nx2  = k.ncols();  
  const Index  nri2 = k_names.nelem();
  Index        iri;

  // Make copy of Kx data.
  Matrix          ktemp(kx);
  ArrayOfString   ktemp_names(kx_names);
  ArrayOfIndex    ktemp_lengths(kx_lengths);
  Matrix          ktemp_aux(kx_aux);
  
//   cout << "ktemp =\n" << ktemp << "\n\n\n";
//   cout << "ktemp_names =\n" << ktemp_names << "\n\n\n";
//   cout << "ktemp_lengths =\n" << ktemp_lengths << "\n\n\n";
//   cout << "ktemp_aux =\n" << ktemp_aux << "\n\n\n";

  if ( nx1 > 0 )
  {
    if ( ny1 != ny2 )
      throw runtime_error(
            "The two WF matrices have different number of rows." ); 
  }

  // Reallocate the Kx data
  kx.resize(         ny2,       nx1+nx2 );
  kx_names.resize(   nri1+nri2          );
  kx_lengths.resize( nri1+nri2          );
  kx_aux.resize(     nx1+nx2,   2       );

  // Move Ktemp to Kx
  if ( nx1 > 0 )
  {
    kx( Range(joker), Range(0,nx1) )     = ktemp;
    kx_aux( Range(0,nx1), Range(joker) ) = ktemp_aux;
    // For the Array types kx_names and kx_lenghts we cannot use Range.  
    for ( iri=0; iri<nri1; iri++ )
    {
      kx_lengths[iri]  = ktemp_lengths[iri];
      kx_names[iri]    = ktemp_names[iri];
    }
  }

  // Calculate the vector length for each identity in K
  Index l = nx2/nri2;

  // Move K to Kx
  kx( Range(joker), Range(nx1,nx2) )     = k;
  kx_aux( Range(nx1,nx2), Range(joker) ) = k_aux;
  // For the Array types kx_names and kx_lenghts we cannot use Range.  
  for ( iri=0; iri<nri2; iri++ )
  {
    kx_names[nri1+iri]   = k_names[iri];
    kx_lengths[nri1+iri] = l;
  } 
}



//   Help functions for grid conversions


void p2grid(
              Vector&   grid,
        const Vector&   pgrid )
{
  grid.resize( pgrid.nelem() );
  grid = pgrid;                 // Matpack can copy the contents of
                                // vectors like this. The dimensions
                                // must be the same! 
  transform( grid, log, grid ); // Matpack has the same transform
                                // functionality like the old
                                // vecmat. THE ORDER OF THE ARGUMENTS
                                // IS DIFFERENT, THOUGH (output first).
  grid *= -1;                   // Mulitply all elements with -1.
}




void grid2grid_index (
                    Matrix&   is,
              const Vector&   x0,
              const Vector&   xp )
{
  const Index n0 = x0.nelem();        // length of original grid
  const Index np = xp.nelem();        // length if retrieval grid
  Index       i0, ip;                // counter for each grid

  // Resize is and set all values to -1
  is.resize( np, 2 );
  is = -1.0;                    // Matpack can set all elements like this.
 
  // Some checks
  if ( np < 2 )
    throw logic_error("The retrieval grid must have a length > 1.");
  for ( i0=1; i0<n0; i0++ )
  {
    if ( x0[i0-1] >= x0[i0] )
      throw logic_error("The grids must be sorted with increasing values.");
  }
  for ( ip=1; ip<np; ip++ )
  {
    if ( xp[ip-1] >= xp[ip] )
      throw logic_error("The grids must be sorted with increasing values.");
  }

  // Do something only if at least one point of X0 is inside XP
  if ( !( (x0[0]>xp[np-1]) || (x0[n0-1]<xp[0]) ) )
  {

    // First point of XP
    i0 = 0;
    for ( ; x0[i0] < xp[0]; i0++ ) {}
    if ( x0[i0] < xp[1] )
    {
      is(0,0) = (Numeric) i0;
      for ( ; (i0<n0) && (x0[i0]<xp[1]); i0++ ) {}
      is(0,1) = (Numeric) (i0 - 1);
    }

    // Points inside XP
    for ( ip=1; ip<(np-1); ip++ )
    {
      i0 = 0;
      for ( ; (i0<n0) && (x0[i0]<=xp[ip-1]); i0++ ) {}
      if ( (i0<n0) && (x0[i0]<xp[ip+1]) )
      {
        is(ip,0) = (Numeric) i0;
        for ( ; (i0<n0) && (x0[i0]<xp[ip+1]); i0++ ) {}
        is(ip,1) = (Numeric) (i0 - 1);
      }
    }

    // Last point of XP
    i0 = 0;
    for ( ; (i0<n0) && (x0[i0]<=xp[np-2]); i0++ ) {}
    if ( (i0<n0) && (x0[i0]<=xp[np-1]) )
    {
      is(np-1,0) = (Numeric) i0;
      for ( ; (i0<n0) && (x0[i0]<=xp[np-1]); i0++ ) {}
      is(np-1,1) = (Numeric) (i0 - 1);
    }
  }
}




void grid2grid_weights (
                    Vector&   w,
              const Vector&   x0,
              const Index&   i1,
              const Index&   i2,
              const Vector&   xp,
              const Index&   ip )
{
  const Index   np = xp.nelem();        // number of retrieval points
  const Index   nw = i2-i1+1;          // number of LOS points affected
  Index         i;

  // Reallocate w
  w.resize(nw);

  // First point of the retrieval grid
  if ( ip == 0 )
  {
    for ( i=0; i<nw; i++ )
      w[i] = ( xp[1] - x0[i1+i] ) / ( xp[1] - xp[0] );
  }
  // Points inside the retrieval grid
  else if ( ip < (np-1) )  
  {
    for ( i=0; i<nw; i++ )
    {
      if ( x0[i1+i] <= xp[ip] )
        w[i] = ( x0[i1+i] - xp[ip-1] ) / ( xp[ip] - xp[ip-1] );
      else
        w[i] = ( xp[ip+1] - x0[i1+i] ) / ( xp[ip+1] - xp[ip] );
    }
  }
  // Last point of the retrieval grid
  else
  {
    for ( i=0; i<nw; i++ )
      w[i] = ( x0[i1+i] - xp[np-2] ) / ( xp [np-1] - xp[np-2] );
  }
}


void grid2grid_weights_total (
                          ArrayOfMatrix&     W,
                    const  Los&               los, 
                    const  Vector&            k_grid) 
         
{
  // Main sizes
  const Index  nza = los.start.nelem();     // number of zenith angles  
  const Index  np  = k_grid.nelem();        // number of pressure grid
  // -log(p) is used as altitude variable. The minus is included to get
  // increasing values, a demand for the grid functions. 
  //  const Vector  lgrid=-1.0*log(k_grid);
  Vector  lgrid;
  p2grid( lgrid, k_grid );
  Vector  lplos;    
  Vector w1;
  W.resize(nza);
  // Indices
  // IP0 and IF0 are the index off-sets for the total K matrix
  Index  iza;                       // zenith angle index
  Index  ip;                        // pressure index
  Matrix is;                       // matrix for storing LOS index

  for ( iza=0; iza<nza; iza++ ) 
    {
      W[iza].resize(los.p[iza].nelem(), np);
      W[iza]=0.0;
      if ( (iza%20)==0 )
        {
          out3 << "\n      ";
          out3 << " " << iza; cout.flush();
        }
      w1=0;
      // Do something only if LOS has any points
      if ( los.p[iza].nelem() > 0 )
        {
          // Get the LOS points affected by each retrieval point
          //      lplos = -1.0 * log(los.p[iza]);
          p2grid( lplos, los.p[iza] );
          grid2grid_index (is, lplos, lgrid);
          
          for ( ip=0; ip<np; ip++ ) 
            {
              // Check if there is anything to do
              
              if ( is(ip,0) >= 0 )
                {
                  // Get the weights for the LOS points
                  grid2grid_weights( w1, lplos, (Index)is(ip,0),  (Index)is(ip,1),
                                     lgrid, ip );
                  // fill the corresponding part of the Matrix W with the weights w1

                  W[iza](Range((Index)is(ip,0),(Index)(is(ip,1)-is(ip,0)+1)), ip)=w1;
                }
            }
        }
    }
}

//   Help functions for absloswfsCalc to calculate absorption LOS WFs


void absloswfs_rte_1pass (
                    Matrix&   k,
                    Vector    y,
              const Index&   start_index,
              const Index&   stop_index,   // this variable is used by
              const Numeric&  lstep,        // absloswfs_down
              const Matrix&   tr,
              const Matrix&   s,
              const Index&    ground,
              const Vector&   e_ground,
              const Vector&   y_ground )
{
  const Index   nf = tr.nrows();     // number of frequencies
  Index         iv;                  // frequency index
  Index         q;                   // LOS index (same notation as in AUG)
  Vector        t(nf);               // transmission to the sensor
  
  t = 1.0;              // Set all elements of t to 1.

  if ( ground && ((ground-1==stop_index) || (ground-1==start_index)) )
    throw logic_error("The ground cannot be at one of the end points."); 

  // Resize K
  k.resize( nf, start_index+1 );

  // We start here at the LOS point closest to the sensor, that is,
  // reversed order compared to RTE_ITERATE  

  // The LOS point closest to the sensor 
  q  = stop_index;
  for ( iv=0; iv<nf; iv++ )    
  {
    t[iv]   *= tr(iv,q);
    y[iv]    = (y[iv]-s(iv,q)*(1.0-tr(iv,q)))/tr(iv,q);
    k(iv,q) = (-lstep/2)*(y[iv]-s(iv,q))*t[iv];
  }

  // Points inside the LOS
  for ( q=stop_index+1; q<start_index; q++ )
  {
    // Not a ground point
    if ( q != (ground-1) )
    {
      for ( iv=0; iv<nf; iv++ )    
      {
        y[iv]    = (y[iv]-s(iv,q)*(1.0-tr(iv,q)))/tr(iv,q);
        k(iv,q) = (-lstep/2)*( 2*(y[iv]-s(iv,q))*tr(iv,q) + 
                                             s(iv,q) - s(iv,q-1) ) *t[iv];
        t[iv]   *= tr(iv,q); 
      }
    }
    // A ground point
    else
    {
      out1 << "WARNING: The function absloswfs_1pass not tested for "
              "ground reflections\n";
      for ( iv=0; iv<nf; iv++ )    
      {
        y[iv]    = ( y[iv] - e_ground[iv]*y_ground[iv] - 
                            s(iv,q)*(1.0-tr(iv,q))*(1-e_ground[iv]) ) / 
                                             ( tr(iv,q) * (1-e_ground[iv]) );
        k(iv,q) = (-lstep/2)*( 2*(y[iv]-s(iv,q))*tr(iv,q)*(1-e_ground[iv])+ 
                                s(iv,q)*(1-e_ground[iv]) + 
                             e_ground[iv]*y_ground[iv] - s(iv,q-1) ) * t[iv];
        t[iv]   *= tr(iv,q) * (1-e_ground[iv]); 
      }
    }
  }

  // The LOS point furthest away from the sensor
  q = start_index;
  for ( iv=0; iv<nf; iv++ )    
    k(iv,q)  = (-lstep/2)*(y[iv]-s(iv,q-1))*t[iv];

  // To check the function: Y shall now be equal to the radiation entering 
  // the atmosphere and T the total transmission
}




void absloswfs_rte_limb (
                    Matrix&   k,
                    Vector    y,              // = y_q^q
              const Vector&   y_space,             
              const Index&   start_index,
              const Numeric&  lstep,
              const Matrix&   tr,
              const Matrix&   s,
              const Index&    ground,
              const Vector&   e_ground )
{
  const Index nf = tr.nrows();     // number of frequencies
  Index       iv;                  // frequency index
  Vector      t1q;                 // transmission tangent point - q squared
  Vector      tqn(nf,1.0);         // transmission q - sensor
  Index       q;                   // LOS index (same notation as in AUG)
  Numeric     tv, tv1 = 0.;             // transmission value for q and q-1
  Vector      yn(y_space);         // = y_space

  // Matpack can initialize a new Vector from another Vector. See how
  // yn is initialized from y_space.

  // Resize K
  k.resize( nf, start_index+1 );

  // Calculate the total transmission
  bl( t1q, start_index, start_index, tr, ground, e_ground );

  // We start at the outermost point
  q  = start_index;       
  for ( iv=0; iv<nf; iv++ )    
  {
    tv1      = tr(iv,q-1);
    t1q[iv] /= tv1*tv1;
    tqn[iv] *= tv1;
    y[iv]    = ( y[iv] - s(iv,q-1)*(1-tv1)*(1+t1q[iv]*tv1) ) / tv1;
    k(iv,q)  = (-lstep/2)*( ( 2*yn[iv] + s(iv,q-1)*(1-2*tv1) ) * 
                              t1q[iv]*tv1 + y[iv] - s(iv,q-1) )*tv1;
  }

  // Points inside the LOS
  for ( q=start_index-1; q>0; q-- )
  {
    for ( iv=0; iv<nf; iv++ )    
    {
      tv1      = tr(iv,q-1);    
      tv       = tr(iv,q);
      t1q[iv] /= tv1*tv1;
      y[iv]    = ( y[iv] - s(iv,q-1)*(1-tv1)*(1+t1q[iv]*tv1) ) / tv1;
      k(iv,q) = (-lstep/2) * ( 
           ( 4*(yn[iv]-s(iv,q))*tv1*tv + 3*(s(iv,q)-s(iv,q-1))*tv1 + 
                                        2*s(iv,q-1) ) * t1q[iv]*tv1 + 
             2*(y[iv]-s(iv,q-1))*tv1 + s(iv,q-1) - s(iv,q) ) * tqn[iv];
      tqn[iv] *= tv1;
      yn[iv]   = yn[iv]*tv + s(iv,q)*(1-tv);
    } 
  }

  // The tangent or ground point
  for ( iv=0; iv<nf; iv++ )    
    k(iv,0)  = (-lstep/2)*( (2*yn[iv]*tv1+s(iv,0)*(1-2*tv1))*t1q[iv] +
                       y[iv] - s(iv,0) ) * tqn[iv];

  // To check the function
  // Without ground reflection: T1Q=1 and Y=0
  // With ground reflection: T1Q=(1-e) and Y=eB
  // print_vector( t1q );
  // print_vector( y );
}




void absloswfs_rte_down (
                    Matrix&   k,
              const Vector&   y,
              const Vector&   y_space,
              const Index&   start_index,
              const Index&   stop_index,
              const Numeric&  lstep,
              const Matrix&   tr,
              const Matrix&   s,
              const Index&    ground,
              const Vector&   e_ground,
              const Vector&   y_ground )
{
  const Index   nf = tr.nrows(); // number of frequencies
        Index   iv;              // frequency index
        Index   q;               // LOS index (same notation as in AUG)
        Vector   y0(nf);          // see below
        Matrix   k2;              // matrix for calling other LOS WFs functions
        Vector   tr0;             // see below

  // Resize K
  k.resize( nf, start_index+1 );

  // Calculate Y0, the intensity reaching the platform altitude at the far
  // end of LOS, that is, from above.
  y0 = y_space;                 // Matpack can copy the contents of
                                // vectors like this. The dimensions
                                // must be the same! 
  rte_iterate( y0, start_index-1, stop_index, tr, s, nf );

  // Calculate TR0, the transmission from the platform altitude down to the
  // tangent point or the ground, and up to the platform again.
  bl( tr0, stop_index, stop_index, tr, ground, e_ground );

  // The indeces below STOP_INDEX are handled by the limb sounding function.
  // The limb function is given Y0 instead of cosmic radiation 
  absloswfs_rte_limb( k2, y, y0, stop_index, lstep, tr, s, ground, e_ground );
  for ( iv=0; iv<nf; iv++ )
  {
    for ( q=0; q<stop_index; q++ )
      k(iv,q) = k2(iv,q);
  }

  // The indeces above STOP_INDEX are handled by the 1pass function.
  // The transmission below STOP_INDEX must here be considered.
  absloswfs_rte_1pass( k2, y0, start_index, stop_index, lstep, tr, s, 
                     ground, e_ground, y_ground );
  for ( iv=0; iv<nf; iv++ )
  {
    for ( q=stop_index+1; q<=start_index; q++ )
      k(iv,q) = k2(iv,q)*tr0[iv];
  }

  // The platform altitude must be treated seperately
  //
  // Calculate the intensity generated below point q-1, YQQ
  Vector yqq(nf);
  rte( yqq, stop_index-1, stop_index-1, tr, s, Vector(nf,0.0), 
                                            ground, e_ground, y_ground );
  //
  // Y0 is moved one step upwards and TR0 one step downwards
  q = stop_index; 
  for ( iv=0; iv<nf; iv++ )
  {
    tr0[iv] /= tr(iv,q-1)*tr(iv,q-1);    
    y0[iv]   = ( y0[iv] - s(iv,q)*(1-tr(iv,q)) ) / tr(iv,q);
    k(iv,q) = (-lstep/2)*( (3*(y0[iv]-s(iv,q))*tr(iv,q-1)*tr(iv,q) + 
                  2*(s(iv,q)-s(iv,q-1))*tr(iv,q-1) + s(iv,q-1) )*tr0[iv] + 
                                         yqq[iv] - s(iv,q-1) )*tr(iv,q-1);
  }
}




void absloswfs_rte (
                    ArrayOfMatrix&   absloswfs,
              const Los&             los,   
              const ArrayOfMatrix&   source,
              const ArrayOfMatrix&   trans,
              const Vector&          y,
              const Vector&          y_space,
              const Vector&          f_mono,
              const Vector&          e_ground,
              const Numeric&         t_ground )
{
  const Index  nza = los.start.nelem();   // number of zenith angles  
  const Index  nf  = f_mono.nelem();      // number of frequencies  
        Vector  yp(nf);                   // part of Y
        Index  iy0=0;                    // y index

  // Set up vector for ground blackbody radiation
  Vector   y_ground(f_mono.nelem()); 
  if ( any_ground(los.ground) )
    planck( y_ground, f_mono, t_ground );

  // Resize the LOS WFs array
  absloswfs.resize( nza );

  // Loop zenith angles
  out3 << "    Zenith angle nr:\n      ";
  for (Index i=0; i<nza; i++ ) 
  {
    if ( ((i+1)%20)==0 )
      out3 << "\n      ";
    out3 << " " << i; cout.flush();
    
    // Do something only if LOS has any points
    if ( los.p[i].nelem() > 0 )
    {

      // Pick out the part of Y corresponding to the present zenith angle
      yp = y[Range(iy0,nf)];

      // The calculations are performed in 3 sub-functions
      //
      // Upward looking (=single pass)
      if ( los.stop[i]==0 )
        absloswfs_rte_1pass( absloswfs[i], yp, los.start[i], 0, los.l_step[i], 
                     trans[i], source[i], los.ground[i], e_ground, y_ground );

      //
      // 1D limb sounding
      else if ( los.start[i] == los.stop[i] )
        absloswfs_rte_limb( absloswfs[i], yp, y_space, los.start[i], 
                los.l_step[i], trans[i], source[i], los.ground[i], e_ground );

      //
      // 1D downward looking
      else 
        absloswfs_rte_down( absloswfs[i], yp, y_space, los.start[i], 
                            los.stop[i], los.l_step[i], trans[i], source[i], 
                                           los.ground[i], e_ground, y_ground );
    }

    iy0 += nf;        
  }

  out3 << "\n";
}




void absloswfs_tau (
                    ArrayOfMatrix&   absloswfs,
              const Los&             los,
              const Vector&          f_mono )
{
  const Index  nza = los.start.nelem();   // number of zenith angles  
  const Index  nf  = f_mono.nelem();      // number of frequencies  
  Numeric      kw, kw2;
  Index        row, col, np;

  // Resize the LOS WFs array
  absloswfs.resize( nza );

  // Loop zenith angles
  out3 << "    Zenith angle nr:\n      ";
  for (Index i=0; i<nza; i++ ) 
  {
    if ( ((i+1)%20)==0 )
      out3 << "\n      ";
    out3 << " " << i; cout.flush();

    np = los.start[i] + 1;

    absloswfs[i].resize( nf, np );
    
    // Do something only if LOS has any points
    if ( los.p[i].nelem() > 0 )
    {

      // Upward looking (=single pass)
      if ( los.stop[i]==0 )
      {
        kw  = los.l_step[i] / 2.0;
        kw2 = los.l_step[i];
        for( row=0; row<nf; row++ )
        {
          absloswfs[i](row,0) = kw;
          absloswfs[i](row,np-1) = kw;
          for( col=1; col<np-1; col++ )
            absloswfs[i](row,col) = kw2;
        }
      }

      // 1D limb sounding
      else if ( los.start[i] == los.stop[i] )
      {
        kw  = los.l_step[i];
        kw2 = los.l_step[i] * 2.0;
        for( row=0; row<nf; row++ )
        {
          absloswfs[i](row,0) = kw;
          absloswfs[i](row,np-1) = kw;
          for( col=1; col<np-1; col++ )
            absloswfs[i](row,col) = kw2;
        }
      }

      // 1D downward looking
      else 
      {
        kw  = los.l_step[i];
        kw2 = los.l_step[i] * 2.0;
        for( row=0; row<nf; row++ )
        {
          absloswfs[i](row,0) = kw;
          absloswfs[i](row,np-1) = kw / 2.0;
          absloswfs[i](row,los.stop[i]) = kw * 1.5;
          for( col=1; col<los.stop[i]; col++ )
            absloswfs[i](row,col) = kw2;
          for( col=los.stop[i]+1; col<np-1; col++ )
            absloswfs[i](row,col) = kw;
        }
      }
    }
  }
  out3 << "\n";
}



//   Functions to calculate source function LOS WFs


void sourceloswfs_1pass (
                    Matrix&   k,
              const Index&   start_index,
              const Index&   stop_index,   // this variable is used by 1D down
              const Matrix&   tr,
              const Index&    ground,
              const Vector&   e_ground )
{
  const Index   nf = tr.nrows();     // number of frequencies
        Index   iv;                  // frequency index
        Vector   t(nf,1.0);           // transmission to the sensor
        Index   q;                   // LOS index (same notation as in AUG) 

  if ( ground && ((ground-1==stop_index) || (ground-1==start_index)) )
    throw logic_error("The ground cannot be at one of the end points."); 

  // Resize K
  k.resize( nf, start_index+1 );

  // We start here at the LOS point closest to the sensor, that is,
  // reversed order compared to RTE_ITERATE  

  // The LOS point closest to the sensor 
  q  = stop_index;
  for ( iv=0; iv<nf; iv++ )    
    k(iv,q)  = ( 1 - tr(iv,q) ) / 2;

  // Points inside the LOS
  for ( q=stop_index+1; q<start_index; q++ )
  {
    // Not a ground point
    if ( q != ground )
    {
      for ( iv=0; iv<nf; iv++ )    
      {
        k(iv,q) = ( 1 - tr(iv,q-1)*tr(iv,q) ) * t[iv] / 2;
        t[iv]  *= tr(iv,q); 
      }
    }
    // A ground point
    else
    {
      out1 << "WARNING: The function sourceloswfs_1pass not tested "
              "for ground reflections\n";

      for ( iv=0; iv<nf; iv++ )    
      {
        k(iv,q) = ( (1-tr(iv,q))*(1-e_ground[iv])*tr(iv,q-1) + 1 - 
                                                     tr(iv,q-1) ) * t[iv] / 2;
        t[iv]  *= tr(iv,q)*(1-e_ground[iv]); 
      }
    }
  }

  // The LOS point furthest away from the sensor
  q = start_index;
  for ( iv=0; iv<nf; iv++ )    
    k(iv,q)  = ( 1 - tr(iv,q-1) ) * t[iv] / 2;
}




void sourceloswfs_limb (
                    Matrix&   k,
              const Index&   start_index,
              const Matrix&   tr,
              const Index&    ground,
              const Vector&   e_ground )
{
  const Index   nf = tr.nrows();      // number of frequencies
        Index   iv;                  // frequency index
        Vector   t1q;                 // transmission tangent point - q squared
        Vector   tqn(nf,1);           // transmission q - sensor
        Index   q;                   // LOS index (same notation as in AUG)

  // Calculate the total transmission
  bl( t1q, start_index, start_index, tr, ground, e_ground );

  // Resize K
  k.resize( nf, start_index+1 );

  // We start at the outermost point
  q  = start_index;       
  for ( iv=0; iv<nf; iv++ )    
  {
    t1q[iv] /= tr(iv,q-1) * tr(iv,q-1);
    k(iv,q)  = ( (1-tr(iv,q-1))*t1q[iv]*tr(iv,q-1) + 1 - tr(iv,q-1) )/2;
  }

  // Points inside the LOS
  for ( q=start_index-1; q>0; q-- )
  {
    for ( iv=0; iv<nf; iv++ )    
    {
      t1q[iv]  /= tr(iv,q-1) * tr(iv,q-1);
      k(iv,q)  = ( (1-tr(iv,q-1)*tr(iv,q))*t1q[iv]*tr(iv,q-1)*
         tr(iv,q) + (1-tr(iv,q-1))*tr(iv,q) + 1 - tr(iv,q) ) * tqn[iv] / 2;
      tqn[iv]  *= tr(iv,q-1);
    } 
  }

  // The tangent or ground point
  for ( iv=0; iv<nf; iv++ )    
    k(iv,0)  = ( (1-tr(iv,0))*(1+t1q[iv]*tr(iv,0)) ) * tqn[iv] / 2;
}




void sourceloswfs_down (
                    Matrix&   k,
              const Index&   start_index,
              const Index&   stop_index,
              const Matrix&   tr,
              const Index&    ground,
              const Vector&   e_ground )
{
  const Index   nf = tr.nrows(); // number of frequencies
        Index   iv;             // frequency index
        Index   q;              // LOS index (same notation as in AUG)
        Matrix   k2;             // matrix for calling other LOS WFs functions
        Vector   tr0;            // see below

  // Resize K
  k.resize( nf, start_index+1 );

  // Calculate TR0, the transmission from the platform altitude down to the
  // tangent point or the ground, and up to the platform again.
  bl( tr0, stop_index, stop_index, tr, ground, e_ground );

  // The indeces below STOP_INDEX are handled by the limb sounding function.
  sourceloswfs_limb( k2, stop_index, tr, ground, e_ground );  
  for ( iv=0; iv<nf; iv++ )
  {
    for ( q=0; q<stop_index; q++ )
      k(iv,q) = k2(iv,q);
  }

  // The indecies above STOP_INDEX are handled by the 1pass function.
  // The transmission below STOP_INDEX must here be considered.
  sourceloswfs_1pass( k2, start_index, stop_index, tr, ground, e_ground );
  for ( iv=0; iv<nf; iv++ )
  {
    for ( q=stop_index+1; q<=start_index; q++ )
      k(iv,q) = k2(iv,q)*tr0[iv];
  }

  // The platform altitude must be treated seperately
  //
  // TR0 is moved one step downwards
  q = stop_index; 
  for ( iv=0; iv<nf; iv++ )
  {
    tr0[iv] /= tr(iv,q-1)*tr(iv,q-1);    
    k(iv,q)  = ( (1-tr(iv,q-1)*tr(iv,q))*tr0[iv]*tr0[iv]*tr(iv,q-1) + 1 - 
                 tr(iv,q-1) ) / 2;
  }
}




void sourceloswfs (
                    ArrayOfMatrix&   sourceloswfs,
              const Los&             los,   
              const ArrayOfMatrix&   trans,
              const Vector&          /* f_mono */,
              const Vector&          e_ground )
{
  const Index  nza = los.start.nelem();   // number of zenith angles  

  // Resize the LOS WFs array
  sourceloswfs.resize(nza);

  // Loop zenith angles
  out3 << "    Zenith angle nr:      ";
  for (Index i=0; i<nza; i++ ) 
  {
    if ( (i%20)==0 )
      out3 << "\n      ";
    out3 << " " << i; cout.flush();
    
    // Do something only if LOS has any points
    if ( los.p[i].nelem() > 0 )
    {
      // The calculations are performed in 3 sub-functions
      //
      // Upward looking (=single pass)
      if ( los.stop[i]==0 )
        sourceloswfs_1pass( sourceloswfs[i], los.start[i], 0, trans[i], 
                                                     los.ground[i], e_ground );
      //
      // 1D limb sounding
      else if ( los.start[i] == los.stop[i] )
        sourceloswfs_limb( sourceloswfs[i], los.start[i], trans[i], 
                                                     los.ground[i], e_ground );
      //
      // 1D downward looking
      else 
        sourceloswfs_down( sourceloswfs[i], los.start[i], los.stop[i], 
                                           trans[i], los.ground[i], e_ground );
    }
  }
  out3 << "\n";
}



//   Core functions for analytical WFs.
//     Analytical expressions are used for species, temperature without 
//     hydrostatic eq. and continuum absorption.
//     These functions are very similar and a change in one of the functions
//     should most likely be included in the other functions.


void k_species (
                    Matrix&          k,
                    ArrayOfString&   k_names,
                    Matrix&          k_aux,
              const Los&             los,           
              const ArrayOfMatrix&   absloswfs,
              const Vector&          p_abs,
              const Vector&          t_abs,             
              const TagGroups&       tags,
              const ArrayOfMatrix&   abs_per_tg,
              const Matrix&          vmrs,
              const Vector&          k_grid,
              const ArrayOfIndex&    tg_nr,
              const String&          unit )
{
  check_lengths( p_abs, "p_abs", t_abs, "t_abs" );  
  check_length_ncol( p_abs, "p_abs", vmrs, "vmrs" );
  if ( tags.nelem() != abs_per_tg.nelem() )
    throw runtime_error(
                       "Lengths of *wfs_tgs* and *abs_per_tg* do not match." );
  if ( los.p.nelem() != absloswfs.nelem() )
    throw runtime_error(
     "The number of zenith angles is not the same in *los* and *absloswfs*." );

  // Main sizes
  const Index  nza = los.start.nelem();      // number of zenith angles  
  const Index  nv  = abs_per_tg[0].nrows();  // number of frequencies
  const Index  ntg = tg_nr.nelem();          // number of retrieval tags to do
  const Index  np  = k_grid.nelem();         // number of retrieval altitudes

  // -log(p) is used as altitude variable. The minus is included to get
  // increasing values, a demand for the grid functions. 
  Vector  lgrid;                      // -log of the retrieval pressures
  p2grid( lgrid, k_grid );
  Vector  lplos;                      // -log of the LOS pressures

  // Indices
  // IP0 and IF0 are the index off-sets for the total K matrix
        Index  itg;                       // Tag index
        Index  iza;                       // Zenith angle index
        Index  ip, ip0=0;                 // Retrieval point indices
        Index  iv, iv0;                   // Frequency indices
        Index  i1, iw;                    // weight indices

  // Other variables
        Matrix  abs;                       // absorption at k_grid
        Matrix  is;                        // matrix for storing LOS index
        Vector  w;                         // weights for LOS WFs
        Vector  a(nv);                     // temporary vector
        Index  tg;                        // present tag nr
        Vector  vmr, p, t;                 // for conversion to VMR and ND 
        Vector  nd;                        // number density

  // Set up K and additional data. Set all values of K to 0
  k.resize(nza*nv,ntg*np);
  k = 0.0;                      // Matpack can set all elements like this.
  k_names.resize(ntg);
  k_aux.resize(ntg*np,2);

  // The calculations
  // Loop order:
  //   1 tags
  //   2 zenith angle
  //   3 retrieval altitudes 
  //   4 frequencies
  for ( itg=0; itg<ntg; itg++ ) 
  {
    // Present tag nr
    tg = tg_nr[itg];

    out2 << "  Doing " << tags[tg][0].Name() << "\n";

    // Get the absorption of the species at the retrieval points
    abs.resize( nv, np );
    interpp( abs, p_abs, abs_per_tg[tg], k_grid );

    // Fill K_NAMES
    {
      ostringstream os;
      os << "Species: " << tags[tg][0].Name();
      k_names[itg] = os.str();
    }

    // Fill column 0 of K_AUX
    for ( ip=0; ip<np; ip++ )
       k_aux(ip0+ip,0) = k_grid[ip];

    // Fill column 1 of K_AUX
    //   frac : fractions of the linearisation profile
    //   vmr  : VMR
    //   nd   : number density
    //
    vmr.resize( np );
    interpp( vmr, p_abs, vmrs(tg,Range(joker)), k_grid ); 
    // The Range(joker) expression selects the entire row.

    if ( unit == "frac" )
    {
      for ( ip=0; ip<np; ip++ )
        k_aux(ip0+ip,1) = 1.0;
    }
    else if ( unit == "vmr" )
    {
      for ( ip=0; ip<np; ip++ )
      {
        for ( iv=0; iv<nv; iv++ )
          abs(iv,ip)   /= vmr[ip];
        k_aux(ip0+ip,1) = vmr[ip];
      }
    }  
    else if ( unit == "nd" )
    {
      nd.resize( np );
      interpp(  nd, p_abs, number_density(p_abs,t_abs), k_grid );
      nd *= vmr;                // Matpack can do element-vise
                                // multiplication like this.
      for ( ip=0; ip<np; ip++ )
      {
        for ( iv=0; iv<nv; iv++ )
          abs(iv,ip)   /= nd[ip];
        k_aux(ip0+ip,1) = nd[ip];
      }
    }
    else
      throw runtime_error(
        "Allowed species retrieval units are \"frac\", \"vmr\" and \"nd\"."); 

    // Set frequency zenith angle index off-set to 0
    iv0 = 0;                 

    // Loop zenith angles
    out3 << "    Zenith angle nr:\n      ";
    for ( iza=0; iza<nza; iza++ ) 
    {
      if ( ((iza+1)%20)==0 )
        out3 << "\n      ";
      out3 << " " << iza; cout.flush();
      
      // Do something only if LOS has any points
      if ( los.p[iza].nelem() > 0 )
      {
        // Get the LOS points affected by each retrieval point
        //        lplos = -1.0 * log(los.p[iza]);
        p2grid( lplos, los.p[iza] );
        grid2grid_index( is, lplos, lgrid );

        // Loop retrieval points
        for ( ip=0; ip<np; ip++ ) 
        {
          // Check if there is anything to do
          if ( is(ip,0) >= 0 )
          {
            // Get the weights for the LOS points
            grid2grid_weights( w, lplos, (Index)is(ip,0), (Index)is(ip,1), 
                                                                 lgrid, ip );

            // Calculate the WFs.
            // A is di/dkappa*dkappa/dkp in a compact form.
            // This is possible as the columns of dkappa/dkp are identical.  

            a = 0.0;            // Matpack can set all elements like this.

            i1 = (Index)is(ip,0);       // first LOS point to consider
            for ( iv=0; iv<nv; iv++ )
            {
              for ( iw=i1; iw<=(Index)is(ip,1); iw++ )
                a[iv] += absloswfs[iza](iv,iw) * w[iw-i1];
              k(iv0+iv,ip0+ip) = a[iv] * abs(iv,ip);                    
            }
          }            
        }
      }

      // Update the frequency index offset
      iv0 += nv;
    }  
    out3 << "\n";

    // Increase retrieval altitude index off-set 
    ip0 += np;
  }  
}




void k_contabs (
                    Matrix&          k,
                    ArrayOfString&   k_names,
                    Matrix&          k_aux,
              const Los&             los,           
              const ArrayOfMatrix&   absloswfs,
              const Vector&          f_mono,
              const Vector&          k_grid,
              const Index&           order,
              const Numeric&         flow,
              const Numeric&         fhigh,
              const String&          lunit )
{

  if ( los.p.nelem() != absloswfs.nelem() )
    {
    throw runtime_error(
                    "The number of zenith angles is not the same in *los* and *absloswfs*." ); 
    check_length_nrow( los.p[0], "los points", absloswfs[0], 
                                                 "the matrices of absloswfs" );
    }
  
  // Main sizes
  const Index  nza = los.start.nelem();     // number of zenith angles  
  const Index  np  = k_grid.nelem();        // number of retrieval altitudes
  const Index  npoints = order+1;           // number of off-set points
  const Index  nv  = f_mono.nelem();        // number of frequencies

  // Check given frequency limits
  if ( flow >= fhigh )
    throw runtime_error(
             "The lower frequency limit equals or is above the upper limit." );
  if ( flow < f_mono[0] )
    throw runtime_error(
                  "The lower frequency limit is below all values of f_mono." );
  if ( fhigh > f_mono[nv-1] )
    throw runtime_error(
                  "The upper frequency limit is above all values of f_mono." );

  // Check length unit
  //
  Numeric scfac;
  //
  if ( lunit == "m" )
    scfac = 1.0;
  else if ( lunit == "km" )
    scfac = 0.001;
  else
    throw runtime_error("Allowed length units are \"m\" and \"km\".");

  // -log(p) is used as altitude variable. The minus is included to get
  // increasing values, a demand for the grid functions. 
  Vector  lgrid;
  p2grid( lgrid, k_grid );
  Vector  lplos;                      // -log of the LOS pressures

  // Indices
  // IP0 and IF0 are the index off-sets for the total K matrix
        Index  ipoint;                    // Off-set point index
        Index  iza;                       // Zenith angle index
        Index  ip, ip0=0;                 // Retrieval point indices
        Index  iv, iv0;                   // Frequency indices
        Index  i1, iw;                    // weight indices

  // Determine first and last frequency index inside given limits
  Index   ilow, ihigh;
  for( ilow=0; ilow<(nv-1) && f_mono[ilow] < flow; ilow++ )
    {}
  for( ihigh=ilow; ihigh<(nv-1) && f_mono[ihigh+1] <= fhigh; ihigh++ )
    {}

  // Other variables
        Vector  fpoints;                   // frequencies of the off-set points
        Vector  b(nv);                     // fit base function
        Matrix  is;                        // matrix for storing LOS index
        Vector  w;                         // weights for LOS WFs
        Vector  a(nv);                     // temporary vector

  // Check that the selected polynomial order
  if ( order < 0 )
    throw runtime_error("The polynomial order must be >= 0."); 

  // Set up K and additional data. Set all values of K to 0
  k.resize(nza*nv,npoints*np);
  k = 0.0;                      // Matpack can set all elements like this.
  k_names.resize(npoints);
  k_aux.resize(npoints*np,2);

  // Calculate the frequencies of the off-set points
  if ( npoints > 1 )
    nlinspace( fpoints, f_mono[ilow], f_mono[ihigh], npoints );
  else
  {
    fpoints.resize( 1 );
    fpoints[0] = ( f_mono[ilow] + f_mono[ihigh] ) / 2.0;
  }  

  // The calculations
  // Loop order:
  //   1 polynomial order
  //   2 zenith angle
  //   3 retrieval altitudes 
  //   4 frequencies
  //
  // (Another loop order should be somewhat more efficient, but to keep this
  // function consistent with k_species, this loop order was selected.)
  //
  out2 << "  You have selected " << npoints << " off-set fit points.\n";
  out2 << "  Length unit is " << lunit << "\n";
  for ( ipoint=0; ipoint<npoints; ipoint++ ) 
  {
    out2 << "  Doing point " << ipoint << "\n";

    // Fill K_NAMES and K_AUX
    {
      ostringstream os;
      os << "Continuum: " << fpoints[ipoint]/1e9 << " GHz, Point " 
         << ipoint+1 << "/" << npoints;
      k_names[ipoint] = os.str();
    }
    for ( ip=0; ip<np; ip++ )
    {
       k_aux(ip0+ip,0) = k_grid[ip];
       k_aux(ip0+ip,1) = 0.0;
    }

    // Set-up base vector for the present fit point 
    b = 1.0;                    // Matpack can set all elements like this.
    if ( npoints > 1 )
    {
      for ( ip=0; ip<npoints; ip++ )
      {
        if ( ip != ipoint )
        { 
          for ( iv=ilow; iv<=ihigh; iv++ )
            b[iv] *= (f_mono[iv]-fpoints[ip]) / ( fpoints[ipoint]-fpoints[ip]);
        }
      }
    }

    // Set frequency zenith angle index off-set to 0
    iv0 = 0;                 

    // Loop zenith angles
    out3 << "    Zenith angle nr:\n      ";
    for ( iza=0; iza<nza; iza++ ) 
    {
      if ( ((iza+1)%20)==0 )
        out3 << "\n      ";
      out3 << " " << iza; cout.flush();
      
      // Do something only if LOS has any points
      if ( los.p[iza].nelem() > 0 )
      {
        // Get the LOS points affected by each retrieval point
        //        lplos = -1.0 * log(los.p[iza]);
        p2grid( lplos, los.p[iza] );
        grid2grid_index( is, lplos, lgrid );

        // Loop retrieval points
        for ( ip=0; ip<np; ip++ ) 
        {
          // Check if there is anything to do
          if ( is(ip,0) >= 0 )
          {
            // Get the weights for the LOS points
            grid2grid_weights( w, lplos, (Index)is(ip,0), (Index) is(ip,1),
                                                                 lgrid, ip );

            // Calculate the WFs.
            // A is di/dkappa*dkappa/dkp in a compact form.
            // This is possible as the columns of dkappa/dkp are identical.  

            a = 0.0;            // Matpack can set all elements like this.

            i1 = (Index)is(ip,0);       // first LOS point to consider
            for ( iv=ilow; iv<=ihigh; iv++ )
            {
              for ( iw=i1; iw<=(Index)is(ip,1); iw++ )
                a[iv] += absloswfs[iza](iv,iw) * w[iw-i1];
              k(iv0+iv,ip0+ip) = scfac * a[iv] * b[iv];                    
            }
          }            
        }
      }

      // Update the frequency index offset
      iv0 += nv;
    }  
    out3 << "\n";

    // Increase retrieval altitude index off-set 
    ip0 += np;
  }  
}




void k_temp_nohydro (
              Matrix&                     k,
              ArrayOfString&              k_names,
              Matrix&                     k_aux,
        const TagGroups&                  tag_groups,
        const Los&                        los,           
        const ArrayOfMatrix&              absloswfs,
        const Vector&                     f_mono,
        const Vector&                     p_abs,
        const Vector&                     t_abs,
        const Vector&                     n2_abs,          
        const Vector&                     h2o_abs,         
        const Matrix&                     vmrs,
        const ArrayOfArrayOfLineRecord&   lines_per_tg,
        const ArrayOfLineshapeSpec&       lineshape,
        const Matrix&                     abs,            
        const ArrayOfMatrix&              trans,
        const Vector&                     e_ground,
        const Vector&                     k_grid,
        const ArrayOfString&              cont_description_names,
        const ArrayOfVector&              cont_description_parameters,
        const ArrayOfString&              cont_description_models )
{
  // Main sizes
  const Index  nza = los.start.nelem();     // number of zenith angles  
  const Index  nv  = f_mono.nelem();        // number of frequencies
  const Index  np  = k_grid.nelem();        // number of retrieval altitudes

  // -log(p) is used as altitude variable. The minus is included to get
  // increasing values, a demand for the grid functions. 
  //  const Vector  lgrid=-1.0*log(k_grid);

  Vector  lgrid;
  p2grid( lgrid, k_grid );
  Vector  lplos;                     // -log of the LOS pressures
  ArrayOfMatrix W;
  Index ip;
  //
  // Indices
  // IP0 and IF0 are the index off-sets for the total K matrix
        Index  iza;                       // zenith angle index
        //Index  ip;                        // retrieval point index
        Index  iv, iv0=0;                 // frequency indices
        Index  i1, iw;                    // weight indices

  // Other variables
        Vector  t(k_grid.nelem());          // temperature at retrieval points
        Matrix  abs1k;                     // absorption for t_abs+1K
        Matrix  dabs_dt;                   // see below
 ArrayOfMatrix  abs_dummy;                 // dummy absorption array
 ArrayOfMatrix  sloswfs;                   // source LOS WFs
        Matrix  is;                        // matrix for storing LOS index
        Vector  w;                         // weights for LOS WFs
        Vector  a(nv), b(nv), pl(f_mono.nelem());  // temporary vectors

 // The scalars are declared to be double to avoid possible numerical problems
 // when using float
        double  c,d;                       // temporary values


  // Set up K and additional data. Set all values of K to 0
  k.resize(nza*nv,np);
  k = 0.0;                      // Matpack can set all elements like this.
  k_names.resize(1);
  k_names[0] = "Temperature: no hydrostatic eq.";
  k_aux.resize(np,2);
  interpp( t, p_abs, t_abs, k_grid ); 
  for ( ip=0; ip<np; ip++ )
  {
     k_aux(ip,0) = k_grid[ip];
     k_aux(ip,1) = t[ip];
  }

  // Calculate absorption for t_abs + 1K to estimate the temperature derivative
  // dabs/dt, the temperature derivative of the absorption at k_grid
  out2 << "  Calculating absorption for t_abs + 1K\n";
  out2 << "  ----- Messages from absCalc: -----\n";
  //
  {
    // Dummy should hold t_abs + 1:
    Vector dummy(t_abs);        // Matpack can initialize a Vector
                                // from another Vector. 
    dummy += 1;                 // Matpack can add element-vise like this.

    absCalc( abs1k, abs_dummy, tag_groups, f_mono, p_abs, dummy, n2_abs, 
             h2o_abs, vmrs, lines_per_tg, lineshape, cont_description_names, 
                        cont_description_models, cont_description_parameters);
  }
  abs_dummy.resize(0);
  //
  out2 << "  ----- Back from absCalc ----------\n";
  //
  // Compute abs1k = abs1k - abs:
  abs1k -= abs;                 // Matpack can subtract element-vise like this.

  dabs_dt.resize( abs1k.nrows(), k_grid.nelem() );
  interpp( dabs_dt, p_abs, abs1k, k_grid );
  abs1k.resize(0,0);

  // Calculate source LOS WFs
  out2 << "  Calculating source LOS WFs\n";
  sourceloswfs( sloswfs, los, trans, f_mono, e_ground );

  // Determine the temperatures at the retrieval points
  out2 << "  Calculating temperature at retrieval points\n";
  interpp( t, p_abs, t_abs, k_grid );

  // The calculations
  // Loop order:
  //   1 zenith angle
  //   2 retrieval altitudes 
  //   3 frequencies
  //
  out2 << "  Calculating the weighting functions\n";
  out3 << "    Zenith angle nr:      ";
  for ( iza=0; iza<nza; iza++ ) 
  {
    if ( (iza%20)==0 )
      out3 << "\n      ";
    out3 << " " << iza; cout.flush();
    
    // Do something only if LOS has any points
    if ( los.p[iza].nelem() > 0 )
    {
      // Get the LOS points affected by each retrieval point
      //      lplos = -1.0 * log(los.p[iza]);
      p2grid( lplos, los.p[iza] );
      grid2grid_index( is, lplos, lgrid );

      // Loop retrieval points
      for ( ip=0; ip<np; ip++ ) 
      {
        // Check if there is anything to do
        if ( is(ip,0) >= 0 )
         {
          // Get the weights for the LOS points
          grid2grid_weights( w, lplos, (Index)is(ip,0), (Index) is(ip,1), 
                                                                 lgrid, ip );

          // Calculate the WFs.
          // A is di/dsigma*dsigma/dSp in a compact form.
          // B is di/dkappa*dkappa/dkp in a compact form.
          // This is possible as the columns of dkappa/dkp are identical and  
          // that dkappa/dkp = dsigma/dSp
          // C is just a temporary value
          // PL is the Planck function for the present temperature value
          //
          a = 0.0;              // Matpack can set all elements like this.
          b = 0.0;                    
          c  = PLANCK_CONST / BOLTZMAN_CONST / t[ip];
          planck( pl, f_mono, t[ip] );
          i1 = (Index)is(ip,0);       // first LOS point to consider
          //
          for ( iv=0; iv<nv; iv++ )
          {
            for ( iw=i1; iw<=(Index)is(ip,1); iw++ )
            {
              a[iv] += sloswfs[iza](iv,iw) * w[iw-i1];
              b[iv] += absloswfs[iza](iv,iw) * w[iw-i1];
            }
            d = c * f_mono[iv];
            k(iv0+iv,ip) = a[iv] * d/t[ip] / (exp(d)-1) * pl[iv] +
                                                      b[iv] * dabs_dt(iv,ip);
            
            // The result for very small values can be NaN. So NaNs (and Infs)
            // are set to 0.
#ifdef HPUX
            if ( !isfinite(k(iv0+iv,ip)) )
#else
            if ( !finite(k(iv0+iv,ip)) )
#endif
              k(iv0+iv,ip) = 0.0;

          }
        }            
      }
    }

    // Update the frequency index offset
    iv0 += nv;
  }  
   out3 << "\n";
}


void kb_spectro(
                Matrix&                             kb,
                ArrayOfString&                      kb_names,
                Matrix&                             /* kb_aux */,
                Matrix&                             S_S,
          const TagGroups&                          wfss_tgs,
          const TagGroups&                          tgs,
          const Vector&                             f_mono,
                const Vector&                             p_abs,
                const Vector&                             t_abs,
                const Vector&                             h2o_abs,
                const Matrix&                             vmrs,
                const ArrayOfArrayOfLineRecord&     lines_per_tg,
                const ArrayOfLineshapeSpec&   lineshape,
                const Los&                                 los,           
          const ArrayOfMatrix&                absweight,
          const Index&                              IndexPar,
          const String&                             StrPar
)
 { 
 // Main sizes
  const Index  nza = los.start.nelem();     // number of zenith angles  
  const Index  nv  = f_mono.nelem();        // number of frequencies
  const Index  np  = p_abs.nelem();         // number of retrieval altitudes
  //----------------------------------------------------------------------
  // Index:
  Index itg, iv, nr_line_total, iza, ip, nr_line; 
  //----------------------------------------------------------------------
 // Other variables: 
  Matrix        abs_line;
  Matrix        abs_line_changed;
  Matrix        dabs;
  dabs.resize(np,nv);
  Matrix        dabs1;
  dabs1.resize(nv,np);

  //dabs.resize(nv, np);
  ArrayOfMatrix abs_dummy1;
  ArrayOfMatrix abs_dummy2;
  Index itg_wfss;
  ArrayOfIndex tag_index;
  tag_index.resize(wfss_tgs.nelem()); 
  tag_index = 0;
  //Index  dummy_line=0;

 // Go through the spectroscopic parameters weighting function
  for (itg_wfss=0; itg_wfss<wfss_tgs.nelem(); ++itg_wfss)
    {
      Index n;
      get_tag_group_index_for_tag_group( n, tgs, wfss_tgs[itg_wfss] );
      tag_index[itg_wfss] =n;
    }

 // Get size for Kb
  nr_line_total=0;
  nr_line=0;
  if (tag_index.nelem()==0)
    {
      ostringstream os;
      os << "No species has been set:   "<<"\n";        
      throw runtime_error(os.str());
    }
  else
    {  
      for ( itg_wfss=0; itg_wfss<tag_index.nelem(); ++itg_wfss )
        {
          itg=tag_index[itg_wfss];
          nr_line_total+=lines_per_tg[itg].nelem();
        }  
    }
  kb.resize(nza*nv,nr_line_total);
  kb=0.0;
  kb_names.resize(nr_line_total);
  S_S.resize(nr_line_total,2);
  S_S=0.0;
  
   // loop:
   // 1: over the tags
   // 2: over the lines
   // 3: over the zenith angles
   // 4: over the frequencies

     for ( itg_wfss=0; itg_wfss<tag_index.nelem(); ++itg_wfss )
       {  
        out2<< " \n"; 
        out2<<"==================================: " <<" \n"; 
        out2<<"==Calculating weighting function for species==: " <<" \n"; 
        out2<< wfss_tgs[itg_wfss]<< "\n";  
        itg=tag_index[itg_wfss];
        out2<<"lines found: " <<lines_per_tg[itg].nelem()<< " \n"; 
        Index nr_line_dummy=0;
         if ( lines_per_tg[itg].nelem()>0)
           {
             Vector gamma;
             gamma.resize(lines_per_tg[itg].nelem());
             gamma=0.0;
             for ( Index li=0; li<lines_per_tg[itg].nelem(); ++li )
               { 
                 ArrayOfArrayOfLineRecord dummy_lines_per_tg; 
                 dummy_lines_per_tg.resize( lines_per_tg.nelem() ); 
                 dummy_lines_per_tg[itg].resize(1);
                 out3<<"==Calculating weighting function for line==: " <<" \n"; 
                 out3<< lines_per_tg[itg][li]<< "\n"; 
                 dummy_lines_per_tg[itg][0] = lines_per_tg[itg][li];
                 // calculate the absorption for unperturbed line; 
                 xsec_per_tgInit( abs_dummy1, tgs, f_mono, p_abs );
                 xsec_per_tgAddLines( abs_dummy1, tgs, f_mono, p_abs, t_abs,
                                      h2o_abs, vmrs, dummy_lines_per_tg, 
                                      lineshape );               
                 absCalcFromXsec(abs_line, abs_dummy2, abs_dummy1 , vmrs );
                // define some local variables      
                 Numeric delta_s;   // perturbation applied
                 Numeric parameter; // spectro parameter
                 Numeric dummy;     // perturb spectro parameter
                 Numeric ER_dummy;  // accuracy from spectral  database
                 Vector ER_dummy2;  // accuracy in absolute and relative units;
                 ER_dummy2.resize(2);
                 ER_dummy2=0.0;
                // apply a perturbation to one  parameter
                 switch ( IndexPar )
                   {
                   case 1: 
                     {
                       delta_s = 10E-25; //perturbation in intensity units;
                       // read the intensity and its uncertainty;
                       parameter  = dummy_lines_per_tg[itg][0].I0();
                       ER_dummy   = dummy_lines_per_tg[itg][0].dI0();
                       // Set to a default value if there is no field for uncertainty;
                       if (ER_dummy == -1) //
                         {ER_dummy = 0.02;}
                       // get uncertainty in (1) absolute values and (2) relative values;
                       ER_dummy2[0] = ER_dummy*parameter;
                       ER_dummy2[1] = ER_dummy*100;
                       // apply uncertainty and update the value of the parameter;
                       dummy = parameter +delta_s;       
                       dummy_lines_per_tg[itg][0].setI0(dummy);
                       break;
                     }
                   case 2: 
                     {
                       delta_s=10; // perturbation in HZ
                                  // read the position and and its uncertainty; 
                       parameter=dummy_lines_per_tg[itg][0].F();
                       ER_dummy = dummy_lines_per_tg[itg][0].dF();
                                 // Set to a default value if there is no field for uncertainty;
                       if (ER_dummy==-1)
                         {ER_dummy=3*1E+9;}
                                    // get uncertainty in absolute values only; 
                       ER_dummy2[0] = ER_dummy;
                       ER_dummy2[1] = ER_dummy;
                                  // apply uncertainty and update the value of the parameter;
                       dummy = parameter +delta_s;       
                       dummy_lines_per_tg[itg][0].setF(dummy);
                       break;
                     }
                   case 3: 
                     {  
                       delta_s=0.001; // perturbation in fraction
                                  // read the agam and and its uncertainty; 
                       parameter = dummy_lines_per_tg[itg][0].Agam();                 
                       ER_dummy = dummy_lines_per_tg[itg][0].dAgam();
                                    // Set to a default value if there is no field for uncertainty;
                       if (ER_dummy==-1)
                         {ER_dummy=2;}
                                   // get uncertainty both in (1) absolute values and (2) relative values;
                       ER_dummy2[0] = parameter*ER_dummy;
                       ER_dummy2[1] = ER_dummy*100;
                                   // apply uncertainty and update the value of the parameter;
                       dummy = parameter +delta_s;       
                       dummy = parameter + parameter *delta_s;   
                       dummy_lines_per_tg[itg][0].setAgam(dummy);
                       break;
                     }
                   case 4: 
                     {
                       delta_s=0.001; // perturbation in fraction
                                   // read the sgam and and its uncertainty; 
                       parameter=dummy_lines_per_tg[itg][0].Sgam();
                       ER_dummy = dummy_lines_per_tg[itg][0].dSgam();
                                   // Set to a default value if there is no field for uncertainty;
                       if (ER_dummy==-1)
                         {ER_dummy=2;}
                                   // get uncertainty both in (1) absolute values and (2) relative values; 
                       ER_dummy2[0] = parameter*ER_dummy;
                       ER_dummy2[1] = ER_dummy*100;
                                  // apply uncertainty and update the value of the parameter;
                       dummy    = parameter +parameter*delta_s;          
                       dummy_lines_per_tg[itg][0].setSgam(dummy);
                       break;
                     }
                   case 5: 
                     {
                       delta_s=0.001; // perturbation in fraction
                                 // read the Nair and and its uncertainty;  
                       parameter=dummy_lines_per_tg[itg][0].Nair();
                       ER_dummy = dummy_lines_per_tg[itg][0].dNair();
                                  // Set to a default value (200%) if there is no field for uncertainty;
                       if (ER_dummy==-1)
                         {ER_dummy=2;}
                       // get uncertainty both in (1) absolute values and (2) relative values; 
                       ER_dummy2[0] = parameter*ER_dummy;
                       ER_dummy2[1] = ER_dummy*100;
                                   // apply uncertainty and update the value of the parameter;
                       dummy = parameter + parameter*delta_s;            
                       dummy_lines_per_tg[itg][0].setNair(dummy);
                       break; 
                     }
                   case 6: 
                     { 
                       delta_s=0.01; // perturbation in fraction
                                  // read the Nair and and its uncertainty; 
                       parameter=dummy_lines_per_tg[itg][0].Nself();
                       ER_dummy = dummy_lines_per_tg[itg][0].dNself();
                                   // Set to a default value (200%) if there is no field for uncertainty;
                       if (ER_dummy==-1)
                         {ER_dummy=2;}
                                   // get uncertainty both in (1) absolute values and (2) relative values; 
                       ER_dummy2[0] = parameter*ER_dummy;
                       ER_dummy2[1] = ER_dummy*100;
                                    // apply uncertainty and update the value of the parameter;
                       dummy = parameter +parameter*delta_s;             
                       dummy_lines_per_tg[itg][0].setNself(dummy);
                       break;
                     }
                   case 7: 
                     {
                       extern const Numeric TORR2PA;
                       delta_s=1;// perturbation i
                                  // read the Nair and and its uncertainty; 
                       parameter=dummy_lines_per_tg[itg][0].Psf();
                       ER_dummy = dummy_lines_per_tg[itg][0].dPsf();
                                    // Set to a default value (1MHz/Torr) if there is no field for uncertainty;
                       if ((ER_dummy==-1)|| (parameter==0.0))
                         {ER_dummy = 1E6 / TORR2PA;}
                       else 
                         {ER_dummy=ER_dummy*parameter;}
                       // apply uncertainty and update the value of the parameter;
                       ER_dummy2[0] = ER_dummy;
                       ER_dummy2[1] = ER_dummy;
                       dummy = parameter + delta_s                 ;     
                       dummy_lines_per_tg[itg][0].setPsf(dummy);         
                       break; 
                     }
                   default: 
                     {
                       ostringstream os;
                       os << "The sectroscopic parameter:   "  << StrPar<<"\n"
                          << "does not exists"; 
                       throw runtime_error(os.str());
                     }
                   }
                 //calculate the absorption coefficient for the perturbed line
                 xsec_per_tgInit( abs_dummy1, tgs, f_mono, p_abs );
                 xsec_per_tgAddLines( abs_dummy1,
                                      tgs, f_mono, p_abs, t_abs, h2o_abs, vmrs,
                                      dummy_lines_per_tg, lineshape );  
                 absCalcFromXsec(abs_line_changed, abs_dummy2, abs_dummy1, vmrs ); 
                 abs_line_changed-=abs_line; 
                 dabs1=abs_line_changed;
                 dabs1/=dummy-parameter;
                 dabs=transpose(dabs1);
                 Index iv0=0;
                 for ( iza=0; iza<nza; iza++ ) 
                   { 
                     if ( (iza%20)==0 )
                       out3 << "\n      ";
                     out3 << " " << iza; cout.flush();
                     // Do something only if LOS has any points
                     if ( los.p[iza].nelem() > 0 ) 
                       {
                         for ( iv=0; iv<nv; iv++ ) 
                           {         
                             for ( ip=0; ip<np; ip++ )
                               { 
                                 kb(iv0+iv, nr_line)+= absweight[iza](iv,ip)* dabs(ip, iv);
                               }
                             // The result for very small values can be NaN. So NaNs (and Infs)
                             // are set to 0.
#ifdef HPUX
                             if ( !isfinite(kb(iv0+iv, nr_line)) )
#else
                               if ( !finite(kb(iv0+iv, nr_line)) )
#endif
                                 kb(iv0+iv,nr_line) = 0.0;
                           }
                         iv0+=nv;
                       }
                   }        
                 Numeric nr = dummy_lines_per_tg[itg][0].F() *1e-9;
                 ostringstream os;
                 // create the sting for the name parameter (written in kb_name)
                 // line (species and center frequency)+ uncertinity + name parameter
                 if (IndexPar == 2)
                   {
                     os <<tgs[itg]<<"@" <<nr<<"GHz"<<" / "<<StrPar<<": " << ER_dummy2[1] <<"Hz";
                   }
                 else if (IndexPar == 7)
                   {
                     os <<tgs[itg]<<"@" <<nr<<"GHz"<<" / "<<StrPar<<": "  << ER_dummy2[1] <<"Hz/Pa";
                   }

                 else
                   {
                     os <<tgs[itg]<<"@" <<nr<<"GHz"<<" / "<<StrPar<<": "  << ER_dummy2[1]<< "%";
                   }
                 kb_names[nr_line]= os.str();
                 S_S(nr_line, Range(joker))= ER_dummy2;
                 ++nr_line;
                 ++nr_line_dummy;
               } 
       }
 }
 }

void kSpectro (
                    Matrix&                               k,
                    ArrayOfString&                   k_names,
                    Matrix&                               /* kb_aux */,
                    Matrix&                               S_S,
          const TagGroups&                          wfss_tgs,
          const TagGroups&                          tgs,
          const Vector&                               f_mono,
          const Vector&                               p_abs,
          const Vector&                               t_abs,
          const Vector&                               z_abs,
          const Vector&                               h2o_abs,
          const Matrix&                               vmrs,
          const ArrayOfArrayOfLineRecord&             lines_per_tg,
          const ArrayOfLineshapeSpec&                    lineshape,
          const Los&                                  los,           
          const ArrayOfMatrix&                 absloswfs,
          // Keywords
          const  Index&                               kw_intens,
          const  Index&                               kw_position,
                const  Index&                               kw_agam,
                const  Index&                               kw_sgam,
                const  Index&                               kw_nair,
          const  Index&                               kw_nself,
                const  Index&                               kw_pSift)
{

  check_if_bool( kw_intens, "do_intens keyword" );
  check_if_bool( kw_position, "do_position keyword" );
  check_if_bool( kw_agam, "do_agam keyword" );
  check_if_bool( kw_sgam, "do_sgam keyword" );
  check_if_bool( kw_nair, "do_nair keyword" );
  check_if_bool( kw_nself, "do_nself keyword" );
  check_if_bool( kw_pSift,  "do_pSift keyword" );
  check_lengths( p_abs, "p_abs", t_abs, "t_abs" );  
  check_lengths( p_abs, "p_abs", z_abs, "z_abs" );  
  check_lengths( p_abs, "p_abs", h2o_abs, "h2o_abs" );  
  check_length_ncol( p_abs, "p_abs", vmrs, "vmrs" );

 // Main sizes
  const Index  nza = los.start.nelem();     // number of zenith angles  
  const Index  nv  = f_mono.nelem();        // number of frequencies
  //
  // Indices
  // IP0 and IF0 are the index off-sets for the total K matrix
  Index  iza;                       // zenith angle index         
  Index itg;
  
  // Other variables
  ArrayOfMatrix W; // keeps  dkappa/dk;
  ArrayOfMatrix absweight; // the product di/dkappa * dkappa/dk;
  absweight.resize(nza);
  ArrayOfIndex  k_lengths;
  Matrix        kb;
  ArrayOfString kb_names;
  Matrix        kb_aux;
  Matrix        abs_line;
  Matrix        abs_line_changed;
  Matrix        ER;
  ArrayOfMatrix abs_dummy1;
  Index         IndexPar; // spectro parameters to be investigated
  String        StrPar;
  Index         ipar = 0;              // index for the parameter

  // Calculate weights for conversion between LOS grids and pressure grid: dkappa/dk;
  grid2grid_weights_total(W, los, p_abs);
  // Calculate the product di/dkappa * dkappa/dk; 
  for ( iza=0; iza<nza; iza++ ) 
    {
      absweight[iza].resize(absloswfs[iza].nrows(), W[iza].ncols());  
      // Do something only if LOS has any points
      if ( los.p[iza].nelem() > 0 )
        {
          mult(absweight[iza], absloswfs[iza], W[iza]); 
        }
    }

// get size of k
  if (kw_intens) 
    {  ++ipar;}
  if (kw_position) 
    {  ++ipar;}
  if (kw_agam) 
    {  ++ipar;}
  if (kw_sgam) 
    {  ++ipar;}
  if (kw_nair) 
    {  ++ipar;}
  if (kw_nself) 
    {  ++ipar;}
  if (kw_pSift) 
    {  ++ipar;}
  Index nr_line_total=0;
  Index itg_wfss;
  ArrayOfIndex tag_index;
  tag_index.resize(wfss_tgs.nelem()); 
  tag_index = 0;


  // Go through the spectroscopic parameters weighting function
  for (itg_wfss=0; itg_wfss<wfss_tgs.nelem(); ++itg_wfss)
    {
      Index n;
      get_tag_group_index_for_tag_group( n, tgs, wfss_tgs[itg_wfss] );
      tag_index[itg_wfss] =n;
    }

  // Get size for Kb
  nr_line_total=0;

  // Calculate the total number of lines for which the weighting
  // function is calculated. In caz no speciaes are specified in wtgss
  // then an error messege is given.
 if (tag_index.nelem()==0)
    {
      ostringstream os;
      os << "No species has been set:   "<<"\n";        
      throw runtime_error(os.str());
    }
  else
    {  
      for ( itg_wfss=0; itg_wfss<tag_index.nelem(); ++itg_wfss )
        {
          itg=tag_index[itg_wfss];
          nr_line_total+=lines_per_tg[itg].nelem();
        }  
    }

 // Check whether the spectroscopic parameters for which to calculate the
 // weighting function are set or lines in the spectral data are found.
  if (ipar == 0)
    {
      ostringstream os;
      os << "No sectroscopic parameters have been set"; 
      throw runtime_error(os.str());
    }
  if (nr_line_total == 0)
    {
      ostringstream os;
      os << " No line for which to calculate weighting function has been found. Catalog empty?";        
      throw runtime_error(os.str());
    }
  
  k.resize(nza*nv, nr_line_total*ipar);
  k=0.0;
  k_names.resize(nr_line_total*ipar);
  S_S.resize(nr_line_total*ipar,2);
  ER.resize(nr_line_total,2);

  // Go through the spectroscopic parameters and calculate the weighting function.
  Index nr_line=0;
 
  // for intemsity
  if (kw_intens) 
    {  
      IndexPar = 1;
      StrPar = "S";
      out1<< " \n"; 
      out1<<"==================================: " <<" \n"; 
      out1<<" ******* Calculating Wfs for "<< StrPar<<" ******\n"; 
      kb_spectro( kb, kb_names, kb_aux, ER, wfss_tgs, tgs, 
                  f_mono, p_abs, t_abs, 
                  h2o_abs, vmrs, lines_per_tg, lineshape, los, absweight, 
                  IndexPar, StrPar); 
      k(Range(joker),Range(nr_line,  kb.ncols()))= kb;
     // ER.resize(kb_names.nelem(),2);
      for ( Index iri=0; iri<kb_names.nelem(); iri++)
        {
          k_names[nr_line+iri]= kb_names[iri]; 
          S_S(nr_line+iri, Range(joker))=ER(iri, Range(joker));
        } 
        nr_line+=kb.ncols();
    }
  // for line position  
  if (kw_position)
    { 
      IndexPar  = 2;
      StrPar = "f0";
      out1 <<" ******* Calculating Wfs for "<< StrPar<<" ******\n"; 
      kb_spectro( kb, kb_names, kb_aux, ER, wfss_tgs, tgs, 
                  f_mono, p_abs, t_abs, 
                  h2o_abs, vmrs, lines_per_tg, 
                  lineshape, los, absweight, IndexPar, StrPar);
      k(Range(joker), Range(nr_line,  kb.ncols()))= kb;
      // ER.resize(kb_names.nelem(),2);
      for ( Index iri=0; iri<kb_names.nelem(); iri++)
        {
          k_names[nr_line+iri]= kb_names[iri];
          S_S(nr_line+iri, Range(joker))=ER(iri, Range(joker));
        } 
              nr_line+=kb.ncols();
        }
  // for air broadening parameter agam
  if (kw_agam)
    { 
      IndexPar = 3;
      StrPar = "agam";
      out1<<" ******* Calculating Wfs for "<< StrPar<<" ******\n"; 
      kb_spectro( kb, kb_names, kb_aux, ER, wfss_tgs,
                  tgs, f_mono, p_abs, t_abs,
                  h2o_abs, vmrs, lines_per_tg, 
                  lineshape, los, absweight, IndexPar, StrPar);
      k(Range(joker),Range(nr_line,  kb.ncols()))= kb;
      //ER.resize(kb_names.nelem(),2);
      for ( Index iri=0; iri<kb_names.nelem(); iri++)
        {
          k_names[nr_line+iri]= kb_names[iri]; 
          S_S(nr_line+iri, Range(joker))=ER(iri, Range(joker));
        } 
             nr_line+=kb.ncols();
        }
// for self broadening parameter sgam
  if (kw_sgam)
    { 
      IndexPar = 4;
      StrPar = "sgam";
      out1 <<" ******* Calculating Wfs for "<< StrPar<<" ******\n"; 
      kb_spectro( kb, kb_names, kb_aux, ER,  wfss_tgs,
                  tgs, f_mono, p_abs, t_abs, 
                  h2o_abs, vmrs, lines_per_tg, 
                  lineshape, los, absweight, IndexPar, StrPar);
      k(Range(joker),Range(nr_line,  kb.ncols()))= kb;
      //ER.resize(kb_names.nelem(),2);
      for ( Index iri=0; iri<kb_names.nelem(); iri++)
        {
          k_names[nr_line+iri]= kb_names[iri]; 
          S_S(nr_line+iri, Range(joker))=ER(iri, Range(joker));
        } 
      nr_line+=kb.ncols();
    }
 // for air broadening parameter nair
  if (kw_nair)
    { 
      IndexPar = 5;
      StrPar = "na";
      out1 <<" ******* Calculating Wfs for "<< StrPar<<" ******\n"; 
      kb_spectro( kb, kb_names, kb_aux, ER, wfss_tgs, 
                  tgs, f_mono, p_abs, t_abs, 
                  h2o_abs, vmrs, lines_per_tg, 
                  lineshape, los, absweight, IndexPar, StrPar);
      k(Range(joker),Range(nr_line,  kb.ncols()))= kb;
      //ER.resize(kb_names.nelem(),2);
      for ( Index iri=0; iri<kb_names.nelem(); iri++)
        {
          k_names[nr_line+iri]= kb_names[iri];
          S_S(nr_line+iri, Range(joker)) = ER(iri, Range(joker));
        } 
      nr_line+=kb.ncols();      
        }
// for self broadening parameter nself
  if (kw_nself)
    { 
      IndexPar = 6;
      StrPar = "ns";
      out1 <<" ******* Calculating Wfs for "<< StrPar<<" ******\n"; 
      kb_spectro( kb, kb_names, kb_aux, ER, wfss_tgs, 
                  tgs, f_mono, p_abs, t_abs, 
                  h2o_abs, vmrs, lines_per_tg, lineshape, los, absweight, 
                  IndexPar, StrPar);
      k(Range(joker),Range(nr_line,  kb.ncols()))= kb;
      // ER.resize(kb_names.nelem(),2); 
      for ( Index iri=0; iri<kb_names.nelem(); iri++)
        {
          k_names[nr_line+iri]= kb_names[iri];
          S_S(nr_line+iri, Range(joker))=ER(iri, Range(joker)); 
        } 
      nr_line+=kb.ncols();        
        }
  // for pressure shift
  if (kw_pSift)
    { 
      IndexPar = 7;
      StrPar = "pSf";
      out1 <<" ******* Calculating Wfs for "<< StrPar<<" ******\n"; 
      kb_spectro( kb, kb_names, kb_aux, ER, wfss_tgs, 
                  tgs, f_mono, p_abs, t_abs,
                  h2o_abs, vmrs, lines_per_tg, lineshape, los, absweight,
                  IndexPar, StrPar);
      k(Range(joker),Range(nr_line,  kb.ncols()))= kb;
      ER.resize(kb_names.nelem(),2);
      for ( Index iri=0; iri<kb_names.nelem(); iri++)
        {
          k_names[nr_line+iri]= kb_names[iri]; 
          S_S(nr_line+iri, Range(joker)) = ER(iri, Range(joker));
        } 
      nr_line+=kb.ncols();
        }
  if (ipar == 0)
     {
       ostringstream os;
       os << "No sectroscopic parameters have been set";        
       throw runtime_error(os.str());
     }
} 
 
//   Workspace methods


void wfs_tgsDefine(// WS Output:
                   TagGroups& wfs_tag_groups,
                   // Control Parameters:
                   const ArrayOfString& tags)
{
  wfs_tag_groups = TagGroups(tags.nelem());

  // Each element of the array of Strings tags defines one tag
  // group. Let's work through them one by one.
  for ( Index i=0; i<tags.nelem(); ++i )
    {
      // There can be a comma separated list of tag definitions, so we
      // need to break the String apart at the commas.
      ArrayOfString tag_def;

      bool go_on = true;
      String these_tags = tags[i];
      while (go_on)
        {
          Index n = these_tags.find(',');
          if ( n == these_tags.npos ) // npos indicates `not found'
            {
              // There are no more commas.
              tag_def.push_back(these_tags);
              go_on = false;
            }
          else
            {
              tag_def.push_back(these_tags.substr(0,n));
              these_tags.erase(0,n+1);
            }
        }

      // Unfortunately, MTL conatains a bug that leads to all elements of
      // the outer Array of an Array<Array>> pointing to the same data
      // after creation. So we need to fix this explicitly:
      wfs_tag_groups[i] = Array<OneTag>();

      for ( Index s=0; s<tag_def.nelem(); ++s )
        {
          // Remove leading whitespace, if there is any:
          while ( ' '  == tag_def[s][0] ||
                  '\t' == tag_def[s][0]    )    tag_def[s].erase(0,1);

          OneTag this_tag(tag_def[s]);

          // Safety check: For s>0 check that the tags belong 
          // to the same species.
          if (s>0)
            if ( wfs_tag_groups[i][0].Species() != this_tag.Species() )
              throw runtime_error(
                       "Tags in a tag group must belong to the same species.");

          wfs_tag_groups[i].push_back(this_tag);
        }
    }

  // Print list of tag groups to the most verbose output stream:
  out3 << "  Defined weighting function tag groups:";
  for ( Index i=0; i<wfs_tag_groups.nelem(); ++i )
    {
      out3 << "\n  " << i << ":";
      for ( Index s=0; s<wfs_tag_groups[i].nelem(); ++s )
        {
          out3 << " " << wfs_tag_groups[i][s].Name();
        }
    }
  out3 << '\n';
}

void wfss_tgsDefine(// WS Output:
                   TagGroups& wfss_tgs,
                   // Control Parameters:
                   const ArrayOfString& tags)
{
  wfss_tgs = TagGroups(tags.nelem());

  // Each element of the array of Strings tags defines one tag
  // group. Let's work through them one by one.
  for ( Index i=0; i<tags.nelem(); ++i )
    {
      // There can be a comma separated list of tag definitions, so we
      // need to break the String apart at the commas.
      ArrayOfString tag_def;

      bool go_on = true;
      String these_tags = tags[i];
      while (go_on)
        {
          Index n = these_tags.find(',');
          if ( n == these_tags.npos ) // npos indicates `not found'
            {
              // There are no more commas.
              tag_def.push_back(these_tags);
              go_on = false;
            }
          else
            {
              tag_def.push_back(these_tags.substr(0,n));
              these_tags.erase(0,n+1);
            }
        }

      // Unfortunately, MTL conatains a bug that leads to all elements of
      // the outer Array of an Array<Array>> pointing to the same data
      // after creation. So we need to fix this explicitly:
      wfss_tgs[i] = Array<OneTag>();

      for ( Index s=0; s<tag_def.nelem(); ++s )
        {
          // Remove leading whitespace, if there is any:
          while ( ' '  == tag_def[s][0] ||
                  '\t' == tag_def[s][0]    )    tag_def[s].erase(0,1);

          OneTag this_tag(tag_def[s]);

          // Safety check: For s>0 check that the tags belong 
          // to the same species.
          if (s>0)
            if ( wfss_tgs[i][0].Species() != this_tag.Species() )
              throw runtime_error(
                       "Tags in a tag group must belong to the same species.");

          wfss_tgs[i].push_back(this_tag);
        }
    }

  // Print list of tag groups to the most verbose output stream:
  out3 << "  Defined weighting function tag groups:";
  for ( Index i=0; i<wfss_tgs.nelem(); ++i )
    {
      out3 << "\n  " << i << ":";
      for ( Index s=0; s<wfss_tgs[i].nelem(); ++s )
        {
          out3 << " " << wfss_tgs[i][s].Name();
        }
    }
  out3 << '\n';
}

void absloswfsCalc (
                    ArrayOfMatrix&   absloswfs,
              const Index&             emission,
              const Los&             los,   
              const ArrayOfMatrix&   source,
              const ArrayOfMatrix&   trans,
              const Vector&          y,
              const Vector&          y_space,
              const Vector&          f_mono,
              const Vector&          e_ground,
              const Numeric&         t_ground )
{
  check_if_bool( emission, "emission" );                                      

  if ( emission == 0 )
    absloswfs_tau ( absloswfs, los, f_mono );
  else
    absloswfs_rte ( absloswfs, los, source, trans, y, y_space, f_mono,     
                                                          e_ground, t_ground );

  // The equations used can produce NaNs and Infs when the optical thickness
  // is very high. This corresponds to extremly low values for the ABS LOS WFs
  // so we set NaNs and Infs to zero.
  Index     irow, icol, ncol;
  Numeric   w;
  for( Index im=0; im<absloswfs.nelem(); im++ )
  {
    for( irow=0; irow<absloswfs[im].nrows(); irow++ )
    {
      ncol = absloswfs[im].ncols();
      for( icol=0; icol<ncol; icol++ )
      {
        w = absloswfs[im](irow,icol);
#ifdef HPUX
        if ( !isfinite(w) )
#else
        if ( !finite(w) )
#endif
          absloswfs[im](irow,icol) = 0.0;
      }
    }
  }
}



void kSpecies (
                    Matrix&          k,
                    ArrayOfString&   k_names,
                    Matrix&          k_aux,
              const Los&             los,           
              const ArrayOfMatrix&   absloswfs,
              const Vector&          p_abs,
              const Vector&          t_abs,             
              const TagGroups&       wfs_tgs,
              const ArrayOfMatrix&   abs_per_tg,
              const Matrix&          vmrs,
              const Vector&          k_grid,
              const String&          unit )
{
  // Check of input is performed in k_species

  const Index   ntg = wfs_tgs.nelem();     // number of retrieval tags
  ArrayOfIndex  tg_nr(ntg);

  for ( Index i=0; i<ntg; i++ )
    tg_nr[i] = i;

  k_species( k, k_names, k_aux, los, absloswfs, p_abs, t_abs, 
             wfs_tgs, abs_per_tg, vmrs, k_grid, tg_nr, unit );
}



void kSpeciesSingle (
                    Matrix&          k,
                    ArrayOfString&   k_names,
                    Matrix&          k_aux,
              const Los&             los,           
              const ArrayOfMatrix&   absloswfs,
              const Vector&          p_abs,
              const Vector&          t_abs,             
              const TagGroups&       wfs_tgs,
              const ArrayOfMatrix&   abs_per_tg,
              const Matrix&          vmrs,
              const Vector&          k_grid,
              const String&          tg,
              const String&          unit )
{
  // Check of input is performed in k_species

  ArrayOfString  tg_name(1);
  tg_name[0] = tg;

  ArrayOfIndex   tg_nr; 
  get_tagindex_for_Strings( tg_nr, wfs_tgs, tg_name );
  
  k_species( k, k_names, k_aux, los, absloswfs, p_abs, t_abs, 
                             wfs_tgs, abs_per_tg, vmrs, k_grid, tg_nr, unit );
}



void kContAbs (
                    Matrix&          k,
                    ArrayOfString&   k_names,
                    Matrix&          k_aux,
              const Los&             los,           
              const ArrayOfMatrix&   absloswfs,
              const Vector&          f_mono,
              const Vector&          k_grid,
              const Index&           order,
              const Numeric&         f_low,
              const Numeric&         f_high,
              const String&          l_unit )
{
  // Input is checked in k_contabs  

  Numeric f1=f_low, f2=f_high;

  if ( f1 < 0 )
    f1 = first( f_mono );
  if ( f2 < 0 )
    f2 = last( f_mono );

  k_contabs( k, k_names, k_aux, los, absloswfs, f_mono, k_grid, order, f1, f2,
                                                                      l_unit );
}



void kTemp (
                Matrix&                      k,
                ArrayOfString&               k_names,
                Matrix&                      k_aux,
          const TagGroups&                   tgs,
          const Vector&                      f_mono,
          const Vector&                      p_abs,
          const Vector&                      t_abs,
          const Vector&                      n2_abs,
          const Vector&                      h2o_abs,
          const Matrix&                      vmrs,
          const Matrix&                      abs0,
          const ArrayOfArrayOfLineRecord&    lines_per_tg,
          const ArrayOfLineshapeSpec&        lineshape,
          const Vector&                      e_ground,
          const Index&                       emission,
          const Vector&                      k_grid,
          const ArrayOfString&               cont_description_names,
          const ArrayOfVector&               cont_description_parameters,
          const ArrayOfString&               cont_description_models,
          const Los&                         los,           
          const ArrayOfMatrix&               absloswfs,
          const ArrayOfMatrix&               trans,
          const Numeric&                     z_plat,
          const Vector&                      za,
          const Numeric&                     l_step,
          const Vector&                      z_abs,
          const Index&                       refr,
          const Index&                       refr_lfac,
          const Vector&                      refr_index,
          const String&                      refr_model,
          const Numeric&                     z_ground,
          const Numeric&                     t_ground,
          const Vector&                      y_space,
          const Numeric&                     r_geoid,
          const Vector&                      hse,
          // Keywords
          const Index&                       kw_hse,
          const Index&                       kw_fast )
{
  check_if_bool( kw_hse, "hse keyword" );
  check_if_bool( kw_fast, "fast keyword" );      
  check_if_bool( emission, "emission" );
  check_if_bool( refr, "refr" );
  check_lengths( p_abs, "p_abs", t_abs, "t_abs" );  
  check_lengths( p_abs, "p_abs", z_abs, "z_abs" );  
  check_lengths( p_abs, "p_abs", h2o_abs, "h2o_abs" );  
  check_lengths( p_abs, "p_abs", n2_abs, "n2_abs" );  
  check_length_ncol( p_abs, "p_abs", abs0, "abs" );
  check_length_ncol( p_abs, "p_abs", vmrs, "vmrs" );
  check_length_nrow( f_mono, "f_mono", abs0, "abs" );
  if ( los.p.nelem() != za.nelem() )
    throw runtime_error(
               "The number of zenith angles of *za* and *los* are different.");
  //
  if( refr ) 
    check_lengths( p_abs, "p_abs", refr_index, "refr_index" );  
  //
  if ( !kw_hse )
  {
    if ( los.p.nelem() != trans.nelem() )
      throw runtime_error(
            "The number of zenith angles of *los* and *trans* are different.");

    if ( los.p.nelem() != absloswfs.nelem() )
      throw runtime_error(
      "The number of zenith angles is not the same in *los* and *absloswfs*.");
  }
  else
  { 
    if ( hse[0]==0 )
      throw runtime_error( "Hydrostatic eq. must be considered generally" 
                           "when calculating WFs with hydrostatic eq.");
  }
  //
  if ( any_ground(los.ground) )  
  {
    if ( t_ground <= 0 )
      throw runtime_error(
          "There are intersections with the ground, but the ground\n"
          "temperature is set to be <=0 (are dummy values used?).");
    if ( e_ground.nelem() != f_mono.nelem() )
      throw runtime_error(
          "There are intersections with the ground, but the frequency and\n"
          "ground emission vectors have different lengths (are dummy values\n"
          "used?).");
  }
  if ( emission ) 
    check_lengths( f_mono, "f_mono", y_space, "y_space" );

  //
  // Three options:
  // 1. no hydrostatic eq., use analytical expressions
  // 2. hydrostatic eq., fast version
  // 3. hydrostatic eq., accurate version
  //


  // No hydrostatic eq., use analytical expressions
  //---------------------------------------------------------------------------
  if ( !kw_hse )
  {
    if ( !emission )
      throw runtime_error(
          "Analytical expressions for temperature and no emission have not\n"
          "yet been implemented. Sorry!");

    k_temp_nohydro( k, k_names, k_aux, tgs, los, absloswfs, f_mono, 
     p_abs, t_abs, n2_abs, h2o_abs, vmrs, lines_per_tg, lineshape, abs0, trans,
     e_ground, k_grid, cont_description_names, cont_description_parameters,
     cont_description_models);
  }


  // Hydrostatic eq., fast version
  //---------------------------------------------------------------------------
  else if ( kw_fast )
  {
    // Main sizes
    const Index  nza  = za.nelem();          // number of zenith angles  
    const Index  nv   = f_mono.nelem();      // number of frequencies
    const Index  np   = k_grid.nelem();      // number of retrieval altitudes
    const Index  nabs = p_abs.nelem();       // number of absorption altitudes
  
    // Vectors for the reference state
    Vector z0(nabs), y0, t0(np);
  
    // Local copy of hse
    Vector hse_local( hse );
  
    // Calculate reference z_abs with a high number of iterations
    hse_local[4] = 5;
    z0 = z_abs;
    hseCalc( z0, p_abs, t_abs, h2o_abs, r_geoid,  hse_local );
    hse_local[4] = hse[4];
  
    // Calculate absorption for + 1K
    Matrix         abs1k, abs(nv,nabs);
    ArrayOfMatrix  abs_dummy;
    //
    {
      Vector  t(t_abs);
      t += 1.;
      //
      out1 << "  Calculating absorption for t_abs + 1K \n";
      out2 << "  ----- Messages from absCalc: --------\n";
      absCalc( abs1k, abs_dummy, tgs, f_mono, p_abs, t, n2_abs, h2o_abs, vmrs, 
                 lines_per_tg, lineshape, cont_description_names, 
                         cont_description_models, cont_description_parameters);
    }
    // Calculate reference spectrum
    out1 << "  Calculating reference spectrum\n";
    out2 << "  ----- Messages from losCalc: --------\n";
    Los    los;
    Vector z_tan;
    losCalc( los, z_tan, z_plat, za, l_step, p_abs, z0, refr, refr_lfac, 
                                               refr_index, z_ground, r_geoid );
    out2 << "  -------------------------------------\n";
    out2 << "  ----- Messages from sourceCalc: -----\n";
    ArrayOfMatrix source, trans;
    sourceCalc( source, emission, los, p_abs, t_abs, f_mono );
    out2 << "  -------------------------------------\n";
    out2 << "  ----- Messages from transCalc: ------\n";
    transCalc( trans, los, p_abs, abs0 );
    out2 << "  -------------------------------------\n";
    out2 << "  ----- Messages from yRte: -----------\n";
    yCalc( y0, emission, los, f_mono, y_space, source, trans, 
                                                          e_ground, t_ground );
    out2 << "  -------------------------------------\n";
  
    // Allocate K and fill aux. variables
    k.resize(nza*nv,np);
    k_names.resize(1);
    k_names[0] = "Temperature: with hydrostatic eq.";
    k_aux.resize(np,2);
    interpp( t0, p_abs, t_abs, k_grid ); 
    for ( Index ip=0; ip<np; ip++ )
    {
       k_aux(ip,0) = k_grid[ip];
       k_aux(ip,1) = t0[ip];
    }
  
    // Determine conversion between grids        
    Matrix is;
    Vector lpabs, lgrid;
    p2grid( lpabs, p_abs );
    p2grid( lgrid, k_grid );
    grid2grid_index( is, lpabs, lgrid );
  
    // Loop retrieval altitudes and calculate new spectra
    //
    Vector y, t(nabs), w, refr4t, z4t(nabs);
    Index  i1, iw, iv;
    //
    for ( Index ip=0; ip<np; ip++ )
    {
      out1 << "  Doing altitude " << ip+1 << "/" << np << "\n";   
  
      // Create absorption matrix and temperature profile corresponding to 
      // temperature disturbance
      grid2grid_weights( w, lpabs, Index(is(ip,0)), Index(is(ip,1)), 
                                                                   lgrid, ip );
      i1 = Index( is(ip,0) );    // first p_abs point to consider

      abs = abs0;       
      t   = t_abs;      
      for ( iw=i1; iw<=Index(is(ip,1)); iw++ )
      {
        t[iw] += w[iw-i1];
        for ( iv=0; iv<nv; iv++ )
          abs(iv,iw) = (1-w[iw-i1])*abs0(iv,iw) + w[iw-i1]*abs1k(iv,iw);
      }
  
      out2 << "  ----- Doing new refractive index ----\n";
      refrCalc ( refr4t, p_abs, t, h2o_abs, refr, refr_model );
      out2 << "  ----- Doing new hydrostatic eq. -----\n";
      z4t = z0;
      hseCalc( z4t, p_abs, t, h2o_abs, r_geoid,  hse );
      out2 << "  ----- Messages from losCalc: --------\n";
      losCalc( los, z_tan, z_plat, za, l_step, p_abs, z4t, refr, refr_lfac, 
                                               refr_index, z_ground, r_geoid );
      out2 << "  -------------------------------------\n";
      out2 << "  ----- Messages from sourceCalc: -----\n";
      ArrayOfMatrix source, trans;
      sourceCalc( source, emission, los, p_abs, t, f_mono );
      out2 << "  -------------------------------------\n";
      out2 << "  ----- Messages from transCalc: ------\n";
      transCalc( trans, los, p_abs, abs );
      out2 << "  -------------------------------------\n";
      out2 << "  ----- Messages from yRte: -----------\n";
      yCalc( y, emission, los, f_mono, y_space, source, trans, 
                                                          e_ground, t_ground );
      out2 << "  -------------------------------------\n";
  
      // Fill K
      for ( iv=0; iv<nza*nv; iv++ )
        k(iv,ip) = y[iv] - y0[iv];
    }
  }


  // Hydrostatic eq., accurate version
  //---------------------------------------------------------------------------
  else
  {
    // Main sizes
    const Index  nza  = za.nelem();          // number of zenith angles  
    const Index  nv   = f_mono.nelem();      // number of frequencies
    const Index  np   = k_grid.nelem();      // number of retrieval altitudes
    const Index  nabs = p_abs.nelem();       // number of absorption altitudes
  
    // Vectors for the reference state
    Vector z0(nabs), y0, t0(np);
  
    // Local copy of hse
    Vector hse_local( hse );
  
    // Calculate reference z_abs with a high number of iterations
    hse_local[4] = 5;
    z0 = z_abs; 
    hseCalc( z0, p_abs, t_abs, h2o_abs, r_geoid,  hse_local );
    hse_local[4] = hse[4];
  
    // Calculate reference spectrum
    out1 << "  Calculating reference spectrum\n";
    out2 << "  ----- Messages from losCalc: --------\n";
    Los    los;
    Vector z_tan;
    losCalc( los, z_tan, z_plat, za, l_step, p_abs, z0, refr, refr_lfac, 
                                               refr_index, z_ground, r_geoid );
    out2 << "  -------------------------------------\n";
    out2 << "  ----- Messages from sourceCalc: -----\n";
    ArrayOfMatrix source, trans;
    sourceCalc( source, emission, los, p_abs, t_abs, f_mono );
    out2 << "  -------------------------------------\n";
    out2 << "  ----- Messages from transCalc: ------\n";
    transCalc( trans, los, p_abs, abs0 );
    out2 << "  -------------------------------------\n";
    out2 << "  ----- Messages from yRte: -----------\n";
    yCalc( y0, emission, los, f_mono, y_space, source, trans, 
                                                          e_ground, t_ground );
    out2 << "  -------------------------------------\n";
  
    // Allocate K and fill aux. variables
    k.resize(nza*nv,np);
    k_names.resize(1);
    k_names[0] = "Temperature: with hydrostatic eq.";
    k_aux.resize(np,2);
    interpp( t0, p_abs, t_abs, k_grid ); 
    for ( Index ip=0; ip<np; ip++ )
    {
       k_aux(ip,0) = k_grid[ip];
       k_aux(ip,1) = t0[ip];
    }
  
    // Determine conversion between grids        
    Matrix is;
    Vector lpabs, lgrid;
    p2grid( lpabs, p_abs );
    p2grid( lgrid, k_grid );
    grid2grid_index( is, lpabs, lgrid );
  
    // Loop retrieval altitudes and calculate new spectra
    //
    Matrix         abs;
    ArrayOfMatrix  abs_dummy;
    Vector y, t(nabs), w, refr4t, z4t(nabs);
    Index  i1, iw, iv;
    //
    for ( Index ip=0; ip<np; ip++ )
    {
      out1 << "  Doing altitude " << ip+1 << "/" << np << "\n";   
  
      // Create disturbed temperature profile
      grid2grid_weights( w, lpabs, Index(is(ip,0)), Index(is(ip,1)), 
                                                                   lgrid, ip );
      i1 = Index( is(ip,0) );    // first p_abs point to consider

      t = t_abs;        
      for ( iw=i1; iw<=Index(is(ip,1)); iw++ )
        t[iw] += w[iw-i1];
  
      out2 << "  ----- Messages from absCalc: --------\n";
      absCalc( abs, abs_dummy, tgs, f_mono, p_abs, t, n2_abs, h2o_abs, vmrs, 
                    lines_per_tg, lineshape, cont_description_names, 
                         cont_description_models, cont_description_parameters);
      out2 << "  ----- Doing new refractive index ----\n";
      refrCalc ( refr4t, p_abs, t, h2o_abs, refr, refr_model );
      out2 << "  ----- Doing new hydrostatic eq. -----\n";
      z4t = z0;
      hseCalc( z4t, p_abs, t, h2o_abs, r_geoid,  hse );
      out2 << "  ----- Messages from losCalc: --------\n";
      losCalc( los, z_tan, z_plat, za, l_step, p_abs, z4t, refr, refr_lfac, 
                                                   refr4t, z_ground, r_geoid );
      out2 << "  -------------------------------------\n";
      out2 << "  ----- Messages from sourceCalc: -----\n";
      ArrayOfMatrix source, trans;
      sourceCalc( source, emission, los, p_abs, t, f_mono );
      out2 << "  -------------------------------------\n";
      out2 << "  ----- Messages from transCalc: ------\n";
      transCalc( trans, los, p_abs, abs );
      out2 << "  -------------------------------------\n";
      out2 << "  ----- Messages from yRte: -----------\n";
      yCalc( y, emission, los, f_mono, y_space, source, trans, 
                                                          e_ground, t_ground );
      out2 << "  -------------------------------------\n";
  
      // Fill K
      for ( iv=0; iv<nza*nv; iv++ )
        k(iv,ip) = y[iv] - y0[iv];
    }
  }
}



void kFrequencyOffSet (
                Matrix&                      k,
                ArrayOfString&               k_names,
                Matrix&                      k_aux,
          const TagGroups&                   tgs,
          const Vector&                      f_mono,
          const Vector&                      p_abs,
          const Vector&                      t_abs,
          const Vector&                      n2_abs,
          const Vector&                      h2o_abs,
          const Matrix&                      vmrs,
          const ArrayOfArrayOfLineRecord&    lines_per_tg,
          const ArrayOfLineshapeSpec&        lineshape,
          const Vector&                      e_ground,
          const Index&                       emission,
          const ArrayOfString&               cont_description_names,
          const ArrayOfVector&               cont_description_parameters,
          const ArrayOfString&               cont_description_models,
          const Los&                         los,           
          const Numeric&                     t_ground,
          const Vector&                      y_space,
          const Vector&                      y,
          // Keywords
          const Numeric&                     delta,
          const String&                      f_unit )
{
  check_if_bool( emission, "emission" );
  check_lengths( p_abs, "p_abs", t_abs, "t_abs" );  
  check_lengths( p_abs, "p_abs", h2o_abs, "h2o_abs" );  
  check_lengths( p_abs, "p_abs", n2_abs, "n2_abs" );  
  check_length_ncol( p_abs, "p_abs", vmrs, "vmrs" );
  //
  if ( any_ground(los.ground) )  
  {
    if ( t_ground <= 0 )
      throw runtime_error(
          "There are intersections with the ground, but the ground\n"
          "temperature is set to be <=0 (are dummy values used?).");
    if ( e_ground.nelem() != f_mono.nelem() )
      throw runtime_error(
          "There are intersections with the ground, but the frequency and\n"
          "ground emission vectors have different lengths (are dummy values\n"
          "used?).");
  }
  if ( emission ) 
    check_lengths( f_mono, "f_mono", y_space, "y_space" );


  // Check frequency unit
  //
  Numeric scfac;
  //
  if ( f_unit == "Hz" )
    scfac = 1.0;
  else if ( f_unit == "kHz" )
    scfac = 1000;
  else if ( f_unit == "MHz" )
    scfac = 1e6;
  else
    throw runtime_error("Allowed frequency units are \"Hz\", \"kHz\" and "
                                                                   "\"MHz\".");


  // Create the disturbed f_mono
  //
  Index    nv     = f_mono.nelem();      // number of frequencies
  Numeric  fdelta = scfac * delta;
  //
  Vector f_dist(nv);
  //
  for( Index ii=0; ii<nv; ii++ )
    { f_dist[ii] = f_mono[ii] + fdelta; }


  // Calculate new spectra 
  Matrix         abs;
  ArrayOfMatrix  abs_dummy;
  out2 << "  ----- Messages from absCalc: --------\n";
  absCalc( abs, abs_dummy, tgs, f_dist, p_abs, t_abs, n2_abs, h2o_abs, vmrs, 
                    lines_per_tg, lineshape, cont_description_names, 
                         cont_description_models, cont_description_parameters);
  ArrayOfMatrix source, trans;
  out2 << "  ----- Messages from sourceCalc: -----\n";
  sourceCalc( source, emission, los, p_abs, t_abs, f_mono );
  out2 << "  -------------------------------------\n";
  out2 << "  ----- Messages from transCalc: ------\n";
  transCalc( trans, los, p_abs, abs );
  out2 << "  -------------------------------------\n";
  Vector y_new;
  out2 << "  ----- Messages from yRte: -----------\n";
  yCalc( y_new, emission, los, f_mono, y_space, source, trans, 
                                                         e_ground, t_ground );
  out2 << "  -------------------------------------\n";

  // Make k one-column matrix of the right size:
  nv = y.nelem();
  k.resize( nv, 1 );

  // k = (y_new - y) / delta
  for( Index ii=0; ii<nv; ii++ )
    { k(ii,0) = ( y_new[ii] - y[ii] ) / delta; }

  k_names.resize( 1 );
  k_names[0] = "Frequency: off-set";
  k_aux.resize( 1, 2 );
  k_aux = 0.0;                  // Matpack can set all elements like this.
}



void kPointingOffSet(
                    Matrix&          k,
                    ArrayOfString&   k_names,
                    Matrix&          k_aux,
              const Numeric&         z_plat,
              const Vector&          za_pencil,
              const Numeric&         l_step,
              const Vector&          p_abs,
              const Vector&          z_abs,
              const Vector&          t_abs,
              const Vector&          f_mono,
              const Index&           refr,
              const Index&           refr_lfac,
              const Vector&          refr_index,
              const Numeric&         z_ground,
              const Numeric&         r_geoid,
              const Matrix&          abs,
              const Index&           emission,
              const Vector&          y_space,
              const Vector&          e_ground,
              const Numeric&         t_ground,
              const Vector&          y,
              const Numeric&         delta )
{
  check_if_bool( emission, "emission" );
  check_if_bool( refr, "refr" );
  check_lengths( p_abs, "p_abs", t_abs, "t_abs" );  
  check_lengths( p_abs, "p_abs", z_abs, "z_abs" );  
  check_length_ncol( p_abs, "p_abs", abs, "abs" );
  check_length_nrow( f_mono, "f_mono", abs, "abs" );
  if ( emission ) 
    check_lengths( f_mono, "f_mono", y_space, "y_space" );
  if ( refr )
    check_lengths( p_abs, "p_abs", refr_index, "refr_index" );  


  // Create new zenith angle grid
  //  const Index  nza = za_pencil.nelem();
  Vector za_new( za_pencil );   // Matpack can initialize a
                                // new Vector from another
                                // Vector.
  za_new += delta;              // With Matpack you can add delta to all
                                // Vector elements like this.

  out2 << "  ----- Messages from losCalc: --------\n";
  Los    los;
  Vector z_tan;
  losCalc( los, z_tan, z_plat, za_new, l_step, p_abs, z_abs, refr, refr_lfac, 
                                              refr_index, z_ground, r_geoid );
  out2 << "  -------------------------------------\n";

  ArrayOfMatrix source, trans;
  out2 << "  ----- Messages from sourceCalc: -----\n";
  sourceCalc( source, emission, los, p_abs, t_abs, f_mono );
  out2 << "  -------------------------------------\n";
  out2 << "  ----- Messages from transCalc: ------\n";
  transCalc( trans, los, p_abs, abs );
  out2 << "  -------------------------------------\n";
  Vector y_new;
  out2 << "  ----- Messages from yRte: -----------\n";
  yCalc( y_new, emission, los, f_mono, y_space, source, trans, 
                                                         e_ground, t_ground );
  out2 << "  -------------------------------------\n";

  // Make k one-column matrix of the right size:
  const Index   nv = y.nelem();
  k.resize( nv, 1 );

  // k = (y_new - y) / delta
  // k = (y_new - y) / delta
  for( Index ii=0; ii<nv; ii++ )
    { k(ii,0) = ( y_new[ii] - y[ii] ) / delta; }
  //k = y_new;                  
  //k -= y;
  //k /= delta;

  k_names.resize( 1 );
  k_names[0] = "Pointing: off-set";
  k_aux.resize( 1, 2 );
  k_aux = 0.0;
}



void kEground(
                    Matrix&          k,
                    ArrayOfString&   k_names,
                    Matrix&          k_aux,
              const Vector&          za_pencil,
              const Vector&          f_mono,
              const Index&           emission,
              const Vector&          y_space,
              const Vector&          e_ground,
              const Numeric&         t_ground,
              const Los&             los,           
              const ArrayOfMatrix&   source,
              const ArrayOfMatrix&   trans,
              const Index&           single_e )
{
  check_if_bool( emission, "emission" );                                      
  check_if_bool( emission, "single_e" );                                      
  check_lengths( f_mono, "f_mono", e_ground, "e_ground" );
  if ( los.p.nelem() != za_pencil.nelem() )
    throw runtime_error(
               "The number of zenith angles of *za* and *los* are different.");
  //
  if ( los.p.nelem() != trans.nelem() )
    throw runtime_error(
            "The number of zenith angles of *los* and *trans* are different.");
  //
  if ( emission )
  {
    if ( los.p.nelem() != source.nelem() )
      throw runtime_error(
           "The number of zenith angles of *los* and *source* are different.");
    check_length_nrow( f_mono, "f_mono", source[0], 
                                            "the source matrices (in source)");
  }
  //
  if ( any_ground(los.ground) )  
  {
    if ( t_ground <= 0 )
      throw runtime_error(
          "There are intersections with the ground, but the ground\n"
          "temperature is set to be <=0 (are dummy values used?).");
    if ( e_ground.nelem() != f_mono.nelem() )
      throw runtime_error(
          "There are intersections with the ground, but the frequency and\n"
          "ground emission vectors have different lengths (are dummy values\n"
          "used?).");
  }
  //
  if ( single_e ) 
  {
    for ( INDEX iv=1; iv<f_mono.nelem(); iv++ )
    {
      if ( e_ground[iv] != e_ground[0] )
        throw runtime_error(
          "A single ground emission value is assumed, but all values of\n"
          "*e_ground* are not the same.");
    }    
  }


  // Main sizes
  const Index  nza  = za_pencil.nelem();   // number of zenith angles  
  const Index  nv   = f_mono.nelem();      // number of frequencies

  // Loop variables
  Index iv, iza;

  // Make calculations assuming a single e
  Vector ksingle( nza*nv );

  // Emission
  if ( emission )
  {
    // Set all values to zero (k is zero if no ground intersection)
    ksingle = 0.0;

    // Local vectors
    Vector bground(nv), y_in(nv), tr(nv);

    for ( iza=0; iza<nza; iza++ )
    {
      if ( los.ground[iza] )
      {
        // Calculate blackbody radiation of the ground
        planck( bground, f_mono, t_ground );

        // Calculate the intensity reflected by the ground
        rte( y_in, los.start[iza], los.ground[iza]-1, trans[iza], source[iza], 
                                               y_space, 0, e_ground, bground );

        // Calculate transmission from the ground to the sensor
        tr = 1.0;
        bl_iterate( tr, los.ground[iza]-1, los.stop[iza]-1, trans[iza], nv );

        for ( iv=0; iv<nv; iv++ )
          ksingle[iza*nv+iv] = ( bground[iv] - y_in[iv] ) * tr[iv];      
      }
    }
  }

  // Optical thicknesses
  else
  {
    Numeric a;
    for ( iv=0; iv<nv; iv++ )
    {
      a = 1 / ( 1 - e_ground[iv] );
      for ( iza=0; iza<nza; iza++ )
        ksingle[iza*nv+iv] = a;
    }
  }

  // A single value for the ground emission, k is a column vector
  if ( single_e )
  {
    k_names.resize( 1 );
    k_names[0] = "Ground emission (single)";
    //
    k.resize( nza*nv, 1 );
    k = ksingle;
    //
    k_aux.resize( 1, 2 );
    k_aux(0,0) = 0;
    k_aux(0,1) = e_ground[0];    
  }

  // An e for each monochromatic frequency, k has nv columns
  else
  {
    k_names.resize( 1 );
    k_names[0] = "Ground emission";
    //
    k.resize( nza*nv, nv );
    k = 0.0;
    k_aux.resize( nv, 2 );
    //
    for ( iv=0; iv<nv; iv++ )
    {
      for ( iza=0; iza<nza; iza++ )
        k(iza*nv+iv,iv) = ksingle[iza*nv+iv];
      k_aux(iv,0) = 0;
      k_aux(iv,1) = e_ground[iv];    
    }
  }
}



void kCalibration(
                    Matrix&          k,
                    ArrayOfString&   k_names,
                    Matrix&          k_aux,
              const Vector&          y,
              const Vector&          f_mono,
              const Vector&          y0,
              const String&          /* name */ )
{
  check_lengths( f_mono, "f_mono", y0, "y0" );

  const Index   ny = y.nelem();
  const Index   nf = f_mono.nelem();
  const Index   nza = ny/nf;

  // Make k one-column matrix of the right size:
  k.resize( ny, 1 );

  // k = y - y0
  Index j,i,i0;
  for ( j=0; j<nza; j++ )    
  {
    i0 = j*nf;
    for ( i=0; i<nf; i++ )
      k(i0+i,0) = y[i0+i] - y0[i];
  }

  k_names.resize( 1 );
  k_names[0] = "Calibration: scaling";
  k_aux.resize( 1, 2 );
  k_aux = 0.0;          
}



void kManual(
                    Matrix&          k,
                    ArrayOfString&   k_names,
                    Matrix&          k_aux,
              const Vector&          y0,      
              const Vector&          y,
              const String&          name,
              const Numeric&         delta,
              const Numeric&         grid,
              const Numeric&         apriori )
{
  check_lengths( y, "y", y0, "y0" );

  // Make k one-column matrix of the right size:
  k.resize( y.nelem(), 1 );

  // Compute (y-y0)
  k = y;
  k -= y0;
  k /= delta;

  k_names.resize(1);
  k_names[0] = name;
  k_aux.resize(1,2);
  k_aux(0,0) = grid;
  k_aux(0,1) = apriori;
}



void kxInit (
                    Matrix&          kx,
                    ArrayOfString&   kx_names,
                    ArrayOfIndex&    kx_lengths,
                    Matrix&          kx_aux )
{
  kx.resize(         0, 0 );
  kx_names.resize(   0    );
  kx_lengths.resize( 0    );
  kx_aux.resize(     0, 0 );
}



void kbInit (
                    Matrix&          kb,
                    ArrayOfString&   kb_names,
                    ArrayOfIndex&    kb_lengths,
                    Matrix&          kb_aux )
{
  kxInit( kb, kb_names, kb_lengths, kb_aux );
}



void kxAppend (
                    Matrix&          kx,
                    ArrayOfString&   kx_names,
                    ArrayOfIndex&    kx_lengths,
                    Matrix&          kx_aux,
              const Matrix&          k,
              const ArrayOfString&   k_names,
              const Matrix&          k_aux )
{
  k_append( kx, kx_names, kx_lengths, kx_aux, k, k_names, k_aux );
}



void kbAppend (
                    Matrix&          kb,
                    ArrayOfString&   kb_names,
                    ArrayOfIndex&    kb_lengths,
                    Matrix&          kb_aux,
              const Matrix&          k,
              const ArrayOfString&   k_names,
              const Matrix&          k_aux )
{
  k_append( kb, kb_names, kb_lengths, kb_aux, k, k_names, k_aux );
}



void kxAllocate (
                    Matrix&          kx,
                    ArrayOfString&   kx_names,
                    ArrayOfIndex&    kx_lengths,
                    Matrix&          kx_aux,
              const Vector&          y,
              const String&          /* y_name */,
              const Index&             ni,
              const Index&             nx )
{
  const Index  ny = y.nelem();

  out2 << "  Allocates memory for a weighting function matrix having:\n" 
       << "    " << ny << " frequency points\n"
       << "    " << nx << " state variables\n"
       << "    " << ni << " retrieval identities\n";

  kx.resize(         ny, nx );
  kx_names.resize(   ni     );
  kx_lengths.resize( ni     );
  kx_aux.resize(     nx, 2  );

  for ( Index i=0; i<ni; i++ )
    kx_lengths[i] = 0;
}



void kbAllocate (
                    Matrix&          kb,
                    ArrayOfString&   kb_names,
                    ArrayOfIndex&    kb_lengths,
                    Matrix&          kb_aux,
              const Vector&          y,
              const String&          y_name,
              const Index&             ni,
              const Index&             nx )
{
  kxAllocate( kb, kb_names, kb_lengths, kb_aux, y, y_name, ni, nx );
}



void kxPutInK (
                    Matrix&          kx,
                    ArrayOfString&   kx_names,
                    ArrayOfIndex&    kx_lengths,
                    Matrix&          kx_aux,
              const Matrix&          k,
              const ArrayOfString&   k_names,
              const Matrix&          k_aux )
{
  const Index  ny  = kx.nrows();
  const Index  nx  = kx.ncols();
  const Index  ni  = kx_names.nelem();
  const Index  nx2 = k.ncols();
  const Index  ni2 = k_names.nelem();

  if ( ny != k.nrows() )
    throw runtime_error("The two WF matrices have different number of rows.");

  // Determine old length of x and number of identities
  Index  ni1, nx1=0;
  for( ni1=0; ni1<ni && kx_lengths[ni1]>0; ni1++ )
    nx1 += kx_lengths[ni1];
  if ( ni1 == ni )
    throw runtime_error("All retrieval/error identity positions used.");

  // Check if new data fit
  if ( (ni1+ni2) > ni )
  {
    ostringstream os;
    os << "There is not room for so many retrieval/error identities.\n" 
       << (ni1+ni2)-ni << " positions are missing";
      throw runtime_error(os.str());
  }
  if ( (nx1+nx2) > nx )
  {
    ostringstream os;
    os << "The k-matrix does not fit in kx.\n" 
       << (nx1+nx2)-nx << " columns are missing";
      throw runtime_error(os.str());
  }

  out2 << "  Putting k in kx as\n" 
       << "    identity " << ni1 << " - " << ni1+ni2-1 << "\n"
       << "      column " << nx1 << " - " << nx1+nx2-1 << "\n";
 
  // Calculate the vector length for each identity in K
  Index l = nx2/ni2;

  // Move K to Kx
  kx(     Range(joker),   Range(nx1,nx2) ) = k;
  kx_aux( Range(nx1,nx2), Range(joker)   ) = k_aux;
  // For the Array types kx_names and kx_lenghts we cannot use Range.  
  for ( Index iri=0; iri<ni2; iri++ )
  {
    kx_names[ni1+iri]   = k_names[iri];
    kx_lengths[ni1+iri] = l;
  }   
}



void kbPutInK (
                    Matrix&          kb,
                    ArrayOfString&   kb_names,
                    ArrayOfIndex&    kb_lengths,
                    Matrix&          kb_aux,
              const Matrix&          k,
              const ArrayOfString&   k_names,
              const Matrix&          k_aux )
{
  kxPutInK( kb, kb_names, kb_lengths, kb_aux, k, k_names, k_aux );
}

---