import pandas as pd
import numpy as np
from scipy.interpolate import interp1d
from loguru import logger
from importlib import resources
from omegaconf import DictConfig, ListConfig
# Physical constants for converting wavelength to angular frequency
C_METERS_PER_SEC = 2.99792458e8 # Speed of light in m/s
HBAR_EV_S = 6.582119569e-16 # Reduced Planck constant in eV*s
[docs]
class DataProvider:
"""
Handles loading and providing material permittivity. It is optimized to create
interpolation functions only once and performs interpolation in frequency
space for better numerical accuracy.
"""
[docs]
def __init__(self, material_config):
self.config = material_config
self.source_type_raw = str(self.config.source_type)
self.source_type = self._normalize_source_type(self.source_type_raw)
logger.info(
f"Initializing DataProvider with source type: '{self.source_type_raw}' "
f"(normalized='{self.source_type}')"
)
if self.source_type == 'excel':
self._setup_interpolator_from_excel()
elif self.source_type == 'constant':
self._setup_constant_epsilon()
elif self.source_type == 'drude':
self._setup_drude_model()
elif self.source_type == 'drude_lorentz':
self._setup_drude_lorentz_model()
elif self.source_type == 'rakic_lorentz_drude':
self._setup_rakic_lorentz_drude_model()
else:
raise ValueError(
f"Unknown material source_type: '{self.source_type_raw}'. "
"Use one of: excel, constant, Drude, Drude-Lorentz, Rakic-Lorentz-Drude."
)
[docs]
@staticmethod
def _normalize_source_type(source_type):
"""Normalize source type values from YAML into internal enum-like strings."""
text = str(source_type).strip().lower().replace("_", "-")
if text in {"excel", "constant", "drude"}:
return text
if text in {"drude-lorentz", "drudelorentz", "drude lorentz"}:
return "drude_lorentz"
if text.replace(" ", "-") == "rakic-lorentz-drude":
return "rakic_lorentz_drude"
return text
[docs]
@staticmethod
def _to_omega_from_ev(value_eV):
"""Convert an energy in eV to angular frequency in rad/s."""
return float(value_eV) / HBAR_EV_S
@staticmethod
def _require_finite(value, key_name):
numeric = float(value)
if not np.isfinite(numeric):
raise ValueError(f"Parameter '{key_name}' must be finite.")
return numeric
@staticmethod
def _require_nonnegative(value, key_name):
numeric = DataProvider._require_finite(value, key_name)
if numeric < 0:
raise ValueError(f"Parameter '{key_name}' must be >= 0.")
return numeric
@staticmethod
def _require_positive(value, key_name):
numeric = DataProvider._require_finite(value, key_name)
if numeric <= 0:
raise ValueError(f"Parameter '{key_name}' must be > 0.")
return numeric
[docs]
def _read_omega_parameter(self, cfg, param_name):
"""Read frequency-like model parameters from either eV or rad/s YAML keys.
Supported keys for a parameter named `omega_p` are:
- `omega_p_eV` (preferred for material-model fits in literature)
- `omega_p_rad_s`
- `omega_p` (legacy fallback interpreted as rad/s)
"""
ev_key = f"{param_name}_eV"
rad_key = f"{param_name}_rad_s"
legacy_key = param_name
has_ev = ev_key in cfg
has_rad = rad_key in cfg
has_legacy = legacy_key in cfg
if sum([has_ev, has_rad, has_legacy]) > 1:
raise ValueError(
f"Specify only one of '{ev_key}', '{rad_key}', or '{legacy_key}' "
"in material config."
)
if has_ev:
value_eV = self._require_nonnegative(cfg[ev_key], ev_key)
return self._to_omega_from_ev(value_eV)
if has_rad:
return self._require_nonnegative(cfg[rad_key], rad_key)
if has_legacy:
return self._require_nonnegative(cfg[legacy_key], legacy_key)
raise ValueError(
f"Missing '{param_name}' in material config. "
f"Provide '{ev_key}', '{rad_key}', or legacy '{legacy_key}'."
)
[docs]
def _setup_interpolator_from_excel(self):
"""Loads data, converts to frequency space, and creates interpolators."""
try:
excel_cfg = self.config.excel_config
filepath = resources.files('mqed') / excel_cfg.filepath
logger.info(f"Loading dispersive material data from: {filepath}")
logger.info(f"Excel sheet name: {excel_cfg.sheet_name}")
lambda_nm, epsilon_real, epsilon_imag = self._read_excel_material_table(
filepath,
excel_cfg.sheet_name,
)
epsilon_complex = epsilon_real + 1j * epsilon_imag
# Convert wavelength to angular frequency (omega)
omega_data = 2 * np.pi * C_METERS_PER_SEC / (lambda_nm * 1e-9)
# Sort data by omega, as interp1d requires monotonically increasing x-values
sort_indices = np.argsort(omega_data)
omega_sorted = omega_data[sort_indices]
epsilon_sorted = epsilon_complex[sort_indices]
self.omega_min = float(omega_sorted[0])
self.omega_max = float(omega_sorted[-1])
# Create and store interpolation functions once
self.interp_real = interp1d(
omega_sorted, epsilon_sorted.real, kind='cubic',
bounds_error=True,
)
self.interp_imag = interp1d(
omega_sorted, epsilon_sorted.imag, kind='cubic',
bounds_error=True,
)
logger.success("Successfully initialized interpolator from Excel data.")
except Exception as e:
logger.error(f"Failed to load or process Excel data: {e}")
raise
@staticmethod
def _read_excel_material_table(filepath, sheet_name):
dataframe = pd.read_excel(filepath, sheet_name=sheet_name, header=None)
if dataframe.shape[1] < 3:
raise ValueError("Excel material sheet must contain at least three columns.")
if dataframe.empty:
raise ValueError("Excel material sheet must contain data rows.")
material_columns = dataframe.iloc[:, :3].copy()
first_row = pd.to_numeric(material_columns.iloc[0], errors="coerce")
if first_row.isna().all():
material_columns = material_columns.iloc[1:].reset_index(drop=True)
elif first_row.isna().any():
raise ValueError(
"Excel material sheet first row must be either three labels or three numeric values."
)
numeric_columns = material_columns.apply(pd.to_numeric, errors="coerce")
if numeric_columns.empty or numeric_columns.isna().any().any():
raise ValueError(
"Excel material wavelength, epsilon_real, and epsilon_imag columns "
"must contain only numeric data."
)
lambda_nm = numeric_columns.iloc[:, 0].to_numpy(dtype=float)
epsilon_real = numeric_columns.iloc[:, 1].to_numpy(dtype=float)
epsilon_imag = numeric_columns.iloc[:, 2].to_numpy(dtype=float)
if not np.all(np.isfinite(lambda_nm)) or np.any(lambda_nm <= 0.0):
raise ValueError("Excel material wavelength values must be finite and positive.")
if np.unique(lambda_nm).size != lambda_nm.size:
raise ValueError("Excel material sheet contains duplicate wavelength values.")
if not np.all(np.isfinite(epsilon_real)) or not np.all(np.isfinite(epsilon_imag)):
raise ValueError("Excel material epsilon values must be finite.")
if lambda_nm.size < 4:
raise ValueError("Excel material sheet needs at least four rows for cubic interpolation.")
return lambda_nm, epsilon_real, epsilon_imag
[docs]
def _setup_constant_epsilon(self):
"""Sets up the provider for a non-dispersive material."""
self.constant_epsilon = complex(self.config.constant_value)
logger.info(f"Using constant, non-dispersive epsilon = {self.constant_epsilon}")
[docs]
def _setup_drude_model(self):
"""Set up single-pole Drude parameters.
Model form used in this code:
epsilon(omega) = eps_inf - omega_p^2 / (omega^2 + i*gamma*omega)
"""
cfg = self.config.drude_config
self.drude_eps_inf = self._require_finite(cfg.eps_inf, "drude_config.eps_inf")
self.drude_omega_p = self._read_omega_parameter(cfg, "omega_p")
self.drude_gamma = self._read_omega_parameter(cfg, "gamma")
self.drude_omega_p = self._require_positive(
self.drude_omega_p,
"drude_config.omega_p",
)
logger.info(
"Using Drude model with eps_inf={}, omega_p={} rad/s, gamma={} rad/s",
self.drude_eps_inf,
self.drude_omega_p,
self.drude_gamma,
)
[docs]
def _setup_drude_lorentz_model(self):
"""Set up validated Drude-Lorentz material parameters."""
cfg = self.config.drude_lorentz_config
self.dl_eps_inf = self._require_finite(cfg.eps_inf, "drude_lorentz_config.eps_inf")
self.dl_omega_p = self._read_omega_parameter(cfg, "omega_p")
self.dl_gamma = self._read_omega_parameter(cfg, "gamma")
self.dl_omega_p = self._require_positive(
self.dl_omega_p,
"drude_lorentz_config.omega_p",
)
oscillators = cfg.oscillators
if not isinstance(oscillators, (list, ListConfig)) or len(oscillators) == 0:
raise ValueError("drude_lorentz_config.oscillators must be a non-empty list.")
self.dl_oscillators = []
for idx, osc in enumerate(oscillators):
if not isinstance(osc, (dict, DictConfig)):
raise ValueError(
f"Oscillator #{idx} must be a dict with strength, omega_0_*, and gamma_*."
)
parsed = {
"strength": self._require_finite(
osc.strength,
f"drude_lorentz_config.oscillators[{idx}].strength",
),
"omega_0": self._read_omega_parameter(osc, "omega_0"),
"gamma": self._read_omega_parameter(osc, "gamma"),
}
parsed["omega_0"] = self._require_positive(
parsed["omega_0"],
f"drude_lorentz_config.oscillators[{idx}].omega_0",
)
self.dl_oscillators.append(parsed)
logger.info(
"Using Drude-Lorentz model with {} oscillator(s)",
len(self.dl_oscillators),
)
def _setup_rakic_lorentz_drude_model(self):
config_name = "rakic_lorentz_drude_config"
if config_name not in self.config:
raise ValueError(f"Missing '{config_name}' in material config.")
cfg = self.config[config_name]
if not isinstance(cfg, (dict, DictConfig)):
raise ValueError(f"'{config_name}' must be a dict.")
self.rakic_eps_inf = self._require_finite(cfg.get("eps_inf", 1), f"{config_name}.eps_inf")
self.rakic_omega_p = self._require_positive(
self._read_omega_parameter(cfg, "omega_p"), f"{config_name}.omega_p"
)
self.rakic_gamma = self._require_nonnegative(
self._read_omega_parameter(cfg, "gamma"), f"{config_name}.gamma"
)
if "drude_strength" not in cfg:
raise ValueError(f"Missing '{config_name}.drude_strength' in material config.")
self.rakic_drude_strength = self._require_nonnegative(
cfg["drude_strength"], f"{config_name}.drude_strength"
)
oscillators = cfg.get("oscillators")
if not isinstance(oscillators, (list, ListConfig)) or not oscillators:
raise ValueError(f"{config_name}.oscillators must be a non-empty list.")
self.rakic_oscillators = []
for idx, osc in enumerate(oscillators):
prefix = f"{config_name}.oscillators[{idx}]"
if not isinstance(osc, (dict, DictConfig)):
raise ValueError(f"{prefix} must be a dict with strength, omega_0_*, and gamma_*.")
if "strength" not in osc:
raise ValueError(f"Missing '{prefix}.strength' in material config.")
self.rakic_oscillators.append({
"strength": self._require_nonnegative(osc["strength"], f"{prefix}.strength"),
"omega_0": self._require_positive(
self._read_omega_parameter(osc, "omega_0"), f"{prefix}.omega_0"
),
"gamma": self._require_nonnegative(
self._read_omega_parameter(osc, "gamma"), f"{prefix}.gamma"
),
})
logger.info(
"Using Rakić Lorentz-Drude model with {} oscillator(s)",
len(self.rakic_oscillators),
)
[docs]
def _epsilon_drude(self, omega):
"""Evaluate single-pole Drude permittivity."""
denominator = omega**2 + 1j * self.drude_gamma * omega
return self.drude_eps_inf - (self.drude_omega_p**2) / denominator
[docs]
def _epsilon_drude_lorentz(self, omega):
"""Evaluate Drude-Lorentz permittivity with configured oscillators."""
epsilon = self.dl_eps_inf - (self.dl_omega_p**2) / (
omega**2 + 1j * self.dl_gamma * omega
)
for osc in self.dl_oscillators:
numerator = osc["strength"] * (osc["omega_0"]**2)
denominator = (osc["omega_0"]**2) - omega**2 - 1j * osc["gamma"] * omega
epsilon = epsilon + numerator / denominator
return epsilon
def _epsilon_rakic_lorentz_drude(self, omega):
plasma_squared = self.rakic_omega_p**2
epsilon = self.rakic_eps_inf - self.rakic_drude_strength * plasma_squared / (
omega**2 + 1j * self.rakic_gamma * omega
)
for osc in self.rakic_oscillators:
epsilon += osc["strength"] * plasma_squared / (
osc["omega_0"]**2 - omega**2 - 1j * osc["gamma"] * omega
)
return epsilon
[docs]
def get_epsilon(self, omega):
"""
Returns the complex relative permittivity (ε_r) for a given angular frequency.
Args:
omega (float): The target angular frequency in rad/s.
Returns:
complex: The complex relative permittivity ε_r(ω).
"""
if self.source_type == 'excel':
omega_eval = float(np.clip(float(omega), self.omega_min, self.omega_max))
real_part = self.interp_real(omega_eval)
imag_part = self.interp_imag(omega_eval)
return complex(real_part, imag_part)
elif self.source_type == 'constant':
return self.constant_epsilon
elif self.source_type == 'drude':
return complex(self._epsilon_drude(float(omega)))
elif self.source_type == 'drude_lorentz':
return complex(self._epsilon_drude_lorentz(float(omega)))
elif self.source_type == 'rakic_lorentz_drude':
return complex(self._epsilon_rakic_lorentz_drude(float(omega)))
raise ValueError(f"Unsupported source type at evaluation time: '{self.source_type_raw}'")