R = 100.0; a = 5.0; r = Sqrt[a (2*R - a)];
n1 = n3 = 1; n2 = 1.52;
light = {}; \[Theta] = 0 \[Degree];
Do[
p1 = {-0.5 *R, y0}; k = n1*{Cos[\[Theta]], -Sin[\[Theta]]};
s = FindRoot[{k[[1]]*(y - p1[[2]]) == k[[2]] *(x - p1[[1]]),
R^2 - (y^2 + (x - R + a)^2) == 0}, {{x, -a}, {y, y0}}]；
p2 = {x, y} /. s;
f = R^2 - (y^2 + (x - R + a)^2);
\[CapitalOmega] = D[f, {{x, y}}] /. Thread[{x, y} -> p2];
\[CapitalOmega] = \[CapitalOmega]/Norm[\[CapitalOmega]];
G1 = k.\[CapitalOmega]; G2 = Sqrt[n2^2 - n1^2 + G1^2];
k = k + (G1 - G2) \[CapitalOmega];
s = FindRoot[{k[[1]] (y - p2[[2]]) == k[[2]] (x - p2[[1]]),
y^2 + (x + R - a)^2 - R^2 == 0}, {{x, a}, {y, y0}}];
p3 = {x, y} /. s;
f = y^2 + (x + R - a)^2 - R^2;
\[CapitalOmega] = D[f, {{x, y}}] /. Thread[{x, y} -> p3];
\[CapitalOmega] = \[CapitalOmega]/Norm[\[CapitalOmega]];
G1 = k.\[CapitalOmega]; G2 = Sqrt[n3^2 - n2^2 + G1^2];
k = k + (G1 - G2) \[CapitalOmega];
P4 = {x, k[[2]]/k[[1]] *(x - p3[[1]] + p3[[2]])} /. x -> 2 R;
AppendTo[light, Graphics[Line[{p1, p2, p3, p4}]]],
{y0, -0.1 r, 0.1 r, 0.05 r}]
\[Alpha] = ArcSin[r/R];
Show[light,
Epilog -> {Think,
Circle[{R - a, 0}, R, {\[Pi] - \[Alpha], \[Pi] + \[Alpha]}],
Circle[{-R + a, 0}, R, {-\[Alpha], \[Alpha]}]},
PlotRange -> {{-0.5 R, 1.5 R}, {-r, r}}，
AxesLabel -> {"x", "y"}]
Clear["Global`*"]