diff --git a/meanas/fdfd/functional.py b/meanas/fdfd/functional.py
index 655d9b8..fe11e51 100644
--- a/meanas/fdfd/functional.py
+++ b/meanas/fdfd/functional.py
@@ -187,3 +187,20 @@ def m2j(omega: complex,
         return m2j_mu
 
 
+def e_tfsf_source(TF_region: field_t,
+                  omega: complex,
+                  dxes: dx_lists_t,
+                  epsilon: field_t,
+                  mu: field_t = None,
+                  ) -> functional_matrix:
+    """
+    Operator that turuns an E-field distribution into a total-field/scattered-field
+    (TFSF) source.
+    """
+    # TODO documentation
+    A = e_full(omega, dxes, epsilon, mu)
+
+    def op(e):
+        neg_iwj = A(TF_region * e) - TF_region * A(e)
+        return neg_iwj / (-1j * omega)
+
diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py
index 774c3d9..3042af4 100644
--- a/meanas/fdfd/operators.py
+++ b/meanas/fdfd/operators.py
@@ -501,3 +501,21 @@ def poynting_h_cross(h: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
          [ bx @ Hz @ fy @ dagx,  zero,                -bz @ Hx @ fy @ dagz],
          [-bx @ Hy @ fz @ dagx,  by @ Hx @ fz @ dagy,  zero]])
     return P
+
+
+def e_tfsf_source(TF_region: vfield_t,
+                  omega: complex,
+                  dxes: dx_lists_t,
+                  epsilon: vfield_t,
+                  mu: vfield_t = None,
+                  ) -> sparse.spmatrix:
+    """
+    Operator that turns an E-field distribution into a total-field/scattered-field
+    (TFSF) source.
+    """
+    # TODO documentation
+    A = e_full(omega, dxes, epsilon, mu)
+    Q = sparse.diags(TF_region)
+    return (A @ Q - Q @ A) / (-1j * omega)
+
+