31 ArrayOfMatrix& abs_xsec_per_species,
32 ArrayOfArrayOfMatrix& dabs_xsec_per_species_dx,
40 const Matrix& abs_vmrs,
43 const Numeric& T_extrapolfac,
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();
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")
80 dabs_t.resize(abs_t.nelem());
92 const Index n_p = abs_p.nelem();
93 const Index n_t = abs_t.nelem();
94 const Index nc_vmrs = abs_vmrs.ncols();
97 "The following variables must all have the same dimension:\n"
100 "abs_vmrs.ncols: ", nc_vmrs)
106 Vector xsec_temp(f_grid.nelem());
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());
118 for (Index ii = 0; ii < abs_species_active.
nelem(); ii++) {
119 const Index i = abs_species_active[ii];
121 for (Index s = 0; s < abs_species[i].
nelem(); s++) {
122 const SpeciesTag& this_species = abs_species[i][s];
125 if (this_species.Type() != Species::TagType::Cia)
continue;
130 abs_cia_data, this_species.Spec(), this_species.cia_2nd_species);
134 const CIARecord& this_cia = abs_cia_data[this_cia_index];
135 Matrix& this_xsec = abs_xsec_per_species[i];
137 if (out2.sufficient_priority()) {
139 out2 <<
" CIA Species found: " + this_species.Name() +
"\n";
145 "Wrong dimension of abs_xsec_per_species.nrows for species ", i,
147 "should match f_grid (", f_grid.nelem(),
") but is ",
148 this_xsec.nrows(),
".")
150 "Wrong dimension of abs_xsec_per_species.ncols for species ", i,
152 "should match abs_p (", abs_p.nelem(),
") but is ",
153 this_xsec.ncols(),
".")
161 "VMR profile for second species in CIA species pair does not exist.\n"
162 "Tag ", this_species.Name(),
" needs a VMR profile of ",
166 for (Index ip = 0; ip < abs_p.nelem(); ip++) {
176 this_cia.
Extract(dxsec_temp_dF,
183 this_cia.
Extract(dxsec_temp_dT,
189 }
catch (
const std::runtime_error& e) {
191 this_species.Name(),
":\n", e.what())
206 this_xsec(joker, ip) += xsec_temp;
208 const Numeric dn_dT =
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();
215 const auto& deriv = jacobian_quantities[iq];
217 if (not deriv.propmattype())
continue;
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;
227 dabs_xsec_per_species_dx[i][iq](iv, ip) +=
239 PropagationMatrix& propmat_clearsky,
240 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
245 const Vector& f_grid,
246 const Numeric& rtp_pressure,
247 const Numeric& rtp_temperature,
248 const Vector& rtp_vmr,
251 const Numeric& T_extrapolfac,
252 const Index& ignore_errors,
258 const Index nf = f_grid.nelem();
259 const Index nq = jacobian_quantities.
nelem();
260 const Index ns = abs_species.
nelem();
265 "*rtp_vmr* must match *abs_species*")
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).")
277 "Must be sorted and increasing.")
279 "Negative VMR (at least one value).")
288 jacobian_quantities);
293 Vector dabs_t{rtp_temperature + dt};
296 dfreq.resize(f_grid.nelem());
297 for (Index iv = 0; iv < f_grid.nelem(); iv++) dfreq[iv] = f_grid[iv] + df;
300 Vector dxsec_temp_dF, dxsec_temp_dT;
305 !std::isnormal(df),
"df must be >0 and not NaN or Inf: ", df)
306 dxsec_temp_dF.resize(f_grid.nelem());
310 !std::isnormal(dt),
"dt must be >0 and not NaN or Inf: ", dt)
311 dxsec_temp_dT.resize(f_grid.nelem());
322 Vector xsec_temp(f_grid.nelem());
327 for (Index ispecies = 0; ispecies < ns; ispecies++) {
328 if (select_abs_species.
nelem() and
329 select_abs_species not_eq abs_species[ispecies])
333 if (not abs_species[ispecies].nelem() or abs_species[ispecies].
Zeeman())
338 for (Index s = 0; s < abs_species[ispecies].
nelem(); ++s) {
339 const SpeciesTag& this_species = abs_species[ispecies][s];
342 if (this_species.Type() != Species::TagType::Cia)
continue;
347 abs_cia_data, this_species.Spec(), this_species.cia_2nd_species);
351 const CIARecord& this_cia = abs_cia_data[this_cia_index];
353 if (out2.sufficient_priority()) {
355 out2 <<
" CIA Species found: " + this_species.Name() +
"\n";
365 "VMR profile for second species in CIA species pair does not exist.\n"
368 " needs a VMR profile of ",
376 const Numeric nd_sec =
388 this_cia.
Extract(dxsec_temp_dF,
396 this_cia.
Extract(dxsec_temp_dT,
398 rtp_temperature + dt,
403 }
catch (
const std::runtime_error& e) {
405 "Problem with CIA species ", this_species.Name(),
":\n", e.what())
409 xsec_temp *= nd_sec *
number_density(rtp_pressure, rtp_temperature) *
411 propmat_clearsky.Kjj() += xsec_temp;
414 const Numeric dnd_dt =
416 const Numeric dnd_dt_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];
424 if (not deriv.propmattype())
continue;
427 dpropmat_clearsky_dx[iq].Kjj()[iv] +=
428 nd_sec * (dxsec_temp_dF[iv] - xsec_temp[iv]) / df * nd *
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) *
435 xsec_temp[iv] * nd_sec * dnd_dt) *
437 }
else if (deriv == abs_species[ispecies]) {
438 dpropmat_clearsky_dx[iq].Kjj()[iv] +=
439 nd_sec * xsec_temp[iv] * nd;
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
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
462 const String& species_tag,
469 "Invalid species tag ", species_tag,
".\n"
470 "This is not recognized as a CIA type.\n")
472 cia_record.
SetSpecies(species.Spec(), species.cia_2nd_species);
481 const Index& clobber,
487 abs_cia_data.push_back(cia_record);
489 abs_cia_data[cia_index] = cia_record;
491 abs_cia_data[cia_index].AppendDataset(cia_record);
499 const String& catalogpath,
502 subfolders.push_back(
"Main-Folder/");
503 subfolders.push_back(
"Alternate-Folder/");
505 abs_cia_data.resize(0);
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;
517 abs_species[sp][iso].Spec(),
518 abs_species[sp][iso].cia_2nd_species);
521 if (cia_index != -1)
continue;
523 cia_names.push_back(
String(Species::toShortName(abs_species[sp][iso].Spec())) +
525 String(Species::toShortName(abs_species[sp][iso].cia_2nd_species)));
528 String(Species::toShortName(abs_species[sp][iso].cia_2nd_species)) +
530 String(Species::toShortName(abs_species[sp][iso].Spec())));
535 for (Index fname = 0; !found && fname < cia_names.
nelem(); fname++) {
536 String cia_name = cia_names[fname];
538 for (Index dir = 0; !found && dir < subfolders.
nelem(); dir++) {
540 checked_dirs.push_back(catalogpath +
"/" + subfolders[dir] +
544 }
catch (
const std::runtime_error& e) {
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);
558 String catfile = *(checked_dirs.end() - 1) + files[0];
561 abs_species[sp][iso].cia_2nd_species);
564 abs_cia_data.push_back(ciar);
570 "Error: No data file found for CIA species ", cia_names[0],
572 "Looked in directories: ", checked_dirs)
589 vector<String> missing_tags;
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;
599 abs_species[sp][iso].Spec(),
600 abs_species[sp][iso].cia_2nd_species);
603 if (cia_index == -1) {
604 missing_tags.push_back(
605 String(Species::toShortName(abs_species[sp][iso].Spec())) +
607 String(Species::toShortName(abs_species[sp][iso].cia_2nd_species)));
612 if (missing_tags.size()) {
616 os <<
"Error: The following CIA tag(s) are missing in input file: ";
617 for (
size_t i = 0; i < missing_tags.size(); i++) {
622 os << missing_tags[i];
630 const String& catalogpath,
637 for (Index i = 0; i < cia_tags.
nelem(); i++) {
642 cia_tags[i].split(species_names,
"-");
645 "ERROR: Cannot parse CIA tag: ", cia_tags[i])
647 this_species_tag.push_back(
648 SpeciesTag(species_names[0] +
"-CIA-" + species_names[1] +
"-0"));
650 species_tags.push_back(this_species_tag);
657 Print(cia_data, 1, verbosity);
667 for (
auto& spec : abs_species) {
668 for (
auto& tag : spec) {
669 if (tag.type == Species::TagType::Cia) {
673 Species::toShortName(tag.cia_2nd_species)));
678 names.erase(std::unique(names.begin(), names.end()), names.end());
680 const std::filesystem::path inpath{basename.c_str()};
681 const bool is_dir = basename.back() ==
'/' or std::filesystem::is_directory(inpath);
683 for (
auto& name : names) {
686 fil /= name +
".xml";
688 fil +=
"." + name +
".xml";
694 "Cannot find any data for ", std::quoted(name),
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.
Header file for work with HITRAN collision induced absorption (CIA).
Index nelem() const ARTS_NOEXCEPT
CIA data for a single pair of molecules.
void SetSpecies(const Species::Species first, const Species::Species second)
Set CIA species.
Species::Species Species(const Index i) const
Return CIA species index.
void ReadFromCIA(const String &filename, const Verbosity &verbosity)
Read CIA catalog file.
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.
String MoleculeName(const Index i) const
Return each molecule name (as a string) that is associated with this CIARecord.
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
std::string var_string(Args &&... args)
#define ARTS_USER_ERROR_IF(condition,...)
ArrayOfString list_directory(const std::string_view dirname)
Return list of files in directory.
This file contains basic functions to handle ASCII files.
bool supports_CIA(const ArrayOfRetrievalQuantity &js)
Returns if the array supports CIA derivatives.
bool is_wind_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a wind parameter in propagation matrix calculations.
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a frequency derivative.
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.
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
bool do_wind_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a wind-based frequency derivative derivative.
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the temperature perturbation if it exists.
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
void abs_cia_dataAddCIARecord(ArrayOfCIARecord &abs_cia_data, const CIARecord &cia_record, const Index &clobber, const Verbosity &)
WORKSPACE METHOD: abs_cia_dataAddCIARecord.
constexpr Numeric SPEED_OF_LIGHT
void abs_cia_dataReadFromXML(ArrayOfCIARecord &abs_cia_data, const ArrayOfArrayOfSpeciesTag &abs_species, const String &filename, const Verbosity &verbosity)
WORKSPACE METHOD: abs_cia_dataReadFromXML.
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.
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)
void abs_cia_dataReadFromCIA(ArrayOfCIARecord &abs_cia_data, const ArrayOfArrayOfSpeciesTag &abs_species, const String &catalogpath, const Verbosity &verbosity)
WORKSPACE METHOD: abs_cia_dataReadFromCIA.
void CIAInfo(const String &catalogpath, const ArrayOfString &cia_tags, const Verbosity &verbosity)
WORKSPACE METHOD: CIAInfo.
void CIARecordReadFromFile(CIARecord &cia_record, const String &species_tag, const String &filename, const Verbosity &verbosity)
WORKSPACE METHOD: CIARecordReadFromFile.
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.
void Print(Workspace &ws, const Agenda &x, const Index &level, const Verbosity &verbosity)
constexpr bool any_negative(const MatpackType &var) noexcept
Checks for negative values.
Declarations having to do with the four output streams.
my_basic_string< char > String
The String type for ARTS.
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
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.