%@AUTEUR: Jean-Michel Sarlat
 %@DATE: 14 avril 2007
 
 % === Préparation du fichier postscript embarquant la séquence. 
 input fichierps;
 InitFichierPS(-305,-185,305,185);
 color cwatt; cwatt = (1,0.9,0.6);
 
 % === Arcsin (en degré)
 def asind(expr t) = angle(1+-+t,t) enddef;
 
 % === Pointer les points avec une couleur cerclée de noir.
 def pointe(expr p,c) =
     fill fullcircle scaled 5 shifted p withcolor black;
     fill fullcircle scaled 3 shifted p withcolor c;
 enddef;
 
 % === L'unité.
 u := 1cm;
 
 % === Les paramètres de la courbe.
 a = 5;
 b = 4;
 c = 3;
 am := asind(abs(c-a)/b);
 if am < 0.1: am := 0.1 fi;
 % show am;
 
 
 % === Les fonctions nécessaires à la construction.
 def f(expr t) = ((a*a-c*c)/t + b*b*t) / (2*a*b) enddef;
 def g(expr t) = asind(f(sind(t))) - t enddef;
 def h(expr t) = asind(f(sind(t))) + t enddef;
 def PointB(expr t) = 
     right scaled b rotated (g(t)) shifted (-a,0) 
 enddef;
 def PointC(expr t) = 
     right scaled b rotated (h(t)) shifted (a,0)
 enddef;
 
 % === Les points.
 pair A, B, C, D, M, N;
 A := (-a,0) scaled u;
 D := (a,0)  scaled u;
 
 % === Les courbes.
 path w,cercle;
 cercle := fullcircle scaled (2*b*u);
 
 % === La construction !
 vardef Construction(expr i, t) =
     beginfig(i+1);
 	FondImage(cwatt);
 	CadreImage(0.6cwatt,3);
 	B := PointB(t) scaled u;
 	C := PointC(t) scaled u; 
 	M := 0.5[B,C];
 	if i=0: N := M; w := N; else: w := w--N fi;
 	N := M;
 	draw cercle shifted (-a*u,0)
 	    withcolor 0.8white;
 	draw cercle shifted (a*u,0)
 	    withcolor 0.8white;
 	draw A--B--C--D
 	    withpen pencircle scaled 2
 	    withcolor 0.4white;
 	pointe(A,black); pointe(D,black);
 	pointe(B,white); pointe(C,white);
 	pointe(M,red);
 	drawarrow w 
 	    withpen pencircle scaled 1.5
 	    withcolor red;
     endfig;
 enddef;
 
 
 % === La boucle !
 nb = 100;
 for i=0 upto 2nb:
     if i<= nb:
 	t := 90  - (90 - am) * cosd(180 * i / nb);
     else:
 	t := -90 + (90 - am) * cosd(180 * (i - nb) / nb);
     fi;
     Construction(i,t);
 endfor;
 
 end