//Description:
// Realistic growth simulation
// uses diffusion limited aggregation 

//Author: Ben Scheele
//Date: mid Nov. --> 12-20-02

#version 3.5;

#declare rad = 1;
 
global_settings{ 
  #if(rad = true)
    radiosity{
    pretrace_start 0.08
    pretrace_end .01
    count 55

    nearest_count 5
    error_bound 1.7
    recursion_limit 2

    low_error_factor 0.75
    gray_threshold 0.0
    minimum_reuse 0.015
    brightness 1

    adc_bailout 0.01/2
  } 
  
  ambient_light off
   #end
}              
 
#declare Size = 350;   //should be an odd integer, needs to be over about 30 to see effect
#declare nrows = Size;
#declare ncols = Size; 
#declare H = Size/2;


camera {       
  location <.5*nrows,-1.0*ncols,-Size*12>
  up y
  right x 
  angle 5.1  
  look_at <.5*nrows,.5*ncols>
}

light_source {
  0*x                 // light's position (translated below)
  color rgb 2.0       // light's color
  area_light
  <15, 0, 0> <0, 0, 15> // lights spread out across this distance (x * z)
  6, 6                // total number of lights in grid (4x*4z = 16 lights)
  adaptive 1          // 0,1,2,3...
  jitter              // adds random softening of light
  circular            // make the shape of the light circular
  orient              // orient light
  translate <-10,-10,-50>*100  // <x y z> position of light
} 

background { rgb 1 }

plane{ z, .5
  pigment{ rgb .2 }
  finish{ phong .5 reflection .01 diffuse .8 ambient 0 
  metallic .2 specular .1 roughness .1 brilliance .8 }
  normal{ waves 0.1 scale .6*Size/100  turbulence .5 translate Size/2*(x+y) }

}


//declare arrays and variables ----------------
 
#declare grow = array[nrows][ncols]

#declare num_parts = 1;     
        
#declare block_1 = 
  sphere{ 0, .4 pigment{ rgb <.7,.5,.1> } 
        finish{ ambient .4 phong .2 reflection .2 }
  }          
  
#declare block_2 = 
  sphere{ 0, .5
    pigment{ rgb <.2,.8,.3> }
        finish{
          ambient .2 diffuse .8 phong 1.5
          reflection{ .1 .2 metallic .15 fresnel on }
          metallic .3 specular .2 roughness .06 brilliance .8
        }
        normal{ bumps .5 scale .001 }
     
      interior{ ior 2.1 }       
  }          
  

// initialize grid, set seed  --------
   
#declare i = 0;
#while(i < nrows) 
  #declare j = 0;
  #while(j < ncols)           
     #declare grow[i][j] = 0;        
  #declare j = j + 1;
  #end
#declare i = i + 1;
#end 
  
#declare mid_row = div(nrows,2);
#declare mid_col = div(ncols,2);

#if( mod(nrows,2) = 1 & mod(ncols,2) = 1) //odd # of rows or columns
  #declare grow[mid_row][mid_col] = 2;
  #declare grow[mid_row-1][mid_col] = 2;
  #declare grow[mid_row+1][mid_col] = 2;
  #declare grow[mid_row][mid_col-1] = 2;
  #declare grow[mid_row][mid_col+1] = 2;
#else
  #declare grow[mid_row][mid_col] = 2;
  #declare grow[mid_row-1][mid_col] = 2;
  #declare grow[mid_row][mid_col-1] = 2;
  #declare grow[mid_row-1][mid_col-1] = 2;  
#end   
   
         
// growing ----------------------------------------------------             

#declare moved_1 = 1; 
#declare moved_2 = 1;
#declare moved_3 = 1;
#declare moved_4 = 1;
#declare Ct = 0;                    

 
#while( (moved_1 = 1 | moved_2 = 1 | moved_3 = 1 | moved_4 = 1) ) //& num_parts < 5000 )  //limits amount of growth
  #declare S1 = seed(400+Ct+clock*40);
  #declare S2 = seed(450+Ct+clock*40);
  #declare S3 = seed(300+Ct+clock*40);
  #declare moved = 0;
 
  
//choose starting position of new particle
//-----------------   
#declare start = 0;
#while( start = 0)       
  
  #if( rand(S1) < .25 & moved_1 = 1)
     #declare side = 1;  //top
  #else #if(rand(S1) < .5 & moved_2 = 1 )
     #declare side = 2;  //right
  #else #if( rand(S1) < .75 & moved_3 = 1)
     #declare side = 3;  //bottom                                           
  #else #if( rand(S1) < 1 & moved_4 = 1)
     #declare side = 4;  //left                                          
  #end
  #end
  #end  
  #end
  
//-------  
 #if(side = 1 | side = 3)  
  #if( side = 1 )
    #declare h = 0; 
    #declare a = int(rand(S2)*ncols); 
  #end   
  #if( side = 3 ) 
    #declare h = nrows-1; 
    #declare a = int(rand(S2)*ncols); 
  #end 
  #if( a != 0 & a != ncols-1)
    #declare grow[h][a] = 1;
    #declare start = 1;  
  #end  
 #end
//-------
 #if(side = 2 | side = 4)
  #if( side = 2 )              
    #declare a = ncols-1;
    #declare h = int(rand(S2)*nrows);     
  #end 
  #if( side = 4 )  
    #declare a = 0;
    #declare h = int(rand(S2)*nrows);
  #end  
  #if( h != 0 & h != nrows-1 )
    #declare grow[h][a] = 1;
    #declare start = 1; 
  #end 
 #end 
//------
  
#end  
//-----------------    

#declare l = .33;
#declare m = .66;
   
  //check neighborhood and move(4 different ways).  ----> try changing motion! 
  
  #if( side = 1 )
   #declare moved_1 = 0;     
   #declare done = 0;
   #while( (h < nrows-1 & a != 0 & a != ncols-1 ) & done = 0)
     #if( grow[h][a-1]=0 & grow[h][a+1]=0 & grow[h+1][a-1]=0 & grow[h+1][a]=0 & grow[h+1][a+1]=0 )
       #declare grow[h][a] = 0;
       #declare h = h + 1; 
         #if( rand(S3) < l)
            #declare a = a - 1;
         #else #if( rand(S3) < m)  
            #declare a = a;  
         #else 
            #declare a = a + 1;
         #end
         #end 
       #declare moved_1 = 1;  
     #else       
       #declare grow[h][a] = 2;
       #declare num_parts = num_parts + 1;
       #declare done = 1;
     #end       
   #end  
  #end  
//---------------------------------------------       
  #if( side = 3)
   #declare moved_3 = 0;     
   #declare done = 0;
   #while( (h > 0 & a != 0 & a != ncols-1) & done = 0)
     #if( grow[h][a-1]=0 & grow[h][a+1]=0 & grow[h-1][a-1]=0 & grow[h-1][a]=0 & grow[h-1][a+1]=0 )
       #declare grow[h][a] = 0;
       #declare h = h - 1; 
         #if( rand(S3) < l)
            #declare a = a - 1;
         #else #if( rand(S3) < m)  
            #declare a = a;  
         #else 
            #declare a = a + 1;
         #end
         #end 
       #declare moved_3 = 1;  
     #else       
       #declare grow[h][a] = 2;
       #declare num_parts = num_parts + 1;
       #declare done = 1;
     #end       
   #end
  #end    
//--------------------------------------------- 
  #if( side = 2)
   #declare moved_2 = 0;     
   #declare done = 0;
   #while( (a > 0 & h != 0 & h != nrows-1) & done = 0)
     #if( grow[h-1][a]=0 & grow[h+1][a]=0 & grow[h-1][a-1]=0 & grow[h][a-1]=0 & grow[h+1][a-1]=0 )
       #declare grow[h][a] = 0;
       #declare a = a - 1; 
         #if( rand(S3) < l)
            #declare h = h - 1;
         #else #if( rand(S3) < m)  
            #declare h = h;  
         #else 
            #declare h = h + 1;
         #end
         #end 
       #declare moved_2 = 1;  
     #else       
       #declare grow[h][a] = 2;
       #declare num_parts = num_parts + 1;
       #declare done = 1;
     #end       
   #end
  #end   

//--------------------------------------------- 
  #if( side = 4)
   #declare moved_4 = 0;     
   #declare done = 0;
   #while( (a < ncols - 1 & h != 0 & h != nrows-1) & done = 0)
     #if( grow[h-1][a]=0 & grow[h+1][a]=0 & grow[h-1][a+1]=0 & grow[h][a+1]=0 & grow[h+1][a+1]=0 )
       #declare grow[h][a] = 0;
       #declare a = a + 1; 
         #if( rand(S3) < l)
            #declare h = h - 1;
         #else #if( rand(S3) < m)
            #declare h = h;  
         #else 
            #declare h = h + 1;
         #end
         #end 
       #declare moved_4 = 1;  
     #else       
       #declare grow[h][a] = 2;
       #declare num_parts = num_parts + 1;
       #declare done = 1;
     #end       
   #end
  #end   
//--------------------------------------------- 

#declare Ct = Ct + 1; //number of particles sent 
#end 
//done growing -----------------------------------------------------
 

//display objects  -----------------

union{
#declare q = 0;
#while(q < nrows)
  #declare r = 0;
  #while(r < ncols)
   
   #declare p = (Size-1.99*sqrt(pow(q-H,2)+pow(r-H,2)))/Size; //radial scaling, offset due to origin being at the corner   

   #if(grow[q][r] = 2 )//& 2.1*sqrt(pow(q-H,2)+pow(r-H,2)) < Size )//limits to a circle
      object{ block_2 
        
        scale <1.3+.5*p,1.3+.5*p,80*p>//*Size/80>
        translate <q+.5,r+.5,0>
      }
 #if(grow[q][r] = 1)  
      object{ block_1 
        translate <q+.5,r+.5,0>
      }
    #end
    #end
    
  #declare r = r + 1;
  #end
#declare q = q + 1;
#end 
             
translate <-nrows/2,-ncols/2,0>
rotate -90*z                    
translate <nrows/2,ncols/2,0>
}