19#include "matpack_concepts.h"
41 const ConstVectorView& f_grid,
42 const Numeric& temperature,
44 const Numeric& T_extrapolfac,
49 const Index nf = f_grid.nelem();
58 if (out3.sufficient_priority()) {
61 os <<
" f_grid: " << f_grid[0] <<
" - " << f_grid[nf - 1] <<
" Hz\n"
62 <<
" data_f_grid: " << data_f_grid[0] <<
" - "
63 << data_f_grid[data_f_grid.nelem() - 1] <<
" Hz\n"
64 <<
" temperature: " << temperature <<
" K\n"
65 <<
" data_T_grid: " << data_T_grid[0] <<
" - "
66 << data_T_grid[data_T_grid.nelem() - 1] <<
" K\n";
77 Index i_fstart, i_fstop;
79 for (i_fstart = 0; i_fstart < nf; ++i_fstart)
80 if (f_grid[i_fstart] >= data_f_grid[0])
break;
83 if (i_fstart == nf)
return;
85 for (i_fstop = nf - 1; i_fstop >= 0; --i_fstop)
86 if (f_grid[i_fstop] <= data_f_grid[data_f_grid.nelem() - 1])
break;
89 if (i_fstop == -1)
return;
92 const Index f_extent = i_fstop - i_fstart + 1;
94 if (out3.sufficient_priority()) {
96 os <<
" " << f_extent <<
" frequency extraction points starting at "
97 <<
"frequency index " << i_fstart <<
".\n";
104 if (f_extent < 1)
return;
107 ConstVectorView f_grid_active = f_grid[Range(i_fstart, f_extent)];
110 VectorView result_active = result[Range(i_fstart, f_extent)];
113 constexpr Index f_order = 3;
117 if (data_f_grid.nelem() < f_order + 1) {
119 os <<
"Not enough frequency grid points in CIA data.\n"
120 <<
"You have only " << data_f_grid.nelem() <<
" grid points.\n"
121 <<
"But need at least " << f_order + 1 <<
".";
122 throw runtime_error(os.str());
128 switch (data_T_grid.nelem()) {
157 }
catch (
const std::runtime_error& e) {
165 throw runtime_error(e.what());
171 const auto f_lag = my_interp::lagrange_interpolation_list<FixedLagrangeInterpolation<f_order>>(f_grid_active, data_f_grid);
179 const auto Tnew = matpack::matpack_constant_data<Numeric, 1>{temperature};
181 const auto T_lag = my_interp::lagrange_interpolation_list<FixedLagrangeInterpolation<1>>(Tnew, data_T_grid);
182 result_active = reinterp(cia_data.
data,
interpweights(f_lag, T_lag), f_lag, T_lag).reduce_rank<0>();
183 }
else if (T_order == 2) {
184 const auto T_lag = my_interp::lagrange_interpolation_list<FixedLagrangeInterpolation<2>>(Tnew, data_T_grid);
185 result_active = reinterp(cia_data.
data,
interpweights(f_lag, T_lag), f_lag, T_lag).reduce_rank<0>();
186 }
else if (T_order == 3) {
187 const auto T_lag = my_interp::lagrange_interpolation_list<FixedLagrangeInterpolation<3>>(Tnew, data_T_grid);
188 result_active = reinterp(cia_data.
data,
interpweights(f_lag, T_lag), f_lag, T_lag).reduce_rank<0>();
190 throw std::runtime_error(
"Cannot have this T_order, you must update the code...");
198 for (Index i = 0; i < result_active.nelem(); ++i)
199 if (result_active[i] < 0) result_active[i] = 0;
213 const Species::Species sp1,
214 const Species::Species sp2) {
215 for (Index i = 0; i < cia_data.
nelem(); i++)
216 if ((cia_data[i].
Species(0) == sp1 && cia_data[i].
Species(1) == sp2) ||
225 const Numeric &temperature,
226 const Numeric &T_extrapolfac,
const Index &robust,
230 Vector result(res.nelem());
231 for (
auto &this_cia :
mdata) {
246 return Species::toShortName(
mspecies[i]);
256 Species::Species spec_ind = Species::fromShortName(name);
261 "Species does not exist in ARTS: ", name)
281 out2 <<
" Reading file: " << filename <<
"\n";
288 Numeric wave_min = -1.;
289 Numeric wave_max = -1.;
299 std::vector<Numeric> temp;
315 if (is.eof())
continue;
317 if (line.
nelem() < 100) {
319 os <<
"Error in line " << nline <<
" reading CIA catalog file "
321 <<
"Header line unexpectedly short: " << endl
324 throw runtime_error(os.str());
329 os <<
"Error in line " << nline <<
" reading CIA catalog file "
332 throw runtime_error(os.str());
339 Numeric set_temp = NAN;
340 Numeric set_wave_min = NAN;
341 Numeric set_wave_max = NAN;
349 if (!istr || std::isnan(set_temp) || std::isnan(set_wave_min) ||
350 std::isnan(set_wave_max)) {
352 os <<
"Error in line " << nline <<
" reading CIA catalog file "
355 throw runtime_error(os.str());
360 if (npoints == -1 || wave_min != set_wave_min || wave_max != set_wave_max) {
363 npoints = set_npoints;
365 wave_min = set_wave_min;
366 wave_max = set_wave_max;
368 freq.resize(set_npoints);
372 if (npoints != set_npoints) {
374 os <<
"Error in line " << nline <<
" reading CIA catalog file "
376 <<
"Inconsistent number of data points. Expected " << npoints
377 <<
", got " << set_npoints;
379 throw runtime_error(os.str());
382 temp.push_back(set_temp);
383 cia.push_back(Vector(npoints));
387 for (Index i = 0; i < npoints; i++) {
398 if (std::isnan(
w) || std::isnan(
c) || is.bad() || istr.bad()) {
400 os <<
"Error in line " << nline <<
" reading CIA catalog file "
401 << filename <<
":" << endl
404 throw runtime_error(os.str());
412 cia[nset][i] =
c / 1e10;
420 os <<
"Error in line " << nline <<
" reading CIA catalog file " << filename
423 throw runtime_error(os.str());
440 const ArrayOfVector& cia) {
442 dataset.
resize(freq.nelem(), temp.nelem());
451 for (Index t = 0; t < temp.nelem(); t++) dataset.
data(joker, t) = cia[t];
452 mdata.push_back(dataset);
471 os <<
"CIARecord output operator not yet implemented." << endl;
Declarations required for the calculation of absorption coefficients.
Constants of physical expressions as constexpr.
Index cia_get_index(const ArrayOfCIARecord &cia_data, const Species::Species sp1, const Species::Species sp2)
Get the index in cia_data for the two given species.
constexpr Numeric SPEED_OF_LIGHT
void cia_interpolation(VectorView result, const ConstVectorView &f_grid, const Numeric &temperature, const GriddedField2 &cia_data, const Numeric &T_extrapolfac, const Index &robust, const Verbosity &verbosity)
Interpolate CIA data.
ostream & operator<<(ostream &os, const CIARecord &)
Output operator for CIARecord.
Header file for work with HITRAN collision induced absorption (CIA).
void cia_interpolation(VectorView result, const ConstVectorView &frequency, const Numeric &temperature, const GriddedField2 &cia_data, const Numeric &T_extrapolfac, const Index &robust, const Verbosity &verbosity)
Interpolate CIA data.
This can be used to make arrays out of anything.
Index nelem() const ARTS_NOEXCEPT
CIA data for a single pair of molecules.
std::array< Species::Species, 2 > mspecies
The pair of molecules associated with these CIA data.
ArrayOfGriddedField2 mdata
The data itself, directly from the HITRAN file.
void AppendDataset(const CIARecord &c2)
Append other CIARecord to this.
void ReadFromCIA(const String &filename, const Verbosity &verbosity)
Read CIA catalog file.
void SetMoleculeName(const Index i, const String &name)
Set each molecule name (from a string) that is associated with this CIARecord.
void Extract(VectorView result, const ConstVectorView &f_grid, const Numeric &temperature, const Numeric &T_extrapolfac, const Index &robust, const Verbosity &verbosity) const
Vector version of extract.
Index DatasetCount() const
Return number of datasets in this record.
String MoleculeName(const Index i) const
Return each molecule name (as a string) that is associated with this CIARecord.
const GriddedField2 & Dataset(Index dataset) const
Return CIA dataset.
void resize(const GriddedField2 &gf)
Make this GriddedField2 the same size as the given one.
void set_grid_name(Index i, const String &s)
Set grid name.
void set_grid(Index i, const Vector &g)
Set a numeric grid.
const Vector & get_numeric_grid(Index i) const
Get a numeric grid.
Input manipulator class for doubles to enable nan and inf parsing.
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR_IF(condition,...)
constexpr bool good_enum(EnumType x) noexcept
Checks if the enum number is good.
void open_input_file(ifstream &file, const std::string_view name)
Open a file for reading.
This file contains basic functions to handle ASCII files.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...