// hypertetraedre

#include "colors.inc"
#include "textures.inc"
#include "shapes.inc"
#include "stones.inc"
#include "woods.inc"                                                
#include "glass.inc"


// ***************************
// Installation de la caméra *
// ***************************
// On installe la caméra et on définit son orientation.
// L'axe des x est vers la droite, l'axe des y vers le haut, l'axe des z vers le fond.

#declare cc=0.1; // pour une animation en 400 images, cc=clock/400

camera {
#declare DistanceObservateur = 9;
#declare HauteurObservateur = 0;
#declare Angle = 2;  // < 0 pour l'oeil gauche (G.POV) . L'image devra être placée à DROITE en vision croisée
                     // > 0 pour l'oeil droit (D.POV) (en degré), l'image devra être placée à GAUCHE. 
                     // En général 3 à 5° degré d'écart entre les deux yeux. Angle est en degré. 

//Définir l'angle pour choisir l'oeil
location  <DistanceObservateur*cos(Angle/180*3.142), HauteurObservateur, DistanceObservateur*sin(Angle/180*3.142)> 
look_at  <0, 0,  0>
rotate -5*z
rotate 40*y
}

// On définit le fond (par défaut noir)

background { color Black}  


// ********
// Objets *
// ********
// On installe la caméra et on définit son orientation.
// L'axe des x est vers la droite, l'axe des y vers le haut, l'axe des z vers le fond.
                  
#declare Rayon = 0.075; 
#declare GrandRayon = 0.3;

// Définition des cinq sommets de l'hypertétraèdre en dimension 4, centré en 0
#declare ir5 = 1/sqrt(5);
#declare A = array[4] {4*ir5,4*ir5,4*ir5,4*ir5};
#declare B = array[4] {1-ir5,1-ir5,1-ir5,-3-ir5};
#declare C = array[4] {1-ir5,1-ir5,-3-ir5,1-ir5};
#declare D = array[4] {1-ir5,-3-ir5,1-ir5,1-ir5};
#declare E = array[4] {-3-ir5,1-ir5,1-ir5,1-ir5};
#declare Z = array[5] {A,B,C,D,E};

/* Définition 
   d'un vecteur N orthogonal à l'hyperplan de projection
   d'une base orthonormée (E1, E2, E3) de cet hyperplan
   d'un vecteur V, direction d'une projection cavalière */
#declare N = array[4] {1,2,-3,0};
#declare E1 = array[4] {2/sqrt(5),-1/sqrt(5),0,0};
#declare E2 = array[4] {3/sqrt(70),6/sqrt(70),5/sqrt(70),0};
#declare E3 = array[4] {0,0,0,1};
#declare V = array[4] {1,2,-3,0};

// Définition du produit scalaire en dimension 4
#macro ps(V1,V2)
   V1[0]*V2[0]+V1[1]*V2[1]+V1[2]*V2[2]+V1[3]*V2[3]
#end  

// Définition de l'opérateur de projection sur l'hyperplan orthogonal à N parallèlement à V
#macro P(M)
   #local psMN = ps(M,N);
   #local psVN = ps(V,N);
   #local pM = array[4] {M[0]-psMN/psVN*V[0], M[1]-psMN/psVN*V[1], M[2]-psMN/psVN*V[2], M[3]-psMN/psVN*V[3]};
   <ps(pM,E1),ps(pM,E2),ps(pM,E3)>
#end

// Definition du produit de deux matrices 4x4 données
#macro prodMat(MatA,MatB)
   #local prodAB = array[4][4] {{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
   #local i=0;
   #while (i<4)
      #local j=0;
      #while (j<4)
         #local k=0;
         #while (k<4)
            #local prodAB[i][j]=prodAB[i][j]+MatA[i][k]*MatB[k][j];
            #local k=k+1;
         #end
         #local j=j+1;
      #end
      #local i=i+1;
   #end
   prodAB
#end

/* Définition 
   d'une matrice orthogonale Q
   de sa transposée TQ (ou inverse) en dimension 4. Celle-ci peut être calculée avec Xcas dans le fichier O4.xws 
   d'une matrice de rotation Rot4 en dimension 4 en fonction de la variable cc, d'angles 2*cc*pi et 4*cc*pi
   du produit Q*Rot4*TQ
   */
#declare Q = array[4][4] {
{-0.737991,-0.532751,-0.174672,-0.3755459},
{0.462882,-0.807860,0.358079,0.069869},
{0.489083,-0.061135,-0.659389,-0.567686},
{0.043668,0.244541,0.637555,-0.729258}};

#declare TQ = array[4][4] {
{-0.737991,0.462882,0.489083,0.043668},
{-0.532751,-0.807860,-0.061135,0.244541},
{-0.174672,0.358079,-0.659389,0.637555},
{-0.3755459,0.069869,-0.567686,-0.729258}};

#declare Rot4 = array[4][4] {
{cos(cc*2*pi),-sin(cc*2*pi),0,0},
{sin(cc*2*pi),cos(cc*2*pi),0,0},
{0,0,cos(cc*4*pi),-sin(cc*4*pi)},
{0,0,sin(cc*4*pi),cos(cc*4*pi)}};

#declare QR = prodMat(Q,Rot4);
#declare QRTQ = prodMat(QR,TQ);

// Définition de l'image d'un vecteur-tableau de 4 éléments M par la rotation QRTQ
#macro Rot(M)
   #local RotV = array[4] {0,0,0,0};
   #local i=0;
   #while (i<4)
     #local j=0;
     #while (j<4)
        #local RotV[i]=RotV[i]+QRTQ[i][j]*M[j];
        #local j=j+1;
     #end
     #local i=i+1;
   #end
   RotV
#end

/* définition des cinq sommets de l'hypertétraèdre, ayant subi la rotation QRTQ, puis la projection P,
    donnés comme éléments d'un espace de dimension 3. Ce sont ces projetés qui sont dessinés */
#declare S = array[5];
#declare i=0;
#while (i<5)
  #declare S[i] = P(Rot(Z[i]));
  #declare i=i+1;
#end  

union {        
#declare i=0;
#while (i<5)
   sphere{S[i], GrandRayon texture{Gold_Metal} pigment{color Yellow}}
   #declare j=i+1;
   #while (j<5)
      cylinder{S[i], S[j], Rayon texture{Silver_Metal} pigment{color White}}
      #declare j=j+1;
   #end
   #declare i=i+1;
#end

/* On peut mettre des faces triangulaires transparentes */
mesh{
#declare i=0;
#while (i<3)
   #declare j=0;
   #while (j<4)
      #declare k=0;
      #while (k<5)
         triangle{S[i], S[j], S[k]}
      #declare k=k+1;
      #end
   #declare j=j+1;
   #end
   #declare i=i+1;
#end

pigment {color Col_Aquamarine_01 transmit 0.98}
finish {irid {0.15 thickness .1 turbulence 1} reflection  0.1}
}
                         

rotate x*20
//rotate -y*90
}

light_source { <0, -3, 0> color White shadowless}
light_source { <0, 0, -2> color White shadowless}
light_source { <3, 0, 0> color White shadowless}
light_source { <3, 0, 0> color Gray50 shadowless}