Moved some files and added coefficients

pynam/dielectric/__init__.py

@@ -0,0 +1 @@
+from pynam.dielectric.nam_dielectric_coefficient_approximator import get_nam_dielectric_coefficients

pynam/dielectric/low_k_nam.py

@@ -0,0 +1,55 @@
+import numpy as np
+from numpy.lib.scimath import sqrt as csqrt
+
+import pynam.util.complex_quad
+
+
+def g(w, wp):
+	return ((wp * (w + wp)) + 1) / (csqrt(wp ** 2 - 1) * csqrt((w + wp) ** 2 - 1))
+
+
+def f(k, e, v):
+	return ((4 / 3) * 1 / (e - 1j * v)) + (4 / 15) * (1 / ((e - 1j * v) ** 3)) * k ** 2
+
+
+def i1(w, wp, k, v):
+	gv = g(w, wp)
+	e1 = csqrt((w + wp) ** 2 - 1)
+	e2 = csqrt(wp ** 2 - 1)
+
+	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):
+	gv = g(w, wp)
+	e1 = csqrt((w + wp) ** 2 - 1)
+	e2 = csqrt(wp ** 2 - 1)
+
+	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):
+	return pynam.util.complex_quad.complex_quadrature(
+		lambda wp: np.tanh((w + wp) / (2 * t)) * (i1(w, wp, k, v)),
+		1 - w, 1
+	)[0]
+
+
+def b_int(wp, w, k, v, t):
+	return (np.tanh((w + wp) / (2 * t)) * i1(w, wp, k, v)) - (np.tanh(wp / (2 * t)) * i2(w, wp, k, v))
+
+
+def b(w, k, v, t, b_max=np.inf):
+	return pynam.util.complex_quad.complex_quadrature(
+		lambda wp: b_int(wp, w, k, v, t), 1, b_max
+	)[0]
+
+
+def sigma_nam_alk(w, k, v, t):
+	return -1j * (3 / 4) * (v / w) * (-a(w, k, v, t) + b(w, k, v, t))

pynam/dielectric/nam_dielectric_coefficient_approximator.py

@@ -0,0 +1,105 @@
+import numpy as np
+import pynam.dielectric.sigma_nam
+import pynam.dielectric.low_k_nam
+
+from typing import Tuple
+
+FIXED_LARGE_MOMENTUM = 1e8
+
+
+class DedimensionalisedParameters(object):
+	def __init__(
+			self,
+			omega: float,
+			sigma_n: float,
+			tau: float,
+			v_f: float,
+			temp: float,
+			critical_temp: float,
+			c_light: float):
+		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)
+		self.t = temp / gap
+		self.a = omega * v_f / (c_light * gap)
+		self.b = sigma_n / omega
+
+
+class NamDielectricCoefficients(object):
+	def __init__(self, a: float, b: float, c: float, d: float):
+		self.a = a
+		self.b = b
+		self.c = c
+		self.d = d
+		self.u_l = np.real((-self.c + 1j * self.d) / (-self.a + 1j * self.b))
+
+	def eps(self):
+
+		def piecewise_eps(u: float, u_c: float):
+			# todo add check for u_c vs u_l
+			if u < self.u_l:
+				return -self.a + 1j * self.b
+			elif self.u_l < u < u_c:
+				return 1 + (-self.c + 1j * self.d) / u
+			else:
+				return 1
+
+		return piecewise_eps
+
+
+def get_dedimensionalised_parameters(
+		omega: float,
+		sigma_n: float,
+		tau: float,
+		v_f: float,
+		temp: float,
+		critical_temp: float,
+		c_light: float) -> DedimensionalisedParameters:
+	return DedimensionalisedParameters(omega, sigma_n, tau, v_f, temp, critical_temp, c_light)
+
+
+def get_small_momentum_coefficients(dedim_params: DedimensionalisedParameters) -> Tuple[float, float]:
+	prefactor = 4j * np.pi * dedim_params.b
+	s = pynam.dielectric.low_k_nam.sigma_nam_alk(dedim_params.xi, 0, dedim_params.nu, dedim_params.t)
+	conductivity = prefactor * s
+	return -np.real(conductivity), np.imag(conductivity)
+
+
+def get_big_momentum_coefficients(dedim_params: DedimensionalisedParameters) -> Tuple[float, float]:
+	prefactor = 4j * np.pi * dedim_params.b * FIXED_LARGE_MOMENTUM / dedim_params.a
+	s = pynam.dielectric.sigma_nam.sigma_nam(dedim_params.xi,
+											 FIXED_LARGE_MOMENTUM,
+											 dedim_params.nu,
+											 dedim_params.t)
+	conductivity = prefactor * s
+	return -np.real(conductivity), np.imag(conductivity)
+
+
+def get_nam_dielectric_coefficients(
+		omega: float,
+		sigma_n: float,
+		tau: float,
+		v_f: float,
+		temp: float,
+		crit_temp: float,
+		c_light: float) -> NamDielectricCoefficients:
+	"""Gets a NamDielectricCoefficients object, using SI unit parameters
+
+	:param omega: frequency
+	:param sigma_n: normal state conductivity
+	:param tau: tau in Hz
+	:param v_f: Fermi velocity, in m/s
+	:param temp: temperature in Hz
+	:param crit_temp: critical temperature, in Hz
+	:param c_light: speed of light, meters per second
+	:return:
+	"""
+
+	dedim = get_dedimensionalised_parameters(omega, sigma_n, tau, v_f, temp, crit_temp, c_light)
+	a, b = get_small_momentum_coefficients(dedim)
+	c, d = get_big_momentum_coefficients(dedim)
+
+	return NamDielectricCoefficients(a, b, c, d)

pynam/dielectric/sigma_nam.py

@@ -0,0 +1,64 @@
+import numpy as np
+from numpy.lib.scimath import sqrt as csqrt
+
+import pynam.util
+
+
+def g(w, wp):
+	return ((wp * (w + wp)) + 1) / (csqrt(wp ** 2 - 1) * csqrt((w + wp) ** 2 - 1))
+
+
+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):
+	gv = g(w, wp)
+	e1 = csqrt((w + wp) ** 2 - 1)
+	e2 = csqrt(wp ** 2 - 1)
+
+	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):
+	gv = g(w, wp)
+	e1 = csqrt((w + wp) ** 2 - 1)
+	e2 = csqrt(wp ** 2 - 1)
+
+	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):
+	result = pynam.util.complex_quad.complex_quadrature(
+		lambda wp: np.tanh((w + wp) / (2 * t)) * (i1(w, wp, k, v)),
+		1 - w, 1,
+		epsabs=1e-10
+	)
+
+	return result[0]
+
+
+def b_int(wp, w, k, v, t):
+	return (np.tanh((w + wp) / (2 * t)) * i1(w, wp, k, v)) - (np.tanh(wp / (2 * t)) * i2(w, wp, k, v))
+
+
+def b(w, k, v, t, b_max=np.inf):
+	return pynam.util.complex_quadrature(
+		lambda wp: b_int(wp, w, k, v, t), 1, b_max
+	)[0]
+
+
+def sigma_nam(w, k, v, t):
+	return -1j * (3 / 4) * (v / w) * (-a(w, k, v, t) + b(w, k, v, t))