refraction.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2008 Patrick Eriksson <patrick.eriksson@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
  ===========================================================================*/
 
#include <cmath>
#include "auto_md.h"
#include "interpolation.h"
#include "refraction.h"
#include "special_interp.h"
 
extern const Numeric DEG2RAD;
extern const Numeric RAD2DEG;
 
 
 
/*===========================================================================
  === The functions (in alphabetical order)
  ===========================================================================*/
 
 
 
void get_refr_index_1d(
              Workspace&  ws,
              Numeric&    refr_index,
              Numeric&    a_pressure,
              Numeric&    a_temperature,
              Vector&     a_vmr_list,
        const Agenda&     refr_index_agenda,
        ConstVectorView   p_grid,
        const Numeric&    r_geoid,
        ConstVectorView   z_field,
        ConstVectorView   t_field,
        ConstMatrixView   vmr_field,
        const Numeric&    r )
{ 
  // Altitude (equal to pressure) grid position
  ArrayOfGridPos   gp(1);
  gridpos( gp, z_field, Vector( 1, r - r_geoid ) );
 
  // Altitude interpolation weights
  Matrix   itw(1,2);
  interpweights( itw, gp );
 
  // Pressure
  Vector   dummy(1);
  itw2p( dummy, p_grid, gp, itw );
  a_pressure = dummy[0];
 
  // Temperature
  interp( dummy, itw, t_field, gp );
  a_temperature = dummy[0];
 
  // VMR
  const Index   ns = vmr_field.nrows();
  //
  a_vmr_list.resize(ns);
  //
  for( Index is=0; is<ns; is++ )
    {
      interp( dummy, itw, vmr_field(is,joker), gp );
      a_vmr_list[is] = dummy[0];
    }
 
  refr_index_agendaExecute( ws,
                            refr_index, a_pressure, a_temperature, a_vmr_list,
                            refr_index_agenda );
}
 
 
 
 
void get_refr_index_2d(
              Workspace&  ws,
              Numeric&    refr_index,
              Numeric&    a_pressure,
              Numeric&    a_temperature,
              Vector&     a_vmr_list,
        const Agenda&     refr_index_agenda,
        ConstVectorView   p_grid,
        ConstVectorView   lat_grid,
        ConstVectorView   r_geoid,
        ConstMatrixView   z_field,
        ConstMatrixView   t_field,
        ConstTensor3View  vmr_field,
        const Numeric&    r,
        const Numeric&    lat )
{ 
  // Determine the geometric altitudes at *lat*
  const Index      np = p_grid.nelem();
  Vector           z_grid(np);
  ArrayOfGridPos   gp_lat(1);
  //
  gridpos( gp_lat, lat_grid, lat );
  z_at_lat_2d( z_grid, p_grid, lat_grid, z_field, gp_lat[0] );
 
  // Determine the geoid radius at *lat*
  Matrix   itw(1,2);
  Vector   dummy(1);
  interpweights( itw, gp_lat );
  interp( dummy, itw, r_geoid, gp_lat );
  const Numeric   rgeoid = dummy[0];
 
  // Altitude (equal to pressure) grid position
  ArrayOfGridPos   gp_p(1);
  gridpos( gp_p, z_grid, Vector( 1, r - rgeoid ) );
 
  // Altitude interpolation weights
  interpweights( itw, gp_p );
 
  // Pressure
  itw2p( dummy, p_grid, gp_p, itw );
  a_pressure = dummy[0];
 
  // Temperature
  itw.resize(1,4);
  interpweights( itw, gp_p, gp_lat );
  interp( dummy, itw, t_field, gp_p, gp_lat );
  a_temperature = dummy[0];
 
  // VMR
  const Index   ns = vmr_field.npages();
  //
  a_vmr_list.resize(ns);
  //
  for( Index is=0; is<ns; is++ )
    {
      interp( dummy, itw, vmr_field(is,joker,joker), gp_p, gp_lat );
      a_vmr_list[is] = dummy[0];
    }
 
  refr_index_agendaExecute( ws,
                            refr_index, a_pressure, a_temperature, a_vmr_list,
                            refr_index_agenda );
}
 
 
 
void get_refr_index_3d(
              Workspace&  ws,
              Numeric&    refr_index,
              Numeric&    a_pressure,
              Numeric&    a_temperature,
              Vector&     a_vmr_list,
        const Agenda&     refr_index_agenda,
        ConstVectorView   p_grid,
        ConstVectorView   lat_grid,
        ConstVectorView   lon_grid,
        ConstMatrixView   r_geoid,
        ConstTensor3View  z_field,
        ConstTensor3View  t_field,
        ConstTensor4View  vmr_field,
        const Numeric&    r,
        const Numeric&    lat,
        const Numeric&    lon )
{ 
  // Determine the geometric altitudes at *lat* and *lon*
  const Index      np = p_grid.nelem();
  Vector           z_grid(np);
  ArrayOfGridPos   gp_lat(1), gp_lon(1);
  //
  gridpos( gp_lat, lat_grid, lat );
  gridpos( gp_lon, lon_grid, lon );
  z_at_latlon( z_grid, p_grid, lat_grid, lon_grid, z_field, 
                                                        gp_lat[0], gp_lon[0] );
  
  // Determine the geoid radius at *lat*
  Matrix   itw(1,4);
  Vector   dummy(1);
  interpweights( itw, gp_lat, gp_lon );
  interp( dummy, itw, r_geoid, gp_lat, gp_lon );
  const Numeric   rgeoid = dummy[0];
 
  // Altitude (equal to pressure) grid position
  ArrayOfGridPos   gp_p(1);
  gridpos( gp_p, z_grid, Vector( 1, r - rgeoid ) );
 
  // Altitude interpolation weights
  itw.resize(1,2);
  interpweights( itw, gp_p );
 
  // Pressure
  itw2p( dummy, p_grid, gp_p, itw );
  a_pressure = dummy[0];
 
  // Temperature
  itw.resize(1,8);
  interpweights( itw, gp_p, gp_lat, gp_lon );
  interp( dummy, itw, t_field, gp_p, gp_lat, gp_lon );
  a_temperature = dummy[0];
 
  // VMR
  const Index   ns = vmr_field.nbooks();
  //
  a_vmr_list.resize(ns);
  //
  for( Index is=0; is<ns; is++ )
    {
      interp( dummy, itw, vmr_field(is,joker,joker,joker), gp_p, gp_lat, 
                                                                      gp_lon );
      a_vmr_list[is] = dummy[0];
    }
 
  refr_index_agendaExecute( ws,
                            refr_index, a_pressure, a_temperature, a_vmr_list,
                            refr_index_agenda );
}
 
 
 
 
void refr_gradients_1d(
              Workspace&  ws,
              Numeric&    refr_index,
              Numeric&    dndr,
              Numeric&    a_pressure,
              Numeric&    a_temperature,
              Vector&     a_vmr_list,
        const Agenda&     refr_index_agenda,
        ConstVectorView   p_grid,
        const Numeric&    r_geoid,
        ConstVectorView   z_field,
        ConstVectorView   t_field,
        ConstMatrixView   vmr_field,
        const Numeric&    r )
{ 
   get_refr_index_1d( ws, refr_index, a_pressure,  a_temperature, a_vmr_list, 
                      refr_index_agenda, p_grid, 
                      r_geoid, z_field, t_field, vmr_field, r );
 
   const Numeric   n0 = refr_index;
 
   get_refr_index_1d( ws, refr_index, a_pressure, a_temperature, a_vmr_list, 
                      refr_index_agenda, p_grid, 
                      r_geoid, z_field, t_field, vmr_field, r+1 );
 
   dndr = refr_index - n0;
 
   refr_index = n0;
}
 
 
 
 
void refr_gradients_2d(
              Workspace&  ws,
              Numeric&    refr_index,
              Numeric&    dndr,
              Numeric&    dndlat,
              Numeric&    a_pressure,
              Numeric&    a_temperature,
              Vector&     a_vmr_list,
        const Agenda&     refr_index_agenda,
        ConstVectorView   p_grid,
        ConstVectorView   lat_grid,
        ConstVectorView   r_geoid,
        ConstMatrixView   z_field,
        ConstMatrixView   t_field,
        ConstTensor3View  vmr_field,
        const Numeric&    r,
        const Numeric&    lat )
{ 
   get_refr_index_2d( ws, refr_index, a_pressure,  a_temperature, a_vmr_list, 
                      refr_index_agenda, p_grid, lat_grid,
                      r_geoid, z_field, t_field, vmr_field, r, lat );
 
   const Numeric   n0 = refr_index;
 
   get_refr_index_2d( ws, refr_index, a_pressure, a_temperature, a_vmr_list, 
                      refr_index_agenda, p_grid, lat_grid, r_geoid, 
                      z_field, t_field, vmr_field, r+1, lat );
 
   dndr = refr_index - n0;
 
   const Numeric   dlat = 1e-4;
 
   get_refr_index_2d( ws, refr_index, a_pressure, a_temperature, a_vmr_list, 
                      refr_index_agenda, p_grid, lat_grid, r_geoid, 
                      z_field, t_field, vmr_field, r, lat+dlat );
 
   dndlat = ( refr_index - n0 ) / ( DEG2RAD * dlat * r ); 
 
   refr_index = n0;
}
 
 
 
 
void refr_gradients_3d(
              Workspace&  ws,
              Numeric&    refr_index,
              Numeric&    dndr,
              Numeric&    dndlat,
              Numeric&    dndlon,
              Numeric&    a_pressure,
              Numeric&    a_temperature,
              Vector&     a_vmr_list,
        const Agenda&     refr_index_agenda,
        ConstVectorView   p_grid,
        ConstVectorView   lat_grid,
        ConstVectorView   lon_grid,
        ConstMatrixView   r_geoid,
        ConstTensor3View  z_field,
        ConstTensor3View  t_field,
        ConstTensor4View  vmr_field,
        const Numeric&    r,
        const Numeric&    lat,
        const Numeric&    lon )
{ 
   get_refr_index_3d( ws, refr_index, a_pressure, a_temperature, a_vmr_list, 
               refr_index_agenda, p_grid, lat_grid, 
               lon_grid, r_geoid, z_field, t_field, vmr_field, r, lat, lon );
 
   const Numeric   n0 = refr_index;
 
   get_refr_index_3d( ws, refr_index, a_pressure, a_temperature, a_vmr_list, 
               refr_index_agenda, p_grid, lat_grid, 
               lon_grid, r_geoid, z_field, t_field, vmr_field, r+1, lat, lon );
 
   dndr = refr_index - n0;
 
   const Numeric   dlat = 1e-4;
 
   get_refr_index_3d( ws, refr_index, a_pressure, a_temperature, a_vmr_list, 
                      refr_index_agenda, p_grid, lat_grid, 
                      lon_grid, r_geoid, z_field, t_field, vmr_field, 
                      r, lat+dlat, lon );
 
   dndlat = ( refr_index - n0 ) / ( DEG2RAD * dlat * r ); 
 
   const Numeric   dlon = 1e-4;
 
   get_refr_index_3d( ws, refr_index, a_pressure, a_temperature, a_vmr_list, 
                      refr_index_agenda, p_grid, lat_grid, 
                      lon_grid, r_geoid, z_field, t_field, vmr_field, 
                      r, lat, lon+dlon);
 
   dndlon = ( refr_index - n0 ) / ( DEG2RAD * dlon * r * cos( DEG2RAD*lat ) ); 
 
   refr_index = n0;
}
 
 
 
 
void refr_index_thayer_1974(
              Numeric&   refr_index,
        const Numeric&   p,
        const Numeric&   t,
        const Numeric&   h2o_vmr )
{
  const Numeric   e = p * h2o_vmr;
 
  refr_index = 1 + ( 77.6e-8 * ( p - e ) +
                                          ( 64.8e-8 + 3.776e-3 / t ) * e ) / t;
}
 
 
void refr_index_ir(
              Numeric&   refr_index,
        const Numeric&   p,
        const Numeric&   t )
{
  const Numeric bn0 = 1.000272620045304;
  const Numeric bk = 288.16 * (pow(bn0,Numeric(2.0))-1.0)/
                                        (1013.25*(pow(bn0,Numeric(2.0))+2.0));
 
  // Pa -> HPa
  refr_index = sqrt((2.0*bk*p/100.0+t)/(t-bk*p/100.0));
}