ARTS  2.4.0(git:4fb77825)
m_refraction.cc
Go to the documentation of this file.
1 /* Copyright (C) 2003-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 
33 /*===========================================================================
34  === External declarations
35  ===========================================================================*/
36 
37 #include <cmath>
38 #include "abs_species_tags.h"
39 #include "absorption.h"
40 #include "arts.h"
41 #include "check_input.h"
42 #include "math_funcs.h"
43 #include "matpackI.h"
44 #include "messages.h"
45 #include "physics_funcs.h"
46 #include "refraction.h"
47 #include "special_interp.h"
48 
49 extern const Numeric ELECTRON_CHARGE;
50 extern const Numeric ELECTRON_MASS;
51 extern const Numeric PI;
52 extern const Numeric VACUUM_PERMITTIVITY;
53 extern const Numeric TORR2PA;
54 
55 /*===========================================================================
56  === WSMs for refr_index_air
57  ===========================================================================*/
58 
59 /* Workspace method: Doxygen documentation will be auto-generated */
62  const Vector& f_grid,
64  const Vector& rtp_vmr,
65  const Index& demand_vmr_value,
66  const Verbosity&) {
67  // The expression used is found in many textbooks, e.g. Rybicki and Lightman
68  // (1979). Note that the refractive index corresponds to the phase velocity.
69 
70  static const Numeric k = ELECTRON_CHARGE * ELECTRON_CHARGE /
72 
73  Numeric edensity = 0;
74 
75  Index ife = -1;
76  for (Index sp = 0; sp < abs_species.nelem() && ife < 0; sp++) {
78  ife = sp;
79  }
80  }
81 
82  if (ife < 0) {
83  if (demand_vmr_value) {
84  throw runtime_error(
85  "Free electrons not found in *abs_species* and "
86  "contribution to refractive index can not be calculated.");
87  }
88  } else {
89  edensity = rtp_vmr[ife];
90 
91  if (edensity > 0) {
92  // Check that lowest frequency not too low
93  // Limit at 100 GHz follows Hartmann and Leitinger, Range errors due
94  // to ionospheric and tropospheric effects for signal frequencies
95  // above 100 HMHz, Bull. Goed., 1984.
96  if (f_grid[0] < 100e6) {
97  throw runtime_error(
98  "All frequencies must be >= 100 MHz, but "
99  "this is not the case.");
100  }
101  if (edensity * k / (f_grid[0] * f_grid[0]) > 0.25) {
102  ostringstream os;
103  os << "All frequencies must at least be twice the plasma frequency.\n"
104  << "For this particular point, the plasma frequency is: "
105  << sqrt(edensity * k) / 1e6 << " MHz.";
106  throw runtime_error(os.str());
107  }
108 
109  const Numeric f = (f_grid[0] + last(f_grid)) / 2.0;
110  const Numeric a = edensity * k / (f * f);
111  const Numeric n = sqrt(1 - a);
112 
113  refr_index_air += n - 1;
114  refr_index_air_group += 1 / n - 1;
115  }
116  }
117 }
118 
119 /* Workspace method: Doxygen documentation will be auto-generated */
122  const Numeric& rtp_pressure,
123  const Numeric& rtp_temperature,
124  const Verbosity&) {
125  static const Numeric bn0 = 1.000272620045304;
126  static const Numeric bn02 = bn0 * bn0;
127  static const Numeric bk = 288.16 * (bn02 - 1.0) / (1013.25 * (bn02 + 2.0));
128 
129  // Pa -> hPa
130  const Numeric n = sqrt((2.0 * bk * rtp_pressure / 100.0 + rtp_temperature) /
131  (rtp_temperature - bk * rtp_pressure / 100.0)) -
132  1;
133 
134  refr_index_air += n;
136 }
137 
138 /* Workspace method: Doxygen documentation will be auto-generated */
141  const Numeric& rtp_pressure,
142  const Numeric& rtp_temperature,
143  const Vector& rtp_vmr,
145  const Numeric& k1,
146  const Numeric& k2,
147  const Numeric& k3,
148  const Verbosity&) {
149  if (abs_species.nelem() != rtp_vmr.nelem())
150  throw runtime_error(
151  "The number of tag groups differ between "
152  "*rtp_vmr* and *abs_species*.");
153 
154  Index firstH2O = find_first_species_tg(
156 
157  Numeric e;
158  if (firstH2O < 0)
159  //throw runtime_error(
160  // "Water vapour is a required (must be a tag group in *abs_species*)." );
161  e = 0.;
162  else
163  e = rtp_pressure * rtp_vmr[firstH2O];
164 
165  const Numeric n =
166  (k1 * (rtp_pressure - e) + (k2 + k3 / rtp_temperature) * e) /
168 
169  refr_index_air += n;
171 }
172 
173 /* Workspace method: Doxygen documentation will be auto-generated */
177  const Numeric& rtp_pressure,
178  const Numeric& rtp_temperature,
179  const Vector& rtp_vmr,
181  const Verbosity&) {
182  //FIXME: Shall n be rescaled for sum(VMW)=1? Doing so now, but is it correct?
183  // Short sensitivity test (tropical, dry air, ~11km tanh) shows that
184  // vmr-normalized n fits significantly better (dtanh~0.5m) than
185  // non-normalized one (dtanh~5m).
186 
187  /*
188  for now, hard-coding the reference refindices and refT/p. could make that
189  some re-setabel parameters (like iso ratios... also regarding storing them in
190  file/data struct per species)
191 */
192  const Numeric p0 = 760. * TORR2PA;
193  const Numeric T0 = 273.15;
194 
195  // Number of refractive species:
196  const Index nrs = 6;
197 
198  // This is hardwired here and quite primitive, but should do the job.
199  // Set refractive index species names.
200  ArrayOfString ref_spec_names(nrs);
201  ref_spec_names[0] = "N2";
202  ref_spec_names[1] = "O2";
203  ref_spec_names[2] = "CO2";
204  ref_spec_names[3] = "H2";
205  ref_spec_names[4] = "He";
206  ref_spec_names[5] = "H2O";
207 
208  // Set reference refractive indices
209  // Values from Newell and Baird, 1965
210  Vector ref_n(nrs);
211  ref_n[0] = 293.81e-6;
212  ref_n[1] = 266.95e-6;
213  ref_n[2] = 495.16e-6;
214  ref_n[3] = 135.77e-6;
215  ref_n[4] = 34.51e-6;
216  // value for H2O from H2O contribution according to refr_index_airMicrowavesEarth
217  // at reference conditions
218  // that is: n_H2O = p_ref/T_ref * (k2 + k3/Tref)
219  ref_n[5] = 5338.89e-6;
220 
221  // Checks
222  if (abs_species.nelem() != rtp_vmr.nelem())
223  throw runtime_error(
224  "The number of tag groups differ between "
225  "*rtp_vmr* and *abs_species*.");
226  /*
227  further checks:
228  ? non-neg T
229  ?
230 */
231 
232  // Data management
233  // Find the location of all refractive species in abs_species. Set to -1 if
234  // not found. The length of array ref_spec_locations is the number of
235  // considered refractive species (in this method: N2, O2, CO2, H2, He).
236  // The value means:
237  // -1 = not in abs_species
238  // N = species is number N in abs_species
239 
240  //Can't use this one as it inside gets the broadening species names and
241  //number. Also, we would have to make a workaround for this_species. So,
242  //instead we use a modified version of this function directly included here.
243  /*find_broad_spec_locations(ref_spec_locations,
244  abs_species,
245  this_species);*/
246 
247  ArrayOfIndex ref_spec_locations(nrs);
248 
249  // Loop over all broadening species and see if we can find them in abs_species.
250  for (Index i = 0; i < nrs; ++i) {
251  // Find associated internal species index (we do the lookup by index, not by name).
252  const Index isi = species_index_from_species_name(ref_spec_names[i]);
253 
254  // Find position of broadening species isi in abs_species. The called
255  // function returns -1 if not found, which is already the correct
256  // treatment for this case that we also want here.
257  ref_spec_locations[i] = find_first_species_tg(abs_species, isi);
258  }
259 
260  // The actual calculation
261  // N_tot = sum (Nref_i * p_i/p_0 * T0/T)
262  // = sum (Nref_i * vmr_i*p/p_0 * T0/T)
263  // = p/p_0 * T0/T * sum ( Nref_i * vmr_i)
264 
265  const Numeric ratioT = T0 / rtp_temperature;
266  const Numeric ratiop = rtp_pressure / p0;
267 
268  Numeric ref_spec_vmr_sum = 0.;
269  Numeric n = 0.;
270 
271  // Add up refractive species, where available:
272  for (Index i = 0; i < nrs; ++i) {
273  if (ref_spec_locations[i] >= 0) {
274  // Add to VMR sum:
275  ref_spec_vmr_sum += rtp_vmr[ref_spec_locations[i]];
276 
277  // refraction contribution (excluding the constant factor p/p_0 * T0/T):
278  n += ref_n[i] * rtp_vmr[ref_spec_locations[i]];
279  }
280  }
281 
282  /*
283  if ( abs(ref_spec_vmr_sum-1) > 0.1 )
284  {
285  ostringstream os;
286  os << "Error: The total VMR of all your defined refractive\n"
287  << "species is " << ref_spec_vmr_sum
288  << ", more than 10% " << "different from 1.\n";
289  throw runtime_error(os.str());
290  }
291  */
292 
293  // normalize refractive index with the considered total VMR:
294  if (ref_spec_vmr_sum != 0) n /= ref_spec_vmr_sum;
295 
296  // now applying the constant factor p/p_0 * T0/T:
297  n *= (ratioT * ratiop);
298 
299  refr_index_air += n;
301 }
302 
303 /*===========================================================================
304  === WSMs for complex_refr_index
305  ===========================================================================*/
306 
307 /* Workspace method: Doxygen documentation will be auto-generated */
309  const Numeric& refr_index_real,
310  const Numeric& refr_index_imag,
311  const Verbosity&) {
312  complex_refr_index.resize(1, 1, 2);
313  complex_refr_index.set_grid_name(0, "Frequency");
314  complex_refr_index.set_grid(0, Vector(1, 0));
315  complex_refr_index.set_grid_name(1, "Temperature");
316  complex_refr_index.set_grid(1, Vector(1, 0));
317  complex_refr_index.set_grid_name(2, "Complex");
318  complex_refr_index.set_grid(2, {"real", "imaginary"});
319 
320  complex_refr_index.data(joker, joker, 0) = refr_index_real;
321  complex_refr_index.data(joker, joker, 1) = refr_index_imag;
322 }
323 
324 /* Workspace method: Doxygen documentation will be auto-generated */
326  const Vector& f_grid,
327  const Vector& t_grid,
328  const Verbosity& verbosity) {
329  if (min(t_grid) < 250) {
330  CREATE_OUT1;
331  out1 << "WARNING! The minimum chosen temperature is " << min(t_grid)
332  << ". Temperatures below 250 K may lead to incorrect values of"
333  " *complex_refr_index*.\n";
334  }
335 
336  const Index nf = f_grid.nelem();
337  const Index nt = t_grid.nelem();
338 
339  complex_refr_index.resize(nf, nt, 2);
340  complex_refr_index.set_grid_name(0, "Frequency");
341  complex_refr_index.set_grid(0, f_grid);
342  complex_refr_index.set_grid_name(1, "Temperature");
343  complex_refr_index.set_grid(1, t_grid);
344  complex_refr_index.set_grid_name(2, "Complex");
345  complex_refr_index.set_grid(2, {"real", "imaginary"});
346 
347  Matrix complex_n;
348  for (Index t = 0; t < nt; ++t) {
349  complex_n_water_liebe93(complex_n, f_grid, t_grid[t]);
350  complex_refr_index.data(joker, t, joker) = complex_n;
351  }
352 }
353 
354 /* Workspace method: Doxygen documentation will be auto-generated */
356  const Vector& f_grid,
357  const Vector& t_grid,
358  const Verbosity&) {
359  const Index nf = f_grid.nelem();
360  const Index nt = t_grid.nelem();
361 
362  // Frequency must be between 10MHz and 3THz
363  const Numeric f_min = 10e6;
364  const Numeric f_max = 3e12;
366  "min of complex_refr_index f_grid", min(f_grid), f_min, f_max);
368  "max of complex_refr_index f_grid", max(f_grid), f_min, f_max);
369 
370  // Temperature must be between 213.16 to 272.16 K
371  const Numeric t_min = 20.;
372  const Numeric t_max = 280.;
374  "min of complex_refr_index t_grid", min(t_grid), t_min, t_max);
376  "max of complex_refr_index t_grid", max(t_grid), t_min, t_max);
377 
378  complex_refr_index.resize(nf, nt, 2);
379  complex_refr_index.set_grid_name(0, "Frequency");
380  complex_refr_index.set_grid(0, f_grid);
381  complex_refr_index.set_grid_name(1, "Temperature");
382  complex_refr_index.set_grid(1, t_grid);
383  complex_refr_index.set_grid_name(2, "Complex");
384  complex_refr_index.set_grid(2, {"real", "imaginary"});
385 
386  Matrix complex_n;
387  for (Index i_t = 0; i_t < nt; ++i_t) {
388  complex_n_ice_matzler06(complex_n, f_grid, t_grid[i_t]);
389  complex_refr_index.data(joker, i_t, joker) = complex_n;
390  }
391 }
392 
393 #ifdef ENABLE_REFICE
394 /* Workspace method: Doxygen documentation will be auto-generated */
396  const Vector& f_grid,
397  const Vector& t_grid,
398  const Verbosity&) {
399  extern const Numeric SPEED_OF_LIGHT;
400  const Index nf = f_grid.nelem();
401  const Index nt = t_grid.nelem();
402 
403  // Frequency must be between 0.0443 to 8.600E+06 microns
404  const Numeric f_max = 1e6 * SPEED_OF_LIGHT / 0.0443;
405  const Numeric f_min = 1e6 * SPEED_OF_LIGHT / 8.6e6;
406  chk_if_in_range("min of scat_f_grid", min(f_grid), f_min, f_max);
407  chk_if_in_range("max of scat_f_grid", max(f_grid), f_min, f_max);
408 
409  // Temperature must be between 213.16 to 272.16 K
410  const Numeric t_min = 213.16;
411  const Numeric t_max = 272.16;
412  chk_if_in_range("min of scat_t_grid", min(t_grid), t_min, t_max);
413  chk_if_in_range("max of scat_t_grid", max(t_grid), t_min, t_max);
414 
415  complex_refr_index.resize(nf, nt, 2);
416  complex_refr_index.set_grid_name(0, "Frequency");
417  complex_refr_index.set_grid(0, f_grid);
418  complex_refr_index.set_grid_name(1, "Temperature");
419  complex_refr_index.set_grid(1, t_grid);
420  complex_refr_index.set_grid_name(2, "Complex");
421  complex_refr_index.set_grid(2, {"real", "imaginary"});
422 
423  Complex n;
424  /*#pragma omp parallel for \
425  if (!arts_omp_in_parallel() && nf > 1) \
426  private(n)*/
427  for (Index f = 0; f < nf; ++f)
428  for (Index t = 0; t < nt; ++t) {
429  n = refice_(1e6 * SPEED_OF_LIGHT / f_grid[f], t_grid[t]);
430  complex_refr_index.data(f, t, 0) = n.real();
431  complex_refr_index.data(f, t, 1) = n.imag();
432  }
433 }
434 
435 #else
436 
437 /* Workspace method: Doxygen documentation will be auto-generated */
438 void complex_refr_indexIceWarren84( // Generic output:
439  GriddedField3&,
440  // Generic input:
441  const Vector&,
442  const Vector&,
443  const Verbosity&) {
444  throw std::runtime_error("ARTS was not compiled with Fortran support.");
445 }
446 
447 #endif
Matrix
The Matrix class.
Definition: matpackI.h:1193
refr_index_airMicrowavesGeneral
void refr_index_airMicrowavesGeneral(Numeric &refr_index_air, Numeric &refr_index_air_group, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Vector &rtp_vmr, const ArrayOfArrayOfSpeciesTag &abs_species, const Verbosity &)
WORKSPACE METHOD: refr_index_airMicrowavesGeneral.
Definition: m_refraction.cc:174
absorption.h
Declarations required for the calculation of absorption coefficients.
refr_index_airInfraredEarth
void refr_index_airInfraredEarth(Numeric &refr_index_air, Numeric &refr_index_air_group, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Verbosity &)
WORKSPACE METHOD: refr_index_airInfraredEarth.
Definition: m_refraction.cc:120
joker
const Joker joker
refr_index_airFreeElectrons
void refr_index_airFreeElectrons(Numeric &refr_index_air, Numeric &refr_index_air_group, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const Vector &rtp_vmr, const Index &demand_vmr_value, const Verbosity &)
WORKSPACE METHOD: refr_index_airFreeElectrons.
Definition: m_refraction.cc:60
last
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:165
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
VACUUM_PERMITTIVITY
const Numeric VACUUM_PERMITTIVITY
ARTS::Var::refr_index_air
Numeric refr_index_air(Workspace &ws) noexcept
Definition: autoarts.h:5532
complex_refr_indexWaterLiebe93
void complex_refr_indexWaterLiebe93(GriddedField3 &complex_refr_index, const Vector &f_grid, const Vector &t_grid, const Verbosity &verbosity)
WORKSPACE METHOD: complex_refr_indexWaterLiebe93.
Definition: m_refraction.cc:325
ARTS::Var::rtp_vmr
Vector rtp_vmr(Workspace &ws) noexcept
Definition: autoarts.h:5826
PI
const Numeric PI
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
complex_refr_indexIceMatzler06
void complex_refr_indexIceMatzler06(GriddedField3 &complex_refr_index, const Vector &f_grid, const Vector &t_grid, const Verbosity &)
WORKSPACE METHOD: complex_refr_indexIceMatzler06.
Definition: m_refraction.cc:355
ELECTRON_MASS
const Numeric ELECTRON_MASS
CREATE_OUT1
#define CREATE_OUT1
Definition: messages.h:205
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
GriddedField3
Definition: gridded_fields.h:284
ARTS::Var::complex_refr_index
GriddedField3 complex_refr_index(Workspace &ws) noexcept
Definition: autoarts.h:2829
matpackI.h
Implementation of Matrix, Vector, and such stuff.
Array
This can be used to make arrays out of anything.
Definition: array.h:108
species_index_from_species_name
Index species_index_from_species_name(String name)
Return species index for given species name.
Definition: absorption.cc:531
min
#define min(a, b)
Definition: legacy_continua.cc:20628
messages.h
Declarations having to do with the four output streams.
complex_refr_indexIceWarren84
void complex_refr_indexIceWarren84(GriddedField3 &, const Vector &, const Vector &, const Verbosity &)
WORKSPACE METHOD: complex_refr_indexIceWarren84.
Definition: m_refraction.cc:438
ELECTRON_CHARGE
const Numeric ELECTRON_CHARGE
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
find_first_species_tg
Index find_first_species_tg(const ArrayOfArrayOfSpeciesTag &tgs, const Index &spec)
Find first occurrence of species in tag groups.
Definition: abs_species_tags.cc:532
physics_funcs.h
This file contains declerations of functions of physical character.
ARTS::Var::abs_species
ArrayOfArrayOfSpeciesTag abs_species(Workspace &ws) noexcept
Definition: autoarts.h:2157
SpeciesTag::TYPE_FREE_ELECTRONS
@ TYPE_FREE_ELECTRONS
Definition: abs_species_tags.h:152
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:49
ARTS::Var::refr_index_air_group
Numeric refr_index_air_group(Workspace &ws) noexcept
Definition: autoarts.h:5553
math_funcs.h
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
max
#define max(a, b)
Definition: legacy_continua.cc:20629
abs_species_tags.h
Header file for stuff related to absorption species tags.
LineShape::Type
Type
Type of line shape to compute.
Definition: lineshapemodel.h:788
ARTS::Var::rtp_pressure
Numeric rtp_pressure(Workspace &ws) noexcept
Definition: autoarts.h:5790
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.
refr_index_airMicrowavesEarth
void refr_index_airMicrowavesEarth(Numeric &refr_index_air, Numeric &refr_index_air_group, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Vector &rtp_vmr, const ArrayOfArrayOfSpeciesTag &abs_species, const Numeric &k1, const Numeric &k2, const Numeric &k3, const Verbosity &)
WORKSPACE METHOD: refr_index_airMicrowavesEarth.
Definition: m_refraction.cc:139
special_interp.h
Header file for special_interp.cc.
SPEED_OF_LIGHT
const Numeric SPEED_OF_LIGHT
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
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
check_input.h
Vector
The Vector class.
Definition: matpackI.h:860
complex_refr_indexConstant
void complex_refr_indexConstant(GriddedField3 &complex_refr_index, const Numeric &refr_index_real, const Numeric &refr_index_imag, const Verbosity &)
WORKSPACE METHOD: complex_refr_indexConstant.
Definition: m_refraction.cc:308
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
refice_
Complex refice_(const Numeric &wavlen, const Numeric &temp)
Calculates complex refractive index of Ice 1H.
TORR2PA
const Numeric TORR2PA
arts.h
The global header file for ARTS.