ARTS 2.5.11 (git: 6827797f)
cia.cc
Go to the documentation of this file.
1
14#include "arts_constants.h"
15#include "check_input.h"
16#include "cia.h"
17#include <cmath>
18#include <numeric>
19#include "matpack_concepts.h"
20#include "species_tags.h"
21#include "absorption.h"
22#include "file.h"
23#include "interp.h"
24
26
40void cia_interpolation(VectorView result,
41 const ConstVectorView& f_grid,
42 const Numeric& temperature,
43 const GriddedField2& cia_data,
44 const Numeric& T_extrapolfac,
45 const Index& robust,
46 const Verbosity& verbosity) {
48
49 const Index nf = f_grid.nelem();
50
51 // Assert that result vector has right size:
52 ARTS_ASSERT(result.nelem() == nf);
53
54 // Get data grids:
55 ConstVectorView data_f_grid = cia_data.get_numeric_grid(0);
56 ConstVectorView data_T_grid = cia_data.get_numeric_grid(1);
57
58 if (out3.sufficient_priority()) {
59 // Some detailed information to the most verbose output stream:
60 ostringstream os;
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";
67 out3 << os.str();
68 }
69
70 // Initialize result to zero (important for those frequencies outside the data grid).
71 result = 0;
72
73 // We want to return result zero for all f_grid points that are outside the
74 // data_f_grid, because some CIA datasets are defined only where the absorption
75 // is not zero. So, we have to find out which part of f_grid is inside
76 // data_f_grid.
77 Index i_fstart, i_fstop;
78
79 for (i_fstart = 0; i_fstart < nf; ++i_fstart)
80 if (f_grid[i_fstart] >= data_f_grid[0]) break;
81
82 // Return directly if all frequencies are below data_f_grid:
83 if (i_fstart == nf) return;
84
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;
87
88 // Return directly if all frequencies are above data_f_grid:
89 if (i_fstop == -1) return;
90
91 // Extent for active frequency vector:
92 const Index f_extent = i_fstop - i_fstart + 1;
93
94 if (out3.sufficient_priority()) {
95 ostringstream os;
96 os << " " << f_extent << " frequency extraction points starting at "
97 << "frequency index " << i_fstart << ".\n";
98 out3 << os.str();
99 }
100
101 // If f_extent is less than one, then the entire data_f_grid is between two
102 // grid points of f_grid. (So that we do not have any f_grid points inside
103 // data_f_grid.) Return also in this case.
104 if (f_extent < 1) return;
105
106 // This is the part of f_grid for which we have to do the interpolation.
107 ConstVectorView f_grid_active = f_grid[Range(i_fstart, f_extent)];
108
109 // We have to create a matching view on the result vector:
110 VectorView result_active = result[Range(i_fstart, f_extent)];
111
112 // Decide on interpolation orders:
113 constexpr Index f_order = 3;
114
115 // The frequency grid has to have enough points for this interpolation
116 // order, otherwise throw a runtime error.
117 if (data_f_grid.nelem() < f_order + 1) {
118 ostringstream os;
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());
123 }
124
125 // For T we have to be adaptive, since sometimes there is only one T in
126 // the data
127 Index T_order;
128 switch (data_T_grid.nelem()) {
129 case 1:
130 T_order = 0;
131 break;
132 case 2:
133 T_order = 1;
134 break;
135 case 3:
136 T_order = 2;
137 break;
138 default:
139 T_order = 3;
140 break;
141 }
142
143 // Check if frequency is inside the range covered by the data:
144 chk_interpolation_grids("Frequency interpolation for CIA continuum",
145 data_f_grid,
146 f_grid_active,
147 f_order);
148
149 // Check if temperature is inside the range covered by the data:
150 if (T_order > 0) {
151 try {
152 chk_interpolation_grids("Temperature interpolation for CIA continuum",
153 data_T_grid,
154 temperature,
155 T_order,
156 T_extrapolfac);
157 } catch (const std::runtime_error& e) {
158 // cout << "Gotcha!\n";
159 if (robust) {
160 // Just return NANs, but continue.
161 result_active = NAN;
162 return;
163 } else {
164 // Re-throw the exception.
165 throw runtime_error(e.what());
166 }
167 }
168 }
169
170 // Find frequency grid positions:
171 const auto f_lag = my_interp::lagrange_interpolation_list<FixedLagrangeInterpolation<f_order>>(f_grid_active, data_f_grid);
172
173 // Do the rest of the interpolation.
174 if (T_order == 0) {
175 // No temperature interpolation in this case, just a frequency interpolation.
176 result_active = reinterp(cia_data.data(joker, 0), interpweights(f_lag), f_lag);
177 } else {
178 // Temperature and frequency interpolation.
179 const auto Tnew = matpack::matpack_constant_data<Numeric, 1>{temperature};
180 if (T_order == 1) {
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>();
189 } else {
190 throw std::runtime_error("Cannot have this T_order, you must update the code...");
191 }
192 }
193
194 // cout << "result_active before: " << result_active << endl;
195
196 // Set negative values to zero. (These could happen due to overshooting
197 // of the higher order interpolation.)
198 for (Index i = 0; i < result_active.nelem(); ++i)
199 if (result_active[i] < 0) result_active[i] = 0;
200
201 // cout << "result_active after: " << result_active << endl;
202}
203
212Index cia_get_index(const ArrayOfCIARecord& cia_data,
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) ||
217 (cia_data[i].Species(0) == sp2 && cia_data[i].Species(1) == sp1))
218 return i;
219
220 return -1;
221}
222
223// Documentation in header file.
224void CIARecord::Extract(VectorView res, const ConstVectorView &f_grid,
225 const Numeric &temperature,
226 const Numeric &T_extrapolfac, const Index &robust,
227 const Verbosity &verbosity) const {
228 res = 0;
229
230 Vector result(res.nelem());
231 for (auto &this_cia : mdata) {
232 cia_interpolation(result, f_grid, temperature, this_cia, T_extrapolfac,
233 robust, verbosity);
234 res += result;
235 }
236}
237
238// Documentation in header file.
239String CIARecord::MoleculeName(const Index i) const {
240 // Assert that i is 0 or 1:
241 ARTS_ASSERT(i >= 0);
242 ARTS_ASSERT(i <= 1);
243
244 // The function species_name_from_species_index internally does an assertion
245 // that the species with this index really exists.
246 return Species::toShortName(mspecies[i]);
247}
248
249// Documentation in header file.
250void CIARecord::SetMoleculeName(const Index i, const String& name) {
251 // Assert that i is 0 or 1:
252 ARTS_ASSERT(i >= 0);
253 ARTS_ASSERT(i <= 1);
254
255 // Find out the species index for name:
256 Species::Species spec_ind = Species::fromShortName(name);
257
258 // Function species_index_from_species_name returns -1 if the species does
259 // not exist. Check this:
260 ARTS_USER_ERROR_IF(not good_enum(spec_ind),
261 "Species does not exist in ARTS: ", name)
262
263 // Assign species:
264 mspecies[i] = spec_ind;
265}
266
268
275void CIARecord::ReadFromCIA(const String& filename,
276 const Verbosity& verbosity) {
278
279 ifstream is;
280
281 out2 << " Reading file: " << filename << "\n";
282 open_input_file(is, filename);
283
284 // Number of points for spectral range in current dataset
285 Index npoints = -1;
286
287 // Min/max wave numbers for this dataset
288 Numeric wave_min = -1.;
289 Numeric wave_max = -1.;
290
291 // Current dataset index
292 Index ndataset = -1;
293
294 // Current set in dataset
295 Index nset = -1;
296
297 // Frequency, temp and cia values for current dataset
298 Vector freq;
299 std::vector<Numeric> temp;
300 ArrayOfVector cia;
301
302 // Keep track of current line in file
303 Index nline = 0;
304
305 mdata.resize(0);
306 istringstream istr;
307
308 while (is) {
309 String line;
310
311 // Extract needed information from dataset header line
313 nline++;
314 getline(is, line);
315 if (is.eof()) continue;
316
317 if (line.nelem() < 100) {
318 ostringstream os;
319 os << "Error in line " << nline << " reading CIA catalog file "
320 << filename << endl
321 << "Header line unexpectedly short: " << endl
322 << line;
323
324 throw runtime_error(os.str());
325 }
326
327 if (is.bad()) {
328 ostringstream os;
329 os << "Error in line " << nline << " reading CIA catalog file "
330 << filename << endl;
331
332 throw runtime_error(os.str());
333 }
334
335 line.erase(0, 20);
336
337 // Data for current set
338 Index set_npoints;
339 Numeric set_temp = NAN;
340 Numeric set_wave_min = NAN;
341 Numeric set_wave_max = NAN;
342
343 istr.str(line);
344 istr.clear();
345 istr >> double_imanip() >> set_wave_min >> set_wave_max;
346 istr >> set_npoints;
347 istr >> double_imanip() >> set_temp;
348
349 if (!istr || std::isnan(set_temp) || std::isnan(set_wave_min) ||
350 std::isnan(set_wave_max)) {
351 ostringstream os;
352 os << "Error in line " << nline << " reading CIA catalog file "
353 << filename << endl;
354
355 throw runtime_error(os.str());
356 }
357
358 // If the min/max wave numbers of this set are different from the
359 // previous one, a new dataset starts
360 if (npoints == -1 || wave_min != set_wave_min || wave_max != set_wave_max) {
361 if (ndataset != -1) AppendDataset(freq, Vector{temp}, cia);
362
363 npoints = set_npoints;
364 ndataset++;
365 wave_min = set_wave_min;
366 wave_max = set_wave_max;
367 nset = 0;
368 freq.resize(set_npoints);
369 temp.resize(0);
370 cia.resize(0);
371 }
372 if (npoints != set_npoints) {
373 ostringstream os;
374 os << "Error in line " << nline << " reading CIA catalog file "
375 << filename << endl
376 << "Inconsistent number of data points. Expected " << npoints
377 << ", got " << set_npoints;
378
379 throw runtime_error(os.str());
380 }
381
382 temp.push_back(set_temp);
383 cia.push_back(Vector(npoints));
384
385 // Read dataset
387 for (Index i = 0; i < npoints; i++) {
388 Numeric w, c;
389
390 nline++;
391 getline(is, line);
392
393 w = c = NAN;
394 istr.str(line);
395 istr.clear();
396 istr >> w >> c;
397
398 if (std::isnan(w) || std::isnan(c) || is.bad() || istr.bad()) {
399 ostringstream os;
400 os << "Error in line " << nline << " reading CIA catalog file "
401 << filename << ":" << endl
402 << line;
403
404 throw runtime_error(os.str());
405 }
406
407 // Convert wavenumbers to Herz:
408 freq[i] = 100. * w * SPEED_OF_LIGHT;
409
410 // Convert binary absorption cross-sections from
411 // cm^5 molec^(-2) to m^5 molec^(-2):
412 cia[nset][i] = c / 1e10;
413 }
414
415 nset++;
416 }
417
418 if (is.bad()) {
419 ostringstream os;
420 os << "Error in line " << nline << " reading CIA catalog file " << filename
421 << endl;
422
423 throw runtime_error(os.str());
424 }
425
426 AppendDataset(freq, Vector{temp}, cia);
427
428 // // For debugging
429 // for(Index i = 0; i < mdata.nelem(); i++)
430 // {
431 // cout << i << " ";
432 // cout << mdata[i].get_numeric_grid(0).nelem() << " ";
433 // cout << mdata[i].get_numeric_grid(1).nelem() << endl;
434 // }
435}
436
438void CIARecord::AppendDataset(const Vector& freq,
439 const Vector& temp,
440 const ArrayOfVector& cia) {
441 GriddedField2 dataset;
442 dataset.resize(freq.nelem(), temp.nelem());
443 dataset.set_grid(0, freq);
444 dataset.set_grid_name(0, "Frequency");
445
446 Vector temp_t;
447 temp_t = temp;
448 dataset.set_grid(1, temp_t);
449 dataset.set_grid_name(1, "Temperature");
450
451 for (Index t = 0; t < temp.nelem(); t++) dataset.data(joker, t) = cia[t];
452 mdata.push_back(dataset);
453}
454
457 for (Index ii = 0; ii < c2.DatasetCount(); ii++) {
458 mdata.push_back(c2.Dataset(ii));
459 }
460}
461
463
470ostream& operator<<(ostream& os, const CIARecord& /* cr */) {
471 os << "CIARecord output operator not yet implemented." << endl;
472 return os;
473}
Declarations required for the calculation of absorption coefficients.
Constants of physical expressions as constexpr.
void chk_interpolation_grids(const String &which_interpolation, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac, const bool islog)
Check interpolation grids.
Definition: check_input.cc:887
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.
Definition: cia.cc:212
constexpr Numeric SPEED_OF_LIGHT
Definition: cia.cc:25
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.
Definition: cia.cc:40
ostream & operator<<(ostream &os, const CIARecord &)
Output operator for CIARecord.
Definition: cia.cc:470
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.
Definition: cia.cc:40
This can be used to make arrays out of anything.
Definition: array.h:31
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:75
CIA data for a single pair of molecules.
Definition: cia.h:51
std::array< Species::Species, 2 > mspecies
The pair of molecules associated with these CIA data.
Definition: cia.h:237
ArrayOfGriddedField2 mdata
The data itself, directly from the HITRAN file.
Definition: cia.h:227
void AppendDataset(const CIARecord &c2)
Append other CIARecord to this.
Definition: cia.cc:456
void ReadFromCIA(const String &filename, const Verbosity &verbosity)
Read CIA catalog file.
Definition: cia.cc:275
void SetMoleculeName(const Index i, const String &name)
Set each molecule name (from a string) that is associated with this CIARecord.
Definition: cia.cc:250
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.
Definition: cia.cc:224
Index DatasetCount() const
Return number of datasets in this record.
Definition: cia.h:90
String MoleculeName(const Index i) const
Return each molecule name (as a string) that is associated with this CIARecord.
Definition: cia.cc:239
const GriddedField2 & Dataset(Index dataset) const
Return CIA dataset.
Definition: cia.h:112
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.
Definition: double_imanip.h:25
Index nelem() const
Definition: mystring.h:172
#define ARTS_ASSERT(condition,...)
Definition: debug.h:84
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:135
constexpr bool good_enum(EnumType x) noexcept
Checks if the enum number is good.
Definition: enums.h:21
void open_input_file(ifstream &file, const std::string_view name)
Open a file for reading.
Definition: file.cc:111
This file contains basic functions to handle ASCII files.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
#define CREATE_OUT2
Definition: messages.h:188
#define CREATE_OUTS
Definition: messages.h:191
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...
#define w
#define c