/*



Spline Tree Include File v1.3 (c) 2000 Andrew Clinton


This file requires any version of Pov-Ray 3.1.  A patched version of
Pov-Ray is no longer required with version 1.2 of this file!

Refer to "readme.txt" for documentation.

Version History
---------------
1.3 - 10/21/00
- Modified leaf image maps (richer, more realistic colours)
- Modified branch texture
- Added roots
- Fixed leaf alignment (default leafs should now connect properly to branches)
- Fixed backward protrusion of blob branches
- Removed branchquality variable (it is now automatically calculated)
- Added ability to introduce noise into leaf rotations
- Added the ability to scale different leaves separately
- Added global transformation option
- Added tendancy variable
- Improved feedback system
- Changed variable declaration system
- Modified some variable names (this will not be backwards compatible!)

---------------
1.2 - 04/02/00
- Added cubic spline functions
- Now works with Pov-Ray official 3.1
- Memory usage improved significantly (about 10MB are used for examples 1 through 6)

---------------
1.1 - 03/26/00
- Removed probability restriction for leaves
- Rearanged code so each object receives a custom transformation
- Fixed bug where repeated warning messages were generated
- Trace speed doubled
- Memory usage and parsing time improved
- Added example scenes
- Changed simplify variable properties
- Removed simplify segments variable (now included in simplify array)
- Removed branchnoise variable
- Removed splinepoints variable
- rdev now behaves more realistically
- Other minor technical improvements

1.0 - 02/15/00
- First version posted
---------------

Questions contact:
Andrew Clinton clintona@ibm.net
http://povplace.addr.com/

*/

// If you are using megapov, uncomment the line below
//#version unofficial MegaPov 0.5;

//////////////////////////////////////////////////////////////////////////
//
// Initialization macros
//
//////////////////////////////////////////////////////////////////////////

// Sets all array elements to initializer
#macro init_array(dim_size,initializer)
	#local valuearray = array[dim_size]
	#local i = 0;
	#while (i < dim_size)
		#local valuearray[i] = initializer;
		#local i = i+1;
	#end
	valuearray
#end

// Sets start and end values in array, then interpolates inner values exponentially
#macro fill_array (num,initialvalue,finalvalue,base)
	
	#local i = 0;
	#local value_array = array[num]
	
	#while (i < num)
		#if (base = 1)
			#local value_array[i] = initialvalue + (finalvalue-initialvalue)*i/(num-1);
		#else
			#local value_array[i] = initialvalue + (finalvalue-initialvalue)*(pow(base,i)-1)/(pow(base,num-1)-1);
		#end
		#local i = i+1;
	#end

	value_array	
	
#end

#macro init_rand(seed1,seed2,seed3) #declare rv = array[3] {seed(seed1),seed(seed2),seed(seed3)} #end

// Variable wrappers for fill_array
#macro init_divisions(start_val,end_val,power) #declare divisions = fill_array(recursionmax,start_val,end_val,power) #end
#macro init_tdev(start_val,end_val,power) #declare tdev = fill_array(recursionmax,start_val,end_val,power) #end
#macro init_rdev(start_val,end_val,power) #declare rdev = fill_array(recursionmax,start_val,end_val,power) #end
#macro init_branchl(start_val,end_val,power) #declare branchl = fill_array(recursionmax+1,start_val,end_val,power) #end
#macro init_branchw(start_val,end_val,power) #declare branchw = fill_array(recursionmax+2,start_val,end_val,power) #end
#macro init_twigprob(start_val,end_val,power) #declare twigprob = fill_array(recursionmax+1,start_val,end_val,power) #end

// Assigns simplify levels
#macro init_simplify(levels,segments)
	#declare simplify = init_array(recursionmax+1,no)
	#declare i = 0;
	#while (i < levels)
		#declare simplify[recursionmax-i] = segments;
		#declare i = i+1;
	#end
#end

#macro init_splinedev (start_val,mid_val)
	#declare splinedev = array[recursionmax+1]
	#declare i = 0;
	#while (i < (recursionmax+1)*0.5)
		#declare splinedev[i] = start_val + (mid_val - start_val)*i/int((recursionmax+1)*0.5);
		#declare splinedev[recursionmax-i] = start_val + (mid_val - start_val)*i/int((recursionmax+1)*0.5);
		#declare i = i+1;
	#end
#end

#macro init_roots(root_recurse)
	#declare rootrecursionmax = root_recurse;
	#declare rootdivisions = fill_array(root_recurse,8,3,0.5)
	#declare rootl = fill_array(root_recurse+1,branchl[0]*0.7,branchl[0]*0.2,0.5)
	#declare rootw = fill_array(root_recurse+2,branchw[0]*0.75,branchw[0]*0.1,0.75)
#end


//////////////////////////////////////////////////////////////////////////
//
// Random Number Macros
//
//////////////////////////////////////////////////////////////////////////

#macro rn (index,absrange)
	(rand(rv[index])*absrange)
#end

#macro rr (index,absrange)
	((rand(rv[index])-0.5)*absrange)
#end

//////////////////////////////////////////////////////////////////////////
//
// Macro point_to
// Returns the rotation to align from one axis to another
//
//////////////////////////////////////////////////////////////////////////

#macro point_to(p)
	#if (p.x=0 & p.y=0 & p.z=0)
		#local rotx = 0;
	#else
		#local rotx = degrees(atan2(sqrt(pow(p.x,2)+pow(p.z,2)),p.y));
	#end
	#if (p.x=0 & p.z=0)
		#local roty = 0;
	#else
		#local roty = degrees(atan2(p.x,p.z));
	#end
	<rotx,roty,0>
#end

//////////////////////////////////////////////////////////////////////////
//
// Spline functions
// create_cubic - returns an array of coefficients for the piecewise
// spline function f(x) = ax^3 + bx^2 + cx + d
// d_cubic - returns a 3d position on the spline
// v_cubic - returns the tangent on the spline
//
//////////////////////////////////////////////////////////////////////////

#macro create_cubic (points)
	#local co = array[dimension_size(points,1)-3][4]
	#local i = 0;
	#while (i < dimension_size(co,1))
		#local co[i][0] = points[i+1];
		#local co[i][1] = (points[i+2]-points[i])*0.5;
		#local co[i][3] = (points[i+3]-points[i+1])*0.5+co[i][1]+2*(points[i+1]-points[i+2]);
		#local co[i][2] = points[i+2]-points[i+1]-co[i][1]-co[i][3];
		#local i = i+1;
	#end
	co
#end

#macro d_cubic (co,pos)
#local sp = (pos/1.0001)*dimension_size(co,1);
co[int(sp)][0] + (sp-int(sp))*(co[int(sp)][1] + (sp-int(sp))*(co[int(sp)][2] + (sp-int(sp))*co[int(sp)][3]))
#end

#macro v_cubic (co,pos)
#local sp = (pos/1.0001)*dimension_size(co,1);
co[int(sp)][1] + (sp-int(sp))*(2*co[int(sp)][2] + (sp-int(sp))*3*co[int(sp)][3])
#end


//////////////////////////////////////////////////////////////////////////
//
// Splinetree
// calls the tree and root generation macros
//
//////////////////////////////////////////////////////////////////////////

#macro splinetree (level, leadslope)
	
	tree (level, leadslope, (<0,0,0>), (<0,0,0>))
	
	#local i = 0;
	#while (i < rootdivisions[0])
		#local xrotate = (40+(rn(1,30)))/(level+1);
		#local yrotate = ((i*360)/rootdivisions[level])+rr(1,(60/rootdivisions[level]));
		root (level+1, leadslope, <xrotate,yrotate,0>, (<0,0,0>))
		#local i = i+1;
	#end
	
#end

//////////////////////////////////////////////////////////////////////////
//
// tree - Creates the tree structure.  There are 3 steps in this procedure:
//  1. Create the branch spline
//  2. Place objects on the spline (blobs or cone/sphere)
//  3. Spawn new branches
//
//////////////////////////////////////////////////////////////////////////

#macro tree (level, leadslope, currentrot, currentpos)
		
	// Give some feedback
	
	#ifdef (feedback)
		#if (level <= feedback)
			
			#local i = 0;
		      #debug "\n"
		      #while (i <= level)
		      	#debug " "
		      	#local i = i+1;
		      #end
		      #debug str(level,0,0)
			
			#switch (level)
			#case (0)
				#debug " - Building tree..."
				#break;
			#case (1)
				#debug " - Creating major branch..."
				#break;
			#case (2)
				#debug " - Creating minor branch..."
				#break;
			#end
		
		#end
	#end
	
	// Generate a spline to be used for the branch
		
	#local projection = vrotate(vrotate(tendancy,<0,-currentrot.y,0>),<-currentrot.x,0,0>);
	#local dir = branchl[level] * vnormalize(<projection.x,1,projection.z>);
	
	#local dev = splinedev[level]*branchl[level];
	#local branchspline = create_cubic (
		array[5] {
			<0, -0.5*branchl[level], tan(radians(leadslope))*0.5*branchl[level]>
			<0, 0, 0>
			<(rr(0,dev))+dir.x*0.3, 0.5*dir.y, (rr(0,dev))+dir.z*0.3>
			<(rr(0,dev))+dir.x, dir.y, (rr(0,dev))+dir.z>
			<(rr(0,dev))+dir.x*1.5, 1.5*dir.y, (rr(0,dev))+dir.z*1.5>
			
		}
	)
	
	// Create the branch using the above spline
	
	
	#if (!simplify[level])
		blob {
			threshold 0.6
			#local splineposition = 0;
			#while (splineposition <= 1)
				
				#local instwidth = branchw[level]-((branchw[level]-branchw[level+1])*splineposition);
				#local local_transform = transform {
					scale <instwidth, min(level+2,5)*instwidth, instwidth>
					rotate point_to(vrotate(vnormalize(v_cubic(branchspline,splineposition)),currentrot))
					translate vrotate(d_cubic(branchspline,splineposition),currentrot)+currentpos
				}
				sphere {
					0,1,1
					transform local_transform
				}
				#local splineposition = splineposition + ((min(level+2,5)/2)*instwidth)/branchl[level];
			#end
			#if (leadslope != 0)
				clipped_by {
					plane { -y,0 rotate currentrot translate currentpos }
				}
			#end
			texture { branchtexture rotate currentrot }
			transform global_transform
		}
	#else
		union {
			#local splineposition = 0;
			#local increment = 1/simplify[level];
			#while (splineposition < 0.999)
				#local instwidth = branchw[level]-((branchw[level]-branchw[level+1])*splineposition);
				#local instwidth2 = branchw[level]-((branchw[level]-branchw[level+1])*(splineposition+increment));
				#local newpos = vrotate(d_cubic(branchspline,splineposition),currentrot)+currentpos;
				#local newpos2 = vrotate(d_cubic(branchspline,splineposition+increment),currentrot)+currentpos;
				sphere {
					newpos, instwidth*0.6
				}
				cone {
					newpos, instwidth*0.6
					newpos2, instwidth2*0.6
				}
				#local splineposition = splineposition + increment;
			#end
			texture { branchtexture rotate currentrot }
			transform global_transform
		}
	#end
	
	// Here we go... Down a recursion level...
	
	#if (level < recursionmax) //If not a terminal branch, split the branch.
		
		// Seamlessly continue the parent branch
		
		#local newpos = vrotate(d_cubic(branchspline,1),currentrot)+currentpos;
		#local newrot = point_to(vrotate(v_cubic(branchspline,1),currentrot));
		tree(level+1,0,newrot,newpos)
		
		// Divide according to the divisions variable
		
		#local division = 1;
		#local staticrand = rn(0,360);
		#while (division < divisions[level])
			#local xrotate = (rdev[level]*0.5)+(rn(0,(rdev[level]*0.5)));
			#local yrotate = ((division*360)/(divisions[level]-1))+rr(0,360/(divisions[level]*6))+staticrand;			
			#local newpos = vrotate((d_cubic(branchspline,1-(tdev[level]*rand(rv[0])))),currentrot)+currentpos;
			#local newrot = point_to(vrotate(vrotate(<0,1,0>,<xrotate,yrotate,0>),currentrot));
			tree(level+1,xrotate,newrot,newpos)
			#local division = division + 1;
		#end
		
		#while(rand(rv[0])<twigprob[level]) //Put some twigs on the branch
			#local xrotate = (rdev[level]*0.5)+(rn(0,(rdev[level]*0.5)));
			#local yrotate = (rn(0,360));
			#local newpos = vrotate((d_cubic(branchspline,1-(rand(rv[0])*0.8))),currentrot)+currentpos;
			#local newrot = point_to(vrotate(vrotate(<0,1,0>,<xrotate,yrotate,0>),currentrot));
			tree(recursionmax,xrotate,newrot,newpos)
		#end
		
	#else // If we are on a terminal branch, create leaves.
		
		#ifdef(leafarray)
			#ifndef(totalprob)
				#local counter = 0;
				#local totalprob = 0;
				#while (counter < dimension_size(leafprob,1))
				#declare totalprob = totalprob + leafprob[counter];
				#local counter = counter + 1;
				#end
			#end
			#local randomnumber = rand(rv[2]);
			#local currentprob = 0;
			#local intcount = 0;
			#while(currentprob < totalprob)
				#if((currentprob<=randomnumber)&(randomnumber<=(currentprob+leafprob[intcount])))
				#local local_transform = transform {
					scale leafscale[intcount]+rr(2,leafscale[intcount]*0.3)
					#if (leafnoise[intcount] = -1)
						rotate y*rr(2,180)
					#else
						rotate <rn(2,leafnoise[intcount]*90), rr(2,leafnoise[intcount]*90), 0>
						rotate currentrot
					#end
					translate vrotate(d_cubic(branchspline,0.98),currentrot)+currentpos
				}
				object {
					leafarray[intcount]
					transform local_transform
					transform global_transform
				}
				#end
				#local currentprob = currentprob + leafprob[intcount];
				#local intcount = intcount + 1;
			#end // while loop
		#end // ifdef leafarray
		
	#end // End #if (down a level)

#end

//////////////////////////////////////////////////////////////////////////
//
// root - Creates underground root structure
//  This is simply the above macro with some modifications
//
//////////////////////////////////////////////////////////////////////////

#macro root (level, leadslope, currentrot, currentpos)
		
	// Generate a spline to be used for the root
	#local actualdeviance = 0.2*rootl[level];
	#local rootspline = create_cubic (
		array[5] {
		<0, -0.5*rootl[level], -0.5*tan(radians(leadslope))*rootl[level]>
		<0, 0, 0>
		<rr(1,actualdeviance), 0.5*rootl[level], rr(1,actualdeviance)>
		<rr(1,actualdeviance), rootl[level], rr(1,actualdeviance)>
		<rr(1,actualdeviance), 1.5*rootl[level], rr(1,actualdeviance)>
		}
	)
	
	// Create the root using the above spline
	blob {
		threshold 0.6
		#local splineposition = 0;
		#while (splineposition <= 1)
			#local instwidth = rootw[level]-((rootw[level]-rootw[level+1])*splineposition);
			#local local_transform = transform {
				scale <instwidth, min(level+2,5)*instwidth, instwidth>
				rotate point_to(vrotate(vnormalize(v_cubic(rootspline,splineposition)),currentrot))
				translate vrotate(d_cubic(rootspline,splineposition),currentrot)+currentpos
			}
			sphere {
				0,1,1
				transform local_transform
				scale <1,-1,1>
			}
			#local splineposition = splineposition + ((min(level+2,5)/1.5)*instwidth)/rootl[level];
		#end
		clipped_by {
			plane { -y,0 rotate currentrot translate currentpos scale <1,-1,1> }
		}
		texture { branchtexture rotate currentrot }
		transform global_transform
	}
	
	// Spawn new roots	
	#if (level < rootrecursionmax)
				
		#local division = 0;
		#local staticrand = rn(1,360);
		#while (division < rootdivisions[level])
			#local xrotate = (30+(rn(1,30)))/(level+1);
			#local yrotate = ((division*360)/rootdivisions[level])+rr(1,(120/rootdivisions[level]))+staticrand;
			#local newpos = vrotate((d_cubic(rootspline,1-rn(1,0.5))),currentrot)+currentpos;
			#local newrot = point_to(vrotate(vrotate(<0,1,0>,<xrotate,yrotate,0>),currentrot));
			root(level+1,xrotate,newrot,newpos)
			#local division = division + 1;
		#end
		
	#end

#end //Macro root
