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/utilities/misc.py

@@ -221,7 +221,8 @@ def yokoya_elliptic(x_radius, y_radius):
             small_semiaxis = x_radius
             large_semiaxis = y_radius
 
-        # Our equations are approximately valid only for qr (ratio) values less than or equal to 0.8.
+        # Our equations are approximately valid only for qr (ratio) values
+        # less than or equal to 0.8.
         qr_th = 0.8
         F = np.sqrt(large_semiaxis**2 - small_semiaxis**2)
         mu_b = np.arccosh(large_semiaxis/F)
@@ -305,14 +306,15 @@ def yokoya_elliptic(x_radius, y_radius):
             ff_values = np.array(function_ff(small_semiaxis_th, F_th, mu_b_th, ip, il))
             L_values = np.array(function_L(mu_b_th, ip, il))
 
-            coeff_long = np.where((ip == 0) & (il == 0), 0.25, np.where((ip == 0) | (il == 0), 0.5, 1.0))
+            coeff_long = np.where( (ip == 0) & (il == 0), 0.25,
+                                    np.where((ip == 0) | (il == 0), 0.5, 1.0) )
             coeff_quad = np.where(il == 0, 0.5, 1.0)
 
-            yoklong = np.sum(coeff_long * ff_values[0] * L_values[0])
-            yokxdip = np.sum(ff_values[1] * L_values[1])
-            yokydip = np.sum(ff_values[2] * L_values[2])
-            yokxquad = -np.sum(coeff_quad * ff_values[3] * L_values[0])
-            yokyquad = -yokxquad
+            yoklong_th = np.sum(coeff_long * ff_values[0] * L_values[0])
+            yokxdip_th = np.sum(ff_values[1] * L_values[1])
+            yokydip_th = np.sum(ff_values[2] * L_values[2])
+            yokxquad_th = -np.sum(coeff_quad * ff_values[3] * L_values[0])
+            yokyquad_th = -yokxquad_th
             
             if y_radius > x_radius:
                 yokxdip_th, yokydip_th = yokydip_th, yokxdip_th
@@ -344,7 +346,8 @@ def yokoya_elliptic(x_radius, y_radius):
             ff_values = np.array(function_ff(small_semiaxis, F, mu_b, ip, il))
             L_values = np.array(function_L(mu_b, ip, il))
 
-            coeff_long = np.where((ip == 0) & (il == 0), 0.25, np.where((ip == 0) | (il == 0), 0.5, 1.0))
+            coeff_long = np.where( (ip == 0) & (il == 0), 0.25,
+                                    np.where((ip == 0) | (il == 0), 0.5, 1.0) )
             coeff_quad = np.where(il == 0, 0.5, 1.0)
 
             yoklong = np.sum(coeff_long * ff_values[0] * L_values[0])