diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c
index ef29643..1af486f 100644
--- a/qpms/scatsystem.c
+++ b/qpms/scatsystem.c
@@ -476,7 +476,58 @@ 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_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.
@@ -1050,7 +1101,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
 
 
 
-complex double *qpms_scatsys_at_omega_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_at_omega_t *ssw
diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c
index 96516ba..c60f146 100644
--- a/qpms/tmatrices.c
+++ b/qpms/tmatrices.c
@@ -1140,7 +1140,25 @@ void qpms_tmatrix_operation_compose_chain_init(qpms_tmatrix_operation_t *dest,
   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) {
diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h
index 8b4d2e6..117e3f9 100644
--- a/qpms/tmatrices.h
+++ b/qpms/tmatrices.h
@@ -27,6 +27,13 @@ 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);
 
@@ -623,6 +630,10 @@ qpms_tmatrix_t *qpms_tmatrix_apply_operation(const qpms_tmatrix_operation_t *op,
 /// 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.