tests/test_vswf_translations_array.c

   1 //c99 -o test_vswf_translations_array -ggdb -I .. test_vswf_translations_array.c -lqpms -lgsl -lm -lblas
   2 #include <qpms/translations.h>
   3 #include <qpms/vswf.h>
   4 #include <stdio.h>
   5 #include <gsl/gsl_rng.h>
   6 #include <gsl/gsl_math.h>
   7 #include <gsl/gsl_randist.h>
   8 #include <assert.h>
   9 #include <qpms/vectors.h>
  10 #include <string.h>
  11 #include <qpms/indexing.h>
  12
  13 char *normstr(qpms_normalisation_t norm) {
  14         //int csphase = qpms_normalisation_t_csphase(norm);
  15         norm = norm & QPMS_NORMALISATION_NORM_BITS;
  16         switch (norm) {
  17                 case QPMS_NORMALISATION_NORM_NONE:
  18                         return "none";
  19                 case QPMS_NORMALISATION_NORM_SPHARM:
  20                         return "spharm";
  21                 case QPMS_NORMALISATION_NORM_POWER:
  22                         return "power";
  23                 default:
  24                         return "!!!undef!!!";
  25         }
  26 }
  27
  28 int test_sphwave_translation(const qpms_trans_calculator *c, qpms_bessel_t wavetype,
  29                 cart3_t o2minuso1, int npoints, cart3_t *o1points);
  30 //int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points);
  31 //int test_planewave_decomposition_silent(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points, double relerrthreshold, double *relerrs);
  32
  33 int main() {
  34         gsl_rng *rng =  gsl_rng_alloc(gsl_rng_ranlxs0);
  35         gsl_rng_set(rng, 666);
  36
  37         qpms_l_t lMax = 8;
  38         //qpms_l_t viewlMax = 2;
  39         int npoints = 10;
  40         double sigma = 4;
  41   //double shiftsigma = 2.;
  42
  43         cart3_t o2minuso1;
  44         o2minuso1.x = 1; //gsl_ran_gaussian(rng, shiftsigma);
  45         o2minuso1.y = 2; //gsl_ran_gaussian(rng, shiftsigma);
  46         o2minuso1.z = 5; //gsl_ran_gaussian(rng, shiftsigma);
  47
  48         cart3_t points[npoints];
  49         double relerrs[npoints];
  50         memset(points, 0, npoints * sizeof(cart3_t));
  51         points[0].x = points[1].y = points[2].z = sigma;
  52         double relerrthreshold = 1e-11;
  53         for (unsigned i = 3; i < npoints; ++i) {
  54                 cart3_t *w = points+i;
  55                 w->x = gsl_ran_gaussian(rng, sigma);
  56                 w->y = gsl_ran_gaussian(rng, sigma);
  57                 w->z = gsl_ran_gaussian(rng, sigma);
  58         }
  59
  60   for(int use_csbit = 0; use_csbit <= 1; ++use_csbit) {
  61     for(int i = 0; i < 3; ++i){
  62       qpms_normalisation_t norm = ((qpms_normalisation_t[])
  63           { QPMS_NORMALISATION_NORM_SPHARM,
  64             QPMS_NORMALISATION_NORM_POWER,
  65             QPMS_NORMALISATION_NORM_NONE
  66             })[i]
  67         | (use_csbit ? QPMS_NORMALISATION_CSPHASE : 0);
  68       qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, norm);
  69       for(int J = 1; J <= 4; ++J)
  70         test_sphwave_translation(c, J, o2minuso1, npoints, points);
  71       qpms_trans_calculator_free(c);
  72     }
  73
  74   }
  75         gsl_rng_free(rng);
  76 }
  77
  78 int test_sphwave_translation(const qpms_trans_calculator *c, qpms_bessel_t wavetype,
  79                 cart3_t sc, int npoints, cart3_t *points) {
  80         puts("==============================================================");
  81         printf("Test translation o2-o1 = %fx̂ + %fŷ + %fẑ", sc.x, sc.y, sc.z);
  82         sph_t ss = cart2sph(sc);
  83         printf(" = %fr̂ @ o1, θ = %f, φ = %f\n", ss.r, ss.theta, ss.phi);
  84         printf("lMax = %d, norm: %s, csphase = %d\n",
  85                         (int)c->lMax, normstr(c->normalisation), qpms_normalisation_t_csphase(c->normalisation));
  86         printf("wave type J = %d\n", wavetype);
  87
  88         qpms_l_t lMax = c->lMax;
  89         qpms_y_t nelem = c->nelem;
  90         csphvec_t N1[nelem], /* N2[nelem], */ M1[nelem] /*, M2[nelem]*/;
  91
  92         for (int i = 0; i < npoints; i++) {
  93                 printf("-------- Point %d --------\n", i);
  94                 cart3_t w1c = points[i];
  95                 cart3_t w2c = cart3_add(w1c, cart3_scale(-1, sc));
  96                 sph_t w1s = cart2sph(w1c);
  97                 sph_t w2s = cart2sph(w2c);
  98                 printf(" = %fx̂ + %fŷ + %fẑ @o1\n", w1c.x, w1c.y, w1c.z);
  99                 printf(" = %fx̂ + %fŷ + %fẑ @o2\n", w2c.x, w2c.y, w2c.z);
 100                 printf(" = %fr̂ @ o1, θ = %f, φ = %f\n", w1s.r, w1s.theta, w1s.phi);
 101                 printf(" = %fr̂ @ o2, θ = %f, φ = %f\n", w2s.r, w2s.theta, w2s.phi);
 102                 printf("Outside the sphere centered in o1 intersecting o2: %s; by %f\n", (w1s.r > ss.r) ? "true" : "false",
 103                                 w1s.r - ss.r);
 104                 if(QPMS_SUCCESS != qpms_vswf_fill(NULL, M1, N1, lMax, w1s, wavetype, c->normalisation))
 105                         abort(); // original wave set
 106
 107     complex double A_whole[nelem][nelem], B_whole[nelem][nelem];
 108     if (qpms_trans_calculator_get_AB_arrays(c,(complex double *) A_whole, (complex double *) B_whole,
 109           1, nelem, sph2csph(ss), (w1s.r > ss.r), wavetype)) abort();
 110
 111
 112                 for(qpms_y_t y1 = 0; y1 < nelem; ++y1) { //index of the wave originating in o1 that will be reconstructed in o2
 113                         qpms_m_t m1;
 114                         qpms_l_t l1;
 115                         qpms_y2mn_p(y1, &m1, &l1);
 116                         printf("*** wave l = %d, m = %d ***\n", l1, m1);
 117
 118                         //complex double A[nelem], B[nelem];
 119       /*
 120                         for(qpms_y_t y2 = 0; y2 < nelem; ++y2){
 121                                 qpms_m_t m2; qpms_l_t l2;
 122                                 qpms_y2mn_p(y2, &m2, &l2);
 123                                 if(qpms_trans_calculator_get_AB_p(c, &(A[y2]), &(B[y2]), m2, l2, m1, l1, ss, (w1s.r > ss.r) , wavetype))
 124                                         abort();
 125                         }
 126       */
 127       // array function instead above
 128       complex double *A = A_whole[y1]; // CHECKME right index ??
 129       complex double *B = B_whole[y1];
 130
 131
 132                         printf("M = ");
 133                         print_csphvec(M1[y1]);
 134                         printf(" @ o1\n  = ");
 135                         ccart3_t M1c = csphvec2ccart(M1[y1], w1s);
 136                         print_ccart3(M1c);
 137                         printf("\n  = ");
 138                         csphvec_t M1s2 = ccart2csphvec(M1c, w2s);
 139                         print_csphvec(M1s2);
 140                         printf(" @ o2\n");
 141                         csphvec_t M2s2 = qpms_eval_vswf(w2s, NULL, A, B, lMax, wavetype, c->normalisation);
 142                         printf("Mr= ");
 143                         print_csphvec(M2s2);
 144                         printf(" @ o2\n");
 145
 146                         printf("N = ");
 147                         print_csphvec(N1[y1]);
 148                         printf(" @ o1\n  = ");
 149                         ccart3_t N1c = csphvec2ccart(N1[y1], w1s);
 150                         print_ccart3(N1c);
 151                         printf("\n  = ");
 152                         csphvec_t N1s2 = ccart2csphvec(N1c, w2s);
 153                         print_csphvec(N1s2);
 154                         printf(" @o2\nNr= ");
 155                         csphvec_t N2s2 = qpms_eval_vswf(w2s, NULL, B, A, lMax, wavetype, c->normalisation);
 156                         print_csphvec(N2s2);
 157                         printf(" @o2\n");
 158                 }
 159         }
 160
 161         return 0; // FIXME something more meaningful here...
 162 }
 163
 164
 165
 166
 167
 168 #if 0
 169
 170 int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points){
 171         qpms_y_t nelem = qpms_lMax2nelem(lMax);
 172         complex double lc[nelem], mc[nelem], ec[nelem];
 173         if (QPMS_SUCCESS !=
 174                 qpms_planewave2vswf_fill_cart(k, E, lc, mc, ec, lMax, norm)) {
 175                 printf("Error\n");
 176                 return -1;
 177         }
 178         printf("==============================================================\n");
 179         printf("Test wave k = %fx̂ + %fŷ + %fẑ", k.x, k.y, k.z);
 180         sph_t k_sph = cart2sph(k);
 181         printf(" = %fr̂ @ θ = %f, φ = %f\n", k_sph.r, k_sph.theta, k_sph.phi);
 182         printf("        E_0 = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ",
 183                 creal(E.x),cimag(E.x),
 184                 creal(E.y),cimag(E.y),
 185                 creal(E.z),cimag(E.z));
 186         csphvec_t E_s = ccart2csphvec(E, k_sph);
 187         printf(" = (%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂  @ k\n",
 188                 creal(E_s.rc), cimag(E_s.rc), creal(E_s.thetac), cimag(E_s.thetac),
 189                 creal(E_s.phic), cimag(E_s.phic));
 190         printf("        lMax = %d, norm: %s, csphase = %d\n",
 191                         (int)lMax, normstr(norm), qpms_normalisation_t_csphase(norm));
 192         printf("a_L: ");
 193         for(qpms_y_t y = 0; y < nelem; ++y) printf("%g+%gj ", creal(lc[y]), cimag(lc[y]));
 194         printf("\na_M: ");
 195         for(qpms_y_t y = 0; y < nelem; ++y) printf("%g+%gj ", creal(mc[y]), cimag(mc[y]));
 196         printf("\na_N: ");
 197         for(qpms_y_t y = 0; y < nelem; ++y) printf("%g+%gj ", creal(ec[y]), cimag(ec[y]));
 198         printf("\n");
 199         for (int i = 0; i < npoints; i++) {
 200                 cart3_t w = points[i];
 201                 sph_t w_sph = cart2sph(w);
 202                 printf("Point %d: x = %f, y = %f, z = %f,\n", i, w.x, w.y, w.z);
 203                 printf("        |r| = %f, θ = %f, φ = %f:\n", w_sph.r, w_sph.theta, w_sph.phi);
 204                 double phase = cart3_dot(k,w);
 205                 printf("  k.r = %f\n", phase);
 206                 complex double phfac = cexp(phase * I);
 207                 ccart3_t Ew = ccart3_scale(phfac, E);
 208                 printf("  pw  E(r) = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ",
 209                         creal(Ew.x),cimag(Ew.x),
 210                         creal(Ew.y),cimag(Ew.y),
 211                         creal(Ew.z),cimag(Ew.z));
 212                 csphvec_t Ew_s = ccart2csphvec(Ew, w_sph);
 213                 printf(" = (%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂  @ r\n",
 214                         creal(Ew_s.rc), cimag(Ew_s.rc),
 215                         creal(Ew_s.thetac), cimag(Ew_s.thetac),
 216                         creal(Ew_s.phic), cimag(Ew_s.phic));
 217                 w_sph.r *= k_sph.r; /// NEVER FORGET THIS!!!
 218                 csphvec_t Ew_s_recomp = qpms_eval_vswf(w_sph,
 219                         lc, mc, ec, lMax, QPMS_BESSEL_REGULAR, norm);
 220                 ccart3_t Ew_recomp = csphvec2ccart(Ew_s_recomp, w_sph);
 221                 printf(" rec E(r) = (%f+%fj)x̂ + (%f+%fj)ŷ + (%f+%fj)ẑ",
 222                         creal(Ew_recomp.x),cimag(Ew_recomp.x),
 223                         creal(Ew_recomp.y),cimag(Ew_recomp.y),
 224                         creal(Ew_recomp.z),cimag(Ew_recomp.z));
 225                 printf(" = (%f+%fj)r̂ + (%f+%fj)θ̂ + (%f+%fj)φ̂  @ r\n",
 226                         creal(Ew_s_recomp.rc), cimag(Ew_s_recomp.rc),
 227                         creal(Ew_s_recomp.thetac), cimag(Ew_s_recomp.thetac),
 228                         creal(Ew_s_recomp.phic), cimag(Ew_s_recomp.phic));
 229                 double relerrfac = 2/(cabs(Ew_s_recomp.rc) + cabs(Ew_s.rc)
 230                                 +cabs(Ew_s_recomp.thetac) + cabs(Ew_s.thetac)
 231                                 +cabs(Ew_s_recomp.phic) + cabs(Ew_s.phic));
 232                 printf(" rel. err. magnitude: %g @ r̂, %g @ θ̂, %g @ φ̂\n",
 233                         cabs(Ew_s_recomp.rc - Ew_s.rc) * relerrfac,
 234                         cabs(Ew_s_recomp.thetac - Ew_s.thetac) * relerrfac,
 235                         cabs(Ew_s_recomp.phic - Ew_s.phic) * relerrfac
 236                       );
 237         }
 238         return 0;
 239 }
 240
 241 int test_planewave_decomposition_silent(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points, double relerrthreshold, double *relerrs) {
 242         qpms_y_t nelem = qpms_lMax2nelem(lMax);
 243         int failcount = 0;
 244         complex double lc[nelem], mc[nelem], ec[nelem];
 245         if (QPMS_SUCCESS !=
 246                 qpms_planewave2vswf_fill_cart(k, E, lc, mc, ec, lMax, norm)) {
 247                 printf("Error\n");
 248                 return -1;
 249         }
 250         sph_t k_sph = cart2sph(k);
 251         csphvec_t E_s = ccart2csphvec(E, k_sph);
 252         for (int i = 0; i < npoints; i++) {
 253                 cart3_t w = points[i];
 254                 sph_t w_sph = cart2sph(w);
 255                 w_sph.r *= k_sph.r;
 256                 double phase = cart3_dot(k,w);
 257                 complex double phfac = cexp(phase * I);
 258                 ccart3_t Ew = ccart3_scale(phfac, E);
 259                 csphvec_t Ew_s = ccart2csphvec(Ew, w_sph);
 260                 csphvec_t Ew_s_recomp = qpms_eval_vswf(w_sph,
 261                         lc, mc, ec, lMax, QPMS_BESSEL_REGULAR, norm);
 262                 ccart3_t Ew_recomp = csphvec2ccart(Ew_s_recomp, w_sph);
 263                 double relerrfac = 2/(cabs(Ew_s_recomp.rc) + cabs(Ew_s.rc)
 264                                 +cabs(Ew_s_recomp.thetac) + cabs(Ew_s.thetac)
 265                                 +cabs(Ew_s_recomp.phic) + cabs(Ew_s.phic));
 266
 267                 double relerr = (cabs(Ew_s_recomp.rc - Ew_s.rc)
 268                         +       cabs(Ew_s_recomp.thetac - Ew_s.thetac)
 269                         +       cabs(Ew_s_recomp.phic - Ew_s.phic)
 270                         ) * relerrfac;
 271                 if(relerrs) relerrs[i] = relerr;
 272                 if(relerr > relerrthreshold) ++failcount;
 273         }
 274         return failcount;
 275 }
 276 #endif