lineshapemodel.h

Go to the documentation of this file.

/* Copyright (C) 2018
 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 lineshapemodel_h
#define lineshapemodel_h
 
#include "constants.h"
#include "enums.h"
#include "file.h"
#include "jacobian.h"
#include "species_tags.h"
#include <algorithm>
#include <numeric>
#include <utility>
 
Jacobian::Line select_derivativeLineShape(const String& var,
                                          const String& coeff);
 
namespace LineShape {
 
ENUMCLASS(TemperatureModel, char,
  None,   // 0
  T0,     // Constant, X0
  T1,     // Standard, X0 * (T0/T) ^ X1
  T2,     // X0 * (T0/T) ^ X1 * (1 + X2 * log(T/T0));
  T3,     // X0 + X1 * (T - T0)
  T4,     // (X0 + X1 * (T0/T - 1)) * (T0/T)^X2;
  T5,     // X0 * (T0/T)^(0.25 + 1.5*X1)
  LM_AER, // X(200) = X0; X(250) = X1; X(298) = X2; X(340) = X3;  Linear interpolation in between
  DPL,    // X0 * (T0/T) ^ X1 + X2 * (T0/T) ^ X3
  POLY    // X0 + X1 * T + X2 * T ^ 2 + X3 * T ^ 3
)
 
 
ENUMCLASS(Variable, char,
  G0,   // Pressure broadening speed-independent
  D0,   // Pressure f-shifting speed-dependent
  G2,   // Pressure broadening speed-dependent
  D2,   // Pressure f-shifting speed-independent
  FVC,  // Frequency of velocity-changing collisions
  ETA,  // Correlation
  Y,    // First order line mixing coefficient
  G,    // Second order line mixing coefficient
  DV    // Second order line mixing f-shifting
)
 
struct ModelParameters {
  static constexpr Index N = 4;
  static_assert(Index(Options::LineShapeCoeff::FINAL) == N, "Must update either LineShapeCoeff options or ModelParameters");
  TemperatureModel type;
  Numeric X0;
  Numeric X1;
  Numeric X2;
  Numeric X3;
  
  constexpr ModelParameters(TemperatureModel intype=TemperatureModel::None,
                            Numeric inX0=std::numeric_limits<Numeric>::quiet_NaN(),
                            Numeric inX1=std::numeric_limits<Numeric>::quiet_NaN(),
                            Numeric inX2=std::numeric_limits<Numeric>::quiet_NaN(),
                            Numeric inX3=std::numeric_limits<Numeric>::quiet_NaN())
  noexcept : type(intype), X0(inX0), X1(inX1), X2(inX2), X3(inX3) {}
  
  template <typename VectorType> constexpr
  ModelParameters(TemperatureModel intype, VectorType&& v) ARTS_NOEXCEPT :
  ModelParameters(intype) {
    const auto n = std::size(v);
    ARTS_ASSERT(n <= N, "Must have at most ", N, " inputs, got: ", n)
    switch (n) {
      case 4: X3 = v[3]; [[fallthrough]];
      case 3: X2 = v[2]; [[fallthrough]];
      case 2: X1 = v[1]; [[fallthrough]];
      case 1: X0 = v[0];
    }
  }
  
  [[nodiscard]] constexpr Numeric special_linemixing_aer(Numeric T) const noexcept {
    if (T < 250.0)
      return X0 + (T - 200.0) * (X1 - X0) / (250.0 - 200.0);
    if (T > 296.0)
      return X2 + (T - 296.0) * (X3 - X2) / (340.0 - 296.0);
    return X1 + (T - 250.0) * (X2 - X1) / (296.0 - 250.0);
  }
  
  [[nodiscard]] constexpr Numeric special_linemixing_aer_dT(Numeric T) const noexcept {
    if (T < 250.0)
      return (X1 - X0) / (250.0 - 200.0);
    if (T > 296.0)
      return (X3 - X2) / (340.0 - 296.0);
    return (X2 - X1) / (296.0 - 250.0);
  }
  
  static constexpr Numeric special_linemixing_aer_dX0(Numeric T) noexcept {
    if (T < 250.0)
      return 1 - (T - 200.0) / (250.0 - 200.0);
    return 0;
  }
  
  static constexpr Numeric special_linemixing_aer_dX1(Numeric T) noexcept {
    if (T < 250.0)
      return (T - 200.0) / (250.0 - 200.0);
    if (T > 296.0)
      return 0;
    return 1 - (T - 250.0) / (296.0 - 250.0);
  }
  
  static constexpr Numeric special_linemixing_aer_dX2(Numeric T) noexcept {
    if (T < 250.0)
      return 0;
    if (T > 296.0)
      return 1 - (T - 296.0)  / (340.0 - 296.0);
    return (T - 250.0) / (296.0 - 250.0);
  }
  
  static constexpr Numeric special_linemixing_aer_dX3(Numeric T) noexcept {
    if (T > 296.0)
      return (T - 296.0) / (340.0 - 296.0);
    return 0;
  }
  
  [[nodiscard]] Numeric at(Numeric T, Numeric T0) const noexcept;
  
  [[nodiscard]] [[nodiscard]] Numeric dX0(Numeric T, Numeric T0) const noexcept;
  
  [[nodiscard]] Numeric dX1(Numeric T, Numeric T0) const noexcept;
  
  [[nodiscard]] Numeric dX2(Numeric T, Numeric T0) const noexcept;
  
  [[nodiscard]] Numeric dX3(Numeric T, Numeric T0) const noexcept;
  
  [[nodiscard]] Numeric dT(Numeric T, Numeric T0) const noexcept;
  
  [[nodiscard]] Numeric dT0(Numeric T, Numeric T0) const noexcept;
 
  friend std::ostream& operator<<(std::ostream& os, const ModelParameters& mp);
 
  friend std::istream& operator>>(std::istream& is, ModelParameters& mp);
};
 
String modelparameters2metadata(const ModelParameters mp, const Numeric T0);
 
Numeric& SingleModelParameter(ModelParameters& mp, const String& type);
 
constexpr ModelParameters modelparameterGetEmpty(const TemperatureModel t) noexcept {
  switch(t) {
    case TemperatureModel::None:   // 0
      return {TemperatureModel::None, 0, 0, 0, 0};
    case TemperatureModel::T0:     // Constant, X0
      return {TemperatureModel::T0, 0, 0, 0, 0};
    case TemperatureModel::T1:     // Standard, X0 * (T0/T) ^ X1
      return {TemperatureModel::T1, 0, 0, 0, 0};
    case TemperatureModel::T2:     // X0 * (T0/T) ^ X1 * (1 + X2 * log(T/T0));
      return {TemperatureModel::T2, 0, 0, 0, 0};
    case TemperatureModel::T3:     // X0 + X1 * (T - T0)
      return {TemperatureModel::T3, 0, 0, 0, 0};
    case TemperatureModel::T4:     // (X0 + X1 * (T0/T - 1)) * (T0/T)^X2;
      return {TemperatureModel::T4, 0, 0, 0, 0};
    case TemperatureModel::T5:     // X0 * (T0/T)^(0.25 + 1.5*X1)
      return {TemperatureModel::T5, 0, 0, 0, 0};
    case TemperatureModel::LM_AER: // X(200) = X0; X(250) = X1; X(298) = X2; X(340) = X3;  Linear interpolation in between
      return {TemperatureModel::LM_AER, 0, 0, 0, 0};
    case TemperatureModel::DPL:    // X0 * (T0/T) ^ X1 + X2 * (T0/T) ^ X3
      return {TemperatureModel::DPL, 0, 0, 0, 0};
    case TemperatureModel::POLY:
      return {TemperatureModel::POLY, 0, 0, 0, 0};
    case TemperatureModel::FINAL:
      return {TemperatureModel::None, 0, 0, 0, 0};
  }
  return {TemperatureModel::None, 0, 0, 0, 0};
}
 
constexpr bool modelparameterEmpty(const ModelParameters mp) noexcept {
  switch(mp.type) {
    case TemperatureModel::None:   // 0
      return true;
    case TemperatureModel::T0:     // Constant, X0
      return (mp.X0 == 0);
    case TemperatureModel::T1:     // Standard, X0 * (T0/T) ^ X1
      return (mp.X0 == 0);
    case TemperatureModel::T2:     // X0 * (T0/T) ^ X1 * (1 + X2 * log(T/T0));
      return (mp.X0 == 0);
    case TemperatureModel::T3:     // X0 + X1 * (T - T0)
      return (mp.X0 == 0 and mp.X1 == 0);
    case TemperatureModel::T4:     // (X0 + X1 * (T0/T - 1)) * (T0/T)^X2;
      return (mp.X0 == 0 and mp.X1 == 0);
    case TemperatureModel::T5:     // X0 * (T0/T)^(0.25 + 1.5*X1)
      return (mp.X0 == 0);
    case TemperatureModel::LM_AER: // X(200) = X0; X(250) = X1; X(298) = X2; X(340) = X3;  Linear interpolation in between
      return (mp.X0 == 0 and mp.X1 == 0 and mp.X2 == 0 and mp.X3 == 0);
    case TemperatureModel::DPL:    // X0 * (T0/T) ^ X1 + X2 * (T0/T) ^ X3
      return (mp.X0 == 0 and mp.X2 == 0);
    case TemperatureModel::POLY:
      return (mp.X0 == 0 and mp.X1 == 0 and mp.X2 == 0 and mp.X3 == 0);
    case TemperatureModel::FINAL:
      return true;
  }
  return true;
}
 
constexpr Numeric modelparameterFirstExponent(const ModelParameters mp) noexcept {
  switch(mp.type) {
    case TemperatureModel::None:   // 0
      return 0;
    case TemperatureModel::T0:     // Constant, X0
      return 0;
    case TemperatureModel::T1:     // Standard, X0 * (T0/T) ^ X1
      return mp.X1;
    case TemperatureModel::T2:     // X0 * (T0/T) ^ X1 * (1 + X2 * log(T/T0));
      return mp.X1;
    case TemperatureModel::T3:     // X0 + X1 * (T - T0)
      return 0;
    case TemperatureModel::T4:     // (X0 + X1 * (T0/T - 1)) * (T0/T)^X2;
      return mp.X2;
    case TemperatureModel::T5:     // X0 * (T0/T)^(0.25 + 1.5*X1)
      return (0.25 + 1.5*mp.X1);
    case TemperatureModel::LM_AER: // X(200) = X0; X(250) = X1; X(298) = X2; X(340) = X3;  Linear interpolation in between
      return 0;
    case TemperatureModel::DPL:    // X0 * (T0/T) ^ X1 + X2 * (T0/T) ^ X3
      return mp.X1;
    case TemperatureModel::POLY:
      return 0;
    case TemperatureModel::FINAL:
      return std::numeric_limits<Numeric>::quiet_NaN();
  }
  return std::numeric_limits<Numeric>::quiet_NaN();
}
 
constexpr Index nVars = Index(Variable::FINAL);
 
class SingleSpeciesModel {
 private:
  std::array<ModelParameters, nVars> X;
 
 public:
  constexpr SingleSpeciesModel(
    ModelParameters G0=ModelParameters{},
    ModelParameters D0=ModelParameters{},
    ModelParameters G2=ModelParameters{},
    ModelParameters D2=ModelParameters{},
    ModelParameters FVC=ModelParameters{},
    ModelParameters ETA=ModelParameters{},
    ModelParameters Y=ModelParameters{},
    ModelParameters G=ModelParameters{},
    ModelParameters DV=ModelParameters{})
      : X({G0, D0, G2, D2, FVC, ETA, Y, G, DV}) {}
 
#define ACCESS_INTERNAL(VARPOS)                                             \
  constexpr ModelParameters& VARPOS() noexcept { return std::get<Index(Variable::VARPOS)>(X); } \
  constexpr ModelParameters VARPOS() const noexcept { return std::get<Index(Variable::VARPOS)>(X); }
  ACCESS_INTERNAL(G0);
  ACCESS_INTERNAL(D0);
  ACCESS_INTERNAL(G2);
  ACCESS_INTERNAL(D2);
  ACCESS_INTERNAL(FVC);
  ACCESS_INTERNAL(ETA);
  ACCESS_INTERNAL(Y);
  ACCESS_INTERNAL(G);
  ACCESS_INTERNAL(DV);
#undef ACCESS_INTERNAL
 
  constexpr std::array<ModelParameters, nVars>& Data() noexcept { return X; }
  
  [[nodiscard]] constexpr const std::array<ModelParameters, nVars>& Data() const noexcept { return X; }
 
  constexpr void Set(Variable var, const ModelParameters& x) noexcept {
#define MODELPARAMCASESETTER(X) \
  case Variable::X:             \
    X() = x;                    \
    break
    switch (var) {
      MODELPARAMCASESETTER(G0);
      MODELPARAMCASESETTER(D0);
      MODELPARAMCASESETTER(G2);
      MODELPARAMCASESETTER(D2);
      MODELPARAMCASESETTER(FVC);
      MODELPARAMCASESETTER(ETA);
      MODELPARAMCASESETTER(Y);
      MODELPARAMCASESETTER(G);
      MODELPARAMCASESETTER(DV);
      case Variable::FINAL: break;
    }
#undef MODELPARAMCASESETTER
  }
 
  [[nodiscard]] constexpr ModelParameters Get(Variable var) const noexcept {
  #define MODELPARAMCASEGETTER(X) case Variable::X: out = X(); break;  
  ModelParameters out{};
  switch (var) {
    MODELPARAMCASEGETTER(G0);
    MODELPARAMCASEGETTER(D0);
    MODELPARAMCASEGETTER(G2);
    MODELPARAMCASEGETTER(D2);
    MODELPARAMCASEGETTER(FVC);
    MODELPARAMCASEGETTER(ETA);
    MODELPARAMCASEGETTER(Y);
    MODELPARAMCASEGETTER(G);
    MODELPARAMCASEGETTER(DV);
    case Variable::FINAL: break;
  }
  return out;
  #undef MODELPARAMCASEGETTER
  }
 
  bifstream& read(bifstream& bif);
 
  bofstream& write(bofstream& bof) const;
  
  [[nodiscard]] std::pair<bool, bool> MatchTypes(const SingleSpeciesModel& other) const noexcept;
 
  friend std::ostream& operator<<(std::ostream& os, const SingleSpeciesModel& ssm);
 
  friend std::istream& operator>>(std::istream& is, SingleSpeciesModel& ssm);
};
 
ENUMCLASS(Type, char,
  DP,    // Doppler
  LP,    // Lorentz
  VP,    // Voigt
  SDVP,  // Speed-dependent Voigt
  HTP    // Hartmann-Tran
)
 
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wreturn-type"
constexpr std::string_view shapetype2metadatastring(Type type) noexcept {
  switch (type) {
    case Type::DP:
      return "The line shape type is the Doppler profile\n";
    case Type::LP:
      return "The line shape type is the Lorentz profile.\n";
    case Type::VP:
      return "The line shape type is the Voigt profile.\n";
    case Type::SDVP:
      return "The line shape type is the speed-dependent Voigt profile.\n";
    case Type::HTP:
      return "The line shape type is the Hartmann-Tran profile.\n";
    case Type::FINAL:
      return "There's an error.\n";
  }
}
#pragma GCC diagnostic pop
 
struct Output {
  Numeric G0{0},  // Pressure broadening speed-independent
      D0{0},      // Pressure f-shifting speed-independent
      G2{0},      // Pressure broadening speed-dependent
      D2{0},      // Pressure f-shifting speed-dependent
      FVC{0},     // Frequency of velocity-changing collisions
      ETA{0},     // Correlation
      Y{0},       // First order line mixing coefficient
      G{0},       // Second order line mixing coefficient
      DV{0};      // Second order line mixing f-shifting
  constexpr Output() noexcept = default;
  constexpr Output(Numeric g0, Numeric d0, Numeric g2,
                   Numeric d2, Numeric fvc, Numeric eta,
                   Numeric y, Numeric g, Numeric dv) noexcept :
    G0(g0), D0(d0), G2(g2), D2(d2), FVC(fvc), ETA(eta), Y(y), G(g), DV(dv) {}
  constexpr Output(const Output&) noexcept = default;
  constexpr Output(Output&&) noexcept = default;
  constexpr Output& operator=(const Output&) noexcept = default;
  constexpr Output& operator=(Output&&) noexcept = default;
 
  friend std::ostream& operator<<(std::ostream& os, Output x);
};
 
constexpr Output mirroredOutput(Output x) noexcept {
  return {x.G0, -x.D0, x.G2, -x.D2, x.FVC, x.ETA, x.Y, x.G, -x.DV};
}
 
constexpr Output negativeOutput(Output x) noexcept {
  return {-x.G0, -x.D0, -x.G2, -x.D2, -x.FVC, -x.ETA, -x.Y, -x.G, -x.DV};
}
 
constexpr Output si2cgs(Output x) noexcept {
  using Conversion::freq2kaycm;
  return {freq2kaycm(x.G0),
          freq2kaycm(x.D0),
          freq2kaycm(x.D2),
          freq2kaycm(x.G2),
          freq2kaycm(x.FVC),
          x.ETA,
          x.Y,
          x.G,
          freq2kaycm(x.DV)};
}
 
constexpr Output differenceOutput(Output y, Output x) noexcept {
  return {y.G0 - x.G0,
          y.D0 - x.D0,
          y.G2 - x.G2,
          y.D2 - x.D2,
          y.FVC - x.FVC,
          y.ETA - x.ETA,
          y.Y - x.Y,
          y.G - x.G,
          y.DV - x.DV};
}
 
Vector vmrs(const ConstVectorView& atmospheric_vmrs,
            const ArrayOfArrayOfSpeciesTag& atmospheric_species,
            const ArrayOfSpecies& lineshape_species) ARTS_NOEXCEPT;
 
Vector mass(const ConstVectorView& atmospheric_vmrs,
            const ArrayOfArrayOfSpeciesTag& atmospheric_species,
            const ArrayOfSpecies& lineshape_species,
            const SpeciesIsotopologueRatios& ir) ARTS_NOEXCEPT;
 
static constexpr std::string_view bath_broadening = "AIR";
 
static constexpr std::string_view self_broadening = "SELF";
 
class Model {
 private:
   std::vector<SingleSpeciesModel> mdata;
 
 public:
  Model(Index n=0) noexcept : mdata(n) {}
  
  explicit Model(const std::vector<SingleSpeciesModel>& assm) noexcept : mdata(assm) {}
  
  Model(const Model& m) noexcept : Model(m.mdata) {}
  
  explicit Model(std::vector<SingleSpeciesModel>&& assm) noexcept : mdata(std::move(assm)) {}
  
  Model(Model&& m) noexcept : Model(std::move(m.mdata)) {}
  
  Model& operator=(const Model& m) = default;
  
  Model& operator=(Model&& m) = default;
  
  Model(Numeric sgam,
        Numeric nself,
        Numeric agam,
        Numeric nair,
        Numeric psf,
        std::array<Numeric, 12> aer_interp = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}) noexcept;
  
  [[nodiscard]] bool OK(Type type, bool self, bool bath,
          const std::size_t nspecies) const noexcept;
 
#define LSPC(XVAR, PVAR)                                                       \
  Numeric XVAR(                                                                \
      Numeric T, Numeric T0, Numeric P [[maybe_unused]], const ConstVectorView& vmrs) \
      const noexcept;
  LSPC(G0, P)
  LSPC(D0, P)
  LSPC(G2, P)
  LSPC(D2, P)
  LSPC(FVC, P)
  LSPC(ETA, 1)
  LSPC(Y, P)
  LSPC(G, P* P)
  LSPC(DV, P* P)
#undef LSPC
 
#define LSPCV(XVAR, PVAR)                                            \
  Numeric d##XVAR##_dVMR(Numeric T,                                  \
                         Numeric T0,                                 \
                         Numeric P [[maybe_unused]],                 \
                         const Index deriv_pos) const noexcept;
  LSPCV(G0, P)
  LSPCV(D0, P)
  LSPCV(G2, P)
  LSPCV(D2, P)
  LSPCV(FVC, P)
  LSPCV(ETA, 1)
  LSPCV(Y, P)
  LSPCV(G, P* P)
  LSPCV(DV, P* P)
#undef LSPCV
 
#define LSPCT(XVAR, PVAR)                                                       \
  Numeric d##XVAR##_dT(                                                         \
      Numeric T, Numeric T0, Numeric P [[maybe_unused]], const ConstVectorView& vmrs)  \
      const noexcept;
  LSPCT(G0, P)
  LSPCT(D0, P)
  LSPCT(G2, P)
  LSPCT(D2, P)
  LSPCT(FVC, P)
  LSPCT(ETA, 1)
  LSPCT(Y, P)
  LSPCT(G, P* P)
  LSPCT(DV, P* P)
#undef LSPCT
 
// All shape model derivatives
#define LSPDC(XVAR, DERIV, PVAR)                                     \
  Numeric d##XVAR##DERIV(Numeric T,                                  \
                         Numeric T0,                                 \
                         Numeric P [[maybe_unused]],                 \
                         Index deriv_pos,                            \
                         const ConstVectorView& vmrs) const noexcept;
  LSPDC(G0, _dT0, P)
  LSPDC(G0, _dX0, P)
  LSPDC(G0, _dX1, P)
  LSPDC(G0, _dX2, P)
  LSPDC(G0, _dX3, P)
  LSPDC(D0, _dT0, P)
  LSPDC(D0, _dX0, P)
  LSPDC(D0, _dX1, P)
  LSPDC(D0, _dX2, P)
  LSPDC(D0, _dX3, P)
  LSPDC(G2, _dT0, P)
  LSPDC(G2, _dX0, P)
  LSPDC(G2, _dX1, P)
  LSPDC(G2, _dX2, P)
  LSPDC(G2, _dX3, P)
  LSPDC(D2, _dT0, P)
  LSPDC(D2, _dX0, P)
  LSPDC(D2, _dX1, P)
  LSPDC(D2, _dX2, P)
  LSPDC(D2, _dX3, P)
  LSPDC(FVC, _dT0, P)
  LSPDC(FVC, _dX0, P)
  LSPDC(FVC, _dX1, P)
  LSPDC(FVC, _dX2, P)
  LSPDC(FVC, _dX3, P)
  LSPDC(ETA, _dT0, 1)
  LSPDC(ETA, _dX0, 1)
  LSPDC(ETA, _dX1, 1)
  LSPDC(ETA, _dX2, 1)
  LSPDC(ETA, _dX3, 1)
  LSPDC(Y, _dT0, P)
  LSPDC(Y, _dX0, P)
  LSPDC(Y, _dX1, P)
  LSPDC(Y, _dX2, P)
  LSPDC(Y, _dX3, P)
  LSPDC(G, _dT0, P* P)
  LSPDC(G, _dX0, P* P)
  LSPDC(G, _dX1, P* P)
  LSPDC(G, _dX2, P* P)
  LSPDC(G, _dX3, P* P)
  LSPDC(DV, _dT0, P* P)
  LSPDC(DV, _dX0, P* P)
  LSPDC(DV, _dX1, P* P)
  LSPDC(DV, _dX2, P* P)
  LSPDC(DV, _dX3, P* P)
#undef LSPDC
 
  [[nodiscard]] Output GetParams(Numeric T,
                   Numeric T0,
                   Numeric P,
                   const ConstVectorView& vmrs) const noexcept;
 
  [[nodiscard]] Output GetParams(Numeric T,
                   Numeric T0,
                   Numeric P,
                   size_t k) const noexcept;
 
  [[nodiscard]] Output GetTemperatureDerivs(Numeric T,
                              Numeric T0,
                              Numeric P,
                              ConstVectorView vmrs) const noexcept;
 
  [[nodiscard]] Output GetVMRDerivs(Numeric T, Numeric T0, Numeric P, const Index pos) const noexcept;
  
  [[nodiscard]] Numeric GetInternalDeriv(Numeric T,
                           Numeric T0,
                           Numeric P,
                           Index pos,
                           const Vector& vmrs,
                           Jacobian::Line deriv) const noexcept;
 
  [[nodiscard]] Index nelem() const { return Index(mdata.size()); }
  
  [[nodiscard]] Index size() const { return Index(mdata.size()); }
  
  void resize(Index n) {mdata.resize(n);}
  
  void reserve(Index n) {mdata.reserve(n);}
  
  SingleSpeciesModel& operator[](Index i) {return mdata[i];}
  
  const SingleSpeciesModel& operator[](Index i) const {return mdata[i];}
  
  auto begin() { return mdata.begin(); }
  auto end() { return mdata.end(); }
 
  [[nodiscard]] const std::vector<SingleSpeciesModel>& Data() const noexcept { return mdata; }
  
  std::vector<SingleSpeciesModel>& Data() noexcept { return mdata; }
 
  void Remove(Index i, ArrayOfSpeciesTag& specs);
  void Remove(Index i, ArrayOfSpecies& specs);
 
  void SetLineMixingModel(SingleSpeciesModel x);
  
  [[nodiscard]] std::pair<bool, bool> Match(const Model& other) const noexcept;
  
  friend
  std::istream& from_linefunctiondata(std::istream& data,
                                      Type& type,
                                      bool& self,
                                      bool& bath,
                                      Model& m,
                                      ArrayOfSpecies& species);
  
  friend
  std::istream& from_artscat4(std::istream& is,
                              Type& type,
                              bool& self,
                              bool& bath,
                              Model& m,
                              ArrayOfSpecies& species,
                              const QuantumIdentifier& qid);
 
  friend std::ostream& operator<<(std::ostream&, const Model&);
 
  friend std::istream& operator>>(std::istream&, Model&);
  
  bifstream& read(bifstream& bif);
  
  bofstream& write(bofstream& bof) const;
};  // Model;
 
 
Model hitran_model(Numeric sgam,
                   Numeric nself,
                   Numeric agam,
                   Numeric nair,
                   Numeric psf);
 
Model lblrtm_model(Numeric sgam,
                   Numeric nself,
                   Numeric agam,
                   Numeric nair,
                   Numeric psf,
                   std::array<Numeric, 12> aer_interp = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0});
 
String ModelShape2MetaData(const Model& m);
 
Model MetaData2ModelShape(const String& s);
 
ArrayOfString ModelMetaDataArray(const Model& m, const bool self, const ArrayOfSpecies& sts, const Numeric T0);
 
std::istream& from_artscat4(std::istream& is,
                            Type& type,
                            bool& self,
                            bool& bath,
                            Model& m,
                            ArrayOfSpecies& species,
                            const QuantumIdentifier& qid);
 
std::istream& from_linefunctiondata(std::istream& data,
                                    Type& type,
                                    bool& self,
                                    bool& bath,
                                    Model& m,
                                    ArrayOfSpecies& species);
 
std::istream& from_linemixingdata(std::istream& data, Model& lsc);
 
std::istream& from_pressurebroadeningdata(std::istream& data,
                                          LineShape::Type& type,
                                          bool& self,
                                          bool& bath,
                                          Model& m,
                                          ArrayOfSpecies& species,
                                          const QuantumIdentifier& qid);
 
namespace LegacyLineFunctionData {
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wreturn-type"
constexpr Index temperaturemodel2legacynelem(TemperatureModel type) noexcept {
  switch (type) {
    case TemperatureModel::None:
      return 0;
    case TemperatureModel::T0:
      return 1;
    case TemperatureModel::T1:
      return 2;
    case TemperatureModel::T2:
      return 3;
    case TemperatureModel::T3:
      return 2;
    case TemperatureModel::T4:
      return 3;
    case TemperatureModel::T5:
      return 2;
    case TemperatureModel::LM_AER:
      return 12;
    case TemperatureModel::DPL:
      return 4;
    case TemperatureModel::POLY:
      return 4;
    case TemperatureModel::FINAL: break;
  }
}
#pragma GCC diagnostic pop
 
std::vector<Variable> lineshapetag2variablesvector(String type);
 
std::vector<Variable> linemixingtag2variablesvector(String type);
};  // namespace LegacyLineFunctionData
 
namespace LegacyLineMixingData {
enum class TypeLM {
  LM_NONE,                  // Reserved for no line mixing
  LM_LBLRTM,                // Reserved for LBLRTM line mixing
  LM_LBLRTM_O2NonResonant,  // Reserved for the non-resonant O2 line in LBLRTM
  LM_1STORDER,  // Reserved for Tretyakov et al. 2005 1st order of line mixing
  LM_2NDORDER,  // Reserved for Makarov et al. 2011 second order of line mixing
  LM_BYBAND  // Reserved for Paris data of relaxation matrix line mixing for band
};
 
LegacyLineMixingData::TypeLM string2typelm(String type);
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wreturn-type"
constexpr Index typelm2nelem(LegacyLineMixingData::TypeLM type) {
  switch (type) {
    case TypeLM::LM_NONE:  // The standard case
      return 0;
    case TypeLM::LM_LBLRTM:  // The LBLRTM case
      return 12;
    case TypeLM::LM_LBLRTM_O2NonResonant:  // Nonresonant is just a tag
      return 1;
    case TypeLM::LM_2NDORDER:  // The 2nd order case
      return 10;
    case TypeLM::LM_1STORDER:  // The 2nd order case
      return 3;
    case TypeLM::LM_BYBAND:  // The band class
      return 1;
  }
}
#pragma GCC diagnostic pop
 
Model vector2modellm(Vector x, LegacyLineMixingData::TypeLM type);
};  // namespace LegacyLineMixingData
 
namespace LegacyPressureBroadeningData {
enum class TypePB {
  PB_NONE,                      // No pressure broadening
  PB_AIR_BROADENING,            // Air broadening and self broadening only
  PB_AIR_AND_WATER_BROADENING,  // Air, water, and self broadening
  PB_PLANETARY_BROADENING,  // Gas broadening as done for solar system planets
  // PB_SD_AIR_VOLUME,                 // HTP in air for SD limit NOT SUPPORTED
  // PB_HTP_AIR_VOLUME,                // HTP in air NOT SUPPORTED
  // PB_VOIGT_TEST_WATER,              // Voigt parameters for testing NOT SUPPORTED
  // PB_SD_TEST_WATER,                 // SD parameters for testing NOT SUPPORTED
  // PB_PURELY_FOR_TESTING             // Testing tag for new input structures --- can be changed by anyone... NOT SUPPORTED
};
 
LegacyPressureBroadeningData::TypePB string2typepb(String type);
 
Index self_listed(const QuantumIdentifier& qid,
                  LegacyPressureBroadeningData::TypePB t);
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wreturn-type"
constexpr Index typepb2nelem(LegacyPressureBroadeningData::TypePB type)  {
  switch (type) {
    case TypePB::PB_NONE:
      return 0;
    case TypePB::PB_AIR_BROADENING:
      return 10;
    case TypePB::PB_AIR_AND_WATER_BROADENING:
      return 9;
    case TypePB::PB_PLANETARY_BROADENING:
      return 20;
  }
}
#pragma GCC diagnostic pop
 
void vector2modelpb(LineShape::Type& mtype,
                    bool& self,
                    bool& bath,
                    Model& m,
                    ArrayOfSpecies& species,
                    Vector x,
                    LegacyPressureBroadeningData::TypePB type,
                    bool self_in_list,
                    Species::Species self_spec);
};  // namespace LegacyPressureBroadeningData
};  // namespace LineShape
 
using LineShapeModelParameters = LineShape::ModelParameters;
using LineShapeModel = LineShape::Model;
using LineShapeSingleSpeciesModel = LineShape::SingleSpeciesModel;
 
using LineShapeType = LineShape::Type;
using LineShapeVariable = LineShape::Variable;
using LineShapeTemperatureModel = LineShape::TemperatureModel;
 
#endif  // lineshapemodel_h