absorptionlines.h

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. */
 
#ifndef absorptionlines_h
#define absorptionlines_h
 
#include "bifstream.h"
#include "bofstream.h"
#include "enums.h"
#include "jacobian.h"
#include "lineshapemodel.h"
#include "matpack.h"
#include "quantum_numbers.h"
#include "species_tags.h"
#include "zeemandata.h"
 
#include <utility>
#include <vector>
 
 
namespace Absorption {
ENUMCLASS(MirroringType, char,
  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
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wreturn-type"
constexpr std::string_view mirroringtype2metadatastring(MirroringType in) noexcept {
  switch (in) {
    case MirroringType::None:
      return "These lines are not mirrored at 0 Hz.\n";
    case MirroringType::Lorentz:
      return "These lines are mirrored around 0 Hz using the Lorentz line shape.\n";
    case MirroringType::SameAsLineShape:
      return "These line are mirrored around 0 Hz using the original line shape.\n";
    case MirroringType::Manual:
      return "There are manual line entries in the catalog to mirror this line.\n";
    case MirroringType::FINAL: break;
  }
}
#pragma GCC diagnostic pop
 
ENUMCLASS(NormalizationType, char,
  None,                // Do not renormalize the line shape
  VVH,                 // Renormalize with Van Vleck and Huber specifications
  VVW,                 // Renormalize with Van Vleck and Weiskopf specifications
  RQ,                  // Renormalize using Rosenkranz's quadratic specifications
  SFS                  // Renormalize using simple frequency scaling of the line strength
)  // NormalizationType
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wreturn-type"
constexpr std::string_view normalizationtype2metadatastring(NormalizationType in) {
  switch (in) {
    case NormalizationType::None:
      return "No re-normalization in the far wing will be applied.\n";
    case 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";
    case NormalizationType::VVW:
      return "van Vleck and Weisskopf far-wing renormalization will be applied, "
        "i.e. F ~ (f/f0)^2\n";
    case NormalizationType::RQ:
      return "Rosenkranz quadratic far-wing renormalization will be applied, "
        "i.e. F ~ hf0/2kT sinh(hf0/2kT) (f/f0)^2\n";
    case NormalizationType::SFS:
      return "Simple frequency scaling of the far-wings will be applied, "
        "i.e. F ~ (f / f0) * ((1 - exp(- hf / kT)) / (1 - exp(- hf0 / kT)))\n";
    case NormalizationType::FINAL: break;
  }
}
#pragma GCC diagnostic pop
 
ENUMCLASS(PopulationType, char,
  LTE,                            // Assume band is in LTE
  NLTE,                           // Assume band is in NLTE and the upper-to-lower ratio is known
  VibTemps,                       // Assume band is in NLTE described by vibrational temperatures and LTE at other levels
  ByHITRANRosenkranzRelmat,       // Assume band needs to compute relaxation matrix to derive HITRAN Y-coefficients
  ByHITRANFullRelmat,             // Assume band needs to compute and directly use the relaxation matrix according to HITRAN
  ByMakarovFullRelmat,            // Assume band needs to compute and directly use the relaxation matrix according to Makarov et al 2020
  ByRovibLinearDipoleLineMixing   // Assume band needs to compute and directly use the relaxation matrix according to Hartmann, Boulet, Robert, 2008, 1st edition
)  // PopulationType
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wreturn-type"
constexpr std::string_view populationtype2metadatastring(PopulationType in) {
  switch (in) {
    case PopulationType::LTE:
      return "The lines are considered as in pure LTE.\n";
    case PopulationType::ByMakarovFullRelmat:
      return "The lines requires relaxation matrix calculations in LTE - Makarov et al 2020 full method.\n";
    case PopulationType::ByRovibLinearDipoleLineMixing:
      return "The lines requires relaxation matrix calculations in LTE - Hartmann, Boulet, Robert, 2008, 1st edition method.\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::VibTemps:
      return "The lines are considered as in NLTE by vibrational temperatures.\n";
    case PopulationType::NLTE:
      return "The lines are considered as in pure NLTE.\n";
    case PopulationType::FINAL: return "There's an error";
  }
}
#pragma GCC diagnostic pop
 
constexpr bool relaxationtype_relmat(PopulationType in) noexcept {
  return in == PopulationType::ByHITRANFullRelmat or
         in == PopulationType::ByMakarovFullRelmat or
         in == PopulationType::ByHITRANRosenkranzRelmat or
         in == PopulationType::ByRovibLinearDipoleLineMixing;
}
 
ENUMCLASS(CutoffType, char,
  None,                             // No cutoff frequency at all
  ByLine                            // The cutoff frequency is at SingleLine::F0 plus the cutoff frequency plus the speed independent pressure shift
)  // CutoffType
 
String cutofftype2metadatastring(CutoffType in, Numeric cutoff);
 
struct SingleLine {
  Numeric F0{};
  
  Numeric I0{};
  
  Numeric E0{};
  
  Numeric glow{};
  
  Numeric gupp{};
  
  Numeric A{};
  
  Zeeman::Model zeeman{};
  
  LineShape::Model lineshape{};
  
  Quantum::Number::LocalState localquanta{};
 
  SingleLine(Numeric F0_=0,
             Numeric I0_=0,
             Numeric E0_=0,
             Numeric glow_=0,
             Numeric gupp_=0,
             Numeric A_=0,
             Zeeman::Model zeeman_=Zeeman::Model(),
             LineShape::Model lineshape_=LineShape::Model(),
             Quantum::Number::LocalState localquanta_={}) :
             F0(F0_),
             I0(I0_),
             E0(E0_),
             glow(glow_),
             gupp(gupp_),
             A(A_),
             zeeman(zeeman_),
             lineshape(std::move(lineshape_)),
             localquanta(std::move(localquanta_)) {}
  
  SingleLine(Quantum::Number::LocalState metaquanta, LineShape::Model metamodel) :
  lineshape(std::move(metamodel)), localquanta(std::move(metaquanta)) {}
  
  
  [[nodiscard]] Index LineShapeElems() const noexcept {return lineshape.nelem();}
  
  [[nodiscard]] Index LocalQuantumElems() const noexcept {return localquanta.val.nelem();}
  
  
  void SetAutomaticZeeman(QuantumIdentifier qid);
  
  void SetLineMixing2SecondOrderData(const Vector& d);
  
  void SetLineMixing2AER(const Vector& d);
  
  bifstream& read(bifstream& bif);
  
  bofstream& write(bofstream& bof) const;
 
  friend std::ostream& operator<<(std::ostream&, const SingleLine&);
 
  friend std::istream& operator>>(std::istream&, SingleLine&);
};  // SingleLine
 
struct SingleLineExternal {
  bool bad=true;
  bool selfbroadening=false;
  bool bathbroadening=false;
  CutoffType cutoff=CutoffType::None;
  MirroringType mirroring=MirroringType::None;
  PopulationType population=PopulationType::LTE;
  NormalizationType normalization=NormalizationType::None;
  LineShape::Type lineshapetype=LineShape::Type::DP;
  Numeric T0=0;
  Numeric cutofffreq=0;
  Numeric linemixinglimit=-1;
  QuantumIdentifier quantumidentity;
  ArrayOfSpecies species;
  SingleLine line;
};
 
struct Lines {
  static constexpr Index version = 2;
 
  bool selfbroadening;
  
  bool bathbroadening;
  
  CutoffType cutoff;
  
  MirroringType mirroring;
  
  PopulationType population;
  
  NormalizationType normalization;
 
  LineShape::Type lineshapetype;
  
  Numeric T0;
  
  Numeric cutofffreq;
  
  Numeric linemixinglimit;
  
  QuantumIdentifier quantumidentity;
  
  ArrayOfSpecies broadeningspecies;
  
  Array<SingleLine> lines;
  
  Lines(bool selfbroadening_=false,
        bool bathbroadening_=false,
        CutoffType cutoff_=CutoffType::None,
        MirroringType mirroring_=MirroringType::None,
        PopulationType population_=PopulationType::LTE,
        NormalizationType normalization_=NormalizationType::None,
        LineShape::Type lineshapetype_=LineShape::Type::DP,
        Numeric T0_=296,
        Numeric cutofffreq_=-1,
        Numeric linemixinglimit_=-1,
        QuantumIdentifier quantumidentity_=QuantumIdentifier(),
        ArrayOfSpecies broadeningspecies_={},
        Array<SingleLine> lines_={}) :
        selfbroadening(selfbroadening_),
        bathbroadening(bathbroadening_),
        cutoff(cutoff_),
        mirroring(mirroring_),
        population(population_),
        normalization(normalization_),
        lineshapetype(lineshapetype_),
        T0(T0_),
        cutofffreq(cutofffreq_),
        linemixinglimit(linemixinglimit_),
        quantumidentity(std::move(quantumidentity_)),
        broadeningspecies(std::move(broadeningspecies_)),
        lines(std::move(lines_)) {
    if (selfbroadening) broadeningspecies.front() = quantumidentity.Species();
    if (bathbroadening) broadeningspecies.back() = Species::Species::Bath;
  }
  
  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_,
        QuantumIdentifier  quantumidentity_,
        ArrayOfSpecies  broadeningspecies_,
        const Quantum::Number::LocalState& metalocalquanta,
        const LineShape::Model& metamodel) :
        selfbroadening(selfbroadening_),
        bathbroadening(bathbroadening_),
        cutoff(cutoff_),
        mirroring(mirroring_),
        population(population_),
        normalization(normalization_),
        lineshapetype(lineshapetype_),
        T0(T0_),
        cutofffreq(cutofffreq_),
        linemixinglimit(linemixinglimit_),
        quantumidentity(std::move(quantumidentity_)),
        broadeningspecies(std::move(broadeningspecies_)),
        lines(nlines, SingleLine(metalocalquanta, metamodel)) {
    if (selfbroadening) broadeningspecies.front() = quantumidentity.Species();
    if (bathbroadening) broadeningspecies.back() = Species::Species::Bath;
  }
  
  void AppendSingleLine(SingleLine&& sl);
  
  void AppendSingleLine(const SingleLine& sl);
  
  [[nodiscard]] bool MatchWithExternal(const SingleLineExternal& sle, const QuantumIdentifier& quantumidentity) const ARTS_NOEXCEPT;
  
  [[nodiscard]] std::pair<bool, bool> Match(const Lines& l) const noexcept;
  
  void sort_by_frequency();
  
  void sort_by_einstein();
  
  [[nodiscard]] String SpeciesName() const noexcept;
  
  [[nodiscard]] String LineShapeMetaData() const noexcept;
  
  [[nodiscard]] Species::Species Species() const noexcept {return quantumidentity.Species();}
  
  [[nodiscard]] Species::IsotopeRecord Isotopologue() const noexcept {return quantumidentity.Isotopologue();}
  
  [[nodiscard]] Index NumLines() const noexcept {return Index(lines.size());}
  
  void MakeLineShapeModelCommon();
  
  [[nodiscard]] Index NumBroadeners() const noexcept {return Index(broadeningspecies.nelem());}
  
  [[nodiscard]] Index NumLocalQuanta() const noexcept {
    return lines.size() ? lines.front().localquanta.val.nelem() : 0;}
  
  [[nodiscard]] Index ZeemanCount(size_t k, Zeeman::Polarization type) const noexcept;
  
  [[nodiscard]] Numeric ZeemanStrength(size_t k, Zeeman::Polarization type, Index i) const noexcept;
  
  [[nodiscard]] Numeric ZeemanSplitting(size_t k, Zeeman::Polarization type, Index i) const noexcept;
  
  void SetAutomaticZeeman() noexcept;
  
  [[nodiscard]] Numeric F_mean(Numeric T=0) const noexcept;
  
  [[nodiscard]] Numeric F_mean(const ConstVectorView& wgts) const noexcept;
  
  [[nodiscard]] bool OnTheFlyLineMixing() const noexcept {
    return population == PopulationType::ByMakarovFullRelmat or
           population == PopulationType::ByRovibLinearDipoleLineMixing;
  }
  
  [[nodiscard]] bool DoLineMixing(Numeric P) const noexcept {
    return linemixinglimit < 0 ? true : linemixinglimit > P;
  }
 
  [[nodiscard]] bool DoVmrDerivative(const QuantumIdentifier& qid) const noexcept {
    return qid.Isotopologue() == quantumidentity.Isotopologue() or 
    (qid.Isotopologue().joker() and qid.Species() == quantumidentity.Species()) or
    std::any_of(broadeningspecies.begin(), broadeningspecies.end(), [s=qid.Species()](auto& a){return a == s;});
  }
 
  [[nodiscard]] bool AnyLinemixing() const noexcept {
    for (auto& line : lines) {
      for (auto& shape : line.lineshape.Data()) {
        if (shape.Y().type not_eq LineShape::TemperatureModel::None or
            shape.G().type not_eq LineShape::TemperatureModel::None or
            shape.DV().type not_eq LineShape::TemperatureModel::None) {
          return true;
        }
      }
    }
    return false;
  }
 
  [[nodiscard]] LineShape::Output ShapeParameters(size_t k, Numeric T, Numeric P, const Vector& vmrs) const noexcept;
  
  [[nodiscard]] LineShape::Output ShapeParameters(size_t k, Numeric T, Numeric P, size_t m) const noexcept;
  
  [[nodiscard]] LineShape::Output ShapeParameters_dT(size_t k, Numeric T, Numeric P, const Vector& vmrs) const noexcept;
  
  [[nodiscard]] Index LineShapePos(const Species::Species spec) const noexcept;
  
  [[nodiscard]] Index LineShapePos(const QuantumIdentifier& qid) const noexcept {
    return LineShapePos(qid.Species());
  }
  
  [[nodiscard]] LineShape::Output ShapeParameters_dVMR(size_t k, Numeric T, Numeric P,
                                         const QuantumIdentifier& vmr_qid) const noexcept;
  
  [[nodiscard]] Numeric ShapeParameter_dInternal(size_t k, Numeric T, Numeric P,
                                   const Vector& vmrs,
                                   const RetrievalQuantity& derivative) const noexcept;
  
  [[nodiscard]] Numeric CutoffFreq(size_t k, Numeric shift=0) const noexcept;
  
  [[nodiscard]] Numeric CutoffFreqMinus(size_t k, Numeric shift=0) const noexcept;
  
  [[nodiscard]] Index BroadeningSpeciesPosition(Species::Species spec) const noexcept {
    if (auto ptr = std::find(broadeningspecies.cbegin(),
      broadeningspecies.cend(), spec); ptr not_eq broadeningspecies.cend())
      return std::distance(broadeningspecies.cbegin(), ptr);
    return -1;
  }
  
  [[nodiscard]] String MetaData() const;
  
  void RemoveLine(Index) noexcept;
  
  SingleLine PopLine(Index) noexcept;
  
  void ReverseLines() noexcept;
  
  [[nodiscard]] Numeric SpeciesMass() const noexcept;
  
  [[nodiscard]] Vector BroadeningSpeciesVMR(const ConstVectorView&, const ArrayOfArrayOfSpeciesTag&) const;
  
  [[nodiscard]] Vector BroadeningSpeciesMass(const ConstVectorView&, const ArrayOfArrayOfSpeciesTag&, const SpeciesIsotopologueRatios&, const Numeric& bath_mass=0) const;
  
  [[nodiscard]] Numeric SelfVMR(const ConstVectorView&, const ArrayOfArrayOfSpeciesTag&) const;
  
  bifstream& read(bifstream& is);
  
  bofstream& write(bofstream& os) const;
  
  [[nodiscard]] bool OK() const noexcept;
  
  [[nodiscard]] Numeric DopplerConstant(Numeric T) const noexcept;
 
  [[nodiscard]] QuantumIdentifier QuantumIdentityOfLine(Index k) const noexcept;
 
  [[nodiscard]] Rational max(QuantumNumberType) const;
 
  friend std::ostream& operator<<(std::ostream&, const Lines&);
 
  friend std::istream& operator>>(std::istream&, Lines&);
};  // 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 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={});
 
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);
 
[[nodiscard]] bool any_cutoff(const Array<Array<Lines>>& abs_lines_per_species);
} // namespace Absorption
 
using AbsorptionSingleLine = Absorption::SingleLine;
using ArrayOfAbsorptionSingleLine = Array<AbsorptionSingleLine>;
using AbsorptionLines = Absorption::Lines;
using ArrayOfAbsorptionLines = Array<AbsorptionLines>;
using ArrayOfArrayOfAbsorptionLines = Array<ArrayOfAbsorptionLines>;
 
using AbsorptionNormalizationType = Absorption::NormalizationType;
using AbsorptionPopulationType = Absorption::PopulationType;
using AbsorptionMirroringType = Absorption::MirroringType;
using AbsorptionCutoffType = Absorption::CutoffType;
 
struct AbsorptionMirroringTagTypeStatus {
  bool None{false}, Lorentz{false}, SameAsLineShape{false}, Manual{false};
  AbsorptionMirroringTagTypeStatus(const ArrayOfArrayOfAbsorptionLines&);
  friend std::ostream& operator<<(std::ostream&, AbsorptionMirroringTagTypeStatus);
};
 
struct AbsorptionNormalizationTagTypeStatus {
  bool None{false}, VVH{false}, VVW{false}, RQ{false}, SFS{false};
  AbsorptionNormalizationTagTypeStatus(const ArrayOfArrayOfAbsorptionLines&);
  friend std::ostream& operator<<(std::ostream&, AbsorptionNormalizationTagTypeStatus);
};
 
struct AbsorptionPopulationTagTypeStatus {
  bool LTE{false}, NLTE{false}, VibTemps{false},
      ByHITRANRosenkranzRelmat{false}, ByHITRANFullRelmat{false},
      ByMakarovFullRelmat{false}, ByRovibLinearDipoleLineMixing{false};
  AbsorptionPopulationTagTypeStatus(const ArrayOfArrayOfAbsorptionLines&);
  friend std::ostream& operator<<(std::ostream&, AbsorptionPopulationTagTypeStatus);
};
 
struct AbsorptionCutoffTagTypeStatus {
  bool None{false}, ByLine{false};
  AbsorptionCutoffTagTypeStatus(const ArrayOfArrayOfAbsorptionLines&);
  friend std::ostream& operator<<(std::ostream&, AbsorptionCutoffTagTypeStatus);
};
 
struct AbsorptionLineShapeTagTypeStatus {
  bool DP{false}, LP{false}, VP{false}, SDVP{false}, HTP{false};
  AbsorptionLineShapeTagTypeStatus(const ArrayOfArrayOfAbsorptionLines&);
  friend std::ostream& operator<<(std::ostream&, AbsorptionLineShapeTagTypeStatus);
};
 
struct AbsorptionTagTypesStatus {
  AbsorptionMirroringTagTypeStatus mirroring;
  AbsorptionNormalizationTagTypeStatus normalization;
  AbsorptionPopulationTagTypeStatus population;
  AbsorptionCutoffTagTypeStatus cutoff;
  AbsorptionLineShapeTagTypeStatus lineshapetype;
 
  AbsorptionTagTypesStatus(const ArrayOfArrayOfAbsorptionLines& lines)
      : mirroring(lines),
        normalization(lines),
        population(lines),
        cutoff(lines),
        lineshapetype(lines) {}
  friend std::ostream& operator<<(std::ostream&, AbsorptionTagTypesStatus);
};
 
struct AbsorptionSpeciesBandIndex {
  Index ispecies;
 
  Index iband;
};
 
AbsorptionSpeciesBandIndex flat_index(
    Index i,
    const ArrayOfArrayOfSpeciesTag& abs_species,
    const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species);
 
#endif  // absorptionlines_h