\[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/15; \[CapitalDelta] = 30*\[Lambda]0; (*Parameters for the grid*)

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*)
\[Sigma] = 10 \[Lambda]0;
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]}];
randomstructure[n_] := Table[RotateRight[ DiskMatrix[2, dim], {RandomInteger[{-Round[dim/2], Round[dim/2] }], RandomInteger[{ -Round[dim/2], Round[dim/2] - 10}]}], {n}];
\[Phi] = Table[ren = Clip[Total[randomstructure[100]] , {0, 1}] + 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*)
   Partition[LinearSolve[M, b], dim], {600}];
\[Phi]meanempty = Mean[\[Phi]];

averages = 200; nscatterers = 100;
scatterers = Table[{RandomInteger[{-Round[dim/2], Round[dim/2] }], RandomInteger[{ -Round[dim/2], Round[dim/2] - 10}]}, {averages}, {nscatterers}];
\[Phi] = Table[ ren = Clip[ Total[Join[ RotateRight[DiskMatrix[2, dim], #] & /@ scatterers[[j]], {2*RotateRight[DiskMatrix[4, dim], {0, -80}]}]  ] , {0, 1}] + 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*)
   Partition[LinearSolve[M, b], dim], {j, 1, averages}];

frames = Table[\[Phi]tmp = Mean[\[Phi][[1 ;; j]]] - \[Phi]meanempty;
   Grid[{
     {Style[StringForm["\!\(\*SubscriptBox[\(E\), \(``\)]\)", NumberForm[j, {3, 0}]], Black, Bold, FontSize -> 16], Style["<\!\(\*SubscriptBox[\(E\), \(tot\)]\)>", Black, Bold, FontSize -> 16], Style["<\!\(\*SubscriptBox[\(E\), \(scatter\)]\)>", Black, Bold, FontSize -> 16]},
     {
      ImageAdd[
       ArrayPlot[
        Transpose[(Re@\[Phi][[j]]/
           Max[Abs@\[Phi][[j]][[(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, FontSize -> 14}, 
        ColorFunctionScaling -> True, ColorFunction -> ReMapC, 
        ImageSize -> 300]
       ,
       ArrayPlot[
        Transpose[Clip[Total[Join[RotateRight[DiskMatrix[2, dim], #] & /@ scatterers[[j]], {2*RotateRight[DiskMatrix[4, dim], {0, -80}]}]  ] , {0, 1}]/5] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False, ImageSize -> 300, ImagePadding -> 0.001]],
      ArrayPlot[ Transpose[(Re@Mean[\[Phi][[1 ;; j]]]/Max[Abs@Mean[\[Phi][[1 ;; j]]][[(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, FontSize -> 14}, ColorFunctionScaling -> True, ColorFunction -> ReMapC,        ImageSize -> 300],
      ArrayPlot[ Transpose[(Re@\[Phi]tmp/Max[Abs@\[Phi]tmp[[(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, FontSize -> 14}, ColorFunctionScaling -> True, ColorFunction -> ReMapC, ImageSize -> 300]
      }}]
   , {j, 1, 100, 1}];
ListAnimate[Join[
  Flatten@Table[Table[frames[[j]], {11 - j}], {j, 1, 10}]
  ,
  frames[[11 ;; 100]]
  ]
 ]