fix: fixes stuff with nam integration

pyewjn/baskets.py

@@ -32,7 +32,7 @@ class CalculationParams(object):
 		tau: float,
 		v_f: float,
 		t_rel: float = 0.8,
-		t_c: float = 1e-11,
+		t_c: float = 1e11,
 		dipole_moment: float = 1,
 	):
 		"""Creates parameter object, SI units

pyewjn/dielectric/__init__.py

@@ -1,4 +1,11 @@
-from pyewjn.dielectric.nam_dielectric_coefficient_approximator import get_nam_dielectric
+from pyewjn.dielectric.nam_dielectric_coefficient_approximator import (
+	get_nam_dielectric,
+	get_unapproximated_nam_dielectric,
+)
 from pyewjn.dielectric.lindhard_dielectric import get_lindhard_dielectric
 
-__all__ = ["get_nam_dielectric", "get_lindhard_dielectric"]
+__all__ = [
+	"get_nam_dielectric",
+	"get_lindhard_dielectric",
+	"get_unapproximated_nam_dielectric",
+]

pyewjn/dielectric/nam_dielectric_coefficient_approximator.py

@@ -12,6 +12,7 @@ _logger = logging.getLogger(__name__)
 
 
 FIXED_LARGE_MOMENTUM = 1e8
+WRT_TC_THRESHOLD = 0.98
 
 
 class DedimensionalisedParameters(object):
@@ -24,16 +25,24 @@ class DedimensionalisedParameters(object):
 		temp: float,
 		critical_temp: float,
 		c_light: float,
+		wrt_tc: bool = False,
 	):
 		gap = 0
 		if temp < critical_temp:
 			# else, problems will happen
 			gap = 3.06 * np.sqrt(critical_temp * (critical_temp - temp))
-		self.xi = omega / gap
+
-		self.nu = 1 / (tau * gap)
+		if wrt_tc:
-		self.t = temp / gap
+			scale = critical_temp
-		self.a = omega * v_f / (c_light * gap)
+		else:
+			scale = gap
+
+		self.xi = omega / scale
+		self.nu = 1 / (tau * scale)
+		self.t = temp / scale
+		self.a = omega * v_f / (c_light * scale)
 		self.b = sigma_n / omega
+		self.delta = gap / scale
 
 
 class NamDielectricCoefficients(object):
@@ -65,9 +74,11 @@ def get_dedimensionalised_parameters(
 	temp: float,
 	critical_temp: float,
 	c_light: float,
+	wrt_tc: bool = False,
 ) -> DedimensionalisedParameters:
+
 	return DedimensionalisedParameters(
-		omega, sigma_n, tau, v_f, temp, critical_temp, c_light
+		omega, sigma_n, tau, v_f, temp, critical_temp, c_light, wrt_tc
 	)
 
 
@@ -157,6 +168,7 @@ def get_unapproximated_nam_dielectric(
 		params.t_rel * params.t_c,
 		params.t_c,
 		constants.c_light,
+		wrt_tc=params.t_rel > WRT_TC_THRESHOLD,
 	)
 	prefactor = 4j * np.pi * dedim.b
 

pyewjn/dielectric/sigma_nam_keep_gap.py

@@ -0,0 +1,67 @@
+import numpy as np
+from numpy.lib.scimath import sqrt as csqrt
+
+import pyewjn.util
+
+
+def g(w, wp, d):
+	return ((wp * (w + wp)) + d**2) / (
+		csqrt(wp**2 - d**2) * csqrt((w + wp) ** 2 - d**2)
+	)
+
+
+def s(k, e, v):
+	return (e - 1j * v) / k
+
+
+def f(k, e, v):
+	sv = s(k, e, v)
+	logv = np.log(np.real_if_close((sv + 1) / (sv - 1)) + 0j)
+	return (1 / k) * (2 * sv + ((1 - sv**2) * logv))
+
+
+def i1(w, wp, k, v, d):
+	gv = g(w, wp, d)
+	e1 = csqrt((w + wp) ** 2 - d**2)
+	e2 = csqrt(wp**2 - d**2)
+
+	f_upper = f(k, np.real(e1 - e2), np.imag(e1 + e2) + 2 * v) * (gv + 1)
+	f_lower = f(k, np.real(-e1 - e2), np.imag(e1 + e2) + 2 * v) * (gv - 1)
+
+	return f_upper + f_lower
+
+
+def i2(w, wp, k, v, d):
+	gv = g(w, wp, d)
+	e1 = csqrt((w + wp) ** 2 - d**2)
+	e2 = csqrt(wp**2 - d**2)
+
+	f_upper = f(k, np.real(e1 - e2), np.imag(e1 + e2) + 2 * v) * (gv + 1)
+	f_lower = f(k, np.real(e1 + e2), np.imag(e1 + e2) + 2 * v) * (gv - 1)
+
+	return f_upper + f_lower
+
+
+def a(w, k, v, t, d):
+	result = pyewjn.util.complex_quad(
+		lambda wp: np.tanh((w + wp) / (2 * t)) * (i1(w, wp, k, v, d)),
+		1 - w,
+		1,
+		epsabs=1e-10,
+	)
+
+	return result[0]
+
+
+def b_int(wp, w, k, v, t, d):
+	return (np.tanh((w + wp) / (2 * t)) * i1(w, wp, k, v, d)) - (
+		np.tanh(wp / (2 * t)) * i2(w, wp, k, v, d)
+	)
+
+
+def b(w, k, v, t, d, b_max=np.inf):
+	return pyewjn.util.complex_quad(lambda wp: b_int(wp, w, k, v, t, d), 1, b_max)[0]
+
+
+def sigma_nam_keep_gap(w, k, v, t, d):
+	return -1j * (3 / 4) * (v / w) * (-a(w, k, v, t, d) + b(w, k, v, t, d))

tests/noise/test_chi_nam.py

@@ -0,0 +1,58 @@
+import numpy as np
+import pytest
+
+import pyewjn.dielectric
+import pyewjn.noise.chi
+from pyewjn.baskets import CalculationParams
+
+
+cutoff_to_use = 5.4596e9
+
+
+@pytest.fixture
+def chi_zz_e_nam():
+	params = CalculationParams(
+		omega=1e9, v_f=2e6, omega_p=3.544907701811032e15, tau=1e-14, t_rel=0.99999
+	)
+	eps_l = pyewjn.dielectric.get_nam_dielectric(cutoff_to_use, params)
+	return pyewjn.noise.chi.get_chi_zz_e(eps_l)
+
+
+@pytest.mark.parametrize(
+	"test_input,expected",
+	[
+		# z   chi_zz_e_nam(z)
+		(1e-5, 4.07695673649665e6),
+		(1e-6, 4.095895777068543e9),
+		# (1e-7, 5.012885033150058e12), commenting this one out because it seems numerically too unstable
+		# (1e-8, 1.441261982619894e16), commenting this one out because it seems numerically too unstable
+	],
+)
+def test_chi_zz_e_nam(chi_zz_e_nam, test_input, expected):
+	actual = chi_zz_e_nam(test_input)
+
+	np.testing.assert_allclose(
+		actual,
+		expected,
+		rtol=0.05,
+		err_msg="chi_zz_e is inaccurate for nam case",
+		verbose=True,
+	)
+
+
+@pytest.mark.parametrize(
+	"test_input,expected",
+	[
+		# z   chi_zz_e_nam(z)
+		(1e-6, 4.095895777068543e9),
+	],
+)
+def test_chi_zz_e_nam_benchmark(benchmark, chi_zz_e_nam, test_input, expected):
+	actual = benchmark(chi_zz_e_nam, test_input)
+	np.testing.assert_allclose(
+		actual,
+		expected,
+		rtol=0.05,
+		err_msg="chi_zz_e is inaccurate for nam case",
+		verbose=True,
+	)