qpms/vswf.c

   1 #include <math.h>
   2 #include <gsl/gsl_math.h>
   3 #include "assert_cython_workaround.h"
   4 #include "vswf.h"
   5 #include "indexing.h"
   6 #include "translations.h" // TODO move qpms_sph_bessel_fill elsewhere
   7 #include "qpms_specfunc.h"
   8 #include <stdlib.h>
   9 #include <string.h>
  10 #include "qpms_error.h"
  11 #include "normalisation.h"
  12
  13
  14 qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() {
  15   qpms_vswf_set_spec_t *s = calloc(1,sizeof(qpms_vswf_set_spec_t));
  16   if (s == NULL) return NULL; // TODO warn
  17   // The rest are zeroes automatically because of calloc:
  18   s->lMax_L = -1;
  19   return s;
  20 }
  21
  22 #define MAX(x,y) (((x)<(y))?(y):(x))
  23
  24 qpms_errno_t qpms_vswf_set_spec_append(qpms_vswf_set_spec_t *s, const qpms_uvswfi_t u) {
  25   qpms_l_t l;
  26   qpms_m_t m;
  27   qpms_vswf_type_t t;
  28   if (qpms_uvswfi2tmn(u, &t, &m, &l)!=QPMS_SUCCESS) return QPMS_ERROR; // TODO WARN
  29   if (s->n + 1 > s->capacity) {
  30     size_t newcap = (s->capacity < 32) ? 32 : 2*s->capacity;
  31     qpms_uvswfi_t *newmem = realloc(s->ilist, newcap * sizeof(qpms_uvswfi_t));
  32     if (newmem == NULL) return QPMS_ENOMEM; // TODO WARN
  33     s->capacity = newcap;
  34     s->ilist = newmem;
  35   }
  36   s->ilist[s->n] = u;
  37   ++s->n;
  38   switch(t) {
  39     case QPMS_VSWF_ELECTRIC:
  40       s->lMax_N = MAX(s->lMax_N, l);
  41       break;
  42     case QPMS_VSWF_MAGNETIC:
  43       s->lMax_M = MAX(s->lMax_M, l);
  44       break;
  45     case QPMS_VSWF_LONGITUDINAL:
  46       s->lMax_L = MAX(s->lMax_L, l);
  47       break;
  48     default:
  49       QPMS_WTF;
  50   }
  51   s->lMax = MAX(s->lMax, l);
  52   return QPMS_SUCCESS;
  53 }
  54
  55 bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
  56                     const qpms_vswf_set_spec_t *b) {
  57   if (a == b) return true;
  58   if (a->norm != b->norm) return false;
  59   if (a->n != b->n) return false;
  60   for (size_t i = 0; i < a->n; ++i)
  61     if (a->ilist[i] != b->ilist[i])
  62       return false;
  63   return true;
  64 }
  65
  66 qpms_vswf_set_spec_t *qpms_vswf_set_spec_copy(const qpms_vswf_set_spec_t *or){
  67   qpms_vswf_set_spec_t *c;
  68   QPMS_CRASHING_MALLOC(c, sizeof(*c));
  69   *c = *or;
  70   c->ilist = malloc(sizeof(qpms_uvswfi_t) * c->n);
  71   memcpy(c->ilist, or->ilist, sizeof(qpms_uvswfi_t)*c->n);
  72   c->capacity = c->n;
  73   return c;
  74 }
  75
  76 qpms_vswf_set_spec_t *qpms_vswf_set_spec_from_lMax(qpms_l_t lMax,
  77     qpms_normalisation_t norm) {
  78   qpms_vswf_set_spec_t *c;
  79   QPMS_CRASHING_MALLOC(c, sizeof(*c));
  80   c->n = c->capacity = 2 * qpms_lMax2nelem(lMax);
  81   c->ilist = malloc(sizeof(qpms_uvswfi_t) * c->capacity);
  82   size_t i = 0;
  83   for (int it = 0; it < 2; ++it)
  84     for (qpms_l_t n = 1; n <= lMax; ++n)
  85       for (qpms_m_t m = -n; m <= n; ++m)
  86         c->ilist[i++] =
  87           qpms_tmn2uvswfi(it ? QPMS_VSWF_MAGNETIC : QPMS_VSWF_ELECTRIC, m, n);
  88   c->norm = norm;
  89   c->lMax = c->lMax_M = c->lMax_N = lMax;
  90   c->lMax_L = -1;
  91   return c;
  92 }
  93
  94 void qpms_vswf_set_spec_free(qpms_vswf_set_spec_t *s) {
  95   if(s) free(s->ilist);
  96   free(s);
  97 }
  98
  99 struct bspec_reindex_pair {
 100   qpms_uvswfi_t ui;
 101   size_t i_orig;
 102 };
 103
 104 static int cmp_bspec_reindex_pair(const void *aa, const void *bb) {
 105   const struct bspec_reindex_pair *a = aa, *b = bb;
 106   if (a->ui < b->ui) return -1;
 107   else if (a->ui == b->ui) return 0;
 108   else return 1;
 109 }
 110
 111 size_t *qpms_vswf_set_reindex(const qpms_vswf_set_spec_t *small, const qpms_vswf_set_spec_t *big) {
 112   QPMS_UNTESTED;
 113   struct bspec_reindex_pair *small_pairs, *big_pairs;
 114   size_t *r;
 115   QPMS_CRASHING_MALLOC(small_pairs, sizeof(struct bspec_reindex_pair) * small->n);
 116   QPMS_CRASHING_MALLOC(big_pairs, sizeof(struct bspec_reindex_pair) * big->n);
 117   QPMS_CRASHING_MALLOC(r, sizeof(size_t) * small->n);
 118   for(size_t i = 0; i < small->n; ++i) {
 119     small_pairs[i].ui = small->ilist[i];
 120     small_pairs[i].i_orig = i;
 121   }
 122   for(size_t i = 0 ; i < big->n; ++i) {
 123     big_pairs[i].ui = big->ilist[i];
 124     big_pairs[i].i_orig = i;
 125   }
 126   qsort(small_pairs, small->n, sizeof(struct bspec_reindex_pair), cmp_bspec_reindex_pair);
 127   qsort(big_pairs, big->n, sizeof(struct bspec_reindex_pair), cmp_bspec_reindex_pair);
 128
 129   size_t bi = 0;
 130   for(size_t si = 0; si < small->n; ++si) {
 131     while(big_pairs[bi].ui < small_pairs[si].ui)
 132       ++bi;
 133     if(big_pairs[bi].ui == small_pairs[si].ui)
 134       r[small_pairs[si].i_orig] = big_pairs[bi].i_orig;
 135     else
 136       r[small_pairs[si].i_orig] = ~(size_t)0;
 137   }
 138
 139   free(small_pairs);
 140   free(big_pairs);
 141   return r;
 142 }
 143
 144 csphvec_t qpms_vswf_single_el_csph(qpms_m_t m, qpms_l_t l, csph_t kdlj,
 145     qpms_bessel_t btyp, qpms_normalisation_t norm) {
 146   lmcheck(l,m);
 147   csphvec_t N;
 148   complex double *bessel;
 149   QPMS_CRASHING_MALLOC(bessel,(l+1)*sizeof(complex double));
 150   QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel));
 151   qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, qpms_normalisation_t_csphase(norm));
 152
 153   complex double eimf = qpms_spharm_azimuthal_part(norm, m, kdlj.phi);
 154   complex double d_eimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, kdlj.phi);
 155   qpms_y_t y = qpms_mn2y(m,l);
 156
 157   N.rc = l*(l+1) * pt.leg[y] * bessel[l] / kdlj.r * eimf;
 158   complex double besselfac = bessel[l-1] - l * bessel[l] / kdlj.r;
 159   N.thetac = pt.tau[y] * besselfac * eimf;
 160   N.phic = pt.pi[y] * besselfac * d_eimf_dmf;
 161
 162   N = csphvec_scale(qpms_normalisation_factor_N_noCS(norm, l, m), N);
 163
 164   qpms_pitau_free(pt);
 165   free(bessel);
 166   return N;
 167 }
 168
 169 csphvec_t qpms_vswf_single_mg_csph(qpms_m_t m, qpms_l_t l, csph_t kdlj,
 170     qpms_bessel_t btyp, qpms_normalisation_t norm) {
 171   lmcheck(l,m);
 172   csphvec_t M;
 173   complex double *bessel;
 174   QPMS_CRASHING_MALLOC(bessel,(l+1)*sizeof(complex double));
 175   QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel));
 176   qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, qpms_normalisation_t_csphase(norm));
 177   complex double eimf = qpms_spharm_azimuthal_part(norm, m, kdlj.phi);
 178   complex double d_eimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, kdlj.phi);
 179   qpms_y_t y = qpms_mn2y(m,l);
 180
 181   M.rc = 0.;
 182   M.thetac = pt.pi[y] * bessel[l] * d_eimf_dmf;
 183   M.phic = -pt.tau[y] * bessel[l] * eimf;
 184
 185   M = csphvec_scale(qpms_normalisation_factor_M_noCS(norm, l, m), M);
 186
 187   qpms_pitau_free(pt);
 188   free(bessel);
 189   return M;
 190 }
 191
 192 csphvec_t qpms_vswf_single_el(qpms_m_t m, qpms_l_t l, sph_t kdlj,
 193     qpms_bessel_t btyp, qpms_normalisation_t norm) {
 194   return qpms_vswf_single_el_csph(m, l, sph2csph(kdlj), btyp, norm);
 195 }
 196
 197 csphvec_t qpms_vswf_single_mg(qpms_m_t m, qpms_l_t l, sph_t kdlj,
 198     qpms_bessel_t btyp, qpms_normalisation_t norm) {
 199   return qpms_vswf_single_mg_csph(m, l, sph2csph(kdlj), btyp, norm);
 200 }
 201
 202 qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
 203     qpms_bessel_t btyp, qpms_normalisation_t norm) {
 204   qpms_vswfset_sph_t *res = malloc(sizeof(qpms_vswfset_sph_t));
 205   res->lMax = lMax;
 206   qpms_y_t nelem = qpms_lMax2nelem(lMax);
 207   QPMS_CRASHING_MALLOC(res->el, sizeof(csphvec_t)*nelem);
 208   QPMS_CRASHING_MALLOC(res->mg, sizeof(csphvec_t)*nelem);
 209   QPMS_ENSURE_SUCCESS(qpms_vswf_fill(NULL, res->mg, res->el, lMax, kdlj, btyp, norm));
 210   return res;
 211 }
 212
 213 void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *w) {
 214   assert(NULL != w && NULL != w->el && NULL != w->mg);
 215   free(w->el);
 216   free(w->mg);
 217   free(w);
 218 }
 219
 220 csphvec_t qpms_vswf_L00(csph_t kr, qpms_bessel_t btyp,
 221     qpms_normalisation_t norm) {
 222   QPMS_UNTESTED;
 223   // CHECKME Is it OK to ignore norm?? (Is L_0^0 the same for all conventions?)
 224   complex double bessel0;
 225   QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, 0, kr.r, &bessel0));
 226   csphvec_t result = {0.25 * M_2_SQRTPI * bessel0, 0, 0};
 227   return result;
 228 }
 229
 230 qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget,
 231     csphvec_t * const mgtarget, csphvec_t * const eltarget, qpms_l_t lMax,
 232     csph_t kr, qpms_bessel_t btyp, const qpms_normalisation_t norm) {
 233   assert(lMax >= 1);
 234   complex double *bessel = malloc((lMax+1)*sizeof(complex double));
 235   qpms_errno_t bessel_retval = qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel);
 236   QPMS_ENSURE_SUCCESS_OR(bessel_retval, QPMS_ESING);
 237   qpms_pitau_t pt = qpms_pitau_get(kr.theta, lMax, qpms_normalisation_t_csphase(norm));
 238   complex double const *pbes = bessel + 1; // starting from l = 1
 239   double const *pleg = pt.leg;
 240   double const *ppi = pt.pi;
 241   double const *ptau = pt.tau;
 242   csphvec_t *plong = longtarget, *pmg = mgtarget, *pel = eltarget;
 243   for(qpms_l_t l = 1; l <= lMax; ++l) {
 244     complex double besfac;
 245     complex double besderfac;
 246     if (kr.r) {
 247       besfac = *pbes / kr.r;
 248     } else {
 249       besfac = (1 == l) ? 1/3. : 0;
 250     }
 251     besderfac = *(pbes-1) - l * besfac;
 252     for(qpms_m_t m = -l; m <= l; ++m) {
 253       complex double eimf = qpms_spharm_azimuthal_part(norm, m, kr.phi);
 254       complex double d_eimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, kr.phi);
 255       if (longtarget) { QPMS_UNTESTED;
 256         double longfac = sqrt(l*(l+1));
 257         plong->rc = // FATAL FIXME: I get wrong result here for plane wave re-expansion
 258           // whenever kr.r > 0 (for waves with longitudinal component, ofcoz)
 259           /*(*(pbes-1) - (l+1)/kr.r* *pbes)*/
 260           (besderfac-besfac)
 261           * (*pleg) * longfac * eimf;
 262         plong->thetac = *ptau * besfac * longfac * eimf;
 263         plong->phic = *ppi * besfac * longfac * d_eimf_dmf;
 264         *plong = csphvec_scale(qpms_normalisation_factor_L_noCS(norm, l, m), *plong);
 265         ++plong;
 266       }
 267       if (eltarget) {
 268         pel->rc = l*(l+1) * (*pleg) * besfac * eimf;
 269         pel->thetac = *ptau * besderfac * eimf;
 270         pel->phic = *ppi * besderfac * d_eimf_dmf;
 271         *pel = csphvec_scale(qpms_normalisation_factor_N_noCS(norm, l, m), *pel);
 272         ++pel;
 273       }
 274       if (mgtarget) {
 275         pmg->rc = 0.;
 276         pmg->thetac = *ppi * (*pbes) * d_eimf_dmf;
 277         pmg->phic = - *ptau * (*pbes) * eimf;
 278         *pmg = csphvec_scale(qpms_normalisation_factor_M_noCS(norm, l, m), *pmg);
 279         ++pmg;
 280       }
 281       ++pleg; ++ppi; ++ptau;
 282     }
 283     ++pbes;
 284   }
 285   free(bessel);
 286   qpms_pitau_free(pt);
 287   return bessel_retval;
 288 }
 289
 290 qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget,
 291     csphvec_t * const mgtarget, csphvec_t * const eltarget, qpms_l_t lMax,
 292     sph_t kr, qpms_bessel_t btyp, qpms_normalisation_t norm) {
 293   csph_t krc = {kr.r, kr.theta, kr.phi};
 294   return qpms_vswf_fill_csph(longtarget, mgtarget, eltarget, lMax,
 295       krc, btyp, norm);
 296 }
 297
 298 // consistency check: this should give the same results as the above function (up to rounding errors)
 299 qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t * const mgtarget, csphvec_t * const eltarget,
 300     qpms_l_t lMax, sph_t kr,
 301     qpms_bessel_t btyp, qpms_normalisation_t norm) {
 302   assert(lMax >= 1);
 303   complex double *bessel = malloc((lMax+1)*sizeof(complex double));
 304   QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel));
 305   complex double const *pbes = bessel + 1; // starting from l = 1
 306
 307   qpms_y_t nelem = qpms_lMax2nelem(lMax);
 308   csphvec_t *a;
 309   QPMS_CRASHING_MALLOC(a, 3*nelem*sizeof(csphvec_t))
 310   csphvec_t * const a1 = a, * const a2 = a1 + nelem, * const a3 = a2 + 2 * nelem;
 311   QPMS_ENSURE_SUCCESS(qpms_vecspharm_fill(a1, a2, a3, lMax, kr, norm));
 312   const csphvec_t *p1 = a1;
 313   const csphvec_t *p2 = a2;
 314   const csphvec_t *p3 = a3;
 315
 316   csphvec_t *plong = longtarget, *pmg = mgtarget, *pel = eltarget;
 317   for(qpms_l_t l = 1; l <= lMax; ++l) {
 318     complex double besfac = *pbes / kr.r;
 319     complex double besderfac = *(pbes-1) - l * besfac;
 320     double sqrtlfac = sqrt(l*(l+1));
 321     for(qpms_m_t m = -l; m <= l; ++m) {
 322       if (longtarget) {
 323         complex double L2Nfac = qpms_normalisation_factor_L_noCS(norm, l, m)
 324           / qpms_normalisation_factor_N_noCS(norm, l, m);
 325         *plong = csphvec_add(csphvec_scale(besderfac-besfac, *p3),
 326             csphvec_scale(sqrtlfac * besfac, *p2));
 327         *plong = csphvec_scale(L2Nfac, *plong);
 328         ++plong;
 329       }
 330       if (eltarget) {
 331         *pel = csphvec_add(csphvec_scale(besderfac, *p2),
 332             csphvec_scale(sqrtlfac * besfac, *p3));
 333         ++pel;
 334       }
 335       if (mgtarget) {
 336         *pmg = csphvec_scale(*pbes, *p1);
 337         ++pmg;
 338       }
 339       ++p1; ++p2; ++p3;
 340     }
 341     ++pbes;
 342   }
 343   free(a);
 344   free(bessel);
 345   return QPMS_SUCCESS;
 346 }
 347
 348 qpms_errno_t qpms_vecspharm_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target,
 349     qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm) {
 350   assert(lMax >= 1);
 351   qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, qpms_normalisation_t_csphase(norm));
 352   double const *pleg = pt.leg;
 353   double const *ppi = pt.pi;
 354   double const *ptau = pt.tau;
 355   csphvec_t *p1 = a1target, *p2 = a2target, *p3 = a3target;
 356   for (qpms_l_t l = 1; l <= lMax; ++l) {
 357     for(qpms_m_t m = -l; m <= l; ++m) {
 358       const complex double Mfac = qpms_normalisation_factor_M_noCS(norm, l, m);
 359       const complex double Nfac = qpms_normalisation_factor_N_noCS(norm, l, m);
 360       const complex double eimf = qpms_spharm_azimuthal_part(norm, m, dir.phi);
 361       const complex double deimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, dir.phi);
 362       if (a1target) {
 363         p1->rc = 0;
 364         p1->thetac = *ppi * deimf_dmf * Mfac;
 365         p1->phic = -*ptau * eimf * Mfac;
 366         ++p1;
 367       }
 368       if (a2target) {
 369         p2->rc = 0;
 370         p2->thetac = *ptau * eimf * Nfac;
 371         p2->phic = *ppi * deimf_dmf * Nfac;
 372         ++p2;
 373       }
 374       if (a3target) {
 375         p3->rc = sqrt(l*(l+1)) * (*pleg) * eimf * Nfac;
 376         p3->thetac = 0;
 377         p3->phic = 0;
 378         ++p3;
 379       }
 380       ++pleg; ++ppi; ++ptau;
 381     }
 382   }
 383   qpms_pitau_free(pt);
 384   return QPMS_SUCCESS;
 385 }
 386
 387 qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target,
 388     qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm) {
 389 #if 1
 390   return qpms_vecspharm_fill(a1target, a2target, a3target, lMax, dir,
 391       qpms_normalisation_dual(norm));
 392 #else
 393   assert(lMax >= 1);
 394   qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, norm);
 395   double const *pleg = pt.leg;
 396   double const *ppi = pt.pi;
 397   double const *ptau = pt.tau;
 398   csphvec_t *p1 = a1target, *p2 = a2target, *p3 = a3target;
 399   for(qpms_l_t l = 1; l <= lMax; ++l) {
 400     for(qpms_m_t m = -l; m <= l; ++m) {
 401       double normfac = 1./qpms_normalisation_t_factor_abssquare(norm, l, m); // factor w.r.t. Kristensson
 402       complex double eimf = cexp(m * dir.phi * I);
 403       if (a1target) {
 404         p1->rc = 0;
 405         p1->thetac = conj(*ppi * normfac * I * eimf);
 406         p1->phic = conj(-*ptau * normfac * eimf);
 407         ++p1;
 408       }
 409       if (a2target) {
 410         p2->rc = 0;
 411         p2->thetac = conj(*ptau * normfac * eimf);
 412         p2->phic = conj(*ppi * normfac * I * eimf);
 413         ++p2;
 414       }
 415       if (a3target) {
 416         p3->rc = conj(sqrt(l*(l+1)) * (*pleg) * normfac * eimf);
 417         p3->thetac = 0;
 418         p3->phic = 0;
 419         ++p3;
 420       }
 421       ++pleg; ++ppi; ++ptau;
 422     }
 423   }
 424   qpms_pitau_free(pt);
 425   return QPMS_SUCCESS;
 426 #endif
 427 }
 428
 429
 430 static inline complex double ipowl(qpms_l_t l) {
 431   switch(l % 4) {
 432     case 0: return 1;
 433             break;
 434     case 1: return I;
 435             break;
 436     case 2: return -1;
 437             break;
 438     case 3: return -I;
 439             break;
 440     default: QPMS_WTF;
 441   }
 442   QPMS_WTF;
 443 }
 444
 445 qpms_errno_t qpms_planewave2vswf_fill_sph(sph_t wavedir, csphvec_t amplitude,
 446     complex double *target_longcoeff, complex double *target_mgcoeff,
 447     complex double *target_elcoeff, qpms_l_t lMax, qpms_normalisation_t norm) {
 448   qpms_y_t nelem = qpms_lMax2nelem(lMax);
 449   csphvec_t * const dual_A1 = malloc(3*nelem*sizeof(csphvec_t)), *const dual_A2 = dual_A1 + nelem,
 450             * const dual_A3 = dual_A2 + nelem;
 451   QPMS_ENSURE_SUCCESS(qpms_vecspharm_dual_fill(dual_A1, dual_A2, dual_A3, lMax, wavedir, norm))
 452   const csphvec_t *pA1 = dual_A1, *pA2 = dual_A2, *pA3 = dual_A3;
 453   complex double *plong = target_longcoeff, *pmg = target_mgcoeff, *pel = target_elcoeff;
 454   for (qpms_l_t l = 1; l <= lMax; ++l) {
 455     complex double prefac1 = 4 * M_PI * ipowl(l);
 456     complex double prefac23 = - 4 * M_PI * ipowl(l+1);
 457     for (qpms_m_t m = -l; m <= l; ++m) {
 458       if(target_longcoeff) *plong = prefac23 * csphvec_dotnc(*pA3, amplitude);
 459       if(target_mgcoeff) *pmg = prefac1 * csphvec_dotnc(*pA1, amplitude);
 460       if(target_elcoeff) *pel = prefac23 * csphvec_dotnc(*pA2, amplitude);
 461       ++pA1; ++pA2; ++pA3; ++plong; ++pmg; ++pel;
 462     }
 463
 464   }
 465   free(dual_A1);
 466   return QPMS_SUCCESS;
 467 }
 468
 469 qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir_cart /*allow complex k?*/, ccart3_t amplitude_cart,
 470     complex double * const longcoeff, complex double * const mgcoeff,
 471     complex double * const elcoeff, qpms_l_t lMax, qpms_normalisation_t norm)
 472 {
 473
 474   sph_t wavedir_sph = cart2sph(wavedir_cart);
 475   csphvec_t amplitude_sphvec = ccart2csphvec(amplitude_cart, wavedir_sph);
 476   return qpms_planewave2vswf_fill_sph(wavedir_sph, amplitude_sphvec,
 477       longcoeff, mgcoeff, elcoeff, lMax, norm);
 478 }
 479
 480 qpms_errno_t qpms_incfield_planewave(complex double *target, const qpms_vswf_set_spec_t *bspec,
 481     const cart3_t evalpoint, const void *args, bool add) {
 482   QPMS_UNTESTED;
 483   const qpms_incfield_planewave_params_t *p = args;
 484
 485   const ccart3_t k_cart = p->use_cartesian ? p->k.cart : csph2ccart(p->k.sph);
 486   const complex double phase = ccart3_dotnc(k_cart, cart32ccart3(evalpoint));
 487   if(cimag(phase))
 488     QPMS_INCOMPLETE_IMPLEMENTATION("Complex-valued wave vector not implemented correctly; cf. docs.");
 489   const complex double phasefac = cexp(I*phase);
 490
 491   // Throw away the imaginary component; TODO handle it correctly
 492   const sph_t k_sph = csph2sph(p->use_cartesian ? ccart2csph(p->k.cart) : p->k.sph);
 493   const csphvec_t E_sph = csphvec_scale(phasefac,
 494       p->use_cartesian ? ccart2csphvec(p->E.cart, k_sph) : p->E.sph);
 495
 496   complex double *lc = NULL, *mc = NULL, *nc = NULL;
 497   const qpms_y_t nelem = qpms_lMax2nelem(bspec->lMax);
 498   if (bspec->lMax_L > 0)  QPMS_CRASHING_MALLOC(lc, nelem * sizeof(complex double));
 499   if (bspec->lMax_M > 0)  QPMS_CRASHING_MALLOC(mc, nelem * sizeof(complex double));
 500   if (bspec->lMax_N > 0)  QPMS_CRASHING_MALLOC(nc, nelem * sizeof(complex double));
 501
 502   qpms_errno_t retval = qpms_planewave2vswf_fill_sph(k_sph, E_sph, lc, mc, nc,
 503       bspec->lMax, bspec->norm);
 504
 505   if (!add) memset(target, 0, bspec->n * sizeof(complex double));
 506   for (size_t i = 0; i < bspec->n; ++i) {
 507     const qpms_uvswfi_t ui = bspec->ilist[i];
 508     if (ui == QPMS_UI_L00) // for l = 0 the coefficient is zero due to symmetry (for real wave vector)
 509       target[i] = 0;
 510     else {
 511       qpms_vswf_type_t t; qpms_y_t y;
 512       QPMS_ENSURE_SUCCESS(qpms_uvswfi2ty_l(ui, &t, &y));
 513       switch(t) {
 514         case QPMS_VSWF_ELECTRIC:
 515           target[i] += nc[y];
 516           break;
 517         case QPMS_VSWF_MAGNETIC:
 518           target[i] += mc[y];
 519           break;
 520         case QPMS_VSWF_LONGITUDINAL:
 521           target[i] += lc[y];
 522           break;
 523         default:
 524           QPMS_WTF;
 525       }
 526     }
 527   }
 528   free(lc); free(mc); free(nc);
 529   return retval;
 530 }
 531
 532
 533 csphvec_t qpms_eval_vswf_csph(csph_t kr,
 534     complex double * const lc, complex double *const mc, complex double *const ec,
 535     qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm)
 536 {
 537   qpms_y_t nelem = qpms_lMax2nelem(lMax);
 538   csphvec_t lsum, msum, esum, lcomp, mcomp, ecomp;
 539   csphvec_kahaninit(&lsum, &lcomp);
 540   csphvec_kahaninit(&msum, &mcomp);
 541   csphvec_kahaninit(&esum, &ecomp);
 542   csphvec_t *lset = NULL, *mset = NULL, *eset = NULL;
 543   if(lc) lset = malloc(nelem * sizeof(csphvec_t));
 544   if(mc) mset = malloc(nelem * sizeof(csphvec_t));
 545   if(ec) eset = malloc(nelem * sizeof(csphvec_t));
 546   qpms_vswf_fill_csph(lset, mset, eset, lMax, kr, btyp, norm);
 547   if(lc) for(qpms_y_t y = 0; y < nelem; ++y)
 548     csphvec_kahanadd(&lsum, &lcomp, csphvec_scale(lc[y], lset[y]));
 549   if(mc) for(qpms_y_t y = 0; y < nelem; ++y)
 550     csphvec_kahanadd(&msum, &mcomp, csphvec_scale(mc[y], mset[y]));
 551   if(ec) for(qpms_y_t y = 0; y < nelem; ++y)
 552     csphvec_kahanadd(&esum, &ecomp, csphvec_scale(ec[y], eset[y]));
 553   if(lc) free(lset);
 554   if(mc) free(mset);
 555   if(ec) free(eset);
 556   //return csphvec_add(esum, csphvec_add(msum, lsum));
 557   csphvec_kahanadd(&esum, &ecomp, msum);
 558   csphvec_kahanadd(&esum, &ecomp, lsum);
 559   return esum;
 560 }
 561
 562 csphvec_t qpms_eval_vswf(sph_t kr,
 563     complex double * const lc, complex double *const mc, complex double *const ec,
 564     qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm) {
 565   csph_t krc = {kr.r, kr.theta, kr.phi};
 566   return qpms_eval_vswf_csph(krc, lc, mc, ec, lMax, btyp, norm);
 567 }
 568
 569 qpms_errno_t qpms_uvswf_fill(csphvec_t *const target, const qpms_vswf_set_spec_t *bspec,
 570     csph_t kr, qpms_bessel_t btyp) {
 571   QPMS_UNTESTED;
 572   QPMS_ENSURE(target, "Target array pointer must not be NULL.");
 573   csphvec_t *el = NULL, *mg = NULL, *lg = NULL;
 574   const qpms_y_t nelem = qpms_lMax2nelem(bspec->lMax);
 575   if (bspec->lMax_L > 0)
 576     QPMS_CRASHING_MALLOC(lg, nelem * sizeof(csphvec_t));
 577   if (bspec->lMax_M > 0)
 578     QPMS_CRASHING_MALLOC(mg, nelem * sizeof(csphvec_t));
 579   if (bspec->lMax_N > 0)
 580     QPMS_CRASHING_MALLOC(el, nelem * sizeof(csphvec_t));
 581   qpms_errno_t retval = qpms_vswf_fill_csph(lg, mg, el, bspec->lMax, kr, btyp, bspec->norm);
 582   for (size_t i = 0; i < bspec->n; ++i) {
 583     const qpms_uvswfi_t ui = bspec->ilist[i];
 584     if (ui == QPMS_UI_L00) // l = 0 longitudinal wave must be calculated separately.
 585       target[i] = qpms_vswf_L00(kr, btyp, bspec->norm);
 586     else {
 587       qpms_vswf_type_t t; qpms_y_t y;
 588       QPMS_ENSURE_SUCCESS(qpms_uvswfi2ty_l(ui, &t, &y));
 589       switch(t) {
 590         case QPMS_VSWF_ELECTRIC:
 591           target[i] = el[y];
 592           break;
 593         case QPMS_VSWF_MAGNETIC:
 594           target[i] = mg[y];
 595           break;
 596         case QPMS_VSWF_LONGITUDINAL:
 597           target[i] = lg[y];
 598           break;
 599         default:
 600           QPMS_WTF;
 601       }
 602     }
 603   }
 604   free(lg);
 605   free(mg);
 606   free(el);
 607   return retval;
 608 }
 609
 610
 611 csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *bspec,
 612     const complex double *coeffs, const csph_t kr,
 613     const qpms_bessel_t btyp) {
 614   QPMS_UNTESTED;
 615   complex double *cM = NULL, *cN = NULL, *cL = NULL, cL00 = 0;
 616   if (bspec->lMax_L > 0)
 617     QPMS_CRASHING_CALLOC(cL, bspec->n, sizeof(complex double));
 618   if (bspec->lMax_M > 0)
 619     QPMS_CRASHING_CALLOC(cM, bspec->n, sizeof(complex double));
 620   if (bspec->lMax_N > 0)
 621     QPMS_CRASHING_CALLOC(cN, bspec->n, sizeof(complex double));
 622   for (size_t i = 0; i < bspec->n; ++i) {
 623     if (bspec->ilist[i] == 0) // L00, needs special care
 624       cL00 = coeffs[i];
 625     else {
 626       qpms_vswf_type_t t;
 627       qpms_y_t y;
 628       QPMS_ENSURE_SUCCESS(qpms_uvswfi2ty_l(bspec->ilist[i], &t, &y));
 629       switch(t) {
 630         case QPMS_VSWF_LONGITUDINAL:
 631           QPMS_ASSERT(cL);
 632           cL[y] = coeffs[i];
 633           break;
 634         case QPMS_VSWF_MAGNETIC:
 635           QPMS_ASSERT(cM);
 636           cM[y] = coeffs[i];
 637           break;
 638         case QPMS_VSWF_ELECTRIC:
 639           QPMS_ASSERT(cN);
 640           cN[y] = coeffs[i];
 641           break;
 642         default:
 643           QPMS_WTF;
 644       }
 645     }
 646   }
 647   csphvec_t result = qpms_eval_vswf_csph(kr, cL, cM, cN, bspec->lMax, btyp, bspec->norm);
 648   free(cM); free(cN); free(cL);
 649   if(cL00)
 650     result = csphvec_add(result,
 651        csphvec_scale(cL00, qpms_vswf_L00(kr, btyp, bspec->norm)));
 652   return result;
 653 }
 654