16#include "matpack_data.h"
33 ConstVectorView old_grid,
34 ConstVectorView new_grid,
38 const Index n_new_grid = new_grid.nelem();
39 const Index n_old_grid = old_grid.nelem();
48 for (Index i = 0; i < n_new_grid; ++i) {
53 while (abs(new_grid[i] - old_grid[j]) >
54 max(abs(new_grid[i]), abs(old_grid[j])) * DBL_EPSILON) {
56 if (j >= n_old_grid) {
58 os <<
"Cannot find new frequency " << i <<
" (" << new_grid[i]
59 <<
"Hz) in the lookup table frequency grid.";
60 throw runtime_error(os.str());
65 out3 <<
" " << new_grid[i] <<
" found, index = " << pos[i] <<
".\n";
103 ConstVectorView current_f_grid,
109 const Index n_current_species = current_species.
nelem();
110 const Index n_current_f_grid = current_f_grid.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();
118 out2 <<
" Original table: " << n_species <<
" species, " << n_f_grid
120 <<
" Adapt to: " << n_current_species <<
" species, "
121 << n_current_f_grid <<
" frequencies.\n";
124 out2 <<
" Table contains no nonlinear species.\n";
129 for (Index s = 0; s < n_nls; ++s) {
133 if (0 ==
t_pert.nelem()) {
134 out2 <<
" Table contains no temperature perturbations.\n";
146 if (0 == n_species) {
148 os <<
"The lookup table should have at least one species.";
149 throw runtime_error(os.str());
156 os <<
"The table must not have duplicate nonlinear species.\n"
158 throw runtime_error(os.str());
162 for (Index i = 0; i < n_nls; ++i) {
164 os <<
"nonlinear_species[" << i <<
"]";
190 if (0 == n_nls_pert) {
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());
206 if (0 ==
t_pert.nelem()) {
213 chk_size(
"xsec",
xsec, 1, n_species, n_f_grid, n_p_grid);
231 Index
b = n_species + n_nls * (n_nls_pert - 1);
242 for (Index i = 0, sp = 0; i < n_species; ++i) {
243 original_spec_pos_in_xsec[i] = sp;
253 if (0 == n_current_species) {
255 os <<
"The list of current species should not be empty.";
256 throw runtime_error(os.str());
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() <<
": ";
271 i_current_species[i] =
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) {
281 i_current_species[i] = -1;
284 os <<
"Species " << current_species[i].Name()
285 <<
" is missing in absorption lookup table.";
286 throw runtime_error(os.str());
290 out3 <<
"found (or trivial).\n";
294 Index n_current_nonlinear_species = 0;
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)
305 if (non_linear[i_current_species[i]]) {
306 out3 <<
" " << current_species[i] <<
"\n";
308 current_non_linear[i] = 1;
309 ++n_current_nonlinear_species;
322 out3 <<
" Looking for Frequencies in lookup table:\n";
328 i_current_f_grid,
f_grid, current_f_grid, verbosity);
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) {
345 new_table.
species[i] = current_species[i];
350 new_table.
f_grid.resize(n_current_f_grid);
351 for (Index i = 0; i < n_current_f_grid; ++i) {
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));
368 new_table.
vmrs_ref(i, Range(joker)) = NAN;
388 new_table.
xsec.resize(
390 n_current_species + n_current_nonlinear_species * (n_nls_pert - 1),
399 for (Index i_s = 0, sp = 0; i_s < n_current_species; ++i_s) {
402 if (current_non_linear[i_s])
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)) =
415 Range(original_spec_pos_in_xsec[i_current_species[i_s]], n_v),
416 i_current_f_grid[i_f],
421 new_table.
xsec(Range(joker), Range(sp, n_v), i_f, Range(joker)) = NAN;
500 const Index& p_interp_order,
501 const Index& t_interp_order,
502 const Index& h2o_interp_order,
503 const Index& f_interp_order,
506 ConstVectorView abs_vmrs,
507 ConstVectorView new_f_grid,
508 const Numeric& extpolfac)
const {
518 const Index n_f_grid =
f_grid.nelem();
521 const Index n_p_grid =
p_grid.nelem();
524 const Index n_t_pert =
t_pert.nelem();
527 const Index n_nls_pert =
nls_pert.nelem();
531 const Index n_new_f_grid = new_f_grid.nelem();
542 Index h2o_index = -1;
550 if (h2o_index == -1) {
552 os <<
"With nonlinear species, at least one species must be a H2O species.";
553 throw runtime_error(os.str());
570 b = n_species + n_nls * (n_nls_pert - 1);
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());
600 if ((n_p_grid < p_interp_order + 1)) {
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());
608 if ((n_nls_pert != 0) && (n_nls_pert < h2o_interp_order + 1)) {
610 os <<
"The number of humidity perturbation grid points in the table ("
612 <<
") is not enough for the desired order of interpolation ("
613 << h2o_interp_order <<
").";
614 throw runtime_error(os.str());
617 if ((n_t_pert != 0) && (n_t_pert < t_interp_order + 1)) {
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());
625 if ((n_f_grid < f_interp_order + 1)) {
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());
636 if (!is_size(abs_vmrs, n_species)) {
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());
651 const ArrayOfLagrangeInterpolation* flag;
652 ArrayOfLagrangeInterpolation flag_local;
658 if (f_interp_order == 0) {
665 const Numeric allowed_f_margin = 0.09;
667 if (n_new_f_grid == n_f_grid) {
673 if (abs(
f_grid[0] - new_f_grid[0]) > allowed_f_margin)
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());
683 if (abs(
f_grid[n_f_grid - 1] - new_f_grid[n_new_f_grid - 1]) > allowed_f_margin)
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]
689 <<
"new_f_grid[n_new_f_grid-1] = " << new_f_grid[n_new_f_grid - 1]
691 throw runtime_error(os.str());
693 }
else if (n_new_f_grid == 1) {
695 flag_local = my_interp::lagrange_interpolation_list<LagrangeInterpolation>(new_f_grid,
f_grid, 0);
698 if (abs(
f_grid[flag_local[0].pos] - new_f_grid[0]) > allowed_f_margin)
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());
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());
715 const Numeric f_max =
f_grid[n_f_grid - 1] +
717 if (new_f_grid[0] < f_min) {
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());
725 if (new_f_grid[n_new_f_grid - 1] > f_max) {
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]
731 <<
"The allowed range is " << f_min <<
" to " << f_max <<
" Hz.";
732 throw runtime_error(os.str());
737 flag_local = my_interp::lagrange_interpolation_list<LagrangeInterpolation>(new_f_grid,
f_grid, f_interp_order);
744 const Index do_T = n_t_pert;
748 for (Index s = 0; s < n_nls; ++s) {
762 const Numeric p_min =
p_grid[n_p_grid - 1] -
764 if ((p > p_max) || (p < p_min)) {
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());
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);
791 const ArrayOfLagrangeInterpolation lag_trivial(1);
794 ArrayOfLagrangeInterpolation tlag_withT(1);
795 const ArrayOfLagrangeInterpolation* tlag;
808 const ArrayOfLagrangeInterpolation* vlag;
809 ArrayOfLagrangeInterpolation vlag_h2o(1);
829 Tensor5 xsec_pre_interpolated;
830 xsec_pre_interpolated.resize(
831 p_interp_order + 1, n_species, 1, 1, n_new_f_grid);
836 Tensor6 itw_withH2O{}, itw_noH2O{};
839 for (Index pi = 0; pi < p_interp_order + 1; ++pi) {
858 const Index this_p_grid_index = plag[0].pos + pi;
890 const Numeric effective_T_ref =
t_ref[this_p_grid_index];
893 const Numeric T_offset = T - effective_T_ref;
900 const Numeric t_max =
902 extpolfac * (
t_pert[n_t_pert - 1] -
t_pert[n_t_pert - 2]);
903 if ((T_offset > t_max) || (T_offset < t_min)) {
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
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 "
913 <<
"We allow a bit of extrapolation, but NOT SO MUCH!";
914 throw runtime_error(os.str());
918 tlag_withT[0] = LagrangeInterpolation(0, T_offset,
t_pert, t_interp_order);
944 const Numeric effective_vmr_ref =
vmrs_ref(h2o_index, this_p_grid_index);
947 const Numeric VMR_frac = abs_vmrs[h2o_index] / effective_vmr_ref;
952 const Numeric x_min =
954 const Numeric x_max =
958 if ((VMR_frac > x_max) || (VMR_frac < x_min)) {
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 "
965 <<
"The reference VMR value there is " << effective_vmr_ref <<
"\n"
966 <<
"The fractional VMR relative to the reference value is "
968 <<
"The allowed range is " << x_min <<
" to " << x_max <<
".\n"
969 <<
"The fractional VMR perturbation grid range in the table is "
971 <<
"We allow a bit of extrapolation, but NOT SO MUCH!";
972 throw runtime_error(os.str());
977 vlag_h2o[0] = LagrangeInterpolation(0, VMR_frac,
nls_pert, h2o_interp_order);
981 if (n_nls < n_species) {
995 for (Index si = 0; si < n_species; ++si) {
998 const Index do_VMR = non_linear[si];
1003 Tensor3View res(xsec_pre_interpolated(
1004 pi, si, Range(joker), Range(joker), Range(joker)));
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());
1022 Index this_h2o_extent;
1025 this_h2o_extent = n_nls_pert;
1028 vlag = &lag_trivial;
1029 this_h2o_extent = 1;
1034 ConstTensor3View this_xsec =
1036 Range(fpi, this_h2o_extent),
1070 sga.resize(n_species, n_new_f_grid);
1072 for (Index pi = 0; pi < p_interp_order + 1; ++pi) {
1080 xsec_pre_interpolated(
1081 pi, Range(joker), Range(joker), Range(joker), Range(joker)) *= pitw[pi];
1085 sga += xsec_pre_interpolated(pi, Range(joker), 0, 0, Range(joker));
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]);
1097 sga(si, Range(joker)) = 0.;
1099 sga(si, Range(joker)) *= (n * abs_vmrs[si]);
1112 os <<
"GasAbsLookup: Output operator not implemented";
base max(const Array< base > &x)
Max function.
Index nelem() const ARTS_NOEXCEPT
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 ¤t_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.
#define ARTS_ASSERT(condition,...)
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.
Implements Zeeman modeling.
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density