argproc.py: enable manual specification of BaseSpec

Remove old (not valid anymore) comments from vswf.h

Doxygenise translations_dbg.h

BUGS.rst

@@ -23,28 +23,6 @@ Plane wave decompositions gives wrong value on the longitudinal part.
 The implementation of the L coefficients OR the longitudinal waves
 is thus probably wrong.
 
-Scattering result asymmetry
----------------------------
-The lattice scattering code (such as finitesqlatzsym-scattery.py) produces
-asymmetric results where one should not get them, due to the system symmetry.
-
-It seems that the asymmetry appears mostly in the y-direction (i.e.
-for example the scattering/absorption cross section at k = (kx, ky, kz)
-is not exactly the same as k = (kx, -ky, kz).
-
-What has been checked (hopefully):
- - The flip operators for electric waves
- - Some weird kind of rounding or other numerical error depending on
-   the position ordering of the matrix (randomized indices give
-   the same asymmetry).
-
-What has not been checked (so well):
- - The x,y,z-flip operators for magnetic waves (i.e. are they really 
-   supposet to bejust the
-   same as the operators for electric waves, just with opposite sign?) 
- - zplane_pq_y
- - the translation operators
-
 
 Singular value asymmetry
 ------------------------
@@ -60,6 +38,3 @@ Moreover, the non-normalized legendre functions that are used in translations.c
 are likely to overflow for even smaller values of l. 
 Therefore TODO: modify translations.c in a way that uses normalised legendre functions everywhere.
 
-Memory management
------------------
-When compiled with optimizations, I get stack smashing error in  qpms_trans_calculator_get_2DFT_longrange_AB_arrays().

CMakeLists.txt

@@ -48,6 +48,7 @@ add_subdirectory(faddeeva)
 
 add_subdirectory (qpms)
 
+add_subdirectory (tests/catch EXCLUDE_FROM_ALL)
 
 enable_testing()
 add_subdirectory (tests EXCLUDE_FROM_ALL)

examples/rectangular/modes/01_compute_modes.sh

@@ -9,10 +9,10 @@ ${MISCDIR}/lat2d_modes.py \
 	-b s${A1X_nm}e-9 s${A1Y_nm}e-9 \
 	-b s${A2X_nm}e-9 s${A2Y_nm}e-9 \
 	-p s${P1X_nm}e-9 s${P1Y_nm}e-9 \
-	-p s${P2X_nm}e-9 s${P2Y_nm}e-9 \
 	-L 3 -m $METAL -r ${RADIUS_nm}e-9 -H ${HEIGHT_nm}e-9 \
-	-k s${KPOINTX_nmi}e9 s${KPOINTY_nmi}e9 \
+	-k 0 0 \
 	-d -3 \
 	-t 0.01 \
 	-c 250 \
 	-P
+	#-k s${B1X_nmi}e9 s${B1Y_nmi}e9 \

misc/finiterectlat-constant-driving.py

@@ -264,6 +264,7 @@ np.savez(outfile, meta=vars(a),
                  #ccd_size = ccd_size if not math.isnan(a.ccd_distance) else None,
                  ccd_points = ccd_points if not math.isnan(a.ccd_distance) else None,
                  ccd_fields = ccd_fields if not math.isnan(a.ccd_distance) else None,
+                 fullvec_poffsets = ss.fullvec_poffsets,
        )
 logging.info("Saved to %s" % outfile)
 

notes/VSWF_from_SSWF.md

@@ -0,0 +1,33 @@
+VSWF expansions in terms of SSWF
+================================
+
+From
+\cite necada_multiple-scattering_2021, eq. (2.19)
+\f[
+	\wfkcout_{\tau lm}\left(\kappa (\vect r - \vect r_1) \right) =
+		\sum_{\tau'l'm'} \tropSr{\kappa(\vect r_2 - \vect r_1)}_{\tau l m;\tau'l'm} \wfkcreg_{\tau'l'm'}(\vect r -\vect r_2),
+	\qquad |\vect r -\vect r_2| < |\vect r_1 - \vect r_2|,
+\f]
+setting \f$ \vect r = \vect r_2\f$ and considering that 
+ \f$ \wfkcreg_{\tau'l'm'}(\vect 0) \ne \vect 0 \f$ only for electric dipole waves (\f$ \tau = \mathrm{E}, l=1 \f$),
+ we have
+\f[
+	\wfkcout_{\tau lm}\left(\kappa (\vect r - \vect r_1) \right) =
+		\sum_{m'} \tropSr{\kappa(\vect r - \vect r_1)}_{\tau l m;\mathrm{E}1m} \wfkcreg_{\mathrm{E}1m'}(\vect 0),
+	\qquad \vect r \ne \vect r_1 .
+\f]
+Combining this with
+\cite necada_multiple-scattering_2021, eq. (2.25)
+\f[
+	\tropSr{\vect d}_{\tau l m; \tau' l' m'} = \sum_{\lambda =|l-l'|+|\tau-\tau'|}^{l+l'} 
+		C^{\lambda}_{\tau l m;\tau' l'm'} \underbrace{ \spharm{\lambda}{m-m'}(\uvec d) h_\lambda^{(1)}(d)}_{\sswfout_\lambda^{m-m'}(\vect d)},
+\f]
+we get 
+\f[
+	\wfkcout_{\tau lm}(\vect d) = \sum_{m'=-1}^1 \wfkcreg_{\mathrm{E}1m'}(\vect 0) 
+		\sum_{\lambda=l-1+|\tau-\tau'|}^{l+1} 
+		C^\lambda_{\tau l m;\mathrm{E}1m'} \sswfout_\lambda^{m-m'}(\vect d).
+\f]
+Note that the VSWF components in this expression are given in global "cartesian" basis, 
+*not* the local orthonormal basis derived from spherical coordinates.
+(This is mostly desirable, because in lattices we need to work with flat coordinates anyway.)

notes/conventions.md

@@ -1,6 +1,10 @@
 VSWF conventions {#vswf_conventions}
 ====================================
 
+*This page provides reference about the VSWF conventions used in the literature.
+For VSWF convention specification in QPMS API, see 
+[SWF Conventions in QPMS](@ref swf_conventions_qpms).*
+
 In general, the (transversal) VSWFs can be defined using (some) vector spherical harmonics
 as follows: \f[
 	\wfm\pr{k\vect r}_{lm} = \sphbes_l(kr) \vshrot_{lm} (\uvec r),\\
@@ -50,7 +54,7 @@ where the connection to negative orders is
 	\dlmfFer{\nu}{m}(x) = (-1)^m \frac{\Gamma\pr{\nu-m+1}}{\Gamma\pr{\nu+m+1}}\dlmfFer{\nu}{m}(x),\\
 	%\dlmfLeg{\nu}{m}(x) =        \frac{\Gamma\pr{\nu-m+1}}{\Gamma\pr{\nu+m+1}}\dlmfLeg{\nu}{m}(x).\\
 \f]
-Note that there are called "Ferrers" functions in DLMF, while the "Legendre" functions have slightly
+Note that they are called "Ferrers" functions in DLMF, while the "Legendre" functions have slightly
 different meaning / conventions (Ferrers functions being defined for \f$ \abs{x} \le 1 \f$, whereas
 Legendre for \f$ \abs{x} \ge 1 \f$. We will not use the DLMF "Legendre" functions here.
 

notes/conventions_qpms.md

@@ -0,0 +1,114 @@
+SWF conventions in QPMS {#swf_conventions_qpms}
+=================================================
+
+*This page describes how (V)SWF conventions are specified
+internally and in QPMS API. For a general overview of VSWF
+conventions in the literature, see [VSWF Conventions](@ref vswf_conventions).*
+
+Convention enumerator
+---------------------
+
+Most of the meaningful phase and normalisation conventions for spherical waves
+can be specified by the enum type @ref qpms_normalisation_t.
+
+The type can be also used to specify conventions that are currently not fully
+supported in QPMS (such as those based on real spherical harmonics).
+
+As an enum type, it does not cover all the conventions possibly imaginable,
+but it does cover the most meaningful ones and most of those that can be found
+in the literature. 
+
+(Most notably, it does not cover the “anti-normalisation”
+that does appear in certain expressions in some literature where the spherical
+harmonics contain unnormalised Legendre functions, so that the basis set of 
+of spherical harmonics has different norms for different signs of *m*
+for the same *l*. This is a bad idea overall and an absolutely atrocious
+approach for numerics. Do not use that.)
+
+
+VSWF evaluation
+---------------
+
+Evaluation of VSWFs using qpms_vswf_fill(),  qpms_eval_vswf(),
+qpms_uvswf_fill() and other functions from vswf.h are evaluated as follows.
+These fuctions take a @ref qpms_normalisation_t as an argument.
+The Ferrers-Legendre functions and the π, τ functions are evaluated 
+by qpms_pitau_get(), which internally uses gsl_sf_legendre_deriv_array_e(). 
+
+Note only the information about the Condon-Shortley
+phase is passed to qpms_pitau_get() – the result of this function
+uses always the GSL_SF_LEGENDRE_SPHARM normalisation,
+and possible normalisation and other phase factors are evaluated afterwards
+using the inline functions
+qpms_normalisation_factor_L_noCS(),
+qpms_normalisation_factor_M_noCS(),
+qpms_normalisation_factor_N_noCS().
+
+Evaluation of vector spherical harmonics only with qpms_vecspharm_fill()
+works similarly but TODO.
+
+TODO reference to pi, tau.
+
+VSWF translation operator evaluation
+------------------------------------
+
+In practice, translation operators are calculated by first creating
+an instance of the qpms_trans_calculator structure, which contains 
+a table of constant normalisation factors for a given phase/normalisation
+convention (it is assumed that the phase/normalisation conventions do not
+change with the translation), and then calling 
+qpms_trans_calculator_get_AB_arrays()
+(or others).
+
+The precomputed factor table in qpms_trans_calculator_t contains a CS phase
+related factor (via qpms_trans_normfac()). 
+
+Function qpms_trans_calculator_get_AB_arrays_buf() then calculates
+the unnormalised (GSL_SF_LEGENDRE_NONE) 
+associated Legendre functions always with CS phase -1;
+and the A, B arrays are filled (via qpms_trans_calculator_get_AB_arrays_precalcbuf())
+by individual calls of qpms_trans_calculator_get_A_precalcbuf()
+and qpms_trans_calculator_B_precalcbuf(), which basically just multiply
+and sum the precalculated constant factors with the radial (Bessel),
+polar (Legendre) and azimuthal (exponential/trigonometric) functions.
+**This means that the normalisation and phase convention is fully embedded
+in the constant factor tables, and nothing is calculated during "runtime".**
+
+The "higher-level" qpms_trans_calculator_get_trans_array() currently
+just calls qpms_trans_calculator_get_AB_arrays() and then reorders
+the elements (using qpms_trans_array_from_AB()), asserting that
+the normalisation conventions remain the same.
+
+
+There seems to be an inconsistency between 
+qpms_trans_calculator_get_B_buf() and 
+qpms_trans_calculator_get_A_buf() on one hand, and
+qpms_trans_calculator_get_AB_buf_p() and
+qpms_trans_calculator_get_AB_arrays_buf() on the other.
+While the latter two functions use always -1 as the CS phase,
+the former two take it from the normalisation enumerator.
+**Although the former two are probably used nowhere in the production,
+this needs to be fixed.**
+
+
+
+Lattice sums
+------------
+
+### Scalar SWFs
+
+### Translation operators
+
+Function qpms_trans_calculator_get_AB_arrays_e32_e()
+first compute the scalar lattice sums (using ewald3_sigma_short(),
+ewald3_sigma_long() and ewald3_sigma0() calls).
+
+These are then transformed into the VSWF translation operator
+elements in a similar manner as in 
+qpms_trans_calculator_get_A_precalcbuf() and
+qpms_trans_calculator_get_B_precalcbuf(), although there some optical
+differences (CHECK!).
+
+### VSWFs
+
+

notes/mathjax_newcommands.js

@@ -31,8 +31,15 @@ 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
+	    tropS: "{\\mathcal{S}}", // Translation operator singular
+	    tropR: "{\\mathcal{R}}", // Translation operator regular
+	    tropSr: ["{{\\mathcal{S}^\\mathrm{#1}}\\pr{{#2}}}", 2, ""], // Translation operator singular
+	    tropRr: ["{{\\mathcal{R}^\\mathrm{#1}}\\pr{{#2}}}", 2, ""], // Translation operator regular
+	    tropSrr: ["{{\\mathcal{S}^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator singular
+	    tropRrr: ["{{\\mathcal{R}^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator regular
+            sswfout: "{\\psi}", // outgoing SSWF
+            sswfreg: "{\\phi}", // regular SSWF
+
 
             // Kristensson's VSWFs, complex version (2014 notes)
 	    wfkc: "{\\vect{y}}", // any wave

qpms/CMakeLists.txt

@@ -16,6 +16,7 @@ add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c e
 	lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
 	bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c
 	lll.c beyn.c trivialgroup.c
+	translations_dbg.c scatsystem_dbg.c
 	)
 use_c99()
 

qpms/__init__.py

@@ -19,7 +19,7 @@ except ImportError as ex:
 from .cyquaternions import CQuat, IRot3
 from .cybspec import VSWFNorm, BaseSpec, default_bspec
 from .cytmatrices import CTMatrix, TMatrixInterpolator, TMatrixGenerator
-from .cytranslations import trans_calculator
+#from .cytranslations import trans_calculator
 from .cymaterials import MaterialInterpolator, EpsMu, LorentzDrudeModel, lorentz_drude, EpsMuGenerator
 from .cycommon import dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType, VSWFType
 from .cywaves import vswf_single

qpms/analysis.py

@@ -0,0 +1,150 @@
+'''
+Auxiliary functions to analyse results obtained using scripts in
+the ../misc directory.
+
+This needs to be kept in sync with the scripts and qpms/argproc.py.
+'''
+import math
+import numpy as np
+import qpms
+from qpms.qpms_c import Particle
+from qpms.cytmatrices import CTMatrix, TMatrixGenerator
+from qpms.cybspec import BaseSpec
+from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude
+from qpms.symmetries import point_group_info
+from qpms import FinitePointGroup, ScatteringSystem, eV, hbar
+eh = eV / hbar
+
+def cleanarray(a, atol=1e-8, copy=True, fullNaN=True):
+    a = np.array(a, copy=copy)
+    sieve = abs(a.real) < atol
+    a[sieve] = 1j * a[sieve].imag
+    sieve = abs(a.imag) < atol
+    a[sieve] = a[sieve].real
+    if fullNaN:
+        if np.all(a == 0):
+            a[...] = 0
+    return a
+
+def nicerot(a, atol=1e-10, copy=True): 
+    '''
+    Gives an array a "nice" phase.
+    '''
+    a = np.array(a, copy=copy)
+    i = np.argmax(abs(a))
+    a = a / a[i] * abs(a[i])
+    return a
+
+def _meta_matspec2emg(matspec):
+    if isinstance(matspec, (float, complex)): # number is interpreted as refractive index
+        return EpsMuGenerator(EpsMu(matspec**2)) 
+    else:
+        return EpsMuGenerator(lorentz_drude[matspec])
+
+
+class FiniteRectLatAnalysis:
+    '''Recreates the scattering system structure from finiterectlat-modes.py or 
+    finiterectlat-constant-driving.py'''
+    def __init__(self, data):
+        meta = data['meta'][()]
+        scatter = 'symmetry_adapted' in meta.keys()# finiterectlat-constant-driving ; TODO unify
+        if scatter: 
+            #thegroup = meta['symmetry_adapted'] # always D2h, this option does something else!
+            thegroup = 'D2h'
+        else: # finiterectlat-modes
+            thegroup = 'D2h' if meta['D2'] else 'D4h'
+        Nx, Ny = meta['size']
+        px, py = meta['period']
+        orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
+        orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
+        orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
+        bspec = BaseSpec(lMax = meta['lMax'])
+        nelem = len(bspec)
+        bg = _meta_matspec2emg(meta['background'])
+        fg = _meta_matspec2emg(meta['material'])
+        if meta.get('height', None) is None:
+            tmgen = TMatrixGenerator.sphere(bg, fg, meta['radius'])
+        else:
+            tmgen = TMatrixGenerator.cylinder(bg, fg, meta['radius'], meta['height'], lMax_extend=meta['lMax_extend'])
+        particles = [Particle(orig_xy[i], tmgen, bspec) for i in np.ndindex(orig_xy.shape[:-1])]
+        sym = FinitePointGroup(point_group_info[thegroup])
+        omega = data['omega'][()] if scatter else meta['centre'] * eh 
+        self.ss, self.ssw = ScatteringSystem.create(particles, bg, omega, sym=sym)
+        ss1, ssw1 = ScatteringSystem.create([Particle([0,0,0], tmgen, bspec)], bg, omega, sym=sym)
+        self.ss1, self.ssw1 = ss1, ssw1
+        
+        # per-particle projectors/transformation to symmetry-adapted bases
+        fvcs1 = np.empty((nelem, nelem), dtype=complex)
+        y = 0
+        iris1 = []
+        for iri1 in range(ss1.nirreps):
+            for j in range(ss1.saecv_sizes[iri1]):
+                pvc1 = np.zeros((ss1.saecv_sizes[iri1],), dtype=complex)
+                pvc1[j] = 1
+                fvcs1[y] = ss1.unpack_vector(pvc1, iri1)
+                fvcs1[y] = cleanarray(nicerot(fvcs1[y], copy=False), copy=False)
+                y += 1
+                iris1.append(iri1)
+        
+        # Mapping between ss particles and grid positions
+        positions = self.ss.positions.astype(np.float32) 
+        # round to get rid of duplicities due to rounding errors
+        positions = positions.round(7-int(math.log(np.amax(abs(positions)),10)))
+        xpositions = np.unique(positions[:,0]) 
+        ypositions = np.unique(positions[:,1])
+        assert(len(xpositions) == Nx)
+        assert(len(ypositions == Ny))
+        # particle positions as integer indices
+        posmap = np.empty((positions.shape[0],2), dtype=int)
+        #invposmap = np.empty((Nx, Ny), dtype=int)
+        for i, pos in enumerate(positions):
+            posmap[i,0] = np.searchsorted(xpositions, positions[i,0])
+            posmap[i,1] = np.searchsorted(ypositions, positions[i,1])
+            #invposmap[posmap[i,0], posmap[i, 1]] = i
+
+        self.meta = meta
+        self.npzfile = data
+        self.fvcs1 = fvcs1
+        self.iris1 = np.array(iris1)
+        self.bspec = bspec
+        self.nelem = nelem
+        self.Nx, self.Ny, self.px, self.py = Nx, Ny, px, py
+        self.posmap = posmap
+        self.fullvec_poffsets = nelem * np.arange(Nx*Ny)
+        self.group = thegroup
+        if scatter:
+            self.typ = 'scatter'
+            self.fvcs1_fromdata = data['fvcs1']
+            self.iris1_fromdata = data['iris1']
+            self.positions_fromdata = data['positions']
+            self.driving_symmetry_adapted = data['symmetry_adapted']
+        else: # modes
+            self.typ = 'modes'
+            self.eigval = data['eigval']
+            self.neig = len(self.eigval)
+            self.eigvec = data['eigvec']
+            self.residuals = data['residuals']
+            self.ranktest_SV = data['ranktest_SV']
+            self.iri = data['iri'][()]
+
+
+    def fullvec2grid(self, fullvec, swapxy=False):
+        arr = np.empty((self.Nx, self.Ny, self.nelem), dtype=complex)
+        for pi, offset in enumerate(self.fullvec_poffsets):
+            ix, iy = self.posmap[pi]
+            arr[ix, iy] = fullvec[offset:offset+self.nelem]
+        return np.swapaxes(arr, 0, 1) if swapxy else arr
+    
+    def saecv2grid(self, saecv, swapxy=False):
+        fullvec = self.ss.unpack_vector(saecv, iri=self.iri)
+        return self.fullvec2grid(fullvec, swapxy)
+    
+    def eigvec_asgrid(self, i, swapxy = False):
+        '''
+        Returns i-th eigenmode as a (Nx, Ny, nelem)-shaped grid of excitation
+        coefficients.
+        '''
+        if self.iri is not None:
+            return self.saecv2grid(self.eigvec[i], swapxy)
+        else:
+            return self.fullvec2grid(self.eigvec[i], swapxy)

qpms/argproc.py

@@ -5,6 +5,7 @@ Common snippets for argument processing in command line scripts.
 import argparse
 import sys
 import warnings
+import ast
 
 def flatten(S):
     if S == []:
@@ -192,6 +193,13 @@ def material_spec(string):
                 raise argparse.ArgumentTypeError("Material specification must be a supported material name %s, or a number" % (str(lorentz_drude.keys()),)) from ve
         return lemat
 
+def string2bspec(string):
+    """
+    Converts string representation of list to BaseSpec.
+    """
+    from .cybspec import BaseSpec
+    return BaseSpec(ast.literal_eval(string))
+
 class ArgParser:
     ''' Common argument parsing engine for QPMS python CLI scripts. '''
     
@@ -248,6 +256,13 @@ class ArgParser:
         mpgrp.add_argument("+H", "++height", nargs=2,
                 action=make_dict_action(argtype=float, postaction='store', first_is_key=True,),
                 help='particle radius (cylinder; labeled)')
+        # Alternatively, add a constant T-matrix
+        mpgrp.add_argmuent("-w", "--vswf-set", nargs=1, default={},
+                action=make_dict_action(argtype=string2basespec, postaction='store', first_is_key=False),
+                help='Manual specification of VSWF set codes (format as a python list of integers); see docs on qpms_uvswfi_t for valid codes or simply use --lMax instead. Overrides --lMax.')
+        mpgrp.add_argmuent("+w", "++vswf-set", nargs=2, default={},
+                action=make_dict_action(argtype=string2basespec, postaction='store', first_is_key=True),
+                help='Manual specification of VSWF set codes (format as a python list of integers); see docs on qpms_uvswfi_t for valid codes or simply use ++lMax instead. Overrides ++lMax and --lMax.')
 
     atomic_arguments = {
             'rectlattice2d_periods': lambda ap: ap.add_argument("-p", "--period", type=float, nargs='+', required=True, help='square/rectangular lattice periods', metavar=('px','[py]')),
@@ -337,6 +352,14 @@ class ArgParser:
             return tmgen
 
     def _add_bspec(self, key):
+        """
+        Transforms a given key into a BaseSpec and registers
+        the BaseSpec with the key.
+        Always returns a BaseSpec (unless an exception occurs).
+        If an equivalent instance of BaseSpec is registered, returns that one
+        (therefore, all equivalent BaseSpecs in the register values should
+        be the same instance).
+        """
         if key in self._bspec_register.keys():
             return self._bspec_register[key]
         else:
@@ -345,7 +368,9 @@ class ArgParser:
                 bspec = key
             elif isinstance(key, int):
                 bspec = self._add_bspec(BaseSpec(lMax=key))
-            else: raise TypeError("Can't register this as a BaseSpec")
+            #else: raise TypeError("Can't register this as a BaseSpec")
+            else: 
+                bspec = self._add_bspec(BaseSpec(key))
             self._bspec_register[key] = bspec
             return bspec
 
@@ -488,8 +513,12 @@ class ArgParser:
             warnings.warn("No particle position (-p or +p) specified, assuming single particle in the origin / single particle per unit cell!")
             a.position[None] = [(0.,0.,0.)]
         for poslabel in a.position.keys():
+            # TODO HERE GOES THE CODE TRYING TO LOAD CONSTANT T-MATRIX
             try:
-                lMax = a.lMax.get(poslabel, False) or a.lMax[None]
+                lMax_or_bspec = ( a.vswf_set.get(poslabel, False) 
+                        or a.lMax.get(poslabel, False)
+                        or a.vswf_set.get(None, False)
+                        or a.lMax[None] )
                 radius = a.radius.get(poslabel, False) or a.radius[None]
                 # Height is "inherited" only together with radius
                 height = a.height.get(poslabel, None) if poslabel in a.radius.keys() else a.height.get(None, None)
@@ -508,7 +537,7 @@ class ArgParser:
             tmspec = (a.background, material, (radius, height, lMax_extend)) 
             self.tmspecs[poslabel] = tmspec
             self.tmgens[poslabel] = self._add_tmg(tmspec)
-            self.bspecs[poslabel] = self._add_bspec(lMax)
+            self.bspecs[poslabel] = self._add_bspec(lMax_or_bspec)
             poslist_cured = []
             for pos in a.position[poslabel]:
                 if len(pos) == 1:

qpms/beyn.h

@@ -3,18 +3,20 @@
  */
 #ifndef BEYN_H
 #define BEYN_H
+#ifdef __cplusplus
+extern "C" {
+#endif
 
 #include <stddef.h>
-#include <complex.h>
 
 /// User-supplied function that provides the (row-major) m × m matrix M(z) whose "roots" are to be found.
 /** Pure C array version */
-typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex double z, void *params);
+typedef int (*beyn_function_M_t)(_Complex double *target_M, size_t m, _Complex double z, void *params);
 
 /// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$.
 /** Pure C array version */
-typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l,
-	       	const complex double *Vhat, complex double z, void *params); 
+typedef int (*beyn_function_M_inv_Vhat_t)(_Complex double *target_M_inv_Vhat, size_t m, size_t l,
+	       	const _Complex double *Vhat, _Complex double z, void *params); 
 
 /// Complex plane integration contour structure.
 typedef struct beyn_contour_t {
@@ -26,15 +28,15 @@ typedef struct beyn_contour_t {
 	 * It does not have to be a centre in some strictly defined sense,
 	 * but it should be "somewhere around" where the contour is.
 	 */
-	complex double centre;
+	_Complex double centre;
 	/// Function testing that a point \a z lies inside the contour (optional).
-	_Bool (*inside_test)(struct beyn_contour_t *, complex double z);
-	complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points.
+	_Bool (*inside_test)(struct beyn_contour_t *, _Complex double z);
+	_Complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points.
 } beyn_contour_t;
 
 /// Complex plane elliptic integration contour with axes parallel to the real, imaginary axes.
 /** Free using free(). */
-beyn_contour_t *beyn_contour_ellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints);
+beyn_contour_t *beyn_contour_ellipse(_Complex double centre, double halfax_re, double halfax_im, size_t npoints);
 
 
 typedef enum {
@@ -47,23 +49,23 @@ typedef enum {
 
 /// Complex plane "half-elliptic" integration contour with axes parallel to the real, imaginary axes.
 /** Free using free(). */
-beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints,
-		beyn_contour_halfellipse_orientation or);
+beyn_contour_t *beyn_contour_halfellipse(_Complex double centre, double halfax_re, double halfax_im, size_t npoints,
+		beyn_contour_halfellipse_orientation o);
 
 /// Similar to halfellipse but with rounded corners.
-beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, double halfax_im, 
+beyn_contour_t *beyn_contour_kidney(_Complex double centre, double halfax_re, double halfax_im, 
 	double rounding, ///< Must be in interval [0, 0.5)
-	size_t n, beyn_contour_halfellipse_orientation or);
+	size_t n, beyn_contour_halfellipse_orientation o);
 
 
 /// Beyn algorithm result structure (pure C array version).
 typedef struct beyn_result_t {
 	size_t neig; ///< Number of eigenvalues found.
 	size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec).
-	complex double *eigval;
-	complex double *eigval_err;
+	_Complex double *eigval;
+	_Complex double *eigval_err;
 	double *residuals;
-	complex double *eigvec; // Rows are the eigenvectors
+	_Complex double *eigvec; // Rows are the eigenvectors
 	double *ranktest_SV;
 
 } beyn_result_t;
@@ -83,4 +85,7 @@ beyn_result_t *beyn_solve(
 		double res_tol ///< (default: `0.0`) TODO DOC.
 	      );
 
+#ifdef __cplusplus
+}
+#endif
 #endif // BEYN_H

qpms/cybspec.pyx

@@ -1,4 +1,5 @@
 import numpy as np
+cimport numpy as np
 import enum
 from .cycommon import get_mn_y, tlm2uvswfi
 __VSWF_norm_dict = {
@@ -35,6 +36,13 @@ cdef class BaseSpec:
     whenever used in other wrapper classes that need the pointer
     to qpms_vswf_set_spec_t, remember to set a (private, probably immutable) reference to qpms.basespec to ensure
     correct reference counting and garbage collection.
+
+    BaseSpec can be constructed either using a single keyword argument `lMax` to include
+    all the transversal vector spherical wavefunction to a specific degree, or
+    using an explicit list of integer basis VSWF codes, \see qpms_uvswfi_t
+    BaseSpec(lMax = <positive_integer>)
+        or
+    BaseSpec([ uvswfi0, ... ])
     '''
     #cdef qpms_vswf_set_spec_t s # in pxd
     #cdef np.ndarray __ilist     # in pxd
@@ -132,6 +140,14 @@ cdef class BaseSpec:
         # TODO raise correct errors (TypeError on bad type of key, IndexError on exceeding index)
         return self.__ilist[key]
 
+    def __repr__(self):
+        cdef qpms_l_t lMax_candidate
+        lMax_candidate = qpms_nelem2lMax(len(self.__ilist) // 2)
+        if lMax_candidate >= 0:
+            if np.array_equal(self.__ilist, BaseSpec(lMax=lMax_candidate).ilist):
+                return 'BaseSpec(lMax=%d)' % lMax_candidate
+        return 'BaseSpec([' + ','.join('%d' % i  for i in self.__ilist) + '])'
+
     property ilist:
         def __get__(self):
             return self.__ilist

qpms/cycommon.pyx

@@ -20,6 +20,12 @@ class BesselType(enum.IntEnum):
     HANKEL_PLUS = QPMS_HANKEL_PLUS
     HANKEL_MINUS = QPMS_HANKEL_MINUS
 
+class EwaldPart(enum.IntEnum):
+    LONG_RANGE = QPMS_EWALD_LONG_RANGE
+    SHORT_RANGE = QPMS_EWALD_SHORT_RANGE
+    FULL = QPMS_EWALD_FULL
+    ZEROTERM = QPMS_EWALD_0TERM
+
 class PointGroupClass(enum.IntEnum):
     CN    = QPMS_PGS_CN
     S2N   = QPMS_PGS_S2N
@@ -138,6 +144,46 @@ def get_y_mn_unsigned(int nmax):
             i = i + 1
     return(ymn_plus, ymn_minus)
 
+def get_nelem_scalar(lMax):
+    # TODO DOC, exceptions etc.
+    #return qpms_lMax2nelem(lMax)
+    return lMax * (lMax + 2) + 1 # = (lMax + 1)**2
+
+## Auxillary function for retrieving the "meshgrid-like" indices; inc. nmax
+@cython.boundscheck(False)
+def get_mn_y_scalar(int nmax):
+    """
+    Auxillary function for retreiving the 'meshgrid-like' indices from the flat indexing; 
+    inc. nmax.
+    ('y to mn' conversion)
+    
+    Parameters
+    ----------
+
+    nmax : int
+        The maximum order to which the SSWFs / Legendre functions etc. will be evaluated.
+        
+    Returns
+    -------
+    
+    output : (m, n)
+        Tuple of two arrays of type np.array(shape=(nmax*nmax + 2*nmax + 1), dtype=np.int),
+        where [(m[y],n[y]) for y in range(nmax*nmax + 2*nmax + 1)] covers all possible 
+        integer pairs n >= 0, -n <= m <= n.
+    """
+    cdef Py_ssize_t nelems = (nmax + 1)**2
+    cdef np.ndarray[np.int_t,ndim=1] m_arr = np.empty([nelems], dtype=np.int)
+    cdef np.ndarray[np.int_t,ndim=1] n_arr = np.empty([nelems], dtype=np.int)
+    cdef Py_ssize_t i = 0
+    cdef np.int_t n, m
+    for n in range(0,nmax+1):
+        for m in range(-n,n+1):
+            m_arr[i] = m
+            n_arr[i] = n
+            i = i + 1
+    return (m_arr, n_arr)
+
+
 
 def tlm2uvswfi(t, l, m):
     ''' TODO doc

qpms/cyewaldtest.pyx

@@ -1,12 +1,36 @@
 from .qpms_cdefs cimport *
+from .qpms_c cimport ScatteringSystem, _ScatteringSystemAtOmega, _ScatteringSystemAtOmegaK
+from .cycommon import EwaldPart, BesselType
 from libc.stdlib cimport malloc, free, calloc
 import numpy as np
+from cython.parallel cimport prange, parallel
+from cython cimport boundscheck, wraparound
 
 cdef extern from "ewald.h":
     void ewald3_2_sigma_long_Delta(cdouble *target, double *err,  int maxn, cdouble x, int xbranch, cdouble z)
     void ewald3_2_sigma_long_Delta_series(cdouble *target, double *err,  int maxn, cdouble x, int xbranch, cdouble z)
     void ewald3_2_sigma_long_Delta_recurrent(cdouble *target, double *err,  int maxn, cdouble x, int xbranch, cdouble z, bint bigimz)
     int complex_gamma_inc_e(double a, cdouble x, int xbranch, qpms_csf_result *result)
+    extern int ewald_factor_ipow_l
+    extern int ewald_factor_ipow_m
+
+cdef extern from "translations_dbg.h":
+    int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
+        cdouble *dest, const csph_t kdlj, const  qpms_bessel_t J) nogil
+    int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c,
+        cdouble *dest, const cdouble k, const cart3_t r, const qpms_bessel_t J) nogil
+    int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
+        cdouble * const sigmas_total, double * serr_total,
+        const double eta, const cdouble k,
+        const cart2_t b1, const cart2_t b2,
+        const cart2_t beta,
+        const cart3_t particle_shift,
+        double maxR, double maxK,
+        const qpms_ewald_part parts) nogil
+
+cdef extern from "scatsystem_dbg.h":
+    int qpms_scatsyswk_test_sswf_basis_e(cdouble *target, const qpms_scatsys_at_omega_k_t *sswk,
+            qpms_bessel_t typ, cart3_t evalpoint, qpms_ewald_part parts) nogil
 
 def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, method='auto'):
     cdef np.ndarray[double, ndim=1] err_np
@@ -29,9 +53,117 @@ def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, met
     else:
         return target_np
 
+def scatsyswk_sswf_basis(_ScatteringSystemAtOmegaK sswk, evalpos, btyp=QPMS_HANKEL_PLUS, ewaldparts=EwaldPart.FULL):
+    """Evaluate Ewald-summed scalar SWFs that are used to compose the tranlation operators.
+
+    Parameters
+    ----------
+    evalpos: array_like of floats and shape (..., 3)
+        Evaluation points in cartesian coordinates.
+
+    Returns
+    -------
+    ndarray of complex, with the shape `evalpos.shape[:-1] + (particle_count, pnelem)`
+    """        
+    if(btyp != QPMS_HANKEL_PLUS):
+        raise NotImplementedError("Only first kind Bessel function-based fields are supported")
+    cdef qpms_bessel_t btyp_c = BesselType(btyp)
+    cdef qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
+    evalpos = np.array(evalpos, dtype=float, copy=False)
+    if evalpos.shape[-1] != 3:
+        raise ValueError("Last dimension of evalpos has to be 3")
+
+    cdef ScatteringSystem ss = sswk.ssw_pyref.ss_pyref
+    cdef qpms_scatsys_t *ss_c = ss.s
+    cdef qpms_scatsys_at_omega_k_t *sswk_c = sswk.rawpointer()
+    cdef const qpms_trans_calculator *c = ss_c[0].c
+    cdef qpms_l_t sswf_lMax = 2*c[0].lMax+1
+    cdef qpms_y_t pnelem = qpms_lMax2nelem_sc(sswf_lMax)
+    cdef qpms_ss_pi_t particle_count = sswk.ssw_pyref.ss_pyref.particle_count
+
+    cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
+    cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], particle_count, pnelem), dtype=complex)
+    cdef cdouble *res
+    cdef cart3_t pos
+    cdef Py_ssize_t i, pi, j
+    with nogil, wraparound(False), parallel():
+        res = <cdouble *> malloc(particle_count * pnelem * sizeof(cdouble))
+        for i in prange(evalpos_a.shape[0]):
+            pos.x = evalpos_a[i,0]
+            pos.y = evalpos_a[i,1]
+            pos.z = evalpos_a[i,2]
+            qpms_scatsyswk_test_sswf_basis_e(res, sswk_c, btyp_c, pos, ewaldparts_c)
+            for pi in range(particle_count):
+                for j in range(pnelem):
+                    results[i,pi,j] = res[pi * pnelem + j]
+        free(res)
+    return results.reshape(evalpos.shape[:-1] + (particle_count, pnelem))
+
+def scatsysw_eval_sswf(_ScatteringSystemAtOmega ssw, evalpos, btyp=QPMS_HANKEL_PLUS):
+    """Evaluate sswfs by internal qpms_trans_calculator at given locations.
+
+    Parameters
+    ----------
+    ssw: _ScatteringSystemAtOmega
+        We take the wavenumber and qpms_trans_calculator from this, nothing else.
+        (In a way, this is an overkill, but there is no up-to-date standalone cython
+        wrapper for qpms_trans_calculator.)
+    evalpos: array_like of floats and shape (..., 3)
+        Evaluation points in cartesian coordinates.
+    btyp: BesselType
+        Kind of Bessel functions to be used
+
+    Returns
+    -------
+    ndarray of complex, with the shape `evalpos.shape[:-1] + (pnelem,)`
+    where `pnelem` is TODO
+    """
+    cdef qpms_bessel_t btyp_c = BesselType(btyp)
+    evalpos = np.array(evalpos, dtype=float, copy=False)
+    if evalpos.shape[-1] != 3:
+        raise ValueError("Last dimension of evalpos has to be 3")
+
+    cdef ScatteringSystem ss = ssw.ss_pyref
+    cdef qpms_scatsys_t *ss_c = ss.s
+    cdef const qpms_trans_calculator *c = ss_c[0].c
+    cdef qpms_l_t sswf_lMax = 2*c[0].lMax+1
+    cdef qpms_y_t pnelem = qpms_lMax2nelem_sc(sswf_lMax)
+
+    cdef cdouble k = ssw.wavenumber
+
+    cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
+    cdef np.ndarray[complex,ndim=2] results = np.empty((evalpos_a.shape[0], pnelem), dtype=complex)
+    cdef cdouble *res
+    cdef cart3_t pos
+    cdef csph_t kdlj
+    cdef Py_ssize_t i, j
+    with nogil, wraparound(False), parallel():
+        res = <cdouble *> malloc(pnelem * sizeof(cdouble))
+        for i in prange(evalpos_a.shape[0]):
+            pos.x = evalpos_a[i,0]
+            pos.y = evalpos_a[i,1]
+            pos.z = evalpos_a[i,2]
+            kdlj = cart2csph(pos)
+            kdlj.r *= k
+            qpms_trans_calculator_test_sswf(c, res, kdlj, btyp_c)
+            for j in range(pnelem):
+                results[i,j] = res[j]
+        free(res)
+    return results.reshape(evalpos.shape[:-1] + (pnelem,))
+
 def gamma_inc(double a, cdouble x, int xbranch=0):
     cdef qpms_csf_result res
     complex_gamma_inc_e(a, x, xbranch, &res)
     return res.val
 
+def ipow_l_factor(n = None):
+    global ewald_factor_ipow_l
+    if n is not None:
+        ewald_factor_ipow_l = n
+    return ewald_factor_ipow_l
 
+def ipow_m_factor(n = None):
+    global ewald_factor_ipow_m
+    if n is not None:
+        ewald_factor_ipow_m = n
+    return ewald_factor_ipow_m

qpms/cytmatrices.pyx

@@ -286,6 +286,10 @@ cdef class TMatrixGenerator:
             self.holder = what
             self.g.function = qpms_tmatrix_generator_interpolator
             self.g.params = <void*>(<TMatrixInterpolator?>self.holder).rawpointer()
+        elif what == 0:
+            self.holder = 0
+            self.g.function = qpms_tmatrix_generator_zero
+            self.g.params = <void*> 0
         elif callable(what):
             warnings.warn("Custom python T-matrix generators are an experimental feature. Also expect it to be slow.")
             self.holder = what
@@ -403,3 +407,15 @@ cdef class TMatrixGenerator:
             EpsMuGenerator(outside), EpsMuGenerator(inside),
             ArcFunction(__ArcCylinder(r, h)), *args, **kwargs))
 
+    @staticmethod
+    def dummy():
+        """Returns a dummy T-matrix generator (returns a zero T-matrix).
+
+        Returns
+        -------
+        tmgen_dummy : TMatrixGenerator
+        """
+        return tmgen_dummy
+
+# pre-generate a dummy TMatrixGenerator (which returns zero T-matrix)
+tmgen_dummy = TMatrixGenerator(0)

qpms/ewald.c

@@ -7,11 +7,23 @@
 #include <complex.h>
 #include "tiny_inlines.h"
 #include "qpms_error.h"
+#include "lattices.h"
 #include <gsl/gsl_integration.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_sf_legendre.h>
 #include <gsl/gsl_sf_expint.h>
 
+/*
+ * Control additional phase factors
+ * that could possibly affect the resulting VSWF
+ * shape.
+ * Currently not used in ewald3_1_z_sigma_long()
+ */
+
+int ewald_factor_ipow_l = 0;
+int ewald_factor_ipow_m = 0;
+
+
 // parameters for the quadrature of integral in (4.6)
 #ifndef INTEGRATION_WORKSPACE_LIMIT
 #define INTEGRATION_WORKSPACE_LIMIT 30000
@@ -125,7 +137,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
 
   c->s1_constfacs_base = malloc(s1_constfacs_sz * sizeof(complex double));
   size_t s1_constfacs_sz_cumsum = 0;
-  for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
+  for (qpms_y_sc_t y = 0; y < c->nelem_sc; ++y) {
     qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
     if ((m + n) % 2 == 0) {
       c->s1_constfacs[y] = c->s1_constfacs_base + s1_constfacs_sz_cumsum;
@@ -161,7 +173,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
   c->S1_constfacs = malloc((1+c->nelem_sc) * sizeof(complex double *));
   //determine sizes
   size_t S1_constfacs_sz = 0;
-  for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
+  for (qpms_y_sc_t y = 0; y < c->nelem_sc; ++y) {
     qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
     const qpms_l_t L_M = n - abs(m);
     for(qpms_l_t j = 0; j <= L_M; ++j) { // outer sum
@@ -176,7 +188,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
 
   c->S1_constfacs_base = malloc(S1_constfacs_sz * sizeof(complex double));
   size_t S1_constfacs_sz_cumsum = 0; // second count
-  for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
+  for (qpms_y_sc_t y = 0; y < c->nelem_sc; ++y) {
     qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
     const complex double yfactor = -2 * ipow(n+1) * M_SQRTPI  
               * factorial((n-m)/2) * factorial((n+m)/2);
@@ -303,7 +315,7 @@ int ewald3_21_xy_sigma_long (
   const bool k_is_real = (cimag(k) == 0);
   assert((latdim & LAT_XYONLY) && (latdim & SPACE3D));
   assert((latdim & LAT1D) || (latdim & LAT2D));
-  const qpms_y_t nelem_sc = c->nelem_sc;
+  const qpms_y_sc_t nelem_sc = c->nelem_sc;
   assert(nelem_sc > 0);
   const qpms_l_t lMax = c->lMax;
   
@@ -435,7 +447,7 @@ int ewald3_21_xy_sigma_long (
       for(qpms_m_t m = -n; m <= n; ++m) {
         if((particle_shift.z == 0) && ((m+n) % 2 != 0)) // odd coefficients are zero.
           continue;
-        const qpms_y_t y = qpms_mn2y_sc(m, n);
+        const qpms_y_sc_t y = qpms_mn2y_sc(m, n);
         size_t constidx = 0; // constants offset
         const complex double e_imalpha_pq = cexp(I*m*arg_pq);
           complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
@@ -444,7 +456,7 @@ int ewald3_21_xy_sigma_long (
           assert((n-abs(m))/2 == c->s1_jMaxes[y]);
           for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
             complex double summand = rbeta_pq_div_k_pow[n-2*j]
-              * e_imalpha_pq  * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop
+              * e_imalpha_pq  * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow(m) // This line can actually go outside j-loop
               * gamma_pq_powm1[2*j]// * Gamma_pq[j] bellow (GGG) after error computation
               * c->s1_constfacs[y][j];
             if(err) {
@@ -455,6 +467,7 @@ int ewald3_21_xy_sigma_long (
             ckahanadd(&jsum, &jsum_c, summand);
           }
           jsum *= phasefac * factor1d; // PFC
+          jsum *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m);
           ckahanadd(target + y, target_c + y, jsum);
 #ifdef EWALD_AUTO_CUTOFF
           kahanadd(&lsum, &lsum_c, cabs(jsum));
@@ -483,6 +496,7 @@ int ewald3_21_xy_sigma_long (
             ckahanadd(&jsum, &jsum_c, jsummand);
           }
           jsum *= phasefac; // factor1d not here, off-axis sums not implemented/allowed.
+          jsum *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m);
           ckahanadd(target + y, target_c + y, jsum);
 #ifdef EWALD_AUTO_CUTOFF
           kahanadd(&lsum, &lsum_c, cabs(jsum));
@@ -496,7 +510,7 @@ int ewald3_21_xy_sigma_long (
 #ifdef EWALD_AUTO_CUTOFF
     ++li;
     double cursum_min = INFINITY;
-    for (qpms_y_t y = 0; y < nelem_sc; ++y) {
+    for (qpms_y_sc_t y = 0; y < nelem_sc; ++y) {
       const double a = cabs(target[y]);
       if (a) // Skip possible symmetry-induced zeros; TODO this might need to be done in completely different way, perhaps on per-y basis
         cursum_min = MIN(cursum_min, a);
@@ -516,10 +530,10 @@ ewald3_21_xy_sigma_long_end_point_loop:
  
   free(err_c);
   free(target_c);
-  for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
+  for(qpms_y_sc_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
     target[y] *= commonfac;
   if(err)
-    for(qpms_y_t y = 0; y < nelem_sc; ++y)
+    for(qpms_y_sc_t y = 0; y < nelem_sc; ++y)
       err[y] *= commonfac;
   return 0;
 }
@@ -549,7 +563,7 @@ int ewald3_1_z_sigma_long (
   assert(particle_shift.x == 0 && particle_shift.y == 0);
   const double beta_z = beta.z;
   const double particle_shift_z = particle_shift_z;
-  const qpms_y_t nelem_sc = c->nelem_sc;
+  const qpms_y_sc_t nelem_sc = c->nelem_sc;
   const qpms_l_t lMax = c->lMax;
   
   // Manual init of the ewald summation targets
@@ -600,7 +614,7 @@ int ewald3_1_z_sigma_long (
     // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array
     // and just fetched for each n
     for(qpms_l_t n = 0; n <= lMax; ++n) {
-        const qpms_y_t y = qpms_mn2y_sc(0, n);
+        const qpms_y_sc_t y = qpms_mn2y_sc(0, n);
         complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
         double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
         for(qpms_l_t j = 0; j <= n/2; ++j) {
@@ -634,10 +648,10 @@ int ewald3_1_z_sigma_long (
  
   free(err_c);
   free(target_c);
-  for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
+  for(qpms_y_sc_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
     target[y] *= commonfac;
   if(err)
-    for(qpms_y_t y = 0; y < nelem_sc; ++y)
+    for(qpms_y_sc_t y = 0; y < nelem_sc; ++y)
       err[y] *= commonfac;
   return 0;
 }
@@ -780,7 +794,7 @@ int ewald3_sigma_short(
 {
   const bool k_is_real = (cimag(k) == 0); // TODO check how the compiler optimises the loops
   const double kreal = creal(k);
-  const qpms_y_t nelem_sc = c->nelem_sc;
+  const qpms_y_sc_t nelem_sc = c->nelem_sc;
   const qpms_l_t lMax = c->lMax;
   gsl_integration_workspace *workspace = 
     gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT);
@@ -902,7 +916,8 @@ int ewald3_sigma_short(
         if(err)
           kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown
               * interr[n])); // TODO include also other errors
-        complex double thesummand = prefacn * R_pq_pown * leg * cintres[n] * e_beta_Rpq * e_imf * min1pow_m_neg(m);
+        complex double thesummand = prefacn * R_pq_pown * leg * cintres[n] * e_beta_Rpq * e_imf * min1pow(m);
+        thesummand *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m);
         ckahanadd(target + y, target_c + y, thesummand);
 #ifdef EWALD_AUTO_CUTOFF
         kahanadd(&lsum, &lsum_c, cabs(thesummand));
@@ -916,7 +931,7 @@ int ewald3_sigma_short(
 #ifdef EWALD_AUTO_CUTOFF
     ++li;
     double cursum_min = INFINITY;
-    for (qpms_y_t y = 0; y < nelem_sc; ++y) {
+    for (qpms_y_sc_t y = 0; y < nelem_sc; ++y) {
       const double a = cabs(target[y]);
       if (a) // Skip possible symmetry-induced zeros; TODO this might need to be done in completely different way, perhaps on per-y basis
         cursum_min = MIN(cursum_min, a);

qpms/ewald.h

@@ -26,22 +26,18 @@
 
 #ifndef EWALD_H
 #define EWALD_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #include <gsl/gsl_sf_result.h>
 #include <stdlib.h>
 #include <gsl/gsl_sf_legendre.h>
 #include <gsl/gsl_errno.h>
 #include <math.h> // for inlined lilgamma
-#include <complex.h>
 #include "qpms_types.h"
-#include "lattices.h"
-
-
-typedef enum {
-	QPMS_EWALD_LONG_RANGE = 1,
-	QPMS_EWALD_SHORT_RANGE = 2,
-	QPMS_EWALD_0TERM = 4,
-	QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
-} qpms_ewald_part;
+#include "lattices_types.h"
+#include <complex.h>
 
 
 /// Use this handler to ignore underflows of incomplete gamma.
@@ -55,11 +51,11 @@ gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;
  */
 typedef struct qpms_ewald3_constants_t { 
 	qpms_l_t lMax;
-	qpms_y_t nelem_sc;
+	qpms_y_sc_t nelem_sc;
 	/// The values of maximum \a j's in the long-range part summation, `[(l-|m|/2)]`.
 	qpms_l_t *s1_jMaxes;
 	/// The constant factors for the long range part of a 2D Ewald sum.
-	complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
+	_Complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
 	/* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
 	 * for m + n EVEN:
 	 *
@@ -73,7 +69,7 @@ typedef struct qpms_ewald3_constants_t {
 	 *
 	 * s1_constfacs[y(m,n)][j] = 0
 	 */
-	complex double *s1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
+	_Complex double *s1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
 	// similarly for the 1D z-axis aligned case; now the indices are [n][j] (as m == 0)
 	/// The constant factors for the long range part of a 1D Ewald sum along the \a z axis.
 	/** If the summation points lie along a different direction, use the formula for
@@ -86,7 +82,7 @@ typedef struct qpms_ewald3_constants_t {
 	// TODO indexing mechanisms
 	
 	/// The constant factors for the long range part of a 2D Ewald sum.
-	complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
+	_Complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
 	/* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
 	 * for m + n EVEN:
 	 *
@@ -100,7 +96,7 @@ typedef struct qpms_ewald3_constants_t {
 	 *
 	 * S1_constfacs[y(m,n)][j] = 0
 	 */
-	complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
+	_Complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
 	/// The constant factors for the long range part of a 1D Ewald sum along the \a z axis.
 	/** If the summation points lie along a different direction, use the formula for
 	 * 2D sum with additional factor of 
@@ -108,7 +104,7 @@ typedef struct qpms_ewald3_constants_t {
 	 */
 
 
-	complex double **s1_constfacs_1Dz; 
+	_Complex double **s1_constfacs_1Dz; 
 	/* These are the actual numbers now:
 	 * s1_constfacs_1Dz[n][j] =
 	 *   
@@ -116,7 +112,7 @@ typedef struct qpms_ewald3_constants_t {
 	 *   --------------------------
 	 *   j! * 2**(2*j) * (n - 2*j)!
 	 */
-	complex double *s1_constfacs_1Dz_base; ///<Internal pointer holding memory for the 1D Ewald sum constant factors.
+	_Complex double *s1_constfacs_1Dz_base; ///<Internal pointer holding memory for the 1D Ewald sum constant factors.
 
 	double *legendre0; /* now with GSL_SF_LEGENDRE_NONE normalisation, because this is what is
 			    * what the multipliers from translations.c count with.
@@ -139,13 +135,13 @@ void qpms_ewald3_constants_free(qpms_ewald3_constants_t *);
 /// Structure for holding complex-valued result of computation and an error estimate.
 /** Similar to gsl_sf_result, but with complex val. */
 typedef struct qpms_csf_result { 
-  complex double val; ///< Calculation result.
+  _Complex double val; ///< Calculation result.
   double err; ///< Error estimate.
 } qpms_csf_result;
 
 
 // [1, (A.9)]
-static inline complex double lilgamma(double t) {
+static inline _Complex double lilgamma(double t) {
   t = fabs(t);
   if (t >= 1) 
     return sqrt(t*t - 1);
@@ -154,8 +150,8 @@ static inline complex double lilgamma(double t) {
 }
 
 // [1, (A.8)], complex version of lilgamma()
-static inline complex double clilgamma(complex double z) {
-	complex double a1 = z - 1, a2 = z + 1;
+static inline _Complex double clilgamma(_Complex double z) {
+	_Complex double a1 = z - 1, a2 = z + 1;
 	// ensure  -pi/2 < arg(z + 1) < 3*pi/2
 	if (creal(a2) < 0 && cimag(a2) <= 0) 
 		a2 = -csqrt(a2);
@@ -180,7 +176,7 @@ static inline complex double clilgamma(complex double z) {
  * even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
  * The side of the branch cut can be determined using `signbit(creal(z))`.
  */
-int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result);
+int cx_gamma_inc_series_e(double a, _Complex double z, qpms_csf_result * result);
 
 /// Incomplete Gamma function as continued fractions.
 /** 
@@ -192,7 +188,7 @@ int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result);
  * even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
  * The side of the branch cut can be determined using `signbit(creal(z))`.
  */
-int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result);
+int cx_gamma_inc_CF_e(double a, _Complex double z, qpms_csf_result * result);
 
 /// Incomplete gamma for complex second argument.
 /** 
@@ -209,7 +205,7 @@ int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result);
  * Another than principal branch can be selected using non-zero \a m 
  * argument.
  */
-int complex_gamma_inc_e(double a, complex double x,
+int complex_gamma_inc_e(double a, _Complex double x,
 	/// Branch index.
 	/** If zero, the principal value is calculated.	 
 	 * Other branches might be chosen using non-zero \a m.
@@ -226,7 +222,7 @@ int complex_gamma_inc_e(double a, complex double x,
 
 /// Exponential integral for complex second argument.
 /** If x is (almost) positive real, it just uses gsl_sf_expint_En_e(). */
-int complex_expint_n_e(int n, complex double x, qpms_csf_result *result);
+int complex_expint_n_e(int n, _Complex double x, qpms_csf_result *result);
 
 
 /// Hypergeometric 2F2, used to calculate some errors.
@@ -245,15 +241,15 @@ int ewald32_sr_integral(double r, double k, double n, double eta, double *result
  * unsuitable especially for big values of \a maxn.
  * 
  */
-void ewald3_2_sigma_long_Delta(complex double *target, double *target_err, int maxn, complex double x,
-	       int xbranch, complex double z);
+void ewald3_2_sigma_long_Delta(_Complex double *target, double *target_err, int maxn, _Complex double x,
+	       int xbranch, _Complex double z);
 
 /// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum.
 /** This function always uses Kambe's (corrected) recurrent formula.
  * For production, use ewald3_2_sigma_long_Delta() instead.
  */
-void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_err, int maxn, complex double x,
-	       int xbranch, complex double z, _Bool bigimz);
+void ewald3_2_sigma_long_Delta_recurrent(_Complex double *target, double *target_err, int maxn, _Complex double x,
+	       int xbranch, _Complex double z, _Bool bigimz);
 
 /// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum.
 /** This function always uses Taylor expansion in \a z.
@@ -263,26 +259,26 @@ void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_
  * parameters maxn = 40, z = 0.5, x = -3. This might be related to the exponential growth
  * of the error.
  */
-void ewald3_2_sigma_long_Delta_series(complex double *target, double *target_err, int maxn, complex double x,
-	       int xbranch, complex double z);
+void ewald3_2_sigma_long_Delta_series(_Complex double *target, double *target_err, int maxn, _Complex double x,
+	       int xbranch, _Complex double z);
 
 // General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space
 
 /// The Ewald sum "self-interaction" term that appears in the lattice sums with zero (direct-space) Bravais lattice displacement.
-int ewald3_sigma0(complex double *result, ///< Pointer to save the result (single complex double).
+int ewald3_sigma0(_Complex double *result, ///< Pointer to save the result (single _Complex double).
 	       	double *err, ///< Pointer to save the result error estimate (single double).
 		const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
 		double eta, ///< Ewald parameter.
-		complex double wavenumber ///< Wavenumber of the background medium.
+		_Complex double wavenumber ///< Wavenumber of the background medium.
 );
 
 /// Short-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{S}(\vect k,\vect s)\f$.
 int ewald3_sigma_short(
-		complex double *target_sigmasr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{S} \f$, must be `c->nelem_sc` long.
+		_Complex double *target_sigmasr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{S} \f$, must be `c->nelem_sc` long.
 		double *target_sigmasr_y_err, ///< Target array for error estimates, must be `c->nelem_sc` long or `NULL`.
 		const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
 		double eta, ///< Ewald parameter.
-	       	complex double wavenumber, ///< Wavenumber of the background medium. 
+	       	_Complex double wavenumber, ///< Wavenumber of the background medium. 
 		/// Lattice dimensionality. 
 		/** Ignored apart from asserts and possible optimisations, as the SR formula stays the same. */
 		LatticeDimensionality latdim,
@@ -293,7 +289,7 @@ int ewald3_sigma_short(
 		 *  In such case, it is the responsibility of the caller to deallocate
 		 *  the generator.
 		 */
-		PGen *pgen_R,
+		struct PGen *pgen_R,
 		/// Indicates whether pgen_R already generates shifted points.
 		/** If false, the behaviour corresponds to the old ewald32_sigma_short_points_and_shift(),
 		 * so the function assumes that the generated points correspond to the unshifted Bravais lattice,
@@ -310,11 +306,11 @@ int ewald3_sigma_short(
 
 /// Long-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{L}(\vect k,\vect s)\f$.
 int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, depending on latdim
-		complex double *target_sigmalr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{L} \f$, must be `c->nelem_sc` long.
+		_Complex double *target_sigmalr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{L} \f$, must be `c->nelem_sc` long.
 		double *target_sigmalr_y_err, ///< Target array for error estimates, must be `c->nelem_sc` long or `NULL`.
 		const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
 		double eta, ///< Ewald parameter.
-	       	complex double wavenumber, ///< Wavenumber of the background medium. 
+	       	_Complex double wavenumber, ///< Wavenumber of the background medium. 
 	        double unitcell_volume, ///< Volume of the (direct lattice) unit cell (with dimension corresponding to the lattice dimensionality).
 		/// Lattice dimensionality. 
 		LatticeDimensionality latdim,
@@ -325,7 +321,7 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep
 		 *  In such case, it is the responsibility of the caller to deallocate
 		 *  the generator.
 		 */
-		PGen *pgen_K, 
+		struct PGen *pgen_K, 
 		/// Indicates whether pgen_K already generates shifted points.
 		/** If false, the behaviour corresponds to the old ewald32_sigma_long_points_and_shift(),
 		 * so the function assumes that the generated points correspond to the unshifted reciprocal Bravais lattice,
@@ -339,4 +335,13 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep
 		cart3_t particle_shift
 		);
 		
+
+// If nonzero, adds an additional factor \f$ i^{nl} \f$ to the Ewald sum result (for debugging).
+extern int ewald_factor_ipow_l;
+// If nonzero, adds an additional factor \f$ i^{nm} \f$ to the Ewald sum result (for debubbing).
+extern int ewald_factor_ipow_m;
+
+#ifdef __cplusplus
+}
+#endif
 #endif //EWALD_H

qpms/ewaldsf.c

@@ -14,6 +14,7 @@
 #include <stdbool.h>
 #include <Faddeeva.h>
 #include "tiny_inlines.h"
+#include "qpms_error.h"
 
 // Some magic constants
 

qpms/gaunt.h

@@ -3,6 +3,10 @@
  */
 #ifndef GAUNT_H
 #define GAUNT_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #include <stdlib.h>
 
 #define _GAUNT_H_MIN(x,y) (((x) > (y)) ? (y) : (x))
@@ -30,4 +34,7 @@ double const * gaunt_table_retrieve_allq(int m, int n, int mu, int nu);
 int gaunt_table_or_xu_fill(double *target, int m, int n, int mu, int nu);
 #endif //GAUNT_PRECOMPILED
 
+#ifdef __cplusplus
+}
+#endif
 #endif //GAUNT_H

qpms/groups.h

@@ -23,6 +23,9 @@
  */
 #ifndef QPMS_GROUPS_H
 #define QPMS_GROUPS_H
+#ifdef __cplusplus
+extern "C" {
+#endif
 
 #include "qpms_types.h"
 #include <assert.h>
@@ -35,7 +38,7 @@ struct qpms_finite_group_irrep_t {
 	/** The r-th row, c-th column of the representation of the i'th element is retrieved as
 	 *  m[i * dim * dim + r * dim + c]
 	 */
-	complex double *m; 
+	_Complex double *m; 
 };
 
 /// A point group with its irreducible representations and some metadata.
@@ -93,4 +96,7 @@ qpms_iri_t qpms_finite_group_find_irrep_by_name(qpms_finite_group_t *G, char *na
 extern const qpms_finite_group_t QPMS_FINITE_GROUP_TRIVIAL;
 extern const qpms_finite_group_t QPMS_FINITE_GROUP_TRIVIAL_G;
 
+#ifdef __cplusplus
+}
+#endif
 #endif // QPMS_GROUPS_H

qpms/indexing.h

@@ -1,57 +1,106 @@
 /*! \file indexing.h
  * \brief Various index conversion functions.
+ *
+ * Index conventions in QPMS
+ * -------------------------
+ *
+ * Depending on circumstances (e.g. whether scalar or vector spherical wavefunctions
+ * are being used), QPMS uses two mappings between the multipole degree/order index 
+ * pair (\a l, \a m), due to historical reasons ocassionally also noted as (\a n, \a m),
+ * and a flat non-negative integer index.
+ * 
+ * Both mappings map the (\a l, \a m) pair in lexicographic order (with 
+ * the basic constraint \f$ |m| \le l \f$), the difference is that the "scalar"
+ * mapping starts from \a l = 0, whereas the "vector" mapping starts from \a l = 1, 
+ * so that for the zeroth elements we have
+ *  * (0, 0) ↔ 0, (1, -1) ↔ 1 in the "scalar" mapping,
+ *  * (1, -1) ↔ 0 in the "vector" mapping.
+ *
+ * The vector SWF flat index is typically denoted by `y` and its type is \ref qpms_y_t, 
+ * whereas  * the scalar SWF flat index is typically denoted by `y_sc` and its 
+ * type is \ref qpms_y_sc_t.
+ *
+ * Moreover, there is another mapping that encodes also the electric/magnetic/scalar
+ * WSWF type index, see \ref qpms_uvswfi_t.
+ *
  */
 #ifndef QPMS_INDEXING_H
 #define QPMS_INDEXING_H
+#ifdef __cplusplus
+extern "C" {
+#endif
 
 #include "qpms_types.h"
 #include <math.h>
 #include "optim.h"
 
-
-
-static inline qpms_y_t qpms_mn2y(qpms_m_t m, qpms_l_t n) {
+/// Mapping from multipole degree, order to flat index (degree ≥ 1)
+static inline qpms_y_t qpms_mn2y(
+    qpms_m_t m, ///< multipole order
+    qpms_l_t n ///< multipole degree (a.k.a. \a l); must be greater than 0
+    ) {
         return n * (n + 1) + m - 1;
 }
 
+/// Mapping from flat index to multipole degree (degree ≥ 1)
 static inline qpms_lm_t qpms_y2n(qpms_y_t y) {
         //return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
         return sqrt(y+1);
 }
 
+/// Mapping from flat index, multipole degree (≥ 1) to multipole order
 static inline qpms_m_t qpms_yn2m(qpms_y_t y, qpms_l_t n) {
         return y-qpms_mn2y(0,n);
 }
 
+/// Mapping from flat index to multipole degree, order (degree ≥ 1)
 static inline void qpms_y2mn_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n){
         *m=qpms_yn2m(y,*n=qpms_y2n(y));
 }
 
+/// Number of spherical multipoles with `lmax` ≥ degree ≥ 1
 static inline qpms_y_t qpms_lMax2nelem(qpms_l_t lmax){
 	return lmax * ((qpms_y_t)lmax + 2);
 }
 
+/// Inverse of qpms_lMax2nelem(); returns -1 if there is no exact match.
+static inline qpms_l_t qpms_nelem2lMax(qpms_y_t nelem) {
+  qpms_l_t l = qpms_y2n(nelem);
+  if (l == -qpms_yn2m(nelem, l))
+    return l - 1;
+  else
+    return -1;
+}
+
 // Scalar versions: they have a place for the 0, 0 term in the beginning
 
-static inline qpms_y_t qpms_mn2y_sc(qpms_m_t m, qpms_l_t n) {
+/// Mapping from multipole degree, order to flat index (degree ≥ 0)
+static inline qpms_y_sc_t qpms_mn2y_sc(
+    qpms_m_t m, ///< multipole order
+    qpms_l_t n ///< multipole degree (a.k.a. \a l); muste be greater or equal than 0
+    ) {
         return n * (n + 1) + m;
 }
 
-static inline qpms_lm_t qpms_y2n_sc(qpms_y_t y) {
+/// Mapping from flat index to multipole degree (degree ≥ 0)
+static inline qpms_lm_t qpms_y2n_sc(qpms_y_sc_t y) {
         //return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
         return sqrt(y);
 }
 
-static inline qpms_m_t qpms_yn2m_sc(qpms_y_t y, qpms_l_t n) {
+/// Mapping from flat index, multipole degree (≥ 0) to multipole order
+static inline qpms_m_t qpms_yn2m_sc(qpms_y_sc_t y, qpms_l_t n) {
         return y-qpms_mn2y_sc(0,n);
 }
 
-static inline void qpms_y2mn_sc_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n){
+/// Mapping from flat index to multipole degree, order (degree ≥ 0)
+static inline void qpms_y2mn_sc_p(qpms_y_sc_t y, qpms_m_t *m, qpms_l_t *n){
         *m=qpms_yn2m_sc(y,*n=qpms_y2n_sc(y));
 }
 
-static inline qpms_y_t qpms_lMax2nelem_sc(qpms_l_t lmax){
-	return lmax * ((qpms_y_t)lmax + 2) + 1;
+/// Number of spherical multipoles with `lmax` ≥ degree ≥ 0
+static inline qpms_y_sc_t qpms_lMax2nelem_sc(qpms_l_t lmax){
+	return lmax * ((qpms_y_sc_t)lmax + 2) + 1;
 }
 
 // TODO maybe enable crashing / validity control by macro definitions...
@@ -69,7 +118,7 @@ static const qpms_uvswfi_t QPMS_UI_L00 = 0;
 /** Returns a non-zero value if the u value is invalid. */
 static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
 		qpms_vswf_type_t *t, qpms_m_t *m, qpms_l_t *n) {
-	*t = u & 3;
+	*t = (qpms_vswf_type_t)(u & 3);
 	qpms_y_sc_t y_sc = u / 4;
 	qpms_y2mn_sc_p(y_sc, m, n);
 	// Test validity
@@ -82,7 +131,7 @@ static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
 /** Does *not* allow for longitudinal waves. */
 static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u,
 		qpms_vswf_type_t *t, qpms_y_t *y) {
-	*t = u & 3;
+	*t = (qpms_vswf_type_t)(u & 3);
 	*y = u / 4 - 1;
 	if (QPMS_UNLIKELY(*t == 0 || *t == 3)) return QPMS_ERROR;
 	if (QPMS_UNLIKELY(*y < 0)) return QPMS_ERROR;
@@ -95,7 +144,7 @@ static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u,
  */
 static inline qpms_errno_t qpms_uvswfi2ty_l(qpms_uvswfi_t u,
 		qpms_vswf_type_t *t, qpms_y_t *y) {
-	*t = u & 3;
+	*t = (qpms_vswf_type_t)(u & 3);
 	*y = u / 4 - 1;
 	if (QPMS_UNLIKELY(*t == 3)) return QPMS_ERROR;
 	if (QPMS_UNLIKELY(*y < 0)) return QPMS_ERROR;
@@ -110,4 +159,7 @@ static inline qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) {
 }
 
 
+#ifdef __cplusplus
+}
+#endif
 #endif //QPMS_INDEXING_H

qpms/kahansum.h

@@ -3,8 +3,9 @@
  */
 #ifndef KAHANSUM_H
 #define KAHANSUM_H
-
-#include <complex.h>
+#ifdef __cplusplus
+extern "C" {
+#endif
 
 static inline void kahaninit(double * const sum, double * const compensation) {
 	*sum = 0;
@@ -19,16 +20,19 @@ static inline void kahanadd(double *sum, double *compensation, double input) {
 }
 
 
-static inline void ckahaninit(complex double * const sum, complex double * const compensation) {
+static inline void ckahaninit(_Complex double * const sum, _Complex double * const compensation) {
 	*sum = 0;
 	*compensation = 0;
 }
 
-static inline void ckahanadd(complex double *sum, complex double *compensation, complex double input) {
-	complex double compensated_input = input - *compensation;
-	complex double nsum = *sum + compensated_input;
+static inline void ckahanadd(_Complex double *sum, _Complex double *compensation, _Complex double input) {
+	_Complex double compensated_input = input - *compensation;
+	_Complex double nsum = *sum + compensated_input;
 	*compensation = (nsum - *sum) - compensated_input;
 	*sum = nsum;
 }
 
+#ifdef __cplusplus
+}
+#endif
 #endif //KAHANSUM_H

qpms/lattices.h

@@ -1,9 +1,15 @@
 /*! \file lattices.h
  * \brief Lattice point generators and lattice vector analysis / transformation.
  *
+ * \bug Header file not C++ compatible.
  */
 #ifndef LATTICES_H
 #define LATTICES_H
+#ifdef __cplusplus // FIXME Not C++ compatible. Include "lattices_types.h" for minimal necessary enum decls.
+extern "C" {
+#endif
+
+#include "lattices_types.h"
 
 #include <math.h>
 #include <stdbool.h>
@@ -30,22 +36,6 @@
  *
  */
 
-typedef enum LatticeDimensionality {
-	LAT1D = 1,
-	LAT2D = 2,
-	LAT3D = 4,
-	SPACE1D = 8,
-	SPACE2D = 16,
-	SPACE3D = 32,
-	LAT_1D_IN_3D = 33,
-	LAT_2D_IN_3D = 34,
-	LAT_3D_IN_3D = 40,
-	// special coordinate arrangements (indicating possible optimisations)
-	LAT_ZONLY = 64,
-	LAT_XYONLY = 128,
-	LAT_1D_IN_3D_ZONLY = 97, // LAT1D | SPACE3D | 64
-	LAT_2D_IN_3D_XYONLY = 162 // LAT2D | SPACE3D | 128
-} LatticeDimensionality;
 
 inline static bool LatticeDimensionality_checkflags(
 		LatticeDimensionality a, LatticeDimensionality flags_a_has_to_contain) {
@@ -945,4 +935,7 @@ int honeycomb_lattice_gen_extend_to_steps(honeycomb_lattice_gen_t *g, int maxste
 int honeycomb_lattice_gen_extend_to_r(honeycomb_lattice_gen_t *g, double r);
 void honeycomb_lattice_gen_free(honeycomb_lattice_gen_t *g);
 
+#ifdef __cplusplus
+}
+#endif
 #endif // LATTICES_H

qpms/materials.c

@@ -5,6 +5,7 @@
 #include <string.h>
 #include "materials.h"
 #include "qpms_error.h"
+#include <stdbool.h>
 
 #define SQ(x) ((x)*(x))
 

qpms/materials.h

@@ -3,7 +3,12 @@
  */
 #ifndef QPMS_MATERIALS_H
 #define QPMS_MATERIALS_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #include "qpms_types.h"
+#include <complex.h>
 #include <gsl/gsl_spline.h>
 
 #ifndef SPEED_OF_LIGHT
@@ -20,36 +25,36 @@ typedef struct qpms_epsmu_generator_t {
 	 *  qpms_permittivity_interpolator_epsmu_g(),
 	 *  qpms_lorentzdrude_epsmu_g().
 	 */
-	qpms_epsmu_t (*function) (complex double omega, const void *params);
+	qpms_epsmu_t (*function) (_Complex double omega, const void *params);
 	const void *params;
 } qpms_epsmu_generator_t;
 
 /// Convenience function for generating material properties at given frequency.
 static inline qpms_epsmu_t qpms_epsmu_generator_eval(
-		qpms_epsmu_generator_t gen, complex double omega) {
+		qpms_epsmu_generator_t gen, _Complex double omega) {
 	return gen.function(omega, gen.params);
 }
 
 /// Constant optical property "generator" for qpms_epsmu_generator_t.
-qpms_epsmu_t qpms_epsmu_const_g(complex double omega, ///< Frequency ignored.
+qpms_epsmu_t qpms_epsmu_const_g(_Complex double omega, ///< Frequency ignored.
 	const void *epsmu ///< Points to the qpms_epsmu_t to be returned.
 	);
 
 /// Gets refractive index of a material from its permeability and permittivity.
 /** \f[ n = \sqrt{\mu_r \varepsilon_r} \f] */
-static inline complex double qpms_refindex(qpms_epsmu_t em) {
+static inline _Complex double qpms_refindex(qpms_epsmu_t em) {
 	return csqrt(em.eps * em.mu);
 }
 
 /// Gets wave number \a k from angular frequency and material permeability and permittivity.
 /** \f[ k = \frac{n\omega}{c_0}  = \frac{\omega\sqrt{\mu_r \varepsilon_r}}{c_0} \f] */
-static inline complex double qpms_wavenumber(complex double omega, qpms_epsmu_t em) {
+static inline _Complex double qpms_wavenumber(_Complex double omega, qpms_epsmu_t em) {
 	return qpms_refindex(em)*omega/SPEED_OF_LIGHT;
 }
 
 /// Gets (relative) wave impedance \f$ \eta_r \f$ from material permeability and permittivity.
 /** \eta_r = \sqrt{\mu_r  / \varepsilon_r} \f] */
-static inline complex double qpms_waveimpedance(qpms_epsmu_t em) {
+static inline _Complex double qpms_waveimpedance(qpms_epsmu_t em) {
 	return csqrt(em.mu / em.eps);
 }
 
@@ -67,7 +72,7 @@ typedef struct qpms_ldparams_triple_t {
  *  \f]
  */
 typedef struct qpms_ldparams_t {
-	complex double eps_inf; ///< Permittivity at infinity.
+	_Complex double eps_inf; ///< Permittivity at infinity.
 	double omega_p; ///< Plasma frequency.
 	size_t n; ///< Number of "oscillators".
 	qpms_ldparams_triple_t data[]; ///< "Oscillator" parameters.
@@ -86,14 +91,14 @@ extern const qpms_ldparams_t *const QPMS_LDPARAMS_PT; ///< Lorentz-Drude paramet
 extern const qpms_ldparams_t *const QPMS_LDPARAMS_W ; ///< Lorentz-Drude parameters from \cite rakic_optical_1998 for tungsten.
 
 /// Lorentz-Drude permittivity.
-complex double qpms_lorentzdrude_eps(complex double omega, const qpms_ldparams_t *);
+_Complex double qpms_lorentzdrude_eps(_Complex double omega, const qpms_ldparams_t *);
 
 /// Lorentz-Drude optical properties, with relative permeability set always to one.
-qpms_epsmu_t qpms_lorentzdrude_epsmu(complex double omega, const qpms_ldparams_t *);
+qpms_epsmu_t qpms_lorentzdrude_epsmu(_Complex double omega, const qpms_ldparams_t *);
 
 /// Lorentz-Drude optical properties, with relative permeability set always to one, compatible with qpms_epsmu_generator_t.
 qpms_epsmu_t qpms_lorentzdrude_epsmu_g(
-		complex double omega, 
+		_Complex double omega, 
 		const void *ldparams ///< Lorentz-Drude parameters, in reality const qpms_ldparams_t *.
 		);
 
@@ -125,14 +130,14 @@ qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(
                 );
 
 /// Evaluates interpolated material permittivity at a given angular frequency.
-complex double qpms_permittivity_interpolator_eps_at_omega(
+_Complex double qpms_permittivity_interpolator_eps_at_omega(
                 const qpms_permittivity_interpolator_t *interp, double omega_SI);
 
 /// Evaluates interpolated material permittivity at a given angular frequency, qpms_epsmu_generator_t compatible version.
 /** Permeability is always set to one. Imaginary part of omega is discarded.
  */
 qpms_epsmu_t qpms_permittivity_interpolator_epsmu_g(
-		complex double omega, ///< Angular frequency. The imaginary part is ignored!
+		_Complex double omega, ///< Angular frequency. The imaginary part is ignored!
                 const void * interpolator ///< Interpolator of type qpms_permittivity_interpolator_t
 		);
 
@@ -148,14 +153,17 @@ double qpms_permittivity_interpolator_omega_max(
 void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp);
 
 /// Relative permittivity from the Drude model.
-static inline complex double qpms_drude_epsilon(
-                complex double eps_inf, ///< Relative permittivity "at infinity".
-                complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material.
-                complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material.
-                complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated.
+static inline _Complex double qpms_drude_epsilon(
+                _Complex double eps_inf, ///< Relative permittivity "at infinity".
+                _Complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material.
+                _Complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material.
+                _Complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated.
                 ) {
         return eps_inf - omega_p*omega_p/(omega*(omega+I*gamma_p));
 }
 
 
+#ifdef __cplusplus
+}
+#endif
 #endif //QPMS_MATERIALS_H

qpms/normalisation.h

@@ -2,14 +2,18 @@
  * \brief Convention-dependent coefficients for VSWFs.
  *
  * See also @ref qpms_normalisation_t and @ref vswf_conventions.
+ *
+ * \bug Header file not C++ compatible.
  */
 #ifndef NORMALISATION_H
 #define NORMALISATION_H
+#ifdef __cplusplus //FIXME not C++ compatible yet (enum bit operations)
+extern "C" {
+#endif
 
 #include "qpms_types.h"
 #include "qpms_error.h"
 #include <math.h>
-#include <complex.h>
 #include "indexing.h"
 #include "optim.h"
 
@@ -36,8 +40,8 @@ static inline double qpms_normalisation_normfactor(qpms_normalisation_t norm, qp
  * 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);
+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 (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_MINUS)) fac *= -1;
 	if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_I)) fac *= I;
 	if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@@ -51,8 +55,8 @@ static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation
  * 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);
+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;
 }
 
@@ -62,8 +66,8 @@ static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t no
  * 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);
+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 (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_MINUS)) fac *= -1;
 	if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_I)) fac *= I;
 	if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@@ -77,14 +81,14 @@ static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation
  * 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);
+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) {
+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);
 }
@@ -95,8 +99,8 @@ static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t
  * 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);
+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 (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_MINUS)) fac *= -1;
 	if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_I)) fac *= I;
 	if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@@ -109,13 +113,13 @@ static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation
  * 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);
+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) {
+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) {
@@ -156,7 +160,7 @@ static inline qpms_normalisation_t qpms_normalisation_dual(qpms_normalisation_t
  * 	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) {
+static inline _Complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) {
 	switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT) 
 			& (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
 		case 0:
@@ -194,7 +198,7 @@ static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t nor
  *
  *
  */
-static inline complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) {
+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(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT) 
 			& (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
@@ -283,4 +287,7 @@ static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t
 }
 #endif
 
+#ifdef __cplusplus
+}
+#endif
 #endif //NORMALISATION_H

qpms/parsing.h

@@ -3,6 +3,9 @@
  */
 #ifndef QPMS_PARSING_H
 #define QPMS_PARSING_H
+#ifdef __cplusplus
+extern "C" {
+#endif
 
 #include <stddef.h>
 
@@ -57,4 +60,7 @@ size_t qpms_parse_doubles_fromfile(
   const char *filepath //< File to read from, or NULL, "", "-" to read from stdin.
 );
 
+#ifdef __cplusplus
+}
+#endif
 #endif // QPMS_PARSING_H

qpms/pointgroups.c

@@ -1,6 +1,7 @@
 #include "pointgroups.h"
 #include <search.h>
 #include <stdlib.h>
+#include <stdbool.h>
 #include <assert.h>
 
 double qpms_pg_quat_cmp_atol = QPMS_QUAT_ATOL;

qpms/pointgroups.h

@@ -3,6 +3,9 @@
  */
 #ifndef POINTGROUPS_H
 #define POINTGROUPS_H
+#ifdef __cplusplus
+extern "C" {
+#endif
 
 #include "qpms_error.h"
 #include "quaternions.h"
@@ -18,9 +21,9 @@ static inline _Bool qpms_pg_is_finite_axial(qpms_pointgroup_class cls) {
         case QPMS_PGS_DN: 
         case QPMS_PGS_DND: 
         case QPMS_PGS_DNH: 
-		return true;
+		return 1;
 	default:
-		return false;
+		return 0;
 	}
 }
 
@@ -81,4 +84,7 @@ _Bool qpms_pg_is_subgroup(qpms_pointgroup_t a, qpms_pointgroup_t b);
 
 
 
+#ifdef __cplusplus
+}
+#endif
 #endif //POINTGROUPS_H

qpms/qpms_c.pxd

@@ -1,4 +1,6 @@
-from .qpms_cdefs cimport qpms_finite_group_t
+from .qpms_cdefs cimport qpms_finite_group_t, qpms_scatsys_t, qpms_iri_t,\
+     qpms_scatsys_at_omega_t, qpms_scatsys_at_omega_k_t
+from .cymaterials cimport EpsMuGenerator
 
 cdef class FinitePointGroup:
     cdef readonly bint owns_data
@@ -6,3 +8,30 @@ cdef class FinitePointGroup:
 
     cdef inline qpms_finite_group_t *rawpointer(self):
         return self.G
+
+
+cdef class ScatteringSystem:
+    cdef list tmgobjs # here we keep the references to occuring TMatrixFunctions (and hence BaseSpecs and TMatrixGenerators)
+    #cdef list Tmatrices # Here we keep the references to occuring T-matrices
+    cdef EpsMuGenerator medium_holder # Here we keep the reference to medium generator
+    cdef qpms_scatsys_t *s
+    cdef FinitePointGroup sym
+
+    cdef qpms_iri_t iri_py2c(self, iri, allow_None=?)
+
+
+cdef class _ScatteringSystemAtOmegaK:
+    cdef qpms_scatsys_at_omega_k_t sswk
+    cdef _ScatteringSystemAtOmega ssw_pyref
+
+    cdef inline qpms_scatsys_at_omega_k_t *rawpointer(self):
+        return &self.sswk
+
+
+cdef class _ScatteringSystemAtOmega:
+    cdef qpms_scatsys_at_omega_t *ssw
+    cdef ScatteringSystem ss_pyref
+    
+    cdef inline qpms_scatsys_at_omega_t *rawpointer(self):
+        return self.ssw
+ 

qpms/qpms_c.pyx

@@ -11,14 +11,14 @@ from .qpms_cdefs cimport *
 from .cyquaternions cimport IRot3, CQuat
 from .cybspec cimport BaseSpec
 from .cycommon cimport make_c_string
-from .cycommon import string_c2py, PointGroupClass, BesselType
+from .cycommon import string_c2py, PointGroupClass, BesselType, EwaldPart
 from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator, TMatrixInterpolator
 from .cymaterials cimport EpsMuGenerator, EpsMu
 from libc.stdlib cimport malloc, free, calloc
 import warnings
 
-from cython.parallel cimport prange, parallel
-from cython cimport boundscheck, wraparound
+from cython.parallel import prange, parallel
+from cython import boundscheck, wraparound
 
 # Set custom GSL error handler. N.B. this is obviously not thread-safe.
 cdef char *pgsl_err_reason
@@ -369,12 +369,6 @@ cdef class ScatteringSystem:
     Currently, it does not have a standard constructor. Use the
     ScatteringSystem.create() method instead.
     '''
-    cdef list tmgobjs # here we keep the references to occuring TMatrixFunctions (and hence BaseSpecs and TMatrixGenerators)
-    #cdef list Tmatrices # Here we keep the references to occuring T-matrices
-    cdef EpsMuGenerator medium_holder # Here we keep the reference to medium generator
-    cdef qpms_scatsys_t *s
-    cdef FinitePointGroup sym
-
     cdef qpms_iri_t iri_py2c(self, iri, allow_None = True):
         if iri is None and allow_None:
             return QPMS_NO_IRREP
@@ -699,7 +693,7 @@ cdef class ScatteringSystem:
                 self.s, iri, 0)
         return target_np
 
-    def pack_matrix(self, fullmatrix, iri):
+    def pack_matrix(self, fullmatrix, iri, version='normal'):
         """Converts (projects) a matrix into an irrep subspace.
 
         Parameters
@@ -731,10 +725,14 @@ cdef class ScatteringSystem:
         cdef np.ndarray[np.complex_t, ndim=2] target_np = np.empty(
                 (rlen, rlen), dtype=complex, order='C')
         cdef cdouble[:,::1] target_view = target_np
-        qpms_scatsys_irrep_pack_matrix(&target_view[0][0], &fullmatrix_view[0][0],
+        if version == 'stupid':
+            qpms_scatsys_irrep_pack_matrix_stupid(&target_view[0][0], &fullmatrix_view[0][0],
+                self.s, iri)
+        else:
+            qpms_scatsys_irrep_pack_matrix(&target_view[0][0], &fullmatrix_view[0][0],
                 self.s, iri)
         return target_np
-    def unpack_matrix(self, packedmatrix, iri):
+    def unpack_matrix(self, packedmatrix, iri, version='normal'):
         """Unpacks an "irrep-packed" excitation coefficient vector to full coordinates.
 
         Parameters
@@ -767,7 +765,11 @@ cdef class ScatteringSystem:
         cdef np.ndarray[np.complex_t, ndim=2] target_np = np.empty(
                 (flen, flen), dtype=complex, order='C')
         cdef cdouble[:,::1] target_view = target_np
-        qpms_scatsys_irrep_unpack_matrix(&target_view[0][0], &packedmatrix_view[0][0],
+        if version == 'stupid':
+            qpms_scatsys_irrep_unpack_matrix_stupid(&target_view[0][0], &packedmatrix_view[0][0],
+                self.s, iri, 0)
+        else:
+            qpms_scatsys_irrep_unpack_matrix(&target_view[0][0], &packedmatrix_view[0][0],
                 self.s, iri, 0)
         return target_np
 
@@ -819,6 +821,16 @@ cdef class ScatteringSystem:
                     qpms_scatsys_periodic_build_translation_matrix_full(&target_view[0][0], self.s, wavenumber, &blochvector_c, eta)
         return target
 
+    def irrep_transform_matrix(self, qpms_iri_t iri):
+        self.check_s()
+        cdef size_t rlen = self.saecv_sizes[iri]
+        cdef size_t fullen = self.fecv_size
+        cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
+                (rlen, fullen), dtype=complex, order='C')
+        cdef cdouble[:,::1] target_view = target
+        qpms_scatsys_irrep_transform_matrix(&target_view[0][0], self.s, iri)
+        return target
+
     def translation_matrix_packed(self, cdouble wavenumber, qpms_iri_t iri, J = QPMS_HANKEL_PLUS):
         self.check_s()
         cdef size_t rlen = self.saecv_sizes[iri]
@@ -962,8 +974,8 @@ cdef class ScatteringSystem:
         cdef ccart3_t res
         cdef cart3_t pos
         cdef Py_ssize_t i
-        with nogil, parallel(), boundscheck(False), wraparound(False):
-            for i in prange(evalpos_a.shape[0]):
+        with nogil:
+            for i in range(evalpos_a.shape[0]):
                 pos.x = evalpos_a[i,0]
                 pos.y = evalpos_a[i,1]
                 pos.z = evalpos_a[i,2]
@@ -1003,18 +1015,13 @@ def empty_lattice_modes_xy(EpsMu epsmu, reciprocal_basis, wavevector, double max
     free(omegas_c)
     return omegas
 
+
 cdef class _ScatteringSystemAtOmegaK:
     '''
     Wrapper over the C qpms_scatsys_at_omega_k_t structure
 
     This represents an infinite periodic system with a given (fixed) frequency and Bloch vector.
     '''
-    cdef qpms_scatsys_at_omega_k_t sswk
-    cdef _ScatteringSystemAtOmega ssw_pyref
-
-    cdef qpms_scatsys_at_omega_k_t *rawpointer(self):
-        return &self.sswk
-    
     property eta:
         """Ewald parameter η"""
         def __get__(self):
@@ -1023,7 +1030,8 @@ cdef class _ScatteringSystemAtOmegaK:
             self.sswk.eta = eta
     
     @boundscheck(False)
-    def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS):
+    def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS,
+            ewaldparts = EwaldPart.FULL):
         """Evaluate electric field for a given excitation coefficient vector (periodic system)
 
         Parameters
@@ -1042,6 +1050,7 @@ cdef class _ScatteringSystemAtOmegaK:
         if(btyp != QPMS_HANKEL_PLUS):
             raise NotImplementedError("Only first kind Bessel function-based fields are supported")
         cdef qpms_bessel_t btyp_c = BesselType(btyp)
+        cdef qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
         evalpos = np.array(evalpos, dtype=float, copy=False)
         if evalpos.shape[-1] != 3:
             raise ValueError("Last dimension of evalpos has to be 3")
@@ -1052,18 +1061,19 @@ cdef class _ScatteringSystemAtOmegaK:
         cdef ccart3_t res
         cdef cart3_t pos
         cdef Py_ssize_t i
-        with nogil, wraparound(False), parallel():
-            for i in prange(evalpos_a.shape[0]):
+        with nogil:
+            for i in range(evalpos_a.shape[0]):
                 pos.x = evalpos_a[i,0]
                 pos.y = evalpos_a[i,1]
                 pos.z = evalpos_a[i,2]
-                res = qpms_scatsyswk_scattered_E(&self.sswk, btyp_c, &scv_view[0], pos)
+                res = qpms_scatsyswk_scattered_E_e(&self.sswk, btyp_c, &scv_view[0], pos, ewaldparts_c)
                 results[i,0] = res.x
                 results[i,1] = res.y
                 results[i,2] = res.z
         return results.reshape(evalpos.shape)
 
-    def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS):
+    @boundscheck(False)
+    def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS, ewaldparts=EwaldPart.FULL):
         # TODO examples
         """Evaluate scattered field "basis" (periodic system)
 
@@ -1085,6 +1095,7 @@ cdef class _ScatteringSystemAtOmegaK:
         if(btyp != QPMS_HANKEL_PLUS):
             raise NotImplementedError("Only first kind Bessel function-based fields are supported")
         cdef qpms_bessel_t btyp_c = BesselType(btyp)
+        cdef qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
         cdef Py_ssize_t fecv_size = self.fecv_size
         evalpos = np.array(evalpos, dtype=float, copy=False)
         if evalpos.shape[-1] != 3:
@@ -1092,20 +1103,20 @@ cdef class _ScatteringSystemAtOmegaK:
         cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
         cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], fecv_size, 3), dtype=complex)
         cdef ccart3_t *res
-        res = <ccart3_t *> malloc(fecv_size*sizeof(ccart3_t))
         cdef cart3_t pos
         cdef Py_ssize_t i, j
-        with nogil, wraparound(False), parallel():
-            for i in prange(evalpos_a.shape[0]):
+        with nogil:
+            res = <ccart3_t *> malloc(fecv_size*sizeof(ccart3_t))
+            for i in range(evalpos_a.shape[0]):
                 pos.x = evalpos_a[i,0]
                 pos.y = evalpos_a[i,1]
                 pos.z = evalpos_a[i,2]
-                qpms_scatsyswk_scattered_field_basis(res, &self.sswk, btyp_c, pos)
+                qpms_scatsyswk_scattered_field_basis_e(res, &self.sswk, btyp_c, pos, ewaldparts_c)
                 for j in range(fecv_size):
                     results[i,j,0] = res[j].x
                     results[i,j,1] = res[j].y
                     results[i,j,2] = res[j].z
-        free(res)
+            free(res)
         return results.reshape(evalpos.shape[:-1] + (self.fecv_size, 3))
     
     property fecv_size: 
@@ -1117,9 +1128,6 @@ cdef class _ScatteringSystemAtOmega:
     that keeps the T-matrix and background data evaluated
     at specific frequency.
     '''
-    cdef qpms_scatsys_at_omega_t *ssw
-    cdef ScatteringSystem ss_pyref
-    
     def check(self): # cdef instead?
         if not self.ssw:
             raise ValueError("_ScatteringSystemAtOmega's ssw-pointer not set. You must not use the default constructor; ScatteringSystem.create() instead")
@@ -1148,9 +1156,6 @@ cdef class _ScatteringSystemAtOmega:
         qpms_scatsysw_apply_Tmatrices_full(&target_view[0], &a_view[0], self.ssw)
         return target_np
     
-    cdef qpms_scatsys_at_omega_t *rawpointer(self):
-        return self.ssw
-
     def scatter_solver(self, iri=None, k=None):
         self.check()
         cdef _ScatteringSystemAtOmegaK sswk # used only for periodic systems
@@ -1234,7 +1239,6 @@ cdef class _ScatteringSystemAtOmega:
 
     @boundscheck(False)
     def scattered_E(self, scatcoeffvector_full, evalpos, blochvector=None, btyp=QPMS_HANKEL_PLUS, bint alt=False): 
-        # FIXME TODO this obviously does not work for periodic systems
         """Evaluate electric field for a given excitation coefficient vector
 
         Parameters
@@ -1271,19 +1275,85 @@ cdef class _ScatteringSystemAtOmega:
         cdef ccart3_t res
         cdef cart3_t pos
         cdef Py_ssize_t i
-        with wraparound(False), nogil, parallel():
-            for i in prange(evalpos_a.shape[0]):
+        for i in range(evalpos_a.shape[0]):
+            pos.x = evalpos_a[i,0]
+            pos.y = evalpos_a[i,1]
+            pos.z = evalpos_a[i,2]
+            if alt:
+                res = qpms_scatsysw_scattered_E__alt(self.ssw, btyp_c, &scv_view[0], pos)
+            else:
+                res = qpms_scatsysw_scattered_E(self.ssw, btyp_c, &scv_view[0], pos)
+            results[i,0] = res.x
+            results[i,1] = res.y
+            results[i,2] = res.z
+        return results.reshape(evalpos.shape)
+
+    @boundscheck(False)
+    def scattered_field_basis(self, evalpos, blochvector=None, particle=None, btyp=QPMS_HANKEL_PLUS):
+        # TODO examples
+        # FIXME periodic case not implemented
+        """Evaluate scattered field "basis" 
+
+        This function enables the evaluation of "scattered" fields
+        generated by the system for many different excitation 
+        coefficients vectors, without the expensive re-evaluation of the respective
+        translation operators for each excitation coefficient vector.
+
+        Parameters
+        ----------
+        evalpos: array_like of floats and shape (..., 3)
+            Evaluation points in cartesian coordinates.
+        blochvector: array_like or None
+            Bloch vector, must be supplied (non-None) for periodic systems, else None.
+        particle_index: int or None (default), optional
+            A valid particle index; if specified, only the contribution of the given particle
+            is evaluated. Not applicable for periodic arrays.
+        btyp: BesselType, optional
+            Kind of the waves. Defaults to BesselType.HANKEL_PLUS.
+
+        Returns
+        ------- 
+        ndarray of complex, with the shape `evalpos.shape[:-1] + (n, 3)`
+            "Basis" fields at the positions given in `evalpos`, in cartesian coordinates.
+            `n` here stays either for `self.fecv_size` (if `particle==None`)
+            or `len(self.bspec_pi(particle))`.
+        """
+        #if(btyp != QPMS_HANKEL_PLUS): # #TODO, IIRC not supported only for periodic systems
+        #    raise NotImplementedError("Only first kind Bessel function-based fields are supported")
+        cdef qpms_bessel_t btyp_c = BesselType(btyp)
+        cdef Py_ssize_t basissize = self.fecv_size if particle is None else len(self.bspec_pi(particle))
+        cdef qpms_ss_pi_t pi = particle
+        evalpos = np.array(evalpos, dtype=float, copy=False)
+        if evalpos.shape[-1] != 3:
+            raise ValueError("Last dimension of evalpos has to be 3")
+        cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
+        cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], basissize, 3), dtype=complex)
+        cdef ccart3_t *res
+        cdef cart3_t pos
+        cdef Py_ssize_t i, j
+        with nogil:
+            res = <ccart3_t *> malloc(basissize*sizeof(ccart3_t)) # thread-local
+            for i in range(evalpos_a.shape[0]):
                 pos.x = evalpos_a[i,0]
                 pos.y = evalpos_a[i,1]
                 pos.z = evalpos_a[i,2]
-                if alt:
-                    res = qpms_scatsysw_scattered_E__alt(self.ssw, btyp_c, &scv_view[0], pos)
+                if particle is None:
+                    qpms_scatsysw_scattered_field_basis(res, self.ssw, btyp_c, pos)
                 else:
-                    res = qpms_scatsysw_scattered_E(self.ssw, btyp_c, &scv_view[0], pos)
-                results[i,0] = res.x
-                results[i,1] = res.y
-                results[i,2] = res.z
-        return results.reshape(evalpos.shape)
+                    qpms_scatsysw_scattered_field_basis_pi(res, self.ssw, pi, btyp_c, pos)
+                for j in range(basissize):
+                    results[i,j,0] = res[j].x
+                    results[i,j,1] = res[j].y
+                    results[i,j,2] = res[j].z
+            free(res)
+        return results.reshape(evalpos.shape[:-1] + (basissize, 3))
+
+    def bspec_pi(self, pi):
+        return self.ss_pyref.bspec_pi(pi)
+
+    property bspecs:
+        def __get__(self):
+            return self.ss_pyref.bspecs
 
 
 cdef class ScatteringMatrix:

qpms/qpms_cdefs.pxd

@@ -147,6 +147,11 @@ cdef extern from "qpms_types.h":
         cdouble mu
     ctypedef enum qpms_coord_system_t:
         pass
+    ctypedef enum qpms_ewald_part:
+        QPMS_EWALD_LONG_RANGE
+        QPMS_EWALD_SHORT_RANGE
+        QPMS_EWALD_FULL
+        QPMS_EWALD_0TERM
     # maybe more if needed
 
 cdef extern from "qpms_error.h":
@@ -185,10 +190,13 @@ cdef extern from "vswf.h":
     csphvec_t qpms_vswf_single_mg_csph(qpms_m_t m, qpms_l_t n, csph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm)
 
 cdef extern from "indexing.h":
-    qpms_y_t qpms_lMax2nelem(qpms_l_t lMax)
-    qpms_uvswfi_t qpms_tmn2uvswfi(qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n)
-    qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t* t, qpms_m_t* m, qpms_l_t* n)
-    qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u)
+    qpms_y_t qpms_lMax2nelem(qpms_l_t lMax) nogil
+    qpms_uvswfi_t qpms_tmn2uvswfi(qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n) nogil
+    qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t* t, qpms_m_t* m, qpms_l_t* n) nogil
+    qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) nogil
+    qpms_y_t qpms_lMax2nelem_sc(qpms_l_t lmax) nogil
+    void qpms_y2mn_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n)
+    qpms_l_t qpms_nelem2lMax(qpms_y_t nelem)
     # maybe more if needed
 
 # Point generators from lattices.h
@@ -251,95 +259,95 @@ cdef extern from "lattices.h":
 
 
 cdef extern from "vectors.h":
-    cart2_t cart2_add(const cart2_t a, const cart2_t b)
-    cart2_t cart2_substract(const cart2_t a, const cart2_t b)
-    cart2_t cart2_scale(const double c, const cart2_t v)
-    double cart2_dot(const cart2_t a, const cart2_t b)
-    double cart2_normsq(const cart2_t a)
-    double cart2norm(const cart2_t v)
-    pol_t cart2pol(const cart2_t cart)
-    sph_t pol2sph_equator(const pol_t pol)
-    sph_t cart22sph(const cart2_t cart)
-    sph_t cart12sph_zaxis(double z)
-    cart3_t cart12cart3z(double z)
-    cart3_t cart22cart3xy(const cart2_t a)
-    cart2_t cart3xy2cart2(const cart3_t a)
-    double cart3_dot(const cart3_t a, const cart3_t b)
-    cart3_t cart3_vectorprod(const cart3_t a, const cart3_t b)
-    double cart3_tripleprod(const cart3_t a, const cart3_t b, const cart3_t c)
-    double cart3_normsq(const cart3_t a)
-    double cart3norm(const cart3_t v)
-    sph_t cart2sph(const cart3_t cart)
-    cart3_t sph2cart(const sph_t sph)
-    cart2_t pol2cart(const pol_t pol)
-    cart3_t pol2cart3_equator(const pol_t pol)
-    cart3_t cart3_add(const cart3_t a, const cart3_t b)
-    cart3_t cart3_substract(const cart3_t a, const cart3_t b)
-    cart3_t cart3_scale(const double c, const cart3_t v)
-    cart3_t cart3_divscale( const cart3_t v, const double c)
-    double cart3_dist(const cart3_t a, const cart3_t b)
-    bint cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol)
-    ccart3_t ccart3_scale(const cdouble c, const ccart3_t v)
-    ccart3_t ccart3_add(const ccart3_t a, const ccart3_t b)
-    ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b)
-    cdouble ccart3_dotnc(const ccart3_t a, const ccart3_t b)
-    ccart3_t cart32ccart3(cart3_t c)
-    csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b)
-    csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b)
-    csphvec_t csphvec_scale(cdouble c, const csphvec_t v)
-    cdouble csphvec_dotnc(const csphvec_t a, const csphvec_t b)
-    sph_t sph_scale(double c, const sph_t s)
-    csph_t sph_cscale(cdouble c, const sph_t s)
-    ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at)
-    ccart3_t csphvec2ccart_csph(const csphvec_t sphvec, const csph_t at)
-    csphvec_t ccart2csphvec(const ccart3_t cartvec, const sph_t at)
-    csph_t sph2csph(sph_t s)
-    sph_t csph2sph(csph_t s)
-    csph_t ccart2csph(const ccart3_t cart)
-    csph_t cart2csph(const cart3_t cart)
-    ccart3_t csph2ccart(const csph_t sph)
-    void ccart3_kahaninit(ccart3_t *sum, ccart3_t *compensation)
-    void csphvec_kahaninit(csphvec_t *sum, csphvec_t *compensation)
-    void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const ccart3_t input)
-    void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t input)
-    double csphvec_norm(const csphvec_t a)
-    double csphvec_reldiff_abstol(const csphvec_t a, const csphvec_t b, double tolerance)
-    double csphvec_reldiff(const csphvec_t a, const csphvec_t b)
-    sph_t anycoord2sph(anycoord_point_t p, qpms_coord_system_t t)
-    cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t)
-    double anycoord_norm(anycoord_point_t p, qpms_coord_system_t t)
-    pol_t anycoord2pol(anycoord_point_t p, qpms_coord_system_t t)
-    cart2_t anycoord2cart2(anycoord_point_t p, qpms_coord_system_t t)
-    double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t)
+    cart2_t cart2_add(const cart2_t a, const cart2_t b) nogil
+    cart2_t cart2_substract(const cart2_t a, const cart2_t b) nogil
+    cart2_t cart2_scale(const double c, const cart2_t v) nogil
+    double cart2_dot(const cart2_t a, const cart2_t b) nogil
+    double cart2_normsq(const cart2_t a) nogil
+    double cart2norm(const cart2_t v) nogil
+    pol_t cart2pol(const cart2_t cart) nogil
+    sph_t pol2sph_equator(const pol_t pol) nogil
+    sph_t cart22sph(const cart2_t cart) nogil
+    sph_t cart12sph_zaxis(double z) nogil
+    cart3_t cart12cart3z(double z) nogil
+    cart3_t cart22cart3xy(const cart2_t a) nogil
+    cart2_t cart3xy2cart2(const cart3_t a) nogil
+    double cart3_dot(const cart3_t a, const cart3_t b) nogil
+    cart3_t cart3_vectorprod(const cart3_t a, const cart3_t b) nogil
+    double cart3_tripleprod(const cart3_t a, const cart3_t b, const cart3_t c) nogil
+    double cart3_normsq(const cart3_t a) nogil
+    double cart3norm(const cart3_t v) nogil
+    sph_t cart2sph(const cart3_t cart) nogil
+    cart3_t sph2cart(const sph_t sph) nogil
+    cart2_t pol2cart(const pol_t pol) nogil
+    cart3_t pol2cart3_equator(const pol_t pol) nogil
+    cart3_t cart3_add(const cart3_t a, const cart3_t b) nogil
+    cart3_t cart3_substract(const cart3_t a, const cart3_t b) nogil
+    cart3_t cart3_scale(const double c, const cart3_t v) nogil
+    cart3_t cart3_divscale( const cart3_t v, const double c) nogil
+    double cart3_dist(const cart3_t a, const cart3_t b) nogil
+    bint cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) nogil
+    ccart3_t ccart3_scale(const cdouble c, const ccart3_t v) nogil
+    ccart3_t ccart3_add(const ccart3_t a, const ccart3_t b) nogil
+    ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) nogil
+    cdouble ccart3_dotnc(const ccart3_t a, const ccart3_t b) nogil
+    ccart3_t cart32ccart3(cart3_t c) nogil
+    csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b) nogil
+    csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b) nogil
+    csphvec_t csphvec_scale(cdouble c, const csphvec_t v) nogil
+    cdouble csphvec_dotnc(const csphvec_t a, const csphvec_t b) nogil
+    sph_t sph_scale(double c, const sph_t s) nogil
+    csph_t sph_cscale(cdouble c, const sph_t s) nogil
+    ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at) nogil
+    ccart3_t csphvec2ccart_csph(const csphvec_t sphvec, const csph_t at) nogil
+    csphvec_t ccart2csphvec(const ccart3_t cartvec, const sph_t at) nogil
+    csph_t sph2csph(sph_t s) nogil
+    sph_t csph2sph(csph_t s) nogil
+    csph_t ccart2csph(const ccart3_t cart) nogil
+    csph_t cart2csph(const cart3_t cart) nogil
+    ccart3_t csph2ccart(const csph_t sph) nogil
+    void ccart3_kahaninit(ccart3_t *sum, ccart3_t *compensation) nogil
+    void csphvec_kahaninit(csphvec_t *sum, csphvec_t *compensation) nogil
+    void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const ccart3_t input) nogil
+    void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t input) nogil
+    double csphvec_norm(const csphvec_t a) nogil
+    double csphvec_reldiff_abstol(const csphvec_t a, const csphvec_t b, double tolerance) nogil
+    double csphvec_reldiff(const csphvec_t a, const csphvec_t b) nogil
+    sph_t anycoord2sph(anycoord_point_t p, qpms_coord_system_t t) nogil
+    cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t) nogil
+    double anycoord_norm(anycoord_point_t p, qpms_coord_system_t t) nogil
+    pol_t anycoord2pol(anycoord_point_t p, qpms_coord_system_t t) nogil
+    cart2_t anycoord2cart2(anycoord_point_t p, qpms_coord_system_t t) nogil
+    double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t) nogil
     void qpms_array_coord_transform(void *dest, qpms_coord_system_t tdest,
-            const void *src, qpms_coord_system_t tsrc, size_t nmemb)
+            const void *src, qpms_coord_system_t tsrc, size_t nmemb) nogil
     void anycoord_arr2something(void *dest, qpms_coord_system_t tdest,
-            const void *src, qpms_coord_system_t tsrc, size_t nmemb)
-    void cart3_to_double_array(double *a, cart3_t b)
-    cart3_t cart3_from_double_array(const double *a)
-    void cart2_to_double_array(double *a, cart2_t b)
-    cart2_t cart2_from_double_array(const double *a)
+            const void *src, qpms_coord_system_t tsrc, size_t nmemb) nogil
+    void cart3_to_double_array(double *a, cart3_t b) nogil
+    cart3_t cart3_from_double_array(const double *a) nogil
+    void cart2_to_double_array(double *a, cart2_t b) nogil
+    cart2_t cart2_from_double_array(const double *a) nogil
     const cart2_t CART2_ZERO
 
 
 cdef extern from "quaternions.h":
-    qpms_quat_t qpms_quat_2c_from_4d(qpms_quat4d_t q)
-    qpms_quat4d_t qpms_quat_4d_from_2c(qpms_quat_t q)
-    qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q)
-    qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q)
-    qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q) 
-    qpms_quat_t qpms_quat_conj(qpms_quat_t q) 
-    double qpms_quat_norm(qpms_quat_t q) 
-    double qpms_quat_imnorm(qpms_quat_t q)
-    qpms_quat_t qpms_quat_normalise(qpms_quat_t q) 
-    qpms_quat_t qpms_quat_log(qpms_quat_t q)
-    qpms_quat_t qpms_quat_exp(qpms_quat_t q)
-    qpms_quat_t qpms_quat_pow(qpms_quat_t q, double exponent)
+    qpms_quat_t qpms_quat_2c_from_4d(qpms_quat4d_t q) nogil
+    qpms_quat4d_t qpms_quat_4d_from_2c(qpms_quat_t q) nogil
+    qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q) nogil
+    qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q) nogil
+    qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q) nogil 
+    qpms_quat_t qpms_quat_conj(qpms_quat_t q) nogil 
+    double qpms_quat_norm(qpms_quat_t q) nogil 
+    double qpms_quat_imnorm(qpms_quat_t q) nogil
+    qpms_quat_t qpms_quat_normalise(qpms_quat_t q) nogil 
+    qpms_quat_t qpms_quat_log(qpms_quat_t q) nogil
+    qpms_quat_t qpms_quat_exp(qpms_quat_t q) nogil
+    qpms_quat_t qpms_quat_pow(qpms_quat_t q, double exponent) nogil
     cdouble qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
-                           qpms_m_t mp, qpms_m_t m)
-    qpms_irot3_t qpms_irot3_mult(qpms_irot3_t p, qpms_irot3_t q)
-    qpms_irot3_t qpms_irot3_pow(qpms_irot3_t p, int n)
-    qpms_quat_t qpms_quat_from_rotvector(cart3_t v)
+                           qpms_m_t mp, qpms_m_t m) nogil
+    qpms_irot3_t qpms_irot3_mult(qpms_irot3_t p, qpms_irot3_t q) nogil
+    qpms_irot3_t qpms_irot3_pow(qpms_irot3_t p, int n) nogil
+    qpms_quat_t qpms_quat_from_rotvector(cart3_t v) nogil
 
 cdef extern from "groups.h":
     struct qpms_finite_group_irrep_t:
@@ -487,6 +495,7 @@ cdef extern from "tmatrices.h":
     qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, cdouble omega, const void *params)
     qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t, cdouble omega, const void *params)
     qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, cdouble omega, const void *params)
+    qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t, cdouble omega, const void *params)
     struct qpms_tmatrix_generator_sphere_param_t:
         qpms_epsmu_generator_t outside
         qpms_epsmu_generator_t inside
@@ -620,6 +629,7 @@ cdef extern from "scatsystem.h":
         size_t *saecv_sizes
         const qpms_finite_group_t *sym
         qpms_scatsys_periodic_info_t per
+        qpms_trans_calculator *c
 
         # per[] and other stuff not currently needed in cython
     void qpms_scatsys_free(qpms_scatsys_t *s)
@@ -635,6 +645,11 @@ cdef extern from "scatsystem.h":
             cdouble omega, const qpms_tolerance_spec_t *tol)
     qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, cdouble omega)
     void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw)
+    cdouble *qpms_scatsys_irrep_transform_matrix(cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri)
+    cdouble *qpms_scatsys_irrep_pack_matrix_stupid(cdouble *target_packed,
+            const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
+    cdouble *qpms_scatsys_irrep_unpack_matrix_stupid(cdouble *target_full, 
+            const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add)
     cdouble *qpms_scatsys_irrep_pack_matrix(cdouble *target_packed,
             const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
     cdouble *qpms_scatsys_irrep_unpack_matrix(cdouble *target_full, 
@@ -705,8 +720,20 @@ cdef extern from "scatsystem.h":
             const cdouble *f_excitation_vector_full, cart3_t where) nogil
     ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp,
             const cdouble *f_excitation_vector_full, cart3_t where) nogil
+    ccart3_t qpms_scatsyswk_scattered_E_e(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp,
+            const cdouble *f_excitation_vector_full, cart3_t where, qpms_ewald_part parts) nogil
+    qpms_errno_t qpms_scatsys_scattered_field_basis(ccart3_t *target, const qpms_scatsys_t *ss,
+            qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil
+    qpms_errno_t qpms_scatsysw_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
+            qpms_bessel_t btyp, cart3_t where) nogil
+    qpms_errno_t qpms_scatsys_scattered_field_basis_pi(ccart3_t *target, const qpms_scatsys_t *ss,
+            qpms_ss_pi_t pi, qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil
+    qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
+            qpms_ss_pi_t pi, qpms_bessel_t btyp, cart3_t where) nogil
     qpms_errno_t qpms_scatsyswk_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
             qpms_bessel_t btyp, cart3_t where) nogil
+    qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
+            qpms_bessel_t btyp, cart3_t where, qpms_ewald_part parts) nogil
     double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector) nogil
 
 
@@ -715,12 +742,6 @@ cdef extern from "ewald.h":
         cdouble val
         double err
 
-    ctypedef enum qpms_ewald_part:
-        QPMS_EWALD_LONG_RANGE
-        QPMS_EWALD_SHORT_RANGE
-        QPMS_EWALD_FULL
-        QPMS_EWALD_0TERM
-
     struct qpms_ewald3_constants_t:
         qpms_l_t lMax
         qpms_y_t nelem_sc

qpms/qpms_error.h

@@ -5,6 +5,9 @@
  */
 #ifndef QPMS_ERROR_H
 #define QPMS_ERROR_H
+#ifdef __cplusplus
+extern "C" {
+#endif
 
 #include "optim.h"
 
@@ -266,4 +269,7 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
 	}\
 }
 
+#ifdef __cplusplus
+}
+#endif
 #endif 

qpms/qpms_specfunc.h

@@ -3,6 +3,10 @@
  */
 #ifndef QPMS_SPECFUNC_H
 #define QPMS_SPECFUNC_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #include "qpms_types.h"
 #include <gsl/gsl_sf_legendre.h>
 
@@ -11,25 +15,25 @@
  ******************************************************************************/
 
 // TODO unify types
-qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex double x, complex double *result_array);
+qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, _Complex double x, _Complex double *result_array);
 
 
 
 typedef struct {
         qpms_l_t lMax;
         double *akn; // coefficients as in DLMF 10.49.1
-	//complex double *bkn; // coefficients of the derivatives
+	//_Complex double *bkn; // coefficients of the derivatives
 } qpms_sbessel_calculator_t;
 
 qpms_sbessel_calculator_t *qpms_sbessel_calculator_init(void);
 void qpms_sbessel_calculator_pfree(qpms_sbessel_calculator_t *c);
 
 qpms_errno_t qpms_sbessel_calc_fill(qpms_sbessel_calculator_t *c, qpms_bessel_t typ, qpms_l_t lmax,
-	double x, complex double *result_array);
+	double x, _Complex double *result_array);
 
-complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, complex double x);
+_Complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, _Complex double x);
 qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t *c, qpms_l_t lmax,
-		complex double x, complex double *result_array);
+		_Complex double x, _Complex double *result_array);
 
 
 /******************************************************************************
@@ -107,4 +111,7 @@ double qpms_legendre0(int m, int n);
 // Associated Legendre polynomial derivative at zero argument (DLMF 14.5.2)
 double qpms_legendred0(int m, int n);
 
+#ifdef __cplusplus
+}
+#endif
 #endif // QPMS_SPECFUNC_H

qpms/qpms_types.h

@@ -3,10 +3,13 @@
  */
 #ifndef QPMS_TYPES_H
 #define QPMS_TYPES_H
-#include <complex.h>
-#include <stdbool.h>
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #include <stddef.h>
 #include <stdint.h>
+#include <stdbool.h>
 
 #ifndef M_PI_2
 #define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487L)
@@ -99,7 +102,9 @@ typedef enum {
  *  These bit flags are used by the functions declared in normalisation.h
  *  that return the appropriate convention-dependent factors.
  *
- *  See @ref vswf_conventions for comparison of the various conventions used.
+ *  \see @ref vswf_conventions for comparison of the various conventions used.
+ *  \see @ref swf_conventions_qpms for description how the conventions are used internally and in the QPMS API.
+ *  
  */
 typedef enum {
 	QPMS_NORMALISATION_UNDEF = 0, ///< Convention undefined. This should not happen.
@@ -189,7 +194,7 @@ typedef struct cart3_t {
 
 /// 3D complex (actually 6D) coordinates. See also vectors.h.
 typedef struct ccart3_t {
-	complex double x, y, z;
+	_Complex double x, y, z;
 } ccart3_t;
 
 /// 3D complex vector pair (represents the E, H fields).
@@ -210,13 +215,13 @@ typedef struct sph_t {
 
 /// Spherical coordinates with complex radial component. See also vectors.h.
 typedef struct csph_t { // Do I really need this???
-	complex double r;
+	_Complex double r;
 	double	theta, phi;
 } csph_t;
 
 /// 3D complex vector components in local spherical basis. See also vectors.h.
 typedef struct csphvec_t {
-	complex double rc, thetac, phic; 
+	_Complex double rc, thetac, phic; 
 } csphvec_t;
 
 /// 2D polar coordinates. See also vectors.h.
@@ -259,7 +264,7 @@ typedef enum {
  * See quaternions.h for "methods".
  */
 typedef struct qpms_quat_t {
-        complex double a, b;
+        _Complex double a, b;
 } qpms_quat_t;
 
 /// Quaternion type as four doubles.
@@ -344,7 +349,7 @@ typedef struct qpms_tmatrix_t {
          *  there is also one with order \a −m.
          */
         const qpms_vswf_set_spec_t *spec;
-        complex double *m; ///< Matrix elements in row-major order.
+        _Complex double *m; ///< Matrix elements in row-major order.
         bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free()
 } qpms_tmatrix_t;
 
@@ -420,11 +425,22 @@ typedef struct qpms_abstract_tmatrix_t {
 
 /// A type holding electric permittivity and magnetic permeability of a material.
 typedef struct qpms_epsmu_t {
-	complex double eps; ///< Relative permittivity.
-	complex double mu; ///< Relative permeability.
+	_Complex double eps; ///< Relative permittivity.
+	_Complex double mu; ///< Relative permeability.
 } qpms_epsmu_t;
 
 struct qpms_tolerance_spec_t; // See tolerances.h
 
+typedef enum {
+	QPMS_EWALD_LONG_RANGE = 1,
+	QPMS_EWALD_SHORT_RANGE = 2,
+	QPMS_EWALD_0TERM = 4,
+	QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
+} qpms_ewald_part;
+
 #define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
+
+#ifdef __cplusplus
+}
+#endif
 #endif // QPMS_TYPES

qpms/qpmsblas.h

@@ -14,6 +14,10 @@
  */
 #ifndef QPMSBLAS_H
 #define QPMSBLAS_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #define QPMS_BLAS_INDEX_T long long int
 
 #ifndef CBLAS_H
@@ -31,4 +35,7 @@ void qpms_zgemm(CBLAS_LAYOUT Order, CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE Tran
 		const _Complex double *B, const QPMS_BLAS_INDEX_T ldb,
 		const _Complex double *beta, _Complex double *C, const QPMS_BLAS_INDEX_T ldc);
 
+#ifdef __cplusplus
+}
+#endif
 #endif //QPMSBLAS_H

qpms/quaternions.h

@@ -3,6 +3,9 @@
  */
 #ifndef QPMS_WIGNER_H
 #define QPMS_WIGNER_H
+#ifdef __cplusplus
+extern "C" {
+#endif
 
 #include "qpms_types.h"
 #include "vectors.h"
@@ -93,7 +96,7 @@ static inline double qpms_quat_norm(const qpms_quat_t q) {
 }
 
 /// Test approximate equality of quaternions.
-static inline bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, double atol) {
+static inline _Bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, double atol) {
 	return qpms_quat_norm(qpms_quat_sub(p,q)) <= atol;
 }
 
@@ -119,7 +122,7 @@ static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p, double atol) {
 }		
 
 /// Test approximate equality of "standardised" quaternions, so that \f$-q\f$ is considered equal to \f$q\f$.
-static inline bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) {
+static inline _Bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) {
 	return qpms_quat_norm(qpms_quat_sub(
 				qpms_quat_standardise(p, atol),
 				qpms_quat_standardise(q, atol))) <= atol;
@@ -202,17 +205,17 @@ static inline qpms_quat_t qpms_quat_from_rotvector(cart3_t v) {
  * The D matrix are calculated using formulae (3), (4), (6), (7) from
  * http://moble.github.io/spherical_functions/WignerDMatrices.html
  */
-complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
+_Complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
 	       qpms_m_t	mp, qpms_m_t m);
 
 /// A VSWF representation element of the O(3) group.
 /**
  * TODO more doc.
  */
-complex double qpms_vswf_irot_elem_from_irot3(
+_Complex double qpms_vswf_irot_elem_from_irot3(
 		const qpms_irot3_t q, ///< The O(3) element in the quaternion representation.
 		qpms_l_t l, qpms_m_t mp, qpms_m_t m,
-		bool pseudo ///< Determines the sign of improper rotations. True for magnetic waves, false otherwise.
+		_Bool pseudo ///< Determines the sign of improper rotations. True for magnetic waves, false otherwise.
 		);
 
 
@@ -251,7 +254,7 @@ static inline qpms_irot3_t qpms_irot3_pow(const qpms_irot3_t p, int n) {
 }
 
 /// Test approximate equality of irot3.
-static inline bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) {
+static inline _Bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) {
 	return qpms_quat_isclose2(p.rot, q.rot, atol) && p.det == q.det;
 }
 
@@ -303,4 +306,7 @@ static inline qpms_irot3_t qpms_irot3_zrot_Nk(double N, double k) {
 	return qpms_irot3_zrot_angle(2 * M_PI * k / N);
 }
 
+#ifdef __cplusplus
+}
+#endif
 #endif //QPMS_WIGNER_H

qpms/scatsystem.c

@@ -24,6 +24,7 @@
 #include "tolerances.h"
 #include "beyn.h"
 #include "tiny_inlines.h"
+#include "lattices.h"
 
 #ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS
 #include "qpmsblas.h"
@@ -1295,7 +1296,6 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f
     QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
   complex double *tmp; // translation matrix, S or W
   QPMS_CRASHING_MALLOC(tmp, SQ(ss->max_bspecn) * sizeof(complex double));
-  memset(target, 0, SQ(full_len) * sizeof(complex double));
   const complex double zero = 0, minusone = -1;
   { // Non-diagonal part; M[piR, piC] = -T[piR] S(piR<-piC)
     size_t fullvec_offsetR = 0;
@@ -1321,7 +1321,14 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f
                 QPMS_EWALD_FULL, eta));
         }
 
-        cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+        if (k == NULL && piC == piR) { 
+          // non-periodic case diagonal block: no "self-interaction", just
+          // fill with zeroes (the ones on the diagonal are added in the very end)
+          QPMS_ENSURE_SUCCESS(LAPACKE_zlaset(CblasRowMajor, 'x',
+                bspecR->n /*m*/, bspecC->n/*n*/, 0 /*alfa: offdiag*/, 0 /*beta: diag*/, 
+                target + fullvec_offsetR*full_len + fullvec_offsetC /*a*/,
+                full_len /*lda*/));
+        } else cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
             bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, 
             &minusone/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/,
             tmp/*b*/, bspecC->n/*ldb*/, &zero/*beta*/,
@@ -2043,6 +2050,61 @@ ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw,
       cvf, where);
 }
 
+qpms_errno_t qpms_scatsys_scattered_field_basis_pi(
+    ccart3_t *target,
+    const qpms_scatsys_t *ss,
+    const qpms_ss_pi_t pi,
+    const qpms_bessel_t btyp,
+    const complex double k,
+    const cart3_t where
+    ) {
+  qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsyswk_scattered_field_basis()");
+  QPMS_UNTESTED;
+  const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi);
+  csphvec_t *vswfs_sph; //Single particle contributions in spherical coordinates
+  QPMS_CRASHING_MALLOC(vswfs_sph, bspec->n * sizeof(*vswfs_sph));
+  const cart3_t particle_pos = ss->p[pi].pos;
+  const csph_t kr = sph_cscale(k, cart2sph(
+      cart3_substract(where, particle_pos)));
+  QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(vswfs_sph, bspec, kr, btyp));
+  for(size_t i = 0; i < bspec->n; ++i)
+      target[i] = csphvec2ccart_csph(vswfs_sph[i], kr);
+  
+  free(vswfs_sph);
+  return QPMS_SUCCESS;
+}
+
+qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(
+    ccart3_t *target, const qpms_scatsys_at_omega_t *ssw, const qpms_ss_pi_t pi,
+    const qpms_bessel_t btyp, const cart3_t where) {
+  return qpms_scatsys_scattered_field_basis_pi(target, ssw->ss, pi, btyp,
+      ssw->wavenumber, where);
+}
+
+qpms_errno_t qpms_scatsys_scattered_field_basis(
+    ccart3_t *target,
+    const qpms_scatsys_t *ss,
+    const qpms_bessel_t btyp,
+    const complex double k,
+    const cart3_t where
+    ) {
+  qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsyswk_scattered_field_basis()");
+  QPMS_UNTESTED;
+  for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) 
+    QPMS_ENSURE_SUCCESS(qpms_scatsys_scattered_field_basis_pi(
+          target + ss->fecv_pstarts[pi], ss, pi, btyp, k, where));
+  return QPMS_SUCCESS;
+}
+
+qpms_errno_t qpms_scatsysw_scattered_field_basis(
+    ccart3_t *target, const qpms_scatsys_at_omega_t *ssw, 
+    const qpms_bessel_t btyp, const cart3_t where) {
+  return qpms_scatsys_scattered_field_basis(target, ssw->ss, btyp, 
+        ssw->wavenumber, where);
+}
+
+
+const qpms_uvswfi_t ELDIPILIST[] = {6, 10, 14};
 
 #define DIPSPECN 3 // We have three basis vectors
 // Evaluates the regular electric dipole waves in the origin. The returned
@@ -2052,11 +2114,15 @@ static inline const qpms_vswf_set_spec_t qpms_fill_regdipoles_0(
   // bspec containing only electric dipoles
   const qpms_vswf_set_spec_t dipspec = {
     .n = DIPSPECN,
+#if  0
     .ilist = (qpms_uvswfi_t[]){
       qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, -1, 1),
       qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC,  0, 1),
       qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1),
     },
+#else 
+    .ilist = ELDIPILIST,
+#endif
     .lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1,
     .capacity=0,
     .norm = normalisation,
@@ -2121,10 +2187,11 @@ ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw,
 
 // For periodic lattices, we use directly the "alternative" implementation, 
 // using translation operator and regular dipole waves at zero
-ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
+ccart3_t qpms_scatsyswk_scattered_E_e(const qpms_scatsys_at_omega_k_t *sswk,
     qpms_bessel_t btyp,
     const complex double *cvf, 
-    const cart3_t where
+    const cart3_t where,
+    const qpms_ewald_part parts
     ) {
   QPMS_UNTESTED;
   if (btyp != QPMS_HANKEL_PLUS) 
@@ -2156,13 +2223,13 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
     const double maxR = sqrt(ss->per.unitcell_volume) * 64;
     const double maxK = 2048 * 2 * M_PI / maxR;
 
-    QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
+    QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32_e(
           ss->c, s, NULL, 
           &dipspec, 1, bspec, dipspec.n,
           sswk->eta, sswk->ssw->wavenumber,
           cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
           cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/, 
-          maxR, maxK));
+          maxR, maxK, parts));
 
     for(size_t i = 0; i < bspec->n; ++i)
       for(size_t j = 0; j < dipspec.n; ++j){
@@ -2176,11 +2243,17 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
   return res;
 }
 
-qpms_errno_t qpms_scatsyswk_scattered_field_basis(
+ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
+    qpms_bessel_t btyp, const complex double *cvf, const cart3_t where) {
+  return qpms_scatsyswk_scattered_E_e(sswk, btyp, cvf, where, QPMS_EWALD_FULL);
+}
+
+qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
     ccart3_t *target,
     const qpms_scatsys_at_omega_k_t *sswk,
     const qpms_bessel_t btyp,
-    const cart3_t where
+    const cart3_t where,
+    const qpms_ewald_part parts
     ) {
   QPMS_UNTESTED;
   if (btyp != QPMS_HANKEL_PLUS) 
@@ -2214,13 +2287,13 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
     const double maxR = sqrt(ss->per.unitcell_volume) * 64;
     const double maxK = 2048 * 2 * M_PI / maxR;
 
-    QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
+    QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32_e(
           ss->c, s, NULL, 
           &dipspec, 1, bspec, dipspec.n,
           sswk->eta, sswk->ssw->wavenumber,
           cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
           cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/, 
-          maxR, maxK));
+          maxR, maxK, parts));
 
     for(size_t i = 0; i < bspec->n; ++i)
       for(size_t j = 0; j < dipspec.n; ++j){
@@ -2236,6 +2309,13 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
   return QPMS_SUCCESS;
 }
 
+qpms_errno_t qpms_scatsyswk_scattered_field_basis(
+    ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
+    const qpms_bessel_t btyp, const cart3_t where) {
+  return qpms_scatsyswk_scattered_field_basis_e(
+      target, sswk, btyp, where, QPMS_EWALD_FULL);
+}
+
 #if 0
 ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
    qpms_iri_t iri, const complex double *cvr, cart3_t where) {

qpms/scatsystem.h

@@ -10,6 +10,10 @@
  */
 #ifndef QPMS_SCATSYSTEM_H
 #define QPMS_SCATSYSTEM_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #include "qpms_types.h"
 #include "vswf.h"
 #include "tmatrices.h"
@@ -107,7 +111,7 @@ typedef struct qpms_ss_orbit_type_t {
 	 *
 	 * TODO doc.
 	 */
-	complex double *irbases;
+	_Complex double *irbases;
 	/// TODO doc.
 	size_t instance_count;
 	/// Cumulative sum of the preceding ot->siza * ot->instance_count;
@@ -243,9 +247,9 @@ typedef struct qpms_scatsys_at_omega_t {
 	 * i.e in the order corresponding to \a ss->tm.
 	 */
 	qpms_tmatrix_t **tm;
-	complex double omega; ///< Angular frequency
+	_Complex double omega; ///< Angular frequency
 	qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency
-	complex double wavenumber; ///< Background medium wave number
+	_Complex double wavenumber; ///< Background medium wave number
 } qpms_scatsys_at_omega_t;
 
 
@@ -291,7 +295,7 @@ typedef struct qpms_scatsys_at_omega_t {
  *
  */
 qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym,
-		complex double omega, const struct qpms_tolerance_spec_t *tol);
+		_Complex double omega, const struct qpms_tolerance_spec_t *tol);
 
 /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load.
 void qpms_scatsys_free(qpms_scatsys_t *s);
@@ -303,50 +307,50 @@ void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw);
 /// Evaluates scattering system T-matrices at a given frequency.
 /** Free the result using qpms_scatsys_at_omega_free() when done. */
 qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss,
-		complex double omega);
+		_Complex double omega);
 
 /// Creates a "full" transformation matrix U that takes a full vector and projects it onto an symmetry adapted basis.
 /** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient.
  * 
  * TODO doc about shape etc.
  */
-complex double *qpms_scatsys_irrep_transform_matrix(complex double *target_U,
+_Complex double *qpms_scatsys_irrep_transform_matrix(_Complex double *target_U,
 		const qpms_scatsys_t *ss, qpms_iri_t iri);
 
 /// Projects a "big" matrix onto an irrep (slow reference implementation).
 /** TODO doc */
-complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed,
-		const complex double *orig_full, const qpms_scatsys_t *ss,
+_Complex double *qpms_scatsys_irrep_pack_matrix_stupid(_Complex double *target_packed,
+		const _Complex double *orig_full, const qpms_scatsys_t *ss,
 		qpms_iri_t iri);
 
 /// Transforms a big "packed" matrix into the full basis (slow reference implementation).
 /** TODO doc */
-complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full,
-		const complex double *orig_packed, const qpms_scatsys_t *ss,
+_Complex double *qpms_scatsys_irrep_unpack_matrix_stupid(_Complex double *target_full,
+		const _Complex double *orig_packed, const qpms_scatsys_t *ss,
 		qpms_iri_t iri, bool add);
 
 /// Projects a "big" matrix onto an irrep.
 /** TODO doc */
-complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
-		const complex double *orig_full, const qpms_scatsys_t *ss,
+_Complex double *qpms_scatsys_irrep_pack_matrix(_Complex double *target_packed,
+		const _Complex double *orig_full, const qpms_scatsys_t *ss,
 		qpms_iri_t iri);
 
 /// Transforms a big "packed" matrix into the full basis.
 /** TODO doc */
-complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
-		const complex double *orig_packed, const qpms_scatsys_t *ss,
+_Complex double *qpms_scatsys_irrep_unpack_matrix(_Complex double *target_full,
+		const _Complex double *orig_packed, const qpms_scatsys_t *ss,
 		qpms_iri_t iri, bool add);
 
 /// Projects a "big" vector onto an irrep.
 /** TODO doc */
-complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
-		const complex double *orig_full, const qpms_scatsys_t *ss,
+_Complex double *qpms_scatsys_irrep_pack_vector(_Complex double *target_packed,
+		const _Complex double *orig_full, const qpms_scatsys_t *ss,
 		qpms_iri_t iri);
 
 /// Transforms a big "packed" vector into the full basis.
 /** TODO doc */
-complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
-		const complex double *orig_packed, const qpms_scatsys_t *ss,
+_Complex double *qpms_scatsys_irrep_unpack_vector(_Complex double *target_full,
+		const _Complex double *orig_packed, const qpms_scatsys_t *ss,
 		qpms_iri_t iri, bool add);
 
 /// Global translation matrix.
@@ -354,31 +358,31 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
  * The diagonal (particle self-) block are filled with zeros (even for regular Bessel waves).
  * This may change in the future.
  */
-complex double *qpms_scatsys_build_translation_matrix_full(
+_Complex double *qpms_scatsys_build_translation_matrix_full(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
-		complex double *target,
+		_Complex double *target,
 		const qpms_scatsys_t *ss,
-		complex double k ///< Wave number to use in the translation matrix.
+		_Complex double k ///< Wave number to use in the translation matrix.
 		);
 
 /// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system.
 /**
  * \returns \a target on success, NULL on error.
  */
-complex double *qpms_scatsysw_build_modeproblem_matrix_full(
+_Complex double *qpms_scatsysw_build_modeproblem_matrix_full(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
-		complex double *target,
+		_Complex double *target,
 		const qpms_scatsys_at_omega_t *ssw
 		);
 
 /// As qpms_scatsys_build_translation_full() but with choice of Bessel function type.
 /** Might be useful for evaluation of cross sections and testing.
  */
-complex double *qpms_scatsys_build_translation_matrix_e_full(
+_Complex double *qpms_scatsys_build_translation_matrix_e_full(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
-		complex double *target,
+		_Complex double *target,
 		const qpms_scatsys_t *ss,
-		complex double k, ///< Wave number to use in the translation matrix.
+		_Complex double k, ///< Wave number to use in the translation matrix.
 		qpms_bessel_t J
 		);
 
@@ -387,12 +391,12 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
  * The diagonal (particle self-) blocks are currently filled with zeros.
  * This may change in the future.
  */
-complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
+_Complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
-		complex double *target,
+		_Complex double *target,
 		const qpms_scatsys_t *ss,
 		qpms_iri_t iri,
-		complex double k, ///< Wave number to use in the translation matrix.
+		_Complex double k, ///< Wave number to use in the translation matrix.
 		qpms_bessel_t J
 		);
 
@@ -400,9 +404,9 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
 /**
  * \returns \a target on success, NULL on error.
  */
-complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
+_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
-		complex double *target,
+		_Complex double *target,
 		const qpms_scatsys_at_omega_t *ssw,
 		qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
 		);
@@ -410,9 +414,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
 /**
  * \returns \a target on success, NULL on error.
  */
-complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
+_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
-		complex double *target,
+		_Complex double *target,
 		const qpms_scatsys_at_omega_t *ssw,
 		qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
 		);
@@ -420,9 +424,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
 /**
  * \returns \a target on success, NULL on error.
  */
-complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
+_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
-		complex double *target,
+		_Complex double *target,
 		const qpms_scatsys_at_omega_t *ssw,
 		qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
 		);
@@ -435,7 +439,7 @@ typedef struct qpms_ss_LU {
 	bool full; ///< true if full matrix; false if irrep-packed.
 	qpms_iri_t iri; ///< Irrep index if `full == false`.
 	/// LU decomposition array.
-	complex double *a;
+	_Complex double *a;
 	/// Pivot index array, size at least max(1,min(m, n)).
 	int *ipiv; 
 } qpms_ss_LU;
@@ -443,14 +447,14 @@ void qpms_ss_LU_free(qpms_ss_LU);
 
 /// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Nonperiodic systems only.
 qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(
-		complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
+		_Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
 		int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
 		const qpms_scatsys_at_omega_t *ssw
 		);
 
 /// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch.
 qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
-		complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
+		_Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
 		int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
 		const qpms_scatsys_at_omega_t *ssw,
 		qpms_iri_t iri
@@ -458,7 +462,7 @@ qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
 
 /// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents.
 qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(
-		complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
+		_Complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
 		int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
 		const qpms_scatsys_at_omega_t *ssw, ///< Must be filled for non-periodic systems.
 		const struct qpms_scatsys_at_omega_k_t *sswk ///< Must be filled for periodic systems, otherwise must be NULL.
@@ -466,16 +470,16 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(
 
 /// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents.
 qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(
-		complex double *modeproblem_matrix_irrep_packed, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
+		_Complex double *modeproblem_matrix_irrep_packed, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
 		int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
 		const qpms_scatsys_at_omega_t *ssw,
 		qpms_iri_t iri
 		);
 
 /// Solves a (possibly partial, irrep-packed) scattering problem \f$ (I-TS)f = Ta_\mathrm{inc} \f$ using a pre-factorised \f$ (I-TS) \f$.
-complex double *qpms_scatsys_scatter_solve(
-		complex double *target_f, ///< Target (full or irrep-packed, depending on `ludata.full`) array for \a f. If NULL, a new one is allocated.
-		const complex double *a_inc, ///< Incident field expansion coefficient vector \a a (full or irrep-packed, depending on `ludata.full`).
+_Complex double *qpms_scatsys_scatter_solve(
+		_Complex double *target_f, ///< Target (full or irrep-packed, depending on `ludata.full`) array for \a f. If NULL, a new one is allocated.
+		const _Complex double *a_inc, ///< Incident field expansion coefficient vector \a a (full or irrep-packed, depending on `ludata.full`).
 		qpms_ss_LU ludata ///< Pre-factorised \f$ I - TS \f$ matrix data.
 		);
 
@@ -495,33 +499,33 @@ typedef struct qpms_scatsys_at_omega_k_t {
 /**
  * \returns \a target on success, NULL on error.
  */
-complex double *qpms_scatsyswk_build_modeproblem_matrix_full(
+_Complex double *qpms_scatsyswk_build_modeproblem_matrix_full(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
-		complex double *target,
+		_Complex double *target,
 		const qpms_scatsys_at_omega_k_t *sswk
 		);
 
 /// Global translation matrix.
-complex double *qpms_scatsys_periodic_build_translation_matrix_full(
+_Complex double *qpms_scatsys_periodic_build_translation_matrix_full(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
-		complex double *target,
+		_Complex double *target,
 		const qpms_scatsys_t *ss,
-		complex double wavenumber, ///< Wave number to use in the translation matrix.
+		_Complex double wavenumber, ///< Wave number to use in the translation matrix.
 		const cart3_t *wavevector, ///< Wavevector / pseudomomentum in cartesian coordinates.
 		double eta ///< Ewald parameter eta. Pass 0 or NaN to use the default value in \a ss.
 		);
 
 /// Global translation matrix. 
-complex double *qpms_scatsyswk_build_translation_matrix_full(
+_Complex double *qpms_scatsyswk_build_translation_matrix_full(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
-		complex double *target,
+		_Complex double *target,
 		const qpms_scatsys_at_omega_k_t *sswk
 		);
 
 
 /// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Periodic systems only.
 qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(
-		complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
+		_Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
 		int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
 		const qpms_scatsys_at_omega_k_t *sswk
 		);
@@ -539,7 +543,7 @@ struct beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
 		const qpms_scatsys_t *ss, 
 		/// Wavevector in cartesian coordinates (must lie in the lattice plane).
 		const double k[3],
-		complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
+		_Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
 		double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
 		double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
 		size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@@ -561,12 +565,12 @@ struct qpms_finite_group_t;
 
 /// Constructs a "full matrix action" of a point group element for an orbit type.
 /** TODO detailed doc */
-complex double *qpms_orbit_action_matrix(
+_Complex double *qpms_orbit_action_matrix(
 		/// Target array. If NULL, a new one is allocated.
 		/** The size of the array is (orbit->size * bspec->n)**2
 		 * (it makes sense to assume all the T-matrices share their spec).
 		 */
-		complex double *target,
+		_Complex double *target,
 		/// The orbit (type).
 		const qpms_ss_orbit_type_t *orbit,
 		/// Base spec of the t-matrices (we don't know it from orbit, as it has 
@@ -579,12 +583,12 @@ complex double *qpms_orbit_action_matrix(
 
 /// Constructs a dense matrix representation of a irrep projector for an orbit type.
 /** TODO detailed doc */
-complex double *qpms_orbit_irrep_projector_matrix(
+_Complex double *qpms_orbit_irrep_projector_matrix(
 		/// Target array. If NULL, a new one is allocated.
 		/** The size of the array is (orbit->size * bspec->n)**2
 		 * (it makes sense to assume all the T-matrices share their spec).
 		 */
-		complex double *target,
+		_Complex double *target,
 		/// The orbit (type).
 		const qpms_ss_orbit_type_t *orbit,
 		/// Base spec of the t-matrices (we don't know it from orbit, as it has 
@@ -596,14 +600,14 @@ complex double *qpms_orbit_irrep_projector_matrix(
 		const qpms_iri_t iri);
 
 /// TODO DOC!!!!!
-complex double *qpms_orbit_irrep_basis(
+_Complex double *qpms_orbit_irrep_basis(
 		/// Here theh size of theh basis shall be saved,
 		size_t *basis_size,
 		/// Target array. If NULL, a new one is allocated.
 		/** The size of the array is basis_size * (orbit->size * bspec->n)
 		 * (it makes sense to assume all the T-matrices share their spec).
 		 */
-		complex double *target,
+		_Complex double *target,
 		/// The orbit (type).
 		const qpms_ss_orbit_type_t *orbit,
 		/// Base spec of the t-matrices (we don't know it from orbit, as it has 
@@ -618,10 +622,10 @@ complex double *qpms_orbit_irrep_basis(
 /// Creates an incident field vector in the full basis, given a function that evaluates the field expansions at points.
 /** TODO detailed doc!
  * \returns target_full if target_full was not NULL, otherwise the newly allocated array. */
-complex double *qpms_scatsys_incident_field_vector_full(
+_Complex double *qpms_scatsys_incident_field_vector_full(
 		/// Target array. If NULL, a new one is allocated.
 		/** The length of the array is ss->fecv_size. */
-		complex double *target_full,
+		_Complex double *target_full,
 		const qpms_scatsys_t *ss,
 		qpms_incfield_t field_at_point,
 		const void *args, ///< Pointer passed as the last argument to (*field_at_point)()
@@ -629,9 +633,9 @@ complex double *qpms_scatsys_incident_field_vector_full(
 		);
 
 /// Applies T-matrices onto an incident field vector in the full basis.
-complex double *qpms_scatsysw_apply_Tmatrices_full(
-		complex double *target_full, /// Target vector array. If NULL, a new one is allocated.
-		const complex double *inc_full, /// Incident field coefficient vector. Must not be NULL.
+_Complex double *qpms_scatsysw_apply_Tmatrices_full(
+		_Complex double *target_full, /// Target vector array. If NULL, a new one is allocated.
+		const _Complex double *inc_full, /// Incident field coefficient vector. Must not be NULL.
 		const qpms_scatsys_at_omega_t *ssw
 		);
 
@@ -650,7 +654,7 @@ struct beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
 		const qpms_scatsys_t *ss, 
 		/// A valid irrep index to search only in one irrep, or QPMS_NO_IRREP for solving the full system.
 		qpms_iri_t iri,
-		complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
+		_Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
 		double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
 		double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
 		size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@@ -673,7 +677,7 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes(
 		const qpms_scatsys_t *ss,
 		double eta, ///< Ewald sum parameter
 		const double *beta_, ///< k-vector of corresponding dimensionality, NULL/ignored for finite system.
-		complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
+		_Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
 		double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
 		double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
 		size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@@ -687,10 +691,10 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes(
 #if 0
 /// Creates a (partial) incident field vector in the symmetry-adapted basis, given a function that evaluates the field expansions at points.
 /** TODO detailed doc! */
-complex double *qpms_scatsys_incident_field_vector_irrep_packed(
+_Complex double *qpms_scatsys_incident_field_vector_irrep_packed(
 		/// Target array. If NULL, a new one is allocated.
 		/** The length of the array is ss->fecv_size. */
-		complex double *target_full,
+		_Complex double *target_full,
 		const qpms_scatsys_t *ss,
 		const qpms_iri_t iri, ///< The index of given irreducible representation of ss->sym.
 		qpms_incfield_t field_at_point,
@@ -706,7 +710,7 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
  *
  * \return Complex electric field at the point defined by \a evalpoint.
  * 
- * \see qpms_scatsysw_scattered_E()
+ * \see qpms_scatsysw_scattered_E(), qpms_scatsys_scattered_field_basis()
  *
  * \see qpms_scatsyswk_scattered_E() for periodic systems.
  *
@@ -715,8 +719,8 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
 ccart3_t qpms_scatsys_scattered_E(
 		const qpms_scatsys_t *ss,
 		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
-		complex double wavenumber, ///< Wavenumber of the background medium.
-		const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
+		_Complex double wavenumber, ///< Wavenumber of the background medium.
+		const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
 		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
 		);
 
@@ -727,14 +731,96 @@ ccart3_t qpms_scatsys_scattered_E(
  *
  * \return Complex electric field at the point defined by \a evalpoint.
  * 
- * \see qpms_scatsys_scattered_E()
+ * \see qpms_scatsys_scattered_E(), qpms_scatsysw_scattered_field_basis()
  *
  * \see qpms_scatsyswk_scattered_E() for periodic systems.
  */
 ccart3_t qpms_scatsysw_scattered_E(
 		const qpms_scatsys_at_omega_t *ssw,
 		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
-		const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
+		const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
+		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
+		);
+
+/// Evaluates a "basis" for electric field at a given point.
+/**
+ * This function evaluates all the included VSWFs from the particles in the system, evaluated 
+ * at a given point. Taking a linear combination of these with the coefficients \a scattcoeff_full[]
+ * would be equivalent to the result of qpms_scatsys_scattered_E().
+ *
+ * \see qpms_scatsysw_scattered_field_basis()
+ *
+ * \see qpms_scatsyswk_scattered_field_basis() for periodic systems.
+ *
+ */
+qpms_errno_t qpms_scatsys_scattered_field_basis(
+		ccart3_t *target, ///< Target array of length \a ss->fecv_size
+		const qpms_scatsys_t *ss,
+		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
+		_Complex double wavenumber, ///< Wavenumber of the background medium.
+		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
+		);
+
+/// Evaluates a "basis" for electric field of a given particle at a given point.
+/**
+ * This function evaluates all the included VSWFs from a particle in the system, evaluated 
+ * at a given point. 
+ *
+ * Finite systems only.
+ *
+ * \see qpms_scatsys_scattered_field_basis() for the whole-array equivalent.
+ *
+ * \see qpms_scatsysw_scattered_field_basis_pi()
+ */
+qpms_errno_t qpms_scatsys_scattered_field_basis_pi(
+		ccart3_t *target, ///< Target array of length at least `qpms_ss_bspec_pi(ssw->ss, pi)->n`
+		const qpms_scatsys_t *ss,
+		qpms_ss_pi_t pi, ///< Particle index
+		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
+		_Complex double wavenumber, ///< Wavenumber of the background medium
+		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
+		);
+
+/// Evaluates a "basis" for electric field at a given point.
+/**
+ * This function evaluates all the included VSWFs from the particles in the system, evaluated 
+ * at a given point. Taking a linear combination of these with the coefficients \a scattcoeff_full[]
+ * would be equivalent to the result of qpms_scatsysw_scattered_E().
+ *
+ * Note that this might require a relatively large amount of memory. Depending
+ * on the application, qpms_scatsysw_scattered_field_basis_pi() might be a better
+ * alternative.
+ *
+ * \see qpms_scatsysw_scattered_field_basis_pi() for the single-particle equivalent.
+ *
+ * \see qpms_scatsys_scattered_field_basis()
+ *
+ * \see qpms_scatsyswk_scattered_field_basis() for periodic systems.
+ *
+ */
+qpms_errno_t qpms_scatsysw_scattered_field_basis(
+		ccart3_t *target, ///< Target array of length \a ss->fecv_size
+		const qpms_scatsys_at_omega_t *ssw,
+		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
+		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
+		);
+
+/// Evaluates a "basis" for electric field of a given particle at a given point.
+/**
+ * This function evaluates all the included VSWFs from a particle in the system, evaluated 
+ * at a given point. 
+ *
+ * Finite systems only.
+ *
+ * \see qpms_scatsysw_scattered_field_basis() for the whole-array equivalent.
+ *
+ * \see qpms_scatsys_scattered_field_basis_pi()
+ */
+qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(
+		ccart3_t *target, ///< Target array of length at least `qpms_ss_bspec_pi(ssw->ss, pi)->n`
+		const qpms_scatsys_at_omega_t *ssw,
+		qpms_ss_pi_t pi, ///< Particle index
+		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
 		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
 		);
 
@@ -750,8 +836,8 @@ ccart3_t qpms_scatsysw_scattered_E(
 ccart3_t qpms_scatsys_scattered_E__alt(
 		const qpms_scatsys_t *ss,
 		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
-		complex double wavenumber, ///< Wavenumber of the background medium.
-		const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
+		_Complex double wavenumber, ///< Wavenumber of the background medium.
+		const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
 		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
 		);
 
@@ -767,7 +853,7 @@ ccart3_t qpms_scatsys_scattered_E__alt(
 ccart3_t qpms_scatsysw_scattered_E__alt(
 		const qpms_scatsys_at_omega_t *ssw,
 		qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
-		const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
+		const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
 		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
 		);
 
@@ -786,10 +872,18 @@ ccart3_t qpms_scatsysw_scattered_E__alt(
 ccart3_t qpms_scatsyswk_scattered_E(
 		const qpms_scatsys_at_omega_k_t *sswk,
 		qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported).
-		const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
+		const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
 		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
 		);
 
+ccart3_t qpms_scatsyswk_scattered_E_e(
+		const qpms_scatsys_at_omega_k_t *sswk,
+		qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported).
+		const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
+		cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the field is evaluated.
+		qpms_ewald_part parts
+		);
+
 /// Evaluates "scattered" field basis functions in a periodic system.
 /**
  *
@@ -802,10 +896,18 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
 		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the basis is evaluated.
 		);
 
+qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
+		ccart3_t *target, ///< Target array of size sswk->ssw->ss->fecv_size
+		const qpms_scatsys_at_omega_k_t *sswk,
+		qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
+		cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
+		qpms_ewald_part parts
+		);
+
 
 /// Adjusted Ewadl parameter to avoid high-frequency breakdown.
 // TODO DOC
-double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double wavevector[3]);
+double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, _Complex double wavenumber, const double wavevector[3]);
 
 #if 0
 /** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector)
@@ -815,10 +917,13 @@ double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber,
  */
 ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
 		qpms_iri_t iri, ///< Irreducible representation
-		const complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri.
+		const _Complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri.
 		cart3_t where, ///< Evaluation point.
-		complex double k ///< Wave number.
+		_Complex double k ///< Wave number.
 		);
 #endif
 
+#ifdef __cplusplus
+}
+#endif
 #endif //QPMS_SCATSYSTEM_H

qpms/scatsystem_dbg.c

@@ -0,0 +1,48 @@
+#include "scatsystem_dbg.h"
+#include "translations_dbg.h"
+#include "indexing.h"
+
+
+// analogue of qpms_scatsyswk_scattered_field_basis_e()
+qpms_errno_t qpms_scatsyswk_test_sswf_basis_e(
+    complex double *target, ///< Target array of size (2 * sswk->ssw->ss->c->lMax + 1) * sswk->ssw->ss->p_count
+    const qpms_scatsys_at_omega_k_t *sswk,
+    qpms_bessel_t btyp, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
+    cart3_t where, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
+    qpms_ewald_part parts
+    ) {
+  QPMS_UNTESTED;
+  if (btyp != QPMS_HANKEL_PLUS)
+    QPMS_NOT_IMPLEMENTED("Only scattered field with first kind Hankel functions currently implemented.");
+  const qpms_scatsys_t *ss = sswk->ssw->ss;
+  if (ss->lattice_dimension != 2)
+    QPMS_NOT_IMPLEMENTED("Only 2D-periodic lattices implemented");
+  //ccart3_t res = {0,0,0};
+  //ccart3_t res_kc = {0,0,0}; // kahan sum compensation
+  
+  const size_t particle_nelem = qpms_lMax2nelem_sc(2*ss->c->lMax + 1);
+
+
+  for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
+    const cart3_t particle_pos = ss->p[pi].pos;
+
+    const cart3_t origin_cart = {0, 0, 0};
+
+    QPMS_ASSERT(sswk->k[2] == 0); // At least not implemented now
+    QPMS_ASSERT(ss->per.lattice_basis[0].z == 0);
+    QPMS_ASSERT(ss->per.lattice_basis[1].z == 0);
+
+    // Same choices as in qpms_ss_ppair_W32xy; TODO make it more dynamic
+    const double maxR = sqrt(ss->per.unitcell_volume) * 64;
+    const double maxK = 2048 * 2 * M_PI / maxR;
+
+		QPMS_ENSURE_SUCCESS(qpms_trans_calculator_test_sswf_e32(ss->c,
+          target + pi * particle_nelem, NULL,
+          sswk->eta, sswk->ssw->wavenumber,
+          cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
+          cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
+          maxR, maxK, parts));
+  }
+
+  return QPMS_SUCCESS;
+}

qpms/scatsystem_dbg.h

@@ -0,0 +1,10 @@
+#include "scatsystem.h"
+
+qpms_errno_t qpms_scatsyswk_test_sswf_basis_e(
+                _Complex double *target, ///< Target array of size sswk->ssw->ss->fecv_size
+                const qpms_scatsys_at_omega_k_t *sswk,
+                qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
+                cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
+		qpms_ewald_part parts
+                );
+

qpms/symmetries.h

@@ -21,40 +21,44 @@
  */
 #ifndef SYMMETRIES_H
 #define SYMMETRIES_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #include "qpms_types.h"
 #include <cblas.h>
 
 /// Dense matrix representation of the z coordinate sign flip operation (xy-plane mirroring).
-complex double *qpms_zflip_uvswi_dense(
-		complex double *target, ///< If NULL, a new array is allocated.
+_Complex double *qpms_zflip_uvswi_dense(
+		_Complex double *target, ///< If NULL, a new array is allocated.
 		const qpms_vswf_set_spec_t *bspec);
 /// Dense matrix representation of the y coordinate sign flip operation (xz-plane mirroring).
-complex double *qpms_yflip_uvswi_dense(
-		complex double *target, ///< If NULL, a new array is allocated.
+_Complex double *qpms_yflip_uvswi_dense(
+		_Complex double *target, ///< If NULL, a new array is allocated.
 		const qpms_vswf_set_spec_t *bspec);
 /// Dense matrix representation of the x coordinate sign flip operation (yz-plane mirroring).
-complex double *qpms_xflip_uvswi_dense(
-		complex double *target, ///< If NULL, a new array is allocated.
+_Complex double *qpms_xflip_uvswi_dense(
+		_Complex double *target, ///< If NULL, a new array is allocated.
 		const qpms_vswf_set_spec_t *bspec);
 /// Dense matrix representation of a rotation around the z-axis
-complex double *qpms_zrot_uvswi_dense(
-		complex double *target, ///< If NULL, a new array is allocated.
+_Complex double *qpms_zrot_uvswi_dense(
+		_Complex double *target, ///< If NULL, a new array is allocated.
 		const qpms_vswf_set_spec_t *bspec,
 		double phi ///< Rotation angle
 		);
 /// Dense matrix representation of a "rational" rotation around the z-axis
 /** Just for convenience. Corresponds to the angle \f$ \phi = 2\piw/N \f$.
  */
-complex double *qpms_zrot_rational_uvswi_dense(
-		complex double *target, ///< If NULL, a new array is allocated.
+_Complex double *qpms_zrot_rational_uvswi_dense(
+		_Complex double *target, ///< If NULL, a new array is allocated.
 		const qpms_vswf_set_spec_t *bspec,
 		int N,
 		int w
 		);
 
 /// Dense matrix (uvswi-indexed) representation of any O(3) transformation.
-complex double *qpms_irot3_uvswfi_dense(
-		complex double *target, ///< If NULL, a new array is allocated.
+_Complex double *qpms_irot3_uvswfi_dense(
+		_Complex double *target, ///< If NULL, a new array is allocated.
 		const qpms_vswf_set_spec_t *bspec,
 		const qpms_irot3_t transf);
 
@@ -74,5 +78,9 @@ size_t qpms_zero_roundoff_clean(double *arr, size_t nmemb, double atol);
  * Works on real and imaginary parts separately.
  * TODO doc.
  */
-size_t qpms_czero_roundoff_clean(complex double *arr, size_t nmemb, double atol);
+size_t qpms_czero_roundoff_clean(_Complex double *arr, size_t nmemb, double atol);
+
+#ifdef __cplusplus
+}
+#endif
 #endif // SYMMETRIES_H

qpms/tiny_inlines.h

@@ -4,6 +4,7 @@
 #ifndef TINY_INLINES_H
 #define TINY_INLINES_H
 #include <stdlib.h>
+#include <complex.h>
 
 static inline int min1pow(int pow) { return (pow % 2) ? -1 : 1; }
 
@@ -19,8 +20,8 @@ static inline int min1pow_m_neg(int m) {
 
 #if 0
 #ifdef __GSL_SF_LEGENDRE_H__
-static inline complex double
-spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m, double P_n_abs_m, complex double exp_imf) {
+static inline _Complex double
+spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m, double P_n_abs_m, _Complex double exp_imf) {
 
 	return;
 }
@@ -28,9 +29,9 @@ spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m,
 #endif
 
 // this has shitty precision:
-// static inline complex double ipow(int x) { return cpow(I, x); }
+// static inline _Complex double ipow(int x) { return cpow(I, x); }
 
-static inline complex double ipow(int x) {
+static inline _Complex double ipow(int x) {
   x = ((x % 4) + 4) % 4;
   switch(x) {
     case 0:

qpms/tmatrices.c

@@ -988,6 +988,12 @@ qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t,
   return QPMS_SUCCESS;
 }
 
+qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t,
+    complex double omega_ignored, const void *params_ignored) {
+  memset(t->m, 0, SQ(t->spec->n)*sizeof(*(t->m)));
+  return QPMS_SUCCESS;
+}
+
 void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) {
   switch(f->typ) {
     case QPMS_TMATRIX_OPERATION_NOOP:

qpms/tmatrices.h

@@ -3,6 +3,10 @@
  */
 #ifndef TMATRICES_H
 #define TMATRICES_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 // #include "qpms_types.h" // included via materials.h
 // #include <gsl/gsl_spline.h> // included via materials.h
 #include "materials.h"
@@ -14,7 +18,7 @@ struct qpms_finite_group_t;
 typedef struct qpms_finite_group_t qpms_finite_group_t;
 
 /// Returns a pointer to the beginning of the T-matrix row \a rowno.
-static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
+static inline _Complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
 	return t->m + (t->spec->n * rowno);
 }
 
@@ -42,7 +46,7 @@ void qpms_tmatrix_free(qpms_tmatrix_t *t);
  *  This function actually checks for identical vswf specs.
  * TODO define constants with "default" atol, rtol for this function.
  */
-bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
+_Bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
 		const double rtol, const double atol);
 
 /// Creates a T-matrix from another matrix and a symmetry operation.
@@ -53,7 +57,7 @@ bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
  */
 qpms_tmatrix_t *qpms_tmatrix_apply_symop(
 		const qpms_tmatrix_t *T, ///< the original T-matrix
-		const complex double *M ///< the symmetry op matrix in row-major format
+		const _Complex double *M ///< the symmetry op matrix in row-major format
 		);
 
 /// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data.
@@ -64,7 +68,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop(
  */
 qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
 		qpms_tmatrix_t *T, ///< the original T-matrix
-		const complex double *M ///< the symmetry op matrix in row-major format
+		const _Complex double *M ///< the symmetry op matrix in row-major format
 		);
 
 /// Symmetrizes a T-matrix with an involution symmetry operation.
@@ -75,7 +79,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
  */
 qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
 		const qpms_tmatrix_t *T, ///< the original T-matrix
-		const complex double *M ///< the symmetry op matrix in row-major format
+		const _Complex double *M ///< the symmetry op matrix in row-major format
 		);
 
 /// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix.
@@ -105,7 +109,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(
  */
 qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
 		qpms_tmatrix_t *T, ///< the original T-matrix
-		const complex double *M ///< the symmetry op matrix in row-major format
+		const _Complex double *M ///< the symmetry op matrix in row-major format
 		);
 
 /// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix, rewriting the original one.
@@ -149,7 +153,7 @@ qpms_errno_t qpms_load_scuff_tmatrix(
 		double **freqs_su, ///< Frequencies in SCUFF units (optional).
 		/// The resulting T-matrices (optional).
 		qpms_tmatrix_t **tmatrices_array, 
-		complex double **tmdata ///< The T-matrices raw contents
+		_Complex double **tmdata ///< The T-matrices raw contents
 		);
 
 /// Tells whether qpms_load_scuff_tmatrix should crash if fopen() fails.
@@ -161,7 +165,7 @@ qpms_errno_t qpms_load_scuff_tmatrix(
  * This is desirable e.g. when used in Python (so that proper exception can
  * be thrown).
  */
-extern bool qpms_load_scuff_tmatrix_crash_on_failure;
+extern _Bool qpms_load_scuff_tmatrix_crash_on_failure;
 
 /// Loads a scuff-tmatrix generated file.
 /** A simple wrapper over qpms_read_scuff_tmatrix() that needs a 
@@ -189,7 +193,7 @@ qpms_errno_t qpms_read_scuff_tmatrix(
 		 * is accessed via
 		 * (*tmdata)[bspec->n*bspec->n*fi + desti*bspec->n + srci].
 		 */
-		complex double ** tmdata 
+		_Complex double ** tmdata 
 		);
 
 /// In-place application of point group elements on raw T-matrix data.
@@ -200,7 +204,7 @@ qpms_errno_t qpms_read_scuff_tmatrix(
  * TODO more doc.
  */
 qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
-		complex double *tmdata, const size_t tmcount,
+		_Complex double *tmdata, const size_t tmcount,
 		const qpms_vswf_set_spec_t *bspec,
 		size_t n_symops,
 		const qpms_irot3_t *symops
@@ -213,7 +217,7 @@ qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
  * TODO more doc.
  */
 qpms_errno_t qpms_symmetrise_tmdata_finite_group(
-		complex double *tmdata, const size_t tmcount,
+		_Complex double *tmdata, const size_t tmcount,
 		const qpms_vswf_set_spec_t *bspec,
 		const qpms_finite_group_t *pointgroup
 		);
@@ -242,9 +246,9 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
 		);
 
 /// Application of T-matrix on a vector of incident field coefficients, \f$ f = Ta \f$.
-complex double *qpms_apply_tmatrix(
-		complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated.
-		const complex double *a, ///< Incident field coefficient array of size T->spec->n.
+_Complex double *qpms_apply_tmatrix(
+		_Complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated.
+		const _Complex double *a, ///< Incident field coefficient array of size T->spec->n.
 		const qpms_tmatrix_t *T ///< T-matrix \a T to apply.
 		);
 
@@ -254,10 +258,11 @@ complex double *qpms_apply_tmatrix(
  * qpms_tmatrix_generator_interpolator()
  * qpms_tmatrix_generator_sphere()
  * qpms_tmatrix_generator_constant()
+ * qpms_tmatrix_generator_zero()
  */
 typedef struct qpms_tmatrix_generator_t {
 	qpms_errno_t (*function) (qpms_tmatrix_t *t, ///< T-matrix to fill.
-			complex double omega, ///< Angular frequency.
+			_Complex double omega, ///< Angular frequency.
 			const void *params ///< Implementation dependent parameters.
 			);
 	const void *params; ///< Parameter pointer passed to the function.
@@ -267,7 +272,7 @@ typedef struct qpms_tmatrix_generator_t {
 qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
 		const qpms_vswf_set_spec_t *bspec,
 		qpms_tmatrix_generator_t gen,
-		complex double omega);
+		_Complex double omega);
 
 
 /// Implementation of qpms_matrix_generator_t that just copies a constant matrix.
@@ -275,11 +280,20 @@ qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
  * the same base spec.
  */
 qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t,
-		complex double omega,
+		_Complex double omega_ignored,
 		/// Source T-matrix, real type is (const qpms_tmatrix_t*).
 		const void *tmatrix_orig
 		);
 
+/// Dummy implementation of qpms_tmatrix_generator_t, yielding a zero matrix.
+/**
+ * This effectively represents a part of the background medium without
+ * any scatterer present.
+ */
+qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t,
+    _Complex double omega_ignored,
+    const void *params_ignored);
+
 /* Fuck this, include the whole <gsl/gsl_spline.h>
 typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct
 typedef struct gsl_interp_type gsl_interp_type;
@@ -295,7 +309,7 @@ extern const gsl_interp_type * gsl_interp_steffen;
 // struct gsl_interp_accel; // use if lookup proves to be too slow
 typedef struct qpms_tmatrix_interpolator_t {
 	const qpms_vswf_set_spec_t *bspec;
-	//bool owns_bspec;
+	//_Bool owns_bspec;
 	gsl_spline **splines_real; ///< There will be a spline object for each nonzero element
 	gsl_spline **splines_imag; ///< There will be a spline object for each nonzero element
 	// gsl_interp_accel **accel_real;
@@ -317,7 +331,7 @@ qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t
 qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Number of freqs and T-matrices provided.
 		const double *freqs, const qpms_tmatrix_t *tmatrices_array, ///< N.B. array of qpms_tmatrix_t, not pointers!
 		const gsl_interp_type *iptype
-	       //, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
+	       //, _Bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
        	       /*, ...? */);
 
 
@@ -325,7 +339,7 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num
 /** As in qpms_tmatrix_interpolator_eval(), the imaginary part of frequency is discarded!
  */
 qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, ///< T-matrix to fill.
-			complex double omega, ///< Angular frequency.
+			_Complex double omega, ///< Angular frequency.
 			const void *interpolator ///< Parameter of type qpms_tmatrix_interpolator_t *.
 );
 
@@ -342,14 +356,14 @@ qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, ///< T-matri
  *
  * TODO better doc.
  */
-complex double *qpms_mie_coefficients_reflection(
-                complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated.
+_Complex double *qpms_mie_coefficients_reflection(
+                _Complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated.
                 const qpms_vswf_set_spec_t *bspec, ///< Defines which of the coefficients are calculated.
                 double a, ///< Radius of the sphere.
-                complex double k_i, ///< Wave number of the internal material of the sphere.
-                complex double k_e, ///< Wave number of the surrounding medium.
-                complex double mu_i, ///< Relative permeability of the sphere material.
-                complex double mu_e, ///< Relative permeability of the surrounding medium.
+                _Complex double k_i, ///< Wave number of the internal material of the sphere.
+                _Complex double k_e, ///< Wave number of the surrounding medium.
+                _Complex double mu_i, ///< Relative permeability of the sphere material.
+                _Complex double mu_e, ///< Relative permeability of the surrounding medium.
                 qpms_bessel_t J_ext, ///< Kind of the "incoming" waves. Most likely QPMS_BESSEL_REGULAR.
                 qpms_bessel_t J_scat ///< Kind of the "scattered" waves. Most likely QPMS_HANKEL_PLUS.
                 );
@@ -357,10 +371,10 @@ complex double *qpms_mie_coefficients_reflection(
 /// Replaces the contents of an existing T-matrix with that of a spherical nanoparticle calculated using the Lorentz-mie theory.
 qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
 		double a, ///< Radius of the sphere.
-		complex double k_i, ///< Wave number of the internal material of the sphere.
-		complex double k_e, ///< Wave number of the surrounding medium.
-		complex double mu_i, ///< Relative permeability of the sphere material.
-		complex double mu_e ///< Relative permeability of the surrounding medium.
+		_Complex double k_i, ///< Wave number of the internal material of the sphere.
+		_Complex double k_e, ///< Wave number of the surrounding medium.
+		_Complex double mu_i, ///< Relative permeability of the sphere material.
+		_Complex double mu_e ///< Relative permeability of the surrounding medium.
 		);
 
 /// Parameter structure for qpms_tmatrix_generator_sphere().
@@ -372,7 +386,7 @@ typedef struct qpms_tmatrix_generator_sphere_param_t {
 
 /// T-matrix generator for spherical particles using Lorentz-Mie solution.
 qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t,
-		complex double omega,
+		_Complex double omega,
 		const void *params ///< Of type qpms_tmatrix_generator_sphere_param_t.
 		);
 
@@ -380,10 +394,10 @@ qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t,
 static inline qpms_tmatrix_t *qpms_tmatrix_spherical(
 		const qpms_vswf_set_spec_t *bspec,
 		double a, ///< Radius of the sphere.
-		complex double k_i, ///< Wave number of the internal material of the sphere.
-		complex double k_e, ///< Wave number of the surrounding medium.
-		complex double mu_i, ///< Relative permeability of the sphere material.
-		complex double mu_e ///< Relative permeability of the surrounding medium.
+		_Complex double k_i, ///< Wave number of the internal material of the sphere.
+		_Complex double k_e, ///< Wave number of the surrounding medium.
+		_Complex double mu_i, ///< Relative permeability of the sphere material.
+		_Complex double mu_e ///< Relative permeability of the surrounding medium.
 		) {
 	qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
 	qpms_tmatrix_spherical_fill(t, a, k_i, k_e, mu_i, mu_e);
@@ -395,8 +409,8 @@ qpms_errno_t qpms_tmatrix_spherical_mu0_fill(
 		qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
 		double a, ///< Radius of the sphere.
 		double omega, ///< Angular frequency.
-		complex double epsilon_fg, ///< Relative permittivity of the sphere.
-		complex double epsilon_bg ///< Relative permittivity of the background medium.
+		_Complex double epsilon_fg, ///< Relative permittivity of the sphere.
+		_Complex double epsilon_bg ///< Relative permittivity of the background medium.
 		); 
 
 /// Convenience function to calculate T-matrix of a non-magnetic spherical particle using the permittivity values.
@@ -404,8 +418,8 @@ static inline qpms_tmatrix_t *qpms_tmatrix_spherical_mu0(
 		const qpms_vswf_set_spec_t *bspec,
 		double a, ///< Radius of the sphere.
 		double omega, ///< Angular frequency.
-		complex double epsilon_fg, ///< Relative permittivity of the sphere.
-		complex double epsilon_bg ///< Relative permittivity of the background medium.
+		_Complex double epsilon_fg, ///< Relative permittivity of the sphere.
+		_Complex double epsilon_bg ///< Relative permittivity of the background medium.
 		) {
 	qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
 	qpms_tmatrix_spherical_mu0_fill(t, a, omega, epsilon_fg, epsilon_bg);
@@ -460,7 +474,7 @@ qpms_arc_function_retval_t qpms_arc_sphere(double theta,
  */
 qpms_errno_t qpms_tmatrix_axialsym_fill(
 		qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
-		complex double omega, ///< Angular frequency.
+		_Complex double omega, ///< Angular frequency.
 		qpms_epsmu_t outside, ///< Optical properties of the outside medium.
 		qpms_epsmu_t inside, ///< Optical properties of the particle's material.
 		qpms_arc_function_t shape, ///< Particle surface parametrisation.
@@ -474,7 +488,7 @@ qpms_errno_t qpms_tmatrix_axialsym_fill(
 /// Creates a new T-matrix of a particle with \f$ C_\infty \f$ symmetry.
 static inline qpms_tmatrix_t *qpms_tmatrix_axialsym(
 		const qpms_vswf_set_spec_t *bspec,
-		complex double omega, ///< Angular frequency.
+		_Complex double omega, ///< Angular frequency.
 		qpms_epsmu_t outside, ///< Optical properties of the outside medium.
 		qpms_epsmu_t inside, ///< Optical properties of the particle's material.
 		qpms_arc_function_t shape, ///< Particle surface parametrisation.
@@ -501,22 +515,22 @@ typedef struct qpms_tmatrix_generator_axialsym_param_t {
 
 /// qpms_tmatrix_axialsym for qpms_tmatrix_generator_t
 qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, ///< T-matrix to fill.
-			complex double omega, ///< Angular frequency.
+			_Complex double omega, ///< Angular frequency.
 			const void *params ///< Parameters of type qpms_tmatrix_generator_axialsym_param_t.
 );
 
 
 /// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful for debugging).
-qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(complex double *target,
-		complex double omega,
+qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(_Complex double *target,
+		_Complex double omega,
 		const qpms_tmatrix_generator_axialsym_param_t *param,
 		qpms_normalisation_t norm,
 		qpms_bessel_t J
 		);
 
 /// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful mostly for debugging).
-qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(complex double *target,
-		complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside,
+qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(_Complex double *target,
+		_Complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside,
 		qpms_arc_function_t shape, qpms_l_t lMaxQR, qpms_normalisation_t norm,
 		qpms_bessel_t J ///< Use QPMS_BESSEL_REGULAR to calculate \f$ R^T\f$ or QPMS_HANKEL_PLUS to calculate \f$ Q^T\f$.
 		);
@@ -537,7 +551,7 @@ typedef struct qpms_tmatrix_function_t {
 
 /// Convenience function to create a new T-matrix from qpms_tmatrix_function_t.
 // FIXME the name is not very intuitive.
-static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, complex double omega) {
+static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, _Complex double omega) {
 	return qpms_tmatrix_init_from_generator(f.spec, *f.gen, omega);
 }
 
@@ -558,9 +572,9 @@ typedef enum {
 struct qpms_tmatrix_operation_lrmatrix {
 	/// Raw matrix data of \a M in row-major order.
 	/** The matrix must be taylored for the given bspec! */
-	complex double *m; 
-        size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double).
-	bool owns_m; ///< Whether \a m is owned by this;
+	_Complex double *m; 
+        size_t m_size; ///< Total size of \a m matrix in terms of sizeof(_Complex double).
+	_Bool owns_m; ///< Whether \a m is owned by this;
 };
 
 struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations.
@@ -591,9 +605,9 @@ struct qpms_tmatrix_operation_compose_chain {
 struct qpms_tmatrix_operation_scmulz {
 	/// Raw matrix data of \a M in row-major order.
 	/** The matrix must be taylored for the given bspec! */
-	complex double *m;
-        size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double).
-	bool owns_m; ///< Whether \a m is owned by this.
+	_Complex double *m;
+        size_t m_size; ///< Total size of \a m matrix in terms of sizeof(_Complex double).
+	_Bool owns_m; ///< Whether \a m is owned by this.
 };
 
 /// Specifies a symmetrisation using a set of rotoreflections (with equal weights) for qpms_tmatrix_operation_t.
@@ -602,7 +616,7 @@ struct qpms_tmatrix_operation_scmulz {
 struct qpms_tmatrix_operation_irot3arr {
 	size_t n; ///< Number of rotoreflections;
 	qpms_irot3_t *ops; ///< Rotoreflection array of size \a n.
-	bool owns_ops; ///< Whether \a ops array is owned by this.
+	_Bool owns_ops; ///< Whether \a ops array is owned by this.
 };
 
 /// A generic T-matrix transformation operator.
@@ -689,4 +703,7 @@ typedef struct qpms_abstract_particle_t {
 typedef qpms_particle_tid_t qpms_abstract_particle_tid_t;
 #endif // 0
 
+#ifdef __cplusplus
+}
+#endif
 #endif //TMATRICES_H

qpms/tolerances.h

@@ -1,6 +1,9 @@
 /*! \file tolerances.h */
 #ifndef QPMS_TOLERANCES_H
 #define QPMS_TOLERANCES_H
+#ifdef __cplusplus
+extern "C" {
+#endif
 
 // TODO DOC
 
@@ -12,4 +15,7 @@ typedef struct qpms_tolerance_spec_t {
 /// A rather arbitrary default tolerance.
 static const qpms_tolerance_spec_t QPMS_TOLERANCE_DEFAULT = {.atol = 1e-9, .rtol =  1e-8};
 
+#ifdef __cplusplus
+}
+#endif
 #endif // QPMS_TOLERANCES_H

qpms/translations.c

@@ -1,5 +1,6 @@
 #include <math.h>
 #include "qpms_types.h"
+#include <complex.h>
 #include "qpms_specfunc.h"
 #include "gaunt.h"
 #include "translations.h"
@@ -15,6 +16,10 @@
 #include "normalisation.h"
 #include "translations_inlines.h"
 
+#if defined LATTICESUMS31 || defined LATTICESUMS32
+#include "lattices.h"
+#endif
+
 /*
  * Define macros with additional factors that "should not be there" according
  * to the "original" formulae but are needed to work with my vswfs.
@@ -825,7 +830,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
     )
 {
 
-  const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
+  const qpms_y_sc_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
   //const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
   const bool doerr = Aerr || Berr;
   const bool do_sigma0 = (particle_shift == 0)//DIFF21((particle_shift.x == 0) && (particle_shift.y == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
@@ -846,15 +851,15 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
   QPMS_ENSURE_SUCCESS(ewald31z_sigma_short_points_and_shift(sigmas_short, serr_short, //DIFF21
       c->e3c, eta, k, nRpoints, Rpoints, beta, particle_shift));
 
-  for(qpms_y_t y = 0; y < nelem2_sc; ++y)
+  for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
     sigmas_total[y] = sigmas_short[y] + sigmas_long[y];
-  if (doerr) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
+  if (doerr) for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
     serr_total[y]  = serr_short[y] + serr_long[y];
 
   complex double sigma0 = 0; double sigma0_err = 0;
   if (do_sigma0) {
     QPMS_ENSURE_SUCCESS(ewald31z_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
-    const qpms_l_t y = qpms_mn2y_sc(0,0);
+    const qpms_y_sc_t y = qpms_mn2y_sc(0,0);
     sigmas_total[y] += sigma0;
     if(doerr) serr_total[y] += sigma0_err;
   }
@@ -872,7 +877,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
             // TODO skip if ... (N.B. skip will be different for 31z and 32)
             for(qpms_l_t q = 0; q <= qmax; ++q) {
               const qpms_l_t p = n + nu - 2*q;
-              const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p); 
+              const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p); 
               const complex double multiplier = c->A_multipliers[i][q];
               complex double sigma = sigmas_total[y_sc];
               ckahanadd(&Asum, &Asumc, multiplier * sigma);
@@ -887,7 +892,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
             double Bsumerr, Bsumerrc; if(Berr) kahaninit(&Bsumerr, &Bsumerrc);
              for(qpms_l_t q = 0; q <= qmax; ++q) {
               const qpms_l_t p_ = n + nu - 2*q + 1;
-              const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p_); 
+              const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p_); 
               const complex double multiplier = c->B_multipliers[i][q-BQ_OFFSET];
               complex double sigma = sigmas_total[y_sc];
               ckahanadd(&Bsum, &Bsumc, multiplier * sigma);
@@ -937,7 +942,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
     )
 {
 
-  const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
+  const qpms_y_sc_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
   //const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
   const bool doerr = Aerr || Berr;
   const bool do_sigma0 = ((particle_shift.x == 0) && (particle_shift.y == 0) && (particle_shift.z == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
@@ -1006,17 +1011,17 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
       PGen_destroy(&Rgen);
   }
 
-  for(qpms_y_t y = 0; y < nelem2_sc; ++y)
+  for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
     sigmas_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? sigmas_short[y] : 0)
                     + ((parts & QPMS_EWALD_LONG_RANGE)  ? sigmas_long[y] : 0);
-  if (doerr) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
+  if (doerr) for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
     serr_total[y]  = ((parts & QPMS_EWALD_SHORT_RANGE) ? serr_short[y] : 0)
                     + ((parts & QPMS_EWALD_LONG_RANGE)  ? serr_long[y] : 0);
 
   complex double sigma0 = 0; double sigma0_err = 0;
   if (do_sigma0 && (parts & QPMS_EWALD_0TERM)) {
     QPMS_ENSURE_SUCCESS(ewald3_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
-    const qpms_l_t y = qpms_mn2y_sc(0,0);
+    const qpms_y_sc_t y = qpms_mn2y_sc(0,0);
     sigmas_total[y] += sigma0;
     if(doerr) serr_total[y] += sigma0_err;
   }
@@ -1034,7 +1039,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
             // TODO skip if ...
             for(qpms_l_t q = 0; q <= qmax; ++q) {
               const qpms_l_t p = n + nu - 2*q;
-              const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p); 
+              const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p); 
               const complex double multiplier = c->A_multipliers[i][q];
               complex double sigma = sigmas_total[y_sc];
               ckahanadd(&Asum, &Asumc, multiplier * sigma);
@@ -1049,7 +1054,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
             double Bsumerr, Bsumerrc; if(Berr) kahaninit(&Bsumerr, &Bsumerrc);
              for(qpms_l_t q = 0; q <= qmax; ++q) {
               const qpms_l_t p_ = n + nu - 2*q + 1;
-              const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p_); 
+              const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p_); 
               const complex double multiplier = c->B_multipliers[i][q-BQ_OFFSET];
               complex double sigma = sigmas_total[y_sc];
               ckahanadd(&Bsum, &Bsumc, multiplier * sigma);

qpms/translations.h

@@ -26,9 +26,12 @@
  */
 #ifndef QPMS_TRANSLATIONS_H
 #define QPMS_TRANSLATIONS_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #include "vectors.h"
 #include "qpms_types.h"
-#include <complex.h>
 #include <stdbool.h>
 #include <stddef.h>
 
@@ -37,10 +40,10 @@
 #endif
 
 
-complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
+_Complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
 		bool r_ge_d, qpms_bessel_t J);
 
-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, csph_t kdlj,
+_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, csph_t kdlj,
 		bool r_ge_d, qpms_bessel_t J);
 
 /// Structure holding the constant factors in normalisation operators.
@@ -65,8 +68,8 @@ typedef struct qpms_trans_calculator {
 	qpms_normalisation_t normalisation;
 	qpms_l_t lMax;
 	qpms_y_t nelem;
-	complex double **A_multipliers;
-	complex double **B_multipliers;
+	_Complex double **A_multipliers;
+	_Complex double **B_multipliers;
 #if 0
 	// Normalised values of the Legendre functions and derivatives
 	// for θ == π/2, i.e. for the 2D case.
@@ -90,44 +93,44 @@ qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, ///< Truncation
 /// 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,
+_Complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c,
 		qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
 		bool r_ge_d, qpms_bessel_t J);
-complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c,
+_Complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c,
 		qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
 		bool r_ge_d, qpms_bessel_t J);
 int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c,
-		complex double *Adest, complex double *Bdest,
+		_Complex double *Adest, _Complex double *Bdest,
 		qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
 		bool r_ge_d, qpms_bessel_t J);
 int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c,
-		complex double *Adest, complex double *Bdest,
+		_Complex double *Adest, _Complex double *Bdest,
 		size_t deststride, size_t srcstride,
 		csph_t kdlj, bool r_ge_d, qpms_bessel_t J); 
 
 // TODO update the types later
-complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
-		int m, int n, int mu, int nu, complex double kdlj_r,
+_Complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
+		int m, int n, int mu, int nu, _Complex double kdlj_r,
 		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 
-complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c,
-		int m, int n, int mu, int nu, complex double kdlj_r,
+_Complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c,
+		int m, int n, int mu, int nu, _Complex double kdlj_r,
 		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 
 int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c,
-		complex double *Adest, complex double *Bdest,
-		int m, int n, int mu, int nu, complex double kdlj_r,
+		_Complex double *Adest, _Complex double *Bdest,
+		int m, int n, int mu, int nu, _Complex double kdlj_r,
 		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 
 int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c,
-		complex double *Adest, complex double *Bdest,
+		_Complex double *Adest, _Complex double *Bdest,
 		size_t deststride, size_t srcstride,
-		complex double kdlj_r, double kdlj_theta, double kdlj_phi,
+		_Complex double kdlj_r, double kdlj_theta, double kdlj_phi,
 		int r_ge_d, int J);
 
 // Convenience functions using VSWF base specs
 qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *c,
-		complex double *target, 
+		_Complex double *target, 
 		/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. 
 		const qpms_vswf_set_spec_t *destspec, size_t deststride,
 		/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. 
@@ -138,12 +141,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
 /// and with automatic \a r_ge_d = `false`.
 qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
 		const qpms_trans_calculator *c,
-		complex double *target, 
+		_Complex double *target, 
 		/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. 
 		const qpms_vswf_set_spec_t *destspec, size_t deststride,
 		/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. 
 		const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
-		complex double k, cart3_t destpos, cart3_t srcpos,
+		_Complex double k, cart3_t destpos, cart3_t srcpos,
 		qpms_bessel_t J
 		/// Workspace has to be at least 2 * c->neleme**2 long
 		);
@@ -156,10 +159,10 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
 
 
 int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
-		complex double *Adest, double *Aerr,
-		complex double *Bdest, double *Berr,
+		_Complex double *Adest, double *Aerr,
+		_Complex double *Bdest, double *Berr,
 		const ptrdiff_t deststride, const ptrdiff_t srcstride,
-		const double eta, const complex double wavenumber,
+		const double eta, const _Complex double wavenumber,
 		cart2_t b1, cart2_t b2,
 		const cart2_t beta,
 		const cart3_t particle_shift,
@@ -167,10 +170,10 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
 		);
 
 int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
-		complex double *Adest, double *Aerr,
-		complex double *Bdest, double *Berr,
+		_Complex double *Adest, double *Aerr,
+		_Complex double *Bdest, double *Berr,
 		const ptrdiff_t deststride, const ptrdiff_t srcstride,
-		const double eta, const complex double wavenumber,
+		const double eta, const _Complex double wavenumber,
 		cart2_t b1, cart2_t b2,
 		const cart2_t beta,
 		const cart3_t particle_shift,
@@ -180,12 +183,12 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
 
 // Convenience functions using VSWF base specs
 qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c,
-		complex double *target, double *err,
+		_Complex double *target, double *err,
 		/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. 
 		const qpms_vswf_set_spec_t *destspec, size_t deststride,
 		/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. 
 		const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
-		const double eta, const complex double wavenumber,
+		const double eta, const _Complex double wavenumber,
 		cart2_t b1, cart2_t b2,
 		const cart2_t beta,
 		const cart3_t particle_shift,
@@ -193,12 +196,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat
 		);
 
 qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calculator *c,
-		complex double *target, double *err,
+		_Complex double *target, double *err,
 		/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. 
 		const qpms_vswf_set_spec_t *destspec, size_t deststride,
 		/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. 
 		const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
-		const double eta, const complex double wavenumber,
+		const double eta, const _Complex double wavenumber,
 		cart2_t b1, cart2_t b2,
 		const cart2_t beta,
 		const cart3_t particle_shift,
@@ -212,8 +215,8 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calcul
 // e31z means that the particles are positioned along the z-axis;
 // their positions and K-values are then denoted by a single z-coordinate
 int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_trans_calculator *c,
-		complex double *Adest, double *Aerr,
-		complex double *Bdest, double *Berr,
+		_Complex double *Adest, double *Aerr,
+		_Complex double *Bdest, double *Berr,
 		const ptrdiff_t deststride, const ptrdiff_t srcstride,
 		const double eta, const double k,
 		const double unitcell_area, // just lattice period
@@ -244,4 +247,7 @@ int qpms_cython_trans_calculator_get_AB_arrays_loop(
 #endif //QPMS_COMPILE_PYTHON_EXTENSIONS
 
 
+#ifdef __cplusplus
+}
+#endif
 #endif // QPMS_TRANSLATIONS_H

qpms/translations_dbg.c

@@ -0,0 +1,149 @@
+#include "qpms_types.h"
+#include "qpms_specfunc.h"
+#include "translations_dbg.h"
+#include "indexing.h"
+#include "tiny_inlines.h"
+#include "qpms_error.h"
+#include "normalisation.h"
+#include "lattices.h"
+
+
+int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
+    complex double *dest, const csph_t kdlj, const  qpms_bessel_t J) {
+  if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) {
+    for (qpms_l_t l = 0; l <= 2*c->lMax+1; ++l)
+      for (qpms_m_t m = 0; m <= l; ++m) {
+          const qpms_y_t y_sc = qpms_mn2y_sc(m, l);
+          dest[y_sc] = NAN+I*NAN;
+      }
+    // TODO warn? different return value?
+    return 0;
+  }
+  double legendre_buf[gsl_sf_legendre_array_n(c->lMax+c->lMax+1)];
+  complex double bessel_buf[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth.
+  double costheta = cos(kdlj.theta);
+  QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1,
+    costheta,-1,legendre_buf));
+  QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf));
+  for (qpms_l_t l = 0; l <= 2*c->lMax+1; ++l) 
+      for (qpms_m_t m = -l; m <= l; ++m) {
+          const qpms_y_t y_sc = qpms_mn2y_sc(m, l);
+          dest[y_sc] = bessel_buf[l]
+                    * legendre_buf[gsl_sf_legendre_array_index(l, abs(m))]
+                    * cexp(I * m * kdlj.phi);
+      }
+  return 0;
+}
+
+int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c,
+    complex double *dest, const complex double k, const cart3_t r, const qpms_bessel_t J) {
+  csph_t kdlj = cart2csph(r);
+  kdlj.r *= k;
+  return qpms_trans_calculator_test_sswf(c, dest, kdlj, J);
+}
+
+#ifdef LATTICESUMS32
+int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
+    complex double * const sigmas_total, double * serr_total,
+    const double eta, const complex double k,
+    const cart2_t b1, const cart2_t b2,
+    const cart2_t beta,
+    const cart3_t particle_shift,
+    double maxR, double maxK,
+    const qpms_ewald_part parts) 
+{
+  const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
+  //const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
+  const bool doerr = !!serr_total;
+  const bool do_sigma0 = ((particle_shift.x == 0) && (particle_shift.y == 0) && (particle_shift.z == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
+
+  complex double *sigmas_short = malloc(sizeof(complex double)*nelem2_sc);
+  complex double *sigmas_long = malloc(sizeof(complex double)*nelem2_sc);
+  //complex double *sigmas_total = malloc(sizeof(complex double)*nelem2_sc);
+  double *serr_short, *serr_long;//, *serr_total;
+  if(doerr) {
+    serr_short = malloc(sizeof(double)*nelem2_sc);
+    serr_long = malloc(sizeof(double)*nelem2_sc);
+    //serr_total = malloc(sizeof(double)*nelem2_sc);
+  } else serr_short = serr_long = serr_total = NULL;
+
+  const double unitcell_area = l2d_unitcell_area(b1, b2);
+  cart2_t rb1, rb2; // reciprocal basis
+  QPMS_ENSURE_SUCCESS(l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2));
+
+  if (parts & QPMS_EWALD_LONG_RANGE) {
+    PGen Kgen = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL,
+#ifdef GEN_KSHIFTEDPOINTS
+        beta,
+#else
+        CART2_ZERO,
+#endif
+        0, true, maxK, false);
+
+    QPMS_ENSURE_SUCCESS(ewald3_sigma_long(sigmas_long, serr_long, c->e3c, eta, k, 
+          unitcell_area, LAT_2D_IN_3D_XYONLY, &Kgen,
+#ifdef GEN_KSHIFTEDPOINTS
+          true,
+#else
+          false,
+#endif
+          cart22cart3xy(beta), particle_shift));
+    if(Kgen.stateData) // PGen not consumed entirely (converged earlier)
+      PGen_destroy(&Kgen);
+  }
+
+  if (parts & QPMS_EWALD_SHORT_RANGE) {
+    PGen Rgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL, 
+#ifdef GEN_RSHIFTEDPOINTS
+        cart2_scale(-1 /*CHECKSIGN*/, cart3xy2cart2(particle_shift)),
+#else
+        CART2_ZERO,
+#endif
+        0, !do_sigma0, maxR, false);
+#ifdef GEN_RSHIFTEDPOINTS // rather ugly hacks, LPTODO cleanup
+    if (particle_shift.z != 0) {
+      const cart3_t zshift = {0, 0, -particle_shift.z /*CHECKSIGN*/}; 
+      Rgen = Pgen_shifted_new(Rgen, zshift);
+    }
+#endif
+
+
+    QPMS_ENSURE_SUCCESS(ewald3_sigma_short(sigmas_short, serr_short, c->e3c, eta, k,
+          particle_shift.z ? LAT_2D_IN_3D : LAT_2D_IN_3D_XYONLY, &Rgen, 
+#ifdef GEN_RSHIFTEDPOINTS
+          true,
+#else
+          false,
+#endif
+          cart22cart3xy(beta), particle_shift));
+
+    if(Rgen.stateData) // PGen not consumed entirely (converged earlier)
+      PGen_destroy(&Rgen);
+  }
+
+  for(qpms_y_t y = 0; y < nelem2_sc; ++y)
+    sigmas_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? sigmas_short[y] : 0)
+                    + ((parts & QPMS_EWALD_LONG_RANGE)  ? sigmas_long[y] : 0);
+  if (doerr) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
+    serr_total[y]  = ((parts & QPMS_EWALD_SHORT_RANGE) ? serr_short[y] : 0)
+                    + ((parts & QPMS_EWALD_LONG_RANGE)  ? serr_long[y] : 0);
+
+  complex double sigma0 = 0; double sigma0_err = 0;
+  if (do_sigma0 && (parts & QPMS_EWALD_0TERM)) {
+    QPMS_ENSURE_SUCCESS(ewald3_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
+    const qpms_l_t y = qpms_mn2y_sc(0,0);
+    sigmas_total[y] += sigma0;
+    if(doerr) serr_total[y] += sigma0_err;
+  }
+
+  free(sigmas_short);
+  free(sigmas_long);
+  //free(sigmas_total);
+  if(doerr) {
+    free(serr_short);
+    free(serr_long);
+    //free(serr_total);
+  }
+  return 0;
+}
+#endif // LATTICESUMS32

qpms/translations_dbg.h

@@ -0,0 +1,52 @@
+/*! \file translations_dbg.h
+ * \brief Auxiliary functions for debugging and testing (complementary to translations.h)
+ * 
+ * The following functions evaluate the SSWFs the same way as
+ * certain functions in translations.h, but they leave them
+ * as they are, without rearranging them into the translation
+ * operators.
+ */
+#ifndef TRANSLATIONS_DBG_H
+
+#include "translations.h"
+
+/// Evaluates scalar spherical wavefunctions at a point (spherical coordinates)
+/**
+ * \f[
+ *   \phi_l^m (kr, \theta, \varphi) = h_l(kr) \rawLeg{l}{m}(\cos \theta) e^{im\varphi} 
+ * \f]
+ *
+ * The wavefunctions are evaluated up to the degree `2 * c->lMax + 1`.
+ *
+ * Currently, THIS ALWAYS USES NON-NORMALISED LEGENDRE FUNCTIONS (TODO FIXME?).
+ */
+int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
+    _Complex double *dest, ///< Result array, length `qpms_lMax2nelem_sc(c->lMax)`
+    const csph_t kdlj, ///< Evaluation point (incl. wavenumber)
+    const  qpms_bessel_t J);
+
+/// Evaluates scalar spherical wavefunctions at a point (cartesian coordinates)
+/**
+ * Currently, THIS ALWAYS USES NON-NORMALISED LEGENDRE FUNCTIONS (TODO FIXME?).
+ *
+ * \see qpms_trans_calculator_test_sswf()
+ */
+int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c,
+    _Complex double *dest, ///< Result array, length `qpms_lMax2nelem_sc(c->lMax)`
+    const _Complex double k, ///< Wavenumber
+    const cart3_t r, ///< Evaluation point (excl. wavenumber)
+    const qpms_bessel_t J);
+
+
+#ifdef LATTICESUMS32
+int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
+    _Complex double * const sigmas_total, double * serr_total,
+    const double eta, const _Complex double k,
+    const cart2_t b1, const cart2_t b2,
+    const cart2_t beta,
+    const cart3_t particle_shift,
+    double maxR, double maxK,
+    const qpms_ewald_part parts);
+#endif // LATTICESUMS32
+
+#endif // TRANSLATIONS_DBG_H

qpms/translations_inlines.h

@@ -6,12 +6,12 @@
 /// Rearranges the default-ordered "A,B" array elements into "bspec"-defined matrix.
 // TODO DOC
 static inline void qpms_trans_array_from_AB(
-		complex double *t,
+		_Complex double *t,
 		const qpms_vswf_set_spec_t *const t_destspec,
 		const size_t t_deststride,
 		const qpms_vswf_set_spec_t *const t_srcspec,
 		const size_t t_srcstride,
-		const complex double *const A, const complex double *const B,
+		const _Complex double *const A, const _Complex double *const B,
 		/// A and B matrices' lMax.
 		/** This also determines their size and stride: they are assumed to
 		 * be square matrices of size `nelem * nelem` where

qpms/trivialgroup.c

@@ -1,4 +1,5 @@
 #include "groups.h"
+#include <complex.h>
 
 /// Trivial group, with one (reduntant) generator.
 /**

qpms/vectors.h

@@ -4,6 +4,7 @@
 #ifndef VECTORS_H
 #define VECTORS_H
 #include <math.h>
+#include <complex.h>
 #ifndef M_PI_2
 #define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487)
 #endif
@@ -198,12 +199,12 @@ static inline double cart3_dist(const cart3_t a, const cart3_t b) {
 	return cart3norm(cart3_substract(a,b));
 }
 
-static inline bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) {
+static inline _Bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) {
 	return cart3_dist(a,b) <= atol + rtol * (cart3norm(b) + cart3norm(a)) * .5;
 }
 
 /// Complex 3D vector scaling.
-static inline ccart3_t ccart3_scale(const complex double c, const ccart3_t v) {
+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;
 }
@@ -221,7 +222,7 @@ static inline ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) {
 }
 
 /// Complex 3D cartesian vector "dot product" without conjugation.
-static inline complex double ccart3_dotnc(const ccart3_t a, const ccart3_t b) {
+static inline _Complex double ccart3_dotnc(const ccart3_t a, const ccart3_t b) {
 	return a.x * b.x + a.y * b.y + a.z * b.z;
 }
 
@@ -244,13 +245,13 @@ static inline csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b)
 }
 
 /// Complex 3D vector (geographic coordinates) scaling.
-static inline csphvec_t csphvec_scale(complex double c, const csphvec_t v) {
+static inline csphvec_t csphvec_scale(_Complex double c, const csphvec_t v) {
 	csphvec_t res = {c * v.rc, c * v.thetac, c * v.phic};
 	return res;
 }
 
 /// Complex 3D vector (geographic coordinates) "dot product" without conjugation.
-static inline complex double csphvec_dotnc(const csphvec_t a, const 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;
 }
@@ -262,7 +263,7 @@ static inline sph_t sph_scale(double c, const sph_t s) {
 }
 
 /// "Complex spherical" coordinate system scaling.
-static inline csph_t sph_cscale(complex double c, const sph_t s) {
+static inline csph_t sph_cscale(_Complex double c, const sph_t s) {
 	csph_t res = {c * s.r, s.theta, s.phi};
 	return res;
 }
@@ -358,8 +359,9 @@ static inline csph_t ccart2csph(const ccart3_t cart) {
 /// Real 3D cartesian to spherical (complex r) coordinates conversion. See @ref coord_conversions.
 static inline csph_t cart2csph(const cart3_t cart) {
 	csph_t sph;
-	sph.r = cart3norm(cart);
-	sph.theta = sph.r ? acos(cart.z / sph.r) : M_PI_2;
+  const double r = cart3norm(cart);
+	sph.r = r;
+	sph.theta = sph.r ? acos(cart.z / r) : M_PI_2;
 	sph.phi = atan2(cart.y, cart.x);
 	return sph;
 }
@@ -578,7 +580,6 @@ static inline double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t) {
 				" to 1D not allowed.");
 }
 
-
 /// Coordinate conversion of point arrays (something to something).
 /** The dest and src arrays must not overlap */
 static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t tdest,
@@ -589,31 +590,31 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
 				sph_t *d = (sph_t *) dest;
 				switch (tsrc & QPMS_COORDS_BITRANGE) {
 					case QPMS_COORDS_SPH: {
-						const sph_t *s = src;
+						const sph_t *s = (const sph_t*) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = s[i];
 						return; 
 					} break;
 					case QPMS_COORDS_CART3: {
-						const cart3_t *s = src;
+						const cart3_t *s = (const cart3_t*) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = cart2sph(s[i]);
 						return; 
 					} break;
 					case QPMS_COORDS_POL: {
-						const pol_t *s = src;
+						const pol_t *s = (const pol_t*) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = pol2sph_equator(s[i]);
 						return;
 					} break;
 					case QPMS_COORDS_CART2: {
-						const cart2_t *s = src;
+						const cart2_t *s = (const cart2_t*) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = cart22sph(s[i]);
 						return;
 					} break;
 					case QPMS_COORDS_CART1: {
-						const double *s = src;
+						const double *s = (const double *) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = cart12sph_zaxis(s[i]);
 						return; 
@@ -627,31 +628,31 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
 				cart3_t *d = (cart3_t *) dest;
 				switch (tsrc & QPMS_COORDS_BITRANGE) {
 					case QPMS_COORDS_SPH: {
-						const sph_t *s = src;
+						const sph_t *s = (const sph_t *) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = sph2cart(s[i]);
 						return;
 					} break;
 					case QPMS_COORDS_CART3: {
-						const cart3_t *s = src;
+						const cart3_t *s = (const cart3_t *) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = s[i];
 						return;
 					} break;
 					case QPMS_COORDS_POL: {
-						const pol_t *s = src;
+						const pol_t *s = (const pol_t *) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = pol2cart3_equator(s[i]);
 						return; 
 					} break;
 					case QPMS_COORDS_CART2: {
-						const cart2_t *s = src;
+						const cart2_t *s = (const cart2_t *) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = cart22cart3xy(s[i]);
 						return;
 					} break;
 					case QPMS_COORDS_CART1: {
-						const double *s = src;
+						const double *s = (const double *)src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = cart12cart3z(s[i]);
 						return; 
@@ -669,13 +670,13 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
 						QPMS_PR_ERROR("Implicit conversion from 3D to 2D coordinates not allowed"); 
 						break;
 					case QPMS_COORDS_POL: {
-						const pol_t *s = src;
+						const pol_t *s = (const pol_t *) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = s[i];
 						return;
 					} break;
 					case QPMS_COORDS_CART2: {
-						const cart2_t *s = src;
+						const cart2_t *s = (const cart2_t *) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = cart2pol(s[i]);
 						return;
@@ -696,13 +697,13 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
 						QPMS_PR_ERROR("Implicit conversion from 3D to 2D coordinates not allowed"); 
 						break;
 					case QPMS_COORDS_POL: {
-						const pol_t *s = src;
+						const pol_t *s = (const pol_t *) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = pol2cart(s[i]);
 						return; 
 					} break;
 					case QPMS_COORDS_CART2: {
-						const cart2_t *s = src;
+						const cart2_t *s = (const cart2_t *) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = s[i];
 						return;
@@ -727,7 +728,7 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
 						QPMS_PR_ERROR("Implicit conversion from 3D to 1D coordinates not allowed"); 
 						break;
 					case QPMS_COORDS_CART1: {
-						const double *s = src;
+						const double *s = (const double *) src;
 						for(size_t i = 0; i < nmemb; ++i)
 							d[i] = s[i];
 						return; 
@@ -899,6 +900,6 @@ static inline cart2_t cart2_from_double_array(const double a[]) {
 
 typedef double matrix3d[3][3];
 typedef double matrix2d[2][2];
-typedef complex double cmatrix3d[3][3];
-typedef complex double cmatrix2d[2][2];
+typedef _Complex double cmatrix3d[3][3];
+typedef _Complex double cmatrix2d[2][2];
 #endif //VECTORS_H

qpms/vswf.c

@@ -9,6 +9,7 @@
 #include <string.h>
 #include "qpms_error.h"
 #include "normalisation.h"
+#include <stdbool.h>
 
 
 qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() {
@@ -52,7 +53,7 @@ qpms_errno_t qpms_vswf_set_spec_append(qpms_vswf_set_spec_t *s, const qpms_uvswf
   return QPMS_SUCCESS;
 }
 
-bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
+_Bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
                     const qpms_vswf_set_spec_t *b) {
   if (a == b) return true;
   if (a->norm != b->norm) return false;
@@ -307,7 +308,7 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
   qpms_y_t nelem = qpms_lMax2nelem(lMax);
   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;
+  csphvec_t * const a1 = a, * const a2 = a1 + nelem, * const a3 = a1 + 2 * nelem;
   QPMS_ENSURE_SUCCESS(qpms_vecspharm_fill(a1, a2, a3, lMax, kr, norm));
   const csphvec_t *p1 = a1; 
   const csphvec_t *p2 = a2;
@@ -478,7 +479,7 @@ qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir_cart /*allow complex
 }
 
 qpms_errno_t qpms_incfield_planewave(complex double *target, const qpms_vswf_set_spec_t *bspec,
-    const cart3_t evalpoint, const void *args, bool add) {
+    const cart3_t evalpoint, const void *args, _Bool add) {
   QPMS_UNTESTED;
   const qpms_incfield_planewave_params_t *p = args;
 

qpms/vswf.h

@@ -7,6 +7,10 @@
  */
 #ifndef QPMS_VSWF_H
 #define QPMS_VSWF_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #include <unistd.h> // ssize_t 
 #include "qpms_types.h"
 #include <gsl/gsl_sf_legendre.h>
@@ -16,7 +20,7 @@
 /// Calculates the (regular VSWF) expansion coefficients of an external incident field.
 typedef qpms_errno_t (*qpms_incfield_t)(
         /// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n.
-        complex double *target,
+        _Complex double *target,
         const qpms_vswf_set_spec_t *bspec,
         const cart3_t evalpoint, ///< Point at which the VSWF expansion is made.
         const void *args, ///< Pointer to additional function-specific arguments.
@@ -82,7 +86,7 @@ qpms_errno_t qpms_uvswf_fill(
 /** SVWF coefficients in \a coeffs must be ordered according to \a setspec->ilist.
  */
 csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec,
-		const complex double *coeffs, ///< SVWF coefficient vector of size setspec->n.
+		const _Complex double *coeffs, ///< SVWF coefficient vector of size setspec->n.
 		csph_t kr, ///< Evaluation point.
 		qpms_bessel_t btyp);
 
@@ -118,7 +122,7 @@ typedef struct qpms_incfield_planewave_params_t {
  */
 qpms_errno_t qpms_incfield_planewave(
         /// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n.
-        complex double *target,
+        _Complex double *target,
         const qpms_vswf_set_spec_t *bspec,
         const cart3_t evalpoint, ///< Point at which the VSWF expansion is made.
         const void *args, ///< Pointer to additional function-specific arguments (converted to (const qpms_incfield_planewave_params_t *)).
@@ -227,23 +231,26 @@ qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *cons
 		qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm);
 
 qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir, ccart3_t amplitude,
-		complex double *targt_longcoeff, complex double *target_mgcoeff, complex double *target_elcoeff,
+		_Complex double *targt_longcoeff, _Complex double *target_mgcoeff, _Complex double *target_elcoeff,
 		qpms_l_t lMax, qpms_normalisation_t norm);
 qpms_errno_t qpms_planewave2vswf_fill_sph(sph_t wavedir, csphvec_t amplitude,
-		complex double *targt_longcoeff, complex double *target_mgcoeff, complex double *target_elcoeff,
+		_Complex double *targt_longcoeff, _Complex double *target_mgcoeff, _Complex double *target_elcoeff,
 		qpms_l_t lMax, qpms_normalisation_t norm);
 
 
 csphvec_t qpms_eval_vswf(sph_t where,
-		complex double *longcoeffs, complex double *mgcoeffs, complex double *elcoeffs,
+		_Complex double *longcoeffs, _Complex double *mgcoeffs, _Complex double *elcoeffs,
 		qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm);
 
 csphvec_t qpms_eval_vswf_csph(csph_t where,
-		complex double *longcoeffs, complex double *mgcoeffs, complex double *elcoeffs,
+		_Complex double *longcoeffs, _Complex double *mgcoeffs, _Complex double *elcoeffs,
 		qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm);
 
 qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
-		qpms_bessel_t btyp, qpms_normalisation_t norm);//NI
-void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *);//NI
+		qpms_bessel_t btyp, qpms_normalisation_t norm);
+void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *);
 
+#ifdef __cplusplus
+}
+#endif
 #endif // QPMS_VSWF_H

qpms/wigner.c

@@ -52,7 +52,7 @@ complex double qpms_wignerD_elem(const qpms_quat_t R,
 }
 
 complex double qpms_vswf_irot_elem_from_irot3(const qpms_irot3_t q,
-    const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m, bool pseudo) {
+    const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m, _Bool pseudo) {
 #ifndef NDEBUG
   qpms_irot3_checkdet(q);
 #endif

setup.py

@@ -111,8 +111,8 @@ cywaves = Extension('qpms.cywaves',
         libraries = common_libs,
         include_dirs = [numpy_includes], 
         )
-cytranslations = Extension('qpms.cytranslations',
-        sources = ['qpms/cytranslations.pyx',
+cytranslations = Extension('qpms.cytranslations_legacy',
+        sources = ['qpms/cytranslations_legacy.pyx',
             'qpms/translations_python.c',
             ],
         extra_compile_args=['-std=c99',
@@ -150,7 +150,7 @@ qpms_c = Extension('qpms.qpms_c',
         sources = ['qpms/qpms_c.pyx',],
         libraries = common_libs,
         include_dirs=['amos', 'qpms', numpy_includes],
-        extra_link_args=['-fopenmp'],
+        extra_link_args=['-fopenmp'], # needed?
         extra_compile_args=['-fopenmp'],
         )
 

tests/catch/CMakeLists.txt

@@ -0,0 +1,7 @@
+find_package(Catch2 REQUIRED)
+add_executable(catchtest all_includes.C t_vswf.C catch_aux.C)
+# TODO ensure that the tests are linked against the "working tree" version of qpms,
+# not the installed one
+target_link_libraries(catchtest PRIVATE Catch2::Catch2WithMain qpms lapacke)
+
+

tests/catch/all_includes.C

@@ -0,0 +1,30 @@
+// Just to test that all headers are C++ compatible
+#include <catch2/catch.hpp>
+#include <vswf.h>
+#include <beyn.h>
+#include <gaunt.h>
+#include <groups.h>
+#include <kahansum.h>
+#include <materials.h>
+#include <optim.h>
+#include <oshacks.h>
+#include <parsing.h>
+#include <qpms_error.h>
+#include <qpms_specfunc.h>
+#include <qpms_types.h>
+#include <scatsystem.h>
+#include <symmetries.h>
+#include <tiny_inlines.h>
+#include <tmatrices.h>
+#include <tolerances.h>
+#include <translations.h>
+#include <vectors.h>
+#include <pointgroups.h>
+#include <lattices_types.h>
+#include <indexing.h>
+#include <ewald.h>
+// C++ type conversion issues:
+#include <quaternions.h>
+//#include <lattices.h>
+//#include <normalisation.h>
+//#include <qpmsblas.h>

tests/catch/t_vswf.C

@@ -0,0 +1,44 @@
+#include <catch2/catch.hpp>
+#include <qpms_error.h>
+#include <vswf.h>
+#include "catch_aux.h"
+using namespace qpmstest;
+
+constexpr qpms_l_t lMax = 6;
+
+static std::vector<double> vswfvals_dvector(double r, double theta, double phi, 
+    qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm) 
+{
+  std::vector<double> vals(2 /*el and mg for now */ * (lMax * (lMax+2)) 
+      * 3 /* 3D vector components */ * 2 /*for real and imag*/);
+  qpms_vswf_fill( 
+      nullptr, (csphvec_t *) &(vals[0]), (csphvec_t *) &(vals[(lMax*(lMax+2)) * 3 * 2]),
+      lMax, sph_t{r, theta, phi},
+      btyp, norm);
+  return vals;
+}
+
+static std::vector<double> vswfvals_dvector_alternative(double r, double theta, double phi, 
+    qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm) 
+{
+  std::vector<double> vals(2 /*el and mg for now */ * (lMax * (lMax+2)) 
+      * 3 /* 3D vector components */ * 2 /*for real and imag*/);
+  qpms_vswf_fill_alternative( 
+      nullptr, (csphvec_t *) &(vals[0]), (csphvec_t *) &(vals[(lMax*(lMax+2)) * 3 * 2]),
+      lMax, sph_t{r, theta, phi},
+      btyp, norm);
+  return vals;
+}
+
+TEST_CASE("VSWF alternative implementation") {
+  using Catch::Matchers::Approx;
+  auto norms = GENERATE(QPMS_NORMALISATION_DEFAULT, QPMS_NORMALISATION_NORM_POWER,
+      QPMS_NORMALISATION_NORM_SPHARM, QPMS_NORMALISATION_CONVENTION_SCUFF);
+  auto thetas = GENERATE(0., M_PI_2, M_PI_2/2, M_PI, 3*M_PI_2, 0.1, -0.2, 4., 12.666);
+  auto phis = GENERATE(0., M_PI_2, M_PI_2/2, M_PI, 3*M_PI_2, 0.1, -0.2, 4., 12.666);
+  auto rs_positive = GENERATE(0.0001, 0.001, 0.01, 0.1, 1., 10., 100.);
+  CHECK_THAT( 
+        vswfvals_dvector(rs_positive, thetas, phis, lMax, QPMS_HANKEL_PLUS, norms),
+        Approx(vswfvals_dvector_alternative(rs_positive, thetas, phis, lMax, QPMS_HANKEL_PLUS, norms))
+      ) ;
+}

tests/test_ewald_consistency.py

@@ -0,0 +1,65 @@
+#!/usr/bin/env python3
+'''
+Test 2D Ewald summation consistency for basis VSWFs, using varying values of eta.
+'''
+import sys
+# Test parameters, rather arbitrary
+
+npoints = 10
+
+etafactors = [0.9, 1.1]
+
+zpointmax = 0. if len(sys.argv) == 1 else float(sys.argv[1])
+lattice_basis=[(580e-9,0),(0,580e-9)]
+positions = [(0,0)] # particle positions
+wavevector = [0.,0.,0.]
+bgparam = 1.52
+omega_eh = 2.1
+lMax = 2
+
+
+import numpy as np
+from qpms import Particle, CTMatrix, lorentz_drude, EpsMuGenerator, TMatrixGenerator, EwaldPart, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, EpsMu, dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType,eV, hbar, c
+from qpms.symmetries import point_group_info
+eh = eV/hbar
+np.random.seed(666)
+
+omega = omega_eh * eh
+
+bspec = BaseSpec(lMax=lMax)
+medium = EpsMuGenerator(bgparam)
+emg_dummy = EpsMuGenerator(10)
+tmg_dummy = TMatrixGenerator.sphere(medium, emg_dummy, 50e-9)
+
+particles = [Particle(pos, tmg_dummy, bspec) for pos in positions]
+
+ss, ssw = ScatteringSystem.create(particles, medium, omega, latticebasis=lattice_basis)
+wavevector = np.array(wavevector)
+
+sswk = ssw._sswk(wavevector)
+eta_orig = sswk.eta
+
+print("default eta: %g" % (eta_orig,))
+print("lattice basis:\n", ss.lattice_basis)
+
+points =  (np.random.rand(npoints) - .5)[:,None] * ss.lattice_basis[0] + (np.random.rand(npoints) - .5)[:,None] * ss.lattice_basis[1]
+points += (np.random.rand(npoints) - .5)[:,None] * np.array([0.,0.,zpointmax])
+
+
+fields_reference = sswk.scattered_field_basis(points)
+
+fields = np.empty((len(etafactors),) + fields_reference.shape, dtype=complex)
+
+fails = 0
+
+for n, etafactor in enumerate(etafactors):
+    sswk.eta = eta_orig * etafactor
+    fields[n] = sswk.scattered_field_basis(points)
+    print(fields[n])
+    print(fields[n]-fields_reference)
+    if np.sum(1-np.isclose(fields[n], fields_reference)) > 0: 
+        fails += 1
+
+print(fails)
+
+sys.exit(int(fails))