ContinuumState.hpp

#pragma once
#include "AdamsMoulton.hpp"
#include <utility>
#include <vector>
class DiracSpinor;
class Grid;
 
namespace DiracODE {
 
 
void solveContinuum(DiracSpinor &Fa, double en, const std::vector<double> &v,
                    double alpha, const DiracSpinor *const VxFa = nullptr,
                    const DiracSpinor *const Fa0 = nullptr);
 
double analytic_f_amplitude(double en, double alpha);
 
 
std::pair<double, double> numerical_f_amplitude(double en, int kappa,
                                                double alpha, double Zeff,
                                                double f_final, double g_final,
                                                double r_final, double dr);
 
struct DiracContinuumDerivative
    : AdamsMoulton::DerivativeMatrix<double, double> {
 
  DiracContinuumDerivative(double in_Zeff, const int in_kappa,
                           const double in_en, const double in_alpha)
      : Zeff(in_Zeff),
        kappa(in_kappa),
        en(in_en),
        alpha(in_alpha),
        cc(1.0 / alpha) {}
 
  double Zeff = 1.0;
  int kappa;
  double en, alpha, cc;
 
  double a(double r) const final { return double(-kappa) / r; }
  double b(double r) const final {
    return (alpha * en + 2.0 * cc + Zeff * alpha / r);
  }
  double c(double r) const final { return -alpha * (Zeff / r + en); }
  double d(double r) const final { return -a(r); }
};
 
double fitQuadratic(double x1, double x2, double x3, double y1, double y2,
                    double y3);
 
} // namespace DiracODE