Fix single translation coefficients.

Former-commit-id: 1c045d994bb46ed16350d294bacc9203854cd2ea

CMakeLists.txt

@@ -30,6 +30,8 @@ cmake_add_fortran_subdirectory (amos
 	NO_EXTERNAL_INSTALL)
 add_subdirectory (qpms)
 
-add_subdirectory (tests)
+
+enable_testing()
+add_subdirectory (tests EXCLUDE_FROM_ALL)
 
 #add_subdirectory (apps/transop-ewald)

qpms/translations.c

@@ -183,10 +183,8 @@ complex double qpms_trans_single_A(qpms_normalisation_t norm,
   double a1q0 = a1q[0];
   if (err) abort();
 
-  int csphase = qpms_normalisation_t_csphase(norm);
-
   double leg[gsl_sf_legendre_array_n(n+nu)];
-  if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,csphase,leg)) abort();
+  if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort();
   complex double bes[n+nu+1];
   if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort();
   complex double sum = 0;
@@ -345,9 +343,6 @@ complex double qpms_trans_single_B(qpms_normalisation_t norm,
     int m, int n, int mu, int nu, sph_t kdlj,
     bool r_ge_d, qpms_bessel_t J) { 
   TROPS_ONLY_EIMF_IMPLEMENTED(norm);
-#ifndef USE_BROKEN_SINGLETC
-  assert(0); // FIXME probably gives wrong values, do not use.
-#endif
   if(r_ge_d) J = QPMS_BESSEL_REGULAR;
   double costheta = cos(kdlj.theta);
 
@@ -368,10 +363,8 @@ complex double qpms_trans_single_B(qpms_normalisation_t norm,
   gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err); if (err) abort();
   a3q0 = a3q[0];
   
-  int csphase = qpms_normalisation_t_csphase(norm);
-
   double leg[gsl_sf_legendre_array_n(n+nu+1)];
-  if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,csphase,leg)) abort();
+  if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort();
   complex double bes[n+nu+2];
   if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort();
 

tests/CMakeLists.txt

@@ -6,3 +6,11 @@ target_include_directories ( test_vswf_translations_array PRIVATE .. )
 add_executable( test_vswf_translations test_vswf_translations.c )
 target_link_libraries( test_vswf_translations qpms gsl blas lapacke amos m )
 target_include_directories ( test_vswf_translations PRIVATE .. )
+
+add_executable( test_single_translations_vs_calc single_translations_vs_calc.c )
+target_link_libraries( test_single_translations_vs_calc qpms gsl lapacke amos m )
+target_include_directories( test_single_translations_vs_calc PRIVATE .. )
+
+add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array )
+
+add_test( NAME single_vs_array_translation_coeffs COMMAND test_single_translations_vs_calc )

tests/single_translations_vs_calc.c

@@ -0,0 +1,96 @@
+// TODO complex kr version
+
+#include <stdbool.h>
+#include <complex.h>
+#include <string.h>
+#include <qpms/translations.h>
+#include <qpms/qpms_error.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <qpms/indexing.h>
+
+#define RTOL (1e-8)
+#define ATOL (1e-14)
+
+int isclose_cmplx(complex double a, complex double b) {
+  return cabs(a-b) <= ATOL + RTOL * .5 * (cabs(a) + cabs(b));
+}
+
+static inline size_t ssq(size_t s) { return s * s; }
+
+int test_AB_single_vs_array(const qpms_trans_calculator *c, qpms_bessel_t wavetype,
+    cart3_t kd_cart) 
+{
+  int fails = 0;
+  sph_t kd_sph = cart2sph(kd_cart);
+
+  complex double A[ssq(c->nelem)], B[ssq(c->nelem)];
+  QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_AB_arrays(c, A, B, c->nelem, 1, kd_sph, false, wavetype));
+  
+  for (qpms_y_t ydest = 0; ydest < c->nelem; ++ydest) {
+    qpms_l_t ldest; qpms_m_t mdest; qpms_y2mn_p(ydest, &mdest, &ldest);
+    for (qpms_y_t ysrc = 0; ysrc < c->nelem; ++ysrc) {
+      qpms_l_t lsrc; qpms_m_t msrc; qpms_y2mn_p(ysrc, &msrc, &lsrc);
+      complex double Asingle = qpms_trans_single_A(c->normalisation, mdest, ldest, msrc, lsrc, kd_sph, false, wavetype);
+      complex double Aarr = A[ysrc + c->nelem * ydest];
+      if (!isclose_cmplx(Asingle, Aarr)) {
+        ++fails;
+        fprintf(stderr, "l=%d,m=%+d <- l=%d,m=%+d: A_single=%.16g%+.16gj,\tA_arr=%.16g%+.16gj,\tdiff=%.g\t(norm=%x)\n",
+            (int)ldest, (int)mdest, (int)lsrc, (int)msrc, creal(Asingle), cimag(Asingle), 
+            creal(Aarr), cimag(Aarr), cabs(Aarr-Asingle), (unsigned int)(c->normalisation));
+      }
+      complex double Bsingle = qpms_trans_single_B(c->normalisation, mdest, ldest, msrc, lsrc, kd_sph, false, wavetype);
+      complex double Barr = B[ysrc + c->nelem * ydest];
+      if (!isclose_cmplx(Bsingle, Barr)) {
+        ++fails;
+        fprintf(stderr, "l=%d,m=%+d <- l=%d,m=%+d: B_single=%.16g%+.16gj,\tB_arr=%.16g%+.16gj,\tdiff=%.g\t(norm=%x)\n",
+            (int)ldest, (int)mdest, (int)lsrc, (int)msrc, creal(Bsingle), cimag(Bsingle), 
+            creal(Barr), cimag(Barr), cabs(Barr-Bsingle), (unsigned int)(c->normalisation));
+      }
+    }
+  }
+  return fails;
+}  
+
+int main() {
+  gsl_rng *rng = gsl_rng_alloc(gsl_rng_ranlxs0);
+  gsl_rng_set(rng, 666);
+
+	qpms_l_t lMax = 3;
+	int npoints = 10;
+	double sigma = 12;
+
+  cart3_t points[npoints];
+  double relerrs[npoints];
+  memset(points, 0, npoints * sizeof(cart3_t));
+  points[0].x = points[1].y = points[2].z = sigma;
+  for (unsigned i = 3; i < npoints; ++i) {
+    cart3_t *w = points+i;
+    w->x = gsl_ran_gaussian(rng, sigma);
+    w->y = gsl_ran_gaussian(rng, sigma);
+    w->z = gsl_ran_gaussian(rng, sigma);
+  }
+
+  int fails = 0;
+
+  for(int use_csbit = 0; use_csbit <= 1; ++use_csbit) {
+    for(int i = 0; i < 3; ++i){
+      qpms_normalisation_t norm = ((qpms_normalisation_t[])
+          { QPMS_NORMALISATION_NORM_SPHARM,
+            QPMS_NORMALISATION_NORM_POWER,
+            QPMS_NORMALISATION_NORM_NONE
+            })[i]
+        | (use_csbit ? QPMS_NORMALISATION_CSPHASE : 0);
+      qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, norm);
+      for(int J = 1; J <= 4; ++J)
+        for(int p = 0; p < npoints; ++p)
+          fails += test_AB_single_vs_array(c, J, points[p]);
+      qpms_trans_calculator_free(c);
+    }
+  }  
+
+  gsl_rng_free(rng);
+  
+  return fails;
+}
+