refraction.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2012 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 "refraction.h"
#include <cmath>
#include "auto_md.h"
#include "complex.h"
#include "geodetic.h"
#include "interpolation.h"
#include "special_interp.h"

extern const Numeric DEG2RAD;
extern const Numeric RAD2DEG;
extern const Numeric TEMP_0_C;

/*===========================================================================
  === The functions (in alphabetical order)
  ===========================================================================*/


void complex_n_water_liebe93(Matrix& complex_n,
                             const Vector& f_grid,
                             const Numeric& t) {
  chk_if_in_range("t", t, TEMP_0_C - 40, TEMP_0_C + 100);
  chk_if_in_range("min of f_grid", min(f_grid), 10e9, 1000e9);
  chk_if_in_range("max of f_grid", max(f_grid), 10e9, 1000e9);

  const Index nf = f_grid.nelem();

  complex_n.resize(nf, 2);

  // Implementation following epswater93.m (by C. Mätzler), part of Atmlab,
  // but numeric values strictly following the paper version (146, not 146.4)
  const Numeric theta = 1 - 300 / t;
  const Numeric e0 = 77.66 - 103.3 * theta;
  const Numeric e1 = 0.0671 * e0;
  const Numeric f1 = 20.2 + 146 * theta + 316 * theta * theta;
  const Numeric e2 = 3.52;
  const Numeric f2 = 39.8 * f1;

  for (Index iv = 0; iv < nf; iv++) {
    const Complex ifGHz(0.0, f_grid[iv] / 1e9);

    Complex n = sqrt(e2 + (e1 - e2) / (Numeric(1.0) - ifGHz / f2) +
                     (e0 - e1) / (Numeric(1.0) - ifGHz / f1));

    complex_n(iv, 0) = n.real();
    complex_n(iv, 1) = n.imag();
  }
}


void complex_n_ice_matzler06(Matrix& complex_n,
                             const Vector& f_grid,
                             const Numeric& t) {
  chk_if_in_range("t", t, 20., 280.);
  chk_if_in_range("min of f_grid", min(f_grid), 10e6, 3000e9);
  chk_if_in_range("max of f_grid", max(f_grid), 10e6, 3000e9);

  const Index nf = f_grid.nelem();

  complex_n.resize(nf, 2);

  // some parametrization constants
  const Numeric B1 = 0.0207;
  const Numeric B2 = 1.16e-11;
  const Numeric b = 335.;

  const Numeric deltabeta = exp(-9.963 + 0.0372 * (t - 273));
  const Numeric ebdt = exp(b / t);
  const Numeric betam = (B1 / t) * ebdt / ((ebdt - 1.) * (ebdt - 1.));

  const Numeric theta = 300. / t - 1;
  const Numeric alfa = (0.00504 + 0.0062 * theta) * exp(-22.1 * theta);
  const Numeric reps = 3.1884 + 9.1e-4 * (t - 273);

  for (Index iv = 0; iv < nf; iv++) {
    Numeric f = f_grid[iv] / 1e9;
    Numeric beta = betam + B2 * f * f + deltabeta;
    Numeric ieps = alfa / f + beta * f;

    Complex eps(reps, ieps);
    Complex n = sqrt(eps);
    complex_n(iv, 0) = n.real();
    complex_n(iv, 1) = n.imag();
  }
}


void get_refr_index_1d(Workspace& ws,
                       Numeric& refr_index_air,
                       Numeric& refr_index_air_group,
                       const Agenda& refr_index_air_agenda,
                       ConstVectorView p_grid,
                       ConstVectorView refellipsoid,
                       ConstTensor3View z_field,
                       ConstTensor3View t_field,
                       ConstTensor4View vmr_field,
                       ConstVectorView f_grid,
                       const Numeric& r) {
  Numeric rtp_pressure, rtp_temperature;
  Vector rtp_vmr;

  // Pressure grid position
  ArrayOfGridPos gp(1);
  gridpos(gp, z_field(joker, 0, 0), Vector(1, r - refellipsoid[0]));

  // Altitude interpolation weights
  Matrix itw(1, 2);
  interpweights(itw, gp);

  // Pressure
  Vector dummy(1);
  itw2p(dummy, p_grid, gp, itw);
  rtp_pressure = dummy[0];

  // Temperature
  interp(dummy, itw, t_field(joker, 0, 0), gp);
  rtp_temperature = dummy[0];

  // VMR
  const Index ns = vmr_field.nbooks();
  //
  rtp_vmr.resize(ns);
  //
  for (Index is = 0; is < ns; is++) {
    interp(dummy, itw, vmr_field(is, joker, 0, 0), gp);
    rtp_vmr[is] = dummy[0];
  }

  refr_index_air_agendaExecute(ws,
                               refr_index_air,
                               refr_index_air_group,
                               rtp_pressure,
                               rtp_temperature,
                               rtp_vmr,
                               f_grid,
                               refr_index_air_agenda);
}


void get_refr_index_2d(Workspace& ws,
                       Numeric& refr_index_air,
                       Numeric& refr_index_air_group,
                       const Agenda& refr_index_air_agenda,
                       ConstVectorView p_grid,
                       ConstVectorView lat_grid,
                       ConstVectorView refellipsoid,
                       ConstTensor3View z_field,
                       ConstTensor3View t_field,
                       ConstTensor4View vmr_field,
                       ConstVectorView f_grid,
                       const Numeric& r,
                       const Numeric& lat) {
  Numeric rtp_pressure, rtp_temperature;
  Vector rtp_vmr;

  // 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(joker, joker, 0), gp_lat[0]);

  // Determine the ellipsoid radius at *lat*
  const Numeric rellips = refell2d(refellipsoid, lat_grid, gp_lat[0]);

  // Altitude (equal to pressure) grid position
  ArrayOfGridPos gp_p(1);
  gridpos(gp_p, z_grid, Vector(1, r - rellips));

  // Altitude interpolation weights
  Matrix itw(1, 2);
  Vector dummy(1);
  interpweights(itw, gp_p);

  // Pressure
  itw2p(dummy, p_grid, gp_p, itw);
  rtp_pressure = dummy[0];

  // Temperature
  itw.resize(1, 4);
  interpweights(itw, gp_p, gp_lat);
  interp(dummy, itw, t_field(joker, joker, 0), gp_p, gp_lat);
  rtp_temperature = dummy[0];

  // VMR
  const Index ns = vmr_field.nbooks();
  //
  rtp_vmr.resize(ns);
  //
  for (Index is = 0; is < ns; is++) {
    interp(dummy, itw, vmr_field(is, joker, joker, 0), gp_p, gp_lat);
    rtp_vmr[is] = dummy[0];
  }

  refr_index_air_agendaExecute(ws,
                               refr_index_air,
                               refr_index_air_group,
                               rtp_pressure,
                               rtp_temperature,
                               rtp_vmr,
                               f_grid,
                               refr_index_air_agenda);
}

void get_refr_index_3d(Workspace& ws,
                       Numeric& refr_index_air,
                       Numeric& refr_index_air_group,
                       const Agenda& refr_index_air_agenda,
                       ConstVectorView p_grid,
                       ConstVectorView lat_grid,
                       ConstVectorView lon_grid,
                       ConstVectorView refellipsoid,
                       ConstTensor3View z_field,
                       ConstTensor3View t_field,
                       ConstTensor4View vmr_field,
                       ConstVectorView f_grid,
                       const Numeric& r,
                       const Numeric& lat,
                       const Numeric& lon) {
  Numeric rtp_pressure, rtp_temperature;
  Vector rtp_vmr;

  // 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 elipsoid radius at *lat*
  const Numeric rellips = refell2d(refellipsoid, lat_grid, gp_lat[0]);

  // Altitude (equal to pressure) grid position
  ArrayOfGridPos gp_p(1);
  gridpos(gp_p, z_grid, Vector(1, r - rellips));

  // Altitude interpolation weights
  Matrix itw(1, 2);
  Vector dummy(1);
  interpweights(itw, gp_p);

  // Pressure
  itw2p(dummy, p_grid, gp_p, itw);
  rtp_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);
  rtp_temperature = dummy[0];

  // VMR
  const Index ns = vmr_field.nbooks();
  //
  rtp_vmr.resize(ns);
  //
  for (Index is = 0; is < ns; is++) {
    interp(
        dummy, itw, vmr_field(is, joker, joker, joker), gp_p, gp_lat, gp_lon);
    rtp_vmr[is] = dummy[0];
  }

  refr_index_air_agendaExecute(ws,
                               refr_index_air,
                               refr_index_air_group,
                               rtp_pressure,
                               rtp_temperature,
                               rtp_vmr,
                               f_grid,
                               refr_index_air_agenda);
}


void refr_gradients_1d(Workspace& ws,
                       Numeric& refr_index_air,
                       Numeric& refr_index_air_group,
                       Numeric& dndr,
                       const Agenda& refr_index_air_agenda,
                       ConstVectorView p_grid,
                       ConstVectorView refellipsoid,
                       ConstTensor3View z_field,
                       ConstTensor3View t_field,
                       ConstTensor4View vmr_field,
                       ConstVectorView f_grid,
                       const Numeric& r) {
  get_refr_index_1d(ws,
                    refr_index_air,
                    refr_index_air_group,
                    refr_index_air_agenda,
                    p_grid,
                    refellipsoid,
                    z_field,
                    t_field,
                    vmr_field,
                    f_grid,
                    r);

  const Numeric n0 = refr_index_air;
  Numeric dummy;

  get_refr_index_1d(ws,
                    refr_index_air,
                    dummy,
                    refr_index_air_agenda,
                    p_grid,
                    refellipsoid,
                    z_field,
                    t_field,
                    vmr_field,
                    f_grid,
                    r + 1);

  dndr = refr_index_air - n0;

  refr_index_air = n0;
}


void refr_gradients_2d(Workspace& ws,
                       Numeric& refr_index_air,
                       Numeric& refr_index_air_group,
                       Numeric& dndr,
                       Numeric& dndlat,
                       const Agenda& refr_index_air_agenda,
                       ConstVectorView p_grid,
                       ConstVectorView lat_grid,
                       ConstVectorView refellipsoid,
                       ConstTensor3View z_field,
                       ConstTensor3View t_field,
                       ConstTensor4View vmr_field,
                       ConstVectorView f_grid,
                       const Numeric& r,
                       const Numeric& lat) {
  get_refr_index_2d(ws,
                    refr_index_air,
                    refr_index_air_group,
                    refr_index_air_agenda,
                    p_grid,
                    lat_grid,
                    refellipsoid,
                    z_field,
                    t_field,
                    vmr_field,
                    f_grid,
                    r,
                    lat);

  const Numeric n0 = refr_index_air;
  Numeric dummy;

  get_refr_index_2d(ws,
                    refr_index_air,
                    dummy,
                    refr_index_air_agenda,
                    p_grid,
                    lat_grid,
                    refellipsoid,
                    z_field,
                    t_field,
                    vmr_field,
                    f_grid,
                    r + 1,
                    lat);

  dndr = refr_index_air - n0;

  const Numeric dlat = 1e-4;

  get_refr_index_2d(ws,
                    refr_index_air,
                    dummy,
                    refr_index_air_agenda,
                    p_grid,
                    lat_grid,
                    refellipsoid,
                    z_field,
                    t_field,
                    vmr_field,
                    f_grid,
                    r,
                    lat + dlat);

  dndlat = (refr_index_air - n0) / (DEG2RAD * dlat * r);

  refr_index_air = n0;
}


void refr_gradients_3d(Workspace& ws,
                       Numeric& refr_index_air,
                       Numeric& refr_index_air_group,
                       Numeric& dndr,
                       Numeric& dndlat,
                       Numeric& dndlon,
                       const Agenda& refr_index_air_agenda,
                       ConstVectorView p_grid,
                       ConstVectorView lat_grid,
                       ConstVectorView lon_grid,
                       ConstVectorView refellipsoid,
                       ConstTensor3View z_field,
                       ConstTensor3View t_field,
                       ConstTensor4View vmr_field,
                       ConstVectorView f_grid,
                       const Numeric& r,
                       const Numeric& lat,
                       const Numeric& lon) {
  get_refr_index_3d(ws,
                    refr_index_air,
                    refr_index_air_group,
                    refr_index_air_agenda,
                    p_grid,
                    lat_grid,
                    lon_grid,
                    refellipsoid,
                    z_field,
                    t_field,
                    vmr_field,
                    f_grid,
                    r,
                    lat,
                    lon);

  const Numeric n0 = refr_index_air;
  Numeric dummy;

  get_refr_index_3d(ws,
                    refr_index_air,
                    dummy,
                    refr_index_air_agenda,
                    p_grid,
                    lat_grid,
                    lon_grid,
                    refellipsoid,
                    z_field,
                    t_field,
                    vmr_field,
                    f_grid,
                    r + 1,
                    lat,
                    lon);

  dndr = refr_index_air - n0;

  const Numeric dlat = 1e-4;

  get_refr_index_3d(ws,
                    refr_index_air,
                    dummy,
                    refr_index_air_agenda,
                    p_grid,
                    lat_grid,
                    lon_grid,
                    refellipsoid,
                    z_field,
                    t_field,
                    vmr_field,
                    f_grid,
                    r,
                    lat + dlat,
                    lon);

  dndlat = (refr_index_air - n0) / (DEG2RAD * dlat * r);

  const Numeric dlon = 1e-4;

  get_refr_index_3d(ws,
                    refr_index_air,
                    dummy,
                    refr_index_air_agenda,
                    p_grid,
                    lat_grid,
                    lon_grid,
                    refellipsoid,
                    z_field,
                    t_field,
                    vmr_field,
                    f_grid,
                    r,
                    lat,
                    lon + dlon);

  dndlon = (refr_index_air - n0) / (DEG2RAD * dlon * r * cos(DEG2RAD * lat));

  refr_index_air = n0;
}