ARTS 2.5.0 (git: 9ee3ac6c)
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
29#include "absorption.h"
30#include "arts.h"
31#include "messages.h"
32
33#include "auto_md.h"
34#include "hitran_xsec.h"
35#include "physics_funcs.h"
36#include "xml_io.h"
37
38extern const Numeric SPEED_OF_LIGHT;
39
40/* Workspace method: Doxygen documentation will be auto-generated */
42 ArrayOfMatrix& abs_xsec_per_species,
43 ArrayOfArrayOfMatrix& dabs_xsec_per_species_dx,
44 // WS Input:
45 const ArrayOfArrayOfSpeciesTag& abs_species,
46 const ArrayOfRetrievalQuantity& jacobian_quantities,
47 const ArrayOfIndex& abs_species_active,
48 const Vector& f_grid,
49 const Vector& abs_p,
50 const Vector& abs_t,
51 const ArrayOfXsecRecord& hitran_xsec_data,
52 const Index& apply_tfit,
53 const Numeric& force_p,
54 const Numeric& force_t,
55 // Verbosity object:
56 const Verbosity& verbosity) {
58
59 {
60 // Check that all parameters that should have the number of tag
61 // groups as a dimension are consistent:
62 const Index n_tgs = abs_species.nelem();
63 const Index n_xsec = abs_xsec_per_species.nelem();
64
65 ARTS_USER_ERROR_IF (n_tgs != n_xsec,
66 "The following variables must all have the same dimension:\n"
67 "abs_species: ", n_tgs, "\n"
68 "abs_xsec_per_species: ", n_xsec)
69 }
70
71 // Jacobian overhead START
72 /* NOTE: The calculations below are inefficient and could
73 be made much better by using interp in Extract to
74 return the derivatives as well. */
75 const bool do_jac = supports_hitran_xsec(jacobian_quantities);
76 const bool do_freq_jac = do_frequency_jacobian(jacobian_quantities);
77 Vector dfreq, dabs_t;
78 const Numeric df = frequency_perturbation(jacobian_quantities);
79 if (do_freq_jac) {
80 dfreq.resize(f_grid.nelem());
81 dfreq = f_grid;
82 dfreq += df;
83 }
84 // if (do_temp_jac)
85 // {
86 // dabs_t.resize(abs_t.nelem());
87 // dabs_t = abs_t;
88 // dabs_t += dt;
89 // }
90 // Jacobian overhead END
91
92 // Useful if there is no Jacobian to calculate
93 ArrayOfMatrix empty;
94
95 {
96 // Check that all parameters that should have the the dimension of p_grid
97 // are consistent:
98 const Index n_p = abs_p.nelem();
99 const Index n_t = abs_t.nelem();
100
101 ARTS_USER_ERROR_IF (n_p != n_t,
102 "The following variables must all have the same dimension:\n"
103 "abs_p: ", n_p, "\n"
104 "abs_t: ", n_t)
105 }
106
107 // Allocate a vector with dimension frequencies for constructing our
108 // cross-sections before adding them (more efficient to allocate this here
109 // outside of the loops)
110 Vector xsec_temp(f_grid.nelem(), 0.);
111
112 // Jacobian vectors START
113 // Vector dxsec_temp_dT;
114 Vector dxsec_temp_dF;
115 if (do_freq_jac) dxsec_temp_dF.resize(f_grid.nelem());
116 // if (do_temp_jac)
117 // dxsec_temp_dT.resize(f_grid.nelem());
118 // Jacobian vectors END
119
120 ArrayOfString fail_msg;
121 bool do_abort = false;
122
123 // Loop over Xsec data sets.
124 // Index ii loops through the outer array (different tag groups),
125 // index s through the inner array (different tags within each goup).
126 for (Index ii = 0; ii < abs_species_active.nelem(); ii++) {
127 const Index i = abs_species_active[ii];
128
129 for (Index s = 0; s < abs_species[i].nelem(); s++) {
130 const SpeciesTag& this_species = abs_species[i][s];
131
132 // Check if this is a HITRAN cross section tag
133 if (this_species.Type() != Species::TagType::HitranXsec) continue;
134
135#ifndef ENABLE_FFTW
136 out0 << "HITRAN XSEC Warning: No FFTW library support enabled, "
137 << "convolution will be extremely slow\n";
138#endif
139
140 Index this_xdata_index =
141 hitran_xsec_get_index(hitran_xsec_data, this_species.Spec());
142 ARTS_USER_ERROR_IF (this_xdata_index < 0,
143 "Cross-section species ", this_species.Name(),
144 " not found in *hitran_xsec_data*.")
145 const XsecRecord& this_xdata = hitran_xsec_data[this_xdata_index];
146 Matrix& this_xsec = abs_xsec_per_species[i];
147 ArrayOfMatrix& this_dxsec = do_jac ? dabs_xsec_per_species_dx[i] : empty;
148
149 // Loop over pressure:
150#pragma omp parallel for if (!arts_omp_in_parallel() && abs_p.nelem() >= 1) \
151 firstprivate(xsec_temp, dxsec_temp_dF)
152 for (Index ip = 0; ip < abs_p.nelem(); ip++) {
153 if (do_abort) continue;
154 const Numeric current_p = force_p < 0 ? abs_p[ip] : force_p;
155 const Numeric current_t = force_t < 0 ? abs_t[ip] : force_t;
156
157 // Get the absorption cross sections from the HITRAN data:
158 try {
159 this_xdata.Extract(
160 xsec_temp, f_grid, current_p, current_t, apply_tfit, verbosity);
161 if (do_freq_jac)
162 this_xdata.Extract(dxsec_temp_dF,
163 dfreq,
164 current_p,
165 current_t,
166 apply_tfit,
167 verbosity);
168 // FIXME: Temperature is not yet taken into account
169 // if(do_temp_jac)
170 // this_xdata.Extract(dxsec_temp_dT, f_grid, dabs_t[ip],
171 // verbosity);
172 } catch (runtime_error& e) {
173 ostringstream os;
174 os << "Problem with HITRAN cross section species "
175 << this_species.Name() << " at pressure level " << ip << " ("
176 << abs_p[ip] / 100. << " hPa):\n"
177 << e.what();
178#pragma omp critical(abs_xsec_per_speciesAddHitranXsec)
179 {
180 do_abort = true;
181 fail_msg.push_back(os.str());
182 }
183 }
184
185 if (!do_jac) {
186 // Add to result variable:
187 this_xsec(joker, ip) += xsec_temp;
188 } else {
189 for (Index iv = 0; iv < xsec_temp.nelem(); iv++) {
190 this_xsec(iv, ip) += xsec_temp[iv];
191 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
192 const auto& deriv = jacobian_quantities[ii];
193
194 if (not deriv.propmattype()) continue;
195
196 if (is_frequency_parameter(deriv))
197 this_dxsec[iq](iv, ip) +=
198 (dxsec_temp_dF[iv] - xsec_temp[iv]) / df;
199 else if (deriv == Jacobian::Line::VMR) {
200 if (species_match(deriv, abs_species[i])) {
201 this_dxsec[iq](iv, ip) += xsec_temp[iv];
202 }
203 }
204 }
205 }
206 }
207 }
208 }
209 }
210
211 if (do_abort) {
212 std::ostringstream os;
213 os << "Error messages from failures:\n";
214 for (const auto& msg : fail_msg) {
215 os << msg << '\n';
216 }
217 ARTS_USER_ERROR (os.str());
218 }
219}
Declarations required for the calculation of absorption coefficients.
The global header file for ARTS.
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
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The Matrix class.
Definition: matpackI.h:1225
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
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
#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:331
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:1168
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 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 frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
Definition: jacobian.cc:1184
const Numeric SPEED_OF_LIGHT
void abs_xsec_per_speciesAddHitranXsec(ArrayOfMatrix &abs_xsec_per_species, ArrayOfArrayOfMatrix &dabs_xsec_per_species_dx, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfIndex &abs_species_active, const Vector &f_grid, const Vector &abs_p, const Vector &abs_t, const ArrayOfXsecRecord &hitran_xsec_data, const Index &apply_tfit, const Numeric &force_p, const Numeric &force_t, const Verbosity &verbosity)
WORKSPACE METHOD: abs_xsec_per_speciesAddHitranXsec.
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
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUTS
Definition: messages.h:209
This file contains declerations of functions of physical character.
Species::Tag SpeciesTag
Definition: species_tags.h:99
This file contains basic functions to handle XML data files.