ARTS  2.4.0(git:4fb77825)
refraction.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Patrick Eriksson <patrick.eriksson@chalmers.se>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
18 /*===========================================================================
19  === File description
20  ===========================================================================*/
21 
30 /*===========================================================================
31  === External declarations
32  ===========================================================================*/
33 
34 #include "refraction.h"
35 #include <cmath>
36 #include "auto_md.h"
37 #include "complex.h"
38 #include "geodetic.h"
39 #include "interpolation.h"
40 #include "special_interp.h"
41 
42 extern const Numeric DEG2RAD;
43 extern const Numeric RAD2DEG;
44 extern const Numeric TEMP_0_C;
45 
46 /*===========================================================================
47  === The functions (in alphabetical order)
48  ===========================================================================*/
49 
51 
72  const Vector& f_grid,
73  const Numeric& t) {
74  chk_if_in_range("t", t, TEMP_0_C - 40, TEMP_0_C + 100);
75  chk_if_in_range("min of f_grid", min(f_grid), 10e9, 1000e9);
76  chk_if_in_range("max of f_grid", max(f_grid), 10e9, 1000e9);
77 
78  const Index nf = f_grid.nelem();
79 
80  complex_n.resize(nf, 2);
81 
82  // Implementation following epswater93.m (by C. Mätzler), part of Atmlab,
83  // but numeric values strictly following the paper version (146, not 146.4)
84  const Numeric theta = 1 - 300 / t;
85  const Numeric e0 = 77.66 - 103.3 * theta;
86  const Numeric e1 = 0.0671 * e0;
87  const Numeric f1 = 20.2 + 146 * theta + 316 * theta * theta;
88  const Numeric e2 = 3.52;
89  const Numeric f2 = 39.8 * f1;
90 
91  for (Index iv = 0; iv < nf; iv++) {
92  const Complex ifGHz(0.0, f_grid[iv] / 1e9);
93 
94  Complex n = sqrt(e2 + (e1 - e2) / (Numeric(1.0) - ifGHz / f2) +
95  (e0 - e1) / (Numeric(1.0) - ifGHz / f1));
96 
97  complex_n(iv, 0) = n.real();
98  complex_n(iv, 1) = n.imag();
99  }
100 }
101 
103 
122  const Vector& f_grid,
123  const Numeric& t) {
124  chk_if_in_range("t", t, 20., 280.);
125  chk_if_in_range("min of f_grid", min(f_grid), 10e6, 3000e9);
126  chk_if_in_range("max of f_grid", max(f_grid), 10e6, 3000e9);
127 
128  const Index nf = f_grid.nelem();
129 
130  complex_n.resize(nf, 2);
131 
132  // some parametrization constants
133  const Numeric B1 = 0.0207;
134  const Numeric B2 = 1.16e-11;
135  const Numeric b = 335.;
136 
137  const Numeric deltabeta = exp(-9.963 + 0.0372 * (t - 273));
138  const Numeric ebdt = exp(b / t);
139  const Numeric betam = (B1 / t) * ebdt / ((ebdt - 1.) * (ebdt - 1.));
140 
141  const Numeric theta = 300. / t - 1;
142  const Numeric alfa = (0.00504 + 0.0062 * theta) * exp(-22.1 * theta);
143  const Numeric reps = 3.1884 + 9.1e-4 * (t - 273);
144 
145  for (Index iv = 0; iv < nf; iv++) {
146  Numeric f = f_grid[iv] / 1e9;
147  Numeric beta = betam + B2 * f * f + deltabeta;
148  Numeric ieps = alfa / f + beta * f;
149 
150  Complex eps(reps, ieps);
151  Complex n = sqrt(eps);
152  complex_n(iv, 0) = n.real();
153  complex_n(iv, 1) = n.imag();
154  }
155 }
156 
158 
194  const Numeric& r) {
196  Vector rtp_vmr;
197 
198  // Pressure grid position
199  ArrayOfGridPos gp(1);
200  gridpos(gp, z_field(joker, 0, 0), Vector(1, r - refellipsoid[0]));
201 
202  // Altitude interpolation weights
203  Matrix itw(1, 2);
204  interpweights(itw, gp);
205 
206  // Pressure
207  Vector dummy(1);
208  itw2p(dummy, p_grid, gp, itw);
209  rtp_pressure = dummy[0];
210 
211  // Temperature
212  interp(dummy, itw, t_field(joker, 0, 0), gp);
213  rtp_temperature = dummy[0];
214 
215  // VMR
216  const Index ns = vmr_field.nbooks();
217  //
218  rtp_vmr.resize(ns);
219  //
220  for (Index is = 0; is < ns; is++) {
221  interp(dummy, itw, vmr_field(is, joker, 0, 0), gp);
222  rtp_vmr[is] = dummy[0];
223  }
224 
228  rtp_pressure,
230  rtp_vmr,
231  f_grid,
233 }
234 
236 
275  const Numeric& r,
276  const Numeric& lat) {
278  Vector rtp_vmr;
279 
280  // Determine the geometric altitudes at *lat*
281  const Index np = p_grid.nelem();
282  Vector z_grid(np);
283  ArrayOfGridPos gp_lat(1);
284  //
285  gridpos(gp_lat, lat_grid, lat);
286  z_at_lat_2d(z_grid, p_grid, lat_grid, z_field(joker, joker, 0), gp_lat[0]);
287 
288  // Determine the ellipsoid radius at *lat*
289  const Numeric rellips = refell2d(refellipsoid, lat_grid, gp_lat[0]);
290 
291  // Altitude (equal to pressure) grid position
292  ArrayOfGridPos gp_p(1);
293  gridpos(gp_p, z_grid, Vector(1, r - rellips));
294 
295  // Altitude interpolation weights
296  Matrix itw(1, 2);
297  Vector dummy(1);
298  interpweights(itw, gp_p);
299 
300  // Pressure
301  itw2p(dummy, p_grid, gp_p, itw);
302  rtp_pressure = dummy[0];
303 
304  // Temperature
305  itw.resize(1, 4);
306  interpweights(itw, gp_p, gp_lat);
307  interp(dummy, itw, t_field(joker, joker, 0), gp_p, gp_lat);
308  rtp_temperature = dummy[0];
309 
310  // VMR
311  const Index ns = vmr_field.nbooks();
312  //
313  rtp_vmr.resize(ns);
314  //
315  for (Index is = 0; is < ns; is++) {
316  interp(dummy, itw, vmr_field(is, joker, joker, 0), gp_p, gp_lat);
317  rtp_vmr[is] = dummy[0];
318  }
319 
323  rtp_pressure,
325  rtp_vmr,
326  f_grid,
328 }
329 
369  const Numeric& r,
370  const Numeric& lat,
371  const Numeric& lon) {
373  Vector rtp_vmr;
374 
375  // Determine the geometric altitudes at *lat* and *lon*
376  const Index np = p_grid.nelem();
377  Vector z_grid(np);
378  ArrayOfGridPos gp_lat(1), gp_lon(1);
379  //
380  gridpos(gp_lat, lat_grid, lat);
381  gridpos(gp_lon, lon_grid, lon);
382  z_at_latlon(
383  z_grid, p_grid, lat_grid, lon_grid, z_field, gp_lat[0], gp_lon[0]);
384 
385  // Determine the elipsoid radius at *lat*
386  const Numeric rellips = refell2d(refellipsoid, lat_grid, gp_lat[0]);
387 
388  // Altitude (equal to pressure) grid position
389  ArrayOfGridPos gp_p(1);
390  gridpos(gp_p, z_grid, Vector(1, r - rellips));
391 
392  // Altitude interpolation weights
393  Matrix itw(1, 2);
394  Vector dummy(1);
395  interpweights(itw, gp_p);
396 
397  // Pressure
398  itw2p(dummy, p_grid, gp_p, itw);
399  rtp_pressure = dummy[0];
400 
401  // Temperature
402  itw.resize(1, 8);
403  interpweights(itw, gp_p, gp_lat, gp_lon);
404  interp(dummy, itw, t_field, gp_p, gp_lat, gp_lon);
405  rtp_temperature = dummy[0];
406 
407  // VMR
408  const Index ns = vmr_field.nbooks();
409  //
410  rtp_vmr.resize(ns);
411  //
412  for (Index is = 0; is < ns; is++) {
413  interp(
414  dummy, itw, vmr_field(is, joker, joker, joker), gp_p, gp_lat, gp_lon);
415  rtp_vmr[is] = dummy[0];
416  }
417 
421  rtp_pressure,
423  rtp_vmr,
424  f_grid,
426 }
427 
429 
457  Numeric& dndr,
465  const Numeric& r) {
470  p_grid,
471  refellipsoid,
472  z_field,
473  t_field,
474  vmr_field,
475  f_grid,
476  r);
477 
478  const Numeric n0 = refr_index_air;
479  Numeric dummy;
480 
483  dummy,
485  p_grid,
486  refellipsoid,
487  z_field,
488  t_field,
489  vmr_field,
490  f_grid,
491  r + 1);
492 
493  dndr = refr_index_air - n0;
494 
495  refr_index_air = n0;
496 }
497 
499 
534  Numeric& dndr,
535  Numeric& dndlat,
544  const Numeric& r,
545  const Numeric& lat) {
550  p_grid,
551  lat_grid,
552  refellipsoid,
553  z_field,
554  t_field,
555  vmr_field,
556  f_grid,
557  r,
558  lat);
559 
560  const Numeric n0 = refr_index_air;
561  Numeric dummy;
562 
565  dummy,
567  p_grid,
568  lat_grid,
569  refellipsoid,
570  z_field,
571  t_field,
572  vmr_field,
573  f_grid,
574  r + 1,
575  lat);
576 
577  dndr = refr_index_air - n0;
578 
579  const Numeric dlat = 1e-4;
580 
583  dummy,
585  p_grid,
586  lat_grid,
587  refellipsoid,
588  z_field,
589  t_field,
590  vmr_field,
591  f_grid,
592  r,
593  lat + dlat);
594 
595  dndlat = (refr_index_air - n0) / (DEG2RAD * dlat * r);
596 
597  refr_index_air = n0;
598 }
599 
601 
640  Numeric& dndr,
641  Numeric& dndlat,
642  Numeric& dndlon,
652  const Numeric& r,
653  const Numeric& lat,
654  const Numeric& lon) {
659  p_grid,
660  lat_grid,
661  lon_grid,
662  refellipsoid,
663  z_field,
664  t_field,
665  vmr_field,
666  f_grid,
667  r,
668  lat,
669  lon);
670 
671  const Numeric n0 = refr_index_air;
672  Numeric dummy;
673 
676  dummy,
678  p_grid,
679  lat_grid,
680  lon_grid,
681  refellipsoid,
682  z_field,
683  t_field,
684  vmr_field,
685  f_grid,
686  r + 1,
687  lat,
688  lon);
689 
690  dndr = refr_index_air - n0;
691 
692  const Numeric dlat = 1e-4;
693 
696  dummy,
698  p_grid,
699  lat_grid,
700  lon_grid,
701  refellipsoid,
702  z_field,
703  t_field,
704  vmr_field,
705  f_grid,
706  r,
707  lat + dlat,
708  lon);
709 
710  dndlat = (refr_index_air - n0) / (DEG2RAD * dlat * r);
711 
712  const Numeric dlon = 1e-4;
713 
716  dummy,
718  p_grid,
719  lat_grid,
720  lon_grid,
721  refellipsoid,
722  z_field,
723  t_field,
724  vmr_field,
725  f_grid,
726  r,
727  lat,
728  lon + dlon);
729 
730  dndlon = (refr_index_air - n0) / (DEG2RAD * dlon * r * cos(DEG2RAD * lat));
731 
732  refr_index_air = n0;
733 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::z_field
Tensor3 z_field(Workspace &ws) noexcept
Definition: autoarts.h:7690
ARTS::Var::lat
Numeric lat(Workspace &ws) noexcept
Definition: autoarts.h:3933
auto_md.h
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
itw2p
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
Definition: special_interp.cc:718
refr_index_air_agendaExecute
void refr_index_air_agendaExecute(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Numeric rtp_pressure, const Numeric rtp_temperature, const Vector &rtp_vmr, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:24874
ARTS::Var::lat_grid
Vector lat_grid(Workspace &ws) noexcept
Definition: autoarts.h:3962
joker
const Joker joker
interpolation.h
Header file for interpolation.cc.
refell2d
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition: geodetic.cc:1174
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
ARTS::Var::refr_index_air_agenda
Agenda refr_index_air_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5539
z_at_latlon
void z_at_latlon(VectorView z, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, const GridPos &gp_lat, const GridPos &gp_lon)
Returns the geomtrical altitudes of p_grid for one latitude and one longitude.
Definition: special_interp.cc:888
RAD2DEG
const Numeric RAD2DEG
ARTS::Var::refr_index_air
Numeric refr_index_air(Workspace &ws) noexcept
Definition: autoarts.h:5532
get_refr_index_1d
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)
get_refr_index_1d
Definition: refraction.cc:184
ARTS::Var::rtp_vmr
Vector rtp_vmr(Workspace &ws) noexcept
Definition: autoarts.h:5826
get_refr_index_2d
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)
get_refr_index_2d
Definition: refraction.cc:264
complex_n_water_liebe93
void complex_n_water_liebe93(Matrix &complex_n, const Vector &f_grid, const Numeric &t)
complex_n_water_liebe93
Definition: refraction.cc:71
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
z_at_lat_2d
void z_at_lat_2d(VectorView z, ConstVectorView p_grid, ConstVectorView lat_grid, ConstMatrixView z_field, const GridPos &gp_lat)
Returns the geomtrical altitudes of p_grid for one latitude.
Definition: special_interp.cc:860
Agenda
The Agenda class.
Definition: agenda_class.h:44
ARTS::Var::lon
Numeric lon(Workspace &ws) noexcept
Definition: autoarts.h:4064
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
ConstTensor4View
A constant view of a Tensor4.
Definition: matpackIV.h:133
complex.h
A class implementing complex numbers for ARTS.
Array
This can be used to make arrays out of anything.
Definition: array.h:108
min
#define min(a, b)
Definition: legacy_continua.cc:20628
ARTS::Var::vmr_field
Tensor4 vmr_field(Workspace &ws) noexcept
Definition: autoarts.h:7130
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
refr_gradients_1d
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)
refr_gradients_1d
Definition: refraction.cc:454
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
ARTS::Var::refr_index_air_group
Numeric refr_index_air_group(Workspace &ws) noexcept
Definition: autoarts.h:5553
refr_gradients_2d
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)
refr_gradients_2d
Definition: refraction.cc:531
geodetic.h
ARTS::Var::refellipsoid
Vector refellipsoid(Workspace &ws) noexcept
Definition: autoarts.h:5519
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
max
#define max(a, b)
Definition: legacy_continua.cc:20629
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
ARTS::Var::rtp_pressure
Numeric rtp_pressure(Workspace &ws) noexcept
Definition: autoarts.h:5790
ARTS::Var::lon_grid
Vector lon_grid(Workspace &ws) noexcept
Definition: autoarts.h:4090
TEMP_0_C
const Numeric TEMP_0_C
oem::Vector
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
ARTS::Var::rtp_temperature
Numeric rtp_temperature(Workspace &ws) noexcept
Definition: autoarts.h:5807
refraction.h
Refraction functions.
ARTS::Var::t_field
Tensor3 t_field(Workspace &ws) noexcept
Definition: autoarts.h:6947
Workspace
Workspace class.
Definition: workspace_ng.h:40
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
special_interp.h
Header file for special_interp.cc.
get_refr_index_3d
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)
Definition: refraction.cc:357
refr_gradients_3d
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)
refr_gradients_3d
Definition: refraction.cc:637
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
DEG2RAD
const Numeric DEG2RAD
chk_if_in_range
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:89
Vector
The Vector class.
Definition: matpackI.h:860
beta
#define beta
Definition: legacy_continua.cc:21435
complex_n_ice_matzler06
void complex_n_ice_matzler06(Matrix &complex_n, const Vector &f_grid, const Numeric &t)
complex_n_ice_matzler06
Definition: refraction.cc:121
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
ns
#define ns
Definition: legacy_continua.cc:22183