construct medium from a transfer function describing the differential…

… equation relating polarization current density to electric field

tests/test_components/test_medium.py

@@ -776,3 +776,20 @@ def create_mediums(n_dataset):
         n_data2 = np.vstack((n_data[0, :, :, :].reshape(1, Ny, Nz, Nf), n_data))
         n_dataset2 = td.ScalarFieldDataArray(n_data2, coords=dict(x=X2, y=Y, z=Z, f=freqs))
         create_mediums(n_dataset=n_dataset2)
+
+
+def test_medium_from_transfer_function():
+    freqs = np.linspace(0.01, 1, 1001)
+    m_DR = td.Drude(eps_inf=1.0, coeffs=[(1.5, 3)])
+    # test from transfer function using Drude model
+    f1 = m_DR.coeffs[0][0]
+    d1 = m_DR.coeffs[0][1]
+    # admitance function in Laplace domain (a/b) modeling Drude differential equation is
+    a = np.array([0, td.EPSILON_0 * (2 * np.pi) ** 2 * f1**2, 0])
+    b = np.array([1, 2 * np.pi * d1, 0])
+
+    m_transfer = td.PoleResidue.from_transfer_function(a, b)
+    assert np.allclose(
+        m_transfer.eps_model(freqs),
+        m_DR.eps_model(freqs),
+    )

tidy3d/components/medium.py

@@ -15,6 +15,7 @@
 import numpy as npo
 import pydantic.v1 as pd
 import xarray as xr
+from scipy import signal
 
 from ..constants import (
     C_0,
@@ -3188,6 +3189,95 @@ def compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldM
 
         return derivative_map
 
+    @classmethod
+    def from_transfer_function(
+        cls, a: np.ndarray, b: np.ndarray, eps_inf: pd.PositiveFloat = 1
+    ) -> PoleResidue:
+        """Construct a PoleResidue model from a material with an admittance function defining the
+        relationship between the electric field and current density in the Laplace domain.
+
+        Parameters
+        ----------
+        eps_inf: pd.PositiveFloat
+            The relative permittivity at infinite frequency.
+        a : np.ndarray
+            Coefficients of the numerator polynomial in decreasing order.
+        b : np.ndarray
+            Coefficients of the denominator polynomial in decreasing order.
+
+        Returns
+        -------
+        :class:`.PoleResidue`
+            The pole residue equivalent.
+
+        Notes
+        -----
+
+            The supplied admittance function relates the electric field to the polarization current density
+            in the Laplace domain and is equivalent to a frequency-dependent complex conductivity
+            :math:`\\sigma(\\omega)`.
+
+            .. math::
+                J_p(s) = Y(s)E(s)
+
+            .. math::
+                Y(s) = \\frac{a_0 + a_1 s + \\dots + a_M s^M}{b_0 + b_1 s + \\dots + b_N s^N}
+
+            An equivalent :class:`.PoleResidue` medium is constructed using an equivalent frequency-dependent
+            complex permittivity defined as
+
+            .. math::
+                \\epsilon(\\omega) = \\epsilon_\\infty - \\frac{1}{\\epsilon_0 s}
+                \\frac{a_0 + a_1 s + \\dots + a_M s^M}{b_0 + b_1 s + \\dots + b_N s^N}.
+        """
+
+        assert not np.any(np.iscomplex(a))
+        assert not np.any(np.iscomplex(b))
+        assert a.ndim == 1
+        assert b.ndim == 1
+
+        # Modify the transfer function defining a complex conductivity to match the complex
+        # frequency-dependent portion of the pole residue model
+        # Divide by -j*omega*epsilon (s*epsilon)
+        b = np.append(b * EPSILON_0, 0)
+        # Compute residues and poles using scipy
+        # TODO can we find a better tolerance?
+        pole_tol = 1  # Frequency range is usually quite large GHz-THz, so tolerance of 1 is relatively accurate.
+        (r, p, k) = signal.residue(a, b, tol=pole_tol, rtype="avg")
+
+        if not len(k) == 0:
+            raise ValidationError("Transfer function is invalid. Direct polynomial term detected.")
+        # Assuming real coefficients for the polynomials, the poles should be real or come as complex conjugate pairs
+        r_filtered = []
+        p_filtered = []
+        for res, pole in zip(list(r), list(p)):
+            # Residue equal to zero interpreted as rational expression was not
+            # in simplest form. So skip this pole.
+            if res == 0:
+                continue
+            # Causal and stability check
+            if np.real(pole) > 0:
+                raise ValidationError("Transfer function is invalid. It is non-causal.")
+            # Check for higher order pole
+            if len(np.argwhere(np.array(p_filtered) == pole)) != 0:
+                raise ValidationError(
+                    "Transfer function is invalid. A higher order pole was detected."
+                )
+            if np.imag(pole) == 0:
+                r_filtered.append(res / 2)
+                p_filtered.append(pole)
+            else:
+                pair_found = len(np.argwhere(np.array(p_filtered) == np.conj(pole))) == 1
+                if not pair_found:
+                    raise ValueError("Failed to find complex-conjugate of pole.")
+                r_filtered.append(res)
+                p_filtered.append(pole)
+
+        pole_residue_from_transfer = PoleResidue(
+            eps_inf=eps_inf, poles=list(zip(p_filtered, r_filtered))
+        )
+        return pole_residue_from_transfer
+
 
 class CustomPoleResidue(CustomDispersiveMedium, PoleResidue):
     """A spatially varying dispersive medium described by the pole-residue pair model.