PROGRAM Rotations;


{
Epreuve informatique de l'Ecole Polytechnique           93.246
--------------------------------------------------------------


                                                    p 
1/   Soit A une matrice carrée d'ordre n. On note A  le pro-
     duit matriciel de A par lui-même p fois et on considère
     la somme :
                     2    3              p
                A + A + A   +   ....  + A           p
          B =   --------------------------- =  ( õ A )/p
           p               p
     
     Discuter les conditions que A doit satisfaire pour que Bp
     ait une limite finie B non nulle lorsque l'entier p tend
     vers l'infini, et écrire un programme Pascal qui calcule
     cette limite.


2/   Dans un repère orthonormé de dimension 3, écrire la ma-
     trice correspondant à une rotation d'angle Ú autour d'un
     axe passant par l'origine O et défini par un vecteur
     directeur Û.

     Ecrire un programme Pascal qui visualise les positions
     successives d'un tétraèdre ABCD tournant autour d'un axe,
     pour des valeurs croissantes de Ú.



3/   Montrer que l'enchaînement de deux rotations R1 et R2
     successives d'axes Û1 et Û2 et d'angles Ú1 et Ú2 est une
     rotation R3. Montrer que la matrice associée à R3 satis-
     fait à la condition de la première question. Ecrire un
     programme Pascal qui calcule l'axe et l'angle de la rota-
     tion R3.


     Application numérique :
          Pour R1 : Û1=(1,2,3)   Ú1= Ò/5
          et   R2 : Û2=(1,3,2)   Ú2= 1 radian,

          Calculer les rotations : R3=R2°R1 et R4=R1°R2


+------------------------------------------------------------+
¦              Imprimer tous les résultats                   ¦
¦     en indiquant chaque fois à quoi ils correspondent      ¦
+------------------------------------------------------------+

                            -=-=-=-
 }



uses
  crt, modubase;



TYPE
    Angle = real;
    Point = record
              x,y,z :real;
              t  : string[1]
            end;

    Vecteur = record
                 x,y,z :real
              end;

    Arete = Record
                p1,p2 : ^point
            end ;

    Tetraedre = record
                      Sommet : Array[1..4] of point;
                end;

    Rotation = record
                     omega : vecteur;
                     theta : angle;
               end;

Const
     reduction = 0.3;
     coefy=0.15;
     coefz=0.32;
     distobs=30;
     origine    :point=(x:0;y:0;z:0);
     observateur:point=(x:distobs;y:distobs*coefy;z:distobs*coefz);

     n=3;
VAR
     Car         : Char ;
     Question    : Char;
     W,Px,Py,Pz  : Point;
     Sommet      : Array[1..20] of Point;
     Xmin, Xmax,Ymin, yv : real;

     
Function ProjX(p:point):real;
begin
     ProjX := (p.y-p.x*coefy)*reduction;
end;

Function ProjY(p:point):real;
begin
     ProjY := (p.z-p.x*coefz)*reduction;
end;
Procedure NewPoint(VAR p:point;x,y,z:real;t:string);
begin
     p.x := x;
     p.y := y;
     p.z := z;
     p.t := t;
end;


Procedure Vectorise(VAR A:Vecteur;P1,P2:Point);
BEGIN
     With A DO
          BEGIN
               x:=P2.x-P1.x;
               y:=P2.y-P1.y;
               z:=P2.z-P1.z;
          END;
END;

Function Module(v:Vecteur):Real;
BEGIN
     With V do
          Module := SQRT(SQR(x)+SQR(y)+SQR(z));

END;

Function Distance(A,B:Point):Real;
VAR
   v:vecteur;
BEGIN
     Vectorise(V,A,B);
     Distance := Module(v);

END;


PROCEDURE Normalise(VAR v:vecteur);
VAR
   Module : Real;
begin
     With v do
          begin
               Module := SQRT(SQR(x)+SQR(y)+SQR(z));
               x := x/Module;
               y := y/Module;
               z := z/Module;
          end;
end;



Procedure TracePoint(p:point);
VAR
   x,y : real;
begin
     x := ProjX(p);
     y := ProjY(p);
     croix(x,y);
     Ecris(p.t);
end ;


Procedure TraceTrait(p1,p2:Point);
BEGIN
      Deplace(ProjX(p1),ProjY(p1));
      Trace(ProjX(p2),ProjY(p2));
END;


Procedure TraceAxe(x,y,z:real;t:string);
VAR
   p : point;
begin
     NewPoint(p,x,y,z,t);
     TraceTrait(W,p);
     Croix(ProjX(p),ProjY(p));
     Ecris(t);
end;




Procedure Vectoriel(V,W:Vecteur;VAR U:Vecteur);
begin
     With U do
          begin
               x := V.y*W.z-V.z*W.y;
               y := V.z*W.x-V.x*W.z;
               z := V.x*W.y-V.y*W.x;
          end;
end;

Function Scalaire(u,v:Vecteur):Real;
BEGIN
     Scalaire := u.x*v.x+u.y*v.y+u.z*v.z;
END;

Function mixte(u,v,w:Vecteur):Real;
var
   uv:vecteur;
BEGIN
      Vectoriel(u,v,uv);
      mixte:=Scalaire(uv,w);
END;
Function determinant(m:matrice):real;
begin
     determinant:=
                  m[1,1]*(m[2,2]*m[3,3]-m[2,3]*m[3,2])
                + m[1,2]*(m[2,3]*m[3,1]-m[2,1]*m[3,3])
                + m[1,3]*(m[2,1]*m[3,2]-m[2,2]*m[3,1]);
end;

Procedure Matrice_nulle(VAR m:matrice);
var
   i,j:integer;
begin
     for i:=1 to matmax do
         for j:=1 to n do
             m[i,j]:=0;
end;

Procedure Matrice_Identite(VAR m:matrice);
var
   i:integer;
begin
     Matrice_nulle(m);
     for i:=1 to n do
         m[i,i]:=1;
end;

Procedure  Add_Matrice(m1,m2:matrice;var m:matrice);
var
   i,j:integer;
begin
     for i:=1 to n do
         for j:=1 to n do
             m[i,j]:=m1[i,j]+m2[i,j];
end;


Procedure Mult_Matrice(m1,m2:matrice;var m:matrice);
var
   i,j,k:integer;
   s:real;
begin
     for i:=1 to n do
         for j:=1 to n do
             begin
                  s:=0;
                  for k:=1 to n do
                       s:=s+m1[i,k]*m2[k,j];
                  m[i,j]:=s;
             end;
end;


Procedure  Mult_scal_Matrice(m1:matrice;k:real;var m:matrice);
var
   i,j:integer;
begin
     for i:=1 to matmax do
         for j:=1 to n do
             m[i,j]:=k*m1[i,j];
end;

Procedure Affiche_rotation(r1:rotation;s:string);
var
   r:rotation;
begin
     r:=r1;
     with r do
      begin
          normalise(omega);
          with omega do
                    WriteLN(s,
                          ' : Omega=(',x:12:6,',',y:12:6,',',z:12:6,
                          ')   Theta=',theta:12:6);
          WriteLN('1+2cos(theta)=',1+2*cos(theta):12:6);
      end;
end;




Procedure Mat_Vec(m:matrice;u:vecteur;Var v:vecteur);
begin
     with u do
       begin
            v.x:=m[1,1]*x+m[1,2]*y+m[1,3]*z;
            v.y:=m[2,1]*x+m[2,2]*y+m[2,3]*z;
            v.z:=m[3,1]*x+m[3,2]*y+m[3,3]*z;
       end;

end;


Procedure Alea_Matrice(VAR m:matrice);
Var
   i,j:integer;
   x:real;
begin
     Randomize;
     for i:=1 to n do
      for j:=1 to n do
          m[i,j]:=2*random-1;
end;

Procedure Affiche_Matrice(m:matrice);
Var
   i,j:integer;
   x:real;
begin
     for i:=1 to n do
     begin                           
      for j:=1 to n do
       begin
           write('  m[',i,',',j,']=',m[i,j]:12:5);
       end;
      WriteLN;
      end;
      WriteLN('Déterminant =', determinant(m):15:8,
              ' Trace=',m[1,1]+m[2,2]+m[3,3]:15:8);

end;

procedure init(titre:string);
begin
     Modegraphique;
     Efface;
     Xmin := -1.8;
     Xmax := +1.7;
     Ymin := -1.4;
     IsoFenetre(Xmin,Xmax,Ymin);
     yv := 1.2;
     Deplace(0,Ymin+0.1);
     Ecris(Titre);


     Couleur(Rouge);
     NewPoint(W,0,0,0,'W');
     TracePoint(W);
     TraceAxe(1.2,0,0,'x');
     TraceAxe(0,1.2,0,'y');
     TraceAxe(0,0,1.2,'z');
     
     Couleur(-Brillant);
end;                     (* init *)

procedure matrice_rotation(r:rotation;VAR m:matrice);
begin
 with r do
  with omega do
   begin
     Normalise(omega);
     m[1,1]:=  cos(theta)+sqr(x)*(1-cos(theta));
     m[2,2]:=  cos(theta)+sqr(y)*(1-cos(theta));
     m[3,3]:=  cos(theta)+sqr(z)*(1-cos(theta));

     m[1,2]:= x*y*(1-cos(theta))+z*sin(theta);
     m[2,1]:= x*y*(1-cos(theta))-z*sin(theta);

     m[1,3]:= x*z*(1-cos(theta))-y*sin(theta);
     m[3,1]:= x*z*(1-cos(theta))+y*sin(theta);

     m[2,3]:= y*z*(1-cos(theta))+x*sin(theta);
     m[3,2]:= y*z*(1-cos(theta))-x*sin(theta);

  end ; {with r }
end;



procedure Analyse_Matrice(m:matrice);
Var
   an,am,mn,moy:matrice;
   n:integer;
   r:rotation;
   u,v,w:vecteur;
begin
     Matrice_nulle(an);
     Matrice_Identite(mn);
     n:=0;
     repeat
           n:=n+1;
           Mult_Matrice(m,mn,mn);
           Add_Matrice(an,mn,an);
           Mult_scal_Matrice(an,1/n,moy);
           if (n mod 40)=0 then
              write('*');

     until (n>200) and (determinant(moy)<0.00001);
     WriteLN;

{     WriteLN('Matrice absorbée=');
     Affiche_matrice(moy); }
     with r.omega do
      begin
           x:=moy[1,1];
           y:=moy[1,2];
           z:=moy[1,3];
           normalise(r.omega);
          { Affiche_rotation(r,'R3'); }
           u.x:=1;
           u.y:=0;
           u.z:=0;
           vectoriel(u,r.omega,u);
           normalise(u);
           Mat_Vec(m,u,v);
           vectoriel(v,u,w);
           r.theta:=asin(module(w));
           r.omega:=w;
           normalise(r.omega);
           Affiche_rotation(r,'Rotation=');

      end;  {with r.omega}


end;


function signe(x:real):integer;
begin
     if x=0 then
        signe:=0
     else
         if x>0 then
            signe:=1
         else
            signe:=-1;
end;

function sens_tetraedre(P1,P2,P3,P4:Point):integer;
var
   v1,v2,v3:vecteur;

begin
     vectorise(v1,P1,P2);
     vectorise(v2,P1,P3);
     vectorise(v3,P1,P4);
     sens_tetraedre:=signe(mixte(v1,v2,v3));
end;


Procedure Dessine_Tetraedre(t:tetraedre);
VAR
   j : integer;
   nbcache :integer;
procedure TraceTrait_si_Visible(i,j:integer);
var
   k,nbcache :integer;
begin
     with t do
      begin
       nbcache:=0;
       for k:=1 to 4 do
        if (k<>i) and (k<>j) then
         if sens_tetraedre(Sommet[i],Sommet[j],Sommet[k],sommet[10-i-j-k])
               <> sens_tetraedre(Sommet[i],Sommet[j],Sommet[k],observateur)
               then
                Inc(nbcache);
       if nbcache<2 then
              TraceTrait(Sommet[i],Sommet[j]);
      end; { with }
end;

begin
         For j:=1 to 3 do
         begin
              TraceTrait_si_Visible(1,1+j);
              TraceTrait_si_Visible(2+ j mod 3, 2+((j+1) mod 3));
         end;
end ;

Procedure Rotation_Tetraedre(r:rotation;t1:tetraedre;VAR t2:tetraedre);
VAR
   j : integer;
   m:matrice;
begin
   Matrice_rotation(r,m);


   For j:=1 to 4 do
       with t1.sommet[j] do
            begin
              t2.sommet[j].x:=m[1,1]*x+m[1,2]*y+m[1,3]*z;
              t2.sommet[j].y:=m[2,1]*x+m[2,2]*y+m[2,3]*z;
              t2.sommet[j].z:=m[3,1]*x+m[3,2]*y+m[3,3]*z;
            end ;
end ;

Procedure Question1;
Var
   an,m,mn,moy:matrice;
   n,i:integer;
   d:real;
   rot:rotation;
Const
   initr:rotation = (omega: (x:1;y:1;z:1);
                     theta: 0.2)             ;

procedure saisie_matrice(var a:matrice);
var i,j:integer;
begin
    { for i:=1 to n do
     for j:=1 to n do
     begin
          write('Elément ',i,',',j,': '); readln(a[i,j]);
     end; }

     a[1,1]:= sqrt(3)/2;   a[1,2]:=-1/2;        a[1,3]:=0          ;
     a[2,1]:= sqrt(2)/4;   a[2,2]:= sqrt(6)/4;  a[2,3]:=sqrt(2)/2  ;
     a[3,1]:=-sqrt(2)/4;   a[3,2]:=-sqrt(6)/4;  a[3,3]:=sqrt(2)/2  ;


end;
begin
{     rot:=initr;
         Randomize;
         rot.omega.x:=2*random-1;
         rot.omega.y:=2*random-1;
         rot.omega.z:=2*random-1;


     Matrice_rotation(rot,m);}
     saisie_matrice(m);

     Affiche_matrice(m);
     PAuse;
     Matrice_nulle(an);
     Matrice_Identite(mn);
     n:=0;
     repeat
           n:=n+1;
           Efface;
 {          Normalise(rot.omega);
           WriteLN('omega=',rot.omega.x:12:8,rot.omega.y:12:8,rot.omega.z:12:8);}
           WriteLN('n=',n);
           Mult_Matrice(m,mn,mn);
           Add_Matrice(an,mn,an);
           Mult_scal_Matrice(an,1/n,moy);
           WriteLN('mn=m^n');
           Affiche_matrice(mn);
           WriteLN('moy=(SIGMA m^n)/n');
           Affiche_matrice(moy);
           PAuse;
           Delay(100);
     until n>1000;
end;

Procedure Question2;
Var
   a,b : matrice;
   s:real;
   r:rotation;
   t1,t2 : tetraedre;
Const
   init1:tetraedre =( Sommet:((x:1;y:3;z:0;t:'A'),
                              (x:0;y:2;z:1;t:'B'),
                              (x:1;y:1;z:-1;t:'C'),
                              (x:4;y:2;z:1;t:'D')));

   initr:rotation = (omega: (x:0.5;y:0.5;z:0.5);
                     theta: 0.02)             ;


begin
     Randomize;
     r := initr;
     with r do
      begin

      end;
     Normalise(r.omega);
     t1:=init1;
     t2:=t1;
     Init('Tétraedre');
     Couleur(-Brillant);
     TraceTrait(Origine,Observateur);
     Dessine_Tetraedre(t1);
     repeat
           Rotation_Tetraedre(r,t1,t2);
           Dessine_Tetraedre(t2);
           Dessine_Tetraedre(t1);         
           Delay(50);
           t1:=t2;
     until KeyPressed;
end;

Procedure Question3;
var
   m:matrice;
begin
     m[1,1]:=sqrt(3)/2;
     m[1,2]:=-1/2;
     m[1,3]:=0;
     m[2,1]:=sqrt(2)/4;
     m[2,2]:=sqrt(6)/4;
     m[2,3]:=sqrt(2)/2;
     m[3,1]:=-sqrt(2)/4;
     m[3,2]:=-sqrt(6)/4;
     m[3,3]:=sqrt(2)/2;
     Affiche_matrice(m);
     Analyse_matrice(m);
     PAuse;
end;

Procedure Question4;
Var
   r,r1,r2,r3:rotation;
   m,m1,m2,m3:matrice;
   {          Pour R1 : Û1=(1,2,3)   Ú1= Ò/5
          et   R2 : Û2=(1,3,2)   Ú2= 1 radian,
   }
Const
   initr1:rotation = (omega: (x:1;y:2;z:3);
                     theta: pi/5)             ;
   initr2:rotation = (omega: (x:1;y:3;z:2);
                     theta: 1)             ;
   initr3:rotation = (omega: (x:0;y:0;z:2);
                     theta: pi/6)             ;


begin
     Matrice_rotation(initr1,m1);
     WriteLN('m1=');
     Affiche_matrice(m1);
     Matrice_rotation(initr2,m2);
     WriteLN('m2=');
     Affiche_matrice(m2);
     Matrice_rotation(initr3,m3);
     WriteLN('m3=');
     Affiche_matrice(m3);

     Affiche_rotation(initr1,'R1:');
     Affiche_rotation(initr2,'R2:');
     Affiche_rotation(initr3,'R3:');
     PAuse;

     Mult_Matrice(m1,m2,m);

     WriteLN('R1°R2');
     Affiche_matrice(m);
     Analyse_matrice(m);

     Mult_Matrice(m2,m1,m);
     WriteLN('R2°R1');
     Affiche_matrice(m);
     Analyse_matrice(m);
     PAuse;

     Mult_Matrice(m1,m2,m);
     Mult_Matrice(m,m3,m);
     WriteLN('R1°R2°R3');
     Affiche_matrice(m);
     Analyse_matrice(m);
     Mult_Matrice(m2,m3,m);
     Mult_Matrice(m,m1,m);
     WriteLN('R2°R3°R1');
     Affiche_matrice(m);
     Analyse_matrice(m);
     Mult_Matrice(m1,m3,m);
     Mult_Matrice(m,m2,m);
     WriteLN('R1°R3°R2');
     Affiche_matrice(m);
     Analyse_matrice(m);
     PAuse;


end;

Procedure Question5;
begin
     ClrScr;

     WriteLN('+------------------------------------------------------------+');
     WriteLN('¦ Question n°5 :  Extension à une dimension > 3              ¦ ');
     WriteLN('¦                                                            ¦');
     WriteLN('¦ Pour que la suite Bp converge, il faut que A admette       ¦');
     WriteLN('¦ au moins une valeur propre égale à 1                       ¦');
     WriteLN('¦ et que toutes les autres aient un module <=1               ¦');
     WriteLN('¦                                                            ¦');
     WriteLN('¦ Alors, Bx est la projection de x sur le sous-espace des    ¦');
     WriteLN('¦ vecteurs propres correspondant à la valeur propre 1.       ¦');
     WriteLN('¦                                                            ¦');
     WriteLN('¦ Si toutes les valeurs propres ont un module égal à 1       ¦');
     WriteLN('¦ elles vont nécessairement par paires de complexes conjugués¦');
     WriteLN('¦ de la forme e(it) et e(-it) où t est un angle de rotation  ¦');
     WriteLN('¦                                                            ¦');
     WriteLN('¦ Il faut alors chercher une décomposition de l''espace       ¦');
     WriteLN('¦ en sous-espaces de dimension 2 dans chacun desquels A      ¦');
     WriteLN('¦ est une rotation et un sous-espace résiduel invariant par A¦');
     WriteLN('+------------------------------------------------------------+');
     PAuse;
end;     


Procedure Presentation;
begin
    ModeTexte;
    Efface;
    WriteLN('Rotations dans l''espace E3             ');
    WriteLN('                                       ');
    WriteLN('                                       ');
    WriteLN(' (1)  Limite de An/n                   ');
    WriteLN(' (2)  Rotation d''un tétraèdre         ');
    WriteLN(' (3)  Recherche de omega et theta      ');
    WriteLN(' (4)  Composition R2°R1 et R1°R2       ');
    WriteLN(' (5)  Extension à une dimension > 3    ');
    WriteLN('                                       ');
end;



begin

  Initgraphique;
  while true do
  begin
     Presentation;
     write('Question choisie ? ');
     readln(question);
     Efface;
     writeln('======== Question n°',question,'===============');
     writeln('');

     case question of
          '1':  Question1;
          '2':  Question2;
          '3' : Question3;
          '4' : Question4;
          '5' : Question5;
     end;
     Pause;
  end;

end.