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