Adding tests of computed medium epsilon vs expected

tests/test_components.py

@@ -1,3 +1,4 @@
+from typing import Dict
 import pytest
 import numpy as np
 import pydantic
@@ -255,6 +256,74 @@ def test_medium_dispersion_create():
         struct = Structure(geometry=Box(size=(1, 1, 1)), medium=medium)
 
 
+def eps_compare(medium: Medium, expected: Dict, tol: float = 1e-5):
+
+    for freq, val in expected.items():
+        # print(f"Expected: {medium.eps_model(freq)}, got: {val}")
+        assert np.abs(medium.eps_model(freq) - val) < tol
+
+
+def test_epsilon_eval():
+    """Compare epsilon evaluated from a dispersive various models to expected."""
+
+    # Dispersive silver model
+    poles_silver = [
+        (a / HBAR, c / HBAR)
+        for (a, c) in [
+            (-2.502e-2 - 8.626e-3j, 5.987e-1 + 4.195e3j),
+            (-2.021e-1 - 9.407e-1j, -2.211e-1 + 2.680e-1j),
+            (-1.467e1 - 1.338e0j, -4.240e0 + 7.324e2j),
+            (-2.997e-1 - 4.034e0j, 6.391e-1 - 7.186e-2j),
+            (-1.896e0 - 4.808e0j, 1.806e0 + 4.563e0j),
+            (-9.396e0 - 6.477e0j, 1.443e0 - 8.219e1j),
+        ]
+    ]
+
+    material = PoleResidue(poles=poles_silver)
+    expected = {
+        2e14: (-102.18389652032306 + 9.22771912188222j),
+        5e14: (-13.517709933590542 + 0.9384819052893092j),
+    }
+    eps_compare(material, expected)
+
+    # Constant and eps, zero sigma
+    material = Medium(permittivity=1.5 ** 2)
+    expected = {
+        2e14: 2.25,
+        5e14: 2.25,
+    }
+    eps_compare(material, expected)
+
+    # Constant eps and sigma
+    material = Medium(permittivity=1.5 ** 2, conductivity=0.1)
+    expected = {
+        2e14: 2.25 + 8.987552009401353j,
+        5e14: 2.25 + 3.5950208037605416j,
+    }
+    eps_compare(material, expected)
+
+    # Constant n and k at a given frequency
+    material = Medium.from_nk(n=1.5, k=0.1, freq=C_0 / 0.8)
+    expected = {
+        2e14: 2.24 + 0.5621108598392753j,
+        5e14: 2.24 + 0.22484434393571015j,
+    }
+    eps_compare(material, expected)
+
+    # Anisotropic material
+    eps = (1.5, 2.0, 2.3)
+    sig = (0.01, 0.03, 0.015)
+    mediums = [Medium(permittivity=eps[i], conductivity=sig[i]) for i in range(3)]
+    material = AnisotropicMedium(xx=mediums[0], yy=mediums[1], zz=mediums[2])
+
+    eps_diag_2 = material.eps_diagonal(2e14)
+    eps_diag_5 = material.eps_diagonal(5e14)
+    assert np.all(np.array(eps_diag_2) == np.array([medium.eps_model(2e14) for medium in mediums]))
+
+    expected = {2e14: np.mean(eps_diag_2), 5e14: np.mean(eps_diag_5)}
+    eps_compare(material, expected)
+
+
 """ modes """
 
 

tidy3d/__init__.py

@@ -42,7 +42,7 @@
 from .components import DATA_TYPE_MAP, ScalarFieldData, ScalarFieldTimeData
 
 # constants imported as `C_0 = td.C_0` or `td.constants.C_0`
-from .constants import inf, C_0, ETA_0
+from .constants import inf, C_0, ETA_0, HBAR
 
 # plugins typically imported as `from tidy3d.plugins import DispersionFitter`
 from . import plugins

tidy3d/components/medium.py

@@ -11,10 +11,31 @@
 from .viz import add_ax_if_none
 from .validators import validate_name_str
 
-from ..constants import C_0, pec_val
+from ..constants import C_0, pec_val, EPSILON_0
 from ..log import log
 
 
+def ensure_freq_in_range(eps_model: Callable[[float], complex]) -> Callable[[float], complex]:
+    """Decorate ``eps_model`` to log warning if frequency supplied is out of bounds."""
+
+    def _eps_model(self, frequency: float) -> complex:
+        """New eps_model function."""
+
+        # if frequency is none, don't check, return original function
+        if frequency is None or self.frequency_range is None:
+            return eps_model(self, frequency)
+
+        fmin, fmax = self.frequency_range
+        if np.any(frequency < fmin) or np.any(frequency > fmax):
+            log.warning(
+                "frequency passed to `Medium.eps_model()`"
+                f"is outside of `Medium.frequency_range` = {self.frequency_range}"
+            )
+        return eps_model(self, frequency)
+
+    return _eps_model
+
+
 """ Medium Definitions """
 
 
@@ -51,6 +72,25 @@ def eps_model(self, frequency: float) -> complex:
             Complex-valued relative permittivity evaluated at ``frequency``.
         """
 
+    @ensure_freq_in_range
+    def eps_diagonal(self, frequency: float) -> Tuple[complex, complex, complex]:
+        """Main diagonal of the complex-valued permittivity tensor as a function of frequency.
+
+        Parameters
+        ----------
+        frequency : float
+            Frequency to evaluate permittivity at (Hz).
+
+        Returns
+        -------
+        complex
+            The diagonal elements of the relative permittivity tensor evaluated at ``frequency``.
+        """
+
+        # This only needs to be overwritten for anisotropic materials
+        eps = self.eps_model(frequency)
+        return (eps, eps, eps)
+
     @add_ax_if_none
     def plot(self, freqs: float, ax: Ax = None) -> Ax:  # pylint: disable=invalid-name
         """Plot n, k of a :class:`Medium` as a function of frequency.
@@ -141,7 +181,7 @@ def nk_to_eps_sigma(n: float, k: float, freq: float) -> Tuple[float, float]:
         eps_complex = AbstractMedium.nk_to_eps_complex(n, k)
         eps_real, eps_imag = eps_complex.real, eps_complex.imag
         omega = 2 * np.pi * freq
-        sigma = omega * eps_imag
+        sigma = omega * eps_imag * EPSILON_0
         return eps_real, sigma
 
     @staticmethod
@@ -166,28 +206,8 @@ def eps_sigma_to_eps_complex(eps_real: float, sigma: float, freq: float) -> comp
         if freq is None:
             return eps_real
         omega = 2 * np.pi * freq
-        return eps_real + 1j * sigma / omega
-
 
-def ensure_freq_in_range(eps_model: Callable[[float], complex]) -> Callable[[float], complex]:
-    """Decorate ``eps_model`` to log warning if frequency supplied is out of bounds."""
-
-    def _eps_model(self, frequency: float) -> complex:
-        """New eps_model function."""
-
-        # if frequency is none, don't check, return original function
-        if frequency is None or self.frequency_range is None:
-            return eps_model(self, frequency)
-
-        fmin, fmax = self.frequency_range
-        if np.any(frequency < fmin) or np.any(frequency > fmax):
-            log.warning(
-                "frequency passed to `Medium.eps_model()`"
-                f"is outside of `Medium.frequency_range` = {self.frequency_range}"
-            )
-        return eps_model(self, frequency)
-
-    return _eps_model
+        return eps_real + 1j * sigma / omega / EPSILON_0
 
 
 """ Dispersionless Medium """
@@ -337,6 +357,26 @@ def eps_model(self, frequency: float) -> complex:
         Tuple[complex, complex, complex]
             Complex-valued relative permittivity for each component evaluated at ``frequency``.
         """
+        eps_xx = self.xx.eps_model(frequency)
+        eps_yy = self.yy.eps_model(frequency)
+        eps_zz = self.zz.eps_model(frequency)
+        return np.mean((eps_xx, eps_yy, eps_zz))
+
+    @ensure_freq_in_range
+    def eps_diagonal(self, frequency: float) -> Tuple[complex, complex, complex]:
+        """Main diagonal of the complex-valued permittivity tensor as a function of frequency.
+
+        Parameters
+        ----------
+        frequency : float
+            Frequency to evaluate permittivity at (Hz).
+
+        Returns
+        -------
+        complex
+            The diagonal elements of the relative permittivity tensor evaluated at ``frequency``.
+        """
+
         eps_xx = self.xx.eps_model(frequency)
         eps_yy = self.yy.eps_model(frequency)
         eps_zz = self.zz.eps_model(frequency)