diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx
index 7ea3b9c..0b5c6be 100644
--- a/qpms/qpms_c.pyx
+++ b/qpms/qpms_c.pyx
@@ -947,13 +947,16 @@ cdef class ScatteringSystem:
         return retdict
 
     @boundscheck(False)
-    def scattered_E(self, cdouble wavenumber, scatcoeffvector_full, evalpos, bint alt=False, btyp=BesselType.HANKEL_PLUS):
+    def scattered_E(self, cdouble wavenumber, scatcoeffvector_full, evalpos, btyp=BesselType.HANKEL_PLUS, bint alt=False):
+        # TODO DOC, periodic systems
+        if self.lattice_dimension != 0: # TODO enable also periodic systems
+            raise NotImplementedError("For periodic system, use _ScatteringSystemAtOmegaK.scattered_E")
         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 np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
-        cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, copy=False)
+        cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, copy=False, dtype=complex)
         cdef cdouble[::1] scv_view = scv
         cdef np.ndarray[complex, ndim=2] results = np.empty((evalpos_a.shape[0],3), dtype=complex)
         cdef ccart3_t res
@@ -1003,6 +1006,8 @@ def empty_lattice_modes_xy(EpsMu epsmu, reciprocal_basis, wavevector, double max
 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
@@ -1018,7 +1023,22 @@ cdef class _ScatteringSystemAtOmegaK:
             self.sswk.eta = eta
     
     @boundscheck(False)
-    def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS): # TODO DOC!!!
+    def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS):
+        """Evaluate electric field for a given excitation coefficient vector (periodic system)
+
+        Parameters
+        ----------
+        scatcoeffvector_full: array_like of shape (self.fecv_size,)
+            Onedimensional excitation coefficient vector, describing SVWFs leaving
+            the scatterers inside the representative unit cell.
+        evalpos: array_like of floats and shape (..., 3)
+            Evaluation points in cartesian coordinates.
+
+        Returns
+        -------
+        array_like of complex, with the same shape as `evalpos`
+            Electric field at the positions given in `evalpos`.
+        """
         if(btyp != QPMS_HANKEL_PLUS):
             raise NotImplementedError("Only first kind Bessel function-based fields are supported")
         cdef qpms_bessel_t btyp_c = BesselType(btyp)
@@ -1026,7 +1046,7 @@ cdef class _ScatteringSystemAtOmegaK:
         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[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, copy=False)
+        cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, dtype=complex, copy=False)
         cdef cdouble[::1] scv_view = scv
         cdef np.ndarray[complex, ndim=2] results = np.empty((evalpos_a.shape[0],3), dtype=complex)
         cdef ccart3_t res
@@ -1042,6 +1062,9 @@ cdef class _ScatteringSystemAtOmegaK:
                 results[i,1] = res.y
                 results[i,2] = res.z
         return results.reshape(evalpos.shape)
+    
+    property fecv_size: 
+        def __get__(self): return self.ssw_pyref.ss_pyref.fecv_size
 
 cdef class _ScatteringSystemAtOmega:
     '''
@@ -1097,6 +1120,7 @@ cdef class _ScatteringSystemAtOmega:
     def _sswk(self, k):
         cdef _ScatteringSystemAtOmegaK sswk = _ScatteringSystemAtOmegaK()
         sswk.sswk.ssw = self.ssw
+        sswk.ssw_pyref=self
         sswk.sswk.k[0] = k[0]
         sswk.sswk.k[1] = k[1]
         sswk.sswk.k[2] = k[2]
@@ -1164,13 +1188,39 @@ cdef class _ScatteringSystemAtOmega:
         return self.ss_pyref.translation_matrix_packed(wavenumber=self.wavenumber, iri=iri, J=J)
 
     @boundscheck(False)
-    def scattered_E(self, scatcoeffvector_full, evalpos, bint alt=False, btyp=QPMS_HANKEL_PLUS):
+    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
+        ----------
+        scatcoeffvector_full: array_like of shape (self.fecv_size,)
+            Onedimensional excitation coefficient vector, describing SVWFs leaving
+            the scatterers.
+        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.
+
+        Returns
+        -------
+        array_like of complex, with the same shape as `evalpos`
+            Electric field at the positions given in `evalpos`.
+
+        See Also
+        --------
+        _ScatteringSystemAtOmegaK.scattered_E: Evaluate scattered field in an infinite periodic system.
+        """
+        if(self.ssw[0].ss[0].lattice_dimension != 0):
+            if blochvector is None:
+                raise ValueError("For a periodic system, blochvector must be specified.")
+            return self._sswk(blochvector).scattered_E(scatcoeffvector_full, evalpos=evalpos, btyp=btyp)
         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 np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
-        cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, copy=False)
+        cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, dtype=complex, copy=False)
         cdef cdouble[::1] scv_view = scv
         cdef np.ndarray[complex, ndim=2] results = np.empty((evalpos_a.shape[0],3), dtype=complex)
         cdef ccart3_t res