add multi-pole Drude medium

docs/api.rst

@@ -117,6 +117,7 @@ Mediums
    Sellmeier
    Debye
    Lorentz
+   Drude
 
 Methods
 -------

tests/test_components.py

@@ -218,10 +218,11 @@ def test_medium_dispersion():
     m_PR = PoleResidue(eps_inf=1.0, poles=[(1 + 2j, 1 + 3j), (2 + 4j, 1 + 5j)])
     m_SM = Sellmeier(coeffs=[(2, 3), (2, 4)])
     m_LZ = Lorentz(eps_inf=1.0, coeffs=[(1, 3, 2), (2, 4, 1)])
+    m_DR = Drude(eps_inf=1.0, coeffs=[(1, 3), (2, 4)])
     m_DB = Debye(eps_inf=1.0, coeffs=[(1, 3), (2, 4)])
 
     freqs = np.linspace(0.01, 1, 1001)
-    for medium in [m_PR, m_SM, m_LZ, m_DB]:
+    for medium in [m_PR, m_SM, m_LZ, m_DR, m_DB]:
         eps_c = medium.eps_model(freqs)
 
 
@@ -230,15 +231,28 @@ def test_medium_dispersion_conversion():
     m_PR = PoleResidue(eps_inf=1.0, poles=[((1 + 2j), (1 + 3j)), ((2 + 4j), (1 + 5j))])
     m_SM = Sellmeier(coeffs=[(2, 3), (2, 4)])
     m_LZ = Lorentz(eps_inf=1.0, coeffs=[(1, 3, 2), (2, 4, 1)])
+    m_DR = Drude(eps_inf=1.0, coeffs=[(1, 3), (2, 4)])
     m_DB = Debye(eps_inf=1.0, coeffs=[(1, 3), (2, 4)])
 
     freqs = np.linspace(0.01, 1, 1001)
-    for medium in [m_PR, m_SM, m_DB, m_LZ]:  # , m_DB]:
+    for medium in [m_PR, m_SM, m_DB, m_LZ, m_DR]:  # , m_DB]:
         eps_model = medium.eps_model(freqs)
         eps_pr = medium.pole_residue.eps_model(freqs)
         np.testing.assert_allclose(eps_model, eps_pr)
 
 
+def test_medium_dispersion_create():
+
+    m_PR = PoleResidue(eps_inf=1.0, poles=[((1 + 2j), (1 + 3j)), ((2 + 4j), (1 + 5j))])
+    m_SM = Sellmeier(coeffs=[(2, 3), (2, 4)])
+    m_LZ = Lorentz(eps_inf=1.0, coeffs=[(1, 3, 2), (2, 4, 1)])
+    m_DR = Drude(eps_inf=1.0, coeffs=[(1, 3), (2, 4)])
+    m_DB = Debye(eps_inf=1.0, coeffs=[(1, 3), (2, 4)])
+
+    for medium in [m_PR, m_SM, m_DB, m_LZ, m_DR]:
+        struct = Structure(geometry=Box(size=(1, 1, 1)), medium=medium)
+
+
 """ modes """
 
 

tidy3d/__init__.py

@@ -17,7 +17,8 @@
 from .components import Box, Sphere, Cylinder, PolySlab
 
 # medium
-from .components import Medium, PoleResidue, Sellmeier, Debye, Lorentz, AnisotropicMedium, PEC
+from .components import Medium, PoleResidue, AnisotropicMedium, PEC
+from .components import Sellmeier, Debye, Drude, Lorentz
 from .components import nk_to_eps_complex, nk_to_eps_sigma, eps_complex_to_nk
 from .components import nk_to_medium, eps_sigma_to_eps_complex
 

tidy3d/components/__init__.py

@@ -13,7 +13,7 @@
 from .geometry import Geometry
 
 # medium
-from .medium import Medium, PoleResidue, Sellmeier, Debye, Lorentz, AnisotropicMedium, PEC
+from .medium import Medium, PoleResidue, Sellmeier, Debye, Drude, Lorentz, AnisotropicMedium, PEC
 from .medium import nk_to_eps_complex, nk_to_eps_sigma, eps_complex_to_nk
 from .medium import nk_to_medium, eps_sigma_to_eps_complex
 

tidy3d/components/medium.py

@@ -511,6 +511,84 @@ def pole_residue(self):
         )
 
 
+class Drude(DispersiveMedium):
+    """A dispersive medium described by the Drude model.
+    The frequency-dependence of the complex-valued permittivity is described by:
+
+    .. math::
+        \\epsilon(f) = \\epsilon_\\infty - \\sum_i
+        \\frac{ f_i^2}{f^2 - jf\\delta_i}
+
+    where :math:`f, f_i, \\delta_i` are in Hz.
+
+    Parameters
+    ----------
+    eps_inf : float = 1.0
+        Relative permittivity at infinite frequency (:math:`\\epsilon_\\infty`).
+    coeffs : List[Tuple[float, float]]
+        List of (:math:`f_i, \\delta_i`) values for model.
+    frequeuncy_range : Tuple[float, float] = (-inf, inf)
+        Range of validity for the medium in Hz.
+        If simulation or plotting functions use frequency out of this range, a warning is thrown.
+    name : str = None
+        Optional name for the medium.
+
+    Example
+    -------
+    >>> drude_medium = Drude(eps_inf=2.0, coeffs=[(1,2), (3,4)])
+    >>> eps = drude_medium.eps_model(200e12)
+    """
+
+    eps_inf: float = 1.0
+    coeffs: List[Tuple[float, float]]
+    type: Literal["Drude"] = "Drude"
+
+    @ensure_freq_in_range
+    def eps_model(self, frequency: float) -> complex:
+        """Complex-valued permittivity as a function of frequency.
+
+        Parameters
+        ----------
+        frequency : float
+            Frequency to evaluate permittivity at (Hz).
+
+        Returns
+        -------
+        complex
+            Complex-valued relative permittivity evaluated at the frequency.
+        """
+        eps = self.eps_inf + 0.0j
+        for (f, delta) in self.coeffs:
+            eps -= (f ** 2) / (frequency ** 2 - 1j * frequency * delta)
+        return eps
+
+    @property
+    def pole_residue(self):
+        """Representation of Medium as a pole-residue model."""
+
+        poles = []
+        for (f, delta) in self.coeffs:
+
+            w = 2 * np.pi * f
+            d = 2 * np.pi * delta
+
+            c0 = -(w ** 2) / 2 / d + 0j
+            a0 = 0j
+
+            c1 = w ** 2 / 2 / d + 0j
+            a1 = d + 0j
+
+            poles.append((a0, c0))
+            poles.append((a1, c1))
+
+        return PoleResidue(
+            eps_inf=self.eps_inf,
+            poles=poles,
+            frequency_range=self.frequency_range,
+            name=self.name,
+        )
+
+
 class Debye(DispersiveMedium):
     """A dispersive medium described by the Debye model.
     The frequency-dependence of the complex-valued permittivity is described by:
@@ -581,7 +659,9 @@ def pole_residue(self):
 
 
 # types of mediums that can be used in Simulation and Structures
-MediumType = Union[Literal[PEC], Medium, AnisotropicMedium, PoleResidue, Sellmeier, Lorentz, Debye]
+MediumType = Union[
+    Literal[PEC], Medium, AnisotropicMedium, PoleResidue, Sellmeier, Lorentz, Debye, Drude
+]
 
 """ Conversion helper functions """