ARTS  2.2.66
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 
20 /*===========================================================================
21  === File description
22  ===========================================================================*/
23 
34 /*===========================================================================
35  === External declarations
36  ===========================================================================*/
37 
38 #include <cmath>
39 #include "auto_md.h"
40 #include "complex.h"
41 #include "interpolation.h"
42 #include "geodetic.h"
43 #include "refraction.h"
44 #include "special_interp.h"
45 
46 extern const Numeric DEG2RAD;
47 extern const Numeric RAD2DEG;
48 extern const Numeric TEMP_0_C;
49 
50 
51 
52 /*===========================================================================
53  === The functions (in alphabetical order)
54  ===========================================================================*/
55 
56 
58 
79  Matrix& complex_n,
80  const Vector& f_grid,
81  const Numeric& t )
82 {
83  chk_if_in_range( "t", t, TEMP_0_C-40, TEMP_0_C+100 );
84  chk_if_in_range( "min of f_grid", min(f_grid), 10e9, 1000e9 );
85  chk_if_in_range( "max of f_grid", max(f_grid), 10e9, 1000e9 );
86 
87  const Index nf = f_grid.nelem();
88 
89  complex_n.resize( nf, 2 );
90 
91  // Implementation following epswater93.m (by C. Mätzler), part of Atmlab,
92  // but numeric values strictly following the paper version (146, not 146.4)
93  const Numeric theta = 1 - 300 / t;
94  const Numeric e0 = 77.66 - 103.3 * theta;
95  const Numeric e1 = 0.0671 * e0;
96  const Numeric f1 = 20.2 + 146 * theta + 316 * theta * theta;
97  const Numeric e2 = 3.52;
98  const Numeric f2 = 39.8 * f1;
99 
100  for( Index iv=0; iv<nf; iv++ )
101  {
102  const Complex ifGHz( 0.0, f_grid[iv]/1e9 );
103 
104  Complex n = sqrt( e2 + (e1-e2) / (Numeric(1.0)-ifGHz/f2) +
105  (e0-e1) / (Numeric(1.0)-ifGHz/f1) );
106 
107  complex_n(iv,0) = n.real();
108  complex_n(iv,1) = n.imag();
109  }
110 }
111 
112 
113 
115 
142  Workspace& ws,
143  Numeric& refr_index_air,
144  Numeric& refr_index_air_group,
145  const Agenda& refr_index_air_agenda,
146  ConstVectorView p_grid,
147  ConstVectorView refellipsoid,
148  ConstTensor3View z_field,
149  ConstTensor3View t_field,
150  ConstTensor4View vmr_field,
151  ConstVectorView f_grid,
152  const Numeric& r )
153 {
154  Numeric rtp_pressure, rtp_temperature;
155  Vector rtp_vmr;
156 
157  // Pressure grid position
158  ArrayOfGridPos gp(1);
159  gridpos( gp, z_field(joker,0,0), Vector( 1, r - refellipsoid[0] ) );
160 
161  // Altitude interpolation weights
162  Matrix itw(1,2);
163  interpweights( itw, gp );
164 
165  // Pressure
166  Vector dummy(1);
167  itw2p( dummy, p_grid, gp, itw );
168  rtp_pressure = dummy[0];
169 
170  // Temperature
171  interp( dummy, itw, t_field(joker,0,0), gp );
172  rtp_temperature = dummy[0];
173 
174  // VMR
175  const Index ns = vmr_field.nbooks();
176  //
177  rtp_vmr.resize(ns);
178  //
179  for( Index is=0; is<ns; is++ )
180  {
181  interp( dummy, itw, vmr_field(is,joker,0,0), gp );
182  rtp_vmr[is] = dummy[0];
183  }
184 
185  refr_index_air_agendaExecute( ws, refr_index_air, refr_index_air_group,
186  rtp_pressure, rtp_temperature, rtp_vmr,
187  f_grid, refr_index_air_agenda );
188 }
189 
190 
191 
193 
222  Workspace& ws,
223  Numeric& refr_index_air,
224  Numeric& refr_index_air_group,
225  const Agenda& refr_index_air_agenda,
226  ConstVectorView p_grid,
227  ConstVectorView lat_grid,
228  ConstVectorView refellipsoid,
229  ConstTensor3View z_field,
230  ConstTensor3View t_field,
231  ConstTensor4View vmr_field,
232  ConstVectorView f_grid,
233  const Numeric& r,
234  const Numeric& lat )
235 {
236  Numeric rtp_pressure, rtp_temperature;
237  Vector rtp_vmr;
238 
239  // Determine the geometric altitudes at *lat*
240  const Index np = p_grid.nelem();
241  Vector z_grid(np);
242  ArrayOfGridPos gp_lat(1);
243  //
244  gridpos( gp_lat, lat_grid, lat );
245  z_at_lat_2d( z_grid, p_grid, lat_grid, z_field(joker,joker,0), gp_lat[0] );
246 
247  // Determine the ellipsoid radius at *lat*
248  const Numeric rellips = refell2d( refellipsoid, lat_grid, gp_lat[0] );
249 
250  // Altitude (equal to pressure) grid position
251  ArrayOfGridPos gp_p(1);
252  gridpos( gp_p, z_grid, Vector( 1, r - rellips ) );
253 
254  // Altitude interpolation weights
255  Matrix itw(1,2);
256  Vector dummy(1);
257  interpweights( itw, gp_p );
258 
259  // Pressure
260  itw2p( dummy, p_grid, gp_p, itw );
261  rtp_pressure = dummy[0];
262 
263  // Temperature
264  itw.resize(1,4);
265  interpweights( itw, gp_p, gp_lat );
266  interp( dummy, itw, t_field(joker,joker,0), gp_p, gp_lat );
267  rtp_temperature = dummy[0];
268 
269  // VMR
270  const Index ns = vmr_field.nbooks();
271  //
272  rtp_vmr.resize(ns);
273  //
274  for( Index is=0; is<ns; is++ )
275  {
276  interp( dummy, itw, vmr_field(is,joker,joker,0), gp_p, gp_lat );
277  rtp_vmr[is] = dummy[0];
278  }
279 
280 
281  refr_index_air_agendaExecute( ws, refr_index_air, refr_index_air_group,
282  rtp_pressure, rtp_temperature, rtp_vmr,
283  f_grid, refr_index_air_agenda );
284 }
285 
286 
287 
316  Workspace& ws,
317  Numeric& refr_index_air,
318  Numeric& refr_index_air_group,
319  const Agenda& refr_index_air_agenda,
320  ConstVectorView p_grid,
321  ConstVectorView lat_grid,
322  ConstVectorView lon_grid,
323  ConstVectorView refellipsoid,
324  ConstTensor3View z_field,
325  ConstTensor3View t_field,
326  ConstTensor4View vmr_field,
327  ConstVectorView f_grid,
328  const Numeric& r,
329  const Numeric& lat,
330  const Numeric& lon )
331 {
332  Numeric rtp_pressure, rtp_temperature;
333  Vector rtp_vmr;
334 
335  // Determine the geometric altitudes at *lat* and *lon*
336  const Index np = p_grid.nelem();
337  Vector z_grid(np);
338  ArrayOfGridPos gp_lat(1), gp_lon(1);
339  //
340  gridpos( gp_lat, lat_grid, lat );
341  gridpos( gp_lon, lon_grid, lon );
342  z_at_latlon( z_grid, p_grid, lat_grid, lon_grid, z_field,
343  gp_lat[0], gp_lon[0] );
344 
345  // Determine the elipsoid radius at *lat*
346  const Numeric rellips = refell2d( refellipsoid, lat_grid, gp_lat[0] );
347 
348  // Altitude (equal to pressure) grid position
349  ArrayOfGridPos gp_p(1);
350  gridpos( gp_p, z_grid, Vector( 1, r - rellips ) );
351 
352  // Altitude interpolation weights
353  Matrix itw(1,2);
354  Vector dummy(1);
355  interpweights( itw, gp_p );
356 
357  // Pressure
358  itw2p( dummy, p_grid, gp_p, itw );
359  rtp_pressure = dummy[0];
360 
361  // Temperature
362  itw.resize(1,8);
363  interpweights( itw, gp_p, gp_lat, gp_lon );
364  interp( dummy, itw, t_field, gp_p, gp_lat, gp_lon );
365  rtp_temperature = dummy[0];
366 
367  // VMR
368  const Index ns = vmr_field.nbooks();
369  //
370  rtp_vmr.resize(ns);
371  //
372  for( Index is=0; is<ns; is++ )
373  {
374  interp( dummy, itw, vmr_field(is,joker,joker,joker), gp_p, gp_lat,
375  gp_lon );
376  rtp_vmr[is] = dummy[0];
377  }
378 
379  refr_index_air_agendaExecute( ws, refr_index_air, refr_index_air_group,
380  rtp_pressure, rtp_temperature, rtp_vmr,
381  f_grid, refr_index_air_agenda );
382 }
383 
384 
385 
387 
413  Workspace& ws,
414  Numeric& refr_index_air,
415  Numeric& refr_index_air_group,
416  Numeric& dndr,
417  const Agenda& refr_index_air_agenda,
418  ConstVectorView p_grid,
419  ConstVectorView refellipsoid,
420  ConstTensor3View z_field,
421  ConstTensor3View t_field,
422  ConstTensor4View vmr_field,
423  ConstVectorView f_grid,
424  const Numeric& r )
425 {
426  get_refr_index_1d( ws, refr_index_air, refr_index_air_group,
427  refr_index_air_agenda, p_grid, refellipsoid,
428  z_field, t_field, vmr_field, f_grid, r );
429 
430  const Numeric n0 = refr_index_air;
431  Numeric dummy;
432 
433  get_refr_index_1d( ws, refr_index_air, dummy,
434  refr_index_air_agenda, p_grid, refellipsoid,
435  z_field, t_field, vmr_field, f_grid, r+1 );
436 
437  dndr = refr_index_air - n0;
438 
439  refr_index_air = n0;
440 }
441 
442 
443 
445 
478  Workspace& ws,
479  Numeric& refr_index_air,
480  Numeric& refr_index_air_group,
481  Numeric& dndr,
482  Numeric& dndlat,
483  const Agenda& refr_index_air_agenda,
484  ConstVectorView p_grid,
485  ConstVectorView lat_grid,
486  ConstVectorView refellipsoid,
487  ConstTensor3View z_field,
488  ConstTensor3View t_field,
489  ConstTensor4View vmr_field,
490  ConstVectorView f_grid,
491  const Numeric& r,
492  const Numeric& lat )
493 {
494  get_refr_index_2d( ws, refr_index_air, refr_index_air_group,
495  refr_index_air_agenda, p_grid, lat_grid, refellipsoid,
496  z_field, t_field, vmr_field, f_grid, r, lat );
497 
498  const Numeric n0 = refr_index_air;
499  Numeric dummy;
500 
501  get_refr_index_2d( ws, refr_index_air, dummy, refr_index_air_agenda, p_grid,
502  lat_grid, refellipsoid, z_field, t_field, vmr_field,
503  f_grid, r+1, lat );
504 
505  dndr = refr_index_air - n0;
506 
507  const Numeric dlat = 1e-4;
508 
509  get_refr_index_2d( ws, refr_index_air, dummy, refr_index_air_agenda, p_grid,
510  lat_grid, refellipsoid, z_field, t_field, vmr_field,
511  f_grid, r, lat+dlat );
512 
513  dndlat = ( refr_index_air - n0 ) / ( DEG2RAD * dlat * r );
514 
515  refr_index_air = n0;
516 }
517 
518 
519 
521 
558  Workspace& ws,
559  Numeric& refr_index_air,
560  Numeric& refr_index_air_group,
561  Numeric& dndr,
562  Numeric& dndlat,
563  Numeric& dndlon,
564  const Agenda& refr_index_air_agenda,
565  ConstVectorView p_grid,
566  ConstVectorView lat_grid,
567  ConstVectorView lon_grid,
568  ConstVectorView refellipsoid,
569  ConstTensor3View z_field,
570  ConstTensor3View t_field,
571  ConstTensor4View vmr_field,
572  ConstVectorView f_grid,
573  const Numeric& r,
574  const Numeric& lat,
575  const Numeric& lon )
576 {
577  get_refr_index_3d( ws, refr_index_air, refr_index_air_group,
578  refr_index_air_agenda, p_grid, lat_grid, lon_grid,
579  refellipsoid, z_field, t_field, vmr_field, f_grid,
580  r, lat, lon );
581 
582  const Numeric n0 = refr_index_air;
583  Numeric dummy;
584 
585  get_refr_index_3d( ws, refr_index_air, dummy, refr_index_air_agenda, p_grid,
586  lat_grid, lon_grid, refellipsoid, z_field, t_field,
587  vmr_field, f_grid, r+1, lat, lon );
588 
589  dndr = refr_index_air - n0;
590 
591  const Numeric dlat = 1e-4;
592 
593  get_refr_index_3d( ws, refr_index_air, dummy, refr_index_air_agenda, p_grid,
594  lat_grid, lon_grid, refellipsoid, z_field, t_field,
595  vmr_field, f_grid, r, lat+dlat, lon );
596 
597  dndlat = ( refr_index_air - n0 ) / ( DEG2RAD * dlat * r );
598 
599  const Numeric dlon = 1e-4;
600 
601  get_refr_index_3d( ws, refr_index_air, dummy, refr_index_air_agenda, p_grid,
602  lat_grid, lon_grid, refellipsoid, z_field, t_field,
603  vmr_field, f_grid, r, lat, lon+dlon);
604 
605  dndlon = ( refr_index_air - n0 ) /
606  ( DEG2RAD * dlon * r * cos( DEG2RAD*lat ) );
607 
608  refr_index_air = n0;
609 }
610 
611 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
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:221
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:315
const Numeric DEG2RAD
A class implementing complex numbers for ARTS.
The Agenda class.
Definition: agenda_class.h:44
The Vector class.
Definition: matpackI.h:556
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Header file for interpolation.cc.
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition: geodetic.cc:1244
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:141
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
A constant view of a Tensor4.
Definition: matpackIV.h:141
void complex_n_water_liebe93(Matrix &complex_n, const Vector &f_grid, const Numeric &t)
complex_n_water_liebe93
Definition: refraction.cc:78
const Numeric TEMP_0_C
std::complex< Numeric > Complex
Definition: complex.h:32
#define max(a, b)
Definition: continua.cc:20461
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
const Joker joker
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
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:14800
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:557
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)
z_at_latlon
void z_at_lat_2d(VectorView z, ConstVectorView p_grid, ConstVectorView lat_grid, ConstMatrixView z_field, const GridPos &gp_lat)
z_at_lat_2d
Header file for special_interp.cc.
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:101
const Numeric RAD2DEG
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
#define ns
Definition: continua.cc:21931
#define min(a, b)
Definition: continua.cc:20460
A constant view of a Tensor3.
Definition: matpackIII.h:139
A constant view of a Vector.
Definition: matpackI.h:292
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:412
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:477
Workspace class.
Definition: workspace_ng.h:47
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
Refraction functions.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
itw2p
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580