From 8b8a06490aebf818cfeff4b959ef2db6fb163e17 Mon Sep 17 00:00:00 2001 From: Deepak Date: Tue, 14 Jul 2020 08:22:07 -0500 Subject: [PATCH] Moving complex quad file to comlpex integrate --- pynam/dielectric/lindhard_dielectric.py | 2 +- pynam/dielectric/low_k_nam.py | 4 +- pynam/dielectric/sigma_nam.py | 4 +- pynam/noise/__init__.py | 0 pynam/noise/zeta.py | 45 +++++++++++++++++++ pynam/util/__init__.py | 2 +- pynam/util/complex_integrate.py | 29 ++++++++++++ pynam/util/complex_quad.py | 16 ------- tests/dielectric/test_lindhard_dielectric.py | 22 +++++++++ tests/noise/__init__.py | 0 tests/noise/test_zeta.py | 28 ++++++++++++ ...plex_quad.py => test_complex_integrate.py} | 4 +- 12 files changed, 132 insertions(+), 24 deletions(-) create mode 100644 pynam/noise/__init__.py create mode 100644 pynam/noise/zeta.py create mode 100644 pynam/util/complex_integrate.py delete mode 100644 pynam/util/complex_quad.py create mode 100644 tests/noise/__init__.py create mode 100644 tests/noise/test_zeta.py rename tests/util/{test_complex_quad.py => test_complex_integrate.py} (65%) diff --git a/pynam/dielectric/lindhard_dielectric.py b/pynam/dielectric/lindhard_dielectric.py index 6ace91d..b38b3c7 100644 --- a/pynam/dielectric/lindhard_dielectric.py +++ b/pynam/dielectric/lindhard_dielectric.py @@ -31,7 +31,7 @@ class LindhardDielectric(object): # converts u from inverse vacuum wavelength to inverse mean free path u = u_inverse_wavelength * self.v_f / self.c_light - if u < LINDHARD_SERIES_THRESHOLD * self.c_light / self.omega: + if u < LINDHARD_SERIES_THRESHOLD * self.v_f / self.omega: return eps_series(u) else: return eps_full_lindhard(u) diff --git a/pynam/dielectric/low_k_nam.py b/pynam/dielectric/low_k_nam.py index 2363303..bb2b5ff 100644 --- a/pynam/dielectric/low_k_nam.py +++ b/pynam/dielectric/low_k_nam.py @@ -35,7 +35,7 @@ def i2(w, wp, k, v): def a(w, k, v, t): - return pynam.util.complex_quad.complex_quadrature( + return pynam.util.complex_quad.complex_quad( lambda wp: np.tanh((w + wp) / (2 * t)) * (i1(w, wp, k, v)), 1 - w, 1 )[0] @@ -46,7 +46,7 @@ def b_int(wp, w, k, v, t): def b(w, k, v, t, b_max=np.inf): - return pynam.util.complex_quad.complex_quadrature( + return pynam.util.complex_quad.complex_quad( lambda wp: b_int(wp, w, k, v, t), 1, b_max )[0] diff --git a/pynam/dielectric/sigma_nam.py b/pynam/dielectric/sigma_nam.py index ec83275..c710282 100644 --- a/pynam/dielectric/sigma_nam.py +++ b/pynam/dielectric/sigma_nam.py @@ -41,7 +41,7 @@ def i2(w, wp, k, v): def a(w, k, v, t): - result = pynam.util.complex_quad.complex_quadrature( + result = pynam.util.complex_quad.complex_quad( lambda wp: np.tanh((w + wp) / (2 * t)) * (i1(w, wp, k, v)), 1 - w, 1, epsabs=1e-10 @@ -55,7 +55,7 @@ def b_int(wp, w, k, v, t): def b(w, k, v, t, b_max=np.inf): - return pynam.util.complex_quadrature( + return pynam.util.complex_quad( lambda wp: b_int(wp, w, k, v, t), 1, b_max )[0] diff --git a/pynam/noise/__init__.py b/pynam/noise/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/pynam/noise/zeta.py b/pynam/noise/zeta.py new file mode 100644 index 0000000..6599600 --- /dev/null +++ b/pynam/noise/zeta.py @@ -0,0 +1,45 @@ +import pynam.util + +from typing import Callable +import numpy as np + + +def get_zeta_p_integrand(eps: Callable[[float], complex]) -> Callable[[float, float], complex]: + """ Gets the integrand function zeta_p_integrand(u, y). + + Returns zeta_p_integrand(u, y), a complex valued function of two momenta in units of vacuum wavelength. + + :param eps: + :return: + """ + def zeta_p_integrand(u: float, y: float) -> complex: + """ + Here y and u are in units of vacuum wavelength, coming from Ford-Weber / from the EWJN noise expressions. + :param u: + :param y: + :return: + """ + u2 = u ** 2 + y2 = y ** 2 + k2 = u2 + y2 + k = np.sqrt(k2) + eps_value = eps(k) + term_1 = y2 / (eps_value - k2) + term_2 = u2 / eps_value + return (term_1 + term_2) / k2 + + return zeta_p_integrand + + +# def get_zeta_p_function(eps: Callable[[float], complex]): +# def zeta_p(u: float) -> complex: +# zeta_p_integrand = get_zeta_integrand(eps) +# +# integral_result = pynam.util.complex_quad(zeta_p_integrand, 0, np.inf) +# +# print(integral_result) +# integral = integral_result[0] +# +# return integral * 2j +# +# return zeta_p diff --git a/pynam/util/__init__.py b/pynam/util/__init__.py index ef592eb..cfb6416 100644 --- a/pynam/util/__init__.py +++ b/pynam/util/__init__.py @@ -1 +1 @@ -from pynam.util.complex_quad import complex_quadrature +from pynam.util.complex_quad import complex_quad, complex_quadrature diff --git a/pynam/util/complex_integrate.py b/pynam/util/complex_integrate.py new file mode 100644 index 0000000..111e109 --- /dev/null +++ b/pynam/util/complex_integrate.py @@ -0,0 +1,29 @@ +import numpy as np +from scipy.integrate import quad, quadrature + + +def complex_quad(func, a, b, **kwargs): + + def real_func(x): + return np.real(func(x)) + + def imag_func(x): + return np.imag(func(x)) + + real_integral = quad(real_func, a, b, **kwargs) + imag_integral = quad(imag_func, a, b, **kwargs) + + return real_integral[0] + 1j * imag_integral[0], real_integral[1:], imag_integral[1:] + +def complex_quadrature(func, a, b, **kwargs): + + def real_func(x): + return np.real(func(x)) + + def imag_func(x): + return np.imag(func(x)) + + real_integral = quadrature(real_func, a, b, **kwargs) + imag_integral = quadrature(imag_func, a, b, **kwargs) + + return real_integral[0] + 1j * imag_integral[0], real_integral[1:], imag_integral[1:] diff --git a/pynam/util/complex_quad.py b/pynam/util/complex_quad.py deleted file mode 100644 index fbc0189..0000000 --- a/pynam/util/complex_quad.py +++ /dev/null @@ -1,16 +0,0 @@ -import numpy as np -from scipy.integrate import quad - - -def complex_quadrature(func, a, b, **kwargs): - - def real_func(x): - return np.real(func(x)) - - def imag_func(x): - return np.imag(func(x)) - - real_integral = quad(real_func, a, b, **kwargs) - imag_integral = quad(imag_func, a, b, **kwargs) - - return real_integral[0] + 1j * imag_integral[0], real_integral[1:], imag_integral[1:] diff --git a/tests/dielectric/test_lindhard_dielectric.py b/tests/dielectric/test_lindhard_dielectric.py index 8b6bf1d..39a7e97 100644 --- a/tests/dielectric/test_lindhard_dielectric.py +++ b/tests/dielectric/test_lindhard_dielectric.py @@ -22,3 +22,25 @@ def test_lindhard_dielectric(test_input, expected): eps_to_test(test_input), expected, decimal=6, err_msg='b function is off' ) + + +@pytest.mark.parametrize("test_input,expected", [ + ((100, 100), -883.3001542404703 + 1.2566370613549341e8j), + ((100, 1e5), 5.827225842825694e7 + 3.933446612656656e7j), + ((100, 1e10), 1.0084823001646925 + 2.0013975538629039e-10j), + ((100, 1e7), 8483.300121667038 + 0.6340397839154446) +]) +def test_zeta_pi_lindhard_dielectric(zeta_p_i_epsilon, test_input, expected): + u, y = test_input + actual = zeta_p_i_epsilon(np.sqrt(u**2 + y**2)) + + np.testing.assert_allclose( + actual, expected, + rtol=10**3.8, err_msg='lindhard dielectric differs from Mathematica' + ) + + +@pytest.fixture +def zeta_p_i_epsilon(): + params = CalculationParams(omega=1e9, omega_p=3.544907701811032e15, tau=1e-14, v_f=2e6) + return pynam.dielectric.get_lindhard_dielectric(params) diff --git a/tests/noise/__init__.py b/tests/noise/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/noise/test_zeta.py b/tests/noise/test_zeta.py new file mode 100644 index 0000000..58e1fd9 --- /dev/null +++ b/tests/noise/test_zeta.py @@ -0,0 +1,28 @@ +import numpy as np +import pytest + +import pynam.dielectric +import pynam.noise.zeta +from pynam.baskets import CalculationParams + + +@pytest.fixture +def zeta_p_integrand_lindhard(): + params = CalculationParams(omega=1e9, v_f=2e6, omega_p=3.544907701811032e15, tau=1e-14) + eps_l = pynam.dielectric.get_lindhard_dielectric(params) + return pynam.noise.zeta.get_zeta_p_integrand(eps_l) + + +@pytest.mark.parametrize("test_input,expected", [ + # u y zeta_p_i(u, y) + ((100, 100), -6.891930153028566e-13 - 7.957747045025948e-9j), + ((100, 1e5), -1.0057257267146669e-10 - 4.0591966623027983e-13j), + ((1e5, 100), 1.1789175285399862e-8 - 7.957833322596519e-9j) +]) +def test_zeta_p_integrand_lindhard(zeta_p_integrand_lindhard, test_input, expected): + actual = zeta_p_integrand_lindhard(*test_input) + + np.testing.assert_allclose( + actual, expected, + rtol=1e-7, err_msg='Zeta_p is inaccurate for Lindhard case', verbose=True + ) diff --git a/tests/util/test_complex_quad.py b/tests/util/test_complex_integrate.py similarity index 65% rename from tests/util/test_complex_quad.py rename to tests/util/test_complex_integrate.py index fb69d04..0df4dc7 100644 --- a/tests/util/test_complex_quad.py +++ b/tests/util/test_complex_integrate.py @@ -1,9 +1,9 @@ import numpy as np -import pynam.util.complex_quad +import pynam.util.complex_integrate def test_complex_quad(): - actual = pynam.util.complex_quad.complex_quadrature(lambda x: x ** 2 + 1j * x ** 3, 0, 6)[0] + actual = pynam.util.complex_integrate.complex_quad(lambda x: x ** 2 + 1j * x ** 3, 0, 6)[0] # int_1^6 dx x^2 + i x^3 should equal (1/3)6^3 + (i/4)6^4 np.testing.assert_almost_equal( actual, (6**3)/3 + 1j*(6**4)/4,