ARTS 2.5.4 (git: 31ce4f0e)
xsec_fit.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 "xsec_fit.h"
27
28#include <algorithm>
29#include <numeric>
30
31#include "absorption.h"
32#include "arts.h"
33#include "check_input.h"
34#include "debug.h"
35#include "interpolation.h"
36#include "physics_funcs.h"
37
39 // The function species_name_from_species_index internally does an assertion
40 // that the species with this index really exists.
41 return Species::toShortName(mspecies);
42}
43
44void XsecRecord::SetVersion(const Index version) {
45 if (version != 2) {
46 ARTS_USER_ERROR("Invalid version, only 2 supported")
47 }
48
49 mversion = version;
50}
51
53 const Vector& f_grid,
54 const Numeric pressure,
55 const Numeric temperature,
56 const Verbosity& verbosity) const {
58
59 const Index nf = f_grid.nelem();
60
61 ARTS_ASSERT(result.nelem() == nf)
62
63 result = 0.;
64
65 const Index ndatasets = mfitcoeffs.nelem();
66 for (Index this_dataset_i = 0; this_dataset_i < ndatasets; this_dataset_i++) {
67 const Vector& data_f_grid = mfitcoeffs[this_dataset_i].get_numeric_grid(0);
68 const Numeric data_fmin = data_f_grid[0];
69 const Numeric data_fmax = data_f_grid[data_f_grid.nelem() - 1];
70
71 if (out3.sufficient_priority()) {
72 ostringstream os;
73 os << " f_grid: " << f_grid[0] << " - " << f_grid[nf - 1]
74 << " Hz\n"
75 << " data_f_grid: " << data_fmin << " - " << data_fmax << " Hz\n"
76 << " pressure: " << pressure << " K\n";
77 out3 << os.str();
78 }
79
80 // We want to return result zero for all f_grid points that are outside the
81 // data_f_grid, because xsec datasets are defined only where the absorption
82 // was measured. So, we have to find out which part of f_grid is inside
83 // data_f_grid.
84 const Numeric* f_grid_begin = f_grid.get_c_array();
85 const Numeric* f_grid_end = f_grid_begin + f_grid.nelem();
86 const Index i_fstart = std::distance(
87 f_grid_begin, std::lower_bound(f_grid_begin, f_grid_end, data_fmin));
88 const Index i_fstop =
89 std::distance(
90 f_grid_begin,
91 std::upper_bound(f_grid_begin + i_fstart, f_grid_end, data_fmax)) -
92 1;
93
94 // Ignore band if all frequencies are below or above data_f_grid:
95 if (i_fstart == nf || i_fstop == -1) continue;
96
97 const Index f_extent = i_fstop - i_fstart + 1;
98
99 if (out3.sufficient_priority()) {
100 ostringstream os;
101 os << " " << f_extent << " frequency extraction points starting at "
102 << "frequency index " << i_fstart << ".\n";
103 out3 << os.str();
104 }
105
106 // If f_extent is less than one, then the entire data_f_grid is between two
107 // grid points of f_grid. (So that we do not have any f_grid points inside
108 // data_f_grid.) Return also in this case.
109 if (f_extent < 1) continue;
110
111 // This is the part of f_grid that lies inside the xsec band
112 const ConstVectorView f_grid_active = f_grid[Range(i_fstart, f_extent)];
113
114 // Determine the index of the first grid points in the xsec band that lie
115 // outside or on the boundary of the part of f_grid that lies inside the
116 // xsec band. We can ignore the remaining data.
117 const Numeric f_grid_fmin = f_grid[i_fstart];
118 const Numeric f_grid_fmax = f_grid[i_fstop];
119
120 const Numeric* data_f_grid_begin = data_f_grid.get_c_array();
121 const Numeric* data_f_grid_end = data_f_grid_begin + data_f_grid.size() - 1;
122 const Index i_data_fstart =
123 std::distance(
124 data_f_grid_begin,
125 std::upper_bound(data_f_grid_begin, data_f_grid_end, f_grid_fmin)) -
126 1;
127 const Index i_data_fstop = std::distance(
128 data_f_grid_begin,
129 std::upper_bound(
130 data_f_grid_begin + i_data_fstart, data_f_grid_end, f_grid_fmax));
131
132 ARTS_ASSERT(i_data_fstart >= 0)
133 ARTS_ASSERT(f_grid[i_fstart] > data_f_grid[i_data_fstart])
134 ARTS_ASSERT(f_grid[i_fstart] < data_f_grid[i_data_fstart + 1])
135 ARTS_ASSERT(f_grid[i_fstop] < data_f_grid[i_data_fstop])
136 ARTS_ASSERT(f_grid[i_fstop] > data_f_grid[i_data_fstop - 1])
137
138 // Extent for active data frequency vector:
139 const Index data_f_extent = i_data_fstop - i_data_fstart + 1;
140
141 // This is the part of the xsec dataset for which we have to do the
142 // interpolation.
143 const Range active_range(i_data_fstart, data_f_extent);
144 const ConstVectorView data_f_grid_active = data_f_grid[active_range];
145
146 Vector fit_result(data_f_grid.nelem());
147 VectorView fit_result_active = fit_result[active_range];
148
149 // We have to create a matching view on the result vector:
150 VectorView result_active = result[Range(i_fstart, f_extent)];
151 Vector xsec_interp(f_extent);
152
153 CalcXsec(fit_result, this_dataset_i, pressure, temperature);
154
155 RemoveNegativeXsec(fit_result);
156
157 // Check if frequency is inside the range covered by the data:
158 chk_interpolation_grids("Frequency interpolation for cross sections",
159 data_f_grid,
160 f_grid_active);
161
162 {
163 // Find frequency grid positions:
164 ArrayOfGridPos f_gp(f_grid_active.nelem());
165 gridpos(f_gp, data_f_grid_active, f_grid_active);
166
167 Matrix itw(f_gp.nelem(), 2);
168 interpweights(itw, f_gp);
169 interp(xsec_interp, itw, fit_result_active, f_gp);
170 }
171
172 result_active += xsec_interp;
173 }
174}
175
177 const Index dataset,
178 const Numeric pressure,
179 const Numeric temperature) const {
180 for (Index i = 0; i < xsec.nelem(); i++) {
181 const ConstVectorView coeffs = mfitcoeffs[dataset].data(i, joker);
182 xsec[i] = coeffs[P00] + coeffs[P10] * temperature + coeffs[P01] * pressure +
183 coeffs[P20] * temperature * temperature;
184 }
185}
186
187// void XsecRecord::CalcDT(VectorView& xsec_dt,
188// const Index dataset,
189// const Numeric temperature) const {
190// for (Index i = 0; i < xsec_dt.nelem(); i++) {
191// const ConstVectorView coeffs = mfitcoeffs[dataset].data(i, joker);
192// xsec_dt[i] = coeffs[P10] + 2. * coeffs[P20] * temperature;
193// }
194// }
195
196// void XsecRecord::CalcDP(VectorView& xsec_dp,
197// const Index dataset,
198// const Numeric pressure) const {
199// for (Index i = 0; i < xsec_dp.nelem(); i++) {
200// const ConstVectorView coeffs = mfitcoeffs[dataset].data(i, joker);
201// xsec_dp[i] = coeffs[P01] + 2. * coeffs[P02] * pressure;
202// }
203// }
204
206 Numeric sum_xsec{};
207 Numeric sum_xsec_non_negative{};
208 for (auto& x : xsec) {
209 sum_xsec += x;
210 if (x < 0.) x = 0.;
211 sum_xsec_non_negative += x;
212 }
213
214 if (sum_xsec > 0. && sum_xsec != sum_xsec_non_negative) {
215 xsec *= sum_xsec / sum_xsec_non_negative;
216 }
217}
218
227 const Species::Species species) {
228 for (Index i = 0; i < xsec_data.nelem(); i++)
229 if (xsec_data[i].Species() == species) return i;
230
231 return -1;
232}
233
234std::ostream& operator<<(std::ostream& os, const XsecRecord& xd) {
235 os << "Species: " << xd.Species() << std::endl;
236 return os;
237}
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
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Vector.
Definition: matpackI.h:512
Numeric * get_c_array() const noexcept
Conversion to plain C-array, const-version.
Definition: matpackI.h:622
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
Index size() const noexcept
Definition: matpackI.h:539
The Matrix class.
Definition: matpackI.h:1261
The range class.
Definition: matpackI.h:159
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
Hitran crosssection class.
Definition: xsec_fit.h:43
Index mversion
Definition: xsec_fit.h:145
static const Index P00
Definition: xsec_fit.h:140
void CalcXsec(VectorView &xsec, Index dataset, Numeric pressure, Numeric temperature) const
Calculate crosssections.
Definition: xsec_fit.cc:176
Species::Species mspecies
Definition: xsec_fit.h:146
ArrayOfGriddedField2 mfitcoeffs
Definition: xsec_fit.h:152
String SpeciesName() const
Return species name.
Definition: xsec_fit.cc:38
void RemoveNegativeXsec(VectorView &xsec) const
Remove negative cross sections and adjust integral.
Definition: xsec_fit.cc:205
void SetVersion(Index version)
Set species name.
Definition: xsec_fit.cc:44
static const Index P20
Definition: xsec_fit.h:143
static const Index P10
Definition: xsec_fit.h:141
void Extract(VectorView result, const Vector &f_grid, Numeric pressure, Numeric temperature, const Verbosity &verbosity) const
Calculate hitran cross section data.
Definition: xsec_fit.cc:52
static const Index P01
Definition: xsec_fit.h:142
const Species::Species & Species() const
Return species index.
Definition: xsec_fit.h:46
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Header file for interpolation.cc.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
#define CREATE_OUTS
Definition: messages.h:208
This file contains declerations of functions of physical character.
std::ostream & operator<<(std::ostream &os, const XsecRecord &xd)
Definition: xsec_fit.cc:234
Index hitran_xsec_get_index(const ArrayOfXsecRecord &xsec_data, const Species::Species species)
Get the index in xsec_fit_data for the given species.
Definition: xsec_fit.cc:226
Methods and classes for HITRAN absorption cross section data.