From ab61ba5713adbff9942efad8ed2e38735965c140 Mon Sep 17 00:00:00 2001 From: Deepak Date: Tue, 14 Jul 2020 08:22:18 -0500 Subject: [PATCH] adding zeta stuff --- pynam/dielectric/low_k_nam.py | 6 +++--- pynam/dielectric/sigma_nam.py | 2 +- pynam/noise/zeta.py | 26 +++++++++++++------------- pynam/util/__init__.py | 2 +- pynam/util/complex_integrate.py | 1 + tests/noise/test_zeta.py | 26 +++++++++++++++++++++++--- tests/util/test_complex_integrate.py | 9 +++++++++ 7 files changed, 51 insertions(+), 21 deletions(-) diff --git a/pynam/dielectric/low_k_nam.py b/pynam/dielectric/low_k_nam.py index bb2b5ff..02c4ba8 100644 --- a/pynam/dielectric/low_k_nam.py +++ b/pynam/dielectric/low_k_nam.py @@ -1,7 +1,7 @@ import numpy as np from numpy.lib.scimath import sqrt as csqrt -import pynam.util.complex_quad +import pynam.util def g(w, wp): @@ -35,7 +35,7 @@ def i2(w, wp, k, v): def a(w, k, v, t): - return pynam.util.complex_quad.complex_quad( + return pynam.util.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_quad( + return pynam.util.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 c710282..2b4d4cb 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_quad( + result = pynam.util.complex_quad( lambda wp: np.tanh((w + wp) / (2 * t)) * (i1(w, wp, k, v)), 1 - w, 1, epsabs=1e-10 diff --git a/pynam/noise/zeta.py b/pynam/noise/zeta.py index 6599600..125c11e 100644 --- a/pynam/noise/zeta.py +++ b/pynam/noise/zeta.py @@ -12,7 +12,7 @@ def get_zeta_p_integrand(eps: Callable[[float], complex]) -> Callable[[float, fl :param eps: :return: """ - def zeta_p_integrand(u: float, y: float) -> complex: + def zeta_p_integrand(y: float, u: float) -> complex: """ Here y and u are in units of vacuum wavelength, coming from Ford-Weber / from the EWJN noise expressions. :param u: @@ -31,15 +31,15 @@ def get_zeta_p_integrand(eps: Callable[[float], complex]) -> Callable[[float, fl 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 +def get_zeta_p_function(eps: Callable[[float], complex]): + def zeta_p(u: float) -> complex: + zeta_p_integrand = get_zeta_p_integrand(eps) + + integral_result = pynam.util.complex_quad(lambda y: zeta_p_integrand(y, u), 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 cfb6416..4a4cfcb 100644 --- a/pynam/util/__init__.py +++ b/pynam/util/__init__.py @@ -1 +1 @@ -from pynam.util.complex_quad import complex_quad, complex_quadrature +from pynam.util.complex_integrate import complex_quad, complex_quadrature diff --git a/pynam/util/complex_integrate.py b/pynam/util/complex_integrate.py index 111e109..cf1540c 100644 --- a/pynam/util/complex_integrate.py +++ b/pynam/util/complex_integrate.py @@ -15,6 +15,7 @@ def complex_quad(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): diff --git a/tests/noise/test_zeta.py b/tests/noise/test_zeta.py index 58e1fd9..0c73ce4 100644 --- a/tests/noise/test_zeta.py +++ b/tests/noise/test_zeta.py @@ -14,10 +14,10 @@ def zeta_p_integrand_lindhard(): @pytest.mark.parametrize("test_input,expected", [ - # u y zeta_p_i(u, y) + # y u 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) + ((1e5, 100), -1.0057257267146669e-10 - 4.0591966623027983e-13j), + ((100, 1e5), 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) @@ -26,3 +26,23 @@ def test_zeta_p_integrand_lindhard(zeta_p_integrand_lindhard, test_input, expect actual, expected, rtol=1e-7, err_msg='Zeta_p is inaccurate for Lindhard case', verbose=True ) + + +@pytest.fixture +def zeta_p_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_function(eps_l) + + +@pytest.mark.parametrize("test_input,expected", [ + # u zeta_p(u) + (1, 0.000199609 - 0.000199608j), +]) +def test_zeta_p(zeta_p_lindhard, test_input, expected): + actual = zeta_p_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_integrate.py b/tests/util/test_complex_integrate.py index 0df4dc7..6088896 100644 --- a/tests/util/test_complex_integrate.py +++ b/tests/util/test_complex_integrate.py @@ -9,3 +9,12 @@ def test_complex_quad(): actual, (6**3)/3 + 1j*(6**4)/4, decimal=7, err_msg='complex quadrature is broken', verbose=True ) + + +def test_complex_quadrature(): + actual = pynam.util.complex_integrate.complex_quadrature(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, + decimal=7, err_msg='complex quadrature is broken', verbose=True + )