ARTS 2.5.4 (git: 31ce4f0e)
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 <algorithm>
30#include <filesystem>
31#include "absorption.h"
32#include "arts.h"
33#include "arts_constants.h"
34#include "auto_md.h"
35#include "cia.h"
36#include "debug.h"
37#include "file.h"
38#include "messages.h"
39#include "physics_funcs.h"
40#include "species.h"
41#include "species_tags.h"
42#include "xml_io.h"
43
45
46/* Workspace method: Doxygen documentation will be auto-generated */
47void abs_xsec_per_speciesAddCIA( // WS Output:
48 ArrayOfMatrix& abs_xsec_per_species,
49 ArrayOfArrayOfMatrix& dabs_xsec_per_species_dx,
50 // WS Input:
51 const ArrayOfArrayOfSpeciesTag& abs_species,
52 const ArrayOfRetrievalQuantity& jacobian_quantities,
53 const ArrayOfIndex& abs_species_active,
54 const Vector& f_grid,
55 const Vector& abs_p,
56 const Vector& abs_t,
57 const Matrix& abs_vmrs,
58 const ArrayOfCIARecord& abs_cia_data,
59 // WS Generic Input:
60 const Numeric& T_extrapolfac,
61 const Index& robust,
62 // Verbosity object:
63 const Verbosity& verbosity) {
65
66 {
67 // Check that all parameters that should have the number of tag
68 // groups as a dimension are consistent:
69 const Index n_tgs = abs_species.nelem();
70 const Index n_xsec = abs_xsec_per_species.nelem();
71 const Index nr_vmrs = abs_vmrs.nrows();
72 // const Index n_cia = abs_cia_data.nelem();
73
74 ARTS_USER_ERROR_IF (n_tgs != n_xsec || n_tgs != nr_vmrs,
75 "The following variables must all have the same dimension:\n"
76 "abs_species: ", n_tgs, "\n"
77 "abs_xsec_per_species: ", n_xsec, "\n"
78 "abs_vmrs.nrows: ", nr_vmrs, "\n")
79 }
80
81 // Jacobian overhead START
82 /* NOTE: The calculations below are inefficient and could
83 be made much better by using interp in Extract to
84 return the derivatives as well. */
85 const bool do_jac = supports_CIA(jacobian_quantities);
86 const bool do_freq_jac = do_frequency_jacobian(jacobian_quantities);
87 const bool do_temp_jac = do_temperature_jacobian(jacobian_quantities);
88 Vector dfreq, dabs_t;
89 const Numeric df = frequency_perturbation(jacobian_quantities);
90 const Numeric dt = temperature_perturbation(jacobian_quantities);
91
92 if (do_freq_jac) {
93 dfreq = f_grid;
94 dfreq += df;
95 }
96 if (do_temp_jac) {
97 dabs_t.resize(abs_t.nelem());
98 dabs_t = abs_t;
99 dabs_t += dt;
100 }
101 // Jacobian overhead END
102
103 // Useful if there is no Jacobian to calculate
104 ArrayOfMatrix empty;
105
106 {
107 // Check that all parameters that should have the the dimension of p_grid
108 // are consistent:
109 const Index n_p = abs_p.nelem();
110 const Index n_t = abs_t.nelem();
111 const Index nc_vmrs = abs_vmrs.ncols();
112
113 ARTS_USER_ERROR_IF (n_p != n_t || n_p != nc_vmrs,
114 "The following variables must all have the same dimension:\n"
115 "abs_p: ", n_p, "\n"
116 "abs_t: ", n_t, "\n"
117 "abs_vmrs.ncols: ", nc_vmrs)
118 }
119
120 // Allocate a vector with dimension frequencies for constructing our
121 // cross-sections before adding them (more efficient to allocate this here
122 // outside of the loops)
123 Vector xsec_temp(f_grid.nelem());
124
125 // Jacobian vectors START
126 Vector dxsec_temp_dT;
127 Vector dxsec_temp_dF;
128 if (do_freq_jac) dxsec_temp_dF.resize(f_grid.nelem());
129 if (do_temp_jac) dxsec_temp_dT.resize(f_grid.nelem());
130 // Jacobian vectors END
131
132 // Loop over CIA data sets.
133 // Index ii loops through the outer array (different tag groups),
134 // index s through the inner array (different tags within each goup).
135 for (Index ii = 0; ii < abs_species_active.nelem(); ii++) {
136 const Index i = abs_species_active[ii];
137
138 for (Index s = 0; s < abs_species[i].nelem(); s++) {
139 const SpeciesTag& this_species = abs_species[i][s];
140
141 // Check if this is a CIA tag
142 if (this_species.Type() != Species::TagType::Cia) continue;
143
144 // Get convenient references of this CIA data record and this
145 // absorption cross-section record:
146 Index this_cia_index = cia_get_index(
147 abs_cia_data, this_species.Spec(), this_species.cia_2nd_species);
148
149 ARTS_ASSERT(this_cia_index != -1);
150
151 const CIARecord& this_cia = abs_cia_data[this_cia_index];
152 Matrix& this_xsec = abs_xsec_per_species[i];
153
154 if (out2.sufficient_priority()) {
155 // Some nice output to out2:
156 out2 << " CIA Species found: " + this_species.Name() + "\n";
157 }
158
159 // Check that the dimension of this_xsec is
160 // consistent with abs_p and f_grid.
161 ARTS_USER_ERROR_IF (this_xsec.nrows() != f_grid.nelem(),
162 "Wrong dimension of abs_xsec_per_species.nrows for species ", i,
163 ":\n"
164 "should match f_grid (", f_grid.nelem(), ") but is ",
165 this_xsec.nrows(), ".")
166 ARTS_USER_ERROR_IF (this_xsec.ncols() != abs_p.nelem(),
167 "Wrong dimension of abs_xsec_per_species.ncols for species ", i,
168 ":\n"
169 "should match abs_p (", abs_p.nelem(), ") but is ",
170 this_xsec.ncols(), ".")
171
172 // Find out index of VMR for the second CIA species.
173 // (The index for the first species is simply i.)
174 Index i_sec = find_first_species(abs_species, this_cia.Species(1));
175
176 // Catch the case that the VMR for the second species does not exist:
177 ARTS_USER_ERROR_IF (i_sec < 0,
178 "VMR profile for second species in CIA species pair does not exist.\n"
179 "Tag ", this_species.Name(), " needs a VMR profile of ",
180 this_cia.MoleculeName(1), "!")
181
182 // Loop over pressure:
183 for (Index ip = 0; ip < abs_p.nelem(); ip++) {
184 // Get the binary absorption cross sections from the CIA data:
185 try {
186 this_cia.Extract(xsec_temp,
187 f_grid,
188 abs_t[ip],
189 this_species.cia_dataset_index,
190 T_extrapolfac,
191 robust,
192 verbosity);
193 if (do_freq_jac)
194 this_cia.Extract(dxsec_temp_dF,
195 dfreq,
196 abs_t[ip],
197 this_species.cia_dataset_index,
198 T_extrapolfac,
199 robust,
200 verbosity);
201 if (do_temp_jac)
202 this_cia.Extract(dxsec_temp_dT,
203 f_grid,
204 dabs_t[ip],
205 this_species.cia_dataset_index,
206 T_extrapolfac,
207 robust,
208 verbosity);
209 } catch (const std::runtime_error& e) {
210 ARTS_USER_ERROR ("Problem with CIA species ",
211 this_species.Name(), ":\n", e.what())
212 }
213
214 // We have to multiply with the number density of the second CIA species.
215 // We do not have to multiply with the first, since we still
216 // want to return a (unary) absorption cross-section, not an
217 // absorption coefficient.
218
219 // Calculate number density from pressure and temperature.
220 const Numeric n =
221 abs_vmrs(i_sec, ip) * number_density(abs_p[ip], abs_t[ip]);
222
223 if (!do_jac) {
224 xsec_temp *= n;
225 // Add to result variable:
226 this_xsec(joker, ip) += xsec_temp;
227 } else {
228 const Numeric dn_dT =
229 abs_vmrs(i_sec, ip) * dnumber_density_dt(abs_p[ip], abs_t[ip]);
230
231 for (Index iv = 0; iv < xsec_temp.nelem(); iv++) {
232 this_xsec(iv, ip) += n * xsec_temp[iv];
233 for (Index iq = 0; iq < jacobian_quantities.nelem();
234 iq++) {
235 const auto& deriv = jacobian_quantities[iq];
236
237 if (not deriv.propmattype()) continue;
238
239 if (is_frequency_parameter(deriv))
240 dabs_xsec_per_species_dx[i][iq](iv, ip) +=
241 n * (dxsec_temp_dF[iv] - xsec_temp[iv]) / df;
242 else if (deriv == Jacobian::Atm::Temperature)
243 dabs_xsec_per_species_dx[i][iq](iv, ip) +=
244 n * (dxsec_temp_dT[iv] - xsec_temp[iv]) / dt +
245 xsec_temp[iv] * dn_dT;
246 else if (species_match(deriv, this_species.cia_2nd_species))
247 dabs_xsec_per_species_dx[i][iq](iv, ip) +=
248 number_density(abs_p[ip], abs_t[ip]) * xsec_temp[iv];
249 }
250 }
251 }
252 }
253 }
254 }
255}
256
257/* Workspace method: Doxygen documentation will be auto-generated */
258void propmat_clearskyAddCIA( // WS Output:
259 PropagationMatrix& propmat_clearsky,
260 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
261 // WS Input:
262 const ArrayOfArrayOfSpeciesTag& abs_species,
263 const ArrayOfSpeciesTag& select_abs_species,
264 const ArrayOfRetrievalQuantity& jacobian_quantities,
265 const Vector& f_grid,
266 const Numeric& rtp_pressure,
267 const Numeric& rtp_temperature,
268 const Vector& rtp_vmr,
269 const ArrayOfCIARecord& abs_cia_data,
270 // WS Generic Input:
271 const Numeric& T_extrapolfac,
272 const Index& ignore_errors,
273 // Verbosity object:
274 const Verbosity& verbosity) {
276
277 // Size of problem
278 const Index nf = f_grid.nelem();
279 const Index nq = jacobian_quantities.nelem();
280 const Index ns = abs_species.nelem();
281
282 // Possible things that can go wrong in this code (excluding line parameters)
283 check_abs_species(abs_species);
284 ARTS_USER_ERROR_IF(rtp_vmr.nelem() not_eq abs_species.nelem(),
285 "*rtp_vmr* must match *abs_species*")
286 ARTS_USER_ERROR_IF(propmat_clearsky.NumberOfFrequencies() not_eq nf,
287 "*f_grid* must match *propmat_clearsky*")
289 not nq and (nq not_eq dpropmat_clearsky_dx.nelem()),
290 "*dpropmat_clearsky_dx* must match derived form of *jacobian_quantities*")
292 not nq and bad_propmat(dpropmat_clearsky_dx, f_grid),
293 "*dpropmat_clearsky_dx* must have frequency dim same as *f_grid*")
295 "Negative frequency (at least one value).")
297 "Must be sorted and increasing.")
299 "Negative VMR (at least one value).")
300 ARTS_USER_ERROR_IF(rtp_temperature <= 0, "Non-positive temperature")
301 ARTS_USER_ERROR_IF(rtp_pressure <= 0, "Non-positive pressure")
302
303 // Jacobian overhead START
304 const Numeric df = frequency_perturbation(jacobian_quantities);
305 const Numeric dt = temperature_perturbation(jacobian_quantities);
306
307 const bool do_jac = supports_CIA(
308 jacobian_quantities); // Throws runtime error if line parameters are wanted since we cannot know if the line is in the Continuum...
309 const bool do_wind_jac = do_wind_jacobian(jacobian_quantities);
310 const bool do_temp_jac = do_temperature_jacobian(jacobian_quantities);
311
312 Vector dfreq;
313 Vector dabs_t{rtp_temperature + dt};
314
315 if (do_wind_jac) {
316 dfreq.resize(f_grid.nelem());
317 for (Index iv = 0; iv < f_grid.nelem(); iv++) dfreq[iv] = f_grid[iv] + df;
318 }
319
320 Vector dxsec_temp_dF, dxsec_temp_dT;
321
322 if (do_jac) {
323 if (do_wind_jac) {
325 !std::isnormal(df), "df must be >0 and not NaN or Inf: ", df)
326 dxsec_temp_dF.resize(f_grid.nelem());
327 }
328 if (do_temp_jac) {
330 !std::isnormal(dt), "dt must be >0 and not NaN or Inf: ", dt)
331 dxsec_temp_dT.resize(f_grid.nelem());
332 }
333 }
334 // Jacobian overhead END
335
336 // Useful if there is no Jacobian to calculate
337 ArrayOfMatrix empty;
338
339 // Allocate a vector with dimension frequencies for constructing our
340 // cross-sections before adding them (more efficient to allocate this here
341 // outside of the loops)
342 Vector xsec_temp(f_grid.nelem());
343
344 // Loop over CIA data sets.
345 // Index ii loops through the outer array (different tag groups),
346 // index s through the inner array (different tags within each goup).
347 for (Index ispecies = 0; ispecies < ns; ispecies++) {
348 if (select_abs_species.nelem() and
349 select_abs_species not_eq abs_species[ispecies])
350 continue;
351
352 // Skip it if there are no species or there is Zeeman requested
353 if (not abs_species[ispecies].nelem() or abs_species[ispecies].Zeeman())
354 continue;
355
356 // Go through the tags in the current tag group to see if they
357 // are continuum tags:
358 for (Index s = 0; s < abs_species[ispecies].nelem(); ++s) {
359 const SpeciesTag& this_species = abs_species[ispecies][s];
360
361 // Check if this is a CIA tag
362 if (this_species.Type() != Species::TagType::Cia) continue;
363
364 // Get convenient references of this CIA data record and this
365 // absorption cross-section record:
366 Index this_cia_index = cia_get_index(
367 abs_cia_data, this_species.Spec(), this_species.cia_2nd_species);
368
369 ARTS_ASSERT(this_cia_index != -1);
370
371 const CIARecord& this_cia = abs_cia_data[this_cia_index];
372
373 if (out2.sufficient_priority()) {
374 // Some nice output to out2:
375 out2 << " CIA Species found: " + this_species.Name() + "\n";
376 }
377
378 // Find out index of VMR for the second CIA species.
379 // (The index for the first species is simply i.)
380 Index i_sec = find_first_species(abs_species, this_cia.Species(1));
381
382 // Catch the case that the VMR for the second species does not exist:
384 i_sec < 0,
385 "VMR profile for second species in CIA species pair does not exist.\n"
386 "Tag ",
387 this_species.Name(),
388 " needs a VMR profile of ",
389 this_cia.MoleculeName(1),
390 "!")
391
392 // We have to multiply with the number density of the second CIA species.
393 // We do not have to multiply with the first, since we still
394 // want to return a (unary) absorption cross-section, not an
395 // absorption coefficient.
396 const Numeric nd_sec =
397 number_density(rtp_pressure, rtp_temperature) * rtp_vmr[i_sec];
398 // Get the binary absorption cross sections from the CIA data:
399
400 try {
401 this_cia.Extract(xsec_temp,
402 f_grid,
403 rtp_temperature,
404 this_species.cia_dataset_index,
405 T_extrapolfac,
406 ignore_errors,
407 verbosity);
408 if (do_wind_jac) {
409 this_cia.Extract(dxsec_temp_dF,
410 dfreq,
411 rtp_temperature,
412 this_species.cia_dataset_index,
413 T_extrapolfac,
414 ignore_errors,
415 verbosity);
416 }
417 if (do_temp_jac) {
418 this_cia.Extract(dxsec_temp_dT,
419 f_grid,
420 rtp_temperature + dt,
421 this_species.cia_dataset_index,
422 T_extrapolfac,
423 ignore_errors,
424 verbosity);
425 }
426 } catch (const std::runtime_error& e) {
428 "Problem with CIA species ", this_species.Name(), ":\n", e.what())
429 }
430
431 if (!do_jac) {
432 xsec_temp *= nd_sec * number_density(rtp_pressure, rtp_temperature) *
433 rtp_vmr[ispecies];
434 propmat_clearsky.Kjj() += xsec_temp;
435 } else { // The Jacobian block
436 const Numeric nd = number_density(rtp_pressure, rtp_temperature);
437 const Numeric dnd_dt =
438 dnumber_density_dt(rtp_pressure, rtp_temperature);
439 const Numeric dnd_dt_sec =
440 dnumber_density_dt(rtp_pressure, rtp_temperature) * rtp_vmr[i_sec];
441 for (Index iv = 0; iv < f_grid.nelem(); iv++) {
442 propmat_clearsky.Kjj()[iv] +=
443 nd_sec * xsec_temp[iv] * nd * rtp_vmr[ispecies];
444 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
445 const auto& deriv = jacobian_quantities[iq];
446
447 if (not deriv.propmattype()) continue;
448
449 if (is_wind_parameter(deriv)) {
450 dpropmat_clearsky_dx[iq].Kjj()[iv] +=
451 nd_sec * (dxsec_temp_dF[iv] - xsec_temp[iv]) / df * nd *
452 rtp_vmr[ispecies];
453 } else if (deriv == Jacobian::Atm::Temperature) {
454 dpropmat_clearsky_dx[iq].Kjj()[iv] +=
455 ((nd_sec * (dxsec_temp_dT[iv] - xsec_temp[iv]) / dt +
456 xsec_temp[iv] * dnd_dt_sec) *
457 nd +
458 xsec_temp[iv] * nd_sec * dnd_dt) *
459 rtp_vmr[ispecies];
460 } else if (deriv == abs_species[ispecies]) {
461 dpropmat_clearsky_dx[iq].Kjj()[iv] +=
462 nd_sec * xsec_temp[iv] * nd;
463 } else if (species_match(deriv, this_species.Spec())) {
464 dpropmat_clearsky_dx[iq].Kjj()[iv] +=
465 nd * nd_sec * xsec_temp[iv] *
466 (this_species.cia_2nd_species == this_species.Spec() ? 2.0
467 : 1.0);
468 } else if (species_match(deriv, this_species.cia_2nd_species)) {
469 dpropmat_clearsky_dx[iq].Kjj()[iv] +=
470 nd * nd * xsec_temp[iv] * rtp_vmr[ispecies] *
471 (this_species.cia_2nd_species == this_species.Spec() ? 2.0
472 : 1.0);
473 }
474 }
475 }
476 }
477 }
478 }
479}
480
481/* Workspace method: Doxygen documentation will be auto-generated */
482void CIARecordReadFromFile( // WS GOutput:
483 CIARecord& cia_record,
484 // WS Generic Input:
485 const String& species_tag,
486 const String& filename,
487 // Verbosity object:
488 const Verbosity& verbosity) {
489 SpeciesTag species(species_tag);
490
491 ARTS_USER_ERROR_IF (species.Type() != Species::TagType::Cia,
492 "Invalid species tag ", species_tag, ".\n"
493 "This is not recognized as a CIA type.\n")
494
495 cia_record.SetSpecies(species.Spec(), species.cia_2nd_species);
496 cia_record.ReadFromCIA(filename, verbosity);
497}
498
499/* Workspace method: Doxygen documentation will be auto-generated */
500void abs_cia_dataAddCIARecord( // WS Output:
501 ArrayOfCIARecord& abs_cia_data,
502 // WS GInput:
503 const CIARecord& cia_record,
504 const Index& clobber,
505 // WS Input:
506 const Verbosity&) {
507 Index cia_index =
508 cia_get_index(abs_cia_data, cia_record.Species(0), cia_record.Species(1));
509 if (cia_index == -1)
510 abs_cia_data.push_back(cia_record);
511 else if (clobber)
512 abs_cia_data[cia_index] = cia_record;
513 else
514 abs_cia_data[cia_index].AppendDataset(cia_record);
515}
516
517/* Workspace method: Doxygen documentation will be auto-generated */
518void abs_cia_dataReadFromCIA( // WS Output:
519 ArrayOfCIARecord& abs_cia_data,
520 // WS Input:
521 const ArrayOfArrayOfSpeciesTag& abs_species,
522 const String& catalogpath,
523 const Verbosity& verbosity) {
524 ArrayOfString subfolders;
525 subfolders.push_back("Main-Folder/");
526 subfolders.push_back("Alternate-Folder/");
527
528 abs_cia_data.resize(0);
529
530 // Loop species tag groups to find CIA tags.
531 // Index sp loops through the tag groups, index iso through the tags within
532 // each group. Despite the name, iso does not denote the isotope!
533 for (Index sp = 0; sp < abs_species.nelem(); sp++) {
534 for (Index iso = 0; iso < abs_species[sp].nelem(); iso++) {
535 if (abs_species[sp][iso].Type() != Species::TagType::Cia) continue;
536
537 ArrayOfString cia_names;
538
539 Index cia_index = cia_get_index(abs_cia_data,
540 abs_species[sp][iso].Spec(),
541 abs_species[sp][iso].cia_2nd_species);
542
543 // If cia_index is not -1, we have already read this datafile earlier
544 if (cia_index != -1) continue;
545
546 cia_names.push_back(String(Species::toShortName(abs_species[sp][iso].Spec())) +
547 "-" +
548 String(Species::toShortName(abs_species[sp][iso].cia_2nd_species)));
549
550 cia_names.push_back(
551 String(Species::toShortName(abs_species[sp][iso].cia_2nd_species)) +
552 "-" +
553 String(Species::toShortName(abs_species[sp][iso].Spec())));
554
555 ArrayOfString checked_dirs;
556
557 bool found = false;
558 for (Index fname = 0; !found && fname < cia_names.nelem(); fname++) {
559 String cia_name = cia_names[fname];
560
561 for (Index dir = 0; !found && dir < subfolders.nelem(); dir++) {
562 ArrayOfString files;
563 checked_dirs.push_back(catalogpath + "/" + subfolders[dir] +
564 cia_name + "/");
565 try {
566 files = list_directory(*(checked_dirs.end() - 1));
567 } catch (const std::runtime_error& e) {
568 continue;
569 }
570
571 for (Index i = files.nelem() - 1; i >= 0; i--) {
572 if (files[i].find(cia_name) != 0 ||
573 files[i].rfind(".cia") != files[i].length() - 4) {
574 files.erase(files.begin() + i);
575 }
576 }
577 if (files.nelem()) {
578 CIARecord ciar;
579
580 found = true;
581 String catfile = *(checked_dirs.end() - 1) + files[0];
582
583 ciar.SetSpecies(abs_species[sp][iso].Spec(),
584 abs_species[sp][iso].cia_2nd_species);
585 ciar.ReadFromCIA(catfile, verbosity);
586
587 abs_cia_data.push_back(ciar);
588 }
589 }
590 }
591
592 ARTS_USER_ERROR_IF (!found,
593 "Error: No data file found for CIA species ", cia_names[0],
594 "\n"
595 "Looked in directories: ", checked_dirs)
596 }
597 }
598}
599
600/* Workspace method: Doxygen documentation will be auto-generated */
601void abs_cia_dataReadFromXML( // WS Output:
602 ArrayOfCIARecord& abs_cia_data,
603 // WS Input:
604 const ArrayOfArrayOfSpeciesTag& abs_species,
605 const String& filename,
606 const Verbosity& verbosity) {
607 xml_read_from_file(filename, abs_cia_data, verbosity);
608
609 // Check that all CIA tags from abs_species are present in the
610 // XML file
611
612 vector<String> missing_tags;
613
614 // Loop species tag groups to find CIA tags.
615 // Index sp loops through the tag groups, index iso through the tags within
616 // each group. Despite the name, iso does not denote the isotope!
617 for (Index sp = 0; sp < abs_species.nelem(); sp++) {
618 for (Index iso = 0; iso < abs_species[sp].nelem(); iso++) {
619 if (abs_species[sp][iso].Type() != Species::TagType::Cia) continue;
620
621 Index cia_index = cia_get_index(abs_cia_data,
622 abs_species[sp][iso].Spec(),
623 abs_species[sp][iso].cia_2nd_species);
624
625 // If cia_index is -1, this CIA tag was not present in the input file
626 if (cia_index == -1) {
627 missing_tags.push_back(
628 String(Species::toShortName(abs_species[sp][iso].Spec())) +
629 "-" +
630 String(Species::toShortName(abs_species[sp][iso].cia_2nd_species)));
631 }
632 }
633 }
634
635 if (missing_tags.size()) {
636 ostringstream os;
637 bool first = true;
638
639 os << "Error: The following CIA tag(s) are missing in input file: ";
640 for (size_t i = 0; i < missing_tags.size(); i++) {
641 if (!first)
642 os << ", ";
643 else
644 first = false;
645 os << missing_tags[i];
646 }
647 ARTS_USER_ERROR (os.str());
648 }
649}
650
651/* Workspace method: Doxygen documentation will be auto-generated */
652void CIAInfo( // Generic Input:
653 const String& catalogpath,
654 const ArrayOfString& cia_tags,
655 const Verbosity& verbosity) {
657
658 ArrayOfArrayOfSpeciesTag species_tags;
659
660 for (Index i = 0; i < cia_tags.nelem(); i++) {
661 ArrayOfSpeciesTag this_species_tag;
662
663 ArrayOfString species_names;
664
665 cia_tags[i].split(species_names, "-");
666
667 ARTS_USER_ERROR_IF (species_names.nelem() != 2,
668 "ERROR: Cannot parse CIA tag: ", cia_tags[i])
669
670 this_species_tag.push_back(
671 SpeciesTag(species_names[0] + "-CIA-" + species_names[1] + "-0"));
672
673 species_tags.push_back(this_species_tag);
674 }
675
676 ArrayOfCIARecord cia_data;
677
678 abs_cia_dataReadFromCIA(cia_data, species_tags, catalogpath, verbosity);
679
680 Print(cia_data, 1, verbosity);
681}
682
684 ArrayOfCIARecord& abs_cia_data,
685 const ArrayOfArrayOfSpeciesTag& abs_species,
686 const String& basename,
687 const Index& robust,
688 const Verbosity& verbosity) {
689 ArrayOfString names{};
690 for (auto& spec : abs_species) {
691 for (auto& tag : spec) {
692 if (tag.type == Species::TagType::Cia) {
693 names.emplace_back(
694 var_string(Species::toShortName(tag.Spec()),
695 "-CIA-",
696 Species::toShortName(tag.cia_2nd_species)));
697 }
698 }
699 }
700
701 names.erase(std::unique(names.begin(), names.end()), names.end());
702
703 const std::filesystem::path inpath{basename.c_str()};
704 const bool is_dir = std::filesystem::is_directory(inpath);
705
706 for (auto& name : names) {
707 auto fil{inpath};
708 if (is_dir) {
709 fil /= name + ".xml";
710 } else {
711 fil += "." + name + ".xml";
712 }
713
714 if (std::filesystem::is_regular_file(fil)) {
715 xml_read_from_file(fil.c_str(), abs_cia_data.emplace_back(), verbosity);
716 } else {
717 ARTS_USER_ERROR_IF(not robust,
718 "Cannot find ",
719 fil,
720 "\nIt is required to be there for ",
721 name,
722 " is found in a tag")
723 }
724 }
725}
Declarations required for the calculation of absorption coefficients.
The global header file for ARTS.
Constants of physical expressions as constexpr.
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:227
Header file for work with HITRAN collision induced absorption (CIA).
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
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:303
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:239
String MoleculeName(const Index i) const
Return each molecule name (as a string) that is associated with this CIARecord.
Definition: cia.cc:267
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The Matrix class.
Definition: matpackI.h:1261
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:899
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
std::string var_string(Args &&... args)
Definition: debug.h:36
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
ArrayOfString list_directory(const std::string_view dirname)
Return list of files in directory.
Definition: file.cc:505
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:1089
bool is_wind_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a wind parameter in propagation matrix calculations.
Definition: jacobian.cc:1002
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a frequency derivative.
Definition: jacobian.cc:1175
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:1119
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1150
bool do_wind_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a wind-based frequency derivative derivative.
Definition: jacobian.cc:1171
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Definition: jacobian.cc:1006
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the temperature perturbation if it exists.
Definition: jacobian.cc:1183
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
Definition: jacobian.cc:1191
#define ns
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:218
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:500
constexpr Numeric SPEED_OF_LIGHT
Definition: m_cia.cc:44
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:601
void propmat_clearskyAddCIA(PropagationMatrix &propmat_clearsky, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfSpeciesTag &select_abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &f_grid, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Vector &rtp_vmr, const ArrayOfCIARecord &abs_cia_data, const Numeric &T_extrapolfac, const Index &ignore_errors, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddCIA.
Definition: m_cia.cc:258
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)
Definition: m_cia.cc:47
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:518
void CIAInfo(const String &catalogpath, const ArrayOfString &cia_tags, const Verbosity &verbosity)
WORKSPACE METHOD: CIAInfo.
Definition: m_cia.cc:652
void CIARecordReadFromFile(CIARecord &cia_record, const String &species_tag, const String &filename, const Verbosity &verbosity)
WORKSPACE METHOD: CIARecordReadFromFile.
Definition: m_cia.cc:482
void abs_cia_dataReadSpeciesSplitCatalog(ArrayOfCIARecord &abs_cia_data, const ArrayOfArrayOfSpeciesTag &abs_species, const String &basename, const Index &robust, const Verbosity &verbosity)
WORKSPACE METHOD: abs_cia_dataReadSpeciesSplitCatalog.
Definition: m_cia.cc:683
void Print(Workspace &ws, const Agenda &x, const Index &level, const Verbosity &verbosity)
Definition: m_general.cc:76
constexpr bool any_negative(const MatpackType &var) noexcept
Checks for negative values.
Definition: math_funcs.h:147
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_OUT1
Definition: messages.h:204
#define CREATE_OUTS
Definition: messages.h:208
my_basic_string< char > String
The String type for ARTS.
Definition: mystring.h:216
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...
constexpr Numeric e
Elementary charge convenience name [C].
Temperature
Definition: jacobian.h:59
Implements Zeeman modeling.
Definition: zeemandata.cc:333
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
bool bad_propmat(const Array< T > &main, const Vector &f_grid, const Index sd=1) noexcept
Checks if a Propagation Matrix or something similar has good grids.
Index find_first_species(const ArrayOfArrayOfSpeciesTag &specs, Species::Species spec) noexcept
void check_abs_species(const ArrayOfArrayOfSpeciesTag &abs_species)
Species::Tag SpeciesTag
Definition: species_tags.h:101
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:151