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