diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt
index 2a63a3f..66036a8 100644
--- a/qpms/CMakeLists.txt
+++ b/qpms/CMakeLists.txt
@@ -3,6 +3,9 @@ find_package(GSL 2.0 REQUIRED)
 find_package(BLAS REQUIRED)
 find_package(LAPACK REQUIRED)
 
+# disable an annoying warning that gives false positives probably due to a bug in gcc
+# and other not very relevant warnings
+set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wno-int-in-bool-context -Wno-comment")
 
 #includes
 set (DIRS ${GSL_INCLUDE_DIRS} ${GSLCBLAS_INCLUDE_DIRS})
@@ -26,7 +29,7 @@ target_link_libraries (qpms
 	)
 target_include_directories (qpms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 
-target_compile_options(qpms PRIVATE -Wall -Wno-return-type)
+target_compile_options(qpms PRIVATE -Wall -Wno-return-type -Wno-unused-variable -Wno-unused-function -Wno-unused-but-set-variable -Wno-unused-label)
 target_compile_definitions(qpms PRIVATE LATTICESUMS32 QPMS_VECTORS_NICE_TRANSFORMATIONS
 	EWALD_AUTO_CUTOFF
 	)
diff --git a/qpms/cytmatrices.pxd b/qpms/cytmatrices.pxd
index 51641d9..3d39d65 100644
--- a/qpms/cytmatrices.pxd
+++ b/qpms/cytmatrices.pxd
@@ -1,5 +1,5 @@
 cimport numpy as np
-from .qpms_cdefs cimport qpms_tmatrix_t, cdouble, qpms_tmatrix_interpolator_t
+from .qpms_cdefs cimport qpms_tmatrix_t, cdouble, qpms_tmatrix_interpolator_t, qpms_tmatrix_generator_t, qpms_tmatrix_function_t
 from .cybspec cimport BaseSpec
 
 cdef class TMatrixInterpolator:
@@ -14,6 +14,27 @@ cdef class TMatrixInterpolator:
     cdef inline qpms_tmatrix_interpolator_t *rawpointer(self):
         return self.interp
 
+cdef class TMatrixGenerator:
+    cdef qpms_tmatrix_generator_t g
+    cdef object holder
+    cdef inline qpms_tmatrix_generator_t raw(self):
+        return self.g
+    cdef inline qpms_tmatrix_generator_t *rawpointer(self):
+        return &(self.g)
+
+cdef class TMatrixFunction:
+    cdef qpms_tmatrix_function_t f
+    cdef readonly TMatrixGenerator generator # reference holder
+    cdef readonly BaseSpec spec # reference holder
+    cdef inline qpms_tmatrix_function_t raw(self):
+        return self.f
+    cdef inline qpms_tmatrix_function_t *rawpointer(self):
+        return &self.f
+
+
+cdef class TMatrixGeneratorTransformed:
+    pass
+
 cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py!
     cdef readonly np.ndarray m # Numpy array holding the matrix data
     cdef readonly BaseSpec spec # Here we hold the base spec for the correct reference counting; TODO check if it gets copied
diff --git a/qpms/cytmatrices.pyx b/qpms/cytmatrices.pyx
index 63bbc14..940cd9e 100644
--- a/qpms/cytmatrices.pyx
+++ b/qpms/cytmatrices.pyx
@@ -228,13 +228,31 @@ cdef class __AxialSymParams:
         qpms_tmatrix_generator_axialsym_RQ_transposed_fill(&arrview[0][0], omega, &self.p, norm, QPMS_BESSEL_REGULAR)
         return arr
 
+cdef class TMatrixFunction:
+    '''
+    Wrapper over qpms_tmatrix_function_t. The main functional difference between this
+    and TMatrixGenerator is that this remembers a specific BaseSpec
+    and its __call__ method takes only one mandatory argument (in addition to self).
+    '''
+    def __init__(self, TMatrixGenerator tmg, BaseSpec spec):
+        self.generator = tmg
+        self.spec = spec
+        self.f.gen = self.generator.rawpointer()
+        self.f.spec = self.spec.rawpointer()
+
+    def __call__(self, cdouble omega, fill = None):
+        cdef CTMatrix tm
+        if fill is None: # make a new CTMatrix
+            tm = CTMatrix(self.spec, None)
+        else: # TODO check whether fill has the same bspec as self?
+            tm = fill
+        if self.f.gen.function(tm.rawpointer(), omega, self.f.gen.params) != 0:
+            raise ValueError("Something went wrong")
+        else:
+            return tm
 
 
 cdef class TMatrixGenerator:
-    cdef qpms_tmatrix_generator_t g
-    cdef object holder
-    cdef qpms_tmatrix_generator_t raw(self):
-        return self.g
     def __init__(self, what):
         if isinstance(what, __MieParams):
             self.holder = what
diff --git a/qpms/materials.h b/qpms/materials.h
index 0db7372..ce27d76 100644
--- a/qpms/materials.h
+++ b/qpms/materials.h
@@ -24,6 +24,12 @@ typedef struct qpms_epsmu_generator_t {
 	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) {
+	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.
 	const void *epsmu ///< Points to the qpms_epsmu_t to be returned.
diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx
index c78e764..5a6e203 100644
--- a/qpms/qpms_c.pyx
+++ b/qpms/qpms_c.pyx
@@ -12,7 +12,8 @@ from .cyquaternions cimport IRot3, CQuat
 from .cybspec cimport BaseSpec
 from .cycommon cimport make_c_string
 from .cycommon import string_c2py, PointGroupClass
-from .cytmatrices cimport CTMatrix
+from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator
+from .cymaterials cimport EpsMuGenerator
 from libc.stdlib cimport malloc, free, calloc
 import warnings
 
@@ -256,17 +257,37 @@ cdef class Particle:
     Wrapper over the qpms_particle_t structure.
     '''
     cdef qpms_particle_t p
-    cdef readonly CTMatrix t # We hold the reference to the T-matrix to ensure correct reference counting
+    cdef readonly TMatrixFunction f # Reference to ensure correct reference counting
 
-    def __cinit__(Particle self, pos, CTMatrix t):
+
+    def __cinit__(Particle self, pos, t, bspec = None):
+        cdef TMatrixGenerator tgen
+        cdef BaseSpec spec
         if(len(pos)>=2 and len(pos) < 4):
             self.p.pos.x = pos[0]
             self.p.pos.y = pos[1]
             self.p.pos.z = pos[2] if len(pos)==3 else 0
         else:
             raise ValueError("Position argument has to contain 3 or 2 cartesian coordinates")
-        self.t = t
-        self.p.tmatrix = self.t.rawpointer()
+        if isinstance(t, CTMatrix):
+            tgen = TMatrixGenerator(t)
+        elif isinstance(t, TMatrixGenerator):
+            tgen = <TMatrixGenerator>t
+        else: raise TypeError('t must be either CTMatrix or TMatrixGenerator, was %s' % str(type(t)))
+        if bspec is not None:
+            spec = bspec
+        else:
+            if isinstance(tgen.holder, CTMatrix):
+                spec = (<CTMatrix>tgen.holder).spec
+            else:
+                raise ValueError("bspec argument must be specified separately for str(type(t))")
+        self.f = TMatrixFunction(tgen, spec)
+        self.p.tmg = self.f.rawpointer()
+        # TODO non-trivial transformations later; if modified, do not forget to update ScatteringSystem constructor
+        self.p.op = qpms_tmatrix_operation_noop
+
+    def __dealloc__(self):
+        qpms_tmatrix_operation_clear(&self.p.op)
 
     cdef qpms_particle_t *rawpointer(Particle self):
         '''Pointer to the qpms_particle_p structure.
@@ -310,56 +331,109 @@ cpdef void scatsystem_set_nthreads(long n):
     qpms_scatsystem_set_nthreads(n)
     return
 
+
 cdef class ScatteringSystem:
     '''
     Wrapper over the C qpms_scatsys_t structure.
     '''
-    cdef list basespecs # Here we keep the references to occuring basespecs
+    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
 
-    def __cinit__(self, particles, FinitePointGroup sym):
-        '''TODO doc.
-        Takes the particles (which have to be a sequence of instances of Particle),
-        fills them together with their t-matrices to the "proto-qpms_scatsys_t"
-        orig and calls qpms_scatsys_apply_symmetry
-        (and then cleans orig)
-        '''
+    def check_s(self): # cdef instead?
+        if self.s == <qpms_scatsys_t *>NULL:
+            raise ValueError("ScatteringSystem's s-pointer not set. You must not use the default constructor; use the create() method instead")
+        #TODO is there a way to disable the constructor outside this module?
+
+    @staticmethod # We don't have any "standard" constructor for this right now
+    def create(particles, medium, FinitePointGroup sym, cdouble omega): # TODO tolerances
+        # These we are going to construct
+        cdef ScatteringSystem self
+        cdef _ScatteringSystemAtOmega pyssw
+        
         cdef qpms_scatsys_t orig # This should be automatically init'd to 0 (CHECKME)
-        cdef qpms_ss_pi_t p_count = len(particles)
-        cdef qpms_ss_tmi_t tm_count = 0
+        cdef qpms_ss_pi_t pi, p_count = len(particles)
+        cdef qpms_ss_tmi_t tmi, tm_count = 0
+        cdef qpms_ss_tmgi_t tmgi, tmg_count = 0
+
+        cdef qpms_scatsys_at_omega_t *ssw
+        cdef qpms_scatsys_t *ss
+
+        cdef Particle p
+        
+        tmgindices = dict()
+        tmgobjs = list()
         tmindices = dict()
-        tmobjs = list()
-        self.basespecs=list()
-        for p in particles: # find and enumerate unique t-matrices
-            if id(p.t) not in tmindices:
-                tmindices[id(p.t)] = tm_count
-                tmobjs.append(p.t)
+        tmlist = list()
+        for p in particles: # find and enumerate unique t-matrix generators
+            if p.p.op.typ != QPMS_TMATRIX_OPERATION_NOOP:
+                raise NotImplementedError("currently, only no-op T-matrix operations are allowed in ScatteringSystem constructor")
+            tmg_key = id(p.f)
+            if tmg_key not in tmgindices:
+                tmgindices[tmg_key] = tmg_count
+                tmgobjs.append(p.f) # Save the references on BaseSpecs and TMatrixGenerators (via TMatrixFunctions)
+                tmg_count += 1
+            # Following lines have to be adjusted when nontrivial operations allowed:
+            tm_derived_key = (tmg_key, None) # TODO unique representation of p.p.op instead of None
+            if tm_derived_key not in tmindices:
+                tmindices[tm_derived_key] = tm_count
+                tmlist.append(tm_derived_key)
                 tm_count += 1
+        cdef EpsMuGenerator mediumgen = EpsMuGenerator(medium)
+        orig.medium = mediumgen.g
+        orig.tmg_count = tmg_count
         orig.tm_count = tm_count
         orig.p_count = p_count
-        for tm in tmobjs: # create references to BaseSpec objects
-            self.basespecs.append(tm.spec)
         try:
-            orig.tm = <qpms_tmatrix_t **>malloc(orig.tm_count * sizeof(orig.tm[0]))
+            orig.tmg = <qpms_tmatrix_function_t *>malloc(orig.tmg_count * sizeof(orig.tmg[0]))
+            if not orig.tmg: raise MemoryError
+            orig.tm = <qpms_ss_derived_tmatrix_t *>malloc(orig.tm_count * sizeof(orig.tm[0]))
             if not orig.tm: raise MemoryError
             orig.p = <qpms_particle_tid_t *>malloc(orig.p_count * sizeof(orig.p[0]))
             if not orig.p: raise MemoryError
+            for tmgi in range(orig.tmg_count):
+                orig.tmg[tmgi] = (<TMatrixFunction?>tmgobjs[tmgi]).raw()
             for tmi in range(tm_count):
-                orig.tm[tmi] = (<CTMatrix?>(tmobjs[tmi])).rawpointer()
+                tm_derived_key = tmlist[tmi]
+                tmgi = tmgindices[tm_derived_key[0]]
+                orig.tm[tmi].tmgi = tmgi
+                orig.tm[tmi].op = qpms_tmatrix_operation_noop # TODO adjust when notrivial operations allowed
             for pi in range(p_count):
-                orig.p[pi].pos = (<Particle?>(particles[pi])).cval().pos
-                orig.p[pi].tmatrix_id = tmindices[id(particles[pi].t)]
-            self.s = qpms_scatsys_apply_symmetry(&orig, sym.rawpointer())
+                p = particles[pi]
+                tmg_key = id(p.f)
+                tm_derived_key = (tmg_key, None) # TODO unique representation of p.p.op instead of None
+                orig.p[pi].pos = p.cval().pos
+                orig.p[pi].tmatrix_id = tmindices[tm_derived_key]
+            ssw = qpms_scatsys_apply_symmetry(&orig, sym.rawpointer(), omega, &QPMS_TOLERANCE_DEFAULT)
+            ss = ssw[0].ss
         finally:
+            free(orig.tmg)
             free(orig.tm)
             free(orig.p)
+        self = ScatteringSystem()
+        self.medium_holder = mediumgen
+        self.s = ss
+        self.tmgobjs = tmgobjs
+        pyssw = _ScatteringSystemAtOmega()
+        pyssw.ssw = ssw
+        pyssw.ss_pyref = self
+        return self, pyssw
+
+    def __call__(self, cdouble omega):
+        self.check_s()
+        cdef _ScatteringSystemAtOmega pyssw = _ScatteringSystemAtOmega()
+        pyssw.ssw = qpms_scatsys_at_omega(self.s, omega)
+        pyssw.ss_pyref = self
+        return pyssw
 
     def __dealloc__(self):
-        qpms_scatsys_free(self.s)
+        if(self.s):
+            qpms_scatsys_free(self.s)
 
     property particles_tmi:
       def __get__(self):
+        self.check_s()
         r = list()
         cdef qpms_ss_pi_t pi
         for pi in range(self.s[0].p_count):
@@ -367,20 +441,27 @@ cdef class ScatteringSystem:
         return r
 
     property fecv_size: 
-        def __get__(self): return self.s[0].fecv_size
+        def __get__(self): 
+            self.check_s()
+            return self.s[0].fecv_size
     property saecv_sizes: 
         def __get__(self): 
+            self.check_s()
             return [self.s[0].saecv_sizes[i] 
                 for i in range(self.s[0].sym[0].nirreps)]
     property irrep_names: 
         def __get__(self): 
+            self.check_s()
             return [string_c2py(self.s[0].sym[0].irreps[iri].name) 
                     if (self.s[0].sym[0].irreps[iri].name) else None
                 for iri in range(self.s[0].sym[0].nirreps)]
     property nirreps: 
-        def __get__(self): return self.s[0].sym[0].nirreps
+        def __get__(self): 
+            self.check_s()
+            return self.s[0].sym[0].nirreps
 
     def pack_vector(self, vect, iri):
+        self.check_s()
         if len(vect) != self.fecv_size: 
             raise ValueError("Length of a full vector has to be %d, not %d" 
                     % (self.fecv_size, len(vect)))
@@ -392,6 +473,7 @@ cdef class ScatteringSystem:
         qpms_scatsys_irrep_pack_vector(&target_view[0], &vect_view[0], self.s, iri)
         return target_np
     def unpack_vector(self, packed, iri):
+        self.check_s()
         if len(packed) != self.saecv_sizes[iri]: 
             raise ValueError("Length of %d. irrep-packed vector has to be %d, not %d"
                     % (iri, self.saecv_sizes, len(packed)))
@@ -404,6 +486,7 @@ cdef class ScatteringSystem:
                 self.s, iri, 0)
         return target_np
     def pack_matrix(self, fullmatrix, iri):
+        self.check_s()
         cdef size_t flen = self.s[0].fecv_size
         cdef size_t rlen = self.saecv_sizes[iri]
         fullmatrix = np.array(fullmatrix, dtype=complex, copy=False, order='C')
@@ -418,6 +501,7 @@ cdef class ScatteringSystem:
                 self.s, iri)
         return target_np
     def unpack_matrix(self, packedmatrix, iri):
+        self.check_s()
         cdef size_t flen = self.s[0].fecv_size
         cdef size_t rlen = self.saecv_sizes[iri]
         packedmatrix = np.array(packedmatrix, dtype=complex, copy=False, order='C')
@@ -432,29 +516,8 @@ cdef class ScatteringSystem:
                 self.s, iri, 0)
         return target_np
 
-    def modeproblem_matrix_full(self, double k):
-        cdef size_t flen = self.s[0].fecv_size
-        cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
-                (flen,flen),dtype=complex, order='C')
-        cdef cdouble[:,::1] target_view = target
-        qpms_scatsys_build_modeproblem_matrix_full(&target_view[0][0], self.s, k)
-        return target
-
-    def modeproblem_matrix_packed(self, double k, qpms_iri_t iri, version='pR'):
-        cdef size_t rlen = self.saecv_sizes[iri]
-        cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
-                (rlen,rlen),dtype=complex, order='C')
-        cdef cdouble[:,::1] target_view = target
-        if (version == 'R'):
-            qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(&target_view[0][0], self.s, iri, k)
-        elif (version == 'pR'):
-          with nogil:
-            qpms_scatsys_build_modeproblem_matrix_irrep_packed(&target_view[0][0], self.s, iri, k)
-        else:
-            qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(&target_view[0][0], self.s, iri, k)
-        return target
-
     def translation_matrix_full(self, double k, J = QPMS_HANKEL_PLUS):
+        self.check_s()
         cdef size_t flen = self.s[0].fecv_size
         cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
                 (flen,flen),dtype=complex, order='C')
@@ -463,6 +526,7 @@ cdef class ScatteringSystem:
         return target
 
     def translation_matrix_packed(self, double k, qpms_iri_t iri, J = QPMS_HANKEL_PLUS):
+        self.check_s()
         cdef size_t rlen = self.saecv_sizes[iri]
         cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
                 (rlen,rlen),dtype=complex, order='C')
@@ -473,6 +537,7 @@ cdef class ScatteringSystem:
     
     property fullvec_psizes:
       def __get__(self):
+        self.check_s()
         cdef np.ndarray[int32_t, ndim=1] ar = np.empty((self.s[0].p_count,), dtype=np.int32)
         cdef int32_t[::1] ar_view = ar
         for pi in range(self.s[0].p_count):
@@ -482,6 +547,7 @@ cdef class ScatteringSystem:
 
     property fullvec_poffsets:
       def __get__(self):
+        self.check_s()
         cdef np.ndarray[intptr_t, ndim=1] ar = np.empty((self.s[0].p_count,), dtype=np.intp)
         cdef intptr_t[::1] ar_view = ar
         cdef intptr_t offset = 0
@@ -492,6 +558,7 @@ cdef class ScatteringSystem:
 
     property positions:
       def __get__(self):
+        self.check_s()
         cdef np.ndarray[np.double_t, ndim=2] ar = np.empty((self.s[0].p_count, 3), dtype=float)
         cdef np.double_t[:,::1] ar_view = ar
         for pi in range(self.s[0].p_count):
@@ -501,6 +568,7 @@ cdef class ScatteringSystem:
         return ar
    
     def planewave_full(self, k_cart, E_cart):
+        self.check_s()
         k_cart = np.array(k_cart)
         E_cart = np.array(E_cart)
         if k_cart.shape != (3,) or E_cart.shape != (3,):
@@ -520,7 +588,27 @@ cdef class ScatteringSystem:
                 self.s, qpms_incfield_planewave, <void *>&p, 0)
         return target_np
 
+cdef class _ScatteringSystemAtOmega:
+    '''
+    Wrapper over the C qpms_scatsys_at_omega_t structure
+    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")
+        self.ss_pyref.check_s()
+        #TODO is there a way to disable the constructor outside this module?
+
+    def __dealloc__(self):
+        if (self.ssw):
+            qpms_scatsys_at_omega_free(self.ssw)
+
     def apply_Tmatrices_full(self, a):
+        self.check()
         if len(a) != self.fecv_size: 
             raise ValueError("Length of a full vector has to be %d, not %d" 
                     % (self.fecv_size, len(a)))
@@ -529,31 +617,67 @@ cdef class ScatteringSystem:
         cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty(
                 (self.fecv_size,), dtype=complex, order='C')
         cdef cdouble[::1] target_view = target_np
-        qpms_scatsys_apply_Tmatrices_full(&target_view[0], &a_view[0], self.s)
+        qpms_scatsysw_apply_Tmatrices_full(&target_view[0], &a_view[0], self.ssw)
         return target_np
     
-    cdef qpms_scatsys_t *rawpointer(self):
-        return self.s
+    cdef qpms_scatsys_at_omega_t *rawpointer(self):
+        return self.ssw
+
+    def scatter_solver(self, iri=None):
+        self.check()
+        return ScatteringMatrix(self, iri)
+
+    property fecv_size: 
+        def __get__(self): return self.ss_pyref.fecv_size
+    property saecv_sizes: 
+        def __get__(self): return self.ss_pyref.saecv_sizes
+    property irrep_names: 
+        def __get__(self): return self.ss_pyref.irrep_names
+    property nirreps: 
+        def __get__(self): return self.ss_pyref.nirreps
+
+    def modeproblem_matrix_full(self):
+        self.check()
+        cdef size_t flen = self.ss_pyref.s[0].fecv_size
+        cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
+                (flen,flen),dtype=complex, order='C')
+        cdef cdouble[:,::1] target_view = target
+        qpms_scatsysw_build_modeproblem_matrix_full(&target_view[0][0], self.ssw)
+        return target
+
+    def modeproblem_matrix_packed(self, qpms_iri_t iri, version='pR'):
+        self.check()
+        cdef size_t rlen = self.saecv_sizes[iri]
+        cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
+                (rlen,rlen),dtype=complex, order='C')
+        cdef cdouble[:,::1] target_view = target
+        if (version == 'R'):
+            qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(&target_view[0][0], self.ssw, iri)
+        elif (version == 'pR'):
+          with nogil:
+            qpms_scatsysw_build_modeproblem_matrix_irrep_packed(&target_view[0][0], self.ssw, iri)
+        else:
+            qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(&target_view[0][0], self.ssw, iri)
+        return target
 
-    def scatter_solver(self, double k, iri=None):
-        return ScatteringMatrix(self, k, iri)
 
 cdef class ScatteringMatrix:
     '''
     Wrapper over the C qpms_ss_LU structure that keeps the factorised mode problem matrix.
     '''
-    cdef ScatteringSystem ss # Here we keep the reference to the parent scattering system
+    cdef _ScatteringSystemAtOmega ssw # Here we keep the reference to the parent scattering system
     cdef qpms_ss_LU lu
 
-    def __cinit__(self, ScatteringSystem ss, double k, iri=None):
-        self.ss = ss
+    def __cinit__(self, _ScatteringSystemAtOmega ssw, iri=None):
+        ssw.check()
+        self.ssw = ssw
         # TODO? pre-allocate the matrix with numpy to make it transparent?
         if iri is None:
-            self.lu = qpms_scatsys_build_modeproblem_matrix_full_LU(
-                    NULL, NULL, ss.rawpointer(), k)
+            self.lu = qpms_scatsysw_build_modeproblem_matrix_full_LU(
+                    NULL, NULL, ssw.rawpointer())
         else:
-            self.lu = qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU(
-                    NULL, NULL, ss.rawpointer(), iri, k)
+            self.lu = qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
+                    NULL, NULL, ssw.rawpointer(), iri)
 
     def __dealloc__(self):
         qpms_ss_LU_free(self.lu)
@@ -566,13 +690,13 @@ cdef class ScatteringMatrix:
         cdef size_t vlen
         cdef qpms_iri_t iri = -1;
         if self.lu.full:
-            vlen = self.lu.ss[0].fecv_size
+            vlen = self.lu.ssw[0].ss[0].fecv_size
             if len(a_inc) != vlen:
                 raise ValueError("Length of a full coefficient vector has to be %d, not %d"
                         % (vlen, len(a_inc)))
         else:
             iri = self.lu.iri
-            vlen = self.lu.ss[0].saecv_sizes[iri]
+            vlen = self.lu.ssw[0].ss[0].saecv_sizes[iri]
             if len(a_inc) != vlen:
                 raise ValueError("Length of a %d. irrep packed coefficient vector has to be %d, not %d"
                         % (iri, vlen, len(a_inc)))
diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd
index 7e9bcd7..6df0435 100644
--- a/qpms/qpms_cdefs.pxd
+++ b/qpms/qpms_cdefs.pxd
@@ -102,6 +102,7 @@ cdef extern from "qpms_types.h":
         QPMS_VSWF_ELECTRIC
         QPMS_VSWF_MAGNETIC
         QPMS_VSWF_LONGITUDINAL
+    ctypedef int32_t qpms_ss_tmgi_t
     ctypedef int32_t qpms_ss_tmi_t
     ctypedef int32_t qpms_ss_pi_t
     ctypedef int qpms_gmi_t
@@ -150,6 +151,11 @@ cdef extern from "qpms_error.h":
     qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types)
     qpms_dbgmsg_flags qpms_dbgmsg_disable(qpms_dbgmsg_flags types)
 
+cdef extern from "tolerances.h":
+    struct qpms_tolerance_spec_t:
+        pass # TODO
+    const qpms_tolerance_spec_t QPMS_TOLERANCE_DEFAULT
+
 
 # This is due to the fact that cython apparently cannot nest the unnamed struct/unions in an obvious way
 ctypedef union qpms_incfield_planewave_params_k:
@@ -476,6 +482,24 @@ cdef extern from "tmatrices.h":
             double omega, cdouble epsilon_fg, cdouble epsilon_bg)
     qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(cdouble *target, cdouble omega,
             const qpms_tmatrix_generator_axialsym_param_t *param, qpms_normalisation_t norm, qpms_bessel_t J)
+    struct qpms_tmatrix_function_t:
+        const qpms_vswf_set_spec_t *spec
+        const qpms_tmatrix_generator_t *gen
+    ctypedef enum qpms_tmatrix_operation_kind_t:
+        QPMS_TMATRIX_OPERATION_NOOP
+        QPMS_TMATRIX_OPERATION_LRMATRIX
+        QPMS_TMATRIX_OPERATION_IROT3
+        QPMS_TMATRIX_OPERATION_IROT3ARR
+        QPMS_TMATRIX_OPERATION_COMPOSE_SUM
+        QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN
+        QPMS_TMATRIX_OPERATION_SCMULZ
+        QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE
+    
+    struct qpms_tmatrix_operation_t:
+        qpms_tmatrix_operation_kind_t typ
+        pass # TODO add the op union later if needed
+    const qpms_tmatrix_operation_t qpms_tmatrix_operation_noop
+    void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *)
 
 cdef extern from "pointgroups.h":
     bint qpms_pg_is_finite_axial(qpms_pointgroup_class cls)
@@ -496,12 +520,19 @@ cdef extern from "scatsystem.h":
     void qpms_scatsystem_set_nthreads(long n)
     struct qpms_particle_t:
         cart3_t pos
-        const qpms_tmatrix_t *tmatrix
+        const qpms_tmatrix_function_t *tmg
+        qpms_tmatrix_operation_t op
     struct qpms_particle_tid_t:
         cart3_t pos
         qpms_ss_tmi_t tmatrix_id
+    struct qpms_ss_derived_tmatrix_t:
+        qpms_ss_tmgi_t tmgi
+        qpms_tmatrix_operation_t op
     struct qpms_scatsys_t:
-        qpms_tmatrix_t **tm
+        qpms_epsmu_generator_t medium
+        qpms_tmatrix_function_t *tmg
+        qpms_ss_tmgi_t tmg_count
+        qpms_ss_derived_tmatrix_t *tm
         qpms_ss_tmi_t tm_count
         qpms_particle_tid_t *p
         qpms_ss_pi_t p_count
@@ -509,10 +540,19 @@ cdef extern from "scatsystem.h":
         size_t fecv_size
         size_t *saecv_sizes
         const qpms_finite_group_t *sym
-    qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym)
     void qpms_scatsys_free(qpms_scatsys_t *s)
     qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path) #NI
     qpms_scatsys_t *qpms_scatsys_load(char *path) #NI
+    struct qpms_scatsys_at_omega_t:
+        const qpms_scatsys_t *ss
+        qpms_tmatrix_t **tm,
+        cdouble omega
+        qpms_epsmu_t medium
+        cdouble wavenumber
+    qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym,
+            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_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, 
@@ -521,41 +561,43 @@ cdef extern from "scatsystem.h":
             const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
     cdouble *qpms_scatsys_irrep_unpack_vector(cdouble *target_full,
             const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add)
-    cdouble *qpms_scatsys_build_modeproblem_matrix_full(cdouble *target,
-            const qpms_scatsys_t *ss, cdouble k)
+    cdouble *qpms_scatsysw_build_modeproblem_matrix_full(cdouble *target,
+            const qpms_scatsys_at_omega_t *ssw)
     cdouble *qpms_scatsys_build_translation_matrix_full(cdouble *target,
             const qpms_scatsys_t *ss, cdouble k)
     cdouble *qpms_scatsys_build_translation_matrix_e_full(cdouble *target,
             const qpms_scatsys_t *ss, cdouble k, qpms_bessel_t J)
-    cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed(cdouble *target,
-            const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) nogil
+    cdouble *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(cdouble *target,
+            const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) nogil
     cdouble *qpms_scatsys_build_translation_matrix_e_irrep_packed(cdouble *target,
             const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k, qpms_bessel_t J) nogil
-    cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
-            cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) nogil
-    cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(
-            cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) nogil
+    cdouble *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
+            cdouble *target, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) nogil
+    cdouble *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
+            cdouble *target, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) nogil
     cdouble *qpms_scatsys_incident_field_vector_full(cdouble *target_full,
             const qpms_scatsys_t *ss, qpms_incfield_t field_at_point, 
             const void *args, bint add)
-    cdouble *qpms_scatsys_apply_Tmatrices_full(cdouble *target_full, const cdouble *inc_full,
-            const qpms_scatsys_t *ss)
+    cdouble *qpms_scatsysw_apply_Tmatrices_full(cdouble *target_full, const cdouble *inc_full,
+            const qpms_scatsys_at_omega_t *ssw)
     struct qpms_ss_LU:
-        const qpms_scatsys_t *ss
+        const qpms_scatsys_at_omega_t *ssw
         bint full
         qpms_iri_t iri
         cdouble *a
         int *ipiv
     void qpms_ss_LU_free(qpms_ss_LU lu)
-    qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU(cdouble *target,
-            int *target_piv, const qpms_scatsys_t *ss, cdouble k)
-    qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU(cdouble *target,
-            int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k)
-    qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(cdouble *modeproblem_matrix_full,
-            int *target_piv, const qpms_scatsys_t *ss)
-    qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(cdouble *modeproblem_matrix_full,
-            int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri)
+    qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(cdouble *target,
+            int *target_piv, const qpms_scatsys_at_omega_t *ssw)
+    qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(cdouble *target,
+            int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri)
+    qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(cdouble *modeproblem_matrix_full,
+            int *target_piv, const qpms_scatsys_at_omega_t *ssw)
+    qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(cdouble *modeproblem_matrix_full,
+            int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri)
     cdouble *qpms_scatsys_scatter_solve(cdouble *target_f, const cdouble *a_inc, qpms_ss_LU ludata)
+    const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(const qpms_scatsys_t *ss, qpms_ss_tmi_t tmi)
+    const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi)
 
 cdef extern from "ewald.h":
     struct qpms_csf_result:
diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h
index 10fb7b8..393157c 100644
--- a/qpms/qpms_error.h
+++ b/qpms/qpms_error.h
@@ -56,6 +56,15 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
 			(size_t) (size));\
 }
 
+#define QPMS_CRASHING_MALLOCPY(dest, src, size) {\
+	(dest) = malloc(size);\
+       	if(QPMS_UNLIKELY(!(dest) && (size)))\
+       		qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__,\
+			"Allocation of %zd bytes for " #dest " failed.",\
+			(size_t) (size));\
+	memcpy((dest), (src), (size));\
+}
+
 #define QPMS_CRASHING_CALLOC(pointer, nmemb, size) {\
 	(pointer) = calloc((nmemb), (size));\
        	if(QPMS_UNLIKELY(!(pointer) && (nmemb) && (size)))\
diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h
index fd73ce0..1f708de 100644
--- a/qpms/qpms_types.h
+++ b/qpms/qpms_types.h
@@ -297,6 +297,9 @@ typedef struct qpms_vswf_set_spec_t {
 /// T-matrix index used in qpms_scatsys_t and related structures.
 typedef int32_t qpms_ss_tmi_t;
 
+/// T-matrix generator index used in qpms_scatsys_t and related structures.
+typedef int32_t qpms_ss_tmgi_t;
+
 /// Particle index used in qpms_scatsys_t and related structures.
 typedef int32_t qpms_ss_pi_t;
 
@@ -407,6 +410,7 @@ typedef struct qpms_epsmu_t {
 	complex double mu; ///< Relative permeability.
 } qpms_epsmu_t;
 
+struct qpms_tolerance_spec_t; // See tolerances.h
 
 #define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
 #endif // QPMS_TYPES
diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c
index a697d0f..febef46 100644
--- a/qpms/scatsystem.c
+++ b/qpms/scatsystem.c
@@ -21,6 +21,7 @@
 #include "tmatrices.h"
 #include <pthread.h>
 #include "kahansum.h"
+#include "tolerances.h"
 
 #ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS
 #include "qpmsblas.h"
@@ -80,7 +81,7 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu
   qpms_ss_orbit_type_t * const ot_new = & (ss->orbit_types[ss->orbit_type_count]);
   ot_new->size = ot_current->size;
 
-  const qpms_vswf_set_spec_t *bspec = ss->tm[ot_current->tmatrices[0]]->spec;
+  const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_tmi(ss, ot_current->tmatrices[0]);
   const size_t bspecn = bspec->n;
   ot_new->bspecn = bspecn;
   
@@ -143,7 +144,9 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu
 
 
 // Almost 200 lines. This whole thing deserves a rewrite!
-qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym) {
+qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig,
+    const qpms_finite_group_t *sym, complex double omega, 
+    const qpms_tolerance_spec_t *tol) {
   // TODO check data sanity
 
   qpms_l_t lMax = 0; // the overall lMax of all base specs.
@@ -171,49 +174,82 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
     for (qpms_ss_pi_t j = 0; j < i; ++j)
       assert(!cart3_isclose(orig->p[i].pos, orig->p[j].pos, 0, QPMS_SCATSYS_LEN_RTOL * lenscale));
 
+
   // Allocate T-matrix, particle and particle orbit info arrays
-  qpms_scatsys_t *ss = malloc(sizeof(qpms_scatsys_t));
+  qpms_scatsys_t *ss;
+  QPMS_CRASHING_MALLOC(ss, sizeof(qpms_scatsys_t));     
   ss->lenscale = lenscale;
   ss->sym = sym;
+  ss->medium = orig->medium;
+
+  // Copy the qpms_tmatrix_fuction_t from orig
+  ss->tmg_count = orig->tmg_count;
+  QPMS_CRASHING_MALLOC(ss->tmg, ss->tmg_count * sizeof(*(ss->tmg)));
+  memcpy(ss->tmg, orig->tmg, ss->tmg_count * sizeof(*(ss->tmg)));
 
   ss->tm_capacity = sym->order * orig->tm_count;
-  ss->tm = malloc(ss->tm_capacity * sizeof(qpms_tmatrix_t *));
+  QPMS_CRASHING_MALLOC(ss->tm, ss->tm_capacity * sizeof(*(ss->tm)));
 
   ss->p_capacity = sym->order * orig->p_count;
-  ss->p = malloc(ss->p_capacity * sizeof(qpms_particle_tid_t));
-  ss->p_orbitinfo = malloc(ss->p_capacity * sizeof(qpms_ss_particle_orbitinfo_t));
+  QPMS_CRASHING_MALLOC(ss->p, ss->p_capacity * sizeof(qpms_particle_tid_t));
+  QPMS_CRASHING_MALLOC(ss->p_orbitinfo, ss->p_capacity * sizeof(qpms_ss_particle_orbitinfo_t));
   for (qpms_ss_pi_t pi = 0; pi < ss->p_capacity; ++pi) {
     ss->p_orbitinfo[pi].t = QPMS_SS_P_ORBITINFO_UNDEF;
     ss->p_orbitinfo[pi].p = QPMS_SS_P_ORBITINFO_UNDEF;
   }
 
-  // Copy T-matrices; checking for duplicities
+  // Evaluate the original T-matrices at omega
+  qpms_tmatrix_t **tm_orig_omega; 
+  QPMS_CRASHING_MALLOC(tm_orig_omega, orig->tmg_count * sizeof(*tm_orig_omega));
+  for(qpms_ss_tmgi_t i = 0; i < orig->tmg_count; ++i) 
+    tm_orig_omega[i] = qpms_tmatrix_init_from_function(orig->tmg[i], omega);
+
+  // Evaluate the medium and derived T-matrices at omega.
+  qpms_scatsys_at_omega_t *ssw;
+  QPMS_CRASHING_MALLOC(ssw, sizeof(*ssw)); // returned
+  ssw->ss = ss;
+  ssw->omega = omega;
+  ssw->medium = qpms_epsmu_generator_eval(ss->medium, omega);
+  ssw->wavenumber = qpms_wavenumber(omega, ssw->medium);
+  // we will be using ss->tm_capacity also for ssw->tm
+  QPMS_CRASHING_MALLOC(ssw->tm, ss->tm_capacity * sizeof(*(ssw->tm))); // returned
+
+  // Evaluate T-matrices at omega; checking for duplicities
   
   ss->max_bspecn = 0; // We'll need it later.for memory alloc estimates.
 
-  qpms_ss_tmi_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities
+  qpms_ss_tmi_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities; VLA!
   ss->tm_count = 0;
   for (qpms_ss_tmi_t i = 0; i < orig->tm_count; ++i) {
+    qpms_tmatrix_t *ti = qpms_tmatrix_apply_operation(&(orig->tm[i].op), tm_orig_omega[orig->tm[i].tmgi]);
     qpms_ss_tmi_t j;
     for (j = 0; j < ss->tm_count; ++j) 
-      if (qpms_tmatrix_isclose(orig->tm[i], ss->tm[j], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL)) {
+      if (qpms_tmatrix_isclose(ti, ssw->tm[j], tol->rtol, tol->atol)) {
         break;
       }
-    if (j == ss->tm_count) { // duplicity not found, copy the t-matrix
-      ss->tm[j] = qpms_tmatrix_copy(orig->tm[i]);
-      ss->max_bspecn = MAX(ss->tm[j]->spec->n, ss->max_bspecn);
-      lMax = MAX(lMax, ss->tm[j]->spec->lMax);
+    if (j == ss->tm_count) { // duplicity not found, save both the "abstract" and "at omega" T-matrices
+      qpms_tmatrix_operation_copy(&ss->tm[j].op, &orig->tm[i].op);
+      ss->tm[j].tmgi = orig->tm[i].tmgi; // T-matrix functions are preserved.
+      ssw->tm[j] = ti;
+      ss->max_bspecn = MAX(ssw->tm[j]->spec->n, ss->max_bspecn);
+      lMax = MAX(lMax, ssw->tm[j]->spec->lMax);
       ++(ss->tm_count);
     } 
+    else qpms_tmatrix_free(ti);
     tm_dupl_remap[i] = j;
     if (normalisation == QPMS_NORMALISATION_UNDEF)
-      normalisation = ss->tm[i]->spec->norm;
+      normalisation = ssw->tm[i]->spec->norm;
     // We expect all bspec norms to be the same.
-    else QPMS_ENSURE(normalisation == ss->tm[j]->spec->norm,
+    else QPMS_ENSURE(normalisation == ssw->tm[j]->spec->norm,
         "Normalisation convention must be the same for all T-matrices."
-        " %d != %d\n", normalisation, ss->tm[j]->spec->norm);
+        " %d != %d\n", normalisation, ssw->tm[j]->spec->norm);
   }
 
+  // Free the original T-matrices at omega
+  for(qpms_ss_tmgi_t i = 0; i < orig->tmg_count; ++i)
+    qpms_tmatrix_free(tm_orig_omega[i]);
+  free(tm_orig_omega);
+
   // Copy particles, remapping the t-matrix indices
   for (qpms_ss_pi_t i = 0; i < orig->p_count; ++(i)) {
     ss->p[i] = orig->p[i];
@@ -222,44 +258,57 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
   ss->p_count = orig->p_count;
 
   // allocate t-matrix symmetry map
-  ss->tm_sym_map = malloc(sizeof(qpms_ss_tmi_t) * sym->order * sym->order * ss->tm_count);
+  QPMS_CRASHING_MALLOC(ss->tm_sym_map, sizeof(qpms_ss_tmi_t) * sym->order * sym->order * ss->tm_count);
 
   // Extend the T-matrices list by the symmetry operations
   for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi) 
     for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi){
-      const size_t d = ss->tm[tmi]->spec->n;
-      complex double M[d][d]; // transformation matrix
-      qpms_irot3_uvswfi_dense(M[0], ss->tm[tmi]->spec, sym->rep3d[gmi]);
-      qpms_tmatrix_t *transformed = qpms_tmatrix_apply_symop(ss->tm[tmi], M[0]);
+      const size_t d = ssw->tm[tmi]->spec->n;
+      complex double *m;
+      QPMS_CRASHING_MALLOC(m, d*d*sizeof(complex double)); // ownership passed to ss->tm[ss->tm_count].op
+      qpms_irot3_uvswfi_dense(m, ssw->tm[tmi]->spec, sym->rep3d[gmi]);
+      qpms_tmatrix_t *transformed = qpms_tmatrix_apply_symop(ssw->tm[tmi], m);
       qpms_ss_tmi_t tmj;
       for (tmj = 0; tmj < ss->tm_count; ++tmj)
-        if (qpms_tmatrix_isclose(transformed, ss->tm[tmj], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL))
+        if (qpms_tmatrix_isclose(transformed, ssw->tm[tmj], tol->rtol, tol->atol))
           break;
       if (tmj < ss->tm_count) { // HIT, transformed T-matrix already exists
+        //TODO some "rounding error cleanup" symmetrisation could be performed here?
         qpms_tmatrix_free(transformed);
-      } else { // MISS, save the matrix and increment the count
-        ss->tm[ss->tm_count] = transformed;
+      } else { // MISS, save the matrix (also the "abstract" one)
+        ssw->tm[ss->tm_count] = transformed;
+	ss->tm[ss->tm_count].tmgi = ss->tm[tmi].tmgi;
+        qpms_tmatrix_operation_compose_chain_init(&(ss->tm[ss->tm_count].op), 2, 1);
+        struct qpms_tmatrix_operation_compose_chain * const o = &(ss->tm[ss->tm_count].op.op.compose_chain);
+        o->ops[0] = & ss->tm[tmi].op; // Let's just borrow this
+        o->ops_owned[0] = false;
+        o->opmem[0].typ = QPMS_TMATRIX_OPERATION_LRMATRIX;
+        o->opmem[0].op.lrmatrix.m = m;
+        o->opmem[0].op.lrmatrix.m_size = d * d;
+        o->ops[1] = o->opmem;
+        o->ops_owned[1] = true;
         ++(ss->tm_count);
       }
       ss->tm_sym_map[gmi + tmi * sym->order] = tmj; // In any case, tmj now indexes the correct transformed matrix
     }
   // Possibly free some space using the new ss->tm_count instead of (old) ss->tm_count*sym->order
-  ss->tm_sym_map = realloc(ss->tm_sym_map, sizeof(qpms_ss_tmi_t) * sym->order * ss->tm_count);
+  QPMS_CRASHING_REALLOC(ss->tm_sym_map, sizeof(qpms_ss_tmi_t) * sym->order * ss->tm_count);
   // tm could be realloc'd as well, but those are just pointers, not likely many.
  
 
   // allocate particle symmetry map
-  ss->p_sym_map = malloc(sizeof(qpms_ss_pi_t) * sym->order * sym->order * ss->p_count);
+  QPMS_CRASHING_MALLOC(ss->p_sym_map, sizeof(qpms_ss_pi_t) * sym->order * sym->order * ss->p_count);
   // allocate orbit type array (TODO realloc later if too long)
   ss->orbit_type_count = 0;
-  ss->orbit_types = calloc(ss->p_count, sizeof(qpms_ss_orbit_type_t));
+  QPMS_CRASHING_CALLOC(ss->orbit_types, ss->p_count, sizeof(qpms_ss_orbit_type_t));
 
-  ss->otspace_end = ss->otspace = malloc( // reallocate later
+  QPMS_CRASHING_MALLOC(ss->otspace, // reallocate later
       (sizeof(qpms_ss_orbit_pi_t) * sym->order * sym->order
        + sizeof(qpms_ss_tmi_t) * sym->order
        + 3*sizeof(size_t) * sym->nirreps
        + sizeof(complex double) * SQ(sym->order * ss->max_bspecn)) * ss->p_count
        );
+  ss->otspace_end = ss->otspace;
   
   // Workspace for the orbit type determination
   qpms_ss_orbit_type_t ot_current;
@@ -348,15 +397,15 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
     }
   }
   // Possibly free some space using the new ss->p_count instead of (old) ss->p_count*sym->order
-  ss->p_sym_map = realloc(ss->p_sym_map, sizeof(qpms_ss_pi_t) * sym->order * ss->p_count);
-  ss->p = realloc(ss->p, sizeof(qpms_particle_tid_t) * ss->p_count); 
-  ss->p_orbitinfo = realloc(ss->p_orbitinfo, sizeof(qpms_ss_particle_orbitinfo_t)*ss->p_count);
+  QPMS_CRASHING_REALLOC(ss->p_sym_map, sizeof(qpms_ss_pi_t) * sym->order * ss->p_count);
+  QPMS_CRASHING_REALLOC(ss->p, sizeof(qpms_particle_tid_t) * ss->p_count); 
+  QPMS_CRASHING_REALLOC(ss->p_orbitinfo, sizeof(qpms_ss_particle_orbitinfo_t)*ss->p_count);
   ss->p_capacity = ss->p_count;
 
   {  // Reallocate the orbit type data space and update the pointers if needed.
     size_t otspace_sz = ss->otspace_end - ss->otspace;
     char *old_otspace = ss->otspace;
-    ss->otspace = realloc(ss->otspace, otspace_sz);
+    QPMS_CRASHING_REALLOC(ss->otspace, otspace_sz);
     ptrdiff_t shift = ss->otspace - old_otspace;
     if(shift) {
       for (size_t oi = 0; oi < ss->orbit_type_count; ++oi) {
@@ -373,15 +422,14 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
 
   // Set ss->fecv_size and ss->fecv_pstarts
   ss->fecv_size = 0;
-  ss->fecv_pstarts = malloc(ss->p_count * sizeof(size_t));
+  QPMS_CRASHING_MALLOC(ss->fecv_pstarts, ss->p_count * sizeof(size_t));
   for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
     ss->fecv_pstarts[pi] = ss->fecv_size;
-    ss->fecv_size += ss->tm[ss->p[pi].tmatrix_id]->spec->n; // That's a lot of dereferencing!
+    ss->fecv_size += ssw->tm[ss->p[pi].tmatrix_id]->spec->n; // That's a lot of dereferencing!
   }
 
-  ss->saecv_sizes = malloc(sizeof(size_t) * sym->nirreps); if (!ss->saecv_sizes) abort();
-  ss->saecv_ot_offsets = malloc(sizeof(size_t) * sym->nirreps * ss->orbit_type_count);
-  if (!ss->saecv_ot_offsets) abort();
+  QPMS_CRASHING_MALLOC(ss->saecv_sizes, sizeof(size_t) * sym->nirreps);
+  QPMS_CRASHING_MALLOC(ss->saecv_ot_offsets, sizeof(size_t) * sym->nirreps * ss->orbit_type_count);
   for(qpms_iri_t iri = 0; iri < sym->nirreps; ++iri) {
     ss->saecv_sizes[iri] = 0;
     for(qpms_ss_oti_t oti = 0; oti < ss->orbit_type_count; ++oti) {
@@ -408,13 +456,14 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
   }
   
   ss->c = qpms_trans_calculator_init(lMax, normalisation);
-  return ss;
+  return ssw;
 }
 
 
 void qpms_scatsys_free(qpms_scatsys_t *ss) {
   if(ss) {
     free(ss->tm);
+    free(ss->tmg);
     free(ss->p);
     free(ss->fecv_pstarts);
     free(ss->tm_sym_map);
@@ -429,7 +478,59 @@ void qpms_scatsys_free(qpms_scatsys_t *ss) {
   free(ss);
 }
 
+void qpms_scatsys_at_omega_refill(qpms_scatsys_at_omega_t *ssw, 
+    complex double omega) {
+  const qpms_scatsys_t * const ss = ssw->ss;
+  ssw->omega = omega;
+  ssw->medium = qpms_epsmu_generator_eval(ss->medium, omega);
+  ssw->wavenumber = qpms_wavenumber(omega, ssw->medium);
+  qpms_tmatrix_t **tmatrices_preop;
+  QPMS_CRASHING_CALLOC(tmatrices_preop, ss->tmg_count, sizeof(*tmatrices_preop));
+  for (qpms_ss_tmgi_t tmgi = 0; tmgi < ss->tmg_count; ++tmgi) 
+    tmatrices_preop[tmgi] = qpms_tmatrix_init_from_function(ss->tmg[tmgi], omega);
+  for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi) 
+    qpms_tmatrix_apply_operation_replace(ssw->tm[tmi], &ss->tm[tmi].op,
+        tmatrices_preop[ss->tm[tmi].tmgi]);
+  for (qpms_ss_tmgi_t tmgi = 0; tmgi < ss->tmg_count; ++tmgi)
+    qpms_tmatrix_free(tmatrices_preop[tmgi]);
+  free(tmatrices_preop);
+}
 
+qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss,
+    complex double omega) {
+  // TODO
+  qpms_scatsys_at_omega_t *ssw;
+  QPMS_CRASHING_MALLOC(ssw, sizeof(*ssw));
+  ssw->omega = omega;
+  ssw->ss = ss;
+  ssw->medium = qpms_epsmu_generator_eval(ss->medium, omega);
+  ssw->wavenumber = qpms_wavenumber(omega, ssw->medium);
+  QPMS_CRASHING_CALLOC(ssw->tm, ss->tm_count, sizeof(*ssw->tm));
+  qpms_tmatrix_t **tmatrices_preop;
+  QPMS_CRASHING_CALLOC(tmatrices_preop, ss->tmg_count, sizeof(*tmatrices_preop));
+  for (qpms_ss_tmgi_t tmgi = 0; tmgi < ss->tmg_count; ++tmgi) 
+    tmatrices_preop[tmgi] = qpms_tmatrix_init_from_function(ss->tmg[tmgi], omega);
+  for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi) {
+    ssw->tm[tmi] = qpms_tmatrix_apply_operation(&ss->tm[tmi].op,
+        tmatrices_preop[ss->tm[tmi].tmgi]); //<- main difference to .._refill()
+    QPMS_ENSURE(ssw->tm[tmi], 
+        "Got NULL pointer from qpms_tmatrix_apply_operation");
+  }
+  for (qpms_ss_tmgi_t tmgi = 0; tmgi < ss->tmg_count; ++tmgi)
+    qpms_tmatrix_free(tmatrices_preop[tmgi]);
+  free(tmatrices_preop);
+  return ssw;
+}
+
+void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw) {
+  if (ssw) {
+    if(ssw->tm) 
+     for(qpms_ss_tmi_t i = 0; i < ssw->ss->tm_count; ++i) 
+      qpms_tmatrix_free(ssw->tm[i]); 
+    free(ssw->tm);
+  }
+  free(ssw);
+}
 
 // (copypasta from symmetries.c)
 // TODO at some point, maybe support also other norms.
@@ -442,7 +543,7 @@ static inline void check_norm_compat(const qpms_vswf_set_spec_t *s)
     case QPMS_NORMALISATION_NORM_SPHARM:
       break;
     default:
-      abort(); // Only SPHARM and POWER norms are supported right now.
+      QPMS_WTF; // Only SPHARM and POWER norms are supported right now.
   }
 }
 
@@ -455,9 +556,8 @@ complex double *qpms_orbit_action_matrix(complex double *target,
   // check_norm_compat(bspec); not needed here, the qpms_irot3_uvswfi_dense should complain if necessary
   const size_t n = bspec->n;
   const qpms_gmi_t N = ot->size;
-  if (target == NULL)
-    target = malloc(n*n*N*N*sizeof(complex double));
-  if (target == NULL) abort();
+  if (target == NULL) 
+    QPMS_CRASHING_MALLOC(target, n*n*N*N*sizeof(complex double));
   memset(target, 0, n*n*N*N*sizeof(complex double));
   complex double tmp[n][n]; // this is the 'per-particle' action
   qpms_irot3_uvswfi_dense(tmp[0], bspec, sym->rep3d[g]); 
@@ -500,8 +600,7 @@ complex double *qpms_orbit_irrep_projector_matrix(complex double *target,
   const size_t n = bspec->n;
   const qpms_gmi_t N = ot->size;
   if (target == NULL)
-    target = malloc(n*n*N*N*sizeof(complex double));
-  if (target == NULL) abort();
+    QPMS_CRASHING_MALLOC(target, n*n*N*N*sizeof(complex double));
   memset(target, 0, n*n*N*N*sizeof(complex double));
   // Workspace for the orbit group action matrices
   complex double *tmp = malloc(n*n*N*N*sizeof(complex double));
@@ -556,7 +655,6 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size,
   const bool newtarget = (target == NULL);
   if (newtarget)
     QPMS_CRASHING_MALLOC(target,n*n*N*N*sizeof(complex double));
-  if (target == NULL) abort();
   memset(target, 0, n*n*N*N*sizeof(complex double));
 
   // Get the projector (also workspace for right sg. vect.)
@@ -564,12 +662,10 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size,
     ot, bspec, sym, iri);
   if(!projector) abort();
   // Workspace for the right singular vectors.
-  complex double *V_H = malloc(n*n*N*N*sizeof(complex double));
-  if(!V_H) abort();
+  complex double *V_H; QPMS_CRASHING_MALLOC(V_H, n*n*N*N*sizeof(complex double));
   // THIS SHOULD NOT BE NECESSARY
-  complex double *U = malloc(n*n*N*N*sizeof(complex double));
-  if(!U) abort();
-  double *s = malloc(n*N*sizeof(double)); if(!s) abort();
+  complex double *U; QPMS_CRASHING_MALLOC(U, n*n*N*N*sizeof(complex double));
+  double *s; QPMS_CRASHING_MALLOC(s, n*N*sizeof(double));
 
   int info = LAPACKE_zgesdd(LAPACK_ROW_MAJOR,
       'A', // jobz; overwrite projector with left sg.vec. and write right into V_H
@@ -646,8 +742,7 @@ complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_pac
     return target_packed; 
   const size_t full_len = ss->fecv_size;
   if (target_packed == NULL)
-    target_packed = malloc(SQ(packedlen)*sizeof(complex double));
-  if (target_packed == NULL) abort();
+    QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double));
   memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
 
   // Workspace for the intermediate matrix
@@ -684,8 +779,7 @@ complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_f
   const size_t packedlen = ss->saecv_sizes[iri];
   const size_t full_len = ss->fecv_size;
   if (target_full == NULL)
-    target_full = malloc(SQ(full_len)*sizeof(complex double));
-  if (target_full == NULL) abort();
+    QPMS_CRASHING_MALLOC(target_full, SQ(full_len)*sizeof(complex double));
   if(!add) memset(target_full, 0, SQ(full_len)*sizeof(complex double));
 
   if(!packedlen) return target_full; // Empty irrep, do nothing.
@@ -724,13 +818,12 @@ complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
     return target_packed; 
   const size_t full_len = ss->fecv_size;
   if (target_packed == NULL)
-    target_packed = malloc(SQ(packedlen)*sizeof(complex double));
-  if (target_packed == NULL) abort();
+    QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double));
   memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
 
   // Workspace for the intermediate particle-orbit matrix result
-  complex double *tmp = malloc(sizeof(complex double) * SQ(ss->max_bspecn)
-      * ss->sym->order); if (!tmp) abort();
+  complex double *tmp;
+  QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order);
 
   const complex double one = 1, zero = 0;
 
@@ -806,15 +899,14 @@ complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
   const size_t packedlen = ss->saecv_sizes[iri];
   const size_t full_len = ss->fecv_size;
   if (target_full == NULL)
-    target_full = malloc(SQ(full_len)*sizeof(complex double));
-  if (target_full == NULL) abort();
+    QPMS_CRASHING_MALLOC(target_full, SQ(full_len)*sizeof(complex double));
   if(!add) memset(target_full, 0, SQ(full_len)*sizeof(complex double));
 
   if(!packedlen) return target_full; // Empty irrep, do nothing.
 
   // Workspace for the intermediate particle-orbit matrix result
-  complex double *tmp = malloc(sizeof(complex double) * SQ(ss->max_bspecn)
-      * ss->sym->order); if (!tmp) abort();
+  complex double *tmp;
+  QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order);
 
   const complex double one = 1, zero = 0;
 
@@ -889,8 +981,7 @@ complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
   const size_t packedlen = ss->saecv_sizes[iri];
   if (!packedlen) return target_packed; // Empty irrep
   if (target_packed == NULL)
-    target_packed = malloc(packedlen*sizeof(complex double));
-  if (target_packed == NULL) abort();
+    QPMS_CRASHING_MALLOC(target_packed, packedlen*sizeof(complex double));
   memset(target_packed, 0, packedlen*sizeof(complex double));
 
   const complex double one = 1;
@@ -928,8 +1019,7 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
 		const qpms_iri_t iri, bool add) {
   const size_t full_len = ss->fecv_size;
   if (target_full == NULL)
-    target_full = malloc(full_len*sizeof(complex double));
-  if (target_full == NULL) abort();
+    QPMS_CRASHING_MALLOC(target_full, full_len*sizeof(complex double));
   if (!add) memset(target_full, 0, full_len*sizeof(complex double));
 
   const complex double one = 1;
@@ -989,11 +1079,11 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
   { // Non-diagonal part; M[piR, piC] = T[piR] S(piR<-piC)
     size_t fullvec_offsetR = 0;
     for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) {
-      const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[piR].tmatrix_id]->spec;
+      const qpms_vswf_set_spec_t *bspecR = qpms_ss_bspec_pi(ss, piR); 
       const cart3_t posR = ss->p[piR].pos;
       size_t fullvec_offsetC = 0;
       for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) {
-        const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec;
+        const qpms_vswf_set_spec_t *bspecC = qpms_ss_bspec_pi(ss, piC);
         if(piC != piR)  { // The diagonal will be dealt with later.
           const cart3_t posC = ss->p[piC].pos;
           QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
@@ -1003,7 +1093,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
         }
         fullvec_offsetC += bspecC->n;
       }
-      assert(fullvec_offsetC = full_len);
+      assert(fullvec_offsetC == full_len);
       fullvec_offsetR += bspecR->n;
     }
     assert(fullvec_offsetR == full_len);
@@ -1014,13 +1104,14 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
 
 
 
-complex double *qpms_scatsys_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,
-    const qpms_scatsys_t *ss,
-    complex double k ///< Wave number to use in the translation matrix.
+    const qpms_scatsys_at_omega_t *ssw
     )
 {
+  const complex double k = ssw->wavenumber;
+  const qpms_scatsys_t *ss = ssw->ss;
   const size_t full_len = ss->fecv_size;
   if(!target)
     QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
@@ -1031,13 +1122,13 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full(
   { // Non-diagonal part; M[piR, piC] = -T[piR] S(piR<-piC)
     size_t fullvec_offsetR = 0;
     for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) {
-      const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[piR].tmatrix_id]->spec;
+      const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[piR].tmatrix_id]->spec;
       const cart3_t posR = ss->p[piR].pos;
       size_t fullvec_offsetC = 0;
       // dest particle T-matrix
-      const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m;
+      const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m;
       for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) {
-        const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec;
+        const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec;
         if(piC != piR) { // The diagonal will be dealt with later.
           const cart3_t posC = ss->p[piC].pos;
           QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
@@ -1064,19 +1155,20 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full(
 }
 
 // Serial reference implementation.
-complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(
+complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
     /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
     complex double *target_packed,
-    const qpms_scatsys_t *ss, qpms_iri_t iri,
-    complex double k ///< Wave number to use in the translation matrix.
+    const qpms_scatsys_at_omega_t *ssw,
+    qpms_iri_t iri
     )
 {
+  const qpms_scatsys_t *ss = ssw->ss;
+  const complex double k = ssw->wavenumber;
   const size_t packedlen = ss->saecv_sizes[iri];
   if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
     return target_packed; 
   if (target_packed == NULL)
-    target_packed = malloc(SQ(packedlen)*sizeof(complex double));
-  if (target_packed == NULL) abort();
+    QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double));
   memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
 
   // some of the following workspaces are probably redundant; TODO optimize later.
@@ -1102,7 +1194,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(
     const size_t packed_orbit_offsetR =
       ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiR] 
       + osnR * otR->irbase_sizes[iri];
-    const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[piR].tmatrix_id]->spec;
+    const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[piR].tmatrix_id]->spec;
     // Orbit coeff vector's full size:
     const size_t orbit_fullsizeR = otR->size * otR->bspecn;
     const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n
@@ -1112,7 +1204,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(
     const cart3_t posR = ss->p[piR].pos;
     if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep
       // dest particle T-matrix
-      const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m;
+      const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m;
       for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { //Column loop
         const qpms_ss_oti_t otiC = ss->p_orbitinfo[piC].t;
         const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC;
@@ -1122,7 +1214,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(
         const size_t packed_orbit_offsetC = 
           ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC]
           + osnC * otC->irbase_sizes[iri];
-        const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec;
+        const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec;
         // Orbit coeff vector's full size:
         const size_t orbit_fullsizeC = otC->size * otC->bspecn;
         const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n
@@ -1174,19 +1266,19 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(
   return target_packed;
 }
 
-complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
+complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
     /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
     complex double *target_packed,
-    const qpms_scatsys_t *ss, qpms_iri_t iri,
-    complex double k ///< Wave number to use in the translation matrix.
+    const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri
     )
 {
+  const qpms_scatsys_t *ss = ssw->ss;
+  const complex double k = ssw->wavenumber;
   const size_t packedlen = ss->saecv_sizes[iri];
   if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
     return target_packed; 
   if (target_packed == NULL)
-    target_packed = malloc(SQ(packedlen)*sizeof(complex double));
-  if (target_packed == NULL) abort();
+    QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double));
   memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
 
   // some of the following workspaces are probably redundant; TODO optimize later.
@@ -1215,7 +1307,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
 
     if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep
       const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n
-      const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[orbitstartpiR].tmatrix_id]->spec;
+      const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[orbitstartpiR].tmatrix_id]->spec;
       // This is the orbit-level matrix projecting the whole orbit onto the irrep.
       const complex double *omR = otR->irbases + otR->irbase_offsets[iri];
       // Orbit coeff vector's full size:
@@ -1231,7 +1323,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
         assert(ss->p_orbitinfo[piR].osn == osnR);
         const cart3_t posR = ss->p[piR].pos;
         // dest particle T-matrix
-        const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m;
+        const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m;
         for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { //Column loop
           const qpms_ss_oti_t otiC = ss->p_orbitinfo[piC].t;
           const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC;
@@ -1241,7 +1333,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
           const size_t packed_orbit_offsetC = 
             ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC]
             + osnC * otC->irbase_sizes[iri];
-          const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec;
+          const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec;
           // Orbit coeff vector's full size:
           const size_t orbit_fullsizeC = otC->size * otC->bspecn;
           const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n
@@ -1294,20 +1386,21 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
   return target_packed;
 }
 
-struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{
-  const qpms_scatsys_t *ss;
+struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{
+  const qpms_scatsys_at_omega_t *ssw;
   qpms_ss_pi_t *opistartR_ptr;
   pthread_mutex_t *opistartR_mutex;
   qpms_iri_t iri;
   complex double *target_packed;
-  complex double k;
 };
 
-static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread(void *arg)
+static void *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread(void *arg)
 {
-  const struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg 
+  const struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg 
     *a = arg;
-  const qpms_scatsys_t *ss = a->ss;
+  const qpms_scatsys_at_omega_t *ssw = a->ssw;
+  const complex double k = ssw->wavenumber;
+  const qpms_scatsys_t *ss = ssw->ss;
   const qpms_iri_t iri = a->iri;
   const size_t packedlen = ss->saecv_sizes[iri];
 
@@ -1347,7 +1440,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread
 
     if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep
       const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n
-      const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[orbitstartpiR].tmatrix_id]->spec;
+      const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[orbitstartpiR].tmatrix_id]->spec;
       // This is the orbit-level matrix projecting the whole orbit onto the irrep.
       const complex double *omR = otR->irbases + otR->irbase_offsets[iri];
       // Orbit coeff vector's full size:
@@ -1363,7 +1456,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread
         assert(ss->p_orbitinfo[piR].osn == osnR);
         const cart3_t posR = ss->p[piR].pos;
         // dest particle T-matrix
-        const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m;
+        const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m;
         for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { //Column loop
           const qpms_ss_oti_t otiC = ss->p_orbitinfo[piC].t;
           const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC;
@@ -1373,7 +1466,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread
           const size_t packed_orbit_offsetC = 
             ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC]
             + osnC * otC->irbase_sizes[iri];
-          const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec;
+          const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec;
           // Orbit coeff vector's full size:
           const size_t orbit_fullsizeC = otC->size * otC->bspecn;
           const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n
@@ -1387,7 +1480,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread
               QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
                     Sblock, // Sblock is S(piR->piC)
                     bspecR, bspecC->n, bspecC, 1,
-                    a->k, posR, posC, QPMS_HANKEL_PLUS));
+                    k, posR, posC, QPMS_HANKEL_PLUS));
 
               SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                   bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/,
@@ -1484,7 +1577,7 @@ static void *qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thre
 
     if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep
       const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n
-      const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[orbitstartpiR].tmatrix_id]->spec;
+      const qpms_vswf_set_spec_t *bspecR = qpms_ss_bspec_pi(ss, orbitstartpiR);
       // This is the orbit-level matrix projecting the whole orbit onto the irrep.
       const complex double *omR = otR->irbases + otR->irbase_offsets[iri];
       // Orbit coeff vector's full size:
@@ -1508,7 +1601,7 @@ static void *qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thre
           const size_t packed_orbit_offsetC = 
             ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC]
             + osnC * otC->irbase_sizes[iri];
-          const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec;
+          const qpms_vswf_set_spec_t *bspecC = qpms_ss_bspec_pi(ss, piC);
           // Orbit coeff vector's full size:
           const size_t orbit_fullsizeC = otC->size * otC->bspecn;
           const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n
@@ -1567,7 +1660,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
 		complex double *target_packed,
 		const qpms_scatsys_t *ss,
 		qpms_iri_t iri,
-		complex double k, ///< Wave number to use in the translation matrix.
+    const complex double k,
 		qpms_bessel_t J
 		)
 {
@@ -1576,8 +1669,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
   if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
     return target_packed; 
   if (target_packed == NULL)
-    target_packed = malloc(SQ(packedlen)*sizeof(complex double));
-  if (target_packed == NULL) abort();
+    QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double));
   memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
 
   qpms_ss_pi_t opistartR = 0;
@@ -1618,26 +1710,24 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
 
 
 // Parallel implementation, now default
-complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed(
+complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
     /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
     complex double *target_packed,
-    const qpms_scatsys_t *ss, qpms_iri_t iri,
-    complex double k ///< Wave number to use in the translation matrix.
+    const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri
     )
 {
-  const size_t packedlen = ss->saecv_sizes[iri];
+  const size_t packedlen = ssw->ss->saecv_sizes[iri];
   if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
     return target_packed; 
   if (target_packed == NULL)
-    target_packed = malloc(SQ(packedlen)*sizeof(complex double));
-  if (target_packed == NULL) abort();
+    QPMS_CRASHING_MALLOC(target_packed,SQ(packedlen)*sizeof(complex double));
   memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
 
   qpms_ss_pi_t opistartR = 0;
   pthread_mutex_t opistartR_mutex;
   QPMS_ENSURE_SUCCESS(pthread_mutex_init(&opistartR_mutex, NULL));
-  const struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg
-    arg = {ss, &opistartR, &opistartR_mutex, iri, target_packed, k};
+  const struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg
+    arg = {ssw, &opistartR, &opistartR_mutex, iri, target_packed};
 
   // FIXME THIS IS NOT PORTABLE:
   long nthreads;
@@ -1658,7 +1748,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed(
   pthread_t thread_ids[nthreads];
   for(long thi = 0; thi < nthreads; ++thi)
     QPMS_ENSURE_SUCCESS(pthread_create(thread_ids + thi, NULL,
-      qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread,
+      qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread,
       (void *) &arg));
   for(long thi = 0; thi < nthreads; ++thi) {
     void *retval;
@@ -1696,18 +1786,18 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
 #endif
 
 
-complex double *qpms_scatsys_apply_Tmatrices_full(
+complex double *qpms_scatsysw_apply_Tmatrices_full(
 		complex double *target_full, const complex double *inc_full, 
-		const qpms_scatsys_t *ss) {
+		const qpms_scatsys_at_omega_t *ssw) {
   QPMS_UNTESTED;
+  const qpms_scatsys_t *ss = ssw->ss;
   if (!target_full) QPMS_CRASHING_CALLOC(target_full, ss->fecv_size,
       sizeof(complex double));
   for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
     complex double *ptarget = target_full + ss->fecv_pstarts[pi];
     const complex double *psrc = inc_full + ss->fecv_pstarts[pi];
-    const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi);
     // TODO check whether T-matrix is non-virtual after virtual t-matrices are implemented.
-    const qpms_tmatrix_t *T = ss->tm[ss->p[pi].tmatrix_id];
+    const qpms_tmatrix_t *T = ssw->tm[ss->p[pi].tmatrix_id];
     qpms_apply_tmatrix(ptarget, psrc, T);
   }
   return target_full;
@@ -1751,8 +1841,9 @@ void qpms_ss_LU_free(qpms_ss_LU lu) {
   free(lu.ipiv);
 }
 
-qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(complex double *mpmatrix_full,
-    int *target_piv, const qpms_scatsys_t *ss) {
+qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(complex double *mpmatrix_full,
+    int *target_piv, const qpms_scatsys_at_omega_t *ssw) {
+  const qpms_scatsys_t *ss = ssw->ss;
   QPMS_ENSURE(mpmatrix_full, "A non-NULL pointer to the pre-calculated mode matrix is required");
   if (!target_piv) QPMS_CRASHING_MALLOC(target_piv, ss->fecv_size * sizeof(int));
   QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, ss->fecv_size, ss->fecv_size,
@@ -1760,45 +1851,45 @@ qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(complex double *mpmatr
   qpms_ss_LU lu;
   lu.a = mpmatrix_full;
   lu.ipiv = target_piv;
-  lu.ss = ss;
+  lu.ssw = ssw;
   lu.full = true;
   lu.iri = -1;
   return lu;
 }
 
-qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(complex double *mpmatrix_packed,
-    int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri) {
+qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(complex double *mpmatrix_packed,
+    int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) {
   QPMS_ENSURE(mpmatrix_packed, "A non-NULL pointer to the pre-calculated mode matrix is required");
-  size_t n = ss->saecv_sizes[iri];
+  size_t n = ssw->ss->saecv_sizes[iri];
   if (!target_piv) QPMS_CRASHING_MALLOC(target_piv, n * sizeof(int));
   QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n,
         mpmatrix_packed, n, target_piv));
   qpms_ss_LU lu;
   lu.a = mpmatrix_packed;
   lu.ipiv = target_piv;
-  lu.ss = ss;
+  lu.ssw = ssw;
   lu.full = false;
   lu.iri = iri;
   return lu;
 }
 
-qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU(
+qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(
     complex double *target, int *target_piv,
-    const qpms_scatsys_t *ss, complex double k){
-  target = qpms_scatsys_build_modeproblem_matrix_full(target, ss, k);
-  return qpms_scatsys_modeproblem_matrix_full_factorise(target, target_piv, ss);
+    const qpms_scatsys_at_omega_t *ssw){
+  target = qpms_scatsysw_build_modeproblem_matrix_full(target, ssw);
+  return qpms_scatsysw_modeproblem_matrix_full_factorise(target, target_piv, ssw);
 }
 
-qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU(
+qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
     complex double *target, int *target_piv,
-    const qpms_scatsys_t *ss, qpms_iri_t iri, complex double k){
-  target = qpms_scatsys_build_modeproblem_matrix_irrep_packed(target, ss, iri, k);
-  return qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(target, target_piv, ss, iri);
+    const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri){
+  target = qpms_scatsysw_build_modeproblem_matrix_irrep_packed(target, ssw, iri);
+  return qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(target, target_piv, ssw, iri);
 }
 
 complex double *qpms_scatsys_scatter_solve(
     complex double *f, const complex double *a_inc, qpms_ss_LU lu) {
-  const size_t n = lu.full ? lu.ss->fecv_size : lu.ss->saecv_sizes[lu.iri];
+  const size_t n = lu.full ? lu.ssw->ss->fecv_size : lu.ssw->ss->saecv_sizes[lu.iri];
   if (!f) QPMS_CRASHING_MALLOC(f, n * sizeof(complex double));
   memcpy(f, a_inc, n*sizeof(complex double)); // It will be rewritten by zgetrs
   QPMS_ENSURE_SUCCESS(LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N' /*trans*/,  n /*n*/, 1 /*nrhs number of right hand sides*/,
diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h
index 16578a6..49bc9b8 100644
--- a/qpms/scatsystem.h
+++ b/qpms/scatsystem.h
@@ -12,6 +12,7 @@
 #define QPMS_SCATSYSTEM_H
 #include "qpms_types.h"
 #include "vswf.h"
+#include "tmatrices.h"
 #include <stdbool.h>
 
 /// Overrides the number of threads spawned by the paralellized functions.
@@ -19,10 +20,13 @@
 void qpms_scatsystem_set_nthreads(long n);
 
 /// A particle, defined by its T-matrix and position.
+/** This is rather only an auxillary intermediate structure to ultimately 
+ * build an qpms_scatsys_t instance */
 typedef struct qpms_particle_t {
 	// Does it make sense to ever use other than cartesian coords for this?
 	cart3_t pos; ///< Particle position in cartesian coordinates.
-	const qpms_tmatrix_t *tmatrix; ///< T-matrix; not owned by qpms_particle_t.
+	const qpms_tmatrix_function_t *tmg; ///< T-matrix function; not owned by qpms_particle_t.
+	qpms_tmatrix_operation_t op; ///< T-matrix transformation operation w.r.t. \a tmg.
 } qpms_particle_t;
 
 struct qpms_finite_group_t;
@@ -121,13 +125,32 @@ typedef struct qpms_ss_particle_orbitinfo {
 	qpms_ss_orbit_pi_t p; ///< Order (sija, ei rankki) of the particle inside that orbit type.
 } qpms_ss_particle_orbitinfo_t;
 
+/// Auxillary type used in qpms_scatsys_t: A recepy to create another T-matrices by symmetry operations.
+typedef struct qpms_ss_derived_tmatrix_t {
+	qpms_ss_tmgi_t tmgi; ///< Index of the corresponding qpms_scatsys_t::tm element.
+	struct qpms_tmatrix_operation_t op; ///< Operation to derive this particular T-matrix.
+} qpms_ss_derived_tmatrix_t;
+
+
 struct qpms_trans_calculator;
+struct qpms_epsmu_generator_t;
 
 typedef struct qpms_scatsys_t {
-	// TODO does bspec belong here?
-	qpms_tmatrix_t **tm; ///< T-matrices in the system
-	qpms_ss_tmi_t tm_count; ///< Number of all different T-matrices
+	struct qpms_epsmu_generator_t medium; ///< Optical properties of the background medium.
+	
+	/// (Template) T-matrix functions in the system.
+	/** The qpms_abstract_tmatrix_t objects (onto which this array member point)
+	 *  are NOT owned by this and must be kept alive for the whole lifetime
+	 *  of all qpms_scatsys_t objects that are built upon them.
+	 */
+	qpms_tmatrix_function_t *tmg; 
+	qpms_ss_tmgi_t tmg_count; ///< Number of all different original T-matrix generators in the system.
+
+	/// All the different T-matrix functions in the system, including those derived from \a tmg elements by symmetries.
+	qpms_ss_derived_tmatrix_t *tm; 
+	qpms_ss_tmi_t tm_count; ///< Number of all different T-matrices in the system (length of tm[]).
 	qpms_ss_tmi_t tm_capacity; ///< Capacity of tm[].
+	
 	qpms_particle_tid_t *p; ///< Particles.
 	qpms_ss_pi_t p_count; ///< Size of particles array.
 	qpms_ss_pi_t p_capacity; ///< Capacity of p[].
@@ -167,23 +190,77 @@ typedef struct qpms_scatsys_t {
 	struct qpms_trans_calculator *c;
 } qpms_scatsys_t;
 
-/// Convenience function to access pi'th particle's bspec.
-static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) {
-	return ss->tm[ss->p[pi].tmatrix_id]->spec;
+/// Retrieve the bspec of \a tmi'th element of \a ss->tm.
+static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(const qpms_scatsys_t *ss, qpms_ss_tmi_t tmi) {
+	return ss->tmg[ss->tm[tmi].tmgi].spec;
 }
 
-/// Creates a new scatsys by applying a symmetry group, copying particles if needed.
-/** In fact, it copies everything except the vswf set specs, so keep them alive until scatsys is destroyed.
- * The following fields must be filled:
- * orig->tm
- * orig->tm_count
- * orig->p
- * orig->p_count
+/// Retrieve the bspec of \a pi'th particle in \a ss->p.
+static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) {
+	return ss->tmg[ss->tm[ss->p[pi].tmatrix_id].tmgi].spec;
+}
+
+typedef struct qpms_scatsys_at_omega_t {
+	const qpms_scatsys_t *ss; ///< Parent scattering system.
+	/// T-matrices from \a ss, evaluated at \a omega.
+	/** The T-matrices are in the same order as in \a ss,
+	 * i.e in the order corresponding to \a ss->tm.
+	 */
+	qpms_tmatrix_t **tm;
+	complex double omega; ///< Angular frequency
+	qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency
+	complex double wavenumber; ///< Background medium wavenumber
+} qpms_scatsys_at_omega_t;
+
+/// Creates a new scatsys by applying a symmetry group onto a "proto-scatsys", copying particles if needed.
+/** In fact, it copies everything except the vswf set specs and qpms_abstract_tmatrix_t instances,
+ * so keep them alive until scatsys is destroyed.
+ *  
+ *  The following fields must be filled in the "proto- scattering system" \a orig:
+ *  * orig->medium – The pointers are copied to the new qpms_scatsys_t instance; 
+ *    the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting 
+ *    qpms_scatsys_t instances are properly destroyed.
+ *  * orig->tmg – The pointers are copied to the new qpms_scatsys_t instance; 
+ *    the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting 
+ *    qpms_scatsys_t instances are properly destroyed. The pointers from orig->tmg, however, are copied.
+ *  * orig->tmg_count
+ *  * orig->tm – Must be filled, although the operations will typically be identities
+ *    (QPMS_TMATRIX_OPERATION_NOOP). N.B. these NOOPs might be replaced with some symmetrisation operation
+ *    in the resulting "full" qpms_scatsys_t instance.
+ *  * orig->tm_count
+ *  * orig->p
+ *  * orig->p_count
+ *
+ *  The resulting qpms_scatsys_t is obtained by actually evaluating the T-matrices
+ *  at the given frequency \a omega and where applicable, these are compared
+ *  by their values with given tolerances. The T-matrix generators are expected
+ *  to preserve the point group symmetries for all frequencies.
+ *
+ *  This particular function will likely change in the future.
+ *
+ *  \returns An instance \a sso of qpms_scatsys_omega_t. Note that \a sso->ss
+ *  must be saved by the caller before destroying \a sso 
+ *  (with qpms_scatsys_at_omega_free(), and destroyed only afterwards with 
+ *  qpms_scatsys_free() when not needed anymore.
+ *  \a sso->ss can be reused for different frequency by a
+ *  qpms_scatsys_at_omega() call.
+ *
  */
-qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym);
+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);
+
 /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load.
 void qpms_scatsys_free(qpms_scatsys_t *s);
 
+/// Destroys a qpms_scatsys_at_omega_t.
+/** Used on results of qpms_scatsys_apply_symmetry() and qpms_scatsys_at_omega(). */
+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);
+
 /// 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.
  * 
@@ -266,37 +343,36 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
 		);
 
 /// Creates the full \f$ (I - TS) \f$ matrix of the scattering system.
-complex double *qpms_scatsys_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,
-		const qpms_scatsys_t *ss,
-		complex double k ///< Wave number to use in the translation matrix.
+		const qpms_scatsys_at_omega_t *ssw
 		);
 /// Creates the mode problem matrix \f$ (I - TS) \f$ directly in the irrep-packed form.
-complex double *qpms_scatsys_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,
-		const qpms_scatsys_t *ss, qpms_iri_t iri,
-		complex double k ///< Wave number to use in the translation matrix.
+		const qpms_scatsys_at_omega_t *ssw,
+		qpms_iri_t iri
 		);
-/// Alternative implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed().
-complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
+/// Alternative implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed().
+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,
-		const qpms_scatsys_t *ss, qpms_iri_t iri,
-		complex double k ///< Wave number to use in the translation matrix.
+		const qpms_scatsys_at_omega_t *ssw,
+		qpms_iri_t iri
 		);
-/// Alternative (serial reference) implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed().
-complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorder_serial(
+/// Alternative (serial reference) implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed().
+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,
-		const qpms_scatsys_t *ss, qpms_iri_t iri,
-		complex double k ///< Wave number to use in the translation matrix.
+		const qpms_scatsys_at_omega_t *ssw,
+		qpms_iri_t iri
 		);
 
 /// LU factorisation (LAPACKE_zgetrf) result holder.
 typedef struct qpms_ss_LU {
-	const qpms_scatsys_t *ss;
+	const qpms_scatsys_at_omega_t *ssw;
 	bool full; ///< true if full matrix; false if irrep-packed.
 	qpms_iri_t iri; ///< Irrep index if `full == false`.
 	/// LU decomposition array.
@@ -307,33 +383,33 @@ typedef struct qpms_ss_LU {
 void qpms_ss_LU_free(qpms_ss_LU);
 
 /// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch.
-qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU(
+qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(
 		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_t *ss,
-		complex double k ///< Wave number to use in the translation matrix.
+		const qpms_scatsys_at_omega_t *ssw
 		);
 
 /// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch.
-qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU(
+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).
 		int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
-		const qpms_scatsys_t *ss, qpms_iri_t iri,
-		complex double k ///< Wave number to use in the translation matrix.
+		const qpms_scatsys_at_omega_t *ssw,
+		qpms_iri_t iri
 		);
 
 /// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents.
-qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(
+qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(
 		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_t *ss
+		const qpms_scatsys_at_omega_t *ssw
 		);
 
 /// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents.
-qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(
+qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(
 		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_t *ss, qpms_iri_t iri
+		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$.
@@ -422,10 +498,10 @@ complex double *qpms_scatsys_incident_field_vector_full(
 		);
 
 /// Applies T-matrices onto an incident field vector in the full basis.
-complex double *qpms_scatsys_apply_Tmatrices_full(
+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_t *ss
+		const qpms_scatsys_at_omega_t *ssw
 		);
 
 
diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c
index 1ce435c..7d21a5c 100644
--- a/qpms/tmatrices.c
+++ b/qpms/tmatrices.c
@@ -53,6 +53,15 @@ qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
   return t;
 }
 
+qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
+                const qpms_vswf_set_spec_t *bspec,
+                qpms_tmatrix_generator_t gen,
+                complex double omega) {
+        qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
+        QPMS_ENSURE_SUCCESS(gen.function(t, omega, gen.params));
+        return t;
+}
+
 qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T) {
   qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
   size_t n = T->spec->n;
@@ -991,3 +1000,224 @@ qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t,
   return QPMS_SUCCESS;
 }
 
+void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) {
+  switch(f->typ) {
+    case QPMS_TMATRIX_OPERATION_NOOP:
+      break;
+    case QPMS_TMATRIX_OPERATION_LRMATRIX:
+      if(f->op.lrmatrix.owns_m)
+        free(f->op.lrmatrix.m);
+      break;
+    case QPMS_TMATRIX_OPERATION_SCMULZ:
+      if(f->op.scmulz.owns_m)
+        free(f->op.scmulz.m);
+      break;
+    case QPMS_TMATRIX_OPERATION_IROT3:
+      break;
+    case QPMS_TMATRIX_OPERATION_IROT3ARR:
+      if(f->op.irot3arr.owns_ops)
+        free(f->op.irot3arr.ops);
+      break;
+    case QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE: // Group never owned
+      break;
+    case QPMS_TMATRIX_OPERATION_COMPOSE_SUM:
+      {
+        struct qpms_tmatrix_operation_compose_sum * const o = 
+          &(f->op.compose_sum);
+        if(o->opmem) {
+          for(size_t i = 0; i < o->n; ++i) 
+            qpms_tmatrix_operation_clear(&(o->opmem[i]));
+          free(o->opmem);
+        }
+        free(o->ops);
+      }
+      break;
+    case QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN:
+      {
+        struct qpms_tmatrix_operation_compose_chain * const o = 
+          &(f->op.compose_chain);
+        if(o->opmem) {
+          for(size_t i = 0; i < o->n; ++i)
+            if(o->ops_owned[i]) {
+              // A bit complicated yet imperfect sanity/bound check,
+              // but this way we at least get rid of the const discard warning.
+              ptrdiff_t d = o->ops[i] - o->opmem;
+              QPMS_ENSURE(d >= 0 && d < o->opmem_size, 
+                  "Bound check failed; is there a bug in the"
+                  " construction of compose_chain operation?\n"
+                  "opmem_size = %zu, o->ops[%zu] - o->opmem = %td.",
+                  o->opmem_size, i, d);
+              qpms_tmatrix_operation_clear(o->opmem + d);
+            }
+          free(o->opmem);
+          free(o->ops_owned);
+        }
+        free(o->ops);
+      }
+      break;
+    default:
+      QPMS_WTF;
+  }
+}
+
+void qpms_tmatrix_operation_copy(qpms_tmatrix_operation_t *dest, const qpms_tmatrix_operation_t *src) {
+  memcpy(dest, src, sizeof(qpms_tmatrix_operation_t));
+  switch(dest->typ) {
+    case QPMS_TMATRIX_OPERATION_NOOP:
+      break;
+    case QPMS_TMATRIX_OPERATION_LRMATRIX:
+      QPMS_CRASHING_MALLOCPY(dest->op.lrmatrix.m, src->op.lrmatrix.m, 
+          sizeof(complex double) * src->op.lrmatrix.m_size);
+      dest->op.lrmatrix.owns_m = true;
+      break;
+    case QPMS_TMATRIX_OPERATION_SCMULZ:
+      QPMS_CRASHING_MALLOCPY(dest->op.scmulz.m, src->op.scmulz.m, 
+          sizeof(complex double) * src->op.scmulz.m_size);
+      dest->op.scmulz.owns_m = true;
+      break;
+    case QPMS_TMATRIX_OPERATION_IROT3:
+      break;
+    case QPMS_TMATRIX_OPERATION_IROT3ARR:
+      QPMS_CRASHING_MALLOCPY(dest->op.irot3arr.ops, src->op.irot3arr.ops,
+          src->op.irot3arr.n * sizeof(qpms_irot3_t));
+      dest->op.irot3arr.owns_ops = true;
+      break;
+    case QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE: // Group never owned
+      break;
+    case QPMS_TMATRIX_OPERATION_COMPOSE_SUM:
+      {
+        struct qpms_tmatrix_operation_compose_sum * const o = 
+          &(dest->op.compose_sum);
+        QPMS_CRASHING_MALLOC(o->ops, o->n * sizeof(*(o->ops)));
+        QPMS_CRASHING_MALLOC(o->opmem, o->n * sizeof(*(o->opmem)));
+        for(size_t i = 0; i < o->n; ++i) {
+          qpms_tmatrix_operation_copy(o->opmem + i, src->op.compose_sum.ops[i]);
+          o->ops[i] = o->opmem + i;
+        }
+      }
+      break;
+    case QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN:
+      {
+        struct qpms_tmatrix_operation_compose_chain * const o = 
+          &(dest->op.compose_chain);
+        QPMS_CRASHING_MALLOC(o->ops, o->n * sizeof(*(o->ops)));
+        QPMS_CRASHING_MALLOC(o->ops_owned, o->n * sizeof(*(o->ops_owned)));
+        QPMS_CRASHING_MALLOC(o->opmem, o->n * sizeof(*(o->opmem)));
+        for(size_t i = 0; i < o->n; ++i) {
+          qpms_tmatrix_operation_copy(o->opmem + i, src->op.compose_chain.ops[i]);
+          o->ops[i] = o->opmem + i;
+          o->ops_owned[i] = true;
+        }
+      }
+      break;
+    default:
+      QPMS_WTF;
+  }
+  dest->typ = src->typ;
+}
+
+void qpms_tmatrix_operation_compose_chain_init(qpms_tmatrix_operation_t *dest,
+    size_t nops, size_t opmem_size) {
+  if (nops == 0) {
+    QPMS_WARN("Tried to create a composed (chain) operation of zero size;"
+        "setting to no-op instead.");
+    *dest = qpms_tmatrix_operation_noop;
+  }
+  if (nops < opmem_size)
+    QPMS_WARN("Allocating buffer for %zu operations, in a chained operation of"
+     " only %zu elemens, that does not seem to make sense.", opmem_size, nops);
+  dest->typ = QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN;
+  struct qpms_tmatrix_operation_compose_chain * const o = 
+          &(dest->op.compose_chain);
+  o->n = nops;
+  QPMS_CRASHING_MALLOC(o->ops, nops * sizeof(*(o->ops)));
+  if (opmem_size != 0) {
+    QPMS_CRASHING_CALLOC(o->ops_owned, o->n, sizeof(_Bool));
+    QPMS_CRASHING_MALLOC(o->opmem, opmem_size * sizeof(*(o->opmem)));
+  } else {
+    o->ops_owned = NULL;
+    o->opmem = NULL;
+  }
+  o->opmem_size = opmem_size;
+}
+
+qpms_tmatrix_t *qpms_tmatrix_mv(qpms_tmatrix_t *dest,
+    const qpms_tmatrix_t *orig) {
+  QPMS_ENSURE(qpms_vswf_set_spec_isidentical(dest->spec, orig->spec),
+     "Basis specifications must be identical!");
+  memcpy(dest->m, orig->m, SQ(orig->spec->n)*sizeof(complex double));
+  return dest;
+}
+
+qpms_tmatrix_t *qpms_tmatrix_apply_operation_replace(qpms_tmatrix_t *dest,
+    const qpms_tmatrix_operation_t *op, const qpms_tmatrix_t *orig) {
+  //QPMS_ENSURE(qpms_vswf_set_spec_isidentical(dest->spec, orig->spec),
+  //    "Basis specifications must be identical!");
+  // TODO should I check also for dest->owns_m?
+  // FIXME this is highly inoptimal; the hierarchy should be such
+  // that _operation() and operation_inplace() call this, and not the
+  // other way around
+  qpms_tmatrix_mv(dest, orig);
+  return qpms_tmatrix_apply_operation_inplace(op, dest);
+}
+
+qpms_tmatrix_t *qpms_tmatrix_apply_operation(
+    const qpms_tmatrix_operation_t *f, const qpms_tmatrix_t *orig) {
+  // Certain operations could be optimized, but the effect would be marginal.
+  qpms_tmatrix_t *res = qpms_tmatrix_copy(orig); 
+  return qpms_tmatrix_apply_operation_inplace(f, res);
+}
+
+static qpms_tmatrix_t *qtao_compose_sum_inplace(qpms_tmatrix_t *T, 
+    const struct qpms_tmatrix_operation_compose_sum *cs) {
+  qpms_tmatrix_t *tmp_target = qpms_tmatrix_init(T->spec);
+  qpms_tmatrix_t *sum = qpms_tmatrix_init(T->spec);
+  for (size_t i = 0; i < cs->n; ++i) {
+    memcpy(tmp_target->m, T->m, SQ(T->spec->n) * sizeof(complex double));
+    QPMS_ENSURE(qpms_tmatrix_apply_operation_inplace(cs->ops[i] , tmp_target),
+        "Got NULL pointer from qpms_tmatrix_apply_operation_inplace, hupsis!");
+    for (size_t j = 0; j < SQ(T->spec->n); ++j)
+      sum->m[j] += tmp_target->m[j];
+  }
+  for(size_t j = 0; j < SQ(T->spec->n); ++j)
+    T->m[j] = sum->m[j] * cs->factor;
+  qpms_tmatrix_free(sum);
+  qpms_tmatrix_free(tmp_target);
+  return T;
+}
+
+static qpms_tmatrix_t *qtao_compose_chain_inplace(qpms_tmatrix_t *T,
+    const struct qpms_tmatrix_operation_compose_chain *cc) {
+  for(size_t i = 0; i < cc->n; ++i)
+    qpms_tmatrix_apply_operation_inplace(cc->ops[i], T);
+  return T;
+}
+
+static qpms_tmatrix_t *qtao_scmulz_inplace(qpms_tmatrix_t *T,
+    const struct qpms_tmatrix_operation_scmulz *s) {
+  for(size_t i = 0; i < SQ(T->spec->n); ++i)
+    T->m[i] *= s->m[i];
+  return T;
+}
+
+qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace(
+    const qpms_tmatrix_operation_t *f, qpms_tmatrix_t *T) {
+  switch(f->typ) {
+    case QPMS_TMATRIX_OPERATION_NOOP:
+      return T;
+    case QPMS_TMATRIX_OPERATION_LRMATRIX:
+      return qpms_tmatrix_apply_symop_inplace(T, f->op.lrmatrix.m);
+    case QPMS_TMATRIX_OPERATION_IROT3:
+      return qpms_tmatrix_symmetrise_irot3arr_inplace(T, 1, &(f->op.irot3));
+    case QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE:
+      return qpms_tmatrix_symmetrise_finite_group_inplace(T, f->op.finite_group);
+    case QPMS_TMATRIX_OPERATION_COMPOSE_SUM:
+      return qtao_compose_sum_inplace(T, &(f->op.compose_sum));
+    case QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN:
+      return qtao_compose_chain_inplace(T, &(f->op.compose_chain));
+    case QPMS_TMATRIX_OPERATION_SCMULZ:
+      return qtao_scmulz_inplace(T, &(f->op.scmulz));
+    default:
+      QPMS_WTF;
+  }
+}
diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h
index c208aa0..117e3f9 100644
--- a/qpms/tmatrices.h
+++ b/qpms/tmatrices.h
@@ -8,6 +8,8 @@
 #include "materials.h"
 #include <stdio.h>
 
+
+
 struct qpms_finite_group_t;
 typedef struct qpms_finite_group_t qpms_finite_group_t;
 
@@ -17,11 +19,21 @@ static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
 }
 
 /// Initialises a zero T-matrix.
+/** \sa qpms_tmatrix_init_from_generator() 
+ *  \sa qpms_tmatrix_init_from_function() */
 qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec);
 
+
 /// Copies a T-matrix, allocating new array for the T-matrix data.
 qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T);
 
+/// Copies a T-matrix contents between two pre-allocated T-matrices.
+/** orig->spec and dest->spec must be identical.
+ *
+ * \returns \a dest 
+ */ 
+qpms_tmatrix_t *qpms_tmatrix_mv(qpms_tmatrix_t *dest, const qpms_tmatrix_t *orig);
+
 /// Destroys a T-matrix.
 void qpms_tmatrix_free(qpms_tmatrix_t *t);
 
@@ -233,7 +245,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
 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
+		const qpms_tmatrix_t *T ///< T-matrix \a T to apply.
 		);
 
 /// Generic T-matrix generator function that fills a pre-initialised qpms_tmatrix_t given a frequency.
@@ -251,6 +263,13 @@ typedef struct qpms_tmatrix_generator_t {
 	const void *params; ///< Parameter pointer passed to the function.
 } qpms_tmatrix_generator_t;
 
+/// Initialises and evaluates a new T-matrix using a generator.
+qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
+		const qpms_vswf_set_spec_t *bspec,
+		qpms_tmatrix_generator_t gen,
+		complex double omega);
+
+
 /// Implementation of qpms_matrix_generator_t that just copies a constant matrix.
 /** N.B. this does almost no checks at all, so use it only for t-matrices with
  * the same base spec.
@@ -427,7 +446,7 @@ qpms_arc_function_retval_t qpms_arc_cylinder(double theta,
 		);
 
 /// Arc parametrisation of spherical particle; for qpms_arc_function_t.
-/** Useful mostly only for benchmarks, as one can use the Mie-Lorentz solution. */
+/** Useful mostly only for benchmarks or debugging, as one can use the Mie-Lorentz solution. */
 qpms_arc_function_retval_t qpms_arc_sphere(double theta,
 		const void *R; ///< Points to double containing particle's radius.
 		);
@@ -503,6 +522,137 @@ qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(complex double *target,
 		);
 
 
+/// An "abstract" T-matrix, contains a T-matrix generator instead of actual data.
+typedef struct qpms_tmatrix_function_t {
+        /** \brief VSWF basis specification, NOT owned by qpms_tmatrix_t by default.
+         *
+         *  Usually not checked for meaningfulness by the functions (methods),
+         *  so the caller should take care that \a spec->ilist does not
+         *  contain any duplicities and that for each wave with order \a m
+         *  there is also one with order \a −m.
+         */
+        const qpms_vswf_set_spec_t *spec;
+        const qpms_tmatrix_generator_t *gen; ///< A T-matrix generator function.
+} 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) {
+	return qpms_tmatrix_init_from_generator(f.spec, *f.gen, omega);
+}
+
+/// Specifies different kinds of operations done on T-matrices
+typedef enum {
+	QPMS_TMATRIX_OPERATION_NOOP, ///< Identity operation.
+	QPMS_TMATRIX_OPERATION_LRMATRIX, ///< General matrix transformation \f$ T' = MTM^\dagger \f$; see @ref qpms_tmatrix_operation_lrmatrix.
+	QPMS_TMATRIX_OPERATION_IROT3, ///< Single rotoreflection specified by a qpms_irot3_t.
+	QPMS_TMATRIX_OPERATION_IROT3ARR, ///< Symmetrise using an array of rotoreflection; see @ref qpms_tmatrix_operation_irot3arr.
+	QPMS_TMATRIX_OPERATION_COMPOSE_SUM, ///< Apply several transformations and sum the results, see @ref qpms_tmatrix_operation_compose_sum.
+	QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN, ///< Chain several different transformations; see @ref qpms_tmatrix_operation_compose_chain.
+	QPMS_TMATRIX_OPERATION_SCMULZ, ///< Elementwise scalar multiplication with a complex matrix; see @ref qpms_tmatrix_operation_scmulz.
+	//QPMS_TMATRIX_OPERATION_POINTGROUP, ///< TODO the new point group in pointgroup.h
+	QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE ///< Legacy qpms_finite_group_t with filled rep3d; see @ref qpms_tmatrix_operation_finite_group.
+} qpms_tmatrix_operation_kind_t;
+
+/// General matrix transformation \f[ T' = MTM^\dagger \f] spec for qpms_tmatrix_operation_t.
+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;
+};
+
+struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations.
+
+/// Specifies a composed operation of type \f$ T' = c\sum_i f_i'(T) \f$ for qpms_tmatrix_operation_t.
+struct qpms_tmatrix_operation_compose_sum {
+	size_t n; ///< Number of operations in ops;
+	const struct qpms_tmatrix_operation_t **ops; ///< Operations array. (Pointers array \a ops[] always owned by this.)
+	double factor; ///< Overall factor \a c.
+	/// (Optional) operations buffer into which elements of \a ops point.
+	/** Owned by this or NULL. If not NULL, all \a ops members are assumed to point into 
+	 * the memory held by \a opmem and to be properly initialised.
+	 * Each \a ops member has to point to _different_ elements of \a opmem.
+	 */ 
+	struct qpms_tmatrix_operation_t *opmem; 
+};
+
+/// Specifies a composed operation of type \f$ T' =  f_{n-1}(f_{n-2}(\dots f_0(T)\dots))) \f$ for qpms_tmatrix_operation_t.
+struct qpms_tmatrix_operation_compose_chain {
+	size_t n; ///< Number of operations in ops;
+	const struct qpms_tmatrix_operation_t **ops; ///< Operations array. (Pointers owned by this.)
+	struct qpms_tmatrix_operation_t *opmem; ///< (Optional) operations buffer into which elements of \a ops point. (Owned by this or NULL.)
+	size_t opmem_size; ///< Length of the opmem array.
+	_Bool *ops_owned; ///< True for all sub operations owned by this and saved in opmem. NULL if opmem is NULL.
+};
+
+/// Specifies an elementwise complex multiplication of type \f$ T'_{ij} = M_{ij}T_{ij} \f$ for qpms_tmatrix_operation_t.
+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.
+};
+
+/// Specifies a symmetrisation using a set of rotoreflections (with equal weights) for qpms_tmatrix_operation_t.
+/** Internally, this is evaluated using a call to qpms_symmetrise_tmdata_irot3arr()
+ * or qpms_symmetrise_tmdata_irot3arr_inplace(). */
+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.
+};
+
+/// A generic T-matrix transformation operator.
+typedef struct qpms_tmatrix_operation_t {
+	/// Specifies the operation kind to be performed and which type \op actually contains.
+	qpms_tmatrix_operation_kind_t typ;
+	union {
+		struct qpms_tmatrix_operation_lrmatrix lrmatrix;
+		struct qpms_tmatrix_operation_scmulz scmulz;
+		qpms_irot3_t irot3; ///< Single rotoreflection; \a typ = QPMS_TMATRIX_OPERATION_IROT3
+		struct qpms_tmatrix_operation_irot3arr irot3arr;
+		struct qpms_tmatrix_operation_compose_sum compose_sum;
+		struct qpms_tmatrix_operation_compose_chain compose_chain;
+		/// Finite group for QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE.
+		/** Not owned by this; \a rep3d must be filled. */
+		const qpms_finite_group_t *finite_group; 
+	} op; ///< Operation data; actual type is determined by \a typ.
+} qpms_tmatrix_operation_t;
+
+static const qpms_tmatrix_operation_t qpms_tmatrix_operation_noop = {.typ = QPMS_TMATRIX_OPERATION_NOOP};
+
+/// Apply an operation on a T-matrix, returning a newly allocated result.
+qpms_tmatrix_t *qpms_tmatrix_apply_operation(const qpms_tmatrix_operation_t *op, const qpms_tmatrix_t *orig);
+
+/// Apply an operation on a T-matrix and replace it with the result.
+qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace(const qpms_tmatrix_operation_t *op, qpms_tmatrix_t *orig);
+
+/// Apply an operation on a T-matrix and replace another one it with the result.
+qpms_tmatrix_t *qpms_tmatrix_apply_operation_replace(qpms_tmatrix_t *dest, 
+		const qpms_tmatrix_operation_t *op, const qpms_tmatrix_t *orig);
+
+/// (Recursively) deallocates internal data of qpms_tmatrix_operation_t.
+/** It does NOT clear the memory pointed to by op it self, but only
+ * heap-allocated data of certain operations, if labeled as owned by it.
+ * In case of compose operations, the recursion stops if the children are
+ * not owned by them (the opmem pointer is NULL).
+ */
+void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f);
+
+/// (Recursively) copies an qpms_tmatrix_operation_t.
+/** Makes copies of all the internal data and takes ownership over them if needed */
+void qpms_tmatrix_operation_copy(qpms_tmatrix_operation_t *target, const qpms_tmatrix_operation_t *src);
+
+/// Inits a new "chain" of composed operations, some of which might be owned.
+void qpms_tmatrix_operation_compose_chain_init(
+		qpms_tmatrix_operation_t *target, ///< The operation structure that will be set to the chain.
+		size_t nops, ///< Number of chained operations (length of the \a ops array)
+		size_t opmem_size ///< Size of the own operations buffer (length of the \a opmem array)
+		);
+
 #if 0
 // Abstract types that describe T-matrix/particle/scatsystem symmetries
 // To be implemented later. See also the thoughts in the beginning of groups.h.
diff --git a/qpms/tolerances.h b/qpms/tolerances.h
new file mode 100644
index 0000000..48e0a13
--- /dev/null
+++ b/qpms/tolerances.h
@@ -0,0 +1,15 @@
+/*! \file tolerances.h */
+#ifndef QPMS_TOLERANCES_H
+#define QPMS_TOLERANCES_H
+
+// TODO DOC
+
+typedef struct qpms_tolerance_spec_t {
+	double atol;
+	double rtol;
+} qpms_tolerance_spec_t;
+
+/// A rather arbitrary default tolerance.
+static const qpms_tolerance_spec_t QPMS_TOLERANCE_DEFAULT = {.atol = 1e-15, .rtol =  1e-8};
+
+#endif // QPMS_TOLERANCES_H