ARTS 2.5.4 (git: 4c0d3b4d)
species_tags.cc
Go to the documentation of this file.
1#include "species_tags.h"
2
3#include <cfloat>
4#include <iterator>
5#include <string_view>
6
7#include "debug.h"
8
9namespace Species {
15constexpr bool is_modern_predefined(const IsotopeRecord& isot) {
16 constexpr std::array modern{
17 find_species_index(Species::Oxygen, "MPM2020"),
18 find_species_index(Species::Water, "ForeignContCKDMT350"),
19 find_species_index(Species::Water, "SelfContCKDMT350"),
20 };
21 const Index self = find_species_index(isot);
22
23 for (auto& x : modern)
24 if (x == self) return true;
25 return false;
26}
27
29 Index i = 0;
30 for (auto& x : Isotopologues) {
32 }
33 return i;
34}
35
51static_assert(legacy_predefined_count() == 58,
52 "2021-11-18: No more legacy models. "
53 "Please see src/predefined to add a modern model. "
54 "Never increment the count of legacy models. "
55 "If you remove any legacy model(s), decrement the count.");
56
57Tag::Tag(String def) : type(TagType::Plain) {
58 // Save input string for error messages:
59 String def_original = def;
60
61 // Name of species and isotopologue (aux variables):
62 String name, isoname;
63 // Aux index:
64 Index n;
65
66 // We cannot set a default value for the isotopologue, because the
67 // default should be `ALL' and the value for `ALL' depends on the
68 // species.
69
70 // Extract the species name:
71 n = def.find('-'); // find the '-'
72 if (n != def.npos) {
73 name = def.substr(0, n); // Extract before '-'
74 def.erase(0, n + 1); // Remove from def
75 } else {
76 // n==def.npos means that def does not contain a '-'. In that
77 // case we assume that it contains just a species name and
78 // nothing else
79 name = def;
80 def = "";
81 }
82
83 // Remove whitespace
84 name.trim();
85
86 // Obtain species index from species name.
87 // (This will also remove possible whitespace.)
88 const Species spec = fromShortName(name);
90 "Bad species name: ",
91 name,
92 " extracted from: ",
93 def_original)
94
95 // Check if species name contains the special tag for
96 // Faraday Rotation
97 if (spec == Species::free_electrons) {
98 constexpr Index ind =
99 find_species_index(IsotopeRecord(Species::free_electrons));
100 spec_ind = ind;
101 type = TagType::FreeElectrons;
102 return;
103 }
104
105 // Check if species name contains the special tag for
106 // Particles
107 if (spec == Species::particles) {
108 constexpr Index ind = find_species_index(IsotopeRecord(Species::particles));
109 spec_ind = ind;
110 type = TagType::Particles;
111 return;
112 }
113
114 // Set "all" species per default by leaving the joker in
115 spec_ind = find_species_index(IsotopeRecord(spec));
117 spec_ind < 0, "Bad species extracted: ", spec, " from ", def_original)
118 if (0 == def.nelem()) {
119 return;
120 }
121
122 // Extract the isotopologue name/Zeeman flag:
123 n = def.find('-'); // find the '-'
124 if (n != def.npos) {
125 isoname = def.substr(0, n); // Extract before '-'
126 def.erase(0, n + 1); // Remove from def
127
128 if ("Z" == isoname) {
129 type = TagType::Zeeman;
130 // Zeeman flag was present, now extract the isotopologue name:
131 n = def.find('-'); // find the '-'
132 if (n != def.npos) {
133 isoname = def.substr(0, n); // Extract before '-'
134 def.erase(0, n + 1); // Remove from def
135 } else {
136 // n==def.npos means that def does not contain a '-'. In that
137 // case we assume that it contains just the isotopologue name and
138 // nothing else.
139 isoname = def;
140 def = "";
141 }
142 }
143
144 if ("HXSEC" == isoname) {
145 type = TagType::HitranXsec;
146 // Hitran Xsec flag was present, now extract the isotopologue name:
147 n = def.find('-'); // find the '-'
148 if (n != def.npos) {
149 isoname = def.substr(0, n); // Extract before '-'
150 def.erase(0, n + 1); // Remove from def
151 } else {
152 // n==def.npos means that def does not contain a '-'. In that
153 // case we assume that it contains just the isotopologue name and
154 // nothing else.
155 isoname = def;
156 def = "";
157 }
158 }
159 } else {
160 // n==def.npos means that def does not contain a '-'. In that
161 // case we assume that it contains just the isotopologue name or
162 // Zeeman flag and nothing else.
163 isoname = def;
164 def = "";
165 if ("Z" == isoname) {
166 type = TagType::Zeeman;
167 // This means that there is nothing else to parse. Apparently
168 // the user wants all isotopologues and no frequency limits.
169 return;
170 }
171 if ("HXSEC" == isoname) {
172 type = TagType::HitranXsec;
173 // This means that there is nothing else to parse. Apparently
174 // the user wants all isotopologues and no frequency limits.
175 return;
176 }
177 }
178
179 if (Joker == isoname) {
180 // The user wants all isotopologues. This already set
181 } else if ("CIA" == isoname) // Check for "cia":
182 {
183 // The user wants this to use the CIA catalog:
184 type = TagType::Cia;
185
186 // We have to read in the second species, and the dataset index
187 n = def.find('-'); // find the '-'
188
190 n == def.npos,
191 "Invalid species tag ",
192 def_original,
193 ".\n"
194 "I am missing a minus sign (and a dataset index after that.)")
195
196 String otherspec = def.substr(0, n); // Extract before '-'
197 def.erase(0, n + 1); // Remove from def
198
199 cia_2nd_species = fromShortName(otherspec);
200 ARTS_USER_ERROR_IF(not good_enum(cia_2nd_species),
201 "Bad species name: ",
202 otherspec,
203 " extracted from: ",
204 def_original)
205
206 // Convert remaining def to dataset index.
207
208 // Check that everything remaining is just numbers.
209 for (Index i = 0; i < def.nelem(); ++i)
211 "Invalid species tag ",
212 def_original,
213 ".\n"
214 "The tag should end with a dataset index")
215
216 // Do the conversion from string to index:
217 std::istringstream is(def);
218 is >> cia_dataset_index;
219
220 def = "";
221 } else {
222 spec_ind = find_species_index(IsotopeRecord(spec, isoname));
223
224 // Check if we found a matching isotopologue:
225 ARTS_USER_ERROR_IF(spec_ind < 0,
226 "Isotopologue ",
227 isoname,
228 " is not a valid isotopologue or "
229 "absorption model for species ",
230 name,
231 ".\n"
232 "Valid options are:\n",
234
235 // Check if the found isotopologue represents a predefined model
236 // (continuum or full absorption model) and set the type accordingly:
237 if (is_predefined_model(Isotopologue()))
238 type = is_modern_predefined(Isotopologue())
239 ? TagType::PredefinedModern
240 : TagType::PredefinedLegacy;
241 }
242
243 if (0 == def.nelem()) {
244 // This means that there is nothing else to parse. Apparently
245 // the user wants no frequency limits. Frequency defaults are
246 // already set, so we can return directly.
247
248 return;
249 }
250
251 ARTS_USER_ERROR_IF(def[0] != Joker[0] && !isdigit(def[0]),
252 "Expected frequency limits, but got \"",
253 def,
254 "\"")
255
256 // Look for the two frequency limits:
257
258 // Extract first frequency
259 n = def.find('-'); // find the '-'
260 if (n != def.npos) {
261 // Frequency as a String:
262 String fname;
263 fname = def.substr(0, n); // Extract before '-'
264 def.erase(0, n + 1); // Remove from def
265
266 // Check for joker:
267 if (Joker == fname) {
268 // The default for lower_freq is already -1, meaning `ALL'.
269 // So there is nothing to do here.
270 } else {
271 ARTS_USER_ERROR_IF(!isdigit(fname[0]),
272 "Expected frequency limit, but got \"",
273 fname,
274 "\"")
275 // Convert to Numeric:
276 char* endptr;
277 lower_freq = strtod(fname.c_str(), &endptr);
278 ARTS_USER_ERROR_IF(endptr != fname.c_str() + fname.nelem(),
279 "Error parsing frequency limit \"",
280 fname,
281 "\"")
282 }
283 } else {
284 // n==def.npos means that def does not contain a '-'. In this
285 // case that is not allowed!
287 "You must either specify both frequency limits\n"
288 "(at least with jokers), or none.");
289 }
290
291 // Now there should only be the upper frequency left in def.
292 // Check for joker:
293 if (Joker== def) {
294 // The default for upper_freq is already -1, meaning `ALL'.
295 // So there is nothing to do here.
296 } else {
298 !isdigit(def[0]), "Expected frequency limit, but got \"", def, "\"")
299 // Convert to Numeric:
300 char* endptr;
301 upper_freq = strtod(def.c_str(), &endptr);
302 ARTS_USER_ERROR_IF(endptr != def.c_str() + def.nelem(),
303 "Error parsing frequency limit \"",
304 def,
305 "\"")
306 }
307}
308
310
322String Tag::Name() const {
323 std::ostringstream os;
324
325 // First the species name:
326 os << toShortName(Isotopologue().spec) << "-";
327
328 // Is this a CIA tag?
329 if (type == TagType::Cia) {
330 os << "CIA-" << toShortName(cia_2nd_species) << "-" << cia_dataset_index;
331
332 } else if (type == TagType::FreeElectrons || type == TagType::Particles) {
333 os << toShortName(Isotopologue().spec);
334 }
335 // Hitran Xsec flag.
336 else if (type == TagType::HitranXsec) {
337 os << "HXSEC";
338 } else {
339 // Zeeman flag.
340 if (type == TagType::Zeeman) os << "Z-";
341
342 // Now the isotopologue. Can be a single isotopologue or ALL.
343 os << Isotopologue().isotname << '-';
344
345 // Now the frequency limits, if there are any. For this we first
346 // need to determine the floating point precision.
347
348 // Determine the precision, depending on whether Numeric is double
349 // or float:
350 int precision;
351#ifdef USE_FLOAT
352 precision = FLT_DIG;
353#else
354#ifdef USE_DOUBLE
355 precision = DBL_DIG;
356#else
357#error Numeric must be double or float
358#endif
359#endif
360
361 if (0 > lower_freq) {
362 // lower_freq < 0 means no lower limit.
363 os << Joker << '-';
364 } else {
365 os << std::setprecision(precision);
366 os << lower_freq << "-";
367 }
368
369 if (0 > upper_freq) {
370 // upper_freq < 0 means no upper limit.
371 os << Joker;
372 } else {
373 os << std::setprecision(precision);
374 os << upper_freq;
375 }
376 }
377 return os.str();
378}
379} // namespace Species
380
382 // There can be a comma separated list of tag definitions, so we
383 // need to break the String apart at the commas.
384 ArrayOfString tag_def;
385
386 bool go_on = true;
387
388 while (go_on) {
389 // Index n = find_first( these_names, ',' );
390 Index n = these_names.find(',');
391 if (n == these_names.npos) // Value npos indicates not found.
392 {
393 // There are no more commas.
394 // cout << "these_names: (" << these_names << ")\n";
395 tag_def.push_back(these_names);
396 go_on = false;
397 } else {
398 tag_def.push_back(these_names.substr(0, n));
399 these_names.erase(0, n + 1);
400 }
401 }
402 // tag_def now holds the different tag Strings for this group.
403
404 // Set size to zero, in case the method has been called before.
405 resize(0);
406
407 for (Index s = 0; s < tag_def.nelem(); ++s) {
408 Species::Tag this_tag(tag_def[s]);
409
410 // Safety checks:
411 if (s > 0) {
412 // Tags inside a group must belong to the same species.
414 front().Isotopologue().spec != this_tag.Isotopologue().spec,
415 "Tags in a tag group must belong to the same species.");
416
417 // Zeeman tags and plain line by line tags must not be mixed. (Because
418 // there can be only one line list per tag group.)
419 ARTS_USER_ERROR_IF(((front().type == Species::TagType::Zeeman) &&
420 (this_tag.type == Species::TagType::Plain)) ||
421 ((front().type == Species::TagType::Plain) &&
422 (this_tag.type == Species::TagType::Zeeman)),
423 "Zeeman tags and plain line-by-line tags must "
424 "not be mixed in the same tag group.");
425 }
426
427 push_back(this_tag);
428 }
429}
430
432 Species::Species spec,
433 Index i) noexcept {
434 const Index n = specs.nelem();
435 for (; i < n; i++)
436 if (specs[i].Species() == spec) return i;
437 return -1;
438}
439
441 Species::Species spec) noexcept {
442 return find_next_species(specs, spec, 0);
443}
444
445std::pair<Index, Index> find_first_species_tag(
446 const ArrayOfArrayOfSpeciesTag& specs, const SpeciesTag& tag) noexcept {
447 for (Index i = 0; i < specs.nelem(); i++) {
448 if (auto ptr = std::find(specs[i].cbegin(), specs[i].cend(), tag);
449 ptr not_eq specs[i].cend())
450 return {i, std::distance(specs[i].cbegin(), ptr)};
451 }
452 return {-1, -1};
453}
454
455std::pair<Index, Index> find_first_isotologue(
456 const ArrayOfArrayOfSpeciesTag& specs,
457 const SpeciesIsotopeRecord& isot) noexcept {
458 for (Index i = 0; i < specs.nelem(); i++) {
459 if (auto ptr =
460 std::find_if(specs[i].cbegin(),
461 specs[i].cend(),
462 [&](auto& tag) { return tag.Isotopologue() == isot; });
463 ptr not_eq specs[i].cend())
464 return {i, std::distance(specs[i].cbegin(), ptr)};
465 }
466 return {-1, -1};
467}
468
470 Index num_free_electrons = 0;
471 for (Index i = 0; i < abs_species.nelem(); ++i) {
472 bool has_free_electrons = false;
473 bool has_particles = false;
474 bool has_hitran_xsec = false;
475 for (Index s = 0; s < abs_species[i].nelem(); ++s) {
476 if (abs_species[i][s].Type() == Species::TagType::FreeElectrons) {
477 num_free_electrons++;
478 has_free_electrons = true;
479 }
480
481 if (abs_species[i][s].Type() == Species::TagType::Particles) {
482 has_particles = true;
483 }
484
485 if (abs_species[i][s].Type() == Species::TagType::HitranXsec) {
486 has_hitran_xsec = true;
487 }
488 }
489
490 ARTS_USER_ERROR_IF(abs_species[i].nelem() > 1 && has_free_electrons,
491 "'free_electrons' must not be combined "
492 "with other tags in the same group.");
493 ARTS_USER_ERROR_IF(num_free_electrons > 1,
494 "'free_electrons' must not be defined "
495 "more than once.");
496
497 ARTS_USER_ERROR_IF(abs_species[i].nelem() > 1 && has_particles,
498 "'particles' must not be combined "
499 "with other tags in the same group.");
500
501 ARTS_USER_ERROR_IF(abs_species[i].nelem() > 1 && has_hitran_xsec,
502 "'hitran_xsec' must not be combined "
503 "with other tags in the same group.");
504 }
505}
506
507String ArrayOfSpeciesTag::Name() const {
508 String out = "";
509 bool first = true;
510 for (auto& x : *this) {
511 if (not first) out += ", ";
512 out += x.Name();
513 first = false;
514 }
515 return out;
516}
517
518std::set<Species::Species> lbl_species(
519 const ArrayOfArrayOfSpeciesTag& abs_species) noexcept {
520 std::set<Species::Species> unique_species;
521 for (auto& specs : abs_species) {
522 if (specs.RequireLines()) unique_species.insert(specs.front().Spec());
523 }
524 return unique_species;
525}
526
528 const Vector& rtp_vmr,
529 const Species spec) ARTS_NOEXCEPT {
530 ARTS_ASSERT(abs_species.nelem() == rtp_vmr.nelem())
531
532 auto pos =
533 std::find_if(abs_species.begin(),
534 abs_species.end(),
535 [spec](const ArrayOfSpeciesTag& tag_group) {
536 return tag_group.nelem() and tag_group.Species() == spec;
537 });
538 return pos == abs_species.end()
539 ? 0
540 : rtp_vmr[std::distance(abs_species.begin(), pos)];
541}
ArrayOfSpeciesTag() noexcept
Definition: species_tags.h:105
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
The Joker class.
Definition: matpackI.h:132
The Vector class.
Definition: matpackI.h:908
Index nelem() const
Number of elements.
Definition: mystring.h:253
void trim()
Trim leading and trailing whitespace.
Definition: mystring.h:232
static const Index npos
Define npos:
Definition: mystring.h:107
Helper macros for debugging.
#define ARTS_NOEXCEPT
Definition: debug.h:80
#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
constexpr bool good_enum(EnumType x) noexcept
Checks if the enum number is good.
Definition: enums.h:22
#define precision
Definition: logic.cc:43
void Error(const String &msg, const Verbosity &verbosity)
WORKSPACE METHOD: Error.
Definition: m_general.cc:397
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
Index nelem(const Lines &l)
Number of lines.
VectorView std(VectorView std, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the standard deviation of the ranged ys.
Definition: raw.cc:205
String isotopologues_names(Species spec)
constexpr bool is_predefined_model(const IsotopeRecord &ir) noexcept
constexpr std::array Isotopologues
A list of all ARTS isotopologues, note how the species enum class input HAS to be sorted.
Definition: isotopologues.h:57
Numeric first_vmr(const ArrayOfArrayOfSpeciesTag &abs_species, const Vector &rtp_vmr, const Species spec) ARTS_NOEXCEPT
constexpr bool is_modern_predefined(const IsotopeRecord &isot)
Checks if the isotname of an isotopologue match a modern model.
Definition: species_tags.cc:15
constexpr Index legacy_predefined_count()
Definition: species_tags.cc:28
constexpr Index find_species_index(const Species spec, const std::string_view isot) noexcept
constexpr std::string_view Joker
Definition: isotopologues.h:13
constexpr Rational end(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the largest M for a polarization type of this transition.
Definition: zeemandata.h:113
constexpr int isdigit(int ch) noexcept
Definition: nonstd.h:24
#define a
Index find_first_species(const ArrayOfArrayOfSpeciesTag &specs, Species::Species spec) noexcept
Index find_next_species(const ArrayOfArrayOfSpeciesTag &specs, Species::Species spec, Index i) noexcept
std::pair< Index, Index > find_first_isotologue(const ArrayOfArrayOfSpeciesTag &specs, const SpeciesIsotopeRecord &isot) noexcept
std::set< Species::Species > lbl_species(const ArrayOfArrayOfSpeciesTag &abs_species) noexcept
void check_abs_species(const ArrayOfArrayOfSpeciesTag &abs_species)
std::pair< Index, Index > find_first_species_tag(const ArrayOfArrayOfSpeciesTag &specs, const SpeciesTag &tag) noexcept
Species::Tag SpeciesTag
Definition: species_tags.h:101
Struct containing all information needed about one isotope.
Definition: isotopologues.h:16