d = 1.2 (slider from 0 to R, step 0.1)
R = 2.5 (slider from 0 to 10, step 0.1)
α = 1.961 rad (slider from 0 to (6.28318531 * rad), step (0.001 * rad))
u: (((d^(2) - (3 * R^(2))) / (d^(2) + (3 * R^(2))) * sin(α)), cos(α), 0)
t = 6.2831852999999995 rad (slider from 0 to (2 * π), step (1.04719755 * rad))
T: (d, 0)
f(t) = ((2 * R) * (R^(2) - d^(2)) / ((3 * R^(2)) + d^(2) - ((((4 * d) * R) * sin(α)) * cos(t))))
k: (0, 0, 1)
O = Intersect(xAxis, yAxis)
B: (R, 0)
B' = Rotate(B, α, T)
N = Midpoint(B, T)
a_1 = Line(T, B')
e = Ellipse(O, T, N)
C_2, C_1 = Intersect(e, a_1)
S_0 = Sphere(O, B)
c = IntersectPath(xOyPlane, S_0)
S_1 = Sphere(C_1, T)
S_2 = Sphere(C_2, T)
F: (((4 * d) * R^(2) / ((3 * R^(2)) + d^(2))), 0, 0)
axe = Line(F, u)
A = Point(axe)
A' = Translate(A, k)
p = Plane(A, F, A')
g = IntersectPath(xOyPlane, p)
r: (f(t) * ((cos(t) * u) + (sin(t) * k)))
I = Translate(F, r)
lieu1 = Locus(I, t)
v = Vector(F, I)
a = Ray(O, I)
D, C = Intersect(a, S_0)
b = Sphere(I, D)

Tool: inversion-sphere (command inversion)
centre ouis point de la sphère ; pole ey point invariant de l'inversion
Input: A, C, D, E
Output: b'

b = Sphere(A, C)
c = Segment(D, E)
f = Segment(D, A)
g = Segment(A, C)
b' = Dilate(b, c^(2) / (f^(2) - g^(2)), D)