ARTS 2.5.4 (git: 4c0d3b4d)
absorptionlines.cc
Go to the documentation of this file.
1/* Copyright (C) 2019
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
32#include "absorptionlines.h"
33
34#include "absorption.h"
35#include "constants.h"
36#include "debug.h"
37#include "enums.h"
38#include "file.h"
39#include "hitran_species.h"
40#include "jpl_species.h"
41#include "linescaling.h"
42#include "quantum_numbers.h"
43
44LineShape::Output Absorption::Lines::ShapeParameters(size_t k, Numeric T, Numeric P, const Vector& vmrs) const noexcept {
45 auto x = lines[k].lineshape.GetParams(T, T0, P, vmrs);
46
47 if (not DoLineMixing(P)) x.Y = x.G = x.DV = 0;
48
49 return x;
50}
51
52LineShape::Output Absorption::Lines::ShapeParameters(size_t k, Numeric T, Numeric P, size_t m) const noexcept {
53 return lines[k].lineshape.GetParams(T, T0, P, m);
54}
55
56LineShape::Output Absorption::Lines::ShapeParameters_dT(size_t k, Numeric T, Numeric P, const Vector& vmrs) const noexcept {
57 auto x = lines[k].lineshape.GetTemperatureDerivs(T, T0, P, vmrs);
58
59 if (not DoLineMixing(P)) x.Y = x.G = x.DV = 0;
60
61 return x;
62}
63
64Index Absorption::Lines::LineShapePos(const Species::Species spec) const noexcept {
65 // Is always first if this is self and self broadening exists
66 if(selfbroadening and spec == quantumidentity.Species()) {
67 return 0;
68 }
69
70 // First and last might be artificial so they should not be checked
71 const Index s = selfbroadening;
72 const Index e = broadeningspecies.nelem() - bathbroadening;
73 for(Index i=s; i<e; i++) {
74 if(spec == broadeningspecies[i]) {
75 return i;
76 }
77 }
78
79 // At this point, the ID is not explicitly among the broadeners, but bath broadening means its VMR still might matter
80 if(bathbroadening)
81 return broadeningspecies.nelem()-1;
82 return -1;
83}
84
85LineShape::Output Absorption::Lines::ShapeParameters_dVMR(size_t k, Numeric T, Numeric P, const QuantumIdentifier& vmr_qid) const noexcept {
86 const Index pos=LineShapePos(vmr_qid);
87 if (pos >= 0)
88 return lines[k].lineshape.GetVMRDerivs(T, T0, P, pos);
89 return LineShape::Output{};
90}
91
92Numeric Absorption::Lines::ShapeParameter_dInternal(size_t k, Numeric T, Numeric P, const Vector& vmrs, const RetrievalQuantity& derivative) const noexcept {
93 const auto self = derivative.Mode() == LineShape::self_broadening;
94 const auto bath = derivative.Mode() == LineShape::bath_broadening;
95 const auto& ls = lines[k].lineshape;
96
97 if(derivative.QuantumIdentity().Species() not_eq Species() or
98 derivative.QuantumIdentity().Isotopologue() not_eq Isotopologue())
99 return 0;
100 if(self and selfbroadening)
101 return ls.GetInternalDeriv(
102 T, T0, P, 0, vmrs, derivative.LineType());
103 if(self)
104 return ls.GetInternalDeriv(
105 T, T0, P, LineShapePos(derivative.QuantumIdentity().Species()), vmrs, derivative.LineType());
106 if(bath and bathbroadening)
107 return ls.GetInternalDeriv(
108 T, T0, P, ls.nelem() - 1, vmrs, derivative.LineType());
109 if(bath)
110 return 0;
111 return ls.GetInternalDeriv(
112 T, T0, P, LineShapePos(derivative.QuantumIdentity().Species()), vmrs, derivative.LineType());
113}
114
116 // Default data and values for this type
118 data.selfbroadening = true;
119 data.bathbroadening = true;
120 data.lineshapetype = LineShape::Type::VP;
121 data.species.resize(2);
122
123 // This always contains the rest of the line to parse. At the
124 // beginning the entire line. Line gets shorter and shorter as we
125 // continue to extract stuff from the beginning.
126 String line;
127
128 // Look for more comments?
129 bool comment = true;
130
131 while (comment) {
132 // Return true if eof is reached:
133 if (is.eof()) return data;
134
135 // Throw runtime_error if stream is bad:
136 ARTS_USER_ERROR_IF (!is, "Stream bad.");
137
138 // Read line from file into linebuffer:
139 getline(is, line);
140
141 // It is possible that we were exactly at the end of the file before
142 // calling getline. In that case the previous eof() was still false
143 // because eof() evaluates only to true if one tries to read after the
144 // end of the file. The following check catches this.
145 if (line.nelem() == 0 && is.eof()) return data;
146
147 // @ as first character marks catalogue entry
148 char c;
149 extract(c, line, 1);
150
151 // check for empty line
152 if (c == '@') {
153 comment = false;
154 }
155 }
156
157 // read the arts identifier String
158 istringstream icecream(line);
159
160 String artsid;
161 icecream >> artsid;
162
163 if (artsid.length() != 0) {
164 // Set the species
165 const auto isotopologue = Species::Tag(Species::update_isot_name(artsid));
166 ARTS_USER_ERROR_IF (isotopologue.is_joker() or isotopologue.type not_eq Species::TagType::Plain,
167 "A line catalog species can only be of the form \"Plain\", meaning it\nhas the form SPECIES-ISONUM.\n"
168 "Your input contains: ", artsid, ". which we cannot interpret as a plain species")
169 data.quantumidentity = QuantumIdentifier{Species::find_species_index(isotopologue.Isotopologue())};
170
171 // Extract center frequency:
172 icecream >> double_imanip() >> data.line.F0;
173
174 Numeric psf;
175 // Extract pressure shift:
176 icecream >> double_imanip() >> psf;
177
178 // Extract intensity:
179 icecream >> double_imanip() >> data.line.I0;
180
181 // Extract reference temperature for Intensity in K:
182 icecream >> double_imanip() >> data.T0;
183
184 // Extract lower state energy:
185 icecream >> data.line.E0;
186
187 // Extract air broadening parameters:
188 Numeric agam, sgam;
189 icecream >> double_imanip() >> agam;
190 icecream >> double_imanip() >> sgam;
191
192 // Extract temperature coefficient of broadening parameters:
193 Numeric nair, nself;
194 icecream >> double_imanip() >> nair;
195 icecream >> double_imanip() >> nself;
196
197 // Extract reference temperature for broadening parameter in K:
198 Numeric tgam;
199 icecream >> double_imanip() >> tgam;
200
201 // Extract the aux parameters:
202 Index naux;
203 icecream >> naux;
204
205 // resize the aux array and read it
206 ArrayOfNumeric maux;
207 maux.resize(naux);
208
209 for (Index j = 0; j < naux; j++) {
210 icecream >> maux[j];
211 //cout << "maux" << j << " = " << maux[j] << "\n";
212 }
213
214 // Extract accuracies:
215 try {
216 Numeric unused_numeric;
217 icecream >> double_imanip() >> unused_numeric /*mdf*/;
218 icecream >> double_imanip() >> unused_numeric /*mdi0*/;
219 icecream >> double_imanip() >> unused_numeric /*dagam*/;
220 icecream >> double_imanip() >> unused_numeric /*dsgam*/;
221 icecream >> double_imanip() >> unused_numeric /*dnair*/;
222 icecream >> double_imanip() >> unused_numeric /*dnself*/;
223 icecream >> double_imanip() >> unused_numeric /*dpsf*/;
224 } catch (const std::runtime_error&) {
225 }
226
227 // Fix if tgam is different from ti0
228 if (tgam != data.T0) {
229 agam = agam * pow(tgam / data.T0, nair);
230 sgam = sgam * pow(tgam / data.T0, nself);
231 psf = psf * pow(tgam / data.T0, (Numeric).25 + (Numeric)1.5 * nair);
232 }
233
234 // Set line shape computer
235 data.line.lineshape = LineShape::Model(sgam, nself, agam, nair, psf);
236 }
237
238 // That's it!
239 data.bad = false;
240 return data;
241}
242
244 // Default data and values for this type
246 data.selfbroadening = true;
247 data.bathbroadening = false;
248 data.lineshapetype = LineShape::Type::VP;
249
250 // This always contains the rest of the line to parse. At the
251 // beginning the entire line. Line gets shorter and shorter as we
252 // continue to extract stuff from the beginning.
253 String line;
254
255 // Look for more comments?
256 bool comment = true;
257
258 while (comment) {
259 // Return true if eof is reached:
260 if (is.eof()) return data;
261
262 // Throw runtime_error if stream is bad:
263 ARTS_USER_ERROR_IF (!is, "Stream bad.");
264
265 // Read line from file into linebuffer:
266 getline(is, line);
267
268 // It is possible that we were exactly at the end of the file before
269 // calling getline. In that case the previous eof() was still false
270 // because eof() evaluates only to true if one tries to read after the
271 // end of the file. The following check catches this.
272 if (line.nelem() == 0 && is.eof()) return data;
273
274 // @ as first character marks catalogue entry
275 char c;
276 extract(c, line, 1);
277
278 // check for empty line
279 if (c == '@') {
280 comment = false;
281 }
282 }
283
284 // read the arts identifier String
285 istringstream icecream(line);
286
287 String artsid;
288 icecream >> artsid;
289
290 if (artsid.length() != 0) {
291
292 // Set line ID
293 const auto isotopologue = Species::Tag(Species::update_isot_name(artsid));
294 ARTS_USER_ERROR_IF (isotopologue.is_joker() or isotopologue.type not_eq Species::TagType::Plain,
295 "A line catalog species can only be of the form \"Plain\", meaning it\nhas the form SPECIES-ISONUM.\n"
296 "Your input contains: ", artsid, ". which we cannot interpret as a plain species")
297 data.quantumidentity = QuantumIdentifier{Species::find_species_index(isotopologue.Isotopologue())};
298
299 // Extract center frequency:
300 icecream >> double_imanip() >> data.line.F0;
301
302 // Extract intensity:
303 icecream >> double_imanip() >> data.line.I0;
304
305 // Extract reference temperature for Intensity in K:
306 icecream >> double_imanip() >> data.T0;
307
308 // Extract lower state energy:
309 icecream >> double_imanip() >> data.line.E0;
310
311 // Extract Einstein A-coefficient:
312 icecream >> double_imanip() >> data.line.A;
313
314 // Extract upper state stat. weight:
315 icecream >> double_imanip() >> data.line.gupp;
316
317 // Extract lower state stat. weight:
318 icecream >> double_imanip() >> data.line.glow;
319
320 LineShape::from_artscat4(icecream,
321 data.lineshapetype,
322 data.selfbroadening,
323 data.bathbroadening,
324 data.line.lineshape,
325 data.species,
326 data.quantumidentity);
327 }
328
329 // That's it!
330 data.bad = false;
331 return data;
332}
333
335 // Default data and values for this type
337
338 LineShape::Model line_mixing_model;
339 bool lmd_found = false;
340
341 // This always contains the rest of the line to parse. At the
342 // beginning the entire line. Line gets shorter and shorter as we
343 // continue to extract stuff from the beginning.
344 String line;
345
346 // Look for more comments?
347 bool comment = true;
348
349 while (comment) {
350 // Return true if eof is reached:
351 if (is.eof()) return data;
352
353 // Throw runtime_error if stream is bad:
354 ARTS_USER_ERROR_IF (!is, "Stream bad.");
355
356 // Read line from file into linebuffer:
357 getline(is, line);
358
359 // It is possible that we were exactly at the end of the file before
360 // calling getline. In that case the previous eof() was still false
361 // because eof() evaluates only to true if one tries to read after the
362 // end of the file. The following check catches this.
363 if (line.nelem() == 0 && is.eof()) return data;
364
365 // @ as first character marks catalogue entry
366 char c;
367 extract(c, line, 1);
368
369 // check for empty line
370 if (c == '@') {
371 comment = false;
372 }
373 }
374
375 // read the arts identifier String
376 istringstream icecream(line);
377
378 try {
379 String artsid;
380 icecream >> artsid;
381
382 if (artsid.length() != 0) {
383 // Set line ID:
384 const auto isotopologue = Species::Tag(Species::update_isot_name(artsid));
385 ARTS_USER_ERROR_IF (isotopologue.is_joker() or isotopologue.type not_eq Species::TagType::Plain,
386 "A line catalog species can only be of the form \"Plain\", meaning it\nhas the form SPECIES-ISONUM.\n"
387 "Your input contains: ", artsid, ". which we cannot interpret as a plain species")
388 data.quantumidentity = QuantumIdentifier{Species::find_species_index(isotopologue.Isotopologue())};
389
390 // Extract center frequency:
391 icecream >> double_imanip() >> data.line.F0;
392
393 // Extract intensity:
394 icecream >> double_imanip() >> data.line.I0;
395
396 // Extract reference temperature for Intensity in K:
397 icecream >> double_imanip() >> data.T0;
398
399 // Extract lower state energy:
400 icecream >> double_imanip() >> data.line.E0;
401
402 // Extract Einstein A-coefficient:
403 icecream >> double_imanip() >> data.line.A;
404
405 // Extract upper state stat. weight:
406 icecream >> double_imanip() >> data.line.gupp;
407
408 // Extract lower state stat. weight:
409 icecream >> double_imanip() >> data.line.glow;
410
411 String token;
412 Index nelem;
413 icecream >> token;
414
415 while (icecream) {
416 // Read pressure broadening (LEGACY)
417 if (token == "PB") {
419 icecream,
420 data.lineshapetype,
421 data.selfbroadening,
422 data.bathbroadening,
423 data.line.lineshape,
424 data.species,
425 data.quantumidentity);
426 icecream >> token;
427 } else if (token == "QN") {
428 // Quantum numbers
429
430 icecream >> token;
431 ARTS_USER_ERROR_IF (token != "UP",
432 "Unknown quantum number tag: ", token)
433
434 icecream >> token;
435 String r;
436 while (icecream and token not_eq "LO") {
437 auto qn = Quantum::Number::toType(token);
438 icecream >> r;
439 icecream >> token;
440
441 if (good_enum(qn)) {
442 if (data.quantumidentity.val.has(qn)) {
444 val.set(r, true);
445 data.quantumidentity.val.set(val);
446 } else {
447 data.quantumidentity.val.add(qn).set(r, true);
448 }
449 }
450 }
451
452 ARTS_USER_ERROR_IF (!is || token != "LO",
453 "Error in catalog. Lower quantum number tag 'LO' not found.")
454
455 icecream >> token;
456 auto qn = Quantum::Number::toType(token);
457 while (icecream and not (token == "LM" or token == "LF" or token == "ZM" or token == "LSM" or token == "PB")) {
458 icecream >> r;
459 icecream >> token;
460
461 if (good_enum(qn)) {
462 if (data.quantumidentity.val.has(qn)) {
464 val.set(r, false);
465 data.quantumidentity.val.set(val);
466 } else {
467 data.quantumidentity.val.add(qn).set(r, false);
468 }
469 }
470
471 qn = Quantum::Number::toType(token);
472 }
473 } else if (token == "LM") {
474 LineShape::from_linemixingdata(icecream, line_mixing_model);
475 icecream >> token;
476 lmd_found = true;
477 } else if (token == "LF") {
479 data.lineshapetype,
480 data.selfbroadening,
481 data.bathbroadening,
482 data.line.lineshape,
483 data.species);
484 // if (data.selfbroadening) data.species[0] = data.quantumidentity.Species();
485 icecream >> token;
486 } else if (token == "ZM") {
487 // Zeeman effect
488 icecream >> data.line.zeeman;
489 icecream >> token;
490 } else if (token == "LSM") {
491 // Line shape modifications
492
493 // Starts with the number of modifications
494 icecream >> nelem;
495 for (Index lsm = 0; lsm < nelem; lsm++) {
496 icecream >> token;
497
498 // cutoff frequency
499 if (token == "CUT") {
500 icecream >> data.cutofffreq;
501 data.cutoff = CutoffType::ByLine;
502 }
503
504 // linemixing pressure limit
505 if (token == "LML") {
506 icecream >> data.linemixinglimit;
507 }
508
509 // mirroring
510 else if (token == "MTM") {
511 icecream >> data.mirroring;
512 }
513
514 // line normalization
515 else if (token == "LNT") {
516 icecream >> data.normalization;
517 }
518
519 else {
520 ARTS_USER_ERROR ("Unknown line modifications given: ", token)
521 }
522 }
523 icecream >> token;
524 } else {
525 ARTS_USER_ERROR ("Unknown line data tag in legacy reading routine: ", token)
526 }
527 }
528 }
529 } catch (const std::runtime_error& e) {
531 "Parse error in catalog line: \n",
532 line, '\n', e.what())
533 }
534
535 if (lmd_found)
536 data.line.lineshape.SetLineMixingModel(line_mixing_model.Data()[0]);
537
538 // That's it!
539 data.bad = false;
540 return data;
541}
542
544 // Default data and values for this type
546 data.selfbroadening = true;
547 data.bathbroadening = true;
548 data.lineshapetype = LineShape::Type::VP;
549 data.species.resize(2);
550 data.species[1] = Species::Species::Bath;
551
552 // This contains the rest of the line to parse. At the beginning the
553 // entire line. Line gets shorter and shorter as we continue to
554 // extract stuff from the beginning.
555 String line;
556
557 // The first item is the molecule number:
558 Index mo;
559
560 // Look for more comments?
561 bool comment = true;
562
563 while (comment) {
564 // Return true if eof is reached:
565 if (is.eof()) return data;
566
567 // Throw runtime_error if stream is bad:
568 ARTS_USER_ERROR_IF (!is, "Stream bad.");
569
570 // Read line from file into linebuffer:
571 getline(is, line);
572
573 // It is possible that we were exactly at the end of the file before
574 // calling getline. In that case the previous eof() was still false
575 // because eof() evaluates only to true if one tries to read after the
576 // end of the file. The following check catches this.
577 if (line.nelem() == 0 && is.eof()) return data;
578
579 // If the catalogue is in dos encoding, throw away the
580 // additional carriage return
581 if (line[line.nelem() - 1] == 13) {
582 line.erase(line.nelem() - 1, 1);
583 }
584
585 mo = 0;
586 // Initialization of mo is important, because mo stays the same
587 // if line is empty.
588 extract(mo, line, 2);
589 comment = false;
590 }
591
592 // Extract isotopologue:
593 char iso;
594 extract(iso, line, 1);
595
596 // Set line data
598 data.species[0] = data.quantumidentity.Species();
599
600 // Position.
601 {
602 // HITRAN position in wavenumbers (cm^-1):
603 Numeric v;
604 // Conversion from wavenumber to Hz. If you multiply a line
605 // position in wavenumber (cm^-1) by this constant, you get the
606 // frequency in Hz.
607 constexpr Numeric w2Hz = Constant::c * 100.;
608
609 // Extract HITRAN postion:
610 extract(v, line, 12);
611
612 // ARTS position in Hz:
613 data.line.F0 = v * w2Hz;
614 }
615
616 // Intensity.
617 {
618 // HITRAN intensity is in cm-1/(molec * cm-2) at 296 Kelvin.
619 // It already includes the isotpic ratio.
620 // The first cm-1 is the frequency unit (it cancels with the
621 // 1/frequency unit of the line shape function).
622 //
623 // We need to do the following:
624 // 1. Convert frequency from wavenumber to Hz (factor 1e2 * c).
625 // 2. Convert [molec * cm-2] to [molec * m-2] (factor 1e-4).
626 // 3. Take out the isotopologue ratio.
627
628 constexpr Numeric hi2arts = 1e-2 * Constant::c;
629
630 Numeric s;
631
632 // Extract HITRAN intensity:
633 extract(s, line, 10);
634 // Convert to ARTS units (Hz / (molec * m-2) ), or shorter: Hz*m^2
635 data.line.I0 = s * hi2arts;
636 // Take out isotopologue ratio:
637 data.line.I0 /= Hitran::ratio_from_lookup(mo, iso);
638 }
639
640 // Einstein coefficient
641 {
642 extract(data.line.A, line, 10);
643 }
644
645 // Air broadening parameters.
646 Numeric agam, sgam;
647 {
648 // HITRAN parameter is in cm-1/atm at 296 Kelvin
649 // All parameters are HWHM (I hope this is true!)
650 Numeric gam;
651 // Conversion from wavenumber to Hz. If you multiply a value in
652 // wavenumber (cm^-1) by this constant, you get the value in Hz.
653 constexpr Numeric w2Hz = Constant::c * 1e2;
654 // Ok, put together the end-to-end conversion that we need:
655 constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
656
657 // Extract HITRAN AGAM value:
658 extract(gam, line, 5);
659
660 // ARTS parameter in Hz/Pa:
661 agam = gam * hi2arts;
662
663 // Extract HITRAN SGAM value:
664 extract(gam, line, 5);
665
666 // ARTS parameter in Hz/Pa:
667 sgam = gam * hi2arts;
668
669 // If zero, set to agam:
670 if (0 == sgam) sgam = agam;
671
672 // cout << "agam, sgam = " << magam << ", " << msgam << endl;
673 }
674
675 // Lower state energy.
676 {
677 // HITRAN parameter is in wavenumbers (cm^-1).
678 // We have to convert this to the ARTS unit Joule.
679
680 // Extract from Catalogue line
681 extract(data.line.E0, line, 10);
682
683 // Convert to Joule:
684 data.line.E0 = wavenumber_to_joule(data.line.E0);
685 }
686
687 // Temperature coefficient of broadening parameters.
688 Numeric nair, nself;
689 {
690 // This is dimensionless, we can also extract directly.
691 extract(nair, line, 4);
692
693 // Set self broadening temperature coefficient to the same value:
694 nself = nair;
695 // cout << "mnair = " << mnair << endl;
696 }
697
698 // Pressure shift.
699 Numeric psf;
700 {
701 // HITRAN value in cm^-1 / atm. So the conversion goes exactly as
702 // for the broadening parameters.
703 Numeric d;
704 // Conversion from wavenumber to Hz. If you multiply a value in
705 // wavenumber (cm^-1) by this constant, you get the value in Hz.
706 constexpr Numeric w2Hz = Constant::c * 1e2;
707 // Ok, put together the end-to-end conversion that we need:
708 constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
709
710 // Extract HITRAN value:
711 extract(d, line, 8);
712
713 // ARTS value in Hz/Pa
714 psf = d * hi2arts;
715 }
716 // Set the accuracies using the definition of HITRAN
717 // indices. If some are missing, they are set to -1.
718
719 // Upper state global quanta
720 {
721 Index eu;
722 extract(eu, line, 15);
723 }
724
725 // Lower state global quanta
726 {
727 Index el;
728 extract(el, line, 15);
729 }
730
731 // Upper state local quanta
732 {
733 Index eul;
734 extract(eul, line, 15);
735 }
736
737 // Lower state local quanta
738 {
739 Index ell;
740 extract(ell, line, 15);
741 }
742
743 // Accuracy index for frequency
744 {
745 Index df;
746 // Extract HITRAN value:
747 extract(df, line, 1);
748 }
749
750 // Accuracy index for intensity
751 {
752 Index di0;
753 // Extract HITRAN value:
754 extract(di0, line, 1);
755 }
756
757 // Accuracy index for air-broadened halfwidth
758 {
759 Index dgam;
760 // Extract HITRAN value:
761 extract(dgam, line, 1);
762 }
763
764 // Accuracy index for self-broadened half-width
765 {
766 Index dgam;
767 // Extract HITRAN value:
768 extract(dgam, line, 1);
769 }
770
771 // Accuracy index for temperature-dependence exponent for agam
772 {
773 Index dn;
774 // Extract HITRAN value:
775 extract(dn, line, 1);
776 }
777
778 // Accuracy index for pressure shift
779 {
780 Index dpsfi;
781 // Extract HITRAN value (given in cm-1):
782 extract(dpsfi, line, 1);
783 }
784
785 // These were all the parameters that we can extract from
786 // HITRAN 2004. However, we still have to set the reference temperatures
787 // to the appropriate value:
788
789 // Reference temperature for Intensity in K.
790 data.T0 = 296.0;
791
792 // Set line shape computer
793 data.line.lineshape = LineShape::hitran_model(sgam, nself, agam, nair, psf);
794 {
795 Index garbage;
796 extract(garbage, line, 13);
797
798 // The statistical weights
799 extract(data.line.gupp, line, 7);
800 extract(data.line.glow, line, 7);
801 }
802
803 // That's it!
804 data.bad = false;
805 return data;
806}
807
809 // Default data and values for this type
811 data.selfbroadening = true;
812 data.bathbroadening = true;
813 data.lineshapetype = LineShape::Type::VP;
814 data.species.resize(2);
815 data.species[1] = Species::Species::Bath;
816
817 // This contains the rest of the line to parse. At the beginning the
818 // entire line. Line gets shorter and shorter as we continue to
819 // extract stuff from the beginning.
820 String line;
821
822 // The first item is the molecule number:
823 Index mo;
824
825 // Look for more comments?
826 bool comment = true;
827
828 while (comment) {
829 // Return true if eof is reached:
830 if (is.eof()) return data;
831
832 // Throw runtime_error if stream is bad:
833 ARTS_USER_ERROR_IF (!is, "Stream bad.");
834
835 // Read line from file into linebuffer:
836 getline(is, line);
837
838 // It is possible that we were exactly at the end of the file before
839 // calling getline. In that case the previous eof() was still false
840 // because eof() evaluates only to true if one tries to read after the
841 // end of the file. The following check catches this.
842 if (line.nelem() == 0 && is.eof()) return data;
843
844 // If the catalogue is in dos encoding, throw away the
845 // additional carriage return
846 if (line[line.nelem() - 1] == 13) {
847 line.erase(line.nelem() - 1, 1);
848 }
849
850 mo = 0;
851 // Initialization of mo is important, because mo stays the same
852 // if line is empty.
853 extract(mo, line, 2);
854 comment = false;
855 }
856
857 // Extract isotopologue:
858 char iso;
859 extract(iso, line, 1);
860
861 // Set line data
863 data.species[0] = data.quantumidentity.Species();
864
865 // Position.
866 {
867 // HITRAN position in wavenumbers (cm^-1):
868 Numeric v;
869 // Conversion from wavenumber to Hz. If you multiply a line
870 // position in wavenumber (cm^-1) by this constant, you get the
871 // frequency in Hz.
872 constexpr Numeric w2Hz = Constant::c * 100.;
873
874 // Extract HITRAN postion:
875 extract(v, line, 12);
876
877 // ARTS position in Hz:
878 data.line.F0 = v * w2Hz;
879 }
880
881 // Intensity.
882 {
883 // HITRAN intensity is in cm-1/(molec * cm-2) at 296 Kelvin.
884 // It already includes the isotpic ratio.
885 // The first cm-1 is the frequency unit (it cancels with the
886 // 1/frequency unit of the line shape function).
887 //
888 // We need to do the following:
889 // 1. Convert frequency from wavenumber to Hz (factor 1e2 * c).
890 // 2. Convert [molec * cm-2] to [molec * m-2] (factor 1e-4).
891 // 3. Take out the isotopologue ratio.
892
893 constexpr Numeric hi2arts = 1e-2 * Constant::c;
894
895 Numeric s;
896
897 // Extract HITRAN intensity:
898 extract(s, line, 10);
899 // Convert to ARTS units (Hz / (molec * m-2) ), or shorter: Hz*m^2
900 data.line.I0 = s * hi2arts;
901 // Take out isotopologue ratio:
902 data.line.I0 /= Hitran::ratio_from_lookup(mo, iso);
903 }
904
905 // Einstein coefficient
906 {
907 extract(data.line.A, line, 10);
908 }
909
910 // Air broadening parameters.
911 Numeric agam, sgam;
912 {
913 // HITRAN parameter is in cm-1/atm at 296 Kelvin
914 // All parameters are HWHM (I hope this is true!)
915 Numeric gam;
916 // Conversion from wavenumber to Hz. If you multiply a value in
917 // wavenumber (cm^-1) by this constant, you get the value in Hz.
918 constexpr Numeric w2Hz = Constant::c * 1e2;
919 // Ok, put together the end-to-end conversion that we need:
920 constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
921
922 // Extract HITRAN AGAM value:
923 extract(gam, line, 5);
924
925 // ARTS parameter in Hz/Pa:
926 agam = gam * hi2arts;
927
928 // Extract HITRAN SGAM value:
929 extract(gam, line, 5);
930
931 // ARTS parameter in Hz/Pa:
932 sgam = gam * hi2arts;
933
934 // If zero, set to agam:
935 if (0 == sgam) sgam = agam;
936
937 // cout << "agam, sgam = " << magam << ", " << msgam << endl;
938 }
939
940 // Lower state energy.
941 {
942 // HITRAN parameter is in wavenumbers (cm^-1).
943 // We have to convert this to the ARTS unit Joule.
944
945 // Extract from Catalogue line
946 extract(data.line.E0, line, 10);
947
948 // Convert to Joule:
949 data.line.E0 = wavenumber_to_joule(data.line.E0);
950 }
951
952 // Temperature coefficient of broadening parameters.
953 Numeric nair, nself;
954 {
955 // This is dimensionless, we can also extract directly.
956 extract(nair, line, 4);
957
958 // Set self broadening temperature coefficient to the same value:
959 nself = nair;
960 // cout << "mnair = " << mnair << endl;
961 }
962
963 // Pressure shift.
964 Numeric psf;
965 {
966 // HITRAN value in cm^-1 / atm. So the conversion goes exactly as
967 // for the broadening parameters.
968 Numeric d;
969 // Conversion from wavenumber to Hz. If you multiply a value in
970 // wavenumber (cm^-1) by this constant, you get the value in Hz.
971 constexpr Numeric w2Hz = Constant::c * 1e2;
972 // Ok, put together the end-to-end conversion that we need:
973 constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
974
975 // Extract HITRAN value:
976 extract(d, line, 8);
977
978 // ARTS value in Hz/Pa
979 psf = d * hi2arts;
980 }
981 // Set the accuracies using the definition of HITRAN
982 // indices. If some are missing, they are set to -1.
983
984 // Upper state global quanta
985 {
986 Index eu;
987 extract(eu, line, 15);
988 }
989
990 // Lower state global quanta
991 {
992 Index el;
993 extract(el, line, 15);
994 }
995
996 // Upper state local quanta
997 {
998 Index eul;
999 extract(eul, line, 15);
1000 }
1001
1002 // Lower state local quanta
1003 {
1004 Index ell;
1005 extract(ell, line, 15);
1006 }
1007
1008 // Accuracy index for frequency
1009 {
1010 Index df;
1011 // Extract HITRAN value:
1012 extract(df, line, 1);
1013 }
1014
1015 // Accuracy index for intensity
1016 {
1017 Index di0;
1018 // Extract HITRAN value:
1019 extract(di0, line, 1);
1020 }
1021
1022 // Accuracy index for air-broadened halfwidth
1023 {
1024 Index dgam;
1025 // Extract HITRAN value:
1026 extract(dgam, line, 1);
1027 }
1028
1029 // Accuracy index for self-broadened half-width
1030 {
1031 Index dgam;
1032 // Extract HITRAN value:
1033 extract(dgam, line, 1);
1034 }
1035
1036 // Accuracy index for temperature-dependence exponent for agam
1037 {
1038 Index dn;
1039 // Extract HITRAN value:
1040 extract(dn, line, 1);
1041 }
1042
1043 // Accuracy index for pressure shift
1044 {
1045 Index dpsfi;
1046 // Extract HITRAN value (given in cm-1):
1047 extract(dpsfi, line, 1);
1048 }
1049
1050 // These were all the parameters that we can extract from
1051 // HITRAN 2004. However, we still have to set the reference temperatures
1052 // to the appropriate value:
1053
1054 // Reference temperature for Intensity in K.
1055 data.T0 = 296.0;
1056
1057 // Set line shape computer
1058 data.line.lineshape = LineShape::hitran_model(sgam, nself, agam, nair, psf);
1059 {
1060 Index garbage;
1061 extract(garbage, line, 13);
1062
1063 // The statistical weights
1064 extract(data.line.gupp, line, 7);
1065 extract(data.line.glow, line, 7);
1066 }
1067
1068 // ADD QUANTUM NUMBER PARSING HERE!
1069 String upper, lower;
1070 std::stringstream ss;
1071 ss.str(line);
1072 ss >> upper >> lower;
1074
1075 // That's it!
1076 data.bad = false;
1077 return data;
1078}
1079
1081 // Default data and values for this type
1082 SingleLineExternal data;
1083 data.selfbroadening = true;
1084 data.bathbroadening = true;
1085 data.lineshapetype = LineShape::Type::VP;
1086 data.species.resize(2);
1087 data.species[1] = Species::Species::Bath;
1088
1089 // This contains the rest of the line to parse. At the beginning the
1090 // entire line. Line gets shorter and shorter as we continue to
1091 // extract stuff from the beginning.
1092 String line;
1093
1094 // The first item is the molecule number:
1095 Index mo;
1096
1097 // Look for more comments?
1098 bool comment = true;
1099
1100 while (comment) {
1101 // Return true if eof is reached:
1102 if (is.eof()) return data;
1103
1104 // Throw runtime_error if stream is bad:
1105 ARTS_USER_ERROR_IF (!is, "Stream bad.");
1106
1107 // Read line from file into linebuffer:
1108 getline(is, line);
1109
1110 // It is possible that we were exactly at the end of the file before
1111 // calling getline. In that case the previous eof() was still false
1112 // because eof() evaluates only to true if one tries to read after the
1113 // end of the file. The following check catches this.
1114 if (line.nelem() == 0 && is.eof()) return data;
1115
1116 // If the catalogue is in dos encoding, throw away the
1117 // additional carriage return
1118 if (line[line.nelem() - 1] == 13) {
1119 line.erase(line.nelem() - 1, 1);
1120 }
1121
1122 mo = 0;
1123 // Initialization of mo is important, because mo stays the same
1124 // if line is empty.
1125 extract(mo, line, 2);
1126 comment = false;
1127 }
1128
1129 // Extract isotopologue:
1130 char iso;
1131 extract(iso, line, 1);
1132
1133 // Set line data
1135 data.species[0] = data.quantumidentity.Species();
1136
1137 // Position.
1138 {
1139 // HITRAN position in wavenumbers (cm^-1):
1140 Numeric v;
1141 // Conversion from wavenumber to Hz. If you multiply a line
1142 // position in wavenumber (cm^-1) by this constant, you get the
1143 // frequency in Hz.
1144 constexpr Numeric w2Hz = Constant::c * 100.;
1145
1146 // Extract HITRAN postion:
1147 extract(v, line, 12);
1148
1149 // ARTS position in Hz:
1150 data.line.F0 = v * w2Hz;
1151 // cout << "mf = " << mf << endl;
1152 }
1153
1154 // Intensity.
1155 {
1156 // HITRAN intensity is in cm-1/(molec * cm-2) at 296 Kelvin.
1157 // It already includes the isotpic ratio.
1158 // The first cm-1 is the frequency unit (it cancels with the
1159 // 1/frequency unit of the line shape function).
1160 //
1161 // We need to do the following:
1162 // 1. Convert frequency from wavenumber to Hz (factor 1e2 * c).
1163 // 2. Convert [molec * cm-2] to [molec * m-2] (factor 1e-4).
1164 // 3. Take out the isotopologue ratio.
1165
1166 constexpr Numeric hi2arts = 1e-2 * Constant::c;
1167
1168 Numeric s;
1169
1170 // Extract HITRAN intensity:
1171 extract(s, line, 10);
1172 // Convert to ARTS units (Hz / (molec * m-2) ), or shorter: Hz*m^2
1173 data.line.I0 = s * hi2arts;
1174 // Take out isotopologue ratio:
1175 data.line.I0 /= Hitran::ratio_from_lookup(mo, iso);
1176 }
1177
1178 // Skip transition probability:
1179 {
1180 Numeric r;
1181 extract(r, line, 10);
1182 }
1183
1184 // Air broadening parameters.
1185 Numeric agam, sgam;
1186 {
1187 // HITRAN parameter is in cm-1/atm at 296 Kelvin
1188 // All parameters are HWHM (I hope this is true!)
1189 Numeric gam;
1190 // Conversion from wavenumber to Hz. If you multiply a value in
1191 // wavenumber (cm^-1) by this constant, you get the value in Hz.
1192 constexpr Numeric w2Hz = Constant::c * 1e2;
1193 // Ok, put together the end-to-end conversion that we need:
1194 constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
1195
1196 // Extract HITRAN AGAM value:
1197 extract(gam, line, 5);
1198
1199 // ARTS parameter in Hz/Pa:
1200 agam = gam * hi2arts;
1201
1202 // Extract HITRAN SGAM value:
1203 extract(gam, line, 5);
1204
1205 // ARTS parameter in Hz/Pa:
1206 sgam = gam * hi2arts;
1207
1208 // If zero, set to agam:
1209 if (0 == sgam) sgam = agam;
1210
1211 // cout << "agam, sgam = " << magam << ", " << msgam << endl;
1212 }
1213
1214 // Lower state energy.
1215 {
1216 // HITRAN parameter is in wavenumbers (cm^-1).
1217 // We have to convert this to the ARTS unit Joule.
1218
1219 // Extract from Catalogue line
1220 extract(data.line.E0, line, 10);
1221
1222 // Convert to Joule:
1223 data.line.E0 = wavenumber_to_joule(data.line.E0);
1224 }
1225
1226 // Temperature coefficient of broadening parameters.
1227 Numeric nair, nself;
1228 {
1229 // This is dimensionless, we can also extract directly.
1230 extract(nair, line, 4);
1231
1232 // Set self broadening temperature coefficient to the same value:
1233 nself = nair;
1234 // cout << "mnair = " << mnair << endl;
1235 }
1236
1237 // Pressure shift.
1238 Numeric psf;
1239 {
1240 // HITRAN value in cm^-1 / atm. So the conversion goes exactly as
1241 // for the broadening parameters.
1242 Numeric d;
1243 // Conversion from wavenumber to Hz. If you multiply a value in
1244 // wavenumber (cm^-1) by this constant, you get the value in Hz.
1245 constexpr Numeric w2Hz = Constant::c * 1e2;
1246 // Ok, put together the end-to-end conversion that we need:
1247 constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
1248
1249 // Extract HITRAN value:
1250 extract(d, line, 8);
1251
1252 // ARTS value in Hz/Pa
1253 psf = d * hi2arts;
1254 }
1255 // Set the accuracies using the definition of HITRAN
1256 // indices. If some are missing, they are set to -1.
1257
1258 //Skip upper state global quanta index
1259 {
1260 Index eu;
1261 extract(eu, line, 3);
1262 }
1263
1264 //Skip lower state global quanta index
1265 {
1266 Index el;
1267 extract(el, line, 3);
1268 }
1269
1270 //Skip upper state local quanta
1271 {
1272 Index eul;
1273 extract(eul, line, 9);
1274 }
1275
1276 //Skip lower state local quanta
1277 {
1278 Index ell;
1279 extract(ell, line, 9);
1280 }
1281
1282 // Accuracy index for frequency reference
1283 {
1284 Index df;
1285 // Extract HITRAN value:
1286 extract(df, line, 1);
1287 }
1288
1289 // Accuracy index for intensity reference
1290 {
1291 Index di0;
1292 // Extract HITRAN value:
1293 extract(di0, line, 1);
1294 }
1295
1296 // Accuracy index for halfwidth reference
1297 {
1298 Index dgam;
1299 // Extract HITRAN value:
1300 extract(dgam, line, 1);
1301 }
1302
1303 // These were all the parameters that we can extract from
1304 // HITRAN. However, we still have to set the reference temperatures
1305 // to the appropriate value:
1306
1307 // Reference temperature for Intensity in K.
1308 data.T0 = 296.0;
1309
1310 // Set line shape computer
1311 data.line.lineshape = LineShape::hitran_model(sgam, nself, agam, nair, psf);
1312
1313 // That's it!
1314 data.bad = false;
1315 return data;
1316}
1317
1319 // Default data and values for this type
1320 SingleLineExternal data;
1321 data.selfbroadening = true;
1322 data.bathbroadening = true;
1323 data.lineshapetype = LineShape::Type::VP;
1324 data.species.resize(2);
1325
1326 // This contains the rest of the line to parse. At the beginning the
1327 // entire line. Line gets shorter and shorter as we continue to
1328 // extract stuff from the beginning.
1329 String line;
1330
1331 // The first item is the molecule number:
1332 Index mo;
1333
1334 // Look for more comments?
1335 bool comment = true;
1336
1337 while (comment) {
1338 // Return true if eof is reached:
1339 if (is.eof()) return data;
1340
1341 // Throw runtime_error if stream is bad:
1342 ARTS_USER_ERROR_IF (!is, "Stream bad.");
1343
1344 // Read line from file into linebuffer:
1345 getline(is, line);
1346 if (line[0] == '>' or line[0] == '%') continue;
1347
1348 // It is possible that we were exactly at the end of the file before
1349 // calling getline. In that case the previous eof() was still false
1350 // because eof() evaluates only to true if one tries to read after the
1351 // end of the file. The following check catches this.
1352 if (line.nelem() == 0 && is.eof()) return data;
1353
1354 // If the catalogue is in dos encoding, throw away the
1355 // additional carriage return
1356 if (line[line.nelem() - 1] == 13) {
1357 line.erase(line.nelem() - 1, 1);
1358 }
1359
1360 mo = 0;
1361 // Initialization of mo is important, because mo stays the same
1362 // if line is empty.
1363 extract(mo, line, 2);
1364 comment = false;
1365 }
1366
1367 // Extract isotopologue:
1368 char iso;
1369 extract(iso, line, 1);
1370
1371 // Set line data
1373
1374 // Position.
1375 {
1376 // HITRAN position in wavenumbers (cm^-1):
1377 Numeric v;
1378 // Conversion from wavenumber to Hz. If you multiply a line
1379 // position in wavenumber (cm^-1) by this constant, you get the
1380 // frequency in Hz.
1381 constexpr Numeric w2Hz = Constant::c * 100.;
1382
1383 // Extract HITRAN postion:
1384 extract(v, line, 12);
1385
1386 // ARTS position in Hz:
1387 data.line.F0 = v * w2Hz;
1388 // cout << "mf = " << mf << endl;
1389 }
1390
1391 // Intensity.
1392 {
1393 // HITRAN intensity is in cm-1/(molec * cm-2) at 296 Kelvin.
1394 // It already includes the isotpic ratio.
1395 // The first cm-1 is the frequency unit (it cancels with the
1396 // 1/frequency unit of the line shape function).
1397 //
1398 // We need to do the following:
1399 // 1. Convert frequency from wavenumber to Hz (factor 1e2 * c).
1400 // 2. Convert [molec * cm-2] to [molec * m-2] (factor 1e-4).
1401 // 3. Take out the isotopologue ratio.
1402
1403 constexpr Numeric hi2arts = 1e-2 * Constant::c;
1404
1405 Numeric s;
1406 if (line[6] == 'D') line[6] = 'E';
1407 // Extract HITRAN intensity:
1408 extract(s,
1409 line,
1410 10); // NOTE: If error shooting, FORTRAN "D" is not read properly.
1411 // Convert to ARTS units (Hz / (molec * m-2) ), or shorter: Hz*m^2
1412 data.line.I0 = s * hi2arts;
1413 // Take out isotopologue ratio:
1414 data.line.I0 /= Hitran::ratio_from_lookup(mo, iso);
1415 }
1416
1417 // Skip transition probability:
1418 {
1419 Numeric r;
1420 extract(r, line, 10);
1421 }
1422
1423 // Air broadening parameters.
1424 Numeric sgam, agam;
1425 {
1426 // HITRAN parameter is in cm-1/atm at 296 Kelvin
1427 // All parameters are HWHM (I hope this is true!)
1428 Numeric gam;
1429 // Conversion from wavenumber to Hz. If you multiply a value in
1430 // wavenumber (cm^-1) by this constant, you get the value in Hz.
1431 constexpr Numeric w2Hz = Constant::c * 1e2;
1432 // Ok, put together the end-to-end conversion that we need:
1433 constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
1434
1435 // Extract HITRAN AGAM value:
1436 extract(gam, line, 5);
1437
1438 // ARTS parameter in Hz/Pa:
1439 agam = gam * hi2arts;
1440
1441 // Extract HITRAN SGAM value:
1442 extract(gam, line, 5);
1443
1444 // ARTS parameter in Hz/Pa:
1445 sgam = gam * hi2arts;
1446
1447 // If zero, set to agam:
1448 if (0 == sgam) sgam = agam;
1449
1450 // cout << "agam, sgam = " << magam << ", " << msgam << endl;
1451 }
1452
1453 // Lower state energy.
1454 {
1455 // HITRAN parameter is in wavenumbers (cm^-1).
1456 // We have to convert this to the ARTS unit Joule.
1457
1458 // Extract from Catalogue line
1459 extract(data.line.E0, line, 10);
1460
1461 // Convert to Joule:
1462 data.line.E0 = wavenumber_to_joule(data.line.E0);
1463 }
1464
1465 // Temperature coefficient of broadening parameters.
1466 Numeric nair, nself;
1467 {
1468 // This is dimensionless, we can also extract directly.
1469 extract(nair, line, 4);
1470
1471 // Set self broadening temperature coefficient to the same value:
1472 nself = nair;
1473 // cout << "mnair = " << mnair << endl;
1474 }
1475
1476 // Pressure shift.
1477 Numeric psf;
1478 {
1479 // HITRAN value in cm^-1 / atm. So the conversion goes exactly as
1480 // for the broadening parameters.
1481 Numeric d;
1482 // Conversion from wavenumber to Hz. If you multiply a value in
1483 // wavenumber (cm^-1) by this constant, you get the value in Hz.
1484 constexpr Numeric w2Hz = Constant::c * 1e2;
1485 // Ok, put together the end-to-end conversion that we need:
1486 constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
1487
1488 // Extract HITRAN value:
1489 extract(d, line, 8);
1490
1491 // ARTS value in Hz/Pa
1492 psf = d * hi2arts;
1493 }
1494 // Set the accuracies using the definition of HITRAN
1495 // indices. If some are missing, they are set to -1.
1496
1497 //Skip upper state global quanta index
1498 {
1499 Index eu;
1500 extract(eu, line, 3);
1501 }
1502
1503 //Skip lower state global quanta index
1504 {
1505 Index el;
1506 extract(el, line, 3);
1507 }
1508
1509 //Skip upper state local quanta
1510 {
1511 Index eul;
1512 extract(eul, line, 9);
1513 }
1514
1515 //Skip lower state local quanta
1516 {
1517 Index ell;
1518 extract(ell, line, 9);
1519 }
1520
1521 // Accuracy index for frequency reference
1522 {
1523 Index df;
1524 // Extract HITRAN value:
1525 extract(df, line, 1);
1526 }
1527
1528 // Accuracy index for intensity reference
1529 {
1530 Index di0;
1531 // Extract HITRAN value:
1532 extract(di0, line, 1);
1533 }
1534
1535 // Accuracy index for halfwidth reference
1536 {
1537 Index dgam;
1538 // Extract HITRAN value:
1539 extract(dgam, line, 1);
1540 }
1541
1542 // These were all the parameters that we can extract from
1543 // HITRAN. However, we still have to set the reference temperatures
1544 // to the appropriate value:
1545
1546 // Reference temperature for Intensity in K.
1547 // (This is fix for HITRAN)
1548 data.T0 = 296.0;
1549
1550 // Skip four
1551 {
1552 Index four;
1553 extract(four, line, 4);
1554 }
1555
1556 // This is the test for the last two characters of the line
1557 {
1558 /*
1559 * 0 is nothing,
1560 * -1 is linemixing on the next line,
1561 * -3 is the non-resonant line
1562 */
1563 Index test;
1564 extract(test, line, 2);
1565 //If the tag is as it should be, then a minus one means that more should be read
1566 if (test == -1 || test == -3)
1567 getline(is, line);
1568 else // the line is done and we are happy to leave
1569 {
1570 data.line.lineshape = LineShape::hitran_model(sgam, nself, agam, nair, psf);
1571
1572 data.bad = false;
1573 return data;
1574 }
1575 }
1576
1577 // In case we are unable to leave, the next line is a line mixing parameter line
1578
1579 // First is the molecular number. This should be the same as above.
1580 {
1581 Index mo2;
1582 extract(mo2, line, 2);
1583 // Skip one
1584
1585 ARTS_USER_ERROR_IF (mo != mo2,
1586 "There is an error in the line mixing\n");
1587 }
1588
1589 Vector Y(4), G(4), T(4);
1590
1591 // These are constants for AER but should be included because we need their grid.
1592 T[0] = 200;
1593 T[1] = 250;
1594 T[2] = 296;
1595 T[3] = 340;
1596
1597 // Next is the Y and G at various temperatures
1598 {
1599 Numeric Y_200K;
1600 extract(Y_200K, line, 13);
1601 Y[0] = Y_200K;
1602 }
1603 {
1604 Numeric G_200K;
1605 extract(G_200K, line, 11);
1606 G[0] = G_200K;
1607 }
1608 {
1609 Numeric Y_250K;
1610 extract(Y_250K, line, 13);
1611 Y[1] = Y_250K;
1612 }
1613 {
1614 Numeric G_250K;
1615 extract(G_250K, line, 11);
1616 G[1] = G_250K;
1617 }
1618 {
1619 Numeric Y_296K;
1620 extract(Y_296K, line, 13);
1621 Y[2] = Y_296K;
1622 }
1623 {
1624 Numeric G_296K;
1625 extract(G_296K, line, 11);
1626 G[2] = G_296K;
1627 }
1628 {
1629 Numeric Y_340K;
1630 extract(Y_340K, line, 13);
1631 Y[3] = Y_340K;
1632 }
1633 {
1634 Numeric G_340K;
1635 extract(G_340K, line, 11);
1636 G[3] = G_340K;
1637 }
1638
1639 // Cpnvert from per Atm and per Atm^2
1640 Y /= Conversion::atm2pa(1);
1642
1643 // ARTS uses (1-iY) as line-mixing factor, LBLRTM uses (1+iY), so we must change sign
1644 Y *= -1;
1645
1646 // Test that this is the end
1647 {
1648 Index test;
1649 extract(test, line, 2);
1650 if (test == -1) {
1652 nself,
1653 agam,
1654 nair,
1655 psf,
1656 {T[0],
1657 T[1],
1658 T[2],
1659 T[3],
1660 Y[0],
1661 Y[1],
1662 Y[2],
1663 Y[3],
1664 G[0],
1665 G[1],
1666 G[2],
1667 G[3]});
1668
1669 data.bad = false;
1670 return data;
1671 }
1672 if (test == -3) {
1674 nself,
1675 agam,
1676 nair,
1677 psf,
1678 {T[0],
1679 T[1],
1680 T[2],
1681 T[3],
1682 Y[0],
1683 Y[1],
1684 Y[2],
1685 Y[3],
1686 G[0],
1687 G[1],
1688 G[2],
1689 G[3]});
1690
1691 data.bad = false;
1692 return data;
1693 }
1694 return data;
1695 }
1696}
1697
1698std::vector<Absorption::Lines> Absorption::split_list_of_external_lines(std::vector<SingleLineExternal>& external_lines,
1699 const std::vector<QuantumNumberType>& localquantas,
1700 const std::vector<QuantumNumberType>& globalquantas)
1701{
1702 std::vector<Lines> lines(0);
1703
1704 // Loop but make copies because we will need to modify some of the data
1705 while(external_lines.size()) {
1706 auto& sle = external_lines.back();
1707
1708 // Adapt broadening to fit with line catalog
1709 if (sle.selfbroadening) sle.species.front() = sle.quantumidentity.Species();
1710 if (sle.bathbroadening) sle.species.back() = Species::Species::Bath;
1711
1712 Quantum::Number::GlobalState global_id{sle.quantumidentity.isotopologue_index};
1713 for (auto qn: globalquantas) if (sle.quantumidentity.val.has(qn)) global_id.val.set(sle.quantumidentity.val[qn]);
1714
1715 Quantum::Number::LocalState local_id{};
1716 for (auto qn: localquantas) if (sle.quantumidentity.val.has(qn)) local_id.val.set(sle.quantumidentity.val[qn]);
1717
1718 // Set the line
1719 auto line = sle.line;
1720 line.localquanta = local_id;
1721
1722 // Either find a line like this in the list of lines or start a new Lines
1723 auto band = std::find_if(lines.begin(), lines.end(), [&](const Lines& li){return li.MatchWithExternal(sle, global_id);});
1724 if (band not_eq lines.end()) {
1725 band -> AppendSingleLine(line);
1726 } else {
1727 lines.push_back(Lines(sle.selfbroadening, sle.bathbroadening, sle.cutoff,
1728 sle.mirroring, sle.population, sle.normalization,
1729 sle.lineshapetype, sle.T0, sle.cutofffreq,
1730 sle.linemixinglimit, global_id, sle.species, {line}));
1731 }
1732 external_lines.pop_back();
1733 }
1734
1735 return lines;
1736}
1737
1738std::ostream& Absorption::operator<<(std::ostream& os, const Absorption::Lines& lines)
1739{
1740 for(auto& line: lines.lines)
1741 os << line << '\n';
1742 return os;
1743}
1744
1745std::istream& Absorption::operator>>(std::istream& is, Lines& lines) {
1746 for(auto& line: lines.lines)
1747 is >> line;
1748 return is;
1749}
1750
1751std::ostream& Absorption::operator<<(std::ostream& os,
1752 const Absorption::SingleLine& line) {
1753 os << line.F0 << ' ' << line.I0 << ' ' << line.E0 << ' ' << line.glow << ' '
1754 << line.gupp << ' ' << line.A << ' ' << line.zeeman << ' '
1755 << line.lineshape;
1756 if (line.localquanta.val.nelem()) os << ' ' << line.localquanta.values();
1757 return os;
1758}
1759
1760std::istream& Absorption::operator>>(std::istream& is, Absorption::SingleLine& line)
1761{
1762 is >> double_imanip() >> line.F0 >> line.I0 >> line.E0 >> line.glow >>
1763 line.gupp >> line.A;
1764 return is >> line.zeeman >> line.lineshape >> line.localquanta;
1765}
1766
1767std::ostream& operator<<(std::ostream& os, const ArrayOfAbsorptionLines& aol)
1768{
1769 for(auto& l: aol)
1770 os << l << '\n';
1771 return os;
1772}
1773
1774std::ostream& operator<<(std::ostream& os, const ArrayOfArrayOfAbsorptionLines& aol)
1775{
1776 for(auto& l: aol)
1777 os << l << '\n';
1778 return os;
1779}
1780
1782{
1784}
1785
1787{
1788 std::ostringstream os;
1789
1790 os << "\nLines meta-data:\n";
1791 os << '\t' << "Species identity:\n";
1792 os << "\t\tSpecies: "<< SpeciesName() << '\n';
1793 os << "\t\tIdentity: "<< quantumidentity << '\n';
1794 os << '\t' << cutofftype2metadatastring(cutoff, cutofffreq);
1795 os << '\t' << populationtype2metadatastring(population);
1796 os << '\t' << normalizationtype2metadatastring(normalization);
1797 os << '\t' << LineShape::shapetype2metadatastring(lineshapetype);
1798 os << '\t' << mirroringtype2metadatastring(mirroring);
1799 os << '\t' << "The reference temperature for all line parameters is "
1800 << T0 << " K.\n";
1801 if(linemixinglimit < 0)
1802 os << '\t' << "If applicable, there is no line mixing limit.\n";
1803 else
1804 os << '\t' << "If applicable, there is a line mixing limit at "
1805 << linemixinglimit << " Pa.\n";
1806
1807 if (not NumLines()) {
1808 os << "\tNo line data is available.\n";
1809 } else {
1810 os << "\tThere are " << NumLines() << " lines available.\n";
1811
1812 auto& line = lines.front();
1813 os << "\tThe front line has:\n";
1814 os << "\t\t" << "f0: " << line.F0 << " Hz\n";
1815 os << "\t\t" << "i0: " << line.I0 << " m^2/Hz\n";
1816 os << "\t\t" << "e0: " << line.E0 << " J\n";
1817 os << "\t\t" << "Lower stat. weight: " << line.glow << " [-]\n";
1818 os << "\t\t" << "Upper stat. weight: " << line.gupp << " [-]\n";
1819 os << "\t\t" << "A: " << line.A << " 1/s\n";
1820 os << "\t\t" << "Zeeman splitting of lower state: " << line.zeeman.gl() << " [-]\n";
1821 os << "\t\t" << "Zeeman splitting of upper state: " << line.zeeman.gu() << " [-]\n";
1822 os << "\t\t" << "Local quantum numbers: " << line.localquanta.val << "\n";
1823
1824 ArrayOfString ls_meta = LineShape::ModelMetaDataArray(line.lineshape,
1825 selfbroadening,
1826 broadeningspecies,
1827 T0);
1828 os << "\t\t" << "Line shape parameters (are normalized by sum(VMR)):\n";
1829 for(auto& ls_form: ls_meta)
1830 os << "\t\t\t" << ls_form << "\n";
1831 }
1832
1833
1834 return os.str();
1835}
1836
1837
1839{
1840 lines.erase(lines.begin() + i);
1841}
1842
1843
1845{
1846 auto line = lines[i];
1847 RemoveLine(i);
1848 return line;
1849}
1850
1852 std::reverse(lines.begin(), lines.end());
1853}
1854
1856 return quantumidentity.Isotopologue().mass;
1857}
1858
1860 const ArrayOfArrayOfSpeciesTag& atm_spec) const
1861{
1862 if (lineshapetype == LineShape::Type::DP)
1863 return Vector(broadeningspecies.nelem(), std::numeric_limits<Numeric>::quiet_NaN());
1864 return LineShape::vmrs(atm_vmrs, atm_spec, broadeningspecies);
1865}
1866
1868 const ArrayOfArrayOfSpeciesTag& atm_spec,
1869 const SpeciesIsotopologueRatios& ir,
1870 const Numeric& bath_mass
1871 ) const
1872{
1873 if (lineshapetype == LineShape::Type::DP)
1874 return Vector(broadeningspecies.nelem(), std::numeric_limits<Numeric>::quiet_NaN());
1875 Vector mass = LineShape::mass(atm_vmrs, atm_spec, broadeningspecies, ir);
1876 if (bathbroadening and bath_mass > 0) mass[mass.nelem()-1] = bath_mass;
1877 return mass;
1878}
1879
1881 const ArrayOfArrayOfSpeciesTag& atm_spec) const
1882{
1883 ARTS_USER_ERROR_IF (atm_vmrs.nelem() not_eq atm_spec.nelem(),
1884 "Bad species and vmr lists");
1885
1886 for (Index i=0; i<atm_spec.nelem(); i++)
1887 if (atm_spec[i].nelem() and atm_spec[i].Species() == Species())
1888 return atm_vmrs[i];
1889 return 0;
1890}
1891
1893 if (not iseven(Jf + lf + 1))
1894 return - sqrt(2 * Jf + 1) * wigner3j(Jf, k, Ji, li, lf - li, -lf);
1895 return + sqrt(2 * Jf + 1) * wigner3j(Jf, k, Ji, li, lf - li, -lf);
1896}
1897
1899 if (not iseven(Jf + N))
1900 return - sqrt(6 * (2 * Jf + 1) * (2 * Ji + 1)) * wigner6j(1, 1, 1, Ji, Jf, N);
1901 return + sqrt(6 * (2 * Jf + 1) * (2 * Ji + 1)) * wigner6j(1, 1, 1, Ji, Jf, N);
1902}
1903
1905{
1906 // Default data and values for this type
1907 SingleLineExternal data;
1908 data.selfbroadening = true;
1909 data.bathbroadening = true;
1910 data.lineshapetype = LineShape::Type::VP;
1911 data.species.resize(2);
1912
1913 // This contains the rest of the line to parse. At the beginning the
1914 // entire line. Line gets shorter and shorter as we continue to
1915 // extract stuff from the beginning.
1916 String line;
1917
1918 // Look for more comments?
1919 bool comment = true;
1920
1921 while (comment) {
1922 // Return true if eof is reached:
1923 if (is.eof()) return data;
1924
1925 // Throw runtime_error if stream is bad:
1926 ARTS_USER_ERROR_IF (!is, "Stream bad.");
1927
1928 // Read line from file into linebuffer:
1929 getline(is, line);
1930
1931 // It is possible that we were exactly at the end of the file before
1932 // calling getline. In that case the previous eof() was still false
1933 // because eof() evaluates only to true if one tries to read after the
1934 // end of the file. The following check catches this.
1935 if (line.nelem() == 0 && is.eof()) return data;
1936
1937 // Because of the fixed FORTRAN format, we need to break up the line
1938 // explicitly in apropriate pieces. Not elegant, but works!
1939
1940 // Extract center frequency:
1941 // Initialization of v is important, because v stays the same
1942 // if line is empty.
1943 // JPL position in MHz:
1944 Numeric v = 0.0;
1945
1946 // Extract JPL position:
1947 extract(v, line, 13);
1948
1949 // check for empty line
1950 if (v != 0.0) {
1951 // ARTS position in Hz:
1952 data.line.F0 = v * 1E6;
1953
1954 comment = false;
1955 }
1956 }
1957
1958 // Accuracy for line position
1959 {
1960 Numeric df;
1961 extract(df, line, 8);
1962 }
1963
1964 // Intensity.
1965 {
1966 // JPL has log (10) of intensity in nm2 MHz at 300 Kelvin.
1967 //
1968 // We need to do the following:
1969 // 1. take 10^intensity
1970 // 2. convert to cm-1/(molecule * cm-2): devide by c * 1e10
1971 // 3. Convert frequency from wavenumber to Hz (factor 1e2 * c)
1972 // 4. Convert [molec * cm-2] to [molec * m-2] (factor 1e-4)
1973
1974 Numeric s;
1975
1976 // Extract JPL intensity:
1977 extract(s, line, 8);
1978
1979 // remove log
1980 s = pow((Numeric)10., s);
1981
1982 // Convert to ARTS units (Hz / (molec * m-2) ), or shorter: Hz*m^2
1983 data.line.I0 = s / 1E12;
1984 }
1985
1986 // Degrees of freedom
1987 {
1988 Index dr;
1989
1990 // Extract degrees of freedom
1991 extract(dr, line, 2);
1992 }
1993
1994 // Lower state energy.
1995 {
1996 // JPL parameter is in wavenumbers (cm^-1).
1997 // We have to convert this to the ARTS unit Joule.
1998
1999 // Extract from Catalogue line
2000 extract(data.line.E0, line, 10);
2001
2002 // Convert to Joule:
2003 data.line.E0 = wavenumber_to_joule(data.line.E0);
2004 }
2005
2006 // Upper state degeneracy
2007 {
2008 Index gup;
2009
2010 // Extract upper state degeneracy
2011 extract(gup, line, 3);
2012 }
2013
2014 // Tag number
2015 Index tag;
2016 {
2017 // Extract Tag number
2018 extract(tag, line, 7);
2019
2020 // make sure tag is not negative (damned jpl cat):
2021 tag = tag > 0 ? tag : -tag;
2022 }
2023
2024 // Set line ID
2026
2027 // Air broadening parameters: unknown to jpl, use old iup forward
2028 // model default values, which is mostly set to 0.0025 GHz/hPa, even
2029 // though for some lines the pressure broadening is given explicitly
2030 // in the program code. The explicitly given values are ignored and
2031 // only the default value is set. Self broadening was in general not
2032 // considered in the old forward model.
2033 Numeric agam, sgam;
2034 {
2035 // ARTS parameter in Hz/Pa:
2036 agam = 2.5E4;
2037
2038 // ARTS parameter in Hz/Pa:
2039 sgam = agam;
2040 }
2041
2042 // Temperature coefficient of broadening parameters. Was set to 0.75
2043 // in old forward model, even though for some lines the parameter is
2044 // given explicitly in the program code. The explicitly given values
2045 // are ignored and only the default value is set. Self broadening
2046 // not considered.
2047 Numeric nair, nself;
2048 {
2049 nair = 0.75;
2050 nself = 0.0;
2051 }
2052
2053 // Reference temperature for broadening parameter in K, was
2054 // generally set to 300 K in old forward model, with the exceptions
2055 // as already mentioned above: //DEPRECEATED but is same as for mti0 so moving on
2056 // {
2057 // mtgam = 300.0;
2058 // }
2059
2060 // Pressure shift: not given in JPL, set to 0
2061 Numeric psf;
2062 { psf = 0.0; }
2063
2064 // These were all the parameters that we can extract from
2065 // JPL. However, we still have to set the reference temperatures
2066 // to the appropriate value:
2067
2068 // Reference temperature for Intensity in K.
2069 data.T0 = 300.0;
2070
2071 // Set line shape computer
2072 data.line.lineshape = LineShape::Model(sgam, nself, agam, nair, psf);
2073
2074 // That's it!
2075 data.bad = false;
2076 return data;
2077}
2078
2079bool Absorption::Lines::OK() const noexcept
2080{
2081 const Index nb = broadeningspecies.nelem();
2082
2083 // Check that the isotopologue is ok
2084 if (not Isotopologue().OK()) return false;
2085
2086 // Test that self and bath is covered by the range if set positive
2087 if (nb < (Index(selfbroadening) + Index(bathbroadening)))
2088 return false;
2089
2090 // Test that the temperature is physical
2091 if (T0 <= 0)
2092 return false;
2093
2094 // Test that all lines have the correct sized line shape model
2095 if (std::any_of(lines.cbegin(), lines.cend(), [nb](auto& line){return line.LineShapeElems() != nb;}))
2096 return false;
2097
2098 // Test that all lines have the correct sized local quantum numbers
2099 if (std::any_of(lines.cbegin(), lines.cend(), [&](auto& line){return line.LocalQuantumElems() != lines.front().LocalQuantumElems();}))
2100 return false;
2101
2102 // Otherwise everything is fine!
2103 return true;
2104}
2105
2107{
2108 return std::sqrt(Constant::doppler_broadening_const_squared * T / SpeciesMass());
2109}
2110
2111
2112
2113namespace Absorption {
2115 std::ostringstream os;
2116 switch (in) {
2117 case CutoffType::None:
2118 os << "No cut-off will be applied.\n"; break;
2119 case CutoffType::ByLine:
2120 os << "The lines will be cut-off " << cutoff << " Hz from the line center + D0.\n"; break;
2121 case CutoffType::FINAL: break;
2122 }
2123 return os.str();
2124}
2125
2127 for (auto& val: localquanta.val) qid.val.set(val); // Copy to same struct
2128
2129 zeeman = Zeeman::Model(qid);
2130}
2131
2136 .Data()[0]);
2137}
2138
2140 const LineShape::ModelParameters Y = {LineShape::TemperatureModel::LM_AER, d[4], d[5], d[6], d[7]};
2141 const LineShape::ModelParameters G = {LineShape::TemperatureModel::LM_AER, d[8], d[9], d[10], d[11]};
2142 for (auto& sm : lineshape.Data()) {
2143 sm.Y() = Y;
2144 sm.G() = G;
2145 }
2146}
2147
2150 bif >> F0 >> I0 >> E0 >> glow >> gupp >> A >> zeeman;
2151
2153 lineshape.read(bif);
2154
2156 for (auto& val: localquanta.val) val.read(bif);
2157
2158 return bif;
2159}
2160
2163 bof << F0 << I0 << E0 << glow << gupp << A << zeeman;
2164
2166 lineshape.write(bof);
2167
2169 for (auto& val: localquanta.val) val.write(bof);
2170
2171 return bof;
2172}
2173
2176 NumLines() not_eq 0 and NumLocalQuanta() not_eq sl.LocalQuantumElems(),
2177 "Error calling appending function, bad size of quantum numbers");
2178
2179 ARTS_USER_ERROR_IF(NumLines() not_eq 0 and
2180 sl.LineShapeElems() not_eq lines[0].LineShapeElems(),
2181 "Error calling appending function, bad size of broadening species");
2182
2183 lines.push_back(std::move(sl));
2184}
2185
2188 NumLines() not_eq 0 and NumLocalQuanta() not_eq sl.LocalQuantumElems(),
2189 "Error calling appending function, bad size of quantum numbers\n"
2190 "Type of quantum numbers in band: ", lines.front().localquanta.val,
2191 "\nType of quantum numbers in new line:", sl.localquanta.val
2192 );
2193
2195 NumLines() not_eq 0 and
2196 sl.LineShapeElems() not_eq lines.front().LineShapeElems(),
2197 "Error calling appending function, bad size of broadening species");
2198
2199 lines.push_back(sl);
2200}
2201
2203 if(sle.bad)
2204 return false;
2205 if(sle.selfbroadening not_eq selfbroadening)
2206 return false;
2207 if(sle.bathbroadening not_eq bathbroadening)
2208 return false;
2209 if(sle.cutoff not_eq cutoff)
2210 return false;
2211 if(sle.mirroring not_eq mirroring)
2212 return false;
2213 if(sle.population not_eq population)
2214 return false;
2215 if(sle.normalization not_eq normalization)
2216 return false;
2217 if(sle.lineshapetype not_eq lineshapetype)
2218 return false;
2219 if(sle.T0 not_eq T0)
2220 return false;
2221 if(sle.cutofffreq not_eq cutofffreq)
2222 return false;
2223 if(sle.linemixinglimit not_eq linemixinglimit)
2224 return false;
2225 if(quantumidentity not_eq qid)
2226 return false;
2227 if(not std::equal(sle.species.cbegin(), sle.species.cend(), broadeningspecies.cbegin(), broadeningspecies.cend()))
2228 return false;
2229 if(NumLines() not_eq 0 and not lines.front().localquanta.same_types_as(sle.line.localquanta))
2230 return false;
2231 if(NumLines() not_eq 0 and not sle.line.lineshape.Match(lines.front().lineshape).first)
2232 return false;
2233 return true;
2234}
2235
2236std::pair<bool, bool> Lines::Match(const Lines& l) const noexcept {
2237 // Note: The pair here is first: matching and second: nullable
2238 if(l.selfbroadening not_eq selfbroadening)
2239 return {false, false};
2240 if(l.bathbroadening not_eq bathbroadening)
2241 return {false, false};
2242 if(l.cutoff not_eq cutoff)
2243 return {false, false};
2244 if(l.mirroring not_eq mirroring)
2245 return {false, false};
2246 if(l.population not_eq population)
2247 return {false, false};
2248 if(l.normalization not_eq normalization)
2249 return {false, false};
2250 if(l.lineshapetype not_eq lineshapetype)
2251 return {false, false};
2252 if(l.T0 not_eq T0)
2253 return {false, false};
2254 if(l.cutofffreq not_eq cutofffreq)
2255 return {false, false};
2256 if(l.linemixinglimit not_eq linemixinglimit)
2257 return {false, false};
2258 if(l.quantumidentity not_eq quantumidentity)
2259 return {false, false};
2260 if(not std::equal(l.broadeningspecies.cbegin(), l.broadeningspecies.cend(), broadeningspecies.cbegin(), broadeningspecies.cend()))
2261 return {false, false};
2262 if(NumLines() not_eq 0 and l.NumLines() not_eq 0 and not lines.front().localquanta.same_types_as(l.lines.front().localquanta))
2263 return {false, false};
2264 if(NumLines() not_eq 0 and l.NumLines() not_eq 0) {
2265 if (auto matchpair = l.lines.front().lineshape.Match(lines.front().lineshape); not matchpair.first) return matchpair;
2266 }
2267
2268 return {true, true};
2269}
2270
2272 std::sort(lines.begin(), lines.end(),
2273 [](const SingleLine& a, const SingleLine& b){return a.F0 < b.F0;});
2274}
2275
2277 std::sort(lines.begin(), lines.end(),
2278 [](const SingleLine& a, const SingleLine& b){return a.A < b.A;});
2279}
2280
2282 return NumLines() ?
2284 "";
2285}
2286
2288 ARTS_ASSERT(qns.val.has(QuantumNumberType::F) or qns.val.has(QuantumNumberType::J))
2289 return qns.val.has(QuantumNumberType::F) ? qns.val[QuantumNumberType::F] : qns.val[QuantumNumberType::J];
2290}
2291
2292Index Lines::ZeemanCount(size_t k, Zeeman::Polarization type) const noexcept {
2293 if (type == Zeeman::Polarization::None) return 1;
2294
2295 // Select F before J but assume one of them exist
2296 auto& val = get(lines[k].localquanta);
2297 return Zeeman::nelem(val.upp(), val.low(), type);
2298}
2299
2302 Index i) const noexcept {
2303 if (type == Zeeman::Polarization::None) return 1.0;
2304
2305 // Select F before J but assume one of them exist
2306 auto& val = get(lines[k].localquanta);
2307 return lines[k].zeeman.Strength(val.upp(), val.low(), type, i);
2308}
2309
2312 Index i) const noexcept {
2313 if (type == Zeeman::Polarization::None) return 0.0;
2314
2315 // Select F before J but assume one of them exist
2316 auto& val = get(lines[k].localquanta);
2317 return lines[k].zeeman.Splitting(val.upp(), val.low(), type, i);
2318}
2319
2321 for(auto& line: lines)
2322 line.SetAutomaticZeeman(quantumidentity);
2323}
2324
2325Numeric Lines::F_mean(const ConstVectorView& wgts) const noexcept {
2326 const Numeric val = std::inner_product(lines.cbegin(), lines.cend(),
2327 wgts.begin(), 0.0, std::plus<>(),
2328 [](const auto& a, const auto& b){return a.F0 * b;});
2329 const Numeric div = wgts.sum();
2330 return val / div;
2331}
2332
2333Numeric Lines::F_mean(Numeric T) const noexcept {
2334 if (T <= 0) T = T0;
2335
2336 const Index n = NumLines();
2337 const Numeric QT = single_partition_function(T, Isotopologue());
2338 const Numeric QT0 = single_partition_function(T0, Isotopologue());
2339 const Numeric ratiopart = QT0 / QT;
2340
2341 Vector wgts(n);
2342 for (Index i=0; i<n; i++) {
2343 const Numeric pop0 = (lines[i].gupp / QT0) * boltzman_factor(T0, lines[i].E0);
2344 const Numeric pop = pop0 * ratiopart * boltzman_ratio(T, T0, lines[i].E0);
2345 const Numeric dip_squared = - lines[i].I0/(pop0 * lines[i].F0 * std::expm1(- (Constant::h * lines[i].F0) / (Constant::k * T0)));
2346 wgts[i] = pop * dip_squared;
2347 }
2348
2349 return F_mean(wgts);
2350}
2351
2352#pragma GCC diagnostic push
2353#pragma GCC diagnostic ignored "-Wreturn-type"
2354Numeric Lines::CutoffFreq(size_t k, Numeric shift) const noexcept {
2355 switch(cutoff) {
2356 case CutoffType::ByLine:
2357 return lines[k].F0 + cutofffreq + (mirroring == MirroringType::Manual ? - shift : shift);
2358 case CutoffType::None:
2360 case CutoffType::FINAL: break;
2361 }
2362}
2363#pragma GCC diagnostic pop
2364
2365#pragma GCC diagnostic push
2366#pragma GCC diagnostic ignored "-Wreturn-type"
2367Numeric Lines::CutoffFreqMinus(size_t k, Numeric shift) const noexcept {
2368 switch(cutoff) {
2369 case CutoffType::ByLine:
2370 return lines[k].F0 - cutofffreq + (mirroring == MirroringType::Manual ? - shift : shift);
2371 case CutoffType::None:
2372 return std::numeric_limits<Numeric>::lowest();
2373 case CutoffType::FINAL: break;
2374 }
2375}
2376#pragma GCC diagnostic pop
2377
2379 for (auto& line: lines)
2380 line.read(is);
2381 return is;
2382}
2383
2385 for (auto& line: lines)
2386 line.write(os);
2387 return os;
2388}
2389
2391 QuantumIdentifier qid = quantumidentity;
2392 for (auto& a: lines[k].localquanta.val) qid.val.set(a);
2393 return qid;
2394}
2395
2397 const Index n = NumLines();
2398 const Index m = NumBroadeners();
2399 if (not n) return;
2400
2401 for (Index j=0; j<m; j++) {
2402 for (Index i=0; i<LineShape::nVars; i++) {
2403
2404 // Find a common type (same or none) or throw if there is none
2405 LineShape::TemperatureModel t=LineShape::TemperatureModel::None;
2406 for (auto& line : lines) {
2407 if (auto& data = line.lineshape[j].Data()[i]; not LineShape::modelparameterEmpty(data)) {
2408 if (t == LineShape::TemperatureModel::None) t = data.type;
2409 ARTS_USER_ERROR_IF(t not_eq data.type, "Cannot make a common line shape model for the band as there are multiple non-empty types: ", data.type, " and ", t)
2410 }
2411 }
2412
2413 // Set the common type
2414 for (auto& line : lines) {
2415 if (auto& data = line.lineshape[j].Data()[i]; data.type not_eq t) {
2417 }
2418 }
2419 }
2420 }
2421}
2422} // namespace Absorption
Numeric wavenumber_to_joule(Numeric e)
A little helper function to convert energy from units of wavenumber (cm^-1) to Joule (J).
Definition: absorption.cc:85
Declarations required for the calculation of absorption coefficients.
std::ostream & operator<<(std::ostream &os, const ArrayOfAbsorptionLines &aol)
Contains the absorption namespace.
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
A constant view of a Vector.
Definition: matpackI.h:517
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:541
Main line shape model class.
void SetLineMixingModel(SingleSpeciesModel x)
Sets the same line mixing model to all species.
bifstream & read(bifstream &bif)
Binary read for Model.
bofstream & write(bofstream &bof) const
Binary write for Model.
const std::vector< SingleSpeciesModel > & Data() const noexcept
The line shape model data.
void set(Value v)
Sets the value if it exists or adds it otherwise.
Index nelem() const
Return number of quantum numbers.
bool has(Types... ts) const ARTS_NOEXCEPT
Returns whether all the Types are part of the list, the types must be sorted.
Value & add(Type t)
Add for manipulation.
Implements rational numbers to work with other ARTS types.
Definition: rational.h:52
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition: jacobian.h:330
The Vector class.
Definition: matpackI.h:908
Main Zeeman Model.
Definition: zeemandata.h:296
Binary output file stream class.
Definition: bifstream.h:42
Binary output file stream class.
Definition: bofstream.h:42
Input manipulator class for doubles to enable nan and inf parsing.
Definition: double_imanip.h:64
Index nelem() const
Number of elements.
Definition: mystring.h:253
Constants of physical expressions as constexpr.
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
This file contains basic functions to handle ASCII files.
#define max(a, b)
Numeric boltzman_factor(Numeric T, Numeric E0)
Computes exp(- E0/kT)
Definition: linescaling.cc:114
Numeric boltzman_ratio(const Numeric &T, const Numeric &T0, const Numeric &E0)
Computes exp(E0/c (T - T0) / (T * T0))
Definition: linescaling.cc:97
Numeric single_partition_function(const Numeric &T, const Species::IsotopeRecord &ir)
Computes the partition function at one temperature.
Definition: linescaling.cc:23
Constains various line scaling functions.
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
void extract(T &x, String &line, Index n)
Extract something from the beginning of a string.
Definition: mystring.h:304
Namespace to contain things required for absorption calculations.
std::ostream & operator<<(std::ostream &, const SingleLine &)
SingleLineExternal ReadFromHitranOnlineStream(istream &is)
Read from HITRAN online.
String cutofftype2metadatastring(CutoffType in, Numeric cutoff)
SingleLineExternal ReadFromJplStream(istream &is)
Read from JPL.
const Quantum::Number::Value & get(const Quantum::Number::LocalState &qns) ARTS_NOEXCEPT
Numeric reduced_rovibrational_dipole(Rational Jf, Rational Ji, Rational lf, Rational li, Rational k=Rational(1))
Compute the reduced rovibrational dipole moment.
SingleLineExternal ReadFromArtscat3Stream(istream &is)
Read from ARTSCAT-3.
Index nelem(const Lines &l)
Number of lines.
SingleLineExternal ReadFromLBLRTMStream(istream &is)
Read from LBLRTM.
Numeric cutoff
std::istream & operator>>(std::istream &, SingleLine &)
SingleLineExternal ReadFromArtscat5Stream(istream &is)
Read from ARTSCAT-5.
Numeric reduced_magnetic_quadrapole(Rational Jf, Rational Ji, Rational N)
Compute the reduced magnetic quadrapole moment.
SingleLineExternal ReadFromHitran2004Stream(istream &is)
Read from newer HITRAN.
SingleLineExternal ReadFromHitran2001Stream(istream &is)
Read from HITRAN before 2004.
SingleLineExternal ReadFromArtscat4Stream(istream &is)
Read from ARTSCAT-4.
std::vector< Lines > split_list_of_external_lines(std::vector< SingleLineExternal > &external_lines, const std::vector< QuantumNumberType > &localquantas={}, const std::vector< QuantumNumberType > &globalquantas={})
Splits a list of lines into proper Lines.
constexpr auto pow2(T x) noexcept -> decltype(x *x)
power of two
Definition: constants.h:66
constexpr auto atm2pa(T x) noexcept -> decltype(x *101 '325.0)
Conversion from Atm to Pa.
Definition: constants.h:446
QuantumIdentifier id_from_lookup(Index mol, char isochar)
Finds the ID of the ARTS species from HITRAN.
Numeric ratio_from_lookup(Index mol, char isochar)
Finds the isotopologue ratio of the species from HITRAN.
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
QuantumIdentifier id_from_lookup(Index tag)
Finds the ID of the ARTS species from JPL.
Definition: jpl_species.cc:204
Model vector2modellm(Vector x, LegacyLineMixingData::TypeLM type)
LineShape::Model from legacy input vector.
ArrayOfString ModelMetaDataArray(const Model &m, const bool self, const ArrayOfSpecies &sts, const Numeric T0)
Model lblrtm_model(Numeric sgam, Numeric nself, Numeric agam, Numeric nair, Numeric psf, std::array< Numeric, 12 > aer_interp)
constexpr ModelParameters modelparameterGetEmpty(const TemperatureModel t) noexcept
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.
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::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.
ValueList from_hitran(std::string_view upp, std::string_view low)
String update_isot_name(const String &old_name)
Updates the name of the isotopologue based on updates of the isotopologues.
constexpr Index find_species_index(const Species spec, const std::string_view isot) noexcept
constexpr Index nelem(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the number of elements of the polarization type of this transition.
Definition: zeemandata.h:142
Polarization
Zeeman polarization selection.
Definition: zeemandata.h:44
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
#define d
#define v
#define a
#define c
#define b
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:747
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
constexpr bool iseven(const Rational r) noexcept
Returns true if even integer.
Definition: rational.h:942
#define N
Definition: rng.cc:164
void sort_by_einstein()
Sort inner line list by Einstein coefficient.
void sort_by_frequency()
Sort inner line list by frequency.
Numeric ZeemanSplitting(size_t k, Zeeman::Polarization type, Index i) const noexcept
Returns the splitting of a Zeeman split line.
SingleLine PopLine(Index) noexcept
Pops a single line.
Index NumLocalQuanta() const noexcept
Number of broadening species.
void RemoveLine(Index) noexcept
Removes a single line.
Numeric F_mean(Numeric T=0) const noexcept
Mean frequency by weight of line strength.
void SetAutomaticZeeman() noexcept
Set Zeeman effect for all lines that have the correct quantum numbers.
Vector BroadeningSpeciesVMR(const ConstVectorView &, const ArrayOfArrayOfSpeciesTag &) const
Returns the VMRs of the broadening species.
LineShape::Output ShapeParameters_dT(size_t k, Numeric T, Numeric P, const Vector &vmrs) const noexcept
Line shape parameters temperature derivatives.
LineShape::Output ShapeParameters_dVMR(size_t k, Numeric T, Numeric P, const QuantumIdentifier &vmr_qid) const noexcept
Line shape parameters vmr derivative.
Array< SingleLine > lines
A list of individual lines.
bool OK() const noexcept
Numeric ShapeParameter_dInternal(size_t k, Numeric T, Numeric P, const Vector &vmrs, const RetrievalQuantity &derivative) const noexcept
Line shape parameter internal derivative.
Index NumLines() const noexcept
Number of lines.
bofstream & write(bofstream &os) const
Binary write for Lines.
Index NumBroadeners() const noexcept
Number of broadening species.
Numeric DopplerConstant(Numeric T) const noexcept
QuantumIdentifier QuantumIdentityOfLine(Index k) const noexcept
bifstream & read(bifstream &is)
Binary read for Lines.
LineShape::Output ShapeParameters(size_t k, Numeric T, Numeric P, const Vector &vmrs) const noexcept
Line shape parameters.
bool MatchWithExternal(const SingleLineExternal &sle, const QuantumIdentifier &quantumidentity) const ARTS_NOEXCEPT
Checks if an external line matches this structure.
String LineShapeMetaData() const noexcept
Meta data for the line shape if it exists.
void MakeLineShapeModelCommon()
Make a common line shape if possible.
String SpeciesName() const noexcept
Species Name.
Numeric CutoffFreq(size_t k, Numeric shift=0) const noexcept
Returns cutoff frequency or maximum value.
void ReverseLines() noexcept
Reverses the order of the internal lines.
std::pair< bool, bool > Match(const Lines &l) const noexcept
Checks if another line list matches this structure.
void AppendSingleLine(SingleLine &&sl)
Appends a single line to the absorption lines.
Numeric SelfVMR(const ConstVectorView &, const ArrayOfArrayOfSpeciesTag &) const
Returns the VMR of the species.
Numeric CutoffFreqMinus(size_t k, Numeric shift=0) const noexcept
Returns negative cutoff frequency or lowest value.
String MetaData() const
Returns a printable statement about the lines.
Index ZeemanCount(size_t k, Zeeman::Polarization type) const noexcept
Returns the number of Zeeman split lines.
Numeric SpeciesMass() const noexcept
Mass of the molecule.
Index LineShapePos(const Species::Species spec) const noexcept
Position among broadening species or -1.
QuantumIdentifier quantumidentity
Catalog ID.
Numeric ZeemanStrength(size_t k, Zeeman::Polarization type, Index i) const noexcept
Returns the strength of a Zeeman split line.
Vector BroadeningSpeciesMass(const ConstVectorView &, const ArrayOfArrayOfSpeciesTag &, const SpeciesIsotopologueRatios &, const Numeric &bath_mass=0) const
Returns the mass of the broadening species.
Single line reading output.
Computations and data for a single absorption line.
Numeric E0
Lower state energy level.
bofstream & write(bofstream &bof) const
Binary write for AbsorptionLines.
bifstream & read(bifstream &bif)
Binary read for AbsorptionLines.
Numeric F0
Central frequency.
Quantum::Number::LocalState localquanta
Local quantum numbers.
LineShape::Model lineshape
Line shape model.
void SetAutomaticZeeman(QuantumIdentifier qid)
Set Zeeman effect by automatic detection.
Zeeman::Model zeeman
Zeeman model.
Numeric A
Einstein spontaneous emission coefficient.
Index LocalQuantumElems() const noexcept
Number of lower quantum numbers.
Numeric glow
Lower level statistical weight.
void SetLineMixing2AER(const Vector &d)
Set the line mixing model to AER kind.
Numeric gupp
Upper level statistical weight.
Index LineShapeElems() const noexcept
Number of lineshape elements.
void SetLineMixing2SecondOrderData(const Vector &d)
Set the line mixing model to 2nd order.
Numeric I0
Reference intensity.
Coefficients and temperature model for SingleSpeciesModel.
A logical struct for global quantum numbers with species identifiers.
Species::Species Species() const noexcept
Species::IsotopeRecord Isotopologue() const noexcept
A logical struct for local quantum numbers.
A complete quantum number value with type information.
constexpr void set(std::string_view s, bool upp)
Set level value.
String FullName() const noexcept
Definition: isotopologues.h:49
Numeric wigner6j(const Rational j1, const Rational j2, const Rational j3, const Rational l1, const Rational l2, const Rational l3)
Wigner 6J symbol.
Numeric wigner3j(const Rational j1, const Rational j2, const Rational j3, const Rational m1, const Rational m2, const Rational m3)
Wigner 3J symbol.