// Diffusion Limited Aggregation / Collapsing Globe '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,120,-160>
  angle 24
  look_at <0, 0, 0>
  }

#declare R1 = seed(2);

#macro Dlaglobe(Nmax, Dia)
  #declare Nlast=0;
  #declare Amax=2.5/Dia*2;
  #declare Cp=array[Nmax][3]
  #declare Cp[0][0]=0;
  #declare Cp[0][1]=0;
  #declare Cp[0][2]=-Dia/2;

  #declare Np=1; #while (Np<Nmax)
    #declare Ao=Amax+2.1/Dia*2;
    #declare Ap=2*pi*rand(R1);
    #declare Cosap=cos(Ap);
    #declare Sinap=sin(Ap);
    #declare Ntmax=Ao*Dia/2;
    #declare Nt=Ntmax; #while (Nt>0)
      #declare Rp=Dia/2*sin(Ao*Nt/Ntmax);
      #declare Cp[Np][0]=Rp*Cosap;
      #declare Cp[Np][1]=Rp*Sinap;
      #declare Cp[Np][2]=-Dia/2*cos(Ao*Nt/Ntmax);
      #declare Ns=Np-1; #while (Ns>-1)
        #if ((Cp[Np][0]-Cp[Ns][0])*(Cp[Np][0]-Cp[Ns][0])+(Cp[Np][1]-Cp[Ns][1])*(Cp[Np][1]-Cp[Ns][1])+(Cp[Np][2]-Cp[Ns][2])*(Cp[Np][2]-Cp[Ns][2])<8) 
          #declare Ntmaxd=(Nt+1)/Ntmax*Ao*10*Dia/2;
          #declare Ntd=Ntmaxd; #while (Ntd>0)
            #declare Ad=(Nt+1)/Ntmax*Ao*Ntd/Ntmaxd;
            #declare Rpd=Dia/2*sin(Ad);
            #declare Cp[Np][0]=Rpd*Cosap;
            #declare Cp[Np][1]=Rpd*Sinap;
            #declare Cp[Np][2]=-Dia/2*cos(Ad);
            #declare Nsd=Np-1; #while (Nsd>-1)
              #if ((Cp[Np][0]-Cp[Nsd][0])*(Cp[Np][0]-Cp[Nsd][0])+(Cp[Np][1]-Cp[Nsd][1])*(Cp[Np][1]-Cp[Nsd][1])+(Cp[Np][2]-Cp[Nsd][2])*(Cp[Np][2]-Cp[Nsd][2])<4)
                #if (Ao*Nt/Ntmax>pi)
                  #declare Nlast=Np;
                  #declare Np=Nmax-1;
                #end
                #declare Nsclose=Nsd;
                #declare Nsd=0;
                #declare Ntd=1;
              #end
            #declare Nsd=Nsd-1; #end
          #declare Ntd=Ntd-1; #end
          #declare Ns=0;
          #declare Nt=1;
        #end
      #declare Ns=Ns-1; #end
    #declare Nt=Nt-1; #end
    #declare Amax=max(Ad,Amax);
  #declare Np=Np+1; #end
  #if (Nlast=0) #declare Nlast=Nmax; #end

  #declare Np=0; #while (Np<min(Nlast,Nmax))
    sphere { <Cp[Np][0],-Cp[Np][2],-Cp[Np][1]> 1
      pigment {color rgbt <0,1-Np/Nlast,Np/Nlast,0.8>}
      finish {
        reflection 1
        phong 2
        brilliance 1
        ambient 0.4
        }
      no_shadow
      }
  #declare Np=Np+1; #end

#end

Dlaglobe(800,60) //(60,30)//(Nmax, Dia)