// VEINS / Parasol '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 {<200, -80, 10> color <1, 1, 1.6>}

camera {
  location <100, -40, -320>
  angle 58
  look_at <5, 140, 0>
  }

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

#declare Jmax=48;
#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) // 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]=Cp[L][3][0]+3; // Radius
    #declare Cp[L][0][1]=Cp[L][3][1]*cos(Cp[L][4][1]); // X
    #declare Cp[L][1][1]=Cp[L][1][0]+800*pow(1/(Nm+6), 2); // 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
union {
#declare I=0; #while (I<60)
  Cpset (I, Jmax)
  Apset (Jmax, I)
  #declare J=0; #while (J<Jmax)
    sphere {<Cp[J][0][1], Cp[J][1][1], Cp[J][2][1]> Rd}
    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}
    #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 }
  no_shadow
  }

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