jwt625/BPM

1import numpy as np2import warnings3 4def slab_mode_source(x, w, n_WG, n0, wavelength, ind_m=0, x0=0):5    """6    Returns the normalized TE mode profile for a symmetric slab waveguide with a lateral shift x0.7    """8    k0 = 2 * np.pi / wavelength9 10    def f_even(beta):11        if beta < n0*k0 or beta > n_WG*k0:12            return None13        inside = n_WG**2 * k0**2 - beta**214        outside = beta**2 - n0**2 * k0**215        if inside <= 0 or outside <= 0:16            return None17        kx = np.sqrt(inside)18        kappa = np.sqrt(outside)19        return kx * np.tan(kx * w / 2) - kappa20 21    def f_odd(beta):22        if beta < n0*k0 or beta > n_WG*k0:23            return None24        inside = n_WG**2 * k0**2 - beta**225        outside = beta**2 - n0**2 * k0**226        if inside <= 0 or outside <= 0:27            return None28        kx = np.sqrt(inside)29        kappa = np.sqrt(outside)30        sin_term = np.sin(kx * w / 2)31        if abs(sin_term) < 1e-12:32            return None33        return - kx * (np.cos(kx * w / 2) / sin_term) - kappa34 35    def valid_even(beta):36        inside = n_WG**2 * k0**2 - beta**237        if inside <= 0:38            return False39        kx = np.sqrt(inside)40        theta = kx * w / 241        m = int(np.floor(2 * theta / np.pi))42        if m % 2 == 0:43            if m == 0 and theta > (np.pi/2 - 0.1):44                return False45            return True46        return False47 48    def valid_odd(beta):49        inside = n_WG**2 * k0**2 - beta**250        if inside <= 0:51            return False52        kx = np.sqrt(inside)53        theta = kx * w / 254        m = int(np.floor(2 * theta / np.pi))55        return (m % 2 == 1)56 57    N = 200058    beta_scan = np.linspace(n0*k0, n_WG*k0, N)59    even_intervals = []60    odd_intervals = []61    f_even_vals = [f_even(b) for b in beta_scan]62    f_odd_vals = [f_odd(b) for b in beta_scan]63    for i in range(N-1):64        if (f_even_vals[i] is not None) and (f_even_vals[i+1] is not None):65            if f_even_vals[i] * f_even_vals[i+1] < 0:66                even_intervals.append((beta_scan[i], beta_scan[i+1]))67        if (f_odd_vals[i] is not None) and (f_odd_vals[i+1] is not None):68            if f_odd_vals[i] * f_odd_vals[i+1] < 0:69                odd_intervals.append((beta_scan[i], beta_scan[i+1]))70                71    def refine_root(f, b_left, b_right):72        for _ in range(50):73            b_mid = 0.5*(b_left+b_right)74            val_mid = f(b_mid)75            if val_mid is None:76                b_right = b_mid77                continue78            if abs(val_mid) < 1e-9:79                return b_mid80            val_left = f(b_left)81            if val_left is None or val_left*val_mid > 0:82                b_left = b_mid83            else:84                b_right = b_mid85        return b_mid86 87    even_roots = []88    for (b_left, b_right) in even_intervals:89        root = refine_root(f_even, b_left, b_right)90        if valid_even(root):91            even_roots.append(root)92    odd_roots = []93    for (b_left, b_right) in odd_intervals:94        root = refine_root(f_odd, b_left, b_right)95        if valid_odd(root):96            odd_roots.append(root)97 98    modes = [("even", r) for r in even_roots] + [("odd", r) for r in odd_roots]99    modes_sorted = sorted(modes, key=lambda tup: tup[1], reverse=True)100    if len(modes_sorted) == 0:101        raise ValueError("No guided slab modes found in [n0*k0, n_WG*k0].")102    if ind_m >= len(modes_sorted):103        warnings.warn(104            f"Requested mode index {ind_m} >= found modes ({len(modes_sorted)}). Using highest mode index {len(modes_sorted)-1}.",105            UserWarning106        )107        ind_m = len(modes_sorted) - 1108 109    parity, beta_chosen = modes_sorted[ind_m]110    inside = n_WG**2 * k0**2 - beta_chosen**2111    outside = beta_chosen**2 - n0**2 * k0**2112    kx = np.sqrt(inside)113    kappa = np.sqrt(outside)114 115    E = np.zeros_like(x, dtype=np.complex128)116    if parity == "even":117        for i, xi in enumerate(x):118            xp = xi - x0119            if abs(xp) <= w/2:120                E[i] = np.cos(kx * xp)121            else:122                E[i] = np.cos(kx * (w/2)) * np.exp(-kappa * (abs(xp)-w/2))123    else:124        for i, xi in enumerate(x):125            xp = xi - x0126            if abs(xp) <= w/2:127                E[i] = np.sin(kx * xp)128            else:129                E[i] = np.sign(xp) * np.sin(kx * (w/2)) * np.exp(-kappa * (abs(xp)-w/2))130    norm = np.sqrt(np.trapz(np.abs(E)**2, x))131    E /= norm132    return E133