File:Grating Diffraction vs Period.webm

\[Lambda]0 = 1.; k0 = 
 N[(2 \[Pi])/\[Lambda]0]; (*The wavelength in vacuum is set to 1, so \
all lengths are now in units of wavelengths*)
\[Delta] = \[Lambda]0/50; \[CapitalDelta] =  20*\[Lambda]0; (*Parameters for the grid*)
\[Sigma] = 20 \[Lambda]0; (*width of the gaussian beam*)

sourcef[x_, y_] := E^(-(x^2/(2 \[Sigma]^2))) E^(-((y + \[CapitalDelta]/            2)^2/(2 (\[Lambda]0/2)^2))) E^(I k0 y);
\[Phi]in = Table[Chop[sourcef[x, y]], {x, -\[CapitalDelta]/2, \[CapitalDelta]/    2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}]; (*Discretized source*)
d = \[Lambda]0/2; (*typical scale of the absorbing layer*)

imn = Table[ Chop[5 (E^-((x + \[CapitalDelta]/2)/d) + E^((x - \[CapitalDelta]/2)/d) + E^-((y + \[CapitalDelta]/2)/d) + E^((y - \[CapitalDelta]/2)/d))], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}]; (*Imaginary part of the refractive index (used to emulate absorbing boundaries)*)
dim = Dimensions[\[Phi]in][[1]];
L = -1/\[Delta]^2*KirchhoffMatrix[GridGraph[{dim, dim}]]; (*Discretized Laplacian*)

stopstep[t_] := 10 (t - 1)^4 + 0.5;
frames = Table[
   period = stopstep[t];
   ren = Table[
     If[-\[CapitalDelta]/2 < y < Cos[(2 \[Pi])/period x] - \[CapitalDelta]/2 + 3, 2, 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
   n = ren + I imn;
   b = -(Flatten[n]^2 - 1) k0^2 Flatten[\[Phi]in]; (*Right-hand side of the equation we want to solve*)
M = L + DiagonalMatrix[SparseArray[Flatten[n]^2 k0^2]]; (*Operator on the left-hand side of the equation we want to solve*)
\[Phi]s = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
   
ImageAdd[
    ArrayPlot[Transpose[(Abs[(\[Phi]in + \[Phi]s)]^2/Max[(Abs[\[Phi]in + \[Phi]s]^2)[[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]]])][[(       4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]],      ColorFunction -> "AvocadoColors", DataReversed -> True, Frame -> False, PlotRange -> {0, 1}]
    ,
ArrayPlot[ Transpose[Re[(n - 1)/5]] [[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], DataReversed -> True ,      ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
    ](*Plot everything*)
   , {t, 0, 1, 1/100}];

ListAnimate[Join[Table[frames[[1]], {2}], frames, Reverse@frames]]