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 <limits>
#include "constants.h"
#include "file.h"
#include "mystring.h"
#include "quantum.h"
#include "wigner_functions.h"

namespace Zeeman {
  
enum class Polarization { 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 -std::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 std::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;
  else 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;
  } else {
    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 == 0)
    return 0.0;
  else {
    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) {}
  
  explicit Model(const QuantumIdentifier& qid) noexcept;

  /* constexpr */ bool empty() const noexcept {
    return std::isnan(mdata.gu) and std::isnan(mdata.gl);
  }

  Numeric& gu() noexcept { return mdata.gu; }
  
  Numeric& gl() noexcept { return mdata.gl; }

  void gu(Numeric x) noexcept { mdata.gu = x; }
  
  void gl(Numeric x) noexcept { mdata.gl = x; }
  
  constexpr Numeric gu() const noexcept { return mdata.gu; }
  
  constexpr Numeric gl() const noexcept { return mdata.gl; }

  Numeric Strength(Rational Ju, Rational Jl, Polarization type, Index n) const {
    using Constant::pow2;

    auto ml = Ml(Ju, Jl, type, n);
    auto mu = Mu(Ju, Jl, type, n);
    auto dm = Rational(dM(type));
    return PolarizationFactor(type) * pow2(wigner3j(Jl, Rational(1), Ju, ml, -dm, -mu));
  }
  
  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 inline std::ostream& operator<<(std::ostream& os, const Model& m);
  
  friend inline std::istream& operator>>(std::istream& is, Model& m);
  
  friend inline std::ostream& operator<<(bofstream& bof, const Model& m);
  
  friend inline std::istream& operator>>(bifstream& bif, Model& m);
};  // Model;

Model GetSimpleModel(const QuantumIdentifier& qid) noexcept;

Model GetAdvancedModel(const QuantumIdentifier& qid) noexcept;

inline std::ostream& operator<<(std::ostream& os, const Model& m) {
  os << m.mdata.gu << ' ' << m.mdata.gl;
  return os;
}

inline std::istream& operator>>(std::istream& is, Model& m) {
  is >> double_imanip() >> m.mdata.gu >> m.mdata.gl;
  return is;
}

inline std::ostream& operator<<(bofstream& bof, const Model& m) {
  bof << m.mdata.gu << m.mdata.gl;
  return bof;
}

inline std::istream& operator>>(bifstream& bif, Model& m) {
  bif >> m.mdata.gu >> m.mdata.gl;
  return bif;
}

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){};

  const Eigen::RowVector4d& attenuation() const noexcept { return att; }
  
  const Eigen::RowVector3d& dispersion() const noexcept { return dis; }
  
  Eigen::RowVector4d& attenuation() noexcept { return att; }
  
  Eigen::RowVector3d& dispersion() noexcept { return dis; }

  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;
};

inline AllPolarizationVectors AllPolarization(Numeric theta,
                                              Numeric eta) noexcept {
  const Numeric ST = std::sin(theta), CT = std::cos(theta), ST2 = ST * ST,
                CT2 = CT * CT, ST2C2E = ST2 * std::cos(2 * eta),
                ST2S2E = ST2 * std::sin(2 * eta);

  AllPolarizationVectors pv;
  pv.sm = PolarizationVector(
      1 + CT2, ST2C2E, ST2S2E, 2 * CT, 4 * CT, 2 * ST2S2E, -2 * ST2C2E);
  pv.pi =
      PolarizationVector(ST2, -ST2C2E, -ST2S2E, 0, 0, -2 * ST2S2E, 2 * ST2C2E);
  pv.sp = PolarizationVector(
      1 + CT2, ST2C2E, ST2S2E, -2 * CT, -4 * CT, 2 * ST2S2E, -2 * ST2C2E);
  return pv;
}

inline AllPolarizationVectors AllPolarization_dtheta(
    Numeric theta, const Numeric eta) noexcept {
  const Numeric ST = std::sin(theta), CT = std::cos(theta),
                C2E = std::cos(2 * eta), S2E = std::sin(2 * eta), dST = CT,
                dST2 = 2 * ST * dST, dCT = -ST, dST2C2E = dST2 * C2E,
                dST2S2E = dST2 * S2E, dCT2 = 2 * CT * dCT;

  AllPolarizationVectors pv;
  pv.sm = PolarizationVector(
      dCT2, dST2C2E, dST2S2E, 2 * dCT, 4 * dCT, 2 * dST2S2E, -2 * dST2C2E);
  pv.pi = PolarizationVector(
      dST2, -dST2C2E, -dST2S2E, 0, 0, -2 * dST2S2E, 2 * dST2C2E);
  pv.sp = PolarizationVector(
      dCT2, dST2C2E, dST2S2E, -2 * dCT, -4 * dCT, 2 * dST2S2E, -2 * dST2C2E);
  return pv;
}

inline AllPolarizationVectors AllPolarization_deta(Numeric theta,
                                                   Numeric eta) noexcept {
  const Numeric ST = std::sin(theta), ST2 = ST * ST, C2E = std::cos(2 * eta),
                S2E = std::sin(2 * eta), dST2C2E = -2 * ST2 * S2E,
                dST2S2E = 2 * ST2 * C2E;

  AllPolarizationVectors pv;
  pv.sm =
      PolarizationVector(0, dST2C2E, dST2S2E, 0, 0, 2 * dST2S2E, -2 * dST2C2E);
  pv.pi = PolarizationVector(
      0, -dST2C2E, -dST2S2E, 0, 0, -2 * dST2S2E, 2 * dST2C2E);
  pv.sp =
      PolarizationVector(0, dST2C2E, dST2S2E, 0, 0, 2 * dST2S2E, -2 * dST2C2E);
  return pv;
}

inline const PolarizationVector& SelectPolarization(
    const AllPolarizationVectors& data, Polarization type) noexcept {
      switch (type) {
    case Polarization::SigmaMinus:
      return data.sm;
    case Polarization::Pi:
      return data.pi;
    case Polarization::SigmaPlus:
      return data.sp;
  }
  std::terminate();
}

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
typedef Zeeman::Model ZeemanModel;

#endif /* zeemandata_h */