Repository: PyEllips/pyElli (language: python)

Issue:
Implement Kramers-Kronig relations
To supply a wider range of optical models it is necessary to use kkr's, because for some dispersions only the real part is analytically defined.
This relates directly to #55.


Discussion / hints:

Fix patch:
diff --git a/docs/dispersions.rst b/docs/dispersions.rst
index c07509d2..5bf41abd 100644
--- a/docs/dispersions.rst
+++ b/docs/dispersions.rst
@@ -141,6 +141,17 @@ Tauc-Lorentz
    fig.update_yaxes(title="Dielectric function")
    fig.update_layout(title="Tauc-Lorentz dispersion with Eg=2 and A=20, E=2.2 and C=1")
 
+Cody-Lorentz
+============
+.. autoclass:: elli.dispersions.CodyLorentz
+
+.. plotly::
+
+   fig = elli.CodyLorentz().get_dielectric_df().plot(backend="plotly")
+   fig.update_xaxes(title="Wavelength (nm")
+   fig.update_yaxes(title="Dielectric function")
+   fig.update_layout(title="Cody-Lorentz dispersion with default values")
+
 
 Gaussian
 ========
diff --git a/docs/kkr.rst b/docs/kkr.rst
new file mode 100644
index 00000000..7d17742c
--- /dev/null
+++ b/docs/kkr.rst
@@ -0,0 +1,5 @@
+Kramers-Kronig relations
+==============================
+
+.. automodule:: elli.kkr.kkr
+    :members:
diff --git a/docs/modules.rst b/docs/modules.rst
index cd99ede1..d17c649b 100644
--- a/docs/modules.rst
+++ b/docs/modules.rst
@@ -16,3 +16,4 @@ API Reference
    plot
    read_write
    helpers
+   kkr
diff --git a/src/elli/dispersions/__init__.py b/src/elli/dispersions/__init__.py
index 8ff0bfe2..2b7ff127 100644
--- a/src/elli/dispersions/__init__.py
+++ b/src/elli/dispersions/__init__.py
@@ -14,3 +14,4 @@
 from .table_spectraray import TableSpectraRay
 from .tanguy import Tanguy
 from .tauc_lorentz import TaucLorentz
+from .cody_lorentz import CodyLorentz
diff --git a/src/elli/dispersions/cody_lorentz.py b/src/elli/dispersions/cody_lorentz.py
new file mode 100644
index 00000000..46376f4e
--- /dev/null
+++ b/src/elli/dispersions/cody_lorentz.py
@@ -0,0 +1,77 @@
+# Encoding: utf-8
+"""Cody-Lorentz dispersion law. Model by Ferlauto et al."""
+from typing import Dict
+import numpy as np
+import numpy.typing as npt
+from scipy.interpolate import interp1d
+
+from ..utils import conversion_wavelength_energy
+from .base_dispersion import Dispersion
+from ..kkr import im2re_reciprocal
+
+
+class CodyLorentz(Dispersion):
+    """Tauc-Lorentz dispersion law. Model by Ferlauto et al.
+
+    Single parameters:
+        :Eg: Bandgap energy (eV). Defaults to 1.6.
+        :A: Amplitude (eV). Defaults to 100.
+        :Et: Energy at which the Urbach tail starts (eV). Defaults to 1.8.
+        :gamma: Broadening (eV). Defaults to 2.4.
+        :Ep:
+            Distance from bandgap for transition from Cody type absorption
+            to Lorentz type absorption (eV). Defaults to 0.8.
+        :E0: Lorentz resonance energy (eV). Defaults to 3.6.
+        :Eu: Exponential decay of the Urbach tail (eV). Defaults to 0.05.
+
+    Repeated parameters:
+        --
+
+    Output:
+        The Cody lorentz dispersion. Please refer to the references for a full formula.
+
+    References:
+        * Ferlauto et al., J. Appl. Phys. 92, 2424 (2002)
+    """
+
+    single_params_template = {
+        "Eg": 1.6,
+        "A": 100,
+        "Et": 1.8,
+        "gamma": 2.4,
+        "Ep": 0.8,
+        "E0": 3.6,
+        "Eu": 0.05,
+    }
+    rep_params_template: Dict[str, float] = {}
+
+    @staticmethod
+    def eps2(E, Eg, A, Et, gamma, Ep, E0, Eu):
+        """The imaginary part of the cody lorentz dispersion"""
+        # pylint: disable=invalid-name
+        def G(E):
+            return (E - Eg) ** 2 / ((E - Eg) ** 2 + Ep**2)
+
+        def L(E):
+            return A * E0 * gamma * E / ((E**2 - E0**2) ** 2 + gamma**2 * E**2)
+
+        E1 = Et * G(Et) * L(Et)
+
+        # fmt: off
+        return (
+            E1 / E * np.exp((E - Et) / Eu) * np.heaviside(Et - E, 1)
+            + G(E) * L(E) * np.heaviside(E - Et, 0)
+        )
+        # fmt: on
+
+    def dielectric_function(self, lbda: npt.ArrayLike) -> npt.NDArray:
+        energy = conversion_wavelength_energy(lbda)
+
+        lbda_broad = np.linspace(50, 10000, 1000)
+        energy_padded = conversion_wavelength_energy(lbda_broad)
+        eps1 = im2re_reciprocal(
+            CodyLorentz.eps2(energy_padded, **self.single_params), lbda_broad
+        )
+        eps1_interp = interp1d(lbda_broad, eps1)(lbda)
+
+        return eps1_interp + 1j * CodyLorentz.eps2(energy, **self.single_params)
diff --git a/src/elli/kkr/__init__.py b/src/elli/kkr/__init__.py
new file mode 100644
index 00000000..94ec84f5
--- /dev/null
+++ b/src/elli/kkr/__init__.py
@@ -0,0 +1,1 @@
+from .kkr import *
diff --git a/src/elli/kkr/kkr.py b/src/elli/kkr/kkr.py
new file mode 100644
index 00000000..a866b182
--- /dev/null
+++ b/src/elli/kkr/kkr.py
@@ -0,0 +1,214 @@
+r"""Calculate Kramers-Kronig relations according to Maclaurin's formula [1]_.
+Here, the differential formulation is used,
+which means that the transformation always diverges to zero at infinity,
+leaving a free infinity offset.
+This offset is typically referred to as :math:`\epsilon(\infty)` in spectroscopy.
+
+Since the Kramers-Kronig relation integrates over the whole spectrum for each point,
+it is very sensitive to sampling changes and non-zero values.
+To obtain best values the y-axis
+should be zero anywhere outside the integration range.
+However, since this is not possible for every dispersion
+relation the calculation should be done on a much wider
+range than the actual interval used.
+Additionally, the discretisation steps of the x-axis should be kept constant.
+
+For the transformation of the real to imaginary part you need to be especially cautios.
+Typically, in spectroscopy the real part of the dielectric
+function is non-zero throughout the whole spectrum
+and especially :math:`\epsilon(\infty) \ne 0`.
+This makes the Kramers-Kronig transformation virtually impossible in these cases
+as it suffers from major uncertanties.
+Hence, this transformation is included in pyElli only for completeness and for the
+special cases it may be applicable.
+
+.. rubric:: References
+
+.. [1] Ohta and Ishida, Appl. Spectroscopy 42, 952 (1988), https://doi.org/10.1366/0003702884430380
+"""
+# pylint: disable=invalid-name
+from typing import Callable
+import numpy as np
+
+
+def _integrate_im(im: np.ndarray, x: np.ndarray, x_i: float) -> np.ndarray:
+    """Calculate the discrete imaginary sum (integral) for the kkr.
+
+    Args:
+        im (numpy.ndarray): The imaginary values from which to calculate.
+        x (numpy.ndarray): The x-axis on which to calculate.
+        x_i (float): The current point around which to integrate.
+
+    Returns:
+        numpy.ndarray: The integral sum
+    """
+    return np.sum(x * im / (x**2 - x_i**2))
+
+
+def _integrate_im_reciprocal(im: np.ndarray, x: np.ndarray, x_i: float) -> np.ndarray:
+    """Calculate the discrete imaginary sum (integral) for the kkr.
+    This formulation uses an 1/x axis to transform a wavelength axis.
+
+    Args:
+        im (numpy.ndarray): The imaginary values from which to calculate.
+        x (numpy.ndarray): The reciprocal x-axis on which to calculate.
+        x_i (float): The current point around which to integrate.
+
+    Returns:
+        numpy.ndarray: The integral sum
+    """
+    return np.sum(im / (x - x**3 / x_i**2))
+
+
+def _integrate_re(re: np.ndarray, x: np.ndarray, x_i: float) -> np.ndarray:
+    """Calculate the discrete real sum (integral) for the kkr.
+
+    Args:
+        re (numpy.ndarray): The real values from which to calculate.
+        x (numpy.ndarray): The x-axis on which to calculate.
+        x_i (float): The current point around which to integrate.
+
+    Returns:
+        numpy.ndarray: The real sum
+    """
+    return np.sum(x_i * re / (x**2 - x_i**2))
+
+
+def _integrate_re_reciprocal(re: np.ndarray, x: np.ndarray, x_i: float) -> np.ndarray:
+    """Calculate the discrete real sum (integral) for the kkr.
+    This formulation uses an 1/x axis to transform a wavelength axis.
+
+    Args:
+        re (numpy.ndarray): The real values from which to calculate.
+        x (numpy.ndarray): The reciprocal x-axis on which to calculate.
+        x_i (float): The current point around which to integrate.
+
+    Returns:
+        numpy.ndarray: The real sum
+    """
+    return np.sum(re / (x_i - x**2 / x_i))
+
+
+def _calc_kkr(
+    t: np.ndarray,
+    x: np.ndarray,
+    trafo: Callable[[np.ndarray, np.ndarray, float], np.ndarray],
+) -> np.ndarray:
+    """Calculates the kramers-kronig relation
+    according to Maclaurin's formula.
+
+    Args:
+        t (np.ndarray): The y-axis on which to transform.
+        x (np.ndarray): The x-axis on which to transform.
+        trafo (Callable[[np.ndarray, np.ndarray, float], np.ndarray]):
+            The transformation function.
+
+    Raises:
+        ValueError: y and x axis must have the same length.
+
+    Returns:
+        np.ndarray: The kkr transformed y-axis
+    """
+    if len(t) != len(x):
+        raise ValueError(
+            "y- and x-axes arrays must have the same length, "
+            f"but have lengths {len(t)} and {len(x)}."
+        )
+
+    integral = np.zeros(len(t))
+    interval = np.diff(x, prepend=x[1] - x[0])
+    odd_y, odd_x = t[1::2], x[1::2]
+    even_y, even_x = t[::2], x[::2]
+    for i, x_i in enumerate(x):
+        if i % 2 == 0:
+            integral[i] = trafo(odd_y, odd_x, x_i)
+        else:
+            integral[i] = trafo(even_y, even_x, x_i)
+
+    return 4 / np.pi * interval * integral
+
+
+def re2im(re: np.ndarray, x: np.ndarray) -> np.ndarray:
+    r"""Calculates the differential Kramers-Kronig relation from the
+    real to imaginary part
+    according to Maclaurin's formula.
+
+    The underlying formula reads:
+
+    .. math::
+        \Delta \Im(x_i) = \Im(x_i) - \Im(\infty) =
+        \frac{2}{\pi} \int_0^\infty \frac{x_i \Re(x)}{x^2 - x_i^2} dx
+
+    Args:
+        re (numpy.ndarray): The real values to transform.
+        x (numpy.ndarray): The axis on which to transform.
+
+    Returns:
+        numpy.ndarray: The transformed imaginary part.
+    """
+    return _calc_kkr(re, x, _integrate_re)
+
+
+def im2re(im: np.ndarray, x: np.ndarray) -> np.ndarray:
+    r"""Calculates the differential Kramers-Kronig relation from the
+    imaginary to real part
+    according to Maclaurin's formula.
+
+    The underlying formula reads:
+
+    .. math::
+        \Delta \Re(x_i) = \Re(x_i) - \Re(\infty) =
+        \frac{2}{\pi} \int_0^\infty \frac{x \Im(x)}{x^2 - x_i^2} dx
+
+    Args:
+        im (numpy.ndarray): The imaginary values to transform.
+        x (numpy.ndarray): The axis on which to transform.
+
+    Returns:
+        numpy.ndarray: The transformed real part.
+    """
+    return _calc_kkr(im, x, _integrate_im)
+
+
+def re2im_reciprocal(re: np.ndarray, x: np.ndarray) -> np.ndarray:
+    r"""Calculates the differential Kramers-Kronig relation from the
+    real to imaginary part
+    according to Maclaurin's formula.
+    This function assumes a reciprocal x-axis, e.g. wavelength in spectroscopy.
+
+    The underlying formula reads:
+
+    .. math::
+        \Delta \Im(x_i) = \Im(x_i) - \Im(\infty) =
+        \frac{2}{\pi} \int_0^\infty \frac{\Re(x)}{x_i - \frac{x^2}{x_i}} dx
+
+    Args:
+        re (numpy.ndarray): The real values to transform.
+        x (numpy.ndarray): The reciprocal axis on which to transform.
+
+    Returns:
+        numpy.ndarray: The transformed imaginary part.
+    """
+    return _calc_kkr(re, x, _integrate_re_reciprocal)
+
+
+def im2re_reciprocal(im: np.ndarray, x: np.ndarray) -> np.ndarray:
+    r"""Calculates the differential Kramers-Kronig relation from the
+    imaginary to real part
+    according to Maclaurin's formula.
+    This function assumes a reciprocal x-axis, e.g. wavelength in spectroscopy.
+
+     The underlying formula reads:
+
+    .. math::
+        \Delta \Re(x_i) = \Re(x_i) - \Re(\infty) =
+        \frac{2}{\pi} \int_0^\infty \frac{x \Im(x)}{1 - \frac{x^2}{x_i^2}} dx
+
+    Args:
+        im (numpy.ndarray): The imaginary values to transform.
+        x (numpy.ndarray): The reciprocal axis on which to transform.
+
+    Returns:
+        numpy.ndarray: The transformed real part.
+    """
+    return _calc_kkr(im, x, _integrate_im_reciprocal)


Test patch:
diff --git a/tests/test_kkr.py b/tests/test_kkr.py
new file mode 100644
index 00000000..c90100ad
--- /dev/null
+++ b/tests/test_kkr.py
@@ -0,0 +1,52 @@
+"""Test kramers kronig relations"""
+import numpy as np
+from numpy.testing import assert_array_almost_equal
+import elli
+from elli.kkr import im2re, im2re_reciprocal
+
+
+def test_tauc_lorentz():
+    """Test whether the kkr reproduces the analytical expression of Tauc-Lorentz"""
+    lbda = np.linspace(1e-2, 2000, 2000)
+    g = elli.TaucLorentz(Eg=5).add(A=20, E=8, C=5)
+    assert_array_almost_equal(
+        im2re_reciprocal(g.get_dielectric(lbda).imag, lbda),
+        g.get_dielectric(lbda).real,
+        decimal=6,
+    )
+
+
+def test_tauc_lorentz_energy():
+    """Test whether the kkr in non reciprocal formulation reproduces the analyitical expression
+    of Tauc-Lorentz"""
+    energy = np.linspace(0, 10, 5000)
+    amp = 20
+    osc_energy = 5
+    gamma = 0.1
+    lorentz = amp / (osc_energy**2 - energy**2 - 1j * gamma * energy)
+
+    assert_array_almost_equal(
+        im2re(lorentz.imag, energy)[10:], lorentz.real[10:], decimal=2
+    )
+
+
+def test_lorentz():
+    """Test whether the kkr reproduces the analyitical expression of a Lorentz oscillator"""
+    lbda = np.linspace(1e-2, 5000, 5000)
+    g = elli.LorentzEnergy().add(A=20, E=5, gamma=5)
+    assert_array_almost_equal(
+        1 + im2re_reciprocal(g.get_dielectric(lbda).imag, lbda)[:-1000],
+        g.get_dielectric(lbda).real[:-1000],
+        decimal=2,
+    )
+
+
+def test_gauss():
+    """KKR reproduces the analytical epxression of gaussian."""
+    lbda = np.linspace(1e-2, 5000, 5000)
+    g = elli.Gaussian().add(A=10, E=8, sigma=5)
+    assert_array_almost_equal(
+        im2re_reciprocal(g.get_dielectric(lbda).imag, lbda)[:-1000],
+        g.get_dielectric(lbda).real[:-1000],
+        decimal=2,
+    )