/*=========================================================
  
  Catena disposta lungo una curva parametrica
      
  Metodi e Modelli Matematici          
  Gian Marco Todesco, novembre 2007
  http://archinet.uniroma3.it/mate/metodi/2007

=========================================================*/

#include "colors.inc"
#include "textures.inc"


//
// camera, sorgente di luce, piano
//

camera {  
  location 0.8*<-2,4,-5>
  right x*image_width/image_height
  look_at <0,0,0>
  angle 40
}

default {
  finish {
    ambient 0.5
    diffuse 0.4
  }
}

light_source { <1,4,-2> color 0.7 }


plane {
  y,-4
  pigment {
    checker White Orange   
    scale 4
  }
} 

//
// Curva (retta sul toro)
//

#macro Curva(T) 
  #local Phi = 4*pi*T;  
  #local Theta = 1.5*Phi;
  #local R = 1+0.2*sin(Theta);
  #local Y = 0.5*cos(Theta);    
  <cos(Phi)*R, Y,sin(Phi)*R>  
#end

/*
#macro Curva(T) 
  #local Phi = 2*pi*T;  
  <cos(Phi)*1, 0,sin(Phi)*1>  
#end
*/

//
// Definisco l'anello della catena
//
#declare Anello = union {
  #local R0 = 0.1; // raggio dei cerchi e semidistanza fra i cilindri 
  #local R1 = 0.05; // semispessore dell'anello
  #local D = 0.1;  // semilunghezza cilindri
  #local SemiCerchio = difference {torus {R0,R1} plane {x,0}} 
  object {SemiCerchio translate D*x}
  object {SemiCerchio rotate 180*z translate -D*x}
  cylinder {<-D,0,R0>,< D,0,R0>, R1}
  cylinder {<-D,0,-R0>,< D,0,-R0>, R1}
  texture {Gold_Metal}
}

// dichiaro una lunghezza un po' inferiore in modo da incatenare gli anelli: vedi sotto
#declare LunghezzaAnello = 0.3;   


//
// Calcolo la lunghezza totale della curva (fra T=0 e T=1)
//
#declare LunghezzaCurva = 0.0;
#declare N = 10000; // approssimo la curva come una spezzata fatta da N segmenti
#declare I = 0;
#while(I<N)
  #local Segmento = Curva(I/N) - Curva((I+1)/N);
  #declare LunghezzaCurva = LunghezzaCurva + vlength(Segmento); 
  #declare I=I+1;
#end
#debug concat("Lunghezza curva = ", str(LunghezzaCurva, 10,5), "\n")


//
// creo la catena
//
union {
  
  // calcolo il numero di anelli. ceil(X) e' il piu' piccolo intero >= X
  #local N = ceil(LunghezzaCurva / LunghezzaAnello); 
  
  // voglio un numero pari di anelli (visto che devono essere messi alternativamente
  // di taglio e di piatto). Se N e' dispari aggiungo 1: mod(N,2) e' il resto della divisione
  // intera fra N e 2 e quindi vale 1 se N e' dispari e 0 se N e' pari
  #local N = N + mod(N,2);
  
  // D e' la distanza reale fra un anello e l'altro 
  // (puo' essere piu' piccola di LunghezzaAnello)
  #local D = LunghezzaCurva/N;
  
  #local T = 0;  
  #local OldP = Curva(0);
  
  // ciclo
  #local I = 0;
  #while(I<N)
    
    // calcolo il punto P = Curva(T) tale che |P - OldP| sia circa uguale a D (poco piu' grande) 
    #local P = Curva(T);
    #while(vlength(P-OldP)<D)
      #local T = T + 0.0005;
      #local P = Curva(T);
    #end
    #local OldP = P;   
    
    // calcolo il triedro E0,E1,E2 sul punto P=Curva(T). E0 e' parallelo alla derivata alla curva    
    #local H = 0.00001;
    #local E0 = vnormalize(Curva(T+H) - Curva(T-H));
    #local E1 = vnormalize(vcross(E0,y));  
    #local E2 = vcross(E0,E1);
    
    // l'anello I-esimo
    object  {
      Anello 
      rotate 90*I*x      
      matrix <
        E0.x,E0.y,E0.z,
        E1.x,E1.y,E1.z,
        E2.x,E2.y,E2.z,
        P.x,P.y,P.z
      >      
    }
    #local I=I+1;
  #end  

}