ARTS 2.5.4 (git: 31ce4f0e)
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 "gas_abs_lookup.h"
27#include <cfloat>
28#include <cmath>
29#include "check_input.h"
30#include "interpolation.h"
31#include "logic.h"
32#include "messages.h"
33#include "physics_funcs.h"
34
36
48 ConstVectorView old_grid,
49 ConstVectorView new_grid,
50 const Verbosity& verbosity) {
52
53 const Index n_new_grid = new_grid.nelem();
54 const Index n_old_grid = old_grid.nelem();
55
56 // Make sure that pos has the right size:
57 ARTS_ASSERT(n_new_grid == pos.nelem());
58
59 // Old grid position:
60 Index j = 0;
61
62 // Loop the new frequencies:
63 for (Index i = 0; i < n_new_grid; ++i) {
64 // We have done runtime checks that both the new and the old
65 // frequency grids are sorted in GasAbsLookup::Adapt, so we can
66 // use the fact here.
67
68 while (abs(new_grid[i] - old_grid[j]) >
69 max(abs(new_grid[i]), abs(old_grid[j])) * DBL_EPSILON) {
70 ++j;
71 if (j >= n_old_grid) {
72 ostringstream os;
73 os << "Cannot find new frequency " << i << " (" << new_grid[i]
74 << "Hz) in the lookup table frequency grid.";
75 throw runtime_error(os.str());
76 }
77 }
78
79 pos[i] = j;
80 out3 << " " << new_grid[i] << " found, index = " << pos[i] << ".\n";
81 }
82}
83
85
118 ConstVectorView current_f_grid,
119 const Verbosity& verbosity) {
122
123 // Some constants we will need:
124 const Index n_current_species = current_species.nelem();
125 const Index n_current_f_grid = current_f_grid.nelem();
126
127 const Index n_species = species.nelem();
128 const Index n_nls = nonlinear_species.nelem();
129 const Index n_nls_pert = nls_pert.nelem();
130 const Index n_f_grid = f_grid.nelem();
131 const Index n_p_grid = p_grid.nelem();
132
133 out2 << " Original table: " << n_species << " species, " << n_f_grid
134 << " frequencies.\n"
135 << " Adapt to: " << n_current_species << " species, "
136 << n_current_f_grid << " frequencies.\n";
137
138 if (0 == n_nls) {
139 out2 << " Table contains no nonlinear species.\n";
140 }
141
142 // Set up a logical array for the nonlinear species
143 ArrayOfIndex non_linear(n_species, 0);
144 for (Index s = 0; s < n_nls; ++s) {
145 non_linear[nonlinear_species[s]] = 1;
146 }
147
148 if (0 == t_pert.nelem()) {
149 out2 << " Table contains no temperature perturbations.\n";
150 }
151
152 // We are constructing a new lookup table, containing just the
153 // species and frequencies that are necessary for the current
154 // calculation. We will build it in this local variable, then copy
155 // it back to *this in the end.
156 GasAbsLookup new_table;
157
158 // First some checks on the lookup table itself:
159
160 // Species:
161 if (0 == n_species) {
162 ostringstream os;
163 os << "The lookup table should have at least one species.";
164 throw runtime_error(os.str());
165 }
166
167 // Nonlinear species:
168 // They should be unique ...
170 ostringstream os;
171 os << "The table must not have duplicate nonlinear species.\n"
172 << "Value of *nonlinear_species*: " << nonlinear_species;
173 throw runtime_error(os.str());
174 }
175
176 // ... and pointing at valid species.
177 for (Index i = 0; i < n_nls; ++i) {
178 ostringstream os;
179 os << "nonlinear_species[" << i << "]";
180 chk_if_in_range(os.str(), nonlinear_species[i], 0, n_species - 1);
181 }
182
183 // Frequency grid:
184 chk_if_increasing("f_grid", f_grid);
185
186 // Pressure grid:
187 chk_if_decreasing("p_grid", p_grid);
188
189 // Reference VMRs:
190 chk_matrix_nrows("vmrs_ref", vmrs_ref, n_species);
191 chk_matrix_ncols("vmrs_ref", vmrs_ref, n_p_grid);
192
193 // Reference temperatur:
194 chk_vector_length("t_ref", t_ref, n_p_grid);
195
196 // Temperature perturbations:
197 // Nothing to check for t_pert, it seems.
198
199 // Perturbations for nonlinear species:
200 // Check that nls_pert is empty if and only if nonlinear_species is
201 // empty:
202 if (0 == n_nls) {
203 chk_vector_length("nls_pert", nls_pert, 0);
204 } else {
205 if (0 == n_nls_pert) {
206 ostringstream os;
207 os << "The vector nls_pert should contain the perturbations\n"
208 << "for the nonlinear species, but it is empty.";
209 throw runtime_error(os.str());
210 }
211 }
212
213 // The table itself, xsec:
214 //
215 // We have to separtely consider the 3 cases described in the
216 // documentation of GasAbsLookup.
217 //
218 // Dimension: [ a, b, c, d ]
219 //
220 if (0 == n_nls) {
221 if (0 == t_pert.nelem()) {
222 // Simplest case (no temperature perturbations,
223 // no vmr perturbations):
224 // a = 1
225 // b = n_species
226 // c = n_f_grid
227 // d = n_p_grid
228 chk_size("xsec", xsec, 1, n_species, n_f_grid, n_p_grid);
229 } else {
230 // Standard case (temperature perturbations,
231 // but no vmr perturbations):
232 // a = n_t_pert
233 // b = n_species
234 // c = n_f_grid
235 // d = n_p_grid
236 chk_size("xsec", xsec, t_pert.nelem(), n_species, n_f_grid, n_p_grid);
237 }
238 } else {
239 // Full case (with temperature perturbations and
240 // vmr perturbations):
241 // a = n_t_pert
242 // b = n_species + n_nonlinear_species * ( n_nls_pert - 1 )
243 // c = n_f_grid
244 // d = n_p_grid
245 Index a = t_pert.nelem();
246 Index b = n_species + n_nls * (n_nls_pert - 1);
247 Index c = n_f_grid;
248 Index d = n_p_grid;
249
250 chk_size("xsec", xsec, a, b, c, d);
251 }
252
253 // We also need indices to the positions of the original species
254 // data in xsec. Nonlinear species take more space, therefor the
255 // position in xsec is not the same as the position in species.
256 ArrayOfIndex original_spec_pos_in_xsec(n_species);
257 for (Index i = 0, sp = 0; i < n_species; ++i) {
258 original_spec_pos_in_xsec[i] = sp;
259 if (non_linear[i])
260 sp += n_nls_pert;
261 else
262 sp += 1;
263 }
264
265 // Now some checks on the input data:
266
267 // The list of current species should not be empty:
268 if (0 == n_current_species) {
269 ostringstream os;
270 os << "The list of current species should not be empty.";
271 throw runtime_error(os.str());
272 }
273
274 // The grid of current frequencies should be monotonically sorted:
275 chk_if_increasing("current_f_grid", current_f_grid);
276
277 // 1. Find and remember the indices of the current species in the
278 // lookup table. At the same time verify that each species is
279 // included in the table exactly once.
280 ArrayOfIndex i_current_species(n_current_species);
281 out3 << " Looking for species in lookup table:\n";
282 for (Index i = 0; i < n_current_species; ++i) {
283 out3 << " " << current_species[i].Name() << ": ";
284
285 try {
286 i_current_species[i] =
287 chk_contains("abs_species", species, current_species[i]);
288
289 } catch (const runtime_error_not_found&) {
290 // Is this one of the trivial species?
291 const auto spec_type = current_species[i].Type();
292 if (spec_type == Species::TagType::Zeeman ||
293 spec_type == Species::TagType::FreeElectrons ||
294 spec_type == Species::TagType::Particles) {
295 // Set to -1 to flag that this needs special handling later on.
296 i_current_species[i] = -1;
297 } else {
298 ostringstream os;
299 os << "Species " << current_species[i].Name()
300 << " is missing in absorption lookup table.";
301 throw runtime_error(os.str());
302 }
303 }
304
305 out3 << "found (or trivial).\n";
306 }
307
308 // 1a. Find out which of the current species are nonlinear species:
309 Index n_current_nonlinear_species = 0; // Number of current nonlinear species
310 ArrayOfIndex current_non_linear(n_current_species, 0); // A logical array to
311 // flag which of the
312 // current species are
313 // nonlinear.
314
315 out3 << " Finding out which of the current species are nonlinear:\n";
316 for (Index i = 0; i < n_current_species; ++i) {
317 if (i_current_species[i] >= 0) // Jump over trivial species here.
318 {
319 // Check if this is a nonlinear species:
320 if (non_linear[i_current_species[i]]) {
321 out3 << " " << current_species[i] << "\n";
322
323 current_non_linear[i] = 1;
324 ++n_current_nonlinear_species;
325 }
326 }
327 }
328
329 // 2. Find and remember the frequencies of the current calculation in
330 // the lookup table. At the same time verify that all frequencies are
331 // included and that no frequency occurs twice.
332
333 // FIXME: This is a bit tricky, because we are comparing
334 // Numerics. Let's see how well this works in practice.
335
336 ArrayOfIndex i_current_f_grid(n_current_f_grid);
337 out3 << " Looking for Frequencies in lookup table:\n";
338
339 // We need no error checking for the next statement, since the
340 // function called throws a runtime error if a frequency
341 // is not found, or if the grids are not ok.
343 i_current_f_grid, f_grid, current_f_grid, verbosity);
344
345 // 3. Use the species and frequency index lists to build the new lookup
346 // table.
347
348 // Species:
349 new_table.species.resize(n_current_species);
350 for (Index i = 0; i < n_current_species; ++i) {
351 if (i_current_species[i] >= 0) {
352 new_table.species[i] = species[i_current_species[i]];
353
354 // Is this a nonlinear species?
355 if (current_non_linear[i]) new_table.nonlinear_species.push_back(i);
356 } else {
357 // Here we handle the case of the trivial species, for which we simply
358 // copy the name:
359
360 new_table.species[i] = current_species[i];
361 }
362 }
363
364 // Frequency grid:
365 new_table.f_grid.resize(n_current_f_grid);
366 for (Index i = 0; i < n_current_f_grid; ++i) {
367 new_table.f_grid[i] = f_grid[i_current_f_grid[i]];
368 }
369
370 // Pressure grid:
371 // new_table.p_grid.resize( n_p_grid );
372 new_table.p_grid = p_grid;
373
374 // Reference VMR profiles:
375 new_table.vmrs_ref.resize(n_current_species, n_p_grid);
376 for (Index i = 0; i < n_current_species; ++i) {
377 if (i_current_species[i] >= 0) {
378 new_table.vmrs_ref(i, Range(joker)) =
379 vmrs_ref(i_current_species[i], Range(joker));
380 } else {
381 // Here we handle the case of the trivial species, for which we set
382 // the reference VMR to NAN.
383 new_table.vmrs_ref(i, Range(joker)) = NAN;
384 }
385 }
386
387 // Reference temperature profile:
388 // new_table.t_ref.resize( t_ref.nelem() );
389 new_table.t_ref = t_ref;
390
391 // Vector of temperature perturbations:
392 // new_table.t_pert.resize( t_pert.nelem() );
393 new_table.t_pert = t_pert;
394
395 // Vector of perturbations for the VMRs of the nonlinear species:
396 // (Should stay empty if we have no nonlinear species)
397 if (0 != new_table.nonlinear_species.nelem()) {
398 // new_table.nls_pert.resize( n_nls_pert );
399 new_table.nls_pert = nls_pert;
400 }
401
402 // Absorption coefficients:
403 new_table.xsec.resize(
404 xsec.nbooks(),
405 n_current_species + n_current_nonlinear_species * (n_nls_pert - 1),
406 n_current_f_grid,
407 xsec.ncols());
408
409 // We have to copy the right species and frequencies from the old to
410 // the new table. Temperature perturbations and pressure grid remain
411 // the same.
412
413 // Do species:
414 for (Index i_s = 0, sp = 0; i_s < n_current_species; ++i_s) {
415 // n_v is the number of VMR perturbations
416 Index n_v;
417 if (current_non_linear[i_s])
418 n_v = n_nls_pert;
419 else
420 n_v = 1;
421
422 // cout << "i_s / sp / n_v = " << i_s << " / " << sp << " / " << n_v << endl;
423 // cout << "orig_pos = " << original_spec_pos_in_xsec[i_current_species[i_s]] << endl;
424
425 // Do frequencies:
426 for (Index i_f = 0; i_f < n_current_f_grid; ++i_f) {
427 if (i_current_species[i_s] >= 0) {
428 new_table.xsec(Range(joker), Range(sp, n_v), i_f, Range(joker)) =
430 Range(original_spec_pos_in_xsec[i_current_species[i_s]], n_v),
431 i_current_f_grid[i_f],
432 Range(joker));
433 } else {
434 // Here we handle the case of the trivial species, which we simply
435 // set to NAN:
436 new_table.xsec(Range(joker), Range(sp, n_v), i_f, Range(joker)) = NAN;
437 }
438
439 // cout << "result: " << xsec( Range(joker),
440 // Range(original_spec_pos_in_xsec[i_current_species[i_s]],n_v),
441 // i_current_f_grid[i_f],
442 // Range(joker) ) << endl;
443 }
444
445 sp += n_v;
446 }
447
448 // 4. Replace original table by the new one.
449 *this = new_table;
450
451 // 5. Initialize log_p_grid.
452 log_p_grid.resize(n_p_grid);
454
455 // 6. Initialize flag_default.
457}
458
460
514 const ArrayOfSpeciesTag& select_abs_species,
515 const Index& p_interp_order,
516 const Index& t_interp_order,
517 const Index& h2o_interp_order,
518 const Index& f_interp_order,
519 const Numeric& p,
520 const Numeric& T,
521 ConstVectorView abs_vmrs,
522 ConstVectorView new_f_grid,
523 const Numeric& extpolfac) const {
524 // 1. Obtain some properties of the lookup table:
525
526 // Number of gas species in the table:
527 const Index n_species = species.nelem();
528
529 // Number of nonlinear species:
530 const Index n_nls = nonlinear_species.nelem();
531
532 // Number of frequencies in the table:
533 const Index n_f_grid = f_grid.nelem();
534
535 // Number of pressure grid points in the table:
536 const Index n_p_grid = p_grid.nelem();
537
538 // Number of temperature perturbations:
539 const Index n_t_pert = t_pert.nelem();
540
541 // Number of nonlinear species perturbations:
542 const Index n_nls_pert = nls_pert.nelem();
543
544 // Number of frequencies in new_f_grid, the frequency grid for which we
545 // want to extract.
546 const Index n_new_f_grid = new_f_grid.nelem();
547
548 // 2. First some checks on the lookup table itself:
549
550 // Most checks here are asserts, because they check the internal
551 // consistency of the lookup table. They should never fail if the
552 // table has been created with ARTS.
553
554 // If there are nonlinear species, then at least one species must be
555 // H2O. We will use that to perturb in the case of nonlinear
556 // species.
557 Index h2o_index = -1;
558 if (n_nls > 0) {
559 h2o_index = find_first_species(species, Species::Species::Water);
560
561 // This is a runtime error, even though it would be more logical
562 // for it to be an assertion, since it is an internal check on
563 // the table. The reason is that it is somewhat awkward to check
564 // for this in other places.
565 if (h2o_index == -1) {
566 ostringstream os;
567 os << "With nonlinear species, at least one species must be a H2O species.";
568 throw runtime_error(os.str());
569 }
570 }
571
572 // Check that the dimension of vmrs_ref is consistent with species and p_grid:
573 ARTS_ASSERT(is_size(vmrs_ref, n_species, n_p_grid));
574
575 // Check dimension of t_ref:
576 ARTS_ASSERT(is_size(t_ref, n_p_grid));
577
578 // Check dimension of xsec:
579 DEBUG_ONLY({
580 Index a, b, c, d;
581 if (0 == n_t_pert)
582 a = 1;
583 else
584 a = n_t_pert;
585 b = n_species + n_nls * (n_nls_pert - 1);
586 c = n_f_grid;
587 d = n_p_grid;
588 // cout << "xsec: "
589 // << xsec.nbooks() << ", "
590 // << xsec.npages() << ", "
591 // << xsec.nrows() << ", "
592 // << xsec.ncols() << "\n";
593 // cout << "a b c d: "
594 // << a << ", "
595 // << b << ", "
596 // << c << ", "
597 // << d << "\n";
598 ARTS_ASSERT(is_size(xsec, a, b, c, d));
599 })
600
601 // Make sure that log_p_grid is initialized:
602 if (log_p_grid.nelem() != n_p_grid) {
603 ostringstream os;
604 os << "The lookup table internal variable log_p_grid is not initialized.\n"
605 << "Use the abs_lookupAdapt method!";
606 throw runtime_error(os.str());
607 }
608
609 // Verify that we have enough pressure, temperature,humdity, and frequency grid points
610 // for the desired interpolation orders. This check is not only
611 // table internal, since abs_nls_interp_order and abs_t_interp_order
612 // are separate WSVs that could have been modified. Hence, these are
613 // runtime errors.
614
615 if ((n_p_grid < p_interp_order + 1)) {
616 ostringstream os;
617 os << "The number of pressure grid points in the table (" << n_p_grid
618 << ") is not enough for the desired order of interpolation ("
619 << p_interp_order << ").";
620 throw runtime_error(os.str());
621 }
622
623 if ((n_nls_pert != 0) && (n_nls_pert < h2o_interp_order + 1)) {
624 ostringstream os;
625 os << "The number of humidity perturbation grid points in the table ("
626 << n_nls_pert
627 << ") is not enough for the desired order of interpolation ("
628 << h2o_interp_order << ").";
629 throw runtime_error(os.str());
630 }
631
632 if ((n_t_pert != 0) && (n_t_pert < t_interp_order + 1)) {
633 ostringstream os;
634 os << "The number of temperature perturbation grid points in the table ("
635 << n_t_pert << ") is not enough for the desired order of interpolation ("
636 << t_interp_order << ").";
637 throw runtime_error(os.str());
638 }
639
640 if ((n_f_grid < f_interp_order + 1)) {
641 ostringstream os;
642 os << "The number of frequency grid points in the table (" << n_f_grid
643 << ") is not enough for the desired order of interpolation ("
644 << f_interp_order << ").";
645 throw runtime_error(os.str());
646 }
647
648 // 3. Checks on the input variables:
649
650 // Check that abs_vmrs has the right dimension:
651 if (!is_size(abs_vmrs, n_species)) {
652 ostringstream os;
653 os << "Number of species in lookup table does not match number\n"
654 << "of species for which you want to extract absorption.\n"
655 << "Have you used abs_lookupAdapt? Or did you miss to add\n"
656 << "some VRM fields (e.g. for free electrons or particles)?\n";
657 throw runtime_error(os.str());
658 }
659
660 // 4. Set up some things we will need later on:
661
662 // 4.a Frequency grid positions
663
664 // Frequency grid positions. The pointer is used to save copying of the
665 // default from the lookup table.
668
669 // With f_interp_order 0 the frequency grid has to have the same size as in the
670 // lookup table, or exactly one element. If it matches the lookup table, we
671 // do no frequency interpolation at all. (We set the frequency grid positions
672 // to the predefined ones that come with the lookup table.)
673 if (f_interp_order == 0) {
674
675 // We do some superficial checks below, to make sure that the
676 // frequency grid is the same as in the lookup table (only first
677 // and last element of f_grid are tested). As margin for
678 // agreement, we pick a value that is just slightly smaller than
679 // the perturbation that is used by the wind jacobian, which is 0.1.
680 const Numeric allowed_f_margin = 0.09;
681
682 if (n_new_f_grid == n_f_grid) {
683 // Use the default flag that is stored in the lookup table itself
684 // (which effectively means no frequency interpolation)
686
687 // Check first f_grid element:
688 if (abs(f_grid[0] - new_f_grid[0]) > allowed_f_margin)
689 {
690 ostringstream os;
691 os << "First frequency in f_grid inconsistent with lookup table.\n"
692 << "f_grid[0] = " << f_grid[0] << "\n"
693 << "new_f_grid[0] = " << new_f_grid[0] << ".";
694 throw runtime_error(os.str());
695 }
696
697 // Check last f_grid element:
698 if (abs(f_grid[n_f_grid - 1] - new_f_grid[n_new_f_grid - 1]) > allowed_f_margin)
699 {
700 ostringstream os;
701 os << "Last frequency in f_grid inconsistent with lookup table.\n"
702 << "f_grid[n_f_grid-1] = " << f_grid[n_f_grid - 1]
703 << "\n"
704 << "new_f_grid[n_new_f_grid-1] = " << new_f_grid[n_new_f_grid - 1]
705 << ".";
706 throw runtime_error(os.str());
707 }
708 } else if (n_new_f_grid == 1) {
709 flag = &flag_local;
710 flag_local = Interpolation::LagrangeVector(new_f_grid, f_grid, 0);
711
712 // Check that we really are on a frequency grid point, for safety's sake.
713 if (abs(f_grid[flag_local[0].pos] - new_f_grid[0]) > allowed_f_margin)
714 {
715 ostringstream os;
716 os << "Cannot find a matching lookup table frequency for frequency "
717 << new_f_grid[0] << ".\n"
718 << "(This check has not been properly tested, so perhaps this is\n"
719 << "a false alarm. Check for this in file gas_abs_lookup.cc.)";
720 throw runtime_error(os.str());
721 }
722 } else {
723 ostringstream os;
724 os << "With f_interp_order 0 the frequency grid has to have the same\n"
725 << "size as in the lookup table, or exactly one element.";
726 throw runtime_error(os.str());
727 }
728 } else {
729 const Numeric f_min = f_grid[0] - 0.5 * (f_grid[1] - f_grid[0]);
730 const Numeric f_max = f_grid[n_f_grid - 1] +
731 0.5 * (f_grid[n_f_grid - 1] - f_grid[n_f_grid - 2]);
732 if (new_f_grid[0] < f_min) {
733 ostringstream os;
734 os << "Problem with gas absorption lookup table.\n"
735 << "At least one frequency is outside the range covered by the lookup table.\n"
736 << "Your new frequency value is " << new_f_grid[0] << " Hz.\n"
737 << "The allowed range is " << f_min << " to " << f_max << " Hz.";
738 throw runtime_error(os.str());
739 }
740 if (new_f_grid[n_new_f_grid - 1] > f_max) {
741 ostringstream os;
742 os << "Problem with gas absorption lookup table.\n"
743 << "At least one frequency is outside the range covered by the lookup table.\n"
744 << "Your new frequency value is " << new_f_grid[n_new_f_grid - 1]
745 << " Hz.\n"
746 << "The allowed range is " << f_min << " to " << f_max << " Hz.";
747 throw runtime_error(os.str());
748 }
749
750 // We do have real frequency interpolation (f_interp_order!=0).
751 flag = &flag_local;
752 flag_local = Interpolation::LagrangeVector(new_f_grid, f_grid, f_interp_order);
753 }
754
755 // 4.b Other stuff
756
757 // Flag for temperature interpolation, if this is not 0 we want
758 // to do T interpolation:
759 const Index do_T = n_t_pert;
760
761 // Set up a logical array for the nonlinear species
762 ArrayOfIndex non_linear(n_species, 0);
763 for (Index s = 0; s < n_nls; ++s) {
764 non_linear[nonlinear_species[s]] = 1;
765 }
766
767 // Calculate the number density for the given pressure and
768 // temperature:
769 // n = n0*T0/p0 * p/T or n = p/kB/t, ideal gas law
770 const Numeric n = number_density(p, T);
771
772 // 5. Determine pressure grid position and interpolation weights:
773
774 // Check that p is inside the grid. (p_grid is sorted in decreasing order.)
775 {
776 const Numeric p_max = p_grid[0] + 0.5 * (p_grid[0] - p_grid[1]);
777 const Numeric p_min = p_grid[n_p_grid - 1] -
778 0.5 * (p_grid[n_p_grid - 2] - p_grid[n_p_grid - 1]);
779 if ((p > p_max) || (p < p_min)) {
780 ostringstream os;
781 os << "Problem with gas absorption lookup table.\n"
782 << "Pressure p is outside the range covered by the lookup table.\n"
783 << "Your p value is " << p << " Pa.\n"
784 << "The allowed range is " << p_min << " to " << p_max << ".\n"
785 << "The pressure grid range in the table is " << p_grid[n_p_grid - 1]
786 << " to " << p_grid[0] << ".\n"
787 << "We allow a bit of extrapolation, but NOT SO MUCH!";
788 throw runtime_error(os.str());
789 }
790 }
791
792 // For sure, we need to store the pressure grid position.
793 // We do the interpolation in log(p). Test have shown that this
794 // gives slightly better accuracy than interpolating in p directly.
795 const auto plag = Interpolation::LagrangeVector(log(p), log_p_grid, p_interp_order);
796
797 // Pressure interpolation weights:
798 const auto pitw = interpweights(plag[0]);
799
800 // Define also other grid positions and interpolation weights here, so that
801 // we do not have to allocate them over and over in the loops below.
802
803 // Define the ArrayOfLagrangeInterpolation that corresponds to "no interpolation at all".
804 const ArrayOfLagrangeInterpolation lag_trivial(1);
805
806 // Temperature grid positions.
807 ArrayOfLagrangeInterpolation tlag_withT(1); // Only a scalar.
808 const ArrayOfLagrangeInterpolation* tlag; // Pointer to either tlag_withT or lag_trivial.
809
810 // Set this_t_interp_order, depending on whether we do T interpolation or not.
811 if (do_T) {
812 tlag = &tlag_withT;
813 } else {
814 // For the !do_T case we simply take the single
815 // temperature that is there, so we point tgp accordingly.
816 tlag = &lag_trivial;
817 }
818
819 // H2O(VMR) grid positions. vlag is what will be used in the interpolation.
820 // Depending on species, it is either pointed to lag_trivial, or to vlag_h2o.
822 ArrayOfLagrangeInterpolation vlag_h2o(1); // only a scalar
823
824 // 6. We do the T and VMR interpolation for the pressure levels
825 // that are used in the pressure interpolation. (How many depends on
826 // p_interp_order.)
827
828 // To store the interpolated result for the p_interp_order+1
829 // pressure levels:
830 // xsec dimensions are:
831 // Temperature
832 // H2O
833 // Frequency
834 // Pressure
835 // Dimensions of pre_interpolated are:
836 // Pressure (interpolation points)
837 // Species
838 // Temperature (always 1)
839 // H2O (always 1)
840 // Frequency
841
842 Tensor5 xsec_pre_interpolated;
843 xsec_pre_interpolated.resize(
844 p_interp_order + 1, n_species, 1, 1, n_new_f_grid);
845
846 // Define variables for interpolation weights outside the loops.
847 // We will make itw point to either the weights with H2O interpolation, or
848 // the ones without.
849 Tensor3OfTensor3 itw_withH2O(0,0,0), itw_noH2O(0,0,0);
850 const Tensor3OfTensor3 *itw;
851
852 for (Index pi = 0; pi < p_interp_order + 1; ++pi) {
853 // Throw a runtime error if one of the reference VMR profiles is zero, but
854 // abs_vmrs is not. (This means that the lookup table was calculated with a
855 // reference profile of zero for that gas.)
856 // for (Index si=0; si<n_species; ++si)
857 // if ( (vmrs_ref(si,pi) == 0) &&
858 // (abs_vmrs[si] != 0) )
859 // {
860 // ostringstream os;
861 // os << "Reference VMR profile is zero, you cannot extract\n"
862 // << "Absorption for this species.\n"
863 // << "Species: " << si
864 // << " (" << get_species_name(species[si]) << ")\n"
865 // << "Lookup table pressure level: " << pi
866 // << " (" << p_grid[pi] << " Pa).";
867 // throw runtime_error( os.str() );
868 // }
869
870 // Index into p_grid:
871 const Index this_p_grid_index = plag[0].pos + pi;
872
873 // Determine temperature grid position. This is only done if we
874 // want temperature interpolation, but the variable tgp has to
875 // be visible also outside for later use:
876 if (do_T) {
877 // Temperature in the atmosphere is altitude
878 // dependent. When we do the interpolation for the pressure level
879 // below and above our point, we should correct the target value of
880 // the interpolation to the altitude (pressure) difference. This
881 // ensures that there is for example no T interpolation if the
882 // desired T is right on the reference profile curve.
883 //
884 // I explicitly compared this with the old option to calculate
885 // the temperature offset relative to the temperature at
886 // this level. The performance in both cases is very
887 // similar. The reason, why I decided to keep this new
888 // version, is that it avoids the problem of needing
889 // oversized temperature perturbations if the pressure
890 // grid is coarse.
891 //
892 // No! The above approach leads to problems when combined with
893 // higher order pressure interpolation. The problem is that
894 // the reference T and VMR profiles may be very
895 // irregular. (For example the H2O profile often has a big
896 // jump near the bottom.) That sometimes leads to negative
897 // effective reference values when the reference profile is
898 // interpolated. I therefore reverted back to the original
899 // version of using the real temperature and humidity, not
900 // the interpolated one.
901
902 // const Numeric effective_T_ref = interp(pitw,t_ref,pgp);
903 const Numeric effective_T_ref = t_ref[this_p_grid_index];
904
905 // Convert temperature to offset from t_ref:
906 const Numeric T_offset = T - effective_T_ref;
907
908 // cout << "T_offset = " << T_offset << endl;
909
910 // Check that temperature offset is inside the allowed range.
911 {
912 const Numeric t_min = t_pert[0] - extpolfac * (t_pert[1] - t_pert[0]);
913 const Numeric t_max =
914 t_pert[n_t_pert - 1] +
915 extpolfac * (t_pert[n_t_pert - 1] - t_pert[n_t_pert - 2]);
916 if ((T_offset > t_max) || (T_offset < t_min)) {
917 ostringstream os;
918 os << "Problem with gas absorption lookup table.\n"
919 << "Temperature T is outside the range covered by the lookup table.\n"
920 << "Your temperature was " << T << " K at a pressure of " << p
921 << " Pa.\n"
922 << "The temperature offset value is " << T_offset << ".\n"
923 << "The allowed range is " << t_min << " to " << t_max << ".\n"
924 << "The temperature perturbation grid range in the table is "
925 << t_pert[0] << " to " << t_pert[n_t_pert - 1] << ".\n"
926 << "We allow a bit of extrapolation, but NOT SO MUCH!";
927 throw runtime_error(os.str());
928 }
929 }
930
931 tlag_withT[0] = LagrangeInterpolation(0, T_offset, t_pert, t_interp_order);
932 }
933
934 // Determine the H2O VMR grid position. We need to do this only
935 // once, since the only species who's VMR is interpolated is
936 // H2O. We do this only if there are nonlinear species, but the
937 // variable has to be visible later.
938 if (n_nls > 0) {
939 // Similar to the T case, we first interpolate the reference
940 // VMR to the pressure of extraction, then compare with
941 // the extraction VMR to determine the offset/fractional
942 // difference for the VMR interpolation.
943 //
944 // No! The above approach leads to problems when combined with
945 // higher order pressure interpolation. The problem is that
946 // the reference T and VMR profiles may be very
947 // irregular. (For example the H2O profile often has a big
948 // jump near the bottom.) That sometimes leads to negative
949 // effective reference values when the reference profile is
950 // interpolated. I therefore reverted back to the original
951 // version of using the real temperature and humidity, not
952 // the interpolated one.
953
954 // const Numeric effective_vmr_ref = interp(pitw,
955 // vmrs_ref(h2o_index, Range(joker)),
956 // pgp);
957 const Numeric effective_vmr_ref = vmrs_ref(h2o_index, this_p_grid_index);
958
959 // Fractional VMR:
960 const Numeric VMR_frac = abs_vmrs[h2o_index] / effective_vmr_ref;
961
962 // Check that VMR_frac is inside the allowed range.
963 {
964 // FIXME: This check depends on how I interpolate VMR.
965 const Numeric x_min =
966 nls_pert[0] - extpolfac * (nls_pert[1] - nls_pert[0]);
967 const Numeric x_max =
968 nls_pert[n_nls_pert - 1] +
969 extpolfac * (nls_pert[n_nls_pert - 1] - nls_pert[n_nls_pert - 2]);
970
971 if ((VMR_frac > x_max) || (VMR_frac < x_min)) {
972 ostringstream os;
973 os << "Problem with gas absorption lookup table.\n"
974 << "VMR for H2O (species " << h2o_index
975 << ") is outside the range covered by the lookup table.\n"
976 << "Your VMR was " << abs_vmrs[h2o_index] << " at a pressure of "
977 << p << " Pa.\n"
978 << "The reference VMR value there is " << effective_vmr_ref << "\n"
979 << "The fractional VMR relative to the reference value is "
980 << VMR_frac << ".\n"
981 << "The allowed range is " << x_min << " to " << x_max << ".\n"
982 << "The fractional VMR perturbation grid range in the table is "
983 << nls_pert[0] << " to " << nls_pert[n_nls_pert - 1] << ".\n"
984 << "We allow a bit of extrapolation, but NOT SO MUCH!";
985 throw runtime_error(os.str());
986 }
987 }
988
989 // For now, do linear interpolation in the fractional VMR.
990 vlag_h2o[0] = LagrangeInterpolation(0, VMR_frac, nls_pert, h2o_interp_order);
991 }
992
993 // Precalculate interpolation weights.
994 if (n_nls < n_species) {
995 // Precalculate weights without H2O interpolation if there are less
996 // nonlinear species than total species. (So at least one species
997 // without H2O interpolation.)
998 itw_noH2O = interpweights(*tlag, lag_trivial, *flag);
999 }
1000 if (n_nls > 0) {
1001 // Precalculate weights with H2O interpolation if there is at least
1002 // one nonlinear species.
1003 itw_withH2O = interpweights(*tlag, vlag_h2o, *flag);
1004 }
1005
1006 // 7. Loop species:
1007 Index fpi = 0;
1008 for (Index si = 0; si < n_species; ++si) {
1009 // Flag for VMR interpolation, if this is not 0 we want to
1010 // do VMR interpolation:
1011 const Index do_VMR = non_linear[si];
1012
1013 // For interpolation result.
1014 // Fixed pressure level and species.
1015 // Free dimension is T, H2O, and frequency.
1016 Tensor3View res(xsec_pre_interpolated(
1017 pi, si, Range(joker), Range(joker), Range(joker)));
1018
1019 // Ignore species such as Zeeman and free_electrons which are not
1020 // stored in the lookup table. For those the result is set to 0.
1021 if (species[si].Zeeman() or species[si].FreeElectrons() or species[si].Particles()) {
1022 if (do_VMR) {
1023 ostringstream os;
1024 os << "Problem with gas absorption lookup table.\n"
1025 << "VMR interpolation is not allowed for species \""
1026 << species[si][0].Name() << "\"";
1027 throw runtime_error(os.str());
1028 }
1029 res = 0.;
1030 fpi++;
1031 continue;
1032 }
1033
1034 // Set h2o related interpolation parameters:
1035 Index this_h2o_extent; // Range of H2O interpolation
1036 if (do_VMR) {
1037 vlag = &vlag_h2o;
1038 this_h2o_extent = n_nls_pert;
1039 itw = &itw_withH2O;
1040 } else {
1041 vlag = &lag_trivial;
1042 this_h2o_extent = 1;
1043 itw = &itw_noH2O;
1044 }
1045
1046 // Get the right view on xsec.
1047 ConstTensor3View this_xsec =
1048 xsec(Range(joker), // Temperature range
1049 Range(fpi, this_h2o_extent), // VMR profile range
1050 Range(joker), // Frequency range
1051 this_p_grid_index); // Pressure index
1052
1053 // Do interpolation.
1054 reinterp(res, // result
1055 this_xsec, // input
1056 *itw, // weights
1057 *tlag,
1058 *vlag,
1059 *flag); // grid positions
1060
1061 // Increase fpi. fpi marks the position of the first profile
1062 // of the current species in xsec. This is needed to find
1063 // the right subsection of xsec in the presence of nonlinear species.
1064 if (do_VMR)
1065 fpi += n_nls_pert;
1066 else
1067 fpi++;
1068
1069 } // End of species loop
1070
1071 // fpi should have reached the end of that dimension of xsec. Check
1072 // this with an assertion:
1073 ARTS_ASSERT(fpi == xsec.npages());
1074
1075 } // End of pressure index loop (below and above gp)
1076
1077 // Now we have to interpolate between the p_interp_order+1 pressure levels
1078
1079 // It is a "red" 1D interpolation case we are talking about here.
1080 // (But for a matrix in frequency and species.) Doing a loop over
1081 // frequency and species with an interp call inside would be
1082 // unefficient, so we do this by hand here.
1083 sga.resize(n_species, n_new_f_grid);
1084 sga = 0;
1085 for (Index pi = 0; pi < p_interp_order + 1; ++pi) {
1086 // Multiply pre interpolated quantities with pressure interpolation weights.
1087 // Dimensions of pre_interpolated are:
1088 // Pressure (interpolation points)
1089 // Species
1090 // Temperature (always 1)
1091 // H2O (always 1)
1092 // Frequency
1093 xsec_pre_interpolated(
1094 pi, Range(joker), Range(joker), Range(joker), Range(joker)) *= pitw[pi];
1095
1096 // Add up in sga.
1097 // Dimensions of sga are (species, frequency)
1098 sga += xsec_pre_interpolated(pi, Range(joker), 0, 0, Range(joker));
1099 }
1100
1101 // Watch out, this is not yet the final result, we
1102 // need to multiply with the number density of the species, i.e.,
1103 // with the total number density n, times the VMR of the
1104 // species:
1105 for (Index si = 0; si < n_species; ++si) {
1106 if (select_abs_species.nelem()) {
1107 if (species[si] == select_abs_species)
1108 sga(si, Range(joker)) *= (n * abs_vmrs[si]);
1109 else
1110 sga(si, Range(joker)) = 0.;
1111 } else {
1112 sga(si, Range(joker)) *= (n * abs_vmrs[si]);
1113 }
1114 }
1115
1116 // That's it, we're done!
1117}
1118
1119const Vector& GasAbsLookup::GetFgrid() const { return f_grid; }
1120
1121const Vector& GasAbsLookup::GetPgrid() const { return p_grid; }
1122
1124ostream& operator<<(ostream& os, const GasAbsLookup& /* gal */) {
1125 os << "GasAbsLookup: Output operator not implemented";
1126 return os;
1127}
void chk_matrix_ncols(const String &x_name, ConstMatrixView x, const Index &l)
chk_matrix_ncols
Definition: check_input.cc:380
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:421
void chk_matrix_nrows(const String &x_name, ConstMatrixView x, const Index &l)
chk_matrix_nrows
Definition: check_input.cc:399
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:86
void chk_if_decreasing(const String &x_name, ConstVectorView x)
chk_if_decreasing
Definition: check_input.cc:325
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:111
void chk_vector_length(const String &x_name, ConstVectorView x, const Index &l)
chk_vector_length
Definition: check_input.cc:257
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:149
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Tensor3.
Definition: matpackIII.h:130
Index ncols() const noexcept
Definition: matpackIV.h:142
Index nbooks() const noexcept
Definition: matpackIV.h:139
Index npages() const noexcept
Definition: matpackIV.h:140
A constant view of a Vector.
Definition: matpackI.h:512
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
An absorption lookup table.
Vector p_grid
The pressure grid for the table [Pa].
Tensor4 xsec
Absorption cross sections.
Vector log_p_grid
The natural log of the pressure grid.
Vector t_pert
The vector of temperature perturbations [K].
const Vector & GetPgrid() const
Vector f_grid
The frequency grid [Hz].
Matrix vmrs_ref
The reference VMR profiles.
const Vector & GetFgrid() const
ArrayOfLagrangeInterpolation flag_default
Frequency grid positions.
Vector nls_pert
The vector of perturbations for the VMRs of the nonlinear species.
void Adapt(const ArrayOfArrayOfSpeciesTag &current_species, ConstVectorView current_f_grid, const Verbosity &verbosity)
Adapt lookup table to current calculation.
Vector t_ref
The reference temperature profile [K].
void Extract(Matrix &sga, const ArrayOfSpeciesTag &select_abs_species, 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.
ArrayOfIndex nonlinear_species
The species tags with non-linear treatment.
ArrayOfArrayOfSpeciesTag species
The species tags for which the table is valid.
Row-major grid creation.
Definition: grids.h:52
The Matrix class.
Definition: matpackI.h:1261
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1010
The range class.
Definition: matpackI.h:159
The Tensor3View class.
Definition: matpackIII.h:246
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1054
The Tensor5 class.
Definition: matpackV.h:516
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:1741
The Vector class.
Definition: matpackI.h:899
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
Subclasses of runtime_error.
Definition: check_input.h:108
#define DEBUG_ONLY(...)
Definition: debug.h:52
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
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.
ostream & operator<<(ostream &os, const GasAbsLookup &)
Output operatior for GasAbsLookup.
Declarations for the gas absorption lookup table.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Header file for interpolation.cc.
Interpolation::Lagrange LagrangeInterpolation
#define abs(x)
long int flag
#define max(a, b)
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:82
bool is_unique(const ArrayOfIndex &x)
Checks if an ArrayOfIndex is unique, i.e., has no duplicate values.
Definition: logic.cc:277
Header file for logic.cc.
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:1431
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition: messages.h:206
#define CREATE_OUT2
Definition: messages.h:205
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
void reinterp(VectorView out, const ConstVectorView &iy, const Grid< Vector, 1 > &iw, const Array< Lagrange > &dim0)
Array< Lagrange > LagrangeVector(const ConstVectorView &xs, const ConstVectorView &xi, const Index polyorder, const Numeric extrapol, const bool do_derivs, const GridType type, const std::pair< Numeric, Numeric > cycle)
Online Transmission UseSurfaceRtprop TransmitterReceiverPath FreeElectrons
Definition: arts_options.h:114
Implements Zeeman modeling.
Definition: zeemandata.cc:333
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition: physics_funcs.h:70
Index find_first_species(const ArrayOfArrayOfSpeciesTag &specs, Species::Species spec) noexcept
#define d
#define a
#define c
#define b