zeemandata.h

Go to the documentation of this file.

/* Copyright (C) 2018 Richard Larsson
 
   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 zeemandata_h
#define zeemandata_h
 
#include "constants.h"
#include "file.h"
#include "mystring.h"
#include "propagationmatrix.h"
#include "quantum.h"
#include "wigner_functions.h"
#include <limits>
 
namespace Zeeman {
  
enum class Polarization : char { SigmaMinus, Pi, SigmaPlus };
 
constexpr Index dM(Polarization type) noexcept {
  switch (type) {
    case Polarization::SigmaMinus:
      return -1;
    case Polarization::Pi:
      return 0;
    case Polarization::SigmaPlus:
      return 1;
  }
  return std::numeric_limits<Index>::max();
}
 
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept {
  switch (type) {
    case Polarization::SigmaMinus:
      if (Ju < Jl)
        return -Ju;
      else if (Ju == Jl)
        return -Ju + 1;
      else
        return -Ju + 2;
    case Polarization::Pi:
      return -min(Ju, Jl);
    case Polarization::SigmaPlus:
      return -Ju;
  }
  return Rational(std::numeric_limits<Index>::max());
}
 
constexpr Rational end(Rational Ju, Rational Jl, Polarization type) noexcept {
  switch (type) {
    case Polarization::SigmaMinus:
      return Ju + 1;
    case Polarization::Pi:
      return min(Ju, Jl);
    case Polarization::SigmaPlus:
      if (Ju < Jl)
        return Ju + 1;
      else if (Ju == Jl)
        return Ju;
      else
        return Jl;
  }
  return Rational(std::numeric_limits<Index>::max());
}
 
constexpr Index nelem(Rational Ju, Rational Jl, Polarization type) noexcept {
  return (end(Ju, Jl, type) - start(Ju, Jl, type)).toIndex() + 1;
}
 
constexpr Rational Mu(Rational Ju,
                      Rational Jl,
                      Polarization type,
                      Index n) noexcept {
  return start(Ju, Jl, type) + n;
}
 
 
constexpr Rational Ml(Rational Ju,
                      Rational Jl,
                      Polarization type,
                      Index n) noexcept {
  return Mu(Ju, Jl, type, n) + dM(type);
}
 
constexpr Numeric PolarizationFactor(Polarization type) noexcept {
  switch (type) {
    case Polarization::SigmaMinus:
      return .75;
    case Polarization::Pi:
      return 1.5;
    case Polarization::SigmaPlus:
      return .75;
  }
  return std::numeric_limits<Numeric>::max();
}
 
constexpr bool GoodHundData(const QuantumNumbers& qns) noexcept {
  if (qns[QuantumNumberType::Hund].isUndefined()) return false;
  switch (Hund(qns[QuantumNumberType::Hund].toIndex())) {
    case Hund::CaseA:
      if (qns[QuantumNumberType::Omega].isUndefined() or
          qns[QuantumNumberType::J].isUndefined() or
          qns[QuantumNumberType::Lambda].isUndefined() or
          qns[QuantumNumberType::S].isUndefined())
        return false;
      break;
    case Hund::CaseB:
      if (qns[QuantumNumberType::N].isUndefined() or
          qns[QuantumNumberType::J].isUndefined() or
          qns[QuantumNumberType::Lambda].isUndefined() or
          qns[QuantumNumberType::S].isUndefined())
        return false;
      break;
    default:
      return false;
  }
  return true;
}
 
constexpr Numeric SimpleGCaseB(Rational N,
                               Rational J,
                               Rational Lambda,
                               Rational S,
                               Numeric GS,
                               Numeric GL) noexcept {
  auto JJ = J * (J + 1);
  auto NN = N * (N + 1);
  auto SS = S * (S + 1);
  auto LL = Lambda * Lambda;
 
  if (JJ == 0)
    return 0.0;
  if (NN not_eq 0) {
    auto T1 = ((JJ + SS - NN) / JJ / 2).toNumeric();
    auto T2 = ((JJ - SS + NN) * LL / NN / JJ / 2).toNumeric();
    return GS * T1 + GL * T2;
  }      
  auto T1 = ((JJ + SS - NN) / JJ / 2).toNumeric();
  return GS * T1;
}
 
constexpr Numeric SimpleGCaseA(Rational Omega,
                               Rational J,
                               Rational Lambda,
                               Rational Sigma,
                               Numeric GS,
                               Numeric GL) noexcept {
  auto JJ = J * (J + 1);
 
  if (JJ == Rational(0))
    return 0.0;
  auto DIV = Omega / JJ;
  auto T1 = (Sigma * DIV).toNumeric();
  auto T2 = (Lambda * DIV).toNumeric();
  return GS * T1 + GL * T2;
 
}
 
constexpr Numeric SimpleG(const QuantumNumbers& qns,
                          const Numeric& GS,
                          const Numeric& GL) noexcept{
  if (not GoodHundData(qns))
    return NAN;
 
  switch (Hund(qns[QuantumNumberType::Hund].toIndex())) {
    case Hund::CaseA:
      return SimpleGCaseA(qns[QuantumNumberType::Omega],
                          qns[QuantumNumberType::J],
                          qns[QuantumNumberType::Lambda],
                          qns[QuantumNumberType::S],
                          GS,
                          GL);
    case Hund::CaseB:
      return SimpleGCaseB(qns[QuantumNumberType::N],
                          qns[QuantumNumberType::J],
                          qns[QuantumNumberType::Lambda],
                          qns[QuantumNumberType::S],
                          GS,
                          GL);
  }
 
  return NAN;
}
 
struct SplittingData {
  Numeric gu, gl;
};
 
class Model {
 private:
  SplittingData mdata;
 
 public:
   constexpr Model(SplittingData gs = {NAN, NAN}) noexcept : mdata(gs) {}
   
   constexpr Model(Numeric gu, Numeric gl) noexcept : Model(SplittingData{gu, gl}) {}
  
  explicit Model(const QuantumIdentifier& qid) noexcept;
 
  [[nodiscard]] /* constexpr */ bool empty() const noexcept {
    return std::isnan(mdata.gu) and std::isnan(mdata.gl);
  }
 
  constexpr Numeric& gu() noexcept { return mdata.gu; }
  
  constexpr Numeric& gl() noexcept { return mdata.gl; }
 
  constexpr void gu(Numeric x) noexcept { mdata.gu = x; }
  
  constexpr void gl(Numeric x) noexcept { mdata.gl = x; }
  
  [[nodiscard]] constexpr Numeric gu() const noexcept { return mdata.gu; }
  
  [[nodiscard]] constexpr Numeric gl() const noexcept { return mdata.gl; }
 
  [[nodiscard]] Numeric Strength(Rational Ju, Rational Jl, Polarization type, Index n) const;
  
  [[nodiscard]] constexpr Numeric Splitting(Rational Ju, Rational Jl, Polarization type, Index n) const
      noexcept {
    using Constant::bohr_magneton;
    using Constant::h;
    constexpr Numeric C = bohr_magneton / h;
 
    return C * (Ml(Ju, Jl, type, n).toNumeric() * gl() -
                Mu(Ju, Jl, type, n).toNumeric() * gu());
  }
 
  friend std::ostream& operator<<(std::ostream& os, const Model& m);
  
  friend std::istream& operator>>(std::istream& is, Model& m);
  
  friend std::ostream& operator<<(bofstream& bof, const Model& m);
  
  friend std::istream& operator>>(bifstream& bif, Model& m);
};  // Model;
 
Model GetSimpleModel(const QuantumIdentifier& qid) ARTS_NOEXCEPT;
 
Model GetAdvancedModel(const QuantumIdentifier& qid) ARTS_NOEXCEPT;
 
std::ostream& operator<<(std::ostream& os, const Model& m);
 
std::istream& operator>>(std::istream& is, Model& m);
 
std::ostream& operator<<(bofstream& bof, const Model& m);
 
std::istream& operator>>(bifstream& bif, Model& m);
 
class PolarizationVector {
 private:
  Eigen::RowVector4d att;  // attenuation vector
  Eigen::RowVector3d dis;  // dispersion vector
 
 public:
  PolarizationVector(Numeric a = 1,
                     Numeric b = 0,
                     Numeric c = 0,
                     Numeric d = 0,
                     Numeric u = 0,
                     Numeric v = 0,
                     Numeric w = 0) noexcept
      : att(a, b, c, d), dis(u, v, w){};
 
  [[nodiscard]] const Eigen::RowVector4d& attenuation() const noexcept { return att; }
  
  [[nodiscard]] const Eigen::RowVector3d& dispersion() const noexcept { return dis; }
  
  Eigen::RowVector4d& attenuation() noexcept { return att; }
  
  Eigen::RowVector3d& dispersion() noexcept { return dis; }
 
  [[nodiscard]] Eigen::Matrix4d matrix() const noexcept {
    return (Eigen::Matrix4d() << att[0],
            att[1],
            att[2],
            att[3],
            att[1],
            att[0],
            dis[0],
            dis[1],
            att[2],
            -dis[0],
            att[0],
            dis[2],
            att[3],
            -dis[1],
            -dis[2],
            att[0])
        .finished();
  }
};
 
struct AllPolarizationVectors {
  PolarizationVector sm, pi, sp;
};
 
AllPolarizationVectors AllPolarization(Numeric theta,
                                       Numeric eta) noexcept;
 
AllPolarizationVectors AllPolarization_dtheta(
    Numeric theta, const Numeric eta) noexcept;
 
AllPolarizationVectors AllPolarization_deta(Numeric theta,
                                            Numeric eta) noexcept;
 
const PolarizationVector& SelectPolarization(
    const AllPolarizationVectors& data, Polarization type) noexcept;
 
void sum(PropagationMatrix& pm, const ComplexVectorView& abs, const PolarizationVector& polvec, const bool do_phase=true) ARTS_NOEXCEPT;
 
void dsum(PropagationMatrix& dpm,
          const ComplexVectorView& abs,
          const ComplexVectorView& dabs,
          const PolarizationVector& polvec,
          const PolarizationVector& dpolvec_dtheta,
          const PolarizationVector& dpolvec_deta,
          const Numeric dH,
          const Numeric dtheta,
          const Numeric deta, const bool do_phase=true) ARTS_NOEXCEPT;
 
struct Derived {
  Numeric H, theta, eta, dH_du, dH_dv, dH_dw, dtheta_du, dtheta_dv, dtheta_dw,
      deta_du, deta_dv, deta_dw;
};
 
Derived FromGrids(Numeric u, Numeric v, Numeric w, Numeric z, Numeric a) noexcept;
 
constexpr Derived FromPreDerived(Numeric H,
                                 Numeric theta,
                                 Numeric eta) noexcept {
  return {H, theta, eta, 0, 0, 0, 0, 0, 0, 0, 0, 0};
}
};  // namespace Zeeman
 
// Typedef to make it easier to use
using ZeemanModel = Zeeman::Model;
 
#endif /* zeemandata_h */