StableRCWA
Header-only C++17 RCWA for 1-D and 2-D periodic gratings that stays accurate at high truncation orders.
RCWA here means rigorous coupled-wave analysis, also called the Fourier modal method.
中文说明
Language note. Program output, error messages, the GUI and the source comments are currently in Chinese.
This README and the C++ API can be used without reading Chinese. The detailed JSON field reference and the deeper
documentation are in Chinese for now. On Windows, run chcp 65001 first so that the UTF-8 output displays correctly.
Why StableRCWA
-
Stable at high orders. Every layer eigenproblem is solved with a unitary Hessenberg-QR iteration, so accuracy does not
degrade as the number of orders grows. Fields propagate by S-matrix recursion. In 1-D TE/TM, modes with exactly zero
propagation constant get an exact treatment; in conical incidence, degenerate layers switch to a Schur basis.
The table below solves the same structure with an elimination-based (LR) eigen solver for comparison:
| Orders M (2M+1 harmonics) | QR (default): Littrow-order efficiency | LR: Littrow-order efficiency | LR: R + T | LR eig_info |
|---|
| 60 | 98.586939 % | 98.586939 % | 100.000000 % | 0 |
| 100 | 98.586941 % | 98.587496 % | 100.000506 % | 0 |
| 150 | 98.586942 % | 86.618778 % | 86.900177 % | 0 |
The LR solver reports success every time. The only warning sign is the eigen-pair residual, which rcwa_demo prints:
it is 1.4e-4, 1.4e-1 and 1.2e-1 for LR. For QR it stays at or below 8.5e-16, and R + T stays at 100.000000 %.
Single-period dielectric grating on a 20-layer mirror, lossless, TE, λ = 1.053 µm, θ = 66.36°.
Command: ./rcwa_demo --case single --M 150 --eig lr. Measured on Apple M1 Pro with Apple clang 17.
The LR numbers change with compiler and CPU; the QR numbers do not.
-
No external dependencies. Include the headers and call a function. The FFT library pocketfft is bundled;
there is nothing to download, build or link first.
-
Complete toolset:
- 1-D gratings, TE/TM;
- conical incidence, with a full 2×2 Jones matrix per order;
- 2-D crossed gratings;
- six-component near fields and phase maps;
- absorbing (complex-index) media.
-
Independently checked:
- against COMSOL finite-element models (2-D for conical incidence, 3-D for crossed gratings), comparing complex amplitudes and not only efficiencies;
- against closed-form solutions and exact physical identities (flux unitarity, reciprocity, mirror, translation and period-repetition equivalence).
-
Three ways to use it:
- C++ headers;
- a JSON-in/JSON-out command-line program, callable from Python, MATLAB or scripts;
- a local web GUI for modelling and viewing fields without writing code.
Quick start (one minute)
You need any C++17 compiler.
git clone https://github.com/aikkkkko/StableRCWA.git
cd StableRCWA
c++ -std=c++17 -O2 -I include -I examples examples/rcwa_demo.cpp -o rcwa_demo
./rcwa_demo --case single --M 20
The first lines of the output:
结构 single，TE，M = 20（N = 41），本征分解用 QR 迭代
第 1 级反射效率 98.586865 %
R = 99.794843 % T = 0.205157 % R + T = 100.000000 %
Two more lines follow: residuals and timing. The structure is lossless, so R + T must be 100 %.
- Linux: add
-pthread to the compile command. Linux has not been tested yet.
- Windows: build with
g++ ... -static, or with cl /std:c++17 /O2 /EHsc /utf-8 /I include /I examples examples\rcwa_demo.cpp.
Use it from C++
#include <cstdio>
#include "stablercwa/rcwa1d.hpp"
using namespace stablercwa;
int main() {
rcwa::Stack stack; // top to bottom: incidence medium, layers, substrate
stack.nI = 1.0; stack.nII = 1.46; stack.period = 0.574; // any length unit, used consistently
// grating layer, 0.38 thick: [0, 0.172] has index 1.46, [0.172, 0.574] is air
stack.layers.push_back(rcwa::Layer::grating(0.38, {0.0, 0.172, 0.574}, {1.46, 1.0}));
stack.layers.push_back(rcwa::Layer::uniform(0.115, 1.96)); // homogeneous film
const int M = 20; // keep orders -M ... M
rcwa::Result r = rcwa::solve(stack, 1.053, 66.36, M, rcwa::Polarization::TE); // wavelength, angle (deg)
for (int m = -1; m <= 1; ++m)
std::printf("m = %+d R = %.6f T = %.6f\n", m, r.DEr[m + M], r.DEt[m + M]);
std::printf("R + T = %.12f\n", r.reflectance() + r.transmittance());
}
c++ -std=c++17 -O2 -I include my_grating.cpp -o my_grating && ./my_grating
m = -1 R = 0.000000 T = 0.000000
m = +0 R = 0.026034 T = 0.497012
m = +1 R = 0.258135 T = 0.218818
R + T = 1.000000000000
To change the setup:
- Absorbing media: write the index as
cplx(n, k) with k > 0.
- Conical incidence:
rcwa::solve_conical in rcwa1d_conical.hpp.
- 2-D crossed gratings:
rcwa::solve2d in rcwa2d.hpp.
- Near fields: the
*_fields.hpp headers.
- Numerical failures (singular system, eigen solver not converged) do not throw. The efficiencies come back as NaN.
In 1-D,
eig_info holds the eigen solver's code; in 2-D, status and message say what went wrong.
- Invalid input (non-increasing edges, |θ| ≥ 90°, …) throws
std::invalid_argument.
Using it from a CMake project
set(STABLERCWA_BUILD_EXAMPLES OFF)
set(STABLERCWA_BUILD_TESTS OFF)
add_subdirectory(StableRCWA)
target_link_libraries(your_app PRIVATE StableRCWA::stablercwa)
JSON in, JSON out
c++ -std=c++17 -O2 -I include -I examples examples/rcwa_json.cpp -o rcwa_json
./rcwa_json request.json --out reply.json
{
"mode": "planar",
"stack": {"nI": 1.0, "nII": 1.46, "period": 0.574,
"layers": [{"d": 0.38, "edges": [0.0, 0.172, 0.574], "n_values": [1.46, 1.0]},
{"d": 0.115, "n": 1.96}]},
"M": 20, "polarization": "TE",
"points": [{"wavelength": 1.053, "theta_deg": 66.36}]
}
The reply holds per-order efficiencies; here R_total = 0.28416918… and T_total = 0.71583081….
- The modes are
planar, conical and periodic2d.
- Listing several
points gives a wavelength/angle sweep. The points run in parallel when options.threads is greater than 1.
- A malformed request returns
{"ok": false, "error": …, "where": …}. A misspelled key also counts as an error.
Graphical interface
The GUI runs as a local web app on 127.0.0.1 and needs no network connection. It requires Python 3 with NumPy, SciPy and Matplotlib.
python3 -m pip install -r gui/requirements.txt
python3 gui/gui.py
- macOS: double-click
gui/启动GUI.command.
- Windows: run
gui\install_windows.bat once (Python 3.11 or 3.12, x64), then gui\start_gui.bat.
At start-up the GUI compiles the C++ solver if a compiler is available. Otherwise it uses the bundled Python solver
and says so on the results page.
Model page
Results page
2-D crossed grating: top view and linked cross-section
What is supported
| C++ headers | Python solver |
|---|
| 1-D gratings, TE/TM | ✓ | ✓ |
| 1-D conical incidence (Jones matrices) | ✓ | ✓ |
| 2-D crossed gratings (rectangular cells, Li factorization) | ✓ | ✓ |
| Absorbing media (complex index) | ✓ | ✓ |
| Six-component near fields and phase | ✓ | ✓ |
| Finite (non-periodic) structures with lateral PML ¹ | – | 1-D, in-plane incidence |
| Uniaxial layers (optic axis along x, y or z) ¹ | – | 1-D, TE/TM |
| FFT-sampled Fourier coefficients, as an alternative to the analytic ones | 1-D, TE/TM and conical | 1-D |
| Magnetic media (μ ≠ 1) | – | experimental, not in the GUI |
| General 6×6 constitutive tensors | – | 1-D, experimental, not in the GUI |
| GPU eigen solver (CUDA cuSOLVER) | optional, experimental | – |
| oneMKL eigen solver as a single-file drop-in for DLL projects (Windows x64) | optional, standalone/ | – |
oneMKL for the RCWA solvers' eigensolves, products and linear solves (mkl::use_mkl, JSON "backend": "mkl") | optional | – |
¹ Inherited from the author's earlier Python package; their validation reports will be added to docs/.
For a sampled FFT layer, use the bundled BSD-3-Clause pocketfft:
auto layer = rcwa::Layer::grating(0.38, {0.0, 0.1722, 0.574}, {1.46, 1.0},
rcwa::FourierMethod::FFT, 65536);
stack.layers.push_back(layer);
The final sample count defaults to 0 (automatic); explicit counts must satisfy N >= 4*M+1, up to 4194304,
and need not be powers of two. Analytic integration remains the default. JSON layers accept
"fourier_method": "fft", "fft_samples": 65536; the GUI also routes eligible 1-D FFT models to C++.
FFT samples ε and 1/ε independently at cell midpoints. Check convergence in both M and N.
See the Fourier API and sampling convention.
Limitations:
- You must check convergence in
M yourself. 2-D converges more slowly than 1-D.
- The incidence medium must be lossless. In conical incidence the exit medium must be lossless too.
- Multi-threaded eigen solves are verified bitwise identical to single-threaded ones only with FMA contraction off (
-ffp-contract=off).
Optional: oneMKL eigensolver drop-in for DLL projects
standalone/ holds MatrixEigenValueAndEigenVector_MKL.cpp, a single source file for C++ DLL
projects on Windows x64. It exports MatrixEigenValueAndEigenVector for real and for complex matrices, returns all
eigenvalues and right eigenvectors sorted by real part, and computes them with Intel oneMKL's LAPACK (dgeev for real
matrices, zgeev for complex ones). It is separate from the header-only library, which still has no dependencies.
- Faster. On a 16-core AMD Ryzen 9 5950X it is 1.5 to 5.5 times faster than
eig_qr with 8 threads for complex
matrices and 3.1 to 12.4 times faster for real ones (n = 300 to 4000). At n = 4000 it takes 36 s for a complex matrix
and 16 s for a real one, against 199 s.
- Needs oneMKL 2025.0.1 or newer. On AMD processors under Windows,
dgeev in oneMKL 2024.2.2 and 2025.0.0 writes
past its workspace and corrupts memory. The file checks the MKL version at run time and does not call dgeev on those
builds.
- One-click runtime.
standalone/install_mkl.bat downloads oneMKL 2026.2.0 from PyPI, checks the SHA-256 digests
and puts the DLLs in standalone/mkl/bin. Nothing is installed system-wide.
The Visual Studio setup, the conventions, the full timing table and the test are in
standalone/README.md.
Optional: oneMKL backend for the RCWA solvers
include/stablercwa/mkl/mkl_backend.hpp hands the eigendecompositions, the dense matrix products and the dense linear
solves of the RCWA solvers to Intel oneMKL. It plugs into the custom_eig, custom_multiply and custom_solve hooks of
the 1-D, conical and 2-D solvers, and mkl::use_mkl sets all three. JSON requests select it with
"options": {"backend": "mkl"}.
#include "stablercwa/mkl/mkl_backend.hpp"
stablercwa::set_num_threads(8); // MKL uses the same number of threads
stablercwa::rcwa::Options opt; // rcwa::Options2D for 2-D
stablercwa::mkl::Settings s;
s.real_path = true; // only if no layer absorbs (all indices real)
stablercwa::mkl::use_mkl(opt, s); // with absorbing layers: mkl::use_mkl(opt)
- Faster. On a 16-core Ryzen 9 5950X with 8 threads, a whole solve is 1.6 to 3.2 times faster for structures
with complex layer matrices and 3.1 to 5.9 times faster for symmetric lossless ones (orders 338 to 3201).
- Real layer matrices can use
dgeev. For a lossless structure that is symmetric about the origin or the cell
centre, the layer matrices are real up to rounding. With real_path on, the backend drops the rounding-level
imaginary part and calls dgeev, which is 2.2 to 2.5 times as fast as zgeev. Other structures use zgeev.
real_path is off by default. Absorption weaker than about k = 1e-13 cannot be told apart from rounding in the matrix,
yet near a high-Q resonance it still matters: in a test case, treating k = 1e-13 as zero changed the efficiencies by
3.5e-6. So turn it on only when every layer index is real. JSON requests decide it from their refractive indices.
- Same results to rounding, not bitwise. Efficiencies and near fields agree with the built-in solvers to about 1e-12,
all 24 frozen 2-D Python references pass, and singular interface systems fail the same way as with the built-in LU. The largest difference, 2.4e-11 with
dgeev (4.6e-12 with zgeev), is in
an ill-conditioned high-contrast TM case. There, dropping only the rounding-level imaginary part already moves QR's
own result by 1.2e-11. MKL's last digits change with the thread count and the processor.
- Optional. The CMake option
STABLERCWA_WITH_MKL (off by default) provides the target stablercwa_mkl. It needs
oneMKL 2025.0.1 or newer, which standalone/install_mkl.bat downloads. The rest of the library still has no
dependencies.
Setup, thread rules and the timing table are in the user guide (in Chinese).
Validation
| Check | Result |
|---|
| 1-D conical, against COMSOL 2-D FEM (8 cases, from a bare interface to a grating with multilayer) | per-channel efficiency ≤ 3.3e-6, complex Jones matrix ≤ 1.1e-5, phase ≤ 0.013° |
| 2-D crossed gratings, against COMSOL 3-D periodic FEM (rectangular pillars) | complex amplitudes 1.2e-4 to 4.1e-4 (shrinking with mesh refinement), efficiency ≤ 2.1e-4 |
| Exact zero propagation constant (1-D), against COMSOL | amplitude 2.1e-9, efficiency 8.3e-11 |
| C++ 2-D solver, against the Python implementation it was ported from | far field ≤ 3.1e-13 (21 cases), near-field slices ≤ 1.1e-12 |
| Unit and regression tests | 21 C++ test programs plus 4 CTest checks of the JSON mode switches (25 in all), all passing on macOS arm64 (Apple clang 17); the 21 test programs also pass on Windows with MinGW g++ 13.1 and MSVC 19.34; Linux is covered by CI only |
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --config Release -j
cd build && ctest -C Release
Details: docs/validation.zh-CN.md.
Explainer pages
Two illustrated write-ups of the COMSOL comparisons, in Chinese, with derivations and figures embedded.
GitHub shows .html files as source code, so read them on the project site (GitHub Pages) or open them from a local clone.
Documentation
Citing
If you use StableRCWA in published work, please cite it as described in CITATION.cff.
License
MIT. The eigen solvers re-implement algorithms from EISPACK (see Acknowledgements).
The bundled pocketfft retains its BSD-3-Clause license;
include that notice when redistributing source or binaries.
Acknowledgements
- Eigen solvers: follow the EISPACK routines CORTH/COMQR2 (QR) and COMHES/COMLR2 (LR, kept for comparison only).
- 2-D Fourier factorization: follows L. Li's rules.
- Independent reference: COMSOL Multiphysics® FEM models.
- GUI example
meent_conical: uses a tutorial structure from the Meent project.
