diff --git a/mbtrack2/utilities/misc.py b/mbtrack2/utilities/misc.py
index c31af03d547338323c8ca1d98e886eefd02fe853..6aaabef166b458927b868cfa816dd967fb711cba 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,156 @@ 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] : M. Migliorati, L. Palumbo, C. Zannini, N. Biancacci, and
+ V. G. Vaccaro, "Resistive wall impedance in elliptical multilayer vacuum
+ chambers." Phys. Rev. Accel. Beams 22, 121001 (2019).
+ [2] : R.L. Gluckstern, J. van Zeijts, and B. Zotter, "Coupling impedance
+ of beam pipes of general cross section." Phys. Rev. E 47, 656 (1993).
+
"""
- if y_radius < x_radius:
- small_semiaxis = y_radius
- large_semiaxis = x_radius
- 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"])
+
+ if (x_radius <= 0) or (y_radius <= 0):
+ raise ValueError("Both radii must be non-zero positive 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.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.0
+ yokxdip = np.pi**2/12
+ yokydip = np.pi**2/24
+ yokxquad = np.pi**2/24
+ yokyquad = -np.pi**2/24
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
+
+ qr = (large_semiaxis - small_semiaxis) / (large_semiaxis + small_semiaxis)
+ if qr == 0:
+ yoklong = 1.0
+ yokxdip = 1.0
+ yokydip = 1.0
+ yokxquad = 0.0
+ yokyquad = 0.0
+
+ else:
+ # Define form factor functions
+ def function_ff(small_semiaxis, F, mu_b, ip, il):
+ coeff_fflong = 2*np.sqrt(2)*small_semiaxis/(np.pi*F)
+ coeff_fftrans = np.sqrt(2)*small_semiaxis**3/(np.pi*F**3)
+
+ fflong = ( coeff_fflong*(-1)**ip/np.cosh(2*ip*mu_b)
+ * (-1)**il/np.cosh(2*il*mu_b) )
+ ffdx = ( coeff_fftrans*(-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 = ( coeff_fftrans*(-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 = ( coeff_fftrans*(-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):
+ common_L = ( 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)) )
+ val_m = ( 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)) )
+ val_d = ( 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)) )
+
+ Lm = common_L + val_m
+ Ldx = common_L + val_d
+ Ldy = common_L - val_d
+
+ return (Lm, Ldx, Ldy)
+
+ ip_range = np.arange(51)
+ il_range = np.arange(51)
+ ip, il = np.meshgrid(ip_range, il_range, indexing="ij")
+
+ 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)
+
+ # Our equations are approximately valid only for qr (ratio) values
+ # less than or equal to 0.8.
+ qr_th = 0.8
+
+ if qr < qr_th:
+ F = np.sqrt(large_semiaxis**2 - small_semiaxis**2)
+ mu_b = np.arccosh(large_semiaxis/F)
+
+ ff_values = np.array(function_ff(small_semiaxis, F, mu_b, ip, il))
+ L_values = np.array(function_L(mu_b, ip, il))
+
+ 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
+
+ # Beyond the threshold (qr >= qr_th), they may not be valid,
+ # but should converge to Yokoya form factors of parallel plates.
+ # Fortunately, beyond this threshold, they should show asymptotic behavior,
+ # so we perform linear interpolation.
+ else:
+ 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
+
+ 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)
+
+ ff_values_th = np.array(function_ff(small_semiaxis_th, F_th, mu_b_th, ip, il))
+ L_values_th = np.array(function_L(mu_b_th, ip, il))
+
+ yoklong_th = np.sum(coeff_long * ff_values_th[0] * L_values_th[0])
+ yokxdip_th = np.sum(ff_values_th[1] * L_values_th[1])
+ yokydip_th = np.sum(ff_values_th[2] * L_values_th[2])
+ yokxquad_th = -np.sum(coeff_quad * ff_values_th[3] * L_values_th[0])
+ yokyquad_th = -yokxquad_th
+
+ if y_radius > x_radius:
+ yokxdip_th, yokydip_th = yokydip_th, yokxdip_th
+ yokxquad_th, yokyquad_th = yokyquad_th, yokxquad_th
+ yokxdip_pp, yokydip_pp = yokydip_pp, yokxdip_pp
+ yokxquad_pp, yokyquad_pp = yokyquad_pp, yokxquad_pp
+
+ 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])
+ 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)
+
return (yoklong, yokxdip, yokydip, yokxquad, yokyquad)