File:Spherical Lens.gif

From Vigyanwiki

Spherical_Lens.gif ‎(543 × 543 pixels, file size: 6.8 MB, MIME type: image/gif, looped, 92 frames, 9.2 s)

Note: Due to technical limitations, thumbnails of high resolution GIF images such as this one will not be animated.

This file is from Wikimedia Commons and may be used by other projects. The description on its file description page there is shown below.

Summary

Description
English: Visualization of light through a spherical lens, as a function of the radii of curvature of the two facets. Notice that we are far from the "thin lens" approximation.
Date
Source https://twitter.com/j_bertolotti/status/1392058658520027138
Author Jacopo Bertolotti
Permission
(Reusing this file)		 https://twitter.com/j_bertolotti/status/1030470604418428929
GIF genesis
InfoField
 This diagram was created with Mathematica

Mathematica 12.0 code

\[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] = 40*\[Lambda]0; (*Parameters for the grid*)
\[Sigma] = 10 \[Lambda]0; (*width of the gaussian beam*)

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]}]; (*Discretized source*)
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[\[Phi]in][[1]];
L = -1/\[Delta]^2*KirchhoffMatrix[GridGraph[{dim, dim}]]; (*Discretized Laplacian*)

ycenter = Map[y0 /. # &, FullSimplify[Solve[(x1)^2 + (y1 - y0)^2 == r^2, {y0}]][[All, 1, All]]  ];

surface2[x_] := Evaluate[Evaluate[((Sqrt[r^2 - (x)^2] + y0) /. {y0 -> ycenter[[1]]}) /. {y1 -> -(\[CapitalDelta]/4), x1 -> \[CapitalDelta]/2, r -> 100 \[CapitalDelta]} ] ];
surface1[x_] := Evaluate[((-Sqrt[r^2 - (x)^2] + y0 - 1) /. {y0 -> ycenter[[2]]}) /. {y1 -> -(\[CapitalDelta]/4), x1 -> \[CapitalDelta]/2, r -> 100 \[CapitalDelta]}];
frames1 = Table[
  ren = Table[ If[y < Re@Evaluate[surface2[x]] && y > Re@surface1[x], n0, 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
  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*)
  \[Phi]s = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
  
  ImageAdd[
   ArrayPlot[ Transpose[(Abs[\[Phi]in + \[Phi]s]/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, 1}],
   ArrayPlot[Transpose@Re[(n - 1)/5] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
   ](*Plot everything*)
  , {n0, 1, 2.5, 0.24}];

surface2[x_] := Evaluate[Evaluate[((Sqrt[r^2 - (x)^2] + y0) /. {y0 -> ycenter[[1]]}) /. {y1 -> -(\[CapitalDelta]/4), x1 -> \[CapitalDelta]/2, r -> (c^2 + (c \[CapitalDelta])/2 + (5 \[CapitalDelta]^2)/16)/(
        2 (c + \[CapitalDelta]/4))} ] ];
surface1[x_] := Evaluate[((-Sqrt[r^2 - (x)^2] + y0 - 1) /. {y0 -> ycenter[[2]]}) /. {y1 -> -(\[CapitalDelta]/4), x1 -> \[CapitalDelta]/2, r -> 100 \[CapitalDelta]}];
frames2 = Table[
  ren = Table[ If[y < Re@Evaluate[surface2[x]] && y > Re@surface1[x], 2.5, 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
  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*)
  \[Phi]s = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
  
  ImageAdd[
   ArrayPlot[
    Transpose[(Abs[\[Phi]in + \[Phi]s]/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, 1}],ArrayPlot[Transpose@Re[(n - 1)/5] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
   ](*Plot everything*)
  , {c, -(\[CapitalDelta]/4) + 0.01, 0, \[CapitalDelta]/(20*3)}];

surface2[x_] := Evaluate[Evaluate[((Sqrt[r^2 - (x)^2] + y0) /. {y0 -> ycenter[[1]]}) /. {y1 -> -(\[CapitalDelta]/4), x1 -> \[CapitalDelta]/2, r -> 5/8 \[CapitalDelta]} ] ];
surface1[x_] := Evaluate[((-Sqrt[r^2 - (x)^2] + y0 - 1) /. {y0 -> ycenter[[2]]}) /. {y1 -> -(\[CapitalDelta]/4), x1 -> \[CapitalDelta]/2, r -> (c^2 + (c \[CapitalDelta])/2 + (5 \[CapitalDelta]^2)/16)/(2 (c + \[CapitalDelta]/4))}];
frames3 = Table[
  ren = Table[If[y < Re@Evaluate[surface2[x]] && y > Re@surface1[x], 2.5, 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
  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*)
  \[Phi]s = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
  
  ImageAdd[
   ArrayPlot[Transpose[(Abs[\[Phi]in + \[Phi]s]/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, 1}],
   ArrayPlot[Transpose@Re[(n - 1)/5] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
   ](*Plot everything*)
  , {c, -(\[CapitalDelta]/4) + 0.01, -\[CapitalDelta]/10, \[CapitalDelta]/(20*3)}];

surface2[x_] := Evaluate[Evaluate[((Sqrt[r^2 - (x)^2] + y0) /. {y0 -> ycenter[[1]]}) /. {y1 -> -(\[CapitalDelta]/4), x1 -> \[CapitalDelta]/2, r -> (c^2 + (c \[CapitalDelta])/2 + (5 \[CapitalDelta]^2)/16)/(2 (c + \[CapitalDelta]/4))} ] ];
surface1[x_] := Evaluate[((-Sqrt[r^2 - (x)^2] + y0 - 1) /. {y0 -> ycenter[[2]]}) /. {y1 -> -(\[CapitalDelta]/4), x1 -> \[CapitalDelta]/2, r -> 109/120 \[CapitalDelta]}];
frames4 = Table[
  ren = Table[If[y < Re@Evaluate[surface2[x]] && y > Re@surface1[x], 2.5, 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
  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*)
  \[Phi]s = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
  
  ImageAdd[
   ArrayPlot[Transpose[(Abs[\[Phi]in + \[Phi]s]/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, 1}], ArrayPlot[Transpose@Re[(n - 1)/5] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
   ](*Plot everything*)
  , {c, 0, -(\[CapitalDelta]/4) + 0.01, -(\[CapitalDelta]/(20*3))}];

surface2[x_] := Evaluate[Evaluate[((Sqrt[r^2 - (x)^2] + y0) /. {y0 -> ycenter[[1]]}) /. {y1 -> -(\[CapitalDelta]/4), x1 -> \[CapitalDelta]/2, r -> 100 \[CapitalDelta]} ] ];
surface1[x_] := Evaluate[((-Sqrt[r^2 - (x)^2] + y0 - 1) /. {y0 -> ycenter[[2]]}) /. {y1 -> -(\[CapitalDelta]/4), x1 -> \[CapitalDelta]/2, r -> (c^2 + (c \[CapitalDelta])/2 + (5 \[CapitalDelta]^2)/16)/(2 (c + \[CapitalDelta]/4))}];
frames5 = Table[
  ren = Table[If[y < Re@Evaluate[surface2[x]] && y > Re@surface1[x], 2.5, 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
  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*)
  \[Phi]s = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
  
  ImageAdd[
   ArrayPlot[Transpose[(Abs[\[Phi]in + \[Phi]s]/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, 1}], ArrayPlot[Transpose@Re[(n - 1)/5] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
   ](*Plot everything*)
  , {c, -\[CapitalDelta]/10, -(\[CapitalDelta]/4) + 0.01, -(\[CapitalDelta]/(20*3))}];

ListAnimate[ Flatten[ Join[Table[frames1[[1]],