Add a scattering-matrix backend for numerical stability in thick/absorbing stacks

[garrekstemo]

Summary

Add a scattering-matrix (S-matrix) backend as a numerically stable alternative to the transfer-matrix product for thick, strongly absorbing, or high-birefringence stacks, where the current formulation can overflow/underflow.

Context

The stack is assembled by multiplying interface matrices L = D_{i-1} \ D_i and diagonal propagation matrices P containing exp(-iωqz/c) (src/general_TMM.jl, propagation_matrix in src/matrix_constructors.jl). For evanescent modes or thick absorbing layers, growing and decaying exponentials appear in the same matrix; their product across many layers becomes ill-conditioned, losing precision or producing Inf/NaN. This is the classic transfer-matrix instability, and it will be exacerbated by the planned graded/chiral feature (many thin slices → long products).

Why

The transfer-matrix method is known to become numerically unstable for thick systems and high birefringence; the standard, well-established remedy is to reformulate propagation in terms of scattering matrices, which keep only bounded (decaying) exponentials and remain stable for arbitrarily thick stacks. Mature anisotropic-multilayer toolkits (e.g. PyLlama, arXiv:2012.05945) adopt the S-matrix approach specifically for this reason.

Proposed approach

1. Implement S-matrix combination (the Redheboum-style "star product") layer-by-layer as an alternative to forming the global Γ, reusing the existing per-layer eigenmodes/D/q from layer_matrices.
2. Extract R/T from the assembled S-matrix; verify it reproduces the transfer-matrix results on the existing test suite within tolerance.
3. Add a stress test that demonstrates the transfer-matrix path losing accuracy (thick absorbing / high-contrast stack) while the S-matrix path stays well-conditioned and energy-conserving.
4. Decide policy: keep transfer-matrix as default and offer S-matrix opt-in, or auto-select based on conditioning / total optical thickness.

Acceptance criteria

• S-matrix backend reproduces the transfer-matrix R/T on existing tests.
• A documented thick/absorbing case where the S-matrix path is stable and the transfer path is not.

Caveats / risks

• Independent of, but synergistic with, the matrix-exponential and graded/chiral issues — all three concern robustness/generality of layer propagation. Sequencing: this one is most valuable once graded/chiral introduces long layer stacks.
• Get the star-product convention consistent with this package's field ordering and exp(-iωt) time convention.

References

• PyLlama: a stable and versatile Python toolkit for EM modeling of multilayered anisotropic media — https://arxiv.org/abs/2012.05945
• Standard references on S-matrix / Redheffer star-product stabilization of layered-media propagation.