ARTS 2.5.4 (git: 31ce4f0e)
cia.cc
Go to the documentation of this file.
1/* Copyright (C) 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
31#include "arts_constants.h"
32#include "check_input.h"
33#include "cia.h"
34#include <cmath>
35#include "species_tags.h"
36#include "absorption.h"
37#include "file.h"
39
41
56 ConstVectorView f_grid,
57 const Numeric& temperature,
58 const GriddedField2& cia_data,
59 const Numeric& T_extrapolfac,
60 const Index& robust,
61 const Verbosity& verbosity) {
63
64 const Index nf = f_grid.nelem();
65
66 // Assert that result vector has right size:
67 ARTS_ASSERT(result.nelem() == nf);
68
69 // Get data grids:
70 ConstVectorView data_f_grid = cia_data.get_numeric_grid(0);
71 ConstVectorView data_T_grid = cia_data.get_numeric_grid(1);
72
73 if (out3.sufficient_priority()) {
74 // Some detailed information to the most verbose output stream:
75 ostringstream os;
76 os << " f_grid: " << f_grid[0] << " - " << f_grid[nf - 1] << " Hz\n"
77 << " data_f_grid: " << data_f_grid[0] << " - "
78 << data_f_grid[data_f_grid.nelem() - 1] << " Hz\n"
79 << " temperature: " << temperature << " K\n"
80 << " data_T_grid: " << data_T_grid[0] << " - "
81 << data_T_grid[data_T_grid.nelem() - 1] << " K\n";
82 out3 << os.str();
83 }
84
85 // Initialize result to zero (important for those frequencies outside the data grid).
86 result = 0;
87
88 // We want to return result zero for all f_grid points that are outside the
89 // data_f_grid, because some CIA datasets are defined only where the absorption
90 // is not zero. So, we have to find out which part of f_grid is inside
91 // data_f_grid.
92 Index i_fstart, i_fstop;
93
94 for (i_fstart = 0; i_fstart < nf; ++i_fstart)
95 if (f_grid[i_fstart] >= data_f_grid[0]) break;
96
97 // Return directly if all frequencies are below data_f_grid:
98 if (i_fstart == nf) return;
99
100 for (i_fstop = nf - 1; i_fstop >= 0; --i_fstop)
101 if (f_grid[i_fstop] <= data_f_grid[data_f_grid.nelem() - 1]) break;
102
103 // Return directly if all frequencies are above data_f_grid:
104 if (i_fstop == -1) return;
105
106 // Extent for active frequency vector:
107 const Index f_extent = i_fstop - i_fstart + 1;
108
109 if (out3.sufficient_priority()) {
110 ostringstream os;
111 os << " " << f_extent << " frequency extraction points starting at "
112 << "frequency index " << i_fstart << ".\n";
113 out3 << os.str();
114 }
115
116 // If f_extent is less than one, then the entire data_f_grid is between two
117 // grid points of f_grid. (So that we do not have any f_grid points inside
118 // data_f_grid.) Return also in this case.
119 if (f_extent < 1) return;
120
121 // This is the part of f_grid for which we have to do the interpolation.
122 ConstVectorView f_grid_active = f_grid[Range(i_fstart, f_extent)];
123
124 // We have to create a matching view on the result vector:
125 VectorView result_active = result[Range(i_fstart, f_extent)];
126
127 // Decide on interpolation orders:
128 constexpr Index f_order = 3;
129
130 // The frequency grid has to have enough points for this interpolation
131 // order, otherwise throw a runtime error.
132 if (data_f_grid.nelem() < f_order + 1) {
133 ostringstream os;
134 os << "Not enough frequency grid points in CIA data.\n"
135 << "You have only " << data_f_grid.nelem() << " grid points.\n"
136 << "But need at least " << f_order + 1 << ".";
137 throw runtime_error(os.str());
138 }
139
140 // For T we have to be adaptive, since sometimes there is only one T in
141 // the data
142 Index T_order;
143 switch (data_T_grid.nelem()) {
144 case 1:
145 T_order = 0;
146 break;
147 case 2:
148 T_order = 1;
149 break;
150 case 3:
151 T_order = 2;
152 break;
153 default:
154 T_order = 3;
155 break;
156 }
157
158 // Check if frequency is inside the range covered by the data:
159 chk_interpolation_grids("Frequency interpolation for CIA continuum",
160 data_f_grid,
161 f_grid_active,
162 f_order);
163
164 // Check if temperature is inside the range covered by the data:
165 if (T_order > 0) {
166 try {
167 chk_interpolation_grids("Temperature interpolation for CIA continuum",
168 data_T_grid,
169 temperature,
170 T_order,
171 T_extrapolfac);
172 } catch (const std::runtime_error& e) {
173 // cout << "Gotcha!\n";
174 if (robust) {
175 // Just return NANs, but continue.
176 result_active = NAN;
177 return;
178 } else {
179 // Re-throw the exception.
180 throw runtime_error(e.what());
181 }
182 }
183 }
184
185 // Find frequency grid positions:
186 const auto f_lag = Interpolation::FixedLagrangeVector<f_order>(f_grid_active, data_f_grid);
187
188 // Do the rest of the interpolation.
189 if (T_order == 0) {
190 // No temperature interpolation in this case, just a frequency interpolation.
191 result_active = reinterp(cia_data.data(joker, 0), interpweights(f_lag), f_lag);
192 } else {
193 // Temperature and frequency interpolation.
194 const auto Tnew = std::array<double, 1>{temperature};
195 if (T_order == 1) {
196 const auto T_lag = Interpolation::FixedLagrangeVector<1>(Tnew, data_T_grid);
197 result_active = reinterp(cia_data.data, interpweights(f_lag, T_lag), f_lag, T_lag).reduce_rank<0>();
198 } else if (T_order == 2) {
199 const auto T_lag = Interpolation::FixedLagrangeVector<2>(Tnew, data_T_grid);
200 result_active = reinterp(cia_data.data, interpweights(f_lag, T_lag), f_lag, T_lag).reduce_rank<0>();
201 } else if (T_order == 3) {
202 const auto T_lag = Interpolation::FixedLagrangeVector<3>(Tnew, data_T_grid);
203 result_active = reinterp(cia_data.data, interpweights(f_lag, T_lag), f_lag, T_lag).reduce_rank<0>();
204 } else {
205 throw std::runtime_error("Cannot have this T_order, you must update the code...");
206 }
207 }
208
209 // cout << "result_active before: " << result_active << endl;
210
211 // Set negative values to zero. (These could happen due to overshooting
212 // of the higher order interpolation.)
213 for (Index i = 0; i < result_active.nelem(); ++i)
214 if (result_active[i] < 0) result_active[i] = 0;
215
216 // cout << "result_active after: " << result_active << endl;
217}
218
228 const Species::Species sp1,
229 const Species::Species sp2) {
230 for (Index i = 0; i < cia_data.nelem(); i++)
231 if ((cia_data[i].Species(0) == sp1 && cia_data[i].Species(1) == sp2) ||
232 (cia_data[i].Species(0) == sp2 && cia_data[i].Species(1) == sp1))
233 return i;
234
235 return -1;
236}
237
238// Documentation in header file.
240 ConstVectorView f_grid,
241 const Numeric& temperature,
242 const Index& dataset,
243 const Numeric& T_extrapolfac,
244 const Index& robust,
245 const Verbosity& verbosity) const {
246 // If there is more than one dataset available for this species pair,
247 // we have to decide on which one to use. The rest is done by helper function
248 // cia_interpolate.
249
250 // Make sure dataset index exists
251 if (dataset >= mdata.nelem()) {
252 ostringstream os;
253 os << "There are only " << mdata.nelem() << " datasets in this CIA file.\n"
254 << "But you are trying to use dataset " << dataset
255 << ". (Zero-based indexing.)";
256 throw runtime_error(os.str());
257 }
258
259 // Get a handle on this dataset:
260 const GriddedField2& this_cia = mdata[dataset];
261
263 result, f_grid, temperature, this_cia, T_extrapolfac, robust, verbosity);
264}
265
266// Documentation in header file.
268 // Assert that i is 0 or 1:
269 ARTS_ASSERT(i >= 0);
270 ARTS_ASSERT(i <= 1);
271
272 // The function species_name_from_species_index internally does an assertion
273 // that the species with this index really exists.
274 return Species::toShortName(mspecies[i]);
275}
276
277// Documentation in header file.
278void CIARecord::SetMoleculeName(const Index i, const String& name) {
279 // Assert that i is 0 or 1:
280 ARTS_ASSERT(i >= 0);
281 ARTS_ASSERT(i <= 1);
282
283 // Find out the species index for name:
284 Species::Species spec_ind = Species::fromShortName(name);
285
286 // Function species_index_from_species_name returns -1 if the species does
287 // not exist. Check this:
288 ARTS_USER_ERROR_IF(not good_enum(spec_ind),
289 "Species does not exist in ARTS: ", name)
290
291 // Assign species:
292 mspecies[i] = spec_ind;
293}
294
296
303void CIARecord::ReadFromCIA(const String& filename,
304 const Verbosity& verbosity) {
306
307 ifstream is;
308
309 out2 << " Reading file: " << filename << "\n";
310 open_input_file(is, filename);
311
312 // Number of points for spectral range in current dataset
313 Index npoints = -1;
314
315 // Min/max wave numbers for this dataset
316 Numeric wave_min = -1.;
317 Numeric wave_max = -1.;
318
319 // Current dataset index
320 Index ndataset = -1;
321
322 // Current set in dataset
323 Index nset = -1;
324
325 // Frequency, temp and cia values for current dataset
326 Vector freq;
328 ArrayOfVector cia;
329
330 // Keep track of current line in file
331 Index nline = 0;
332
333 mdata.resize(0);
334 istringstream istr;
335
336 while (is) {
337 String line;
338
339 // Extract needed information from dataset header line
341 nline++;
342 getline(is, line);
343 if (is.eof()) continue;
344
345 if (line.nelem() < 100) {
346 ostringstream os;
347 os << "Error in line " << nline << " reading CIA catalog file "
348 << filename << endl
349 << "Header line unexpectedly short: " << endl
350 << line;
351
352 throw runtime_error(os.str());
353 }
354
355 if (is.bad()) {
356 ostringstream os;
357 os << "Error in line " << nline << " reading CIA catalog file "
358 << filename << endl;
359
360 throw runtime_error(os.str());
361 }
362
363 line.erase(0, 20);
364
365 // Data for current set
366 Index set_npoints;
367 Numeric set_temp = NAN;
368 Numeric set_wave_min = NAN;
369 Numeric set_wave_max = NAN;
370
371 istr.str(line);
372 istr.clear();
373 istr >> double_imanip() >> set_wave_min >> set_wave_max;
374 istr >> set_npoints;
375 istr >> double_imanip() >> set_temp;
376
377 if (!istr || std::isnan(set_temp) || std::isnan(set_wave_min) ||
378 std::isnan(set_wave_max)) {
379 ostringstream os;
380 os << "Error in line " << nline << " reading CIA catalog file "
381 << filename << endl;
382
383 throw runtime_error(os.str());
384 }
385
386 // If the min/max wave numbers of this set are different from the
387 // previous one, a new dataset starts
388 if (npoints == -1 || wave_min != set_wave_min || wave_max != set_wave_max) {
389 if (ndataset != -1) AppendDataset(freq, temp, cia);
390
391 npoints = set_npoints;
392 ndataset++;
393 wave_min = set_wave_min;
394 wave_max = set_wave_max;
395 nset = 0;
396 freq.resize(set_npoints);
397 temp.resize(0);
398 cia.resize(0);
399 }
400 if (npoints != set_npoints) {
401 ostringstream os;
402 os << "Error in line " << nline << " reading CIA catalog file "
403 << filename << endl
404 << "Inconsistent number of data points. Expected " << npoints
405 << ", got " << set_npoints;
406
407 throw runtime_error(os.str());
408 }
409
410 temp.push_back(set_temp);
411 cia.push_back(Vector(npoints));
412
413 // Read dataset
415 for (Index i = 0; i < npoints; i++) {
416 Numeric w, c;
417
418 nline++;
419 getline(is, line);
420
421 w = c = NAN;
422 istr.str(line);
423 istr.clear();
424 istr >> w >> c;
425
426 if (std::isnan(w) || std::isnan(c) || is.bad() || istr.bad()) {
427 ostringstream os;
428 os << "Error in line " << nline << " reading CIA catalog file "
429 << filename << ":" << endl
430 << line;
431
432 throw runtime_error(os.str());
433 }
434
435 // Convert wavenumbers to Herz:
436 freq[i] = 100. * w * SPEED_OF_LIGHT;
437
438 // Convert binary absorption cross-sections from
439 // cm^5 molec^(-2) to m^5 molec^(-2):
440 cia[nset][i] = c / 1e10;
441 }
442
443 nset++;
444 }
445
446 if (is.bad()) {
447 ostringstream os;
448 os << "Error in line " << nline << " reading CIA catalog file " << filename
449 << endl;
450
451 throw runtime_error(os.str());
452 }
453
454 AppendDataset(freq, temp, cia);
455
456 // // For debugging
457 // for(Index i = 0; i < mdata.nelem(); i++)
458 // {
459 // cout << i << " ";
460 // cout << mdata[i].get_numeric_grid(0).nelem() << " ";
461 // cout << mdata[i].get_numeric_grid(1).nelem() << endl;
462 // }
463}
464
467 const ArrayOfNumeric& temp,
468 const ArrayOfVector& cia) {
469 GriddedField2 dataset;
470 dataset.resize(freq.nelem(), temp.nelem());
471 dataset.set_grid(0, freq);
472 dataset.set_grid_name(0, "Frequency");
473
474 Vector temp_t;
475 temp_t = temp;
476 dataset.set_grid(1, temp_t);
477 dataset.set_grid_name(1, "Temperature");
478
479 for (Index t = 0; t < temp.nelem(); t++) dataset.data(joker, t) = cia[t];
480 mdata.push_back(dataset);
481}
482
485 for (Index ii = 0; ii < c2.DatasetCount(); ii++) {
486 mdata.push_back(c2.Dataset(ii));
487 }
488}
489
491
498ostream& operator<<(ostream& os, const CIARecord& /* cr */) {
499 os << "CIARecord output operator not yet implemented." << endl;
500 return os;
501}
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:906
void cia_interpolation(VectorView result, 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:55
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:227
constexpr Numeric SPEED_OF_LIGHT
Definition: cia.cc:40
ostream & operator<<(ostream &os, const CIARecord &)
Output operator for CIARecord.
Definition: cia.cc:498
Header file for work with HITRAN collision induced absorption (CIA).
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
CIA data for a single pair of molecules.
Definition: cia.h:68
std::array< Species::Species, 2 > mspecies
The pair of molecules associated with these CIA data.
Definition: cia.h:257
ArrayOfGriddedField2 mdata
The data itself, directly from the HITRAN file.
Definition: cia.h:247
void AppendDataset(const CIARecord &c2)
Append other CIARecord to this.
Definition: cia.cc:484
void ReadFromCIA(const String &filename, const Verbosity &verbosity)
Read CIA catalog file.
Definition: cia.cc:303
void SetMoleculeName(const Index i, const String &name)
Set each molecule name (from a string) that is associated with this CIARecord.
Definition: cia.cc:278
Index DatasetCount() const
Return number of datasets in this record.
Definition: cia.h:107
void Extract(VectorView result, ConstVectorView f_grid, const Numeric &temperature, const Index &dataset, const Numeric &T_extrapolfac, const Index &robust, const Verbosity &verbosity) const
Vector version of extract.
Definition: cia.cc:239
String MoleculeName(const Index i) const
Return each molecule name (as a string) that is associated with this CIARecord.
Definition: cia.cc:267
const GriddedField2 & Dataset(Index dataset) const
Return CIA dataset.
Definition: cia.h:129
A constant view of a Vector.
Definition: matpackI.h:512
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
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.
The range class.
Definition: matpackI.h:159
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
Input manipulator class for doubles to enable nan and inf parsing.
Definition: double_imanip.h:42
Index nelem() const
Definition: mystring.h:189
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
constexpr bool good_enum(EnumType x) noexcept
Checks if the enum number is good.
Definition: enums.h:22
void open_input_file(ifstream &file, const std::string_view name)
Open a file for reading.
Definition: file.cc:128
This file contains basic functions to handle ASCII files.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
#define temp
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
#define CREATE_OUT2
Definition: messages.h:205
#define CREATE_OUTS
Definition: messages.h:208
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...
constexpr Numeric e
Elementary charge convenience name [C].
void reinterp(VectorView out, const ConstVectorView &iy, const Grid< Vector, 1 > &iw, const Array< Lagrange > &dim0)
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:53
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
#define w
#define c