Merge branch 'master' into magnoise

pyewjn/dielectric/__init__.py

@@ -3,9 +3,13 @@ from pyewjn.dielectric.nam_dielectric_coefficient_approximator import (
 	get_unapproximated_nam_dielectric,
 )
 from pyewjn.dielectric.lindhard_dielectric import get_lindhard_dielectric
+from pyewjn.dielectric.lindhard_dielectric_transverse import (
+	get_lindhard_dielectric_transverse,
+)
 
 __all__ = [
 	"get_nam_dielectric",
 	"get_lindhard_dielectric",
+	"get_lindhard_dielectric_transverse",
 	"get_unapproximated_nam_dielectric",
 ]

pyewjn/dielectric/lindhard_dielectric_transverse.py

@@ -0,0 +1,70 @@
+import numpy as np
+from pyewjn.baskets import CalculationConstants, CalculationParams
+
+TRANSVERSE_THRESHOLD = 1e4
+
+
+class LindhardDielectricTransverse(object):
+	def __init__(
+		self,
+		params: CalculationParams,
+		constants: CalculationConstants = CalculationConstants(),
+		thres=TRANSVERSE_THRESHOLD,
+	):
+
+		self.series_threshold = thres
+		self.omega = params.omega
+		self.v_f = params.v_f
+		self.omega_p = params.omega_p
+		self.tau = params.tau
+		self.c_light = constants.c_light
+
+		self.s = 1 / (self.tau * self.omega)
+		self.prefactor = 3 * (self.omega_p**2) / (self.omega**2)
+
+	def get_eps(self):
+		def eps_lindhard(u_inverse_wavelength: float) -> complex:
+			"""the lindhard dielectric function
+
+			:param u_inverse_wavelength: u is in units of the reciprocal vacuum wavelength (omega / c_light)
+			:return: returns the value of epsilon, dimensionless
+			"""
+
+			# converts u from inverse vacuum wavelength to inverse mean free path
+			# want to convert to q = vf k, where k is wavevector in SI units
+			q = u_inverse_wavelength * (self.v_f * self.omega) / (self.c_light)
+
+			# if u_inverse_wavelength < self.series_threshold * self.v_f / self.omega:
+			# 	return eps_series(q)
+			# else:
+			# 	return eps_full_lindhard(q)
+			return eps_full_lindhard(q)
+
+		def eps_series(q: float) -> complex:
+			pass
+
+		def eps_full_lindhard(q: float) -> complex:
+
+			return internal_eps_t_full(q, 1 / self.tau, self.omega_p, self.omega)
+
+		return eps_lindhard
+
+
+def internal_eps_t_full(
+	q: float,
+	nu: float,
+	wp: float,
+	w: float,
+):
+	s = nu / w
+	qtw = q / w
+
+	log_val = np.log((1 + qtw + 1j * s) / (1 - qtw + 1j * s))
+	parens = 1 + 1j * s - (((1 + 1j * s) ** 2 - qtw**2) / (2 * qtw)) * log_val
+	return 1 - (3 / 2) * ((wp) / (w)) ** 2 * (1 / (qtw**2)) * parens
+
+
+def get_lindhard_dielectric_transverse(
+	params: CalculationParams, constants: CalculationConstants = CalculationConstants()
+):
+	return LindhardDielectricTransverse(params, constants).get_eps()

pyewjn/noise/chi.py

@@ -19,3 +19,18 @@ def get_chi_zz_e(eps: Callable[[float], complex]) -> Callable[[float], float]:
 		return integral[0] / (z**3)
 
 	return chi_zz_e
+
+
+def get_chi_zz_b(eps_t: Callable[[float], complex]) -> Callable[[float], float]:
+	im_ref_s = pyewjn.noise.im_ref.get_im_ref_s(eps_t)
+
+	def chi_zz_b(z: float) -> float:
+		def integrand(y: float) -> float:
+			return (y**2) * im_ref_s(y / z) * np.exp(-2 * y)
+
+		integral = scipy.integrate.quad(
+			integrand, 0, np.inf, epsabs=1e-10, epsrel=1e-10
+		)
+		return integral[0] / (z**3)
+
+	return chi_zz_b

pyewjn/noise/im_ref.py

@@ -3,6 +3,7 @@ from typing import Callable
 import numpy as np
 
 import pyewjn.noise.zeta
+import pyewjn.util
 
 
 def get_im_ref_p(eps: Callable[[float], complex]) -> Callable[[float], float]:
@@ -13,3 +14,19 @@ def get_im_ref_p(eps: Callable[[float], complex]) -> Callable[[float], float]:
 		return np.imag((np.pi * 1j * u - zeta_p_val) / (np.pi * 1j * u + zeta_p_val))
 
 	return im_ref_p
+
+
+def get_im_ref_s(eps_t: Callable[[float], complex]) -> Callable[[float], float]:
+	def integrand(kappa: float, u: float) -> complex:
+		k_eff2 = kappa**2 + u**2
+		k_eff4 = k_eff2**2
+		k_eff = k_eff2 ** (1 / 2)
+		return eps_t(k_eff) / k_eff4
+
+	def im_ref_s(u: float) -> float:
+		integral = pyewjn.util.complex_quad(
+			lambda kappa: integrand(kappa, u), 0, np.inf, epsabs=1e-12
+		)
+		return (1 / (4 * u**2)) * ((4 * u**3 * integral) - 1)
+
+	return im_ref_s

tests/dielectric/test_lindhard_dielectric_transverse.py

@@ -0,0 +1,27 @@
+import pyewjn.dielectric
+import numpy as np
+import pytest
+from pyewjn.baskets import CalculationParams
+
+
+def get_common_lindhard_dielectric():
+	params = CalculationParams(omega=1e9, omega_p=3.5e15, tau=1e-14, v_f=2e6)
+	return pyewjn.dielectric.get_lindhard_dielectric_transverse(params)
+
+
+@pytest.mark.skip(reason="Not actually correct values")
+@pytest.mark.parametrize(
+	"test_input,expected",
+	[
+		(10, -1222.185185062794 + 1.2249999998777178e8j),
+		(1000, 16924.14814718176 + 1.2250000020552777e8j),
+		(1e8, 83.687499999706 + 0.00022417398943752126j),
+	],
+)
+def test_lindhard_dielectric_transverse(test_input, expected):
+
+	eps_to_test = get_common_lindhard_dielectric()
+
+	np.testing.assert_almost_equal(
+		eps_to_test(test_input), expected, decimal=6, err_msg="b function is off"
+	)