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