From 6ecb74f1b9dc5a2295b3ecbe2aa911826fdc473c Mon Sep 17 00:00:00 2001 From: Deepak Mallubhotla Date: Wed, 8 Jun 2022 09:37:39 -0500 Subject: [PATCH] feat: adds r_s and chi_zz_b functions --- .../lindhard_dielectric_transverse.py | 9 ++++--- pyewjn/noise/chi.py | 15 +++++++++++ pyewjn/noise/im_ref.py | 17 ++++++++++++ .../test_lindhard_dielectric_transverse.py | 27 +++++++++++++++++++ 4 files changed, 64 insertions(+), 4 deletions(-) create mode 100644 tests/dielectric/test_lindhard_dielectric_transverse.py diff --git a/pyewjn/dielectric/lindhard_dielectric_transverse.py b/pyewjn/dielectric/lindhard_dielectric_transverse.py index 2d15a24..9cff203 100644 --- a/pyewjn/dielectric/lindhard_dielectric_transverse.py +++ b/pyewjn/dielectric/lindhard_dielectric_transverse.py @@ -34,10 +34,11 @@ 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: - return eps_series(q) - else: - return eps_full_lindhard(q) + # if u_inverse_wavelength < self.series_threshold * self.v_f / self.omega: + # return eps_series(q) + # else: + # return eps_full_lindhard(q) + return eps_full_lindhard(q) def eps_series(q: float) -> complex: pass diff --git a/pyewjn/noise/chi.py b/pyewjn/noise/chi.py index 290c5af..63d5096 100644 --- a/pyewjn/noise/chi.py +++ b/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 diff --git a/pyewjn/noise/im_ref.py b/pyewjn/noise/im_ref.py index 8574800..5e4abfd 100644 --- a/pyewjn/noise/im_ref.py +++ 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 diff --git a/tests/dielectric/test_lindhard_dielectric_transverse.py b/tests/dielectric/test_lindhard_dielectric_transverse.py new file mode 100644 index 0000000..18afd3d --- /dev/null +++ b/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" + )