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