It would be nice to have an tutorial example showing that an arbitrary incident wave (e.g. a Gaussian source) can be scattered off of a periodic surface (say in 2d) by Brillouin-zone (BZ) integration of unit-cell calculations. (Examples of this can also be found in fmmax, and it's sometimes called the "phased-array method" or something like that).
Given an arbitrary current source $J(x,y)$ above surface with period $a$, we just do a bunch of Bloch-periodic simulations with a current source $\hat{J}(k,x,y) e^{ikx}$ for $k \in (-\pi/a, \pi/a]$ (the BZ), where $\hat{J}(k,x,y)$ is periodic in $k$, using the identity:
$$
J(x,y) = \frac{a}{2\pi} \int_{-\pi/a}^{+\pi/a} dk , \hat{J}(k,x,y) e^{ikx} \Longleftrightarrow \hat{J}(k,x,y) = \sum_{n=-\infty}^{\infty} J(x,y) e^{-kna}
$$
You can just compute this sum numerically for each source point $(x,y)$ if $J(x,y)$ is localized, since then the summand presumably decays rapidly.
Then, by linearity, you can just sum up the resulting fields $\hat{E}(k,x,y)$ from the unit-cell calculations to find the total field $E(x,y) = \frac{a}{2\pi} \int_{-\pi/a}^{+\pi/a} \hat{E}(k,x,y) , dk$ at any desired point. (By the usual orthogonality relations, you can also compute things like the transmitted or reflected powers by just integrating the powers from each unit-cell calculation.)
For example, we could have a Gaussian source $J(x,0)$ above a periodic grating, and compute the transmitted beam profile and the total transmitted power. (One could also use a Gaussian beam: in that case you'd have to go from our analytical Gaussian beam fields to the corresponding equivalent current $J$, not sure if we have a high-level API for this?)
