File:Light scattering - Averaging.gif

\[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/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] = 10 \[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*)
optimizedpoint = Round[0.7 dim];
sourcelist = Table[x0, {x0, -\[CapitalDelta]/2 + 2, \[CapitalDelta]/2 - 2, 1}];
sourcedim = Dimensions[sourcelist][[1]];
ren0 = 2;

\[Phi] = Table[
   ren = ren0*Clip[Total[Table[RotateRight[DiskMatrix[5, dim], {RandomInteger[{-Round[dim/2], Round[dim/2] }], RandomInteger[{ -Round[dim/8], Round[dim/2] - 10}]}], {200}]], {0, 1}] + 1;
   n = ren + I imn;
   phases = Table[
     \[Phi]in = Table[E^(-((x - sourcelist[[j]])^2/(2  (\[Lambda]0/2)^2))) E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]0/2)^2))), {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];(*Discretized source*)
     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*)
     Arg[\[Phi][[Round[dim/2] , optimizedpoint ]] ]
     , {j, 1, sourcedim, 1}];
   \[Phi]in = Total@Table[Table[E^(-((x - sourcelist[[j]])^2/(2  (\[Lambda]0/2)^2))) E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]0/2)^2))) E^(-I phases[[j]]), {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}], {j, 1, sourcedim}];(*Discretized source*)
   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] (*Solve the linear system*)
   , {50}];

\[Phi]0 = Table[ren = 
    ren0*Clip[Total[Table[RotateRight[DiskMatrix[5, dim], {RandomInteger[{-Round[dim/2], Round[dim/2] }], RandomInteger[{ -Round[dim/8], Round[dim/2] - 10}]}], {250}]], {0, 1}] + 1;
   n = ren + I imn;
   \[Phi]in = Total@Table[Table[E^(-((x - sourcelist[[j]])^2/(2  (\[Lambda]0/2)^2))) E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]0/2)^2))) , {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}], {j, 1, sourcedim}];(*Discretized source*)
   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] (*Solve the linear system*)
   , {50}];

frames = Table[GraphicsRow[{
     ArrayPlot[
      Transpose[(Mean[Abs[#]^2 & /@ \[Phi]0[[1 ;; j]]  ]/Max[avg\[Phi]0abs[[(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}, LabelStyle -> {Black, Bold}, ColorFunctionScaling -> False],
     ArrayPlot[
      Transpose[(Mean[Re[#] & /@ \[Phi]0[[1 ;; j]]  ]/Sqrt[Max[avg\[Phi]0abs[[(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 -> ReMapC, DataReversed -> True, Frame -> False, PlotRange -> {-1, 1}, LabelStyle -> {Black, Bold}, ColorFunctionScaling -> True]
     }]
   , {j, 1, 50}];
ListAnimate[frames]