diff --git a/mbtrack2/utilities/misc.py b/mbtrack2/utilities/misc.py
index e58a848e9aed95336ced13ff1efbb71bded8319f..a8cb9704b569edc1f2a8d0469012592668696547 100644
--- a/mbtrack2/utilities/misc.py
+++ b/mbtrack2/utilities/misc.py
@@ -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)