ARTS 2.5.0 (git: 9ee3ac6c)
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"
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
40extern const Numeric PI;
41
42Numeric func_2straights(const Numeric x, const Vector& coeffs) {
43 ARTS_ASSERT(coeffs.nelem() == 3);
44 return (x <= coeffs[0]) ? coeffs[1] * x
45 : coeffs[2] * (x - coeffs[0]) + coeffs[1] * coeffs[0];
46}
47
48Numeric 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.
56 return Species::toShortName(mspecies);
57}
58
59void convolve(Vector& result,
60 const ConstVectorView& xsec,
61 const ConstVectorView& lorentz) {
62 Index n_xsec = xsec.nelem();
63 Index n_lorentz = lorentz.nelem();
64 // ARTS_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
80void 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 {
158
159 const Index nf = f_grid.nelem();
160
161 // Assert that result vector has right size:
162 ARTS_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 constexpr 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 // Find frequency grid positions:
271 const auto f_lag = Interpolation::FixedLagrangeVector<f_order>(f_grid_active, data_f_grid_active);
272
273 if (pressure > mrefpressure[this_dataset_i] &&
274 mrefpressure[this_dataset_i] > 0.) {
275 // Apply pressure dependent broadening and set negative values to zero.
276 // (These could happen due to overshooting of the higher order interpolation.)
277 const Numeric pdiff = pressure - mrefpressure[this_dataset_i];
278 const Numeric fwhm = func_2straights(pdiff, mcoeffs);
279 // std::cout << mcoeffs << " - ";
280 // std::cout << "pdiff: " << pdiff << " - fwhm: " << fwhm << " - fstep: "
281 // << f_grid[i_fstart] + f_grid[i_fstart + 1] << std::endl;
282
283 Vector f_lorentz(data_f_extent);
284 Numeric lsum = 0.;
285 for (Index i = 0; i < data_f_extent; i++) {
286 f_lorentz[i] =
287 lorentz_pdf(data_f_grid[i_data_fstart + i],
288 data_f_grid[i_data_fstart + data_f_extent / 2 - 1],
289 fwhm / 2.);
290 lsum += f_lorentz[i];
291 }
292
293 f_lorentz /= lsum;
294
295 Vector data_result(xsec_active.nelem());
296#ifdef ENABLE_FFTW
297 fftconvolve(
298 data_result,
299 xsec_active,
300 f_lorentz[Range(f_lorentz.nelem() / 4, f_lorentz.nelem() / 2, 1)]);
301#else
302 convolve(
303 data_result,
304 xsec_active,
305 f_lorentz[Range(f_lorentz.nelem() / 4, f_lorentz.nelem() / 2, 1)]);
306#endif /* ENABLE_FFTW */
307
308 // TODO: Add to result_active here
309 // Check if frequency is inside the range covered by the data:
310 chk_interpolation_grids("Frequency interpolation for cross sections",
311 data_f_grid,
312 f_grid_active,
313 f_order);
314
315 xsec_interp = reinterp(data_result, interpweights(f_lag), f_lag);
316 } else {
317 xsec_interp = reinterp(xsec_active, interpweights(f_lag), f_lag);
318 }
319
320 result_active += xsec_interp;
321 }
322}
323
332 const Species::Species species) {
333 for (Index i = 0; i < xsec_data.nelem(); i++)
334 if (xsec_data[i].Species() == species) return i;
335
336 return -1;
337}
338
339std::ostream& operator<<(std::ostream& os, const XsecRecord& xd) {
340 os << "Species: " << xd.Species() << std::endl;
341 return os;
342}
Declarations required for the calculation of absorption coefficients.
The global header file for ARTS.
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
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
A constant view of a Vector.
Definition: matpackI.h:489
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The range class.
Definition: matpackI.h:165
The VectorView class.
Definition: matpackI.h:626
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array, const-version.
Definition: matpackI.cc:271
The Vector class.
Definition: matpackI.h:876
ArrayOfVector mtintersect
Definition: hitran_xsec.h:123
Vector mrefpressure
Definition: hitran_xsec.h:118
Species::Species mspecies
Definition: hitran_xsec.h:116
Vector mcoeffs
Definition: hitran_xsec.h:117
String SpeciesName() const
Return species name.
Definition: hitran_xsec.cc:53
Vector mreftemperature
Definition: hitran_xsec.h:119
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
ArrayOfVector mxsecs
Definition: hitran_xsec.h:121
ArrayOfVector mtslope
Definition: hitran_xsec.h:122
Species::Species Species() const
Return species index.
Definition: hitran_xsec.h:40
ArrayOfVector mfgrids
Definition: hitran_xsec.h:120
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
Numeric func_2straights(const Numeric x, const Vector &coeffs)
Definition: hitran_xsec.cc:42
std::ostream & operator<<(std::ostream &os, const XsecRecord &xd)
Definition: hitran_xsec.cc:339
Index hitran_xsec_get_index(const ArrayOfXsecRecord &xsec_data, const Species::Species species)
Get the index in hitran_xsec_data for the given species.
Definition: hitran_xsec.cc:331
void convolve(Vector &result, const ConstVectorView &xsec, const ConstVectorView &lorentz)
Definition: hitran_xsec.cc:59
const Numeric PI
Numeric lorentz_pdf(const Numeric x, const Numeric x0, const Numeric gamma)
Definition: hitran_xsec.cc:48
Methods and classes for HITRAN absorption cross section data.
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
#define CREATE_OUTS
Definition: messages.h:209
Index nelem(const Lines &l)
Number of lines.
void reinterp(VectorView out, const ConstVectorView &iy, const Grid< Vector, 1 > &iw, const Array< Lagrange > &dim0)
void sum(PropagationMatrix &pm, const ComplexVectorView &abs, const PolarizationVector &polvec, const bool do_phase) ARTS_NOEXCEPT
Sums the Zeeman components into a propagation matrix.
Definition: zeemandata.cc:375