feat: adds r_s and chi_zz_b functions

pyewjn/dielectric/lindhard_dielectric_transverse.py

@@ -34,9 +34,10 @@ class LindhardDielectricTransverse(object):
 			# want to convert to q = vf k, where k is wavevector in SI units
 			q = u_inverse_wavelength * (self.v_f * self.omega) / (self.c_light)
 
-			if u_inverse_wavelength < self.series_threshold * self.v_f / self.omega:
+			# if u_inverse_wavelength < self.series_threshold * self.v_f / self.omega:
-				return eps_series(q)
+			# 	return eps_series(q)
-			else:
+			# else:
+			# 	return eps_full_lindhard(q)
 			return eps_full_lindhard(q)
 
 		def eps_series(q: float) -> complex:

pyewjn/noise/chi.py

@@ -19,3 +19,18 @@ def get_chi_zz_e(eps: Callable[[float], complex]) -> Callable[[float], float]:
 		return integral[0] / (z**3)
 
 	return chi_zz_e
+
+
+def get_chi_zz_b(eps_t: Callable[[float], complex]) -> Callable[[float], float]:
+	im_ref_s = pyewjn.noise.im_ref.get_im_ref_s(eps_t)
+
+	def chi_zz_b(z: float) -> float:
+		def integrand(y: float) -> float:
+			return (y**2) * im_ref_s(y / z) * np.exp(-2 * y)
+
+		integral = scipy.integrate.quad(
+			integrand, 0, np.inf, epsabs=1e-10, epsrel=1e-10
+		)
+		return integral[0] / (z**3)
+
+	return chi_zz_b

pyewjn/noise/im_ref.py

@@ -3,6 +3,7 @@ from typing import Callable
 import numpy as np
 
 import pyewjn.noise.zeta
+import pyewjn.util
 
 
 def get_im_ref_p(eps: Callable[[float], complex]) -> Callable[[float], float]:
@@ -13,3 +14,19 @@ def get_im_ref_p(eps: Callable[[float], complex]) -> Callable[[float], float]:
 		return np.imag((np.pi * 1j * u - zeta_p_val) / (np.pi * 1j * u + zeta_p_val))
 
 	return im_ref_p
+
+
+def get_im_ref_s(eps_t: Callable[[float], complex]) -> Callable[[float], float]:
+	def integrand(kappa: float, u: float) -> complex:
+		k_eff2 = kappa**2 + u**2
+		k_eff4 = k_eff2**2
+		k_eff = k_eff2 ** (1 / 2)
+		return eps_t(k_eff) / k_eff4
+
+	def im_ref_s(u: float) -> float:
+		integral = pyewjn.util.complex_quad(
+			lambda kappa: integrand(kappa, u), 0, np.inf, epsabs=1e-12
+		)
+		return (1 / (4 * u**2)) * ((4 * u**3 * integral) - 1)
+
+	return im_ref_s

tests/dielectric/test_lindhard_dielectric_transverse.py

@@ -0,0 +1,27 @@
+import pyewjn.dielectric
+import numpy as np
+import pytest
+from pyewjn.baskets import CalculationParams
+
+
+def get_common_lindhard_dielectric():
+	params = CalculationParams(omega=1e9, omega_p=3.5e15, tau=1e-14, v_f=2e6)
+	return pyewjn.dielectric.get_lindhard_dielectric_transverse(params)
+
+
+@pytest.mark.skip(reason="Not actually correct values")
+@pytest.mark.parametrize(
+	"test_input,expected",
+	[
+		(10, -1222.185185062794 + 1.2249999998777178e8j),
+		(1000, 16924.14814718176 + 1.2250000020552777e8j),
+		(1e8, 83.687499999706 + 0.00022417398943752126j),
+	],
+)
+def test_lindhard_dielectric_transverse(test_input, expected):
+
+	eps_to_test = get_common_lindhard_dielectric()
+
+	np.testing.assert_almost_equal(
+		eps_to_test(test_input), expected, decimal=6, err_msg="b function is off"
+	)