rte.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2008
   Patrick Eriksson <Patrick.Eriksson@rss.chalmers.se>
   Stefan Buehler   <sbuehler@ltu.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
  ===========================================================================*/
 
#include <cmath>
#include <stdexcept>
#include "auto_md.h"
#include "check_input.h"
#include "logic.h"
#include "math_funcs.h"
#include "montecarlo.h"
#include "physics_funcs.h"
#include "ppath.h"
#include "rte.h"
#include "special_interp.h"
#include "lin_alg.h"
 
extern const Numeric BOLTZMAN_CONST;
extern const Numeric PLANCK_CONST;
extern const Numeric SPEED_OF_LIGHT;
extern const Numeric PI;
extern const Numeric DEG2RAD;
extern const Numeric RAD2DEG;
 
 
 
/*===========================================================================
  === The functions in alphabetical order
  ===========================================================================*/
 
 
void apply_y_unit( 
            MatrixView   iy, 
         const String&   y_unit, 
       ConstVectorView   f_grid,
   const ArrayOfIndex&   i_pol )
{
  // The code is largely identical between the two apply_y_unit functions.
  // If any change here, remember to update the other function.
 
  const Index nf = iy.nrows();
  const Index ns = iy.ncols();
 
  assert( f_grid.nelem() == nf );
  assert( i_pol.nelem() == ns );
 
  if( y_unit == "1" )
    {}
 
  else if( y_unit == "RJBT" )
    {
      for( Index iv=0; iv<nf; iv++ )
        {
          const Numeric scfac = invrayjean( 1, f_grid[iv] );
          for( Index is=0; is<ns; is++ )
            {
              if( i_pol[is] < 5 )           // Stokes components
                { iy(iv,is) *= scfac; }
              else                          // Measuement single pols
                { iy(iv,is) *= 2*scfac; }
            }
        }
    }
 
  else if( y_unit == "PlanckBT" )
    {
      for( Index iv=0; iv<nf; iv++ )
        {
          for( Index is=ns-1; is>=0; is-- ) // Order must here be reversed
            {
              if( i_pol[is] == 1 )
                { iy(iv,is) = invplanck( iy(iv,is), f_grid[iv] ); }
              else if( i_pol[is] < 5 )
                { 
                  assert( i_pol[0] == 1 );
                  iy(iv,is) = 
                    invplanck( 0.5*(iy(iv,0)+iy(iv,is)), f_grid[iv] ) -
                    invplanck( 0.5*(iy(iv,0)-iy(iv,is)), f_grid[iv] );
                }
              else
                { iy(iv,is) = invplanck( 2*iy(iv,is), f_grid[iv] ); }
            }
        }
    }
  
  else if ( y_unit == "W/(m^2 m sr)" )
    {
      for( Index iv=0; iv<nf; iv++ )
        {
          const Numeric scfac = f_grid[iv] * ( f_grid[iv] / SPEED_OF_LIGHT );
          for( Index is=0; is<ns; is++ )
            { iy(iv,is) *= scfac; }
        }
    }
  
  else if ( y_unit == "W/(m^2 m-1 sr)" )
    {
      iy *= SPEED_OF_LIGHT;
    }
 
  else
    {
      ostringstream os;
      os << "Unknown option: y_unit = \"" << y_unit << "\"\n" 
         << "Recognised choices are: \"1\", \"RJBT\", \"PlanckBT\""
         << "\"W/(m^2 m sr)\" and \"W/(m^2 m-1 sr)\""; 
      
      throw runtime_error( os.str() );      
    }
}
 
 
 
 
void apply_y_unit2( 
   Tensor3View           J,
   ConstMatrixView       iy, 
   const String&         y_unit, 
   ConstVectorView       f_grid,
   const ArrayOfIndex&   i_pol )
{
  // The code is largely identical between the two apply_y_unit functions.
  // If any change here, remember to update the other function.
 
  const Index nf = iy.nrows();
  const Index ns = iy.ncols();
  const Index np = J.npages();
 
  assert( J.nrows() == nf );
  assert( J.ncols() == ns );
  assert( f_grid.nelem() == nf );
  assert( i_pol.nelem() == ns );
 
  if( y_unit == "1" )
    {}
 
  else if( y_unit == "RJBT" )
    {
      for( Index iv=0; iv<nf; iv++ )
        {
          const Numeric scfac = invrayjean( 1, f_grid[iv] );
          for( Index is=0; is<ns; is++ )
            {
              if( i_pol[is] < 5 )           // Stokes componenets
                {
                  for( Index ip=0; ip<np; ip++ )
                    { J(ip,iv,is) *= scfac; }
                }
              else                          // Measuement single pols
                {
                  for( Index ip=0; ip<np; ip++ )
                    { J(ip,iv,is) *= 2*scfac; }
                }
            }
        }
    }
 
  else if( y_unit == "PlanckBT" )
    {
      for( Index iv=0; iv<f_grid.nelem(); iv++ )
        {
          for( Index is=ns-1; is>=0; is-- )
            {
              Numeric scfac = 1;
              if( i_pol[is] == 1 )
                { scfac = dinvplanckdI( iy(iv,is), f_grid[iv] ); }
              else if( i_pol[is] < 5 )
                {
                  assert( i_pol[0] == 1 );
                  scfac = 
                    dinvplanckdI( 0.5*(iy(iv,0)+iy(iv,is)), f_grid[iv] ) +
                    dinvplanckdI( 0.5*(iy(iv,0)-iy(iv,is)), f_grid[iv] );
                }
              else
                { scfac = dinvplanckdI( 2*iy(iv,is), f_grid[iv] ); }
              //
              for( Index ip=0; ip<np; ip++ )
                { J(ip,iv,is) *= scfac; }
            }
        }
    }
 
  else if ( y_unit == "W/(m^2 m sr)" )
    {
      for( Index iv=0; iv<nf; iv++ )
        {
          const Numeric scfac = f_grid[iv] * ( f_grid[iv] / SPEED_OF_LIGHT );
          for( Index ip=0; ip<np; ip++ )
            {
              for( Index is=0; is<ns; is++ )
                { J(ip,iv,is) *= scfac; }
            }
        }
    }
  
  else if ( y_unit == "W/(m^2 m-1 sr)" )
    {
      J *= SPEED_OF_LIGHT;
    }
  
  else
    {
      ostringstream os;
      os << "Unknown option: y_unit = \"" << y_unit << "\"\n" 
         << "Recognised choices are: \"1\", \"RJBT\", \"PlanckBT\""
         << "\"W/(m^2 m sr)\" and \"W/(m^2 m-1 sr)\""; 
      
      throw runtime_error( os.str() );      
    }
  
}
 
 
 
 
void ext2trans(//Output and Input:
              MatrixView trans_mat,
              //Input
              ConstMatrixView ext_mat_av,
              const Numeric& l_step )
{
  //Stokes dimension:
  Index stokes_dim = ext_mat_av.nrows();
 
  //Check inputs:
  assert( is_size(trans_mat, stokes_dim, stokes_dim) );
  assert( is_size(ext_mat_av, stokes_dim, stokes_dim) );
  assert( l_step >= 0 );
  assert( !is_singular( ext_mat_av ) );
 
  // Any changes here should also be implemented in rte_step_std.
 
  //--- Scalar case: ---------------------------------------------------------
  if( stokes_dim == 1 )
    {
      trans_mat(0,0) = exp(-ext_mat_av(0,0) * l_step);
    }
 
 
  //--- Vector case: ---------------------------------------------------------
 
  //- Unpolarised
  else if( is_diagonal(ext_mat_av) )
    {
      const Numeric tv = exp(-ext_mat_av(0,0) * l_step);
 
      trans_mat = 0;
 
      for( Index i=0; i<stokes_dim; i++ )
        {
          trans_mat(i,i)  = tv;
        }
    }
 
 
  //- General case
  else
    {
      Matrix ext_mat_ds = ext_mat_av;
      ext_mat_ds *= -l_step; 
 
      Index q = 10;  // index for the precision of the matrix exp function
 
      matrix_exp( trans_mat, ext_mat_ds, q );
    }
}
 
 
 
 
void get_ppath_atmvars( 
        Vector&      ppath_p, 
        Vector&      ppath_t, 
        Matrix&      ppath_vmr, 
        Vector&      ppath_wind_u, 
        Vector&      ppath_wind_v, 
        Vector&      ppath_wind_w, 
  const Ppath&       ppath,
  const Index&       atmosphere_dim,
  ConstVectorView    p_grid,
  ConstTensor3View   t_field,
  ConstTensor4View   vmr_field,
  ConstTensor3View   wind_u_field,
  ConstTensor3View   wind_v_field,
  ConstTensor3View   wind_w_field )
{
  const Index   np  = ppath.np;
 
  // Pressure:
  ppath_p.resize(np);
  Matrix itw_p(np,2);
  interpweights( itw_p, ppath.gp_p );      
  itw2p( ppath_p, p_grid, ppath.gp_p, itw_p );
  
  // Temperature:
  ppath_t.resize(np);
  Matrix   itw_field;
  interp_atmfield_gp2itw( itw_field, atmosphere_dim, 
                          ppath.gp_p, ppath.gp_lat, ppath.gp_lon );
  interp_atmfield_by_itw( ppath_t,  atmosphere_dim, t_field, 
                          ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
 
  // VMR fields:
  const Index ns = vmr_field.nbooks();
  ppath_vmr.resize(ns,np);
  for( Index is=0; is<ns; is++ )
    {
      interp_atmfield_by_itw( ppath_vmr(is, joker), atmosphere_dim,
                              vmr_field( is, joker, joker,  joker ), 
                              ppath.gp_p, ppath.gp_lat, ppath.gp_lon, 
                              itw_field );
    }
 
  // Winds:
  //
  ppath_wind_w.resize(np);
  if( wind_w_field.npages() > 0 ) 
    { 
      interp_atmfield_by_itw( ppath_wind_w,  atmosphere_dim, wind_w_field, 
                          ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
    }
  else
    { ppath_wind_w = 0; }
  //
  ppath_wind_v.resize(np);
  if( wind_v_field.npages() > 0 ) 
    { 
      interp_atmfield_by_itw( ppath_wind_v,  atmosphere_dim, wind_v_field, 
                          ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
    }
  else
    { ppath_wind_v = 0; }
  //
  ppath_wind_u.resize(np);
  if( atmosphere_dim > 2 )
    {
      if( wind_u_field.npages() > 0 ) 
        { 
          interp_atmfield_by_itw( ppath_wind_u,  atmosphere_dim, wind_u_field, 
                          ppath.gp_p, ppath.gp_lat, ppath.gp_lon, itw_field );
        }
      else
        { ppath_wind_u = 0; }
    }
  else
    { ppath_wind_u = 0; }
}
 
 
 
 
void get_ppath_pnd( 
        Matrix&         ppath_pnd, 
  const Ppath&          ppath,
  const Index&          atmosphere_dim,
  const ArrayOfIndex&   cloudbox_limits,
  ConstTensor4View      pnd_field )
{
  ppath_pnd.resize( pnd_field.nbooks(), ppath.np );
 
  Matrix   itw_field;
  ArrayOfGridPos   gpc_p, gpc_lat, gpc_lon;
  //
  interp_cloudfield_gp2itw( itw_field, gpc_p, gpc_lat, gpc_lon, 
                            ppath.gp_p, ppath.gp_lat, ppath.gp_lon, 
                            atmosphere_dim, cloudbox_limits );
  for( Index ip=0; ip<pnd_field.nbooks(); ip++ )
    {
      interp_atmfield_by_itw( ppath_pnd(ip,joker), atmosphere_dim,
                              pnd_field(ip,joker,joker,joker), 
                              gpc_p, gpc_lat, gpc_lon, itw_field );
    }
}
 
 
 
void get_ppath_rtvars( 
        Workspace&   ws,
        Tensor3&     ppath_abs_scalar, 
        Matrix&      ppath_tau,
        Vector&      total_tau,
        Matrix&      ppath_emission, 
  const Agenda&      abs_scalar_gas_agenda,
  const Agenda&      emission_agenda,
  const Ppath&       ppath,
  ConstVectorView    ppath_p, 
  ConstVectorView    ppath_t, 
  ConstMatrixView    ppath_vmr, 
  ConstVectorView    ppath_wind_u, 
  ConstVectorView    ppath_wind_v, 
  ConstVectorView    ppath_wind_w, 
  ConstVectorView    f_grid, 
  const Index&       atmosphere_dim,
  const Index&       emission_do )
{
  // Sizes
  const Index   np   = ppath.np;
  const Index   nabs = ppath_vmr.nrows();
  const Index   nf   = f_grid.nelem();
 
  // Init variables
  ppath_abs_scalar.resize( nf, nabs, np );
  ppath_tau.resize( nf, np-1 );
  total_tau.resize( nf );
  total_tau = 0;
  if( emission_do )
    { ppath_emission.resize( nf, np ); }
  else
    { ppath_emission.resize( 0, 0 ); }
 
  // Mean of extreme frequencies
  const Numeric f0 = (f_grid[0]+f_grid[nf-1])/2.0;
 
  for( Index ip=0; ip<np; ip++ )
    {
      // Doppler shift
      //
      Numeric rte_doppler = 0;
      //
      if( ppath_wind_v[ip]!=0 || ppath_wind_u[ip]!=0 || ppath_wind_w[ip]!=0 )
        {
          // za nd aa of winds and LOS
          const Numeric v = sqrt( pow(ppath_wind_u[ip],2) + 
                           pow(ppath_wind_v[ip],2) + pow(ppath_wind_w[ip],2) );
          Numeric aa_w = 0, aa_p = 0;
          if( atmosphere_dim < 3 )
            {
              if( ppath_wind_v[ip]<0 ) 
                { aa_w = PI; }  // Negative v-winds for 1 and 2D
              if( atmosphere_dim == 2  &&  ppath.los(ip,0) < 0 )
                { aa_p = PI; }  // Negative sensor za for 2D
            }
          else //3D 
            { 
              aa_w = atan2( ppath_wind_u[ip], ppath_wind_v[ip] ); 
              aa_p = DEG2RAD * ppath.los(ip,1);
            }
          const Numeric za_w = acos( ppath_wind_w[ip] / v );
          const Numeric za_p = DEG2RAD * fabs( ppath.los(ip,0) );
          //
          // Actual shift
          const Numeric costerm = cos(za_w) * cos(za_p) + 
                                  sin(za_w) * sin(za_p)*cos(aa_w-aa_p); 
 
          rte_doppler = -v * f0 * costerm / SPEED_OF_LIGHT;
        }
 
      // We must use temporary variables as the agenda input must be
      // free to be resized
      Vector   evector;
      Matrix   sgmatrix;
      //
      abs_scalar_gas_agendaExecute( ws, sgmatrix, -1, rte_doppler, ppath_p[ip], 
                                    ppath_t[ip], ppath_vmr(joker,ip), 
                                    abs_scalar_gas_agenda );
      ppath_abs_scalar(joker,joker,ip) = sgmatrix;
      //
      if( emission_do )
        {
          emission_agendaExecute( ws, evector, ppath_t[ip], emission_agenda );
          ppath_emission(joker,ip) = evector;
        }
 
      // Partial and total tau
      //
      if( ip > 0 )
        {
          for( Index iv=0; iv<nf; iv++ )
            { 
              ppath_tau(iv,ip-1) = 0.5 * ppath.l_step[ip-1] * (
                                         ppath_abs_scalar(iv,joker,ip-1).sum() +
                                         ppath_abs_scalar(iv,joker,ip).sum() );
              total_tau[iv] += ppath_tau(iv,ip-1);
            }
        }
    }
}
 
 
 
 
void get_ppath_cloudrtvars( 
        Workspace&                     ws,
        Tensor3&                       ppath_asp_abs_vec, 
        Tensor4&                       ppath_asp_ext_mat, 
        Tensor3&                       ppath_pnd_abs_vec, 
        Tensor4&                       ppath_pnd_ext_mat, 
        Tensor4&                       ppath_transmission,
        Tensor3&                       total_transmission,
        Matrix&                        ppath_emission, 
  Array<ArrayOfSingleScatteringData>&  scat_data,
  const Agenda&                        abs_scalar_gas_agenda,
  const Agenda&                        emission_agenda,
  const Agenda&                        opt_prop_gas_agenda,
  const Ppath&                         ppath,
  ConstVectorView                      ppath_p, 
  ConstVectorView                      ppath_t, 
  ConstMatrixView                      ppath_vmr, 
  ConstVectorView                      ppath_wind_u, 
  ConstVectorView                      ppath_wind_v, 
  ConstVectorView                      ppath_wind_w, 
  ConstMatrixView                      ppath_pnd, 
  const Index&                         use_mean_scat_data,
  const ArrayOfSingleScatteringData&   scat_data_raw,
  const Index&                         stokes_dim,
  ConstVectorView                      f_grid, 
  const Index&                         atmosphere_dim,
  const Index&                         emission_do,
  const Verbosity&                     verbosity)
{
  // Sizes
  const Index   np   = ppath.np;
  const Index   nf   = f_grid.nelem();
 
  // Init variables
  ppath_asp_abs_vec.resize( nf, stokes_dim, np );
  ppath_asp_ext_mat.resize( nf, stokes_dim, stokes_dim, np );
  ppath_pnd_abs_vec.resize( nf, stokes_dim, np );
  ppath_pnd_ext_mat.resize( nf, stokes_dim, stokes_dim, np );
  ppath_transmission.resize( nf, stokes_dim, stokes_dim, np-1 );
  total_transmission.resize( nf, stokes_dim, stokes_dim );
  //
  if( emission_do )
    { ppath_emission.resize( nf, np ); }
  else
    { ppath_emission.resize( 0, 0 ); }
 
  // Mean of extreme frequencies
  const Numeric f0 = (f_grid[0]+f_grid[nf-1])/2.0;
 
 
  // Particle single scattering properties (are independent of position)
  //
  if( use_mean_scat_data )
    {
      scat_data.resize( 1 );
      scat_data_monoCalc( scat_data[0], scat_data_raw, Vector(1,f0), 0, verbosity );
    }
  else
    {
      scat_data.resize( nf );
      for( Index iv=0; iv<nf; iv++ )
        { scat_data_monoCalc( scat_data[iv], scat_data_raw, f_grid, iv, verbosity ); }
    }
 
 
  // Loop ppath points
  //
  for( Index ip=0; ip<np; ip++ )
    {
      // Doppler shift 
      // 
      assert( ppath_wind_u[ip] == 0 );
      assert( ppath_wind_v[ip] == 0 );
      assert( ppath_wind_w[ip] == 0 );
      //
      const Numeric rte_doppler = 0;
 
      // Emission source term
      //
      if( emission_do )
        {
          Vector   evector;   // Agenda must be free to resize
          emission_agendaExecute( ws, evector, ppath_t[ip], emission_agenda );
          ppath_emission(joker,ip) = evector;
        }
 
      // Absorption species properties 
        {
          Matrix   asgmatrix;   // Agendas must be free to resize
          Matrix   abs_vec( nf, stokes_dim, 0 );
          Tensor3  ext_mat( nf, stokes_dim, stokes_dim, 0 );
          //
          abs_scalar_gas_agendaExecute( ws, asgmatrix, -1, rte_doppler, 
                                        ppath_p[ip], ppath_t[ip], 
                                        ppath_vmr(joker,ip), 
                                        abs_scalar_gas_agenda );
          opt_prop_gas_agendaExecute( ws, ext_mat, abs_vec, -1, asgmatrix, 
                                      opt_prop_gas_agenda );
          ppath_asp_ext_mat(joker,joker,joker,ip) = ext_mat;
          ppath_asp_abs_vec(joker,joker,ip)       = abs_vec;
        }
 
      // Particle properties 
        {
          // Direction of outgoing scattered radiation (which is reversed to
          // LOS). Note that rte_los2 is only used for extracting scattering
          // properties.
          Vector rte_los2;
          mirror_los( rte_los2, ppath.los(ip,joker), atmosphere_dim );
 
          // Extinction and absorption
          if( use_mean_scat_data )
            {
              Vector   abs_vec( stokes_dim, 0 );
              Matrix   ext_mat( stokes_dim, stokes_dim, 0 );
              opt_propCalc( ext_mat, abs_vec,
                            rte_los2[0], rte_los2[1], scat_data[0], 
                            stokes_dim, ppath_pnd(joker,ip), ppath_t[ip],
                            verbosity);
              for( Index iv=0; iv<nf; iv++ )
                { 
                  ppath_pnd_ext_mat(iv,joker,joker,ip) = ext_mat;
                  ppath_pnd_abs_vec(iv,joker,ip)       = abs_vec;
                }
            }
          else
            {
              for( Index iv=0; iv<nf; iv++ )
                { 
                  Vector   abs_vec( stokes_dim, 0 );
                  Matrix   ext_mat( stokes_dim, stokes_dim, 0 );
                  opt_propCalc( ext_mat, abs_vec,
                                rte_los2[0], rte_los2[1], scat_data[iv], 
                                stokes_dim, ppath_pnd(joker,ip), ppath_t[ip],
                                verbosity );
                  ppath_pnd_ext_mat(iv,joker,joker,ip) = ext_mat;
                  ppath_pnd_abs_vec(iv,joker,ip)       = abs_vec;
                }
            }
        }
       
      // Transmission
      //
      if( ip == 0 )   // Set total_transmission to identity matrices
        {
          for( Index iv=0; iv<nf; iv++ )
            { id_mat( total_transmission(iv,joker,joker) ); } 
        }
      else
        {
          for( Index iv=0; iv<nf; iv++ )
            { 
              // Average extinction matrix
              Matrix  ext_mat_av( stokes_dim, stokes_dim,0 );
              for( Index is1=0; is1<stokes_dim; is1++ )
                { 
                  for( Index is2=0; is2<stokes_dim; is2++ )
                    {
                      ext_mat_av(is1,is2) = 0.5 * (
                                          ppath_asp_ext_mat(iv,is1,is2,ip-1) +
                                          ppath_asp_ext_mat(iv,is1,is2,ip)   +
                                          ppath_pnd_ext_mat(iv,is1,is2,ip-1) +
                                          ppath_pnd_ext_mat(iv,is1,is2,ip) );
                    }
                }
              // Transmission
              ext2trans( ppath_transmission(iv,joker,joker,ip-1),
                         ext_mat_av, ppath.l_step[ip-1] );  
              const Matrix tmp = total_transmission(iv,joker,joker);
              mult( total_transmission(iv,joker,joker), tmp,
                    ppath_transmission(iv,joker,joker,ip-1) );
            }
        }
    }
}
 
 
 
 
Range get_rowindex_for_mblock( 
  const Sparse&   sensor_response, 
  const Index&    imblock )
{
  const Index   n1y = sensor_response.nrows();
  return Range( n1y*imblock, n1y );
}
 
 
 
 
void iy_transmission_for_scalar_tau( 
       Tensor3&     iy_transmission,
  const Index&      stokes_dim,
  ConstVectorView   tau )
{
  iy_transmission.resize( tau.nelem(), stokes_dim, stokes_dim );
  iy_transmission = 0;
  for( Index i=0; i<tau.nelem(); i++ )
    { 
      const Numeric t = exp( -tau[i] );
      for( Index is=0; is<stokes_dim; is++ )
        { 
          iy_transmission(i,is,is) = t;
        }
    }
}
 
 
 
 
void iy_transmission_mult( 
       Tensor3&      iy_trans_new,
  ConstTensor3View   iy_transmission,
  ConstTensor3View   trans )
{
  const Index nf = iy_transmission.npages();
  const Index ns = iy_transmission.ncols();
 
  assert( ns == iy_transmission.nrows() );
  assert( nf == trans.npages() );
  assert( ns == trans.nrows() );
  assert( ns == trans.ncols() );
 
  iy_trans_new.resize( nf, ns, ns );
 
  for( Index iv=0; iv<nf; iv++ )
    {
      mult( iy_trans_new(iv,joker,joker), iy_transmission(iv,joker,joker), 
                                                    trans(iv,joker,joker) );
    } 
}
 
 
 
 
void iy_transmission_mult_scalar_tau( 
       Tensor3&      iy_trans_new,
  ConstTensor3View   iy_transmission,
  ConstVectorView    tau )
{
  const Index nf = iy_transmission.npages();
 
  assert( iy_transmission.ncols() == iy_transmission.nrows() );
  assert( nf == tau.nelem() );
 
  iy_trans_new = iy_transmission;
 
  for( Index iv=0; iv<nf; iv++ )
    { iy_trans_new(iv,joker,joker) *= exp( -tau[iv] ); } 
}
 
 
 
 
void mirror_los(
        Vector&     los_mirrored,
  ConstVectorView   los, 
  const Index&      atmosphere_dim )
{
  los_mirrored.resize(2);
  //
  if( atmosphere_dim == 1 )
    { 
      los_mirrored[0] = 180 - los[0]; 
      los_mirrored[1] = 180; 
    }
  else if( atmosphere_dim == 2 )
    {
      los_mirrored[0] = 180 - fabs( los[0] ); 
      if( los[0] >= 0 )
        { los_mirrored[1] = 180; }
      else
        { los_mirrored[1] = 0; }
    }
  else if( atmosphere_dim == 3 )
    { 
      los_mirrored[0] = 180 - los[0]; 
      los_mirrored[1] = los[1] + 180; 
      if( los_mirrored[1] > 180 )
        { los_mirrored[1] -= 360; }
    }
}    
 
 
 
 
void rte_step_std(//Output and Input:
              VectorView stokes_vec,
              MatrixView trans_mat,
              //Input
              ConstMatrixView ext_mat_av,
              ConstVectorView abs_vec_av,
              ConstVectorView sca_vec_av,
              const Numeric& l_step,
              const Numeric& rte_planck_value )
{
  //Stokes dimension:
  Index stokes_dim = stokes_vec.nelem();
 
  //Check inputs:
  assert(is_size(trans_mat, stokes_dim, stokes_dim));
  assert(is_size(ext_mat_av, stokes_dim, stokes_dim));
  assert(is_size(abs_vec_av, stokes_dim));
  assert(is_size(sca_vec_av, stokes_dim));
  assert( rte_planck_value >= 0 );
  assert( l_step >= 0 );
  assert (!is_singular( ext_mat_av ));
 
  // Any changes here associated with the extinction matrix should also be
  // implemented in ext2mat.
 
  // Check, if only the first component of abs_vec is non-zero:
  bool unpol_abs_vec = true;
 
  for (Index i = 1; i < stokes_dim; i++)
    if (abs_vec_av[i] != 0)
      unpol_abs_vec = false;
 
  bool unpol_sca_vec = true;
 
  for (Index i = 1; i < stokes_dim; i++)
    if (sca_vec_av[i] != 0)
      unpol_sca_vec = false;
 
 
  //--- Scalar case: ---------------------------------------------------------
  if( stokes_dim == 1 )
    {
      trans_mat(0,0) = exp(-ext_mat_av(0,0) * l_step);
      stokes_vec[0]  = stokes_vec[0] * trans_mat(0,0) +
                       (abs_vec_av[0] * rte_planck_value + sca_vec_av[0]) /
        ext_mat_av(0,0) * (1 - trans_mat(0,0));
    }
 
 
  //--- Vector case: ---------------------------------------------------------
 
  // We have here two cases, diagonal or non-diagonal ext_mat_gas
  // For diagonal ext_mat_gas, we expect abs_vec_gas to only have a
  // non-zero value in position 1.
 
  //- Unpolarised
  else if( is_diagonal(ext_mat_av) && unpol_abs_vec && unpol_sca_vec )
    {
      trans_mat      = 0;
      trans_mat(0,0) = exp(-ext_mat_av(0,0) * l_step);
 
      // Stokes dim 1
      //   assert( ext_mat_av(0,0) == abs_vec_av[0] );
      //   Numeric transm = exp( -l_step * abs_vec_av[0] );
      stokes_vec[0] = stokes_vec[0] * trans_mat(0,0) +
                      (abs_vec_av[0] * rte_planck_value + sca_vec_av[0]) /
                      ext_mat_av(0,0) * (1 - trans_mat(0,0));
 
      // Stokes dims > 1
      for( Index i=1; i<stokes_dim; i++ )
        {
          //      assert( abs_vec_av[i] == 0.);
          trans_mat(i,i) = trans_mat(0,0);
          stokes_vec[i]  = stokes_vec[i] * trans_mat(i,i) +
                       sca_vec_av[i] / ext_mat_av(i,i)  * (1 - trans_mat(i,i));
        }
    }
 
 
  //- General case
  else
    {
      //Initialize internal variables:
 
      // Matrix LU used for LU decompostion and as dummy variable:
      Matrix LU(stokes_dim, stokes_dim);
      ArrayOfIndex indx(stokes_dim); // index for pivoting information
      Vector b(stokes_dim); // dummy variable
      Vector x(stokes_dim); // solution vector for K^(-1)*b
      Matrix I(stokes_dim, stokes_dim);
 
      Vector B_abs_vec(stokes_dim);
      B_abs_vec = abs_vec_av;
      B_abs_vec *= rte_planck_value;
 
      for (Index i=0; i<stokes_dim; i++)
        b[i] = B_abs_vec[i] + sca_vec_av[i];  // b = abs_vec * B + sca_vec
 
      // solve K^(-1)*b = x
      ludcmp(LU, indx, ext_mat_av);
      lubacksub(x, LU, b, indx);
 
      Matrix ext_mat_ds(stokes_dim, stokes_dim);
      ext_mat_ds = ext_mat_av;
      ext_mat_ds *= -l_step; // ext_mat_ds = -ext_mat*ds
 
      Index q = 10;  // index for the precision of the matrix exp function
      //Matrix exp_ext_mat(stokes_dim, stokes_dim);
      //matrix_exp(exp_ext_mat, ext_mat_ds, q);
      matrix_exp( trans_mat, ext_mat_ds, q);
 
      Vector term1(stokes_dim);
      Vector term2(stokes_dim);
 
      id_mat(I);
      for(Index i=0; i<stokes_dim; i++)
        {
          for(Index j=0; j<stokes_dim; j++)
            LU(i,j) = I(i,j) - trans_mat(i,j); // take LU as dummy variable
          // LU(i,j) = I(i,j) - exp_ext_mat(i,j); // take LU as dummy variable
        }
      mult(term2, LU, x); // term2: second term of the solution of the RTE with
                          //fixed scattered field
 
      // term1: first term of solution of the RTE with fixed scattered field
      //mult(term1, exp_ext_mat, stokes_vec);
      mult( term1, trans_mat, stokes_vec );
 
      for (Index i=0; i<stokes_dim; i++)
        stokes_vec[i] = term1[i] + term2[i];  // Compute the new Stokes Vector
    }
}
 
 
 
 
void surface_calc(
              Matrix&         iy,
        const Tensor3&        I,
        const Matrix&         surface_los,
        const Tensor4&        surface_rmatrix,
        const Matrix&         surface_emission )
{
  // Some sizes
  const Index   nf         = I.nrows();
  const Index   stokes_dim = I.ncols();
  const Index   nlos       = surface_los.nrows();
 
  iy = surface_emission;
  
  // Loop *surface_los*-es. If no such LOS, we are ready.
  if( nlos > 0 )
    {
      for( Index ilos=0; ilos<nlos; ilos++ )
        {
          Vector rtmp(stokes_dim);  // Reflected Stokes vector for 1 frequency
 
          for( Index iv=0; iv<nf; iv++ )
            {
          mult( rtmp, surface_rmatrix(ilos,iv,joker,joker), I(ilos,iv,joker) );
          iy(iv,joker) += rtmp;
            }
        }
    }
}
 
 
 
 
void trans_step_std(//Output and Input:
              VectorView stokes_vec,
              MatrixView trans_mat,
              //Input
              ConstMatrixView ext_mat_av,
              const Numeric& l_step )
{
  // Checks made in *ext2trans*.
  ext2trans( trans_mat, ext_mat_av, l_step );  
 
  Vector tmp = stokes_vec;
 
  mult( stokes_vec, trans_mat, tmp );
}