Fig := proc(t)
              local M1,M2,m,L1,L3,L4,tang1,tang2,C1,P,T1,xM2,yM2;
              global a, b, c,el,F1,F2,OO,cir;
              point(M1, a*cos(t), b*sin(t)); 
              line( L1, [F1, M1]): 
              EQ([a*cos(t),b*sin(t)],[-c,0]):  
              solve(`%`,y):
              m:=coeff(`%`,x):   
              line(L3,y=m*(x-c)): 
              op(solve({x^2/a^2 + y^2/b^2 - 1 =0, y=m*(x-c)},{x,y},explicit)[2])[1]:
              xM2:=rhs(`%`):  
              op(solve({x^2/a^2 + y^2/b^2 - 1 =0, y=m*(x-c)},{x,y},explicit)[2])[2]:
              yM2:=rhs(`%`):  
              point(M2,xM2,yM2): 
              line( tang1,x*xM2/a^2 + y*yM2/b^2 = 1):  
              line( tang2,x*cos(t)/a + y*sin(t)/b = 1): 
              intersection(C1,L1,tang1): 
              midpoint(P,M2,C1): 
             #intersection(P,tang1,tang2): 
              line(L4,[F1,M2]):
              triangle(T1, [C1,F1,M2]): 
 display
  ( [ textplot
      ( [ [ -c, 0, "F1"],
          [ c,  0, "F2"] ,
          [ coordinates(P)[], "P"],                          
          [ coordinates(M1)[], "M1"],                          
          [ coordinates(M2)[], "M2"],                         
          [ coordinates(C1)[], "C1"]
        ],
        align={"above",'right'}
      ),
      draw
      ( [ el(color=red),
          M1(color=black,symbol=solidcircle, symbolsize=16),
          M2(color=black,symbol=solidcircle, symbolsize=16),                      
          L1( color=black),
          L3(color=green),
          L4(color=green),                      
          tang1(color=red),  tang2(color=red),cir(color=black),
          P(color=blue, symbol=solidcircle, symbolsize=16),   
          F1(color=blue, symbol=solidcircle, symbolsize=16),                      
          F2(color=red, symbol=solidcircle, symbolsize=16),
          C1(color=black, symbol=solidcircle, symbolsize=16),
          T1(color=blue, filled=true, transparency=0.95)
        ],
        axes=none
      )
    ]
  );
end:
      
  
 nFig := 180.0:

   Figs := seq(Fig(2*Pi*i/nFig), i = 0 .. nFig-1.0):
   display(Figs, insequence = true);
I eliminated an intersection but there is still a problem with P, How to solve it. Thank you.

 Fig := proc(t)
              local M1,M2,m,L1,L3,L4,tang1,tang2,C1,P,T1,xM2,yM2;
              global a, b, c,el,F1,F2,OO,cir;
              point(M1, a*cos(t), b*sin(t)); 
              line( L1, [F1, M1]): 
              EQ([a*cos(t),b*sin(t)],[-c,0]):  
              solve(`%`,y):
              m:=coeff(`%`,x):   
              line(L3,y=m*(x-c)): 
              op(solve({x^2/a^2 + y^2/b^2 - 1 =0, y=m*(x-c)},{x,y},explicit)[2])[1]:
              xM2:=rhs(`%`):  
              op(solve({x^2/a^2 + y^2/b^2 - 1 =0, y=m*(x-c)},{x,y},explicit)[2])[2]:
              yM2:=rhs(`%`):  
              point(M2,xM2,yM2): 
              line( tang1,x*xM2/a^2 + y*yM2/b^2 = 1):  
              line( tang2,x*cos(t)/a + y*sin(t)/b = 1): 
              intersection(C1,L1,tang1): 
              intersection(P,tang1,tang2): line(L4,[F1,M2]):
              triangle(T1, [C1,F1,M2]): 
      draw
      ( [ el(color=red),
          M1(color=black,symbol=solidcircle, symbolsize=16),
          M2(color=black,symbol=solidcircle, symbolsize=16),                      
          L1( color=black),
          L3(color=green),
          L4(color=green),                      
          tang1(color=red),  tang2(color=red),cir(color=black),
          P(color=blue, symbol=solidcircle, symbolsize=16),   
          F1(color=blue, symbol=solidcircle, symbolsize=16),                      
          F2(color=red, symbol=solidcircle, symbolsize=16),
          C1(color=black, symbol=solidcircle, symbolsize=16),
          T1(color=blue, filled=true, transparency=0.95)
        ],
        axes=none
      ):end:
  
 nFig := 180.0:

   Figs := seq(Fig(2*Pi*i/nFig), i = 0.1 .. nFig-1.0):
   display(Figs, insequence = true);
intersection: two given lines are ParallelLine, no intersection
Error, (in geometry:-draw) unknown geometric object  P
Error, (in plots:-display) expecting plot structure but received: Figs
How to overcome this difficulty? Thank a lot.

I am sorry; The program does not work With messages : 
Error, (in geometry:-draw) improper op or subscript selector
Error, (in plots:-display) expecting plot structure but received: Figs

Is there any way how to suppress the q on p jumps ? Thank you.

Fig := proc(t) local M, l1, q; global a, b, p; point(M, a*cos(t), b*sin(t)); line(l1, y = (Student[Calculus1]):-Tangent(solve(Equation(p), y)[2], x = HorizontalCoord(M))); reflection(q, p, l1); draw([p(color = blue), q(color = blue), l1(color = black), M(color = red, symbol = solidcircle, symbolsize = 16)], axes = normal, scaling = constrained); end proc
nFig := 60;
Figs := seq(Fig(2*Pi*i/nFig), i = 1 .. nFig);
Error, (in Student:-Calculus1:-Tangent) the slope is not defined at the point `x` = -7
display(Figs, insequence = true);
Error, (in plots:-display) expecting plot structure but received: Figs

Iam sorry; How to manage this error. Thank you.
 

I want to trace the tangent in P to p 

a := 7;
b := 3;
t := (2*Pi)/3;
NULL;
ellipse(p, x^2/a^2 + y^2/b^2 - 1, [x, y]);
point(P, a*cos(t), b*sin(t));
NULL;
rotation(q, p, t, 'clockwise', P);
tangentpc(Tg, P, p);
Error, (in geometry:-tangentpc) wrong type of arguments

draw([p(color = blue), P, q(color = blue),Tg(color=red), axes = normal, scaling = constrained);
Why this error message.? Thank you.

how to highlight that the sum of the sides is constant ? Thank you.
 Fig := proc(t) 
              local xD,yD, D, C,Points,tex; 
              global A, B, b, Omega1, EL1, EL2;
              xD:=Omega1[1]+aa*cos(t):
              yD:=bb*sin(t):
              D:= [xD, yD]; C := [xD + b, yD];
              Points := pointplot([A[], B[], C[], C[], D[], E[], Omega1[]], symbol = solidcircle, color = [red],     symbolsize = 6): 
              tex := textplot([[A[], "A"], [B[], "B"], [C[], "C"], [D[], "D"], [E[], "E"], [Omega1[], "Ω1"]], align = ["above", "right"]):
              display([polygonplot([A, B, C, D], color = blue, filled = true, transparency = 0.9), Points,tex,EL1, EL2], axes = normal, scaling = constrained); 
         end proc:
Fig((3*Pi)/4);
display([seq(Fig((2*Pi*i)/40), i = 1 .. 80)], insequence = true):
 

Thank you for your answer , And without geometry package ?

I still have a problem
 restart:
  with(plots):
  with(plottools):
  _EnvHorizontalName:= 'x':
  _EnvVerticalName:= 'y':
  EqBIS:= proc(P, U, V) 
               local a, eq, eq1, M1, t, PU, PV, bissec1:
               uses geometry, LinearAlgebra:
               description "P est le sommet de l'angle dont on chercche la bissectrice" :
               a:= (P - U)/Norm(P - U, 2) + (P - V)/Norm(P - V, 2): 
               M1:= P + a*t:
               eq1:= op(eliminate({x = M1[1], y = M1[2]}, t)): 
               eq:= op(eq1[2]):
               line(bis, eq):
               line(lv, x = 2*sqrt(r1*R)):
               intersection(Omega, bis, lv):
               return coordinates(Omega);
         end proc:
  r1:= 1/2:
  r2:= r1/2:
  R:= r1*(21 - 12*sqrt(3)):
  R:= r1*(21 - 12*sqrt(3)):
  a:= arc([0, 0], 2*r1, Pi/6 .. (5*Pi)/6):
  b:= arc([0, 0], r1, Pi/6 .. (5*Pi)/6):
  co:= EqBIS(<sqrt(3)/2, 1/2>, <0, 0>, <0, 1/2>):
  A:=[cos(Pi/6), sin(Pi/6)]:
  B:=[cos(5*Pi/6), sin(5*Pi/6)]:
  Oo:=[0,0]:
  Op:=[0,1/2]:
  poly:=[A,B,Oo]:
  Omega:= [co[1], co[2]]:
  Omega1:= [-co[1], co[2]]:
  r3:= 3/16:
  EF:= sqrt(r3):
  r:= (150 - 72*sqrt(3))/193*1/2:
  alpha:= -5/3*r + 1/2*1/2:
  p:= sqrt(3)/3*1/2 - sqrt(3)/18*r:
  p2:= textplot([[A[], "A"], [B[], "B"], [Oo[], "O"]], align = ["above", "right"]):
  display( [ a,
             b,
             seq(plot([[1/2*cos(Pi/6+i*1/5), 1/2*sin(Pi/6+i*1/5)],[0,0]],i=1..10),color=black),
             sector([0,0], 1/2, Pi/6.. 5*Pi/6, color =white,transparency=0.3),
             polygonplot(poly, thickness = 3, color = blue,transparency=0.7),
             disk(Omega, R, color = blue),
             disk(Omega1, R, color = blue),
             disk([0, 3/4], 1/4, color = yellow),
             disk([EF, 1/2 + r3], r3, color = green), disk([-EF, 1/2 + r3], r3, color = green),
             circle([-EF, 1/2 + r3], r3, color = black),
             circle([p, 3/4 + alpha], r, color = black),disk([p, 3/4 + alpha], r, color = red),
             circle([-p, 3/4 + alpha], r, color = black),disk([-p, 3/4 + alpha], r, color = red)
          ],
          axes = none,
          scaling = constrained,
          size = [500, 500]
        );
Warning, data could not be converted to float Matrix, I don't understand this message. Thank you for your help.
 

Thank you for your answer; I found : 

restart;
with(plots):
with(geometry):
_EnvHorizontalName := x:
_EnvVerticalName := y:
R := 11:
r := 7:
a := sqrt(R*r):

b := 2:
circle(C1, [point(P1, [0, 0]), R]):
circle(C2, [point(P2, [R + 2*b + r, 0]), r]):
ellipse(p, (x - R - b)^2/b^2 + y^2/a^2 = 1):
draw([C1(color = yellow, filled = true), 
C2(color = red, filled = true), p(color = blue, filled = true), 
C1(color = black), C2(color = black), p(color = black)], 
axes = none, view = [-15 .. 35, -15 .. 15], scaling = constrained):
alpha := arctan((R - r)/(R + 2*b + r));
long := cos(alpha)*(R + 2*b + r);
evalf(%);
circle(C2, [point(P2, [long, r - R]), r]);
rotation(p1, p, alpha, 'clockwise');
detail(p1);
point(A, 0, -R);
point(B, long, -R);
line(L1, [A, B]);
point(cen, [(143*sqrt(5))/25, -(26*sqrt(5))/25]);
line(L, [P1, P2]);
reflection(L2, L1, L);
detail(L2);
draw([C1(color = yellow, filled = true), C2(color = red, filled = true), p1(color = blue, filled = true), C1(color = black), C2(color = black), p1(color = black), L1(color = black), L2(color = black)], axes = none, view = [-15 .. 35, -15 .. 15], scaling = constrained);
The drawing is a sangatku.

restart;
with(plots):
with(geometry):
_EnvHorizontalName := x:
_EnvVerticalName := y:
R := 7:
point(A, [0, R]):
line(L1, y = sqrt(3)*x + R):
line(L2, y = -sqrt(3)*x + R):
line(L3, y = R/3):
intersection(B, L1, L3):
intersection(C, L2, L3):

triangle(ABC, [A, B, C]):
circle(C1, [point(P1, [0, 0]), R]):
circle(C2, [point(P2, [0, R/3 + (2*R)/9]), (2*R)/9]):

center(C2), coordinates(center(C2)):
circle(C3, [reflection(P3, center(C2), C), radius(C2)]):
circle(C4, [reflection(P4, center(C2), B), radius(C2)]):
circle(C5,[reflection(P5, center(C3), L2), radius(C2)]):

detail(C5):
circle(C6, [point(P5, [0, R/3 - (2*R)/9]), radius(C2)]):
circle(C7, [point(P7,[0, R/3 - (2*R)/9-(4*R)/9]), radius(C2)]):
circle(C8, [point(P8,[0, R/3 - (2*R)/9-(8*R)/9]), radius(C2)]):

circle(Cir, [point(Po, [0, -R/3 ]), (2*R)/3]):
draw([L1(color = blue), 
ABC(color = red, transparency = 0.5, filled = true), 
L2(color = blue), L3(color = blue), 
C1(color = blue, thickness = 3), 
C1(color = yellow, transparency = 0.8, filled = true),
C2(color = blue, filled = true),
C3(color = blue, filled = true),
C4(color = blue, filled = true),
C5(color = blue, filled = true),
C6(color = blue, filled = true),
C7(color = blue, filled = true),
C8(color = blue, filled = true),
Cir(color = red,transparency = 0.5,filled=true)], 
axes = normal,
 view = [-R .. R, -R .. R], 
scaling = constrained);
I think a troucle with a symetry in respect to L2

I try differently
rotation(pol2,pol1,Pi/3,clockwise,o):
Error, (in geometry:-rotation) wrong type of arguments   Still that kind of error.

  restart:
  with(plots): with(geometry):
  _EnvHorizontalName := 'x':
  _EnvVerticalName := 'y':
  a := 7:
  point(E, 0, a*sqrt(3)/2):
  point(B, -a/2, 0):
  point(C, a/2, 0):
  point(o, 0, a*sqrt(3)/6):
  point(A, 0, a/2):
  point(H, 0, 0):
  R := (3-sqrt(3))*sqrt(2)*a/12:
  point(J, 0, a*sqrt(3)/6 - R):
  line(L1, -(7*x*sqrt(3))/6 + (7*y)/2 - (49*sqrt(3))/12 = 0):
  reflection(J1, J, L1):coordJ1:=coordinates(J1):
  
  line(L2, -(7*x*sqrt(3))/6 - (7*y)/2 + (49*sqrt(3))/12 = 0):
  reflection(J2, J, L2):coordJ2:=coordinates(J2):
  triangle(Tr1, [E, B, C]):
  triangle(Tr2, [A, B, C]):
  StretchRotation(E1, E, B, Pi/4, clockwise, sqrt(2)/2):
  coordinates(E1):
  StretchRotation(E2, E, C, Pi/4, counterclockwise, sqrt(2)/2):
  coordinates(E2):
  triangle(Tr3, [E, B, E1]):
  triangle(Tr4, [E, C, E2]):
  triangle(Tr5, [B, C, J]):
  circle(cir, [point(P1,[0,a*sqrt(3)/6]), R]):
  poly1 := Matrix( [ [0, a*sqrt(3)/2], 
                     coordJ1,
                    [0, a/2],
                    coordJ2
                  ],
                  datatype = float
                ):
  pol1 := polygonplot(poly1, colour = green, transparency = 0.7, gridlines):
  #poly2 := Matrix([ coordinates(J),coordinates(E1), coordJ2,|-a/2,0]],
                  datatype = float
                ):
  #pol2 := polygonplot(poly2, colour = green, transparency = 0.7, gridlines):

  tex := textplot([[0, a*sqrt(3)/2, "E"],[0,a*sqrt(3)/6 - R,"J"],
                    [-7/4 + (7*sqrt(3))/4, -7/4 + (7*sqrt(3))/4,"E1"],
                    [7/4 - (7*sqrt(3))/4, -7/4 + (7*sqrt(3))/4,"E2"],
                    [-7/4 + (((7*sqrt(3))/6 - 7*(3 - sqrt(3))*sqrt(2)/12)*sqrt(3))/2, 
                    (7*sqrt(3))/6 + 7*(3 - sqrt(3))*   sqrt(2)/24,"J1"],
                    [7/4 - (((7*sqrt(3))/6 - 7*(3 - sqrt(3))*sqrt(2)/12)*sqrt(3))/2, 
                    (7*sqrt(3))/6 + 7*(3 - sqrt(3))*   sqrt(2)/24,"J2"]], 'align' = {'above', 'right'}):
  display( [ draw( [ Tr1(color = cyan),
                     Tr3(color = green),
                     cir(color=blue),
                     Tr2(color = red),
                     Tr4(color = grey),
                     Tr1(color=blue)],'view' = [-5 .. 5, 0 .. 7
                   ], 
                   axes = normal,
                   scaling = constrained,
                   size=[800,800]
                   ),
            tex
           ]
          );

Error, `)` unexpected

  line(L1, [B, o]):
  Equation(L1):
  line(L1, -(7*x*sqrt(3))/6 + (7*y)/2 - (49*sqrt(3))/12 = 0):
  reflection(J1, J, L1):coordinates(J1):
  triangle(Tr6, [B, J1, E]):
  line(L2, [C, o]):
  Equation(L2):
  line(L2, -(7*x*sqrt(3))/6 - (7*y)/2 + (49*sqrt(3))/12 = 0):
  reflection(J2, J, L2):coordinates(J2):
  triangle(Tr7, [C, J2, E]):
  triangle(T1, [E, J1, A]):
  triangle(T2, [E, C, E2]):
  triangle(T3, [B, H, J]):
  triangle(T4, [C, H, J]):

  display( [ draw( [ cir(color = orange, filled = true, transparency = 0.1),
                     Tr6(color = blue, filled = true, transparency = 0.2),
                     Tr5(color = blue, filled = true, transparency = 0.2),
                     Tr7(color = blue, filled = true, transparency = 0.2)
                     
                   ],
                   axes = none,
                   scaling = constrained
                 ),
             pol1
           ]
         );
very difficult to find where I made mistakes; Thank you to help me.

Ihave still errors.
restart:
  with(geometry):
  _EnvHorizontalName := x:
  _EnvVerticalName := y:
  b := a*(1/2 + 1/6*sqrt(45 - 24*sqrt(3)))^2:
  r := b*sqrt(b)/(sqrt(a + b) + sqrt(a)):
  assume(a > 0);
  ellipse(E, x^2/a^2 + y^2/b^2 = 1):
  circle(C3, [point(P3, [(1 + sqrt(3))*r, -r]), r]):
  simplify(Equation(C3)):
  sol:= solve( [Equation(C3), Equation(E), x>0, y<0]):
  point(PP, [eval(x, sol), eval(y,sol)]):
  simplify(coordinates(PP));
 
Warning, solve may be ignoring assumptions on the input variables.
Error, invalid input: eval expects its 2nd argument, eqns, to be of type {integer, equation, set(equation)}, but received {x = ((3/13)*3^(1/2)*(15-8*3^(1/2))^(1/2)-(5/26)*3^(1/2)+(11/26)*(15-8*3^(1/2))^(1/2)+3/13)*a, y = -(1/4)*a*(4124*((3/13)*3^(1/2)*(15-8*3^(1/2))^(1/2)-(5/26)*3^(1/2)+(11/26)*(15-8*3^(1/2))^(1/2)+3/13)^4*3^(1/2)*(45-24*3^(1/2))^(1/2)+38376*((3/13)*3^(1/2)*(15-8*3^(1/2))^(1/2)-(5/26)*3^(1/2)+(11/26)*(15-8*3^(1/2))^(1/2)+3/13)^4*3^(1/2)+1476*((3/13)*3^(1/2)*(15-8*3^(1/2))^(1/2)-(5/26)*3^(1/2)+(11/26)*(15-8*3^(1/2))^(1/2)+3/13)^3*3^(1/2)*(45-24*3^(1/2))^(1/2)-7145*((3/13)*3^(1/2)*(15-8*3^(1/2))^(1/2)-(5/26)*3^(1/2)+(11/26)*(15-8*3^(1/2))^(1/2)+3/13)^4*(45-24*3^(1/2))^(1/2) ... 9/13)*(15-8*3^(1/2))^(1/2)+70542/13+426*(45-24*3^(1/2))^(1/2)), 0 < a}
Error, (in geometry:-coordinates) wrong type of argument

I am trying to know the coordinates of PP with a and b parameters
p and q are the lengths of the horizontal and vertical sides of the triangle T1
more I know that : xPP=a-2*a*p*(b-q)./(2*a*b-p*q): that  yPP=b-2*b*q*(a-p)./(2*a*b-p*q):and that 2*(a-p)*(b-q)=p*q)
Program attempt :
restart:
  with(geometry):
  _EnvHorizontalName := x:
  _EnvVerticalName := y:
  b := a*(1/2 + 1/6*sqrt(45 - 24*sqrt(3)))^2:
  r := b*sqrt(b)/(sqrt(a + b) + sqrt(a)):
  ellipse(E, x^2/a^2 + y^2/b^2 = 1):
  circle(C3, [point(P3, [(1 + sqrt(3))*r, -r]), r]):
  Equation(C3);

 
ellipse: hint: unable to determine if 1/a^4/(1/2+1/6*(45-24*3^(1/2))^(1/2))^4 is zero
Error, (in geometry:-ellipse) the given polynomial/equation is not an algebraic representation of a ellipse
circle: hint: a^(3/2)*(1/2+1/6*(45-24*3^(1/2))^(1/2))^3/((a+a*(1/2+1/6*(45-24*3^(1/2))^(1/2))^2)^(1/2)+a^(1/2)) > 0
Error, (in geometry:-circle) not enough information: the radius might not be positive
Error, (in geometry:-Equation) wrong type of argument
NULL;
Is it possible to find the coordinates of PP with a and b as parameters. Thank you very much.