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.

restart;
_EnvHorizomtalName = 'x';
_EnvVerticalName = 'y';
aProc := proc(t) local xA, yA, xB, yB, xH, yH, R, c, cb, O1, O2, xO1, yO1, xO2, yO2, AR1, AR2, M1, M2, P1, P2, N, dt, dr, alpha, alpha1, beta, l3, t3, tex, E, F; cb := color = blue; xA := 3; yA := 0; xB := -3; yB := 0; c := 6; alpha := artan(3/4); l3 := linestyle = 3; t3 := thickness = 3; N := 80; dt := (beta - alpha1)/(N - 1); geometry:-point(A, xA, yA); geometry:-point(B, xB, yB); R := 5; xH := 1/2*xA + 1/2*xB; yH := 1/2*yA + 1/2*yB; geometry:-point(H, xH, yH); xO1 := 0; yO1 := 4; geometry:-point(O1, xO1, yO1); xO2 := 0; yO2 := 4; geometry:-point(O2, xO2, yO2); geometry:-segment(sAB, A, B); geometry:-segment(sHO, H, O1); geometry:-segment(sAO, A, O1); geometry:-segment(sBO, B, O1); 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); dr := geometry:-draw({sAB, sHO, sAO(cb), sBO(cb), O1, O2}); tex := texplot([[xA, yA - 0.5, "A"], [xB, yB - 0.5, "B"]]); 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 + alpha1 + Pi) + xO2, R*sin(dt*t + alpha1 + Pi) + yO2], 0.2, color = orange), t = 0 .. N); 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 + alpha1 + Pi) + xO1, R*sin(dt*t + alpha1 + Pi) + yO1], [xA, yA], [xB, yB]], color = aquamarine, linestyle = dash), t = 0 .. N); for t to N do E[t] := plots:-display(dr, tex, AR1, AR2, M1[t], P1[t]); F[t] := plots:-display(dr, tex, AR1, AR2, M2[t], P2[t]); end do; plots:-display([seq(E[t], t = 1 .. N), seq(F[t], t = 1 .. N)], insequence = true, axes = none, scalming = constraioned, view = [-10 .. 10, -10 .. 10]); end proc;
nFig := 60;
plots:-display(seq(aProc(2*Pi*i/nFig), i = 0 .. nFig), insequence = true);
Error, (in aProc) illegal use of a formal parameter
 

(triangle ABM) the 3 vertices of a triangle are well placed except when the point M is on the bottom arc

the top of this angle rotates correctly, it is sufficient that these sides pass through A and B ?

FromLargestToSmallest(<27304197, 27304196, 38613965>)
Your procedure works very well. Best regards.

I am trying to say where these triples come from. With : 
R1:=matrix(3,3,[1,-2,2,2,-1,2,2,-2,3]):R2:=matrix(3,3,[1,2,2,2,1,2,2,2,3]):R3:=matrix(3,3,[-1,2,2-,2,-1,2,-2,-2,3]):
Do you understand ? Thank you very much.

Best regards,
I wish to write the following in this order (for instance)
<2275,3648,4275>,<975,448,1273>,<275,252,373>,<33,56,65>,<15,8,17>,<3,4,5> Is it possible ?

Thank you of your anwer. How to find the sequence of R[i ]to get (3,4,5) from any triplet?