Skip to content
Snippets Groups Projects
Commit b74830df authored by Keon Hee KIM's avatar Keon Hee KIM
Browse files

bug fix, refactor nested loops and optimize performance in misc.py

parent 80e7ba6a
No related branches found
No related tags found
2 merge requests!14Faster CircularResistiveWall,!3Yokoya elliptic (reopened)
This commit is part of merge request !14. Comments created here will be created in the context of that merge request.
Hide whitespace changes
Inline Side-by-side
mbtrack2/utilities/misc.py
+ 107
− 319
View file @ b74830df
......@@ -202,13 +202,13 @@ def yokoya_elliptic(x_radius, y_radius):
elif x_radius==np.inf and y_radius==np.inf:
raise ValueError("Both radii have infinite values.")
elif x_radius==np.inf:
yoklong = 1
yoklong = 1.0
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
yoklong = 1.0
yokxdip = np.pi**2/12
yokydip = np.pi**2/24
yokxquad = np.pi**2/24
......@@ -227,185 +227,98 @@ def yokoya_elliptic(x_radius, y_radius):
mu_b = np.arccosh(large_semiaxis/F)
qr = (large_semiaxis - small_semiaxis) / (large_semiaxis + small_semiaxis)
def function_ff(small_semiaxis, F, mu_b, ip, il):
fflong = ( 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) )
ffdx = ( 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) )
ffdy = ( 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) )
ffquad = ( 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) )
return (fflong, ffdx, ffdy, ffquad)
def function_L(mu_b, ip, il):
Lm = ( 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)) )
Ldx = ( 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)) )
Ldy = ( 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)) )
return (Lm, Ldx, Ldy)
if qr == 0:
yoklong = 1
yokxdip = 1
yokydip = 1
yokxquad = 0
yokyquad = 0
yoklong = 1.0
yokxdip = 1.0
yokydip = 1.0
yokxquad = 0.0
yokyquad = 0.0
# Beyond this threshold, they are not valid (they can even diverge),
# Beyond this threshold (qr_th), 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
yoklong_pp = 1.0
yokxdip_pp = np.pi**2/24
yokydip_pp = np.pi**2/12
yokxquad_pp = -np.pi**2/24
yokyquad_pp = np.pi**2/24
if y_radius > x_radius:
yokxdip_pp, yokydip_pp = yokydip_pp, yokxdip_pp
yokxquad_pp, yokyquad_pp = yokyquad_pp, yokxquad_pp
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_th = 0.0
yokxdip_th = 0.0
yokydip_th = 0.0
yokxquad_th = 0.0
ip_range = np.arange(50)
il_range = np.arange(50)
ip, il = np.meshgrid(ip_range, il_range, indexing='ij')
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_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
if y_radius > x_radius:
yokxdip_th, yokydip_th = yokydip_th, yokxdip_th
yokxquad_th, yokyquad_th = yokyquad_th, yokxquad_th
qr_array = np.array([qr_th, 1.0])
yoklong_array = np.array([yoklong_th, yoklong_pp])
yokxdip_array = np.array([yokxdip_th, yokxdip_pp])
yokydip_array = np.array([yokydip_th, yokydip_pp])
......@@ -419,155 +332,30 @@ def yokoya_elliptic(x_radius, y_radius):
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
yoklong = 0.0
yokxdip = 0.0
yokydip = 0.0
yokxquad = 0.0
ip_range = np.arange(50)
il_range = np.arange(50)
ip, il = np.meshgrid(ip_range, il_range, indexing='ij')
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_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
if y_radius > x_radius:
yokxdip, yokydip = yokydip, yokxdip
yokxquad, yokyquad = yokyquad, yokxquad
return (yoklong, yokxdip, yokydip, yokxquad, yokyquad)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

The place to be for your source code