tests/single_translations_vs_calc.c

   1 // TODO complex kr version
   2
   3 #include <stdbool.h>
   4 #include <complex.h>
   5 #include <string.h>
   6 #include <qpms/translations.h>
   7 #include <qpms/qpms_error.h>
   8 #include <gsl/gsl_rng.h>
   9 #include <gsl/gsl_randist.h>
  10 #include <qpms/indexing.h>
  11
  12 #define RTOL (1e-8)
  13 #define ATOL (1e-14)
  14
  15 int isclose_cmplx(complex double a, complex double b) {
  16   return cabs(a-b) <= ATOL + RTOL * .5 * (cabs(a) + cabs(b));
  17 }
  18
  19 static inline size_t ssq(size_t s) { return s * s; }
  20
  21 int test_AB_single_vs_array(const qpms_trans_calculator *c, qpms_bessel_t wavetype,
  22     cart3_t kd_cart)
  23 {
  24   int fails = 0;
  25   csph_t kd_sph = sph2csph(cart2sph(kd_cart));
  26
  27   complex double A[ssq(c->nelem)], B[ssq(c->nelem)];
  28   QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_AB_arrays(c, A, B, c->nelem, 1, kd_sph, false, wavetype));
  29
  30   for (qpms_y_t ydest = 0; ydest < c->nelem; ++ydest) {
  31     qpms_l_t ldest; qpms_m_t mdest; qpms_y2mn_p(ydest, &mdest, &ldest);
  32     for (qpms_y_t ysrc = 0; ysrc < c->nelem; ++ysrc) {
  33       qpms_l_t lsrc; qpms_m_t msrc; qpms_y2mn_p(ysrc, &msrc, &lsrc);
  34       complex double Asingle = qpms_trans_single_A(c->normalisation, mdest, ldest, msrc, lsrc, kd_sph, false, wavetype);
  35       complex double Aarr = A[ysrc + c->nelem * ydest];
  36       if (!isclose_cmplx(Asingle, Aarr)) {
  37         ++fails;
  38         fprintf(stderr, "l=%d,m=%+d <- l=%d,m=%+d: A_single=%.16g%+.16gj,\tA_arr=%.16g%+.16gj,\tdiff=%.g\t(norm=%x)\n",
  39             (int)ldest, (int)mdest, (int)lsrc, (int)msrc, creal(Asingle), cimag(Asingle),
  40             creal(Aarr), cimag(Aarr), cabs(Aarr-Asingle), (unsigned int)(c->normalisation));
  41       }
  42       complex double Bsingle = qpms_trans_single_B(c->normalisation, mdest, ldest, msrc, lsrc, kd_sph, false, wavetype);
  43       complex double Barr = B[ysrc + c->nelem * ydest];
  44       if (!isclose_cmplx(Bsingle, Barr)) {
  45         ++fails;
  46         fprintf(stderr, "l=%d,m=%+d <- l=%d,m=%+d: B_single=%.16g%+.16gj,\tB_arr=%.16g%+.16gj,\tdiff=%.g\t(norm=%x)\n",
  47             (int)ldest, (int)mdest, (int)lsrc, (int)msrc, creal(Bsingle), cimag(Bsingle),
  48             creal(Barr), cimag(Barr), cabs(Barr-Bsingle), (unsigned int)(c->normalisation));
  49       }
  50     }
  51   }
  52   return fails;
  53 }
  54
  55 int main() {
  56   gsl_rng *rng = gsl_rng_alloc(gsl_rng_ranlxs0);
  57   gsl_rng_set(rng, 666);
  58
  59         qpms_l_t lMax = 3;
  60         int npoints = 10;
  61         double sigma = 12;
  62
  63   cart3_t points[npoints];
  64   double relerrs[npoints];
  65   memset(points, 0, npoints * sizeof(cart3_t));
  66   points[0].x = points[1].y = points[2].z = sigma;
  67   for (unsigned i = 3; i < npoints; ++i) {
  68     cart3_t *w = points+i;
  69     w->x = gsl_ran_gaussian(rng, sigma);
  70     w->y = gsl_ran_gaussian(rng, sigma);
  71     w->z = gsl_ran_gaussian(rng, sigma);
  72   }
  73
  74   int fails = 0;
  75
  76   for(int use_csbit = 0; use_csbit <= 1; ++use_csbit) {
  77     for(int i = 0; i < 3; ++i){
  78       qpms_normalisation_t norm = ((qpms_normalisation_t[])
  79           { QPMS_NORMALISATION_NORM_SPHARM,
  80             QPMS_NORMALISATION_NORM_POWER,
  81             QPMS_NORMALISATION_NORM_NONE
  82             })[i]
  83         | (use_csbit ? QPMS_NORMALISATION_CSPHASE : 0);
  84       qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, norm);
  85       for(int J = 1; J <= 4; ++J)
  86         for(int p = 0; p < npoints; ++p)
  87           fails += test_AB_single_vs_array(c, J, points[p]);
  88       qpms_trans_calculator_free(c);
  89     }
  90   }
  91
  92   gsl_rng_free(rng);
  93
  94   return fails;
  95 }
  96