ARTS 2.5.4 (git: 31ce4f0e)
lineshapemodel.cc
Go to the documentation of this file.
1/* Copyright (C) 2018
2 Richard Larsson <ric.larsson@gmail.com>
3
4
5 This program is free software; you can redistribute it and/or modify it
6 under the terms of the GNU General Public License as published by the
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
20 USA. */
21
35#include "arts_conversions.h"
36#include "lineshapemodel.h"
37#include "matpackI.h"
38
40 const String& coeff) {
41 // Test viability of model variables
42 const auto var_type = LineShape::toVariableOrThrow(var);
43
44 // Test viability of model coefficients
45 const auto coeff_type = Options::toLineShapeCoeffOrThrow(coeff);
46
47// Define a repetitive pattern. Update if/when there are more coefficients in the future
48#define ReturnJacPropMatType(ID) \
49 case LineShape::Variable::ID: \
50 switch (coeff_type) { \
51 case Options::LineShapeCoeff::X0: \
52 return Jacobian::Line::Shape##ID##X0; \
53 case Options::LineShapeCoeff::X1: \
54 return Jacobian::Line::Shape##ID##X1; \
55 case Options::LineShapeCoeff::X2: \
56 return Jacobian::Line::Shape##ID##X2; \
57 case Options::LineShapeCoeff::X3: \
58 return Jacobian::Line::Shape##ID##X3; \
59 case Options::LineShapeCoeff::FINAL: \
60 return Jacobian::Line::FINAL; \
61 } break
62
63 switch (var_type) {
73 case LineShape::Variable::FINAL: {
74 /* Leave last */ }
75 }
76
77#undef ReturnJacPropMatType
78
79 return Jacobian::Line::FINAL;
80}
81
82std::istream& LineShape::from_artscat4(std::istream& is,
83 Type& mtype,
84 bool& self,
85 bool& bath,
86 Model& m,
87 ArrayOfSpecies& species,
88 const QuantumIdentifier& qid) {
89 // Set or reset variables
90 mtype = Type::VP;
91 self = true;
92 bath = false;
93 m.mdata = std::vector<SingleSpeciesModel>(7);
94 species = ArrayOfSpecies(7);
95
96 // Set species
97 species[1] = Species::fromShortName("N2");
98 species[2] = Species::fromShortName("O2");
99 species[3] = Species::fromShortName("H2O");
100 species[4] = Species::fromShortName("CO2");
101 species[5] = Species::fromShortName("H2");
102 species[6] = Species::fromShortName("He");
103
104 // Temperature types
105 for (auto& v : m.mdata) {
106 v.G0().type = TemperatureModel::T1;
107 v.D0().type = TemperatureModel::T5;
108 }
109
110 // G0 main coefficient
111 for (auto& v : m.mdata) is >> double_imanip() >> v.G0().X0;
112
113 // G0 exponent is same as D0 exponent
114 for (auto& v : m.mdata) {
115 is >> double_imanip() >> v.G0().X1;
116 v.D0().X1 = v.G0().X1;
117 }
118
119 // D0 coefficient
120 m.mdata.front().D0().X0 = 0;
121 for (int k = 1; k < 7; k++)
122 is >> double_imanip() >> m.mdata[k].D0().X0;
123
124 // Special case when self is part of this list, it needs to be removed
125 for (int k = 1; k < 7; k++) {
126 if (qid.Species() == species[k]) {
127 ARTS_USER_ERROR_IF(m.mdata.front().G0().X0 not_eq m.mdata[k].G0().X0 or
128 m.mdata.front().G0().X1 not_eq m.mdata[k].G0().X1 or
129 m.mdata.front().D0().X1 not_eq m.mdata[k].D0().X1,
130 "Species is ", qid.Isotopologue(), " and this is a broadening species in ARTSCAT-4.\n"
131 "Despite this, values representing self and ", qid.Isotopologue(), " does not match "
132 "in input string\n")
133 m.mdata.front().D0().X0 = m.mdata[k].D0().X0;
134 m.Remove(k, species);
135 break;
136 }
137 }
138
139 return is;
140}
141
142std::istream& LineShape::from_linefunctiondata(std::istream& data,
143 Type& mtype,
144 bool& self,
145 bool& bath,
146 Model& m,
147 ArrayOfSpecies& species) {
148 self = bath = false;
149 Index specs;
150 String s;
151
152 // The first tag should give the line shape scheme
153 data >> mtype;
154 check_enum_error(mtype, "Bad Data");
155
156 // Order of elements for line shape
157 const auto shapeparams =
159
160 // The second tag should give the line mixing scheme
161 data >> s;
162
163 // Order of elements for line mixing
164 const auto mixingparams =
166
167 // The third tag should contain the number of species
168 data >> specs;
169 species.resize(specs);
170 m.mdata.resize(specs);
171
172 ARTS_USER_ERROR_IF (not specs and mtype not_eq Type::DP,
173 "Need at least one species for non-Doppler line shapes");
174
175 // For all species, we need to set the methods to compute them
176 for (Index i = 0; i < specs; i++) {
177 // This should be a species tag or one of the specials, SELF or BATH
178 data >> s;
179 if (s == self_broadening) {
180 // If the species is self, then we need to flag this
181 self = true;
182 ARTS_USER_ERROR_IF (i not_eq 0,
183 "Self broadening must be first, it is not\n");
184 } else if (s == bath_broadening) {
185 // If the species is air, then we need to flag this
186 bath = true;
187 species[i] = Species::Species::Bath;
188 ARTS_USER_ERROR_IF (i not_eq specs - 1,
189 "Air/bath broadening must be last, it is not\n");
190 } else {
191 // Otherwise, we hope we find a species
192 species[i] = Species::fromShortName(s);
193 ARTS_USER_ERROR_IF (not good_enum(species[i]),
194 "Encountered ", s, " in a position where a species should have been "
195 "defined.\nPlease check your pressure broadening data structure and ensure "
196 "that it follows the correct conventions.\n")
197 }
198
199 // For all parameters
200 for (auto& params : {shapeparams, mixingparams}) {
201 for (auto& param : params) {
202 data >> s; // Should contain a temperature tag
203
204 const auto type = toTemperatureModel(s);
205 const Index ntemp =
207
208 m.mdata[i].Data()[Index(param)].type = type;
209 if (ntemp <= ModelParameters::N) {
210 switch (ntemp) {
211 case 1:
212 data >> double_imanip() >> m.mdata[i].Data()[Index(param)].X0;
213 break;
214 case 2:
215 data >> double_imanip() >> m.mdata[i].Data()[Index(param)].X0 >>
216 m.mdata[i].Data()[Index(param)].X1;
217 break;
218 case 3:
219 data >> double_imanip() >> m.mdata[i].Data()[Index(param)].X0 >>
220 m.mdata[i].Data()[Index(param)].X1 >>
221 m.mdata[i].Data()[Index(param)].X2;
222 break;
223 case 0:
224 break;
225 default:
227 "Unknown number of input parameters in Legacy mode.");
228 }
229 } else { // Has to be the only allowed interpolation case
230 ARTS_USER_ERROR_IF (ntemp > 12,
231 "Too many input parameters in interpolation results Legacy mode.");
233 data >> double_imanip() >> temp; // should be 200
234 data >> double_imanip() >> temp; // should be 250
235 data >> double_imanip() >> temp; // should be 296
236 data >> double_imanip() >> temp; // should be 340
237 data >> double_imanip() >> m.mdata[i].Y().X0 >> m.mdata[i].Y().X1 >> m.mdata[i].Y().X2 >> m.mdata[i].Y().X3
238 >> m.mdata[i].G().X0 >> m.mdata[i].G().X1 >> m.mdata[i].G().X2 >> m.mdata[i].G().X3;
239 }
240 }
241 }
242 }
243
244 return data;
245}
246
248 std::istream& data,
249 LineShape::Type& mtype,
250 bool& self,
251 bool& bath,
252 Model& m,
253 ArrayOfSpecies& species,
254 const QuantumIdentifier& qid) {
255 String s;
256 data >> s;
257
260 const auto self_in_list = LegacyPressureBroadeningData::self_listed(qid, type);
261
262 Vector x(n);
263 for (auto& num : x) data >> double_imanip() >> num;
264
266 self,
267 bath,
268 m,
269 species,
270 x,
271 type,
272 self_in_list,
273 qid.Species());
274
275 return data;
276}
277
278std::istream& LineShape::from_linemixingdata(std::istream& data,
279 LineShape::Model& lsc) {
280 String s;
281 data >> s;
282
283 const auto type = LegacyLineMixingData::string2typelm(s);
284 const auto n = LegacyLineMixingData::typelm2nelem(type);
285
286 Vector x(n);
287 for (auto& num : x) data >> double_imanip() >> num;
288
290
291 return data;
292}
293
294#pragma GCC diagnostic push
295#pragma GCC diagnostic ignored "-Wreturn-type"
297 LineShape::Type& mtype,
298 bool& self,
299 bool& bath,
300 Model& m,
301 ArrayOfSpecies& species,
302 Vector x,
304 bool self_in_list,
305 Species::Species self_spec) {
306 switch (type) {
307 case TypePB::PB_NONE:
308 mtype = LineShape::Type::DP;
309 self = bath = false;
310 m = Model();
311 species.resize(0);
312 return;
313 case TypePB::PB_AIR_BROADENING:
314 mtype = LineShape::Type::VP;
315 self = bath = true;
316 m = Model(x[0], x[1], x[2], x[3], x[4]);
317 species.resize(2);
318 species[0] = self_spec;
319 species[1] = Species::Species::Bath;
320 return;
321 case TypePB::PB_AIR_AND_WATER_BROADENING:
322 if (self_in_list) {
323 mtype = LineShape::Type::VP;
324 self = false;
325 bath = true;
326 m.Data().resize(2);
327 m.Data()[0].G0() = {TemperatureModel::T1, x[0], x[1], 0, 0};
328 m.Data()[0].D0() = {TemperatureModel::T5, x[2], x[1], 0, 0};
329 m.Data()[1].G0() = {TemperatureModel::T1, x[3], x[4], 0, 0};
330 m.Data()[1].D0() = {TemperatureModel::T5, x[5], x[4], 0, 0};
331 species.resize(2);
332 species[0] = Species::fromShortName("H2O");
333 species[1] = Species::Species::Bath;
334 return;
335 } else {
336 mtype = LineShape::Type::VP;
337 self = bath = true;
338 m.Data().resize(2);
339 m.Data()[0].G0() = {TemperatureModel::T1, x[0], x[1], 0, 0};
340 m.Data()[0].D0() = {TemperatureModel::T5, x[2], x[1], 0, 0};
341 m.Data()[2].G0() = {TemperatureModel::T1, x[3], x[4], 0, 0};
342 m.Data()[2].D0() = {TemperatureModel::T5, x[5], x[4], 0, 0};
343 m.Data()[1].G0() = {TemperatureModel::T1, x[6], x[7], 0, 0};
344 m.Data()[1].D0() = {TemperatureModel::T5, x[8], x[7], 0, 0};
345 species.resize(3);
346 species[0] = self_spec;
347 species[1] = Species::fromShortName("H2O");
348 species[2] = Species::Species::Bath;
349 return;
350 }
351 case TypePB::PB_PLANETARY_BROADENING:
352 if (self_in_list) {
353 mtype = LineShape::Type::VP;
354 self = bath = false;
355 m.Data().resize(6);
356 m.Data()[0].G0() = {TemperatureModel::T1, x[1], x[8], 0, 0};
357 m.Data()[0].D0() = {TemperatureModel::T5, x[14], x[8], 0, 0};
358 m.Data()[1].G0() = {TemperatureModel::T1, x[2], x[9], 0, 0};
359 m.Data()[1].D0() = {TemperatureModel::T5, x[15], x[9], 0, 0};
360 m.Data()[2].G0() = {TemperatureModel::T1, x[3], x[10], 0, 0};
361 m.Data()[2].D0() = {TemperatureModel::T5, x[16], x[10], 0, 0};
362 m.Data()[3].G0() = {TemperatureModel::T1, x[4], x[11], 0, 0};
363 m.Data()[3].D0() = {TemperatureModel::T5, x[17], x[11], 0, 0};
364 m.Data()[4].G0() = {TemperatureModel::T1, x[5], x[12], 0, 0};
365 m.Data()[4].D0() = {TemperatureModel::T5, x[18], x[12], 0, 0};
366 m.Data()[5].G0() = {TemperatureModel::T1, x[6], x[13], 0, 0};
367 m.Data()[5].D0() = {TemperatureModel::T5, x[19], x[13], 0, 0};
368 species = {Species::fromShortName("N2"),
369 Species::fromShortName("O2"),
370 Species::fromShortName("H2O"),
371 Species::fromShortName("CO2"),
372 Species::fromShortName("H2"),
373 Species::fromShortName("He")};
374 return;
375 } else {
376 mtype = LineShape::Type::VP;
377 self = true;
378 bath = false;
379 m.Data().resize(7);
380 m.Data()[0].G0() = {TemperatureModel::T1, x[0], x[7], 0, 0};
381 // ssm[0].D0() = ...
382 m.Data()[1].G0() = {TemperatureModel::T1, x[1], x[8], 0, 0};
383 m.Data()[1].D0() = {TemperatureModel::T5, x[14], x[8], 0, 0};
384 m.Data()[2].G0() = {TemperatureModel::T1, x[2], x[9], 0, 0};
385 m.Data()[2].D0() = {TemperatureModel::T5, x[15], x[9], 0, 0};
386 m.Data()[3].G0() = {TemperatureModel::T1, x[3], x[10], 0, 0};
387 m.Data()[3].D0() = {TemperatureModel::T5, x[16], x[10], 0, 0};
388 m.Data()[4].G0() = {TemperatureModel::T1, x[4], x[11], 0, 0};
389 m.Data()[4].D0() = {TemperatureModel::T5, x[17], x[11], 0, 0};
390 m.Data()[5].G0() = {TemperatureModel::T1, x[5], x[12], 0, 0};
391 m.Data()[5].D0() = {TemperatureModel::T5, x[18], x[12], 0, 0};
392 m.Data()[6].G0() = {TemperatureModel::T1, x[6], x[13], 0, 0};
393 m.Data()[6].D0() = {TemperatureModel::T5, x[19], x[13], 0, 0};
394 species.resize(7);
395 species[0] = self_spec;
396 species[1] = Species::fromShortName("N2");
397 species[2] = Species::fromShortName("O2");
398 species[3] = Species::fromShortName("H2O");
399 species[4] = Species::fromShortName("CO2");
400 species[5] = Species::fromShortName("H2");
401 species[6] = Species::fromShortName("He");
402 return;
403 }
404 }
405}
406#pragma GCC diagnostic pop
407
410 Model y(1);
411 switch (type) {
412 case TypeLM::LM_NONE:
413 break;
414 case TypeLM::LM_LBLRTM:
415 y.Data().front().Y().type = LineShape::TemperatureModel::LM_AER;
416 y.Data().front().G().type = LineShape::TemperatureModel::LM_AER;
417 y.Data().front().Y().X0 = x[4];
418 y.Data().front().Y().X1 = x[5];
419 y.Data().front().Y().X2 = x[6];
420 y.Data().front().Y().X3 = x[7];
421 y.Data().front().G().X0 = x[8];
422 y.Data().front().G().X1 = x[9];
423 y.Data().front().G().X2 = x[10];
424 y.Data().front().G().X3 = x[11];
425 break;
426 case TypeLM::LM_LBLRTM_O2NonResonant:
427 y.Data().front().G().type = LineShape::TemperatureModel::T0;
428 y.Data().front().G().X0 = x[0];
429 break;
430 case TypeLM::LM_2NDORDER:
431 y.Data().front().Y().type = LineShape::TemperatureModel::T4;
432 y.Data().front().Y().X0 = x[0];
433 y.Data().front().Y().X1 = x[1];
434 y.Data().front().Y().X2 = x[7];
435 y.Data().front().G().type = LineShape::TemperatureModel::T4;
436 y.Data().front().G().X0 = x[2];
437 y.Data().front().G().X1 = x[3];
438 y.Data().front().G().X2 = x[8];
439 y.Data().front().DV().type = LineShape::TemperatureModel::T4;
440 y.Data().front().DV().X0 = x[4];
441 y.Data().front().DV().X1 = x[5];
442 y.Data().front().DV().X2 = x[9];
443 break;
444 case TypeLM::LM_1STORDER:
445 y.Data().front().Y().type = LineShape::TemperatureModel::T1;
446 y.Data().front().Y().X0 = x[1];
447 y.Data().front().Y().X1 = x[2];
448 break;
449 case TypeLM::LM_BYBAND:
450 break;
451 }
452 return y;
453}
454
455
456Vector LineShape::vmrs(const ConstVectorView& atmospheric_vmrs,
457 const ArrayOfArrayOfSpeciesTag& atmospheric_species,
458 const ArrayOfSpecies& lineshape_species) ARTS_NOEXCEPT {
459 ARTS_ASSERT (atmospheric_species.nelem() == atmospheric_vmrs.nelem(), "Bad atmospheric inputs");
460
461 const Index n = lineshape_species.nelem();
462
463 // Initialize list of VMRS to 0
464 Vector line_vmrs(n, 0);
465
466 // We need to know if bath is an actual species
467 const bool bath = n and lineshape_species.back() == Species::Species::Bath;
468
469 // Loop species
470 for (Index i = 0; i < n - bath; i++) {
471 const Species::Species target = lineshape_species[i];
472
473 // Find species in list or do nothing at all
474 Index this_species_index = -1;
475 for (Index j = 0; j < atmospheric_species.nelem(); j++) {
476 if (atmospheric_species[j].Species() == target) {
477 this_species_index = j;
478 }
479 }
480
481 // Set to non-zero in-case species exists
482 if (this_species_index not_eq -1) {
483 line_vmrs[i] = atmospheric_vmrs[this_species_index];
484 }
485 }
486
487 // Renormalize, if bath-species exist this is automatic.
488 if (bath) {
489 line_vmrs[n - 1] = 1.0 - line_vmrs.sum();
490 } else if(line_vmrs.sum() == 0) { // Special case
491 } else {
492 line_vmrs /= line_vmrs.sum();
493 }
494
495 return line_vmrs;
496}
497
498Vector LineShape::mass(const ConstVectorView& atmospheric_vmrs,
499 const ArrayOfArrayOfSpeciesTag& atmospheric_species,
500 const ArrayOfSpecies& lineshape_species,
502 ARTS_ASSERT (atmospheric_species.nelem() == atmospheric_vmrs.nelem(),
503 "Bad atmospheric inputs");
504
505 const Index n = lineshape_species.nelem();
506
507 // Initialize list of VMRS to 0
508 Vector line_vmrs(n, 0);
509 Vector line_mass(n, 0);
510
511 // We need to know if bath is an actual species
512 const bool bath = n and lineshape_species.back() == Species::Species::Bath;
513
514 // Loop species ignoring self and bath
515 for (Index i = 0; i < n - bath; i++) {
516 // Select target in-case this is self-broadening
517 const Species::Species target = lineshape_species[i];
518
519 // Find species in list or do nothing at all
520 Index this_species_index = -1;
521 for (Index j = 0; j < atmospheric_species.nelem(); j++) {
522 if (atmospheric_species[j].Species() == target) {
523 this_species_index = j;
524 }
525 }
526
527 // Set to non-zero in-case species exists
528 if (this_species_index not_eq -1) {
529 line_vmrs[i] = atmospheric_vmrs[this_species_index];
530 line_mass[i] = Species::mean_mass(target, ir);
531 }
532 }
533
534 // Renormalize, if bath-species exist this is automatic.
535 if(line_vmrs.sum() == 0) { // Special case
536 } else if (bath) {
537 line_mass[n - 1] = (line_vmrs * line_mass) / line_vmrs.sum();
538 }
539
540 return line_mass;
541}
542
543namespace LineShape {
544std::ostream& operator<<(std::ostream& os, const Model& m) {
545 for (auto& data : m.Data()) os << data;
546 return os;
547}
548
549std::istream& operator>>(std::istream& is, Model& m) {
550 for (auto& data : m.Data()) is >> data;
551 return is;
552}
553
555 String out = "";
556 const auto& vars = enumtyps::VariableTypes;
557
558 for (auto& var : vars) {
559 if (std::any_of(m.Data().cbegin(), m.Data().cend(), [var](auto& x) {
560 return x.Get(var).type not_eq TemperatureModel::None;
561 })) {
562 out += String(toString(var)) + ' ';
563 for (auto& ssm : m.Data())
564 out += String(toString(ssm.Get(var).type)) + ' ';
565 }
566 }
567
568 if (out.size()) out.pop_back();
569
570 return out;
571}
572
574 if (s.nelem() == 0) return Model();
575
576 const auto& names = enumstrs::VariableNames;
577
578 std::istringstream str(s);
579 String part;
580 Variable var = Variable::ETA;
581 TemperatureModel tm = TemperatureModel::None;
582 Index i = -100000;
583
584 std::vector<SingleSpeciesModel> ssms(0);
585 while (not str.eof()) {
586 str >> part;
587 if (std::any_of(names.cbegin(), names.cend(), [part](auto x) {
588 return part == x;
589 })) {
590 i = -1;
591 var = toVariable(part);
592 } else {
593 i++;
594 tm = toTemperatureModel(part);
595 }
596
597 if (i < 0) continue;
598 if (i < Index(ssms.size()))
599 goto add_var;
600 else {
601 ssms.emplace_back();
602 add_var:
603 auto mp = ssms[i].Get(var);
604 mp.type = tm;
605 ssms[i].Set(var, mp);
606 }
607 }
608
609 return Model(ssms);
610}
611
613 std::ostringstream os;
614 switch (mp.type) {
615 case TemperatureModel::None:
616 os << 0;
617 break;
618 case TemperatureModel::T0:
619 os << mp.X0;
620 break;
621 case TemperatureModel::T1:
622 os << mp.X0 << " * (" << T0 << "/T)^" << mp.X1;
623 break;
624 case TemperatureModel::T2:
625 os << mp.X0 << " * (" << T0 << "/T)^" << mp.X1 << " / (1 + " << mp.X2
626 << " * log(T/" << T0 << "))";
627 break;
628 case TemperatureModel::T3:
629 os << mp.X0 << " + " << mp.X1 << " * (" << T0 << " - T)";
630 break;
631 case TemperatureModel::T4:
632 os << "(" << mp.X0 << " + " << mp.X1 << " * (" << T0 << "/T - 1)) * ("
633 << T0 << "/T)^" << mp.X2;
634 break;
635 case TemperatureModel::T5:
636 os << mp.X0 << " * (" << T0 << "/T)^(0.25 + 1.5 * " << mp.X1 << ")";
637 break;
638 case TemperatureModel::LM_AER:
639 os << '('
640 << "Linear interpolation to y(x) from x-ref = [200, 250, 296, 340] and y-ref = ["
641 << mp.X0 << ", " << mp.X1 << ", " << mp.X2 << ", " << mp.X3 << ']'
642 << ')';
643 break;
644 case TemperatureModel::DPL:
645 os << '(' << mp.X0 << " * (" << T0 << "/T)^" << mp.X1 << " + " << mp.X2
646 << " * (" << T0 << "/T)^" << mp.X3 << ')';
647 break;
648 case TemperatureModel::POLY:
649 os << '(' << mp.X0 << " + " << mp.X1 << " * T + " << mp.X2
650 << " * T * T + " << mp.X3 << " * T * T * T)";
651 case TemperatureModel::FINAL:
652 break;
653 }
654
655 return os.str();
656}
657
659 const bool self,
660 const ArrayOfSpecies& sts,
661 const Numeric T0) {
662 const auto& vars = enumtyps::VariableTypes;
663 ArrayOfString as(0);
664
665 for (Index i = 0; i < Index(Variable::FINAL); i++) {
666 Variable var = vars[i];
667
668 if (std::any_of(m.Data().cbegin(), m.Data().cend(), [var](auto& x) {
669 return x.Get(var).type not_eq TemperatureModel::None;
670 })) {
671 std::ostringstream os;
672 os << var << " ~ ";
673 for (Index j = 0; j < sts.nelem(); j++) {
674 if (j == 0 and self)
675 os << "VMR(" << self_broadening << ") * "
676 << modelparameters2metadata(m.Data().front().Get(var), T0);
677 else
678 os << "VMR(" << toShortName(sts[j]) << ") * "
679 << modelparameters2metadata(m.Data()[j].Get(var), T0);
680
681 if (sts[j] not_eq sts.back()) os << " + ";
682 }
683 as.push_back(os.str());
684 }
685 }
686
687 return as;
688}
689
690#pragma GCC diagnostic push
691#pragma GCC diagnostic ignored "-Wreturn-type"
693 if (type == "X0")
694 return mp.X0;
695 if (type == "X1")
696 return mp.X1;
697 if (type == "X2")
698 return mp.X2;
699 if (type == "X3")
700 return mp.X3;
702 "Type: ", type, ", is not accepted. "
703 "See documentation for accepted types\n")
704}
705#pragma GCC diagnostic pop
706
707std::ostream& operator<<(std::ostream& os, const ModelParameters& mp) {
708 os << mp.type << ' ' << mp.X0 << ' ' << mp.X1 << ' '
709 << mp.X2 << ' ' << mp.X3 << ' ';
710 return os;
711}
712
713std::istream& operator>>(std::istream& is, ModelParameters& mp) {
714 is >> mp.type >> double_imanip() >> mp.X0 >> mp.X1 >> mp.X2 >> mp.X3;
715 return is;
716}
717
719 using std::log;
720 using std::pow;
721
722 switch (type) {
723 case TemperatureModel::None:
724 return 0;
725 case TemperatureModel::T0:
726 return X0;
727 case TemperatureModel::T1:
728 return X0 * pow(T0 / T, X1);
729 case TemperatureModel::T2:
730 return X0 * pow(T0 / T, X1) * (1 + X2 * log(T / T0));
731 case TemperatureModel::T3:
732 return X0 + X1 * (T - T0);
733 case TemperatureModel::T4:
734 return (X0 + X1 * (T0 / T - 1.)) * pow(T0 / T, X2);
735 case TemperatureModel::T5:
736 return X0 * pow(T0 / T, 0.25 + 1.5 * X1);
737 case TemperatureModel::LM_AER:
738 return special_linemixing_aer(T);
739 case TemperatureModel::DPL:
740 return X0 * pow(T0 / T, X1) + X2 * pow(T0 / T, X3);
741 case TemperatureModel::POLY:
742 return X0 + X1 * T + X2 * T * T + X3 * T * T * T;
743 case TemperatureModel::FINAL: { /* Leave last */ }
744 }
745 return std::numeric_limits<Numeric>::quiet_NaN();
746}
747
749 using std::log;
750 using std::pow;
751
752 switch (type) {
753 case TemperatureModel::None:
754 return 0;
755 case TemperatureModel::T0:
756 return 1;
757 case TemperatureModel::T1:
758 return pow(T0 / T, X1);
759 case TemperatureModel::T2:
760 return pow(T0 / T, X1) * (1 + X2 * log(T / T0));
761 case TemperatureModel::T3:
762 return 1;
763 case TemperatureModel::T4:
764 return pow(T0 / T, X2);
765 case TemperatureModel::T5:
766 return pow(T0 / T, 1.5 * X1 + 0.25);
767 case TemperatureModel::LM_AER:
768 return special_linemixing_aer_dX0(T);
769 case TemperatureModel::DPL:
770 return pow(T0 / T, X1);
771 case TemperatureModel::POLY:
772 return 1;
773 case TemperatureModel::FINAL: { /* Leave last */ }
774 }
775 return std::numeric_limits<Numeric>::quiet_NaN();
776}
777
779 using std::log;
780 using std::pow;
781
782 switch (type) {
783 case TemperatureModel::None:
784 return 0;
785 case TemperatureModel::T0:
786 return 0;
787 case TemperatureModel::T1:
788 return X0 * pow(T0 / T, X1) * log(T0 / T);
789 case TemperatureModel::T2:
790 return X0 * pow(T0 / T, X1) * (X2 * log(T / T0) + 1.) * log(T0 / T);
791 case TemperatureModel::T3:
792 return (T - T0);
793 case TemperatureModel::T4:
794 return pow(T0 / T, X2) * (T0 / T - 1.);
795 case TemperatureModel::T5:
796 return 1.5 * X0 * pow(T0 / T, 1.5 * X1 + 0.25) * log(T0 / T);
797 case TemperatureModel::LM_AER:
798 return special_linemixing_aer_dX1(T);
799 case TemperatureModel::DPL:
800 return X0 * pow(T0 / T, X1) * log(T0 / T);
801 case TemperatureModel::POLY:
802 return T;
803 case TemperatureModel::FINAL: {/* Leave last */ }
804 }
805 return std::numeric_limits<Numeric>::quiet_NaN();
806}
807
809 using std::log;
810 using std::pow;
811
812 switch (type) {
813 case TemperatureModel::None:
814 return 0;
815 case TemperatureModel::T0:
816 return 0;
817 case TemperatureModel::T1:
818 return 0;
819 case TemperatureModel::T2:
820 return X0 * pow(T0 / T, X1) * log(T / T0);
821 case TemperatureModel::T3:
822 return 0;
823 case TemperatureModel::T4:
824 return pow(T0 / T, X2) * (X0 + X1 * (T0 / T - 1)) * log(T0 / T);
825 case TemperatureModel::T5:
826 return 0;
827 case TemperatureModel::LM_AER:
828 return special_linemixing_aer_dX2(T);
829 case TemperatureModel::DPL:
830 return pow(T0 / T, X3);
831 case TemperatureModel::POLY:
832 return T * T;
833 case TemperatureModel::FINAL: {/* Leave last */ }
834 }
835 return std::numeric_limits<Numeric>::quiet_NaN();
836}
837
839 using std::log;
840 using std::pow;
841
842 switch (type) {
843 case TemperatureModel::None:
844 return 0;
845 case TemperatureModel::T0:
846 return 0;
847 case TemperatureModel::T1:
848 return 0;
849 case TemperatureModel::T2:
850 return 0;
851 case TemperatureModel::T3:
852 return 0;
853 case TemperatureModel::T4:
854 return 0;
855 case TemperatureModel::T5:
856 return 0;
857 case TemperatureModel::LM_AER:
858 return special_linemixing_aer_dX3(T);
859 case TemperatureModel::DPL:
860 return X2 * pow(T0 / T, X3) * log(T0 / T);
861 case TemperatureModel::POLY:
862 return T * T * T;
863 case TemperatureModel::FINAL: {/* Leave last */ }
864 }
865 return std::numeric_limits<Numeric>::quiet_NaN();
866}
867
869 using std::log;
870 using std::pow;
871
872 switch (type) {
873 case TemperatureModel::None:
874 return 0;
875 case TemperatureModel::T0:
876 return 0;
877 case TemperatureModel::T1:
878 return -X0 * X1 * pow(T0 / T, X1) / T;
879 case TemperatureModel::T2:
880 return -X0 * X1 * pow(T0 / T, X1) * (X2 * log(T / T0) + 1.) / T +
881 X0 * X2 * pow(T0 / T, X1) / T;
882 case TemperatureModel::T3:
883 return X1;
884 case TemperatureModel::T4:
885 return -X2 * pow(T0 / T, X2) * (X0 + X1 * (T0 / T - 1.)) / T -
886 T0 * X1 * pow(T0 / T, X2) / pow(T, 2);
887 case TemperatureModel::T5:
888 return -X0 * pow(T0 / T, 1.5 * X1 + 0.25) * (1.5 * X1 + 0.25) / T;
889 case TemperatureModel::LM_AER:
890 return special_linemixing_aer_dT(T);
891 case TemperatureModel::DPL:
892 return -X0 * X1 * pow(T0 / T, X1) / T + -X2 * X3 * pow(T0 / T, X3) / T;
893 case TemperatureModel::POLY:
894 return X1 + 2 * X2 * T + 3 * X3 * T * T;
895 case TemperatureModel::FINAL: {/* Leave last */ }
896 }
897 return std::numeric_limits<Numeric>::quiet_NaN();
898}
899
901 using std::log;
902 using std::pow;
903
904 switch (type) {
905 case TemperatureModel::None:
906 return 0;
907 case TemperatureModel::T0:
908 return 0;
909 case TemperatureModel::T1:
910 return X0 * X1 * pow(T0 / T, X1) / T0;
911 case TemperatureModel::T2:
912 return X0 * X1 * pow(T0 / T, X1) * (X2 * log(T / T0) + 1.) / T0 -
913 X0 * X2 * pow(T0 / T, X1) / T0;
914 case TemperatureModel::T3:
915 return -X1;
916 case TemperatureModel::T4:
917 return X2 * pow(T0 / T, X2) * (X0 + X1 * (T0 / T - 1.)) / T0 +
918 X1 * pow(T0 / T, X2) / T;
919 case TemperatureModel::T5:
920 return X0 * pow(T0 / T, 1.5 * X1 + 0.25) * (1.5 * X1 + 0.25) / T0;
921 case TemperatureModel::LM_AER:
922 return 0;
923 case TemperatureModel::DPL:
924 return X0 * X1 * pow(T0 / T, X1) / T0 + X2 * X3 * pow(T0 / T, X3) / T0;
925 case TemperatureModel::POLY:
926 return 0;
927 case TemperatureModel::FINAL: {/* Leave last */ }
928 }
929 return std::numeric_limits<Numeric>::quiet_NaN();
930}
931
933 for (auto& data: X) {
934 Index x;
935 bif >> x >> data.X0 >> data.X1 >> data.X2 >> data.X3;
936 data.type = TemperatureModel(x);
937 }
938 return bif;
939}
940
942 for (auto& data: X) {
943 bof << Index(data.type) << data.X0 << data.X1 << data.X2 << data.X3;
944 }
945 return bof;
946}
947
948std::pair<bool, bool> SingleSpeciesModel::MatchTypes(const SingleSpeciesModel& other) const noexcept {
949 bool match=true, nullable=true;
950 for (Index i=0; i<nVars; i++) {
951 match = match and other.X[i].type == X[i].type;
952 nullable = nullable and (modelparameterEmpty(other.X[i]) or modelparameterEmpty(X[i]) or other.X[i].type == X[i].type);
953 }
954 return {match, nullable};
955}
956
957std::ostream& operator<<(std::ostream& os, const SingleSpeciesModel& ssm) {
958 for (const auto& mp : ssm.Data())
959 if (mp.type not_eq TemperatureModel::None)
960 os << mp.X0 << ' ' << mp.X1 << ' ' << mp.X2 << ' ' << mp.X3 << ' ';
961 return os;
962}
963
964std::istream& operator>>(std::istream& is, SingleSpeciesModel& ssm) {
965 for (auto& mp : ssm.Data())
966 if(mp.type not_eq TemperatureModel::None)
967 is >> double_imanip() >> mp.X0 >> mp.X1 >> mp.X2 >> mp.X3;
968 return is;
969}
970
971std::ostream& operator<<(std::ostream& os, Output x) {
972 return os << "G0: " << x.G0 << " D0: " << x.D0 << " G2: " << x.G2
973 << " D2: " << x.D2 << " FVC: " << x.FVC << " ETA: " << x.ETA
974 << " Y: " << x.Y << " G: " << x.G << " DV: " << x.DV;
975}
976
978 Numeric nself,
979 Numeric agam,
980 Numeric nair,
981 Numeric psf,
982 std::array<Numeric, 12> aer_interp) noexcept : mdata(2) {
983 mdata.front().G0() = {TemperatureModel::T1, sgam, nself, 0, 0};
984 mdata.front().D0() = {TemperatureModel::T5, psf, nair, 0, 0};
985
986 mdata.back().G0() = {TemperatureModel::T1, agam, nair, 0, 0};
987 mdata.back().D0() = {TemperatureModel::T5, psf, nair, 0, 0};
988
989 if (std::any_of(aer_interp.cbegin(), aer_interp.cend(), [](auto x){return x not_eq 0;})) {
990 mdata.front().Y().type = TemperatureModel::LM_AER;
991 mdata.front().Y().X0 = aer_interp[4];
992 mdata.front().Y().X1 = aer_interp[5];
993 mdata.front().Y().X2 = aer_interp[6];
994 mdata.front().Y().X3 = aer_interp[7];
995 mdata.front().G().type = TemperatureModel::LM_AER;
996 mdata.front().G().X0 = aer_interp[8];
997 mdata.front().G().X1 = aer_interp[9];
998 mdata.front().G().X2 = aer_interp[10];
999 mdata.front().G().X3 = aer_interp[11];
1000
1001 mdata.back().Y().type = TemperatureModel::LM_AER;
1002 mdata.back().Y().X0 = aer_interp[4];
1003 mdata.back().Y().X1 = aer_interp[5];
1004 mdata.back().Y().X2 = aer_interp[6];
1005 mdata.back().Y().X3 = aer_interp[7];
1006 mdata.back().G().type = TemperatureModel::LM_AER;
1007 mdata.back().G().X0 = aer_interp[8];
1008 mdata.back().G().X1 = aer_interp[9];
1009 mdata.back().G().X2 = aer_interp[10];
1010 mdata.back().G().X3 = aer_interp[11];
1011 }
1012}
1013
1015 Numeric nself,
1016 Numeric agam,
1017 Numeric nair,
1018 Numeric psf) {
1019 Model m(2);
1020
1021 m.Data().front().G0() = {TemperatureModel::T1, sgam, nself, 0, 0};
1022 m.Data().front().D0() = {TemperatureModel::T0, psf, 0, 0, 0};
1023
1024 m.Data().back().G0() = {TemperatureModel::T1, agam, nair, 0, 0};
1025 m.Data().back().D0() = {TemperatureModel::T0, psf, 0, 0, 0};
1026
1027 return m;
1028}
1029
1031 Numeric nself,
1032 Numeric agam,
1033 Numeric nair,
1034 Numeric psf,
1035 std::array<Numeric, 12> aer_interp) {
1036 Model m(2);
1037
1038 m.Data().front().G0() = {TemperatureModel::T1, sgam, nself, 0, 0};
1039 m.Data().front().D0() = {TemperatureModel::T0, psf, 0, 0, 0};
1040
1041 m.Data().back().G0() = {TemperatureModel::T1, agam, nair, 0, 0};
1042 m.Data().back().D0() = {TemperatureModel::T0, psf, 0, 0, 0};
1043
1044 if (std::any_of(aer_interp.cbegin(), aer_interp.cend(), [](auto x){return x not_eq 0;})) {
1045 m.Data().front().Y().type = TemperatureModel::LM_AER;
1046 m.Data().front().Y().X0 = aer_interp[4];
1047 m.Data().front().Y().X1 = aer_interp[5];
1048 m.Data().front().Y().X2 = aer_interp[6];
1049 m.Data().front().Y().X3 = aer_interp[7];
1050 m.Data().front().G().type = TemperatureModel::LM_AER;
1051 m.Data().front().G().X0 = aer_interp[8];
1052 m.Data().front().G().X1 = aer_interp[9];
1053 m.Data().front().G().X2 = aer_interp[10];
1054 m.Data().front().G().X3 = aer_interp[11];
1055
1056 m.Data().back().Y().type = TemperatureModel::LM_AER;
1057 m.Data().back().Y().X0 = aer_interp[4];
1058 m.Data().back().Y().X1 = aer_interp[5];
1059 m.Data().back().Y().X2 = aer_interp[6];
1060 m.Data().back().Y().X3 = aer_interp[7];
1061 m.Data().back().G().type = TemperatureModel::LM_AER;
1062 m.Data().back().G().X0 = aer_interp[8];
1063 m.Data().back().G().X1 = aer_interp[9];
1064 m.Data().back().G().X2 = aer_interp[10];
1065 m.Data().back().G().X3 = aer_interp[11];
1066 }
1067
1068 return m;
1069}
1070
1071bool Model::OK(Type type, bool self, bool bath,
1072 const std::size_t nspecies) const noexcept {
1073 Index n = mdata.size();
1074 Index m = Index(self) + Index(bath);
1075 bool needs_any = type not_eq Type::DP;
1076 return not (n not_eq Index(nspecies) or m > n or (needs_any and not n));
1077}
1078
1079#define LSPC(XVAR, PVAR) \
1080 Numeric Model::XVAR( \
1081 Numeric T, Numeric T0, Numeric P [[maybe_unused]], \
1082 const ConstVectorView& vmrs) \
1083 const noexcept { \
1084 return PVAR * \
1085 std::inner_product( \
1086 mdata.cbegin(), \
1087 mdata.cend(), \
1088 vmrs.begin(), \
1089 0.0, \
1090 std::plus<>(), \
1091 [&](auto& x, auto vmr) -> Numeric { \
1092 return vmr * x.XVAR().at(T, T0); \
1093 }); \
1094 }
1095 LSPC(G0, P)
1096 LSPC(D0, P)
1097 LSPC(G2, P)
1098 LSPC(D2, P)
1099 LSPC(FVC, P)
1100 LSPC(ETA, 1)
1101 LSPC(Y, P)
1102 LSPC(G, P* P)
1103 LSPC(DV, P* P)
1104#undef LSPC
1105
1106#define LSPCV(XVAR, PVAR) \
1107 Numeric Model::d##XVAR##_dVMR(Numeric T, \
1108 Numeric T0, \
1109 Numeric P [[maybe_unused]], \
1110 const Index deriv_pos) const noexcept { \
1111 if (deriv_pos not_eq -1) \
1112 return PVAR * mdata[deriv_pos].XVAR().at(T, T0); \
1113 return 0; \
1114 }
1115 LSPCV(G0, P)
1116 LSPCV(D0, P)
1117 LSPCV(G2, P)
1118 LSPCV(D2, P)
1119 LSPCV(FVC, P)
1120 LSPCV(ETA, 1)
1121 LSPCV(Y, P)
1122 LSPCV(G, P* P)
1123 LSPCV(DV, P* P)
1124#undef LSPCV
1125
1126#define LSPCT(XVAR, PVAR) \
1127 Numeric Model::d##XVAR##_dT( \
1128 Numeric T, Numeric T0, Numeric P [[maybe_unused]], \
1129 const ConstVectorView& vmrs) \
1130 const noexcept { \
1131 return PVAR * \
1132 std::inner_product( \
1133 mdata.cbegin(), \
1134 mdata.cend(), \
1135 vmrs.begin(), \
1136 0.0, \
1137 std::plus<>(), \
1138 [&](auto& x, auto vmr) -> Numeric { \
1139 return vmr * x.XVAR().dT(T, T0); \
1140 }); \
1141 }
1142 LSPCT(G0, P)
1143 LSPCT(D0, P)
1144 LSPCT(G2, P)
1145 LSPCT(D2, P)
1146 LSPCT(FVC, P)
1147 LSPCT(ETA, 1)
1148 LSPCT(Y, P)
1149 LSPCT(G, P* P)
1150 LSPCT(DV, P* P)
1151#undef LSPCT
1152
1153#define LSPDC(XVAR, DERIV, PVAR) \
1154 Numeric Model::d##XVAR##_##DERIV(Numeric T, \
1155 Numeric T0, \
1156 Numeric P [[maybe_unused]], \
1157 Index deriv_pos, \
1158 const ConstVectorView& vmrs) const noexcept { \
1159 if (deriv_pos not_eq -1) \
1160 return vmrs[deriv_pos] * PVAR * \
1161 mdata[deriv_pos].XVAR().DERIV(T, T0); \
1162 return 0; \
1163 }
1164 LSPDC(G0, dT0, P)
1165 LSPDC(G0, dX0, P)
1166 LSPDC(G0, dX1, P)
1167 LSPDC(G0, dX2, P)
1168 LSPDC(G0, dX3, P)
1169 LSPDC(D0, dT0, P)
1170 LSPDC(D0, dX0, P)
1171 LSPDC(D0, dX1, P)
1172 LSPDC(D0, dX2, P)
1173 LSPDC(D0, dX3, P)
1174 LSPDC(G2, dT0, P)
1175 LSPDC(G2, dX0, P)
1176 LSPDC(G2, dX1, P)
1177 LSPDC(G2, dX2, P)
1178 LSPDC(G2, dX3, P)
1179 LSPDC(D2, dT0, P)
1180 LSPDC(D2, dX0, P)
1181 LSPDC(D2, dX1, P)
1182 LSPDC(D2, dX2, P)
1183 LSPDC(D2, dX3, P)
1184 LSPDC(FVC, dT0, P)
1185 LSPDC(FVC, dX0, P)
1186 LSPDC(FVC, dX1, P)
1187 LSPDC(FVC, dX2, P)
1188 LSPDC(FVC, dX3, P)
1189 LSPDC(ETA, dT0, 1)
1190 LSPDC(ETA, dX0, 1)
1191 LSPDC(ETA, dX1, 1)
1192 LSPDC(ETA, dX2, 1)
1193 LSPDC(ETA, dX3, 1)
1194 LSPDC(Y, dT0, P)
1195 LSPDC(Y, dX0, P)
1196 LSPDC(Y, dX1, P)
1197 LSPDC(Y, dX2, P)
1198 LSPDC(Y, dX3, P)
1199 LSPDC(G, dT0, P* P)
1200 LSPDC(G, dX0, P* P)
1201 LSPDC(G, dX1, P* P)
1202 LSPDC(G, dX2, P* P)
1203 LSPDC(G, dX3, P* P)
1204 LSPDC(DV, dT0, P* P)
1205 LSPDC(DV, dX0, P* P)
1206 LSPDC(DV, dX1, P* P)
1207 LSPDC(DV, dX2, P* P)
1208 LSPDC(DV, dX3, P* P)
1209#undef LSPDC
1210
1212 Numeric T0,
1213 Numeric P,
1214 const ConstVectorView& vmrs) const noexcept {
1215 return {G0(T, T0, P, vmrs),
1216 D0(T, T0, P, vmrs),
1217 G2(T, T0, P, vmrs),
1218 D2(T, T0, P, vmrs),
1219 FVC(T, T0, P, vmrs),
1220 ETA(T, T0, P, vmrs),
1221 Y(T, T0, P, vmrs),
1222 G(T, T0, P, vmrs),
1223 DV(T, T0, P, vmrs)};
1224}
1225
1227 Numeric T0,
1228 Numeric P,
1229 size_t k) const noexcept {
1230 return {P * mdata[k].G0().at(T, T0),
1231 P * mdata[k].D0().at(T, T0),
1232 P * mdata[k].G2().at(T, T0),
1233 P * mdata[k].D2().at(T, T0),
1234 P * mdata[k].FVC().at(T, T0),
1235 mdata[k].ETA().at(T, T0),
1236 P * mdata[k].Y().at(T, T0),
1237 P * P * mdata[k].G().at(T, T0),
1238 P * P * mdata[k].DV().at(T, T0)};
1239}
1240
1242 Numeric T0,
1243 Numeric P,
1244 ConstVectorView vmrs) const noexcept {
1245 return {dG0_dT(T, T0, P, vmrs),
1246 dD0_dT(T, T0, P, vmrs),
1247 dG2_dT(T, T0, P, vmrs),
1248 dD2_dT(T, T0, P, vmrs),
1249 dFVC_dT(T, T0, P, vmrs),
1250 dETA_dT(T, T0, P, vmrs),
1251 dY_dT(T, T0, P, vmrs),
1252 dG_dT(T, T0, P, vmrs),
1253 dDV_dT(T, T0, P, vmrs)};
1254}
1255
1256Output Model::GetVMRDerivs(Numeric T, Numeric T0, Numeric P, const Index pos) const noexcept {
1257 return {dG0_dVMR(T, T0, P, pos),
1258 dD0_dVMR(T, T0, P, pos),
1259 dG2_dVMR(T, T0, P, pos),
1260 dD2_dVMR(T, T0, P, pos),
1261 dFVC_dVMR(T, T0, P, pos),
1262 dETA_dVMR(T, T0, P, pos),
1263 dY_dVMR(T, T0, P, pos),
1264 dG_dVMR(T, T0, P, pos),
1265 dDV_dVMR(T, T0, P, pos)};
1266}
1267
1269 Numeric T0,
1270 Numeric P,
1271 Index pos,
1272 const Vector& vmrs,
1273 Jacobian::Line deriv) const noexcept {
1274 if (pos < 0) return 0;
1275
1276#define RETURNINTERNALDERIVATIVE(TYPE) \
1277 case Jacobian::Line::Shape##TYPE##X0: \
1278 return d##TYPE##_dX0(T, T0, P, pos, vmrs); \
1279 case Jacobian::Line::Shape##TYPE##X1: \
1280 return d##TYPE##_dX1(T, T0, P, pos, vmrs); \
1281 case Jacobian::Line::Shape##TYPE##X2: \
1282 return d##TYPE##_dX2(T, T0, P, pos, vmrs); \
1283 case Jacobian::Line::Shape##TYPE##X3: \
1284 return d##TYPE##_dX3(T, T0, P, pos, vmrs)
1285 switch (deriv) {
1295 default:
1296 return std::numeric_limits<Numeric>::quiet_NaN();
1297 }
1298#undef RETURNINTERNALDERIVATIVE
1299}
1300
1302 mdata.erase(mdata.begin() + i);
1303 specs.erase(specs.begin() + i);
1304}
1305
1307 mdata.erase(mdata.begin() + i);
1308 specs.erase(specs.begin() + i);
1309}
1310
1312 for (auto& ssm : mdata) {
1313 ssm.Y() = x.Y();
1314 ssm.G() = x.G();
1315 ssm.DV() = x.DV();
1316 }
1317}
1318
1319std::pair<bool, bool> Model::Match(const Model& other) const noexcept {
1320 const Index n = nelem();
1321 if (other.nelem() not_eq n) return {false, false};
1322
1323 bool match = true, nullable = true;
1324 for (Index i=0; i<n; i++) {
1325 const auto x = mdata[i].MatchTypes(other[i]);
1326 match = match and x.first;
1327 nullable = nullable and x.second;
1328 }
1329
1330 return {match, nullable};
1331}
1332
1334 for (auto& data: mdata)
1335 data.read(bif);
1336 return bif;
1337}
1338
1340 for (auto& data: mdata)
1341 data.write(bof);
1342 return bof;
1343}
1344
1345namespace LegacyLineFunctionData {
1346std::vector<Variable> lineshapetag2variablesvector(String type) {
1347 if (type == String("DP"))
1348 return {};
1349 if (type == String("LP"))
1350 return {Variable::G0, Variable::D0};
1351 if (type == String("VP"))
1352 return {Variable::G0, Variable::D0};
1353 if (type == String("SDVP"))
1355 if (type == String("HTP"))
1356 return {Variable::G0,
1363 "Type: ", type, ", is not accepted. "
1364 "See documentation for accepted types\n")
1365}
1366
1367std::vector<Variable> linemixingtag2variablesvector(String type) {
1368 if (type == "#")
1369 return {};
1370 if (type == "LM1")
1371 return {Variable::Y};
1372 if (type == "LM2")
1373 return {Variable::Y, Variable::G, Variable::DV};
1374 if (type == "INT")
1375 return {};
1376 if (type == "ConstG")
1377 return {Variable::G};
1379 "Type: ", type, ", is not accepted. "
1380 "See documentation for accepted types\n")
1381}
1382} // namespace LegacyLineFunctionData
1383
1384namespace LegacyLineMixingData {
1386 if (type == "NA") // The standard case
1387 return TypeLM::LM_NONE;
1388 if (type == "LL") // The LBLRTM case
1389 return TypeLM::LM_LBLRTM;
1390 if (type == "NR") // The LBLRTM O2 non-resonant case
1392 if (type == "L2") // The 2nd order case
1393 return TypeLM::LM_2NDORDER;
1394 if (type == "L1") // The 2nd order case
1395 return TypeLM::LM_1STORDER;
1396 if (type == "BB") // The band class
1397 return TypeLM::LM_BYBAND;
1399 "Type: ", type, ", is not accepted. "
1400 "See documentation for accepted types\n")
1401}
1402} // namespace LegacyLineMixingData
1403
1404
1405namespace LegacyPressureBroadeningData {
1407 if (type == "NA") // The none case
1408 return TypePB::PB_NONE;
1409 if (type == "N2") // Air Broadening is N2 broadening mostly...
1411 if (type == "WA") // Water and Air Broadening
1413 if (type == "AP") // Planetary broadening
1416 "Type: ", type, ", is not accepted. "
1417 "See documentation for accepted types\n")
1418}
1419
1423 (qid.Species() == Species::fromShortName("N2") or
1424 qid.Species() == Species::fromShortName("O2") or
1425 qid.Species() == Species::fromShortName("H2O") or
1426 qid.Species() == Species::fromShortName("CO2") or
1427 qid.Species() == Species::fromShortName("H2") or
1428 qid.Species() == Species::fromShortName("He")))
1429 return true;
1430 if (t == TypePB::PB_AIR_AND_WATER_BROADENING and qid.Species() == Species::fromShortName("H2O"))
1431 return true;
1432 return false;
1433}
1434} // namespace LegacyPressureBroadeningData
1435} // namespace LineShape
Common ARTS conversions.
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Vector.
Definition: matpackI.h:512
Numeric sum() const ARTS_NOEXCEPT
The sum of all elements of a Vector.
Definition: matpackI.cc:48
Main line shape model class.
void SetLineMixingModel(SingleSpeciesModel x)
Sets the same line mixing model to all species.
Output GetVMRDerivs(Numeric T, Numeric T0, Numeric P, const Index pos) const noexcept
Derivative of GetParams(...) wrt VMR.
void Remove(Index i, ArrayOfSpeciesTag &specs)
Remove species and data at position.
std::vector< SingleSpeciesModel > mdata
Model(Index n=0) noexcept
Default init just sets the size.
Numeric GetInternalDeriv(Numeric T, Numeric T0, Numeric P, Index pos, const Vector &vmrs, Jacobian::Line deriv) const noexcept
Derivative of GetParams(...) wrt Coefficient.
Output GetTemperatureDerivs(Numeric T, Numeric T0, Numeric P, ConstVectorView vmrs) const noexcept
Derivative of GetParams(...) wrt T.
bool OK(Type type, bool self, bool bath, const std::size_t nspecies) const noexcept
The Model is good to use.
bifstream & read(bifstream &bif)
Binary read for Model.
Output GetParams(Numeric T, Numeric T0, Numeric P, const ConstVectorView &vmrs) const noexcept
Compute all shape parameters.
bofstream & write(bofstream &bof) const
Binary write for Model.
const std::vector< SingleSpeciesModel > & Data() const noexcept
The line shape model data.
std::pair< bool, bool > Match(const Model &other) const noexcept
Compute the line shape parameters for a single broadening species.
constexpr std::array< ModelParameters, nVars > & Data() noexcept
Get internal Data reference.
std::array< ModelParameters, nVars > X
bofstream & write(bofstream &bof) const
Binary write for SingleSpeciesModel.
bifstream & read(bifstream &bif)
Binary read for SingleSpeciesModel.
std::pair< bool, bool > MatchTypes(const SingleSpeciesModel &other) const noexcept
The Vector class.
Definition: matpackI.h:899
Binary output file stream class.
Definition: bifstream.h:43
Binary output file stream class.
Definition: bofstream.h:42
Input manipulator class for doubles to enable nan and inf parsing.
Definition: double_imanip.h:42
Index nelem() const
Definition: mystring.h:189
#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 std::string_view toString(EnergyLevelMapType x) noexcept
constexpr void check_enum_error(EnumType type, Messages ... args)
Checks if the enum class type is good and otherwise throws an error message composed by variadic inpu...
Definition: enums.h:85
constexpr bool good_enum(EnumType x) noexcept
Checks if the enum number is good.
Definition: enums.h:22
Array< Species::Species > ArrayOfSpecies
#define temp
#define LSPCV(XVAR, PVAR)
#define ReturnJacPropMatType(ID)
#define RETURNINTERNALDERIVATIVE(TYPE)
#define LSPCT(XVAR, PVAR)
#define LSPC(XVAR, PVAR)
#define LSPDC(XVAR, DERIV, PVAR)
Jacobian::Line select_derivativeLineShape(const String &var, const String &coeff)
Return the derivative type based on string input.
Contains the line shape namespace.
Implementation of Matrix, Vector, and such stuff.
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
my_basic_string< char > String
The String type for ARTS.
Definition: mystring.h:216
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric k
Boltzmann constant convenience name [J/K].
std::vector< Variable > linemixingtag2variablesvector(String type)
Line mixing models from string.
std::vector< Variable > lineshapetag2variablesvector(String type)
Line shape models from string.
constexpr Index temperaturemodel2legacynelem(TemperatureModel type) noexcept
Length per variable for temperature model.
LegacyLineMixingData::TypeLM string2typelm(String type)
Line mixing types from string.
constexpr Index typelm2nelem(LegacyLineMixingData::TypeLM type)
Line mixing types to number.
Model vector2modellm(Vector x, LegacyLineMixingData::TypeLM type)
LineShape::Model from legacy input vector.
TypeLM
Line mixing types that used to exist.
Index self_listed(const QuantumIdentifier &qid, LegacyPressureBroadeningData::TypePB t)
Pressure broadening if self exist.
void vector2modelpb(LineShape::Type &mtype, bool &self, bool &bath, Model &m, ArrayOfSpecies &species, Vector x, LegacyPressureBroadeningData::TypePB type, bool self_in_list, Species::Species self_spec)
LineShape::Model from legacy input vector.
constexpr Index typepb2nelem(LegacyPressureBroadeningData::TypePB type)
Pressure broadening types to number of elements.
LegacyPressureBroadeningData::TypePB string2typepb(String type)
Pressure broadening types from string.
TypePB
Pressure broadening types that used to exist.
Computations of line shape derived parameters.
Definition: lineshape.cc:23
Model lblrtm_model(Numeric sgam, Numeric nself, Numeric agam, Numeric nair, Numeric psf, std::array< Numeric, 12 > aer_interp)
Numeric & SingleModelParameter(ModelParameters &mp, const String &type)
Get a coefficient from ModelParameters by name.
Model MetaData2ModelShape(const String &s)
std::istream & from_linefunctiondata(std::istream &data, Type &type, bool &self, bool &bath, Model &m, ArrayOfSpecies &species)
constexpr bool modelparameterEmpty(const ModelParameters mp) noexcept
Model hitran_model(Numeric sgam, Numeric nself, Numeric agam, Numeric nair, Numeric psf)
Vector mass(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const ArrayOfSpecies &lineshape_species, const SpeciesIsotopologueRatios &ir) ARTS_NOEXCEPT
Returns a mass vector for this model's main calculations.
std::istream & from_pressurebroadeningdata(std::istream &data, LineShape::Type &type, bool &self, bool &bath, Model &m, ArrayOfSpecies &species, const QuantumIdentifier &qid)
Legacy reading of old deprecated PressureBroadeningData class.
ArrayOfString ModelMetaDataArray(const LineShape::Model &m, const bool self, const ArrayOfSpecies &sts, const Numeric T0)
String modelparameters2metadata(const ModelParameters mp, const Numeric T0)
std::istream & operator>>(std::istream &is, Model &m)
String ModelShape2MetaData(const Model &m)
std::istream & from_artscat4(std::istream &is, Type &type, bool &self, bool &bath, Model &m, ArrayOfSpecies &species, const QuantumIdentifier &qid)
std::ostream & operator<<(std::ostream &os, const Model &m)
std::istream & from_linemixingdata(std::istream &data, Model &lsc)
Legacy reading of old deprecated LineMixingData class.
constexpr Index nVars
Current max number of line shape variables.
Vector vmrs(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const ArrayOfSpecies &lineshape_species) ARTS_NOEXCEPT
Returns a VMR vector for this model's main calculations.
VectorView var(VectorView var, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the variance of the ranged ys.
Definition: raw.cc:159
constexpr Numeric mean_mass(Species spec, const IsotopologueRatios &ir) noexcept
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
#define N
Definition: rng.cc:164
Coefficients and temperature model for SingleSpeciesModel.
Numeric dT0(Numeric T, Numeric T0) const noexcept
Numeric dX0(Numeric T, Numeric T0) const noexcept
Numeric dT(Numeric T, Numeric T0) const noexcept
Numeric at(Numeric T, Numeric T0) const noexcept
Numeric dX1(Numeric T, Numeric T0) const noexcept
Numeric dX3(Numeric T, Numeric T0) const noexcept
Numeric dX2(Numeric T, Numeric T0) const noexcept
A logical struct for global quantum numbers with species identifiers.
Species::Species Species() const noexcept
Species::IsotopeRecord Isotopologue() const noexcept
#define v