ARTS  2.4.0(git:4fb77825)
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 "abs_species_tags.h"
34 #include "absorption.h"
35 #include "file.h"
36 #include "interpolation_poly.h"
37 
38 extern const Numeric SPEED_OF_LIGHT;
39 
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  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  const 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  ArrayOfGridPosPoly f_gp(f_grid_active.nelem()), T_gp(1);
185  gridpos_poly(f_gp, data_f_grid, f_grid_active, f_order);
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 
191  Matrix itw(f_gp.nelem(), f_order + 1);
192  interpweights(itw, f_gp);
193  interp(result_active, itw, cia_data.data(joker, 0), f_gp);
194  } else {
195  // Temperature and frequency interpolation.
196 
197  // Find temperature grid position:
198  gridpos_poly(T_gp, data_T_grid, temperature, T_order, T_extrapolfac);
199 
200  // Calculate combined interpolation weights:
201  Tensor3 itw(f_gp.nelem(), T_gp.nelem(), (f_order + 1) * (T_order + 1));
202  interpweights(itw, f_gp, T_gp);
203 
204  // Make a matrix view of the results vector:
205  MatrixView result_matrix(result_active);
206 
207  // cout << "result_matrix r/c: " << result_matrix.nrows() << " / "
208  // << result_matrix.ncols() << endl;
209  // cout << "result_matrix before: " << result_matrix << endl;
210 
211  // cout << "cia_data: " << cia_data.data << endl;
212 
213  // Actual interpolation:
214  interp(result_matrix, itw, cia_data.data, f_gp, T_gp);
215 
216  // cout << "result_matrix after: " << result_matrix << endl;
217  }
218 
219  // cout << "result_active before: " << result_active << endl;
220 
221  // Set negative values to zero. (These could happen due to overshooting
222  // of the higher order interpolation.)
223  for (Index i = 0; i < result_active.nelem(); ++i)
224  if (result_active[i] < 0) result_active[i] = 0;
225 
226  // cout << "result_active after: " << result_active << endl;
227 }
228 
238  const Index sp1,
239  const Index sp2) {
240  for (Index i = 0; i < cia_data.nelem(); i++)
241  if ((cia_data[i].Species(0) == sp1 && cia_data[i].Species(1) == sp2) ||
242  (cia_data[i].Species(0) == sp2 && cia_data[i].Species(1) == sp1))
243  return i;
244 
245  return -1;
246 }
247 
248 // Documentation in header file.
251  const Numeric& temperature,
252  const Index& dataset,
253  const Numeric& T_extrapolfac,
254  const Index& robust,
255  const Verbosity& verbosity) const {
256  // If there is more than one dataset available for this species pair,
257  // we have to decide on which one to use. The rest is done by helper function
258  // cia_interpolate.
259 
260  // Make sure dataset index exists
261  if (dataset >= mdata.nelem()) {
262  ostringstream os;
263  os << "There are only " << mdata.nelem() << " datasets in this CIA file.\n"
264  << "But you are trying to use dataset " << dataset
265  << ". (Zero-based indexing.)";
266  throw runtime_error(os.str());
267  }
268 
269  // Get a handle on this dataset:
270  const GriddedField2& this_cia = mdata[dataset];
271 
273  result, f_grid, temperature, this_cia, T_extrapolfac, robust, verbosity);
274 }
275 
276 // Documentation in header file.
278  // Assert that i is 0 or 1:
279  assert(i >= 0);
280  assert(i <= 1);
281 
282  // The function species_name_from_species_index internally does an assertion
283  // that the species with this index really exists.
285 }
286 
287 // Documentation in header file.
288 void CIARecord::SetMoleculeName(const Index i, const String& name) {
289  // Assert that i is 0 or 1:
290  assert(i >= 0);
291  assert(i <= 1);
292 
293  // Find out the species index for name:
294  Index spec_ind = species_index_from_species_name(name);
295 
296  // Function species_index_from_species_name returns -1 if the species does
297  // not exist. Check this:
298  if (spec_ind < 0) {
299  ostringstream os;
300  os << "Species does not exist in ARTS: " << name;
301  throw runtime_error(os.str());
302  }
303 
304  // Assign species:
305  mspecies[i] = spec_ind;
306 }
307 
309 
316 void CIARecord::ReadFromCIA(const String& filename,
317  const Verbosity& verbosity) {
318  CREATE_OUT2;
319 
320  ifstream is;
321 
322  out2 << " Reading file: " << filename << "\n";
323  open_input_file(is, filename);
324 
325  // Number of points for spectral range in current dataset
326  Index npoints = -1;
327 
328  // Min/max wave numbers for this dataset
329  Numeric wave_min = -1.;
330  Numeric wave_max = -1.;
331 
332  // Current dataset index
333  Index ndataset = -1;
334 
335  // Current set in dataset
336  Index nset = -1;
337 
338  // Frequency, temp and cia values for current dataset
339  Vector freq;
341  ArrayOfVector cia;
342 
343  // Keep track of current line in file
344  Index nline = 0;
345 
346  mdata.resize(0);
347  istringstream istr;
348 
349  while (is) {
350  String line;
351 
352  // Extract needed information from dataset header line
354  nline++;
355  getline(is, line);
356  if (is.eof()) continue;
357 
358  if (line.nelem() < 100) {
359  ostringstream os;
360  os << "Error in line " << nline << " reading CIA catalog file "
361  << filename << endl
362  << "Header line unexpectedly short: " << endl
363  << line;
364 
365  throw runtime_error(os.str());
366  }
367 
368  if (is.bad()) {
369  ostringstream os;
370  os << "Error in line " << nline << " reading CIA catalog file "
371  << filename << endl;
372 
373  throw runtime_error(os.str());
374  }
375 
376  line.erase(0, 20);
377 
378  // Data for current set
379  Index set_npoints;
380  Numeric set_temp = NAN;
381  Numeric set_wave_min = NAN;
382  Numeric set_wave_max = NAN;
383 
384  istr.str(line);
385  istr.clear();
386  istr >> set_wave_min >> set_wave_max >> set_npoints >> set_temp;
387 
388  if (!istr || std::isnan(set_temp) || std::isnan(set_wave_min) ||
389  std::isnan(set_wave_max)) {
390  ostringstream os;
391  os << "Error in line " << nline << " reading CIA catalog file "
392  << filename << endl;
393 
394  throw runtime_error(os.str());
395  }
396 
397  // If the min/max wave numbers of this set are different from the
398  // previous one, a new dataset starts
399  if (npoints == -1 || wave_min != set_wave_min || wave_max != set_wave_max) {
400  if (ndataset != -1) AppendDataset(freq, temp, cia);
401 
402  npoints = set_npoints;
403  ndataset++;
404  wave_min = set_wave_min;
405  wave_max = set_wave_max;
406  nset = 0;
407  freq.resize(set_npoints);
408  temp.resize(0);
409  cia.resize(0);
410  }
411  if (npoints != set_npoints) {
412  ostringstream os;
413  os << "Error in line " << nline << " reading CIA catalog file "
414  << filename << endl
415  << "Inconsistent number of data points. Expected " << npoints
416  << ", got " << set_npoints;
417 
418  throw runtime_error(os.str());
419  }
420 
421  temp.push_back(set_temp);
422  cia.push_back(Vector(npoints));
423 
424  // Read dataset
426  for (Index i = 0; i < npoints; i++) {
427  Numeric w, c;
428 
429  nline++;
430  getline(is, line);
431 
432  w = c = NAN;
433  istr.str(line);
434  istr.clear();
435  istr >> w >> c;
436 
437  if (std::isnan(w) || std::isnan(c) || is.bad() || istr.bad()) {
438  ostringstream os;
439  os << "Error in line " << nline << " reading CIA catalog file "
440  << filename << ":" << endl
441  << line;
442 
443  throw runtime_error(os.str());
444  }
445 
446  // Convert wavenumbers to Herz:
447  freq[i] = 100. * w * SPEED_OF_LIGHT;
448 
449  // Convert binary absorption cross-sections from
450  // cm^5 molec^(-2) to m^5 molec^(-2):
451  cia[nset][i] = c / 1e10;
452  }
453 
454  nset++;
455  }
456 
457  if (is.bad()) {
458  ostringstream os;
459  os << "Error in line " << nline << " reading CIA catalog file " << filename
460  << endl;
461 
462  throw runtime_error(os.str());
463  }
464 
465  AppendDataset(freq, temp, cia);
466 
467  // // For debugging
468  // for(Index i = 0; i < mdata.nelem(); i++)
469  // {
470  // cout << i << " ";
471  // cout << mdata[i].get_numeric_grid(0).nelem() << " ";
472  // cout << mdata[i].get_numeric_grid(1).nelem() << endl;
473  // }
474 }
475 
478  const ArrayOfNumeric& temp,
479  const ArrayOfVector& cia) {
480  GriddedField2 dataset;
481  dataset.resize(freq.nelem(), temp.nelem());
482  dataset.set_grid(0, freq);
483  dataset.set_grid_name(0, "Frequency");
484 
485  Vector temp_t;
486  temp_t = temp;
487  dataset.set_grid(1, temp_t);
488  dataset.set_grid_name(1, "Temperature");
489 
490  for (Index t = 0; t < temp.nelem(); t++) dataset.data(joker, t) = cia[t];
491  mdata.push_back(dataset);
492 }
493 
496  for (Index ii = 0; ii < c2.DatasetCount(); ii++) {
497  mdata.push_back(c2.Dataset(ii));
498  }
499 }
500 
502 
509 ostream& operator<<(ostream& os, const CIARecord& /* cr */) {
510  os << "CIARecord output operator not yet implemented." << endl;
511  return os;
512 }
GriddedField2
Definition: gridded_fields.h:237
Matrix
The Matrix class.
Definition: matpackI.h:1193
GriddedField::set_grid_name
void set_grid_name(Index i, const String &s)
Set grid name.
Definition: gridded_fields.h:165
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
absorption.h
Declarations required for the calculation of absorption coefficients.
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
joker
const Joker joker
temp
#define temp
Definition: legacy_continua.cc:20951
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
cia_interpolation
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
chk_interpolation_grids
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:989
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
operator<<
ostream & operator<<(ostream &os, const CIARecord &)
Output operator for CIARecord.
Definition: cia.cc:509
open_input_file
void open_input_file(ifstream &file, const String &name)
Open a file for reading.
Definition: file.cc:147
Species
QuantumIdentifier::QType Index LowerQuantumNumbers Species
Definition: arts_api_classes.cc:255
CREATE_OUT2
#define CREATE_OUT2
Definition: messages.h:206
CIARecord::ReadFromCIA
void ReadFromCIA(const String &filename, const Verbosity &verbosity)
Read CIA catalog file.
Definition: cia.cc:316
GriddedField2::data
Matrix data
Definition: gridded_fields.h:281
Array< GridPosPoly >
species_index_from_species_name
Index species_index_from_species_name(String name)
Return species index for given species name.
Definition: absorption.cc:531
CIARecord::SetMoleculeName
void SetMoleculeName(const Index i, const String &name)
Set each molecule name (from a string) that is associated with this CIARecord.
Definition: cia.cc:288
GriddedField::get_numeric_grid
const Vector & get_numeric_grid(Index i) const
Get a numeric grid.
Definition: gridded_fields.cc:87
species_name_from_species_index
String species_name_from_species_index(const Index spec_ind)
Return species name for given species index.
Definition: absorption.cc:569
my_basic_string< char >
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
CIARecord::DatasetCount
Index DatasetCount() const
Return number of datasets in this record.
Definition: cia.h:106
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:49
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
abs_species_tags.h
Header file for stuff related to absorption species tags.
CIARecord::mdata
ArrayOfGriddedField2 mdata
The data itself, directly from the HITRAN file.
Definition: cia.h:232
my_basic_string::nelem
Index nelem() const
Number of elements.
Definition: mystring.h:246
interpolation_poly.h
Header file for interpolation_poly.cc.
CIARecord::Dataset
const GriddedField2 & Dataset(Index dataset) const
Return CIA dataset.
Definition: cia.h:128
oem::Vector
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
Range
The range class.
Definition: matpackI.h:160
GriddedField::set_grid
void set_grid(Index i, const Vector &g)
Set a numeric grid.
Definition: gridded_fields.cc:201
CIARecord
CIA data for a single pair of molecules.
Definition: cia.h:67
SPEED_OF_LIGHT
const Numeric SPEED_OF_LIGHT
file.h
This file contains basic functions to handle ASCII files.
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
gridpos_poly
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
The maximum difference from 1 that we allow for a sum check.
Definition: interpolation_poly.cc:120
CIARecord::MoleculeName
String MoleculeName(const Index i) const
Return each molecule name (as a string) that is associated with this CIARecord.
Definition: cia.cc:277
Vector
The Vector class.
Definition: matpackI.h:860
CIARecord::Extract
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:249
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
cia_get_index
Index cia_get_index(const ArrayOfCIARecord &cia_data, const Index sp1, const Index sp2)
Get the index in cia_data for the two given species.
Definition: cia.cc:237
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
cia.h
Header file for work with HITRAN collision induced absorption (CIA).
GriddedField2::resize
void resize(const GriddedField2 &gf)
Make this GriddedField2 the same size as the given one.
Definition: gridded_fields.h:271
CREATE_OUTS
#define CREATE_OUTS
Definition: messages.h:209
CIARecord::mspecies
Index mspecies[2]
The pair of molecules associated with these CIA data.
Definition: cia.h:242
CIARecord::AppendDataset
void AppendDataset(const CIARecord &c2)
Append other CIARecord to this.
Definition: cia.cc:495