n = 2; (*Refractive index of the bottom medium*)
z = 5; (*The bottom of the plot is at zero. This set how high you want to plot*)

h[x_] := 1 + 0.03 Cos[x] + 0.1 Sin[2.1 x] + 0.11 Sin[2.8 x]; (*profile of the refractive medium*)
dh[x_] :=  Evaluate[D[h[x], x]]; (*Derivative of the profile (that locally determines how Snell law works*)
R[x_] := x + (z - h[x]) dh[x] (n Sqrt[1 - (n^2 - 1) dh[x]^2] - 1)/(1 - (n^2 - 1) dh[x]^2); (*Where the ray ends up*)
Show[
 Plot[h[x], {x, 0, 2 \[Pi]}, Background -> Black, PlotStyle -> {Gray}, Filling -> Axis, PlotRange -> {0, z}, Axes -> False, FillingStyle -> Gray],
 Graphics[{
   Yellow, Thick, Opacity[0.1], Table[{Line[{{y, 0}, {y, h[y]}}], Line[{{y, h[y]}, If[Im[R[y]] == 0, {R[y], z}, {y, h[y]}]}]}, {y, 0, 2 \[Pi], 0.02}]
   }]
 , ImageSize -> Large
 ] (*Plot of the caustics*)

(*Generate the animation*)
listy =  RandomSample@Table[y, {y, 0, 2 \[Pi], 0.015}]; (*Randomize the order the rays appear*)
dim =  Dimensions[listy][[1]];
p1 = Table[
  Show[
   Plot[h[x], {x, 0, 2 \[Pi]}, Background -> Black, PlotStyle -> {Gray}, Filling -> Axis, PlotRange -> {0, z},  Axes -> False, FillingStyle -> Gray],
   Graphics[{
     Yellow, Thick, Opacity[0.15], Table[{Line[{{y, 0}, {y, h[y]}}], Line[{{y, h[y]}, If[Im[R[y]] == 0, {R[y], z}, {y, h[y]}]}]}, {y, listy[[1 ;; k]]}]
     }]
   , ImageSize -> Large ]
  , {k, 1, dim, 2}]; (*Plot every two frames, otherwise the animation is too heavy*)
ListAnimate[p1]

我，本作品著作权人，特此采用以下许可协议发表本作品：

	本作品采用知识共享CC0 1.0 通用公有领域贡献许可协议授权。
	采用本宣告发表本作品的人，已在法律允许的范围内，通过在全世界放弃其对本作品拥有的著作权法规定的所有权利（包括所有相关权利），将本作品贡献至公有领域。您可以复制、修改、传播和表演本作品，将其用于商业目的，无需要求授权。 http://creativecommons.org/publicdomain/zero/1.0/deed.enCC0Creative Commons Zero, Public Domain Dedicationfalsefalse