Coherence in periodic multilayers

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.

The CCL specimen presets follow the original SciApp measurement notebooks: 530 nm high, 50% duty-cycle polymer reliefs on a 50 µm substrate. The periods are 10⁶/740, 10⁶/840 and 10⁶/948 nm. Both the relief and substrate use CCL_grid_substrate.csv, exactly as those notebook models do. The browser evaluates the full table instead of taking every eighth row.

The optional coating follows the three-slice geometry in Diplomka-model.ipynb: coating above the ridges, the relief layer, then coating in the valleys. It uses your chosen original coating dispersion table. This is an exploratory combination, not the aluminium material from that notebook or a fitted reproduction of its measured specimen. No aluminium spectrum is substituted for a polymer coating.

Thin layers remain coherent and include evanescent Fourier orders. A homogeneous thick substrate connects two thin subsystems. The page calculates the coherent response of the complete stack, a Gaussian average of coherent powers over substrate thickness, and an incoherent partial-wave sum. The grating can face either side of the substrate.

The 3D model and optical solver share the same ordered thin-layer definitions. The model shows five periods along x, with grooves uniform along y and depth z increasing downward. Select a region in the model or its button in the layer stack to see physical thickness, depth, material identity, and complex optical properties at the selected wavelength. A material edit uses the same shared control as the main experiment and triggers recalculation; selecting a region alone does not.

Physical scale uses a common scale on all axes. Layer detail compresses thick substrates and enlarges very thin slices; the inspector states the displayed height factor while retaining the physical dimensions. Exploded gaps are purely visual, and translucent exterior boxes are finite portions of semi-infinite media. Colours and opacity identify regions; they do not depict absorption or a solved field. The optional coating consists of horizontal ridge and valley slices, with no sidewall coating. The layer list remains usable if WebGL is unavailable.

All seven original material tables are bundled, with all 1,088 rows each. ANGSTROMS is divided by ten to get nanometres. E1 and E2 are the real and imaginary parts of relative permittivity; they are not refractive index and extinction coefficient, and are not squared. Original files can be downloaded in the material inspector.

The original Python uses ε₁ − iε₂. This implementation uses the conjugate convention ε₁ + iε₂ with exp(−iωt), so positive ε₂ produces attenuation. At tabulated wavelengths the values are unchanged; intermediate wavelengths use linear interpolation. The tables cover 193–1700 nm; this playground exposes 400–1700 nm and never extrapolates. Materials are re-evaluated at every wavelength in a sweep.

The original CCL substrate has ε₂ = 0 in this range. Its zero absorption comes from the data. Complex grating and coating spectra retain their absorption. The thick-layer selector offers the two transparent substrate datasets and air. Lossy thick layers, anisotropy and polarization mixing are not implemented here.

Some long-wavelength entries in CCL_grid_lak.csv have negative fitted ε₂. The original file is preserved, but those settings are unsupported in this passive solver. They produce an explicit error or a gap in the spectrum rather than silently clipping ε₂ to zero. The later substrate and coating fits are separate choices; they are never silently substituted for the originals.

The total absorptance is inferred from the exterior flux deficit and is not an independent conservation check. The volume view separately reconstructs local electromagnetic dissipation; no temperature or thermal-transport model is included. The source files are fitted optical spectra supplied with the notebooks, not raw ellipsometer observations.

Chapter 4 of the thesis describes coherence-transforming matrices and recursive partial-wave summation (Ciompa, 2018, ch. 4). The partial-wave approach comes from Postava, Yamaguchi and Kantor’s treatment of thin coherent and thick incoherent layers (Postava et al., 2002) (reading note). For uncoupled TE or TM polarization with fully decorrelated diffraction channels, each amplitude matrix becomes a power-transfer matrix, including the ratio of outgoing to incoming normal admittance. The total downward population is the solution of a geometric series:

A and B are the two thin subsystems; propagation in the transparent thick layer has unit intensity transmission. The modal populations include the complete partial-wave series; the chart displays only its initial contributions. Evanescent modes are retained inside thin subsystems and excluded from transport through the thick layer, as in the 2024 paper (Postava et al., 2024, §3.1).

The 2024 co-authored paper refines this treatment: orders with the same longitudinal propagation constant can retain mutual coherence (Postava et al., 2024, §2.3) (reading note). The default includes those correlations by propagating a reduced covariance matrix with C → SCS†. At normal incidence this retains the ±m pairs. Turn off the symmetric-order option to compare the original fully decorrelated-order model. The full anisotropic Mueller-matrix formalism is outside this scalar TE/TM implementation.

The average samples a Gaussian distribution of thicknesses between d − 3σ and d + 3σ and normalizes its weights. It averages powers, never complex amplitudes. The controls expose σ and sample count separately. At σ = 0 this is exactly the coherent solution.

Finite thickness averaging does not necessarily equal complete decoherence: the spread has to be wide enough to decorrelate the phase (Postava et al., 2024, §1). Near-degenerate diffraction channels may need a wide thickness distribution to decorrelate. Increasing only the sample count improves quadrature but does not increase the physical spread. This is thickness inhomogeneity, not a calibrated temporal coherence-length control.

The coherent calculation includes evanescent coupling across the substrate. The incoherent model assumes that coupling is negligible, and warns when an evanescent channel has appreciable amplitude across it. Exact grazing settings are reported as numerical exceptions. Wavelength samples can miss dense fringes in the original 50 µm substrate; use a narrow spectral range or increase sampling. Material dispersion, diffraction and interference can all affect the spectrum, so a smooth curve is not, by itself, a convergence study.

The dispersive browser solver uses complex Fourier convolution matrices, inverse factorization for TM (Li, 1996), thin-segment matrix exponentiation, and stable scattering composition (Li, 1996). Automated tests compare individual reflected and transmitted orders against ten fixtures generated by executing the original SciApp Python function bodies. They include both polarizations, oblique incidence, absorbing coatings, both substrate orientations, and coherent and incoherent results. Two of these fixtures are TE cases using the original 2024 symmetric-order projection.

Additional checks verify every supported material-table knot, analytic Fresnel and Airy slab limits, zero-spread averaging, positive absorption, power conservation in transparent samples, and explicit errors for invalid or unsupported inputs. These tests verify selected numerical cases. Experimental agreement across the parameter domain has not been established. Increase the Fourier harmonics toward the paper’s 31 orders (Postava et al., 2024, §3.1) before interpreting fine spectral features.

The original Python and reference notebooks also remain available in calculationsApp’s incoherent-paper notebooks. The local fixture generator is scripts/rcwa-incoherent-reference.py. Python and SciPy are used only during development; browser calculations use no remote solver.

import { type RCWASettings } from './rcwa';
import { thinLayerSpecs } from './rcwa-incoherent-structure';
import { materialPermittivity, type Material } from './rcwa-materials';
import { admittances, interfaceScattering, dispersiveStack, diagonal, complexRoot } from './rcwa-dispersive';
import { complex, sub, propagate, norm } from './rcwa-math';
import { zeros, eye, solve, times, plus, block, join, type Dense, type Scattering } from './rcwa-2d-math';

/** Scalar, isotropic specialization of Ciompa (2018), §4.3, and
 * SciApp/python/funkce.py:calculate_incoherent_partial_waves.
 * All dimensions in nm. Each thick-layer traversal removes unequal-kz order coherence by default.
 * Thin grating scattering still includes evanescent harmonics and interference. */
export interface IncoherentSettings extends RCWASettings {
  spacerIndex: number; spacerThickness: number; thicknessSpread: number;
  placement: 'substrate' | 'cap'; averageSamples: number;
  ridgeMaterial: Material; gapMaterial: Material; spacerMaterial: Material; coatingMaterial: Material; coatingThickness: number; preserveSymmetry: boolean;
}
export const DEFAULT_INCOHERENT: IncoherentSettings = {
  wavelength: 633, period: 1000000 / 740, thickness: 530, fill: .5, angle: 6,
  ridgeIndex: 1.5, grooveIndex: 1, incidentIndex: 1, substrateIndex: 1,
  harmonics: 5, polarization: 'TE', spacerIndex: 1.5, spacerThickness: 50000,
  thicknessSpread: 100, placement: 'substrate', averageSamples: 41,
  ridgeMaterial: 'CCL_grid_substrate', gapMaterial: 'air', spacerMaterial: 'CCL_grid_substrate', coatingMaterial: 'CCL_grid_coating', coatingThickness: 0, preserveSymmetry: true,
};
export interface Powers { reflection: number; transmission: number; r: number[]; t: number[]; absorption: number }
export interface IncoherentResult {
  settings: IncoherentSettings; coherent: Powers; incoherent: Powers; averaged: Powers;
  orders: { order: number; angleR: number | null; angleT: number | null; angleInside: number | null }[];
  passes: { reflection: number; transmission: number; remaining: number }[];
  warnings: string[];
}
const value = (a: Dense, i: number, j: number) => complex(a.re[i * a.cols + j], a.im[i * a.cols + j]);
const column = (a: Dense, col: number) => block(a, 0, col, a.rows, 1);
const sum = (a: number[]) => a.reduce((x, y) => x + y, 0);
function powers(r: number[], t: number[]): Powers {
  const reflection = sum(r), transmission = sum(t);
  if (![...r, ...t].every(Number.isFinite)) throw new Error('The optical calculation produced non-finite powers.');
  return { r, t, reflection, transmission, absorption: 1 - reflection - transmission };
}
function powerMatrix(a: Dense, input: number[], output: number[]) {
  const p = zeros(a.rows, a.cols);
  for (let i = 0; i < a.rows; i++) for (let j = 0; j < a.cols; j++) p.re[i * a.cols + j] = input[j] > 0 ? norm(value(a, i, j)) * output[i] / input[j] : 0;
  return p;
}
function intensityS(a: Scattering, left: number[], right: number[]): Scattering {
  return { r: powerMatrix(a.r, left, left), t: powerMatrix(a.t, left, right),
    rBack: powerMatrix(a.rBack, right, right), tBack: powerMatrix(a.tBack, right, left) };
}
export function solveIncoherent(input: IncoherentSettings): IncoherentResult {
  const ridge = materialPermittivity(input.ridgeMaterial, input.wavelength), gap = materialPermittivity(input.gapMaterial, input.wavelength);
  const spacer = materialPermittivity(input.spacerMaterial, input.wavelength), coating = input.coatingThickness > 0 ? materialPermittivity(input.coatingMaterial, input.wavelength) : complex(1);
  // Original specimen substrates have E2=0 over the exposed spectral range.
  if (spacer.im !== 0) throw new Error('Choose one of the transparent substrate spectra for the thick layer. Absorbing spectra are supported in the thin grating and coating.');
  const s = { ...input, ridgeIndex: complexRoot(ridge).re, grooveIndex: complexRoot(gap).re, spacerIndex: complexRoot(spacer).re };
  if (!Number.isInteger(s.harmonics) || s.harmonics < 0 || s.harmonics > 15 || !Number.isFinite(s.angle) || Math.abs(s.angle) >= 80 || !Number.isFinite(s.period) || s.period <= 0 || !Number.isFinite(s.fill) || s.fill < 0 || s.fill > 1 || !Number.isFinite(s.thickness) || s.thickness < 0 || !Number.isFinite(s.coatingThickness) || s.coatingThickness < 0 || s.coatingThickness > s.thickness || ![s.incidentIndex, s.substrateIndex].every(v => Number.isFinite(v) && v > 0)) throw new Error('Invalid geometry, angle, exterior index, or Fourier harmonics.');
  if (!['substrate', 'cap'].includes(s.placement)) throw new Error('Choose a substrate or capping layer.');
  if (![s.spacerIndex, s.spacerThickness].every(v => Number.isFinite(v) && v > 0) || !Number.isFinite(s.thicknessSpread) || s.thicknessSpread < 0 || 3 * s.thicknessSpread >= s.spacerThickness) throw new Error('The layer thickness must exceed three times its thickness spread; indices must be positive.');
  if (!Number.isInteger(s.averageSamples) || s.averageSamples < 3 || s.averageSamples > 201) throw new Error('Use 3–201 thickness samples.');
  const n = 2 * s.harmonics + 1, center = s.harmonics, tm = s.polarization === 'TM';
  const lateral = s.incidentIndex * Math.sin(s.angle * Math.PI / 180);
  if (Math.abs(lateral) >= s.spacerIndex) throw new Error('The central wave must propagate in the thick layer.');
  const kx = Array.from({ length: n }, (_, i) => lateral + (i - center) * s.wavelength / s.period);
  if ([s.incidentIndex, s.spacerIndex, s.substrateIndex].some(index => kx.some(k => Math.abs(index ** 2 - k ** 2) < 1e-10))) throw new Error('An order is exactly at grazing cutoff. Move the wavelength or angle slightly.');
  const left = complex(s.incidentIndex ** 2), right = complex(s.substrateIndex ** 2);
  const yl = admittances(kx, left, tm), ym = admittances(kx, spacer, tm), yr = admittances(kx, right, tm);
  const layers = thinLayerSpecs(s).map(layer => ({ thickness: layer.thickness, fill: layer.fill,
    ridge: materialPermittivity(layer.ridge, s.wavelength), gap: materialPermittivity(layer.gap, s.wavelength) }));
  const g = dispersiveStack(kx, s.wavelength, tm, s.placement === 'cap' ? spacer : left, s.placement === 'cap' ? right : spacer, layers);
  const a = s.placement === 'cap' ? interfaceScattering(yl, ym) : g;
  const b = s.placement === 'cap' ? g : interfaceScattering(ym, yr);
  const fl = yl.map(y => y.re), fm = ym.map(y => y.re), fr = yr.map(y => y.re);
  const fromS = (total: Scattering) => powers(Array.from(column(powerMatrix(total.r, fl, fl), center).re), Array.from(column(powerMatrix(total.t, fl, fr), center).re));
  const coherentAt = (thickness: number) => {
    const p = diagonal(kx.map(k => propagate(complexRoot(complex(s.spacerIndex ** 2 - k ** 2)), 2 * Math.PI * thickness / s.wavelength)));
    return fromS(join(join(a, { r: zeros(n), rBack: zeros(n), t: p, tBack: p }), b));
  };
  const coherent = coherentAt(s.spacerThickness);
  const ai = intensityS(a, fl, fm), bi = intensityS(b, fm, fr);
  const first = column(ai.t, center), loop = times(ai.rBack, bi.r);
  // Exclude unreachable trapped channels (e.g. TIR modes of a uniform slab),
  // whose arbitrary population otherwise makes I-Q singular.
  const reachable = new Set<number>();
  first.re.forEach((v, i) => { if (v > 1e-24) reachable.add(i); });
  for (let pass = 0; pass < n; pass++) for (const j of [...reachable]) for (let i = 0; i < n; i++) if (loop.re[i * n + j] > 1e-24) reachable.add(i);
  const indices = [...reachable], system = eye(indices.length), rhs = zeros(indices.length, 1);
  indices.forEach((i, row) => { rhs.re[row] = first.re[i]; indices.forEach((j, col) => { system.re[row * indices.length + col] -= loop.re[i * n + j]; }); });
  const totalDown = zeros(n, 1);
  if (indices.length) { const reduced = solve(system, rhs); indices.forEach((i, row) => { totalDown.re[i] = reduced.re[row]; }); }
  const reflected = plus(column(ai.r, center), times(ai.tBack, times(bi.r, totalDown)));
  const transmitted = times(bi.t, totalDown);
  let incoherent = powers(Array.from(reflected.re), Array.from(transmitted.re));
  let passes = [{ reflection: sum(Array.from(column(ai.r, center).re)), transmission: 0, remaining: sum(Array.from(first.re)) }];
  let down = first;
  for (let i = 0; i < 8; i++) {
    const up = times(bi.r, down), r = times(ai.tBack, up), t = times(bi.t, down);
    down = times(ai.rBack, up);
    passes.push({ reflection: sum(Array.from(r.re)), transmission: sum(Array.from(t.re)), remaining: sum(Array.from(down.re)) });
  }
  if (s.preserveSymmetry && kx.some((k, i) => kx.some((v, j) => i !== j && fm[i] > 0 && fm[j] > 0 && Math.abs(k*k-v*v) < 1e-12))) {
    const symmetric = symmetricSum(a,b,fl,fm,fr,kx,center); incoherent = symmetric.power; passes = symmetric.passes;
  }
  let averaged = coherent;
  if (s.thicknessSpread > 0) {
    const r = Array(n).fill(0), t = Array(n).fill(0); let weights = 0;
    for (let i = 0; i < s.averageSamples; i++) {
      const z = -3 + 6 * (i + .5) / s.averageSamples, weight = Math.exp(-z * z / 2);
      const p = coherentAt(s.spacerThickness + z * s.thicknessSpread); weights += weight;
      for (let j = 0; j < n; j++) { r[j] += weight * p.r[j]; t[j] += weight * p.t[j]; }
    }
    averaged = powers(r.map(v => v / weights), t.map(v => v / weights));
  }
  const warnings: string[] = [];
  if ([coherent, incoherent, averaged].some(p => p.absorption < -1e-6 || p.absorption > 1 + 1e-6)) warnings.push('The passive power bounds are violated. Treat this setting as numerically unreliable.');
  if (!s.preserveSymmetry && Math.abs(s.angle) < 1e-8) warnings.push('At normal incidence, symmetric ±m channels can retain mutual coherence. This calculation currently uses the original thesis’s fully decorrelated-order model.');
  if ([-center - 1, center + 1].some(m => Math.abs(lateral + m * s.wavelength / s.period) < Math.max(s.incidentIndex, s.spacerIndex, s.substrateIndex))) warnings.push('Some propagating orders are omitted. Increase the Fourier harmonics.');
  const decay = kx.filter(k => k * k > s.spacerIndex ** 2).map(k => Math.exp(-2 * Math.PI * Math.sqrt(k * k - s.spacerIndex ** 2) * s.spacerThickness / s.wavelength));
  if (decay.some(v => v > 1e-4)) warnings.push('An evanescent order can cross this layer. The thick-layer incoherent approximation may be inadequate near this cutoff.');
  if (s.thicknessSpread > 0 && 4 * Math.PI * s.spacerIndex / s.wavelength * (6 * s.thicknessSpread / s.averageSamples) > .5) warnings.push('Thickness averaging may be undersampled. Increase the thickness samples and compare.');
  const angle = (k: number, index: number) => Math.abs(k) < index ? Math.asin(k / index) * 180 / Math.PI : null;
  return { settings: { ...s }, coherent, incoherent, averaged, passes, warnings,
    orders: kx.map((k, i) => ({ order: i - center, angleR: angle(k, s.incidentIndex), angleT: angle(k, s.substrateIndex), angleInside: angle(k, s.spacerIndex) })) };
}

/** Retain equal-kz cross-order coherences (2024 paper, symmetric interactions).
 * The state contains only pairs that survive a thick-layer traversal. Thin
 * subsystems transform C as S C S† before projecting onto those pairs again. */
function symmetricSum(a: Scattering, b: Scattering, fl: number[], fm: number[], fr: number[], kx: number[], center: number) {
  const n = fm.length, pairs: [number, number][] = [];
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) if (fm[i] > 0 && fm[j] > 0 && Math.abs(kx[i] ** 2 - kx[j] ** 2) < 1e-12) pairs.push([i, j]);
  const normalized = (s: Dense, left: number[], right: number[]) => {
    const out = zeros(s.rows, s.cols);
    for (let i = 0; i < s.rows; i++) for (let j = 0; j < s.cols; j++) if (left[j] > 0 && right[i] > 0) {
      const f = Math.sqrt(right[i] / left[j]), p = i * s.cols + j; out.re[p] = s.re[p] * f; out.im[p] = s.im[p] * f;
    }
    return out;
  };
  const rTop = normalized(a.rBack, fm, fm), rBottom = normalized(b.r, fm, fm);
  const tTop = normalized(a.tBack, fm, fl), tBottom = normalized(b.t, fm, fr), initial = normalized(a.t, fl, fm);
  const product = (a: ReturnType<typeof complex>, b: ReturnType<typeof complex>) => complex(a.re * b.re + a.im * b.im, a.im * b.re - a.re * b.im);
  const operator = (s: Dense, outputPairs: [number, number][]) => {
    const op = zeros(outputPairs.length, pairs.length);
    outputPairs.forEach(([i, j], row) => pairs.forEach(([k, l], col) => {
      const v = product(value(s, i, k), value(s, j, l)); op.re[row * op.cols + col] = v.re; op.im[row * op.cols + col] = v.im;
    })); return op;
  };
  const first = zeros(pairs.length, 1);
  pairs.forEach(([i,j], row) => { const v = product(value(initial, i, center), value(initial, j, center)); first.re[row] = v.re; first.im[row] = v.im; });
  const bottom = operator(rBottom, pairs), loop = times(operator(rTop, pairs), bottom);
  const reachable = new Set<number>(); first.re.forEach((v,i) => { if (Math.hypot(v,first.im[i]) > 1e-24) reachable.add(i); });
  for (let pass=0;pass<pairs.length;pass++) for(const j of [...reachable]) for(let i=0;i<pairs.length;i++) if(norm(value(loop,i,j))>1e-48) reachable.add(i);
  const indices=[...reachable], lhs=eye(indices.length), rhs=zeros(indices.length,1);
  indices.forEach((i,row)=>{rhs.re[row]=first.re[i];rhs.im[row]=first.im[i];indices.forEach((j,col)=>{const p=row*indices.length+col;lhs.re[p]-=loop.re[i*loop.cols+j];lhs.im[p]-=loop.im[i*loop.cols+j];});});
  const down=zeros(pairs.length,1);
  if(indices.length){const solved=solve(lhs,rhs);indices.forEach((i,row)=>{down.re[i]=solved.re[row];down.im[i]=solved.im[row];});}
  const exteriorPairs=Array.from({length:n},(_,i)=>[i,i] as [number,number]);
  const outR=times(operator(tTop,exteriorPairs),bottom), outT=operator(tBottom,exteriorPairs);
  const direct=column(powerMatrix(a.r,fl,fl),center), r=plus(direct,times(outR,down)), t=times(outT,down);
  const trace=(v:Dense)=>pairs.reduce((total,[i,j],row)=>total+(i===j?v.re[row]:0),0);
  const passes=[{reflection:sum(Array.from(direct.re)),transmission:0,remaining:trace(first)}];
  let current=first;
  for(let i=0;i<8;i++){const r=times(outR,current),t=times(outT,current);current=times(loop,current);passes.push({reflection:sum(Array.from(r.re)),transmission:sum(Array.from(t.re)),remaining:trace(current)});}
  return {power:powers(Array.from(r.re),Array.from(t.re)),passes};
}

import { complex, mul, sub, add, scale, norm, type Complex } from './rcwa-math';
import { zeros, eye, plus, times, solve, block, put, join, transparent, maxNorm, type Dense, type Scattering } from './rcwa-2d-math';
export const divide = (a: Complex, b: Complex) => scale(mul(a, complex(b.re, -b.im)), 1 / norm(b));
export function complexRoot(z: Complex): Complex {
  const r = Math.hypot(z.re, z.im), re = Math.sqrt(Math.max(0, (r + z.re) / 2));
  return complex(re, re > 1e-15 ? z.im / (2 * re) : Math.sqrt(Math.max(0, (r - z.re) / 2)));
}
export const diagonal = (values: Complex[]) => {
  const a = zeros(values.length); values.forEach((v, i) => { a.re[i * a.cols + i] = v.re; a.im[i * a.cols + i] = v.im; }); return a;
};
export function admittances(kx: number[], eps: Complex, tm: boolean) {
  return kx.map(k => { const q = complexRoot(sub(eps, complex(k * k))); return tm ? divide(q, eps) : q; });
}
export function interfaceScattering(left: Complex[], right: Complex[]): Scattering {
  const r = left.map((y, i) => divide(sub(y, right[i]), add(y, right[i])));
  return { r: diagonal(r), rBack: diagonal(r.map(v => scale(v, -1))),
    t: diagonal(left.map((y, i) => divide(scale(y, 2), add(y, right[i])))),
    tBack: diagonal(right.map((y, i) => divide(scale(y, 2), add(left[i], y)))) };
}
export function complexConvolution(n: number, fill: number, ridge: Complex, gap: Complex) {
  const out = zeros(n), contrast = sub(ridge, gap);
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) {
    const m = i - j, c = m === 0 ? fill : fill === 0 || fill === 1 ? 0 : Math.sin(Math.PI * m * fill) / (Math.PI * m);
    const v = add(scale(contrast, c), i === j ? gap : complex()); out.re[i * n + j] = v.re; out.im[i * n + j] = v.im;
  }
  return out;
}
/** Complex Maxwell propagation in F/G ports: exp(iG k0 d), evaluated on
 * a thin segment and doubled as scattering matrices. Handles absorbing media
 * without a complex eigensolver or exponentially growing full-layer transfer. */
function propagateLayer(g: Dense, length: number): Scattering {
  const n = g.rows / 2;
  if (length === 0) return transparent(n);
  let bound = 0;
  for (let i = 0; i < g.rows; i++) { let row = 0; for (let j = 0; j < g.cols; j++) row += Math.hypot(g.re[i * g.cols + j], g.im[i * g.cols + j]); bound = Math.max(bound, row); }
  const doubling = Math.max(0, Math.ceil(Math.log2(Math.max(1, length * bound / .5))));
  const generator = zeros(2 * n), step = length / 2 ** doubling;
  for (let i = 0; i < g.re.length; i++) { generator.re[i] = -g.im[i] * step; generator.im[i] = g.re[i] * step; }
  let exp = eye(2 * n), term = eye(2 * n);
  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; }
    exp = plus(exp, term); if (maxNorm(term) < 1e-16) break;
    if (k === 24) throw new Error('Complex layer propagation did not converge.');
  }
  const basis = zeros(2 * n), inverse = zeros(2 * n);
  for (let i = 0; i < n; i++) {
    basis.re[i * 2 * n + i] = basis.re[i * 2 * n + n + i] = basis.re[(n + i) * 2 * n + i] = 1;
    basis.re[(n + i) * 2 * n + n + i] = -1;
  }
  for (let i = 0; i < basis.re.length; i++) inverse.re[i] = basis.re[i] / 2;
  const t = times(times(inverse, exp), basis), 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)), r = times(inv, ba);
  for (let i = 0; i < r.re.length; i++) { r.re[i] *= -1; r.im[i] *= -1; }
  let s: Scattering = { r, tBack: inv, t: plus(aa, times(ab, r)), rBack: times(ab, inv) };
  for (let i = 0; i < doubling; i++) s = join(s, s);
  return s;
}
export interface DispersiveLayer { thickness: number; fill: number; ridge: Complex; gap: Complex }
export function dispersiveStack(kx: number[], wavelength: number, tm: boolean, left: Complex, right: Complex, layers: DispersiveLayer[]): Scattering {
  return dispersivePortStack(kx,wavelength,tm,admittances(kx,left,tm),admittances(kx,right,tm),layers);
}
/** Algebraic reference admittances permit stable field reconstruction at an
 * internal split without introducing fictitious grazing physical ports. */
export function dispersivePortStack(kx:number[],wavelength:number,tm:boolean,left:Complex[],right:Complex[],layers:DispersiveLayer[]):Scattering {
  const n = kx.length, unit = Array.from({ length: n }, () => complex(1));
  let result = interfaceScattering(left, unit);
  for (const layer of layers) {
    if (layer.thickness === 0) continue;
    const eps = complexConvolution(n, layer.fill, layer.ridge, layer.gap);
    const upper = tm ? solve(complexConvolution(n, layer.fill, divide(complex(1), layer.ridge), divide(complex(1), layer.gap)), eye(n)) : eye(n);
    const lower = tm ? eye(n) : { ...eps, re: eps.re.slice(), im: eps.im.slice() };
    if (tm) {
      const inv = solve(eps, eye(n));
      for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) { lower.re[i * n + j] -= kx[i] * inv.re[i * n + j] * kx[j]; lower.im[i * n + j] -= kx[i] * inv.im[i * n + j] * kx[j]; }
    } else for (let i = 0; i < n; i++) lower.re[i * n + i] -= kx[i] ** 2;
    const g = zeros(2 * n); put(g, upper, 0, n); put(g, lower, n, 0);
    result = join(result, propagateLayer(g, 2 * Math.PI * layer.thickness / wavelength));
  }
  return join(result, interfaceScattering(unit, right));
}