File:Multiple Scattering02.gif

c = 1; (*speed of light*)
\[Omega]min = 1.; \[Omega]max = 5.;
\[Omega]0 = Mean[{\[Omega]min, \[Omega]max}]; \[Sigma]0 = (\[Omega]max - \[Omega]min)/10;
amplitude[w_] := E^(-(1/(2 \[Sigma]0^2)) ((w - \[Omega]0)^2) );

k0 = N[\[Omega]0/c];  \[Lambda]0 = N[(2 \[Pi])/k0]; d = \[Lambda]0/2; (*typical scale of the absorbing layer*)
\[Delta] = \[Lambda]0/10; \[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]];
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 = 1.5 - 1;
ren = ren0*Clip[Total[ Table[RotateRight[DiskMatrix[3, dim], {RandomInteger[{0, dim}], RandomInteger[{-50, Round[dim/2] - 15}]}], {200}] ], {0, 1}] + 1;
n = ren + I imn;
\[Sigma] = 2 \[Lambda]0;
sourcef[x_, y_, w_] := Sqrt[w/c] E^(-(x^2/(2 \[Sigma]^2))) E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]0/2)^2))) E^(I w/c y);
\[Delta]\[Omega] = (\[Omega]max - \[Omega]min)/200;
\[Phi] = Table[
   \[Phi]in = Table[amplitude[\[Omega]]*sourcef[x, y, \[Omega]] , {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
   b = -(Flatten[n]^2 - 1) (\[Omega]/c)^2 Flatten[\[Phi]in]; (*Right-hand side of the equation we want to solve*)
   M = L + DiagonalMatrix[ SparseArray[Flatten[n]^2 (\[Omega]/c)^2]]; (*Operator on the left-hand side of the equation we want to solve*)
   Partition[LinearSolve[M, b], dim], {\[Omega], \[Omega]min, \[Omega]max, 1*\[Delta]\[Omega]}];
\[Phi]dim = Dimensions[\[Phi]][[1]];

frames = Table[ ImageAdd[
ArrayPlot[ Transpose[((Abs@Total[Table[\[Phi][[j]] E^(I (\[Omega]min + \[Delta]\[Omega] (j - 1) ) t), {j, 1, \[Phi]dim}] ][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]])/3)^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)/5] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
    ]
   , {t, 32, -60, -1}];

ListAnimate[frames]