diff --git a/mbtrack2/utilities/misc.py b/mbtrack2/utilities/misc.py
index c31af03d547338323c8ca1d98e886eefd02fe853..6addd0211d7f9344dcfd0a7d97c7737546b0f863 100644
--- a/mbtrack2/utilities/misc.py
+++ b/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)