uvswf-indexed vswf evaluation

Former-commit-id: adfc6f7b57aafce70f8adfd2abe03432c283320b

qpms/indexing.h

@@ -72,8 +72,8 @@ static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
 	return QPMS_SUCCESS;
 }
 
-/// Conversion from universal VSWF index u to type and the traditional layout index.
+/// Conversion from universal VSWF index u to type and the traditional layout index (\a l > 0).
-/** Does not for allow longitudinal waves. */
+/** Does *not* allow 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;
@@ -83,6 +83,18 @@ static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u,
 	return QPMS_SUCCESS;
 }
 
+/// Conversion from universal VSWF index u to type and the traditional (vector) layout index.
+/** *Does* allow for longitudinal waves, but returns an error if l == 0.
+ * (The only legit VSWF with l = 0 is the longitudinal one; u == 0.) 
+ */
+static inline qpms_errno_t qpms_uvswfi2ty_l(qpms_uvswfi_t u,
+		qpms_vswf_type_t *t, qpms_y_t *y) {
+	*t = u & 3;
+	*y = u / 4 - 1;
+	if (*t == 3) return QPMS_ERROR;
+	if (*y < 0) return QPMS_ERROR;
+	return QPMS_SUCCESS;
+}
 
 /// Extract degree \a m from an universal VSWF index \a u.
 static inline qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) {

qpms/qpms_error.h

@@ -45,7 +45,8 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
 
 #define QPMS_DEBUG(type, msg, ...) qpms_debug_at_flf(__FILE__,__LINE__,__func__,type,msg, ##__VA_ARGS__)
 
-#define QPMS_CRASHING_MALLOC(pointer, size) {(pointer) = malloc(size); if(!pointer && (size)) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Allocation of %zd bytes for " #pointer " failed.", (size_t) (size));}
+#define QPMS_CRASHING_MALLOC(pointer, size) {(pointer) = malloc(size); if(!(pointer) && (size)) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Allocation of %zd bytes for " #pointer " failed.", (size_t) (size));}
+#define QPMS_CRASHING_CALLOC(pointer, nmemb, size) {(pointer) = calloc((nmemb), (size)); if(!(pointer) && (nmemb) && (size)) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Allocation of %zd bytes for " #pointer " failed.", (size_t)((nmemb)*(size)));}
 
 #define QPMS_CRASHING_REALLOC(pointer, size) {(pointer) = realloc((pointer), size); if(!(pointer) && (size)) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Rellocation of %zd bytes for " #pointer " failed.", (size_t) (size));}
 
@@ -65,6 +66,8 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
 
 #define QPMS_ENSURE(x, msg, ...) {if(!(x)) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); }
 
+#define QPMS_ASSERT(x) {if(!(x)) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__, "Unexpected error. This is most certainly a bug.")
+
 #define QPMS_NOT_IMPLEMENTED(msg, ...) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__, \
 		"Not implemented:" msg, ##__VA_ARGS__)
 

qpms/vswf.c

@@ -150,9 +150,18 @@ void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *w) {
   free(w);
 }
 
-qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, csphvec_t * const mgtarget, csphvec_t * const eltarget,
+csphvec_t qpms_vswf_L00(sph_t kdrj, qpms_bessel_t btyp,
-    qpms_l_t lMax, sph_t kr,
+    qpms_normalisation_t norm) { 
-    qpms_bessel_t btyp, qpms_normalisation_t norm) {
+  QPMS_UNTESTED;
+  // CHECKME Is it OK to ignore norm?? (Is L_0^0 the same for all conventions?)
+  complex double bessel0;
+  QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, 0, kr.r, &bessel0));
+  return {0.25 * M_2_SQRTPI * bessel0, 0, 0};
+}
+
+qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, 
+    csphvec_t * const mgtarget, csphvec_t * const eltarget, qpms_l_t lMax,
+    sph_t kr, qpms_bessel_t btyp, qpms_normalisation_t norm) {
   assert(lMax >= 1);
   complex double *bessel = malloc((lMax+1)*sizeof(complex double));
   if(qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel)) abort();
@@ -173,7 +182,7 @@ qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, csphvec_t * const mgtar
     besderfac = *(pbes-1) - l * besfac;
     for(qpms_m_t m = -l; m <= l; ++m) {
       complex double eimf = cexp(m * kr.phi * I);
-      if (longtarget) {
+      if (longtarget) { QPMS_UNTESTED;
         complex double longfac = sqrt(l*(l+1)) * eimf;
         plong->rc = // FATAL FIXME: I get wrong result here for plane wave re-expansion 
           // whenever kr.r > 0 (for waves with longitudinal component, ofcoz)
@@ -406,12 +415,47 @@ csphvec_t qpms_eval_vswf(sph_t kr,
 }
 
 
-#if 0
 csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec,
-    const complex double *coeffs, sph_t evalpoint,
+    const complex double *coeffs, const sph_t kr,
-    qpms_bessel_t btyp) {
+    const qpms_bessel_t btyp) {
+  QPMS_UNTESTED;
   const qpms_l_t lMax = b->lMax;
-  double *M, *N, *L;
+  complex double *cM = NULL, *cN = NULL, *cL = NULL, cL00 = 0;
-
+  if (b->lMax_L > 0)
+    QPMS_CRASHING_CALLOC(cL, lMax, sizeof(complex double));
+  if (b->lMax_M > 0)
+    QPMS_CRASHING_CALLOC(cM, lMax, sizeof(complex double));
+  if (b->lMax_N > 0)
+    QPMS_CRASHING_CALLOC(cN, lMax, sizeof(complex double));
+  for (size_t i = 0; i < bspec->n; ++i) {
+    if (bspec->ilist[i] == 0) // L00, needs special care
+      cL00 = coeffs[i];
+    else {
+      qpms_vswf_type_t t;
+      qpms_y_t y;
+      QPMS_ENSURE_SUCCESS(qpms_uvswfi2ty_l(bspec->ilist[i], &t, &y));
+      switch(t) {
+        case QPMS_VSWF_LONGITUDINAL:
+          QPMS_ASSERT(cL);
+          cL[y] = coeffs[i];
+          break;
+        case QPMS_VSWF_MAGNETIC:
+          QPMS_ASSERT(cM);
+          cM[y] = coeffs[i];
+          break;
+        case QPMS_VSWF_ELECTRIC:
+          QPMS_ASSERT(cN);
+          cN[y] = coeffs[i];
+          break;
+        default:
+          QPMS_WTF;
+      }
+    }
+  }
+  csphvec_t result = qpms_eval_vswf(kr, cL, cM, cN, lMax, btyp, norm);
+  if(cL00)
+    result = csphvec_add(result,
+       scphvec_scale(cL00, qpms_vswf_L00(kr, btyp, norm)));
+  return result;
 }
-#endif 
+

qpms/vswf.h

@@ -51,8 +51,8 @@ qpms_errno_t qpms_uvswf_fill(
 		const qpms_vswf_set_spec_t *setspec,
 		sph_t evaluation_point, qpms_bessel_t btyp);
 
-/// NOT IMPLEMENTED Evaluates field specified by SVWF coefficients at a given point.
+/// Evaluates field specified by SVWF coefficients at a given point.
-/** SVWF coefficients in \a coeffs must be ordered according to \a setspec->ilist
+/** SVWF coefficients in \a coeffs must be ordered according to \a setspec->ilist.
  */
 csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec,
 		const complex double *coeffs, sph_t evaluation_point,
@@ -83,9 +83,29 @@ qpms_errno_t qpms_legendre_deriv_y_get(double **result, double **result_deriv, d
 qpms_errno_t qpms_legendre_deriv_y_fill(double *where, double *where_deriv, double x, 
 		qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase); 
 
-/* some of the result targets may be NULL */
+
-qpms_errno_t qpms_vswf_fill(csphvec_t *resultL, csphvec_t *resultM, csphvec_t *resultN, qpms_l_t lMax, sph_t kdrj,
+/// Evaluate the zeroth-degree longitudinal VSWF \f$ \mathbf{L}_0^0 \f$.
-		qpms_bessel_t btyp, qpms_normalisation_t norm);
+csphvec_t qpms_vswf_L00(
+		sph_t kdrj, //< VSWF evaluation point.
+		qpms_bessel_t btyp,
+		qpms_normalisation_t norm);
+
+/// Evaluate VSWFs at a given point from \a l = 1 up to a given degree \a lMax.
+/**
+ * The target arrays \a resultL, \a resultM, \a resultN have to be large enough to contain
+ * \a lMax * (\a lMax + 2) elements. If NULL is passed instead, the corresponding SVWF type
+ * is not evaluated.
+ * 
+ * Does not evaluate the zeroth-order wave \f$ \mathbf{L}_0^0 \f$. 
+ * If you need that, use qpms_vswf_L00().
+ */
+qpms_errno_t qpms_vswf_fill(
+	csphvec_t *resultL, //< Target array for longitudinal VSWFs.
+       	csphvec_t *resultM, //< Target array for magnetic VSWFs.
+       	csphvec_t *resultN, //< Target array for electric VSWFs.
+	qpms_l_t lMax, //< Maximum multipole degree to be calculated.
+	sph_t kdrj, //< VSWF evaluation point.
+	qpms_bessel_t btyp, qpms_normalisation_t norm);
 // Should give the same results: for consistency checks
 qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *resultL, csphvec_t *resultM, csphvec_t *resultN, qpms_l_t lMax, sph_t kdrj,
 		qpms_bessel_t btyp, qpms_normalisation_t norm);