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

@@ -133,6 +133,11 @@ class CircularResistiveWall(WakeField):
         to the characteristic time t0.
         
         Eq. (11) in [1] is completely identical to Eq. (22) in [2].
+
+        The real parts of the last two terms of Eq. (11) in [1] are the same,
+        and the imaginary parts have the same magnitude but opposite signs.
+        Therefore, the former term was doubled, the latter term was eliminated,
+        and only the real part was taken to speed up the calculation.
         
         The fundamental theorem of beam loading [3] is applied for the exact
         expression of the longitudinal wake function: Wl(0) = Wl(0+)/2.
@@ -188,10 +193,15 @@ class CircularResistiveWall(WakeField):
         Eq. (11) in [1] is completely identical to Eq. (25) in [2].
 
         There are typos in both Eq. (11) in [1] and Eq. (25) in [2].
-        Corrected the typos in the last two terms of exact expression
-        for transverse wake function in Eq. (11), of [1].
+        Corrected the typos in the last two terms of exact expression for
+        transverse wake function in Eq. (11), of [1].
         It is necessary to multiply Eq. (25) in [2] by c*t0.
 
+        The real parts of the last two terms of Eq. (11) in [1] are the same,
+        and the imaginary parts have the same magnitude but opposite signs.
+        Therefore, the former term was doubled, the latter term was eliminated,
+        and only the real part was taken to speed up the calculation.
+
         Parameters
         ----------
         time : array of float
@@ -230,31 +240,23 @@ class CircularResistiveWall(WakeField):
         return wt
 
     def __LongWakeExact(self, t, factor):
-        wl = np.real( factor *
-                     ( 4 * np.exp(-1 * t / self.t0) *
-                      np.cos(np.sqrt(3) * t / self.t0)
-                      + wofz( 1j *
-                             np.sqrt(2 * t / self.t0) )
-                      - wofz( np.exp(1j * np.pi / 6) *
-                              np.sqrt(2 * t / self.t0) )
-                      - wofz( -1 * np.exp(-1j * 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)
+             + np.real( wofz( 1j *
+                        np.sqrt(2 * t / self.t0) ) )
+             - 2 * np.real( wofz( np.exp(1j * np.pi / 6) *
+                            np.sqrt(2 * t / self.t0) ) ) )
         return wl
 
     def __TransWakeExact(self, t, factor):
-        wt = np.real( factor *
-                     ( 2 * np.exp(-1 * t / self.t0) *
-                      ( np.sqrt(3) *
-                        np.sin(np.sqrt(3) * t / self.t0) -
-                        np.cos(np.sqrt(3) * t / self.t0) )
-                      + wofz( 1j *
-                              np.sqrt(2 * t / self.t0) )
-                      + np.exp(-1j * np.pi / 3) *
-                        wofz( np.exp(1j * np.pi / 6) *
-                              np.sqrt(2 * t / self.t0) )
-                      + np.exp(1j * np.pi / 3) *
-                        wofz( -1 * np.exp(-1j * 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)
+               - np.cos(np.sqrt(3) * t / self.t0) )
+             + np.real( wofz( 1j *
+                        np.sqrt(2 * t / self.t0) ) )
+             + 2 * np.real( np.exp(-1j * np.pi / 3) *
+                            wofz( np.exp(1j * np.pi / 6) *
+                            np.sqrt(2 * t / self.t0) ) ) )
         return wt
 
     def __LongWakeApprox(self, t):