\[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] = 30*\[Lambda]0; (*Parameters for the grid*) \[Sigma]s = 4 \[Lambda]0;
ReMapC[x_] := RGBColor[(2 x - 1) UnitStep[x - 0.5], 0, (1 - 2 x) UnitStep[0.5 - x]];
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*)

sourcef[x_, y_] := E^(-(x^2/(2 \[Sigma]s^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]}];
ren0 = 2; nscatterers = 100;
scatterers = Table[{RandomInteger[Round /@ {-dim/2, dim/2}], RandomInteger[Round /@ {-dim/8, dim/8}]}, nscatterers];
basicren = Clip[Sum[RotateRight[DiskMatrix[2, dim], scatterers[[j]] ], {j, nscatterers}], {0, 1}] ;

frames = Table[
   ren = 0.5 GaussianFilter[basicren, \[Sigma]]/Max[GaussianFilter[basicren, \[Sigma]]] + 1;
   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] = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
   
   Grid[{{Style["Re(E)", White, Bold, Large], Style["|E\!\(\*SuperscriptBox[\(|\), \(2\)]\)", White, Bold, 
       Large]}, {
      ImageAdd[
       ArrayPlot[ Transpose[(Re[\[Phi]][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]])/0.075], DataReversed -> True, Frame -> False, PlotRange -> {-1, 1}, LabelStyle -> {Black, Bold}, ColorFunctionScaling -> True, ColorFunction -> ReMapC, ClippingStyle -> {Blue, Red}, ImageSize -> 300, Background -> Black]
       ,
       ArrayPlot[Transpose[(ren - 1)/1] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
       ]
      ,
      ImageAdd[
       ArrayPlot[ Transpose[((Abs[\[Phi]][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]])/0.075)^2], DataReversed -> True, Frame -> False, PlotRange -> {0, 1}, LabelStyle -> {Black, Bold}, ColorFunctionScaling -> True, ColorFunction -> "AvocadoColors", ClippingStyle -> White, Background -> Black, ImageSize -> 300]
       ,
       ArrayPlot[Transpose[(ren - 1)/1] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
       ]
      }}, Background -> Black]
   , {\[Sigma], 1, 50, 2}];

ListAnimate[ Join[Table[frames[[1]], {5}], frames, Table[frames[[-1]], {5}], Reverse[frames]] ]