Procedural OpenSCAD: random trees

I've been playing with OpenSCAD, and it has some interesting limitations.  One of those is that variables are actually constants set at compile time.  This became immediately apparent when I tried to save a value to a variable within a for loop - the variable is not allowed to have it's value changed!  It means that OpenSCAD is well suited for parametric design, where changing some constants will scale the end result appropriately.

As such, I set myself the task of generating random trees to try and find a way around this issue.  I think the solution is probably inefficient, as the dynamic nature of the algorithm is actually stored as constants during compilation, rather than used and thrown away during run time.

Anyway, this code might be of interest.  It will generate a different tree every time you compile.  I have been using the OpenSCAD version 2013.6.  Creating a render from this takes my machine a very long time.  Be careful with those top level constants!

INITIAL_SCALE=4;  // used to multiply the thickness.
DEFAULT_THCK=2;  // used with scale above.

MAX_BRANCHES=12;  // how many possible branches
MAX_BRANCH_LEN=45; // how long they grow


module growTree() {


 // Place a trunk
  cylinder_ep( [0,0,0], [0,10,60], 10,8);
 // fire off the branching routine
 growBranches([0,10,60], INITIAL_SCALE);

/* Warning: this is a recursive function,
 * meaning it calls itself.  It will not
 * call itself when the scaling factor
 * creates cylinders with a thickness of
 * less than 1.  If you set your initial
 * scale high, you will get a lot of 
 * branches, and may choke your 
 * computer.
module growBranches( origin, scale) {

 // We can't make random numbers at runtime
 // so we have to generate them prior.

 // seed the random number generator
 seed = 42;

 // randomise how much to branch
 // use round() to make it a whole number
 num_branches = round( rands( 1,MAX_BRANCHES,1)[0]);
 // we need a length for each branch 
 branch_len = rands( 2,MAX_BRANCH_LEN, num_branches );

 green_shade = rands( 0.2, 0.7, 1)[0];
 // We create a rotation on the z axis
 // to move our x,y coordinates, and a 
 // rotation on the x axis to move the z
 // coordinates (i.e, similar to azimuth)
 x_a = rands( 0,100, num_branches);
 z_a = rands( 0,360, num_branches);

 // This function only runs and calls itself
    // if the default thickness scales to more
 // than or equal to 1
 if( DEFAULT_THCK * scale >= 1 ) {

  for( i=[0:num_branches] ) {
   color( "Brown" )
   cylinder_ep( origin, 
    [ origin[0] + (branch_len[i]*scale) * cos(z_a[i]) - sin(z_a[i]),
      origin[1] + (branch_len[i]*scale) * sin(z_a[i]) + cos(z_a[i]),
      origin[2] + (branch_len[i]*scale) * cos(x_a[i]) 
    , DEFAULT_THCK*scale, 1 );
   // Call this same function again, passing the end point
           // as the origin of the next call.
   // Note, the scale value is decreased by * 0.4
     [ origin[0] + (branch_len[i]*scale) * cos(z_a[i]) - sin(z_a[i]),
       origin[1] + (branch_len[i]*scale) * sin(z_a[i]) + cos(z_a[i]),
         origin[2] + (branch_len[i]*scale) * cos(x_a[i]) ], 
       scale*0.4 );

   if( branch_len[i] > 0 ) {
    // place some green
     [ origin[0] + (branch_len[i]*scale) * cos(z_a[i]) - sin(z_a[i]),
       origin[1] + (branch_len[i]*scale) * sin(z_a[i]) + cos(z_a[i]),
       origin[2] + (branch_len[i]*scale) * cos(x_a[i]) 
     color( [0,green_shade,0])
      %sphere( r=(5*scale));

// Draw a cylinder between two specified 
// vectors p1, p2, with start and end
// radius r1, r2.
// source:
module cylinder_ep(p1, p2, r1, r2) {


 // Create a vector storing the differences in positions
 assign(vector = [p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] ])

 // store the distance
 assign(distance = sqrt( pow(vector[0], 2) + pow(vector[1], 2) + pow(vector[2], 2)))

 // move from the origin by the differences in position vector
 translate(p1 + vector/2)

 //rotation of XoY plane by the Z axis with the angle 
 // of the [p1 p2] line projection with the X axis on 
 //the XoY plane
 rotate([0, 0, atan2(vector[1], vector[0])] ) 
  // rotation of ZoX plane by the y axis with the angle 
  // given by the z coordinate and the sqrt(x^2 + y^2)) 
  // point in the XoY plane
  rotate([0, atan2(sqrt(pow(vector[0], 2)+pow(vector[1], 2)),vector[2]), 0])
  cylinder(h = distance, r1 = r1, r2 = r2,center=true);