absorptionlines.h

Go to the documentation of this file.

/* Copyright (C) 2019
 Richard Larsson <larsson@mps.mpg.de>
 
 
 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. */
 
#ifndef absorptionlines_h
#define absorptionlines_h
 
#include <vector>
#include "bifstream.h"
#include "bofstream.h"
#include "lineshapemodel.h"
#include "matpack.h"
#include "quantum.h"
#include "zeemandata.h"
 
namespace Absorption {
enum class MirroringType : Index {
  None,             // No mirroring
  Lorentz,          // Mirror, but use Lorentz line shape
  SameAsLineShape,  // Mirror using the same line shape
  Manual,           // Mirror by having a line in the array of line record with negative F0
};  // MirroringType
 
inline MirroringType string2mirroringtype(const String& in) {
  if (in == "None")
    return MirroringType::None;
  else if (in == "Lorentz")
    return MirroringType::Lorentz;
  else if (in == "Same")
    return MirroringType::SameAsLineShape;
  else if (in == "Manual")
    return MirroringType::Manual;
  else
    throw std::runtime_error("Cannot recognize the mirroring type");
}
 
inline String mirroringtype2string(MirroringType in) {
  if (in == MirroringType::None)
    return "None";
  else if (in == MirroringType::Lorentz)
    return "Lorentz";
  else if (in == MirroringType::SameAsLineShape)
    return "Same";
  else if (in == MirroringType::Manual)
    return "Manual";
  std::terminate();
}
 
inline String mirroringtype2metadatastring(MirroringType in) {
  if (in == MirroringType::None)
    return "These lines are not mirrored at 0 Hz.\n";
  else if (in == MirroringType::Lorentz)
    return "These lines are mirrored around 0 Hz using the Lorentz line shape.\n";
  else if (in == MirroringType::SameAsLineShape)
    return "These line are mirrored around 0 Hz using the original line shape.\n";
  else if (in == MirroringType::Manual)
    return "There are manual line entries in the catalog to mirror this line.\n";
  std::terminate();
}
 
enum class NormalizationType : Index {
  None,  // Do not renormalize the line shape
  VVH,   // Renormalize with Van Vleck and Huber specifications
  VVW,   // Renormalize with Van Vleck and Weiskopf specifications
  RosenkranzQuadratic,  // Renormalize using Rosenkranz's quadratic specifications
};  // LineNormalizationType
 
inline NormalizationType string2normalizationtype(const String& in) {
  if (in == "None")
    return NormalizationType::None;
  else if (in == "VVH")
    return NormalizationType::VVH;
  else if (in == "VVW")
    return NormalizationType::VVW;
  else if (in == "RQ")
    return NormalizationType::RosenkranzQuadratic;
  else
    throw std::runtime_error("Cannot recognize the normalization type");
}
 
inline String normalizationtype2string(NormalizationType in) {
  if (in == NormalizationType::None)
    return "None";
  else if (in == NormalizationType::VVH)
    return "VVH";
  else if (in == NormalizationType::VVW)
    return "VVW";
  else if (in == NormalizationType::RosenkranzQuadratic)
    return "RQ";
  std::terminate();
}
 
inline String normalizationtype2metadatastring(NormalizationType in) {
  if (in == NormalizationType::None)
    return "No re-normalization in the far wing will be applied.\n";
  else if (in == NormalizationType::VVH)
    return "van Vleck and Huber far-wing renormalization will be applied, "
      "i.e. F ~ (f tanh(hf/2kT))/(f0 tanh(hf0/2kT))\n";
  else if (in == NormalizationType::VVW)
    return "van Vleck and Weisskopf far-wing renormalization will be applied, "
      "i.e. F ~ (f/f0)^2\n";
  else if (in == NormalizationType::RosenkranzQuadratic)
    return "Rosenkranz quadratic far-wing renormalization will be applied, "
      "i.e. F ~ hf0/2kT sinh(hf0/2kT) (f/f0)^2\n";
  std::terminate();
}
 
enum class PopulationType : Index {
  ByLTE,                          // Assume line is in LTE
  ByNLTEVibrationalTemperatures,  // Assume line is in NLTE described by vibrational temperatures
  ByNLTEPopulationDistribution,   // Assume line is in NLTE and the upper-to-lower ratio is known
  ByHITRANRosenkranzRelmat,       // Assume line needs to compute relaxation matrix to derive HITRAN Y-coefficients
  ByHITRANFullRelmat,             // Assume line needs to compute and directly use the relaxation matrix
};  // PopulationType
 
inline PopulationType string2populationtype(const String& in) {
  if (in == "LTE")
    return PopulationType::ByLTE;
  else if (in == "ByHITRANRosenkranzRelmat")
    return PopulationType::ByHITRANRosenkranzRelmat;
  else if (in == "ByHITRANFullRelmat")
    return PopulationType::ByHITRANFullRelmat;
  else if (in == "NLTE-VibrationalTemperatures")
    return PopulationType::ByNLTEVibrationalTemperatures;
  else if (in == "NLTE")
    return PopulationType::ByNLTEPopulationDistribution;
  else
    throw std::runtime_error("Cannot recognize the population type");
}
 
inline String populationtype2string(PopulationType in) {
  switch (in) {
    case PopulationType::ByLTE:
      return "LTE";
    case PopulationType::ByHITRANFullRelmat:
      return "ByHITRANFullRelmat";
    case PopulationType::ByHITRANRosenkranzRelmat:
      return "ByHITRANRosenkranzRelmat";
    case PopulationType::ByNLTEVibrationalTemperatures:
      return "NLTE-VibrationalTemperatures";
    case PopulationType::ByNLTEPopulationDistribution:
      return "NLTE";
  } std::terminate();
}
 
inline String populationtype2metadatastring(PopulationType in) {
  switch (in) {
    case PopulationType::ByLTE:
      return "The lines are considered as in pure LTE.\n";
    case PopulationType::ByHITRANFullRelmat:
      return "The lines requires relaxation matrix calculations in LTE - HITRAN full method.\n";
    case PopulationType::ByHITRANRosenkranzRelmat:
      return "The lines requires Relaxation matrix calculations in LTE - HITRAN Rosenkranz method.\n";
    case PopulationType::ByNLTEVibrationalTemperatures:
      return "The lines are considered as in NLTE by vibrational temperatures.\n";
    case PopulationType::ByNLTEPopulationDistribution:
      return "The lines are considered as in pure NLTE.\n";
  } std::terminate();
}
 
inline bool relaxationtype_relmat(PopulationType in) {
  return in == PopulationType::ByHITRANFullRelmat or
         in == PopulationType::ByHITRANRosenkranzRelmat;
}
 
enum class CutoffType : Index {
  None,                // No cutoff frequency at all
  LineByLineOffset,    // The cutoff frequency is at SingleLine::F0 plus the cutoff frequency
  BandFixedFrequency,  // The curoff frequency is the cutoff frequency for all SingleLine(s)
};  // LineCutoffType
 
inline CutoffType string2cutofftype(const String& in) {
  if (in == "None")
    return CutoffType::None;
  else if (in == "ByLine")
    return CutoffType::LineByLineOffset;
  else if (in == "ByBand")
    return CutoffType::BandFixedFrequency;
  else
    throw std::runtime_error("Cannot recognize the cutoff type");
}
 
inline String cutofftype2string(CutoffType in) {
  if (in == CutoffType::None)
    return "None";
  else if (in == CutoffType::LineByLineOffset)
    return "ByLine";
  else if (in == CutoffType::BandFixedFrequency)
    return "ByBand";
  std::terminate();
}
 
inline String cutofftype2metadatastring(CutoffType in, Numeric cutoff) {
  std::ostringstream os;
  if (in == CutoffType::None)
    os << "No cut-off will be applied.\n";
  else if (in == CutoffType::LineByLineOffset)
    os << "The lines will be cut-off " << cutoff << " Hz from the line center.\n";
  else if (in == CutoffType::BandFixedFrequency)
    os << "All lines are cut-off at " << cutoff << " Hz.\n";
  return os.str();
}
 
class SingleLine {
private:
  Numeric mF0;
  
  Numeric mI0;
  
  Numeric mE0;
  
  Numeric mglow;
  
  Numeric mgupp;
  
  Numeric mA;
  
  Zeeman::Model mzeeman;
  
  LineShape::Model mlineshape;
  
  std::vector<Rational> mlowerquanta;
  
  std::vector<Rational> mupperquanta;
 
public:
  SingleLine(Numeric F0=0,
             Numeric I0=0,
             Numeric E0=0,
             Numeric glow=0,
             Numeric gupp=0,
             Numeric A=0,
             Zeeman::Model zeeman=Zeeman::Model(),
             const LineShape::Model& lineshape=LineShape::Model(),
             const std::vector<Rational>& lowerquanta={},
             const std::vector<Rational>& upperquanta={}) :
             mF0(F0),
             mI0(I0),
             mE0(E0),
             mglow(glow),
             mgupp(gupp),
             mA(A),
             mzeeman(zeeman),
             mlineshape(lineshape),
             mlowerquanta(lowerquanta),
             mupperquanta(upperquanta) {}
  
  SingleLine(size_t nbroadeners, size_t nquanta, const LineShape::Model& metamodel) :
  mlineshape(metamodel), mlowerquanta(nquanta), mupperquanta(nquanta) {
    if(Index(nbroadeners) not_eq mlineshape.nelem())
      throw std::runtime_error("Mismatch between broadeners and model");
  }
  
  
  Index LineShapeElems() const noexcept {return mlineshape.nelem();}
  
  Index LowerQuantumElems() const noexcept {return mlowerquanta.size();}
  
  Index UpperQuantumElems() const noexcept {return mupperquanta.size();}
  
  
  Numeric F0() const noexcept {return mF0;}
  
  Numeric E0() const noexcept {return mE0;}
  
  Numeric I0() const noexcept {return mI0;}
  
  Numeric A() const noexcept {return mA;}
  
  Numeric g_low() const noexcept {return mglow;}
  
  Numeric g_upp() const noexcept {return mgupp;}
  
  Zeeman::Model Zeeman() const noexcept {return mzeeman;}
  
  const LineShape::Model& LineShape() const noexcept {return mlineshape;}
  
  const std::vector<Rational>& LowerQuantumNumbers() const noexcept {return mlowerquanta;}
  
  const std::vector<Rational>& UpperQuantumNumbers() const noexcept {return mupperquanta;}
  
  
  Numeric& F0() noexcept {return mF0;}
  
  Numeric& E0() noexcept {return mE0;}
  
  Numeric& I0() noexcept {return mI0;}
  
  Numeric& A() noexcept {return mA;}
  
  Numeric& g_low() noexcept {return mglow;}
  
  Numeric& g_upp() noexcept {return mgupp;}
  
  Zeeman::Model& Zeeman() noexcept {return mzeeman;}
  
  LineShape::Model& LineShape() noexcept {return mlineshape;}
  
  std::vector<Rational>& LowerQuantumNumbers() noexcept {return mlowerquanta;}
  
  std::vector<Rational>& UpperQuantumNumbers() noexcept {return mupperquanta;}
  
  
  void F0(Numeric x) noexcept {mF0 = x;}
  
  void E0(Numeric x) noexcept {mE0 = x;}
  
  void I0(Numeric x) noexcept {mI0 = x;}
  
  void A(Numeric x) noexcept {mA = x;}
  
  void g_low(Numeric x) noexcept {mglow = x;}
  
  void g_upp(Numeric x) noexcept {mgupp = x;}
  
  
  Rational LowerQuantumNumber(size_t i) const noexcept {return mlowerquanta[i];}
  
  Rational UpperQuantumNumber(size_t i) const noexcept {return mupperquanta[i];}
  
  Rational& LowerQuantumNumber(size_t i) noexcept {return mlowerquanta[i];}
  
  Rational& UpperQuantumNumber(size_t i) noexcept {return mupperquanta[i];}
  
  bool SameQuantumNumbers(const SingleLine& sl) const noexcept;
  
  
  void SetAutomaticZeeman(QuantumIdentifier qid, const std::vector<QuantumNumberType>& keys) {
    for(size_t i=0; i<keys.size(); i++) {
      qid.LowerQuantumNumber(keys[i]) = mlowerquanta[i];
      qid.UpperQuantumNumber(keys[i]) = mupperquanta[i];
    }
    
    mzeeman = Zeeman::Model(qid);
  }
  
  void SetLineMixing2SecondOrderData(const Vector& d) {
    mlineshape.SetLineMixingModel(
      LineShape::LegacyLineMixingData::vector2modellm(
        d, LineShape::LegacyLineMixingData::TypeLM::LM_2NDORDER)
      .Data()[0]);
  }
  
  void 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 : mlineshape.Data()) {
      sm.Y() = Y;
      sm.G() = G;
    }
  }
  
  bifstream& read(bifstream& bif) {
    bif >> mF0 >> mI0 >> mE0 >> mglow >> mgupp >> mA >> mzeeman;
    
    mlineshape.read(bif);
    
    for (auto& rat: mlowerquanta) rat.read(bif);
    
    for (auto& rat: mupperquanta) rat.read(bif);
    
    return bif;
  }
  
  bofstream& write(bofstream& bof) const {
    bof << mF0 << mI0 << mE0 << mglow << mgupp << mA << mzeeman;
    
    mlineshape.write(bof);
    
    for (auto& rat: mlowerquanta) rat.write(bof);
    
    for (auto& rat: mupperquanta) rat.write(bof);
    
    return bof;
  }
};  // SingleLine
 
std::ostream& operator<<(std::ostream&, const SingleLine&);
 
std::istream& operator>>(std::istream&, SingleLine&);
 
struct SingleLineExternal {
  bool bad=true;
  bool selfbroadening=false;
  bool bathbroadening=false;
  CutoffType cutoff=CutoffType::None;
  MirroringType mirroring=MirroringType::None;
  PopulationType population=PopulationType::ByLTE;
  NormalizationType normalization=NormalizationType::None;
  LineShape::Type lineshapetype=LineShape::Type::DP;
  Numeric T0=0;
  Numeric cutofffreq=0;
  Numeric linemixinglimit=-1;
  QuantumIdentifier quantumidentity=QuantumIdentifier(QuantumIdentifier::TRANSITION, -1, -1);
  ArrayOfSpeciesTag species;
  SingleLine line;
};
 
class Lines {
private:
  bool mselfbroadening;
  
  bool mbathbroadening;
  
  CutoffType mcutoff;
  
  MirroringType mmirroring;
  
  PopulationType mpopulation;
  
  NormalizationType mnormalization;
 
  LineShape::Type mlineshapetype;
  
  Numeric mT0;
  
  Numeric mcutofffreq;
  
  Numeric mlinemixinglimit;
  
  QuantumIdentifier mquantumidentity;
  
  std::vector<QuantumNumberType> mlocalquanta;
  
  ArrayOfSpeciesTag mbroadeningspecies;
  
  std::vector<SingleLine> mlines;
  
public:
  Lines(bool selfbroadening=false,
        bool bathbroadening=false,
        CutoffType cutoff=CutoffType::None,
        MirroringType mirroring=MirroringType::None,
        PopulationType population=PopulationType::ByLTE,
        NormalizationType normalization=NormalizationType::None,
        LineShape::Type lineshapetype=LineShape::Type::DP,
        Numeric T0=296,
        Numeric cutofffreq=-1,
        Numeric linemixinglimit=-1,
        const QuantumIdentifier& quantumidentity=QuantumIdentifier(),
        const std::vector<QuantumNumberType>& localquanta={},
        const ArrayOfSpeciesTag& broadeningspecies={},
        const std::vector<SingleLine>& lines={}) :
        mselfbroadening(selfbroadening),
        mbathbroadening(bathbroadening),
        mcutoff(cutoff),
        mmirroring(mirroring),
        mpopulation(population),
        mnormalization(normalization),
        mlineshapetype(lineshapetype),
        mT0(T0),
        mcutofffreq(cutofffreq),
        mlinemixinglimit(linemixinglimit),
        mquantumidentity(quantumidentity),
        mlocalquanta(localquanta),
        mbroadeningspecies(broadeningspecies),
        mlines(lines) {};
  
  Lines(bool selfbroadening,
        bool bathbroadening,
        size_t nlines,
        CutoffType cutoff,
        MirroringType mirroring,
        PopulationType population,
        NormalizationType normalization,
        LineShape::Type lineshapetype,
        Numeric T0,
        Numeric cutofffreq,
        Numeric linemixinglimit,
        const QuantumIdentifier& quantumidentity,
        const std::vector<QuantumNumberType>& localquanta,
        const ArrayOfSpeciesTag& broadeningspecies,
        const LineShape::Model& metamodel) :
        mselfbroadening(selfbroadening),
        mbathbroadening(bathbroadening),
        mcutoff(cutoff),
        mmirroring(mirroring),
        mpopulation(population),
        mnormalization(normalization),
        mlineshapetype(lineshapetype),
        mT0(T0),
        mcutofffreq(cutofffreq),
        mlinemixinglimit(linemixinglimit),
        mquantumidentity(quantumidentity),
        mlocalquanta(localquanta),
        mbroadeningspecies(broadeningspecies),
        mlines(nlines,
               SingleLine(broadeningspecies.size(),
               localquanta.size(), metamodel)) {};
  
  void AppendSingleLine(SingleLine&& sl) {
    if(NumLocalQuanta() not_eq sl.LowerQuantumElems() or
       NumLocalQuanta() not_eq sl.UpperQuantumElems())
      throw std::runtime_error("Error calling appending function, bad size of quantum numbers");
    
    if(NumLines() not_eq 0 and 
       sl.LineShapeElems() not_eq mlines[0].LineShapeElems())
      throw std::runtime_error("Error calling appending function, bad size of broadening species");
    
    mlines.push_back(std::move(sl));
  }
  
  void AppendSingleLine(const SingleLine& sl) {
    if(NumLocalQuanta() not_eq sl.LowerQuantumElems() or
       NumLocalQuanta() not_eq sl.UpperQuantumElems())
      throw std::runtime_error("Error calling appending function, bad size of quantum numbers");
    
    if(NumLines() not_eq 0 and 
       sl.LineShapeElems() not_eq mlines[0].LineShapeElems())
      throw std::runtime_error("Error calling appending function, bad size of broadening species");
    
    mlines.push_back(sl);
  }
  
  bool MatchWithExternal(const SingleLineExternal& sle, const QuantumIdentifier& quantumidentity) const noexcept {
    if(sle.bad)
      return false;
    else if(sle.selfbroadening not_eq mselfbroadening)
      return false;
    else if(sle.bathbroadening not_eq mbathbroadening)
      return false;
    else if(sle.cutoff not_eq mcutoff)
      return false;
    else if(sle.mirroring not_eq mmirroring)
      return false;
    else if(sle.population not_eq mpopulation)
      return false;
    else if(sle.normalization not_eq mnormalization)
      return false;
    else if(sle.lineshapetype not_eq mlineshapetype)
      return false;
    else if(sle.T0 not_eq mT0)
      return false;
    else if(sle.cutofffreq not_eq mcutofffreq)
      return false;
    else if(sle.linemixinglimit not_eq mlinemixinglimit)
      return false;
    else if(quantumidentity not_eq mquantumidentity)
      return false;
    else if(not std::equal(sle.species.cbegin(), sle.species.cend(), mbroadeningspecies.cbegin(), mbroadeningspecies.cend()))
      return false;
    else if(NumLines() not_eq 0 and not sle.line.LineShape().Match(mlines[0].LineShape()))
      return false;
    else
      return true;
  }
  
  bool Match(const Lines& l) const noexcept {
    if(l.mselfbroadening not_eq mselfbroadening)
      return false;
    else if(l.mbathbroadening not_eq mbathbroadening)
      return false;
    else if(l.mcutoff not_eq mcutoff)
      return false;
    else if(l.mmirroring not_eq mmirroring)
      return false;
    else if(l.mpopulation not_eq mpopulation)
      return false;
    else if(l.mnormalization not_eq mnormalization)
      return false;
    else if(l.mlineshapetype not_eq mlineshapetype)
      return false;
    else if(l.mT0 not_eq mT0)
      return false;
    else if(l.mcutofffreq not_eq mcutofffreq)
      return false;
    else if(l.mlinemixinglimit not_eq mlinemixinglimit)
      return false;
    else if(l.mquantumidentity not_eq mquantumidentity)
      return false;
    else if(not std::equal(l.mbroadeningspecies.cbegin(), l.mbroadeningspecies.cend(), mbroadeningspecies.cbegin(), mbroadeningspecies.cend()))
      return false;
    else if(not std::equal(l.mlocalquanta.cbegin(), l.mlocalquanta.cend(), mlocalquanta.cbegin(), mlocalquanta.cend()))
      return false;
    else if(NumLines() not_eq 0 and l.NumLines() not_eq 0 and not l.mlines[0].LineShape().Match(mlines[0].LineShape()))
      return false;
    else
      return true;
  }
  
  void sort_by_frequency() {
    std::sort(mlines.begin(), mlines.end(),
              [](const SingleLine& a, const SingleLine& b){return a.F0() < b.F0();});
  }
  
  void sort_by_einstein() {
    std::sort(mlines.begin(), mlines.end(),
              [](const SingleLine& a, const SingleLine& b){return a.A() < b.A();});
  }
  
  void truncate_global_quantum_numbers() {
    mquantumidentity.SetTransition(QuantumNumbers(), QuantumNumbers());
  }
  
  String SpeciesName() const noexcept;
  
  String UpperQuantumNumbers() const noexcept;
  
  String LowerQuantumNumbers() const noexcept;
  
  String LineShapeMetaData() const noexcept {
    return NumLines() ?
      LineShape::ModelShape2MetaData(mlines[0].LineShape()) :
      "";
  }
  
  Index Species() const noexcept {return mquantumidentity.Species();}
  
  Index Isotopologue() const noexcept {return mquantumidentity.Isotopologue();}
  
  Index NumLines() const noexcept {return Index(mlines.size());}
  
  const std::vector<SingleLine>& AllLines() const noexcept {return mlines;}
  
  std::vector<SingleLine>& AllLines() noexcept {return mlines;}
  
  Index NumBroadeners() const noexcept {return Index(mbroadeningspecies.nelem());}
  
  Index NumLocalQuanta() const noexcept {return Index(mlocalquanta.size());}
  
  void RemoveUnusedLocalQuantums();
  
  void RemoveLocalQuantum(size_t);
  
  Rational LowerQuantumNumber(size_t k, QuantumNumberType qnt) const noexcept;
  
  Rational UpperQuantumNumber(size_t k, QuantumNumberType qnt) const noexcept;
  
  Rational& LowerQuantumNumber(size_t k, QuantumNumberType qnt) noexcept;
  
  Rational& UpperQuantumNumber(size_t k, QuantumNumberType qnt) noexcept;
  
  Index ZeemanCount(size_t k, Zeeman::Polarization type) const noexcept {
    if (UpperQuantumNumber(k, QuantumNumberType::F).isDefined() and LowerQuantumNumber(k, QuantumNumberType::F).isDefined()) {
      return Zeeman::nelem(UpperQuantumNumber(k, QuantumNumberType::F),
                           LowerQuantumNumber(k, QuantumNumberType::F),
                           type);
    } else {
      return Zeeman::nelem(UpperQuantumNumber(k, QuantumNumberType::J),
                           LowerQuantumNumber(k, QuantumNumberType::J),
                           type);
    }
  }
  
  Numeric ZeemanStrength(size_t k, Zeeman::Polarization type, Index i) const noexcept {
    if (UpperQuantumNumber(k, QuantumNumberType::F).isDefined() and LowerQuantumNumber(k, QuantumNumberType::F).isDefined()) {
      return mlines[k].Zeeman().Strength(UpperQuantumNumber(k, QuantumNumberType::F),
                                         LowerQuantumNumber(k, QuantumNumberType::F),
                                         type, i);
    } else {
      return mlines[k].Zeeman().Strength(UpperQuantumNumber(k, QuantumNumberType::J),
                                         LowerQuantumNumber(k, QuantumNumberType::J),
                                         type, i);
    }
  }
  
  Numeric ZeemanSplitting(size_t k, Zeeman::Polarization type, Index i) const noexcept {
    if (UpperQuantumNumber(k, QuantumNumberType::F).isDefined() and LowerQuantumNumber(k, QuantumNumberType::F).isDefined()) {
      return mlines[k].Zeeman().Splitting(UpperQuantumNumber(k, QuantumNumberType::F),
                                          LowerQuantumNumber(k, QuantumNumberType::F),
                                          type, i);
    } else {
      return mlines[k].Zeeman().Splitting(UpperQuantumNumber(k, QuantumNumberType::J),
                                          LowerQuantumNumber(k, QuantumNumberType::J),
                                          type, i);
    }
  }
  
  void SetAutomaticZeeman() noexcept {
    for(auto& line: mlines)
      line.SetAutomaticZeeman(mquantumidentity, mlocalquanta);
  }
  
  Numeric F0(size_t k) const noexcept {return mlines[k].F0();}
  
  Numeric& F0(size_t k) noexcept {return mlines[k].F0();}
  
  Numeric F_mean() const noexcept {
    const Numeric val = std::inner_product(mlines.cbegin(), mlines.cend(),
                                          mlines.cbegin(), 0.0, std::plus<Numeric>(),
                                          [](const auto& a, const auto& b){return a.F0() * b.I0();});
    const Numeric div = std::accumulate(mlines.cbegin(), mlines.cend(), 0.0,
                                        [](const auto& a, const auto& b){return a + b.I0();});
    return  val / div;
  }
  
  Numeric F_mean(const ConstVectorView wgts) const noexcept {
    const Numeric val = std::inner_product(mlines.cbegin(), mlines.cend(),
                                           wgts.begin(), 0.0, std::plus<Numeric>(),
                                           [](const auto& a, const auto& b){return a.F0() * b;});
    const Numeric div = wgts.sum();
    return  val / div;
  }
  
  Numeric E0(size_t k) const noexcept {return mlines[k].E0();}
  
  Numeric& E0(size_t k) noexcept {return mlines[k].E0();}
  
  Numeric I0(size_t k) const noexcept {return mlines[k].I0();}
  
  Numeric& I0(size_t k) noexcept {return mlines[k].I0();}
  
  Numeric A(size_t k) const noexcept {return mlines[k].A();}
  
  Numeric& A(size_t k) noexcept {return mlines[k].A();}
  
  Numeric g_low(size_t k) const noexcept {return mlines[k].g_low();}
  
  Numeric& g_low(size_t k) noexcept {return mlines[k].g_low();}
  
  Numeric g_upp(size_t k) const noexcept {return mlines[k].g_upp();}
  
  Numeric& g_upp(size_t k) noexcept {return mlines[k].g_upp();}
  
  MirroringType Mirroring() const noexcept {return mmirroring;}
  
  void Mirroring(MirroringType x) noexcept {mmirroring = x;}
  
  static bool validIndexForMirroring(Index x) noexcept {
    constexpr auto keys = stdarrayify(Index(MirroringType::None), MirroringType::None, MirroringType::Lorentz, MirroringType::SameAsLineShape, MirroringType::Manual);
    return std::any_of(keys.cbegin(), keys.cend(), [x](auto y){return x == y;});
  }
  
  static MirroringType string2Mirroring(const String& in) noexcept {
    if (in == "None")
      return MirroringType::None;
    else if (in == "Lorentz")
      return MirroringType::Lorentz;
    else if (in == "Same")
      return MirroringType::SameAsLineShape;
    else if (in == "Manual")
      return MirroringType::Manual;
    else
      return MirroringType(-1);
  }
  
  NormalizationType Normalization() const noexcept {return mnormalization;}
  
  void Normalization(NormalizationType x) noexcept {mnormalization = x;}
  
  static bool validIndexForNormalization(Index x) noexcept {
    constexpr auto keys = stdarrayify(Index(NormalizationType::None), NormalizationType::VVH, NormalizationType::VVW, NormalizationType::RosenkranzQuadratic);
    return std::any_of(keys.cbegin(), keys.cend(), [x](auto y){return x == y;});
  }
  
  static NormalizationType string2Normalization(const String& in) noexcept {
    if (in == "None")
      return NormalizationType::None;
    else if (in == "VVH")
      return NormalizationType::VVH;
    else if (in == "VVW")
      return NormalizationType::VVW;
    else if (in == "RQ")
      return NormalizationType::RosenkranzQuadratic;
    else
      return NormalizationType(-1);
  }
  
  CutoffType Cutoff() const noexcept {return mcutoff;}
  
  void Cutoff(CutoffType x) noexcept {mcutoff = x;}
  
  static bool validIndexForCutoff(Index x) noexcept {
    constexpr auto keys = stdarrayify(Index(CutoffType::None), CutoffType::LineByLineOffset, CutoffType::BandFixedFrequency);
    return std::any_of(keys.cbegin(), keys.cend(), [x](auto y){return x == y;});
  }
  
  static CutoffType string2Cutoff(const String& in) noexcept {
    if (in == "None")
      return CutoffType::None;
    else if (in == "ByLine")
      return CutoffType::LineByLineOffset;
    else if (in == "ByBand")
      return CutoffType::BandFixedFrequency;
    else
      return CutoffType(-1);
  }
  
  PopulationType Population() const noexcept {return mpopulation;}
  
  void Population(PopulationType x) noexcept {mpopulation = x;}
  
  static bool validIndexForPopulation(Index x) noexcept {
    constexpr auto keys = stdarrayify(Index(PopulationType::ByLTE), PopulationType::ByHITRANFullRelmat, PopulationType::ByHITRANRosenkranzRelmat, PopulationType::ByNLTEVibrationalTemperatures, PopulationType::ByNLTEPopulationDistribution);
    return std::any_of(keys.cbegin(), keys.cend(), [x](auto y){return x == y;});
  }
  
  static PopulationType string2Population(const String& in) noexcept {
    if (in == "LTE")
      return PopulationType::ByLTE;
    else if (in == "ByHITRANFullRelmat")
      return PopulationType::ByHITRANFullRelmat;
    else if (in == "ByHITRANRosenkranzRelmat")
      return PopulationType::ByHITRANRosenkranzRelmat;
    else if (in == "NLTE-VibrationalTemperatures")
      return PopulationType::ByNLTEVibrationalTemperatures;
    else if (in == "NLTE")
      return PopulationType::ByNLTEPopulationDistribution;
    else
      return PopulationType(-1);
  }
  
  LineShape::Type LineShapeType() const noexcept {return mlineshapetype;}
  
  void LineShapeType(LineShape::Type x) noexcept {mlineshapetype = x;}
  
  static bool validIndexForLineShapeType(Index x) noexcept {
    constexpr auto keys = stdarrayify(Index(LineShape::Type::DP), LineShape::Type::LP, LineShape::Type::VP, LineShape::Type::SDVP, LineShape::Type::HTP);
    return std::any_of(keys.cbegin(), keys.cend(), [x](auto y){return x == y;});
  }
  
  static LineShape::Type string2LineShapeType(const String& type) noexcept {
    if (type == "DP")
      return LineShape::Type::DP;
    else if (type == String("LP"))
      return LineShape::Type::LP;
    else if (type == String("VP"))
      return LineShape::Type::VP;
    else if (type == String("SDVP"))
      return LineShape::Type::SDVP;
    else if (type == String("HTP"))
      return LineShape::Type::HTP;
    else
      return LineShape::Type(-1);
  }
  
  bool DoLineMixing(Numeric P) const noexcept {
    return mlinemixinglimit < 0 ? true : mlinemixinglimit > P;
  }
  
  LineShape::Output ShapeParameters(size_t k, Numeric T, Numeric P, const Vector& vmrs) const noexcept;
  
  LineShape::Output ShapeParameters_dT(size_t k, Numeric T, Numeric P, const Vector& vmrs) const noexcept;
  
  Index LineShapePos(const Index& spec) const noexcept;
  
  Index LineShapePos(const QuantumIdentifier& qid) const noexcept {
    return LineShapePos(qid.Species());
  }
  
  LineShape::Output ShapeParameters_dVMR(size_t k, Numeric T, Numeric P,
                                         const QuantumIdentifier& vmr_qid) const noexcept;
  
  Numeric ShapeParameter_dInternal(size_t k, Numeric T, Numeric P,
                                   const Vector& vmrs,
                                   const RetrievalQuantity& derivative) const noexcept;
  
  Numeric CutoffFreq(size_t k) const noexcept {
    switch(mcutoff) {
      case CutoffType::LineByLineOffset:
        return F0(k) + mcutofffreq;
      case CutoffType::BandFixedFrequency:
        return mcutofffreq;
      case CutoffType::None:
        return std::numeric_limits<Numeric>::max();
    }
    std::terminate();
  }
  
  Numeric CutoffFreqMinus(size_t k, Numeric fmean) const noexcept {
    switch(mcutoff) {
      case CutoffType::LineByLineOffset:
        return F0(k) - mcutofffreq;
      case CutoffType::BandFixedFrequency:
        return mcutofffreq - 2*fmean;
      case CutoffType::None:
        return std::numeric_limits<Numeric>::lowest();
    }
    std::terminate();
  }
  
  Numeric T0() const noexcept {
    return mT0;
  }
  
  void T0(Numeric x) noexcept {
    mT0 = x;
  }
  
  Numeric CutoffFreqValue() const noexcept {
    return mcutofffreq;
  }
  
  void CutoffFreqValue(Numeric x) noexcept {
    mcutofffreq = x;
  }
  
  Numeric LinemixingLimit() const noexcept {
    return mlinemixinglimit;
  }
  
  void LinemixingLimit(Numeric x) noexcept {
    mlinemixinglimit = x;
  }
  
  const std::vector<QuantumNumberType>& LocalQuanta() const noexcept {
    return mlocalquanta;
  }
  
  std::vector<QuantumNumberType>& LocalQuanta() noexcept {
    return mlocalquanta;
  }
  
  const ArrayOfSpeciesTag& BroadeningSpecies() const noexcept {
    return mbroadeningspecies;
  }
  
  ArrayOfSpeciesTag& BroadeningSpecies() noexcept {
    return mbroadeningspecies;
  }
  
  bool Self() const noexcept {
    return mselfbroadening;
  }
  
  void Self(bool x) noexcept {
    mselfbroadening = x;
  }
  
  bool Bath() const noexcept {
    return mbathbroadening;
  }
  
  void Bath(bool x) noexcept {
    mbathbroadening = x;
  }
  
  const QuantumIdentifier& QuantumIdentity() const noexcept {
    return mquantumidentity;
  }
  
  QuantumIdentifier& QuantumIdentity() noexcept {
    return mquantumidentity;
  }
  
  QuantumIdentifier QuantumIdentityOfLine(Index k) const noexcept {
    QuantumIdentifier qid_copy(mquantumidentity);
    for (size_t i=0; i<mlocalquanta.size(); i++) {
      qid_copy.UpperQuantumNumber(mlocalquanta[i]) = mlines[k].UpperQuantumNumber(i);
      qid_copy.LowerQuantumNumber(mlocalquanta[i]) = mlines[k].LowerQuantumNumber(i);
    }
    return qid_copy;
  }
  
  String MetaData() const;
  
  void RemoveLine(Index) noexcept;
  
  SingleLine PopLine(Index) noexcept;
  
  SingleLine& Line(Index) noexcept;
  
  const SingleLine& Line(Index) const noexcept;
  
  void ReverseLines() noexcept;
  
  Numeric SpeciesMass() const noexcept;
  
  Vector BroadeningSpeciesVMR(const ConstVectorView, const ArrayOfArrayOfSpeciesTag&) const;
  
  Numeric SelfVMR(const ConstVectorView, const ArrayOfArrayOfSpeciesTag&) const;
  
  bifstream& read(bifstream& is) {
    for (auto& line: mlines)
      line.read(is);
    return is;
  }
  
  bofstream& write(bofstream& os) const {
    for (auto& line: mlines)
      line.write(os);
    return os;
  }
  
  bool OK() const noexcept;
};  // Lines
 
std::ostream& operator<<(std::ostream&, const Lines&);
std::istream& operator>>(std::istream&, Lines&);
 
SingleLineExternal ReadFromArtscat3Stream(istream& is);
 
SingleLineExternal ReadFromArtscat4Stream(istream& is);
 
SingleLineExternal ReadFromArtscat5Stream(istream& is);
 
SingleLineExternal ReadFromLBLRTMStream(istream& is);
 
SingleLineExternal ReadFromHitran2004Stream(istream& is);
 
SingleLineExternal ReadFromHitranOnlineStream(istream& is);
 
SingleLineExternal ReadFromHitran2001Stream(istream& is);
 
SingleLineExternal ReadFromMytran2Stream(istream& is);
 
SingleLineExternal ReadFromJplStream(istream& is);
 
std::vector<Lines> split_list_of_external_lines(std::vector<SingleLineExternal>& external_lines,
                                                const std::vector<QuantumNumberType>& localquantas={},
                                                const std::vector<QuantumNumberType>& globalquantas={});
 
Lines createEmptyCopy(const Lines& al) noexcept;
 
bool line_in_id(const Lines& band, const QuantumIdentifier& id, size_t line_index);
 
bool line_is_id(const Lines& band, const QuantumIdentifier& id, size_t line_index);
 
bool id_in_line(const Lines& band, const QuantumIdentifier& id, size_t line_index);
 
bool line_upper_in_id(const Lines& band, const QuantumIdentifier& id, size_t line_index);
 
bool line_lower_in_id(const Lines& band, const QuantumIdentifier& id, size_t line_index);
 
bool id_in_line_upper(const Lines& band, const QuantumIdentifier& id, size_t line_index);
 
bool id_in_line_lower(const Lines& band, const QuantumIdentifier& id, size_t line_index);
 
inline Index nelem(const Lines& l) {return l.NumLines();}
 
inline Index nelem(const Array<Lines>& l) {Index n=0; for (auto& x:l) n+=nelem(x); return n;}
 
inline Index nelem(const Array<Array<Lines>>& l) {Index n=0; for (auto& x:l) n+=nelem(x); return n;}
 
Numeric reduced_rovibrational_dipole(Rational Jf, Rational Ji, Rational lf, Rational li, Rational k = Rational(1));
 
Numeric reduced_magnetic_quadrapole(Rational Jf, Rational Ji, Rational N);
};  // Absorption
 
typedef Absorption::SingleLine AbsorptionSingleLine;
typedef Absorption::Lines AbsorptionLines;
typedef Array<AbsorptionLines> ArrayOfAbsorptionLines;
typedef Array<ArrayOfAbsorptionLines> ArrayOfArrayOfAbsorptionLines;
 
std::ostream& operator<<(std::ostream&, const ArrayOfAbsorptionLines&);
 
std::ostream& operator<<(std::ostream&, const ArrayOfArrayOfAbsorptionLines&);
 
#endif  // absorptionlines_h