linerecord.h

Go to the documentation of this file.

/* Copyright (C) 2000-2013
   Stefan Buehler <sbuehler@ltu.se>
   Axel von Engeln <engeln@uni-bremen.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 linerecord_h
#define linerecord_h

#include <cmath>
#include <stdexcept>
#include "array.h"
#include "constants.h"
#include "lineshapemodel.h"
#include "matpackI.h"
#include "messages.h"
#include "mystring.h"
#include "quantum.h"
#include "wigner_functions.h"
#include "zeemandata.h"

/* Forward declaration of classes */
class SpeciesRecord;
class IsotopologueRecord;

enum class MirroringType : Index {
  // Describes the type of mirroring line effects.
  // LineRecord effect:  Adds MTM (String)MirroringType under line modifications when saving ARTSCAT5
  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
};

enum class LineNormalizationType : Index {
  // Describes the type of normalization line effects.
  // LineRecord effect:  Adds LNT (String)LineNormalizationType under line modifications when saving ARTSCAT5
  // 0 is ignored on printing
  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
};

enum class LinePopulationType : Index {
  // Describes the type of population level counter.
  // LineRecord effect:  Adds POP (String)LinePopulationType under line modifications when saving ARTSCAT5
  ByLTE,                      // Assume line is in LTE
  ByVibrationalTemperatures,  // Assume line is in NLTE described by vibrational temperatures
  ByPopulationDistribution,  // Assume line is in NLTE and the upper-to-lower ratio is known
};

class LineRecord {
 public:
  LineRecord()
      : mversion(5),
        mqid(QuantumIdentifier::TRANSITION, 1000000, 1000000),
        mf(0.),
        mi0(0.),
        mti0(0.),
        melow(0.),
        mevlow(-1.0),
        mnlte_lower_index(-1),
        mevupp(-1.0),
        mnlte_upper_index(-1),
        ma(NAN),
        mgupper(NAN),
        mglower(NAN),
        mstandard(true),
        mcutoff(-1.0),
        mlinemixing_limit(-1.0),
        mmirroring(MirroringType::None),
        mlinenorm(LineNormalizationType::None),
        mpopulation(LinePopulationType::ByLTE) { /* Nothing to do here. */
  }

  LineRecord(Index species,
             Index isotopologue,
             Numeric f,
             Numeric psf,
             Numeric i0,
             Numeric ti0,
             Numeric elow,
             Numeric agam,
             Numeric sgam,
             Numeric nair,
             Numeric nself,
             Numeric /*tgam*/,
             const ArrayOfNumeric& /*aux*/,
             Numeric /* df */,
             Numeric /* di0 */,
             Numeric /* dagam */,
             Numeric /* dsgam */,
             Numeric /* dnair */,
             Numeric /* dnself */,
             Numeric /* dpsf */)
      : mversion(5),
        mqid(QuantumIdentifier::TRANSITION, species, isotopologue),
        mf(f),
        mi0(i0),
        mti0(ti0),
        melow(elow),
        mevlow(-1.0),
        mnlte_lower_index(-1),
        mevupp(-1.0),
        mnlte_upper_index(-1),
        ma(NAN),
        mgupper(NAN),
        mglower(NAN),
        mlineshapemodel(sgam, nself, agam, nair, psf),
        mstandard(true),
        mcutoff(-1.0),
        mlinemixing_limit(-1.0),
        mmirroring(MirroringType::None),
        mlinenorm(LineNormalizationType::None),
        mpopulation(LinePopulationType::ByLTE) {}

  String VersionString() const;

  Index Version() const { return mversion; }
  void SetVersion5() { mversion = 5; }

  /* Set to latest compatible version of the catalog to store all data */
  void SetVersionToLatest() { mversion = 5; }

  Index Species() const { return mqid.Species(); }

  Index Isotopologue() const { return mqid.Isotopologue(); }

  String Name() const;

  const SpeciesRecord& SpeciesData() const;

  const IsotopologueRecord& IsotopologueData() const;

  Numeric F() const { return mf; }

  void setF(Numeric new_mf) { mf = new_mf; }

  Numeric Psf() const { return mlineshapemodel.Data().back().D0().X0; }

  Numeric I0() const { return mi0; }

  void setI0(Numeric new_mi0) { mi0 = new_mi0; }

  Numeric Ti0() const { return mti0; }
  //FIXME: really in cm-1 still?
  Numeric Elow() const { return melow; }
  void SetElow(Numeric new_melow) { melow = new_melow; }
  //FIXME: really in cm-1 still?
  Numeric Evlow() const { return mevlow; }
  void SetEvlow(Numeric evlow) { mevlow = evlow; }
  Index NLTELowerIndex() const { return mnlte_lower_index; }
  void SetNLTELowerIndex(Index nlte_lower_index) {
    mnlte_lower_index = nlte_lower_index;
  }
  //FIXME: really in cm-1 still?
  Numeric Evupp() const { return mevupp; }
  void SetEvupp(Numeric evupp) { mevupp = evupp; }
  Index NLTEUpperIndex() const { return mnlte_upper_index; }
  void SetNLTEUpperIndex(Index nlte_upper_index) {
    mnlte_upper_index = nlte_upper_index;
  }

  Numeric Agam() const { return mlineshapemodel.Data().back().G0().X0; }

  Numeric Sgam() const { return mlineshapemodel.Data().front().G0().X0; }

  Numeric Nair() const { return mlineshapemodel.Data().back().G0().X1; }

  Numeric Nself() const { return mlineshapemodel.Data().front().G0().X1; }

  Numeric A() const { return ma; }

  Numeric electric_dipole_moment_squared() const {
    using namespace Constant;
    return 3 * h * epsilon_0 / (16 * pi) * ma * pow3(c / mf);
  }

  Numeric magnetic_quadrapole_moment_squared() const {
    using namespace Constant;
    return 3 * h / (16 * pi * mu_0) * ma * pow3(c / mf);
  }

  Numeric reduced_rovibrational_dipole(Rational k = 1) const {
    const Rational Jf = LowerQuantumNumber(QuantumNumberType::J);
    const Rational Ji = UpperQuantumNumber(QuantumNumberType::J);
    const Rational lf = LowerQuantumNumber(QuantumNumberType::l2);
    const Rational li = UpperQuantumNumber(QuantumNumberType::l2);
    const Numeric val =
        sqrt(2 * Jf + 1) * wigner3j(Jf, k, Ji, li, lf - li, -lf);
    if ((Jf + lf + 1) % 2)
      return -val;
    else
      return +val;
  }

  Numeric reduced_magnetic_quadrapole() const {
    const Rational Jf = LowerQuantumNumber(QuantumNumberType::J);
    const Rational Ji = UpperQuantumNumber(QuantumNumberType::J);
    const Rational N = LowerQuantumNumber(QuantumNumberType::N);
    constexpr Rational one(1, 1);
    const Numeric val = sqrt(6 * (2 * Jf + 1) * (2 * Ji + 1)) *
                        wigner6j(one, one, one, Ji, Jf, N);
    if ((Jf + N) % 2)
      return -val;
    else
      return +val;
  }

  Numeric G_upper() const { return mgupper; }

  Numeric G_lower() const { return mglower; }

  void SetQuantumNumberLower(const Index i, const Rational r) {
    mqid.LowerQuantumNumbers().Set(i, r);
  }
  void SetQuantumNumberLower(const String i, const Rational r) {
    mqid.LowerQuantumNumbers().Set(i, r);
  }
  void SetQuantumNumberLower(const QuantumNumberType i, const Rational r) {
    mqid.LowerQuantumNumbers().Set(i, r);
  }
  void SetQuantumNumberUpper(const QuantumNumberType i, const Rational r) {
    mqid.UpperQuantumNumbers().Set(i, r);
  }
  void SetQuantumNumberUpper(const String i, const Rational r) {
    mqid.UpperQuantumNumbers().Set(i, r);
  }
  void SetQuantumNumberUpper(const Index i, const Rational r) {
    mqid.UpperQuantumNumbers().Set(i, r);
  }

  const QuantumIdentifier& QuantumIdentity() const { return mqid; }
  const QuantumNumbers& LowerQuantumNumbers() const {
    return mqid.LowerQuantumNumbers();
  }
  const QuantumNumbers& UpperQuantumNumbers() const {
    return mqid.UpperQuantumNumbers();
  }
  Rational LowerQuantumNumber(QuantumNumberType X) const {
    return mqid.LowerQuantumNumber(X);
  }
  Rational UpperQuantumNumber(QuantumNumberType X) const {
    return mqid.UpperQuantumNumber(X);
  }
  bool InQuantumID(const QuantumIdentifier& qid) const { return mqid.In(qid); }
  bool UpperStateInQuantumID(const QuantumIdentifier& qid) const {
    return mqid.InUpper(qid);
  }
  bool LowerStateInQuantumID(const QuantumIdentifier& qid) const {
    return mqid.InLower(qid);
  }

  bool do_linemixing(const Numeric& P) const noexcept {
    if (mlinemixing_limit < 0)
      return true;
    else
      return mlinemixing_limit > P;
  }

  void SetLineMixingParameters(
      Numeric& Y,
      Numeric& G,
      Numeric& DV,
      const Numeric& T,
      const Numeric& P,
      const ConstVectorView vmrs,
      const ArrayOfArrayOfSpeciesTag& abs_species) const {
    if (do_linemixing(P)) {
      const auto v = mlineshapemodel.vmrs(vmrs, abs_species, mqid);
      Y = mlineshapemodel.Y(T, mti0, P, v);
      G = mlineshapemodel.G(T, mti0, P, v);
      DV = mlineshapemodel.DV(T, mti0, P, v);
    } else
      Y = G = DV = 0;
  }

  void SetPressureBroadeningParameters(
      Numeric& G0,
      Numeric& G2,
      Numeric& ETA,
      Numeric& D0,
      Numeric& D2,
      Numeric& FVC,
      const Numeric& T,
      const Numeric& P,
      const ConstVectorView vmrs,
      const ArrayOfArrayOfSpeciesTag& abs_species) const {
    const auto v = mlineshapemodel.vmrs(vmrs, abs_species, mqid);
    G0 = mlineshapemodel.G0(T, mti0, P, v);
    D0 = mlineshapemodel.D0(T, mti0, P, v);
    G2 = mlineshapemodel.G2(T, mti0, P, v);
    D2 = mlineshapemodel.D2(T, mti0, P, v);
    FVC = mlineshapemodel.FVC(T, mti0, P, v);
    ETA = mlineshapemodel.ETA(T, mti0, P, v);
  }

  void SetPressureBroadeningParametersTemperatureDerivative(
      Numeric& dG0,
      Numeric& dG2,
      Numeric& dETA,
      Numeric& dD0,
      Numeric& dD2,
      Numeric& dFVC,
      const Numeric& T,
      const Numeric& P,
      const ConstVectorView vmrs,
      const ArrayOfArrayOfSpeciesTag& abs_species) const {
    const auto v = mlineshapemodel.vmrs(vmrs, abs_species, mqid);
    dG0 = mlineshapemodel.dG0_dT(T, mti0, P, v);
    dD0 = mlineshapemodel.dD0_dT(T, mti0, P, v);
    dG2 = mlineshapemodel.dG2_dT(T, mti0, P, v);
    dD2 = mlineshapemodel.dD2_dT(T, mti0, P, v);
    dFVC = mlineshapemodel.dFVC_dT(T, mti0, P, v);
    dETA = mlineshapemodel.dETA_dT(T, mti0, P, v);
  }

  void SetLineMixing2SecondOrderData(ConstVectorView d) {
    mlineshapemodel.SetLineMixingModel(
        LineShape::LegacyLineMixingData::vector2modellm(
            d, LineShape::LegacyLineMixingData::TypeLM::LM_2NDORDER)
            .Data()[0]);
  }

  void SetLineMixing2AER(ConstVectorView d) {
    for (auto& sm : mlineshapemodel.Data())
      for (Index i = 0; i < d.nelem(); i++) sm.Interp()[i] = d[i];
  }

  LineShape::Model& GetLineShapeModel() { return mlineshapemodel; }
  const LineShape::Model& GetLineShapeModel() const { return mlineshapemodel; }

  void SetLineShapeModel(const LineShape::Model& lsm) { mlineshapemodel = lsm; }

  void LineShapeModelOnlyAir() {
    if (mlineshapemodel.Bath())
      mlineshapemodel = LineShape::Model(mlineshapemodel.Data().back(),
                                         mlineshapemodel.ModelType());
    else
      throw std::runtime_error(
          "Cannot keep only pure air since no air broadening exists");
  }

  void LineShapeModelRemoveSelf() {
    if (mlineshapemodel.Self()) mlineshapemodel.Remove(0);
  }

  void SetLineShapeModelParameters(const LineShape::ModelParameters& x,
                                   const String& spec,
                                   const String& var) {
    mlineshapemodel.Set(x, spec, LineShape::string2variable(var));
  }

  void SetLineShapeModelParameter(const Numeric& X,
                                  const String& spec,
                                  const String& var,
                                  const String& coeff) {
    auto x = GetLineShapeModelParameters(spec, var);
    LineShape::SingleModelParameter(x, coeff) = X;
    SetLineShapeModelParameters(x, spec, var);
  }

  LineShape::ModelParameters GetLineShapeModelParameters(
      const String& spec, const String& var) const {
    return mlineshapemodel.Get(spec, LineShape::string2variable(var));
  }

  void SetSpecial() { mstandard = false; }
  void SetStandard() { mstandard = true; }

  bool LineMixingByBand() const { return not mstandard; }

  bool LineShapeModelHasAir() const noexcept { return mlineshapemodel.Bath(); }

  void SetAirPressureBroadening(Numeric& G0,
                                Numeric& D0,
                                const Numeric& T,
                                const Numeric& P,
                                const Numeric& self_vmr) const {
    G0 = P * (mlineshapemodel.Data().front().compute(
                  T, mti0, LineShape::Variable::G0) *
                  self_vmr +
              mlineshapemodel.Data().back().compute(
                  T, mti0, LineShape::Variable::G0) *
                  (1 - self_vmr));
    D0 = P * (mlineshapemodel.Data().front().compute(
                  T, mti0, LineShape::Variable::D0) *
                  self_vmr +
              mlineshapemodel.Data().back().compute(
                  T, mti0, LineShape::Variable::D0) *
                  (1 - self_vmr));
  }

  Numeric PressureBroadeningAirBroadeningNair() const { return Nair(); }
  Numeric PressureBroadeningAirBroadeningPsf() const { return Psf(); }
  Numeric PressureBroadeningAirBroadeningAgam() const { return Agam(); }

  Numeric GetPrepInternalDerivative(const Numeric& T,
                                    const Numeric& P,
                                    const Vector& vmrs,
                                    const RetrievalQuantity& derivative) const
      noexcept {
    const bool self = derivative.Mode() == LineShape::self_broadening;
    const bool bath = derivative.Mode() == LineShape::bath_broadening;
    if (mlineshapemodel.Self() and self)
      return mlineshapemodel.GetInternalDeriv(
          T, mti0, P, 0, vmrs, derivative.PropMatType());
    else if (self) {
      const auto pos = mlineshapemodel.this_species(mqid);
      return mlineshapemodel.GetInternalDeriv(
          T, mti0, P, pos, vmrs, derivative.PropMatType());
    } else if (mlineshapemodel.Bath() and bath)
      return mlineshapemodel.GetInternalDeriv(T,
                                              mti0,
                                              P,
                                              mlineshapemodel.nelem() - 1,
                                              vmrs,
                                              derivative.PropMatType());
    else if (bath)
      return 0;
    else {
      const auto pos =
          mlineshapemodel.this_species(derivative.QuantumIdentity());
      return mlineshapemodel.GetInternalDeriv(
          T, mti0, P, pos, vmrs, derivative.PropMatType());
    }
  }

  Numeric GetInternalDerivative(const Numeric& T,
                                const Numeric& P,
                                const ConstVectorView vmrs,
                                const ArrayOfArrayOfSpeciesTag& abs_species,
                                const RetrievalQuantity& derivative) const {
    auto v = mlineshapemodel.vmrs(vmrs, abs_species, mqid);
    return GetPrepInternalDerivative(T, P, v, derivative);
  }

  LineShape::Output GetPrepShapeParams(const Numeric& T,
                                       const Numeric& P,
                                       const Vector& vmrs) const noexcept {
    auto x = mlineshapemodel.GetParams(T, mti0, P, vmrs);
    if (not do_linemixing(P)) x.Y = x.G = x.DV = 0;
    return x;
  }

  LineShape::Output GetShapeParams(
      const Numeric& T,
      const Numeric& P,
      const ConstVectorView vmrs,
      const ArrayOfArrayOfSpeciesTag& abs_species) const {
    auto v = mlineshapemodel.vmrs(vmrs, abs_species, mqid);
    return GetPrepShapeParams(T, P, v);
  }

  LineShape::Output GetPrepShapeParams_dT(const Numeric& T,
                                          const Numeric& P,
                                          const Vector& vmrs) const noexcept {
    auto x = mlineshapemodel.GetTemperatureDerivs(T, mti0, P, vmrs);
    if (not do_linemixing(P)) x.Y = x.G = x.DV = 0;
    return x;
  }

  LineShape::Output GetShapeParams_dT(
      const Numeric& T,
      const Numeric& P,
      const ConstVectorView vmrs,
      const ArrayOfArrayOfSpeciesTag& abs_species) const {
    auto v = mlineshapemodel.vmrs(vmrs, abs_species, mqid);
    return GetPrepShapeParams_dT(T, P, v);
  }

  LineShape::Output GetShapeParams_dVMR(const Numeric& T,
                                        const Numeric& P,
                                        const QuantumIdentifier& vmr_qi) const
      noexcept {
    const bool self = vmr_qi.Species() == mqid.Species();
    if (mlineshapemodel.Self() and self) {
      auto x = mlineshapemodel.GetVMRDerivs(T, mti0, P, 0);

      if (mlineshapemodel.Bath())
        x = LineShape::differenceOutput(
            x,
            mlineshapemodel.GetVMRDerivs(
                T, mti0, P, mlineshapemodel.nelem() - 1));

      if (not do_linemixing(P)) x.Y = x.G = x.DV = 0;

      return x;
    } else if (mlineshapemodel.Bath() and self)
      return {0, 0, 0, 0, 0, 0, 0, 0, 0};
    else {
      const auto pos = mlineshapemodel.this_species(vmr_qi);
      auto x = mlineshapemodel.GetVMRDerivs(T, mti0, P, pos);

      if (mlineshapemodel.Bath())
        x = LineShape::differenceOutput(
            x,
            mlineshapemodel.GetVMRDerivs(
                T, mti0, P, mlineshapemodel.nelem() - 1));

      if (not do_linemixing(P)) x.Y = x.G = x.DV = 0;
      return x;
    }
  }

  void ZeemanModelInit() { mzeemanmodel = Zeeman::Model(mqid); }
  void ZeemanModelInitZero() { mzeemanmodel = Zeeman::Model(); }
  Numeric ZeemanModelStrength(Zeeman::Polarization type, Index i) const {
    return mzeemanmodel.Strength(UpperQuantumNumber(QuantumNumberType::J),
                                 LowerQuantumNumber(QuantumNumberType::J),
                                 type,
                                 i);
  }
  Numeric ZeemanModelSplitting(Zeeman::Polarization type, Index i) const {
    return mzeemanmodel.Splitting(UpperQuantumNumber(QuantumNumberType::J),
                                  LowerQuantumNumber(QuantumNumberType::J),
                                  type,
                                  i);
  }
  Index ZeemanModelLineCount(Zeeman::Polarization type) const {
    return Zeeman::nelem(UpperQuantumNumber(QuantumNumberType::J),
                         LowerQuantumNumber(QuantumNumberType::J),
                         type);
  }
  Zeeman::Model ZeemanModel() const { return mzeemanmodel; }
  Zeeman::Model& ZeemanModel() { return mzeemanmodel; }

  const Numeric& CutOff() const { return mcutoff; }
  void SetCutOff(const Numeric& cutoff) { mcutoff = cutoff; }

  const Numeric& LineMixingLimit() const { return mlinemixing_limit; }
  void SetLineMixingLimit(const Numeric& limit) { mlinemixing_limit = limit; }

  const MirroringType& GetMirroringType() const { return mmirroring; }
  void SetMirroringType(const MirroringType in) { mmirroring = in; }
  String GetMirroringTypeString() const;

  const LineNormalizationType& GetLineNormalizationType() const {
    return mlinenorm;
  }
  void SetLineNormalizationType(const LineNormalizationType in) {
    mlinenorm = in;
  }
  String GetLineNormalizationTypeString() const;

  LineShape::Type GetLineShapeType() const noexcept {
    return mlineshapemodel.ModelType();
  }

  const LinePopulationType& GetLinePopulationType() const {
    return mpopulation;
  }
  void SetLinePopulationType(const LinePopulationType in) { mpopulation = in; }
  String GetLinePopulationTypeString() const;

  bool IsNotSameSpecIso(const LineRecord& other) const {
    return mqid.Species() not_eq other.mqid.Species() or
           mqid.Isotopologue() not_eq other.mqid.Isotopologue();
  }

  bool ReadFromHitran2001Stream(istream& is, const Verbosity& verbosity);

  bool ReadFromLBLRTMStream(istream& is, const Verbosity& verbosity);

  bool ReadFromHitran2004Stream(istream& is,
                                const Verbosity& verbosity,
                                const Numeric fmin = 0);

  bool ReadFromMytran2Stream(istream& is, const Verbosity& verbosity);

  bool ReadFromJplStream(istream& is, const Verbosity& verbosity);

  bool ReadFromArtscat3Stream(istream& is, const Verbosity& verbosity);

  bool ReadFromArtscat4Stream(istream& is, const Verbosity& verbosity);

  bool ReadFromArtscat5Stream(istream& is, const Verbosity& verbosity);

 private:
  // Version number:
  Index mversion;

  // Identity of line
  QuantumIdentifier mqid;

  // The line center frequency in Hz:
  Numeric mf;

  // The line intensity in m^2/Hz:
  Numeric mi0;

  // Reference temperature for I0 in K:
  Numeric mti0;

  // Lower state energy in Joules:
  Numeric melow;

  // Lower state vibrational energy in Joules:
  Numeric mevlow;  // NOTE: Not stored in binary data

  // Lower state vibrational energy index:
  Index mnlte_lower_index;  // NOTE: Not stored in binary data

  // Upper state vibrational energy in Joules:
  Numeric mevupp;  // NOTE: Not stored in binary data

  // Upper state vibrational energy index:
  Index mnlte_upper_index;  // NOTE: Not stored in binary data

  // Einstein A-coefficient in 1/s:
  Numeric ma;

  // Upper state stat. weight:
  Numeric mgupper;

  // Lower state stat. weight:
  Numeric mglower;

  LineShape::Model mlineshapemodel;

  bool mstandard;

  Zeeman::Model mzeemanmodel;

  Numeric mcutoff;

  Numeric mlinemixing_limit;

  MirroringType mmirroring;

  LineNormalizationType mlinenorm;

  LinePopulationType mpopulation;
};

ostream& operator<<(ostream& os, const LineRecord& lr);

//======================================================================
//         Typedefs for LineRecord Arrays
//======================================================================

typedef Array<LineRecord> ArrayOfLineRecord;

typedef Array<Array<LineRecord> > ArrayOfArrayOfLineRecord;

LinePopulationType LinePopulationTypeFromString(const String& in);

MirroringType MirroringTypeFromString(const String& in);

LineNormalizationType LineNormalizationTypeFromString(const String& in);

#endif  // linerecord_h