diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt
index 838b2e1..9073d04 100644
--- a/qpms/CMakeLists.txt
+++ b/qpms/CMakeLists.txt
@@ -24,7 +24,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 version.c
-	translations_dbg.c
+	translations_dbg.c scatsystem_dbg.c
 	)
 use_c99()
 
diff --git a/qpms/cycommon.pyx b/qpms/cycommon.pyx
index 31abcfa..feee540 100644
--- a/qpms/cycommon.pyx
+++ b/qpms/cycommon.pyx
@@ -144,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
diff --git a/qpms/cyewaldtest.pyx b/qpms/cyewaldtest.pyx
index 1bb7105..410f3da 100644
--- a/qpms/cyewaldtest.pyx
+++ b/qpms/cyewaldtest.pyx
@@ -1,6 +1,10 @@
 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)
@@ -10,6 +14,24 @@ cdef extern from "ewald.h":
     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
     cdef double[::1] err_view
@@ -31,6 +53,104 @@ 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)
diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd
index 45a18e0..7394265 100644
--- a/qpms/qpms_cdefs.pxd
+++ b/qpms/qpms_cdefs.pxd
@@ -193,10 +193,11 @@ 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
     # maybe more if needed
 
 # Point generators from lattices.h
@@ -259,95 +260,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:
@@ -628,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)
diff --git a/qpms/scatsystem_dbg.c b/qpms/scatsystem_dbg.c
new file mode 100644
index 0000000..66db149
--- /dev/null
+++ b/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;
+}
diff --git a/qpms/scatsystem_dbg.h b/qpms/scatsystem_dbg.h
new file mode 100644
index 0000000..a39707d
--- /dev/null
+++ b/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
+                );
+