Faster CircularResistiveWall

Update CircularResistiveWall class in resistive_wall.py to replace numerical integrations with an analytical forms for improved performance.

• Numerical integrations are removed.
• LongWakeExact and TransWakeExact are used even in the case where time > 20*t0.
• "atol" for LongWakeExact is removed.

I believe "atol" is not needed when considering the fundamental theorem of beam loading. This is because Wl(0) = Wl(0+)/2 is comes from the equation Wl(0) = [Wl(0-)*l/2 + Wl(0+)*l/2] / l, which represents the average Wl value of a very narrow bin (as the length of bin l -> 0) at the origin. Defining "atol" for the fundamental theorem of beam loading introduces the possibility that more than one point could be divided by 2, which I believe is mathematically and physically unacceptable. Thus, it is correct to consider the fundamental theorem of beam loading only at the origin. PyHEADTAIL also considers the fundamental theorem of beam loading only at the origin. Please see def function_longitudinal at the link (line 341).

mbtrack2/impedance/resistive_wall.py

@@ -62,10 +62,6 @@ class CircularResistiveWall(WakeField):
     exact : bool, optional
         If False, approxmiated formulas are used for the wake function 
         computations.
-    atol : float, optional
-        Absolute tolerance used to enforce fundamental theorem of beam loading
-        for the exact expression of the longitudinal wake function.
-        Default is 1e-20.
     
     References
     ----------
@@ -88,8 +84,7 @@ class CircularResistiveWall(WakeField):
                  length,
                  rho,
                  radius,
-                 exact=False,
-                 atol=1e-20):
+                 exact=False):
         super().__init__()
 
         self.length = length
@@ -104,7 +99,7 @@ class CircularResistiveWall(WakeField):
         Z2 = c / omega * length * (1 + np.sign(frequency) * 1j) * rho / (
             np.pi * radius**3 * skin_depth(frequency, rho))
 
-        Wl = self.LongitudinalWakeFunction(time, exact, atol)
+        Wl = self.LongitudinalWakeFunction(time, exact)
         Wt = self.TransverseWakeFunction(time, exact)
 
         Zlong = Impedance(variable=frequency,
@@ -127,7 +122,7 @@ class CircularResistiveWall(WakeField):
         super().append_to_model(Wxdip)
         super().append_to_model(Wydip)
 
-    def LongitudinalWakeFunction(self, time, exact=False, atol=1e-20):
+    def LongitudinalWakeFunction(self, time, exact=False):
         """
         Compute the longitudinal wake function of a circular resistive wall 
         using Eq. (11), of [1], or approxmiated expression Eq. (24), of [2].
@@ -136,8 +131,8 @@ class CircularResistiveWall(WakeField):
         
         Eq. (11) in [1] is completely equivalent to Eq. (22) in [2].
         
-        If some time value is smaller than atol, then the fundamental theorem 
-        of beam loading is applied: Wl(0) = Wl(0+)/2.
+        The fundamental theorem  of beam loading is applied for the exact
+        expression of the longitudinal wake function: Wl(0) = Wl(0+)/2.
 
         Parameters
         ----------
@@ -145,10 +140,6 @@ class CircularResistiveWall(WakeField):
             Time points where the wake function is evaluated in [s].
         exact : bool, optional
             If True, the exact expression is used. The default is False.
-        atol : float, optional
-            Absolute tolerance used to enforce fundamental theorem of beam loading
-            for the exact expression of the longitudinal wake function.
-            Default is 1e-20.
 
         Returns
         -------
@@ -171,17 +162,11 @@ class CircularResistiveWall(WakeField):
         if exact == True:
             idx2 = time == 0
             idx3 = np.logical_not(np.logical_or(idx1, idx2))
-            # idx4 = np.isclose(time, 0, atol=atol) & (time >= 0)
             factor = self.Z0 * c / (3 * np.pi * self.radius**2) * self.length
             if np.any(idx2):
                 # fundamental theorem of beam loading
                 wl[idx2] = 3 * factor / 2
             wl[idx3] = self.__LongWakeExact(time[idx3], factor)
-            # if np.any(idx4):
-            #     closest_to_zero_idx = np.argmin(time[idx4])
-            #     # Get the actual index in the original array
-            #     closest_to_zero_idx = np.where(idx4)[0][closest_to_zero_idx]
-            #     wl[closest_to_zero_idx] *= 0.5
         else:
             idx2 = np.logical_not(idx1)
             wl[idx2] = self.__LongWakeApprox(time[idx2])