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

@@ -8,6 +8,7 @@ from pathlib import Path
 import numpy as np
 import pandas as pd
 from scipy.interpolate import interp1d
+from scipy.special import gamma, factorial, hyp2f1
 
 from mbtrack2.impedance.wakefield import Impedance
 from mbtrack2.utilities.spectrum import spectral_density
@@ -168,7 +169,6 @@ def tune_shift_from_effective_impedance(Zeff):
 def yokoya_elliptic(x_radius, y_radius):
     """
     Compute Yokoya factors for an elliptic beam pipe.
-    Function adapted from N. Mounet IW2D.
 
     Parameters
     ----------
@@ -189,38 +189,386 @@ def yokoya_elliptic(x_radius, y_radius):
         Yokoya factor for the quadrupolar horizontal impedance.
     yokyquad : float
         Yokoya factor for the quadrupolar vertical impedance.
+    
+    References
+    ----------
+    [1] : Phys. Rev. Accel. Beams 22, 121001 (2019)
+    [2] : Phys. Rev. E 47, 656 (1993)
+
     """
-    if y_radius < x_radius:
-        small_semiaxis = y_radius
-        large_semiaxis = x_radius
+    
+    if x_radius==0 or y_radius==0:
+        raise ValueError("Both radii must be non-zero values.")
+    elif x_radius==np.inf and y_radius==np.inf:
+        raise ValueError("Both radii have infinite values.")
+    elif x_radius==np.inf:
+        yoklong = 1
+        yokxdip = np.pi**2/24
+        yokydip = np.pi**2/12
+        yokxquad = -np.pi**2/24
+        yokyquad = np.pi**2/24
+    elif y_radius==np.inf:
+        yoklong = 1
+        yokxdip = np.pi**2/12
+        yokydip = np.pi**2/24
+        yokxquad = np.pi**2/24
+        yokyquad = -np.pi**2/24
     else:
-        small_semiaxis = x_radius
-        large_semiaxis = y_radius
-
-    path_to_file = Path(__file__).parent
-    file = path_to_file / "data" / "Yokoya_elliptic_from_Elias_USPAS.csv"
-
-    # read Yokoya factors interpolation file
-    # BEWARE: columns are ratio, dipy, dipx, quady, quadx
-    yokoya_file = pd.read_csv(file)
-    ratio_col = yokoya_file["x"]
-    # compute semi-axes ratio (first column of this file)
-    ratio = (large_semiaxis-small_semiaxis) / (large_semiaxis+small_semiaxis)
-
-    # interpolate Yokoya file at the correct ratio
-    yoklong = 1
-
-    if y_radius < x_radius:
-        yokydip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
-        yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
-        yokyquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
-        yokxquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])
-    else:
-        yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
-        yokydip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
-        yokxquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
-        yokyquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])
+        if y_radius < x_radius:
+            small_semiaxis = y_radius
+            large_semiaxis = x_radius
+        else:
+            small_semiaxis = x_radius
+            large_semiaxis = y_radius
+
+        # Our equations are approximately valid only for qr (ratio) values less than or equal to 0.9.
+        qr_th = 0.9
+        F = np.sqrt(large_semiaxis**2 - small_semiaxis**2)
+        mu_b = np.arccosh(large_semiaxis/F)
+        qr = (large_semiaxis - small_semiaxis) / (large_semiaxis + small_semiaxis)
+
+        if qr == 0:
+            yoklong = 1
+            yokxdip = 1
+            yokydip = 1
+            yokxquad = 0
+            yokyquad = 0
+
+        # Beyond this threshold, they are not valid (they can even diverge),
+        # but should converge to Yokoya form factors of parallel plates.
+        # Fortunately, beyond this threshold, they show asymptotic behavior,
+        # so we perform linear interpolation.
+        elif qr > qr_th:
+            if y_radius < x_radius:
+                yoklong_pp = 1
+                yokxdip_pp = np.pi**2/24
+                yokydip_pp = np.pi**2/12
+                yokxquad_pp = -np.pi**2/24
+                yokyquad_pp = np.pi**2/24
+            else:
+                yoklong_pp = 1
+                yokxdip_pp = np.pi**2/12
+                yokydip_pp = np.pi**2/24
+                yokxquad_pp = np.pi**2/24
+                yokyquad_pp = -np.pi**2/24
+
+            small_semiaxis_th = large_semiaxis*(1-qr_th)/(1+qr_th)
+            F_th = np.sqrt(large_semiaxis**2 - small_semiaxis_th**2)
+            mu_b_th = np.arccosh(large_semiaxis/F_th)
+
+            yoklong_th = 0
+            yokxdip_th = 0
+            yokydip_th = 0
+            yokxquad_th = 0
+        
+            if y_radius < x_radius:
+                for ip in range(50):
+                    for il in range(50):
+                        if (ip==0) and (il==0):
+                            yoklong_th += ( np.sqrt(2)*small_semiaxis_th/(2*np.pi*F_th)
+                                           *(-1)**ip/np.cosh(2*ip*mu_b_th)
+                                           *(-1)**il/np.cosh(2*il*mu_b_th)
+                                           * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                           *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                           *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                           + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b_th)
+                                           *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                           *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b_th)) ) )
+                        elif (ip==0 and il!=0) or (ip!=0 and il==0):
+                            yoklong_th += ( np.sqrt(2)*small_semiaxis_th/(np.pi*F_th)
+                                           *(-1)**ip/np.cosh(2*ip*mu_b_th)
+                                           *(-1)**il/np.cosh(2*il*mu_b_th)
+                                           * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                           *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                           *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                           + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b_th)
+                                           *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                           *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b_th)) ) )
+                        else:
+                            yoklong_th += ( 2*np.sqrt(2)*small_semiaxis_th/(np.pi*F_th)
+                                           *(-1)**ip/np.cosh(2*ip*mu_b_th)
+                                           *(-1)**il/np.cosh(2*il*mu_b_th)
+                                           * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                           *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                           *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                           + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b_th)
+                                           *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                           *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b_th)) ) )
+                        yokxdip_th += ( np.sqrt(2)*small_semiaxis_th**3/(np.pi*F_th**3)
+                                       *(-1)**ip*(2*ip+1)/np.cosh((2*ip+1)*mu_b_th)
+                                       *(-1)**il*(2*il+1)/np.cosh((2*il+1)*mu_b_th)
+                                       * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                       *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                       *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                       + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+3)*mu_b_th)
+                                       *gamma(1.5+ip+il)/(gamma(0.5)*factorial(ip+il+1))
+                                       *hyp2f1(0.5, ip+il+1.5, ip+il+2, np.exp(-4*mu_b_th)) ) )
+                        yokydip_th += ( np.sqrt(2)*small_semiaxis_th**3/(np.pi*F_th**3)
+                                       *(-1)**ip*(2*ip+1)/np.sinh((2*ip+1)*mu_b_th)
+                                       *(-1)**il*(2*il+1)/np.sinh((2*il+1)*mu_b_th)
+                                       * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                       *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                       *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                       - np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+3)*mu_b_th)
+                                       *gamma(1.5+ip+il)/(gamma(0.5)*factorial(ip+il+1))
+                                       *hyp2f1(0.5, ip+il+1.5, ip+il+2, np.exp(-4*mu_b_th)) ) )
+                        if (il==0):
+                            yokxquad_th -= ( np.sqrt(2)*small_semiaxis_th**3/(2*np.pi*F_th**3)
+                                            *(-1)**ip*(2*ip)**2/np.cosh(2*ip*mu_b_th)
+                                            *(-1)**il/np.cosh(2*il*mu_b_th)
+                                            * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                            *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                            *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                            + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b_th)
+                                            *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                            *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b_th)) ) )
+                        else:
+                            yokxquad_th -= ( np.sqrt(2)*small_semiaxis_th**3/(np.pi*F_th**3)
+                                            *(-1)**ip*(2*ip)**2/np.cosh(2*ip*mu_b_th)
+                                            *(-1)**il/np.cosh(2*il*mu_b_th)
+                                            * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                            *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                            *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                            + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b_th)
+                                            *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                            *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b_th)) ) )
+                yokyquad_th = -yokxquad_th
+            else:
+                for ip in range(50):
+                    for il in range(50):
+                        if (ip==0) and (il==0):
+                            yoklong_th += ( np.sqrt(2)*small_semiaxis_th/(2*np.pi*F_th)
+                                           *(-1)**ip/np.cosh(2*ip*mu_b_th)
+                                           *(-1)**il/np.cosh(2*il*mu_b_th)
+                                           * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                           *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                           *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                           + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b_th)
+                                           *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                           *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b_th)) ) )
+                        elif (ip==0 and il!=0) or (ip!=0 and il==0):
+                            yoklong_th += ( np.sqrt(2)*small_semiaxis_th/(np.pi*F_th)
+                                           *(-1)**ip/np.cosh(2*ip*mu_b_th)
+                                           *(-1)**il/np.cosh(2*il*mu_b_th)
+                                           * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                           *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                           *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                           + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b_th)
+                                           *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                           *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b_th)) ) )
+                        else:
+                            yoklong_th += ( 2*np.sqrt(2)*small_semiaxis_th/(np.pi*F_th)
+                                           *(-1)**ip/np.cosh(2*ip*mu_b_th)*(-1)**il/np.cosh(2*il*mu_b_th)
+                                           * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                           *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                           *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                           + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b_th)
+                                           *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                           *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b_th)) ) )
+                        yokydip_th += ( np.sqrt(2)*small_semiaxis_th**3/(np.pi*F_th**3)
+                                       *(-1)**ip*(2*ip+1)/np.cosh((2*ip+1)*mu_b_th)
+                                       *(-1)**il*(2*il+1)/np.cosh((2*il+1)*mu_b_th)
+                                       * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                       *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                       *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                       + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+3)*mu_b_th)
+                                       *gamma(1.5+ip+il)/(gamma(0.5)*factorial(ip+il+1))
+                                       *hyp2f1(0.5, ip+il+1.5, ip+il+2, np.exp(-4*mu_b_th)) ) )
+                        yokxdip_th += ( np.sqrt(2)*small_semiaxis_th**3/(np.pi*F_th**3)
+                                       *(-1)**ip*(2*ip+1)/np.sinh((2*ip+1)*mu_b_th)
+                                       *(-1)**il*(2*il+1)/np.sinh((2*il+1)*mu_b_th)
+                                       * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                       *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                       *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                       - np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+3)*mu_b_th)
+                                       *gamma(1.5+ip+il)/(gamma(0.5)*factorial(ip+il+1))
+                                       *hyp2f1(0.5, ip+il+1.5, ip+il+2, np.exp(-4*mu_b_th)) ) )
+                        if (il==0):
+                            yokyquad_th -= ( np.sqrt(2)*small_semiaxis_th**3/(2*np.pi*F_th**3)
+                                            *(-1)**ip*(2*ip)**2/np.cosh(2*ip*mu_b_th)
+                                            *(-1)**il/np.cosh(2*il*mu_b_th)
+                                            * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                            *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                            *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                            + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b_th)
+                                            *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                            *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b_th)) ) )
+                        else:
+                            yokyquad_th -= ( np.sqrt(2)*small_semiaxis_th**3/(np.pi*F_th**3)
+                                            *(-1)**ip*(2*ip)**2/np.cosh(2*ip*mu_b_th)
+                                            *(-1)**il/np.cosh(2*il*mu_b_th)
+                                            * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b_th)
+                                            *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                            *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b_th))
+                                            + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b_th)
+                                            *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                            *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b_th)) ) )
+                yokxquad_th = -yokyquad_th
+
+            qr_array = np.array([qr_th, 1])
+            yoklong_array = np.array([yoklong_th, yoklong_pp])
+            yokxdip_array = np.array([yokxdip_th, yokxdip_pp])
+            yokydip_array = np.array([yokydip_th, yokydip_pp])
+            yokxquad_array = np.array([yokxquad_th, yokxquad_pp])
+            yokyquad_array = np.array([yokyquad_th, yokyquad_pp])
+
+            yoklong = np.interp(qr, qr_array, yoklong_array)
+            yokxdip = np.interp(qr, qr_array, yokxdip_array)
+            yokydip = np.interp(qr, qr_array, yokydip_array)
+            yokxquad = np.interp(qr, qr_array, yokxquad_array)
+            yokyquad = np.interp(qr, qr_array, yokyquad_array)
 
+        else:
+            yoklong = 0
+            yokxdip = 0
+            yokydip = 0
+            yokxquad = 0
+
+            if y_radius < x_radius:
+                for ip in range(50):
+                    for il in range(50):
+                        if (ip==0) and (il==0):
+                            yoklong += ( np.sqrt(2)*small_semiaxis/(2*np.pi*F)
+                                        *(-1)**ip/np.cosh(2*ip*mu_b)
+                                        *(-1)**il/np.cosh(2*il*mu_b)
+                                        * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                        *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                        *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                        + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b)
+                                        *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                        *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b)) ) )
+                        elif (ip==0 and il!=0) or (ip!=0 and il==0):
+                            yoklong += ( np.sqrt(2)*small_semiaxis/(np.pi*F)
+                                        *(-1)**ip/np.cosh(2*ip*mu_b)
+                                        *(-1)**il/np.cosh(2*il*mu_b)
+                                        * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                        *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                        *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                        + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b)
+                                        *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                        *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b)) ) )
+                        else:
+                            yoklong += ( 2*np.sqrt(2)*small_semiaxis/(np.pi*F)
+                                        *(-1)**ip/np.cosh(2*ip*mu_b)
+                                        *(-1)**il/np.cosh(2*il*mu_b)
+                                        * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                        *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                        *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                        + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b)
+                                        *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                        *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b)) ) )
+                        yokxdip += ( np.sqrt(2)*small_semiaxis**3/(np.pi*F**3)
+                                    *(-1)**ip*(2*ip+1)/np.cosh((2*ip+1)*mu_b)
+                                    *(-1)**il*(2*il+1)/np.cosh((2*il+1)*mu_b)
+                                    * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                    *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                    *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                    + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+3)*mu_b)
+                                    *gamma(1.5+ip+il)/(gamma(0.5)*factorial(ip+il+1))
+                                    *hyp2f1(0.5, ip+il+1.5, ip+il+2, np.exp(-4*mu_b)) ) )
+                        yokydip += ( np.sqrt(2)*small_semiaxis**3/(np.pi*F**3)
+                                    *(-1)**ip*(2*ip+1)/np.sinh((2*ip+1)*mu_b)
+                                    *(-1)**il*(2*il+1)/np.sinh((2*il+1)*mu_b)
+                                    * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                    *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                    *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                    - np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+3)*mu_b)
+                                    *gamma(1.5+ip+il)/(gamma(0.5)*factorial(ip+il+1))
+                                    *hyp2f1(0.5, ip+il+1.5, ip+il+2, np.exp(-4*mu_b)) ) )
+                        if (il==0):
+                            yokxquad -= ( np.sqrt(2)*small_semiaxis**3/(2*np.pi*F**3)
+                                         *(-1)**ip*(2*ip)**2/np.cosh(2*ip*mu_b)
+                                         *(-1)**il/np.cosh(2*il*mu_b)
+                                         * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                         *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                         *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                         + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b)
+                                         *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                         *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b)) ) )
+                        else:
+                            yokxquad -= ( np.sqrt(2)*small_semiaxis**3/(np.pi*F**3)
+                                         *(-1)**ip*(2*ip)**2/np.cosh(2*ip*mu_b)
+                                         *(-1)**il/np.cosh(2*il*mu_b)
+                                         * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                         *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                         *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                         + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b)
+                                         *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                         *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b)) ) )
+                yokyquad = -yokxquad
+            else:
+                for ip in range(50):
+                    for il in range(50):
+                        if (ip==0) and (il==0):
+                            yoklong += ( np.sqrt(2)*small_semiaxis/(2*np.pi*F)
+                                        *(-1)**ip/np.cosh(2*ip*mu_b)
+                                        *(-1)**il/np.cosh(2*il*mu_b)
+                                        * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                        *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                        *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                        + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b)
+                                        *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                        *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b)) ) )
+                        elif (ip==0 and il!=0) or (ip!=0 and il==0):
+                            yoklong += ( np.sqrt(2)*small_semiaxis/(np.pi*F)
+                                        *(-1)**ip/np.cosh(2*ip*mu_b)
+                                        *(-1)**il/np.cosh(2*il*mu_b)
+                                        * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                        *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                        *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                        + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b)
+                                        *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                        *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b)) ) )
+                        else:
+                            yoklong += ( 2*np.sqrt(2)*small_semiaxis/(np.pi*F)
+                                        *(-1)**ip/np.cosh(2*ip*mu_b)
+                                        *(-1)**il/np.cosh(2*il*mu_b)
+                                        * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                        *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                        *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                        + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b)
+                                        *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                        *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b)) ) )
+                        yokydip += ( np.sqrt(2)*small_semiaxis**3/(np.pi*F**3)
+                                    *(-1)**ip*(2*ip+1)/np.cosh((2*ip+1)*mu_b)
+                                    *(-1)**il*(2*il+1)/np.cosh((2*il+1)*mu_b)
+                                    * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                    *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                    *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                    + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+3)*mu_b)
+                                    *gamma(1.5+ip+il)/(gamma(0.5)*factorial(ip+il+1))
+                                    *hyp2f1(0.5, ip+il+1.5, ip+il+2, np.exp(-4*mu_b)) ) )
+                        yokxdip += ( np.sqrt(2)*small_semiaxis**3/(np.pi*F**3)
+                                    *(-1)**ip*(2*ip+1)/np.sinh((2*ip+1)*mu_b)
+                                    *(-1)**il*(2*il+1)/np.sinh((2*il+1)*mu_b)
+                                    * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                    *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                    *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                    - np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+3)*mu_b)
+                                    *gamma(1.5+ip+il)/(gamma(0.5)*factorial(ip+il+1))
+                                    *hyp2f1(0.5, ip+il+1.5, ip+il+2, np.exp(-4*mu_b)) ) )
+                        if (il==0):
+                            yokyquad -= ( np.sqrt(2)*small_semiaxis**3/(2*np.pi*F**3)
+                                         *(-1)**ip*(2*ip)**2/np.cosh(2*ip*mu_b)
+                                         *(-1)**il/np.cosh(2*il*mu_b)
+                                         * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                         *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                         *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                         + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b)
+                                         *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                         *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b)) ) )
+                        else:
+                            yokyquad -= ( np.sqrt(2)*small_semiaxis**3/(np.pi*F**3)
+                                         *(-1)**ip*(2*ip)**2/np.cosh(2*ip*mu_b)
+                                         *(-1)**il/np.cosh(2*il*mu_b)
+                                         * ( np.sqrt(2)*np.pi*np.exp(-(2*abs(ip-il)+1)*mu_b)
+                                         *gamma(0.5+abs(ip-il))/(gamma(0.5)*factorial(abs(ip-il)))
+                                         *hyp2f1(0.5, abs(ip-il)+0.5, abs(ip-il)+1, np.exp(-4*mu_b))
+                                         + np.sqrt(2)*np.pi*np.exp(-(2*ip+2*il+1)*mu_b)
+                                         *gamma(0.5+ip+il)/(gamma(0.5)*factorial(ip+il))
+                                         *hyp2f1(0.5, ip+il+0.5, ip+il+1, np.exp(-4*mu_b)) ) )
+                yokxquad = -yokyquad
+    
     return (yoklong, yokxdip, yokydip, yokxquad, yokyquad)