"""2D FDTD solver with empirical biological tissue parameters."""
from typing import Callable, Optional
import numpy as np
[docs]
class FDTD2D:
"""
2D FDTD solver (TMz polarisation) with empirical tissue dielectric properties.
"""
[docs]
def __init__(
self,
nx: int,
ny: int,
dx: float,
dy: float,
tissue_type: str = "gray_matter",
frequency: float = 1e9,
pml_thickness: int = 10,
pml_strength: float = 0.1,
):
"""
Parameters
----------
nx, ny : int
Grid dimensions (number of Yee cells).
dx, dy : float
Spatial steps (m).
tissue_type : str
'gray_matter' or 'white_matter' (default 'gray_matter').
frequency : float
Operating frequency (Hz) for dielectric properties (default 1e9).
pml_thickness : int
Number of PML cells on each boundary.
pml_strength : float
PML conductivity scaling factor.
"""
self.nx = nx
self.ny = ny
self.dx = dx
self.dy = dy
# Physical constants
self.eps0 = 8.854187817e-12
self.mu0 = 4e-7 * np.pi
self.c = 1.0 / np.sqrt(self.eps0 * self.mu0)
# Empirical tissue properties (at 1 GHz)
if tissue_type == "gray_matter":
self.eps_r = 52.0
self.sigma = 0.98 # S/m
elif tissue_type == "white_matter":
self.eps_r = 38.0
self.sigma = 0.60 # S/m
else:
raise ValueError(f"Unknown tissue type: {tissue_type}")
self.eps = self.eps_r * self.eps0
self.mu = self.mu0
self.sigma_m = 0.0 # magnetic conductivity (assumed zero)
self.pml_thick = pml_thickness
self.pml_strength = pml_strength
# Field arrays (Yee staggered)
self.Ex = np.zeros((nx, ny))
self.Ey = np.zeros((nx, ny))
self.Hz = np.zeros((nx, ny))
self._init_pml()
self.dt = self._courant_dt()
def _courant_dt(self, safety: float = 0.99) -> float:
inv_dx2 = 1.0 / self.dx**2
inv_dy2 = 1.0 / self.dy**2
dt_max = 1.0 / (self.c * np.sqrt(inv_dx2 + inv_dy2))
return safety * dt_max
def _init_pml(self):
self.pml_sigma = np.zeros((self.nx, self.ny))
pml = self.pml_thick
for i in range(pml):
sigma_val = self.pml_strength * ((pml - i) / pml) ** 3
self.pml_sigma[i, :] += sigma_val
self.pml_sigma[self.nx - 1 - i, :] += sigma_val
self.pml_sigma[:, i] += sigma_val
self.pml_sigma[:, self.ny - 1 - i] += sigma_val
[docs]
def step(
self, source: Optional[Callable[[int, int, float], float]] = None, t: float = 0.0
) -> None:
# Update Hz
for i in range(self.nx - 1):
for j in range(self.ny - 1):
dEx_dy = (self.Ex[i, j + 1] - self.Ex[i, j]) / self.dy
dEy_dx = (self.Ey[i + 1, j] - self.Ey[i, j]) / self.dx
self.Hz[i, j] += (
(self.dt / self.mu) * (dEx_dy - dEy_dx) * (1 - self.pml_sigma[i, j] * self.dt)
)
# Update Ex
for i in range(self.nx):
for j in range(self.ny - 1):
dHz_dy = (self.Hz[i, j] - self.Hz[i, j - 1]) / self.dy if j > 0 else 0.0
self.Ex[i, j] += (self.dt / self.eps) * (dHz_dy - self.sigma * self.Ex[i, j])
# Update Ey
for i in range(self.nx - 1):
for j in range(self.ny):
dHz_dx = (self.Hz[i, j] - self.Hz[i - 1, j]) / self.dx if i > 0 else 0.0
self.Ey[i, j] += (self.dt / self.eps) * (-dHz_dx - self.sigma * self.Ey[i, j])
# Add source
if source is not None:
for i in range(self.nx):
for j in range(self.ny):
src = source(i, j, t)
self.Ex[i, j] += src * self.dt / self.eps
self.Ey[i, j] += src * self.dt / self.eps
[docs]
def get_intensity(self) -> np.ndarray:
Ex_grid = self.Ex[:-1, :]
Ey_grid = self.Ey[:, :-1]
return Ex_grid**2 + Ey_grid**2