ARTS 2.5.0 (git: 9ee3ac6c)
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 Index& p_interp_order,
515 const Index& t_interp_order,
516 const Index& h2o_interp_order,
517 const Index& f_interp_order,
518 const Numeric& p,
519 const Numeric& T,
520 ConstVectorView abs_vmrs,
521 ConstVectorView new_f_grid,
522 const Numeric& extpolfac) const {
523 // 1. Obtain some properties of the lookup table:
524
525 // Number of gas species in the table:
526 const Index n_species = species.nelem();
527
528 // Number of nonlinear species:
529 const Index n_nls = nonlinear_species.nelem();
530
531 // Number of frequencies in the table:
532 const Index n_f_grid = f_grid.nelem();
533
534 // Number of pressure grid points in the table:
535 const Index n_p_grid = p_grid.nelem();
536
537 // Number of temperature perturbations:
538 const Index n_t_pert = t_pert.nelem();
539
540 // Number of nonlinear species perturbations:
541 const Index n_nls_pert = nls_pert.nelem();
542
543 // Number of frequencies in new_f_grid, the frequency grid for which we
544 // want to extract.
545 const Index n_new_f_grid = new_f_grid.nelem();
546
547 // 2. First some checks on the lookup table itself:
548
549 // Most checks here are asserts, because they check the internal
550 // consistency of the lookup table. They should never fail if the
551 // table has been created with ARTS.
552
553 // If there are nonlinear species, then at least one species must be
554 // H2O. We will use that to perturb in the case of nonlinear
555 // species.
556 Index h2o_index = -1;
557 if (n_nls > 0) {
558 h2o_index = find_first_species(species, Species::Species::Water);
559
560 // This is a runtime error, even though it would be more logical
561 // for it to be an assertion, since it is an internal check on
562 // the table. The reason is that it is somewhat awkward to check
563 // for this in other places.
564 if (h2o_index == -1) {
565 ostringstream os;
566 os << "With nonlinear species, at least one species must be a H2O species.";
567 throw runtime_error(os.str());
568 }
569 }
570
571 // Check that the dimension of vmrs_ref is consistent with species and p_grid:
572 ARTS_ASSERT(is_size(vmrs_ref, n_species, n_p_grid));
573
574 // Check dimension of t_ref:
575 ARTS_ASSERT(is_size(t_ref, n_p_grid));
576
577 // Check dimension of xsec:
578 DEBUG_ONLY({
579 Index a, b, c, d;
580 if (0 == n_t_pert)
581 a = 1;
582 else
583 a = n_t_pert;
584 b = n_species + n_nls * (n_nls_pert - 1);
585 c = n_f_grid;
586 d = n_p_grid;
587 // cout << "xsec: "
588 // << xsec.nbooks() << ", "
589 // << xsec.npages() << ", "
590 // << xsec.nrows() << ", "
591 // << xsec.ncols() << "\n";
592 // cout << "a b c d: "
593 // << a << ", "
594 // << b << ", "
595 // << c << ", "
596 // << d << "\n";
597 ARTS_ASSERT(is_size(xsec, a, b, c, d));
598 })
599
600 // Make sure that log_p_grid is initialized:
601 if (log_p_grid.nelem() != n_p_grid) {
602 ostringstream os;
603 os << "The lookup table internal variable log_p_grid is not initialized.\n"
604 << "Use the abs_lookupAdapt method!";
605 throw runtime_error(os.str());
606 }
607
608 // Verify that we have enough pressure, temperature,humdity, and frequency grid points
609 // for the desired interpolation orders. This check is not only
610 // table internal, since abs_nls_interp_order and abs_t_interp_order
611 // are separate WSVs that could have been modified. Hence, these are
612 // runtime errors.
613
614 if ((n_p_grid < p_interp_order + 1)) {
615 ostringstream os;
616 os << "The number of pressure grid points in the table (" << n_p_grid
617 << ") is not enough for the desired order of interpolation ("
618 << p_interp_order << ").";
619 throw runtime_error(os.str());
620 }
621
622 if ((n_nls_pert != 0) && (n_nls_pert < h2o_interp_order + 1)) {
623 ostringstream os;
624 os << "The number of humidity perturbation grid points in the table ("
625 << n_nls_pert
626 << ") is not enough for the desired order of interpolation ("
627 << h2o_interp_order << ").";
628 throw runtime_error(os.str());
629 }
630
631 if ((n_t_pert != 0) && (n_t_pert < t_interp_order + 1)) {
632 ostringstream os;
633 os << "The number of temperature perturbation grid points in the table ("
634 << n_t_pert << ") is not enough for the desired order of interpolation ("
635 << t_interp_order << ").";
636 throw runtime_error(os.str());
637 }
638
639 if ((n_f_grid < f_interp_order + 1)) {
640 ostringstream os;
641 os << "The number of frequency grid points in the table (" << n_f_grid
642 << ") is not enough for the desired order of interpolation ("
643 << f_interp_order << ").";
644 throw runtime_error(os.str());
645 }
646
647 // 3. Checks on the input variables:
648
649 // Check that abs_vmrs has the right dimension:
650 if (!is_size(abs_vmrs, n_species)) {
651 ostringstream os;
652 os << "Number of species in lookup table does not match number\n"
653 << "of species for which you want to extract absorption.\n"
654 << "Have you used abs_lookupAdapt? Or did you miss to add\n"
655 << "some VRM fields (e.g. for free electrons or particles)?\n";
656 throw runtime_error(os.str());
657 }
658
659 // 4. Set up some things we will need later on:
660
661 // 4.a Frequency grid positions
662
663 // Frequency grid positions. The pointer is used to save copying of the
664 // default from the lookup table.
667
668 // With f_interp_order 0 the frequency grid has to have the same size as in the
669 // lookup table, or exactly one element. If it matches the lookup table, we
670 // do no frequency interpolation at all. (We set the frequency grid positions
671 // to the predefined ones that come with the lookup table.)
672 if (f_interp_order == 0) {
673
674 // We do some superficial checks below, to make sure that the
675 // frequency grid is the same as in the lookup table (only first
676 // and last element of f_grid are tested). As margin for
677 // agreement, we pick a value that is just slightly smaller than
678 // the perturbation that is used by the wind jacobian, which is 0.1.
679 const Numeric allowed_f_margin = 0.09;
680
681 if (n_new_f_grid == n_f_grid) {
682 // Use the default flag that is stored in the lookup table itself
683 // (which effectively means no frequency interpolation)
685
686 // Check first f_grid element:
687 if (abs(f_grid[0] - new_f_grid[0]) > allowed_f_margin)
688 {
689 ostringstream os;
690 os << "First frequency in f_grid inconsistent with lookup table.\n"
691 << "f_grid[0] = " << f_grid[0] << "\n"
692 << "new_f_grid[0] = " << new_f_grid[0] << ".";
693 throw runtime_error(os.str());
694 }
695
696 // Check last f_grid element:
697 if (abs(f_grid[n_f_grid - 1] - new_f_grid[n_new_f_grid - 1]) > allowed_f_margin)
698 {
699 ostringstream os;
700 os << "Last frequency in f_grid inconsistent with lookup table.\n"
701 << "f_grid[n_f_grid-1] = " << f_grid[n_f_grid - 1]
702 << "\n"
703 << "new_f_grid[n_new_f_grid-1] = " << new_f_grid[n_new_f_grid - 1]
704 << ".";
705 throw runtime_error(os.str());
706 }
707 } else if (n_new_f_grid == 1) {
708 flag = &flag_local;
709 flag_local = Interpolation::LagrangeVector(new_f_grid, f_grid, 0);
710
711 // Check that we really are on a frequency grid point, for safety's sake.
712 if (abs(f_grid[flag_local[0].pos] - new_f_grid[0]) > allowed_f_margin)
713 {
714 ostringstream os;
715 os << "Cannot find a matching lookup table frequency for frequency "
716 << new_f_grid[0] << ".\n"
717 << "(This check has not been properly tested, so perhaps this is\n"
718 << "a false alarm. Check for this in file gas_abs_lookup.cc.)";
719 throw runtime_error(os.str());
720 }
721 } else {
722 ostringstream os;
723 os << "With f_interp_order 0 the frequency grid has to have the same\n"
724 << "size as in the lookup table, or exactly one element.";
725 throw runtime_error(os.str());
726 }
727 } else {
728 const Numeric f_min = f_grid[0] - 0.5 * (f_grid[1] - f_grid[0]);
729 const Numeric f_max = f_grid[n_f_grid - 1] +
730 0.5 * (f_grid[n_f_grid - 1] - f_grid[n_f_grid - 2]);
731 if (new_f_grid[0] < f_min) {
732 ostringstream os;
733 os << "Problem with gas absorption lookup table.\n"
734 << "At least one frequency is outside the range covered by the lookup table.\n"
735 << "Your new frequency value is " << new_f_grid[0] << " Hz.\n"
736 << "The allowed range is " << f_min << " to " << f_max << " Hz.";
737 throw runtime_error(os.str());
738 }
739 if (new_f_grid[n_new_f_grid - 1] > f_max) {
740 ostringstream os;
741 os << "Problem with gas absorption lookup table.\n"
742 << "At least one frequency is outside the range covered by the lookup table.\n"
743 << "Your new frequency value is " << new_f_grid[n_new_f_grid - 1]
744 << " Hz.\n"
745 << "The allowed range is " << f_min << " to " << f_max << " Hz.";
746 throw runtime_error(os.str());
747 }
748
749 // We do have real frequency interpolation (f_interp_order!=0).
750 flag = &flag_local;
751 flag_local = Interpolation::LagrangeVector(new_f_grid, f_grid, f_interp_order);
752 }
753
754 // 4.b Other stuff
755
756 // Flag for temperature interpolation, if this is not 0 we want
757 // to do T interpolation:
758 const Index do_T = n_t_pert;
759
760 // Set up a logical array for the nonlinear species
761 ArrayOfIndex non_linear(n_species, 0);
762 for (Index s = 0; s < n_nls; ++s) {
763 non_linear[nonlinear_species[s]] = 1;
764 }
765
766 // Calculate the number density for the given pressure and
767 // temperature:
768 // n = n0*T0/p0 * p/T or n = p/kB/t, ideal gas law
769 const Numeric n = number_density(p, T);
770
771 // 5. Determine pressure grid position and interpolation weights:
772
773 // Check that p is inside the grid. (p_grid is sorted in decreasing order.)
774 {
775 const Numeric p_max = p_grid[0] + 0.5 * (p_grid[0] - p_grid[1]);
776 const Numeric p_min = p_grid[n_p_grid - 1] -
777 0.5 * (p_grid[n_p_grid - 2] - p_grid[n_p_grid - 1]);
778 if ((p > p_max) || (p < p_min)) {
779 ostringstream os;
780 os << "Problem with gas absorption lookup table.\n"
781 << "Pressure p is outside the range covered by the lookup table.\n"
782 << "Your p value is " << p << " Pa.\n"
783 << "The allowed range is " << p_min << " to " << p_max << ".\n"
784 << "The pressure grid range in the table is " << p_grid[n_p_grid - 1]
785 << " to " << p_grid[0] << ".\n"
786 << "We allow a bit of extrapolation, but NOT SO MUCH!";
787 throw runtime_error(os.str());
788 }
789 }
790
791 // For sure, we need to store the pressure grid position.
792 // We do the interpolation in log(p). Test have shown that this
793 // gives slightly better accuracy than interpolating in p directly.
794 const auto plag = Interpolation::LagrangeVector(log(p), log_p_grid, p_interp_order);
795
796 // Pressure interpolation weights:
797 const auto pitw = interpweights(plag[0]);
798
799 // Define also other grid positions and interpolation weights here, so that
800 // we do not have to allocate them over and over in the loops below.
801
802 // Define the ArrayOfLagrangeInterpolation that corresponds to "no interpolation at all".
803 const ArrayOfLagrangeInterpolation lag_trivial(1);
804
805 // Temperature grid positions.
806 ArrayOfLagrangeInterpolation tlag_withT(1); // Only a scalar.
807 const ArrayOfLagrangeInterpolation* tlag; // Pointer to either tlag_withT or lag_trivial.
808
809 // Set this_t_interp_order, depending on whether we do T interpolation or not.
810 if (do_T) {
811 tlag = &tlag_withT;
812 } else {
813 // For the !do_T case we simply take the single
814 // temperature that is there, so we point tgp accordingly.
815 tlag = &lag_trivial;
816 }
817
818 // H2O(VMR) grid positions. vlag is what will be used in the interpolation.
819 // Depending on species, it is either pointed to lag_trivial, or to vlag_h2o.
821 ArrayOfLagrangeInterpolation vlag_h2o(1); // only a scalar
822
823 // 6. We do the T and VMR interpolation for the pressure levels
824 // that are used in the pressure interpolation. (How many depends on
825 // p_interp_order.)
826
827 // To store the interpolated result for the p_interp_order+1
828 // pressure levels:
829 // xsec dimensions are:
830 // Temperature
831 // H2O
832 // Frequency
833 // Pressure
834 // Dimensions of pre_interpolated are:
835 // Pressure (interpolation points)
836 // Species
837 // Temperature (always 1)
838 // H2O (always 1)
839 // Frequency
840
841 Tensor5 xsec_pre_interpolated;
842 xsec_pre_interpolated.resize(
843 p_interp_order + 1, n_species, 1, 1, n_new_f_grid);
844
845 // Define variables for interpolation weights outside the loops.
846 // We will make itw point to either the weights with H2O interpolation, or
847 // the ones without.
848 Tensor3OfTensor3 itw_withH2O(0,0,0), itw_noH2O(0,0,0);
849 const Tensor3OfTensor3 *itw;
850
851 for (Index pi = 0; pi < p_interp_order + 1; ++pi) {
852 // Throw a runtime error if one of the reference VMR profiles is zero, but
853 // abs_vmrs is not. (This means that the lookup table was calculated with a
854 // reference profile of zero for that gas.)
855 // for (Index si=0; si<n_species; ++si)
856 // if ( (vmrs_ref(si,pi) == 0) &&
857 // (abs_vmrs[si] != 0) )
858 // {
859 // ostringstream os;
860 // os << "Reference VMR profile is zero, you cannot extract\n"
861 // << "Absorption for this species.\n"
862 // << "Species: " << si
863 // << " (" << get_species_name(species[si]) << ")\n"
864 // << "Lookup table pressure level: " << pi
865 // << " (" << p_grid[pi] << " Pa).";
866 // throw runtime_error( os.str() );
867 // }
868
869 // Index into p_grid:
870 const Index this_p_grid_index = plag[0].pos + pi;
871
872 // Determine temperature grid position. This is only done if we
873 // want temperature interpolation, but the variable tgp has to
874 // be visible also outside for later use:
875 if (do_T) {
876 // Temperature in the atmosphere is altitude
877 // dependent. When we do the interpolation for the pressure level
878 // below and above our point, we should correct the target value of
879 // the interpolation to the altitude (pressure) difference. This
880 // ensures that there is for example no T interpolation if the
881 // desired T is right on the reference profile curve.
882 //
883 // I explicitly compared this with the old option to calculate
884 // the temperature offset relative to the temperature at
885 // this level. The performance in both cases is very
886 // similar. The reason, why I decided to keep this new
887 // version, is that it avoids the problem of needing
888 // oversized temperature perturbations if the pressure
889 // grid is coarse.
890 //
891 // No! The above approach leads to problems when combined with
892 // higher order pressure interpolation. The problem is that
893 // the reference T and VMR profiles may be very
894 // irregular. (For example the H2O profile often has a big
895 // jump near the bottom.) That sometimes leads to negative
896 // effective reference values when the reference profile is
897 // interpolated. I therefore reverted back to the original
898 // version of using the real temperature and humidity, not
899 // the interpolated one.
900
901 // const Numeric effective_T_ref = interp(pitw,t_ref,pgp);
902 const Numeric effective_T_ref = t_ref[this_p_grid_index];
903
904 // Convert temperature to offset from t_ref:
905 const Numeric T_offset = T - effective_T_ref;
906
907 // cout << "T_offset = " << T_offset << endl;
908
909 // Check that temperature offset is inside the allowed range.
910 {
911 const Numeric t_min = t_pert[0] - extpolfac * (t_pert[1] - t_pert[0]);
912 const Numeric t_max =
913 t_pert[n_t_pert - 1] +
914 extpolfac * (t_pert[n_t_pert - 1] - t_pert[n_t_pert - 2]);
915 if ((T_offset > t_max) || (T_offset < t_min)) {
916 ostringstream os;
917 os << "Problem with gas absorption lookup table.\n"
918 << "Temperature T is outside the range covered by the lookup table.\n"
919 << "Your temperature was " << T << " K at a pressure of " << p
920 << " Pa.\n"
921 << "The temperature offset value is " << T_offset << ".\n"
922 << "The allowed range is " << t_min << " to " << t_max << ".\n"
923 << "The temperature perturbation grid range in the table is "
924 << t_pert[0] << " to " << t_pert[n_t_pert - 1] << ".\n"
925 << "We allow a bit of extrapolation, but NOT SO MUCH!";
926 throw runtime_error(os.str());
927 }
928 }
929
930 tlag_withT[0] = LagrangeInterpolation(0, T_offset, t_pert, t_interp_order);
931 }
932
933 // Determine the H2O VMR grid position. We need to do this only
934 // once, since the only species who's VMR is interpolated is
935 // H2O. We do this only if there are nonlinear species, but the
936 // variable has to be visible later.
937 if (n_nls > 0) {
938 // Similar to the T case, we first interpolate the reference
939 // VMR to the pressure of extraction, then compare with
940 // the extraction VMR to determine the offset/fractional
941 // difference for the VMR interpolation.
942 //
943 // No! The above approach leads to problems when combined with
944 // higher order pressure interpolation. The problem is that
945 // the reference T and VMR profiles may be very
946 // irregular. (For example the H2O profile often has a big
947 // jump near the bottom.) That sometimes leads to negative
948 // effective reference values when the reference profile is
949 // interpolated. I therefore reverted back to the original
950 // version of using the real temperature and humidity, not
951 // the interpolated one.
952
953 // const Numeric effective_vmr_ref = interp(pitw,
954 // vmrs_ref(h2o_index, Range(joker)),
955 // pgp);
956 const Numeric effective_vmr_ref = vmrs_ref(h2o_index, this_p_grid_index);
957
958 // Fractional VMR:
959 const Numeric VMR_frac = abs_vmrs[h2o_index] / effective_vmr_ref;
960
961 // Check that VMR_frac is inside the allowed range.
962 {
963 // FIXME: This check depends on how I interpolate VMR.
964 const Numeric x_min =
965 nls_pert[0] - extpolfac * (nls_pert[1] - nls_pert[0]);
966 const Numeric x_max =
967 nls_pert[n_nls_pert - 1] +
968 extpolfac * (nls_pert[n_nls_pert - 1] - nls_pert[n_nls_pert - 2]);
969
970 if ((VMR_frac > x_max) || (VMR_frac < x_min)) {
971 ostringstream os;
972 os << "Problem with gas absorption lookup table.\n"
973 << "VMR for H2O (species " << h2o_index
974 << ") is outside the range covered by the lookup table.\n"
975 << "Your VMR was " << abs_vmrs[h2o_index] << " at a pressure of "
976 << p << " Pa.\n"
977 << "The reference VMR value there is " << effective_vmr_ref << "\n"
978 << "The fractional VMR relative to the reference value is "
979 << VMR_frac << ".\n"
980 << "The allowed range is " << x_min << " to " << x_max << ".\n"
981 << "The fractional VMR perturbation grid range in the table is "
982 << nls_pert[0] << " to " << nls_pert[n_nls_pert - 1] << ".\n"
983 << "We allow a bit of extrapolation, but NOT SO MUCH!";
984 throw runtime_error(os.str());
985 }
986 }
987
988 // For now, do linear interpolation in the fractional VMR.
989 vlag_h2o[0] = LagrangeInterpolation(0, VMR_frac, nls_pert, h2o_interp_order);
990 }
991
992 // Precalculate interpolation weights.
993 if (n_nls < n_species) {
994 // Precalculate weights without H2O interpolation if there are less
995 // nonlinear species than total species. (So at least one species
996 // without H2O interpolation.)
997 itw_noH2O = interpweights(*tlag, lag_trivial, *flag);
998 }
999 if (n_nls > 0) {
1000 // Precalculate weights with H2O interpolation if there is at least
1001 // one nonlinear species.
1002 itw_withH2O = interpweights(*tlag, vlag_h2o, *flag);
1003 }
1004
1005 // 7. Loop species:
1006 Index fpi = 0;
1007 for (Index si = 0; si < n_species; ++si) {
1008 // Flag for VMR interpolation, if this is not 0 we want to
1009 // do VMR interpolation:
1010 const Index do_VMR = non_linear[si];
1011
1012 // For interpolation result.
1013 // Fixed pressure level and species.
1014 // Free dimension is T, H2O, and frequency.
1015 Tensor3View res(xsec_pre_interpolated(
1016 pi, si, Range(joker), Range(joker), Range(joker)));
1017
1018 // Ignore species such as Zeeman and free_electrons which are not
1019 // stored in the lookup table. For those the result is set to 0.
1020 if (species[si].Zeeman() or species[si].FreeElectrons() or species[si].Particles()) {
1021 if (do_VMR) {
1022 ostringstream os;
1023 os << "Problem with gas absorption lookup table.\n"
1024 << "VMR interpolation is not allowed for species \""
1025 << species[si][0].Name() << "\"";
1026 throw runtime_error(os.str());
1027 }
1028 res = 0.;
1029 fpi++;
1030 continue;
1031 }
1032
1033 // Set h2o related interpolation parameters:
1034 Index this_h2o_extent; // Range of H2O interpolation
1035 if (do_VMR) {
1036 vlag = &vlag_h2o;
1037 this_h2o_extent = n_nls_pert;
1038 itw = &itw_withH2O;
1039 } else {
1040 vlag = &lag_trivial;
1041 this_h2o_extent = 1;
1042 itw = &itw_noH2O;
1043 }
1044
1045 // Get the right view on xsec.
1046 ConstTensor3View this_xsec =
1047 xsec(Range(joker), // Temperature range
1048 Range(fpi, this_h2o_extent), // VMR profile range
1049 Range(joker), // Frequency range
1050 this_p_grid_index); // Pressure index
1051
1052 // Do interpolation.
1053 reinterp(res, // result
1054 this_xsec, // input
1055 *itw, // weights
1056 *tlag,
1057 *vlag,
1058 *flag); // grid positions
1059
1060 // Increase fpi. fpi marks the position of the first profile
1061 // of the current species in xsec. This is needed to find
1062 // the right subsection of xsec in the presence of nonlinear species.
1063 if (do_VMR)
1064 fpi += n_nls_pert;
1065 else
1066 fpi++;
1067
1068 } // End of species loop
1069
1070 // fpi should have reached the end of that dimension of xsec. Check
1071 // this with an assertion:
1072 ARTS_ASSERT(fpi == xsec.npages());
1073
1074 } // End of pressure index loop (below and above gp)
1075
1076 // Now we have to interpolate between the p_interp_order+1 pressure levels
1077
1078 // It is a "red" 1D interpolation case we are talking about here.
1079 // (But for a matrix in frequency and species.) Doing a loop over
1080 // frequency and species with an interp call inside would be
1081 // unefficient, so we do this by hand here.
1082 sga.resize(n_species, n_new_f_grid);
1083 sga = 0;
1084 for (Index pi = 0; pi < p_interp_order + 1; ++pi) {
1085 // Multiply pre interpolated quantities with pressure interpolation weights.
1086 // Dimensions of pre_interpolated are:
1087 // Pressure (interpolation points)
1088 // Species
1089 // Temperature (always 1)
1090 // H2O (always 1)
1091 // Frequency
1092 xsec_pre_interpolated(
1093 pi, Range(joker), Range(joker), Range(joker), Range(joker)) *= pitw[pi];
1094
1095 // Add up in sga.
1096 // Dimensions of sga are (species, frequency)
1097 sga += xsec_pre_interpolated(pi, Range(joker), 0, 0, Range(joker));
1098 }
1099
1100 // Watch out, this is not yet the final result, we
1101 // need to multiply with the number density of the species, i.e.,
1102 // with the total number density n, times the VMR of the
1103 // species:
1104 for (Index si = 0; si < n_species; ++si)
1105 sga(si, Range(joker)) *= (n * abs_vmrs[si]);
1106
1107 // That's it, we're done!
1108}
1109
1110const Vector& GasAbsLookup::GetFgrid() const { return f_grid; }
1111
1112const Vector& GasAbsLookup::GetPgrid() const { return p_grid; }
1113
1115ostream& operator<<(ostream& os, const GasAbsLookup& /* gal */) {
1116 os << "GasAbsLookup: Output operator not implemented";
1117 return os;
1118}
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
Number of elements.
Definition: array.h:195
A constant view of a Tensor3.
Definition: matpackIII.h:132
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:58
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:55
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:64
A constant view of a Vector.
Definition: matpackI.h:489
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
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.
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.
Vector t_ref
The reference temperature profile [K].
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:1225
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
The range class.
Definition: matpackI.h:165
The Tensor3View class.
Definition: matpackIII.h:239
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1058
The Tensor5 class.
Definition: matpackV.h:506
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:1741
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
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:79
bool is_unique(const ArrayOfIndex &x)
Checks if an ArrayOfIndex is unique, i.e., has no duplicate values.
Definition: logic.cc:274
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:1472
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition: messages.h:207
#define CREATE_OUT2
Definition: messages.h:206
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)
Implements Zeeman modeling.
Definition: zeemandata.cc:281
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
#define d
#define a
#define c
#define b
Index find_first_species(const ArrayOfArrayOfSpeciesTag &specs, Species::Species spec) noexcept