Question: Finding the location of the poles of a normal chord in an ellipse.

How to find the location of the poles of a normal chord in an ellipse ?
Here is my code :

restart; with(geometry); with(plots); `local`(O);
_EnvHorizontalName := x; _EnvVerticalName := y;
corde := a*x/cos(theta)-b*y/sin(theta) = a^2-b^2;
isolate(corde, a/cos(theta));
Error, (in isolate) a*x/cos(theta)-b*y/sin(theta) = a^2-b^2 does not contain a/cos(theta)
eq1 := (a^2-b^2)*X/a^2 = a/cos(theta);
c := solve(eq1, cos(theta));
eq2 := (a^2-b^2)*Y/b^2 = -b/sin(theta);
s := solve(eq2, sin(theta));
lieu := simplify(expand((a^2-b^2)^2*X^2*Y^2*(c^2+s^2 = 1)));
allvalues(eliminate({eq1, eq2}, theta))[1][2];
ell := x^2/a^2+y^2/b^2 = 1;
P := [a*cos(theta), b*sin(theta)];
tgP := x*P[1]/a^2+y*P[2]/b^2 = 1;

sol := solve({corde, ell}, {x, y});
tgP1 := simplify(x*rhs(sol[2][1])/a^2+y*rhs(sol[2][2])/b^2 = 1);

Drawing in a case
a := 5; b := 3; theta := (1/6)*Pi;

line(l1, corde); conic(co, ell);
Pole(P1, l1, co); coordinates(P1);
a := 5; b := 3; theta := (1/6)*Pi;
Ell := implicitplot(ell, x = -a .. a, y = -b .. b, color = red);
Cor := implicitplot(corde, x = -a-1 .. a, y = -b-1 .. b, color = blue);
TgP := implicitplot(tgP, x = 0 .. 10, y = -5 .. 10, color = magenta);
TgP1 := implicitplot(tgP1, x = -5 .. 10, y = -5 .. 10, color = magenta);
lieu := subs(X = x, Y = y, lieu);
subs(x = 125*sqrt(3)*(1/24), y = -27/8, lieu);
Lieu := implicitplot(lieu, x = -a .. a, y = -b .. b, color = green);
dr := draw(P1);
display([Ell, Cor, Lieu, TgP, TgP1, dr], axes = normal, view = [-10 .. 10, -10 .. 10], scaling = constrained);

Why the drawin of the location (lieu) does not appear ? Thank you.