qpms/tests/ray_planewaves.c

   1 // c99 -o ../../tests/pwraytests/raypwtest1 -I .. ray_planewaves.c -lgsl -lm ../legendre.c ../vswf.c ../bessel.c -lblas
   2 #include <gsl/gsl_rng.h>
   3 #include <gsl/gsl_math.h>
   4 #include <gsl/gsl_randist.h>
   5 #include <assert.h>
   6 #include "vswf.h"
   7 #include "vectors.h"
   8 #include "string.h"
   9 #include "indexing.h"
  10
  11 char *normstr(qpms_normalisation_t norm) {
  12         //int csphase = qpms_normalisation_t_csphase(norm);
  13         norm = qpms_normalisation_t_normonly(norm);
  14         switch (norm) {
  15                 case QPMS_NORMALISATION_NONE:
  16                         return "none";
  17                 case QPMS_NORMALISATION_SPHARM:
  18                         return "spharm";
  19                 case QPMS_NORMALISATION_POWER:
  20                         return "power";
  21                 default:
  22                         return "!!!undef!!!";
  23         }
  24 }
  25 #define MIN(x,y) (((x)<(y))?(x):(y))
  26
  27 int main(int argc, char **argv) {
  28         gsl_rng *rng =  gsl_rng_alloc(gsl_rng_ranlxs0);
  29         gsl_rng_set(rng, 666);
  30
  31
  32   int nraysteps = 512;
  33
  34
  35   // parametry: prefix, kx, ky, kz, E0x, E0y, E0z, norma, lmax
  36   // defaults
  37   double kx = 1, ky = 0, kz = 0, E0x = 0, E0y = 0, E0z = 1;
  38         qpms_l_t lMax = 10;
  39   qpms_normalisation_t norm = QPMS_NORMALISATION_NONE_CS;
  40
  41   switch(MIN(argc, 10)) {
  42     case 10: lMax = atoi(argv[9]);
  43     case 9: norm = atoi(argv[8]);
  44     case 8: E0z = strtod(argv[7], NULL);
  45     case 7: E0y = strtod(argv[6], NULL);
  46     case 6: E0x = strtod(argv[5], NULL);
  47     case 5: kz = strtod(argv[4], NULL);
  48     case 4: ky = strtod(argv[3], NULL);
  49     case 3: kx = strtod(argv[2], NULL);
  50     case 2: break;
  51     case 1: fputs("At least one argument needed (filename prefix)\n", stderr);
  52   }
  53
  54   qpms_y_t nelem = qpms_lMax2nelem(lMax);
  55         int nrays = 10;
  56         double sigma = 4;
  57         cart3_t rays[nrays];
  58         double relerrs[nrays];
  59         memset(rays, 0, nrays * sizeof(cart3_t));
  60         rays[0].x = rays[1].y = rays[2].z = sigma;
  61         double relerrthreshold = 1e-11;
  62         for (unsigned i = 3; i < nrays; ++i) {
  63                 cart3_t *w = rays+i;
  64                 w->x = gsl_ran_gaussian(rng, sigma);
  65                 w->y = gsl_ran_gaussian(rng, sigma);
  66                 w->z = gsl_ran_gaussian(rng, sigma);
  67         }
  68
  69
  70
  71   cart3_t k = {kx, ky, kz};
  72   k = cart3_scale(1./cart3norm(k),k); // automatically normalise
  73
  74         ccart3_t E0 = {E0x, E0y, E0z};
  75   sph_t ks = cart2sph(k);
  76   csphvec_t E0s = ccart2csphvec(E0, ks);
  77   complex double Lc[nelem], Mc[nelem], Nc[nelem];
  78   if (qpms_planewave2vswf_fill_sph(ks, E0s, Lc, Mc, Nc, lMax, norm)) abort();
  79
  80
  81   // each ray will have its file
  82   FILE *rfile[nrays];
  83   {
  84     char fname[strlen(argv[1])+20];
  85     for (int p = 0; p < nrays; ++p) {
  86       sprintf(fname, "%s.%d", argv[1], p);
  87       rfile[p] = fopen(fname, "w");
  88       fputs("##\tnorm\tlMax\tnelem\n", rfile[p]);
  89       fprintf(rfile[p], "#\t%d\t%d\t%d\n",
  90           (int)norm, (int)lMax, (int)nelem);
  91       fputs("##\tkx\tky\tkz\tR(E0x)\tI(E0x)\tR(E0y)\tI(E0y)\tR(E0z)\tI(E0z)\n", rfile[p]);
  92       fprintf(rfile[p], "#\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",
  93           k.x, k.y, k.z,
  94           creal(E0.x), cimag(E0.x),
  95           creal(E0.y), cimag(E0.y),
  96           creal(E0.z), cimag(E0.z)
  97       );
  98       sph_t ks = cart2sph(k);
  99       csphvec_t E0s = ccart2csphvec(E0, ks);
 100       fputs("##\tkr\tkθ\tkφ\tR(E0r)\tI(E0r)\tR(E0θ)\tI(E0θ)\tR(E0φ)\tI(E0φ)\n", rfile[p]);
 101       fprintf(rfile[p], "#\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",
 102           ks.r, ks.theta, ks.phi,
 103           creal(E0s.rc), cimag(E0s.rc),
 104           creal(E0s.thetac), cimag(E0s.thetac),
 105           creal(E0s.phic), cimag(E0s.phic)
 106       );
 107       fprintf(rfile[p], "#\tThis ray cart:\t%g\t%g\t%g\n",
 108           rays[p].x, rays[p].y, rays[p].z);
 109       sph_t ps = cart2sph(rays[p]);
 110       fprintf(rfile[p], "#\tThis ray sph:\t%g\t%g\t%g\n",
 111           ps.r, ps.theta, ps.phi);
 112       fputs("## TODO Coefficients of the wave decomposition\n", rfile[p]); // TODO print Lc, Mc, Nc
 113       fputs("\t"
 114           "\t\t\t\t\t\t" "\t"
 115           "\t\t\t\t\t\t"
 116           "\t\t\t\t\t\t"
 117           "\t\t\t\t\t\t" "\t\t\t\t\t\t", rfile[p]);
 118       for(qpms_y_t y = 0; y < nelem; ++y){ // TODO print also the longitudinal waves?
 119         for(int i = 0; i < 3; ++i)
 120           fprintf(rfile[p], "%g\t%g\t", creal(Mc[y]), cimag(Mc[y]));
 121         for(int i = 0; i < 3; ++i)
 122           fprintf(rfile[p], "%g\t%g\t", creal(Nc[y]), cimag(Nc[y]));
 123       }
 124       fputc('\n', rfile[p]);
 125
 126       fputs("#step\t"
 127         "x\ty\tz\tr\tθ\tφ\tphase\t"
 128         "R(Ex)\tI(Ex)\tR(Ey)\tI(Ey)\tR(Ez)\tI(Ez)\t"
 129         "R(Er)\tI(Er)\tR(Eθ)\tI(Eθ)\tR(Eφ)\tI(Eφ)\t"
 130         "R(rEx)\tI(rEx)\tR(rEy)\tI(rEy)\tR(rEz)\tI(rEz)\t"
 131         "R(rEr)\tI(rEr)\tR(rEθ)\tI(rEθ)\tR(rEφ)\tI(rEφ)\t",
 132         rfile[p]
 133       );
 134
 135
 136       for(qpms_y_t y = 0; y < nelem; ++y) {
 137         qpms_l_t l2; qpms_m_t m2;
 138         qpms_y2mn_p(y, &m2, &l2);
 139         fprintf(rfile[p], // TODO print also the longitudinal waves?
 140           "R(M2(%d,%d).r)\tI(M2(%d,%d).r)\t"
 141           "R(M2(%d,%d).θ)\tI(M2(%d,%d).θ)\t"
 142           "R(M2(%d,%d).φ)\tI(M2(%d,%d).φ)\t"
 143           "R(N2(%d,%d).r)\tI(N2(%d,%d).r)\t"
 144           "R(N2(%d,%d).θ)\tI(N2(%d,%d).θ)\t"
 145           "R(N2(%d,%d).φ)\tI(N2(%d,%d).φ)\t",
 146           l2,m2,l2,m2, l2,m2,l2,m2, l2,m2,l2,m2,
 147           l2,m2,l2,m2, l2,m2,l2,m2, l2,m2,l2,m2
 148         );
 149       }
 150       fputc('\n', rfile[p]);
 151     }
 152   }
 153
 154
 155
 156   for (int step = 0; step <= nraysteps; ++step)  {
 157     double stepratio = (double) step / (double) nraysteps;
 158     for(int p = 0; p < nrays; ++p) {
 159       cart3_t cart = cart3_scale(stepratio, rays[p]);
 160       sph_t sph = cart2sph(cart);
 161       double phase = cart3_dot(k, cart);
 162       ccart3_t E = ccart3_scale(cexp(I*phase), E0);
 163       csphvec_t L[nelem], M[nelem], N[nelem], Es, rEs;
 164       Es = ccart2csphvec(E, sph);
 165       rEs = qpms_eval_vswf(sph, Lc, Mc, Nc, lMax, QPMS_BESSEL_REGULAR, norm); // TODO maybe also check this by summing manually Lc * L etc.
 166       ccart3_t rEc = csphvec2ccart(rEs, sph);
 167       if (qpms_vswf_fill(L, M, N, lMax, sph, QPMS_BESSEL_REGULAR, norm)) abort();
 168
 169       fprintf(rfile[p],"%d\t"
 170           "%g\t%g\t%g\t%g\t%g\t%g\t" "%g\t"
 171           "%g\t%g\t%g\t%g\t%g\t%g\t"
 172           "%g\t%g\t%g\t%g\t%g\t%g\t"
 173           "%g\t%g\t%g\t%g\t%g\t%g\t"
 174           "%g\t%g\t%g\t%g\t%g\t%g\t",
 175           step, cart.x, cart.y, cart.z, sph.r, sph.theta, sph.phi, phase,
 176           creal(E.x), cimag(E.x), creal(E.y), cimag(E.y), creal(E.z), cimag(E.z),
 177           creal(Es.rc), cimag(Es.rc), creal(Es.thetac), cimag(Es.thetac), creal(Es.phic), cimag(Es.phic),
 178           creal(rEc.x), cimag(rEc.x), creal(rEc.y), cimag(rEc.y), creal(rEc.z), cimag(rEc.z),
 179           creal(rEs.rc), cimag(rEs.rc), creal(rEs.thetac), cimag(rEs.thetac), creal(rEs.phic), cimag(rEs.phic)
 180       );
 181
 182       for(qpms_y_t y = 0; y < nelem; ++y) {
 183         csphvec_t /*Lsc,*/ Msc, Nsc; // TODO print also the longitudinal waves?
 184         //Lsc = csphvec_scale(Lc[y], L);
 185         Msc = M[y];//csphvec_scale(Mc[y], M[y]);
 186         Nsc = N[y];//csphvec_scale(Nc[y], N[y]);
 187         fprintf(rfile[p], "%g\t%g\t%g\t%g\t%g\t%g\t""%g\t%g\t%g\t%g\t%g\t%g\t",
 188             creal(Msc.rc), cimag(Msc.rc),
 189             creal(Msc.thetac), cimag(Msc.thetac),
 190             creal(Msc.phic), cimag(Msc.phic),
 191             creal(Nsc.rc), cimag(Nsc.rc),
 192             creal(Nsc.thetac), cimag(Nsc.thetac),
 193             creal(Nsc.phic), cimag(Nsc.phic)
 194         );
 195       }
 196       fputc('\n', rfile[p]);
 197     }
 198   }
 199
 200   for (int p = 0; p < nrays; ++p)
 201     fclose(rfile[p]);
 202         gsl_rng_free(rng);
 203 }
 204