203 lines
6.3 KiB
Python
203 lines
6.3 KiB
Python
"""2D Mie scattering analytical solution for a dielectric cylinder (TM polarization).
|
||
|
||
Computes the exact scattered and total fields for a circular dielectric cylinder
|
||
under plane-wave illumination u_inc = exp(i·k0·x).
|
||
|
||
Line-by-line translation of the validated MATLAB reference (result/mie.py).
|
||
"""
|
||
|
||
import numpy as np
|
||
from scipy.special import jv, hankel1
|
||
from typing import Optional, Tuple
|
||
|
||
|
||
def mie_scattered_field(
|
||
points: np.ndarray,
|
||
k0: float,
|
||
eps_r: float,
|
||
radius: float,
|
||
cx: float = 0.5,
|
||
cy: float = 0.5,
|
||
) -> np.ndarray:
|
||
"""Compute the scattered E_z field at arbitrary query points.
|
||
|
||
The scattered field is u_scat = u_total − u_inc, valid both inside and
|
||
outside the cylinder. This matches the FEM scattered-field formulation.
|
||
|
||
Parameters
|
||
----------
|
||
points : (N, 2) np.ndarray — (x, y) coordinates
|
||
k0 : float — vacuum wavenumber
|
||
eps_r : float — relative permittivity
|
||
radius : float — cylinder radius
|
||
cx, cy : float — cylinder centre
|
||
|
||
Returns
|
||
-------
|
||
E_scat : (N,) np.complex128
|
||
"""
|
||
m = np.sqrt(eps_r)
|
||
k1 = k0 * m # wavenumber inside cylinder
|
||
x_size = k0 * radius # size parameter
|
||
|
||
# ── polar coordinates relative to cylinder centre ──
|
||
dx = points[:, 0] - cx
|
||
dy = points[:, 1] - cy
|
||
R = np.sqrt(dx * dx + dy * dy)
|
||
Phi = np.arctan2(dy, dx) # [-π, π], matches MATLAB cart2pol
|
||
|
||
# ── Wiscombe truncation (matches MATLAB round(…)) ──
|
||
N_trunc = int(np.round(x_size + 4.05 * x_size ** (1.0 / 3.0) + 2))
|
||
N_trunc = max(N_trunc, 3)
|
||
|
||
E_scat = np.zeros(len(points), dtype=np.complex128)
|
||
E_int = np.zeros(len(points), dtype=np.complex128)
|
||
|
||
for n in range(-N_trunc, N_trunc + 1):
|
||
# boundary values — matches MATLAB besselj / besselh(…, 1, …)
|
||
J_nx = jv(n, x_size)
|
||
J_nmx = jv(n, k1 * radius)
|
||
H_nx = hankel1(n, x_size)
|
||
|
||
# derivatives via recurrence Z'_n = ½ (Z_{n-1} − Z_{n+1})
|
||
J_nx_p = 0.5 * (jv(n - 1, x_size) - jv(n + 1, x_size))
|
||
J_nmx_p = 0.5 * (jv(n - 1, k1 * radius) - jv(n + 1, k1 * radius))
|
||
H_nx_p = 0.5 * (hankel1(n - 1, x_size) - hankel1(n + 1, x_size))
|
||
|
||
# TM scattering coefficient a_n
|
||
num_a = m * J_nx * J_nmx_p - J_nx_p * J_nmx
|
||
den_a = J_nmx * H_nx_p - m * J_nmx_p * H_nx
|
||
a_n = num_a / den_a
|
||
|
||
# internal coefficient c_n
|
||
num_c = J_nx * H_nx_p - J_nx_p * H_nx # Wronskian (2i/(π x) from theory)
|
||
c_n = num_c / den_a
|
||
|
||
# phase factor iⁿ · exp(i·n·φ)
|
||
phase = (1j) ** n * np.exp(1j * n * Phi)
|
||
|
||
# scattered field (valid outside the cylinder)
|
||
out = R >= radius
|
||
if out.any():
|
||
E_scat[out] += a_n * hankel1(n, k0 * R[out]) * phase[out]
|
||
|
||
# internal total field (valid inside the cylinder)
|
||
inside = R < radius
|
||
if inside.any():
|
||
E_int[inside] += c_n * jv(n, k1 * R[inside]) * phase[inside]
|
||
|
||
# phase reference at cylinder centre (matches MATLAB phase_shift)
|
||
phase_shift = np.exp(1j * k0 * cx)
|
||
E_scat *= phase_shift
|
||
E_int *= phase_shift
|
||
|
||
# ── scattered field inside cylinder = internal total − incident ──
|
||
E_inc = np.exp(1j * k0 * points[:, 0])
|
||
inside = R < radius
|
||
if inside.any():
|
||
E_scat[inside] = E_int[inside] - E_inc[inside]
|
||
|
||
return E_scat
|
||
|
||
|
||
def mie_total_field(
|
||
points: np.ndarray,
|
||
k0: float,
|
||
eps_r: float,
|
||
radius: float,
|
||
cx: float = 0.5,
|
||
cy: float = 0.5,
|
||
) -> np.ndarray:
|
||
"""Compute the total E_z field.
|
||
|
||
Outside: u_inc + u_scat
|
||
Inside: internal field (refracted wave)
|
||
"""
|
||
m = np.sqrt(eps_r)
|
||
k1 = k0 * m
|
||
x_size = k0 * radius
|
||
|
||
dx = points[:, 0] - cx
|
||
dy = points[:, 1] - cy
|
||
R = np.sqrt(dx * dx + dy * dy)
|
||
Phi = np.arctan2(dy, dx)
|
||
|
||
N_trunc = int(np.round(x_size + 4.05 * x_size ** (1.0 / 3.0) + 2))
|
||
N_trunc = max(N_trunc, 3)
|
||
|
||
E_scat = np.zeros(len(points), dtype=np.complex128)
|
||
E_int = np.zeros(len(points), dtype=np.complex128)
|
||
|
||
for n in range(-N_trunc, N_trunc + 1):
|
||
J_nx = jv(n, x_size)
|
||
J_nmx = jv(n, k1 * radius)
|
||
H_nx = hankel1(n, x_size)
|
||
|
||
J_nx_p = 0.5 * (jv(n - 1, x_size) - jv(n + 1, x_size))
|
||
J_nmx_p = 0.5 * (jv(n - 1, k1 * radius) - jv(n + 1, k1 * radius))
|
||
H_nx_p = 0.5 * (hankel1(n - 1, x_size) - hankel1(n + 1, x_size))
|
||
|
||
num_a = m * J_nx * J_nmx_p - J_nx_p * J_nmx
|
||
den_a = J_nmx * H_nx_p - m * J_nmx_p * H_nx
|
||
a_n = num_a / den_a
|
||
|
||
num_c = J_nx * H_nx_p - J_nx_p * H_nx
|
||
c_n = num_c / den_a
|
||
|
||
phase = (1j) ** n * np.exp(1j * n * Phi)
|
||
|
||
out = R >= radius
|
||
if out.any():
|
||
E_scat[out] += a_n * hankel1(n, k0 * R[out]) * phase[out]
|
||
|
||
inside = R < radius
|
||
if inside.any():
|
||
E_int[inside] += c_n * jv(n, k1 * R[inside]) * phase[inside]
|
||
|
||
phase_shift = np.exp(1j * k0 * cx)
|
||
E_scat *= phase_shift
|
||
E_int *= phase_shift
|
||
|
||
E_inc = np.exp(1j * k0 * points[:, 0])
|
||
|
||
E_total = np.zeros(len(points), dtype=np.complex128)
|
||
E_total[R >= radius] = E_inc[R >= radius] + E_scat[R >= radius]
|
||
E_total[R < radius] = E_int[R < radius]
|
||
|
||
return E_total
|
||
|
||
|
||
def mie_grid_solution(
|
||
k0: float,
|
||
eps_r: float,
|
||
radius: float,
|
||
cx: float = 0.5,
|
||
cy: float = 0.5,
|
||
x_range: Tuple[float, float] = (0.0, 1.0),
|
||
y_range: Tuple[float, float] = (0.0, 1.0),
|
||
Nx: int = 400,
|
||
Ny: int = 400,
|
||
) -> dict:
|
||
"""Compute Mie solution on a regular grid (for plotting / visual checks).
|
||
|
||
Returns a dict with keys: X, Y, R, Phi, E_inc, E_scat, E_total.
|
||
"""
|
||
x_vec = np.linspace(x_range[0], x_range[1], Nx)
|
||
y_vec = np.linspace(y_range[0], y_range[1], Ny)
|
||
X, Y = np.meshgrid(x_vec, y_vec)
|
||
|
||
points = np.column_stack([X.ravel(), Y.ravel()])
|
||
dx = points[:, 0] - cx
|
||
dy = points[:, 1] - cy
|
||
R = np.sqrt(dx * dx + dy * dy).reshape(Ny, Nx)
|
||
Phi = np.arctan2(dy, dx).reshape(Ny, Nx)
|
||
|
||
E_inc = np.exp(1j * k0 * X)
|
||
E_scat = mie_scattered_field(points, k0, eps_r, radius, cx, cy).reshape(Ny, Nx)
|
||
E_total = mie_total_field(points, k0, eps_r, radius, cx, cy).reshape(Ny, Nx)
|
||
|
||
return {
|
||
"X": X, "Y": Y, "R": R, "Phi": Phi,
|
||
"E_inc": E_inc, "E_scat": E_scat, "E_total": E_total,
|
||
}
|