O = Intersect(xAxis, yAxis)
xminus = Ray(O, (-1.0, 0.0))
xplus = Ray(O, (1.0, 0.0))
D = Point(yAxis)
D' = Mirror(D, xplus)
T_1 = Point(xminus)
T_2 = Point(xplus)
upleft = Ray(D, (x(D) - 1.0, y(D)))
upright = Ray(D, (x(D) + 1.0, y(D)))
T_3 = Point(upleft)
T_4 = Point(upright)
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 = Parabola(F, l)
plane = Segment(T_3, T'_3)
i = 1.5 (slider from 1 to 3, step 0.05)
p_1 = If(sqrt((x - x(T_2)) * y(T_4)^(2.0) / (x(T_4) - x(T_2))) < y(T_4), sqrt((x - x(T_2)) * y(T_4)^(2.0) / (x(T_4) - x(T_2))))
p_2 = If(-sqrt((x - x(T_2)) * y(T_4)^(2.0) / (x(T_4) - x(T_2))) > -y(T_4), -sqrt((x - x(T_2)) * y(T_4)^(2.0) / (x(T_4) - x(T_2))))
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 = If(y(R_1) ≠ y(R_2), Tangent[I_2, p], Tangent[I'_2, p])
n_2 = If(y(R_1) ≠ y(R_2), OrthogonalLine[I_2, t_2], OrthogonalLine[I'_2, t_2])
c_2 = If(y(R_1) ≠ y(R_2), Circle[I_2, 1.0], Circle[I'_2, 1.0])
In_3, In_4 = Intersect(c_2, n_2)
Ir_3, Ir_4 = Intersect(c_2, raymid)
α_2 = If(y(R_1) ≠ y(R_2), Angle[Ir_3, I_2, In_3], Angle[Ir_3, I'_2, In_3])
In'_4 = If(y(R_1) ≠ y(R_2), Rotate[In_4, -asin(sin(α_2) * i), I_2], Rotate[In_4, -asin(sin(α_2) * i), I'_2])
segmentmid = If(y(R_1) ≠ y(R_2), Segment[I_1, I_2], Segment[I_1, I'_2])
rayout = If(y(R_1) ≠ y(R_2), Ray[I_2, In'_4], Ray[I'_2, In'_4])