absorptionlines.cc

Go to the documentation of this file.

/* Copyright (C) 2019
 Richard Larsson <ric.larsson@gmail.com>
 
 
 This program is free software; you can redistribute it and/or modify it
 under the terms of the GNU General Public License as published by the
 Free Software Foundation; either version 2, or (at your option) any
 later version.
 
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.
 
 
 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
 USA. */
 
#include "absorptionlines.h"
 
#include "absorption.h"
#include "constants.h"
#include "debug.h"
#include "enums.h"
#include "file.h"
#include "hitran_species.h"
#include "jpl_species.h"
#include "linescaling.h"
#include "quantum_numbers.h"
 
LineShape::Output Absorption::Lines::ShapeParameters(size_t k, Numeric T, Numeric P, const Vector& vmrs) const noexcept {
  auto x = lines[k].lineshape.GetParams(T, T0, P, vmrs);
  
  if (not DoLineMixing(P)) x.Y = x.G = x.DV = 0;
  
  return x;
}
 
LineShape::Output Absorption::Lines::ShapeParameters(size_t k, Numeric T, Numeric P, size_t m) const noexcept {
  return lines[k].lineshape.GetParams(T, T0, P, m);
}
 
LineShape::Output Absorption::Lines::ShapeParameters_dT(size_t k, Numeric T, Numeric P, const Vector& vmrs) const noexcept {
  auto x = lines[k].lineshape.GetTemperatureDerivs(T, T0, P, vmrs);
  
  if (not DoLineMixing(P)) x.Y = x.G = x.DV = 0;
  
  return x;
}
 
Index Absorption::Lines::LineShapePos(const Species::Species spec) const noexcept {
  // Is always first if this is self and self broadening exists
  if(selfbroadening and spec == quantumidentity.Species()) {
    return 0;
  }
  
  // First and last might be artificial so they should not be checked
  const Index s = selfbroadening;
  const Index e = broadeningspecies.nelem() - bathbroadening;
  for(Index i=s; i<e; i++) {
    if(spec == broadeningspecies[i]) {
      return i;
    }
  }
  
  // At this point, the ID is not explicitly among the broadeners, but bath broadening means its VMR still might matter
  if(bathbroadening)
    return broadeningspecies.nelem()-1;
  return -1;
}
 
LineShape::Output Absorption::Lines::ShapeParameters_dVMR(size_t k, Numeric T, Numeric P, const QuantumIdentifier& vmr_qid) const noexcept {
  const Index pos=LineShapePos(vmr_qid);
  if (pos >= 0)
    return lines[k].lineshape.GetVMRDerivs(T, T0, P, pos);
  return LineShape::Output{};
}
 
Numeric Absorption::Lines::ShapeParameter_dInternal(size_t k, Numeric T, Numeric P, const Vector& vmrs, const RetrievalQuantity& derivative) const noexcept {
  const auto self = derivative.Mode() == LineShape::self_broadening;
  const auto bath = derivative.Mode() == LineShape::bath_broadening;
  const auto& ls = lines[k].lineshape;
  
  if(derivative.QuantumIdentity().Species() not_eq Species() or
      derivative.QuantumIdentity().Isotopologue() not_eq Isotopologue())
    return 0;
  if(self and selfbroadening)
    return ls.GetInternalDeriv(
      T, T0, P, 0, vmrs, derivative.LineType());
  if(self)
    return ls.GetInternalDeriv(
      T, T0, P, LineShapePos(derivative.QuantumIdentity().Species()), vmrs, derivative.LineType());
  if(bath and bathbroadening)
    return ls.GetInternalDeriv(
      T, T0, P, ls.nelem() - 1, vmrs, derivative.LineType());
  if(bath)
    return 0;
  return ls.GetInternalDeriv(
    T, T0, P, LineShapePos(derivative.QuantumIdentity().Species()), vmrs, derivative.LineType());
}
 
Absorption::SingleLineExternal Absorption::ReadFromArtscat3Stream(istream& is) {
  // Default data and values for this type
  SingleLineExternal data;
  data.selfbroadening = true;
  data.bathbroadening = true;
  data.lineshapetype = LineShape::Type::VP;
  data.species.resize(2);
 
  // This always contains the rest of the line to parse. At the
  // beginning the entire line. Line gets shorter and shorter as we
  // continue to extract stuff from the beginning.
  String line;
 
  // Look for more comments?
  bool comment = true;
 
  while (comment) {
    // Return true if eof is reached:
    if (is.eof()) return data;
 
    // Throw runtime_error if stream is bad:
    ARTS_USER_ERROR_IF (!is, "Stream bad.");
 
    // Read line from file into linebuffer:
    getline(is, line);
 
    // It is possible that we were exactly at the end of the file before
    // calling getline. In that case the previous eof() was still false
    // because eof() evaluates only to true if one tries to read after the
    // end of the file. The following check catches this.
    if (line.nelem() == 0 && is.eof()) return data;
 
    // @ as first character marks catalogue entry
    char c;
    extract(c, line, 1);
 
    // check for empty line
    if (c == '@') {
      comment = false;
    }
  }
 
  // read the arts identifier String
  istringstream icecream(line);
 
  String artsid;
  icecream >> artsid;
 
  if (artsid.length() != 0) {
    // Set the species
    const auto isotopologue = Species::Tag(Species::update_isot_name(artsid));
    ARTS_USER_ERROR_IF (isotopologue.is_joker() or isotopologue.type not_eq Species::TagType::Plain,
                        "A line catalog species can only be of the form \"Plain\", meaning it\nhas the form SPECIES-ISONUM.\n"
                        "Your input contains: ", artsid, ". which we cannot interpret as a plain species")
    data.quantumidentity = QuantumIdentifier{Species::find_species_index(isotopologue.Isotopologue())};
 
    // Extract center frequency:
    icecream >> double_imanip() >> data.line.F0;
 
    Numeric psf;
    // Extract pressure shift:
    icecream >> double_imanip() >> psf;
 
    // Extract intensity:
    icecream >> double_imanip() >> data.line.I0;
 
    // Extract reference temperature for Intensity in K:
    icecream >> double_imanip() >> data.T0;
 
    // Extract lower state energy:
    icecream >> data.line.E0;
 
    // Extract air broadening parameters:
    Numeric agam, sgam;
    icecream >> double_imanip() >> agam;
    icecream >> double_imanip() >> sgam;
 
    // Extract temperature coefficient of broadening parameters:
    Numeric nair, nself;
    icecream >> double_imanip() >> nair;
    icecream >> double_imanip() >> nself;
 
    // Extract reference temperature for broadening parameter in K:
    Numeric tgam;
    icecream >> double_imanip() >> tgam;
 
    // Extract the aux parameters:
    Index naux;
    icecream >> naux;
 
    // resize the aux array and read it
    ArrayOfNumeric maux;
    maux.resize(naux);
 
    for (Index j = 0; j < naux; j++) {
      icecream >> maux[j];
      //cout << "maux" << j << " = " << maux[j] << "\n";
    }
 
    // Extract accuracies:
    try {
      Numeric unused_numeric;
      icecream >> double_imanip() >> unused_numeric /*mdf*/;
      icecream >> double_imanip() >> unused_numeric /*mdi0*/;
      icecream >> double_imanip() >> unused_numeric /*dagam*/;
      icecream >> double_imanip() >> unused_numeric /*dsgam*/;
      icecream >> double_imanip() >> unused_numeric /*dnair*/;
      icecream >> double_imanip() >> unused_numeric /*dnself*/;
      icecream >> double_imanip() >> unused_numeric /*dpsf*/;
    } catch (const std::runtime_error&) {
    }
 
    // Fix if tgam is different from ti0
    if (tgam != data.T0) {
      agam = agam * pow(tgam / data.T0, nair);
      sgam = sgam * pow(tgam / data.T0, nself);
      psf = psf * pow(tgam / data.T0, (Numeric).25 + (Numeric)1.5 * nair);
    }
 
    // Set line shape computer
    data.line.lineshape = LineShape::Model(sgam, nself, agam, nair, psf);
  }
 
  // That's it!
  data.bad = false;
  return data;
}
 
Absorption::SingleLineExternal Absorption::ReadFromArtscat4Stream(istream& is) {
  // Default data and values for this type
  SingleLineExternal data;
  data.selfbroadening = true;
  data.bathbroadening = false;
  data.lineshapetype = LineShape::Type::VP;
 
  // This always contains the rest of the line to parse. At the
  // beginning the entire line. Line gets shorter and shorter as we
  // continue to extract stuff from the beginning.
  String line;
 
  // Look for more comments?
  bool comment = true;
 
  while (comment) {
    // Return true if eof is reached:
    if (is.eof()) return data;
 
    // Throw runtime_error if stream is bad:
    ARTS_USER_ERROR_IF (!is, "Stream bad.");
 
    // Read line from file into linebuffer:
    getline(is, line);
 
    // It is possible that we were exactly at the end of the file before
    // calling getline. In that case the previous eof() was still false
    // because eof() evaluates only to true if one tries to read after the
    // end of the file. The following check catches this.
    if (line.nelem() == 0 && is.eof()) return data;
 
    // @ as first character marks catalogue entry
    char c;
    extract(c, line, 1);
 
    // check for empty line
    if (c == '@') {
      comment = false;
    }
  }
 
  // read the arts identifier String
  istringstream icecream(line);
 
  String artsid;
  icecream >> artsid;
 
  if (artsid.length() != 0) {
 
    // Set line ID
    const auto isotopologue = Species::Tag(Species::update_isot_name(artsid));
    ARTS_USER_ERROR_IF (isotopologue.is_joker() or isotopologue.type not_eq Species::TagType::Plain,
      "A line catalog species can only be of the form \"Plain\", meaning it\nhas the form SPECIES-ISONUM.\n"
      "Your input contains: ", artsid, ". which we cannot interpret as a plain species")
    data.quantumidentity = QuantumIdentifier{Species::find_species_index(isotopologue.Isotopologue())};
 
    // Extract center frequency:
    icecream >> double_imanip() >> data.line.F0;
 
    // Extract intensity:
    icecream >> double_imanip() >> data.line.I0;
 
    // Extract reference temperature for Intensity in K:
    icecream >> double_imanip() >> data.T0;
 
    // Extract lower state energy:
    icecream >> double_imanip() >> data.line.E0;
 
    // Extract Einstein A-coefficient:
    icecream >> double_imanip() >> data.line.A;
 
    // Extract upper state stat. weight:
    icecream >> double_imanip() >> data.line.gupp;
 
    // Extract lower state stat. weight:
    icecream >> double_imanip() >> data.line.glow;
 
    LineShape::from_artscat4(icecream, 
                             data.lineshapetype,
                             data.selfbroadening,
                             data.bathbroadening,
                             data.line.lineshape,
                             data.species,
                             data.quantumidentity);
  }
 
  // That's it!
  data.bad = false;
  return data;
}
 
Absorption::SingleLineExternal Absorption::ReadFromArtscat5Stream(istream& is) {
  // Default data and values for this type
  SingleLineExternal data;
 
  LineShape::Model line_mixing_model;
  bool lmd_found = false;
 
  // This always contains the rest of the line to parse. At the
  // beginning the entire line. Line gets shorter and shorter as we
  // continue to extract stuff from the beginning.
  String line;
 
  // Look for more comments?
  bool comment = true;
  
  while (comment) {
    // Return true if eof is reached:
    if (is.eof()) return data;
 
    // Throw runtime_error if stream is bad:
    ARTS_USER_ERROR_IF (!is, "Stream bad.");
 
    // Read line from file into linebuffer:
    getline(is, line);
 
    // It is possible that we were exactly at the end of the file before
    // calling getline. In that case the previous eof() was still false
    // because eof() evaluates only to true if one tries to read after the
    // end of the file. The following check catches this.
    if (line.nelem() == 0 && is.eof()) return data;
 
    // @ as first character marks catalogue entry
    char c;
    extract(c, line, 1);
 
    // check for empty line
    if (c == '@') {
      comment = false;
    }
  }
 
  // read the arts identifier String
  istringstream icecream(line);
 
  try {
    String artsid;
    icecream >> artsid;
    
    if (artsid.length() != 0) {
      // Set line ID:
      const auto isotopologue = Species::Tag(Species::update_isot_name(artsid));
      ARTS_USER_ERROR_IF (isotopologue.is_joker() or isotopologue.type not_eq Species::TagType::Plain,
        "A line catalog species can only be of the form \"Plain\", meaning it\nhas the form SPECIES-ISONUM.\n"
        "Your input contains: ", artsid, ". which we cannot interpret as a plain species")
      data.quantumidentity = QuantumIdentifier{Species::find_species_index(isotopologue.Isotopologue())};
 
      // Extract center frequency:
      icecream >> double_imanip() >> data.line.F0;
 
      // Extract intensity:
      icecream >> double_imanip() >> data.line.I0;
 
      // Extract reference temperature for Intensity in K:
      icecream >> double_imanip() >> data.T0;
 
      // Extract lower state energy:
      icecream >> double_imanip() >> data.line.E0;
 
      // Extract Einstein A-coefficient:
      icecream >> double_imanip() >> data.line.A;
 
      // Extract upper state stat. weight:
      icecream >> double_imanip() >> data.line.gupp;
 
      // Extract lower state stat. weight:
      icecream >> double_imanip() >> data.line.glow;
 
      String token;
      Index nelem;
      icecream >> token;
 
      while (icecream) {
        // Read pressure broadening (LEGACY)
        if (token == "PB") {
          LineShape::from_pressurebroadeningdata(
              icecream,
              data.lineshapetype,
              data.selfbroadening,
              data.bathbroadening,
              data.line.lineshape,
              data.species,
              data.quantumidentity);
          icecream >> token;
        } else if (token == "QN") {
          // Quantum numbers
 
          icecream >> token;
          ARTS_USER_ERROR_IF (token != "UP",
            "Unknown quantum number tag: ", token)
 
          icecream >> token;
          String r;
          while (icecream and token not_eq "LO") {
            auto qn = Quantum::Number::toType(token);
            icecream >> r;
            icecream >> token;
 
            if (good_enum(qn)) {
              if (data.quantumidentity.val.has(qn)) {
                Quantum::Number::Value val = data.quantumidentity.val[qn];
                val.set(r, true);
                data.quantumidentity.val.set(val);
              } else {
                data.quantumidentity.val.add(qn).set(r, true);
              }
            }
          }
 
          ARTS_USER_ERROR_IF (!is || token != "LO",
            "Error in catalog. Lower quantum number tag 'LO' not found.")
 
          icecream >> token;
          auto qn = Quantum::Number::toType(token);
          while (icecream and not (token == "LM" or token == "LF" or token == "ZM" or token == "LSM" or token == "PB")) {
            icecream >> r;
            icecream >> token;
 
            if (good_enum(qn)) {
              if (data.quantumidentity.val.has(qn)) {
                Quantum::Number::Value val = data.quantumidentity.val[qn];
                val.set(r, false);
                data.quantumidentity.val.set(val);
              } else {
                data.quantumidentity.val.add(qn).set(r, false);
              }
            }
 
            qn = Quantum::Number::toType(token);
          }
        } else if (token == "LM") {
          LineShape::from_linemixingdata(icecream, line_mixing_model);
          icecream >> token;
          lmd_found = true;
        } else if (token == "LF") {
          LineShape::from_linefunctiondata(icecream,
                                           data.lineshapetype,
                                           data.selfbroadening,
                                           data.bathbroadening,
                                           data.line.lineshape,
                                           data.species);
          // if (data.selfbroadening) data.species[0] = data.quantumidentity.Species();
          icecream >> token;
        } else if (token == "ZM") {
          // Zeeman effect
          icecream >> data.line.zeeman;
          icecream >> token;
        } else if (token == "LSM") {
          // Line shape modifications
 
          // Starts with the number of modifications
          icecream >> nelem;
          for (Index lsm = 0; lsm < nelem; lsm++) {
            icecream >> token;
 
            // cutoff frequency
            if (token == "CUT") {
              icecream >> data.cutofffreq;
              data.cutoff = CutoffType::ByLine;
            }
 
            // linemixing pressure limit
            if (token == "LML") {
              icecream >> data.linemixinglimit;
            }
 
            // mirroring
            else if (token == "MTM") {
              icecream >> data.mirroring;
            }
 
            // line normalization
            else if (token == "LNT") {
              icecream >> data.normalization;
            }
 
            else {
              ARTS_USER_ERROR ("Unknown line modifications given: ", token)
            }
          }
          icecream >> token;
        } else {
          ARTS_USER_ERROR ("Unknown line data tag in legacy reading routine: ", token)
        }
      }
    }
  } catch (const std::runtime_error& e) {
    ARTS_USER_ERROR (
                        "Parse error in catalog line: \n",
                        line, '\n', e.what())
  }
 
  if (lmd_found)
    data.line.lineshape.SetLineMixingModel(line_mixing_model.Data()[0]);
 
  // That's it!
  data.bad = false;
  return data;
}
 
Absorption::SingleLineExternal Absorption::ReadFromHitran2004Stream(istream& is) {
  // Default data and values for this type
  SingleLineExternal data;
  data.selfbroadening = true;
  data.bathbroadening = true;
  data.lineshapetype = LineShape::Type::VP;
  data.species.resize(2);
  data.species[1] = Species::Species::Bath;
 
  // This contains the rest of the line to parse. At the beginning the
  // entire line. Line gets shorter and shorter as we continue to
  // extract stuff from the beginning.
  String line;
 
  // The first item is the molecule number:
  Index mo;
  
  // Look for more comments?
  bool comment = true;
  
  while (comment) {
    // Return true if eof is reached:
    if (is.eof()) return data;
    
    // Throw runtime_error if stream is bad:
    ARTS_USER_ERROR_IF (!is, "Stream bad.");
    
    // Read line from file into linebuffer:
    getline(is, line);
    
    // It is possible that we were exactly at the end of the file before
    // calling getline. In that case the previous eof() was still false
    // because eof() evaluates only to true if one tries to read after the
    // end of the file. The following check catches this.
    if (line.nelem() == 0 && is.eof()) return data;
    
    // If the catalogue is in dos encoding, throw away the
    // additional carriage return
    if (line[line.nelem() - 1] == 13) {
      line.erase(line.nelem() - 1, 1);
    }
    
    mo = 0;
    // Initialization of mo is important, because mo stays the same
    // if line is empty.
    extract(mo, line, 2);
    comment = false;
  }
  
  // Extract isotopologue:
  char iso;
  extract(iso, line, 1);
  
  // Set line data
  data.quantumidentity = Hitran::id_from_lookup(mo, iso);
  data.species[0] = data.quantumidentity.Species();
 
  // Position.
  {
    // HITRAN position in wavenumbers (cm^-1):
    Numeric v;
    // Conversion from wavenumber to Hz. If you multiply a line
    // position in wavenumber (cm^-1) by this constant, you get the
    // frequency in Hz.
    constexpr Numeric w2Hz = Constant::c * 100.;
 
    // Extract HITRAN postion:
    extract(v, line, 12);
 
    // ARTS position in Hz:
    data.line.F0 = v * w2Hz;
  }
 
  // Intensity.
  {
    // HITRAN intensity is in cm-1/(molec * cm-2) at 296 Kelvin.
    // It already includes the isotpic ratio.
    // The first cm-1 is the frequency unit (it cancels with the
    // 1/frequency unit of the line shape function).
    //
    // We need to do the following:
    // 1. Convert frequency from wavenumber to Hz (factor 1e2 * c).
    // 2. Convert [molec * cm-2] to [molec * m-2] (factor 1e-4).
    // 3. Take out the isotopologue ratio.
 
    constexpr Numeric hi2arts = 1e-2 * Constant::c;
 
    Numeric s;
 
    // Extract HITRAN intensity:
    extract(s, line, 10);
    // Convert to ARTS units (Hz / (molec * m-2) ), or shorter: Hz*m^2
    data.line.I0 = s * hi2arts;
    // Take out isotopologue ratio:
    data.line.I0 /= Hitran::ratio_from_lookup(mo, iso);
  }
 
  // Einstein coefficient
  {
    extract(data.line.A, line, 10);
  }
 
  // Air broadening parameters.
  Numeric agam, sgam;
  {
    // HITRAN parameter is in cm-1/atm at 296 Kelvin
    // All parameters are HWHM (I hope this is true!)
    Numeric gam;
    // Conversion from wavenumber to Hz. If you multiply a value in
    // wavenumber (cm^-1) by this constant, you get the value in Hz.
    constexpr Numeric w2Hz = Constant::c * 1e2;
    // Ok, put together the end-to-end conversion that we need:
    constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
 
    // Extract HITRAN AGAM value:
    extract(gam, line, 5);
 
    // ARTS parameter in Hz/Pa:
    agam = gam * hi2arts;
 
    // Extract HITRAN SGAM value:
    extract(gam, line, 5);
 
    // ARTS parameter in Hz/Pa:
    sgam = gam * hi2arts;
 
    // If zero, set to agam:
    if (0 == sgam) sgam = agam;
 
    //    cout << "agam, sgam = " << magam << ", " << msgam << endl;
  }
 
  // Lower state energy.
  {
    // HITRAN parameter is in wavenumbers (cm^-1).
    // We have to convert this to the ARTS unit Joule.
 
    // Extract from Catalogue line
    extract(data.line.E0, line, 10);
 
    // Convert to Joule:
    data.line.E0 = wavenumber_to_joule(data.line.E0);
  }
 
  // Temperature coefficient of broadening parameters.
  Numeric nair, nself;
  {
    // This is dimensionless, we can also extract directly.
    extract(nair, line, 4);
 
    // Set self broadening temperature coefficient to the same value:
    nself = nair;
    //    cout << "mnair = " << mnair << endl;
  }
 
  // Pressure shift.
  Numeric psf;
  {
    // HITRAN value in cm^-1 / atm. So the conversion goes exactly as
    // for the broadening parameters.
    Numeric d;
    // Conversion from wavenumber to Hz. If you multiply a value in
    // wavenumber (cm^-1) by this constant, you get the value in Hz.
    constexpr Numeric w2Hz = Constant::c * 1e2;
    // Ok, put together the end-to-end conversion that we need:
    constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
 
    // Extract HITRAN value:
    extract(d, line, 8);
 
    // ARTS value in Hz/Pa
    psf = d * hi2arts;
  }
  // Set the accuracies using the definition of HITRAN
  // indices. If some are missing, they are set to -1.
 
  // Upper state global quanta
  {
    Index eu;
    extract(eu, line, 15);
  }
 
  // Lower state global quanta
  {
    Index el;
    extract(el, line, 15);
  }
 
  // Upper state local quanta
  {
    Index eul;
    extract(eul, line, 15);
  }
 
  // Lower state local quanta
  {
    Index ell;
    extract(ell, line, 15);
  }
 
  // Accuracy index for frequency
  {
    Index df;
    // Extract HITRAN value:
    extract(df, line, 1);
  }
 
  // Accuracy index for intensity
  {
    Index di0;
    // Extract HITRAN value:
    extract(di0, line, 1);
  }
 
  // Accuracy index for air-broadened halfwidth
  {
    Index dgam;
    // Extract HITRAN value:
    extract(dgam, line, 1);
  }
 
  // Accuracy index for self-broadened half-width
  {
    Index dgam;
    // Extract HITRAN value:
    extract(dgam, line, 1);
  }
 
  // Accuracy index for temperature-dependence exponent for agam
  {
    Index dn;
    // Extract HITRAN value:
    extract(dn, line, 1);
  }
 
  // Accuracy index for pressure shift
  {
    Index dpsfi;
    // Extract HITRAN value (given in cm-1):
    extract(dpsfi, line, 1);
  }
 
  // These were all the parameters that we can extract from
  // HITRAN 2004. However, we still have to set the reference temperatures
  // to the appropriate value:
 
  // Reference temperature for Intensity in K.
  data.T0 = 296.0;
 
  // Set line shape computer
  data.line.lineshape = LineShape::hitran_model(sgam, nself, agam, nair, psf);
  {
    Index garbage;
    extract(garbage, line, 13);
 
    // The statistical weights
    extract(data.line.gupp, line, 7);
    extract(data.line.glow, line, 7);
  }
 
  // That's it!
  data.bad = false;
  return data;
}
 
Absorption::SingleLineExternal Absorption::ReadFromHitranOnlineStream(istream& is) {
  // Default data and values for this type
  SingleLineExternal data;
  data.selfbroadening = true;
  data.bathbroadening = true;
  data.lineshapetype = LineShape::Type::VP;
  data.species.resize(2);
  data.species[1] = Species::Species::Bath;
 
  // This contains the rest of the line to parse. At the beginning the
  // entire line. Line gets shorter and shorter as we continue to
  // extract stuff from the beginning.
  String line;
 
  // The first item is the molecule number:
  Index mo;
 
  // Look for more comments?
  bool comment = true;
 
  while (comment) {
    // Return true if eof is reached:
    if (is.eof()) return data;
 
    // Throw runtime_error if stream is bad:
    ARTS_USER_ERROR_IF (!is, "Stream bad.");
 
    // Read line from file into linebuffer:
    getline(is, line);
 
    // It is possible that we were exactly at the end of the file before
    // calling getline. In that case the previous eof() was still false
    // because eof() evaluates only to true if one tries to read after the
    // end of the file. The following check catches this.
    if (line.nelem() == 0 && is.eof()) return data;
 
    // If the catalogue is in dos encoding, throw away the
    // additional carriage return
    if (line[line.nelem() - 1] == 13) {
      line.erase(line.nelem() - 1, 1);
    }
    
    mo = 0;
    // Initialization of mo is important, because mo stays the same
    // if line is empty.
    extract(mo, line, 2);
    comment = false;
  }
 
  // Extract isotopologue:
  char iso;
  extract(iso, line, 1);
  
  // Set line data
  data.quantumidentity = Hitran::id_from_lookup(mo, iso);
  data.species[0] = data.quantumidentity.Species();
 
  // Position.
  {
    // HITRAN position in wavenumbers (cm^-1):
    Numeric v;
    // Conversion from wavenumber to Hz. If you multiply a line
    // position in wavenumber (cm^-1) by this constant, you get the
    // frequency in Hz.
    constexpr Numeric w2Hz = Constant::c * 100.;
 
    // Extract HITRAN postion:
    extract(v, line, 12);
 
    // ARTS position in Hz:
    data.line.F0 = v * w2Hz;
  }
 
  // Intensity.
  {
    // HITRAN intensity is in cm-1/(molec * cm-2) at 296 Kelvin.
    // It already includes the isotpic ratio.
    // The first cm-1 is the frequency unit (it cancels with the
    // 1/frequency unit of the line shape function).
    //
    // We need to do the following:
    // 1. Convert frequency from wavenumber to Hz (factor 1e2 * c).
    // 2. Convert [molec * cm-2] to [molec * m-2] (factor 1e-4).
    // 3. Take out the isotopologue ratio.
 
    constexpr Numeric hi2arts = 1e-2 * Constant::c;
 
    Numeric s;
 
    // Extract HITRAN intensity:
    extract(s, line, 10);
    // Convert to ARTS units (Hz / (molec * m-2) ), or shorter: Hz*m^2
    data.line.I0 = s * hi2arts;
    // Take out isotopologue ratio:
    data.line.I0 /= Hitran::ratio_from_lookup(mo, iso);
  }
 
  // Einstein coefficient
  {
    extract(data.line.A, line, 10);
  }
 
  // Air broadening parameters.
  Numeric agam, sgam;
  {
    // HITRAN parameter is in cm-1/atm at 296 Kelvin
    // All parameters are HWHM (I hope this is true!)
    Numeric gam;
    // Conversion from wavenumber to Hz. If you multiply a value in
    // wavenumber (cm^-1) by this constant, you get the value in Hz.
    constexpr Numeric w2Hz = Constant::c * 1e2;
    // Ok, put together the end-to-end conversion that we need:
    constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
 
    // Extract HITRAN AGAM value:
    extract(gam, line, 5);
 
    // ARTS parameter in Hz/Pa:
    agam = gam * hi2arts;
 
    // Extract HITRAN SGAM value:
    extract(gam, line, 5);
 
    // ARTS parameter in Hz/Pa:
    sgam = gam * hi2arts;
 
    // If zero, set to agam:
    if (0 == sgam) sgam = agam;
 
    //    cout << "agam, sgam = " << magam << ", " << msgam << endl;
  }
 
  // Lower state energy.
  {
    // HITRAN parameter is in wavenumbers (cm^-1).
    // We have to convert this to the ARTS unit Joule.
 
    // Extract from Catalogue line
    extract(data.line.E0, line, 10);
 
    // Convert to Joule:
    data.line.E0 = wavenumber_to_joule(data.line.E0);
  }
 
  // Temperature coefficient of broadening parameters.
  Numeric nair, nself;
  {
    // This is dimensionless, we can also extract directly.
    extract(nair, line, 4);
 
    // Set self broadening temperature coefficient to the same value:
    nself = nair;
    //    cout << "mnair = " << mnair << endl;
  }
 
  // Pressure shift.
  Numeric psf;
  {
    // HITRAN value in cm^-1 / atm. So the conversion goes exactly as
    // for the broadening parameters.
    Numeric d;
    // Conversion from wavenumber to Hz. If you multiply a value in
    // wavenumber (cm^-1) by this constant, you get the value in Hz.
    constexpr Numeric w2Hz = Constant::c * 1e2;
    // Ok, put together the end-to-end conversion that we need:
    constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
 
    // Extract HITRAN value:
    extract(d, line, 8);
 
    // ARTS value in Hz/Pa
    psf = d * hi2arts;
  }
  // Set the accuracies using the definition of HITRAN
  // indices. If some are missing, they are set to -1.
 
  // Upper state global quanta
  {
    Index eu;
    extract(eu, line, 15);
  }
 
  // Lower state global quanta
  {
    Index el;
    extract(el, line, 15);
  }
 
  // Upper state local quanta
  {
    Index eul;
    extract(eul, line, 15);
  }
 
  // Lower state local quanta
  {
    Index ell;
    extract(ell, line, 15);
  }
 
  // Accuracy index for frequency
  {
    Index df;
    // Extract HITRAN value:
    extract(df, line, 1);
  }
 
  // Accuracy index for intensity
  {
    Index di0;
    // Extract HITRAN value:
    extract(di0, line, 1);
  }
 
  // Accuracy index for air-broadened halfwidth
  {
    Index dgam;
    // Extract HITRAN value:
    extract(dgam, line, 1);
  }
 
  // Accuracy index for self-broadened half-width
  {
    Index dgam;
    // Extract HITRAN value:
    extract(dgam, line, 1);
  }
 
  // Accuracy index for temperature-dependence exponent for agam
  {
    Index dn;
    // Extract HITRAN value:
    extract(dn, line, 1);
  }
 
  // Accuracy index for pressure shift
  {
    Index dpsfi;
    // Extract HITRAN value (given in cm-1):
    extract(dpsfi, line, 1);
  }
 
  // These were all the parameters that we can extract from
  // HITRAN 2004. However, we still have to set the reference temperatures
  // to the appropriate value:
 
  // Reference temperature for Intensity in K.
  data.T0 = 296.0;
 
  // Set line shape computer
  data.line.lineshape = LineShape::hitran_model(sgam, nself, agam, nair, psf);
  {
    Index garbage;
    extract(garbage, line, 13);
 
    // The statistical weights
    extract(data.line.gupp, line, 7);
    extract(data.line.glow, line, 7);
  }
  
  // ADD QUANTUM NUMBER PARSING HERE!
  String upper, lower;
  std::stringstream ss;
  ss.str(line);
  ss >> upper >> lower;
  data.quantumidentity.val = Quantum::Number::from_hitran(upper, lower);
 
  // That's it!
  data.bad = false;
  return data;
}
 
Absorption::SingleLineExternal Absorption::ReadFromHitran2001Stream(istream& is) {
  // Default data and values for this type
  SingleLineExternal data;
  data.selfbroadening = true;
  data.bathbroadening = true;
  data.lineshapetype = LineShape::Type::VP;
  data.species.resize(2);
  data.species[1] = Species::Species::Bath;
 
  // This contains the rest of the line to parse. At the beginning the
  // entire line. Line gets shorter and shorter as we continue to
  // extract stuff from the beginning.
  String line;
 
  // The first item is the molecule number:
  Index mo;
  
  // Look for more comments?
  bool comment = true;
  
  while (comment) {
    // Return true if eof is reached:
    if (is.eof()) return data;
    
    // Throw runtime_error if stream is bad:
    ARTS_USER_ERROR_IF (!is, "Stream bad.");
    
    // Read line from file into linebuffer:
    getline(is, line);
    
    // It is possible that we were exactly at the end of the file before
    // calling getline. In that case the previous eof() was still false
    // because eof() evaluates only to true if one tries to read after the
    // end of the file. The following check catches this.
    if (line.nelem() == 0 && is.eof()) return data;
    
    // If the catalogue is in dos encoding, throw away the
    // additional carriage return
    if (line[line.nelem() - 1] == 13) {
      line.erase(line.nelem() - 1, 1);
    }
    
    mo = 0;
    // Initialization of mo is important, because mo stays the same
    // if line is empty.
    extract(mo, line, 2);
    comment = false;
  }
  
  // Extract isotopologue:
  char iso;
  extract(iso, line, 1);
  
  // Set line data
  data.quantumidentity = Hitran::id_from_lookup(mo, iso);
  data.species[0] = data.quantumidentity.Species();
 
  // Position.
  {
    // HITRAN position in wavenumbers (cm^-1):
    Numeric v;
    // Conversion from wavenumber to Hz. If you multiply a line
    // position in wavenumber (cm^-1) by this constant, you get the
    // frequency in Hz.
    constexpr Numeric w2Hz = Constant::c * 100.;
 
    // Extract HITRAN postion:
    extract(v, line, 12);
 
    // ARTS position in Hz:
    data.line.F0 = v * w2Hz;
    //    cout << "mf = " << mf << endl;
  }
 
  // Intensity.
  {
    // HITRAN intensity is in cm-1/(molec * cm-2) at 296 Kelvin.
    // It already includes the isotpic ratio.
    // The first cm-1 is the frequency unit (it cancels with the
    // 1/frequency unit of the line shape function).
    //
    // We need to do the following:
    // 1. Convert frequency from wavenumber to Hz (factor 1e2 * c).
    // 2. Convert [molec * cm-2] to [molec * m-2] (factor 1e-4).
    // 3. Take out the isotopologue ratio.
 
    constexpr Numeric hi2arts = 1e-2 * Constant::c;
 
    Numeric s;
 
    // Extract HITRAN intensity:
    extract(s, line, 10);
    // Convert to ARTS units (Hz / (molec * m-2) ), or shorter: Hz*m^2
    data.line.I0 = s * hi2arts;
    // Take out isotopologue ratio:
    data.line.I0 /= Hitran::ratio_from_lookup(mo, iso);
  }
 
  // Skip transition probability:
  {
    Numeric r;
    extract(r, line, 10);
  }
 
  // Air broadening parameters.
  Numeric agam, sgam;
  {
    // HITRAN parameter is in cm-1/atm at 296 Kelvin
    // All parameters are HWHM (I hope this is true!)
    Numeric gam;
    // Conversion from wavenumber to Hz. If you multiply a value in
    // wavenumber (cm^-1) by this constant, you get the value in Hz.
    constexpr Numeric w2Hz = Constant::c * 1e2;
    // Ok, put together the end-to-end conversion that we need:
    constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
 
    // Extract HITRAN AGAM value:
    extract(gam, line, 5);
 
    // ARTS parameter in Hz/Pa:
    agam = gam * hi2arts;
 
    // Extract HITRAN SGAM value:
    extract(gam, line, 5);
 
    // ARTS parameter in Hz/Pa:
    sgam = gam * hi2arts;
 
    // If zero, set to agam:
    if (0 == sgam) sgam = agam;
 
    //    cout << "agam, sgam = " << magam << ", " << msgam << endl;
  }
 
  // Lower state energy.
  {
    // HITRAN parameter is in wavenumbers (cm^-1).
    // We have to convert this to the ARTS unit Joule.
 
    // Extract from Catalogue line
    extract(data.line.E0, line, 10);
 
    // Convert to Joule:
    data.line.E0 = wavenumber_to_joule(data.line.E0);
  }
 
  // Temperature coefficient of broadening parameters.
  Numeric nair, nself;
  {
    // This is dimensionless, we can also extract directly.
    extract(nair, line, 4);
 
    // Set self broadening temperature coefficient to the same value:
    nself = nair;
    //    cout << "mnair = " << mnair << endl;
  }
 
  // Pressure shift.
  Numeric psf;
  {
    // HITRAN value in cm^-1 / atm. So the conversion goes exactly as
    // for the broadening parameters.
    Numeric d;
    // Conversion from wavenumber to Hz. If you multiply a value in
    // wavenumber (cm^-1) by this constant, you get the value in Hz.
    constexpr Numeric w2Hz = Constant::c * 1e2;
    // Ok, put together the end-to-end conversion that we need:
    constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
 
    // Extract HITRAN value:
    extract(d, line, 8);
 
    // ARTS value in Hz/Pa
    psf = d * hi2arts;
  }
  // Set the accuracies using the definition of HITRAN
  // indices. If some are missing, they are set to -1.
 
  //Skip upper state global quanta index
  {
    Index eu;
    extract(eu, line, 3);
  }
 
  //Skip lower state global quanta index
  {
    Index el;
    extract(el, line, 3);
  }
 
  //Skip upper state local quanta
  {
    Index eul;
    extract(eul, line, 9);
  }
 
  //Skip lower state local quanta
  {
    Index ell;
    extract(ell, line, 9);
  }
 
  // Accuracy index for frequency reference
  {
    Index df;
    // Extract HITRAN value:
    extract(df, line, 1);
  }
 
  // Accuracy index for intensity reference
  {
    Index di0;
    // Extract HITRAN value:
    extract(di0, line, 1);
  }
 
  // Accuracy index for halfwidth reference
  {
    Index dgam;
    // Extract HITRAN value:
    extract(dgam, line, 1);
  }
 
  // These were all the parameters that we can extract from
  // HITRAN. However, we still have to set the reference temperatures
  // to the appropriate value:
 
  // Reference temperature for Intensity in K.
  data.T0 = 296.0;
 
  // Set line shape computer
  data.line.lineshape = LineShape::hitran_model(sgam, nself, agam, nair, psf);
 
  // That's it!
  data.bad = false;
  return data;
}
 
Absorption::SingleLineExternal Absorption::ReadFromLBLRTMStream(istream& is) {
  // Default data and values for this type
  SingleLineExternal data;
  data.selfbroadening = true;
  data.bathbroadening = true;
  data.lineshapetype = LineShape::Type::VP;
  data.species.resize(2);
 
  // This contains the rest of the line to parse. At the beginning the
  // entire line. Line gets shorter and shorter as we continue to
  // extract stuff from the beginning.
  String line;
 
  // The first item is the molecule number:
  Index mo;
  
  // Look for more comments?
  bool comment = true;
  
  while (comment) {
    // Return true if eof is reached:
    if (is.eof()) return data;
    
    // Throw runtime_error if stream is bad:
    ARTS_USER_ERROR_IF (!is, "Stream bad.");
    
    // Read line from file into linebuffer:
    getline(is, line);
    if (line[0] == '>' or line[0] == '%') continue;
    
    // It is possible that we were exactly at the end of the file before
    // calling getline. In that case the previous eof() was still false
    // because eof() evaluates only to true if one tries to read after the
    // end of the file. The following check catches this.
    if (line.nelem() == 0 && is.eof()) return data;
    
    // If the catalogue is in dos encoding, throw away the
    // additional carriage return
    if (line[line.nelem() - 1] == 13) {
      line.erase(line.nelem() - 1, 1);
    }
    
    mo = 0;
    // Initialization of mo is important, because mo stays the same
    // if line is empty.
    extract(mo, line, 2);
    comment = false;
  }
  
  // Extract isotopologue:
  char iso;
  extract(iso, line, 1);
  
  // Set line data
  data.quantumidentity = Hitran::id_from_lookup(mo, iso);
 
  // Position.
  {
    // HITRAN position in wavenumbers (cm^-1):
    Numeric v;
    // Conversion from wavenumber to Hz. If you multiply a line
    // position in wavenumber (cm^-1) by this constant, you get the
    // frequency in Hz.
    constexpr Numeric w2Hz = Constant::c * 100.;
 
    // Extract HITRAN postion:
    extract(v, line, 12);
 
    // ARTS position in Hz:
    data.line.F0 = v * w2Hz;
    //    cout << "mf = " << mf << endl;
  }
 
  // Intensity.
  {
    // HITRAN intensity is in cm-1/(molec * cm-2) at 296 Kelvin.
    // It already includes the isotpic ratio.
    // The first cm-1 is the frequency unit (it cancels with the
    // 1/frequency unit of the line shape function).
    //
    // We need to do the following:
    // 1. Convert frequency from wavenumber to Hz (factor 1e2 * c).
    // 2. Convert [molec * cm-2] to [molec * m-2] (factor 1e-4).
    // 3. Take out the isotopologue ratio.
 
    constexpr Numeric hi2arts = 1e-2 * Constant::c;
 
    Numeric s;
    if (line[6] == 'D') line[6] = 'E';
    // Extract HITRAN intensity:
    extract(s,
            line,
            10);  // NOTE:  If error shooting, FORTRAN "D" is not read properly.
    // Convert to ARTS units (Hz / (molec * m-2) ), or shorter: Hz*m^2
    data.line.I0 = s * hi2arts;
    // Take out isotopologue ratio:
    data.line.I0 /=  Hitran::ratio_from_lookup(mo, iso);
  }
 
  // Skip transition probability:
  {
    Numeric r;
    extract(r, line, 10);
  }
 
  // Air broadening parameters.
  Numeric sgam, agam;
  {
    // HITRAN parameter is in cm-1/atm at 296 Kelvin
    // All parameters are HWHM (I hope this is true!)
    Numeric gam;
    // Conversion from wavenumber to Hz. If you multiply a value in
    // wavenumber (cm^-1) by this constant, you get the value in Hz.
    constexpr Numeric w2Hz = Constant::c * 1e2;
    // Ok, put together the end-to-end conversion that we need:
    constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
 
    // Extract HITRAN AGAM value:
    extract(gam, line, 5);
 
    // ARTS parameter in Hz/Pa:
    agam = gam * hi2arts;
 
    // Extract HITRAN SGAM value:
    extract(gam, line, 5);
 
    // ARTS parameter in Hz/Pa:
    sgam = gam * hi2arts;
 
    // If zero, set to agam:
    if (0 == sgam) sgam = agam;
 
    //    cout << "agam, sgam = " << magam << ", " << msgam << endl;
  }
 
  // Lower state energy.
  {
    // HITRAN parameter is in wavenumbers (cm^-1).
    // We have to convert this to the ARTS unit Joule.
 
    // Extract from Catalogue line
    extract(data.line.E0, line, 10);
 
    // Convert to Joule:
    data.line.E0 = wavenumber_to_joule(data.line.E0);
  }
 
  // Temperature coefficient of broadening parameters.
  Numeric nair, nself;
  {
    // This is dimensionless, we can also extract directly.
    extract(nair, line, 4);
 
    // Set self broadening temperature coefficient to the same value:
    nself = nair;
    //    cout << "mnair = " << mnair << endl;
  }
 
  // Pressure shift.
  Numeric psf;
  {
    // HITRAN value in cm^-1 / atm. So the conversion goes exactly as
    // for the broadening parameters.
    Numeric d;
    // Conversion from wavenumber to Hz. If you multiply a value in
    // wavenumber (cm^-1) by this constant, you get the value in Hz.
    constexpr Numeric w2Hz = Constant::c * 1e2;
    // Ok, put together the end-to-end conversion that we need:
    constexpr Numeric hi2arts = w2Hz / Conversion::atm2pa(1);
 
    // Extract HITRAN value:
    extract(d, line, 8);
 
    // ARTS value in Hz/Pa
    psf = d * hi2arts;
  }
  // Set the accuracies using the definition of HITRAN
  // indices. If some are missing, they are set to -1.
 
  //Skip upper state global quanta index
  {
    Index eu;
    extract(eu, line, 3);
  }
 
  //Skip lower state global quanta index
  {
    Index el;
    extract(el, line, 3);
  }
 
  //Skip upper state local quanta
  {
    Index eul;
    extract(eul, line, 9);
  }
 
  //Skip lower state local quanta
  {
    Index ell;
    extract(ell, line, 9);
  }
 
  // Accuracy index for frequency reference
  {
    Index df;
    // Extract HITRAN value:
    extract(df, line, 1);
  }
 
  // Accuracy index for intensity reference
  {
    Index di0;
    // Extract HITRAN value:
    extract(di0, line, 1);
  }
 
  // Accuracy index for halfwidth reference
  {
    Index dgam;
    // Extract HITRAN value:
    extract(dgam, line, 1);
  }
 
  // These were all the parameters that we can extract from
  // HITRAN. However, we still have to set the reference temperatures
  // to the appropriate value:
 
  // Reference temperature for Intensity in K.
  // (This is fix for HITRAN)
  data.T0 = 296.0;
 
  // Skip four
  {
    Index four;
    extract(four, line, 4);
  }
 
  // This is the test for the last two characters of the line
  {
    /* 
     *   0 is nothing, 
     *  -1 is linemixing on the next line, 
     *  -3 is the non-resonant line 
     */
    Index test;
    extract(test, line, 2);
    //If the tag is as it should be, then a minus one means that more should be read
    if (test == -1 || test == -3)
      getline(is, line);
    else  // the line is done and we are happy to leave
    {
      data.line.lineshape = LineShape::hitran_model(sgam, nself, agam, nair, psf);
      
      data.bad = false;
      return data;
    }
  }
 
  // In case we are unable to leave, the next line is a line mixing parameter line
 
  // First is the molecular number.  This should be the same as above.
  {
    Index mo2;
    extract(mo2, line, 2);
    // Skip one
 
    ARTS_USER_ERROR_IF (mo != mo2,
                        "There is an error in the line mixing\n");
  }
 
  Vector Y(4), G(4), T(4);
 
  // These are constants for AER but should be included because we need their grid.
  T[0] = 200;
  T[1] = 250;
  T[2] = 296;
  T[3] = 340;
 
  // Next is the Y  and G at various temperatures
  {
    Numeric Y_200K;
    extract(Y_200K, line, 13);
    Y[0] = Y_200K;
  }
  {
    Numeric G_200K;
    extract(G_200K, line, 11);
    G[0] = G_200K;
  }
  {
    Numeric Y_250K;
    extract(Y_250K, line, 13);
    Y[1] = Y_250K;
  }
  {
    Numeric G_250K;
    extract(G_250K, line, 11);
    G[1] = G_250K;
  }
  {
    Numeric Y_296K;
    extract(Y_296K, line, 13);
    Y[2] = Y_296K;
  }
  {
    Numeric G_296K;
    extract(G_296K, line, 11);
    G[2] = G_296K;
  }
  {
    Numeric Y_340K;
    extract(Y_340K, line, 13);
    Y[3] = Y_340K;
  }
  {
    Numeric G_340K;
    extract(G_340K, line, 11);
    G[3] = G_340K;
  }
 
  // Cpnvert from per Atm and per Atm^2
  Y /= Conversion::atm2pa(1);
  G /= Constant::pow2(Conversion::atm2pa(1));
  
  // ARTS uses (1-iY) as line-mixing factor, LBLRTM uses (1+iY), so we must change sign
  Y *= -1;
  
  // Test that this is the end
  {
    Index test;
    extract(test, line, 2);
    if (test == -1) {
      data.line.lineshape = LineShape::lblrtm_model(sgam,
                                         nself,
                                         agam,
                                         nair,
                                         psf,
                                         {T[0],
                                          T[1],
                                          T[2],
                                          T[3],
                                          Y[0],
                                          Y[1],
                                          Y[2],
                                          Y[3],
                                          G[0],
                                          G[1],
                                          G[2],
                                          G[3]});
 
      data.bad = false;
      return data;
    }
    if (test == -3) {
      data.line.lineshape = LineShape::lblrtm_model(sgam,
                                         nself,
                                         agam,
                                         nair,
                                         psf,
                                         {T[0],
                                          T[1],
                                          T[2],
                                          T[3],
                                          Y[0],
                                          Y[1],
                                          Y[2],
                                          Y[3],
                                          G[0],
                                          G[1],
                                          G[2],
                                          G[3]});
      
      data.bad = false;
      return data;
    }
    return data;
  }
}
 
std::vector<Absorption::Lines> Absorption::split_list_of_external_lines(std::vector<SingleLineExternal>& external_lines,
                                                                        const std::vector<QuantumNumberType>& localquantas,
                                                                        const std::vector<QuantumNumberType>& globalquantas)
{
  std::vector<Lines> lines(0);
  
  // Loop but make copies because we will need to modify some of the data
  while(external_lines.size()) {
    auto& sle = external_lines.back();
    
    // Adapt broadening to fit with line catalog
    if (sle.selfbroadening) sle.species.front() = sle.quantumidentity.Species();
    if (sle.bathbroadening) sle.species.back() = Species::Species::Bath;
 
    Quantum::Number::GlobalState global_id{sle.quantumidentity.isotopologue_index};
    for (auto qn: globalquantas) if (sle.quantumidentity.val.has(qn)) global_id.val.set(sle.quantumidentity.val[qn]);
 
    Quantum::Number::LocalState local_id{};
    for (auto qn: localquantas) if (sle.quantumidentity.val.has(qn)) local_id.val.set(sle.quantumidentity.val[qn]);
    
    // Set the line
    auto line = sle.line;
    line.localquanta = local_id;
    
    // Either find a line like this in the list of lines or start a new Lines
    auto band = std::find_if(lines.begin(), lines.end(), [&](const Lines& li){return li.MatchWithExternal(sle, global_id);});
    if (band not_eq lines.end()) {
      band -> AppendSingleLine(line);
    } else {
      lines.push_back(Lines(sle.selfbroadening, sle.bathbroadening, sle.cutoff,
                            sle.mirroring, sle.population, sle.normalization,
                            sle.lineshapetype, sle.T0, sle.cutofffreq,
                            sle.linemixinglimit, global_id, sle.species, {line}));
    }
    external_lines.pop_back();
  }
  
  return lines;
}
 
std::ostream& Absorption::operator<<(std::ostream& os, const Absorption::Lines& lines)
{
  for(auto& line: lines.lines)
    os << line << '\n';
  return os;
}
 
std::istream& Absorption::operator>>(std::istream& is, Lines& lines) {
  for(auto& line: lines.lines)
    is >> line;
  return is;
}
 
std::ostream& Absorption::operator<<(std::ostream& os,
                                     const Absorption::SingleLine& line) {
  os << line.F0 << ' ' << line.I0 << ' ' << line.E0 << ' ' << line.glow << ' '
     << line.gupp << ' ' << line.A << ' ' << line.zeeman << ' '
     << line.lineshape;
  if (line.localquanta.val.nelem()) os << ' ' << line.localquanta.values();
  return os;
}
 
std::istream& Absorption::operator>>(std::istream& is, Absorption::SingleLine& line)
{
  is >> double_imanip() >> line.F0 >> line.I0 >> line.E0 >> line.glow >>
      line.gupp >> line.A;
  return is >> line.zeeman >> line.lineshape >> line.localquanta;
}
 
std::ostream& operator<<(std::ostream& os, const ArrayOfAbsorptionLines& aol)
{
  for(auto& l: aol)
    os << l << '\n';
  return os;
}
 
std::ostream& operator<<(std::ostream& os, const ArrayOfArrayOfAbsorptionLines& aol)
{
  for(auto& l: aol)
    os << l << '\n';
  return os;
}
 
String Absorption::Lines::SpeciesName() const noexcept
{
  return quantumidentity.Isotopologue().FullName();
}
 
String Absorption::Lines::MetaData() const
{
  std::ostringstream os;
  
  os << "\nLines meta-data:\n";
  os << '\t' << "Species identity:\n";
  os << "\t\tSpecies: "<< SpeciesName() << '\n';
  os << "\t\tIdentity: "<< quantumidentity << '\n';
  os << '\t' << cutofftype2metadatastring(cutoff, cutofffreq);
  os << '\t' << populationtype2metadatastring(population);
  os << '\t' << normalizationtype2metadatastring(normalization);
  os << '\t' << LineShape::shapetype2metadatastring(lineshapetype);
  os << '\t' << mirroringtype2metadatastring(mirroring);
  os << '\t' << "The reference temperature for all line parameters is "
     << T0 << " K.\n";
  if(linemixinglimit < 0)
    os << '\t' << "If applicable, there is no line mixing limit.\n";
  else
    os << '\t' << "If applicable, there is a line mixing limit at "
       << linemixinglimit << " Pa.\n";
  
  if (not NumLines()) {
    os << "\tNo line data is available.\n";
  } else {
    os << "\tThere are " << NumLines() << " lines available.\n";
    
    auto& line = lines.front();
    os << "\tThe front line has:\n";
    os << "\t\t" << "f0: " << line.F0 << " Hz\n";
    os << "\t\t" << "i0: " << line.I0 << " m^2/Hz\n";
    os << "\t\t" << "e0: " << line.E0 << " J\n";
    os << "\t\t" << "Lower stat. weight: " << line.glow << " [-]\n";
    os << "\t\t" << "Upper stat. weight: " << line.gupp << " [-]\n";
    os << "\t\t" << "A: " << line.A << " 1/s\n";
    os << "\t\t" << "Zeeman splitting of lower state: " << line.zeeman.gl() << " [-]\n";
    os << "\t\t" << "Zeeman splitting of upper state: " << line.zeeman.gu() << " [-]\n";
    os << "\t\t" << "Local quantum numbers: " << line.localquanta.val << "\n";
 
    ArrayOfString ls_meta = LineShape::ModelMetaDataArray(line.lineshape,
                                                          selfbroadening,
                                                          broadeningspecies,
                                                          T0);
    os << "\t\t" << "Line shape parameters (are normalized by sum(VMR)):\n";
    for(auto& ls_form: ls_meta)
      os << "\t\t\t" << ls_form << "\n";
  }
  
  
  return os.str();
}
 
 
void Absorption::Lines::RemoveLine(Index i) noexcept
{
  lines.erase(lines.begin() + i);
}
 
 
Absorption::SingleLine Absorption::Lines::PopLine(Index i) noexcept
{
  auto line = lines[i];
  RemoveLine(i);
  return line;
}
 
void Absorption::Lines::ReverseLines() noexcept {
  std::reverse(lines.begin(), lines.end());
}
 
Numeric Absorption::Lines::SpeciesMass() const noexcept {
  return quantumidentity.Isotopologue().mass;
}
 
Vector Absorption::Lines::BroadeningSpeciesVMR(const ConstVectorView& atm_vmrs,
                                               const ArrayOfArrayOfSpeciesTag& atm_spec) const
{
  if (lineshapetype == LineShape::Type::DP)
    return Vector(broadeningspecies.nelem(), std::numeric_limits<Numeric>::quiet_NaN());
  return LineShape::vmrs(atm_vmrs, atm_spec, broadeningspecies);
}
 
Vector Absorption::Lines::BroadeningSpeciesMass(const ConstVectorView& atm_vmrs,
                                                const ArrayOfArrayOfSpeciesTag& atm_spec,
                                                const SpeciesIsotopologueRatios& ir,
                                                const Numeric& bath_mass
                                               ) const
{
  if (lineshapetype == LineShape::Type::DP)
    return Vector(broadeningspecies.nelem(), std::numeric_limits<Numeric>::quiet_NaN());
  Vector mass = LineShape::mass(atm_vmrs, atm_spec, broadeningspecies, ir);
  if (bathbroadening and bath_mass > 0) mass[mass.nelem()-1] = bath_mass;
  return mass;
}
 
Numeric Absorption::Lines::SelfVMR(const ConstVectorView& atm_vmrs,
                                   const ArrayOfArrayOfSpeciesTag& atm_spec) const
{
  ARTS_USER_ERROR_IF (atm_vmrs.nelem() not_eq atm_spec.nelem(),
                      "Bad species and vmr lists");
    
  for (Index i=0; i<atm_spec.nelem(); i++)
    if (atm_spec[i].nelem() and atm_spec[i].Species() == Species())
      return atm_vmrs[i];
  return 0;
}
 
Numeric Absorption::reduced_rovibrational_dipole(Rational Jf, Rational Ji, Rational lf, Rational li, Rational k) {
  if (not iseven(Jf + lf + 1))
    return - sqrt(2 * Jf + 1) * wigner3j(Jf, k, Ji, li, lf - li, -lf);
  return + sqrt(2 * Jf + 1) * wigner3j(Jf, k, Ji, li, lf - li, -lf);
}
 
Numeric Absorption::reduced_magnetic_quadrapole(Rational Jf, Rational Ji, Rational N) {
  if (not iseven(Jf + N))
    return - sqrt(6 * (2 * Jf + 1) * (2 * Ji + 1)) * wigner6j(1, 1, 1, Ji, Jf, N);
  return + sqrt(6 * (2 * Jf + 1) * (2 * Ji + 1)) * wigner6j(1, 1, 1, Ji, Jf, N);
}
 
Absorption::SingleLineExternal Absorption::ReadFromJplStream(istream& is)
{
  // Default data and values for this type
  SingleLineExternal data;
  data.selfbroadening = true;
  data.bathbroadening = true;
  data.lineshapetype = LineShape::Type::VP;
  data.species.resize(2);
 
  // This contains the rest of the line to parse. At the beginning the
  // entire line. Line gets shorter and shorter as we continue to
  // extract stuff from the beginning.
  String line;
 
  // Look for more comments?
  bool comment = true;
 
  while (comment) {
    // Return true if eof is reached:
    if (is.eof()) return data;
 
    // Throw runtime_error if stream is bad:
    ARTS_USER_ERROR_IF (!is, "Stream bad.");
 
    // Read line from file into linebuffer:
    getline(is, line);
 
    // It is possible that we were exactly at the end of the file before
    // calling getline. In that case the previous eof() was still false
    // because eof() evaluates only to true if one tries to read after the
    // end of the file. The following check catches this.
    if (line.nelem() == 0 && is.eof()) return data;
 
    // Because of the fixed FORTRAN format, we need to break up the line
    // explicitly in apropriate pieces. Not elegant, but works!
 
    // Extract center frequency:
    // Initialization of v is important, because v stays the same
    // if line is empty.
    // JPL position in MHz:
    Numeric v = 0.0;
 
    // Extract JPL position:
    extract(v, line, 13);
 
    // check for empty line
    if (v != 0.0) {
      // ARTS position in Hz:
      data.line.F0 = v * 1E6;
 
      comment = false;
    }
  }
 
  // Accuracy for line position
  {
    Numeric df;
    extract(df, line, 8);
  }
 
  // Intensity.
  {
    // JPL has log (10) of intensity in nm2 MHz at 300 Kelvin.
    //
    // We need to do the following:
    // 1. take 10^intensity
    // 2. convert to cm-1/(molecule * cm-2): devide by c * 1e10
    // 3. Convert frequency from wavenumber to Hz (factor 1e2 * c)
    // 4. Convert [molec * cm-2] to [molec * m-2] (factor 1e-4)
 
    Numeric s;
 
    // Extract JPL intensity:
    extract(s, line, 8);
 
    // remove log
    s = pow((Numeric)10., s);
 
    // Convert to ARTS units (Hz / (molec * m-2) ), or shorter: Hz*m^2
    data.line.I0 = s / 1E12;
  }
 
  // Degrees of freedom
  {
    Index dr;
 
    // Extract degrees of freedom
    extract(dr, line, 2);
  }
 
  // Lower state energy.
  {
    // JPL parameter is in wavenumbers (cm^-1).
    // We have to convert this to the ARTS unit Joule.
 
    // Extract from Catalogue line
    extract(data.line.E0, line, 10);
 
    // Convert to Joule:
    data.line.E0 = wavenumber_to_joule(data.line.E0);
  }
 
  // Upper state degeneracy
  {
    Index gup;
 
    // Extract upper state degeneracy
    extract(gup, line, 3);
  }
 
  // Tag number
  Index tag;
  {
    // Extract Tag number
    extract(tag, line, 7);
 
    // make sure tag is not negative (damned jpl cat):
    tag = tag > 0 ? tag : -tag;
  }
 
  // Set line ID
  data.quantumidentity = Jpl::id_from_lookup(tag);
 
  // Air broadening parameters: unknown to jpl, use old iup forward
  // model default values, which is mostly set to 0.0025 GHz/hPa, even
  // though for some lines the pressure broadening is given explicitly
  // in the program code. The explicitly given values are ignored and
  // only the default value is set. Self broadening was in general not
  // considered in the old forward model.
  Numeric agam, sgam;
  {
    // ARTS parameter in Hz/Pa:
    agam = 2.5E4;
 
    // ARTS parameter in Hz/Pa:
    sgam = agam;
  }
 
  // Temperature coefficient of broadening parameters. Was set to 0.75
  // in old forward model, even though for some lines the parameter is
  // given explicitly in the program code. The explicitly given values
  // are ignored and only the default value is set. Self broadening
  // not considered.
  Numeric nair, nself;
  {
    nair = 0.75;
    nself = 0.0;
  }
 
  // Reference temperature for broadening parameter in K, was
  // generally set to 300 K in old forward model, with the exceptions
  // as already mentioned above: //DEPRECEATED but is same as for mti0 so moving on
  //   {
  //     mtgam = 300.0;
  //   }
 
  // Pressure shift: not given in JPL, set to 0
  Numeric psf;
  { psf = 0.0; }
 
  // These were all the parameters that we can extract from
  // JPL. However, we still have to set the reference temperatures
  // to the appropriate value:
 
  // Reference temperature for Intensity in K.
  data.T0 = 300.0;
 
  // Set line shape computer
  data.line.lineshape = LineShape::Model(sgam, nself, agam, nair, psf);
 
  // That's it!
  data.bad = false;
  return data;
}
 
bool Absorption::Lines::OK() const noexcept
{
  const Index nb = broadeningspecies.nelem();
 
  // Check that the isotopologue is ok
  if (not Isotopologue().OK()) return false;
  
  // Test that self and bath is covered by the range if set positive
  if (nb < (Index(selfbroadening) + Index(bathbroadening)))
    return false;
 
  // Test that the temperature is physical
  if (T0 <= 0)
    return false;
  
  // Test that all lines have the correct sized line shape model
  if (std::any_of(lines.cbegin(), lines.cend(), [nb](auto& line){return line.LineShapeElems() != nb;}))
    return false;
  
  // Test that all lines have the correct sized local quantum numbers
  if (std::any_of(lines.cbegin(), lines.cend(), [&](auto& line){return line.LocalQuantumElems() != lines.front().LocalQuantumElems();}))
    return false;
  
  // Otherwise everything is fine!
  return true;
}
 
Numeric Absorption::Lines::DopplerConstant(Numeric T) const noexcept
{
  return std::sqrt(Constant::doppler_broadening_const_squared * T / SpeciesMass());
}
 
 
 
namespace Absorption {
String cutofftype2metadatastring(CutoffType in, Numeric cutoff) {
  std::ostringstream os;
  switch (in) {
    case CutoffType::None:
      os << "No cut-off will be applied.\n"; break;
    case CutoffType::ByLine:
      os << "The lines will be cut-off " << cutoff << " Hz from the line center + D0.\n"; break;
    case CutoffType::FINAL: break;
  }
  return os.str();
}
 
void SingleLine::SetAutomaticZeeman(QuantumIdentifier qid) {
  for (auto& val: localquanta.val) qid.val.set(val); // Copy to same struct
  
  zeeman = Zeeman::Model(qid);
}
 
void SingleLine::SetLineMixing2SecondOrderData(const Vector& d) {
  lineshape.SetLineMixingModel(
    LineShape::LegacyLineMixingData::vector2modellm(
      d, LineShape::LegacyLineMixingData::TypeLM::LM_2NDORDER)
    .Data()[0]);
}
 
void SingleLine::SetLineMixing2AER(const Vector& d) {
  const LineShape::ModelParameters Y = {LineShape::TemperatureModel::LM_AER, d[4], d[5], d[6], d[7]};
  const LineShape::ModelParameters G = {LineShape::TemperatureModel::LM_AER, d[8], d[9], d[10], d[11]};
  for (auto& sm : lineshape.Data()) {
    sm.Y() = Y;
    sm.G() = G;
  }
}
 
bifstream& SingleLine::read(bifstream& bif) {
  bif >> F0 >> I0 >> E0 >> glow >> gupp >> A >> zeeman;
  
  lineshape.read(bif);
  
  for (auto& val: localquanta.val) val.read(bif);
  
  return bif;
}
 
bofstream& SingleLine::write(bofstream& bof) const {
  bof << F0 << I0 << E0 << glow << gupp << A << zeeman;
  
  lineshape.write(bof);
  
  for (auto& val: localquanta.val) val.write(bof);
  
  return bof;
}
 
void Lines::AppendSingleLine(SingleLine&& sl) {
  ARTS_USER_ERROR_IF(
      NumLines() not_eq 0 and NumLocalQuanta() not_eq sl.LocalQuantumElems(),
      "Error calling appending function, bad size of quantum numbers");
 
  ARTS_USER_ERROR_IF(NumLines() not_eq 0 and 
    sl.LineShapeElems() not_eq lines[0].LineShapeElems(),
                     "Error calling appending function, bad size of broadening species");
  
  lines.push_back(std::move(sl));
}
 
void Lines::AppendSingleLine(const SingleLine& sl) {
  ARTS_USER_ERROR_IF(
      NumLines() not_eq 0 and NumLocalQuanta() not_eq sl.LocalQuantumElems(),
      "Error calling appending function, bad size of quantum numbers\n"
      "Type of quantum numbers in band: ", lines.front().localquanta.val,
      "\nType of quantum numbers in new line:", sl.localquanta.val
      );
 
  ARTS_USER_ERROR_IF(
      NumLines() not_eq 0 and
          sl.LineShapeElems() not_eq lines.front().LineShapeElems(),
      "Error calling appending function, bad size of broadening species");
 
  lines.push_back(sl);
}
 
bool Lines::MatchWithExternal(const SingleLineExternal& sle, const QuantumIdentifier& qid) const ARTS_NOEXCEPT {
  if(sle.bad)
    return false;
  if(sle.selfbroadening not_eq selfbroadening)
    return false;
  if(sle.bathbroadening not_eq bathbroadening)
    return false;
  if(sle.cutoff not_eq cutoff)
    return false;
  if(sle.mirroring not_eq mirroring)
    return false;
  if(sle.population not_eq population)
    return false;
  if(sle.normalization not_eq normalization)
    return false;
  if(sle.lineshapetype not_eq lineshapetype)
    return false;
  if(sle.T0 not_eq T0)
    return false;
  if(sle.cutofffreq not_eq cutofffreq)
    return false;
  if(sle.linemixinglimit not_eq linemixinglimit)
    return false;
  if(quantumidentity not_eq qid)
    return false;
  if(not std::equal(sle.species.cbegin(), sle.species.cend(), broadeningspecies.cbegin(), broadeningspecies.cend()))
    return false;
  if(NumLines() not_eq 0 and not lines.front().localquanta.same_types_as(sle.line.localquanta))
    return false;
  if(NumLines() not_eq 0 and not sle.line.lineshape.Match(lines.front().lineshape).first)
    return false;
  return true;
}
 
std::pair<bool, bool> Lines::Match(const Lines& l) const noexcept {
  // Note: The pair here is first: matching and second: nullable
  if(l.selfbroadening not_eq selfbroadening)
    return {false, false};
  if(l.bathbroadening not_eq bathbroadening)
    return {false, false};
  if(l.cutoff not_eq cutoff)
    return {false, false};
  if(l.mirroring not_eq mirroring)
    return {false, false};
  if(l.population not_eq population)
    return {false, false};
  if(l.normalization not_eq normalization)
    return {false, false};
  if(l.lineshapetype not_eq lineshapetype)
    return {false, false};
  if(l.T0 not_eq T0)
    return {false, false};
  if(l.cutofffreq not_eq cutofffreq)
    return {false, false};
  if(l.linemixinglimit not_eq linemixinglimit)
    return {false, false};
  if(l.quantumidentity not_eq quantumidentity)
    return {false, false};
  if(not std::equal(l.broadeningspecies.cbegin(), l.broadeningspecies.cend(), broadeningspecies.cbegin(), broadeningspecies.cend()))
    return {false, false};
  if(NumLines() not_eq 0 and l.NumLines() not_eq 0 and not lines.front().localquanta.same_types_as(l.lines.front().localquanta))
    return {false, false};
  if(NumLines() not_eq 0 and l.NumLines() not_eq 0) {
    if (auto matchpair = l.lines.front().lineshape.Match(lines.front().lineshape); not matchpair.first) return matchpair;
  }
  
  return {true, true};
}
 
void Lines::sort_by_frequency() {
  std::sort(lines.begin(), lines.end(),
            [](const SingleLine& a, const SingleLine& b){return a.F0 < b.F0;});
}
 
void Lines::sort_by_einstein() {
  std::sort(lines.begin(), lines.end(),
            [](const SingleLine& a, const SingleLine& b){return a.A < b.A;});
}
 
String Lines::LineShapeMetaData() const noexcept {
  return NumLines() ?
  LineShape::ModelShape2MetaData(lines[0].lineshape) :
  "";
}
 
const Quantum::Number::Value& get(const Quantum::Number::LocalState& qns) ARTS_NOEXCEPT {
  ARTS_ASSERT(qns.val.has(QuantumNumberType::F) or qns.val.has(QuantumNumberType::J))
  return qns.val.has(QuantumNumberType::F) ? qns.val[QuantumNumberType::F] : qns.val[QuantumNumberType::J];
}
 
Index Lines::ZeemanCount(size_t k, Zeeman::Polarization type) const noexcept {
  if (type == Zeeman::Polarization::None) return 1;
 
  // Select F before J but assume one of them exist
  auto& val = get(lines[k].localquanta);
  return Zeeman::nelem(val.upp(), val.low(), type);
}
 
Numeric Lines::ZeemanStrength(size_t k,
                              Zeeman::Polarization type,
                              Index i) const noexcept {
  if (type == Zeeman::Polarization::None) return 1.0;
 
  // Select F before J but assume one of them exist
  auto& val = get(lines[k].localquanta);
  return lines[k].zeeman.Strength(val.upp(), val.low(), type, i);
}
 
Numeric Lines::ZeemanSplitting(size_t k,
                               Zeeman::Polarization type,
                               Index i) const noexcept {
  if (type == Zeeman::Polarization::None) return 0.0;
 
  // Select F before J but assume one of them exist
  auto& val = get(lines[k].localquanta);
  return lines[k].zeeman.Splitting(val.upp(), val.low(), type, i);
}
 
void Lines::SetAutomaticZeeman() noexcept {
  for(auto& line: lines)
    line.SetAutomaticZeeman(quantumidentity);
}
 
Numeric Lines::F_mean(const ConstVectorView& wgts) const noexcept {
  const Numeric val = std::inner_product(lines.cbegin(), lines.cend(),
                                         wgts.begin(), 0.0, std::plus<>(),
                                         [](const auto& a, const auto& b){return a.F0 * b;});
  const Numeric div = wgts.sum();
  return  val / div;
}
 
Numeric Lines::F_mean(Numeric T) const noexcept {
  if (T <= 0) T = T0;
 
  const Index n = NumLines();
  const Numeric QT = single_partition_function(T, Isotopologue());
  const Numeric QT0 = single_partition_function(T0, Isotopologue());
  const Numeric ratiopart = QT0 / QT;
 
  Vector wgts(n);
  for (Index i=0; i<n; i++) {
    const Numeric pop0 = (lines[i].gupp / QT0) * boltzman_factor(T0, lines[i].E0);
    const Numeric pop = pop0 * ratiopart * boltzman_ratio(T, T0, lines[i].E0);
    const Numeric dip_squared = - lines[i].I0/(pop0 * lines[i].F0 * std::expm1(- (Constant::h * lines[i].F0) / (Constant::k * T0)));
    wgts[i] = pop * dip_squared;
  }
 
  return F_mean(wgts);
}
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wreturn-type"
Numeric Lines::CutoffFreq(size_t k, Numeric shift) const noexcept {
  switch(cutoff) {
    case CutoffType::ByLine:
      return lines[k].F0 + cutofffreq + (mirroring == MirroringType::Manual ? - shift : shift);
    case CutoffType::None:
      return std::numeric_limits<Numeric>::max();
    case CutoffType::FINAL: break;
  }
}
#pragma GCC diagnostic pop
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wreturn-type"
Numeric Lines::CutoffFreqMinus(size_t k, Numeric shift) const noexcept {
  switch(cutoff) {
    case CutoffType::ByLine:
      return lines[k].F0 - cutofffreq + (mirroring == MirroringType::Manual ? - shift : shift);
    case CutoffType::None:
      return std::numeric_limits<Numeric>::lowest();
    case CutoffType::FINAL: break;
  }
}
#pragma GCC diagnostic pop
 
bifstream& Lines::read(bifstream& is) {
  for (auto& line: lines)
    line.read(is);
  return is;
}
 
bofstream& Lines::write(bofstream& os) const {
  for (auto& line: lines)
    line.write(os);
  return os;
}
 
QuantumIdentifier Lines::QuantumIdentityOfLine(Index k) const noexcept {
  QuantumIdentifier qid = quantumidentity;
  for (auto& a: lines[k].localquanta.val) qid.val.set(a);
  return qid;
}
 
void Lines::MakeLineShapeModelCommon() {
  const Index n = NumLines();
  const Index m = NumBroadeners();
  if (not n) return;
  
  for (Index j=0; j<m; j++) {
    for (Index i=0; i<LineShape::nVars; i++) {
      
      // Find a common type (same or none) or throw if there is none
      LineShape::TemperatureModel t=LineShape::TemperatureModel::None;
      for (auto& line : lines) {
        if (auto& data = line.lineshape[j].Data()[i]; not LineShape::modelparameterEmpty(data)) {
          if (t == LineShape::TemperatureModel::None) t = data.type;
          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)
        }
      }
      
      // Set the common type
      for (auto& line : lines) {
        if (auto& data = line.lineshape[j].Data()[i]; data.type not_eq t) {
          data = LineShape::modelparameterGetEmpty(t);
        }
      }
    }
  }
}
} // namespace Absorption