Implement uvswf plane wave.

Former-commit-id: 591a742362fd311d8d7fbf12b147c8324e077f3c

qpms/indexing.h

@@ -7,6 +7,8 @@
 #include "qpms_types.h"
 #include <math.h>
 
+
+
 static inline qpms_y_t qpms_mn2y(qpms_m_t m, qpms_l_t n) {
         return n * (n + 1) + m - 1;
 }
@@ -59,6 +61,9 @@ static inline qpms_uvswfi_t qpms_tmn2uvswfi(
 	return t + 4 * qpms_mn2y_sc(m, n);
 }
 
+/// Constant denoting the l=0, m=0 longitudinal spherical vector wave.
+static const qpms_uvswfi_t QPMS_UI_L00 = 0;
+
 /// Conversion from universal VSWF index u to type, order and degree.
 /** Returns a non-zero value if the u value is invalid. */
 static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
@@ -73,7 +78,7 @@ static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
 }
 
 /// Conversion from universal VSWF index u to type and the traditional layout index (\a l > 0).
-/** Does *not* allow longitudinal waves. */
+/** Does *not* allow for longitudinal waves. */
 static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u,
 		qpms_vswf_type_t *t, qpms_y_t *y) {
 	*t = u & 3;

qpms/vectors.h

@@ -143,6 +143,7 @@ static inline cart3_t pol2cart3_equator(const pol_t pol) {
 	return cart3;
 }
 
+
 /// 3D vector addition.
 static inline cart3_t cart3_add(const cart3_t a, const cart3_t b) {
 	cart3_t res = {a.x+b.x, a.y+b.y, a.z+b.z};
@@ -188,6 +189,17 @@ static inline ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) {
 	return res;
 }
 
+/// Complex 3D cartesian vector "dot product" without conjugation.
+static inline complex double ccart3_dotnc(const ccart3_t a, const ccart3_t b) {
+	return a.x * b.x + a.y * b.y + a.z * b.z;
+}
+
+/// Convert cart3_t to ccart3_t.
+static inline ccart3_t cart32ccart3(cart3_t c){
+	ccart3_t res = {c.x, c.y, c.z};
+	return res;
+}
+
 /// Complex 3D vector (geographic coordinates) addition.
 static inline csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b) {
 	csphvec_t res = {a.rc + b.rc, a.thetac + b.thetac, a.phic + b.phic};

qpms/vswf.c

@@ -399,9 +399,9 @@ qpms_errno_t qpms_planewave2vswf_fill_sph(sph_t wavedir, csphvec_t amplitude,
     complex double prefac1 = 4 * M_PI * ipowl(l);
     complex double prefac23 = - 4 * M_PI * ipowl(l+1);
     for (qpms_m_t m = -l; m <= l; ++m) {
-      *plong = prefac23 * csphvec_dotnc(*pA3, amplitude);
+      if(target_longcoeff) *plong = prefac23 * csphvec_dotnc(*pA3, amplitude);
-      *pmg = prefac1 * csphvec_dotnc(*pA1, amplitude);
+      if(target_mgcoeff) *pmg = prefac1 * csphvec_dotnc(*pA1, amplitude);
-      *pel = prefac23 * csphvec_dotnc(*pA2, amplitude);
+      if(target_elcoeff) *pel = prefac23 * csphvec_dotnc(*pA2, amplitude);
       ++pA1; ++pA2; ++pA3; ++plong; ++pmg; ++pel;
     }
 
@@ -421,6 +421,59 @@ qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir_cart /*allow complex
       longcoeff, mgcoeff, elcoeff, lMax, norm);
 }
 
+qpms_errno_t qpms_incfield_planewave(complex double *target, const qpms_vswf_set_spec_t *bspec,
+    const cart3_t evalpoint, const void *args, bool add) {
+  QPMS_UNTESTED;
+  const qpms_incfield_planewave_params_t *p = args;
+
+  const ccart3_t k_cart = p->use_cartesian ? p->k.cart : csph2ccart(p->k.sph);
+  const complex double phase = ccart3_dotnc(k_cart, cart32ccart3(evalpoint));
+  if(cimag(phase))
+    QPMS_INCOMPLETE_IMPLEMENTATION("Complex-valued wave vector not implemented correctly; cf. docs.");
+  const complex double phasefac = cexp(I*phase);
+
+  // Throw away the imaginary component; TODO handle it correctly
+  const sph_t k_sph = csph2sph(p->use_cartesian ? ccart2csph(p->k.cart) : p->k.sph);
+  const csphvec_t E_sph = csphvec_scale(phasefac, 
+      p->use_cartesian ? ccart2csphvec(p->E.cart, k_sph) : p->E.sph);
+
+  complex double *lc = NULL, *mc = NULL, *nc = NULL;
+  const qpms_y_t nelem = qpms_lMax2nelem(bspec->lMax);
+  if (bspec->lMax_L > 0)  QPMS_CRASHING_MALLOC(lc, nelem * sizeof(complex double));
+  if (bspec->lMax_M > 0)  QPMS_CRASHING_MALLOC(mc, nelem * sizeof(complex double));
+  if (bspec->lMax_N > 0)  QPMS_CRASHING_MALLOC(nc, nelem * sizeof(complex double));
+
+  qpms_errno_t retval = qpms_planewave2vswf_fill_sph(k_sph, E_sph, lc, mc, nc,
+      bspec->lMax, bspec->norm);
+
+  if (!add) memset(target, 0, bspec->n * sizeof(complex double));
+  for (size_t i = 0; i < bspec->n; ++i) {
+    const qpms_uvswfi_t ui = bspec->ilist[i];
+    if (ui == QPMS_UI_L00) // for l = 0 the coefficient is zero due to symmetry (for real wave vector)
+      target[i] = 0;
+    else {
+      qpms_vswf_type_t t; qpms_y_t y;
+      QPMS_ENSURE_SUCCESS(qpms_uvswfi2ty_l(ui, &t, &y));
+      switch(t) {
+        case QPMS_VSWF_ELECTRIC:
+          target[i] += nc[y];
+          break;
+        case QPMS_VSWF_MAGNETIC:
+          target[i] += mc[y];
+          break;
+        case QPMS_VSWF_LONGITUDINAL:
+          target[i] += lc[y];
+          break;
+        default:
+          QPMS_WTF;
+      }
+    }
+  }
+  free(lc); free(mc); free(nc);
+  return retval;
+}
+  
+
 csphvec_t qpms_eval_vswf_csph(csph_t kr,
     complex double * const lc, complex double *const mc, complex double *const ec,
     qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm)
@@ -457,6 +510,48 @@ csphvec_t qpms_eval_vswf(sph_t kr,
   return qpms_eval_vswf_csph(krc, lc, mc, ec, lMax, btyp, norm);
 }
 
+qpms_errno_t qpms_uvswf_fill(csphvec_t *const target, const qpms_vswf_set_spec_t *bspec,
+    csph_t kr, qpms_bessel_t btyp) {
+  QPMS_UNTESTED;
+  QPMS_ENSURE(target, "Target array pointer must not be NULL."); 
+  csphvec_t *el = NULL, *mg = NULL, *lg = NULL;
+  const qpms_y_t nelem = qpms_lMax2nelem(bspec->lMax);
+  if (bspec->lMax_L > 0)
+    QPMS_CRASHING_MALLOC(lg, nelem * sizeof(csphvec_t));
+  if (bspec->lMax_M > 0)
+    QPMS_CRASHING_MALLOC(mg, nelem * sizeof(csphvec_t));
+  if (bspec->lMax_N > 0)
+    QPMS_CRASHING_MALLOC(el, nelem * sizeof(csphvec_t));
+  qpms_errno_t retval = qpms_vswf_fill_csph(lg, mg, el, bspec->lMax, kr, btyp, bspec->norm);
+  for (size_t i = 0; i < bspec->n; ++i) {
+    const qpms_uvswfi_t ui = bspec->ilist[i];
+    if (ui == QPMS_UI_L00) // l = 0 longitudinal wave must be calculated separately.
+      target[i] = qpms_vswf_L00(kr, btyp, bspec->norm);
+    else {
+      qpms_vswf_type_t t; qpms_y_t y;
+      QPMS_ENSURE_SUCCESS(qpms_uvswfi2ty_l(ui, &t, &y));
+      switch(t) {
+        case QPMS_VSWF_ELECTRIC:
+          target[i] = el[y];
+          break;
+        case QPMS_VSWF_MAGNETIC:
+          target[i] = mg[y];
+          break;
+        case QPMS_VSWF_LONGITUDINAL:
+          target[i] = lg[y];
+          break;
+        default:
+          QPMS_WTF;
+      }
+    }
+  }
+  free(lg);
+  free(mg);
+  free(el);
+  return retval;
+}
+  
+
 csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *bspec,
     const complex double *coeffs, const csph_t kr,
     const qpms_bessel_t btyp) {
@@ -501,3 +596,4 @@ csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *bspec,
        csphvec_scale(cL00, qpms_vswf_L00(kr, btyp, bspec->norm)));
   return result;
 }
+

qpms/vswf.h

@@ -76,9 +76,9 @@ csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec,
 
 // --- qpms_incfield_t instances and their arguments
 
+/// Parameter structure for qpms_incfield_planewave()
 typedef struct qpms_incfield_planewane_params_t {
 	bool use_cartesian; ///< If true, wave direction k and amplitude E are specified in cartesian coordinates (via k.cart, E.cart). If false, k is specified in spherical coordinates and E are specified in the corresponding geographical coordinates (via k.sph, E.sph).
-	bool allow_longitudinal; ///< Whether to include longitudinal part or not (if k and E are not orthogonal). Usually false makes sense.
 	union {
 		ccart3_t cart;
 		csph_t sph;
@@ -91,12 +91,24 @@ typedef struct qpms_incfield_planewane_params_t {
 	
 
 /// Calculates the (regular VSWF) expansion coefficients of a plane wave.
+/** The wave amplitude and wave vector is defined by struct qpms_incfield_planewave_params_t.
+ *
+ *  If the wave vector and amplitude are not orthogonal (i.e. the plane wave is not
+ *  fully transversal) and bspec->lMaxL is non-negative, 
+ *  the corresponding longitudinal components are calculated as well.
+ *  
+ *  For complex k vectors, the implementation is not completely correct right now.
+ *  Locally, it corresponds to decomposition of a plane wave with a real \a k
+ *  (using the real part of the \a k supplied), just the whole decomposition
+ *  is modulated by the origin-dependent factor \f$ \vect E e^{i \vect k \cdot \vect r} \f$
+ *  with \f$ \vect k \f$ complex.
+ */
 qpms_errno_t qpms_incfield_planewave(
         /// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n.
         complex double *target,
         const qpms_vswf_set_spec_t *bspec,
         const cart3_t evalpoint, ///< Point at which the VSWF expansion is made.
-        const void *args, ///< Pointer to additional function-specific arguments.
+        const void *args, ///< Pointer to additional function-specific arguments (converted to (const qpms_incfield_planewave_params_t *)).
         bool add ///< If true, add to target; rewrite target if false.
 );