\[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/20; \[CapitalDelta] = 20*\[Lambda]0; (*Parameters for the grid*)
sourcef[x_, y_] := E^(I k0 y); (*Important! The source MUST be a solution of the Helmholtz equation in vacuum*)
\[Phi]in = Table[sourcef[x, y], {x, -\[CapitalDelta]/2, \[CapitalDelta]/ 2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/ 2, \[Delta]}]; (*Discretized source*)
d = \[Lambda]0/1; (*typical scale of the absorbing layer*)
imn = Table[5 I (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[\[Phi]in][[1]];
L = -1/\[Delta]^2*KirchhoffMatrix[GridGraph[{dim, dim}]]; (*Discretized Laplacian*)
r[t_] := (\[Lambda]0*3 - \[Delta]*5) Sin[\[Pi]/ 2 t]^2 + \[Delta]*5; (*The radius changes with the parameter t*)
frames = Table[
n = Table[
If[y^2 + x^2 <= r[t]^2, 2, 1], {x, -\[CapitalDelta]/ 2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/ 2, \[CapitalDelta]/2, \[Delta]}] + imn; (*Matrix with the refractive index*)
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]s = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
ImageAdd[
ArrayPlot[
Transpose[Abs[\[Phi]in + \[Phi]s]^2/ Max[Abs[\[Phi]in + \[Phi]s]^2]][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], ColorFunction -> "AvocadoColors" , DataReversed -> True,
Frame -> False, PlotRange -> {0, 0.8}], ArrayPlot[Transpose@Re[(n - 1)/10] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
](*Plot everything*)
, {t, 0, 1, 1./100}];
ListAnimate[frames]
You must be logged in to post a comment.