from_nk to support arbitrary nk values

Author: weiliangjin2021

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 - `tidy3d.plugins.design` tool to explore user-defined design spaces.
 - `ModeData.dispersion` and `ModeSolverData.dispersion` are calculated together with the group index.
+- A utility function `td.medium_from_nk()` that automatically constructs a dispersivless medium when permittivity>=1, and a single-pole Lorentz medium when permittivity<1.
 
 ### Changed
 

tests/test_components/test_medium.py

@@ -81,6 +81,60 @@ def test_medium_conversions():
     assert np.isclose(k, k_)
 
 
+def test_lorentz_medium_conversions(log_capture):
+    freq = 3.0
+
+    # lossless, eps_r > 1
+    eps_complex = 2 + 0j
+    n, k = td.Lorentz.eps_complex_to_nk(eps_complex)
+    medium = td.Lorentz.from_nk(n, k, freq)
+    assert_log_level(log_capture, "WARNING")
+    eps_model = medium.eps_model(freq)
+    assert np.isclose(eps_complex, eps_model)
+
+    # lossless, eps_r < 1
+    eps_complex = 0.5 + 0j
+    n, k = td.Lorentz.eps_complex_to_nk(eps_complex)
+    medium = td.Lorentz.from_nk(n, k, freq)
+    eps_model = medium.eps_model(freq)
+    assert np.isclose(eps_complex, eps_model)
+
+    # lossy, eps_r < 1
+    eps_complex = 0.5 + 0.1j
+    n, k = td.Lorentz.eps_complex_to_nk(eps_complex)
+    medium = td.Lorentz.from_nk(n, k, freq)
+    eps_model = medium.eps_model(freq)
+    assert np.isclose(eps_complex, eps_model)
+
+    # lossy, eps_r > 1
+    eps_complex = 1.5 + 2j
+    n, k = td.Lorentz.eps_complex_to_nk(eps_complex)
+    medium = td.Lorentz.from_nk(n, k, freq)
+    assert_log_level(log_capture, "WARNING")
+    eps_model = medium.eps_model(freq)
+    assert np.isclose(eps_complex, eps_model)
+
+
+def test_medium_from_nk():
+    freq = 3.0
+
+    # lossy, eps_r < 1
+    eps_complex = 0.5 + 0.1j
+    n, k = td.AbstractMedium.eps_complex_to_nk(eps_complex)
+    medium = td.medium_from_nk(n, k, freq)
+    eps_model = medium.eps_model(freq)
+    assert np.isclose(eps_complex, eps_model)
+    assert medium.type == "Lorentz"
+
+    # lossy, eps_r > 1
+    eps_complex = 1.5 + 2j
+    n, k = td.AbstractMedium.eps_complex_to_nk(eps_complex)
+    medium = td.medium_from_nk(n, k, freq)
+    eps_model = medium.eps_model(freq)
+    assert np.isclose(eps_complex, eps_model)
+    assert medium.type == "Medium"
+
+
 def test_PEC():
 
     _ = td.Structure(geometry=td.Box(size=(1, 1, 1)), medium=td.PEC)

tidy3d/__init__.py

@@ -12,7 +12,7 @@
 
 # medium
 from .components.medium import Medium, PoleResidue, AnisotropicMedium, PEC, PECMedium
-from .components.medium import Medium2D, PEC2D
+from .components.medium import Medium2D, PEC2D, medium_from_nk
 from .components.medium import Sellmeier, Debye, Drude, Lorentz
 from .components.medium import CustomMedium, CustomPoleResidue
 from .components.medium import CustomSellmeier, FullyAnisotropicMedium
@@ -326,4 +326,5 @@ def set_logging_level(level: str) -> None:
     "IndexedDataArray",
     "TriangularGridDataset",
     "TetrahedralGridDataset",
+    "medium_from_nk",
 ]

tidy3d/components/medium.py

@@ -1138,6 +1138,12 @@ def from_nk(cls, n: float, k: float, freq: float, **kwargs):
             medium containing the corresponding ``permittivity`` and ``conductivity``.
         """
         eps, sigma = AbstractMedium.nk_to_eps_sigma(n, k, freq)
+        if eps < 1:
+            raise ValidationError(
+                "Dispersiveless medium must have 'permittivity>=1`. "
+                "Please use 'Lorentz.from_nk()' to covert to a Lorentz medium, or the utility "
+                "function 'td.medium_from_nk()' to automatically return the proper medium type."
+            )
         return cls(permittivity=eps, conductivity=sigma, **kwargs)
 
 
@@ -2998,6 +3004,59 @@ def _all_larger(coeff_a, coeff_b) -> bool:
             return np.all(coeff_a.values > coeff_b.values)
         return coeff_a > coeff_b
 
+    @classmethod
+    def from_nk(cls, n: float, k: float, freq: float, **kwargs):
+        """Convert ``n`` and ``k`` values at frequency ``freq`` to a single-pole Lorentz
+        medium.
+
+        Parameters
+        ----------
+        n : float
+            Real part of refractive index.
+        k : float = 0
+            Imaginary part of refrative index.
+        freq : float
+            Frequency to evaluate permittivity at (Hz).
+
+        Returns
+        -------
+        :class:`Lorentz`
+            Lorentz medium having refractive index n+ik at frequency ``freq``.
+        """
+        eps_complex = AbstractMedium.nk_to_eps_complex(n, k)
+        eps_r, eps_i = eps_complex.real, eps_complex.imag
+        if eps_r >= 1:
+            log.warning(
+                "For 'permittivity>=1', it is more computationally efficient to "
+                "use a dispersiveless medium constructed from 'Medium.from_nk()'."
+            )
+        # first, lossless medium
+        if isclose(eps_i, 0):
+            if eps_r < 1:
+                fp = np.sqrt((eps_r - 1) / (eps_r - 2)) * freq
+                return cls(
+                    eps_inf=1,
+                    coeffs=[
+                        (1, fp, 0),
+                    ],
+                )
+            return cls(
+                eps_inf=1,
+                coeffs=[
+                    ((eps_r - 1) / 2, np.sqrt(2) * freq, 0),
+                ],
+            )
+        # lossy medium
+        alpha = (eps_r - 1) / eps_i
+        delta_p = freq / 2 / (alpha**2 - alpha + 1)
+        fp = np.sqrt((alpha**2 + 1) / (alpha**2 - alpha + 1)) * freq
+        return cls(
+            eps_inf=1,
+            coeffs=[
+                (eps_i, fp, delta_p),
+            ],
+        )
+
 
 class CustomLorentz(CustomDispersiveMedium, Lorentz):
     """A spatially varying dispersive medium described by the Lorentz model.
@@ -4843,3 +4902,27 @@ def is_comp_pec_2d(self, comp: Axis, axis: Axis):
 # types of mediums that can be used in Simulation and Structures
 
 MediumType = Union[MediumType3D, Medium2D]
+
+# Utility function
+def medium_from_nk(n: float, k: float, freq: float, **kwargs) -> Union[Medium, Lorentz]:
+    """Convert ``n`` and ``k`` values at frequency ``freq`` to :class:`Medium` if ``Re[epsilon]>=1``,
+    or :class:`Lorentz` if if ``Re[epsilon]<1``.
+
+    Parameters
+    ----------
+    n : float
+        Real part of refractive index.
+    k : float = 0
+        Imaginary part of refrative index.
+    freq : float
+        Frequency to evaluate permittivity at (Hz).
+
+    Returns
+    -------
+    Union[:class:`Medium`, :class:`Lorentz`]
+        Dispersionless medium or Lorentz medium having refractive index n+ik at frequency ``freq``.
+    """
+    eps_complex = AbstractMedium.nk_to_eps_complex(n, k)
+    if eps_complex.real >= 1:
+        return Medium.from_nk(n, k, freq, **kwargs)
+    return Lorentz.from_nk(n, k, freq, **kwargs)