Skip to content
Snippets Groups Projects

Yokoya elliptic

Closed Yokoya elliptic
yokoya_elliptic into develop
Closed Keon Hee KIM requested to merge yokoya_elliptic into develop
Compare
and
1 file
+ 378
− 30
Compare changes
  • Side-by-side
  • Inline
mbtrack2/utilities/misc.py
+ 378
− 30
@@ -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 should 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)

The place to be for your source code