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]

我，本作品的著作權持有者，決定用以下授權條款發佈本作品：

	此檔案在創用CC CC0 1.0 通用公有領域貢獻宣告之下分發。
	在此宣告之下分發本作品者，已依據各國著作權法，在全世界放棄其對本作品所擁有的著作權及所有相關相似的法律權利，從而將本作品貢獻至公有領域。您可以複製、修改、分發和演示該作品，用於任何商業用途，所有這些都不需要請求授權。 http://creativecommons.org/publicdomain/zero/1.0/deed.enCC0Creative Commons Zero, Public Domain Dedicationfalsefalse