From f458e94995d9d183da32f01877f54c8098739c49 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Fri, 2 Feb 2018 14:23:15 +0000
Subject: [PATCH] Plane wave decomposition compiles but yields wrong results.

Former-commit-id: 0cfe82f1efe0e68d41a2f4bf315dfac4618cf435
---
 qpms/bessel.c     |  2 +-
 qpms/qpms_types.h |  3 +++
 qpms/vectors.h    | 18 ++++++++++++++++--
 qpms/vswf.c       |  1 -
 4 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/qpms/bessel.c b/qpms/bessel.c
index 5035a56..b5488b3 100644
--- a/qpms/bessel.c
+++ b/qpms/bessel.c
@@ -1,7 +1,7 @@
+#include <assert.h>
 #include "translations.h"
 #include <stdlib.h>
 #include <gsl/gsl_sf_bessel.h>
-#include <assert.h>
 
 int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) {
 	int retval;
diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h
index bc03db0..e4b4712 100644
--- a/qpms/qpms_types.h
+++ b/qpms/qpms_types.h
@@ -7,6 +7,9 @@
 #ifndef M_PI_2
 #define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487)
 #endif
+#ifndef M_PI
+#define M_PI (3.14159265358979323846)
+#endif
 
 // integer index types
 typedef int qpms_lm_t;
diff --git a/qpms/vectors.h b/qpms/vectors.h
index c95285e..afc64d1 100644
--- a/qpms/vectors.h
+++ b/qpms/vectors.h
@@ -34,6 +34,16 @@ static inline cart3_t cart3_add(const cart3_t a, const cart3_t b) {
 	return res;
 }
 
+static inline cart3_t cart3_scale(const double c, const cart3_t v) {
+	cart3_t res = {c * v.x, c * v.y, c * v.z};
+	return res;
+}
+
+static inline ccart3_t ccart3_scale(const complex  double c, const ccart3_t v) {
+	ccart3_t res = {c * v.x, c * v.y, c * v.z};
+	return res;
+}
+
 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};
 	return res;
@@ -44,13 +54,17 @@ static inline csphvec_t csphvec_scale(complex double c, const csphvec_t v) {
 	return res;
 }
 
-static inline complex double csphvec_dotnc(csphvec_t a, csphvec_t b) {
+static inline complex double csphvec_dotnc(const csphvec_t a, const csphvec_t b) {
 	//N.B. no complex conjugation done here
 	return a.rc * b.rc + a.thetac * b.thetac + a.phic * b.phic;
 }
 
+static inline double cart3_dot(const cart3_t a, const cart3_t b) {
+	return a.x * b.x + a.y * b.y + a.z * b.z;
+}
+
 // equivalent to sph_loccart2cart in qpms_p.py
-static inline ccart3_t csphvec2cart(const csphvec_t sphvec, const sph_t at) {
+static inline ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at) {
 	const double st = sin(at.theta);
 	const double ct = cos(at.theta);
 	const double sf = sin(at.phi);
diff --git a/qpms/vswf.c b/qpms/vswf.c
index 5d21ee8..e6e367c 100644
--- a/qpms/vswf.c
+++ b/qpms/vswf.c
@@ -400,7 +400,6 @@ static inline complex double ipowl(qpms_l_t l) {
 qpms_errno_t qpms_planewave2vswf_fill_sph(sph_t wavedir, csphvec_t amplitude,
 		complex double *target_longcoeff, complex double *target_mgcoeff,
 	       	complex double *target_elcoeff, qpms_l_t lMax, qpms_normalisation_t norm) {
-	abort(); //NI
 	qpms_y_t nelem = qpms_lMax2nelem(lMax);
 	csphvec_t * const dual_A1 = malloc(3*nelem*sizeof(csphvec_t)), *const dual_A2 = dual_A1 + nelem,
 		  * const dual_A3 = dual_A2 + nelem;