// hyper-24-cell

#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; // pour une animation en 400 images, cc=clock/400

camera {//orthographic
#declare DistanceObservateur = 4;
#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*pi), HauteurObservateur, DistanceObservateur*sin(Angle/180*pi)> 
look_at  <0, 0,  0>
}

// 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.01; 
#declare GrandRayon = 0.08;

/* Définition 1
   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,1,1,1};
#declare E1 = array[4] {1/sqrt(2),-1/sqrt(2),0,0};  
#declare E2 = array[4] {1/sqrt(6),1/sqrt(6),-2/sqrt(6),0};  
#declare E3 = array[4] {1/sqrt(12),1/sqrt(12),1/sqrt(12),-3/sqrt(12)};
#declare V = array[4] {1,1,1,1}; 

// Définition 2. La mettre en commentaire pour utiliser la définition 1.
/*#declare N = array[4] {1,2,-1,3};
#declare E1 = array[4] {2/sqrt(5),-1/sqrt(5),0,0};  
#declare E2 = array[4] {1/sqrt(30),2/sqrt(30),5/sqrt(30),0};
#declare E3 = array[4] {1/sqrt(10),2/sqrt(10),-1/sqrt(10),-2/sqrt(10)};  
#declare V = array[4] {1,1.5,-1.2,2.7}; */

// 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 d'un vecteur M 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);
//#declare QRTQ = Rot4;

// Définition de l'image d'un vecteur-tableau de 4 éléments VM 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 24 sommets de l'hyper-24-cell. Ils ont pour coordonnées quatre nombres dont deux sont nuls et deux valent 1 ou -1.
   Un sommet est désigné par S[i,j,k,l] où i et j varient de 0 à 3 et indiquent les indices des composantes non nulles, 
   et k et l varient de 0 à 1, codant les valeurs possibles -1 et 1 des composantes non nulles d'indice respectivement i et j du sommet.
   On effectue ensuite la rotation QRTQ et la projection P pour avoir un élément de l'espace de dimension 3.
   Exemple : S[0][1][1][1] est l'image de {1,1,0,0}
   */
   
#declare S = array[4][4][2][2];
#declare i=0;
#while (i<3)
  #declare j=i+1;
  #while (j<4)
    #declare k=0;
    #while (k<2)
      #declare l=0;
      #while (l<2)
        #declare Stemp = array[4] {0,0,0,0};
        #declare Stemp[i] = 2*k-1;
        #declare Stemp[j] = 2*l-1;        
        #declare S[i][j][k][l] = P(Rot(Stemp));
        #declare l=l+1;
      #end
      #declare k=k+1;
    #end
    #declare j=j+1;
  #end
  #declare i=i+1;
#end

/* Définition des arêtes. Les 24 sommets se répartissent en trois hyperoctaèdres. Deux quelconques d'entre eux forment un hypercube.
   Les arêtes du 24-cell sont les arêtes des trois hypercubes possibles.
   On définit d'abord une fonction qui dessine les arêtes d'un hypercube. On donne en paramètre un sommet M de l'hypercube et une
   liste T de quatre vecteurs de translation qui permettent de le compléter, ceux-ci étant donnés après rotation et projection */


#macro aretesHypercube(M,T)
  #local uu = array[4] {0,0,0,0};
  #local i=0; // i est l'indice tel que la tranlation à effectuer soit T[i]
  #while (i<4)
    #local j=0; // j<k<l sont les autres indices, variant de 0 à 3, distincts de i
    #if (j=i)
      #local j=i+1;
    #end
    #local k=j+1;
    #if (k=i)
      #local k=i+1;
    #end
    #local l=k+1;
    #if (l=i)
      #local l=i+1;
    #end
    #local uu[j]=0;
    #while (uu[j]<2)
      #local uu[k]=0;
      #while (uu[k]<2)
        #local uu[l]=0;
        #while (uu[l]<2)
          cylinder{M+uu[j]*T[j]+uu[k]*T[k]+uu[l]*T[l], 
                   M+uu[j]*T[j]+uu[k]*T[k]+uu[l]*T[l]+T[i],
                   Rayon texture{Silver_Metal} pigment{color White}}
          #local uu[l]=uu[l]+1;
        #end
        #local uu[k]=uu[k]+1;
      #end
      #local uu[j]=uu[j]+1;
    #end  
    #local i=i+1;
  #end
#end  

// Si on veut sélectionner la couleur des sommets en fonction des trois hyperoctaèdres.
#macro couleur(i,j)
  #if (((i=0) & (j=1)) | ((i=2) & (j=3)))
    color Yellow
  #else
    #if (((i=0) & (j=2)) | ((i=1) & (j=3)))
      color Blue
    #else
      color Red
    #end
  #end
#end

// On affiche les sommets. Décommenter la ligne donnant la couleur souhaitée.
union {        
#declare i=0;
#while (i<3)
  #declare j=i+1;
  #while (j<4)
    #declare k=0;
    #while (k<2)
      #declare l=0;
      #while (l<2)
        sphere{S[i][j][k][l], GrandRayon texture{Gold_Metal} pigment{couleur(i,j)}}
        //sphere{S[i][j][k][l], GrandRayon texture{Gold_Metal} pigment{color Yellow}}
        #declare l=l+1;
      #end
      #declare k=k+1;
    #end
    #declare j=j+1;
  #end
  #declare i=i+1;
#end

#declare T = array[4];
#declare T[0] = P(Rot(array[4] {0,-1,0,-1}));
#declare T[1] = P(Rot(array[4] {0,-1,0,1}));
#declare T[2] = P(Rot(array[4] {-1,0,-1,0}));
#declare T[3] = P(Rot(array[4] {-1,0,1,0}));
aretesHypercube(S[0][1][1][1], T)

#declare T[0] = P(Rot(array[4] {-1,0,0,-1}));
#declare T[1] = P(Rot(array[4] {-1,0,0,1}));
#declare T[2] = P(Rot(array[4] {0,-1,-1,0}));
#declare T[3] = P(Rot(array[4] {0,-1,1,0}));
aretesHypercube(S[0][1][1][1], T)

#declare T[0] = P(Rot(array[4] {0,0,-1,-1}));
#declare T[1] = P(Rot(array[4] {0,0,-1,1}));
#declare T[2] = P(Rot(array[4] {-1,-1,0,0}));
#declare T[3] = P(Rot(array[4] {1,-1,0,0}));
aretesHypercube(S[1][2][1][1], T)

rotate -35.6*x
rotate z*60
}

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}