PyStratify

Light scattering and emission by multilayered films, cylinders and spheres, in Python. One modal-response core solves all three geometries under plane waves, electric and magnetic point dipoles and focused (Cassegrain) beams, with numerics rebuilt so that every quantity stays finite and accurate: thin metal shells, multipole orders in the thousands, emitters a fraction of a nanometre from a surface.

PyStratify is the engine of multilayer.app, a free browser calculator for the same three geometries, with a theory guide and every comparison with other codes plotted.

Why PyStratify

• Three geometries, one method. The transfer-matrix method for films (after Byrnes's tmm, whose function names it keeps), its cylindrical analogue for coated wires and fibres, and for spheres the recursive transfer-matrix method of Moroz (2005) as used by STRATIFY (Rasskazov, Carney & Moroz, OSA Continuum 3, 2290 (2020)).
• Sources beyond the plane wave. Electric and magnetic dipoles inside or outside the structure, with decay rates split into radiated, guided and absorbed power; focused beams as in FT-IR microspectroscopy.
• Spheres go further. Absorbing, magnetic, gain and chiral (Pasteur) shells; 2D materials and Feibelman d-parameters on the interfaces; electric, magnetic, chiral and quadrupole emitters; angle-resolved far fields, near fields, stored energy, Casimir–Polder potentials and two-point Green's dyadics.
• Accurate where others break down. Special functions as complex logarithms and a scaled two-sided recursion keep 1e-12 agreement with 60-digit references, energy balance to 1e-14 at 0.05 nm from a lossless interface, and absorbing spheres to size parameter 20 000 (docs/NUMERICS.md).
• Thin-film ellipsometry, end to end. Dispersion models (Cauchy, Sellmeier, Lorentz, Gaussian, Tauc-Lorentz, Cody-Lorentz, Drude, poles, Bruggeman and Maxwell Garnett mixtures), Ψ/Δ in the instruments' convention with incoherent backside reflection, Woollam and SENTECH files, and least-squares fits with MSE, 90 % limits and correlations: fits of real measurements agree with pyElli + lmfit to 4e-7.
• Checked in public. Every comparison with another code is a script in benchmarks/ that anyone can rerun (below).
• Light to install. Pure Python on NumPy and SciPy; whole spectra are vectorised.

Install

pip install pystratify

Python ≥ 3.10, NumPy ≥ 1.24, SciPy ≥ 1.10. The development version: pip install git+https://github.com/multilayerapp/PyStratify.git.

Quick start

Lengths share one unit of your choosing (here nanometres), angles are in radians, and indices are n + ik. No optical constants are shipped: take them from a database such as refractiveindex.info.

Films. The planar module keeps the tmm interface: here a quarter-wave Bragg mirror, 8 pairs of TiO₂ / SiO₂ on glass.

import numpy as np
import pystratify as ps
from pystratify import planar

n_list = [1.0] + [2.35, 1.46] * 8 + [1.52]                    # incident medium, layers, substrate
d_list = [np.inf] + [600 / 4 / 2.35, 600 / 4 / 1.46] * 8 + [np.inf]
R = [planar.coh_tmm("s", n_list, d_list, 0.0, wl)["R"] for wl in np.linspace(400, 900, 501)]

Cylinders. A silicon nanowire of radius 80 nm at normal incidence, both polarisations.

wire = ps.solve_cylinder([80.0], [3.6, 1.0], 775.0)            # radii, indices (host last), wavelength
ps.cross_widths(wire, "axial-electric")["efficiency_scattering"]
ps.cross_widths(wire, "axial-magnetic")["efficiency_scattering"]

Spheres. A silica–gold nanoshell in water over a whole spectrum, with its multipole decomposition.

wavelength = np.linspace(450, 1100, 651)
n_gold = 0.2 + 3.5j                                           # use tabulated n, k for real work
n = np.stack([np.full(651, 1.45), np.full(651, n_gold), np.full(651, 1.33)], axis=1)
sol = ps.solve([60.0, 75.0], n, wavelength)                   # radii of every interface, n per shell
cs = ps.cross_sections(sol)
cs.q_ext, cs.q_sca, cs.q_abs                                  # efficiencies over wavelength
cs.ext_by_order[:, ps.TM, 0] / cs.geometric                   # the electric dipole's share

Emitters. Decay rates of a dipole beside a sphere, and above a silver mirror through the unified solve_problem interface, which takes the same arguments for all three geometries.

rates = ps.decay_rates([40.0], [n_gold, 1.0], 600.0, r=40.0 + np.geomspace(0.2, 100, 60))
rates.total, rates.radiative, rates.nonradiative               # (positions, [radial, tangential])

mirror = ps.Problem("films", [np.inf, 200.0, np.inf], [1.0, 0.135 + 3.99j, 1.52], 633.0,
                    ps.PointDipole(-100.0, "electric", layer=0))  # 100 nm above 200 nm of silver on glass
result = ps.solve_problem(mirror)
result["total"], result["escape"], result["absorbed"]          # perpendicular, parallel, isotropic

Ellipsometry. Fit a measured Woollam file: the native oxide on a silicon wafer.

from pystratify import dispersion as dm, ellipsometry as el

data = el.read_woollam("wafer.dat").select(wavelength=(250, 1000))       # Psi, Delta, sigmas
silicon = dm.Table.from_refractiveindex("Si_Green-2008.yml")            # a refractiveindex.info page
oxide = dm.Sellmeier(B=(0.6961663, 0.4079426, 0.8974794), C=(0.0684043**2, 0.1162414**2, 9.896161**2))
params = el.Parameters().add("d", 5.0, min=0, max=50)                  # nm
result = el.fit(lambda p: el.Sample([el.Layer(oxide, p["d"])], substrate=silicon), params, data,
                grid={"d": np.arange(0, 20, 1.0)})
print(result.report())                                                  # d = 2.245 ± 0.006 nm (90 %), MSE

Focused beams. An FT-IR objective (NA 0.5, central obscuration 0.23) on a polymer fibre.

beam = ps.FocusedBeam(0.5, 0.23)                              # NA, obscuration NA
wavenumber = np.linspace(1000, 4000, 61)                      # cm⁻¹
fibre = ps.focused_cylinders([5.0], [1.49 + 0.01j, 1.0], 1e4 / wavenumber, beam)   # lengths in µm here
fibre.apparent_absorbance

Far more is in the user guide: conventions, every source and output, chiral media, 2D sheets and the full sphere toolbox. The runnable scripts in examples/ (multipole spectra, scattering patterns, ellipsometry fits of measured data, near-field maps, decay rates, circular dichroism, quadrupole emitters, graphene-coated spheres…) each save a figure, and crosscheck/ recomputes figures of classic papers from their stated parameters.

Verified against independent codes

Codes that solve the same equations exactly agree with PyStratify to between 1e-15 and 2e-6, the larger differences being the reference's own precision (float32 in tmmax, a contour integral in smuthi). Solvers that discretize space agree to their mesh or grid error, 0.7–4%; for Meep and MNPBEM the difference was shown to shrink as the grid or mesh is refined. Each point is a script in benchmarks/ that writes its own report, and python benchmarks/verify_all.py reruns the whole study in pinned environments (benchmarks/README.md). MATRIX.md maps every check by geometry, source and quantity, and the Verification page plots each one.

Speed

Typical timings on one core:

More, and the accuracy to expect in every regime: docs/NUMERICS.md.

Documentation

PyStratify is before 1.0: the names listed as frozen in the API contract keep their meaning from 1.0 on. Questions and bug reports are welcome in issues.

Tests

pip install -e ".[test]"
pytest

Citing

PyStratify implements the method of STRATIFY; please cite the paper,

I. L. Rasskazov, P. S. Carney and A. Moroz, "STRATIFY: a comprehensive and versatile MATLAB code for a multilayered sphere," OSA Continuum 3, 2290 (2020), doi:10.1364/OSAC.399979,

and the software itself: doi:10.5281/zenodo.23287561 names every version and resolves to the latest; each release also has its own DOI on that Zenodo page, which is the one to cite for the exact code you ran (CITATION.cff).

Licence

GPL-3.0-or-later, as STRATIFY (this is a derivative work). The planar routines derive from Byrnes's tmm and keep its MIT notice (src/pystratify/PLANAR-LICENSE.txt).