// TWIGS3D / Feelers '07.11
// by Tsutomu HIGO  URL:http://www.asahi-net.or.jp/~nj2t-hg/

#version 3.6;
global_settings {
  max_trace_level 5
  assumed_gamma 2.2
  }

light_source {<2000, 10000, -5000> color <1, 1, 1>}
light_source {<-60000, -30000, -40000> color <0, 0, 0>}

camera {
  location <-30, -100, -100>
  angle 4.8
  look_at <0.2, 0.1, 0>
  }

#declare Rndm = seed(29);
#declare R1 = seed(5);
#declare Rdg = 5; //radius of arranging sphere

#macro Dla31 (Nmax, Rds)
  #declare Cp = array[Nmax][6]
  #declare Cp[0][0] = 0; // x of sphere's center
  #declare Cp[0][1] = 0; // z of sphere's center
  #declare Cp[0][2] = 0; // y of sphere's center
  #declare Cp[0][3] = 0; // Distanse of sphere's center
  #declare Cp[0][4] = 2*pi*rand(R1); // Longitude Angle of sphere's center (keido)
  #declare Cp[0][5] = pi/2*(0.5+0.5*rand(R1)); // Latitude Angle of sphere's center (ido)
  sphere {
    <Cp[0][0], Cp[0][2], Cp[0][1]> Rds
    pigment {color rgb <1, 1-0/Nmax*0.8, 0.2>}
    finish {
      reflection 0.3
      phong 2
      brilliance 1
      ambient 0.4
      }
    no_shadow
    }

  #declare Np = 1; #while (Np < Nmax)
    #declare Cp[Np][4] = 2*pi*rand(R1);
    #declare Cp[Np][5] = pi/2*(0.5+0.5*rand(R1));  //pi*rand(R1)-pi/2;
    #if (Np = 1)
          #declare Cp[Np][3] = 2; //Templary size
          #declare Cp[Np][0] = Cp[Np][3]*cos(Cp[Np][5])*cos(Cp[Np][4]);
          #declare Cp[Np][1] = Cp[Np][3]*cos(Cp[Np][5])*sin(Cp[Np][4]);
          #declare Cp[Np][2] = Cp[Np][3]*sin(Cp[Np][5]);
          #declare Nsclose = 0;
    #else
      #declare Ns = Np-1; #while (Ns>0)
          #declare Cp[Np][3] = Cp[Ns][3];
          #declare Cp[Np][0] = Cp[Np][3]*cos(Cp[Np][5])*cos(Cp[Np][4]);
          #declare Cp[Np][1] = Cp[Np][3]*cos(Cp[Np][5])*sin(Cp[Np][4]);
          #declare Cp[Np][2] = Cp[Np][3]*sin(Cp[Np][5]);
          #declare Det = 0;
        #declare Ldis = sqrt(pow(Cp[Np][0]-Cp[Ns][0],2)+pow(Cp[Np][1]-Cp[Ns][1],2)+pow(Cp[Np][2]-Cp[Ns][2],2));
        #if (Ldis < 1) //< 2)
          #declare Apn = asin(Ldis/2/Cp[Ns][3]);
          #declare Cp[Np][3] = Cp[Ns][3]*cos(2*Apn)+sqrt(4-pow(Ldis*cos(Apn),2)); // Real size
          #declare Cp[Np][0] = Cp[Np][3]*cos(Cp[Np][5])*cos(Cp[Np][4]);
          #declare Cp[Np][1] = Cp[Np][3]*cos(Cp[Np][5])*sin(Cp[Np][4]);
          #declare Cp[Np][2] = Cp[Np][3]*sin(Cp[Np][5]);
          #declare Nsclose = Ns;
          #declare Det = 1;
          #declare Ns = 1;
        #end
      #declare Ns = Ns-1; #end
        #if (Det = 0)
          #declare Cp[Np][3] = 2;
          #declare Cp[Np][0] = Cp[Np][3]*cos(Cp[Np][5])*cos(Cp[Np][4]);
          #declare Cp[Np][1] = Cp[Np][3]*cos(Cp[Np][5])*sin(Cp[Np][4]);
          #declare Cp[Np][2] = Cp[Np][3]*sin(Cp[Np][5]);
          #declare Nsclose = 0;
        #end
    #end

  union {
      sphere {<Cp[Np][0], Cp[Np][2], Cp[Np][1]> Rds}
      cylinder {<Cp[Np][0], Cp[Np][2], Cp[Np][1]> <Cp[Nsclose][0], Cp[Nsclose][2], Cp[Nsclose][1]> Rds}
      pigment {color rgb <1,0.7,0.7>}
      finish {
        reflection 0.3
        phong 2
        brilliance 1
        ambient 0.5
        }
      no_shadow
      }

  #declare Np = Np+1; #end

#end

//initialize the coordinates of beads (X=Cp[I][0] Y=Cp[I][1])
#declare Imax=100;
#declare Cd=array[Imax][3]
#declare I=0; #while (I<Imax)
  #declare J=0; #while (J<3)
    #declare Cd[I][J]=0;
  #declare J=J+1; #end
#declare I=I+1; #end

//from cos to angle (0-2*pi)
#macro Acosc (Rc,Dy,Ac)
  #declare Ac=acos(Rc);
  #if (Dy<0) #declare Ac=2*pi-Ac; #end
  #if (Ac=0) #declare Ac=2*pi; #end
#end

//rotate coordinates Nrot at angle Arot
#macro Rota (Nrot,Arot)
  #declare Cdx=Cd[Nrot][0];
  #declare Cdy=Cd[Nrot][1];
  #declare Cd[Nrot][0]=Cdx*cos(Arot)-Cdy*sin(Arot);
  #declare Cd[Nrot][1]=Cdx*sin(Arot)+Cdy*cos(Arot);
#end

//coordinates of starting 7 beads (arranged bead radius = 1)
#declare I=0; #while (I<6)
  #declare Pe=1+I;
  #declare Cd[Pe][0]=Cd[0][0]+2*cos(pi/3*I);
  #declare Cd[Pe][1]=Cd[0][1]+2*sin(pi/3*I);
#declare I=I+1; #end

//initial declarement
#declare Ap=0;
#declare As=0;
#declare Al=0;
#declare Ae=0;
#declare An=0;
#declare Ad=0;

#macro Spculc (N,P,Pe)
  //reduction ratio from plane to arranging shape at point P
  #declare Ao=(sqrt(Cd[P][0]*Cd[P][0]+Cd[P][1]*Cd[P][1]))/2/Rdg;
  #declare Ro=sin(Ao)*cos(Ao)/Ao;
  //rotate to X axis
    Acosc(Cd[N][0]/sqrt(Cd[N][0]*Cd[N][0]+Cd[N][1]*Cd[N][1]),Cd[N][1],An)
    Rota (N,-An)
    Rota (N+1,-An) #declare Cd[N+1][1]=Cd[N+1][1]*Ro;
    Rota (P,-An)   #declare Cd[P][1]=Cd[P][1]*Ro;
  #declare Rs=sqrt(pow(Cd[P][0]-Cd[N][0],2)+pow(Cd[P][1]-Cd[N][1],2))/4;
  Acosc((Cd[P][0]-Cd[N][0])/4/Rs,Cd[P][1]-Cd[N][1],Ap)
  #declare Re=sqrt(pow(Cd[N+1][0]-Cd[N][0],2)+pow(Cd[N+1][1]-Cd[N][1],2))/4;
  Acosc((Cd[N+1][0]-Cd[N][0])/4/Re,Cd[N+1][1]-Cd[N][1],Ae)
  Acosc(Rs,1,As)
  Acosc(Re,1,Al)
  #if (rand(Rndm)>0.3) #declare Ad=0.5*rand(Rndm); #end //confuse the arrangement
  #if (Ae-Ap<0) #declare Ap=Ap-2*pi; #end
  #declare Bo=acos(1/Rdg/tan(pi-2*Ao)); 
  #if (Ao>pi/4) #declare Rb=2*Bo/pi; #else #declare Rb=1; #end
  #if (Ae-Ap-As-Al-Ad>0)
    #declare I=0; #while (I<int((Ae-Ap-As-Al-Ad)*Rb*3/pi)+1)
      #declare Pe=P+1+I;
      #declare Cd[Pe][0]=Cd[N][0]+2*cos(Ap+As+Ad+pi/3*I/Rb);
      #declare Cd[Pe][1]=Cd[N][1]+2*sin(Ap+As+Ad+pi/3*I/Rb)/Ro;
      Rota (Pe,An)
    #declare I=I+1; #end
    #declare Ad=0;
  #else
    #if (Ae-Ap-As-Al>0)
      #declare I=0; #while (I<int((Ae-Ap-As-Al)*Rb*3/pi)+1)
        #declare Pe=P+1+I;
        #declare Cd[Pe][0]=Cd[N][0]+2*cos(Ap+As+pi/3*I/Rb);
        #declare Cd[Pe][1]=Cd[N][1]+2*sin(Ap+As+pi/3*I/Rb)/Ro;
        Rota (Pe,An)
      #declare I=I+1; #end
    #end
  #end
  //rotate to former position
    Rota (N,An)
    #declare Cd[N+1][1]=Cd[N+1][1]/Ro; Rota (N+1,An)
    #declare Cd[P][1]=Cd[P][1]/Ro;     Rota (P,An)
#end

#declare Ps=6;
#declare Jdl=Imax-int(Rdg*pi)*2;
#declare Jd=1; #while (Jd<Jdl)
  Spculc (Jd,Ps,Pe)
  #declare Ps=Pe;
#if (Ao>pi/2-1/Rdg) #declare Jd=Jdl; #end
#declare Jd=Jd+1; #end

#declare I=0; #while (I<Pe-1)
  //reduction ratio from plane to arranging shape at point P
  #if (I>0)
    #declare Ai=sqrt(Cd[I][0]*Cd[I][0]+Cd[I][1]*Cd[I][1])/2/Rdg;
    #declare Ri=sin(Ai)*cos(Ai)/Ai;
  #else
    #declare Ai=0;
    #declare Ri=1;
  #end
  //arranging beads and calculating the coordinates of the holes position
    #declare Cd[I][0]=Ri*Cd[I][0];
    #declare Cd[I][1]=Ri*Cd[I][1];
    #declare Cd[I][2]=Rdg-2*Rdg*sin(Ai)*sin(Ai);
#declare I=I+1; #end
             
//initialize the coordinates for mesh Cp(x,y,z,r)
#declare Vmax=72; //
#declare Hmax=2*Vmax;
#declare Cp=array[Hmax][Vmax+2][4]
#declare V=0; #while (V<Vmax+2)
#declare H=0; #while (H<Hmax)
#declare I=0; #while (I<4)
  #declare Cp[H][V][I]=0;
#declare I=I+1; #end
#declare H=H+1; #end
#declare V=V+1; #end

//triangle of the surface
#macro Trid (Ph,Pv)
  #declare Ph1=mod(Ph-1+Hmax,Hmax);
  triangle {
    <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]>,
    <Cp[Ph1][Pv-1][0],Cp[Ph1][Pv-1][1],Cp[Ph1][Pv-1][2]>,
    <Cp[Ph][Pv-1][0],Cp[Ph][Pv-1][1],Cp[Ph][Pv-1][2]>
    }
  triangle {
    <Cp[Ph1][Pv][0],Cp[Ph1][Pv][1],Cp[Ph1][Pv][2]>,
    <Cp[Ph1][Pv-1][0],Cp[Ph1][Pv-1][1],Cp[Ph1][Pv-1][2]>,
    <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]>
    }
#end

//normal vector of the plar coordinates point
#macro Vd (Ph,Pv)
  <(Cp[mod(Ph-1+Hmax,Hmax)][Pv][1]-Cp[Ph][Pv][1])*(Cp[Ph][Pv-1][2]-Cp[Ph][Pv][2])-(Cp[Ph][Pv-1][1]-Cp[Ph][Pv][1])*(Cp[mod(Ph-1+Hmax,Hmax)][Pv][2]-Cp[Ph][Pv][2]),
  (Cp[Ph][Pv-1][0]-Cp[Ph][Pv][0])*(Cp[mod(Ph-1+Hmax,Hmax)][Pv][2]-Cp[Ph][Pv][2])-(Cp[mod(Ph-1+Hmax,Hmax)][Pv][0]-Cp[Ph][Pv][0])*(Cp[Ph][Pv-1][2]-Cp[Ph][Pv][2]),
  (Cp[mod(Ph-1+Hmax,Hmax)][Pv][0]-Cp[Ph][Pv][0])*(Cp[Ph][Pv-1][1]-Cp[Ph][Pv][1])-(Cp[Ph][Pv-1][0]-Cp[Ph][Pv][0])*(Cp[mod(Ph-1+Hmax,Hmax)][Pv][1]-Cp[Ph][Pv][1])> 
#end

//smooth_triangle of the surface
#macro Trisd (Ph,Pv)
  #declare Ph1=mod(Ph-1+Hmax,Hmax);
    smooth_triangle { 
             <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]> Vd(Ph,Pv)
             <Cp[Ph1][Pv-1][0],Cp[Ph1][Pv-1][1],Cp[Ph1][Pv-1][2]> Vd(Ph1,Pv-1)
             <Cp[Ph][Pv-1][0],Cp[Ph][Pv-1][1],Cp[Ph][Pv-1][2]> Vd(Ph,Pv-1)
             }
    smooth_triangle { 
             <Cp[Ph1][Pv][0],Cp[Ph1][Pv][1],Cp[Ph1][Pv][2]> Vd(Ph1,Pv)
             <Cp[Ph1][Pv-1][0],Cp[Ph1][Pv-1][1],Cp[Ph1][Pv-1][2]> Vd(Ph1,Pv-1)
             <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]> Vd (Ph,Pv)
             } 
#end

//radius caliculation and the translation from polar to xyz coordinsates
#declare V=0; #while (V<Vmax+2)
  #declare H=0; #while (H<Hmax)  
    #if (V=0) #declare Vvmax=0; #else #declare Vvmax=0.999998*(V-1)/Vmax+0.000001; #end
    #declare Cp[H][V][3]=Rdg; //radius
    #declare Cp[H][V][0]=Cp[H][V][3]*sin(pi*Vvmax)*cos(2*pi*H/Hmax);
    #declare Cp[H][V][1]=Cp[H][V][3]*sin(pi*Vvmax)*sin(2*pi*H/Hmax);
    #declare Cp[H][V][2]=-Cp[H][V][3]*cos(pi*Vvmax);
    #declare P=0; #while (P<Pe)
      #declare Cd2=pow(Cp[H][V][0]-Cd[P][0],2)+pow(Cp[H][V][1]-Cd[P][1],2)+pow(Cp[H][V][2]-Cd[P][2],2);
      #declare Cp[H][V][3]=Cp[H][V][3]-1.6*0.1/(Cd2+0.1);
    #declare P=P+1; #end
    #declare Cp[H][V][0]=Cp[H][V][3]*sin(pi*Vvmax)*cos(2*pi*H/Hmax);
    #declare Cp[H][V][1]=Cp[H][V][3]*sin(pi*Vvmax)*sin(2*pi*H/Hmax);
    #declare Cp[H][V][2]=-Cp[H][V][3]*cos(pi*Vvmax);
  #declare H=H+1; #end
#declare V=V+1; #end
       
//main
union {
  mesh {
    #declare V=2; #while (V<Vmax+2)
    #declare H=0; #while (H<Hmax)
      Trisd (H,V)
    #declare H=H+1; #end
    #declare V=V+1; #end
    }      

#declare P = 0; #while (P < Pe-1)
  union {
    Dla31 (100, 0.2)
    scale 0.15
    rotate z*-90
    rotate z*180*atan(Cd[P][1]/sqrt(Cd[P][0]*Cd[P][0]+Cd[P][2]*Cd[P][2]))/pi 
      #if (Cd[P][0] < 0) rotate y*(180-180*atan(Cd[P][2]/Cd[P][0])/pi)
      #else rotate y*-180*atan(Cd[P][2]/Cd[P][0])/pi 
      #end
    translate <0.6*Cd[P][0], 0.6*Cd[P][1], 0.6*Cd[P][2]> 
    } 
  #declare P = P+1; #end
    pigment {color rgbt <1, 0, 0, 0.3>}
    finish {
      reflection 0.7
      phong 2
      brilliance 1
      ambient 0.3
      }
    rotate <30, 200, 60>
    no_shadow
    }         

background {color rgb <0.01, 0, 0.2>}