File:Photonic Band Gap vs Wavelength.webm

\[Lambda]0 = 0.5; 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/10; \[CapitalDelta] = 50*\[Lambda]0; (*Parameters for the grid*)
ReMapC[x_] := RGBColor[(2 x - 1) UnitStep[x - 0.5], 0, (1 - 2 x) UnitStep[0.5 - x]];
\[Sigma] = 7 \[Lambda]0;
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[imn][[1]];
L = -1/\[Delta]^2*KirchhoffMatrix[GridGraph[{dim, dim}]]; (*Discretized Laplacian*)
ren0 = 2;
ren = Total[ Flatten[Table[RotateRight[DiskMatrix[7, dim], {30 i, 30 j}] , {i, -8, 8}, {j, -3, 8}  ], 1] ] + 1;
n = ren + I imn;
\[Lambda] =.;
frames = Table[
    k = N[(2 \[Pi])/\[Lambda]];    
   sourcef[x_, y_] := E^(-(x^2/(2 \[Sigma]^2))) E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]/2)^2))) E^(I k y);
   \[Phi]in = Table[Chop[sourcef[x, y] ], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
   b = -(Flatten[n]^2 - 1) k^2 Flatten[\[Phi]in]; (*Right-hand side of the equation we want to solve*)
   M = L + DiagonalMatrix[SparseArray[Flatten[n]^2 k^2]]; (*Operator on the left-hand side of the equation we want to solve*)
   \[Phi] = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
   ImageAdd[
    ArrayPlot[Transpose[(Re@\[Phi]/Max[Abs@\[Phi][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]]])][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], DataReversed -> True, Frame -> False, PlotRange -> {-1, 1}, LabelStyle -> {Black, Bold}, ColorFunctionScaling -> True, ColorFunction -> ReMapC]
    ,
    ArrayPlot[Transpose[(ren - 1)/5] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
    ]
   , {\[Lambda], 1., 2., 0.005}];
ListAnimate[frames,5]