diff --git a/CMakeLists.txt b/CMakeLists.txt index e3e7eeb..7f80c81 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/qpms/translations.c b/qpms/translations.c index 87cdd9a..de3cdbb 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -148,10 +148,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; @@ -207,9 +205,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); @@ -230,10 +225,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(); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c9210d5..61ae546 100644 --- a/tests/CMakeLists.txt +++ b/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 ) diff --git a/tests/single_translations_vs_calc.c b/tests/single_translations_vs_calc.c new file mode 100644 index 0000000..6545fae --- /dev/null +++ b/tests/single_translations_vs_calc.c @@ -0,0 +1,96 @@ +// TODO complex kr version + +#include +#include +#include +#include +#include +#include +#include +#include + +#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; +} +