ARTS 2.5.4 (git: 4c0d3b4d)
m_hitran_xsec.cc
Go to the documentation of this file.
1/* Copyright (C) 2017 Oliver Lemke <oliver.lemke@uni-hamburg.de> and
2 Stefan Buehler <stefan.buehler@uni-hamburg.de>.
3
4 This program is free software; you can redistribute it and/or modify it
5 under the terms of the GNU General Public License as published by the
6 Free Software Foundation; either version 2, or (at your option) any
7 later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 USA. */
18
28#include "arts.h"
29#include "global_data.h"
30#include "hitran_xsec.h"
31#include "jacobian.h"
32#include "m_xml.h"
33#include "messages.h"
34#include "physics_funcs.h"
35
36/* Workspace method: Doxygen documentation will be auto-generated */
37void ReadXsecData(ArrayOfXsecRecord& hitran_xsec_data,
38 const ArrayOfArrayOfSpeciesTag& abs_species,
39 const String& basename,
40 const Verbosity& verbosity) {
41 // Build a set of species indices. Duplicates are ignored.
42 std::set<Species::Species> unique_species;
43 for (auto& asp : abs_species) {
44 for (auto& sp : asp) {
45 if (sp.Type() == Species::TagType::HitranXsec) {
46 unique_species.insert(sp.Spec());
47 }
48 }
49 }
50
51 String tmpbasename = basename;
52 if (basename.length() && basename[basename.length() - 1] != '/') {
53 tmpbasename += '.';
54 }
55
56 // Read xsec data for all active species and collect them in hitran_xsec_data
57 hitran_xsec_data.clear();
58 for (auto& species_name : unique_species) {
59 XsecRecord xsec_coeffs;
60 const String filename{tmpbasename +
61 String(Species::toShortName(species_name)) + ".xml"};
62
63 try {
64 ReadXML(xsec_coeffs, "", filename, "", verbosity);
65
66 hitran_xsec_data.push_back(xsec_coeffs);
67 } catch (const std::exception& e) {
69 "Error reading coefficients file:\n", filename, "\n", e.what());
70 }
71 }
72}
73
74/* Workspace method: Doxygen documentation will be auto-generated */
76 PropagationMatrix& propmat_clearsky,
77 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
78 // WS Input:
79 const ArrayOfArrayOfSpeciesTag& abs_species,
80 const ArrayOfRetrievalQuantity& jacobian_quantities,
81 const Vector& f_grid,
82 const Numeric& rtp_pressure,
83 const Numeric& rtp_temperature,
84 const Vector& rtp_vmr,
85 const ArrayOfXsecRecord& hitran_xsec_data,
86 const Numeric& force_p,
87 const Numeric& force_t,
88 // Verbosity object:
89 const Verbosity& verbosity) {
91
92 // Forward simulations and their error handling
93 ARTS_USER_ERROR_IF(rtp_vmr.nelem() not_eq abs_species.nelem(),
94 "Mismatch dimensions on species and VMR inputs");
96 propmat_clearsky.NumberOfFrequencies() not_eq f_grid.nelem(),
97 "Mismatch dimensions on internal matrices of xsec and frequency");
98
99 // Derivatives and their error handling
100 if (dpropmat_clearsky_dx.nelem()) {
102 dpropmat_clearsky_dx.nelem() not_eq jacobian_quantities.nelem(),
103 "Mismatch dimensions on xsec derivatives and Jacobian grids");
105 std::any_of(dpropmat_clearsky_dx.cbegin(),
106 dpropmat_clearsky_dx.cend(),
107 [&f_grid](auto& x) {
108 return x.NumberOfFrequencies() not_eq f_grid.nelem();
109 }),
110 "Mismatch dimensions on internal matrices of xsec derivatives and frequency");
111 }
112
113 // Jacobian overhead START
114 /* NOTE: The calculations below are inefficient and could
115 be made much better by using interp in Extract to
116 return the derivatives as well. */
117 // Jacobian vectors START
118 Vector dxsec_temp_dT;
119 Vector dxsec_temp_dF;
120 Vector dfreq;
121 // Jacobian vectors END
122 const bool do_jac = supports_hitran_xsec(jacobian_quantities);
123 const bool do_freq_jac = do_frequency_jacobian(jacobian_quantities);
124 const bool do_temp_jac = do_temperature_jacobian(jacobian_quantities);
125 const Numeric df = frequency_perturbation(jacobian_quantities);
126 const Numeric dt = temperature_perturbation(jacobian_quantities);
127 if (do_freq_jac) {
128 dfreq.resize(f_grid.nelem());
129 dfreq = f_grid;
130 dfreq += df;
131 dxsec_temp_dF.resize(f_grid.nelem());
132 }
133 if (do_temp_jac) {
134 dxsec_temp_dT.resize(f_grid.nelem());
135 }
136 // Jacobian overhead END
137
138 // Useful if there is no Jacobian to calculate
139 ArrayOfMatrix empty;
140
141 // Allocate a vector with dimension frequencies for constructing our
142 // cross-sections before adding them (more efficient to allocate this here
143 // outside of the loops)
144 Vector xsec_temp(f_grid.nelem(), 0.);
145
146 ArrayOfString fail_msg;
147 bool do_abort = false;
148 // Loop over Xsec data sets.
149 // Index ii loops through the outer array (different tag groups),
150 // index s through the inner array (different tags within each goup).
151 for (Index i = 0; i < abs_species.nelem(); i++) {
152 for (Index s = 0; s < abs_species[i].nelem(); s++) {
153 const SpeciesTag& this_species = abs_species[i][s];
154
155 // Check if this is a HITRAN cross section tag
156 if (this_species.Type() != Species::TagType::HitranXsec) continue;
157
158 Index this_xdata_index =
159 hitran_xsec_get_index(hitran_xsec_data, this_species.Spec());
160 ARTS_USER_ERROR_IF(this_xdata_index < 0,
161 "Cross-section species ",
162 this_species.Name(),
163 " not found in *hitran_xsec_data*.")
164 const XsecRecord& this_xdata = hitran_xsec_data[this_xdata_index];
165 // ArrayOfMatrix& this_dxsec = do_jac ? dpropmat_clearsky_dx[i] : empty;
166
167 if (do_abort) continue;
168 const Numeric current_p = force_p < 0 ? rtp_pressure : force_p;
169 const Numeric current_t = force_t < 0 ? rtp_temperature : force_t;
170
171 // Get the absorption cross sections from the HITRAN data:
172 this_xdata.Extract(xsec_temp, f_grid, current_p, current_t, verbosity);
173 if (do_freq_jac) {
174 this_xdata.Extract(
175 dxsec_temp_dF, dfreq, current_p, current_t, verbosity);
176 }
177 if (do_temp_jac) {
178 this_xdata.Extract(
179 dxsec_temp_dT, f_grid, current_p, current_t + dt, verbosity);
180 }
181 }
182
183 // Add to result variable:
184 Numeric nd = number_density(rtp_pressure, rtp_temperature);
185 if (!do_jac) {
186 xsec_temp *= nd * rtp_vmr[i];
187 propmat_clearsky.Kjj() += xsec_temp;
188 } else {
189 Numeric dnd_dt = dnumber_density_dt(rtp_pressure, rtp_temperature);
190 for (Index f = 0; f < f_grid.nelem(); f++) {
191 propmat_clearsky.Kjj()[f] += xsec_temp[f] * nd * rtp_vmr[i];
192 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
193 const auto& deriv = jacobian_quantities[iq];
194
195 if (!deriv.propmattype()) continue;
196
197 if (is_frequency_parameter(deriv)) {
198 dpropmat_clearsky_dx[iq].Kjj()[f] +=
199 (dxsec_temp_dF[f] - xsec_temp[f]) * nd * rtp_vmr[i] / df;
200 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR ||
201 deriv == Jacobian::Line::VMR) {
202 if (species_match(deriv, abs_species[i])) {
203 dpropmat_clearsky_dx[iq].Kjj()[f] +=
204 xsec_temp[f] * nd * rtp_vmr[i];
205 }
206 } else if (deriv == Jacobian::Atm::Temperature) {
207 dpropmat_clearsky_dx[iq].Kjj()[f] +=
208 ((dxsec_temp_dT[f] - xsec_temp[f]) / dt * nd +
209 xsec_temp[f] * dnd_dt) *
210 rtp_vmr[i];
211 }
212 }
213 }
214 }
215 }
216}
217
218/* Workspace method: Doxygen documentation will be auto-generated */
220 ArrayOfMatrix& abs_xsec_per_species,
221 ArrayOfArrayOfMatrix& /* dabs_xsec_per_species_dx */,
222 // WS Input:
223 const ArrayOfArrayOfSpeciesTag& abs_species,
224 const ArrayOfRetrievalQuantity& /* jacobian_quantities */,
225 const ArrayOfIndex& abs_species_active,
226 const Vector& f_grid,
227 const Vector& abs_p,
228 const Vector& abs_t,
229 const ArrayOfXsecRecord& hitran_xsec_data,
230 const Numeric& force_p,
231 const Numeric& force_t,
232 // Verbosity object:
233 const Verbosity& verbosity) {
235
236 {
237 // Check that all parameters that should have the number of tag
238 // groups as a dimension are consistent:
239 const Index n_tgs = abs_species.nelem();
240 const Index n_xsec = abs_xsec_per_species.nelem();
241
243 n_tgs != n_xsec,
244 "The following variables must all have the same dimension:\n"
245 "abs_species: ",
246 n_tgs,
247 "\n"
248 "abs_xsec_per_species: ",
249 n_xsec)
250 }
251
252 // Useful if there is no Jacobian to calculate
253 ArrayOfMatrix empty;
254
255 {
256 // Check that all parameters that should have the the dimension of p_grid
257 // are consistent:
258 const Index n_p = abs_p.nelem();
259 const Index n_t = abs_t.nelem();
260
262 n_p != n_t,
263 "The following variables must all have the same dimension:\n"
264 "abs_p: ",
265 n_p,
266 "\n"
267 "abs_t: ",
268 n_t)
269 }
270
271 // Allocate a vector with dimension frequencies for constructing our
272 // cross-sections before adding them (more efficient to allocate this here
273 // outside of the loops)
274 Vector xsec_temp(f_grid.nelem(), 0.);
275
276 ArrayOfString fail_msg;
277 bool do_abort = false;
278
279 // Loop over Xsec data sets.
280 // Index ii loops through the outer array (different tag groups),
281 // index s through the inner array (different tags within each goup).
282 for (Index ii = 0; ii < abs_species_active.nelem(); ii++) {
283 const Index i = abs_species_active[ii];
284
285 for (Index s = 0; s < abs_species[i].nelem(); s++) {
286 const SpeciesTag& this_species = abs_species[i][s];
287
288 // Check if this is a HITRAN cross section tag
289 if (this_species.Type() != Species::TagType::HitranXsec) continue;
290
291 Index this_xdata_index =
292 hitran_xsec_get_index(hitran_xsec_data, this_species.Spec());
293 ARTS_USER_ERROR_IF(this_xdata_index < 0,
294 "Cross-section species ",
295 this_species.Name(),
296 " not found in *hitran_xsec_data*.")
297 const XsecRecord& this_xdata = hitran_xsec_data[this_xdata_index];
298 Matrix& this_xsec = abs_xsec_per_species[i];
299
300 // Loop over pressure:
301#pragma omp parallel for if (!arts_omp_in_parallel() && abs_p.nelem() >= 1) \
302 firstprivate(xsec_temp)
303 for (Index ip = 0; ip < abs_p.nelem(); ip++) {
304 if (do_abort) continue;
305 const Numeric current_p = force_p < 0 ? abs_p[ip] : force_p;
306 const Numeric current_t = force_t < 0 ? abs_t[ip] : force_t;
307
308 // Get the absorption cross sections from the HITRAN data:
309 try {
310 this_xdata.Extract(xsec_temp,
311 f_grid,
312 current_p,
313 current_t,
314 verbosity);
315 } catch (runtime_error& e) {
316 ostringstream os;
317 os << "Problem with HITRAN cross section species "
318 << this_species.Name() << " at pressure level " << ip << " ("
319 << abs_p[ip] / 100. << " hPa):\n"
320 << e.what() << "\n";
321#pragma omp critical(abs_xsec_per_speciesAddHitranXsec)
322 {
323 do_abort = true;
324 fail_msg.push_back(os.str());
325 }
326 }
327 this_xsec(joker, ip) += xsec_temp;
328 }
329 }
330 }
331
332 if (do_abort) {
333 ARTS_USER_ERROR(fail_msg);
334 }
335}
The global header file for ARTS.
This can be used to make arrays out of anything.
Definition: array.h:108
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:541
The Matrix class.
Definition: matpackI.h:1270
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
The Vector class.
Definition: matpackI.h:908
void resize(Index n)
Resize function.
Definition: matpackI.cc:411
Hitran crosssection class.
Definition: hitran_xsec.h:43
void Extract(VectorView result, const Vector &f_grid, Numeric pressure, Numeric temperature, const Verbosity &verbosity) const
Calculate hitran cross section data.
Definition: hitran_xsec.cc:52
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
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:226
Methods and classes for HITRAN absorption cross section data.
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a frequency derivative.
Definition: jacobian.cc:1166
bool species_match(const RetrievalQuantity &rq, const ArrayOfSpeciesTag &ast)
Returns if the Retrieval quantity is VMR derivative for all the species in the species tags.
Definition: jacobian.cc:1114
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1145
bool supports_hitran_xsec(const ArrayOfRetrievalQuantity &js)
Returns if the array supports HITRAN cross-section derivatives.
Definition: jacobian.cc:1088
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Definition: jacobian.cc:1001
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the temperature perturbation if it exists.
Definition: jacobian.cc:1174
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
Definition: jacobian.cc:1182
Routines for setting up the jacobian.
void abs_xsec_per_speciesAddHitranXsec(ArrayOfMatrix &abs_xsec_per_species, ArrayOfArrayOfMatrix &, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfRetrievalQuantity &, const ArrayOfIndex &abs_species_active, const Vector &f_grid, const Vector &abs_p, const Vector &abs_t, const ArrayOfXsecRecord &hitran_xsec_data, const Numeric &force_p, const Numeric &force_t, const Verbosity &verbosity)
WORKSPACE METHOD: abs_xsec_per_speciesAddHitranXsec.
void propmat_clearskyAddHitranXsec(PropagationMatrix &propmat_clearsky, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &f_grid, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Vector &rtp_vmr, const ArrayOfXsecRecord &hitran_xsec_data, const Numeric &force_p, const Numeric &force_t, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddHitranXsec.
void ReadXsecData(ArrayOfXsecRecord &hitran_xsec_data, const ArrayOfArrayOfSpeciesTag &abs_species, const String &basename, const Verbosity &verbosity)
WORKSPACE METHOD: ReadXsecData.
Workspace methods and template functions for supergeneric XML IO.
void ReadXML(T &v, const String &v_name, const String &f, const String &f_name, const Verbosity &verbosity)
WORKSPACE METHOD: ReadXML.
Definition: m_xml.h:41
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
Declarations having to do with the four output streams.
#define CREATE_OUTS
Definition: messages.h:209
my_basic_string< char > String
The String type for ARTS.
Definition: mystring.h:287
Particulates ENUMCLASS(Line, char, VMR, Strength, Center, ShapeG0X0, ShapeG0X1, ShapeG0X2, ShapeG0X3, ShapeD0X0, ShapeD0X1, ShapeD0X2, ShapeD0X3, ShapeG2X0, ShapeG2X1, ShapeG2X2, ShapeG2X3, ShapeD2X0, ShapeD2X1, ShapeD2X2, ShapeD2X3, ShapeFVCX0, ShapeFVCX1, ShapeFVCX2, ShapeFVCX3, ShapeETAX0, ShapeETAX1, ShapeETAX2, ShapeETAX3, ShapeYX0, ShapeYX1, ShapeYX2, ShapeYX3, ShapeGX0, ShapeGX1, ShapeGX2, ShapeGX3, ShapeDVX0, ShapeDVX1, ShapeDVX2, ShapeDVX3, ECS_SCALINGX0, ECS_SCALINGX1, ECS_SCALINGX2, ECS_SCALINGX3, ECS_BETAX0, ECS_BETAX1, ECS_BETAX2, ECS_BETAX3, ECS_LAMBDAX0, ECS_LAMBDAX1, ECS_LAMBDAX2, ECS_LAMBDAX3, ECS_DCX0, ECS_DCX1, ECS_DCX2, ECS_DCX3, NLTE) static_assert(Index(Line ArrayOfSpeciesTagVMR
Definition: jacobian.h:101
Temperature
Definition: jacobian.h:58
This file contains declerations of functions of physical character.
constexpr Numeric dnumber_density_dt(Numeric p, Numeric t) noexcept
dnumber_density_dT
Definition: physics_funcs.h:84
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition: physics_funcs.h:70
Species::Tag SpeciesTag
Definition: species_tags.h:101