Diffraction from a 1D lamellar grating

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.

One rectangular grating layer between two homogeneous, semi-infinite media. The structure is periodic along x, uniform along the grooves (y), and illuminated in the x–z plane. Materials are isotropic, nonmagnetic, lossless, and have constant refractive indices. TE has its electric field along the grooves; TM has its magnetic field along them.

This version does not model absorbing materials, dispersion, multilayer stacks, conical incidence, or incoherence. A zero-height layer reduces to an interface. Equal ridge and gap indices give a uniform thin film.

The optional phase-resolved slice reconstructs the transverse component: Eᵧ for TE and Hᵧ for TM, normalized to the incident component. Its squared magnitude is not the total electric-field intensity or energy density. Animation advances the phase of a monochromatic steady-state solution; it does not simulate a pulse arriving at the grating.

RCWA, also called the Fourier modal method (Li, 2014, §13.1), expands the material and fields into spatial Fourier harmonics, following the formulation in thesis §3.1 (p. 33). The integer m labels the diffraction order; M sets the truncation, retaining 2M + 1 orders. The ridge is centered in each period, so its Fourier coefficients are real:

For TE, the eigenproblem is (E − K²)W = Wq². For TM, inverse Fourier factorization gives (I − KE⁻¹K)W = BWq². This is the TM reformulation reported in 1996 by Lalanne and Morris (Lalanne & Morris, 1996) and by Granet and Guizal (Granet & Guizal, 1996); Li showed that it converges faster because it satisfies the boundary conditions uniformly (Li, 1996) (reading note). A symmetric transformation by B⁻¹ᐟ² lets both problems use a real symmetric eigensolver. Tangential fields are matched at both interfaces in a complex linear system.

Forward and backward modal amplitudes are referenced to opposite layer faces, so evanescent propagation factors decay. Stable layered-grating algorithms use the same normalization: exponentially growing factors never appear (Li, 1996)(Li, 2014, §13.5.2). Modes at zero longitudinal wavevector use a constant/linear basis instead of shifting the physical wavelength. The convention is exp(−iωt), with z increasing into the substrate.

Y is q for TE and q/n² for the TM magnetic-field amplitudes. Evanescent orders participate in interface matching, even though their far-field normal power is zero. For a broader account of fields, modes, and scattering matrices, see the RCWA mathematics documentation.

This browser solver adapts the centered Fourier coefficients and mode-matching formulation in SciApp’s Python solver. Development also used the NanoEM implementation in the local calculationsApp repository as a numerical reference.

The browser implementation reduces the original vector formulation to scalar TE/TM problems for a single isotropic layer. The isotropic mode basis is discussed in thesis §3.3 (p. 38). Its automated checks compare each diffraction order against eight fixtures generated by running the original NanoEM solver, and check Fresnel interfaces, analytic thin films, power conservation, normal-incidence symmetry, zero-height layers, and grazing modes.

The normalization of order-resolved powers follows the thesis’s diffraction-efficiency formulation (p. 43). The displayed balance error is |1 − R − T|. A small value checks conservation within the truncated calculation; it does not prove that enough Fourier modes or wavelength samples were used. Increase the harmonic count and sample density to check whether the quantities you care about have stabilized. This playground has not been validated across every possible setting.

The beam is a monochromatic sheet of light, localized along x and uniform along the grooves. Its incident amplitude at z = 0 is a Gaussian centered on x = 0, with a phase tilt set by the incidence angle. The radius w is measured along the interface: the field amplitude falls to 1/e at x = ±w (squared magnitude to 1/e²).

Each sampled angle gets its own RCWA solution. Their complex fields are added coherently, then the real part or squared magnitude is displayed. This gives reflection, refraction, spreading, and grating diffraction for the finite beam. The incident peak at the interface is normalized to one. Animation advances optical phase in a continuous beam; it is not a time-domain pulse simulation.

The angular spectrum is sampled uniformly in lateral wavevector using 41, 81 (default), or 121 midpoint samples over ±6/(k₀w), restricted to propagating angles below 88°. A warning appears when this removes appreciable Gaussian amplitude. Increase the angular sample count to check convergence. The calculation uses an infinite periodic structure; the drawing shows a finite observation window.

The power table and spectrum remain diagnostics for the central-angle plane wave, not integrated powers of the finite beam. See TU Delft’s angular-spectrum explanation for the plane-wave decomposition (Kalkman, n.d., §3.6).

import { sampleField, solveRCWA, type RCWASettings, type FieldBounds, type FieldGrid } from './rcwa';

export interface Illumination { kind: 'plane' | 'beam'; radius: number; samples: number }
export const DEFAULT_ILLUMINATION: Illumination = { kind: 'plane', radius: 2, samples: 81 };

/**
 * Angular spectrum of a 2D monochromatic beam, uniform along the grooves.
 * At z=0: F_inc(x) = exp(-(x/w)^2) exp(i k0 Kc x), before spectral truncation.
 * radius is the 1/e amplitude radius along the interface, expressed in vacuum λ.
 * Uniform midpoint quadrature in K=kx/k0, restricted to forward propagating waves.
 */
export function beamSpectrum(s: RCWASettings, radius: number, samples: number) {
  if (!Number.isFinite(radius) || radius < 1 || radius > 4) throw new Error('Beam radius must be between 1 and 4 wavelengths.');
  if (!Number.isInteger(samples) || samples < 21 || samples > 121) throw new Error('Use 21–121 angular samples.');
  const k0w = 2 * Math.PI * radius, center = s.incidentIndex * Math.sin(s.angle * Math.PI / 180);
  const span = 6 / k0w, limit = s.incidentIndex * Math.sin(88 * Math.PI / 180);
  const low = Math.max(-limit, center - span), high = Math.min(limit, center + span);
  if (high <= low) throw new Error('The beam has no incident propagating angular components.');
  const amplitude = (k: number) => Math.exp(-(((k - center) * k0w / 2) ** 2));
  const components = Array.from({ length: samples }, (_, i) => {
    const k = low + (i + .5) * (high - low) / samples;
    return { angle: Math.asin(k / s.incidentIndex) * 180 / Math.PI, weight: amplitude(k) };
  });
  const sum = components.reduce((total, component) => total + component.weight, 0);
  components.forEach(component => { component.weight /= sum; });
  const truncated = (low > center - span && amplitude(low) > .001) || (high < center + span && amplitude(high) > .001);
  return { components, truncated };
}

/** Coherent sum of complex RCWA fields, never a sum of intensities or an envelope mask. */
export function sampleBeam(s: RCWASettings, illumination: Illumination, width = 320, height = 200, bounds?: FieldBounds): FieldGrid {
  const { components, truncated } = beamSpectrum(s, illumination.radius, illumination.samples);
  const radius = illumination.radius * s.wavelength, depth = 3 * s.wavelength;
  const halfWidth = 3 * radius + depth * Math.abs(Math.tan(s.angle * Math.PI / 180));
  const viewport = bounds ?? { xMin: -halfWidth, xMax: halfWidth, zMin: -depth, zMax: s.thickness + depth };
  const real = new Float64Array(width * height), imaginary = new Float64Array(width * height);
  const warnings = new Set<string>();
  if (truncated) warnings.add('The angular spectrum is clipped to propagating waves below 88°. At this angle and width the incident profile departs from an ideal Gaussian.');
  for (const component of components) {
    const result = solveRCWA({ ...s, angle: component.angle });
    if (result.balanceError > 1e-6) throw new Error('An angular component failed power balance. Try a nearby angle or width.');
    result.warnings.forEach(warning => warnings.add(warning));
    const field = sampleField(result, width, height, viewport);
    for (let i = 0; i < real.length; i++) {
      real[i] += component.weight * field.real[i];
      imaginary[i] += component.weight * field.imaginary[i];
    }
  }
  const shortest = s.wavelength / Math.max(s.incidentIndex, s.substrateIndex, s.ridgeIndex, s.grooveIndex);
  if (Math.max((viewport.xMax - viewport.xMin) / (width - 1), (viewport.zMax - viewport.zMin) / (height - 1)) > shortest / 4) {
    warnings.add('The field image has fewer than four pixels per shortest material wavelength. Fine wave fringes may be undersampled; reduce the width or angle.');
  }
  let peak = 0;
  for (let i = 0; i < real.length; i++) peak = Math.max(peak, Math.hypot(real[i], imaginary[i]));
  return { ...viewport, width, height, real, imaginary, peak, beam: {
    samples: components.length, radius, minAngle: components[0].angle,
    maxAngle: components.at(-1)!.angle, warnings: [...warnings],
  } };
}

import { add, sub, mul, scale, norm, complex, root, propagate, matrix, multiply, symmetricEigen, symmetricPower, solveComplex, type Complex } from './rcwa-math';

/**
 * Lossless, isotropic, single binary layer; incidence perpendicular to the grooves.
 * Reference: calculationsApp/python/lib/functions/{rcwa,calc_coherent_point,s_matrix}.py
 * and SciApp/python/funkce.py. Uses their centered rectangular Fourier coefficients
 * and inverse factorization for TM, reduced to scalar TE/TM eigenproblems.
 * Coordinates here: periodic x, grooves y, depth +z; exp(-i omega t).
 */
export interface RCWASettings {
  wavelength: number; period: number; thickness: number; fill: number; angle: number;
  ridgeIndex: number; grooveIndex: number; incidentIndex: number; substrateIndex: number;
  harmonics: number; polarization: 'TE' | 'TM';
}
export const DEFAULT_RCWA: RCWASettings = {
  wavelength: 633, period: 700, thickness: 250, fill: .5, angle: 0,
  ridgeIndex: 2, grooveIndex: 1, incidentIndex: 1, substrateIndex: 1.5,
  harmonics: 7, polarization: 'TE',
};
export interface DiffractionOrder {
  order: number; reflection: number; transmission: number;
  reflectionAngle: number | null; transmissionAngle: number | null;
}
export interface RCWAResult {
  settings: RCWASettings; reflection: number; transmission: number; balanceError: number;
  orders: DiffractionOrder[]; warnings: string[];
  modes: { w: number[][]; q: Complex[]; a: Complex[]; b: Complex[]; linear: boolean[] };
  reflected: Complex[]; transmitted: Complex[]; kx: number[];
}
export interface FieldBounds { xMin: number; xMax: number; zMin: number; zMax: number }
export interface FieldGrid extends FieldBounds {
  width: number; height: number; real: Float64Array; imaginary: Float64Array; peak: number;
  beam?: { samples: number; radius: number; minAngle: number; maxAngle: number; warnings: string[] };
}
export interface SpectrumPoint { wavelength: number; reflection: number | null; transmission: number | null; error?: string }

function validate(s: RCWASettings) {
  for (const key of ['wavelength', 'period', 'ridgeIndex', 'grooveIndex', 'incidentIndex', 'substrateIndex'] as const) {
    if (!Number.isFinite(s[key]) || s[key] <= 0) throw new Error(`${key} must be positive and finite.`);
  }
  if (!Number.isFinite(s.thickness) || s.thickness < 0 || !Number.isFinite(s.fill) || s.fill < 0 || s.fill > 1 || !Number.isFinite(s.angle) || Math.abs(s.angle) >= 89) throw new Error('Invalid thickness, fill fraction, or incidence angle.');
  if (!Number.isInteger(s.harmonics) || s.harmonics < 0 || s.harmonics > 20) throw new Error('Use 0–20 harmonics on each side of the zero order.');
  if (s.polarization !== 'TE' && s.polarization !== 'TM') throw new Error('Choose TE or TM polarization.');
}

/** Analytic Fourier convolution for a ridge centered at x=0 in each period. */
export function convolution(n: number, fill: number, ridge: number, groove: number) {
  if (fill === 0 || fill === 1 || ridge === groove) return matrix(n, (i, j) => i === j ? (fill === 0 ? groove : fill === 1 ? ridge : groove) : 0);
  return matrix(n, (i, j) => i === j ? fill * ridge + (1 - fill) * groove :
    (ridge - groove) * Math.sin(Math.PI * fill * (i - j)) / (Math.PI * (i - j)));
}

export function solveRCWA(settings: RCWASettings): RCWAResult {
  validate(settings);
  const s = { ...settings }, n = 2 * s.harmonics + 1, tm = s.polarization === 'TM';
  const k0 = 2 * Math.PI / s.wavelength, length = k0 * s.thickness;
  const kx = Array.from({ length: n }, (_, i) => s.incidentIndex * Math.sin(s.angle * Math.PI / 180) + (i - s.harmonics) * s.wavelength / s.period);
  const eps = convolution(n, s.fill, s.ridgeIndex ** 2, s.grooveIndex ** 2);
  const inverseEps = convolution(n, s.fill, 1 / s.ridgeIndex ** 2, 1 / s.grooveIndex ** 2);
  let w: number[][], d: number[][], eigenvalues: number[];
  if (tm) {
    // (I - K E^-1 K) H = q² B H, B = convolution(1/epsilon).
    // B^-1/2 makes the eigenproblem real symmetric; D = B W supplies admittance.
    const epsInverse = symmetricPower(eps, -1), transform = symmetricPower(inverseEps, -.5);
    const c = matrix(n, (i, j) => +(i === j) - kx[i] * epsInverse[i][j] * kx[j]);
    const eig = symmetricEigen(multiply(multiply(transform, c), transform));
    eigenvalues = eig.values; w = multiply(transform, eig.vectors); d = multiply(inverseEps, w);
  } else {
    const eig = symmetricEigen(matrix(n, (i, j) => eps[i][j] - (i === j ? kx[i] ** 2 : 0)));
    eigenvalues = eig.values; w = eig.vectors; d = w;
  }
  const q = eigenvalues.map(root), linear = eigenvalues.map(value => Math.abs(value) < 1e-12);
  const phase = q.map(value => propagate(value, length));
  const outerQ = (index: number) => kx.map(value => root(index ** 2 - value ** 2));
  const qr = outerQ(s.incidentIndex), qt = outerQ(s.substrateIndex);
  const yr = qr.map(value => scale(value, tm ? 1 / s.incidentIndex ** 2 : 1));
  const yt = qt.map(value => scale(value, tm ? 1 / s.substrateIndex ** 2 : 1));
  const incidentFlux = yr[s.harmonics].re;
  const system: Complex[][] = Array.from({ length: 2 * n }, () => Array.from({ length: 2 * n }, () => complex()));
  const rhs = Array.from({ length: 2 * n }, (_, i) => complex(i === s.harmonics ? 2 * incidentFlux : 0));
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) {
    const wi = complex(w[i][j]), v = scale(q[j], d[i][j]);
    if (linear[j]) {
      // At q=0 the two exponentials coincide. Use exact constant/linear solutions
      // F=W(a+k0*z*b), G=-i D b instead of perturbing the wavelength or angle.
      system[i][j] = mul(yr[i], wi);
      system[i][n + j] = complex(0, -d[i][j]);
      system[n + i][j] = scale(mul(yt[i], wi), -1);
      system[n + i][n + j] = sub(complex(0, -d[i][j]), scale(mul(yt[i], wi), length));
    } else {
      // Modes are referenced to opposite faces, keeping all evanescent factors ≤1.
      system[i][j] = add(v, mul(yr[i], wi));
      system[i][n + j] = mul(sub(mul(yr[i], wi), v), phase[j]);
      system[n + i][j] = mul(sub(v, mul(yt[i], wi)), phase[j]);
      system[n + i][n + j] = scale(add(v, mul(yt[i], wi)), -1);
    }
  }
  // An unexcited, uncoupled harmonic grazing in identical media has an arbitrary
  // constant field: both interface equations only constrain its derivative.
  // Select zero amplitude for that free solution, leaving the incident mode alone.
  if (s.fill === 0 || s.fill === 1 || s.ridgeIndex === s.grooveIndex) {
    for (let i = 0; i < n; i++) {
      if (i !== s.harmonics && linear[i] && norm(qr[i]) < 1e-14 && norm(qt[i]) < 1e-14) {
        system[i] = Array.from({ length: 2 * n }, (_, j) => complex(j === i ? 1 : 0));
      }
    }
  }
  const coefficients = solveComplex(system, rhs), a = coefficients.slice(0, n), b = coefficients.slice(n);
  const reflected = Array.from({ length: n }, (_, i) => {
    let value = complex(i === s.harmonics ? -1 : 0);
    for (let j = 0; j < n; j++) value = add(value, scale(linear[j] ? a[j] : add(a[j], mul(phase[j], b[j])), w[i][j]));
    return value;
  });
  const transmitted = Array.from({ length: n }, (_, i) => {
    let value = complex();
    for (let j = 0; j < n; j++) value = add(value, scale(linear[j] ? add(a[j], scale(b[j], length)) : add(mul(phase[j], a[j]), b[j]), w[i][j]));
    return value;
  });
  const orders = kx.map((value, i) => ({
    order: i - s.harmonics,
    reflection: norm(reflected[i]) * yr[i].re / incidentFlux,
    transmission: norm(transmitted[i]) * yt[i].re / incidentFlux,
    reflectionAngle: qr[i].re > 0 ? Math.asin(value / s.incidentIndex) * 180 / Math.PI : null,
    transmissionAngle: qt[i].re > 0 ? Math.asin(value / s.substrateIndex) * 180 / Math.PI : null,
  }));
  const reflection = orders.reduce((sum, order) => sum + order.reflection, 0);
  const transmission = orders.reduce((sum, order) => sum + order.transmission, 0);
  if (![reflection, transmission].every(Number.isFinite)) throw new Error('The calculation produced non-finite powers. Try a nearby setting.');
  const balanceError = Math.abs(1 - reflection - transmission), warnings: string[] = [];
  if (balanceError > 1e-6) warnings.push('Power balance is inaccurate at this setting. Treat the result as unreliable.');
  const omitted = [-s.harmonics - 1, s.harmonics + 1].some(m => Math.abs(s.incidentIndex * Math.sin(s.angle * Math.PI / 180) + m * s.wavelength / s.period) < Math.max(s.incidentIndex, s.substrateIndex));
  if (omitted) warnings.push('Some propagating orders are outside the retained harmonics. Increase the harmonic count.');
  if ([...qr, ...qt].some(value => norm(value) < 1e-14)) warnings.push('An exterior order is at grazing cutoff. It carries zero normal power at the threshold.');
  return { settings: s, reflection, transmission, balanceError, orders, warnings, modes: { w, q, a, b, linear }, reflected, transmitted, kx };
}

/** Normalized transverse scalar field: Ey for TE, Hy for TM (not total |E|²). */
export function sampleField(result: RCWAResult, width = 192, height = 128, bounds?: FieldBounds): FieldGrid {
  const s = result.settings, { w, q, a, b, linear } = result.modes, n = result.kx.length;
  const k0 = 2 * Math.PI / s.wavelength;
  const { xMin, xMax, zMin, zMax } = bounds ?? { xMin: -1.5 * s.period, xMax: 1.5 * s.period, zMin: -.8 * s.wavelength, zMax: s.thickness + .8 * s.wavelength };
  const real = new Float64Array(width * height), imaginary = new Float64Array(width * height);
  const lateral = Array.from({ length: width }, (_, x) => result.kx.map(k => propagate(complex(k), k0 * (xMin + (xMax - xMin) * x / (width - 1)))));
  let peak = 0;
  for (let row = 0; row < height; row++) {
    const z = zMin + (zMax - zMin) * row / (height - 1);
    const field = Array.from({ length: n }, () => complex());
    if (z < 0 || z > s.thickness) {
      for (let i = 0; i < n; i++) {
        const outer = root((z < 0 ? s.incidentIndex : s.substrateIndex) ** 2 - result.kx[i] ** 2);
        field[i] = mul(z < 0 ? result.reflected[i] : result.transmitted[i], propagate(outer, k0 * (z < 0 ? -z : z - s.thickness)));
        if (z < 0 && i === s.harmonics) field[i] = add(field[i], propagate(outer, k0 * z));
      }
    } else {
      for (let j = 0; j < n; j++) {
        const modal = linear[j] ? add(a[j], scale(b[j], k0 * z)) : add(mul(a[j], propagate(q[j], k0 * z)), mul(b[j], propagate(q[j], k0 * (s.thickness - z))));
        for (let i = 0; i < n; i++) field[i] = add(field[i], scale(modal, w[i][j]));
      }
    }
    for (let x = 0; x < width; x++) {
      let value = complex();
      for (let i = 0; i < n; i++) value = add(value, mul(field[i], lateral[x][i]));
      real[row * width + x] = value.re; imaginary[row * width + x] = value.im;
      peak = Math.max(peak, Math.sqrt(norm(value)));
    }
  }
  return { width, height, real, imaginary, peak, xMin, xMax, zMin, zMax };
}

/** Small dense linear algebra for the lossless RCWA solver. No runtime dependencies. */
export type Matrix = number[][];
export const matrix = (n: number, f: (i: number, j: number) => number): Matrix =>
  Array.from({ length: n }, (_, i) => Array.from({ length: n }, (_, j) => f(i, j)));
export const identity = (n: number) => matrix(n, (i, j) => +(i === j));
export const multiply = (a: Matrix, b: Matrix) => matrix(a.length, (i, j) =>
  a[i].reduce((sum, value, k) => sum + value * b[k][j], 0));

/** Cyclic Jacobi rotations: orthonormal eigenvectors of a real symmetric matrix. */
export function symmetricEigen(input: Matrix) {
  const n = input.length, a = input.map(row => [...row]), vectors = identity(n);
  const scale = Math.max(1, ...a.flat().map(Math.abs));
  for (let sweep = 0; sweep < 80; sweep++) {
    let largest = 0;
    for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) {
      const off = a[p][q];
      largest = Math.max(largest, Math.abs(off));
      if (Math.abs(off) < 1e-14 * scale) continue;
      const tau = (a[q][q] - a[p][p]) / (2 * off);
      const t = (tau >= 0 ? 1 : -1) / (Math.abs(tau) + Math.hypot(1, tau));
      const c = 1 / Math.hypot(1, t), s = t * c;
      a[p][p] -= t * off; a[q][q] += t * off;
      a[p][q] = a[q][p] = 0;
      for (let k = 0; k < n; k++) {
        if (k !== p && k !== q) {
          const ap = a[k][p], aq = a[k][q];
          a[k][p] = a[p][k] = c * ap - s * aq;
          a[k][q] = a[q][k] = s * ap + c * aq;
        }
        const vp = vectors[k][p], vq = vectors[k][q];
        vectors[k][p] = c * vp - s * vq;
        vectors[k][q] = s * vp + c * vq;
      }
    }
    if (largest < 1e-12 * scale) return { values: a.map((row, i) => row[i]), vectors };
  }
  throw new Error('The mode solver did not converge. Try fewer harmonics or a lower index contrast.');
}

export function symmetricPower(a: Matrix, power: number) {
  const { values, vectors: u } = symmetricEigen(a);
  if (values.some(value => value <= 0)) throw new Error('Expected positive dielectric permittivity.');
  return matrix(a.length, (i, j) => values.reduce((sum, value, k) => sum + u[i][k] * value ** power * u[j][k], 0));
}

export type Complex = { re: number; im: number };
export const complex = (re = 0, im = 0): Complex => ({ re, im });
export const add = (a: Complex, b: Complex) => complex(a.re + b.re, a.im + b.im);
export const sub = (a: Complex, b: Complex) => complex(a.re - b.re, a.im - b.im);
export const mul = (a: Complex, b: Complex) => complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re);
export const scale = (a: Complex, s: number) => complex(a.re * s, a.im * s);
export const norm = (a: Complex) => a.re ** 2 + a.im ** 2;
export const div = (a: Complex, b: Complex) => scale(mul(a, complex(b.re, -b.im)), 1 / norm(b));
export const root = (value: number) => value >= 0 ? complex(Math.sqrt(value)) : complex(0, Math.sqrt(-value));
/** exp(i q distance), with positive imaginary q giving decay, never growth. */
export const propagate = (q: Complex, distance: number) => scale(complex(Math.cos(q.re * distance), Math.sin(q.re * distance)), Math.exp(-q.im * distance));

/** Complex Gaussian elimination with scaled partial pivoting. */
export function solveComplex(input: Complex[][], rhs: Complex[]) {
  const n = rhs.length, a = input.map((row, i) => [...row, rhs[i]]);
  const rowScale = a.map(row => Math.max(...row.slice(0, n).map(value => Math.sqrt(norm(value)))));
  for (let k = 0; k < n; k++) {
    let pivot = k;
    for (let i = k + 1; i < n; i++) {
      if (Math.sqrt(norm(a[i][k])) / rowScale[i] > Math.sqrt(norm(a[pivot][k])) / rowScale[pivot]) pivot = i;
    }
    if (!Number.isFinite(norm(a[pivot][k])) || Math.sqrt(norm(a[pivot][k])) < 1e-13 * rowScale[pivot] || !rowScale[pivot]) {
      throw new Error('Interface matching is singular at this setting. Try a nearby wavelength or angle.');
    }
    [a[k], a[pivot]] = [a[pivot], a[k]];
    [rowScale[k], rowScale[pivot]] = [rowScale[pivot], rowScale[k]];
    for (let i = k + 1; i < n; i++) {
      const factor = div(a[i][k], a[k][k]);
      for (let j = k + 1; j <= n; j++) a[i][j] = sub(a[i][j], mul(factor, a[k][j]));
      a[i][k] = complex();
    }
  }
  const x = Array.from({ length: n }, () => complex());
  for (let i = n - 1; i >= 0; i--) {
    let value = a[i][n];
    for (let j = i + 1; j < n; j++) value = sub(value, mul(a[i][j], x[j]));
    x[i] = div(value, a[i][i]);
  }
  return x;
}