ARTS  2.3.1285(git:92a29ea9-dirty)
hitran_xsec.cc
Go to the documentation of this file.
1 /* Copyright (C) 2018 Oliver Lemke <oliver.lemke@uni-hamburg.de>
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 
26 #include "arts.h"
27 #include "interpolation_poly.h"
28 
29 #ifdef ENABLE_FFTW
30 
31 #include <complex.h>
32 #include <fftw3.h>
33 
34 #endif /* ENABLE_FFTW */
35 
36 #include "absorption.h"
37 #include "check_input.h"
38 #include "hitran_xsec.h"
39 
40 extern const Numeric PI;
41 
42 Numeric func_2straights(const Numeric x, const Vector& coeffs) {
43  assert(coeffs.nelem() == 3);
44  return (x <= coeffs[0]) ? coeffs[1] * x
45  : coeffs[2] * (x - coeffs[0]) + coeffs[1] * coeffs[0];
46 }
47 
48 Numeric lorentz_pdf(const Numeric x, const Numeric x0, const Numeric gamma) {
49  const Numeric xx0 = x - x0;
50  return gamma / PI / (xx0 * xx0 + gamma * gamma);
51 }
52 
54  // The function species_name_from_species_index internally does an assertion
55  // that the species with this index really exists.
57 }
58 
59 void convolve(Vector& result,
60  const ConstVectorView& xsec,
61  const ConstVectorView& lorentz) {
62  Index n_xsec = xsec.nelem();
63  Index n_lorentz = lorentz.nelem();
64  // assert(n_xsec == n_lorentz);
65  Vector temp(n_xsec + n_lorentz - 1);
66 
67  for (Index i = 0; i < n_xsec + n_lorentz - 1; ++i) {
68  Numeric sum = 0.0;
69  for (Index j = 0; j <= i; ++j) {
70  sum += ((j < n_xsec) && (i - j < n_lorentz)) ? xsec[j] * lorentz[i - j]
71  : 0.0;
72  }
73  temp[i] = sum;
74  }
75  result = temp[Range(n_lorentz / 2, n_xsec, 1)];
76 }
77 
78 #ifdef ENABLE_FFTW
79 
80 void fftconvolve(VectorView& result,
81  const Vector& xsec,
82  const Vector& lorentz) {
83  int n_p = (int)(xsec.nelem() + lorentz.nelem() - 1);
84  int n_p_2 = n_p / 2 + 1;
85 
86  double* xsec_in = fftw_alloc_real((size_t)n_p);
87  fftw_complex* xsec_out = fftw_alloc_complex((size_t)n_p_2);
88  memcpy(xsec_in, xsec.get_c_array(), sizeof(double) * xsec.nelem());
89  memset(&xsec_in[xsec.nelem()], 0, sizeof(double) * (n_p - xsec.nelem()));
90 
91  fftw_plan plan;
92 #pragma omp critical(fftw_call)
93  plan = fftw_plan_dft_r2c_1d(n_p, xsec_in, xsec_out, FFTW_ESTIMATE);
94 
95  fftw_execute(plan);
96 
97 #pragma omp critical(fftw_call)
98  fftw_destroy_plan(plan);
99 
100  fftw_free(xsec_in);
101 
102  double* lorentz_in = fftw_alloc_real((size_t)n_p);
103  fftw_complex* lorentz_out = fftw_alloc_complex((size_t)n_p_2);
104  memcpy(lorentz_in, lorentz.get_c_array(), sizeof(double) * lorentz.nelem());
105  memset(&lorentz_in[lorentz.nelem()],
106  0,
107  sizeof(double) * (n_p - lorentz.nelem()));
108 
109 #pragma omp critical(fftw_call)
110  plan = fftw_plan_dft_r2c_1d(n_p, lorentz_in, lorentz_out, FFTW_ESTIMATE);
111 
112  fftw_execute(plan);
113 
114 #pragma omp critical(fftw_call)
115  fftw_destroy_plan(plan);
116 
117  fftw_free(lorentz_in);
118 
119  fftw_complex* fft_in = fftw_alloc_complex((size_t)n_p_2);
120  double* fft_out = fftw_alloc_real((size_t)n_p);
121  memcpy(fft_in, xsec_out, sizeof(fftw_complex) * n_p_2);
122 
123  for (Index i = 0; i < n_p_2; i++) {
124  fft_in[i][0] =
125  xsec_out[i][0] * lorentz_out[i][0] - xsec_out[i][1] * lorentz_out[i][1];
126  fft_in[i][1] =
127  xsec_out[i][0] * lorentz_out[i][1] + xsec_out[i][1] * lorentz_out[i][0];
128  }
129 
130 #pragma omp critical(fftw_call)
131  plan = fftw_plan_dft_c2r_1d(n_p, fft_in, fft_out, FFTW_ESTIMATE);
132 
133  fftw_execute(plan);
134 
135 #pragma omp critical(fftw_call)
136  fftw_destroy_plan(plan);
137 
138  fftw_free(fft_in);
139 
140  for (Index i = 0; i < xsec.nelem(); i++) {
141  result[i] = fft_out[i + (int)lorentz.nelem() / 2] / n_p;
142  }
143 
144  fftw_free(xsec_out);
145  fftw_free(lorentz_out);
146  fftw_free(fft_out);
147 }
148 
149 #endif /* ENABLE_FFTW */
150 
152  ConstVectorView f_grid,
153  const Numeric& pressure,
154  const Numeric& temperature,
155  const Index& apply_tfit,
156  const Verbosity& verbosity) const {
157  CREATE_OUTS;
158 
159  const Index nf = f_grid.nelem();
160 
161  // Assert that result vector has right size:
162  assert(result.nelem() == nf);
163 
164  // Initialize result to zero (important for those frequencies outside the data grid).
165  result = 0.;
166 
167  const Index ndatasets = mxsecs.nelem();
168  for (Index this_dataset_i = 0; this_dataset_i < ndatasets; this_dataset_i++) {
169  const Vector& data_f_grid = mfgrids[this_dataset_i];
170  const Numeric fmin = data_f_grid[0];
171  const Numeric fmax = data_f_grid[mfgrids[this_dataset_i].nelem() - 1];
172  const Index data_nf = mfgrids[this_dataset_i].nelem();
173 
174  if (out3.sufficient_priority()) {
175  // Some detailed information to the most verbose output stream:
176  ostringstream os;
177  os << " f_grid: " << f_grid[0] << " - " << f_grid[nf - 1]
178  << " Hz\n"
179  << " data_f_grid: " << fmin << " - " << fmax << " Hz\n"
180  << " pressure: " << pressure << " K\n";
181  out3 << os.str();
182  }
183 
184  // We want to return result zero for all f_grid points that are outside the
185  // data_f_grid, because xsec datasets are defined only where the absorption
186  // was measured. So, we have to find out which part of f_grid is inside
187  // data_f_grid.
188  Index i_fstart, i_fstop;
189 
190  for (i_fstart = 0; i_fstart < nf; ++i_fstart)
191  if (f_grid[i_fstart] >= fmin) break;
192 
193  // Return directly if all frequencies are below data_f_grid:
194  if (i_fstart == nf) continue;
195 
196  for (i_fstop = nf - 1; i_fstop >= 0; --i_fstop)
197  if (f_grid[i_fstop] <= fmax) break;
198 
199  // Return directly if all frequencies are above data_f_grid:
200  if (i_fstop == -1) continue;
201 
202  // Extent for active frequency vector:
203  const Index f_extent = i_fstop - i_fstart + 1;
204 
205  if (out3.sufficient_priority()) {
206  ostringstream os;
207  os << " " << f_extent << " frequency extraction points starting at "
208  << "frequency index " << i_fstart << ".\n";
209  out3 << os.str();
210  }
211 
212  // If f_extent is less than one, then the entire data_f_grid is between two
213  // grid points of f_grid. (So that we do not have any f_grid points inside
214  // data_f_grid.) Return also in this case.
215  if (f_extent < 3) continue;
216 
217  // This is the part of f_grid for which we have to do the interpolation.
218  ConstVectorView f_grid_active = f_grid[Range(i_fstart, f_extent)];
219 
220  // We also need to determine the range in the xsec dataset that's inside
221  // f_grid. We can ignore the remaining data.
222  Index i_data_fstart, i_data_fstop;
223 
224  for (i_data_fstart = 0; i_data_fstart < data_nf; ++i_data_fstart)
225  if (data_f_grid[i_data_fstart] >= fmin) break;
226 
227  for (i_data_fstop = data_nf - 1; i_data_fstop >= 0; --i_data_fstop)
228  if (data_f_grid[i_data_fstop] <= fmax) break;
229 
230  // Extent for active data frequency vector:
231  const Index data_f_extent = i_data_fstop - i_data_fstart + 1;
232 
233  // This is the part of f_grid for which we have to do the interpolation.
234  ConstVectorView data_f_grid_active =
235  data_f_grid[Range(i_data_fstart, data_f_extent)];
236 
237  // This is the part of the xsec dataset for which we have to do the
238  // interpolation.
239  Range active_range(i_data_fstart, data_f_extent);
240  ConstVectorView xsec_active = mxsecs[this_dataset_i][active_range];
241 
242  Vector xsec_active_tfit;
243 
244  if (apply_tfit != 0 && mtslope[this_dataset_i].nelem() > 1) {
245  xsec_active_tfit = mtslope[this_dataset_i][active_range];
246  xsec_active_tfit *= temperature - mreftemperature[this_dataset_i];
247  xsec_active_tfit += mtintersect[this_dataset_i][active_range];
248  xsec_active_tfit /= 10000;
249  xsec_active_tfit += mxsecs[this_dataset_i][active_range];
250 
251  xsec_active = xsec_active_tfit;
252  }
253  // We have to create a matching view on the result vector:
254  VectorView result_active = result[Range(i_fstart, f_extent)];
255  Vector xsec_interp(f_extent);
256 
257  // Decide on interpolation orders:
258  const Index f_order = 3;
259 
260  // The frequency grid has to have enough points for this interpolation
261  // order, otherwise throw a runtime error.
262  if (data_f_grid.nelem() < f_order + 1) {
263  ostringstream os;
264  os << "Not enough frequency grid points in Hitran Xsec data.\n"
265  << "You have only " << data_f_grid.nelem() << " grid points.\n"
266  << "But need at least " << f_order + 1 << ".";
267  throw runtime_error(os.str());
268  }
269 
270  if (pressure > mrefpressure[this_dataset_i] &&
271  mrefpressure[this_dataset_i] > 0.) {
272  // Apply pressure dependent broadening and set negative values to zero.
273  // (These could happen due to overshooting of the higher order interpolation.)
274  const Numeric pdiff = pressure - mrefpressure[this_dataset_i];
275  const Numeric fwhm = func_2straights(pdiff, mcoeffs);
276  // std::cout << mcoeffs << " - ";
277  // std::cout << "pdiff: " << pdiff << " - fwhm: " << fwhm << " - fstep: "
278  // << f_grid[i_fstart] + f_grid[i_fstart + 1] << std::endl;
279 
280  Vector f_lorentz(data_f_extent);
281  Numeric lsum = 0.;
282  for (Index i = 0; i < data_f_extent; i++) {
283  f_lorentz[i] =
284  lorentz_pdf(data_f_grid[i_data_fstart + i],
285  data_f_grid[i_data_fstart + data_f_extent / 2 - 1],
286  fwhm / 2.);
287  lsum += f_lorentz[i];
288  }
289 
290  f_lorentz /= lsum;
291 
292  Vector data_result(xsec_active.nelem());
293 #ifdef ENABLE_FFTW
294  fftconvolve(
295  data_result,
296  xsec_active,
297  f_lorentz[Range(f_lorentz.nelem() / 4, f_lorentz.nelem() / 2, 1)]);
298 #else
299  convolve(
300  data_result,
301  xsec_active,
302  f_lorentz[Range(f_lorentz.nelem() / 4, f_lorentz.nelem() / 2, 1)]);
303 #endif /* ENABLE_FFTW */
304 
305  // TODO: Add to result_active here
306  // Check if frequency is inside the range covered by the data:
307  chk_interpolation_grids("Frequency interpolation for cross sections",
308  data_f_grid,
309  f_grid_active,
310  f_order);
311 
312  {
313  // Find frequency grid positions:
314  ArrayOfGridPosPoly f_gp(f_grid_active.nelem()), T_gp(1);
315  gridpos_poly(f_gp, data_f_grid_active, f_grid_active, f_order);
316 
317  Matrix itw(f_gp.nelem(), f_order + 1);
318  interpweights(itw, f_gp);
319  interp(xsec_interp, itw, data_result, f_gp);
320  }
321  } else {
322  // Find frequency grid positions:
323  ArrayOfGridPosPoly f_gp(f_grid_active.nelem()), T_gp(1);
324  gridpos_poly(f_gp, data_f_grid_active, f_grid_active, f_order);
325 
326  Matrix itw(f_gp.nelem(), f_order + 1);
327  interpweights(itw, f_gp);
328  interp(xsec_interp, itw, xsec_active, f_gp);
329  }
330 
331  result_active += xsec_interp;
332  }
333 }
334 
343  const Index species) {
344  for (Index i = 0; i < xsec_data.nelem(); i++)
345  if (xsec_data[i].Species() == species) return i;
346 
347  return -1;
348 }
349 
350 std::ostream& operator<<(std::ostream& os, const XsecRecord& xd) {
351  os << "Species: " << xd.Species() << std::endl;
352  return os;
353 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
The VectorView class.
Definition: matpackI.h:610
A class implementing complex numbers for ARTS.
Index nelem() const
Number of elements.
Definition: array.h:195
#define x0
std::ostream & operator<<(std::ostream &os, const XsecRecord &xd)
Definition: hitran_xsec.cc:350
The Vector class.
Definition: matpackI.h:860
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
ArrayOfVector mfgrids
Definition: hitran_xsec.h:120
The range class.
Definition: matpackI.h:160
Numeric lorentz_pdf(const Numeric x, const Numeric x0, const Numeric gamma)
Definition: hitran_xsec.cc:48
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
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
String species_name_from_species_index(const Index spec_ind)
Return species name for given species index.
Definition: absorption.cc:569
ArrayOfVector mtintersect
Definition: hitran_xsec.h:123
const Numeric * get_c_array() const
Conversion to plain C-array, const-version.
Definition: matpackI.cc:272
The global header file for ARTS.
_CS_string_type str() const
Definition: sstream.h:491
#define CREATE_OUTS
Definition: messages.h:209
Methods and classes for HITRAN absorption cross section data.
Index mspecies
Definition: hitran_xsec.h:116
Vector mrefpressure
Definition: hitran_xsec.h:118
#define temp
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Vector mcoeffs
Definition: hitran_xsec.h:117
The Matrix class.
Definition: matpackI.h:1193
Declarations required for the calculation of absorption coefficients.
void Extract(VectorView result, ConstVectorView f_grid, const Numeric &pressure, const Numeric &temperature, const Index &apply_tfit, const Verbosity &verbosity) const
Interpolate cross section data.
Definition: hitran_xsec.cc:151
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.
void convolve(Vector &result, const ConstVectorView &xsec, const ConstVectorView &lorentz)
Definition: hitran_xsec.cc:59
const Numeric PI
ArrayOfVector mxsecs
Definition: hitran_xsec.h:121
Index Species() const
Return species index.
Definition: hitran_xsec.h:40
Header file for interpolation_poly.cc.
Numeric func_2straights(const Numeric x, const Vector &coeffs)
Definition: hitran_xsec.cc:42
Vector mreftemperature
Definition: hitran_xsec.h:119
A constant view of a Vector.
Definition: matpackI.h:476
Index nelem(const Lines &l)
Number of lines.
Index hitran_xsec_get_index(const ArrayOfXsecRecord &xsec_data, const Index species)
Get the index in hitran_xsec_data for the given species.
Definition: hitran_xsec.cc:342
String SpeciesName() const
Return species name.
Definition: hitran_xsec.cc:53
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
ArrayOfVector mtslope
Definition: hitran_xsec.h:122