// Diffusion Limited Aggregation / Coral Crown Ver.1 '00.10
// By Tsutomu HIGO  E-mail: nj2t-hg@asahi-net.or.jp

#version 3.1;
global_settings {
  assumed_gamma 2.2
  }

light_source {<0, 150, -50> color <1,1,1>}

camera {
  location <0,40,-50>
  angle 80
  look_at <0, 10, 0>
  }

#declare R1 = seed(7);

#macro Dlaring(Nmax, Dia)
  #declare Amax=int(pi*Dia/2*3);
  #declare Ca=array[Amax]
  #declare Na=0; #while (Na<Amax)
    #declare Ca[Na]=0;
  #declare Na=Na+1; #end
  #declare Cp=array[Nmax][3]
  #declare Ap=2*pi*rand(R1);
  #declare Cp[0][0]=Dia/2*cos(Ap);
  #declare Cp[0][1]=Dia/2*sin(Ap);
  #declare Cp[0][2]=0;
  sphere {
    <Cp[0][0],Cp[0][1],Cp[0][2]> 1
    pigment {color rgbt <1,1,1,0.8>}
    finish {
      reflection 1
      phong 2
      brilliance 1
      ambient 0.4
      }
    no_shadow
    }

  #declare Np=1; #while (Np<Nmax)
    #declare Ap=2*pi*rand(R1);
    #declare Pa=int(Ap*Dia/4*3);
    #declare Ho=max(
       max(
         max(Cp[Ca[mod(Pa-4+Amax,Amax)]][2],Cp[Ca[mod(Pa-3+Amax,Amax)]][2]),
         max(Cp[Ca[mod(Pa-2+Amax,Amax)]][2],Cp[Ca[mod(Pa-1+Amax,Amax)]][2])
         ),
       max(
         max(max(Cp[Ca[mod(Pa,Amax)]][2],Cp[Ca[mod(Pa+1,Amax)]][2]),Cp[Ca[mod(Pa+2,Amax)]][2]),
         max(Cp[Ca[mod(Pa+3,Amax)]][2],Cp[Ca[mod(Pa+4,Amax)]][2])
         )
       )+2.1;
    #declare Ntmax=Ho*10;
    #declare Nt=Ntmax; #while (Nt>-1)
      #declare Cp[Np][0]=Dia/2*cos(Ap);
      #declare Cp[Np][1]=Dia/2*sin(Ap);
      #declare Cp[Np][2]=Ho*Nt/Ntmax;
      #declare Ns=-4; #while (Ns<5)
          #declare Nns=Ca[mod(Pa+Ns+Amax,Amax)];
        #if (pow(Cp[Np][0]-Cp[Nns][0],2)+pow(Cp[Np][1]-Cp[Nns][1],2)+pow(Cp[Np][2]-Cp[Nns][2],2)<4)
          #declare Nt=0;
        #end
      #declare Ns=Ns+1; #end
    #declare Nt=Nt-1; #end
    sphere { <Cp[Np][0],Cp[Np][1],Cp[Np][2]> 1
      pigment {color rgbt <1,1-0.9*Np/Nmax,1-Np/Nmax,0.8>}
      finish {
        reflection 1
        phong 2
        brilliance 1
        ambient 0.4
        }
      no_shadow
      rotate x*-90 rotate y*-90
      }
  #declare Ca[mod(Pa,Amax)]=Np;
  #declare Np=Np+1; #end

#end

Dlaring(600,50) //(Nmax, Dia)