Updated yokoya_elliptic with an analytical implementation

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,7 +189,34 @@ 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 (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:
         if y_radius < x_radius:
             small_semiaxis = y_radius
             large_semiaxis = x_radius
@@ -197,29 +224,120 @@ def yokoya_elliptic(x_radius, y_radius):
             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)
+        qr = (large_semiaxis - small_semiaxis) / (large_semiaxis + small_semiaxis)
 
-    # interpolate Yokoya file at the correct ratio
-    yoklong = 1
+        if qr == 0:
+            yoklong = 1.0
+            yokxdip = 1.0
+            yokydip = 1.0
+            yokxquad = 0.0
+            yokyquad = 0.0
 
-    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"])
+            # 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)