atm_funcs.cc

Go to the documentation of this file.

/* Copyright (C) 2000, 2001 Patrick Eriksson <patrick@rss.chalmers.se>
                            Stefan Buehler   <sbuehler@uni-bremen.de>

   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 <math.h>
#include <stdexcept>
#include "arts.h"
#include "matpackI.h"
#include "messages.h"          
#include "math_funcs.h"          
#include "make_vector.h"

extern const Numeric DEG2RAD;
extern const Numeric RAD2DEG;
extern const Numeric PLANCK_CONST;
extern const Numeric SPEED_OF_LIGHT;
extern const Numeric BOLTZMAN_CONST;



//   Physical functions

//
void planck (
             MatrixView      B, 
             ConstVectorView f,
             ConstVectorView t )
{
  // Double must be used here (if not, a becomes 0 when using float)
  static const double  a = 2.0*PLANCK_CONST/(SPEED_OF_LIGHT*SPEED_OF_LIGHT);
  static const double  b = PLANCK_CONST/BOLTZMAN_CONST;

  const Index    n_f  = f.nelem();
  const Index    n_t  = t.nelem();
  Index    i_f, i_t;
  Numeric   c, d;

  assert( n_f==B.nrows() );
  assert( n_t==B.ncols() );

  for ( i_f=0; i_f<n_f; i_f++ )
  {
    c = a * f[i_f]*f[i_f]*f[i_f];
    d = b * f[i_f];
    for ( i_t=0; i_t<n_t; i_t++ )
      B(i_f,i_t) = c / (exp(d/t[i_t]) - 1.0);
  }
}



//
void planck (
             VectorView    B,
             ConstVectorView    f,
             Numeric   t )
{
  // Double must be used here (if not, a becomes 0 when using float)
  static const double  a = 2.0*PLANCK_CONST/(SPEED_OF_LIGHT*SPEED_OF_LIGHT);
  static const double  b = PLANCK_CONST/BOLTZMAN_CONST;
         const double  c = b/t; 

  assert( B.nelem()==f.nelem() );

  for ( Index i=0; i<f.nelem(); i++ )
  {
    B[i] = a * f[i]*f[i]*f[i] / ( exp( f[i]*c ) - 1.0 );
  }
}



//
void invplanck (
                   VectorView   y,
              ConstVectorView   f,
              ConstVectorView   za )
{
  const Index   nf  = f.nelem();
  const Index   nza = za.nelem();
  const Index   ny  = y.nelem();
        Index   i0;

  // Use always double to avoid numerical problem (see invrayjean)
  const double   a = PLANCK_CONST/BOLTZMAN_CONST;
  const double   b = 2*PLANCK_CONST/(SPEED_OF_LIGHT*SPEED_OF_LIGHT);
        double   c,d;

  // Check input
  if ( max(y) > 1e-4 )  
    throw runtime_error("The spectrum cannot be in expected intensity unit "
                        "(impossible value detected).");
  //
  if ( nf*nza != ny )  
  {
    ostringstream os;
    os << "The length of *y* does not match *f_mono* and *za_pencil*.\n"
       << "y.nelem():         " << y.nelem() << "\n"
       << "Should be f_mono.nelem()*za_pencil.nelem(): "
       << f.nelem() * za.nelem() << "\n"
       << "f_mono.nelem():  " << f.nelem() << "\n"
       << "za_pencil.nelem(): " << za.nelem();
    throw runtime_error(os.str());
  }

  for ( Index i=0; i<nf; i++ )
  {
    c = a*f[i];
    d = b*f[i]*f[i]*f[i];
    for ( Index j=0; j<nza; j++ )    
    {
      i0 = j*nf + i;
      y[i0] = c / ( log(d/y[i0]+1) );
    }
  }
}



//
void invrayjean (
                   VectorView   y,
              ConstVectorView   f,
              ConstVectorView   za )
{
  const Index   nf  = f.nelem();
  const Index   nza = za.nelem();
  const Index   ny  = y.nelem();
        Index   i0;

  // The function returned NaNs when a and b were set to be Numeric (PE 010404)
  // Integrated this calculation into the for-loop to avoid numerical overflow
  // (OLE)
  //const double   a = SPEED_OF_LIGHT*SPEED_OF_LIGHT/(2*BOLTZMAN_CONST);
        double   b;


 // Commenting out the check of the input. This check is not any
 // longer true when the weighting functions are converted to
 // brightness temperature units. This occurs, e.g., when the line
 // intensity weighting functions are converted to temperature
 // units.

   // Check input
     //if ( max(y) > 1e-4 )  
     // throw runtime_error("The spectrum is not in expected intensity unit "
     //                    "(impossible value detected).");
  //
  if ( nf*nza != ny )  
  {
    ostringstream os;
    os << "The length of *y* does not match *f_mono* and *za_pencil*.\n"
       << "y.nelem():         " << y.nelem() << "\n"
       << "Should be f_mono.nelem()*za_pencil.nelem(): "
       << f.nelem() * za.nelem() << "\n"
       << "f_mono.nelem():  " << f.nelem() << "\n"
       << "za_pencil.nelem(): " << za.nelem();
    throw runtime_error(os.str());
  }

  for ( Index i=0; i<nf; i++ )
  {
    //b = a/(f[i]*f[i]);
    // a splitted and changed order of calculation to avoid
    // numerical instabilities. (OLE 2002-08-26)
    b = SPEED_OF_LIGHT / ((double)f[i] * (double)f[i])
      * SPEED_OF_LIGHT / (2 * BOLTZMAN_CONST);
    for ( Index j=0; j<nza; j++ )    
    {
      i0 = j*nf + i;
      y[i0] = b * (double)y[i0];
    }
  }
}



//
Numeric number_density (
                        Numeric   p,
                        Numeric   t )
{
  assert( 0!=t );
  return  p/t/BOLTZMAN_CONST;
}



//
Vector number_density (
                       ConstVectorView    p,
                       ConstVectorView    t )
{
  assert( p.nelem()==t.nelem() );

  // Calculate p / (t*BOLTZMAN_CONST):

  Vector dummy(p);              // Matpack can initialize a
                                // new Vector from another
                                // Vector.
  dummy /= t;                   // Element-vise divide by t.
  dummy /= BOLTZMAN_CONST;      // Divide all elements by BOLTZMAN_CONST.

  return dummy; 
}



//
Numeric g_of_z (
                Numeric   r_geoid,
                Numeric   g0,
                Numeric   z )
{
  return g0 * pow( r_geoid/(r_geoid+z), 2 );
}

//
Numeric g_of_lat (Numeric   latitude)
{
  extern const Numeric EARTH_EQUATORIAL_RADIUS;
  extern const Numeric EARTH_POLAR_RADIUS;
  extern const Numeric DEG2RAD;

  /* check latitude range */
  check_if_in_range( -90, 90, latitude, "latitude" );

  /* convert from geocentric to geographic */
  Numeric latG = atan( (EARTH_EQUATORIAL_RADIUS*EARTH_EQUATORIAL_RADIUS) /
                       (EARTH_POLAR_RADIUS*EARTH_POLAR_RADIUS) * tan(latitude*DEG2RAD));

  /* calculate in geographic coordinates */
  return 9.7803267714 * ( ( 1.0 + 0.0019318514 * sin(latG) * sin(latG) ) /
                          sqrt( 1.0 - 0.00669438 * sin(latG) * sin(latG) ) );
}





//   Core functions for RTE and BL 

//
void rte_iterate (
                  VectorView        y, 
                  const Index       start_index,
                  const Index       stop_index,
                  ConstMatrixView   tr,
                  ConstMatrixView   s,
                  const Index       n_f )
{
  Index   i_f;        // frequency index
  Index   i_z;        // LOS index
  Index   i_step;     // step order, -1 or 1

  if ( start_index >= stop_index )
    i_step = -1;

  else
    i_step = 1;

  for ( i_z=start_index; i_z!=(stop_index+i_step); i_z+=i_step ) 
  {
    for ( i_f=0; i_f<n_f; i_f++ )    
      y[i_f] = y[i_f]*tr(i_f,i_z) + s(i_f,i_z) * ( 1.0-tr(i_f,i_z) );
  }
}



//
void rte (
          VectorView        y,
          const Index       start_index,
          const Index       stop_index,
          ConstMatrixView   tr,
          ConstMatrixView   s,
          ConstVectorView   y_space,
          const Index       ground,
          ConstVectorView   e_ground,
          ConstVectorView   y_ground )
{
  const Index   n_f = tr.nrows();              // number of frequencies
  Index         i_f;                           // frequency index
  Index         i_break;                       // break index for looping
  Index         i_start;                       // variable for second loop

  // Init Y with Y_SPACE
  y = y_space;                  // Matpack can copy the contents of
                                // vectors like this. The dimensions
                                // must be the same! 

  // Check if LOS inside the atmosphere (if START_Index=0 -> Y=Y_SPACE)
  if ( start_index > 0 )
  {
    // Determine break index for looping, either 1 or the ground
    if ( ground > 0 )
      i_break = ground-1;
    else
      i_break = 0;       

    // Make first loop
    rte_iterate( y, start_index-1, i_break, tr, s, n_f );

    // We are now at the sensor, the ground or the tangent point
    // We are ready only if we are at the sensor.
    // If at sensor, we have that STOP_Index=0 and GROUND=0
    if ( !(stop_index==0 && ground==0) )
    {
      // Set most common values for I_START and I_BREAK
      i_start = 0;
      i_break = stop_index - 1;
      
      // If at the ground, include ground reflection. 
      // The loop can continue both downwards or upwards
      if ( ground > 0 )
      {      
        for ( i_f=0; i_f<n_f; i_f++ )    
          y[i_f] = y[i_f]*(1.0-e_ground[i_f]) + y_ground[i_f]*e_ground[i_f];
        
        if ( ground > 1 )  // 2D case, loop downwards
        {
         i_start = ground - 2;
         i_break = 0;
        }
      }

      // Make second loop
      rte_iterate( y, i_start, i_break, tr, s, n_f );

    } // second part
  } // if any values
}



//
void bl_iterate (
             VectorView   y,
       const Index   start_index,
       const Index   stop_index,
       ConstMatrixView   tr,
       const Index    n_f )
{
  Index   i_f;        // frequency index
  Index   i_z;        // LOS index
     Index   i_step;     // step order, -1 or 1

  if ( start_index >= stop_index )
    i_step = -1;
  else
    i_step = 1;

  for ( i_z=start_index; i_z!=(stop_index+i_step); i_z+=i_step ) 
  {
    for ( i_f=0; i_f<n_f; i_f++ )    
      y[i_f] *= tr(i_f,i_z);
  }
}



//
void bl (
             Vector&   y,
       const Index   start_index,
       const Index   stop_index,
       ConstMatrixView   tr,
       const Index    ground,
       ConstVectorView   e_ground )
{
  if ( start_index < stop_index )
    throw runtime_error("The start index cannot be "
                        "smaller than the stop index." );

  const Index   nf = tr.nrows();      // number of frequencies
  Index         iy;                   // frequency index

  // Init Y
  y.resize( nf );
  y = 1.0;

  // Loop steps passed twice
  if ( stop_index > 0 )
  {
    bl_iterate( y, 0, stop_index-1, tr, nf );
    y *= y;                     // Calculate the square of y element-vise.
  }

  // Loop remaining steps
  if ( start_index != stop_index )
    bl_iterate( y, stop_index, start_index-1, tr, nf );

  // Include effect of ground reflection
  if ( ground > 0 )
  {
    for ( iy=0; iy<nf; iy++ )    
      y[iy] *= ( 1.0 - e_ground[iy] );
  }
}



//   Conversion and interpolation of pressure and altitude grids.

//
void z2p(
         VectorView      p,
         ConstVectorView z0,
         ConstVectorView p0,
         ConstVectorView z )
{
  assert( p.nelem()==z.nelem() );
  if ( z.nelem() > 0 )
  {
    // Local variable to store log pressure grid:
    Vector logp0(p0.nelem());
    transform( logp0, log, p0 );        // This calculates logp0 = log(p0).

    interp_lin_vector( p, z0, logp0, z );
    transform( p, exp, p );             // This calculates p = exp(p).
  }
}



//
void interpp(
             VectorView          x, 
             ConstVectorView     p0,
             ConstVectorView     x0,
             ConstVectorView     p )
{
  assert( x.nelem()==p.nelem() );

  // Local variables to store log pressure grids:
  Vector logp0(p0.nelem());
  Vector logp(p.nelem());
  transform( logp0, log, p0 );  // This calculates logp0 = log(p0).
  transform( logp,  log, p  );  // This calculates logp  = log(p).

  interp_lin_vector( x, logp0, x0, logp );
}

void interpp_cloud(
                   VectorView     x, 
                   ConstVectorView     p0,
                   ConstVectorView     x0,
                   ConstVectorView     p )
{
  assert( x.nelem()==p.nelem() );

  interp_lin_vector( x, p0, x0, p );
}


//
void interpp(
             MatrixView  A,
             ConstVectorView  p0, 
             ConstMatrixView  A0, 
             ConstVectorView  p )
{
  assert( A.nrows()  == A0.nrows() );
  assert( A.ncols()  == p.nelem()  ); 

  // Local variables to store log pressure grids:
  Vector logp0(p0.nelem());
  Vector logp(p.nelem());
  transform( logp0, log, p0 );  // This calculates logp0 = log(p0).
  transform( logp,  log, p );   // This calculates logp0 = log(p0).

  interp_lin_matrix( A, logp0, A0, logp );
}



//
Numeric interpp(
                ConstVectorView     p0,
                ConstVectorView     x0,
                const Numeric       p )
{
  // Local variable to store log pressure grid:
  Vector logp0(p0.nelem());
  transform( logp0, log, p0 );  // This calculates logp0 = log(p0).

  return interp_lin( logp0, x0, log(p) );
}



//
void interpz(
             VectorView     x, 
             ConstVectorView     p0,
             ConstVectorView     z0,
             ConstVectorView     x0,
             ConstVectorView     z )
{
  assert( x.nelem()==z.nelem() ); 
  Vector p(z.nelem());
  z2p( p, z0, p0, z );
  interpp( x, p0, x0, p );
}



//
Numeric interpz(
        ConstVectorView     p0,
        ConstVectorView     z0,
        ConstVectorView     x0,
        const Numeric    z )
{
  Vector x(1);
  MakeVector Z(z);
  interpz( x, p0, z0, x0, Z );
  return x[0];
}



//   Tangent altitudes.

//
Numeric ztan_geom(
        const Numeric   za,
        const Numeric   z_plat,
        const Numeric   r_geoid )
{
  Numeric  z_tan;
  if ( za >= 90 )   
    z_tan = (r_geoid+z_plat)*sin(DEG2RAD*za) - r_geoid; 
  else
    z_tan = 999e3;
  return z_tan;
}




Numeric n_for_z(
        const Numeric      z,
        ConstVectorView       p_abs,
        ConstVectorView       z_abs,
        ConstVectorView       refr_index,
        const Numeric      atm_limit )

{
  if ( z > atm_limit )
    return 1.0;
  else
    return interpz( p_abs, z_abs, refr_index, z );
}



Numeric refr_constant( 
        const Numeric      r_geoid,
        const Numeric      za,
        const Numeric      z_plat,
        ConstVectorView       p_abs,
        ConstVectorView       z_abs,
        const Numeric      atm_limit,
        ConstVectorView       refr_index )
{
  Numeric n_plat = n_for_z( z_plat, p_abs, z_abs, refr_index, atm_limit );

  return (r_geoid + z_plat) * sin(DEG2RAD*za) * n_plat;
}



//
Numeric ztan_refr(
        const Numeric   c,
        const Numeric   za,
        const Numeric   z_plat,
        const Numeric   z_ground,
        ConstVectorView    p_abs,
        ConstVectorView    z_abs,
        ConstVectorView    refr_index,
        const Numeric   r_geoid )
{
  const Numeric atm_limit = last(z_abs);
  if ( za < 90 )   //=== Upward ==========================================
    return ztan_geom( za, z_plat, r_geoid );
  else
  {
    const Index  n = z_abs.nelem();
          Index  i;

    for ( i=(n-1); (i>=0) && (r_geoid+z_abs[i])*refr_index[i]>c; i-- )
    {
      if ( z_abs[i] <= z_ground ) //=== Ground intersection ==============
      {
        Numeric n_ground =  n_for_z(z_ground,p_abs,z_abs,refr_index,atm_limit);
        Numeric theta = RAD2DEG*asin(c/n_ground/(r_geoid+z_ground));
        return ztan_geom( 180-theta, z_ground, r_geoid );
      }
    }
    if ( i == (n-1) )  //=== outside the atmosphere ======================
      return ztan_geom( za, z_plat, r_geoid );
    else               //=== z_tan inside the atmosphere =================
    {
      Vector zs(2), cs(2);
      zs[0] = z_abs[i];
      zs[1] = z_abs[i+1];
      cs[0] = (r_geoid+z_abs[i])*refr_index[i];
      cs[1] = (r_geoid+z_abs[i+1])*refr_index[i+1];
      return interp_lin( cs, zs, c );
    }
  }
}
                           
//
void e_eq_water( VectorView          e_eq,
                  ConstVectorView    t )
{
  Index cur_elem;
  double a, b, c, d, e, t_cur;

  // make sure e_eq and t has the same length
  assert( e_eq.nelem() == t.nelem() );

  // Coefficients for Sontag's formula
  a = -6096.9385;
  b = 21.2409642;
  c = -2.711193e-2;
  d = 1.673952e-5;
  e = 2.433502;
  
  for ( cur_elem = 0; cur_elem < t.nelem(); cur_elem++ )
    {
      t_cur = t[cur_elem];
      e_eq[cur_elem] = exp( a / t_cur + b + c * t_cur + d * t_cur * t_cur + e * log( t_cur ) );
    }
}

//
void e_eq_ice( VectorView          e_eq,
                  ConstVectorView    t )
{
  Index cur_elem;
  double a, b, c, d, e, t_cur;
  
  // make sure e_eq and t has the same length
  assert( e_eq.nelem() == t.nelem() );

  // Coefficients for Sontag's formula
  a = -6024.5282;
  b = 29.32707;
  c = 1.0613868e-2;
  d = -1.3198825e-5;
  e = -0.49382577;

  for ( cur_elem = 0; cur_elem < t.nelem(); cur_elem++ )
    {
      t_cur = t[cur_elem];
      e_eq[cur_elem] = exp( a / t_cur + b + c * t_cur + d * t_cur * t_cur + e * log( t_cur ) );
    }
}

---