qpms/tests/translations_xutable_test.c

   1 // c99 -ggdb -I .. -DUSE_XU_ANTINORMALISATION -DUSE_BROKEN_SINGLETC -o translations_xutable translations_xutable_test.c ../translations.c ../gaunt.c -lgsl -lblas -lm
   2 #include "translations.h"
   3 #include <stdio.h>
   4 //#include <math.h>
   5 #include <complex.h>
   6
   7 typedef struct {
   8         qpms_normalisation_t norm;
   9         int m, n, mu, nu;
  10         sph_t kdlj;
  11         qpms_bessel_t J;
  12         complex double result_A, result_B;
  13 } testcase_single_trans_t;
  14
  15 testcase_single_trans_t testcases_xu[] = {
  16 #include "testcases_translations_Xu"
  17 };
  18
  19 int lMax=20;
  20
  21 int main() {
  22         qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_XU);
  23
  24         for(testcase_single_trans_t *tc = testcases_xu; tc->J != QPMS_BESSEL_UNDEF; tc++) {
  25                 if (!tc->n || !tc->nu || tc->n > lMax || tc->nu > lMax ) continue;
  26
  27                 printf("m=%d, n=%d, mu=%d, nu=%d,\n", tc->m,tc->n,tc->mu,tc->nu);
  28                 complex double A = qpms_trans_single_A(QPMS_NORMALISATION_XU,tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, false, tc->J);
  29                 complex double B = qpms_trans_single_B(QPMS_NORMALISATION_XU,tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, false, tc->J);
  30                 complex double A2 = qpms_trans_calculator_get_A(c, tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, false, tc->J);
  31                 complex double B2 = qpms_trans_calculator_get_B(c, tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, false, tc->J);
  32                 printf("A  = %.16e+%.16ej, relerr=%.16e, J=%d\n",
  33                                 creal(A), cimag(A),  (0 == cabs(tc->result_A - A)) ? 0 :
  34                                 cabs(tc->result_A - A)/((cabs(A) < cabs(tc->result_A)) ? cabs(A) : cabs(tc->result_A)),
  35                                 tc->J);
  36                 printf("A' = %.16e+%.16ej, relerr=%.16e, relerr2=%.3e\n",
  37                                 creal(A2), cimag(A2),  (0 == cabs(tc->result_A - A2)) ? 0 :
  38                                 cabs(tc->result_A - A2)/((cabs(A2) < cabs(tc->result_A)) ? cabs(A2) : cabs(tc->result_A)),
  39                                 (0 == cabs(A - A2)) ? 0 :
  40                                 cabs(A - A2)/((cabs(A2) < cabs(A)) ? cabs(A2) : cabs(A))
  41                       );
  42                 printf("B  = %.16e+%.16ej, relerr=%.16e, J=%d\n",
  43                                 creal(B), cimag(B),  (0 == cabs(tc->result_B - B)) ? 0 :
  44                                 cabs(tc->result_B - B)/((cabs(B) < cabs(tc->result_B)) ? cabs(B) : cabs(tc->result_B)),
  45                                 tc->J);
  46                 printf("B' = %.16e+%.16ej, relerr=%.16e, relerr2=%.3e\n",
  47                                 creal(B2), cimag(B2),  (0 == cabs(tc->result_B - B2)) ? 0 :
  48                                 cabs(tc->result_B - B2)/((cabs(B2) < cabs(tc->result_B)) ? cabs(B2) : cabs(tc->result_B)),
  49                                 (0 == cabs(B - B2)) ? 0 :
  50                                 cabs(B - B2)/((cabs(B2) < cabs(B)) ? cabs(B2) : cabs(B))
  51                       );
  52         }
  53         complex double A,B;
  54         // Test of zero R
  55         sph_t kdlj = {0, 1, 2};
  56         int m = -1, n = 1, mu = -1, nu = 1;
  57         qpms_trans_calculator_get_AB_p(c,&A,&B,m,n,mu,nu,kdlj,false,3);
  58         printf("A = %.6e+%.6ej, B = %.6e+%.6ej\n", creal(A),cimag(A),creal(B),cimag(B));
  59         qpms_trans_calculator_free(c);
  60 }
  61
  62