# The Mirror curves created by Maple
restart;
E2 := solve({C[y] - P[y] = k[2]*(C[x] - P[x]), 
             C[y] - Y(t) = k[1]*(C[x] - X(t))}, {C[x], C[y]}):
assign(E2):
k[1] := diff(Y(t), t)/diff(X(t),t):
k[2] := -1/k[1]:
M[x] := collect(normal(P[x] + 2*(C[x] - P[x])),diff);
M[y] := collect(normal(P[y] + 2*(C[y] - P[y])),diff);
X(t) := t:  Y(t) := t^2/2: 
P[x]:= 0:
m[x] := simplify(M[x]):  m[y] := simplify(M[y]):
SPy := [seq(i/2, i = -6..6)]:
p1 := plot({seq([m[x], m[y], t=-4..4], P[y] = SPy)},
            color = black, thickness = 1):
p2 := plot([X(t), Y(t), t = -4..4], view = [-4 .. 4, -4 .. 4], 
           thickness = 2):
with(plots):
display({p1, p2}, scaling = constrained);
P[y] := 'P[y]':
plot3d(subs(P[y] = Py,[P[y], m[x], m[y]]), Py = -3..3, 
       t = -3..3, axes = boxed, orientation = [15,30], 
       labels = ['Py', 'Mx', 'My'], grid = [40, 40] );

# The Inverse Problem
X(t) := 'X(t)': Y(t) := 'Y(t)': 
P[x] := 'P[x]': P[y] := 'P[y]':
C1[x] := (X1(t) + P[x])/2:  
C1[y] := (Y1(t) + P[y])/2:
C2[x] := (X2(t) + P[x])/2:  
C2[y] := (Y2(t) + P[y])/2:
k[1] := -1*(P[x] - X1(t))/(P[y] - Y1(t)):
k[2] := -1*(P[x] - X2(t))/(P[y] - Y2(t)): 
E1 := solve({Theta - C1[y] = k[1]*(Xi - C1[x]), 
      Theta - C2[y] = k[2]*(Xi - C2[x])}, {Xi, Theta}):
assign(E1):
X1(t) := X(t) - dX(t)*dt: Y1(t) := Y(t) - dY(t)*dt: 
X2(t) := X(t) + dX(t)*dt: Y2(t) := Y(t) + dY(t)*dt: 
Xi := collect(normal(Xi),[dt,dX(t),dY(t),X(t),Y(t)]);  
Theta := collect(normal(Theta),[dt,dX(t),dY(t),Y(t),X(t)]);  
Xi := limit(Xi, dt = 0):  Theta := limit(Theta, dt = 0):
dX(t) := diff(X(t), t):  dY(t) := diff(Y(t), t):
Xi := collect(Xi, diff);
Theta := collect(Theta, diff);

# The Circle as the Mirror Curve
X(t) := cos(t):  Y(t) := sin(t):  P[y] := 0:
xi := simplify(Xi);  theta := simplify(Theta);
X(t) := xi:  Y(t) := theta:
m[x] := simplify(M[x]);   m[y] := simplify(M[y]);
E1 := simplify(subs(t = solve(xi = X, t), 
                Y^2  = theta^2), symbolic):
E2 := collect(4*E1, [X^2, X ,P[x]]):
E3 := factor((E2)):
E4 := normal(map(t->t/(P[x]^2-1), subs(X = (Xs+P[x])/2, E3))):
E5 := subs(Xs = 2*X-P[x] ,map(t->t+Xs^2, -expand(E4)));
X(t) := cos(t):  Y(t) := sin(t):
P[y] := 0:   P[x] := 'P[x]':
SPx:=[seq(i/10,i=0..9)]:
p1 := plot({seq([xi, theta, t = 0..2*Pi], P[x] = SPx)},
            color = black, thickness = 0):
P[x] := 'P[x]': 
SPx := [seq(1+i/10, i = 1..10)]:
PxR := arccos(1/P[x])*99/100:
p2r := plot({seq([xi, theta, t = -PxR..PxR], P[x]=SPx)},
            color = black, linestyle = 3):
P[x] := 'P[x]': PxL := arccos(1/P[x])*101/100: 
p2l := plot({seq([xi, theta, t = PxL..2*Pi - PxL], 
             P[x]=SPx)}, color = black, linestyle = 3):
p3 := plot([X(t), Y(t), t = 0..2*Pi], 
            view = [-1.5 ..2, -1.5.. 1.5], thickness = 2):
display({p1, p2r, p2l, p3}, scaling = constrained);

# The Line as a Mirror Circle
X(t) := t:  Y(t) := 0:  P[x] := 0:  P[y] := 'P[y]':
xi := simplify(Xi);  theta := simplify(Theta);
X(t) := xi: Y(t) := theta:
m[x] := simplify(M[x]);  m[y] := simplify(M[y]);
E1 := subs(xi = X, Y = theta);
E2 := Y = X:  E3 := Y =- X:
solve({E1, E2}, {X,Y});  solve({E1, E3}, {X,Y});
SPy := [seq(i/4, i = 1..16)]:
p1 := plot({seq([xi, theta, t = -6..6], P[y] = SPy)}, 
              color = black, thickness = 1):
p2 := plot([[-6, 6], [0, 0], [6, 6]], thickness = 2):
display({p1, p2}, view = [-6..6, 0..8], scaling = constrained);

# The Continuous Deformation
a:=1: 
b:=5/8:
rho := sqrt(a^2*cos(t)^2 + b^2*sin(t)^2):
X(t) := rho*cos(t):  
Y(t) := rho*sin(t):
p1 := plot([X(t), Y(t), t = -Pi..Pi], thickness = 3):
P[x] := 'P[x]':  
P[y] := 0:
xi := simplify(Xi,symbolic):  
theta := simplify(Theta,symbolic):
SPx := [seq(i/10, i = -9..9)]:
p2 := plot({seq([xi, theta, t = -Pi..Pi], P[x] = SPx)},
            thickness = 1, color = black):
display({p1, p2}, axes = boxed, scaling = constrained);
P[x] := 0:  
P[y] := 'P[y]':
xi := simplify(Xi,symbolic): 
theta := simplify(Theta, symbolic):
SPy := [seq(i/10, i =-6..6)]:
p2 := plot({seq([xi, theta, t= -Pi..Pi], P[y] = SPy)},
            thickness = 1, color = black):
display({p1, p2}, axes = boxed, scaling = constrained);
P[x] := cos(tau)/2: 
P[y]:=sin(tau)/2:
xi := simplify(Xi, symbolic): 
theta := simplify(Theta, symbolic):
plot3d([xi, theta, tau] ,t = -Pi..Pi, tau = -Pi..Pi , 
          grid = [90, 90], axes = none, orientation = [-110, 90]);