diff --git a/mbtrack2/impedance/resistive_wall.py b/mbtrack2/impedance/resistive_wall.py
index 7928c481385b2ece101df4dca5b1503595e784f9..31d2d35cce96338cf5af45341c9eac0ee254ed9e 100644
--- a/mbtrack2/impedance/resistive_wall.py
+++ b/mbtrack2/impedance/resistive_wall.py
@@ -5,9 +5,9 @@ Define resistive wall elements based on the WakeField class.
 
 import numpy as np
 from scipy.constants import c, epsilon_0, mu_0
-from scipy.special import wofz as _scipy_wofz
 
 from mbtrack2.impedance.wakefield import Impedance, WakeField, WakeFunction
+from mbtrack2.tracking.particles_electromagnetic_fields import _wofz
 
 
 def skin_depth(frequency, rho, mu_r=1, epsilon_r=1):
@@ -238,23 +238,13 @@ class CircularResistiveWall(WakeField):
             idx2 = np.logical_not(idx1)
             wt[idx2] = self.__TransWakeApprox(time[idx2])
         return wt
-    
-    def __wofz(self, z):
-        """
-        Compute the Faddeeva function w(z) = exp(-z**2) * erfc(-i*z).
-
-        Returns
-        -------
-        tuple
-            Real and imaginary parts of the Faddeeva function.
-        """
-        res = _scipy_wofz(z)
-        return res.real, res.imag
 
     def __LongWakeExact(self, t, factor):
-        w1re, _ = self.__wofz( 1j * np.sqrt(2 * t / self.t0) )
-        w2re, _ = self.__wofz( np.exp(1j * np.pi / 6) *
-                                np.sqrt(2 * t / self.t0) )
+        w1re, _ = _wofz( 0, np.sqrt(2 * t / self.t0) )
+        w2re, _ = _wofz( np.cos(np.pi/6) *
+                         np.sqrt(2 * t / self.t0),
+                         np.sin(np.pi/6) *
+                         np.sqrt(2 * t / self.t0) )
 
         wl = factor * ( 4 * np.exp(-1 * t / self.t0) *
              np.cos(np.sqrt(3) * t / self.t0)
@@ -262,9 +252,11 @@ class CircularResistiveWall(WakeField):
         return wl
 
     def __TransWakeExact(self, t, factor):
-        w1re, _ = self.__wofz( 1j * np.sqrt(2 * t / self.t0) )
-        w2re, w2im = self.__wofz( np.exp(1j * np.pi / 6) *
-                                    np.sqrt(2 * t / self.t0) )
+        w1re, _ = _wofz( 0, np.sqrt(2 * t / self.t0) )
+        w2re, w2im = _wofz( np.cos(np.pi/6) *
+                            np.sqrt(2 * t / self.t0),
+                            np.sin(np.pi/6) *
+                            np.sqrt(2 * t / self.t0) )
 
         wt = factor * ( 2 * np.exp(-1 * t / self.t0) *
              ( np.sqrt(3) * np.sin(np.sqrt(3) * t / self.t0)