qpms/tests/vswf_errs.c

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