//  
//  pov file of snake cube.(03.10.06)
//  
//  written by h.hamanaka hammer@sci.hyogo-u.ac.jp
//  
//  


#declare N=24;             //number of units

#declare r=0.04;           // roundishness of units
#declare pos=<-3,-3,0>;      // the position of the initial end unit
#declare in=<0,0,1>;      // getting in vector at initial end surface
#declare out=<0,1,0>;      // getting out vector from initial end unit to the next unit

#declare ro=array[24] {    //how to twist at each joint 1:cockwise 2:180 degree 3:anti clockwise
                           //leave the first element 0
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0

/* propeller 
0,1,1,2,3,1,0,3,
1,1,1,2,3,1,0,3,
1,1,1,2,3,1,0,3
*/

};

#declare colA=<0,0,1>;       // one color of alternative colors of units.
#declare colB=<1,1,1>;      // the other color of alternative colors of units.


camera{
  location <15,13,-12>*1.3
  look_at <0,0,0>
  angle 30
}


background{ color <1,1,1>}

light_source{
  <30,30,-30>
  color <1,1,1>
}


light_source{  
  <15,30,-30>
  color <1,1,1>
}

object{
  plane{y,-5}
  pigment{color <1,1,1>}
}     


//----------------------------------------------------------------------------------------------------



// normalize orthogonalize in,out vectors
#declare in=vnormalize(in);
#declare out=vnormalize(out-vdot(out,in)*in);

// set the third vector which, with in and out, forms a basis.
#declare axis=vcross(in,out);


// set in,out vectors as the gettin in,out vectors of (0)-th unit
// (0)-th unit will not appear. its dummy for the convenience of induction.

#declare dummy=in;
#declare in=out;
#declare out=dummy;








// declare the object of one unit.

#declare unit=object{
union{
  polygon{5,
    <0,r,r>,<0,1-r,r>,<0,1-r,1-r/tan(pi/8)>,<0,r,1-r/tan(pi/8)>,<0,r,r>
  }
  polygon{5,
    <-r,r,0>,<-(1-r/tan(pi/8)),r,0>,<-(1-r/tan(pi/8)),1-r,0>,<-r,1-r,0>,<-r,r,0>
  }
  polygon{5,
    <-r/(sqrt(2)*tan(pi/8)),  r,  1-r/(sqrt(2)*tan(pi/8))>,
    <-1+r/(sqrt(2)*tan(pi/8)),r,  r/(sqrt(2)*tan(pi/8))>,
    <-1+r/(sqrt(2)*tan(pi/8)),1-r,r/(sqrt(2)*tan(pi/8))>,
    <-r/(sqrt(2)*tan(pi/8)),  1-r,1-r/(sqrt(2)*tan(pi/8))>,
    <-r/(sqrt(2)*tan(pi/8)),  r,1-r/(sqrt(2)*tan(pi/8))>
  }
  
  polygon{4,
    <-r,1,r>,<-1+r/tan(pi/8),1,r>,<-r,1,1-r/tan(pi/8)>,<-r,1,r>
  }
  
  polygon{4,
    <-r,0,r>,<-1+r/tan(pi/8),0,r>,<-r,0,1-r/tan(pi/8)>,<-r,0,r>
  }
  
  cylinder{
    <-r,r,r>,<-r,1-r,r>,r
  }
  cylinder{
    <-r,r,1-r/tan(pi/8)>,
    <-r,1-r,1-r/tan(pi/8)>,
    r
  }
  cylinder{
    <-1+r/tan(pi/8),r,r>,
    <-1+r/tan(pi/8),1-r,r>,
    r
  }
  cylinder{
    <-1+r/tan(pi/8),r,  r>,
    <-r,            r,1-r/tan(pi/8)>
    r
  }
  cylinder{
    <-1+r/tan(pi/8),r,r>,
    <-r,r,r>
    r
  }
  cylinder{
    <-r,r,1-r/tan(pi/8)>,
    <-r,r,r>
    r
  }
  
  cylinder{
    <-1+r/tan(pi/8),1-r,r>,
    <-r,1-r,1-r/tan(pi/8)>
    r
  }
  cylinder{
    <-1+r/tan(pi/8),1-r,r>,
    <-r,1-r,r>
    r
  }
  cylinder{
    <-r,1-r,1-r/tan(pi/8)>,
    <-r,1-r,r>
    r
  }
  sphere{<-r,1-r,1-r/tan(pi/8)>,r}
  sphere{<-r,r,1-r/tan(pi/8)>,r}
  sphere{<-1+r/tan(pi/8),1-r,r>,r}
  sphere{<-1+r/tan(pi/8),r,r>,r}
  sphere{<-r,1-r,r>,r}
  sphere{<-r,r,r>,r}
  
    
  }
  scale 2
  translate<1,-1,0>
  rotate<0,0,-90>
  rotate<-90,0,0>

#declare r11=vdot(axis,x);
#declare r12=vdot(axis,y);
#declare r13=vdot(axis,z);
    
#declare r21=vdot(in,x);
#declare r22=vdot(in,y);
#declare r23=vdot(in,z);
    
#declare r31=vdot(out,x);
#declare r32=vdot(out,y);
#declare r33=vdot(out,z);
matrix<r11,r12,r13,
       r21,r22,r23,
       r31,r32,r33,
       0,0,0>


  
}



//set the snake cube!

object{
  union{
    
#declare i=0;
#while (i<N)


//decompose x,y,z vectors into the linear combinations of in,out,axis
    
#declare rx1=vdot(x,in);
#declare rx2=vdot(x,out);
#declare rx3=vdot(x,axis);
    
#declare ry1=vdot(y,in);
#declare ry2=vdot(y,out);
#declare ry3=vdot(y,axis);
    
#declare rz1=vdot(z,in);
#declare rz2=vdot(z,out);
#declare rz3=vdot(z,axis);

//make new linear combinations newx,newy,newz by replace the coefficient of in and out
    
#declare newx=rx2*in+rx1*out+rx3*axis;
#declare newy=ry2*in+ry1*out+ry3*axis;
#declare newz=rz2*in+rz1*out+rz3*axis;
    
//make the matrix by docomposing newx,newy,newz
//this matrix maps x,y,z into newx,newy,newz respectively.
//in otherword this matrix maps in,out,axis into out,in,axis respectively

#declare r11=vdot(newx,x);
#declare r12=vdot(newx,y);
#declare r13=vdot(newx,z);
    
#declare r21=vdot(newy,x);
#declare r22=vdot(newy,y);
#declare r23=vdot(newy,z);
    
#declare r31=vdot(newz,x);
#declare r32=vdot(newz,y);
#declare r33=vdot(newz,z);
    
//ang is the angle we twist at the joint.
    
#declare ang=90*ro[i];

//make vectors newx,newy,newz which are vectors obtained by
//rotating x,y,z around 'in'vector (-ang) degrees.
#declare newx=vaxis_rotate(x,in,-ang);
#declare newy=vaxis_rotate(y,in,-ang);
#declare newz=vaxis_rotate(z,in,-ang);

//make another matrix which maps x,y,z in to newx,newy,newz.
//i.e., this matrix is the rotating matrix around 'in'vector (-and) degrees.    
#declare p11=vdot(newx,x);
#declare p12=vdot(newx,y);
#declare p13=vdot(newx,z);
    
#declare p21=vdot(newy,x);
#declare p22=vdot(newy,y);
#declare p23=vdot(newy,z);
    
#declare p31=vdot(newz,x);
#declare p32=vdot(newz,y);
#declare p33=vdot(newz,z);
    
    
#if(mod(i,2)=0)
#declare col=colA;
#else
#declare col=colB;
#end
    
    
    
#declare unit=object{
    unit
    matrix<p11,p12,p13,
    p21,p22,p23,
    p31,p32,p33,
    0,0,0>
    
    matrix<r11,r12,r13,
    r21,r22,r23,
    r31,r32,r33,
    0,0,0>
  } 

  object{
    unit
    pigment{color col}
    translate pos
  }



#declare dummy=in;
#declare in=out;
#declare out=vaxis_rotate(dummy,in,ang);
#declare axis=vcross(out,in);
#declare pos=pos+in+out;    

#declare i=i+1;          
#end
    
  }
}