quantum.h

Go to the documentation of this file.

/* Copyright (C) 2013
   Oliver Lemke <olemke@core-dump.info>

   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 quantum_h
#define quantum_h

#include <array>
#include <iostream>
#include <map>
#include <numeric>
#include <stdexcept>
#include "array.h"
#include "interpolation.h"
#include "matpack.h"
#include "mystring.h"
#include "rational.h"

enum class QuantumNumberType : Index {
  J = 0,   // Total angular momentum
  dJ,      // Delta total angular momentum
  M,       // Projection of J along magnetic field
  N,       // J minus spin
  dN,      // Delta J minus spin
  S,       // Spin angular momentum (from electrons) NOTE: S_global for HITRAN S
  tau,
  n,
  F,       // J + nuclear spin
  Ka,      //(This is a projection of J along one axis)
  Kc,      //(This is a projection of J along another axis)
  Omega,   // This is an absolute projection of J and S
  i,       //(Is related to Omega)
  Lambda,  // This is Sigma or Pi or Lambda states (as seen in literature)
  alpha,   // Alpha from HITRAN
  Sym,     // Symmetry expression
  parity,  // parity value (+/-)
  kronigParity,  // ???
  v1,      // Vibrational mode 1
  v2,      // Vibrational mode 2
  v3,      // Vibrational mode 3
  v4,      // Vibrational mode 4
  v5,      // Vibrational mode 5
  v6,      // Vibrational mode 6
  v7,
  v8,
  v9,
  v10,
  v11,
  v12,
  l1,      // The absolute sum of l_j for v_j
  l2,      // Vibrational angular momentum associated with v2
  l3,
  l4,
  l5,
  l6,
  l7,
  l8,
  l9,
  l10,
  l11,
  l12,
  pm,      // Symmetry type for l=0
  r,       // Rank of the level within a set of the same vibrational symmetry
  S_global,  // Symmetry of the level
  ElectronState,  // Electronic state
  n_global,  // Torosional quanta
  C,         // Another symmetry expression
  Hund,  // Flag for Hund case type.  This flag lets Zeeman know what to expect
  FINAL_ENTRY  // We need this to determine the number of elements in this enum
};

inline QuantumNumberType string2quantumnumbertype(const String& s) {
  #define INPUT_QUANTUM(ID) \
  if (s == #ID) return QuantumNumberType::ID
  INPUT_QUANTUM(J);
  else INPUT_QUANTUM(dJ);
  else INPUT_QUANTUM(M);
  else INPUT_QUANTUM(N);
  else INPUT_QUANTUM(dN);
  else INPUT_QUANTUM(tau);
  else INPUT_QUANTUM(n);
  else INPUT_QUANTUM(S);
  else INPUT_QUANTUM(F);
  else if (s.find("F#") < s.length()) return QuantumNumberType::F;  // HITRAN has many names for F
  else if (s == "K") return QuantumNumberType::Ka;  // HITRAN name
  else INPUT_QUANTUM(Ka);
  else INPUT_QUANTUM(Kc);
  else INPUT_QUANTUM(Omega);
  else INPUT_QUANTUM(i);
  else INPUT_QUANTUM(Lambda);
  else INPUT_QUANTUM(alpha);
  else INPUT_QUANTUM(Sym);
  else INPUT_QUANTUM(parity);
  else INPUT_QUANTUM(kronigParity);
  else if (s == "v") return QuantumNumberType::v1;  // HITRAN name
  else INPUT_QUANTUM(v1);
  else INPUT_QUANTUM(v2);
  else INPUT_QUANTUM(v3);
  else INPUT_QUANTUM(v4);
  else INPUT_QUANTUM(v5);
  else INPUT_QUANTUM(v6);
  else INPUT_QUANTUM(v7);
  else INPUT_QUANTUM(v8);
  else INPUT_QUANTUM(v9);
  else INPUT_QUANTUM(v10);
  else INPUT_QUANTUM(v11);
  else INPUT_QUANTUM(v12);
  else if (s == "l") return QuantumNumberType::l1;  // HITRAN name
  else INPUT_QUANTUM(l1);
  else INPUT_QUANTUM(l2);
  else INPUT_QUANTUM(l3);
  else INPUT_QUANTUM(l4);
  else INPUT_QUANTUM(l5);
  else INPUT_QUANTUM(l6);
  else INPUT_QUANTUM(l7);
  else INPUT_QUANTUM(l8);
  else INPUT_QUANTUM(l9);
  else INPUT_QUANTUM(l10);
  else INPUT_QUANTUM(l11);
  else INPUT_QUANTUM(l12);
  else INPUT_QUANTUM(pm);
  else INPUT_QUANTUM(r);
  else INPUT_QUANTUM(S_global);
  else INPUT_QUANTUM(ElectronState);
  else if (s == "ElecStateLabel") return QuantumNumberType::ElectronState;  // HITRAN name
  else INPUT_QUANTUM(n_global);
  else INPUT_QUANTUM(C);
  else INPUT_QUANTUM(Hund);
  else return QuantumNumberType::FINAL_ENTRY;
  #undef INPUT_QUANTUM
}

inline String quantumnumbertype2string(QuantumNumberType s) {
  #define INPUT_QUANTUM(ID) \
  if (s == QuantumNumberType::ID) return #ID
  INPUT_QUANTUM(J);
  else INPUT_QUANTUM(dJ);
  else INPUT_QUANTUM(M);
  else INPUT_QUANTUM(N);
  else INPUT_QUANTUM(dN);
  else INPUT_QUANTUM(tau);
  else INPUT_QUANTUM(n);
  else INPUT_QUANTUM(S);
  else INPUT_QUANTUM(F);
  else INPUT_QUANTUM(Ka);
  else INPUT_QUANTUM(Kc);
  else INPUT_QUANTUM(Omega);
  else INPUT_QUANTUM(i);
  else INPUT_QUANTUM(Lambda);
  else INPUT_QUANTUM(alpha);
  else INPUT_QUANTUM(Sym);
  else INPUT_QUANTUM(parity);
  else INPUT_QUANTUM(kronigParity);
  else INPUT_QUANTUM(v1);
  else INPUT_QUANTUM(v2);
  else INPUT_QUANTUM(v3);
  else INPUT_QUANTUM(v4);
  else INPUT_QUANTUM(v5);
  else INPUT_QUANTUM(v6);
  else INPUT_QUANTUM(v7);
  else INPUT_QUANTUM(v8);
  else INPUT_QUANTUM(v9);
  else INPUT_QUANTUM(v10);
  else INPUT_QUANTUM(v11);
  else INPUT_QUANTUM(v12);
  else INPUT_QUANTUM(l1);
  else INPUT_QUANTUM(l2);
  else INPUT_QUANTUM(l3);
  else INPUT_QUANTUM(l4);
  else INPUT_QUANTUM(l5);
  else INPUT_QUANTUM(l6);
  else INPUT_QUANTUM(l7);
  else INPUT_QUANTUM(l8);
  else INPUT_QUANTUM(l9);
  else INPUT_QUANTUM(l10);
  else INPUT_QUANTUM(l11);
  else INPUT_QUANTUM(l12);
  else INPUT_QUANTUM(pm);
  else INPUT_QUANTUM(r);
  else INPUT_QUANTUM(S_global);
  else INPUT_QUANTUM(ElectronState);
  else INPUT_QUANTUM(n_global);
  else INPUT_QUANTUM(C);
  else INPUT_QUANTUM(Hund);
  throw std::runtime_error("Bad quantum number type");
  #undef INPUT_QUANTUM
}

enum class Hund : Index { CaseA = int('a'), CaseB = int('b') };

class QuantumNumbers {
 public:
  typedef std::array<Rational, Index(QuantumNumberType::FINAL_ENTRY)>
      QuantumContainer;

  constexpr QuantumNumbers() noexcept
      : mqnumbers({RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 3
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 6
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 9
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 12
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 15
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 18
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 21
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 24
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 27
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 30
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 33
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 36
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 39
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 42
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 45
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 48
                   RATIONAL_UNDEFINED}) {}
  
  // Temporary initialization until there is a more reliable way to setup 
  // a full initialization for all quantum numbers but only choose a few
  // select ones based on the problem at hand
  constexpr QuantumNumbers(Rational J, Rational N, Rational v) noexcept
      : mqnumbers({J, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 3
                   N, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 6
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 9
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 12
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 15
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 18
                   v, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 21
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 24
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 27
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 30
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 33
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 36
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 39
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 42
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 45
                   RATIONAL_UNDEFINED, RATIONAL_UNDEFINED, RATIONAL_UNDEFINED,  // 48
                   RATIONAL_UNDEFINED}) {}

  constexpr Rational operator[](const Index qn) const noexcept {
    return mqnumbers[qn];
  }

  constexpr Rational operator[](const QuantumNumberType qn) const noexcept {
    return mqnumbers[Index(qn)];
  }
  
  Rational& operator[](const Index qn) noexcept {
    return mqnumbers[qn];
  }
  
  Rational& operator[](const QuantumNumberType qn) noexcept {
    return mqnumbers[Index(qn)];
  }

  void Set(Index qn, Rational r) {
    assert(qn < Index(QuantumNumberType::FINAL_ENTRY));
    mqnumbers[qn] = r;
  }

  void Set(QuantumNumberType qn, Rational r) {
    assert(qn != QuantumNumberType::FINAL_ENTRY);
    mqnumbers[Index(qn)] = r;
  }

  void Set(String qn, Rational r) {
    mqnumbers[Index(string2quantumnumbertype(qn))] = r;
  }

  const QuantumContainer& GetNumbers() const { return mqnumbers; }

  Index nNumbers() const {
    return std::accumulate(
        mqnumbers.cbegin(), mqnumbers.cend(), 0, [](Index i, Rational r) {
          return r.isUndefined() ? i : i + 1;
        });
  }

  bool Compare(const QuantumNumbers& qn) const;

 private:
  QuantumContainer mqnumbers;
};

typedef Array<QuantumNumbers> ArrayOfQuantumNumbers;

class QuantumIdentifier {
 public:
  typedef enum : Index { TRANSITION, ENERGY_LEVEL, ALL, NONE } QType;

  constexpr QuantumIdentifier() noexcept
      : mqtype(QuantumIdentifier::NONE), mspecies(-1), miso(-1) {}

  constexpr QuantumIdentifier(const QuantumIdentifier::QType& qt,
                              const Index species,
                              const Index iso) noexcept
      : mqtype(qt), mspecies(species), miso(iso) {}

  constexpr QuantumIdentifier(const Index spec,
                              const Index isot,
                              const QuantumNumbers& upper,
                              const QuantumNumbers& lower) noexcept
      : mqtype(QuantumIdentifier::TRANSITION),
        mspecies(spec),
        miso(isot),
        mqm({upper, lower}) {}

  QuantumIdentifier(const Index spec,
                    const Index isot,
                    const std::vector<QuantumNumberType>& keys,
                    const std::vector<Rational>& upper,
                    const std::vector<Rational>& lower)
      : mqtype(QuantumIdentifier::TRANSITION),
        mspecies(spec),
        miso(isot) {
          for(size_t i=0; i<keys.size(); i++) {
            mqm[TRANSITION_UPPER_INDEX][keys[i]] = upper[i];
            mqm[TRANSITION_LOWER_INDEX][keys[i]] = lower[i];
          }
        }

  constexpr QuantumIdentifier(const Index spec,
                              const Index isot,
                              const QuantumNumbers& qnr) noexcept
      : mqtype(QuantumIdentifier::ENERGY_LEVEL),
        mspecies(spec),
        miso(isot),
        mqm({qnr}) {}

  explicit QuantumIdentifier(String x) { SetFromString(x); }

  static constexpr Index TRANSITION_UPPER_INDEX = 0;

  static constexpr Index TRANSITION_LOWER_INDEX = 1;

  static constexpr Index ENERGY_LEVEL_INDEX = 0;

  void Species(Index sp) { mspecies = sp; }

  void Isotopologue(Index iso) { miso = iso; }

  void SetTransition(const QuantumNumbers& upper, const QuantumNumbers& lower);

  void SetEnergyLevel(const QuantumNumbers& q);

  void SetAll() { mqtype = QuantumIdentifier::ALL; };
  
  void SetNone() { mqtype = QuantumIdentifier::NONE; };

  void SetTransition() { mqtype = QuantumIdentifier::TRANSITION; };

  void SetFromString(String str);

  void SetFromStringForCO2Band(String upper, String lower, String iso);

  constexpr QType Type() const { return mqtype; }
  
  void Type(QType x) { mqtype = x; }
  
  static bool validIndexForType(Index x) noexcept {
    constexpr auto keys = stdarrayify(Index(TRANSITION), ENERGY_LEVEL, ALL, NONE);
    return std::any_of(keys.cbegin(), keys.cend(), [x](auto y){return x == y;});
  }
  
  static QType string2Type(const String& str) noexcept {
    if (std::string("ENERGY_LEVEL") == str) {
      return QuantumIdentifier::ENERGY_LEVEL;
    } else if (std::string("TRANSITION") == str) {
      return QuantumIdentifier::TRANSITION;
    } else if (std::string("ALL") == str) {
      return QuantumIdentifier::ALL;
    } else if (std::string("NONE") == str) {
      return QuantumIdentifier::NONE;
    } else {
      return QType(-1);
    }
  }

  String TypeStr() const;

  String SpeciesName() const;

  constexpr Index Species() const { return mspecies; }

  Index& Species() { return mspecies; }

  constexpr Index Isotopologue() const { return miso; }

  Index& Isotopologue() { return miso; }

  const std::array<QuantumNumbers, 2>& QuantumMatch() const { return mqm; }

  std::array<QuantumNumbers, 2>& QuantumMatch() { return mqm; }

  constexpr QuantumIdentifier UpperQuantumId() const noexcept {
    return QuantumIdentifier(mspecies, miso, mqm[TRANSITION_UPPER_INDEX]);
  };

  constexpr QuantumIdentifier LowerQuantumId() const noexcept {
    return QuantumIdentifier(mspecies, miso, mqm[TRANSITION_LOWER_INDEX]);
  };

  const QuantumNumbers& UpperQuantumNumbers() const noexcept {
    return mqm[TRANSITION_UPPER_INDEX];
  };

  const QuantumNumbers& LowerQuantumNumbers() const noexcept {
    return mqm[TRANSITION_LOWER_INDEX];
  };

  constexpr Rational UpperQuantumNumber(QuantumNumberType X) const noexcept {
    return mqm[TRANSITION_UPPER_INDEX][X];
  };

  constexpr Rational LowerQuantumNumber(QuantumNumberType X) const noexcept {
    return mqm[TRANSITION_LOWER_INDEX][X];
  };
  
  Rational& UpperQuantumNumber(QuantumNumberType X) noexcept {
    return mqm[TRANSITION_UPPER_INDEX][X];
  };
  
  Rational& LowerQuantumNumber(QuantumNumberType X) noexcept {
    return mqm[TRANSITION_LOWER_INDEX][X];
  };

  const QuantumNumbers& EnergyLevelQuantumNumbers() const noexcept {
    return mqm[ENERGY_LEVEL_INDEX];
  }
  
  constexpr Rational EnergyLevelQuantumNumber(QuantumNumberType X) const noexcept {
    return mqm[ENERGY_LEVEL_INDEX][X];
  };

  QuantumNumbers& UpperQuantumNumbers() {
    return mqm[TRANSITION_UPPER_INDEX];
  };

  QuantumNumbers& LowerQuantumNumbers() {
    return mqm[TRANSITION_LOWER_INDEX];
  };

  QuantumNumbers& EnergyLevelQuantumNumbers() {
    return mqm[ENERGY_LEVEL_INDEX];
  }

  bool In(const QuantumIdentifier& other) const;

  bool InLower(const QuantumIdentifier& other) const;

  bool InUpper(const QuantumIdentifier& other) const;

  bool any_quantumnumbers() const;

  bool IsEnergyLevelType() const { return mqtype == ENERGY_LEVEL; }

 private:
  QType mqtype;
  Index mspecies;
  Index miso;
  std::array<QuantumNumbers, 2> mqm;
};

inline bool operator==(const QuantumIdentifier& a, const QuantumIdentifier& b) {
  if (!(a.Isotopologue() == b.Isotopologue() && a.Species() == b.Species() &&
        a.Type() == b.Type()))
    return false;

  if (a.Type() == QuantumIdentifier::ENERGY_LEVEL)
    return a.QuantumMatch()[a.ENERGY_LEVEL_INDEX].Compare(
        b.QuantumMatch()[b.ENERGY_LEVEL_INDEX]);
  else if (a.Type() == QuantumIdentifier::TRANSITION)
    return a.QuantumMatch()[a.TRANSITION_LOWER_INDEX].Compare(
               b.QuantumMatch()[b.TRANSITION_LOWER_INDEX]) &&
           a.QuantumMatch()[a.TRANSITION_UPPER_INDEX].Compare(
               b.QuantumMatch()[b.TRANSITION_UPPER_INDEX]);
  else if (a.Type() == QuantumIdentifier::ALL)
    return true;
  else if (a.Type() == QuantumIdentifier::NONE)
    return false;
  else
    throw std::runtime_error("Programmer error --- added type is missing");
}

inline bool operator!=(const QuantumIdentifier& a, const QuantumIdentifier& b) {
  return not operator==(a, b);
}

inline bool operator==(const QuantumNumbers& a, const QuantumNumbers& b) {
  return a.Compare(b) and b.Compare(a);
}

typedef Array<QuantumIdentifier> ArrayOfQuantumIdentifier;

bool IsValidQuantumNumberName(String name);

void ThrowIfQuantumNumberNameInvalid(String name);

std::istream& operator>>(std::istream& is, QuantumNumbers& qn);

std::ostream& operator<<(std::ostream& os, const QuantumNumbers& qn);

std::ostream& operator<<(std::ostream& os, const QuantumIdentifier& qi);

std::ostream& operator<<(std::ostream& os, QuantumNumberType t);

void update_id(QuantumIdentifier& qid, const std::vector<std::array<String, 2> >& upper_list, const std::vector<std::array<String, 2> >& lower_list);

#endif