Vector diffraction from a 2D pillar lattice

The default result is a three-dimensional volume of reconstructed electromagnetic quantities. Electric and magnetic squared magnitudes, time-averaged Poynting flux, signed normal flux, energy density and dissipated-power density are evaluated separately. Reflected and transmitted mode arrows give exterior diffraction directions and normalized powers; evanescent channels are retained in the field solution but have no propagating arrows.

In the 1D models the volume is uniform along the grooves. The incoherent-layer page reconstructs nominal or finite thickness-averaged fields, with quadratic quantities averaged after each coherent field solution. This field ensemble is distinct from the exact complete-decoherence power limit. Stored energy in transparent dispersive regions uses the Brillouin derivative; absorbing regions are omitted from that observable because the supplied optical tables do not uniquely specify stored energy. The view includes formulas, units, reference normalization, sampling limits and source links.

Unit-cell diagrams are part of the structure editor. Additional complex-field slices and phase diagnostics remain optional result views.

“2D” refers to periodicity in both x and y. A rectangular dielectric pillar sits at the center of each unit cell and extends through one layer of uniform height. Above and below are homogeneous, semi-infinite media. The structure is infinite; the field image shows a three-by-three-cell observation window.

The electric and magnetic fields each have three components. At an arbitrary azimuth, the incident s and p polarizations couple through the structure: a scattered order can contain both. The field image reconstructs Eₓ, Eᵧ, and E_z on a plane outside the patterned layer. The squared-magnitude view adds their squared magnitudes and normalizes by the incident electric-vector magnitude. It is not a map of normal power flux or energy density.

This version uses a monochromatic plane wave, lossless isotropic nonmagnetic materials, constant refractive indices, one patterned layer, and rectangular pillars. The 1D lab also offers finite Gaussian-beam examples; finite beams are not implemented here. The Three.js view includes reconstructed field slices inside the patterned layer.

Polar angle θ is measured from the surface normal. Azimuth φ is measured from +x toward +y. At normal incidence, φ defines the incident polarization basis. Every outgoing order has its own azimuth and polar angle. At exactly normal propagation its azimuth is geometrically undefined; the displayed basis is conventional.

The centered rectangular footprint has analytic two-dimensional Fourier coefficients. This generalizes the one-dimensional Fourier construction in thesis §3.1; the present 2D vector implementation is a browser extension. With normalized sinc(u) = sin(πu)/(πu), each coefficient factors into an x term and a y term:

Retaining m,n = −M…M gives N = (2M + 1)² spatial orders. The convolution matrix couples all retained orders, and the propagation state contains 4N tangential field coefficients:

For the exp(−iωt) convention, a thin segment propagates by exp(iGk₀d). The code evaluates a converged matrix series only after scaling to a thin segment, converts it to a scattering matrix, and repeatedly doubles it using Redheffer composition, the star product and doubling algorithm for identical cascaded segments (Rumpf, 2011, §4.1–4.2). This avoids propagating a thick layer with an exponentially growing transfer matrix, the main criterion for a stable layered-grating recursion (Li, 1996). Exterior s/p waves supply the physical boundary conditions.

The underlying Fourier Maxwell formulation and scattering approach are described in Raymond Rumpf’s RCWA formulation lectures (Rumpf, n.d.) and the RCWA mathematics documentation (rcwa documentation, n.d.). The implementation extends the original Python solver’s Fourier and boundary-matching workflow; its propagation algorithm is a new browser implementation.

The Fourier-convergence discussion in the thesis (p. 41) provides the one-dimensional background. This 2D solver uses direct Fourier products (the Laurent rule). Li derived when products of discontinuous periodic functions need a different factorization (Li, 1996), and later applied the correct rules to crossed gratings, with much faster convergence (Li, 1997); that formulation is not implemented here. Sharp material boundaries, especially at high contrast and for p polarization, can converge slowly. The 1D lab uses an inverse factorization for TM, so a finite-order TM strip calculation is not expected to match that implementation exactly. The TE strip limit and homogeneous-layer limits do match.

The default calculation retains nine spatial orders. Compare 25 and 49 orders before interpreting detailed powers. Neither small power-balance error nor agreement between two coarse truncations establishes convergence across the whole spectrum. The field image has finite pixel resolution and may not resolve every near-field feature.

An exact grazing exterior order makes the boundary basis degenerate. This implementation reports that setting as a numerical exception and asks for a nearby wavelength or angle; it does not silently perturb the physical inputs. Affected spectrum samples appear as gaps and are marked failed in the CSV.

Lines connect sampled wavelengths and can miss narrow resonances. Material indices stay constant during a sweep. The model does not include material absorption, dispersion, multilayer patterning, finite arrays, or fabrication imperfections.

Both views show three by three unit cells. The geometry editor offers physical scale and a rescaled detail mode. The field-results view uses a common physical scale on all axes. Physical coordinates (x, y, z) map to scene coordinates (x, −z, y), so depth increases downward. Translucent boxes depict finite portions of the semi-infinite media, and the gold arrow marks the incident wave direction.

Choose an x–y, x–z, or y–z field plane and move it through the structure. External fields use the reflected and transmitted Fourier coefficients. Inside the layer, the solver splits the scattering system at the requested depth, solves for the forward and backward amplitudes there, and reconstructs all three electric components. The shader displays the sampled complex field and advances its monochromatic phase; it does not run a time-domain pulse simulation.

Two rendering regimes are available. Slices display the complex field on a movable plane, with signed-component phase animation. Volume reconstructs |E|² throughout a 48³ or 72³ grid and ray-marches that grid on the GPU. The threshold hides weaker values, while density and opacity control the visual accumulation along the viewing ray; these are display settings, not physical absorption. The volume intensity is stationary for monochromatic illumination. Colour limits adapt to the current plane or volume, so compare numerical scales when changing settings. Sharp dielectric boundaries and limited Fourier or grid resolution can produce ringing or unresolved detail. Reconstruction runs in separate workers, and the volume is reused while the underlying solution and grid size stay unchanged. The renderer releases its resources on navigation and keeps the 2D plots available when WebGL2 is unavailable.

The automated checks include three independent Python/SciPy references that use a full-layer matrix exponential and direct physical boundary matching instead of scattering doubling. They compare each outgoing order at normal and oblique incidence, including circular polarization and a rectangular lattice.

Additional checks cover Fresnel interfaces at multiple azimuths, both polarizations of uniform thin films, total internal reflection, the 1D TE strip limit, rotation of a square lattice, full electric-vector normalization in uniform space, populated mixed (m,n) orders, and power conservation with thick high-contrast layers. Field-slice checks cover the analytic plane wave inside uniform media, continuity of tangential electric coefficients at both patterned interfaces, agreement with the existing external x–y plots, and the coordinate mapping of all three slice orientations. Volume checks compare voxel values with independently reconstructed x–y slices at multiple depths and verify uniform-space intensity for s, p, and circular polarization.

The reference generator is scripts/rcwa-2d-reference.py, with NumPy and SciPy as development-only dependencies. No Python, scientific package, or remote solver is required by the browser. These checks validate the stated formulation and selected cases; they are not experimental validation of a fabricated structure.

import { zeros, eye, solve, times, plus, put, block, slab, type Dense } from './rcwa-2d-math';

export interface RCWA2DSettings {
  wavelength: number; periodX: number; periodY: number; thickness: number; fillX: number; fillY: number;
  pillarIndex: number; gapIndex: number; incidentIndex: number; substrateIndex: number;
  theta: number; phi: number; polarization: 's' | 'p' | 'diagonal' | 'circular'; harmonics: number;
}
export const DEFAULT_RCWA_2D: RCWA2DSettings = {
  wavelength: 633, periodX: 700, periodY: 700, thickness: 220, fillX: .5, fillY: .5,
  pillarIndex: 2, gapIndex: 1, incidentIndex: 1, substrateIndex: 1.5,
  theta: 0, phi: 0, polarization: 's', harmonics: 1,
};
export interface Order2D {
  m: number; n: number; kx: number; ky: number; reflection: number; transmission: number;
  reflectedS: number; reflectedP: number; transmittedS: number; transmittedP: number;
  thetaR: number | null; thetaT: number | null; azimuth: number;
}
export interface RCWA2DResult {
  settings: RCWA2DSettings; reflection: number; transmission: number; balanceError: number; orders: Order2D[]; warnings: string[];
  reflected: Dense; transmitted: Dense; incident: Dense;
  topForward: Dense; bottomBackward: Dense;
}
const rectCoefficient = (order: number, fill: number) => order === 0 ? fill : fill === 0 || fill === 1 ? 0 : Math.sin(Math.PI * order * fill) / (Math.PI * order);

/** Centered rectangular pillars, using direct (Laurent) Fourier products. */
export function geometry2D(s: RCWA2DSettings) {
  const orders = Array.from({ length: (2 * s.harmonics + 1) ** 2 }, (_, i) => ({ m: i % (2 * s.harmonics + 1) - s.harmonics, n: Math.floor(i / (2 * s.harmonics + 1)) - s.harmonics }));
  const count = orders.length, eps = zeros(count), angle = s.theta * Math.PI / 180, phi = s.phi * Math.PI / 180;
  const kx = orders.map(o => s.incidentIndex * Math.sin(angle) * Math.cos(phi) + o.m * s.wavelength / s.periodX);
  const ky = orders.map(o => s.incidentIndex * Math.sin(angle) * Math.sin(phi) + o.n * s.wavelength / s.periodY);
  for (let i = 0; i < count; i++) for (let j = 0; j < count; j++) eps.re[i * count + j] = (i === j ? s.gapIndex ** 2 : 0) + (s.pillarIndex ** 2 - s.gapIndex ** 2) * rectCoefficient(orders[i].m - orders[j].m, s.fillX) * rectCoefficient(orders[i].n - orders[j].n, s.fillY);
  const inv = solve(eps, eye(count)), g = zeros(4 * count);
  // d[Ex,Ey,Hx,Hy]/d(k0 z) = i G [Ex,Ey,Hx,Hy], exp(-i omega t).
  const set = (r: number, c: number, i: number, j: number, v: number) => { g.re[(r * count + i) * g.cols + c * count + j] = v; };
  for (let i = 0; i < count; i++) for (let j = 0; j < count; j++) {
    const d = +(i === j), e = eps.re[i * count + j], v = inv.re[i * count + j];
    set(0, 2, i, j, kx[i] * v * ky[j]); set(0, 3, i, j, d - kx[i] * v * kx[j]);
    set(1, 2, i, j, ky[i] * v * ky[j] - d); set(1, 3, i, j, -ky[i] * v * kx[j]);
    set(2, 0, i, j, -d * kx[i] * ky[i]); set(2, 1, i, j, d * kx[i] ** 2 - e);
    set(3, 0, i, j, e - d * ky[i] ** 2); set(3, 1, i, j, d * ky[i] * kx[i]);
  }
  return { orders, kx, ky, g, inversePermittivity: inv };
}

/** Physical s/p amplitudes: propagating modes have unit electric-vector magnitude. */
function exterior(kx: number[], ky: number[], index: number, direction: 1 | -1, phi: number) {
  const n = kx.length, e = zeros(2 * n), h = zeros(2 * n), fullE = zeros(3 * n, 2 * n), qReal: number[] = [];
  for (let i = 0; i < n; i++) {
    const rho = Math.hypot(kx[i], ky[i]), q2 = index ** 2 - rho ** 2;
    if (Math.abs(q2) < 1e-11) throw new Error('A diffraction order is exactly at grazing cutoff. Shift the wavelength or angle slightly.');
    const qr = q2 > 0 ? Math.sqrt(q2) : 0, qi = q2 < 0 ? Math.sqrt(-q2) : 0;
    const ux = rho > 1e-12 ? kx[i] / rho : Math.cos(phi), uy = rho > 1e-12 ? ky[i] / rho : Math.sin(phi);
    const set = (a: Dense, r: number, c: number, re: number, im = 0) => { a.re[r * a.cols + c] = re; a.im[r * a.cols + c] = im; };
    set(e, i, i, -uy); set(e, n + i, i, ux);
    set(e, i, n + i, direction * qr * ux / index, direction * qi * ux / index);
    set(e, n + i, n + i, direction * qr * uy / index, direction * qi * uy / index);
    set(h, i, i, -direction * qr * ux, -direction * qi * ux);
    set(h, n + i, i, -direction * qr * uy, -direction * qi * uy);
    set(h, i, n + i, -index * uy); set(h, n + i, n + i, index * ux);
    qReal.push(qr);
    for (let r = 0; r < 2 * n; r++) for (const c of [i, n + i]) set(fullE, r, c, e.re[r * e.cols + c], e.im[r * e.cols + c]);
    set(fullE, 2 * n + i, n + i, -rho / index);
  }
  const a = zeros(2 * n), b = zeros(2 * n);
  for (let i = 0; i < 2 * n; i++) for (let j = 0; j < 2 * n; j++) {
    const other = i < n ? i + n : i - n, sign = i < n ? -1 : 1, pos = i * 2 * n + j, hp = other * 2 * n + j;
    a.re[pos] = (e.re[pos] - sign * h.re[hp]) / 2; a.im[pos] = (e.im[pos] - sign * h.im[hp]) / 2;
    b.re[pos] = (e.re[pos] + sign * h.re[hp]) / 2; b.im[pos] = (e.im[pos] + sign * h.im[hp]) / 2;
  }
  return { a, b, fullE, qReal };
}

export function solveRCWA2D(settings: RCWA2DSettings): RCWA2DResult {
  const s = { ...settings };
  for (const key of ['wavelength', 'periodX', 'periodY', 'pillarIndex', 'gapIndex', 'incidentIndex', 'substrateIndex'] as const) if (!Number.isFinite(s[key]) || s[key] <= 0) throw new Error(`${key} must be positive and finite.`);
  if (!Number.isInteger(s.harmonics) || s.harmonics < 0 || s.harmonics > 3) throw new Error('Use 0–3 harmonics per axis.');
  if (![s.fillX, s.fillY].every(v => Number.isFinite(v) && v >= 0 && v <= 1) || !Number.isFinite(s.thickness) || s.thickness < 0 || !Number.isFinite(s.theta) || Math.abs(s.theta) >= 85 || !Number.isFinite(s.phi)) throw new Error('Invalid geometry or incidence angle.');
  if (!['s', 'p', 'diagonal', 'circular'].includes(s.polarization)) throw new Error('Invalid incident polarization.');
  const { orders: indices, kx, ky, g } = geometry2D(s), count = indices.length, n = 2 * count;
  const phi = s.phi * Math.PI / 180;
  const above = exterior(kx, ky, s.incidentIndex, 1, phi), reflectedBasis = exterior(kx, ky, s.incidentIndex, -1, phi), below = exterior(kx, ky, s.substrateIndex, 1, phi);
  const source = zeros(n, 1), zero = Math.floor(count / 2);
  source.re[zero] = s.polarization === 's' ? 1 : s.polarization === 'p' ? 0 : Math.SQRT1_2;
  source.re[count + zero] = s.polarization === 'p' ? 1 : s.polarization === 'diagonal' ? Math.SQRT1_2 : 0;
  source.im[count + zero] = s.polarization === 'circular' ? Math.SQRT1_2 : 0;
  const layer = slab(g, 2 * Math.PI * s.thickness / s.wavelength);
  const aIn = times(above.a, source), bIn = times(above.b, source), system = zeros(2 * n), rhs = zeros(2 * n, 1);
  put(system, plus(reflectedBasis.b, times(layer.r, reflectedBasis.a), -1), 0, 0);
  put(system, times(layer.tBack, below.b), 0, n, -1);
  put(system, times(layer.t, reflectedBasis.a), n, 0, -1);
  put(system, plus(below.a, times(layer.rBack, below.b), -1), n, n);
  put(rhs, plus(times(layer.r, aIn), bIn, -1), 0, 0); put(rhs, times(layer.t, aIn), n, 0);
  const amplitudes = solve(system, rhs), r = block(amplitudes, 0, 0, n, 1), t = block(amplitudes, n, 0, n, 1);
  const flux = s.incidentIndex * Math.cos(s.theta * Math.PI / 180);
  const power = (a: Dense, i: number, q: number) => (a.re[i] ** 2 + a.im[i] ** 2) * q / flux;
  const orders = indices.map((o, i): Order2D => {
    const reflectedS = power(r, i, reflectedBasis.qReal[i]), reflectedP = power(r, count + i, reflectedBasis.qReal[i]);
    const transmittedS = power(t, i, below.qReal[i]), transmittedP = power(t, count + i, below.qReal[i]);
    return { ...o, kx: kx[i], ky: ky[i], reflectedS, reflectedP, transmittedS, transmittedP, reflection: reflectedS + reflectedP, transmission: transmittedS + transmittedP,
      thetaR: reflectedBasis.qReal[i] > 0 ? Math.asin(Math.hypot(kx[i], ky[i]) / s.incidentIndex) * 180 / Math.PI : null,
      thetaT: below.qReal[i] > 0 ? Math.asin(Math.hypot(kx[i], ky[i]) / s.substrateIndex) * 180 / Math.PI : null,
      azimuth: Math.atan2(ky[i], kx[i]) * 180 / Math.PI };
  });
  const reflection = orders.reduce((sum, o) => sum + o.reflection, 0), transmission = orders.reduce((sum, o) => sum + o.transmission, 0), balanceError = Math.abs(1 - reflection - transmission);
  if (![reflection, transmission].every(Number.isFinite)) throw new Error('The calculation produced non-finite power.');
  const warnings: string[] = [];
  if (balanceError > 1e-6) warnings.push('The power-balance check failed. Treat this result as unreliable.');
  const maxIndex = Math.max(s.incidentIndex, s.substrateIndex), baseX = kx[zero], baseY = ky[zero];
  for (let m = -s.harmonics - 1; m <= s.harmonics + 1; m++) for (let n = -s.harmonics - 1; n <= s.harmonics + 1; n++) {
    if (Math.abs(m) <= s.harmonics && Math.abs(n) <= s.harmonics) continue;
    if (Math.hypot(baseX + m * s.wavelength / s.periodX, baseY + n * s.wavelength / s.periodY) < maxIndex && !warnings.some(w => w.includes('omitted'))) warnings.push('Some propagating orders are omitted. Increase the harmonic count before interpreting total powers.');
  }
  return { settings: s, reflection, transmission, balanceError, orders, warnings,
    topForward: plus(aIn, times(reflectedBasis.a, r)), bottomBackward: times(below.b, t),
    reflected: times(reflectedBasis.fullE, r), transmitted: times(below.fullE, t), incident: times(above.fullE, source) };
}

export interface Field2D { width: number; height: number; real: Float64Array[]; imaginary: Float64Array[]; peak: number; intensityPeak: number }
/** Full electric vector on an x-y plane outside the patterned layer, over 3×3 cells. */
export function sampleField2D(result: RCWA2DResult, side: 'above' | 'below', distance: number, width = 128, height = 128): Field2D {
  const s = result.settings, count = result.orders.length, k0 = 2 * Math.PI / s.wavelength;
  const real = [0, 1, 2].map(() => new Float64Array(width * height)), imaginary = [0, 1, 2].map(() => new Float64Array(width * height));
  const field = side === 'above' ? result.reflected : result.transmitted, index = side === 'above' ? s.incidentIndex : s.substrateIndex;
  for (let i = 0; i < count; i++) {
    const order = result.orders[i], q2 = index ** 2 - order.kx ** 2 - order.ky ** 2;
    const q = Math.sqrt(Math.abs(q2)), damping = q2 < 0 ? Math.exp(-q * k0 * distance) : 1, phaseZ = q2 > 0 ? q * k0 * distance : 0;
    for (let y = 0; y < height; y++) for (let x = 0; x < width; x++) {
      const pos = y * width + x;
      const xy = k0 * (order.kx * s.periodX * (3 * x / (width - 1) - 1.5) + order.ky * s.periodY * (3 * y / (height - 1) - 1.5));
      const cr = damping * Math.cos(xy + phaseZ), ci = damping * Math.sin(xy + phaseZ);
      for (let c = 0; c < 3; c++) {
        const fr = field.re[c * count + i], fi = field.im[c * count + i];
        real[c][pos] += fr * cr - fi * ci; imaginary[c][pos] += fr * ci + fi * cr;
        if (side === 'above' && order.m === 0 && order.n === 0) {
          const ir = result.incident.re[c * count + i], ii = result.incident.im[c * count + i];
          real[c][pos] += ir * Math.cos(xy - phaseZ) - ii * Math.sin(xy - phaseZ);
          imaginary[c][pos] += ir * Math.sin(xy - phaseZ) + ii * Math.cos(xy - phaseZ);
        }
      }
    }
  }
  let peak = 0, intensityPeak = 0;
  for (let i = 0; i < width * height; i++) {
    let intensity = 0;
    for (let c = 0; c < 3; c++) { const v = real[c][i] ** 2 + imaginary[c][i] ** 2; intensity += v; peak = Math.max(peak, Math.sqrt(v)); }
    intensityPeak = Math.max(intensityPeak, intensity);
  }
  return { width, height, real, imaginary, peak, intensityPeak };
}

/** Dense complex matrices. Separate real/imaginary buffers keep the 2D solver fast. */
export interface Dense { rows: number; cols: number; re: Float64Array; im: Float64Array }
export const zeros = (rows: number, cols = rows): Dense => ({ rows, cols, re: new Float64Array(rows * cols), im: new Float64Array(rows * cols) });
export function eye(n: number) { const a = zeros(n); for (let i = 0; i < n; i++) a.re[i * n + i] = 1; return a; }
export function plus(a: Dense, b: Dense, factor = 1) {
  const c = zeros(a.rows, a.cols);
  for (let i = 0; i < c.re.length; i++) { c.re[i] = a.re[i] + factor * b.re[i]; c.im[i] = a.im[i] + factor * b.im[i]; }
  return c;
}
export function times(a: Dense, b: Dense) {
  const c = zeros(a.rows, b.cols);
  for (let i = 0; i < a.rows; i++) for (let k = 0; k < a.cols; k++) {
    const ar = a.re[i * a.cols + k], ai = a.im[i * a.cols + k];
    if (ar === 0 && ai === 0) continue;
    for (let j = 0; j < b.cols; j++) {
      const pos = i * b.cols + j, bp = k * b.cols + j;
      c.re[pos] += ar * b.re[bp] - ai * b.im[bp]; c.im[pos] += ar * b.im[bp] + ai * b.re[bp];
    }
  }
  return c;
}
export function block(a: Dense, row: number, col: number, rows: number, cols = rows) {
  const b = zeros(rows, cols);
  for (let i = 0; i < rows; i++) for (let j = 0; j < cols; j++) {
    b.re[i * cols + j] = a.re[(row + i) * a.cols + col + j]; b.im[i * cols + j] = a.im[(row + i) * a.cols + col + j];
  }
  return b;
}
export function put(a: Dense, b: Dense, row: number, col: number, scale = 1) {
  for (let i = 0; i < b.rows; i++) for (let j = 0; j < b.cols; j++) {
    a.re[(row + i) * a.cols + col + j] = scale * b.re[i * b.cols + j]; a.im[(row + i) * a.cols + col + j] = scale * b.im[i * b.cols + j];
  }
}
export function solve(a0: Dense, b0: Dense) {
  const a = { ...a0, re: a0.re.slice(), im: a0.im.slice() }, b = { ...b0, re: b0.re.slice(), im: b0.im.slice() };
  const n = a.rows, m = b.cols;
  const scales = Array.from({ length: n }, (_, i) => {
    let maximum = 0; for (let j = 0; j < n; j++) maximum = Math.max(maximum, Math.hypot(a.re[i * n + j], a.im[i * n + j])); return maximum;
  });
  for (let k = 0; k < n; k++) {
    let pivot = k;
    for (let i = k + 1; i < n; i++) if (Math.hypot(a.re[i * n + k], a.im[i * n + k]) / scales[i] > Math.hypot(a.re[pivot * n + k], a.im[pivot * n + k]) / scales[pivot]) pivot = i;
    if (!scales[pivot] || Math.hypot(a.re[pivot * n + k], a.im[pivot * n + k]) < 1e-13 * scales[pivot]) throw new Error('A boundary system is singular. Try a nearby wavelength or angle.');
    for (const target of [a, b]) for (let j = 0; j < target.cols; j++) {
      const u = k * target.cols + j, v = pivot * target.cols + j;
      [target.re[u], target.re[v]] = [target.re[v], target.re[u]]; [target.im[u], target.im[v]] = [target.im[v], target.im[u]];
    }
    [scales[k], scales[pivot]] = [scales[pivot], scales[k]];
    const dr = a.re[k * n + k], di = a.im[k * n + k], den = dr * dr + di * di;
    for (let i = k + 1; i < n; i++) {
      const ar = a.re[i * n + k], ai = a.im[i * n + k], fr = (ar * dr + ai * di) / den, fi = (ai * dr - ar * di) / den;
      for (const target of [a, b]) for (let j = target === a ? k + 1 : 0; j < target.cols; j++) {
        const u = i * target.cols + j, v = k * target.cols + j;
        target.re[u] -= fr * target.re[v] - fi * target.im[v]; target.im[u] -= fr * target.im[v] + fi * target.re[v];
      }
      a.re[i * n + k] = a.im[i * n + k] = 0;
    }
  }
  const x = zeros(n, m);
  for (let i = n - 1; i >= 0; i--) for (let j = 0; j < m; j++) {
    let re = b.re[i * m + j], im = b.im[i * m + j];
    for (let k = i + 1; k < n; k++) {
      const ap = i * n + k, xp = k * m + j;
      re -= a.re[ap] * x.re[xp] - a.im[ap] * x.im[xp]; im -= a.re[ap] * x.im[xp] + a.im[ap] * x.re[xp];
    }
    const ar = a.re[i * n + i], ai = a.im[i * n + i], den = ar * ar + ai * ai;
    x.re[i * m + j] = (re * ar + im * ai) / den; x.im[i * m + j] = (im * ar - re * ai) / den;
  }
  return x;
}
export function maxNorm(a: Dense) { let value = 0; for (let i = 0; i < a.re.length; i++) value = Math.max(value, Math.hypot(a.re[i], a.im[i])); return value; }

export interface Scattering { r: Dense; tBack: Dense; t: Dense; rBack: Dense }
export const transparent = (n: number): Scattering => ({ r: zeros(n), tBack: eye(n), t: eye(n), rBack: zeros(n) });

/** Redheffer composition: a lies above b. Inverts no growing propagation factors. */
export function join(a: Scattering, b: Scattering): Scattering {
  const n = a.r.rows;
  const rhs = zeros(n, 2 * n); put(rhs, times(b.r, a.t), 0, 0); put(rhs, b.tBack, 0, n);
  const solved = solve(plus(eye(n), times(b.r, a.rBack), -1), rhs);
  const left = block(solved, 0, 0, n), right = block(solved, 0, n, n);
  return { r: plus(a.r, times(a.tBack, left)), tBack: times(a.tBack, right),
    t: times(b.t, plus(a.t, times(a.rBack, left))), rBack: plus(b.rBack, times(b.t, times(a.rBack, right))) };
}

/** exp(i G length) only for a thin segment, then stable scattering doubling. */
export function slab(g: Dense, length: number): Scattering {
  const size = g.rows, n = size / 2, harmonics = n / 2;
  if (length === 0) return transparent(n);
  let infinityNorm = 0;
  for (let i = 0; i < size; i++) { let sum = 0; for (let j = 0; j < size; j++) sum += Math.abs(g.re[i * size + j]); infinityNorm = Math.max(infinityNorm, sum); }
  const doublings = Math.max(0, Math.ceil(Math.log2(Math.max(1, length * infinityNorm / .5))));
  const step = length / 2 ** doublings, generator = zeros(size);
  for (let i = 0; i < g.re.length; i++) generator.im[i] = g.re[i] * step;
  let propagator = eye(size), term = eye(size);
  for (let k = 1; k <= 24; k++) {
    term = times(term, generator);
    for (let i = 0; i < term.re.length; i++) { term.re[i] /= k; term.im[i] /= k; }
    propagator = plus(propagator, term);
    if (maxNorm(term) < 1e-16) break;
    if (k === 24) throw new Error('Thin-layer propagation did not converge.');
  }
  // Algebraic reference ports: E=a+b, H=J(a-b), J=[0,-I;I,0].
  // These carry a positive/negative normal flux without fictitious grazing modes.
  const c = zeros(size), inverse = zeros(size);
  for (let i = 0; i < n; i++) {
    const j = i < harmonics ? i + harmonics : i - harmonics, sign = i < harmonics ? -1 : 1;
    c.re[i * size + i] = c.re[i * size + n + i] = 1;
    c.re[(n + i) * size + j] = sign; c.re[(n + i) * size + n + j] = -sign;
    inverse.re[i * size + i] = inverse.re[(n + i) * size + i] = .5;
    inverse.re[i * size + n + j] = -.5 * sign; inverse.re[(n + i) * size + n + j] = .5 * sign;
  }
  const t = times(times(inverse, propagator), c);
  const aa = block(t, 0, 0, n), ab = block(t, 0, n, n), ba = block(t, n, 0, n), bb = block(t, n, n, n);
  const inv = solve(bb, eye(n)), reflection = times(inv, ba);
  for (let i = 0; i < reflection.re.length; i++) { reflection.re[i] *= -1; reflection.im[i] *= -1; }
  let s: Scattering = { r: reflection, tBack: inv, t: plus(aa, times(ab, reflection)), rBack: times(ab, inv) };
  for (let i = 0; i < doublings; i++) s = join(s, s);
  return s;
}