i = 1.5 (slider from 1 to 3, step 0.05)
O = Intersect(xAxis, yAxis)
D = Point(yAxis)
upright = Ray(D, (x(D) + 1.0, y(D)))
T_4 = Point(upright)
xplus = Ray(O, (1.0, 0.0))
T_2 = Point(xplus)
a: (x(T_2)^(2.0) * (i^(2.0) - 1.0) + 2.0 * x(T_2) * x(T_4) * (1.0 - i^(2.0)) + x(T_4)^(2.0) * (i^(2.0) - 1.0) - y(T_4)^(2.0)) / (2.0 * (x(T_2) - x(T_4)) * (i^(2.0) - 1.0))
b: sqrt(a^(2.0) * (i^(2.0) - 1.0))
e: sqrt(a^(2.0) + b^(2.0))
F_2: (x(T_2) - a + e, 0.0)
F_1: (x(T_2) - (e + a), 0.0)
xminus = Ray(O, (-1.0, 0.0))
D' = Mirror(D, xplus)
T_1 = Point(xminus)
upleft = Ray(D, (x(D) - 1.0, y(D)))
T_3 = Point(upleft)
T'_3 = Mirror(T_3, xAxis)
T'_4 = Mirror(T_4, xAxis)
up = Segment(T_3, T_4)
down = Segment(T'_3, T'_4)
xF: (y(T_4)^(2.0) / 4.0 - x(T_2)^(2.0) + x(T_2) * x(T_4)) / (x(T_4) - x(T_2))
xL: 2.0 * x(T_2) - xF
F: (xF, 0.0)
L: (xL, 0.0)
l = OrthogonalLine(L, xAxis)
p = CircumcircleArc(T_4, T_2, T'_4)
plane = Segment(T_3, T'_3)
r_1 = Segment((x(T_3) - 2.0, y(T'_3)), (x(T_3) - 2.0, y(T_3)))
r_2 = Segment((x(T_3) - 4.0, y(T'_3)), (x(T_3) - 4.0, y(T_3)))
R_1 = Point(r_1)
R_2 = Point(r_2)
rayin = Ray(R_2, R_1)
I_1 = Intersect(rayin, plane)
n_1 = OrthogonalLine(I_1, plane)
c_1 = Circle(I_1, 1.0)
Ir_1, Ir_2 = Intersect(c_1, rayin)
In_1, In_2 = Intersect(c_1, n_1)
α_1 = Angle(In_1, I_1, Ir_1)
In'_2 = Rotate(In_2, asin(sin(α_1) / i), I_1)
segmentin = Segment(R_2, I_1)
raymid = Ray(I_1, In'_2)
I_2, I'_2 = Intersect(p, raymid)
t_2, c = Tangent(I_2, p)
n_2 = OrthogonalLine(I_2, t_2)
c_2 = Circle(I_2, 1.0)
In_3, In_4 = Intersect(c_2, n_2)
Ir_3, Ir_4 = Intersect(c_2, raymid)
α_2: -(Angle[Ir_3, I_2, In_3])
In'_4 = Rotate(In_4, -asin(sin(α_2) * i), I_2)
xm: (x(T_2)^(2.0) * (i^(2.0) - 1.0) + x(T_4)^(2.0) * (1.0 - i^(2.0)) + y(T_4)^(2.0)) / (2.0 * (x(T_2) - x(T_4)) * (i^(2.0) - 1.0))
p_1 = If(x ≥ x(T_4) ∧ x ≤ x(T_2), sqrt(((x - xm)^(2.0) / a^(2.0) - 1.0) * b^(2.0)))
p_2 = If(x ≥ x(T_4) ∧ x ≤ x(T_2), -sqrt(((x - xm)^(2.0) / a^(2.0) - 1.0) * b^(2.0)))
segmentmid = Segment(I_1, I_2)
rayout = Ray(I_2, In'_4)