Merge branch 'newnorm'

Former-commit-id: 7c1941cb4436ab1d7920c85df4a82ccba24cccb2

notes/Electrodynamics.bib

@@ -353,6 +353,7 @@ Cross-referenced as UMIACS-TR-2001-44},
   abstract = {This thesis addresses optical binding - a new area of interest within the field of optical micromanipulation. It presents, for the first time, a rigorous numerical simulation of some of the key results, along with new experimental findings and also physical interpretations of the results. In an optical trap particles are attracted close to areas of high optical intensities and intensity gradients. So, for example, if two lasers are pointed towards each other (a counter propagating trap) then a single particle is trapped in the centre of the two beams \textendash{} the system is analogous to a particle being held by two springs in a potential well. If one increases the number of particles in the trap then naively one would expect all the particles to collect in the centre of the well. However, the effect of optical binding means that the presence of one particle affects the distribution of light experienced by another particle, resulting in extremely complex interactions that can lead to unusual 1D and 2D structures to form within the trap. Optical binding is not only of theoretical interest but also has applications in micromanipulation and assembly.},
   language = {en},
   publisher = {{Springer Science \& Business Media}},
+  url = {http://gen.lib.rus.ec/book/index.php?md5=3B222192BB23C6D062BAABB76696D0E4},
   author = {Taylor, Jonathan M.},
   month = jul,
   year = {2011},
@@ -1817,6 +1818,7 @@ The anapole is an intriguing example of a nonradiating source useful in the stud
   abstract = {This book is an introduction to some of the most important properties of electromagnetic waves and their interaction with passive materials and scatterers. The main purpose of the book is to give a theoretical treatment of these scattering phenomena, and to illustrate numerical computations of some canonical scattering problems for different geometries and materials. The scattering theory is also important in the theory of passive antennas, and this book gives several examples on this topic. Topics covered include an introduction to the basic equations used in scattering; the Green functions and dyadics; integral representation of fields; introductory scattering theory; scattering in the time domain; approximations and applications; spherical vector waves; scattering by spherical objects; the null-field approach; and propagation in stratified media. The book is organised along two tracks, which can be studied separately or together. Track 1 material is appropriate for a first reading of the textbook, while Track 2 contains more advanced material suited for the second reading and for reference. Exercises are included for each chapter.},
   language = {English},
   publisher = {{Scitech Publishing}},
+  url = {http://gen.lib.rus.ec/book/index.php?md5=00CCB3E221E741ADDB2E236FD4A9F002},
   author = {Kristensson, Gerhard},
   month = jul,
   year = {2016},

notes/conventions.md

@@ -1,5 +1,5 @@
-VSWF conventions
+VSWF conventions {#vswf_conventions}
-================
+====================================
 
 In general, the (transversal) VSWFs can be defined using (some) vector spherical harmonics
 as follows: \f[
@@ -63,6 +63,135 @@ GSL computes \f$ \rawFer{l}{m} \f$ unless the corresponding `csphase` argument i
 but can be tested by running `gsl_sf_legendre_array_e` for some specific arguments and comparing signs.
 
 
+Convention effects on symmetry operators
+----------------------------------------
+
+### Spherical harmonics
+
+Let' have two different (complex) spherical harmonic conventions connected by constant factors:
+\f[
+	\spharm[a]{l}{m} = c^\mathrm{a}_{lm}\spharm{l}{m}.
+\f]
+
+Both sets can be used to describe an angular function \f$ f \f$
+\f[
+	f = \sum_{lm} f^\mathrm{a}_{lm} \spharm[a]{l}{m} 
+	  = \sum_{lm} f^\mathrm{a}_{lm} c^\mathrm{a}_{lm}\spharm{l}{m}
+	  = \sum_{lm} f_{lm} \spharm{l}{m}.
+\f]
+
+If we perform a (symmetry) transformation \f$ g \f$ acting on the \f$ \spharm{l}{m} \f$
+basis via matrix \f$ D(g)_{l,m;l',m'} \f$, i.e.
+\f[
+	g\pr{\spharm{l}{m}} = \sum_{l'm'} D(g)_{l,m;l'm'} \spharm{l'}{m'},
+\f]
+we see
+\f[
+	g(f) = \sum_{lm} f_{lm}\sum_{l'm'} D(g)_{l,m;l'm'} \spharm{l'}{m'} 
+	     = \sum_{lm} f^\mathrm{a}_{lm} c^\mathrm{a}_{lm}\sum_{l'm'} D(g)_{l,m;l'm'} \spharm{l'}{m'}.
+\f]
+Rewriting the transformation action in the second basis
+\f[
+	g\pr{\spharm[a]{l}{m}} = \sum_{l'm'} D(g)^\mathrm{a}_{l,m;l'm'} \spharm[a]{l'}{m'},\\
+	g(f) = \sum_{lm} f^\mathrm{a}_{lm}\sum_{l'm'} D(g)^\mathrm{a}_{l,m;l'm'} \spharm[a]{l'}{m'},
+\f]
+and performing some substitutions,
+\f[ 
+	g(f) = \sum_{lm} \frac{f_{lm}}{c^\mathrm{a}_{lm}}
+		\sum_{l'm'} D(g)^\mathrm{a}_{l,m;l'm'} c^\mathrm{a}_{l'm'}\spharm{l'}{m'},
+\f]
+and comparing, we get
+\f[
+	D(g)^\mathrm{a}_{l,m;l'm'} = \frac{c^\mathrm{a}_{lm}}{c^\mathrm{a}_{l'm'}}D(g)_{l,m;l'm'}.
+\f]
+
+If the  difference between conventions is in particular Condon-Shortley phase, 
+this means a \f$ (-1)^{m-m'} \f$ factor between the transformation matrices.
+This does not affect the matrices for the inversion and
+mirror symmetry operations with
+respect to the \a xy, \a yz and \a xz planes, because they are all diagonal
+or anti-diagonal with respect to \a m (hence \f$ m-m \f$ is either zero
+or anyways even integer).
+It does, however, affect rotations, flipping the sign of the rotations
+along the \a z axis.
+
+Apparently, a constant complex factor independent of \f$ l,m \f$
+does nothing to the form of the tranformation matrix.
+
+These conclusions about transformations of spherical harmonics
+hold also for the VSWFs built on top of them.
+
+
+
+Convention effect on translation operators
+------------------------------------------
+
+Let us declare VSWFs in Kristensson's conventions below, 
+\f$ \wfkc \f$ \cite kristensson_spherical_2014, 
+\f$ \wfkr \f$ \cite kristensson_scattering_2016, as the "canonical"
+spherical waves based on complex and real spherical harmonics, respectively.
+They both have the property that the translation operators \f$ \tropRrr{}{},\tropSrr{}{} \f$ 
+that transform
+the VSWF field expansion coefficients between different origins, e.g.
+\f[
+	\wfkcreg(\vect{r}) = \tropRrr{\vect r}{\vect r'} \wfkcreg(\vect{r'}),
+\f]
+actually consist of two different submatrices $A,B$ for the same-type and different-type
+(in the sense of "electric" versus "magnetic" waves) that repeat themselves once:
+\f[
+	\begin{bmatrix} \wfkcreg_1(\vect{r}) \\ \wfkcreg_2(\vect{r}) \end{bmatrix} 
+	= \begin{bmatrix} A & B \\ B & A \end{bmatrix}(\vect{r} \leftarrow \vect{r'})
+	\begin{bmatrix} \wfkcreg_1(\vect{r'}) \\ \wfkcreg_2(\vect{r'}) \end{bmatrix}.
+\f]
+(This symmetry holds also for singular translation operators \f$ \tropSrr{}{} \f$
+and real spherical harmonics based VSWFs \f$ \wfkr \f$.)
+
+However, the symmetry above will not hold like this in some stupider convention.
+Let's suppose that one uses a different convention with some additional coefficients
+compared to the canonical one,
+\f[ 
+	\wfm_{lm} = \alpha_{\wfm lm} \wfkc_{1lm},\\
+	\wfe_{lm} = \alpha_{\wfe lm} \wfkc_{2lm}.\\
+\f]
+and with field expansion (WLOG assume regular fields only)
+\f[ \vect E = c_{\wfe l m} \wfe_{lm} + c_{\wfm l m } \wfm_{lm}. \f]
+Under translations, the coefficients then transform like
+\f[
+	\begin{bmatrix} \alpha_\wfe(\vect{r}) \\ \alpha_\wfm(\vect{r}) \end{bmatrix} 
+	= \begin{bmatrix} R_{\wfe\wfe} & R_{\wfe\wfm} \\ 
+	 	          R_{\wfm\wfe} & R_{\wfm\wfm} 
+	  \end{bmatrix}(\vect{r} \leftarrow \vect{r'})
+	\begin{bmatrix} \alpha_\wfe(\vect{r'}) \\ \alpha_\wfm(\vect{r'}) \end{bmatrix},
+\f]
+and by substituting and comparing the expressions for canonical waves above, one gets
+\f[
+	R_{\wfe,lm;\wfe,l'm'} = \alpha_{\wfe lm}^{-1} A_{lm,l'm'} \alpha_{\wfe l'm'},\\
+	R_{\wfe,lm;\wfm,l'm'} = \alpha_{\wfe lm}^{-1} B_{lm,l'm'} \alpha_{\wfm l'm'},\\
+	R_{\wfm,lm;\wfe,l'm'} = \alpha_{\wfm lm}^{-1} B_{lm,l'm'} \alpha_{\wfe l'm'},\\
+	R_{\wfm,lm;\wfm,l'm'} = \alpha_{\wfm lm}^{-1} A_{lm,l'm'} \alpha_{\wfm l'm'}.
+\f]
+	
+If the coefficients for magnetic and electric waves are the same,
+\f$ \alpha_{\wfm lm} = \alpha_{\wfe lm} \f$, the translation operator 
+can be written in the same symmetric form as with the canonical convention,
+just the matrices \f$ A, B\f$ will be different inside. 
+
+If the coefficients differ (as in SCUFF-EM convention, where there
+is a relative \a i -factor between electric and magnetic waves),
+the functions such as qpms_trans_calculator_get_AB_arrays() will 
+compute \f$ R_{\wfe\wfe}, R_{\wfe\wfm} \f$ for \f$ A, B \f$ arrays.
+The remaining matrices' elements must then be obtained as
+\f[
+	R_{\wfm,lm;\wfe,l'm'} = \alpha_{\wfm lm}^{-1} \alpha_{\wfe lm} 
+		R_{\wfe,lm;\wfm,l'm'} \alpha_{\wfm l'm'}^{-1} \alpha_{\wfe l'm'}
+		= g_{lm}R_{\wfe,lm;\wfm,l'm'}g_{l'm'}, \\
+	R_{\wfm,lm;\wfm,l'm'} = \alpha_{\wfm lm}^{-1} \alpha_{\wfe lm} 
+		R_{\wfe,lm;\wfe,l'm'} \alpha_{\wfe l'm'}^{-1} \alpha_{\wfm l'm'}
+		= g_{lm}R_{\wfe,lm;\wfe,l'm'}g_{l'm'}^{-1},
+\f]
+where the coefficients \f$ g_{lm} \f$ can be obtained by
+qpms_normalisation_factor_N_M().
+
 Literature convention tables
 ----------------------------
 
@@ -141,7 +270,7 @@ Literature convention tables
 	\vect E = \sum_\alpha \pr{ \wcrreg_\alpha \wfrreg_\alpha + \wcrout_\alpha \wfrout_\alpha }, \\
 	\vect H = \frac{1}{Z_0Z^r} \sum_\alpha \pr{ \wcrreg_\alpha \sigma_\alpha\wfrreg_\overline{\alpha} +
 		 \wcrout_\alpha \sigma_\alpha\wfrout_\overline{\alpha}},
-\f] where \f$ \sigma_{lmM} = +1, \sigma_{lmN}=-1, \overline{lmM}=lmM, \overline{lmN}=lmM, \f$  cf. eq. (6). The notation is not extremely consistent throughout Reid's memo.	| 	| 	|
+\f] where \f$ \sigma_{lmM} = +1, \sigma_{lmN}=-1, \overline{lmM}=lmN, \overline{lmN}=lmM, \f$  cf. eq. (6). The notation is not extremely consistent throughout Reid's memo.	| 	| 	|
 | Taylor \cite taylor_optical_2011	| \f[
 	\wfet_{mn}^{(j)}	=	\frac{n(n+1)}{kr}\sqrt{\frac{2n+1}{4\pi}\frac{\left(n-m\right)!}{\left(n+m\right)!}}\Fer[Taylor]{n}{m}\left(\cos\theta\right)e^{im\phi}z_{n}^{j}\left(kr\right)\uvec{r} \\
 		+\left[\tilde{\tau}_{mn}\left(\cos\theta\right)\uvec{\theta}+i\tilde{\pi}_{mn}\left(\cos\theta\right)\uvec{\phi}\right]e^{im\phi}\frac{1}{kr}\frac{\ud\left(kr\,z_{n}^{j}\left(kr\right)\right)}{\ud(kr)}, \\ 

notes/mathjax_newcommands.js

@@ -31,6 +31,8 @@ MathJax.Hub.Config({
 	    spharm: ["{{Y_{\\mathrm{#1}}}_{#2}^{#3}}", 3, ""], // Spherical harmonics
 	    spharmR: ["{{Y_{\\mathrm{#1}}}_{\\mathrm{#1}{#2}{#3}}", 4, ""], // Spherical harmonics
 	    csphase: "\\mathsf{C_{CS}}", // Condon-Shortley phase
+	    tropSrr: ["{{S^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator singular
+	    tropRrr: ["{{R^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator regular
 
             // Kristensson's VSWFs, complex version (2014 notes)
 	    wfkc: "{\\vect{y}}", // any wave

qpms/CMakeLists.txt

@@ -12,7 +12,7 @@ include_directories(${DIRS})
 
 add_library (qpms translations.c tmatrices.c vecprint.c vswf.c wigner.c
 	lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
-	bessel.c own_zgemm.c parsing.c)
+	bessel.c own_zgemm.c parsing.c scatsystem.c)
 use_c99()
 
 set(LIBS ${LIBS} ${GSL_LIBRARIES} ${GSLCBLAS_LIBRARIES})

qpms/bessel.c

@@ -8,6 +8,7 @@
 #include <complex.h>
 #include "qpms_error.h"
 #include <amos.h>
+#include <math.h>
 
 #ifndef M_LN2
 #define M_LN2   0.69314718055994530942  /* log_e 2 */

qpms/groups.h

@@ -25,6 +25,7 @@
 #define QPMS_GROUPS_H
 
 #include "qpms_types.h"
+#include <assert.h>
 
 /// To be used only in qpms_finite_group_t
 struct qpms_finite_group_irrep_t {

qpms/legendre.c

@@ -5,14 +5,16 @@
 #include <stdlib.h>
 #include "indexing.h"
 #include <string.h>
+#include "qpms_error.h"
 
 // Legendre functions also for negative m, see DLMF 14.9.3
 qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, double x, qpms_l_t lMax,
     gsl_sf_legendre_t lnorm, double csphase)
 {
-  size_t n = gsl_sf_legendre_array_n(lMax);
+  const size_t n = gsl_sf_legendre_array_n(lMax);
-  double *legendre_tmp = malloc(n * sizeof(double));
+  double *legendre_tmp, *legendre_deriv_tmp;
-  double *legendre_deriv_tmp = malloc(n * sizeof(double));
+  QPMS_CRASHING_MALLOC(legendre_tmp, n * sizeof(double));
+  QPMS_CRASHING_MALLOC(legendre_deriv_tmp, n * sizeof(double));
   int gsl_errno = gsl_sf_legendre_deriv_array_e(
       lnorm, (size_t)lMax, x, csphase, legendre_tmp,legendre_deriv_tmp);
   for (qpms_l_t l = 1; l <= lMax; ++l)
@@ -22,153 +24,69 @@ qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, do
       target[y] = legendre_tmp[i];
       target_deriv[y] = legendre_deriv_tmp[i];
     }
-  switch(lnorm) {
+
-    case GSL_SF_LEGENDRE_NONE:
+  // Fill negative m's.
-      for (qpms_l_t l = 1; l <= lMax; ++l)
+  for (qpms_l_t l = 1; l <= lMax; ++l)
-        for (qpms_m_t m = 1; m <= l; ++m) {
+    for (qpms_m_t m = 1; m <= l; ++m) {
-          qpms_y_t y = qpms_mn2y(-m,l);
+      qpms_y_t y = qpms_mn2y(-m,l);
-          size_t i = gsl_sf_legendre_array_index(l,m);
+      size_t i = gsl_sf_legendre_array_index(l,m);
-          // viz DLMF 14.9.3, čert ví, jak je to s cs fasí.
+      // cf. DLMF 14.9.3, but we're normalised.
-          double factor = exp(lgamma(l-m+1)-lgamma(l+m+1))*((m%2)?-1:1);
+      double factor = ((m%2)?-1:1); 
-          target[y] = factor * legendre_tmp[i];
+      target[y] = factor * legendre_tmp[i];
-          target_deriv[y] = factor * legendre_deriv_tmp[i];
+      target_deriv[y] = factor * legendre_deriv_tmp[i];
-        }
+    }
-      break;
+
-    case GSL_SF_LEGENDRE_SCHMIDT:
-    case GSL_SF_LEGENDRE_SPHARM:
-    case GSL_SF_LEGENDRE_FULL:
-      for (qpms_l_t l = 1; l <= lMax; ++l)
-        for (qpms_m_t m = 1; m <= l; ++m) {
-          qpms_y_t y = qpms_mn2y(-m,l);
-          size_t i = gsl_sf_legendre_array_index(l,m);
-          // viz DLMF 14.9.3, čert ví, jak je to s cs fasí.
-          double factor = ((m%2)?-1:1); // this is the difference from the unnormalised case
-          target[y] = factor * legendre_tmp[i];
-          target_deriv[y] = factor * legendre_deriv_tmp[i];
-        }
-      break;
-    default:
-      abort(); //NI
-      break;
-  }
   free(legendre_tmp);
   free(legendre_deriv_tmp);
-  return QPMS_SUCCESS;
+  return gsl_errno;
 }
 
 qpms_errno_t qpms_legendre_deriv_y_get(double **target, double **dtarget, double x, qpms_l_t lMax, gsl_sf_legendre_t lnorm,
     double csphase)
 {
 
-  *target = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
+  const qpms_y_t nelem = qpms_lMax2nelem(lMax);
-  *dtarget = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
+  QPMS_CRASHING_MALLOC(target, nelem * sizeof(double));
+  QPMS_CRASHING_MALLOC(dtarget, nelem * sizeof(double));
   return qpms_legendre_deriv_y_fill(*target, *dtarget, x, lMax, lnorm, csphase);
 }
 
-qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm)
+qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, const double csphase)
 {
-  const double csphase = qpms_normalisation_t_csphase(norm);
+  QPMS_ENSURE(fabs(csphase) == 1, "The csphase argument must be either 1 or -1, not %g.", csphase);
-  norm = qpms_normalisation_t_normonly(norm);
   qpms_pitau_t res;
-  qpms_y_t nelem = qpms_lMax2nelem(lMax);
+  const qpms_y_t nelem = qpms_lMax2nelem(lMax);
-  res.pi = malloc(nelem * sizeof(double));
+  QPMS_CRASHING_MALLOC(res.leg, nelem * sizeof(double));
-  res.tau = malloc(nelem * sizeof(double));
+  QPMS_CRASHING_MALLOC(res.pi, nelem * sizeof(double));
+  QPMS_CRASHING_MALLOC(res.tau, nelem * sizeof(double));
   double ct = cos(theta), st = sin(theta);
   if (1 == fabs(ct)) { // singular case, use DLMF 14.8.2
-    memset(res.pi, 0, nelem*sizeof(double));
+    memset(res.pi, 0, nelem * sizeof(double));
-    memset(res.tau, 0, nelem*sizeof(double));
+    memset(res.tau, 0, nelem * sizeof(double));
-    res.leg = calloc(nelem, sizeof(double));
+    memset(res.leg, 0, nelem * sizeof(double));
-    switch(norm) {
+    for (qpms_l_t l = 1; l <= lMax; ++l) {
-      /* FIXME get rid of multiplicating the five lines */
+      res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)/(4*M_PI *l*(l+1)));
-      case QPMS_NORMALISATION_NONE:
+      double fl = 0.25 * sqrt((2*l+1)*M_1_PI);
-        for (qpms_l_t l = 1; l <= lMax; ++l) {
+      int lpar = (l%2)?-1:1;
-          res.leg[qpms_mn2y(0, l)] = (l%2)?ct:1.;
+      res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
-          double p = l*(l+1)/2;
+      res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
-          const double n = 0.5;
+      res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase;
-          int lpar = (l%2)?-1:1;
+      res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase;
-          res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * p * csphase;
-          res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * n * csphase;
-          res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * p * csphase;
-          res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * n * csphase;
-        }
-        break;
-#ifdef USE_XU_ANTINORMALISATION // broken
-      case QPMS_NORMALISATION_XU: // Rather useless except for testing.
-        for (qpms_l_t l = 1; l <= lMax; ++l) {
-          res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt(0.25*M_1_PI *(2*l+1)/(l*(l+1)));
-          double p = l*(l+1)/2 /* as in _NONE */ 
-            * sqrt(0.25 * M_1_PI * (2*l + 1));
-          double n = 0.5 /* as in _NONE */
-            * sqrt(0.25 * M_1_PI * (2*l + 1)) / (l * (l+1));
-          int lpar = (l%2)?-1:1;
-          res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * p * csphase;
-          res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * n * csphase;
-          res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * p * csphase;
-          res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * n * csphase;
-        }
-        break;
-#endif
-      case QPMS_NORMALISATION_TAYLOR:
-        for (qpms_l_t l = 1; l <= lMax; ++l) {
-          res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)*0.25*M_1_PI);
-          double fl = 0.25 * sqrt((2*l+1)*l*(l+1)*M_1_PI);
-          int lpar = (l%2)?-1:1;
-          res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
-          res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
-          res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase;
-          res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase;
-        }
-        break;
-      case QPMS_NORMALISATION_POWER:
-        for (qpms_l_t l = 1; l <= lMax; ++l) {
-          res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)/(4*M_PI *l*(l+1)));
-          double fl = 0.25 * sqrt((2*l+1)*M_1_PI);
-          int lpar = (l%2)?-1:1;
-          res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
-          res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase;
-          res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase;
-          res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase;
-        }
-        break;
-      default:
-        abort();
     }
   }
   else { // cos(theta) in (-1,1), use normal calculation
-    double *legder = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
+    double *legder;
-    res.leg = malloc(sizeof(double)*qpms_lMax2nelem(lMax));
+    QPMS_CRASHING_MALLOC(legder, nelem * sizeof(double));
-    if (qpms_legendre_deriv_y_fill(res.leg, legder, ct, lMax,
+    QPMS_ENSURE_SUCCESS(qpms_legendre_deriv_y_fill(res.leg, legder, ct, lMax,
-          norm == QPMS_NORMALISATION_NONE ? GSL_SF_LEGENDRE_NONE
+          GSL_SF_LEGENDRE_SPHARM, csphase));
-          : GSL_SF_LEGENDRE_SPHARM, csphase))
+    // Multiply by the "power normalisation" factor
-      abort();
+    for (qpms_l_t l = 1; l <= lMax; ++l) {
-    if (norm == QPMS_NORMALISATION_POWER)
+      double prefac = 1./sqrt(l*(l+1));
-      /* for None (=non-normalized) and Taylor (=sph. harm. normalized)
+      for (qpms_m_t m = -l; m <= l; ++m) {
-       * the correct normalisation is already obtained from gsl_sf_legendre_deriv_array_e().
+        res.leg[qpms_mn2y(m,l)] *= prefac;
-       * However, Kristensson ("power") normalisation differs from Taylor
+        legder[qpms_mn2y(m,l)] *= prefac;
-       * by 1/sqrt(l*(l+1)) factor.
-       */
-      for (qpms_l_t l = 1; l <= lMax; ++l) {
-        double prefac = 1./sqrt(l*(l+1));
-        for (qpms_m_t m = -l; m <= l; ++m) {
-          res.leg[qpms_mn2y(m,l)] *= prefac;
-          legder[qpms_mn2y(m,l)] *= prefac;
-        }
       }
-#ifdef USE_XU_ANTINORMALISATION
+    }
-    else if (norm == QPMS_NORMALISATION_XU)
-      /* for Xu (anti-normalized), we start from spharm-normalized Legendre functions
-       * Do not use this normalisation except for testing
-       */
-      // FIXME PROBABLY BROKEN HERE
-      for (qpms_l_t l = 1; l <= lMax; ++l) {
-        double prefac = (2*l + 1) / sqrt(4*M_PI / (l*(l+1)));
-        for (qpms_m_t m = -l; m <= l; ++m) {
-          double fac = prefac * exp(lgamma(l+m+1) - lgamma(l-m+1));
-          res.leg[qpms_mn2y(m,l)] *= fac;
-          legder[qpms_mn2y(m,l)] *= fac;
-        }
-      }
-#endif
     for (qpms_l_t l = 1; l <= lMax; ++l) {
       for (qpms_m_t m = -l; m <= l; ++m) {
         res.pi [qpms_mn2y(m,l)] = m / st * res.leg[qpms_mn2y(m,l)];

qpms/normalisation.h

@@ -0,0 +1,283 @@
+/*! \file normalisation.h
+ * \brief Convention-dependent coefficients for VSWFs.
+ *
+ * See also @ref qpms_normalisation_t and @ref vswf_conventions.
+ */
+#ifndef NORMALISATION_H
+#define NORMALISATION_H
+
+#include "qpms_types.h"
+#include "qpms_error.h"
+#include <math.h>
+#include <complex.h>
+#include "indexing.h"
+
+/// Returns the (real positive) common norm factor of a given normalisation compared to the reference convention.
+/** Does NOT perform the inversion if QPMS_NORMALISATION_INVERSE is set. */
+static inline double qpms_normalisation_normfactor(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
+	switch (norm & QPMS_NORMALISATION_NORM_BITS) {
+		case QPMS_NORMALISATION_NORM_POWER:
+			return 1;
+		case QPMS_NORMALISATION_NORM_SPHARM:
+			return sqrt(l*(l+1));
+		case QPMS_NORMALISATION_NORM_NONE: // TODO more precision
+			return sqrt(l*(l+1) * 4*M_PI / (2*l+1)) *
+			       	exp(0.5*(lgamma(l+m+1) - lgamma(l-m+1)));
+		default:
+			QPMS_WTF;
+	}
+}
+
+
+
+/// Returns the factors of a magnetic basis VSWF of a given convention compared to the reference convention.
+/**
+ * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
+ * phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
+ */
+static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
+	complex double fac = qpms_normalisation_normfactor(norm, l, m);
+	if (norm & QPMS_NORMALISATION_M_MINUS) fac *= -1;
+	if (norm & QPMS_NORMALISATION_M_I) fac *= I;
+	if (norm & QPMS_NORMALISATION_INVERSE) fac = 1/fac;
+	return fac;
+}
+
+
+/// Returns the factors of a magnetic basis VSWF of a given convention compared to the reference convention.
+/**
+ * This version takes into account the Condon-Shortley phase bit. 
+ * Do not use if the C.-S. has already been taken into account e.g. in
+ * a `gsl_sf_legendre_*_e()` call.
+ */
+static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
+	complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m);
+	return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
+}
+
+
+/// Returns the factors of a electric basis VSWF of a given convention compared to the reference convention.
+/**
+ * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
+ * phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
+ */
+static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
+	complex double fac = qpms_normalisation_normfactor(norm, l, m);
+	if (norm & QPMS_NORMALISATION_N_MINUS) fac *= -1;
+	if (norm & QPMS_NORMALISATION_N_I) fac *= I;
+	if (norm & QPMS_NORMALISATION_INVERSE) fac = 1/fac;
+	return fac;
+}
+
+
+/// Returns the factors of a electric basis VSWF of a given convention compared to the reference convention.
+/**
+ * This version takes into account the Condon-Shortley phase bit. 
+ * Do not use if the C.-S. has already been taken into account e.g. in
+ * a `gsl_sf_legendre_*_e()` call.
+ */
+static inline complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
+	complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m);
+	return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
+}
+
+
+/// Returns the factors of a electric basis VSWF divided by the factor of a magnetic VWFS of a given convention, compared to the reference one.
+static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
+	return qpms_normalisation_factor_N_noCS(norm, l, m) 
+		/ qpms_normalisation_factor_M_noCS(norm, l, m);
+}
+
+
+/// Returns the factors of a longitudinal basis VSWF of a given convention compared to the reference convention.
+/**
+ * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
+ * phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
+ */
+static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
+	complex double fac = qpms_normalisation_normfactor(norm, l, m);
+	if (norm & QPMS_NORMALISATION_L_MINUS) fac *= -1;
+	if (norm & QPMS_NORMALISATION_L_I) fac *= I;
+	if (norm & QPMS_NORMALISATION_INVERSE) fac = 1/fac;
+	return fac;
+}
+
+/// Returns the factors of a longitudinal basis VSWF of a given convention compared to the reference convention.
+/**
+ * This version takes into account the Condon-Shortley phase bit. 
+ * Do not use if the C.-S. has already been taken into account e.g. in
+ * a `gsl_sf_legendre_*_e()` call.
+ */
+static inline complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
+	complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m);
+	return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
+}
+
+/// Returns the factors of a basis VSWF of a given convention compared to the reference convention.
+static inline complex double qpms_normalisation_factor_uvswfi(const qpms_normalisation_t norm, qpms_uvswfi_t ui) {
+	qpms_vswf_type_t t; qpms_m_t m; qpms_l_t l;
+	qpms_uvswfi2tmn(ui, &t, &m, &l);
+	switch(t) {
+		case QPMS_VSWF_MAGNETIC:
+			return qpms_normalisation_factor_M(norm, l, m);
+		case QPMS_VSWF_ELECTRIC:
+			return qpms_normalisation_factor_N(norm, l, m);
+		case QPMS_VSWF_LONGITUDINAL:
+			return qpms_normalisation_factor_L(norm, l, m);
+		default:
+			QPMS_WTF;
+	}
+}
+
+
+/// Returns normalisation flags corresponding to the dual spherical harmonics / waves.
+/**
+ * This reverses the normalisation factors returned by qpms_normalisation_factor_*
+ * and conjugates the asimuthal part for complex spherical harmonics, 
+ * \f$ e^{\pm im\phi} \leftrightarrow e^{\mp im\phi} \f$.
+ */
+static inline qpms_normalisation_t qpms_normalisation_dual(qpms_normalisation_t norm) {
+	norm ^= QPMS_NORMALISATION_INVERSE;
+	if (!(norm & QPMS_NORMALISATION_SPHARM_REAL))
+		norm ^= QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE;
+	return norm;
+}
+
+/// Returns the asimuthal part of a spherical harmonic.
+/** Returns \f[ e^{im\phi} \f] for standard complex spherical harmonics,
+ * \f[ e^{-im\phi \f] for complex spherical harmonics 
+ * and QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE set.
+ *
+ * For real spherical harmonics, this gives
+ * \f[ 
+ * 	\sqrt{2}\cos{m \phi} \quad \mbox{if } m>0, \\
+ * 	\sqrt{2}\sin{m \phi} \quad \mbox{if } m<0, \\
+ * 	0 \quad \mbox{if } m>0. \\
+ * \f]
+ */
+static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) {
+	switch(norm & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
+		case 0:
+			return cexp(I*m*phi);
+		case QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE:
+			return cexp(-I*m*phi);
+		case QPMS_NORMALISATION_SPHARM_REAL:
+			if (m > 0) return M_SQRT2 * cos(m*phi);
+			else if (m < 0) return M_SQRT2 * sin(m*phi);
+			else return 1.;
+		default:
+			QPMS_WTF;
+	}
+}
+
+/// Returns derivative of the asimuthal part of a spherical harmonic divided by \a m.
+/**
+ *
+ * This is used to evaluate the VSWFs together with the \a pi member array of the
+ * qpms_pitau_t structure.
+ *
+ * Returns \f[ i e^{im\phi} \f] for standard complex spherical harmonics,
+ * \f[-i e^{-i\phi \f] for complex spherical harmonics 
+ * and QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE set.
+ *
+ * For real spherical harmonics, this gives
+ * \f[ 
+ * 	-\sqrt{2}\sin{m \phi} \quad \mbox{if } m>0, \\
+ * 	\sqrt{2}\cos{m \phi} \quad \mbox{if } m<0, \\
+ * 	-1 \quad \mbox{if } \mbox{if }m=0. \\
+ * \f]
+ *
+ * (The value returned for \f$ m = 0 \f$ should not actually be used for
+ * anything except for multiplying by zero.)
+ *
+ *
+ */
+static inline complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) {
+	if(m==0) return 0;
+	switch(norm & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
+		case 0:
+			return I*cexp(I*m*phi);
+		case QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE:
+			return -I*cexp(-I*m*phi);
+		case QPMS_NORMALISATION_SPHARM_REAL:
+			if (m > 0) return -M_SQRT2 * sin(m*phi);
+			else if (m < 0) return M_SQRT2 * cos(m*phi);
+			else return -1;
+		default:
+			QPMS_WTF;
+	}
+}
+
+#if 0 // legacy code moved from qpms_types.h. TODO cleanup
+/// Returns the normalisation convention code without the Condon-Shortley phase.
+static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) {
+	return norm & (~QPMS_NORMALISATION_T_CSBIT);
+}
+
+/* Normalisation of the spherical waves is now scattered in at least three different files:
+ * here, we have the norm in terms of radiated power of outgoing wave.
+ * In file legendre.c, function qpms_pitau_get determines the norm used in the vswf.c
+ * spherical vector wave norms. The "dual" waves in vswf.c use the ..._abssquare function below.
+ * In file translations.c, the normalisations are again set by hand using the normfac and lognormfac
+ * functions.
+ */
+#include <math.h>
+#include <assert.h>
+// relative to QPMS_NORMALISATION_KRISTENSSON_CS, i.e.
+// P_l^m[normtype] = P_l^m[Kristensson]
+static inline double qpms_normalisation_t_factor(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
+	int csphase = qpms_normalisation_t_csphase(norm);
+	norm = qpms_normalisation_t_normonly(norm);
+	double factor;
+	switch (norm) {
+		case QPMS_NORMALISATION_KRISTENSSON:
+			factor = 1.;
+			break;
+		case QPMS_NORMALISATION_TAYLOR:
+			factor = sqrt(l*(l+1));
+			break;
+		case QPMS_NORMALISATION_NONE:
+			factor = sqrt(l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)));
+			break;
+#ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c
+		case QPMS_NORMALISATION_XU:
+			factor = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
+			break;
+#endif
+		default:
+			assert(0);
+	}
+	factor *= (m%2)?(-csphase):1;
+	return factor;
+}
+
+
+// TODO move elsewhere
+static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
+	norm = qpms_normalisation_t_normonly(norm);
+	switch (norm) {
+		case QPMS_NORMALISATION_KRISTENSSON:
+			return 1.;
+			break;
+		case QPMS_NORMALISATION_TAYLOR:
+			return l*(l+1);
+			break;
+		case QPMS_NORMALISATION_NONE:
+			return l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
+			break;
+#ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c
+		case QPMS_NORMALISATION_XU:
+			{
+			  double fac = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
+			  return fac * fac;
+			}
+			break;
+#endif 
+		default:
+			assert(0);
+			return NAN;
+	}
+}
+#endif
+
+#endif //NORMALISATION_H

qpms/qpms_c.pyx

@@ -1,3 +1,8 @@
+"""@package qpms_c
+Cythonized parts of QPMS; mostly wrappers over the C data structures
+to make them available in Python.
+"""
+
 # Cythonized parts of QPMS here
 # -----------------------------
 
@@ -20,19 +25,15 @@ class VSWFType(enum.IntEnum):
     L = QPMS_VSWF_LONGITUDINAL
 
 class VSWFNorm(enum.IntEnum):
-    #XU = QPMS_NORMALISATION_XU
+    # TODO try to make this an enum.IntFlag if supported
-    #XU_CS = QPMS_NORMALISATION_XU_CS
+    # TODO add the other flags from qpms_normalisation_t as well
-    NONE = QPMS_NORMALISATION_NONE
+    UNNORM = QPMS_NORMALISATION_NORM_NONE
-    NONE_CS = QPMS_NORMALISATION_NONE_CS
+    UNNORM_CS = QPMS_NORMALISATION_NORM_NONE | QPMS_NORMALISATION_CSPHASE
-    POWER = QPMS_NORMALISATION_POWER
+    POWERNORM = QPMS_NORMALISATION_NORM_POWER
-    POWER_CS = QPMS_NORMALISATION_POWER_CS
+    POWERNORM_CS = QPMS_NORMALISATION_NORM_POWER | QPMS_NORMALISATION_CSPHASE
-    SPHARM = QPMS_NORMALISATION_SPHARM
+    SPHARMNORM = QPMS_NORMALISATION_NORM_SPHARM
-    SPHARM_CS = QPMS_NORMALISATION_SPHARM_CS
+    SPHARMNORM_CS = QPMS_NORMALISATION_NORM_SPHARM | QPMS_NORMALISATION_CSPHASE
     UNDEF = QPMS_NORMALISATION_UNDEF
-    KRISTENSSON = QPMS_NORMALISATION_KRISTENSSON
-    KRISTENSSON_CS = QPMS_NORMALISATION_KRISTENSSON_CS
-    TAYLOR = QPMS_NORMALISATION_TAYLOR
-    TAYLOR_CS = QPMS_NORMALISATION_TAYLOR_CS
 
 try:
     class DebugFlags(enum.IntFlag): # Should be IntFlag if python version >= 3.6
@@ -721,7 +722,7 @@ cdef class BaseSpec:
         if 'norm' in kwargs.keys():
             self.s.norm = kwargs['norm']
         else:
-            self.s.norm = QPMS_NORMALISATION_POWER_CS
+            self.s.norm = <qpms_normalisation_t>(QPMS_NORMALISATION_NORM_POWER | QPMS_NORMALISATION_CSPHASE)
         # set the other metadata
         cdef qpms_l_t l
         self.s.lMax_L = -1

qpms/qpms_cdefs.pxd

@@ -25,19 +25,25 @@ cdef extern from "qpms_types.h":
         cart2_t cart2
         pol_t pol
     ctypedef enum qpms_normalisation_t:
-        QPMS_NORMALISATION_XU
-        QPMS_NORMALISATION_XU_CS
-        QPMS_NORMALISATION_NONE
-        QPMS_NORMALISATION_NONE_CS
-        QPMS_NORMALISATION_KRISTENSSON
-        QPMS_NORMALISATION_KRISTENSSON_CS
-        QPMS_NORMALISATION_POWER
-        QPMS_NORMALISATION_POWER_CS
-        QPMS_NORMALISATION_TAYLOR
-        QPMS_NORMALISATION_TAYLOR_CS
-        QPMS_NORMALISATION_SPHARM
-        QPMS_NORMALISATION_SPHARM_CS
         QPMS_NORMALISATION_UNDEF
+        QPMS_NORMALISATION_INVERSE
+        QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE
+        QPMS_NORMALISATION_SPHARM_REAL
+        QPMS_NORMALISATION_CSPHASE
+        QPMS_NORMALISATION_M_I
+        QPMS_NORMALISATION_M_MINUS
+        QPMS_NORMALISATION_N_I
+        QPMS_NORMALISATION_N_MINUS
+        QPMS_NORMALISATION_L_I
+        QPMS_NORMALISATION_L_MINUS
+        QPMS_NORMALISATION_NORM_BITSTART
+        QPMS_NORMALISATION_NORM_POWER
+        QPMS_NORMALISATION_NORM_SPHARM
+        QPMS_NORMALISATION_NORM_NONE
+        QPMS_NORMALISATION_NORM_BITS
+        QPMS_NORMALISATION_CONVENTION_KRISTENSSON_REAL
+        QPMS_NORMALISATION_CONVENTION_KRISTENSSON
+        QPMS_NORMALISATION_CONVENTION_SCUFF
     ctypedef enum qpms_bessel_t:
         QPMS_BESSEL_REGULAR
         QPMS_BESSEL_SINGULAR

qpms/qpms_specfunc.h

@@ -1,3 +1,6 @@
+/*! \file qpms_specfunc.h
+ * \brief Various special and auxillary functions.
+ */
 #ifndef QPMS_SPECFUNC_H
 #define QPMS_SPECFUNC_H
 #include "qpms_types.h"
@@ -52,17 +55,40 @@ double *qpms_legendre_minus1d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); /
 
 
 
-// array of Legendre and pi, tau auxillary functions (see [1,(37)])
+/// Array of Legendre and and auxillary \f$\pi_{lm}, \tau_{lm} \f$ functions.
-// This should handle correct evaluation for theta -> 0 and theta -> pi
+/**
+ * The leg, pi, tau arrays are indexed using the standard qpms_mn2y() VSWF indexing.
+ */
 typedef struct {
 	//qpms_normalisation_t norm;
 	qpms_l_t lMax;
 	//qpms_y_t nelem;
 	double *leg, *pi, *tau;
 } qpms_pitau_t;
-qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm);
+
-void qpms_pitau_free(qpms_pitau_t);//NI
+/// Returns an array of normalised Legendre and and auxillary \f$\pi_{lm}, \tau_{lm} \f$ functions.
-void qpms_pitau_pfree(qpms_pitau_t*);//NI
+/**
+ * The normalised Legendre function here is defined as
+ * \f[ 
+ * 	\Fer[norm.]{l}{m} = \csphase^{-1} 
+ * 		\sqrt{\frac{1}{l(l+1)}\frac{(l-m)!(2l+1)}{4\pi(l+m)!}},
+ * \f] i.e. obtained using `gsl_sf_legendre_array_e()` with 
+ * `norm = GSL_SF_LEGENDRE_SPHARM` and multiplied by \f$ \sqrt{l(l+1)} \f$.
+ *
+ * The auxillary functions are defined as
+ * \f[
+ * 	\pi_{lm}(\cos \theta) = \frac{m}{\sin \theta} \Fer[norm.]{l}{m}(\cos\theta),\\
+ * 	\tau_{lm}(\cos \theta) = \frac{\ud}{\ud \theta} \Fer[norm.]{l}{m}(\cos\theta)
+ * \f]
+ * with appropriate limit expression used if \f$ \abs{\cos\theta} = 1 \f$.
+ *
+ * When done, don't forget to deallocate the memory using qpms_pitau_free().
+ *
+ */
+qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, double csphase);
+/// Frees the dynamically allocated arrays from qpms_pitau_t.
+void qpms_pitau_free(qpms_pitau_t);
+//void qpms_pitau_pfree(qpms_pitau_t*);//NI
 
 // Associated Legendre polynomial at zero argument (DLMF 14.5.1) DEPRECATED?
 double qpms_legendre0(int m, int n);

qpms/qpms_types.h

@@ -14,6 +14,9 @@
 #ifndef M_PI
 #define M_PI (3.14159265358979323846)
 #endif
+#ifndef M_SQRT2
+#define M_SQRT2 (1.4142135623730950488)
+#endif
 
 // integer index types
 typedef int qpms_lm_t;
@@ -86,120 +89,83 @@ typedef enum {
 
 
 
-/// Vector spherical wavefuction normalisation (and sign) convention codes.
+/// Vector spherical wavefuction normalisation and phase convention codes.
-/** Throughout the literature, various conventions for VSWF bases are used.
+/**
- *  The meaningful ones are the "power" and "spherical harmonic" normalisation
+ *  Throughout the literature, various conventions for VSWF bases are used.
- *  conventions, as the (\a l, \a m) and (\a l, \a −m) waves of the same type have the same
+ *  These bit flags are used by the functions declared in normalisation.h
- *  intensities.
+ *  that return the appropriate convention-dependent factors.
- *  One might also encounter a very inconvenient and messy "antinormalisation"
- *  used in Xu (TODO reference).
  *
- *  Moreover, VSWFs might use various sign convention. Usually they either
+ *  See @ref vswf_conventions for comparison of the various conventions used.
- *  carry the Condon-Shortley phase \f$ (-1)^m \f$ or not, which is also saved here.
- *
- *  TODO references and exact definitions.
  */
 typedef enum {
+	QPMS_NORMALISATION_UNDEF = 0, ///< Convention undefined. This should not happen.
+       	/// Flag indicating that qpms_normalisition_factor_* should actually return values inverse to the default.
+	QPMS_NORMALISATION_INVERSE = 1, 
+	/** Flag indicating inversion of the asimuthal phase for complex spherical harmonics (i.e. \f$ e^{-im\phi} \f$
+	 * instead of \f$ e^{im\phi} \f$.
+	 */
+	QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE = 2,
+	/// Flag indicating use of the real spherical harmonics.
+	/** If QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE is unset, negative \a m
+	 * correspond to sine in the asimuthal factor; if set, undefined behaviour.
+	 */
+	QPMS_NORMALISATION_SPHARM_REAL = 4,
+	/// Flag indicating usage of Condon-Shortley phase.
+	/** If set, the Ferrers functions and everything derived from them 
+	 * (spherical harmonics, VSWFs) will include a \f$ (-1)^m \f$ factor.
+	 *
+	 * On implementation level, this means that the relevant `gsl_sf_legendre_*_e()`
+	 * functions will be called with argument `csphase = -1.` instead of `+1.`.
+	 */
+	QPMS_NORMALISATION_CSPHASE = 8,
+	QPMS_NORMALISATION_M_I = 16, ///< Include an additional \a i -factor into the magnetic waves.
+	QPMS_NORMALISATION_M_MINUS = 32, ///< Include an additional \f$-1\f$ -factor into the magnetic waves.
+	QPMS_NORMALISATION_N_I = 64, ///< Include an additional \a i -factor into the electric waves.
+	QPMS_NORMALISATION_N_MINUS = 128, ///< Include an additional \f$-1\f$ -factor into the magnetic waves.
+	QPMS_NORMALISATION_L_I = 256, ///< Include an additional \a i -factor into the longitudinal waves.
+	QPMS_NORMALISATION_L_MINUS = 512, ///< Include an additional \f$-1\f$ -factor into the longitudinal waves.
+	QPMS_NORMALISATION_NORM_BITSTART = 65536, 
+	/// The VSWFs shall be power-normalised. This is the "default".
+	/**
+	 * Power normalisation is used e.g. in \cite kristensson_spherical_2014 (complex spherical
+	 * harmonics with Condon-Shortley phase) or \cite kristensson_scattering_2016 (real
+	 * spherical harmonics). This is also the reference for all the other normalisation conventions,
+	 * meaning that qpms_normalisation_factor_M() and qpms_normalisation_factor_N() shall
+	 * always return `1. + 0.*I` if `norm == QPMS_NORMALISATION_NORM_POWER`.
+	 */
+	QPMS_NORMALISATION_NORM_POWER = QPMS_NORMALISATION_NORM_BITSTART * 1,
+	/// The VSWFs shall be normalised as in \cite taylor_optical_2011 .
+	/** This includes a \f$ \sqrt{l(l+1)} \f$ factor compared to the power normalisation. */
+	QPMS_NORMALISATION_NORM_SPHARM = QPMS_NORMALISATION_NORM_BITSTART * 3,
+	/// The VSWFs shall be created using spherical harmonics without any normalisation. Do not use.
+	/** This includes a \f[ 
+	 * 	\sqrt{l(l+1)} \left(\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}\right)^{-\frac{1}{2}}
+	 * \f] factor compared to the power normalisation.
+	 *
+	 * Note that this has no sense whatsoever for real spherical harmonics.
+	 * Again, do not use this.
+	 */
+	QPMS_NORMALISATION_NORM_NONE = QPMS_NORMALISATION_NORM_BITSTART * 2,
+	QPMS_NORMALISATION_NORM_BITS = QPMS_NORMALISATION_NORM_POWER 
+		| QPMS_NORMALISATION_NORM_NONE | QPMS_NORMALISATION_NORM_SPHARM,
 
-#define QPMS_NORMALISATION_T_CSBIT 128 ///< A flag used in qpms_normalisation_t indicating the Condon-Shortley phase.
+	/// VSWF convention used in \cite kristensson_scattering_2016
-
+	QPMS_NORMALISATION_CONVENTION_KRISTENSSON_REAL = QPMS_NORMALISATION_NORM_POWER
-#ifdef USE_XU_ANTINORMALISATION
+		| QPMS_NORMALISATION_SPHARM_REAL,
-	// As in TODO
+	/// VSWF convention used in \cite kristensson_spherical_2014
-	QPMS_NORMALISATION_XU = 4, ///< such that the numerical values in Xu's tables match, not recommended to use otherwise
+	QPMS_NORMALISATION_CONVENTION_KRISTENSSON = QPMS_NORMALISATION_NORM_POWER
-	QPMS_NORMALISATION_XU_CS = QPMS_NORMALISATION_XU | QPMS_NORMALISATION_T_CSBIT, 
+		| QPMS_NORMALISATION_CSPHASE,
-#endif
+	/// VSWF convention used in SCUFF-EM \cite reid_electromagnetism_2016
-	QPMS_NORMALISATION_NONE = 3, ///< genuine unnormalised waves (with unnormalised Legendre polynomials)
+	QPMS_NORMALISATION_CONVENTION_SCUFF = QPMS_NORMALISATION_NORM_POWER
-	QPMS_NORMALISATION_KRISTENSSON = 2, ///< As in http://www.eit.lth.se/fileadmin/eit/courses/eit080f/Literature/book.pdf, power-normalised
+		| QPMS_NORMALISATION_CSPHASE | QPMS_NORMALISATION_M_I
-	QPMS_NORMALISATION_POWER = QPMS_NORMALISATION_KRISTENSSON, 
+		| QPMS_NORMALISATION_N_MINUS
-	// as in TODO
-	QPMS_NORMALISATION_TAYLOR = 1,
-	QPMS_NORMALISATION_SPHARM = QPMS_NORMALISATION_TAYLOR,
-	// Variants with Condon-Shortley phase
-	QPMS_NORMALISATION_NONE_CS = QPMS_NORMALISATION_NONE | QPMS_NORMALISATION_T_CSBIT,
-	QPMS_NORMALISATION_KRISTENSSON_CS = QPMS_NORMALISATION_KRISTENSSON | QPMS_NORMALISATION_T_CSBIT, 
-	QPMS_NORMALISATION_POWER_CS = QPMS_NORMALISATION_KRISTENSSON_CS,
-	QPMS_NORMALISATION_TAYLOR_CS = QPMS_NORMALISATION_TAYLOR | QPMS_NORMALISATION_T_CSBIT,
-	QPMS_NORMALISATION_SPHARM_CS = QPMS_NORMALISATION_TAYLOR_CS,
-	QPMS_NORMALISATION_UNDEF = 0
 } qpms_normalisation_t;
 
 /// Determine whether the convention includes Condon-Shortley phase (-1) or not (+1).
 static inline int qpms_normalisation_t_csphase(qpms_normalisation_t norm) {
-	return (norm & QPMS_NORMALISATION_T_CSBIT)? -1 : 1;
+	return (norm & QPMS_NORMALISATION_CSPHASE)? -1 : 1;
 }
 
-/// Returns the normalisation convention code without the Condon-Shortley phase.
-static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) {
-	return norm & (~QPMS_NORMALISATION_T_CSBIT);
-}
-
-// TODO move the inlines elsewhere
-/* Normalisation of the spherical waves is now scattered in at least three different files:
- * here, we have the norm in terms of radiated power of outgoing wave.
- * In file legendre.c, function qpms_pitau_get determines the norm used in the vswf.c
- * spherical vector wave norms. The "dual" waves in vswf.c use the ..._abssquare function below.
- * In file translations.c, the normalisations are again set by hand using the normfac and lognormfac
- * functions.
- */
-#include <math.h>
-#include <assert.h>
-// relative to QPMS_NORMALISATION_KRISTENSSON_CS, i.e.
-// P_l^m[normtype] = P_l^m[Kristensson]
-static inline double qpms_normalisation_t_factor(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
-	int csphase = qpms_normalisation_t_csphase(norm);
-	norm = qpms_normalisation_t_normonly(norm);
-	double factor;
-	switch (norm) {
-		case QPMS_NORMALISATION_KRISTENSSON:
-			factor = 1.;
-			break;
-		case QPMS_NORMALISATION_TAYLOR:
-			factor = sqrt(l*(l+1));
-			break;
-		case QPMS_NORMALISATION_NONE:
-			factor = sqrt(l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)));
-			break;
-#ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c
-		case QPMS_NORMALISATION_XU:
-			factor = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
-			break;
-#endif
-		default:
-			assert(0);
-	}
-	factor *= (m%2)?(-csphase):1;
-	return factor;
-}
-
-
-// TODO move elsewhere
-static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
-	norm = qpms_normalisation_t_normonly(norm);
-	switch (norm) {
-		case QPMS_NORMALISATION_KRISTENSSON:
-			return 1.;
-			break;
-		case QPMS_NORMALISATION_TAYLOR:
-			return l*(l+1);
-			break;
-		case QPMS_NORMALISATION_NONE:
-			return l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
-			break;
-#ifdef USE_XU_ANTINORMALISATION // broken probably in legendre.c
-		case QPMS_NORMALISATION_XU:
-			{
-			  double fac = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
-			  return fac * fac;
-			}
-			break;
-#endif 
-		default:
-			assert(0);
-			return NAN;
-	}
-}
-
-
 /// Bessel function kinds.
 typedef enum {
 	QPMS_BESSEL_REGULAR = 1, ///< regular (spherical) Bessel function \a j (Bessel function of the first kind)

qpms/scatsystem.c

@@ -431,10 +431,10 @@ void qpms_scatsys_free(qpms_scatsys_t *ss) {
 // (perhaps just use qpms_normalisation_t_factor() at the right places)
 static inline void check_norm_compat(const qpms_vswf_set_spec_t *s)
 {
-  switch (qpms_normalisation_t_normonly(s->norm)) {
+  switch (s->norm & QPMS_NORMALISATION_NORM_BITS) {
-    case QPMS_NORMALISATION_POWER:
+    case QPMS_NORMALISATION_NORM_POWER:
       break;
-    case QPMS_NORMALISATION_SPHARM:
+    case QPMS_NORMALISATION_NORM_SPHARM:
       break;
     default:
       abort(); // Only SPHARM and POWER norms are supported right now.
@@ -447,7 +447,7 @@ complex double *qpms_orbit_action_matrix(complex double *target,
   assert(sym); assert(g < sym->order);
   assert(sym->rep3d);
   assert(ot); assert(ot->size > 0);
-  check_norm_compat(bspec);
+  // check_norm_compat(bspec); not needed here, the qpms_irot3_uvswfi_dense should complain if necessary
   const size_t n = bspec->n;
   const qpms_gmi_t N = ot->size;
   if (target == NULL)
@@ -491,7 +491,7 @@ complex double *qpms_orbit_irrep_projector_matrix(complex double *target,
   assert(sym->rep3d);
   assert(ot); assert(ot->size > 0);
   assert(iri < sym->nirreps); assert(sym->irreps);
-  check_norm_compat(bspec);
+  // check_norm_compat(bspec); // probably not needed here, let the called functions complain if necessary, but CHEKME
   const size_t n = bspec->n;
   const qpms_gmi_t N = ot->size;
   if (target == NULL)
@@ -545,7 +545,7 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size,
   assert(sym->rep3d);
   assert(ot); assert(ot->size > 0);
   assert(iri < sym->nirreps); assert(sym->irreps);
-  check_norm_compat(bspec);
+  check_norm_compat(bspec); // Here I'm not sure; CHECKME
   const size_t n = bspec->n;
   const qpms_gmi_t N = ot->size;
   const bool newtarget = (target == NULL);
@@ -1451,7 +1451,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR(
   for(long thi = 0; thi < nthreads; ++thi)
     QPMS_ENSURE_SUCCESS(pthread_create(thread_ids + thi, NULL,
       qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread,
-      &arg));
+      (void *) &arg));
   for(long thi = 0; thi < nthreads; ++thi) {
     void *retval;
     QPMS_ENSURE_SUCCESS(pthread_join(thread_ids[thi], &retval));

qpms/symmetries.c

@@ -2,38 +2,46 @@
 #include "tiny_inlines.h"
 #include "indexing.h"
 #include "wigner.h"
+#include "qpms_error.h"
 
 // TODO at some point, maybe support also other norms.
 // (perhaps just use qpms_normalisation_t_factor() at the right places)
 static inline void check_norm_compat(const qpms_vswf_set_spec_t *s)
 {
-  switch (qpms_normalisation_t_normonly(s->norm)) {
+  switch (s->norm & QPMS_NORMALISATION_NORM_BITS) {
-    case QPMS_NORMALISATION_POWER:
+    case QPMS_NORMALISATION_NORM_POWER:
       break;
-    case QPMS_NORMALISATION_SPHARM:
+    case QPMS_NORMALISATION_NORM_SPHARM:
       break;
     default:
-      abort(); // Only SPHARM and POWER norms are supported right now.
+      QPMS_NOT_IMPLEMENTED("At the moment, only spherical harmonics of spherical harmonics or power normalisations implemented.");
   }
 }
 
+
+static inline void ONLY_EIMF_IMPLEMENTED(const qpms_normalisation_t norm)
+{
+	if (norm & QPMS_NORMALISATION_SPHARM_REAL)
+		QPMS_NOT_IMPLEMENTED("Support for real spherical harmonics not implemented yet.");
+}
+
+
 // Used in the functions below to ensure memory allocation and checks for bspec validity
 static inline complex double *ensure_alloc(complex double *target,
     const qpms_vswf_set_spec_t *bspec) {
   check_norm_compat(bspec);
   const size_t n = bspec->n;
   if (target == NULL)
-    target = malloc(n * n * sizeof(complex double));
+    QPMS_CRASHING_MALLOC(target, n * n * sizeof(complex double));
-  if (target == NULL) abort();
   return target;
 }
 
 
-
 complex double *qpms_zflip_uvswi_dense(
     complex double *target,
     const qpms_vswf_set_spec_t *bspec) 
 {
+  check_norm_compat(bspec);
   target = ensure_alloc(target, bspec);
   const size_t n = bspec->n;
 
@@ -69,6 +77,8 @@ complex double *qpms_yflip_uvswi_dense(
     complex double *target,
     const qpms_vswf_set_spec_t *bspec) 
 {
+  check_norm_compat(bspec);
+  ONLY_EIMF_IMPLEMENTED(bspec->norm);
   target = ensure_alloc(target, bspec);
   const size_t n = bspec->n;
 
@@ -104,6 +114,8 @@ complex double *qpms_xflip_uvswi_dense(
     complex double *target,
     const qpms_vswf_set_spec_t *bspec) 
 {
+  check_norm_compat(bspec);
+  ONLY_EIMF_IMPLEMENTED(bspec->norm);
   target = ensure_alloc(target, bspec);
   const size_t n = bspec->n;
 
@@ -142,6 +154,9 @@ complex double *qpms_zrot_uvswi_dense(
                 double phi ///< Rotation angle
                 )
 {
+  QPMS_UNTESTED; // not sure about the C.-S. phase. Don't forget documenting it as well.
+  check_norm_compat(bspec);
+  ONLY_EIMF_IMPLEMENTED(bspec->norm);
   target = ensure_alloc(target, bspec);
   const size_t n = bspec->n;
 
@@ -182,6 +197,9 @@ complex double *qpms_irot3_uvswfi_dense(
     const qpms_vswf_set_spec_t *bspec,
     const qpms_irot3_t t)
 {
+  QPMS_UNTESTED; // not sure about the C.-S. phase. Don't forget documenting it as well.
+  check_norm_compat(bspec);
+  ONLY_EIMF_IMPLEMENTED(bspec->norm);
   target = ensure_alloc(target, bspec);
   const size_t n = bspec->n;
 

qpms/tests/test_vswf_translations.c

@@ -12,18 +12,14 @@
 
 char *normstr(qpms_normalisation_t norm) {
 	//int csphase = qpms_normalisation_t_csphase(norm);
-	norm = qpms_normalisation_t_normonly(norm);
+	norm = norm & QPMS_NORMALISATION_NORM_BITS;
 	switch (norm) {
-		case QPMS_NORMALISATION_NONE:
+		case QPMS_NORMALISATION_NORM_NONE:
 			return "none";
-		case QPMS_NORMALISATION_SPHARM:
+		case QPMS_NORMALISATION_NORM_SPHARM:
 			return "spharm";
-		case QPMS_NORMALISATION_POWER:
+		case QPMS_NORMALISATION_NORM_POWER:
 			return "power";
-#ifdef USE_XU_ANTINORMALISATION
-    case QPMS_NORMALISATION_XU:
-      return "xu";
-#endif
 		default:
 			return "!!!undef!!!";
 	}
@@ -62,14 +58,13 @@ int main() {
 	}
 	
 	for(int use_csbit = 0; use_csbit <= 1; ++use_csbit) {
-		for(int i = 1; 
+		for(int i = 0; i < 3; ++i){
-#ifdef USE_XU_ANTINORMALISATION
+			qpms_normalisation_t norm = ((qpms_normalisation_t[]) 
-        i <= 4; 
+          { QPMS_NORMALISATION_NORM_SPHARM,
-#else
+            QPMS_NORMALISATION_NORM_POWER,
-        i <= 3; 
+            QPMS_NORMALISATION_NORM_NONE
-#endif
+            })[i]
-        ++i){
+        | (use_csbit ? QPMS_NORMALISATION_CSPHASE : 0);
-			qpms_normalisation_t norm = i | (use_csbit ? QPMS_NORMALISATION_T_CSBIT : 0);
 			qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, norm);
 			for(int J = 1; J <= 4; ++J)
 				test_sphwave_translation(c, J, o2minuso1, npoints, points);

qpms/tmatrices.c

@@ -17,6 +17,7 @@
 #include "qpms_error.h"
 #include "tmatrices.h"
 #include "qpms_specfunc.h"
+#include "normalisation.h"
 
 #define HBAR (1.05457162825e-34)
 #define ELECTRONVOLT (1.602176487e-19)
@@ -352,10 +353,9 @@ qpms_errno_t qpms_read_scuff_tmatrix(
   if (!(freqs && n && tmdata)) 
     qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
       "freqs, n, and tmdata are mandatory arguments and must not be NULL.");
-  if (bs->norm != QPMS_NORMALISATION_POWER_CS) // CHECKME CORRECT?
+  if(bs->norm & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE 
-    qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
+        | QPMS_NORMALISATION_SPHARM_REAL))
-        "Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)"
+    QPMS_NOT_IMPLEMENTED("Sorry, only standard complex-spherical harmonic based waves are supported right now");
-        " norm supported right now.");
   int n_alloc = 128; // First chunk to allocate
   *n = 0;
   *freqs = malloc(n_alloc * sizeof(double));
@@ -413,7 +413,11 @@ qpms_errno_t qpms_read_scuff_tmatrix(
       /* This element has not been requested in bs->ilist. */
       continue;
     else
-      (*tmdata)[(*n-1)*bs->n*bs->n + desti*bs->n + srci] = tr + I*ti;
+      (*tmdata)[(*n-1)*bs->n*bs->n + desti*bs->n + srci] = (tr + I*ti)
+        * qpms_normalisation_factor_uvswfi(bs->norm, srcui)
+        / qpms_normalisation_factor_uvswfi(bs->norm, destui)
+        * qpms_normalisation_factor_uvswfi(QPMS_NORMALISATION_CONVENTION_SCUFF, destui)
+        / qpms_normalisation_factor_uvswfi(QPMS_NORMALISATION_CONVENTION_SCUFF, srcui);
   }
   free(linebuf);
   // free some more memory

qpms/tmatrices.h

@@ -142,6 +142,13 @@ qpms_errno_t qpms_load_scuff_tmatrix(
 /// Loads a scuff-tmatrix generated file.
 /** A simple wrapper over qpms_read_scuff_tmatrix() that needs a 
  * path instead of open FILE.
+ *
+ * The T-matrix is transformed from the VSWF basis defined by
+ * QPMS_NORMALISATION_CONVENTION_SCUFF into the basis defined
+ * by convention bspec->norm.
+ * 
+ * Right now, bspec->norm with real or "reversed complex" spherical
+ * harmonics are not supported.
  */
 qpms_errno_t qpms_read_scuff_tmatrix(
 		FILE *f, ///< An open stream with the T-matrix data.

qpms/translations.c

@@ -13,7 +13,7 @@
 #include <stdlib.h> //abort()
 #include <gsl/gsl_sf_coupling.h>
 #include "qpms_error.h"
-
+#include "normalisation.h"
 
 /*
  * Define macros with additional factors that "should not be there" according
@@ -53,6 +53,12 @@
 #endif
 
 
+// Translation operators for real sph. harm. based waves are not yet implemented...
+static inline void TROPS_ONLY_EIMF_IMPLEMENTED(qpms_normalisation_t norm) {
+  if (norm & (QPMS_NORMALISATION_SPHARM_REAL | QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE)) 
+    QPMS_NOT_IMPLEMENTED("Translation operators for real or inverse complex spherical harmonics based waves are not implemented.");
+}
+
 /*
  * References:
  * [Xu_old] Yu-Lin Xu, Journal of Computational Physics 127, 285–298 (1996)
@@ -127,60 +133,26 @@ int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *
 
 static inline double qpms_trans_normlogfac(qpms_normalisation_t norm,
     int m, int n, int mu, int nu) {
-  //int csphase = qpms_normalisation_t csphase(norm); // probably not needed here
-  norm = qpms_normalisation_t_normonly(norm);
-  switch(norm) {
-    case QPMS_NORMALISATION_KRISTENSSON:
-    case QPMS_NORMALISATION_TAYLOR:
       return -0.5*(lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1));
-      break;
-    case QPMS_NORMALISATION_NONE:
-      return -(lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1));
-      break;
-#ifdef USE_XU_ANTINORMALISATION
-    case QPMS_NORMALISATION_XU:
-      return 0;
-      break;
-#endif
-    default:
-      abort();
-  }
 }
 
 static inline double qpms_trans_normfac(qpms_normalisation_t norm,
     int m, int n, int mu, int nu) {
   int csphase = qpms_normalisation_t_csphase(norm);
-  norm = qpms_normalisation_t_normonly(norm);
   /* Account for csphase here. Alternatively, this could be done by
    * using appropriate csphase in the legendre polynomials when calculating
    * the translation operator.
    */
   double normfac = (1 == csphase) ? min1pow(m-mu) : 1.;
-  switch(norm) {
+  normfac *= sqrt((n*(n+1.))/(nu*(nu+1.)));
-    case QPMS_NORMALISATION_KRISTENSSON:
+  normfac *= sqrt((2.*n+1)/(2.*nu+1));
-      normfac *= sqrt((n*(n+1.))/(nu*(nu+1.)));
-      normfac *= sqrt((2.*n+1)/(2.*nu+1));
-      break;
-    case QPMS_NORMALISATION_TAYLOR:
-      normfac *= sqrt((2.*n+1)/(2.*nu+1));
-      break;
-    case QPMS_NORMALISATION_NONE:
-      normfac *= (2.*n+1)/(2.*nu+1);
-      break;
-#ifdef USE_XU_ANTINORMALISATION
-    case QPMS_NORMALISATION_XU:
-      break;
-#endif
-    default:
-      abort();
-  }
-
   return normfac;
 }
 
 complex double qpms_trans_single_A(qpms_normalisation_t norm,
     int m, int n, int mu, int nu, sph_t kdlj,
     bool r_ge_d, qpms_bessel_t J) {
+  TROPS_ONLY_EIMF_IMPLEMENTED(norm);
   if(r_ge_d) J = QPMS_BESSEL_REGULAR;
 
   double costheta = cos(kdlj.theta);
@@ -221,7 +193,9 @@ complex double qpms_trans_single_A(qpms_normalisation_t norm,
 
   double normlogfac = qpms_trans_normlogfac(norm,m,n,mu,nu);
   double normfac = qpms_trans_normfac(norm,m,n,mu,nu);
-
+  /// N<-N type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^.
+  normfac *=  qpms_normalisation_factor_N_noCS(norm, nu, mu)
+              / qpms_normalisation_factor_N_noCS(norm, n, m);
   // ipow(n-nu) is the difference from the Taylor formula!
   presum *= /*ipow(n-nu) * */
     (normfac * exp(normlogfac))
@@ -352,6 +326,7 @@ complex double qpms_trans_single_B_Xu(int m, int n, int mu, int nu, sph_t kdlj,
 complex double qpms_trans_single_B(qpms_normalisation_t norm,
     int m, int n, int mu, int nu, sph_t kdlj,
     bool r_ge_d, qpms_bessel_t J) { 
+  TROPS_ONLY_EIMF_IMPLEMENTED(norm);
 #ifndef USE_BROKEN_SINGLETC
   assert(0); // FIXME probably gives wrong values, do not use.
 #endif
@@ -408,6 +383,9 @@ complex double qpms_trans_single_B(qpms_normalisation_t norm,
 
   double normlogfac = qpms_trans_normlogfac(norm,m,n,mu,nu);
   double normfac = qpms_trans_normfac(norm,m,n,mu,nu);
+  /// N<-M type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^.
+  normfac *=  qpms_normalisation_factor_M_noCS(norm, nu, mu)
+              / qpms_normalisation_factor_N_noCS(norm, n, m);
 
   // ipow(n-nu) is the difference from the "old Taylor" formula	
   presum *= /*ipow(n-nu) * */(exp(normlogfac) * normfac)
@@ -546,6 +524,9 @@ static void qpms_trans_calculator_multipliers_A_general(
 
   double normlogfac = qpms_trans_normlogfac(norm,m,n,mu,nu);
   double normfac = qpms_trans_normfac(norm,m,n,mu,nu);
+  /// N<-N type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^.
+  normfac *=  qpms_normalisation_factor_N_noCS(norm, nu, mu)
+              / qpms_normalisation_factor_N_noCS(norm, n, m);
 
   normfac *= min1pow(m); //different from old Taylor
 
@@ -610,6 +591,9 @@ void qpms_trans_calculator_multipliers_B_general(
 
   double normlogfac= qpms_trans_normlogfac(norm,m,n,mu,nu);
   double normfac = qpms_trans_normfac(norm,m,n,mu,nu);
+  /// N<-M type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^.
+  normfac *=  qpms_normalisation_factor_M_noCS(norm, nu, mu)
+              / qpms_normalisation_factor_N_noCS(norm, n, m);
 
   double presum = min1pow(1-m) * (2*nu+1)/(2.*(n*(n+1))) 
     * exp(lgamma(n+m+1) - lgamma(n-m+1) + lgamma(nu-mu+1) - lgamma(nu+mu+1)
@@ -680,9 +664,11 @@ void qpms_trans_calculator_multipliers_B_general(
 
 
   int csphase = qpms_normalisation_t_csphase(norm); //TODO FIXME use this
-  norm = qpms_normalisation_t_normonly(norm);
   double normlogfac= qpms_trans_normlogfac(norm,m,n,mu,nu);
   double normfac = qpms_trans_normfac(norm,m,n,mu,nu);
+  /// N<-M type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^.
+  normfac *=  qpms_normalisation_factor_M_noCS(norm, nu, mu)
+              / qpms_normalisation_factor_N_noCS(norm, n, m);
   // TODO use csphase to modify normfac here!!!!
   // normfac = xxx ? -normfac : normfac;
   normfac *= min1pow(m);//different from old taylor
@@ -831,51 +817,18 @@ static void qpms_trans_calculator_multipliers_B_Taylor(
 }
 
 int qpms_trans_calculator_multipliers_A(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) {
-  switch (qpms_normalisation_t_normonly(norm)) {
-    case QPMS_NORMALISATION_TAYLOR:
-#ifdef USE_SEPARATE_TAYLOR
-      qpms_trans_calculator_multipliers_A_Taylor(dest,m,n,mu,nu,qmax);
-      return 0;
-      break;
-#endif
-    case QPMS_NORMALISATION_NONE:
-#ifdef USE_XU_ANTINORMALISATION
-    case QPMS_NORMALISATION_XU:
-#endif   
-    case QPMS_NORMALISATION_KRISTENSSON:
       qpms_trans_calculator_multipliers_A_general(norm, dest, m, n, mu, nu, qmax);
       return 0;
-      break;
-    default:
-      abort();
-  }
-  assert(0);
 }
 
 int qpms_trans_calculator_multipliers_B(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int Qmax) {
-  switch (qpms_normalisation_t_normonly(norm)) {
-    case QPMS_NORMALISATION_TAYLOR:
-#ifdef USE_SEPARATE_TAYLOR
-      qpms_trans_calculator_multipliers_B_Taylor(dest,m,n,mu,nu,Qmax);
-      return 0;
-      break;
-#endif
-    case QPMS_NORMALISATION_NONE:
-#ifdef USE_XU_ANTINORMALISATION
-    case QPMS_NORMALISATION_XU:
-#endif   
-    case QPMS_NORMALISATION_KRISTENSSON:
       qpms_trans_calculator_multipliers_B_general(norm, dest, m, n, mu, nu, Qmax);
       return 0;
-      break;
-    default:
-      abort();
-  }
-  assert(0);
 }
 
 qpms_trans_calculator
 *qpms_trans_calculator_init (const int lMax, const qpms_normalisation_t normalisation) {
+  TROPS_ONLY_EIMF_IMPLEMENTED(normalisation);
   assert(lMax > 0);
   qpms_trans_calculator *c = malloc(sizeof(qpms_trans_calculator));
   c->lMax = lMax;
@@ -951,6 +904,7 @@ static inline complex double qpms_trans_calculator_get_A_precalcbuf(const qpms_t
     int m, int n, int mu, int nu, sph_t kdlj,
     bool r_ge_d, qpms_bessel_t J,
     const complex double *bessel_buf, const double *legendre_buf) {
+  TROPS_ONLY_EIMF_IMPLEMENTED(c->normalisation);
   size_t i = qpms_trans_calculator_index_mnmunu(c, m, n, mu, nu);
   size_t qmax = c->A_multipliers[i+1] - c->A_multipliers[i] - 1;
   assert(qmax == gaunt_q_max(-m,n,mu,nu));
@@ -977,33 +931,20 @@ complex double qpms_trans_calculator_get_A_buf(const qpms_trans_calculator *c,
     // TODO warn? 
     return NAN+I*NAN;
   int csphase = qpms_normalisation_t_csphase(c->normalisation);
-  switch(qpms_normalisation_t_normonly(c->normalisation)) {
+  
-    // TODO use normalised legendre functions for Taylor and Kristensson
+  double costheta = cos(kdlj.theta);
-    case QPMS_NORMALISATION_TAYLOR:
+  if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,
-    case QPMS_NORMALISATION_KRISTENSSON:
+        costheta,csphase,legendre_buf)) abort();
-    case QPMS_NORMALISATION_NONE:
+  if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort();
-#ifdef USE_XU_ANTINORMALISATION
+  return qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
-    case QPMS_NORMALISATION_XU:   
+      kdlj,r_ge_d,J,bessel_buf,legendre_buf);
-#endif
-      {
-        double costheta = cos(kdlj.theta);
-        if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,
-              costheta,csphase,legendre_buf)) abort();
-        if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort();
-        return qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
-            kdlj,r_ge_d,J,bessel_buf,legendre_buf);
-      }
-      break;
-    default:
-      abort();
-  }
-  assert(0);
 }
 
 static inline complex double qpms_trans_calculator_get_B_precalcbuf(const qpms_trans_calculator *c,
     int m, int n, int mu, int nu, sph_t kdlj,
     bool r_ge_d, qpms_bessel_t J,
     const complex double *bessel_buf, const double *legendre_buf) {
+  TROPS_ONLY_EIMF_IMPLEMENTED(c->normalisation);
   size_t i = qpms_trans_calculator_index_mnmunu(c, m, n, mu, nu);
   size_t qmax = c->B_multipliers[i+1] - c->B_multipliers[i] - (1 - BQ_OFFSET);
   assert(qmax == gauntB_Q_max(-m,n,mu,nu));
@@ -1030,26 +971,12 @@ complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator *c,
     // TODO warn? 
     return NAN+I*NAN;
   int csphase = qpms_normalisation_t_csphase(c->normalisation);
-  switch(qpms_normalisation_t_normonly(c->normalisation)) {
+  double costheta = cos(kdlj.theta);
-    case QPMS_NORMALISATION_TAYLOR:
+  if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
-    case QPMS_NORMALISATION_KRISTENSSON:
+        costheta,csphase,legendre_buf)) abort();
-    case QPMS_NORMALISATION_NONE:
+  if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort();
-#ifdef USE_XU_ANTINORMALISATION
+  return qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
-    case QPMS_NORMALISATION_XU:   
+      kdlj,r_ge_d,J,bessel_buf,legendre_buf);
-#endif
-      {
-        double costheta = cos(kdlj.theta);
-        if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
-              costheta,csphase,legendre_buf)) abort();
-        if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort();
-        return qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
-            kdlj,r_ge_d,J,bessel_buf,legendre_buf);
-      }
-      break;
-    default:
-      abort();
-  }
-  assert(0);
 }
 
 int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c,
@@ -1064,29 +991,15 @@ int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c,
     // TODO warn? different return value?
     return 0;
   }
-  switch(qpms_normalisation_t_normonly(c->normalisation)) {
+  double costheta = cos(kdlj.theta);
-    case QPMS_NORMALISATION_TAYLOR:
+  if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
-    case QPMS_NORMALISATION_KRISTENSSON:
+        costheta,-1,legendre_buf)) abort();
-    case QPMS_NORMALISATION_NONE:
+  if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort();
-#ifdef USE_XU_ANTINORMALISATION
+  *Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
-    case QPMS_NORMALISATION_XU:   
+      kdlj,r_ge_d,J,bessel_buf,legendre_buf);
-#endif
+  *Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
-      {
+      kdlj,r_ge_d,J,bessel_buf,legendre_buf);
-        double costheta = cos(kdlj.theta);
+  return 0;
-        if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
-              costheta,-1,legendre_buf)) abort();
-        if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort();
-        *Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
-            kdlj,r_ge_d,J,bessel_buf,legendre_buf);
-        *Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
-            kdlj,r_ge_d,J,bessel_buf,legendre_buf);
-        return 0;
-      }
-      break;
-    default:
-      abort();
-  }
-  assert(0);
 }
 
 
@@ -1105,40 +1018,31 @@ int qpms_trans_calculator_get_AB_arrays_buf(const qpms_trans_calculator *c,
     // TODO warn? different return value?
     return 0;
   }
-  switch(qpms_normalisation_t_normonly(c->normalisation)) {
+  {
-    case QPMS_NORMALISATION_TAYLOR:
+    double costheta = cos(kdlj.theta);
-    case QPMS_NORMALISATION_POWER:
+    if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1,
-    case QPMS_NORMALISATION_NONE:
+          costheta,-1,legendre_buf)) abort();
-      {
+    if (qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf)) abort();
-        double costheta = cos(kdlj.theta);
+    size_t desti = 0, srci = 0;
-        if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1,
+    for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) {
-              costheta,-1,legendre_buf)) abort();
+      for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) {
-        if (qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf)) abort();
-        size_t desti = 0, srci = 0;
-        for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) {
-          for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) {
 #ifndef NDEBUG
-            size_t assertindex = qpms_trans_calculator_index_mnmunu(c,m,n,mu,nu);
+        size_t assertindex = qpms_trans_calculator_index_mnmunu(c,m,n,mu,nu);
 #endif
-            assert(assertindex == desti*c->nelem + srci);
+        assert(assertindex == desti*c->nelem + srci);
-            *(Adest + deststride * desti + srcstride * srci) = 
+        *(Adest + deststride * desti + srcstride * srci) = 
-              qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
+          qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
-                  kdlj,r_ge_d,J,bessel_buf,legendre_buf);
+              kdlj,r_ge_d,J,bessel_buf,legendre_buf);
-            *(Bdest + deststride * desti + srcstride * srci) = 
+        *(Bdest + deststride * desti + srcstride * srci) = 
-              qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
+          qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
-                  kdlj,r_ge_d,J,bessel_buf,legendre_buf);
+              kdlj,r_ge_d,J,bessel_buf,legendre_buf);
-            ++srci;
+        ++srci;
-          }
-          ++desti;
-          srci = 0;
-        }
-        return 0;
       }
-      break;
+      ++desti;
-    default:
+      srci = 0;
-      abort();
+    }
+    return 0;
   }
-  assert(0);
 }
 
 complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c,
@@ -1287,10 +1191,6 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
     if(doerr) serr_total[y] += sigma0_err;
   }
   
-  switch(qpms_normalisation_t_normonly(c->normalisation)) {
-    case QPMS_NORMALISATION_TAYLOR:
-    case QPMS_NORMALISATION_POWER:
-    case QPMS_NORMALISATION_NONE:
       {
         ptrdiff_t desti = 0, srci = 0;
         for (qpms_l_t n = 1; n <= c->lMax; ++n) for (qpms_m_t m = -n; m <= n; ++m) {
@@ -1335,10 +1235,6 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
           srci = 0;
         }
       }
-      break;
-    default:
-      abort();
-  }
 
   free(sigmas_short);
   free(sigmas_long);
@@ -1563,10 +1459,6 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
     if(doerr) serr_total[y] += sigma0_err;
   }
   
-  switch(qpms_normalisation_t_normonly(c->normalisation)) {
-    case QPMS_NORMALISATION_TAYLOR:
-    case QPMS_NORMALISATION_POWER:
-    case QPMS_NORMALISATION_NONE:
       {
         ptrdiff_t desti = 0, srci = 0;
         for (qpms_l_t n = 1; n <= c->lMax; ++n) for (qpms_m_t m = -n; m <= n; ++m) {
@@ -1611,10 +1503,6 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
           srci = 0;
         }
       }
-      break;
-    default:
-      abort();
-  }
 
   free(sigmas_short);
   free(sigmas_long);
@@ -1670,10 +1558,6 @@ int qpms_trans_calculator_get_shortrange_AB_arrays_buf(const qpms_trans_calculat
     // TODO warn? different return value?
     return 0;
   }
-  switch(qpms_normalisation_t_normonly(c->normalisation)) {
-    case QPMS_NORMALISATION_TAYLOR:
-    case QPMS_NORMALISATION_POWER:
-    case QPMS_NORMALISATION_NONE: 
       {
         double costheta = cos(kdlj.theta);
         if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1,
@@ -1698,11 +1582,6 @@ int qpms_trans_calculator_get_shortrange_AB_arrays_buf(const qpms_trans_calculat
         }
         return 0;
       }
-      break;
-    default:
-      abort();
-  }
-  assert(0);
 }
 
 int qpms_trans_calculator_get_shortrange_AB_buf_p(const qpms_trans_calculator *c,
@@ -1718,10 +1597,6 @@ int qpms_trans_calculator_get_shortrange_AB_buf_p(const qpms_trans_calculator *c
     // TODO warn? different return value?
     return 0;
   }
-  switch(qpms_normalisation_t_normonly(c->normalisation)) {
-    case QPMS_NORMALISATION_TAYLOR:
-    case QPMS_NORMALISATION_KRISTENSSON:
-    case QPMS_NORMALISATION_NONE:
       {
         double costheta = cos(kdlj.theta);
         if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
@@ -1735,11 +1610,6 @@ int qpms_trans_calculator_get_shortrange_AB_buf_p(const qpms_trans_calculator *c
             kdlj,false,J,bessel_buf,legendre_buf);
         return 0;
       }
-      break;
-    default:
-      abort();
-  }
-  assert(0);
 }
 
 // Short-range parts of the translation coefficients
@@ -1826,13 +1696,6 @@ int qpms_trans_calculator_get_2DFT_longrange_AB_buf_p(const qpms_trans_calculato
   assert (J == QPMS_HANKEL_PLUS);
   assert(k_sph.theta == M_PI_2);
   
-  switch(qpms_normalisation_t_normonly(c->normalisation)) {
-    case QPMS_NORMALISATION_TAYLOR:
-    case QPMS_NORMALISATION_KRISTENSSON:
-    case QPMS_NORMALISATION_NONE:
-#ifdef USE_XU_ANTINORMALISATION
-    case QPMS_NORMALISATION_XU:   
-#endif
       {
         //double costheta = cos(kdlj.theta);
         //if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
@@ -1846,11 +1709,6 @@ int qpms_trans_calculator_get_2DFT_longrange_AB_buf_p(const qpms_trans_calculato
             k_sph,J,lrhankel_recparts_buf);
         return 0;
       }
-      break;
-    default:
-      abort();
-  }
-  assert(0);
 }
 
 // Fourier transforms of the long-range parts of the translation coefficients
@@ -1884,10 +1742,6 @@ int qpms_trans_calculator_get_2DFT_longrange_AB_arrays_buf(const qpms_trans_calc
     return 0;
   }
 #endif 
-  switch(qpms_normalisation_t_normonly(c->normalisation)) {
-    case QPMS_NORMALISATION_TAYLOR:
-    case QPMS_NORMALISATION_POWER:
-    case QPMS_NORMALISATION_NONE:
       {
         lrhankel_recpart_fill(lrhankel_recparts_buf, 2*c->lMax+1,
             lrk_cutoff, c->hct, kappa, cv, k0, k_sph.r);
@@ -1910,11 +1764,6 @@ int qpms_trans_calculator_get_2DFT_longrange_AB_arrays_buf(const qpms_trans_calc
         }
         return 0;
       }
-      break;
-    default:
-      abort();
-  }
-  assert(0);
 }
 
 // FIXME i get stack smashing error inside the following function if compiled with optimizations

qpms/translations.h

@@ -1,3 +1,29 @@
+/*! \file translations.h
+ *  \brief VSWF translation operator.
+ *
+ * ### Argument conventions
+ *
+ * A single wave with indices mu, nu is re-expanded at kdlj into waves with indices m, n,
+ * i.e. in the following functions, the first arguments over which one sums (multiplied
+ * by the waves with new origin).
+ *
+ * HOWEVER, this means that if a field has an expansion with coeffs a(mu, nu)
+ * at the original origin, with the basis at the new origin, the coeffs will be
+ * a(m, n) = \sum_{mu,nu} A(m, n, mu, nu) a(mu, nu).
+ *
+ * With qpms_trans_calculator_get_AB_arrays_buf (and other functions from *AB_arrays*
+ * family), one can choose the stride. And it seems that the former stride argument (now called
+ * destride) and the latter (now called srcstride) are connected to (m,n) and (mu,nu) indices,
+ * respectively. Seems consistent.
+ *
+ *
+ * #### r_ge_d argument:
+ *
+ * If r_ge_d == true, the translation coefficients are calculated using regular bessel functions, 
+ * regardless of what J argument is.
+ *
+ *
+ */
 #ifndef QPMS_TRANSLATIONS_H
 #define QPMS_TRANSLATIONS_H
 #include "vectors.h"
@@ -14,33 +40,6 @@
 #include "ewald.h"
 #endif
 
-/*
- * Argument conventions:
- *
- * A single wave with indices mu, nu is re-expanded at kdlj into waves with indices m, n,
- * i.e. in the following functions, the first arguments over which one sums (multiplied
- * by the waves with new origin).
- *
- * HOWEVER, this means that if a field has an expansion with coeffs a(mu, nu)
- * at the original origin, with the basis at the new origin, the coeffs will be
- * a(m, n) = \sum_{mu,nu} A(m, n, mu, nu) a(mu, nu).
- *
- * With qpms_trans_calculator_get_AB_arrays_buf (and other functions from *AB_arrays*
- * family), one can choose the stride. And it seems that the former stride argument (now called
- * destride) and the latter (now called srcstride) are connected to (m,n) and (mu,nu) indices,
- * respectively. Seems consistent.
- *
- */
-
-/*
- * r_ge_d argument:
- *
- * If r_ge_d == true, the translation coefficients are calculated using regular bessel functions, 
- * regardless of what J argument is.
- *
- */
-
-
 // TODO replace the xplicit "Taylor" functions with general,
 // taking qpms_normalisation_t argument.
 complex double qpms_trans_single_A_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
@@ -61,6 +60,24 @@ complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l
 complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
 		bool r_ge_d, qpms_bessel_t J);
 
+/// Structure holding the constant factors in normalisation operators.
+/**
+ * The VSWF translation operator elements are rather complicated linear
+ * combinations of Bessel functions and spherical harmonics.
+ * The preceding factors are rather complicated but need to be calculated
+ * (for a single normalisation convention)
+ * only once and then recycled during the operator evaluation for different
+ * translations.
+ *
+ * This structure is initialised with qpms_trans_calculator_t() and
+ * holds all the constant factors up to a truncation
+ * degree \a lMax.
+ *
+ * The destructor function is qpms_trans_calculator_free().
+ *
+ * If Ewald sums are enabled at build, it also holds the constant
+ * factors useful for lattice sums of translation operator.
+ */
 typedef struct qpms_trans_calculator {
 	qpms_normalisation_t normalisation;
 	qpms_l_t lMax;
@@ -86,8 +103,11 @@ typedef struct qpms_trans_calculator {
 	double *legendre0; // Zero-argument Legendre functions – this might go outside #ifdef in the end...
 } qpms_trans_calculator;
 
-
+/// Initialise a qpms_trans_calculator_t instance for a given convention and truncation degree.
-qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, qpms_normalisation_t nt);
+qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, ///< Truncation degree.
+	       	qpms_normalisation_t nt
+		);
+/// Destructor for qpms_trans_calculator_t.
 void qpms_trans_calculator_free(qpms_trans_calculator *);
 
 complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c,

qpms/vswf.c

@@ -8,6 +8,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include "qpms_error.h"
+#include "normalisation.h"
 
 
 qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() {
@@ -98,35 +99,45 @@ csphvec_t qpms_vswf_single_el(qpms_m_t m, qpms_l_t l, sph_t kdlj,
     qpms_bessel_t btyp, qpms_normalisation_t norm) {
   lmcheck(l,m);
   csphvec_t N;
-  complex double *bessel = malloc((l+1)*sizeof(complex double));
+  complex double *bessel;
-  if(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)) abort();
+  QPMS_CRASHING_MALLOC(bessel,(l+1)*sizeof(complex double));
-  qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, norm);
+  QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel));
-  complex double eimf = cexp(m * kdlj.phi * I);
+  qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, qpms_normalisation_t_csphase(norm));
+
+  complex double eimf = qpms_spharm_azimuthal_part(norm, m, kdlj.phi);
+  complex double d_eimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, kdlj.phi);
   qpms_y_t y = qpms_mn2y(m,l);
 
   N.rc = l*(l+1) * pt.leg[y] * bessel[l] / kdlj.r * eimf;
   complex double besselfac = bessel[l-1] - l * bessel[l] / kdlj.r;
   N.thetac = pt.tau[y] * besselfac * eimf;
-  N.phic = pt.pi[y] * besselfac * I * eimf;
+  N.phic = pt.pi[y] * besselfac * d_eimf_dmf;
+
+  N = csphvec_scale(qpms_normalisation_factor_N_noCS(norm, l, m), N);
 
   qpms_pitau_free(pt);
   free(bessel);
   return N;
 }
+
 csphvec_t qpms_vswf_single_mg(qpms_m_t m, qpms_l_t l, sph_t kdlj,
     qpms_bessel_t btyp, qpms_normalisation_t norm) {
   lmcheck(l,m);
   csphvec_t M;
-  complex double *bessel = malloc((l+1)*sizeof(complex double));
+  complex double *bessel;
-  if(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)) abort();
+  QPMS_CRASHING_MALLOC(bessel,(l+1)*sizeof(complex double));
-  qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, norm);
+  QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel));
-  complex double eimf = cexp(m * kdlj.phi * I);
+  qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, qpms_normalisation_t_csphase(norm));
+  complex double eimf = qpms_spharm_azimuthal_part(norm, m, kdlj.phi);
+  complex double d_eimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, kdlj.phi);
   qpms_y_t y = qpms_mn2y(m,l);
 
   M.rc = 0.;
-  M.thetac = pt.pi[y] * bessel[l] * I * eimf;
+  M.thetac = pt.pi[y] * bessel[l] * d_eimf_dmf;
   M.phic = -pt.tau[y] * bessel[l] * eimf;
 
+  M = csphvec_scale(qpms_normalisation_factor_M_noCS(norm, l, m), M);
+
   qpms_pitau_free(pt);
   free(bessel);
   return M;
@@ -137,10 +148,9 @@ qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
   qpms_vswfset_sph_t *res = malloc(sizeof(qpms_vswfset_sph_t));
   res->lMax = lMax;
   qpms_y_t nelem = qpms_lMax2nelem(lMax);
-  res->el = malloc(sizeof(csphvec_t)*nelem);
+  QPMS_CRASHING_MALLOC(res->el, sizeof(csphvec_t)*nelem);
-  res->mg = malloc(sizeof(csphvec_t)*nelem);
+  QPMS_CRASHING_MALLOC(res->mg, sizeof(csphvec_t)*nelem);
-  if(QPMS_SUCCESS != qpms_vswf_fill(NULL, res->mg, res->el, lMax, kdlj, btyp, norm))
+  QPMS_ENSURE_SUCCESS(qpms_vswf_fill(NULL, res->mg, res->el, lMax, kdlj, btyp, norm));
-    abort(); // or return NULL? or rather assert?
   return res;
 }
 
@@ -163,11 +173,11 @@ csphvec_t qpms_vswf_L00(csph_t kr, qpms_bessel_t btyp,
 
 qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget, 
     csphvec_t * const mgtarget, csphvec_t * const eltarget, qpms_l_t lMax,
-    csph_t kr, qpms_bessel_t btyp, qpms_normalisation_t norm) {
+    csph_t kr, qpms_bessel_t btyp, const 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();
-  qpms_pitau_t pt = qpms_pitau_get(kr.theta, lMax, norm);
+  qpms_pitau_t pt = qpms_pitau_get(kr.theta, lMax, qpms_normalisation_t_csphase(norm));
   complex double const *pbes = bessel + 1; // starting from l = 1
   double const *pleg = pt.leg;
   double const *ppi = pt.pi;
@@ -183,28 +193,32 @@ qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget,
     }
     besderfac = *(pbes-1) - l * besfac;
     for(qpms_m_t m = -l; m <= l; ++m) {
-      complex double eimf = cexp(m * kr.phi * I);
+      complex double eimf = qpms_spharm_azimuthal_part(norm, m, kr.phi);
+      complex double d_eimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, kr.phi);
       if (longtarget) { QPMS_UNTESTED;
-        complex double longfac = sqrt(l*(l+1)) * eimf;
+        double longfac = sqrt(l*(l+1));
         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;
+          * (*pleg) * longfac * eimf;
-        plong->thetac = *ptau * besfac * longfac;
+        plong->thetac = *ptau * besfac * longfac * eimf;
-        plong->phic = *ppi * I * besfac * longfac;
+        plong->phic = *ppi * besfac * longfac * d_eimf_dmf;
+        *plong = csphvec_scale(qpms_normalisation_factor_L_noCS(norm, l, m), *plong);
         ++plong;
       }
       if (eltarget) {
         pel->rc = l*(l+1) * (*pleg) * besfac * eimf;
         pel->thetac = *ptau * besderfac * eimf;
-        pel->phic = *ppi * besderfac * I * eimf;
+        pel->phic = *ppi * besderfac * d_eimf_dmf;
+        *pel = csphvec_scale(qpms_normalisation_factor_N_noCS(norm, l, m), *pel);
         ++pel;
       }
       if (mgtarget) {
         pmg->rc = 0.;
-        pmg->thetac = *ppi * (*pbes) * I * eimf;
+        pmg->thetac = *ppi * (*pbes) * d_eimf_dmf;
         pmg->phic = - *ptau * (*pbes) * eimf;
+        *pmg = csphvec_scale(qpms_normalisation_factor_M_noCS(norm, l, m), *pmg);
         ++pmg;
       }
       ++pleg; ++ppi; ++ptau;
@@ -234,7 +248,9 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
   complex double const *pbes = bessel + 1; // starting from l = 1
 
   qpms_y_t nelem = qpms_lMax2nelem(lMax);
-  csphvec_t * const a1 = malloc(3*nelem*sizeof(csphvec_t)), * const a2 = a1 + nelem, * const a3 = a2 + nelem;
+  csphvec_t *a;
+  QPMS_CRASHING_MALLOC(a, 3*nelem*sizeof(csphvec_t))
+  csphvec_t * const a1 = a, * const a2 = a1 + nelem, * const a3 = a2 + 2 * nelem;
   if(qpms_vecspharm_fill(a1, a2, a3, lMax, kr, norm)) abort();
   const csphvec_t *p1 = a1; 
   const csphvec_t *p2 = a2;
@@ -246,10 +262,12 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
     complex double besderfac = *(pbes-1) - l * besfac;
     double sqrtlfac = sqrt(l*(l+1));
     for(qpms_m_t m = -l; m <= l; ++m) {
-      complex double eimf = cexp(m * kr.phi * I); // FIXME unused variable?!!!
       if (longtarget) {
+        complex double L2Nfac = qpms_normalisation_factor_L_noCS(norm, l, m)
+          / qpms_normalisation_factor_N_noCS(norm, l, m);
         *plong = csphvec_add(csphvec_scale(besderfac-besfac, *p3),
             csphvec_scale(sqrtlfac * besfac, *p2));
+        *plong = csphvec_scale(L2Nfac, *plong);
         ++plong;
       }
       if (eltarget) {
@@ -265,7 +283,7 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
     }
     ++pbes;
   }
-  free(a1);
+  free(a);
   free(bessel);
   return QPMS_SUCCESS;
 }
@@ -273,28 +291,31 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
 qpms_errno_t qpms_vecspharm_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target,
     qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm) {
   assert(lMax >= 1);
-  qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, norm);
+  qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, qpms_normalisation_t_csphase(norm));
   double const *pleg = pt.leg;
   double const *ppi = pt.pi;
   double const *ptau = pt.tau;
   csphvec_t *p1 = a1target, *p2 = a2target, *p3 = a3target;
   for (qpms_l_t l = 1; l <= lMax; ++l) {
     for(qpms_m_t m = -l; m <= l; ++m) {
-      complex double eimf = cexp(m * dir.phi * I);
+      const complex double Mfac = qpms_normalisation_factor_M_noCS(norm, l, m);
+      const complex double Nfac = qpms_normalisation_factor_N_noCS(norm, l, m);
+      const complex double eimf = qpms_spharm_azimuthal_part(norm, m, dir.phi);
+      const complex double deimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, dir.phi);
       if (a1target) {
         p1->rc = 0;
-        p1->thetac = *ppi * I * eimf;
+        p1->thetac = *ppi * deimf_dmf * Mfac;
-        p1->phic = -*ptau * eimf;
+        p1->phic = -*ptau * eimf * Mfac;
         ++p1;
       }
       if (a2target) {
         p2->rc = 0;
-        p2->thetac = *ptau * eimf;
+        p2->thetac = *ptau * eimf * Nfac;
-        p2->phic = *ppi * I * eimf;
+        p2->phic = *ppi * deimf_dmf * Nfac;
         ++p2;
       }
       if (a3target) {
-        p3->rc = sqrt(l*(l+1)) * (*pleg) * eimf;
+        p3->rc = sqrt(l*(l+1)) * (*pleg) * eimf * Nfac;
         p3->thetac = 0;
         p3->phic = 0;
         ++p3;
@@ -308,6 +329,10 @@ qpms_errno_t qpms_vecspharm_fill(csphvec_t *const a1target, csphvec_t *const a2t
 
 qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target,
     qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm) {
+#if 1
+  return qpms_vecspharm_fill(a1target, a2target, a3target, lMax, dir, 
+      qpms_normalisation_dual(norm));
+#else
   assert(lMax >= 1);
   qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, norm);
   double const *pleg = pt.leg;
@@ -341,6 +366,7 @@ qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *cons
   }
   qpms_pitau_free(pt);
   return QPMS_SUCCESS;
+#endif
 }
 
 

qpms/vswf.h

@@ -11,7 +11,8 @@
 #include "qpms_types.h"
 #include <gsl/gsl_sf_legendre.h>
 
-// Methods for qpms_vswf_spec_t
+// ---------------Methods for qpms_vswf_spec_t-----------------------
+//
 /// Creates a qpms_vswf_set_spec_t structure with an empty list of wave indices.
 qpms_vswf_set_spec_t *qpms_vswf_set_spec_init(void);
 /// Appends a VSWF index to a \ref qpms_vswf_set_spec_t, also updating metadata.
@@ -60,6 +61,8 @@ csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec,
 		csph_t kr, ///< Evaluation point.
 		qpms_bessel_t btyp);
 
+// -----------------------------------------------------------------------
+
 /// Electric wave N.
 csphvec_t qpms_vswf_single_el(int m, int n, sph_t kdlj,
 		qpms_bessel_t btyp, qpms_normalisation_t norm);
@@ -87,10 +90,14 @@ qpms_errno_t qpms_legendre_deriv_y_fill(double *where, double *where_deriv, doub
 
 
 /// Evaluate the zeroth-degree longitudinal VSWF \f$ \mathbf{L}_0^0 \f$.
+/**
+ * Any `norm` is being ignored right now.
+ */
 csphvec_t qpms_vswf_L00(
 		csph_t kdrj, //< VSWF evaluation point.
 		qpms_bessel_t btyp,
-		qpms_normalisation_t norm);
+		qpms_normalisation_t norm //< Ignored!
+		);
 
 /// Evaluate VSWFs at a given point from \a l = 1 up to a given degree \a lMax.
 /**