// VEINS / Minaret '04.08
// 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 {<0, 1000, -500> color <1, 1, 1.6>}

camera {
  location <0, 70, -120>
  angle 110
  look_at <10, 16, 0>
  }

#declare Imax=1000;
#declare Cp=array[Imax][5][2]
#declare Rd=0.5;

#declare Jmax=36;
#declare I=0; #while (I<Imax)
  #declare J=0; #while (J<4)
    #declare Cp[I][J][0]=0; // X, Y, Z, R
    #declare Cp[I][J][1]=0; // X, Y, Z, R
  #declare J=J+1; #end
  #declare Cp[I][4][0]=mod(2*pi/Jmax*I, 2*pi); // Angle
  #declare Cp[I][4][1]=Cp[I][4][0]; // Angle
#declare I=I+1; #end

#macro Cpset (Ro, Jm) // Rearrange of Cp[0-Jm][0-4][0-1]
  #declare L=0; #while (L<Jm)
    #declare Cp[L][4][1]=Cp[L][4][0]-0.01; // Angle
    #if (mod(Ro, 6)=0)
      #if (mod(L+2*Ro, 12)=0)
        #declare Jm=Jm+1;
          #declare M=Jm-1; #while (M>L)
            #declare N=0; #while (N<5)
              #declare Cp[M][N][0]=Cp[mod(M-1+Jm, Jm)][N][0];
            #declare Cp[M][N][1]=Cp[M][N][0];
          #declare N=N+1; #end
        #declare M=M-1; #end
      #declare Cp[L][4][1]=Cp[L][4][0]-2*pi/(Jm+1)*0.81+0.01; // Angle
      #end
    #end 
  #declare L=L+1; #end
#end

#macro Apset (Jm, Nm, Ne) // declare of Cp[0-Jm][0-4][0-1]
  #declare L=0; #while (L<Jm)
    #if (Cp[L][4][1]-Cp[mod(L-1+Jm, Jm)][4][1]>0) #declare Cpd=Cp[mod(L-1+Jm, Jm)][4][1]; 
      #else #declare Cpd=Cp[mod(L-1+Jm, Jm)][4][1]-2*pi;
    #end
    #if (Cp[mod(L+1, Jm)][4][1]-Cp[L][4][1]>0) #declare Cpu=Cp[mod(L+1, Jm)][4][1]; 
      #else #declare Cpu=Cp[mod(L+1, Jm)][4][1]+2*pi;
    #end
    #declare Cp[L][4][1]=(2*Cpd+Cp[L][4][1]+2*Cpu)/5;
    #declare Cp[L][3][1]=Ne/2*pow(sin(Nm/Ne*pi), 2); // Radius
    #declare Cp[L][0][1]=Cp[L][3][1]*cos(Cp[L][4][1]); // X
    #declare Cp[L][1][1]=I; // Y
    #declare Cp[L][2][1]=Cp[L][3][1]*sin(Cp[L][4][1]); // Z
  #declare L=L+1; #end
#end

// draw spheres and cylinders
#declare In=70;
#declare Sc=3;
union {
#declare I=1; #while (I<In)
  Cpset (I, Jmax)
  Apset (Jmax, I, In)
  #declare J=0; #while (J<Jmax)
    sphere {<Cp[J][0][1], Cp[J][1][1], Cp[J][2][1]> Rd/Sc}
    cylinder {<Cp[J][0][1], Cp[J][1][1], Cp[J][2][1]>, <Cp[J][0][0], Cp[J][1][0], Cp[J][2][0]>, Rd/Sc}
    #declare K=0; #while (K<5)
      #declare Cp[J][K][0]=Cp[J][K][1];
    #declare K=K+1; #end
  #declare J=J+1; #end
#declare I=I+1; #end
  pigment {color rgb <0.5, 0.5, 0.5>}
  finish {reflection 0.5 ambient 0.7 phong 0.3 brilliance 0.5 diffuse 1.2}
  scale Sc
  no_shadow
  translate <0,-In,0> rotate z*-30
  }

background {color rgb <0, 0, 0.3>}