From a7c33fb5df1d7849db398f3470441d701e4a733a Mon Sep 17 00:00:00 2001
From: Deepak <dmallubhotla+github@gmail.com>
Date: Tue, 14 Jul 2020 10:59:33 -0500
Subject: [PATCH] adding tighter error tolerances and adding checks for imrpl

---
 pynam/noise/im_ref.py      | 15 +++++++++++++++
 pynam/noise/zeta.py        |  8 ++++----
 tests/noise/test_im_ref.py | 31 +++++++++++++++++++++++++++++++
 tests/noise/test_zeta.py   |  6 ++++++
 4 files changed, 56 insertions(+), 4 deletions(-)
 create mode 100644 pynam/noise/im_ref.py
 create mode 100644 tests/noise/test_im_ref.py

diff --git a/pynam/noise/im_ref.py b/pynam/noise/im_ref.py
new file mode 100644
index 0000000..0029857
--- /dev/null
+++ b/pynam/noise/im_ref.py
@@ -0,0 +1,15 @@
+from typing import Callable
+
+import numpy as np
+
+import pynam.noise.zeta
+
+
+def get_im_ref_p(eps: Callable[[float], complex]) -> Callable[[float], float]:
+	zeta_p = pynam.noise.zeta.get_zeta_p_function(eps)
+
+	def im_ref_p(u: float) -> float:
+		zeta_p_val = zeta_p(u)
+		return np.imag((np.pi * 1j * u - zeta_p_val) / (np.pi * 1j * u + zeta_p_val))
+
+	return im_ref_p
diff --git a/pynam/noise/zeta.py b/pynam/noise/zeta.py
index 9b44fb1..e9ff7fe 100644
--- a/pynam/noise/zeta.py
+++ b/pynam/noise/zeta.py
@@ -62,10 +62,10 @@ def get_zeta_p_function(eps: Callable[[float], complex]):
 	def zeta_p(u: float) -> complex:
 		zeta_p_integrand = get_zeta_p_integrand(eps)
 
-		i1_small = pynam.util.complex_quad(lambda x: integrand1_small_x(x, u), 0, SMALL_X_BOUNDARY)
-		i1_big = pynam.util.complex_quad(lambda x: integrand1_big_x(x, u), SMALL_X_BOUNDARY, np.inf)
-		i2_small = pynam.util.complex_quad(lambda x: integrand2_small_x(x, u), 0, SMALL_X_BOUNDARY)
-		i2_big = pynam.util.complex_quad(lambda x: integrand2_big_x(x, u), SMALL_X_BOUNDARY, np.inf)
+		i1_small = pynam.util.complex_quad(lambda x: integrand1_small_x(x, u), 0, SMALL_X_BOUNDARY, epsabs=1e-12)
+		i1_big = pynam.util.complex_quad(lambda x: integrand1_big_x(x, u), SMALL_X_BOUNDARY, np.inf, epsabs=1e-12)
+		i2_small = pynam.util.complex_quad(lambda x: integrand2_small_x(x, u), 0, SMALL_X_BOUNDARY, epsabs=1e-12)
+		i2_big = pynam.util.complex_quad(lambda x: integrand2_big_x(x, u), SMALL_X_BOUNDARY, np.inf, epsabs=1e-12)
 
 		integral = sum(term[0] for term in [i1_small, i2_small, i1_big, i2_big])
 
diff --git a/tests/noise/test_im_ref.py b/tests/noise/test_im_ref.py
new file mode 100644
index 0000000..b0f3e82
--- /dev/null
+++ b/tests/noise/test_im_ref.py
@@ -0,0 +1,31 @@
+import numpy as np
+import pytest
+
+import pynam.dielectric
+import pynam.noise.im_ref
+from pynam.baskets import CalculationParams
+
+
+@pytest.fixture
+def im_ref_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.im_ref.get_im_ref_p(eps_l)
+
+
+@pytest.mark.parametrize("test_input,expected", [
+	# u   im_ref_p_l(u)
+	# needs to be close in range around 1/z, so from 1e4 to 1e8
+	(1e4, 1.821722334939806e-8),
+	(1e5, 1.602855764970752e-8),
+	(1e6, 1.704326041013161e-8),
+	(1e7, 2.674124312031195e-8),
+	(1e8, 7.441319151047531e-8),
+])
+def test_im_ref_p_lindhard(im_ref_p_lindhard, test_input, expected):
+	actual = im_ref_p_lindhard(test_input)
+
+	np.testing.assert_allclose(
+		actual, expected,
+		rtol=1e-4, err_msg='imrp is inaccurate for Lindhard case', verbose=True
+	)
diff --git a/tests/noise/test_zeta.py b/tests/noise/test_zeta.py
index 3523678..72a51dd 100644
--- a/tests/noise/test_zeta.py
+++ b/tests/noise/test_zeta.py
@@ -38,7 +38,13 @@ def zeta_p_lindhard():
 @pytest.mark.parametrize("test_input,expected", [
 	# u    zeta_p(u)
 	(1, 0.000199609 - 0.000199608j),
+	# (10, 0.00019960929309663014 - 0.00019927000998506335j),
+	# (100, 0.0001996175250684056 - 0.0001654898843938523j),
+	# (1e3, 0.0002003339895748246 + 0.003212370020888438j),
+	# (1e4, 0.00028616168676982363 + 0.34096962141224463j),
 	(1e5, 0.0025183067257958545 + 34.11087430547122j),
+	(1e6, 0.026829658454640887 + 3411.0870128247902j),
+	(1e7, 0.4292211181081069 + 341088.797211291j),
 	(1e8, 14.348462224076096 + 3.391157983312813e7j)
 ])
 def test_zeta_p(zeta_p_lindhard, test_input, expected):