資料來源

https://commons.wikimedia.org/wiki/File:Bessel_Beam_-_Self_healing.ogg

Methematica代碼：

\[Lambda]0 = 0.25; 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] = 60*\[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] = 20 \[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*)

f = 10;
sourcef1[x_, y_] :=(*\[ExponentialE]^(-(x^2/(2 \[Sigma]^2)))*) E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]0/2)^2))) E^(I k0 y) E^(-I k0/(2 f) x^2);
\[Phi]in1 = Table[Chop[sourcef1[x, y]], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}]; 
f2 = 2;
sourcef2[x_, y_] :=(*\[ExponentialE]^(-(x^2/(2 \[Sigma]^2)))*) E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]0/2)^2))) E^(I k0 y) E^(-I k0/(2 f2) Abs[x]);
\[Phi]in2 = Table[Chop[sourcef2[x, y]], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}]; (*Discretized source*)
frames = Table[
   ren = Table[ If[-(\[CapitalDelta]/50) < x + shift < \[CapitalDelta]/50 && -\[CapitalDelta]/400 < y + \[CapitalDelta]/4 < \[CapitalDelta]/400, 2 I, 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
   n = ren + I imn;
   M = L + DiagonalMatrix[SparseArray[Flatten[n]^2 k0^2]]; (*Operator on the left-hand side of the equation we want to solve*)
   b = -(Flatten[n]^2 - 1) k0^2 Flatten[\[Phi]in1]; (*Right-hand side of the equation we want to solve*)
   \[Phi] = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
   b2 = -(Flatten[n]^2 - 1) k0^2 Flatten[\[Phi]in2]; (*Right-hand side of the equation we want to solve*)
   \[Phi]2 = Partition[LinearSolve[M, b2], dim]; (*Solve the linear system*)
  
    Grid[{{
      Show[ImageAdd[
        ArrayPlot[
         Transpose[(Re@\[Phi]/Max[Abs@\[Phi][[(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 -> {-0.5, 0.5}, LabelStyle -> {Black, Bold}, ColorFunctionScaling -> True, ColorFunction -> ReMapC, ClippingStyle -> {Blue, Red}],
        ArrayPlot[Transpose[Im@ren/1] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
        ], PlotLabel -> "Gaussian beam", LabelStyle -> {Black, Bold, FontSize -> 26}, ImageSize -> Medium]
      ,
      Show[ImageAdd[
       ArrayPlot[
        Transpose[(Re@\[Phi]2/Max[Abs@\[Phi]2[[(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 -> {-0.5, 0.5}, LabelStyle -> {Black, Bold}, ColorFunctionScaling -> True, ColorFunction -> ReMapC, ClippingStyle -> {Blue, Red}]
        ,
        ArrayPlot[Transpose[Im@ren/1] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
        ], PlotLabel -> "Bessel beam", LabelStyle -> {Black, Bold, FontSize -> 26}, ImageSize -> Medium]
      }}]
   , {shift, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[CapitalDelta]/50}];
ListAnimate[frames]