#include <catch2/catch.hpp>
#include <qpms_error.h>
#include <vswf.h>
#include "catch_aux.h"
using namespace qpmstest;

constexpr qpms_l_t lMax = 6;

static std::vector<double> vswfvals_dvector(double r, double theta, double phi, 
    qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm) 
{
  std::vector<double> vals(2 /*el and mg for now */ * (lMax * (lMax+2)) 
      * 3 /* 3D vector components */ * 2 /*for real and imag*/);
  qpms_vswf_fill( 
      nullptr, (csphvec_t *) &(vals[0]), (csphvec_t *) &(vals[(lMax*(lMax+2)) * 3 * 2]),
      lMax, sph_t{r, theta, phi},
      btyp, norm);
  return vals;
}

static std::vector<double> vswfvals_dvector_alternative(double r, double theta, double phi, 
    qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm) 
{
  std::vector<double> vals(2 /*el and mg for now */ * (lMax * (lMax+2)) 
      * 3 /* 3D vector components */ * 2 /*for real and imag*/);
  qpms_vswf_fill_alternative( 
      nullptr, (csphvec_t *) &(vals[0]), (csphvec_t *) &(vals[(lMax*(lMax+2)) * 3 * 2]),
      lMax, sph_t{r, theta, phi},
      btyp, norm);
  return vals;
}

TEST_CASE("VSWF alternative implementation") {
  using Catch::Matchers::Approx;
  auto norms = GENERATE(QPMS_NORMALISATION_DEFAULT, QPMS_NORMALISATION_NORM_POWER,
      QPMS_NORMALISATION_NORM_SPHARM, QPMS_NORMALISATION_CONVENTION_SCUFF);
  auto thetas = GENERATE(0., M_PI_2, M_PI_2/2, M_PI, 3*M_PI_2, 0.1, -0.2, 4., 12.666);
  auto phis = GENERATE(0., M_PI_2, M_PI_2/2, M_PI, 3*M_PI_2, 0.1, -0.2, 4., 12.666);
  auto rs_positive = GENERATE(0.0001, 0.001, 0.01, 0.1, 1., 10., 100.);
  CHECK_THAT( 
        vswfvals_dvector(rs_positive, thetas, phis, lMax, QPMS_HANKEL_PLUS, norms),
        Approx(vswfvals_dvector_alternative(rs_positive, thetas, phis, lMax, QPMS_HANKEL_PLUS, norms))
      ) ;
}