ARTS 2.5.0 (git: 9ee3ac6c)
m_cia.cc
Go to the documentation of this file.
1/* Copyright (C) 2012 Oliver Lemke <olemke@core-dump.info> and Stefan
2 Buehler <sbuehler@ltu.se>.
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 "auto_md.h"
32#include "cia.h"
33#include "file.h"
34#include "messages.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 */
41void abs_xsec_per_speciesAddCIA( // WS Output:
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 Matrix& abs_vmrs,
52 const ArrayOfCIARecord& abs_cia_data,
53 // WS Generic Input:
54 const Numeric& T_extrapolfac,
55 const Index& robust,
56 // Verbosity object:
57 const Verbosity& verbosity) {
59
60 {
61 // Check that all parameters that should have the number of tag
62 // groups as a dimension are consistent:
63 const Index n_tgs = abs_species.nelem();
64 const Index n_xsec = abs_xsec_per_species.nelem();
65 const Index nr_vmrs = abs_vmrs.nrows();
66 // const Index n_cia = abs_cia_data.nelem();
67
68 ARTS_USER_ERROR_IF (n_tgs != n_xsec || n_tgs != nr_vmrs,
69 "The following variables must all have the same dimension:\n"
70 "abs_species: ", n_tgs, "\n"
71 "abs_xsec_per_species: ", n_xsec, "\n"
72 "abs_vmrs.nrows: ", nr_vmrs, "\n")
73 }
74
75 // Jacobian overhead START
76 /* NOTE: The calculations below are inefficient and could
77 be made much better by using interp in Extract to
78 return the derivatives as well. */
79 const bool do_jac = supports_CIA(jacobian_quantities);
80 const bool do_freq_jac = do_frequency_jacobian(jacobian_quantities);
81 const bool do_temp_jac = do_temperature_jacobian(jacobian_quantities);
82 Vector dfreq, dabs_t;
83 const Numeric df = frequency_perturbation(jacobian_quantities);
84 const Numeric dt = temperature_perturbation(jacobian_quantities);
85
86 if (do_freq_jac) {
87 dfreq = f_grid;
88 dfreq += df;
89 }
90 if (do_temp_jac) {
91 dabs_t.resize(abs_t.nelem());
92 dabs_t = abs_t;
93 dabs_t += dt;
94 }
95 // Jacobian overhead END
96
97 // Useful if there is no Jacobian to calculate
98 ArrayOfMatrix empty;
99
100 {
101 // Check that all parameters that should have the the dimension of p_grid
102 // are consistent:
103 const Index n_p = abs_p.nelem();
104 const Index n_t = abs_t.nelem();
105 const Index nc_vmrs = abs_vmrs.ncols();
106
107 ARTS_USER_ERROR_IF (n_p != n_t || n_p != nc_vmrs,
108 "The following variables must all have the same dimension:\n"
109 "abs_p: ", n_p, "\n"
110 "abs_t: ", n_t, "\n"
111 "abs_vmrs.ncols: ", nc_vmrs)
112 }
113
114 // Allocate a vector with dimension frequencies for constructing our
115 // cross-sections before adding them (more efficient to allocate this here
116 // outside of the loops)
117 Vector xsec_temp(f_grid.nelem());
118
119 // Jacobian vectors START
120 Vector dxsec_temp_dT;
121 Vector dxsec_temp_dF;
122 if (do_freq_jac) dxsec_temp_dF.resize(f_grid.nelem());
123 if (do_temp_jac) dxsec_temp_dT.resize(f_grid.nelem());
124 // Jacobian vectors END
125
126 // Loop over CIA data sets.
127 // Index ii loops through the outer array (different tag groups),
128 // index s through the inner array (different tags within each goup).
129 for (Index ii = 0; ii < abs_species_active.nelem(); ii++) {
130 const Index i = abs_species_active[ii];
131
132 for (Index s = 0; s < abs_species[i].nelem(); s++) {
133 const SpeciesTag& this_species = abs_species[i][s];
134
135 // Check if this is a CIA tag
136 if (this_species.Type() != Species::TagType::Cia) continue;
137
138 // Get convenient references of this CIA data record and this
139 // absorption cross-section record:
140 Index this_cia_index = cia_get_index(
141 abs_cia_data, this_species.Spec(), this_species.cia_2nd_species);
142
143 ARTS_ASSERT(this_cia_index != -1);
144
145 const CIARecord& this_cia = abs_cia_data[this_cia_index];
146 Matrix& this_xsec = abs_xsec_per_species[i];
147
148 if (out2.sufficient_priority()) {
149 // Some nice output to out2:
150 out2 << " CIA Species found: " + this_species.Name() + "\n";
151 }
152
153 // Check that the dimension of this_xsec is
154 // consistent with abs_p and f_grid.
155 ARTS_USER_ERROR_IF (this_xsec.nrows() != f_grid.nelem(),
156 "Wrong dimension of abs_xsec_per_species.nrows for species ", i,
157 ":\n"
158 "should match f_grid (", f_grid.nelem(), ") but is ",
159 this_xsec.nrows(), ".")
160 ARTS_USER_ERROR_IF (this_xsec.ncols() != abs_p.nelem(),
161 "Wrong dimension of abs_xsec_per_species.ncols for species ", i,
162 ":\n"
163 "should match abs_p (", abs_p.nelem(), ") but is ",
164 this_xsec.ncols(), ".")
165
166 // Find out index of VMR for the second CIA species.
167 // (The index for the first species is simply i.)
168 Index i_sec = find_first_species(abs_species, this_cia.Species(1));
169
170 // Catch the case that the VMR for the second species does not exist:
171 ARTS_USER_ERROR_IF (i_sec < 0,
172 "VMR profile for second species in CIA species pair does not exist.\n"
173 "Tag ", this_species.Name(), " needs a VMR profile of ",
174 this_cia.MoleculeName(1), "!")
175
176 // Loop over pressure:
177 for (Index ip = 0; ip < abs_p.nelem(); ip++) {
178 // Get the binary absorption cross sections from the CIA data:
179 try {
180 this_cia.Extract(xsec_temp,
181 f_grid,
182 abs_t[ip],
183 this_species.cia_dataset_index,
184 T_extrapolfac,
185 robust,
186 verbosity);
187 if (do_freq_jac)
188 this_cia.Extract(dxsec_temp_dF,
189 dfreq,
190 abs_t[ip],
191 this_species.cia_dataset_index,
192 T_extrapolfac,
193 robust,
194 verbosity);
195 if (do_temp_jac)
196 this_cia.Extract(dxsec_temp_dT,
197 f_grid,
198 dabs_t[ip],
199 this_species.cia_dataset_index,
200 T_extrapolfac,
201 robust,
202 verbosity);
203 } catch (const std::runtime_error& e) {
204 ARTS_USER_ERROR ("Problem with CIA species ",
205 this_species.Name(), ":\n", e.what())
206 }
207
208 // We have to multiply with the number density of the second CIA species.
209 // We do not have to multiply with the first, since we still
210 // want to return a (unary) absorption cross-section, not an
211 // absorption coefficient.
212
213 // Calculate number density from pressure and temperature.
214 const Numeric n =
215 abs_vmrs(i_sec, ip) * number_density(abs_p[ip], abs_t[ip]);
216
217 if (!do_jac) {
218 xsec_temp *= n;
219 // Add to result variable:
220 this_xsec(joker, ip) += xsec_temp;
221 } else {
222 const Numeric dn_dT =
223 abs_vmrs(i_sec, ip) * dnumber_density_dt(abs_p[ip], abs_t[ip]);
224
225 for (Index iv = 0; iv < xsec_temp.nelem(); iv++) {
226 this_xsec(iv, ip) += n * xsec_temp[iv];
227 for (Index iq = 0; iq < jacobian_quantities.nelem();
228 iq++) {
229 const auto& deriv = jacobian_quantities[iq];
230
231 if (not deriv.propmattype()) continue;
232
233 if (is_frequency_parameter(deriv))
234 dabs_xsec_per_species_dx[i][iq](iv, ip) +=
235 n * (dxsec_temp_dF[iv] - xsec_temp[iv]) / df;
236 else if (deriv == Jacobian::Atm::Temperature)
237 dabs_xsec_per_species_dx[i][iq](iv, ip) +=
238 n * (dxsec_temp_dT[iv] - xsec_temp[iv]) / dt +
239 xsec_temp[iv] * dn_dT;
240 else if (species_match(deriv, this_species.cia_2nd_species))
241 dabs_xsec_per_species_dx[i][iq](iv, ip) +=
242 number_density(abs_p[ip], abs_t[ip]) * xsec_temp[iv];
243 }
244 }
245 }
246 }
247 }
248 }
249}
250
251/* Workspace method: Doxygen documentation will be auto-generated */
252void CIARecordReadFromFile( // WS GOutput:
253 CIARecord& cia_record,
254 // WS Generic Input:
255 const String& species_tag,
256 const String& filename,
257 // Verbosity object:
258 const Verbosity& verbosity) {
259 SpeciesTag species(species_tag);
260
261 ARTS_USER_ERROR_IF (species.Type() != Species::TagType::Cia,
262 "Invalid species tag ", species_tag, ".\n"
263 "This is not recognized as a CIA type.\n")
264
265 cia_record.SetSpecies(species.Spec(), species.cia_2nd_species);
266 cia_record.ReadFromCIA(filename, verbosity);
267}
268
269/* Workspace method: Doxygen documentation will be auto-generated */
270void abs_cia_dataAddCIARecord( // WS Output:
271 ArrayOfCIARecord& abs_cia_data,
272 // WS GInput:
273 const CIARecord& cia_record,
274 const Index& clobber,
275 // WS Input:
276 const Verbosity&) {
277 Index cia_index =
278 cia_get_index(abs_cia_data, cia_record.Species(0), cia_record.Species(1));
279 if (cia_index == -1)
280 abs_cia_data.push_back(cia_record);
281 else if (clobber)
282 abs_cia_data[cia_index] = cia_record;
283 else
284 abs_cia_data[cia_index].AppendDataset(cia_record);
285}
286
287/* Workspace method: Doxygen documentation will be auto-generated */
288void abs_cia_dataReadFromCIA( // WS Output:
289 ArrayOfCIARecord& abs_cia_data,
290 // WS Input:
291 const ArrayOfArrayOfSpeciesTag& abs_species,
292 const String& catalogpath,
293 const Verbosity& verbosity) {
294 ArrayOfString subfolders;
295 subfolders.push_back("Main-Folder/");
296 subfolders.push_back("Alternate-Folder/");
297
298 abs_cia_data.resize(0);
299
300 // Loop species tag groups to find CIA tags.
301 // Index sp loops through the tag groups, index iso through the tags within
302 // each group. Despite the name, iso does not denote the isotope!
303 for (Index sp = 0; sp < abs_species.nelem(); sp++) {
304 for (Index iso = 0; iso < abs_species[sp].nelem(); iso++) {
305 if (abs_species[sp][iso].Type() != Species::TagType::Cia) continue;
306
307 ArrayOfString cia_names;
308
309 Index cia_index = cia_get_index(abs_cia_data,
310 abs_species[sp][iso].Spec(),
311 abs_species[sp][iso].cia_2nd_species);
312
313 // If cia_index is not -1, we have already read this datafile earlier
314 if (cia_index != -1) continue;
315
316 cia_names.push_back(String(Species::toShortName(abs_species[sp][iso].Spec())) +
317 "-" +
318 String(Species::toShortName(abs_species[sp][iso].cia_2nd_species)));
319
320 cia_names.push_back(
321 String(Species::toShortName(abs_species[sp][iso].cia_2nd_species)) +
322 "-" +
323 String(Species::toShortName(abs_species[sp][iso].Spec())));
324
325 ArrayOfString checked_dirs;
326
327 bool found = false;
328 for (Index fname = 0; !found && fname < cia_names.nelem(); fname++) {
329 String cia_name = cia_names[fname];
330
331 for (Index dir = 0; !found && dir < subfolders.nelem(); dir++) {
332 ArrayOfString files;
333 checked_dirs.push_back(catalogpath + "/" + subfolders[dir] +
334 cia_name + "/");
335 try {
336 list_directory(files, *(checked_dirs.end() - 1));
337 } catch (const std::runtime_error& e) {
338 continue;
339 }
340
341 for (Index i = files.nelem() - 1; i >= 0; i--) {
342 if (files[i].find(cia_name) != 0 ||
343 files[i].rfind(".cia") != files[i].length() - 4) {
344 files.erase(files.begin() + i);
345 }
346 }
347 if (files.nelem()) {
348 CIARecord ciar;
349
350 found = true;
351 String catfile = *(checked_dirs.end() - 1) + files[0];
352
353 ciar.SetSpecies(abs_species[sp][iso].Spec(),
354 abs_species[sp][iso].cia_2nd_species);
355 ciar.ReadFromCIA(catfile, verbosity);
356
357 abs_cia_data.push_back(ciar);
358 }
359 }
360 }
361
362 ARTS_USER_ERROR_IF (!found,
363 "Error: No data file found for CIA species ", cia_names[0],
364 "\n"
365 "Looked in directories: ", checked_dirs)
366 }
367 }
368}
369
370/* Workspace method: Doxygen documentation will be auto-generated */
371void abs_cia_dataReadFromXML( // WS Output:
372 ArrayOfCIARecord& abs_cia_data,
373 // WS Input:
374 const ArrayOfArrayOfSpeciesTag& abs_species,
375 const String& filename,
376 const Verbosity& verbosity) {
377 xml_read_from_file(filename, abs_cia_data, verbosity);
378
379 // Check that all CIA tags from abs_species are present in the
380 // XML file
381
382 vector<String> missing_tags;
383
384 // Loop species tag groups to find CIA tags.
385 // Index sp loops through the tag groups, index iso through the tags within
386 // each group. Despite the name, iso does not denote the isotope!
387 for (Index sp = 0; sp < abs_species.nelem(); sp++) {
388 for (Index iso = 0; iso < abs_species[sp].nelem(); iso++) {
389 if (abs_species[sp][iso].Type() != Species::TagType::Cia) continue;
390
391 Index cia_index = cia_get_index(abs_cia_data,
392 abs_species[sp][iso].Spec(),
393 abs_species[sp][iso].cia_2nd_species);
394
395 // If cia_index is -1, this CIA tag was not present in the input file
396 if (cia_index == -1) {
397 missing_tags.push_back(
398 String(Species::toShortName(abs_species[sp][iso].Spec())) +
399 "-" +
400 String(Species::toShortName(abs_species[sp][iso].cia_2nd_species)));
401 }
402 }
403 }
404
405 if (missing_tags.size()) {
406 ostringstream os;
407 bool first = true;
408
409 os << "Error: The following CIA tag(s) are missing in input file: ";
410 for (size_t i = 0; i < missing_tags.size(); i++) {
411 if (!first)
412 os << ", ";
413 else
414 first = false;
415 os << missing_tags[i];
416 }
417 ARTS_USER_ERROR (os.str());
418 }
419}
420
421/* Workspace method: Doxygen documentation will be auto-generated */
422void CIAInfo( // Generic Input:
423 const String& catalogpath,
424 const ArrayOfString& cia_tags,
425 const Verbosity& verbosity) {
427
428 ArrayOfArrayOfSpeciesTag species_tags;
429
430 for (Index i = 0; i < cia_tags.nelem(); i++) {
431 ArrayOfSpeciesTag this_species_tag;
432
433 ArrayOfString species_names;
434
435 cia_tags[i].split(species_names, "-");
436
437 ARTS_USER_ERROR_IF (species_names.nelem() != 2,
438 "ERROR: Cannot parse CIA tag: ", cia_tags[i])
439
440 this_species_tag.push_back(
441 SpeciesTag(species_names[0] + "-CIA-" + species_names[1] + "-0"));
442
443 species_tags.push_back(this_species_tag);
444 }
445
446 ArrayOfCIARecord cia_data;
447
448 abs_cia_dataReadFromCIA(cia_data, species_tags, catalogpath, verbosity);
449
450 Print(cia_data, 1, verbosity);
451}
Declarations required for the calculation of absorption coefficients.
The global header file for ARTS.
void Print(const T &in, const Index &level, const Verbosity &verbosity)
WORKSPACE METHOD: Print.
Definition: m_general.h:83
Index cia_get_index(const ArrayOfCIARecord &cia_data, const Species::Species sp1, const Species::Species sp2)
Get the index in cia_data for the two given species.
Definition: cia.cc:225
Header file for work with HITRAN collision induced absorption (CIA).
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
CIA data for a single pair of molecules.
Definition: cia.h:68
void SetSpecies(const Species::Species first, const Species::Species second)
Set CIA species.
Definition: cia.h:148
Species::Species Species(const Index i) const
Return CIA species index.
Definition: cia.h:97
void ReadFromCIA(const String &filename, const Verbosity &verbosity)
Read CIA catalog file.
Definition: cia.cc:301
void Extract(VectorView result, ConstVectorView f_grid, const Numeric &temperature, const Index &dataset, const Numeric &T_extrapolfac, const Index &robust, const Verbosity &verbosity) const
Vector version of extract.
Definition: cia.cc:237
String MoleculeName(const Index i) const
Return each molecule name (as a string) that is associated with this CIARecord.
Definition: cia.cc:265
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
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
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
void list_directory(ArrayOfString &files, String dirname)
Return list of files in directory.
Definition: file.cc:543
This file contains basic functions to handle ASCII files.
bool supports_CIA(const ArrayOfRetrievalQuantity &js)
Returns if the array supports CIA derivatives.
Definition: jacobian.cc:1084
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 do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1147
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:1176
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
Definition: jacobian.cc:1184
void abs_cia_dataAddCIARecord(ArrayOfCIARecord &abs_cia_data, const CIARecord &cia_record, const Index &clobber, const Verbosity &)
WORKSPACE METHOD: abs_cia_dataAddCIARecord.
Definition: m_cia.cc:270
void abs_cia_dataReadFromXML(ArrayOfCIARecord &abs_cia_data, const ArrayOfArrayOfSpeciesTag &abs_species, const String &filename, const Verbosity &verbosity)
WORKSPACE METHOD: abs_cia_dataReadFromXML.
Definition: m_cia.cc:371
void abs_xsec_per_speciesAddCIA(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 Matrix &abs_vmrs, const ArrayOfCIARecord &abs_cia_data, const Numeric &T_extrapolfac, const Index &robust, const Verbosity &verbosity)
WORKSPACE METHOD: abs_xsec_per_speciesAddCIA.
Definition: m_cia.cc:41
void abs_cia_dataReadFromCIA(ArrayOfCIARecord &abs_cia_data, const ArrayOfArrayOfSpeciesTag &abs_species, const String &catalogpath, const Verbosity &verbosity)
WORKSPACE METHOD: abs_cia_dataReadFromCIA.
Definition: m_cia.cc:288
const Numeric SPEED_OF_LIGHT
void CIAInfo(const String &catalogpath, const ArrayOfString &cia_tags, const Verbosity &verbosity)
WORKSPACE METHOD: CIAInfo.
Definition: m_cia.cc:422
void CIARecordReadFromFile(CIARecord &cia_record, const String &species_tag, const String &filename, const Verbosity &verbosity)
WORKSPACE METHOD: CIARecordReadFromFile.
Definition: m_cia.cc:252
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_OUT1
Definition: messages.h:205
#define CREATE_OUTS
Definition: messages.h:209
my_basic_string< char > String
The String type for ARTS.
Definition: mystring.h:287
Temperature
Definition: jacobian.h:57
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
Index find_first_species(const ArrayOfArrayOfSpeciesTag &specs, Species::Species spec) noexcept
Species::Tag SpeciesTag
Definition: species_tags.h:99
This file contains basic functions to handle XML data files.
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.h:160