// BEADS / Cup '99.12 (v.1.2 '00.01)
// By Tsutomu HIGO  E-mail: nj2t-hg@asahi-net.or.jp

#version 3.1;
global_settings {
  max_trace_level 6
  assumed_gamma 2.2
  }

light_source {<40, -100, -50> color <1,1,1>}

camera {
  location <-10,-80,-100>
  angle 60
  look_at <7, 40, 0>
  }

#declare Rndm = seed(89);
#declare Rdg=46;  //radius of arranging sphere

//initialize the coordinates of beads (X=Cp[I][0] Y=Cp[I][1])
#declare Imax=4200;
#declare Cp=array[Imax][2]
#declare I=0; #while (I<Imax)
  #declare J=0; #while (J<2)
    #declare Cp[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 Cpx=Cp[Nrot][0];
  #declare Cpy=Cp[Nrot][1];
  #declare Cp[Nrot][0]=Cpx*cos(Arot)-Cpy*sin(Arot);
  #declare Cp[Nrot][1]=Cpx*sin(Arot)+Cpy*cos(Arot);
#end

//coordinates of starting 7 beads (arranged bead radius = 1)
#declare I=0; #while (I<6)
  #declare Pe=1+I;
  #declare Cp[Pe][0]=Cp[0][0]+2*cos(pi/3*I);
  #declare Cp[Pe][1]=Cp[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(Cp[P][0]*Cp[P][0]+Cp[P][1]*Cp[P][1]))/2/Rdg;
  #declare Ro=sin(Ao)*cos(Ao)/Ao;
  //rotate to X axis
    Acosc(Cp[N][0]/sqrt(Cp[N][0]*Cp[N][0]+Cp[N][1]*Cp[N][1]),Cp[N][1],An)
    Rota (N,-An)
    Rota (N+1,-An) #declare Cp[N+1][1]=Cp[N+1][1]*Ro;
    Rota (P,-An)   #declare Cp[P][1]=Cp[P][1]*Ro;
  #declare Rs=sqrt(pow(Cp[P][0]-Cp[N][0],2)+pow(Cp[P][1]-Cp[N][1],2))/4;
  Acosc((Cp[P][0]-Cp[N][0])/4/Rs,Cp[P][1]-Cp[N][1],Ap)
  #declare Re=sqrt(pow(Cp[N+1][0]-Cp[N][0],2)+pow(Cp[N+1][1]-Cp[N][1],2))/4;
  Acosc((Cp[N+1][0]-Cp[N][0])/4/Re,Cp[N+1][1]-Cp[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 Cp[Pe][0]=Cp[N][0]+2*cos(Ap+As+Ad+pi/3*I/Rb);
      #declare Cp[Pe][1]=Cp[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 Cp[Pe][0]=Cp[N][0]+2*cos(Ap+As+pi/3*I/Rb);
        #declare Cp[Pe][1]=Cp[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 Cp[N+1][1]=Cp[N+1][1]/Ro; Rota (N+1,An)
    #declare Cp[P][1]=Cp[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(Cp[I][0]*Cp[I][0]+Cp[I][1]*Cp[I][1])/2/Rdg;
    #declare Ri=sin(Ai)*cos(Ai)/Ai;
  #else
    #declare Ai=0;
    #declare Ri=1;
  #end
  //arranging beads
  sphere { <Ri*Cp[I][0],Ri*Cp[I][1],-2*Rdg*sin(Ai)*sin(Ai)> 1
    pigment { color rgbt <1-0.7*I/Pe,0.3+0.7*I/Pe,0,0.99> }
    finish { reflection 1.5 ambient 0.2 phong 0.3 brilliance 0.5 diffuse 2 }
      }
#declare I=I+1; #end

background {color rgbt <0.02,0.02,0>}