Despite your animations I can not answer these questions :
1.-Why 2 triangles CMP and ABH are similar ?
2.-At what curve E is the line MP tangent ?
Can you help me ?

It is a corollary of Poncelet’s theorem : The line that joins a focal spot F to a point P common in M and M' to an ellipse is bisector inside the angle MFM"

Do you want to tell me what is the interest of the calculation of alpha1 compared to alpha ? alpha can be drawn;

Problem with the 3 last lines of the code..
(sin(Setof[2]))^2
Error, invalid subscript selector
simpascn:=arcsin(sqrt(factor(  ))); # the t's eliminate
Error, Invalid label reference
 alpha1:=(subs(Setof[2]=simpascn,alpha))*180/Pi
Error, invalid subscript selector
plot(alpha1,t=-10..10,view=[-10..10,0..180])
Warning, expecting only range variable t in expression alpha1 to be plotted but found name alpha1
 

Thank you a lot for your code.
Howeven ir does not work with Maple 2021.

I show that alpha is constant using this procedure :
TANG := proc(t)
local alpha, xM, yM, M, P, Q, TT, PF2, QF2;
xM := a*cos(t);
yM := b*sin(t);
point(M, xM, yM);
line(TT, 1/49*x*xM + 1/25*y*yM - 1 = 0);
intersection(P, tang1, TT);
intersection(Q, tang2, TT);
line(PF2, [P, F2]); line(QF2, [Q, F2]);
alpha := FindAngle(PF2, QF2);
print(evalf(180 - 180*arctan(alpha)/Pi, 5));
end proc;
seq(TANG(i*Pi/10), i = 0 .. 5);

More cleaily how to show that alpha remains constant when t varies ? Thank you.
restart;
with(plots);
with(geometry);
_EnvHorizontalName := 'x';
_EnvVerticalName := 'y';
x0 := 100;
y0 := 40;
a := 7;
b := 5;
c := sqrt(a^2 - b^2);
ellipse(el, x^2/a^2 + y^2/b^2 - 1);
point(F1, -c, 0);
point(F2, c, 0);
eq := simplify((a^2 - x0^2)*(y - y0)^2 + (b^2 - y0^2)*(x - x0)^2 + 2*x0*y0*(x - x0)*(y - y0)) = 0;
sol := solve({eq}, {y});
line(tang1, op(sol[1]));
line(tang2, op(sol[2]));
sol2 := op(solve({op(sol[1]), x^2/a^2 + y^2/b^2 - 1 = 0}, {x, y}));
xM2 := rhs(sol2[1]);
yM2 := rhs(sol2[2]);
point(A, xM2, yM2);
sol3 := op(solve({op(sol[2]), x^2/a^2 + y^2/b^2 - 1 = 0}, {x, y}));
xM3 := rhs(sol3[1]);
yM3 := rhs(sol3[2]);
point(B, xM3, yM3);
line(Pol, [A, B]);
simplify(Equation(Pol));
isolate(%, y);
TANG := proc(t) local xM, yM; xM := a*cos(t); yM := b*sin(t); RETURN(expand(1/49*x*xM + 1/25*y*yM - 1 = 0)); end proc;
t := -0.15;
line(TT, TANG(t));
Equation(TT);
point(M, a*cos(t), b*sin(t));
intersection(P, tang1, TT);
intersection(Q, tang2, TT);
line(PF2, [P, F2]);
line(QF2, [Q, F2]);
alpha := FindAngle(PF2, QF2);
print(evalf(180 - 180*arctan(alpha)/Pi, 4));
display(textplot([[-c, 0, "F1"], [c, 0, "F2"], [coordinates(B)[], "B"], [coordinates(A)[], "A "], [coordinates(M)[], "M "], [coordinates(P)[], "P "], [coordinates(Q)[], "Q "]], align = {"above", 'right'}), draw([el(color = red), A(color = black, symbol = solidcircle, symbolsize = 16), PF2(color = brown), QF2(color = brown), B(color = black, symbol = solidcircle, symbolsize = 16), M(color = black, symbol = solidcircle, symbolsize = 16), P(color = black, symbol = solidcircle, symbolsize = 16), tang1(color = green), tang2(color = green), TT(color = green), F1(color = blue, symbol = solidcircle, symbolsize = 16), Q(color = blue, symbol = solidcircle, symbolsize = 16), F2(color = red, symbol = solidcircle, symbolsize = 16)]), axes = none, view = [-7 .. 18, -7 .. 12])

the second proposal is appropriate. Thank you.

AI helped me out; An equation of polar line is x0*x/a^2+y0*y/b^2=1.

There is a parametric equation for the 2 tangents from an external point (eq); It would be impossible to find a polar equation.

restart;
with(plots);
with(geometry);
_EnvHorizontalName := 'x';
_EnvVerticalName := 'y';
assume(1/(a^2*b^2) <> 0);
ellipse(el, x^2/a^2 + y^2/b^2 = 1);
eq := simplify((a^2 - x0^2)*(y - y0)^2 + (b^2 - y0^2)*(x - x0)^2 + 2*x0*y0*(x - x0)*(y - y0)) = 0;
sol := solve({eq}, {y});
line(tang1, op(sol[1]));
slope(tang1);
line(tang2, op(sol[2]));
slope(tang2);
sol2 := op(solve({op(sol[1]), x^2/a^2 + y^2/b^2 - 1 = 0}, {x, y}));
xM2 := rhs(sol2[1]);
yM2 := rhs(sol2[2]);
point(M2, xM2, yM2);
sol3 := op(solve({op(sol[2]), x^2/a^2 + y^2/b^2 - 1 = 0}, {x, y}));
xM3 := rhs(sol3[1]);
yM3 := rhs(sol3[2]);
point(M3, xM3, yM3);
line(Pol, [M2, M3]);
equation(Pol);#??
Can you help me to find equation(Pol) ? Thank you.

restart;
with(plots);
with(geometry);
_EnvHorizontalName := 'x';
_EnvVerticalName := 'y';
x0 := 10;
y0 := 9;
a := 7;
b := 5;
c := sqrt(a^2 - b^2);
ellipse(el, x^2/a^2 + y^2/b^2 - 1);
point(F1, -c, 0);
point(F2, c, 0);
point(P, x0, y0);
eq := simplify((a^2 - x0^2)*(y - y0)^2 + (b^2 - y0^2)*(x - x0)^2 + 2*x0*y0*(x - x0)*(y - y0)) = 0;
sol := solve({eq}, {y});
line(tang1, op(sol[1]));
slope(tang1);
line(tang2, op(sol[2]));
slope(tang2);
sol2 := op(solve({op(sol[1]), x^2/a^2 + y^2/b^2 - 1 = 0}, {x, y}));
xM2 = rhs(sol2[1]);
yM2 = rhs(sol2[2]);
point(M2, xM2, yM2);
sol3 := op(solve({op(sol[2]), x^2/a^2 + y^2/b^2 - 1 = 0}, {x, y}));
xM3 = rhs(sol3[1]);
yM3 = rhs(sol3[2]);
point(M3, xM3, yM3);
assume(xM2 - xM3 <> 0);
line(Pol, [M2, M3]);
display(textplot([[-c, 0, "F1"], [c, 0, "F2"], [coordinates(P)[], "P"]], align = {"above", 'right'}), draw([el(color = red), P(color = black, symbol = solidcircle, symbolsize = 16), tang1(color = green), tang2(color = green), Pol(color = red, thickness = 3), F1(color = blue, symbol = solidcircle, symbolsize = 16), F2(color = red, symbol = solidcircle, symbolsize = 16)]), axes = none);
Warning, data could not be converted to float Matrix
Do you want to tell me why polar is not tracing ? Thank you very much/

Yes you have ggod memory.
I really have a hard nipple

(y - y0 - m1) + (-x + x0);
(y - y0 - m2) + (-x + x0);
line(tang1, (y - y0 - m1) + (-x + x0) = 0);
line(tang2, (y - y0 - m2) + (-x + x0) = 0);
display(textplot*([-c, 0, "F1"], [c, 0, "F2"], [coordinates(P)[], "P"], align = {"above", 'right'}), draw(el(color = red), P(color = black, symbol = solidcircle, symbolsize = 16), tang1(color = green), tang2(color = green), F1(color = blue, symbol = solidcircle, symbolsize = 16), F2(color = red, symbol = solidcircle, symbolsize = 16), axes = none));
Error, invalid input: plottools:-line expects its 1st argument, x, to be of type list, but received tang1
Error, invalid input: plottools:-line expects its 1st argument, x, to be of type list, but received tang2
Error, (in geometry:-coordinates) wrong type of argument
how to solve problems ?

restart;
with(geometry);
with(plots);
_EnvHorizomtalName = 'x';
_EnvVerticalName = 'y';
_local(D);
line(delta, y = 1/3*x + 1, [x, y]);
line(deltap, y = (-1)/4*x - 2, [x, y]);
line(D, y = 3*x - 5, [x, y]);
point(S, -1, 0);
intersection(P, delta, D);
intersection(Pp, deltap, D);
intersection(N, delta, deltap);
projection(H, S, delta);
projection(K, S, deltap);
projection(M, S, D);
circle(c1, [H, K, M], 'centername' = O1);
coordinates(O1);
omega := FindAngle(delta, deltap);
line(SP, [S, P]);
line(SPp, [S, Pp]);
alpha := FindAngle(SP, SPp);
segment(sMH, M, H);
segment(sMK, M, K);
line(MH, [M, H]);
line(MK, [M, K]);
omega - alpha;
beta := FindAngle(MH, MK);
is(beta + omega - alpha = 0);
point(S1, 2*HorizontalCoord(O1) + 1, 2*VerticalCoord(O1));
coordinates(S1);
line(SS1, [S, S1]);
c := distance(S, S1)/2;
a := radius(c1);
b := sqrt(a^2 - c^2);
ell := x^2/a^2 + y^2/b^2 = 1;
ellip := implicitplot(ell, x = -2 .. 2, y = -2 .. 2, color = blue);
solve({Equation(SS1, [x, y])}, {y});
rotation(elli2, ellip, arccos(-3283/4253), 'clockwise', O1);
display(textplot([[coordinates(S)[], "S"], [coordinates(P)[], "P"], [coordinates(H)[], "H"], [coordinates(N)[], "N"], [coordinates(K)[], "K"], [coordinates(M)[], "M"], [coordinates(Pp)[], "P'"], [coordinates(S1)[], "S1"], [coordinates(O1)[], "O1"]], font = [times, bold, 16], align = [above, right]), draw([delta(color = green), deltap(color = blue), D(color = red), c1(color = black), MH(color = orange), MK(color = orange), SP(color = magenta), SPp(color = magenta), S(color = black, symbol = solidcircle, symbolsize = 18), H(color = black, symbol = solidcircle, symbolsize = 12), K(color = black, symbol = solidcircle, symbolsize = 12), M(color = black, symbol = solidcircle, symbolsize = 12), P(color = black, symbol = solidcircle, symbolsize = 12), S1(color = black, symbol = solidcircle, symbolsize = 12), O1(color = black, symbol = solidcircle, symbolsize = 12)]), ellip, scaling = constrained, axes = normal);
                         [-965  -3283]
                         [----, -----]
                         [5218  5218 ]

                                     /7 \
                      omega := arctan|--|
                                     \11/

                                     /276\
                      alpha := arctan|---|
                                     \61 /

                          /7 \         /276\
                    arctan|--| - arctan|---|
                          \11/         \61 /

                                    /2609\
                      beta := arctan|----|
                                    \2603/

                              true

                               S1

                         [1644  -3283]
                         [----, -----]
                         [2609  2609 ]

                              SS1

                    1             (1/2)        (1/2)
            c := -------- 28866098      6806881     
                 13613762                           

                   1             (1/2)         (1/2)
           a := -------- 34595401      13613762     
                13613762                            

                           2          (1/2)
                     b := ---- 2520294     
                          2609             

                       13613762  2   2609  2    
                ell := -------- x  + ---- y  = 1
                       34595401      3864       

                     /      3283   3283  \ 
                    { y = - ---- - ---- x }
                     \      4253   4253  / 

Error, (in geometry:-rotation) wrong type of arguments

NULL;I corrected the scale error of the previous program
I still haven’t understood my problem with display. Thank you very much.

restart;
with(geometry);
with(plots);
_EnvHorizontalName := 'x';
_EnvVerticalName := 'y';
cb := color = blue;
t3 := thickness = 3;
l3 := linestyle = 3;
xA := 3;
yA := 0;
xB := -3;
yB := 0;
c := 6;
point(A, xA, yA);
point(B, xB, yB);
R := 5;
alpha := arctan(3/4);
evalf(%*180/Pi);
                          36.86989765

xH := (xA + xB)/2;
yH := (yA + yB)/2;
point(H, xH, yH);
xO1 := 0;
yO1 := 4;
point(O1, xO1, yO1);
xO2 := 0;
yO2 := -4;
point(O2, xO2, yO2);
segment(sAB, A, B);
segment(sHO, H, O1);
segment(sAO, A, O1);
segment(sBO, B, O1);
                              sAB

alpha1 := arctan((yO1 - yA)/(xO1 - xA));
beta := Pi + arctan((yO1 - yB)/(xO1 - xB));
AR1 := plottools[arc]([xO1, yO1], R, alpha1 .. beta, l3);
AR2 := plottools[arc]([xO2, yO2], R, -beta .. -alpha1, l3);
N := 80;
dt := (beta - alpha1)/(N - 1);
dr := draw({O1, O2, sAB, sHO, sAO(cb), sBO(cb)});
tex := textplot([[xA, yA - 0.5, "A"], [xB, yB - 0.5, "B"], [xO1, yO1 + 0.5, "O1"], [xO2, yO2 - 0.5, "O2"], [xH, yH - 0.5, "H"]]);
M1 := seq(plottools[disk]([R*cos(dt*t + alpha1) + xO1, R*sin(dt*t + alpha1) + yO1], 0.2, color = orange), t = 0 .. N);
M2 := seq(plottools[disk]([R*cos(dt*t + Pi + alpha1) + xO2, R*sin(dt*t + Pi + alpha1) + yO2], 0.2, color = orange), t = 0 .. N);
point(m2, [R*cos(dt*t + Pi + alpha1) + xO2, R*sin(dt*t + Pi + alpha1) + yO2]);
P1 := seq(plottools[polygon]([[R*cos(dt*t + alpha1) + xO1, R*sin(dt*t + alpha1) + yO1], [xA, yA], [xB, yB]], color = aquamarine, linestyle = dash), t = 0 .. N);
P2 := seq(plottools[polygon]([[R*cos(dt*t + Pi + alpha1) + xO1, R*sin(dt*t + Pi + alpha1) + yO1], [xA, yA], [xB, yB]], color = aquamarine, linestyle = dash), t = 0 .. N);
triangle(ABM, [A, B, m2]);
AreCollinear: hint: could not determine if 24+30*sin(2/79*t*arctan(4/3)+1/79*t*Pi-arctan(4/3)) is zero
Error, (in geometry:-triangle) not enough information: the three points might be AreCollinear
for t to N do
    E[t] := display(dr, tex, AR1, AR2, M1[t], P1[t]);
    F[t] := display(dr, tex, AR1, AR2, M2[t], ABM);
end do;
Error, (in plots:-display) expecting plot structure but received: ABM
NULL;

display([seq(E[t], t = 1 .. N), seq(F[t], t = 1 .. N)], insequence = true, axes = none, scaling = constrained, view = [-10 .. 10, -10 .. 10]);
Error, (in plots:-display) expecting plot structures but received: [E[2], E[3], E[4], E[5], E[6], E[7], E[8], E[9], E[10], E[11], E[12], E[13], E[14], E[15], E[16], E[17], E[18], E[19], E[20], E[21], E[22], E[23], E[24], E[25], E[26], E[27], E[28], E[29], E[30], E[31], E[32], E[33], E[34], E[35], E[36], E[37], E[38], E[39], E[40], E[41], E[42], E[43], E[44], E[45], E[46], E[47], E[48], E[49], E[50], E[51], E[52], E[53], E[54], E[55], E[56], E[57], E[58], E[59], E[60], E[61], E[62], E[63], E[64], E[65], E[66], E[67], E[68], E[69], E[70], E[71], E[72], E[73], E[74], E[75], E[76], E[77], E[78], E[79], E[80], F[1], F[2], F[3], F[4], F[5], ...  F[71], F[72], F[73], F[74], F[75], F[76], F[77], F[78], F[79], F[80]]
NULL;
how to improve this program ? Thank you.