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