From 4dc151d5b978ee0c0b1878386508dd5d5ad5853b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Sun, 4 Feb 2018 15:05:41 +0000
Subject: [PATCH] Correct plane wave decomposition at origin (wrong elsewhere
 for longitudinal waves

Former-commit-id: 4dc70fde23e31828fb959da429d716bd39e84e47
---
 qpms/vswf.c | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/qpms/vswf.c b/qpms/vswf.c
index 2690456..c91cb8c 100644
--- a/qpms/vswf.c
+++ b/qpms/vswf.c
@@ -230,13 +230,23 @@ qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, csphvec_t * const mgtar
 	double const *ptau = pt.tau;
 	csphvec_t *plong = longtarget, *pmg = mgtarget, *pel = eltarget;
 	for(qpms_l_t l = 1; l <= lMax; ++l) {
-		complex double besfac = *pbes / kr.r;
-		complex double besderfac = *(pbes-1) - l * besfac;
+		complex double besfac;
+		complex double besderfac;
+	       	if (kr.r) {
+			besfac = *pbes / kr.r;
+		} else {
+			besfac = (1 == l) ? 1/3. : 0;
+		}
+		besderfac = *(pbes-1) - l * besfac;
 		for(qpms_m_t m = -l; m <= l; ++m) {
 			complex double eimf = cexp(m * kr.phi * I);
 			if (longtarget) {
 				complex double longfac = sqrt(l*(l+1)) * eimf;
-				plong->rc = (besderfac-besfac) * (*pleg) * longfac;
+				plong->rc = // FATAL FIXME: I get wrong result here for plane wave re-expansion 
+					// whenever kr.r > 0 (for waves with longitudinal component, ofcoz)
+					/*(*(pbes-1) - (l+1)/kr.r* *pbes)*/
+					(besderfac-besfac) 
+					* (*pleg) * longfac;
 				plong->thetac = *ptau * besfac * longfac;
 				plong->phic = *ppi * I * besfac * longfac;
 				++plong;