ARTS  2.2.66
gas_abs_lookup.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.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 
26 #include <cmath>
27 #include "gas_abs_lookup.h"
28 #include "interpolation.h"
29 #include "interpolation_poly.h"
30 #include "make_vector.h"
31 #include "logic.h"
32 #include "check_input.h"
33 #include "messages.h"
34 #include "physics_funcs.h"
35 
37 
48  ConstVectorView old_grid,
49  ConstVectorView new_grid,
50  const Verbosity& verbosity )
51 {
53 
54  const Index n_new_grid = new_grid.nelem();
55  const Index n_old_grid = old_grid.nelem();
56 
57  // Make sure that pos has the right size:
58  assert( n_new_grid == pos.nelem() );
59 
60  // The frequency difference in Hz that we are willing to tolerate. 1
61  // Hz seems to be on the safe side.
62  const Numeric tolerance = 1;
63 
64  // Old grid position:
65  Index j = 0;
66 
67  // Loop the new frequencies:
68  for ( Index i=0; i<n_new_grid; ++i )
69  {
70  // We have done runtime checks that both the new and the old
71  // frequency gris are sorted in GasAbsLookup::Adapt, so we can
72  // use the fact here.
73 
74  while ( abs(new_grid[i]-old_grid[j]) > tolerance )
75  {
76  ++j;
77  if ( j>=n_old_grid )
78  {
79  ostringstream os;
80  os << "Cannot find new frequency " << i
81  << " (" << new_grid[i] << "Hz) in the lookup table frequency grid.";
82  throw runtime_error( os.str() );
83  }
84  }
85 
86  pos[i] = j;
87  out3 << " " << new_grid[i] << " found, index = " << pos[i] << ".\n";
88  }
89 }
90 
92 
124 void GasAbsLookup::Adapt( const ArrayOfArrayOfSpeciesTag& current_species,
125  ConstVectorView current_f_grid,
126  const Verbosity& verbosity )
127 {
128  CREATE_OUT2;
129  CREATE_OUT3;
130 
131  // Some constants we will need:
132  const Index n_current_species = current_species.nelem();
133  const Index n_current_f_grid = current_f_grid.nelem();
134 
135  const Index n_species = species.nelem();
136  const Index n_nls = nonlinear_species.nelem();
137  const Index n_nls_pert = nls_pert.nelem();
138  const Index n_f_grid = f_grid.nelem();
139  const Index n_p_grid = p_grid.nelem();
140 
141  out2 << " Original table: " << n_species << " species, "
142  << n_f_grid << " frequencies.\n"
143  << " Adapt to: " << n_current_species << " species, "
144  << n_current_f_grid << " frequencies.\n";
145 
146  if ( 0 == n_nls )
147  {
148  out2 << " Table contains no nonlinear species.\n";
149  }
150 
151  // Set up a logical array for the nonlinear species
152  ArrayOfIndex non_linear(n_species,0);
153  for ( Index s=0; s<n_nls; ++s )
154  {
155  non_linear[nonlinear_species[s]] = 1;
156  }
157 
158  if ( 0 == t_pert.nelem() )
159  {
160  out2 << " Table contains no temperature perturbations.\n";
161  }
162 
163  // We are constructing a new lookup table, containing just the
164  // species and frequencies that are necessary for the current
165  // calculation. We will build it in this local variable, then copy
166  // it back to *this in the end.
167  GasAbsLookup new_table;
168 
169  // First some checks on the lookup table itself:
170 
171  // Species:
172  if ( 0 == n_species )
173  {
174  ostringstream os;
175  os << "The lookup table should have at least one species.";
176  throw runtime_error( os.str() );
177  }
178 
179  // Nonlinear species:
180  // They should be unique ...
182  {
183  ostringstream os;
184  os << "The table must not have duplicate nonlinear species.\n"
185  << "Value of *nonlinear_species*: " << nonlinear_species;
186  throw runtime_error( os.str() );
187  }
188 
189  // ... and pointing at valid species.
190  for ( Index i=0; i<n_nls; ++i )
191  {
192  ostringstream os;
193  os << "nonlinear_species[" << i << "]";
194  chk_if_in_range( os.str(),
196  0,
197  n_species-1 );
198  }
199 
200  // Frequency grid:
201  chk_if_increasing( "f_grid", f_grid );
202 
203  // Pressure grid:
204  chk_if_decreasing( "p_grid", p_grid );
205 
206  // Reference VMRs:
207  chk_matrix_nrows( "vmrs_ref", vmrs_ref, n_species );
208  chk_matrix_ncols( "vmrs_ref", vmrs_ref, n_p_grid );
209 
210  // Reference temperatur:
211  chk_vector_length( "t_ref", t_ref, n_p_grid );
212 
213  // Temperature perturbations:
214  // Nothing to check for t_pert, it seems.
215 
216  // Perturbations for nonlinear species:
217  // Check that nls_pert is empty if and only if nonlinear_species is
218  // empty:
219  if ( 0 == n_nls )
220  {
221  chk_vector_length( "nls_pert", nls_pert, 0 );
222  }
223  else
224  {
225  if ( 0 == n_nls_pert )
226  {
227  ostringstream os;
228  os << "The vector nls_pert should contain the perturbations\n"
229  << "for the nonlinear species, but it is empty.";
230  throw runtime_error( os.str() );
231  }
232  }
233 
234  // The table itself, xsec:
235  //
236  // We have to separtely consider the 3 cases described in the
237  // documentation of GasAbsLookup.
238  //
239  // Dimension: [ a, b, c, d ]
240  //
241  if ( 0==n_nls )
242  {
243  if ( 0==t_pert.nelem() )
244  {
245  // Simplest case (no temperature perturbations,
246  // no vmr perturbations):
247  // a = 1
248  // b = n_species
249  // c = n_f_grid
250  // d = n_p_grid
251  chk_size( "xsec", xsec,
252  1,
253  n_species,
254  n_f_grid,
255  n_p_grid );
256  }
257  else
258  {
259  // Standard case (temperature perturbations,
260  // but no vmr perturbations):
261  // a = n_t_pert
262  // b = n_species
263  // c = n_f_grid
264  // d = n_p_grid
265  chk_size( "xsec", xsec,
266  t_pert.nelem(),
267  n_species,
268  n_f_grid,
269  n_p_grid );
270  }
271  }
272  else
273  {
274  // Full case (with temperature perturbations and
275  // vmr perturbations):
276  // a = n_t_pert
277  // b = n_species + n_nonlinear_species * ( n_nls_pert - 1 )
278  // c = n_f_grid
279  // d = n_p_grid
280  Index a = t_pert.nelem();
281  Index b = n_species
282  + n_nls * ( n_nls_pert - 1 );
283  Index c = n_f_grid;
284  Index d = n_p_grid;
285 
286  chk_size( "xsec", xsec, a, b, c, d );
287  }
288 
289  // We also need indices to the positions of the original species
290  // data in xsec. Nonlinear species take more space, therefor the
291  // position in xsec is not the same as the position in species.
292  ArrayOfIndex original_spec_pos_in_xsec(n_species);
293  for (Index i=0,sp=0; i<n_species; ++i)
294  {
295  original_spec_pos_in_xsec[i] = sp;
296  if (non_linear[i]) sp += n_nls_pert;
297  else sp += 1;
298  }
299 
300 
301 
302  // Now some checks on the input data:
303 
304  // The list of current species should not be empty:
305  if ( 0==n_current_species )
306  {
307  ostringstream os;
308  os << "The list of current species should not be empty.";
309  throw runtime_error( os.str() );
310  }
311 
312  // The grid of current frequencies should be monotonically sorted:
313  chk_if_increasing( "current_f_grid", current_f_grid );
314 
315 
316  // 1. Find and remember the indices of the current species in the
317  // lookup table. At the same time verify that each species is
318  // included in the table exactly once.
319  ArrayOfIndex i_current_species(n_current_species);
320  out3 << " Looking for species in lookup table:\n";
321  for ( Index i=0; i<n_current_species; ++i )
322  {
323  out3 << " " << get_tag_group_name( current_species[i] ) << ": ";
324 
325  try {
326  i_current_species[i] =
327  chk_contains( "abs_species", species, current_species[i] );
328 
329  } catch (runtime_error_not_found) {
330 
331  // Is this one of the trivial species?
332  const Index spec_type = current_species[i][0].Type();
333  if (spec_type == SpeciesTag::TYPE_ZEEMAN ||
334  spec_type == SpeciesTag::TYPE_FREE_ELECTRONS ||
335  spec_type == SpeciesTag::TYPE_PARTICLES)
336  {
337  // Set to -1 to flag that this needs special handling later on.
338  i_current_species[i] = -1;
339  }
340  else
341  {
342  ostringstream os;
343  os << "Species " << get_tag_group_name( current_species[i] )
344  << " is missing in absorption lookup table.";
345  throw runtime_error( os.str() );
346  }
347  }
348 
349  out3 << "found (or trivial).\n";
350  }
351 
352  // 1a. Find out which of the current species are nonlinear species:
353  Index n_current_nonlinear_species = 0; // Number of current nonlinear species
354  ArrayOfIndex current_non_linear(n_current_species,0); // A logical array to
355  // flag which of the
356  // current species are
357  // nonlinear.
358 
359  out3 << " Finding out which of the current species are nonlinear:\n";
360  for ( Index i=0; i<n_current_species; ++i )
361  {
362  if (i_current_species[i]>=0) // Jump over trivial species here.
363  {
364  // Check if this is a nonlinear species:
365  if (non_linear[i_current_species[i]])
366  {
367  out3 << " " << get_tag_group_name( current_species[i] ) << "\n";
368 
369  current_non_linear[i] = 1;
370  ++n_current_nonlinear_species;
371  }
372  }
373  }
374 
375  // 2. Find and remember the frequencies of the current calculation in
376  // the lookup table. At the same time verify that all frequencies are
377  // included and that no frequency occurs twice.
378 
379  // FIXME: This is a bit tricky, because we are comparing
380  // Numerics. Let's see how well this works in practice.
381 
382  ArrayOfIndex i_current_f_grid(n_current_f_grid);
383  out3 << " Looking for Frequencies in lookup table:\n";
384 
385  // We need no error checking for the next statement, since the
386  // function called throws a runtime error if a frequency
387  // is not found, or if the grids are not ok.
388  find_new_grid_in_old_grid( i_current_f_grid,
389  f_grid,
390  current_f_grid,
391  verbosity );
392 
393 
394  // 3. Use the species and frequency index lists to build the new lookup
395  // table.
396 
397  // Species:
398  new_table.species.resize( n_current_species );
399  for ( Index i=0; i<n_current_species; ++i )
400  {
401  if (i_current_species[i]>=0)
402  {
403  new_table.species[i] = species[i_current_species[i]];
404 
405  // Is this a nonlinear species?
406  if (current_non_linear[i])
407  new_table.nonlinear_species.push_back( i );
408  }
409  else
410  {
411  // Here we handle the case of the trivial species, for which we simply
412  // copy the name:
413 
414  new_table.species[i] = current_species[i];
415  }
416  }
417 
418  // Frequency grid:
419  new_table.f_grid.resize( n_current_f_grid );
420  for ( Index i=0; i<n_current_f_grid; ++i )
421  {
422  new_table.f_grid[i] = f_grid[i_current_f_grid[i]];
423  }
424 
425  // Pressure grid:
426  // new_table.p_grid.resize( n_p_grid );
427  new_table.p_grid = p_grid;
428 
429  // Reference VMR profiles:
430  new_table.vmrs_ref.resize( n_current_species,
431  n_p_grid );
432  for ( Index i=0; i<n_current_species; ++i )
433  {
434  if (i_current_species[i]>=0)
435  {
436  new_table.vmrs_ref(i, Range(joker))
437  = vmrs_ref(i_current_species[i], Range(joker));
438  }
439  else
440  {
441  // Here we handle the case of the trivial species, for which we set
442  // the reference VMR to NAN.
443  new_table.vmrs_ref(i, Range(joker)) = NAN;
444  }
445  }
446 
447  // Reference temperature profile:
448  // new_table.t_ref.resize( t_ref.nelem() );
449  new_table.t_ref = t_ref;
450 
451  // Vector of temperature perturbations:
452  // new_table.t_pert.resize( t_pert.nelem() );
453  new_table.t_pert = t_pert;
454 
455  // Vector of perturbations for the VMRs of the nonlinear species:
456  // (Should stay empty if we have no nonlinear species)
457  if ( 0 != new_table.nonlinear_species.nelem() )
458  {
459  // new_table.nls_pert.resize( n_nls_pert );
460  new_table.nls_pert = nls_pert;
461  }
462 
463  // Absorption coefficients:
464  new_table.xsec.resize( xsec.nbooks(),
465  n_current_species+n_current_nonlinear_species*(n_nls_pert-1),
466  n_current_f_grid,
467  xsec.ncols() );
468 
469  // We have to copy the right species and frequencies from the old to
470  // the new table. Temperature perturbations and pressure grid remain
471  // the same.
472 
473  // Do species:
474  for ( Index i_s=0,sp=0; i_s<n_current_species; ++i_s )
475  {
476  // n_v is the number of VMR perturbations
477  Index n_v;
478  if (current_non_linear[i_s])
479  n_v = n_nls_pert;
480  else
481  n_v = 1;
482 
483 // cout << "i_s / sp / n_v = " << i_s << " / " << sp << " / " << n_v << endl;
484 // cout << "orig_pos = " << original_spec_pos_in_xsec[i_current_species[i_s]] << endl;
485 
486  // Do frequencies:
487  for ( Index i_f=0; i_f<n_current_f_grid; ++i_f )
488  {
489  if (i_current_species[i_s]>=0)
490  {
491  new_table.xsec(Range(joker),
492  Range(sp,n_v),
493  i_f,
494  Range(joker) )
495  =
496  xsec(Range(joker),
497  Range(original_spec_pos_in_xsec[i_current_species[i_s]],n_v),
498  i_current_f_grid[i_f],
499  Range(joker) );
500  }
501  else
502  {
503  // Here we handle the case of the trivial species, which we simply
504  // set to NAN:
505  new_table.xsec(Range(joker),
506  Range(sp,n_v),
507  i_f,
508  Range(joker) ) = NAN;
509 
510  }
511 
512 
513 // cout << "result: " << xsec( Range(joker),
514 // Range(original_spec_pos_in_xsec[i_current_species[i_s]],n_v),
515 // i_current_f_grid[i_f],
516 // Range(joker) ) << endl;
517 
518  }
519 
520  sp += n_v;
521  }
522 
523 
524  // 4. Replace original table by the new one.
525  *this = new_table;
526 
527  // 5. Initialize log_p_grid.
528  log_p_grid.resize( n_p_grid );
530  log,
531  p_grid );
532 
533  // 6. Initialize fgp_default.
534  fgp_default.resize(f_grid.nelem());
536 
537 }
538 
540 
594  const Index& p_interp_order,
595  const Index& t_interp_order,
596  const Index& h2o_interp_order,
597  const Index& f_interp_order,
598  const Numeric& p,
599  const Numeric& T,
600  ConstVectorView abs_vmrs,
601  ConstVectorView new_f_grid,
602  const Numeric& extpolfac) const
603 {
604  // 1. Obtain some properties of the lookup table:
605 
606  // Number of gas species in the table:
607  const Index n_species = species.nelem();
608 
609  // Number of nonlinear species:
610  const Index n_nls = nonlinear_species.nelem();
611 
612  // Number of frequencies in the table:
613  const Index n_f_grid = f_grid.nelem();
614 
615  // Number of pressure grid points in the table:
616  const Index n_p_grid = p_grid.nelem();
617 
618  // Number of temperature perturbations:
619  const Index n_t_pert = t_pert.nelem();
620 
621  // Number of nonlinear species perturbations:
622  const Index n_nls_pert = nls_pert.nelem();
623 
624  // Number of frequencies in new_f_grid, the frequency grid for which we
625  // want to extract.
626  const Index n_new_f_grid = new_f_grid.nelem();
627 
628  // 2. First some checks on the lookup table itself:
629 
630  // Most checks here are asserts, because they check the internal
631  // consistency of the lookup table. They should never fail if the
632  // table has been created with ARTS.
633 
634  // If there are nonlinear species, then at least one species must be
635  // H2O. We will use that to perturb in the case of nonlinear
636  // species.
637  Index h2o_index = -1;
638  if (n_nls>0)
639  {
640  h2o_index = find_first_species_tg( species,
642 
643  // This is a runtime error, even though it would be more logical
644  // for it to be an assertion, since it is an internal check on
645  // the table. The reason is that it is somewhat awkward to check
646  // for this in other places.
647  if ( h2o_index == -1 )
648  {
649  ostringstream os;
650  os << "With nonlinear species, at least one species must be a H2O species.";
651  throw runtime_error( os.str() );
652  }
653  }
654 
655  // Check that the dimension of vmrs_ref is consistent with species and p_grid:
656  assert( is_size( vmrs_ref, n_species, n_p_grid) );
657 
658  // Check dimension of t_ref:
659  assert( is_size( t_ref, n_p_grid) );
660 
661  // Check dimension of xsec:
662  DEBUG_ONLY({
663  Index a,b,c,d;
664  if ( 0 == n_t_pert ) a = 1;
665  else a = n_t_pert;
666  b = n_species + n_nls * ( n_nls_pert - 1 );
667  c = n_f_grid;
668  d = n_p_grid;
669 // cout << "xsec: "
670 // << xsec.nbooks() << ", "
671 // << xsec.npages() << ", "
672 // << xsec.nrows() << ", "
673 // << xsec.ncols() << "\n";
674 // cout << "a b c d: "
675 // << a << ", "
676 // << b << ", "
677 // << c << ", "
678 // << d << "\n";
679  assert( is_size( xsec, a, b, c, d ) );
680  })
681 
682  // Make sure that log_p_grid is initialized:
683  if (log_p_grid.nelem() != n_p_grid)
684  {
685  ostringstream os;
686  os << "The lookup table internal variable log_p_grid is not initialized.\n"
687  << "Use the abs_lookupAdapt method!";
688  throw runtime_error( os.str() );
689  }
690 
691  // Verify that we have enough pressure, temperature,humdity, and frequency grid points
692  // for the desired interpolation orders. This check is not only
693  // table internal, since abs_nls_interp_order and abs_t_interp_order
694  // are separate WSVs that could have been modified. Hence, these are
695  // runtime errors.
696 
697  if ( (n_p_grid < p_interp_order+1) )
698  {
699  ostringstream os;
700  os << "The number of pressure grid points in the table ("
701  << n_p_grid << ") is not enough for the desired order of interpolation ("
702  << p_interp_order << ").";
703  throw runtime_error( os.str() );
704  }
705 
706  if ( (n_nls_pert != 0) && (n_nls_pert < h2o_interp_order+1) )
707  {
708  ostringstream os;
709  os << "The number of humidity perturbation grid points in the table ("
710  << n_nls_pert << ") is not enough for the desired order of interpolation ("
711  << h2o_interp_order << ").";
712  throw runtime_error( os.str() );
713  }
714 
715  if ( (n_t_pert != 0) && (n_t_pert < t_interp_order+1) )
716  {
717  ostringstream os;
718  os << "The number of temperature perturbation grid points in the table ("
719  << n_t_pert << ") is not enough for the desired order of interpolation ("
720  << t_interp_order << ").";
721  throw runtime_error( os.str() );
722  }
723 
724  if ( (n_f_grid < f_interp_order+1) )
725  {
726  ostringstream os;
727  os << "The number of frequency grid points in the table ("
728  << n_f_grid << ") is not enough for the desired order of interpolation ("
729  << f_interp_order << ").";
730  throw runtime_error( os.str() );
731  }
732 
733 
734  // 3. Checks on the input variables:
735 
736  // Check that abs_vmrs has the right dimension:
737  if ( !is_size( abs_vmrs, n_species ) )
738  {
739  ostringstream os;
740  os << "Number of species in lookup table does not match number\n"
741  << "of species for which you want to extract absorption.\n"
742  << "Have you used abs_lookupAdapt? Or did you miss to add\n"
743  << "some VRM fields (e.g. for free electrons or particles)?\n";
744  throw runtime_error( os.str() );
745  }
746 
747 
748  // 4. Set up some things we will need later on:
749 
750  // 4.a Frequency grid positions
751 
752  // Frequency grid positions. The pointer is used to save copying of the
753  // default from the lookup table.
754  const ArrayOfGridPosPoly *fgp;
755  ArrayOfGridPosPoly fgp_local;
756 
757  // With f_interp_order 0 the frequency grid has to have the same size as in the
758  // lookup table, or exactly one element. If it matches the lookup table, we
759  // do no frequency interpolation at all. (We set the frequency grid positions
760  // to the predefined ones that come with the lookup table.)
761  if (f_interp_order==0)
762  {
763  if (n_new_f_grid==n_f_grid) {
764 
765  // Use the default fgp that is stored in the lookup table itself
766  // (which effectively means no frequency interpolation)
767  fgp = &fgp_default;
768 
769  // Check identitiy of first and last element for safety's sake (only
770  // till 1Hz identity; checking float for total identity is sometimes
771  // problematic)
772 
773  if (abs(f_grid[0]-new_f_grid[0]) > 1)
774  //if (f_grid[0]!=new_f_grid[0])
775  {
776  ostringstream os;
777  os << "First frequency in f_grid inconsistent with lookup table.\n"
778  << "f_grid[0] = " << f_grid[0] << "\n"
779  << "new_f_grid[0] = " << new_f_grid[0] << ".";
780  throw runtime_error( os.str() );
781  }
782 
783  if (abs(f_grid[n_f_grid-1]-new_f_grid[n_new_f_grid-1]) > 1)
784  //if (f_grid[n_f_grid-1]!=new_f_grid[n_new_f_grid-1])
785  {
786  ostringstream os;
787  os << "Last frequency in f_grid inconsistent with lookup table.\n"
788  << "f_grid[n_f_grid-1] = " << f_grid[n_f_grid-1] << "\n"
789  << "new_f_grid[n_new_f_grid-1] = " << new_f_grid[n_new_f_grid-1] << ".";
790  throw runtime_error( os.str() );
791  }
792  }
793  else if (n_new_f_grid==1) {
794  fgp = &fgp_local;
795  fgp_local.resize(1);
796  gridpos_poly( fgp_local, f_grid, new_f_grid, 0 );
797 
798  // Check that we really are on a frequency grid point, for safety's sake.
799  if (fgp_local[0].w[0]!=1)
800  {
801  ostringstream os;
802  os << "Cannot find a matching lookup table frequency for frequency\n"
803  << new_f_grid[0] << ".";
804  throw runtime_error( os.str() );
805  }
806  }
807  else {
808  ostringstream os;
809  os << "With f_interp_order 0 the frequency grid has to have the same\n"
810  << "size as in the lookup table, or exactly one element.";
811  throw runtime_error( os.str() );
812  }
813  }
814  else
815  {
816  // We do have real frequency interpolation (f_interp_order!=0).
817  fgp = &fgp_local;
818  fgp_local.resize(n_new_f_grid);
819  gridpos_poly( fgp_local, f_grid, new_f_grid, f_interp_order);
820  }
821 
822 
823  // 4.b Other stuff
824 
825  // Flag for temperature interpolation, if this is not 0 we want
826  // to do T interpolation:
827  const Index do_T = n_t_pert;
828 
829 
830  // Set up a logical array for the nonlinear species
831  ArrayOfIndex non_linear(n_species,0);
832  for ( Index s=0; s<n_nls; ++s )
833  {
834  non_linear[nonlinear_species[s]] = 1;
835  }
836 
837 
838  // Calculate the number density for the given pressure and
839  // temperature:
840  // n = n0*T0/p0 * p/T or n = p/kB/t, ideal gas law
841  const Numeric n = number_density( p, T );
842 
843 
844  // 5. Determine pressure grid position and interpolation weights:
845 
846  // Check that p is inside the grid. (p_grid is sorted in decreasing order.)
847  {
848  const Numeric p_max = p_grid[0] + 0.5*(p_grid[0]-p_grid[1]);
849  const Numeric p_min = p_grid[n_p_grid-1] - 0.5*(p_grid[n_p_grid-2]-p_grid[n_p_grid-1]);
850  if ( ( p > p_max ) ||
851  ( p < p_min ) )
852  {
853  ostringstream os;
854  os << "Problem with gas absorption lookup table.\n"
855  << "Pressure p is outside the range covered by the lookup table.\n"
856  << "Your p value is " << p << " Pa.\n"
857  << "The allowed range is " << p_min << " to " << p_max << ".\n"
858  << "The pressure grid range in the table is " << p_grid[n_p_grid-1]
859  << " to " << p_grid[0] << ".\n"
860  << "We allow a bit of extrapolation, but NOT SO MUCH!";
861  throw runtime_error( os.str() );
862  }
863  }
864 
865  // For sure, we need to store the pressure grid position.
866  // We do the interpolation in log(p). Test have shown that this
867  // gives slightly better accuracy than interpolating in p directly.
868  ArrayOfGridPosPoly pgp(1);
869  gridpos_poly( pgp,
870  log_p_grid,
871  log(p),
872  p_interp_order,
873  extpolfac );
874 
875  // Pressure interpolation weights:
876  Vector pitw;
877  pitw.resize(p_interp_order+1);
878  interpweights(pitw,pgp[0]);
879 
880 
881  // Define also other grid positions and interpolation weights here, so that
882  // we do not have to allocate them over and over in the loops below.
883 
884  // Define the ArrayOfGridPosPoly that corresponds to "no interpolation at all".
885  ArrayOfGridPosPoly gp_trivial(1);
886  gp_trivial[0].idx.resize(1);
887  gp_trivial[0].w.resize(1);
888  gp_trivial[0].idx[0] = 0;
889  gp_trivial[0].w[0] = 1;
890 
891  // Temperature grid positions.
892  ArrayOfGridPosPoly tgp_withT(1); // Only a scalar.
893  ArrayOfGridPosPoly *tgp; // Pointer to either tgp_withT or gp_trivial.
894 
895  // Set this_t_interp_order, depending on whether we do T interpolation or not.
896  Index this_t_interp_order; // Local T interpolation order
897  if (do_T)
898  {
899  this_t_interp_order = t_interp_order;
900  tgp = &tgp_withT;
901  }
902  else
903  {
904  // For the !do_T case we simply take the single
905  // temperature that is there, so we point tgp accordingly.
906  this_t_interp_order = 0;
907  tgp = &gp_trivial;
908  }
909 
910 
911  // H2O(VMR) grid positions. vgp is what will be used in the interpolation.
912  // Depending on species, it is either pointed to gp_trivial, or to vgp_h2o.
913  ArrayOfGridPosPoly *vgp;
914  ArrayOfGridPosPoly vgp_h2o(1); // only a scalar
915 
916 
917  // 6. We do the T and VMR interpolation for the pressure levels
918  // that are used in the pressure interpolation. (How many depends on
919  // p_interp_order.)
920 
921  // To store the interpolated result for the p_interp_order+1
922  // pressure levels:
923  // xsec dimensions are:
924  // Temperature
925  // H2O
926  // Frequency
927  // Pressure
928  // Dimensions of pre_interpolated are:
929  // Pressure (interpolation points)
930  // Species
931  // Temperature (always 1)
932  // H2O (always 1)
933  // Frequency
934 
935  Tensor5 xsec_pre_interpolated;
936  xsec_pre_interpolated.resize(p_interp_order+1, n_species,
937  1, 1, n_new_f_grid );
938 
939  // Define variables for interpolation weights outside the loops.
940  // We will make itw point to either the weights with H2O interpolation, or
941  // the ones without.
942  Tensor4 itw_withH2O, itw_noH2O, *itw;
943 
944  for ( Index pi=0; pi<p_interp_order+1; ++pi )
945  {
946  // Throw a runtime error if one of the reference VMR profiles is zero, but
947  // abs_vmrs is not. (This means that the lookup table was calculated with a
948  // reference profile of zero for that gas.)
949 // for (Index si=0; si<n_species; ++si)
950 // if ( (vmrs_ref(si,pi) == 0) &&
951 // (abs_vmrs[si] != 0) )
952 // {
953 // ostringstream os;
954 // os << "Reference VMR profile is zero, you cannot extract\n"
955 // << "Absorption for this species.\n"
956 // << "Species: " << si
957 // << " (" << get_species_name(species[si]) << ")\n"
958 // << "Lookup table pressure level: " << pi
959 // << " (" << p_grid[pi] << " Pa).";
960 // throw runtime_error( os.str() );
961 // }
962 
963  // Index into p_grid:
964  const Index this_p_grid_index = pgp[0].idx[pi];
965 
966  // Determine temperature grid position. This is only done if we
967  // want temperature interpolation, but the variable tgp has to
968  // be visible also outside for later use:
969  if (do_T)
970  {
971 
972  // Temperature in the atmosphere is altitude
973  // dependent. When we do the interpolation for the pressure level
974  // below and above our point, we should correct the target value of
975  // the interpolation to the altitude (pressure) difference. This
976  // ensures that there is for example no T interpolation if the
977  // desired T is right on the reference profile curve.
978  //
979  // I explicitly compared this with the old option to calculate
980  // the temperature offset relative to the temperature at
981  // this level. The performance in both cases is very
982  // similar. The reason, why I decided to keep this new
983  // version, is that it avoids the problem of needing
984  // oversized temperature perturbations if the pressure
985  // grid is coarse.
986  //
987  // No! The above approach leads to problems when combined with
988  // higher order pressure interpolation. The problem is that
989  // the reference T and VMR profiles may be very
990  // irregular. (For example the H2O profile often has a big
991  // jump near the bottom.) That sometimes leads to negative
992  // effective reference values when the reference profile is
993  // interpolated. I therefore reverted back to the original
994  // version of using the real temperature and humidity, not
995  // the interpolated one.
996 
997  // const Numeric effective_T_ref = interp(pitw,t_ref,pgp);
998  const Numeric effective_T_ref = t_ref[this_p_grid_index];
999 
1000  // Convert temperature to offset from t_ref:
1001  const Numeric T_offset = T - effective_T_ref;
1002 
1003  // cout << "T_offset = " << T_offset << endl;
1004 
1005  // Check that temperature offset is inside the allowed range.
1006  {
1007  const Numeric t_min = t_pert[0] - 0.5*(t_pert[1]-t_pert[0]);
1008  const Numeric t_max = t_pert[n_t_pert-1] + 0.5*(t_pert[n_t_pert-1]-t_pert[n_t_pert-2]);
1009  if ( ( T_offset > t_max ) ||
1010  ( T_offset < t_min ) )
1011  {
1012  ostringstream os;
1013  os << "Problem with gas absorption lookup table.\n"
1014  << "Temperature T is outside the range covered by the lookup table.\n"
1015  << "Your temperature was " << T
1016  << " K at a pressure of " << p << " Pa.\n"
1017  << "The temperature offset value is " << T_offset << ".\n"
1018  << "The allowed range is " << t_min << " to " << t_max << ".\n"
1019  << "The temperature perturbation grid range in the table is "
1020  << t_pert[0] << " to " << t_pert[n_t_pert-1] << ".\n"
1021  << "We allow a bit of extrapolation, but NOT SO MUCH!";
1022  throw runtime_error( os.str() );
1023  }
1024  }
1025 
1026  gridpos_poly( tgp_withT, t_pert, T_offset, t_interp_order, extpolfac );
1027  }
1028 
1029  // Determine the H2O VMR grid position. We need to do this only
1030  // once, since the only species who's VMR is interpolated is
1031  // H2O. We do this only if there are nonlinear species, but the
1032  // variable has to be visible later.
1033  if (n_nls>0)
1034  {
1035  // Similar to the T case, we first interpolate the reference
1036  // VMR to the pressure of extraction, then compare with
1037  // the extraction VMR to determine the offset/fractional
1038  // difference for the VMR interpolation.
1039  //
1040  // No! The above approach leads to problems when combined with
1041  // higher order pressure interpolation. The problem is that
1042  // the reference T and VMR profiles may be very
1043  // irregular. (For example the H2O profile often has a big
1044  // jump near the bottom.) That sometimes leads to negative
1045  // effective reference values when the reference profile is
1046  // interpolated. I therefore reverted back to the original
1047  // version of using the real temperature and humidity, not
1048  // the interpolated one.
1049 
1050 // const Numeric effective_vmr_ref = interp(pitw,
1051 // vmrs_ref(h2o_index, Range(joker)),
1052 // pgp);
1053  const Numeric effective_vmr_ref = vmrs_ref(h2o_index, this_p_grid_index);
1054 
1055 
1056  // Fractional VMR:
1057  const Numeric VMR_frac = abs_vmrs[h2o_index] / effective_vmr_ref;
1058 
1059  // Check that VMR_frac is inside the allowed range.
1060  {
1061  // FIXME: This check depends on how I interpolate VMR.
1062  const Numeric x_min = nls_pert[0] - 0.5*(nls_pert[1]-nls_pert[0]);
1063  const Numeric x_max = nls_pert[n_nls_pert-1]
1064  + 0.5*(nls_pert[n_nls_pert-1]-nls_pert[n_nls_pert-2]);
1065 
1066  if ( ( VMR_frac > x_max ) ||
1067  ( VMR_frac < x_min ) )
1068  {
1069  ostringstream os;
1070  os << "Problem with gas absorption lookup table.\n"
1071  << "VMR for H2O (species " << h2o_index
1072  << ") is outside the range covered by the lookup table.\n"
1073  << "Your VMR was " << abs_vmrs[h2o_index]
1074  << " at a pressure of " << p << " Pa.\n"
1075  << "The reference VMR value there is " << effective_vmr_ref << "\n"
1076  << "The fractional VMR relative to the reference value is "
1077  << VMR_frac << ".\n"
1078  << "The allowed range is " << x_min << " to " << x_max << ".\n"
1079  << "The fractional VMR perturbation grid range in the table is "
1080  << nls_pert[0] << " to " << nls_pert[n_nls_pert-1] << ".\n"
1081  << "We allow a bit of extrapolation, but NOT SO MUCH!";
1082  throw runtime_error( os.str() );
1083  }
1084  }
1085 
1086  // For now, do linear interpolation in the fractional VMR.
1087  gridpos_poly( vgp_h2o, nls_pert, VMR_frac, h2o_interp_order, extpolfac );
1088  }
1089 
1090 
1091  // Precalculate interpolation weights.
1092  if (n_nls<n_species) {
1093  // Precalculate weights without H2O interpolation if there are less
1094  // nonlinear species than total species. (So at least one species
1095  // without H2O interpolation.)
1096  itw_noH2O.resize(1, 1, n_new_f_grid,
1097  (this_t_interp_order+1)*
1098  (1)* // H2O dimension
1099  (f_interp_order+1));
1100  interpweights(itw_noH2O, *tgp, gp_trivial, *fgp);
1101  }
1102  if (n_nls>0) {
1103  // Precalculate weights with H2O interpolation if there is at least
1104  // one nonlinear species.
1105  itw_withH2O.resize(1, 1, n_new_f_grid,
1106  (this_t_interp_order+1)*
1107  (h2o_interp_order+1)*
1108  (f_interp_order+1));
1109  interpweights(itw_withH2O, *tgp, vgp_h2o, *fgp);
1110  }
1111 
1112  // 7. Loop species:
1113  Index fpi=0;
1114  for ( Index si=0; si<n_species; ++si )
1115  {
1116  // Flag for VMR interpolation, if this is not 0 we want to
1117  // do VMR interpolation:
1118  const Index do_VMR = non_linear[si];
1119 
1120  // For interpolation result.
1121  // Fixed pressure level and species.
1122  // Free dimension is T, H2O, and frequency.
1123  Tensor3View res(xsec_pre_interpolated(pi,si,
1124  Range(joker),
1125  Range(joker),
1126  Range(joker)));
1127 
1128  // Ignore species such as Zeeman and free_electrons which are not
1129  // stored in the lookup table. For those the result is set to 0.
1130  if (is_zeeman(species[si])
1131  || species[si][0].Type() == SpeciesTag::TYPE_FREE_ELECTRONS
1132  || species[si][0].Type() == SpeciesTag::TYPE_PARTICLES)
1133  {
1134  if (do_VMR)
1135  {
1136  ostringstream os;
1137  os << "Problem with gas absorption lookup table.\n"
1138  << "VMR interpolation is not allowed for species \""
1139  << species[si][0].Name() << "\"";
1140  throw runtime_error(os.str());
1141  }
1142  res = 0.;
1143  fpi++;
1144  continue;
1145  }
1146 
1147  // Set h2o related interpolation parameters:
1148  Index this_h2o_extent; // Range of H2O interpolation
1149  if (do_VMR) {
1150  vgp = &vgp_h2o;
1151  this_h2o_extent = n_nls_pert;
1152  itw = &itw_withH2O;
1153  } else {
1154  vgp = &gp_trivial;
1155  this_h2o_extent = 1;
1156  itw = &itw_noH2O;
1157  }
1158 
1159  // Get the right view on xsec.
1160  ConstTensor3View this_xsec
1161  = xsec(Range(joker), // Temperature range
1162  Range(fpi,this_h2o_extent), // VMR profile range
1163  Range(joker), // Frequency range
1164  this_p_grid_index ); // Pressure index
1165 
1166  // Do interpolation.
1167  interp(res, // result
1168  *itw, // weights
1169  this_xsec, // input
1170  *tgp, *vgp, *fgp); // grid positions
1171 
1172  // Increase fpi. fpi marks the position of the first profile
1173  // of the current species in xsec. This is needed to find
1174  // the right subsection of xsec in the presence of nonlinear species.
1175  if (do_VMR)
1176  fpi += n_nls_pert;
1177  else
1178  fpi++;
1179 
1180  } // End of species loop
1181 
1182  // fpi should have reached the end of that dimension of xsec. Check
1183  // this with an assertion:
1184  assert( fpi==xsec.npages() );
1185 
1186  } // End of pressure index loop (below and above gp)
1187 
1188  // Now we have to interpolate between the p_interp_order+1 pressure levels
1189 
1190  // It is a "red" 1D interpolation case we are talking about here.
1191  // (But for a matrix in frequency and species.) Doing a loop over
1192  // frequency and species with an interp call inside would be
1193  // unefficient, so we do this by hand here.
1194  sga.resize(n_species, n_new_f_grid);
1195  sga = 0;
1196  for ( Index pi=0; pi<p_interp_order+1; ++pi )
1197  {
1198  // Multiply pre interpolated quantities with pressure interpolation weights.
1199  // Dimensions of pre_interpolated are:
1200  // Pressure (interpolation points)
1201  // Species
1202  // Temperature (always 1)
1203  // H2O (always 1)
1204  // Frequency
1205  xsec_pre_interpolated(pi,
1206  Range(joker),
1207  Range(joker),
1208  Range(joker),
1209  Range(joker)) *= pitw[pi];
1210 
1211  // Add up in sga.
1212  // Dimensions of sga are (species, frequency)
1213  sga += xsec_pre_interpolated(pi, Range(joker),
1214  0, 0, Range(joker));
1215  }
1216 
1217  // Watch out, this is not yet the final result, we
1218  // need to multiply with the number density of the species, i.e.,
1219  // with the total number density n, times the VMR of the
1220  // species:
1221  for ( Index si=0; si<n_species; ++si )
1222  sga(si,Range(joker)) *= ( n * abs_vmrs[si] );
1223 
1224 
1225  // That's it, we're done!
1226 }
1227 
1228 
1230 {
1231  return f_grid;
1232 }
1233 
1234 
1236 {
1237  return p_grid;
1238 }
1239 
1240 
1242 ostream& operator<< (ostream& os, const GasAbsLookup& /* gal */)
1243 {
1244  os << "GasAbsLookup: Output operator not implemented";
1245  return os;
1246 }
1247 
1248 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
void find_new_grid_in_old_grid(ArrayOfIndex &pos, ConstVectorView old_grid, ConstVectorView new_grid, const Verbosity &verbosity)
Find positions of new grid points in old grid.
Index nelem() const
Number of elements.
Definition: array.h:176
Vector nls_pert
The vector of perturbations for the VMRs of the nonlinear species.
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:2430
void chk_if_decreasing(const String &x_name, ConstVectorView x)
chk_if_decreasing
Definition: check_input.cc:327
Vector t_pert
The vector of temperature perturbations [K].
Declarations having to do with the four output streams.
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:132
The Vector class.
Definition: matpackI.h:556
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
ostream & operator<<(ostream &os, const GasAbsLookup &)
Output operatior for GasAbsLookup.
The Tensor4 class.
Definition: matpackIV.h:383
The range class.
Definition: matpackI.h:148
String get_tag_group_name(const ArrayOfSpeciesTag &tg)
Return the name of a tag group as a string.
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:69
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:679
bool is_zeeman(const ArrayOfSpeciesTag &tg)
Is this a Zeeman tag group?
ArrayOfIndex nonlinear_species
The species tags with non-linear treatment.
Header file for interpolation.cc.
void Adapt(const ArrayOfArrayOfSpeciesTag &current_species, ConstVectorView current_f_grid, const Verbosity &verbosity)
Adapt lookup table to current calculation.
Vector log_p_grid
The natural log of the pressure grid.
void chk_matrix_ncols(const String &x_name, ConstMatrixView x, const Index &l)
chk_matrix_ncols
Definition: check_input.cc:402
Numeric number_density(const Numeric &p, const Numeric &t)
number_density
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
bool is_unique(const ArrayOfIndex &x)
Checks if an ArrayOfIndex is unique, i.e., has no duplicate values.
Definition: logic.cc:352
Tensor4 xsec
Absorption cross sections.
#define DEBUG_ONLY(...)
Definition: debug.h:37
Matrix vmrs_ref
The reference VMR profiles.
Vector f_grid
The frequency grid [Hz].
Index chk_contains(const String &x_name, const Array< T > &x, const T &what)
Check if an array contains a value.
Definition: check_input.h:164
The Tensor3View class.
Definition: matpackIII.h:232
Index find_first_species_tg(const ArrayOfArrayOfSpeciesTag &tgs, const Index &spec)
Find first occurrence of species in tag groups.
An absorption lookup table.
#define abs(x)
Definition: continua.cc:20458
const Vector & GetFgrid() const
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 for the gas absorption lookup table.
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
Set up grid positions for higher order interpolation.
Index species_index_from_species_name(String name)
Return species index for given species name.
Definition: absorption.cc:1260
Vector t_ref
The reference temperature profile [K].
Header file for logic.cc.
Vector p_grid
The pressure grid for the table [Pa].
This can be used to make arrays out of anything.
Definition: array.h:40
Header file for interpolation_poly.cc.
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:465
void Extract(Matrix &sga, const Index &p_interp_order, const Index &t_interp_order, const Index &h2o_interp_order, const Index &f_interp_order, const Numeric &p, const Numeric &T, ConstVectorView abs_vmrs, ConstVectorView new_f_grid, const Numeric &extpolfac) const
Extract scalar gas absorption coefficients from the lookup table.
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
The class MakeVector is a special kind of Vector that can be initialized explicitly from one or more ...
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
A constant view of a Tensor3.
Definition: matpackIII.h:139
A constant view of a Vector.
Definition: matpackI.h:292
void chk_vector_length(const String &x_name, ConstVectorView x, const Index &l)
chk_vector_length
Definition: check_input.cc:233
const Vector & GetPgrid() const
void chk_matrix_nrows(const String &x_name, ConstMatrixView x, const Index &l)
chk_matrix_nrows
Definition: check_input.cc:431
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1838
ArrayOfGridPosPoly fgp_default
Frequency grid positions.
Subclasses of runtime_error.
Definition: check_input.h:123
#define CREATE_OUT3
Definition: messages.h:214
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:91
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:81
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
#define CREATE_OUT2
Definition: messages.h:213
ArrayOfArrayOfSpeciesTag species
The species tags for which the table is valid.
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1403
The Tensor5 class.
Definition: matpackV.h:451
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580