qpms/tests/vswf_translation_test_rays.c

   1 // c99 -o ../../tests/raytests/n_nu0 -ggdb -I .. vswf_translation_test_rays.c ../translations.c ../vswf.c ../gaunt.c ../legendre.c -lgsl -lm -lblas
   2
   3 #include "translations.h"
   4 #include "vswf.h"
   5 #include <stdio.h>
   6 #include <stdlib.h>
   7 #include <gsl/gsl_rng.h>
   8 #include <gsl/gsl_math.h>
   9 #include <gsl/gsl_randist.h>
  10 #include <assert.h>
  11 #include "vectors.h"
  12 #include "string.h"
  13 #include "indexing.h"
  14
  15 #define MIN(x,y) (((x)<(y))?(x):(y))
  16
  17 char *normstr(qpms_normalisation_t norm) {
  18   //int csphase = qpms_normalisation_t_csphase(norm);
  19   norm = qpms_normalisation_t_normonly(norm);
  20   switch (norm) {
  21     case QPMS_NORMALISATION_NONE:
  22       return "none";
  23     case QPMS_NORMALISATION_SPHARM:
  24       return "spharm";
  25     case QPMS_NORMALISATION_POWER:
  26       return "power";
  27     default:
  28       return "!!!undef!!!";
  29   }
  30 }
  31
  32 int test_sphwave_translation(const qpms_trans_calculator *c, qpms_bessel_t wavetype,
  33     cart3_t o2minuso1, int npoints, cart3_t *o1points);
  34 //int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points);
  35 //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);
  36
  37 int main(int argc, char **argv) {
  38   gsl_rng *rng =  gsl_rng_alloc(gsl_rng_ranlxs0);
  39   gsl_rng_set(rng, 666);
  40
  41   // parametry: prefix, n, m, norma, lmax, J
  42   // defaults
  43   qpms_m_t m1 = 0;
  44   qpms_l_t l1 = 1;
  45   qpms_normalisation_t norm = QPMS_NORMALISATION_NONE_CS;
  46   qpms_l_t lMax = 10;
  47   qpms_bessel_t J = 1;
  48
  49   switch(MIN(argc,7)) {
  50     case 7: J = atoi(argv[6]);
  51     case 6: lMax = atoi(argv[5]);
  52     case 5: norm = atoi(argv[4]);
  53     case 4: m1 = atoi(argv[3]);
  54     case 3: l1 = atoi(argv[2]);
  55     case 2: break; // first argument is the filename prefix
  56     case 1: fputs("At least one argument needed (filename prefix)\n", stderr);
  57   }
  58   qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, norm);
  59   qpms_y_t nelem = c->nelem;
  60   qpms_y_t y1 = qpms_mn2y(m1, l1);
  61
  62   //qpms_l_t viewlMax = 2;
  63   int npoints = 3;
  64   double sigma = 4.;
  65   //double shiftsigma = 0.01;
  66
  67   cart3_t o2minuso1_max;
  68   // o2minuso1.x = gsl_ran_gaussian(rng, shiftsigma);
  69   // o2minuso1.y = gsl_ran_gaussian(rng, shiftsigma);
  70   // o2minuso1.z = gsl_ran_gaussian(rng, shiftsigma);
  71   o2minuso1_max.x = 1;
  72   o2minuso1_max.y = 2;
  73   o2minuso1_max.z = 5;
  74   int nraysteps = 512;
  75   cart3_t o2mo1_step = cart3_scale(1./nraysteps, o2minuso1_max);
  76
  77   // measurement points
  78   cart3_t points[npoints];
  79   double relerrs[npoints];
  80   memset(points, 0, npoints * sizeof(cart3_t));
  81   points[0].x = points[1].y = points[2].z = sigma;
  82   points[3].x = points[3].y = 1./M_SQRT2;
  83   double relerrthreshold = 1e-11;
  84   for (unsigned i = 3; i < npoints; ++i) {
  85     cart3_t *w = points+i;
  86     w->x = gsl_ran_gaussian(rng, sigma);
  87     w->y = gsl_ran_gaussian(rng, sigma);
  88     w->z = gsl_ran_gaussian(rng, sigma);
  89   }
  90   // each point will have its file
  91   FILE *pfile[npoints];
  92   {
  93     char fname[strlen(argv[1])+20];
  94     for (int p = 0; p < npoints; ++p) {
  95       sprintf(fname, "%s.%d", argv[1], p);
  96       pfile[p] = fopen(fname, "w");
  97       fputs("##\tnu\tmu\t%d\tnorm\tlMax\tnelem\tJ\n", pfile[p]);
  98       fprintf(pfile[p], "#\t%d\t%d\t%d\t%d\t%d\t%d\t%d\n",
  99           (int)l1, (int)m1, (int)y1, (int)norm, (int)lMax, (int)nelem, (int)J);
 100       fprintf(pfile[p], "#\tThis point cart:\t%g\t%g\t%g\n",
 101           points[p].x, points[p].y, points[p].z);
 102       sph_t ps = cart2sph(points[p]);
 103       fprintf(pfile[p], "#\tThis point sph:\t%g\t%g\t%g\n",
 104           ps.r, ps.theta, ps.phi);
 105       csphvec_t M1[nelem], N1[nelem]; // mohl bych to počítat jednotlivě...
 106       if(QPMS_SUCCESS != qpms_vswf_fill(NULL, M1, N1, lMax, ps, J, norm)) abort();
 107       fprintf(pfile[p],"##\tM(%d,%d) at this point:\n", l1, m1);
 108       fputs("##\tRe(rc)\tIm(rc)\tRe(tc)\tIm(tc)\tRe(fc)\tIm(fc)\n",pfile[p]);
 109       fprintf(pfile[p],"#\t%g\t%g\t%g\t%g\t%g\t%g\n",
 110           creal(M1[y1].rc),cimag(M1[y1].rc),
 111           creal(M1[y1].thetac),cimag(M1[y1].thetac),
 112           creal(M1[y1].phic),cimag(M1[y1].phic)
 113       );
 114       fprintf(pfile[p],"##\tN(%d,%d) at this point:\n", l1, m1);
 115       fputs("##\tRe(rc)\tIm(rc)\tRe(tc)\tIm(tc)\tRe(fc)\tIm(fc)\n",pfile[p]);
 116       fprintf(pfile[p],"#\t%g\t%g\t%g\t%g\t%g\t%g\n",
 117           creal(N1[y1].rc),cimag(N1[y1].rc),
 118           creal(N1[y1].thetac),cimag(N1[y1].thetac),
 119           creal(N1[y1].phic),cimag(N1[y1].phic)
 120       );
 121       // TODO print column headers here
 122       fputs("## TODO print column headers here\n", pfile[p]);
 123       fputs("#step\t"
 124         "(o2-o1).x\t(o2-o1).y\t(o2-o1).z\t(o2-o1).r\t(o2-o1).θ\t(o2-o1).φ\t"
 125         "(x-o2).x\t(x-o2).y\t(x-o2).z\t(x-o2).r\t(x-o2).θ\t(x-o2).φ\t",
 126         pfile[p]
 127       );
 128
 129       fprintf(pfile[p], //original and reconstructed waves at new origin
 130           "R(M1@2(%d,%d).r)\tI(M1@2(%d,%d).r)\t"
 131           "R(M1@2(%d,%d).θ)\tI(M1@2(%d,%d).θ)\t"
 132           "R(M1@2(%d,%d).φ)\tI(M1@2(%d,%d).φ)\t"
 133           "R(N1@2(%d,%d).r)\tI(N1@2(%d,%d).r)\t"
 134           "R(N1@2(%d,%d).θ)\tI(N1@2(%d,%d).θ)\t"
 135           "R(N1@2(%d,%d).φ)\tI(N1@2(%d,%d).φ)\t"
 136           "R(Mr1@2(%d,%d).r)\tI(Mr1@2(%d,%d).r)\t"
 137           "R(Mr1@2(%d,%d).θ)\tI(Mr1@2(%d,%d).θ)\t"
 138           "R(Mr1@2(%d,%d).φ)\tI(Mr1@2(%d,%d).φ)\t"
 139           "R(Nr1@2(%d,%d).r)\tI(Nr1@2(%d,%d).r)\t"
 140           "R(Nr1@2(%d,%d).θ)\tI(Nr1@2(%d,%d).θ)\t"
 141           "R(Nr1@2(%d,%d).φ)\tI(Nr1@2(%d,%d).φ)\t",
 142           l1,m1,l1,m1, l1,m1,l1,m1, l1,m1,l1,m1,
 143           l1,m1,l1,m1, l1,m1,l1,m1, l1,m1,l1,m1,
 144           l1,m1,l1,m1, l1,m1,l1,m1, l1,m1,l1,m1,
 145           l1,m1,l1,m1, l1,m1,l1,m1, l1,m1,l1,m1
 146       );
 147
 148       for(qpms_y_t y = 0; y < nelem; ++y) {
 149         qpms_l_t l2; qpms_m_t m2;
 150         qpms_y2mn_p(y, &m2, &l2);
 151         fprintf(pfile[p], "R(A(%d,%d))\tI(A(%d,%d))\t"
 152           "R(B(%d,%d))\tI(B(%d,%d))\t"
 153           "R(M2(%d,%d).r)\tI(M2(%d,%d).r)\t"
 154           "R(M2(%d,%d).θ)\tI(M2(%d,%d).θ)\t"
 155           "R(M2(%d,%d).φ)\tI(M2(%d,%d).φ)\t"
 156           "R(N2(%d,%d).r)\tI(N2(%d,%d).r)\t"
 157           "R(N2(%d,%d).θ)\tI(N2(%d,%d).θ)\t"
 158           "R(N2(%d,%d).φ)\tI(N2(%d,%d).φ)\t",
 159           l2,m2,l2,m2, l2,m2,l2,m2, l2,m2,l2,m2, l2,m2,l2,m2, l2,m2,l2,m2,
 160           l2,m2,l2,m2, l2,m2,l2,m2, l2,m2,l2,m2
 161         );
 162       }
 163       fputc('\n', pfile[p]);
 164     }
 165   }
 166
 167   for (int r = 0; r <= nraysteps; ++r) {
 168     cart3_t sc = cart3_scale(r, o2mo1_step);
 169     sph_t ss = cart2sph(sc);
 170     csphvec_t N1[nelem], M1[nelem];
 171     for(int p = 0; p < npoints; p++){
 172       FILE *f = pfile[p];
 173       cart3_t w1c = points[p];
 174       cart3_t w2c = cart3_add(w1c, cart3_scale(-1.,sc));
 175       sph_t w1s = cart2sph(w1c);
 176       sph_t w2s = cart2sph(w2c);
 177       fprintf(f, "%d\t" "%g\t%g\t%g\t%g\t%g\t%g\t" "%g\t%g\t%g\t%g\t%g\t%g\t",
 178           r, sc.x, sc.y, sc.z, ss.r, ss.theta, ss.phi,
 179           w2c.x, w2c.y, w2c.z, w2s.r, w2s.theta, w2s.phi);
 180       complex double A[nelem], B[nelem]; // translation ceofficients
 181       if(QPMS_SUCCESS != qpms_vswf_fill(NULL, M1, N1, lMax, w1s, J, norm)) abort();
 182       csphvec_t M2at2[nelem], N2at2[nelem], M1at2, N1at2, M1Rat2, N1Rat2;
 183       if(QPMS_SUCCESS != qpms_vswf_fill(NULL, M2at2, N2at2, lMax, w2s, J, norm)) abort();
 184       for(qpms_y_t y2 = 0; y2 < nelem; ++y2) {
 185         qpms_l_t l2; qpms_m_t m2;
 186         qpms_y2mn_p(y2, &m2, &l2);
 187         if (QPMS_SUCCESS != qpms_trans_calculator_get_AB_p(c, &(A[y2]), &(B[y2]),
 188 #ifdef REVERSE
 189               m1, l1, m2, l2,
 190 #else
 191               m2, l2, m1, l1, // !!! FIXME mám správné pořadí??? !!!
 192 #endif
 193             ss, true /* FIXME Pro J != 1 */,  J)) abort();
 194       }
 195
 196       M1at2 = ccart2csphvec(csphvec2ccart(M1[y1], w1s), w2s);
 197       N1at2 = ccart2csphvec(csphvec2ccart(N1[y1], w1s), w2s);
 198       fprintf(f, "%g\t%g\t%g\t%g\t%g\t%g\t" "%g\t%g\t%g\t%g\t%g\t%g\t",
 199         creal(M1at2.rc), cimag(M1at2.rc),
 200         creal(M1at2.thetac), cimag(M1at2.thetac),
 201         creal(M1at2.phic), cimag(M1at2.phic),
 202         creal(N1at2.rc), cimag(N1at2.rc),
 203         creal(N1at2.thetac), cimag(N1at2.thetac),
 204         creal(N1at2.phic), cimag(N1at2.phic)
 205       );
 206       M1Rat2 = qpms_eval_vswf(w2s, NULL, A, B, lMax, J, norm);
 207       N1Rat2 = qpms_eval_vswf(w2s, NULL, B, A, lMax, J, norm);
 208       fprintf(f, "%g\t%g\t%g\t%g\t%g\t%g\t" "%g\t%g\t%g\t%g\t%g\t%g\t",
 209         creal(M1Rat2.rc), cimag(M1Rat2.rc),
 210         creal(M1Rat2.thetac), cimag(M1Rat2.thetac),
 211         creal(M1Rat2.phic), cimag(M1Rat2.phic),
 212         creal(N1Rat2.rc), cimag(N1Rat2.rc),
 213         creal(N1Rat2.thetac), cimag(N1Rat2.thetac),
 214         creal(N1Rat2.phic), cimag(N1Rat2.phic)
 215       );
 216
 217       for(qpms_y_t y2 = 0; y2 < nelem; ++y2){
 218         fprintf(f, "%g\t%g\t" "%g\t%g\t" "%g\t%g\t%g\t%g\t%g\t%g\t"
 219           "%g\t%g\t%g\t%g\t%g\t%g\t",
 220           creal(A[y2]), cimag(A[y2]),
 221           creal(B[y2]), cimag(B[y2]),
 222           creal(M2at2[y2].rc), cimag(M2at2[y2].rc),
 223           creal(M2at2[y2].thetac), cimag(M2at2[y2].thetac),
 224           creal(M2at2[y2].phic), cimag(M2at2[y2].phic),
 225           creal(N2at2[y2].rc), cimag(N2at2[y2].rc),
 226           creal(N2at2[y2].thetac), cimag(N2at2[y2].thetac),
 227           creal(N2at2[y2].phic), cimag(N2at2[y2].phic)
 228         );
 229       }
 230       fputc('\n', f);
 231     }
 232   }
 233
 234   for (int p = 0; p < npoints; ++p)
 235     fclose(pfile[p]);
 236   gsl_rng_free(rng);
 237   return 0;
 238 }
 239