ARTS  2.2.66
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 <cmath>
32 #include "cia.h"
33 #include "interpolation_poly.h"
34 #include "absorption.h"
35 #include "abs_species_tags.h"
36 #include "file.h"
37 
38 extern 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)
60 {
62 
63  const Index nf = f_grid.nelem();
64 
65  // Assert that result vector has right size:
66  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  {
74  // Some detailed information to the most verbose output stream:
75  ostringstream os;
76  os << " f_grid: " << f_grid[0] << " - "
77  << f_grid[nf-1] << " Hz\n"
78  << " data_f_grid: " << data_f_grid[0] << " - "
79  << data_f_grid[data_f_grid.nelem()-1] << " Hz\n"
80  << " temperature: " << temperature << " K\n"
81  << " data_T_grid: " << data_T_grid[0] << " - "
82  << data_T_grid[data_T_grid.nelem()-1] << " K\n";
83  out3 << os.str();
84  }
85 
86  // Initialize result to zero (important for those frequencies outside the data grid).
87  result = 0;
88 
89  // We want to return result zero for all f_grid points that are outside the
90  // data_f_grid, because some CIA datasets are defined only where the absorption
91  // is not zero. So, we have to find out which part of f_grid is inside
92  // data_f_grid.
93  Index i_fstart, i_fstop;
94 
95  for (i_fstart=0; i_fstart<nf; ++i_fstart)
96  if (f_grid[i_fstart] >= data_f_grid[0]) break;
97 
98  // Return directly if all frequencies are below data_f_grid:
99  if (i_fstart==nf) return;
100 
101  for (i_fstop=nf-1; i_fstop>=0; --i_fstop)
102  if (f_grid[i_fstop] <= data_f_grid[data_f_grid.nelem()-1]) break;
103 
104  // Return directly if all frequencies are above data_f_grid:
105  if (i_fstop==-1) return;
106 
107  // Extent for active frequency vector:
108  const Index f_extent = i_fstop-i_fstart+1;
109 
110  if (out3.sufficient_priority())
111  {
112  ostringstream os;
113  os << " " << f_extent << " frequency extraction points starting at "
114  << "frequency index " << i_fstart << ".\n";
115  out3 << os.str();
116  }
117 
118  // If f_extent is less than one, then the entire data_f_grid is between two
119  // grid points of f_grid. (So that we do not have any f_grid points inside
120  // data_f_grid.) Return also in this case.
121  if (f_extent < 1) return;
122 
123 
124  // This is the part of f_grid for which we have to do the interpolation.
125  ConstVectorView f_grid_active = f_grid[Range(i_fstart, f_extent)];
126 
127  // We have to create a matching view on the result vector:
128  VectorView result_active = result[Range(i_fstart, f_extent)];
129 
130 
131  // Decide on interpolation orders:
132  const Index f_order = 3;
133 
134  // The frequency grid has to have enough points for this interpolation
135  // order, otherwise throw a runtime error.
136  if ( data_f_grid.nelem() < f_order+1 )
137  {
138  ostringstream os;
139  os << "Not enough frequency grid points in CIA data.\n"
140  << "You have only " << data_f_grid.nelem() << " grid points.\n"
141  << "But need at least " << f_order+1 << ".";
142  throw runtime_error(os.str());
143  }
144 
145 
146  // For T we have to be adaptive, since sometimes there is only one T in
147  // the data
148  Index T_order;
149  switch (data_T_grid.nelem()) {
150  case 1:
151  T_order = 0;
152  break;
153  case 2:
154  T_order = 1;
155  break;
156  case 3:
157  T_order = 2;
158  break;
159  default:
160  T_order = 3;
161  break;
162  }
163 
164  // Check if frequency is inside the range covered by the data:
165  chk_interpolation_grids("Frequency interpolation for CIA continuum",
166  data_f_grid,
167  f_grid_active,
168  f_order);
169 
170 
171  // Check if temperature is inside the range covered by the data:
172  if (T_order > 0) {
173  try {
174  chk_interpolation_grids("Temperature interpolation for CIA continuum",
175  data_T_grid,
176  temperature,
177  T_order,
178  T_extrapolfac);
179  } catch (runtime_error e) {
180  // cout << "Gotcha!\n";
181  if (robust) {
182  // Just return NANs, but continue.
183  result_active = NAN;
184  return;
185  } else {
186  // Re-throw the exception.
187  throw runtime_error(e.what());
188  }
189  }
190  }
191 
192  // Find frequency grid positions:
193  ArrayOfGridPosPoly f_gp(f_grid_active.nelem()), T_gp(1);
194  gridpos_poly(f_gp, data_f_grid, f_grid_active, f_order);
195 
196  // Do the rest of the interpolation.
197  if (T_order == 0)
198  {
199  // No temperature interpolation in this case, just a frequency interpolation.
200 
201  Matrix itw(f_gp.nelem(),f_order+1);
202  interpweights(itw,f_gp);
203  interp(result_active, itw, cia_data.data(joker,0), f_gp);
204  }
205  else
206  {
207  // Temperature and frequency interpolation.
208 
209  // Find temperature grid position:
210  gridpos_poly(T_gp, data_T_grid, temperature, T_order, T_extrapolfac);
211 
212  // Calculate combined interpolation weights:
213  Tensor3 itw(f_gp.nelem(),T_gp.nelem(),(f_order+1)*(T_order+1));
214  interpweights(itw, f_gp, T_gp);
215 
216  // Make a matrix view of the results vector:
217  MatrixView result_matrix(result_active);
218 
219 // cout << "result_matrix r/c: " << result_matrix.nrows() << " / "
220 // << result_matrix.ncols() << endl;
221 // cout << "result_matrix before: " << result_matrix << endl;
222 
223  // cout << "cia_data: " << cia_data.data << endl;
224 
225  // Actual interpolation:
226  interp(result_matrix, itw, cia_data.data, f_gp, T_gp);
227 
228 // cout << "result_matrix after: " << result_matrix << endl;
229  }
230 
231 // cout << "result_active before: " << result_active << endl;
232 
233  // Set negative values to zero. (These could happen due to overshooting
234  // of the higher order interpolation.)
235  for (Index i=0; i<result_active.nelem(); ++i)
236  if (result_active[i]<0)
237  result_active[i] = 0;
238 
239 // cout << "result_active after: " << result_active << endl;
240 
241 }
242 
243 
253  const Index sp1, const Index sp2)
254 {
255  for (Index i = 0; i < cia_data.nelem(); i++)
256  if ((cia_data[i].Species(0) == sp1 && cia_data[i].Species(1) == sp2)
257  || (cia_data[i].Species(0) == sp2 && cia_data[i].Species(1) == sp1))
258  return i;
259 
260  return -1;
261 }
262 
263 
264 // Documentation in header file.
266  ConstVectorView f_grid,
267  const Numeric& temperature,
268  const Index& dataset,
269  const Numeric& T_extrapolfac,
270  const Index& robust,
271  const Verbosity& verbosity) const
272 {
273  // If there is more than one dataset available for this species pair,
274  // we have to decide on which one to use. The rest is done by helper function
275  // cia_interpolate.
276 
277  // Make sure dataset index exists
278  if ( dataset >= mdata.nelem() )
279  {
280  ostringstream os;
281  os << "There are only " << mdata.nelem() << " datasets in this CIA file.\n"
282  << "But you are trying to use dataset " << dataset << ". (Zero-based indexing.)";
283  throw runtime_error(os.str());
284  }
285 
286 
287  // Get a handle on this dataset:
288  const GriddedField2& this_cia = mdata[dataset];
289 
290  cia_interpolation(result,
291  f_grid,
292  temperature,
293  this_cia,
294  T_extrapolfac,
295  robust,
296  verbosity);
297 }
298 
299 
300 // Documentation in header file.
302 {
303  // Assert that i is 0 or 1:
304  assert(i>=0);
305  assert(i<=1);
306 
307  // The function species_name_from_species_index internally does an assertion
308  // that the species with this index really exists.
310 }
311 
312 // Documentation in header file.
314  const String& name)
315 {
316  // Assert that i is 0 or 1:
317  assert(i>=0);
318  assert(i<=1);
319 
320  // Find out the species index for name:
321  Index spec_ind = species_index_from_species_name(name);
322 
323  // Function species_index_from_species_name returns -1 if the species does
324  // not exist. Check this:
325  if ( spec_ind < 0 )
326  {
327  ostringstream os;
328  os << "Species does not exist in ARTS: " << name;
329  throw runtime_error(os.str());
330  }
331 
332  // Assign species:
333  mspecies[i] = spec_ind;
334 }
335 
336 
338 
345 void CIARecord::ReadFromCIA(const String& filename, const Verbosity& verbosity)
346 {
347  CREATE_OUT2;
348 
349  ifstream is;
350 
351  out2 << " Reading file: " << filename << "\n";
352  open_input_file(is, filename);
353 
354  // Number of points for spectral range in current dataset
355  Index npoints = -1;
356 
357  // Min/max wave numbers for this dataset
358  Numeric wave_min = -1.;
359  Numeric wave_max = -1.;
360 
361  // Current dataset index
362  Index ndataset = -1;
363 
364  // Current set in dataset
365  Index nset = -1;
366 
367  // Frequency, temp and cia values for current dataset
368  Vector freq;
370  ArrayOfVector cia;
371 
372  // Keep track of current line in file
373  Index nline = 0;
374 
375  mdata.resize(0);
376  istringstream istr;
377 
378  while (is)
379  {
380  String line;
381 
382  // Extract needed information from dataset header line
384  nline++;
385  getline(is,line);
386  if (is.eof()) continue;
387 
388  if (line.nelem() < 100)
389  {
390  ostringstream os;
391  os << "Error in line " << nline
392  << " reading CIA catalog file " << filename << endl
393  << "Header line unexpectedly short: " << endl << line;
394 
395  throw runtime_error(os.str());
396  }
397 
398  if (is.bad())
399  {
400  ostringstream os;
401  os << "Error in line " << nline
402  << " reading CIA catalog file " << filename << endl;
403 
404  throw runtime_error(os.str());
405  }
406 
407  line.erase(0, 20);
408 
409  // Data for current set
410  Index set_npoints;
411  Numeric set_temp = NAN;
412  Numeric set_wave_min = NAN;
413  Numeric set_wave_max = NAN;
414 
415  istr.str(line);
416  istr.clear();
417  istr >> set_wave_min >> set_wave_max >> set_npoints >> set_temp;
418 
419  if (!istr || isnan(set_temp) || isnan(set_wave_min) || isnan(set_wave_max))
420  {
421  ostringstream os;
422  os << "Error in line " << nline
423  << " reading CIA catalog file " << filename << endl;
424 
425  throw runtime_error(os.str());
426  }
427 
428  // If the min/max wave numbers of this set are different from the
429  // previous one, a new dataset starts
430  if (npoints == -1 || wave_min != set_wave_min || wave_max != set_wave_max)
431  {
432  if (ndataset != -1)
433  AppendDataset(freq, temp, cia);
434 
435  npoints = set_npoints;
436  ndataset++;
437  wave_min = set_wave_min;
438  wave_max = set_wave_max;
439  nset = 0;
440  freq.resize(set_npoints);
441  temp.resize(0);
442  cia.resize(0);
443  }
444  if (npoints != set_npoints)
445  {
446  ostringstream os;
447  os << "Error in line " << nline
448  << " reading CIA catalog file " << filename << endl
449  << "Inconsistent number of data points. Expected " << npoints
450  << ", got " << set_npoints;
451 
452  throw runtime_error(os.str());
453  }
454 
455  temp.push_back(set_temp);
456  cia.push_back(Vector(npoints));
457 
458  // Read dataset
460  for (Index i = 0; i < npoints; i++)
461  {
462  Numeric w, c;
463 
464  nline++;
465  getline(is,line);
466 
467  w = c = NAN;
468  istr.str(line);
469  istr.clear();
470  istr >> w >> c;
471 
472  if (isnan(w) || isnan(c) || is.bad() || istr.bad())
473  {
474  ostringstream os;
475  os << "Error in line " << nline
476  << " reading CIA catalog file " << filename << ":" << endl
477  << line;
478 
479  throw runtime_error(os.str());
480  }
481 
482  // Convert wavenumbers to Herz:
483  freq[i] = 100. * w * SPEED_OF_LIGHT;
484 
485  // Convert binary absorption cross-sections from
486  // cm^5 molec^(-2) to m^5 molec^(-2):
487  cia[nset][i] = c / 1e10;
488  }
489 
490  nset++;
491  }
492 
493  if (is.bad())
494  {
495  ostringstream os;
496  os << "Error in line " << nline
497  << " reading CIA catalog file " << filename << endl;
498 
499  throw runtime_error(os.str());
500  }
501 
502  AppendDataset(freq, temp, cia);
503 
504 // // For debugging
505 // for(Index i = 0; i < mdata.nelem(); i++)
506 // {
507 // cout << i << " ";
508 // cout << mdata[i].get_numeric_grid(0).nelem() << " ";
509 // cout << mdata[i].get_numeric_grid(1).nelem() << endl;
510 // }
511 }
512 
513 
516  const ArrayOfNumeric& temp,
517  const ArrayOfVector& cia)
518 {
519  GriddedField2 dataset;
520  dataset.resize(freq.nelem(), temp.nelem());
521  dataset.set_grid(0, freq);
522  dataset.set_grid_name(0, "Frequency");
523 
524  Vector temp_t;
525  temp_t = temp;
526  dataset.set_grid(1, temp_t);
527  dataset.set_grid_name(1, "Temperature");
528 
529  for (Index t = 0; t < temp.nelem(); t++)
530  dataset.data(joker, t) = cia[t];
531  mdata.push_back(dataset);
532 }
533 
534 
536 
543 ostream& operator<<(ostream& os, const CIARecord& /* cr */)
544 {
545  os << "CIARecord output operator not yet implemented." << endl;
546  return os;
547 }
548 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
Header file for work with HITRAN collision induced absorption (CIA).
The VectorView class.
Definition: matpackI.h:372
Index nelem() const
Number of elements.
Definition: array.h:176
ConstVectorView get_numeric_grid(Index i) const
Get a numeric grid.
The Vector class.
Definition: matpackI.h:556
The MatrixView class.
Definition: matpackI.h:679
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
void ReadFromCIA(const String &filename, const Verbosity &verbosity)
Read CIA catalog file.
Definition: cia.cc:345
ArrayOfGriddedField2 mdata
The data itself, directly from the HITRAN file.
Definition: cia.h:246
The range class.
Definition: matpackI.h:148
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:679
This file contains basic functions to handle ASCII files.
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
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.
String species_name_from_species_index(const Index spec_ind)
Return species name for given species index.
Definition: absorption.cc:1301
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:265
void AppendDataset(const Vector &freq, const ArrayOfNumeric &temp, const ArrayOfVector &cia)
Append dataset to mdata.
Definition: cia.cc:515
The Tensor3 class.
Definition: matpackIII.h:348
Index mspecies[2]
The pair of molecules associated with these CIA data.
Definition: cia.h:256
#define CREATE_OUTS
Definition: messages.h:216
void set_grid(Index i, const Vector &g)
Set a numeric grid.
void resize(const GriddedField2 &gf)
Make this GriddedField2 the same size as the given one.
const Joker joker
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
Index nelem() const
Number of elements.
Definition: mystring.h:278
Declarations required for the calculation of absorption coefficients.
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
Set up grid positions for higher order interpolation.
void SetMoleculeName(const Index i, const String &name)
Set each molecule name (from a string) that is associated with this CIARecord.
Definition: cia.cc:313
Index species_index_from_species_name(String name)
Return species index for given species name.
Definition: absorption.cc:1260
void open_input_file(ifstream &file, const String &name)
Open a file for reading.
Definition: file.cc:166
Header file for interpolation_poly.cc.
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
A constant view of a Vector.
Definition: matpackI.h:292
String MoleculeName(const Index i) const
Return each molecule name (as a string) that is associated with this CIARecord.
Definition: cia.cc:301
Header file for stuff related to absorption species tags.
#define temp
Definition: continua.cc:20773
void set_grid_name(Index i, const String &s)
Set grid name.
CIA data for a single pair of molecules.
Definition: cia.h:68
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:252
ostream & operator<<(ostream &os, const CIARecord &)
Output operator for CIARecord.
Definition: cia.cc:543
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
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
#define CREATE_OUT2
Definition: messages.h:213
const Numeric SPEED_OF_LIGHT