//=========================================
// MakeTree 1.0b
//=========================================
// Leaf & Tree macro for Pov-ray 3.1 & higher
// (c) 1999 Gilles Tran tran@inapg.inra.fr
//=========================================
// This file contains 4 macros to build trees
// MakeTree : "wrapper" macro for the MakeBranch macro
// MakeBranch : recursive macro that creates a tree structure
// MakeRoot : recursive macro that creates a 2-level root branch
// MakeLeaf : creates different kinds of leaves
// and 2 vector utilities : vCone and mAlign
//-----------------------------------------
// Look at the end of this file for definition of parameters
// See extree*.pov files for examples of leaf & tree definitions 
// See the showtree.pov file for a template of how to use the include files generated by these macros
//-----------------------------------------
// The size of the leafless tree (lower and upper corner of the containing box)
// is displayed after parsing, and a line is added at the end of generated file (if dofile=true)
//-----------------------------------------
// Credits go to :
// Daniel Skarda for providing the first algorithm I used for tree recursivity
// Margus Ramst for ideas taken from his worm.inc & fur.inc macro
// Steve Pigeon for parameter ideas from his arbre.inc macro
// Ken Tyler for adding leaves in the previous version of this macro
// John VanSickle for his tutorials on vector calculations
//-----------------------------------------
// Release history
// 1.0b April 10, 1999
//      - fixes bug (leaves with stalks weren't properly centered at <0,0,0> in the files)
//
// 1.0a April 1, 1999
//      - display of container size both on screen and at the end of the file
//
// 1.0  March 26, 1999
//      - initial release

//=========================================
// macro MakeBranch
//-----------------------------------------
// This is the main tree macro.
// This macro can create trees that look realistic but no particular tree species.
// It cannot create cone-like trees such as pine trees.
// Beware that large recursion levels (>6) can lead to large parsing times.

// Basically, the MakeBranch creates a branch at position [pos], in the direction [vb].
// The branch is made of segments. 
// [nb] other branches grow on this branch and then the macro calls itself.
// It goes on until a certain level of recursivity [level0] is attained.

// - Leaves are turned on by setting the [leafproba] parameter > 0->1
// - Roots are turned on by setting the [rootproba] parameter > 0->1
// - Trunk twigs are turned on by setting the [twigproba] parameter > 0->1
// - Gravity, antigravity or wind are turned on by the [vpush] vector.
// - Gnarledness is turned on by setting the [fgnarl] value > 0
// - The [dotexture] flag (true/false) applies the txtTree texture to the
//   branch and root segments according to their position

// Look at the end of this file or in the extree*.pov files for other parameter definitions.

//-----------------------------------------
// Note on creation of include files
//-----------------------------------------
// Trees can take long to parse so once you created one you like
// you may want to save it as an include file for later use.
// For this, set the [dofile] parameter to true
// and give names to the tree file, the foliage file and the leaf file
// See the end of this file or the extree*.pov files for examples

//-----------------------------------------
// Note on texturing
//-----------------------------------------
// You will need to define two textures, one for the leaf (txtLeaf) and one for the tree bark (txtTree)
// How the texture is applied depends on the choice of the [dotexture] parameter during the tree generation
// If dotexture=true, then the tree include file contains already the call for the texture
//                    and the texture is positionned according to its position in the tree. This allows for
//                    more realistic texturing but can be expensive memorywise.
// If dotexture=false then you need to apply it here.
// The leaf texture is applied on the whole foliage (though you can add it manually in the leaf include file
// if you want to have leaves individually textured).
//-----------------------------------------

#macro MakeBranch(level,pos,vb,dotop) 
#if (level>level0-3) #warning concat("level ",str(level,0,0),"\n") #end
#if (level>0)     
//-----------------------------------------
// The branch is made up of segments
// First we create the branch segments
//-----------------------------------------
#local nseg=max(1,nseg0-level0+level);  // number of branch segments
#local lbseg=lb[level]/nseg; // length of branch segment
#local posb=array[nseg+1] // consecutive branch positions
#local posb[0]=pos; // start of the branch

#local i=0;

#local ay= 360*(0.5-rand(rd));
#while (i<nseg)
        //-----------------------------------------
        // Segment calculation and positionning
        //-----------------------------------------
        #local rbseg1=rb[level]+(rb[level-1]-rb[level])*i/nseg; // start radius of segment
        #local rbseg2=rb[level]+(rb[level-1]-rb[level])*(i+1)/nseg;// end radius of segment
        // calculates the force (wind, gravity) to apply to the segment
        // the force is stronger at the end of the branch
        #if (vlength(vpush)=0) #local vp=vpush; #else #local vp=vnormalize(vpush)*i/nseg; #end

        // then we add the gnarledness
        #local vg=<0.5-rand(rd),0.5-rand(rd),0.5-rand(rd)>;
        #if (vlength(vg)>0) #local vg=vnormalize(vg); #end

        // now we calculate the position ; vb is the original direction of the segment
         #local posb[i+1]=posb[i]+lbseg*vnormalize(vb+vp*fpush+vg*fgnarl)*(1+(0.5-rand(rd))*stdlseg);
        // segments must be kept over the ground
        #if (posb[i+1].y<aboveground) 
                #local posb[i+1]=<posb[i+1].x,aboveground+rand(rd)*i/nseg,posb[i+1].z>; 
        #end
        // segments must be kept under the sky !

        #if (posb[i+1].y>belowsky) 
                #local posb[i+1]=<posb[i+1].x,belowsky-rand(rd)*i/nseg,posb[i+1].z>; 
        #end

        //-----------------------------------------
        // a cone and a sphere : this is where the tree is created
        //-----------------------------------------
        //#warning concat("x",str(posb[i+1].x,0,3)," y ",str(posb[i+1].y,0,3)," z ",str(posb[i+1].z,0,3),">\n")
        #if(dofile = true)
                #if (dotexture = true) // texture is applied to the branch segment

                        #write(filehandle,"union{cone{",posb[i],",",rbseg1,",",posb[i+1],",",rbseg2,"}\n")
                        #write(filehandle,"sphere{",posb[i],",",rbseg1,"}\n")
                        #write(filehandle,"texture{txtTree ")
                        mAlign(posb[i],posb[i+1],true)
                        #write(filehandle,"}}\n")
                #else
                        #write(filehandle,"cone{",posb[i],",",rbseg1,",",posb[i+1],",",rbseg2,"}\n")
                        #write(filehandle,"sphere{",posb[i],",",rbseg1,"}\n")
                
                #end
        #end                                                
                union{
                        cone{posb[i],rbseg1,posb[i+1],rbseg2}
                        sphere{posb[i],rbseg1}
                        #if (dotexture= true)
                                texture{txtTree mAlign(posb[i],posb[i+1],false)}
                        #end
                }   
        //-----------------------------------------
        // Leaf calculation and positionning
        // This could be improved a lot (for now the leaves just hang from the branch)
        //-----------------------------------------
        #if (rand(rdl)<=leafproba & level<=leaflevel)
                #local alz=alz0*(0.5-rand(rdl));
                #local alx=alx0+stdalx*(0.5-rand(rdl));
                #local P1=posb[i];
                #local P2=posb[i+1];
                #local aly=degrees(atan2(P2.x-P1.x,P2.z-P1.z+0.0001))-180;

                #if (dofile = true)
                        #write(filehandle2,"object{Leaf scale ",(1+stdlsize*rand(rdl))," rotate z*",alz," rotate x*",alx," rotate y*",aly, " translate <",P2.x,",",P2.y,",",P2.z,">}\n")
                #end                                                                 
                object{Leaf scale (1+stdlsize*rand(rdl)) rotate z*alz rotate x*alx rotate y*aly translate P2}
        
        #end
        #if (posb[i+1].x>xMax) #declare xMax=posb[i+1].x; #end
        #if (posb[i+1].y>yMax) #declare yMax=posb[i+1].y; #end
        #if (posb[i+1].z>zMax) #declare zMax=posb[i+1].z; #end
        #if (posb[i+1].x<xMin) #declare xMin=posb[i+1].x; #end
        #if (posb[i+1].y<yMin) #declare yMin=posb[i+1].y; #end
        #if (posb[i+1].z<zMin) #declare zMin=posb[i+1].z; #end
        #local i=i+1;

#end            
//-----------------------------------------
// All segments are now created for the branch
//-----------------------------------------
// Then we create the new branches
//-----------------------------------------
// One branch must continue the present one (otherwise we would get a stump)
//-----------------------------------------
#local new_level=level-1;                  
#local ax=ab[level]+ stdax*(0.5-rand(rd));
#local ay=stday*(0.5-rand(rd));
#local new_vb=vCone(posb[nseg-1],posb[nseg],ax,ay)
MakeBranch(new_level,posb[nseg],new_vb,false)
//-----------------------------------------
// The nb-1 other branches grow on the present branch
//-----------------------------------------
#local j=1;      
#while (j<nb)
        #if (rand(rd)<=branchproba)
//                #local i=int(jb*nseg)+j*(nseg-int(jb*nseg))/nb; // spiraling branches effect
                #local i=int(jb*nseg)+rand(rd)*(nseg-int(jb*nseg));
                #local ay=j*360/nb + stday*(0.5-rand(rd));
                #local ax=ab[level]+ stdax*(0.5-rand(rd));
                #local new_vb=vCone(posb[i],posb[i+1],ax,ay)
                #local new_pos=posb[i]+(posb[i+1]-posb[i])*rand(rd);
                MakeBranch(new_level,new_pos,new_vb,false)
        #end
        #local j=j+1;
#end

//-----------------------------------------
// Sometimes it's better to add a vertical branch on the trunk (dotop=true)
//-----------------------------------------
#if (dotop=true & level=level0)
        #local ax=stdax*(0.5-rand(rd));
        #local ay=stday*(0.5-rand(rd));
        #local new_vb=vCone(posb[nseg-1],posb[nseg],ax,ay)
        MakeBranch(new_level,posb[nseg],new_vb,true)
#end
//-----------------------------------------
// Do you want some twigs on the trunk ?
//-----------------------------------------
#if (twigproba>0 & level=level0)
        #local i=0;
        #while (i<nseg)
                #if (rand(rd)<=twigproba)
                        #local ay=360*(0.5-rand(rd));
                        #local ax=ab[level]+ stdax*(0.5-rand(rd));
                        #local new_vb=vCone(posb[i],posb[i+1],ax,ay)   
                        #local new_pos=posb[i]+(posb[i+1]-posb[i])*rand(rd);
                        MakeBranch(1,new_pos,new_vb,false)
                                
                #end
                #local i=i+1;
        #end

#end
//-----------------------------------------
// Roots anyone ?
//-----------------------------------------
#if (rootproba>0 & level=level0)                         
        #local new_vb=vroot;
        #local new_pos=yroot;
        #local i=0;
        #while (i<nroot)
                #if (rand(rd)<=rootproba)
                        #local ay=((i*360/nroot)+stday*(0.5-rand(rd)));
                        #if (dofile=true)
                                #write(filehandle,"union{\n")
                        #end
                        union{MakeRoot(new_level,new_pos,new_vb) rotate ay*y}
                        #if (dofile=true)
                                #write(filehandle,"rotate y*",ay,"}\n")
                        #end
                #end
                #local i=i+1;
        #end
#end                        

//-----------------------------------------
#end
#end
//-----------------------------------------
// end of macro MakeBranch
//=========================================

//=========================================
// macro MakeRoot
//-----------------------------------------
// Creates a two-level root
// It has the size of a level0-1 branch
// and starts at yroot, in the vroot direction
// Each root branch is rotated in the MakeBranch macro
//-----------------------------------------
#macro MakeRoot(level,pos,vb) 
//-----------------------------------------
// The root is made up of segments
// First we create the root segments
//-----------------------------------------
#if (level>level0-2)
#warning concat("root level ",str(level,0,0),"\n")
#local nseg=max(1,nseg0-level0+level);  // number of branch root segments
#local lbseg=lb[level]/nseg; // length of branch root segment
#local posb=array[nseg+1] // consecutive branch root positions
#local posb[0]=pos; // start of the root

#local i=0;

#local ay= 360*(0.5-rand(rd));                                               
#while (i<nseg)
//        #warning concat("i ",str(i,0,0),"\n")
        //-----------------------------------------
        // Segment calculation and positionning
        //-----------------------------------------
        #local rbseg1=rb[level]*(1-i/nseg); // start radius of segment
        #local rbseg2=rb[level]*(1-(i+1)/nseg); // start radius of segment
        // calculates the gravity to apply to the segment
        // the force is stronger at the end of the root
        #local vp=-y*i/nseg; #if (vlength(vp)>0) #local vp=vnormalize(vp); #end

        // then we add the gnarledness
        #local vg=<(0.5-rand(rd))*2,0.5-rand(rd),2*(0.5-rand(rd))>;
        #if (vlength(vg)>0) #local vg=vnormalize(vg); #end

        // now we calculate the position ; vb is the original direction of the segment
        #local posb[i+1]=posb[i]+lbseg*vnormalize(vb+vp*0.2+vg*fgnarl)*(1+(0.5-rand(rd))*stdlseg);
        #if (posb[i+1].y>0 & (i/nseg)>=0.3) #local posb[i+1]=<posb[i+1].x,0,posb[i+1].z>; #end
//        #warning concat("i=",str(i,0,0)," x=",str(posb[i+1].x,0,3)," y=",str(posb[i+1].y,0,3)," z=",str(posb[i+1].z,0,3),"\n")

        //-----------------------------------------
        // a cone and a sphere : this is where the root is created
        //-----------------------------------------

//        #warning concat("rbseg1 ",str(vlength(posb[i]),0,3)," rbseg2 ",str(vlength(posb[i]),0,3),"\n")
        #if (dofile = true)
                #if (dotexture = true) // texture is applied to the root segment
                        #write(filehandle,"union{cone{",posb[i],",",rbseg1,",",posb[i+1],",",rbseg2,"}\n")
                        #write(filehandle,"sphere{",posb[i],",",rbseg1,"}\n")
                        #write(filehandle,"texture{txtTree ")
                        mAlign(posb[i],posb[i+1],true)
                        #write(filehandle,"}}\n")
                #else                   // no texture applied (faster parsing)
                        #write(filehandle,"cone{",posb[i],",",rbseg1,",",posb[i+1],",",rbseg2,"}\n")
                        #write(filehandle,"sphere{",posb[i],",",rbseg1,"}\n")
              
                #end
        #end
        union{
                cone{posb[i],rbseg1,posb[i+1],rbseg2}
                sphere{posb[i],rbseg1}
                #if (dotexture = true)
                        texture{txtTree mAlign(posb[i],posb[i+1],false)}
                #end
                }
        #if (posb[i+1].x>xMax) #declare xMax=posb[i+1].x; #end
        #if (posb[i+1].y>yMax) #declare yMax=posb[i+1].y; #end
        #if (posb[i+1].z>zMax) #declare zMax=posb[i+1].z; #end
        #if (posb[i+1].x<xMin) #declare xMin=posb[i+1].x; #end
        #if (posb[i+1].y<yMin) #declare yMin=posb[i+1].y; #end
        #if (posb[i+1].z<zMin) #declare zMin=posb[i+1].z; #end


        #local i=i+1;

#end    


//-----------------------------------------
// All segments are now created for the root branch
//-----------------------------------------
// The [nb] other roots grow on the present root branch
//-----------------------------------------
#local j=0;      
#local new_level=level-1;
#while (j<nb)
        #if (rand(rd)<=branchproba)
                #local i=int(jb*nseg)+j*(nseg-int(jb*nseg))/nb;
                #local ay=360*(0.5-rand(rd));
                #local ax=ab[level]+ stdax*(0.5-rand(rd));
                #local new_vb=vCone(posb[i],posb[i+1],ax,ay)
                #if (new_vb.y>0) #local new_vb=<new_vb.x,-rand(rd)*new_vb.y,new_vb.z>;#end
                #local new_pos=posb[i]+(posb[i+1]-posb[i])*rand(rd);
                MakeRoot(new_level,new_pos,new_vb)
        #end
        #local j=j+1;

#end

#end        
#end
//-----------------------------------------
// end of macro MakeRoot
//=========================================

//=========================================
// macro MakeTree
//-----------------------------------------
// Only a wrapper for the MakeBranch macro
//-----------------------------------------
#macro MakeTree()
// First initialize the container dimension
#declare xMax=pos0.x;
#declare yMax=pos0.y;
#declare zMax=pos0.z;
#declare xMin=pos0.x;
#declare yMin=pos0.y;
#declare zMin=pos0.z;

// Fill the size and angle values for each recursion level
FillTreeArrays(level0,lb0,qlb,rb0,qrb,ab0,qab)
       #if(dofile=true)
            #warning concat(ftname," tree file creation start\n")
            #fopen filehandle ftname write      // creates the tree file (contains the branches)
            #write(filehandle,"union{\n")
            #if (leafproba>0)
                #warning concat(fvname," foliage file creation starts\n")
                #fopen filehandle2 fvname write    // creates the foliage file (contains the leaf meshes)
                #write(filehandle2,"union{\n")
            #end
       #end

// Create the tree
       union{MakeBranch(level0,pos0,v0,dotop)}
       
       #if(dofile = true)
            #if (leafproba>0)
                #write (filehandle2,"}\n")          // close the foliage file
                #fclose filehandle2
                #warning concat(fvname," foliage file created\n")
            #end
            #write (filehandle,"}\n")           // closes the tree file
            #write (filehandle,"// Tree in box{<",xMin,",",yMin,",",zMin,">,<",xMax,",",yMax,",",zMax,"> pigment{Green}}\n")
            #fclose filehandle
            #warning concat(ftname," tree file created\n")
       #end
       #warning concat("Tree goes from <",str(xMin,0,3),",",str(yMin,0,3),",",str(zMin,0,3),"> to <",str(xMax,0,3),",",str(yMax,0,3),",",str(zMax,0,3),">\n")

#end
//-----------------------------------------
// end of MakeTree macro
//=========================================

//=========================================
// macro FillTreeArrays
//-----------------------------------------
// This utility macro fills the branch length, radius and angle arrays
// for each recursion level
//-----------------------------------------
#macro FillTreeArrays(level,lb0,qlb,rb0,qrb,ab0,qab)
#declare lb=array[level+1]
#declare rb=array[level+1]
#declare ab=array[level+1]
#local i=level;
#declare lb[i]=lb0;
#declare rb[i]=rb0;
#declare ab[i]=ab0;
#warning concat("level=",str(i,0,0)," lb=",str(lb[i],0,3)," rb=",str(rb[i],0,3),"\n")
#local i=level-1;  
#while (i>=0)
        #declare lb[i]=lb[i+1]*qlb;
        #declare rb[i]=rb[i+1]*qrb;
        #declare ab[i]=ab[i+1]*qab;
#warning concat("level=",str(i,0,0)," lb=",str(lb[i],0,3)," rb=",str(rb[i],0,3),"\n")
        
#local i=i-1;
#end
#end                        
//-----------------------------------------
// end of macro FillTreeArrays
//=========================================

//=========================================
// macro MakeLeaf
//-----------------------------------------
// Creates a mesh of a curved leaf in the -z direction and that bends downward
// The leaf has a stalk in the z direction and is folded in two
// The leaf is not smoothed
// Note 1 : the final leaf is positionned at <0,0,0>
// Note 2 : in fact the leaf is built in the z direction and rotated y*180 later
//-----------------------------------------
/*
#declare lsize=0.3;  // leaf size
#declare seg=20;     // nb of leaf segments and stalk segments : increase to smooth
#declare ll=5;       // leaf length
#declare wl=0.6;     // leaf width 
#declare fl=0.2;     // depth of leaf fold
#declare lpow=1;     // modifies the leaf shape : lpow=3 makes heart-shaped leaf
#declare al=100;     // leaf bending angle : the leaf bends forward until this angle
#declare apow=1;     // curve power, how slow the leaf bends
#declare ndents=0;   // dents in the leaf (8 = "oak" leaf). May require a high seg >20
#declare nlobes=5;   // number of lobes (individual leaves)
#declare alobes=160; // angle made by all the lobes
#declare qlobes=0.9; // size of the farest lobe (0.9 = size will be 0.9*leaf length)                               
#declare ls=5;       // stalk length (0=no stalk)
#declare ws=0.2;     // width of stalk
#declare as=30;      // stalk bending angle : the leaf bends backward
#declare dofile=false;          // true creates a file ; false otherwise
#declare ffname="gtleaf.inc";   // file name
*/
//-----------------------------------------
// Leaf examples
//-----------------------------------------
// 5-lobes leaf
// #declare Leaf=object{MakeLeaf(0.3,10,5,0.6,0.2,1,100,1,0,5,160,0.9,3,0.12,10,dofile,ffname)}
// Wide, heart-shaped leaf
// #declare Leaf=object{MakeLeaf(0.3,10,5,2,0.4,3,100,1,0,1,0,0,3,0.12,10,dofile,ffname)} 
// "Oak", dented leaf (ndents=8)
// #declare Leaf=object{MakeLeaf(0.3,40,5,1,0.1,1,100,1,8,1,0,0,3,0.12,10,dofile,ffname)}
// Trilobed, narrow, very dented leaf (ndents=20)
//#declare Leaf=object{MakeLeaf(0.3,40,5,0.3,0.1,1,100,1,20,3,140,0.6,3,0.12,10,dofile,ffname)} 
// Wide, star-shaped, leaf
//#declare Leaf=object{MakeLeaf(0.3,10,4,1,0.4,1,100,0.5,1,5,200,0.5,3,0.12,10,dofile,ffname)} 
// Narrow, long, leaf, very pointy (lpow=3)
//#declare Leaf=object{MakeLeaf(0.3,60,15,1,0.3,3,140,1,0,1,0,0,3,0.12,10,dofile,ffname)}
// Narrow, long, leaf, straight borders (lpow=0.3)
//#declare Leaf=object{MakeLeaf(0.3,60,15,1,0.3,0.3,100,1,0,1,0,0,3,0.12,10,dofile,ffname)}
// Twin short leaf
//#declare Leaf=object{MakeLeaf(0.3,40,5,0.7,0.5,1,100,1,0,2,80,0.6,3,0.1,10,dofile,ffname)} 
// 7 lobes leaf (conic when seen from the side)
//#declare Leaf=object{MakeLeaf(0.3,40,5,1,0.5,2,100,1,0,7,280,0.8,3,0.1,10,dofile,ffname)} 
// 7 lobes leaf (flat when seen from the side)
//#declare Leaf=object{MakeLeaf(0.3,40,5,1,0.5,2,100,1,0,7,280,0.8,10,0.1,-40,dofile,ffname)} 
//-----------------------------------------
#macro MakeLeaf(lsize,seg,ll,wl,fl,lpow,al,apow,ndents,nlobes,alobes,qlobes,ls,ws,as,dofile,ffname)
#warning "create leaf\n"
#if(dofile=true)          
        #warning concat(ffname," individual leaf file creation starts\n")

        #fopen filehandle ffname write  // creates the leaf mesh (individual leaf)
        #write(filehandle,"mesh{\n")
#end
mesh{
//-----------------------------------------
// leaf
//-----------------------------------------
// leaf points
//-----------------------------------------
#local lseg=ll/seg;  // length of leaf segment
#local nI=3; #local nJ=seg+1;
#local nP=nI*nJ; #local P=array[nP] // array of leaf points
#local pl=<0,0,0>;   // start of leaf
#local j=0; 
#while (j<nJ)          
        #local tl=j/seg;
        #if (j>0) #local pl=pl+lseg*vaxis_rotate(z,x,pow(tl,apow)*al);#end
        
                #local P[j*nI]=pl-x*wl*(pow(sin(tl*pi),lpow)*(2+(pow(sin(tl*pi*ndents),2))))/3;
                #local P[j*nI+1]=pl-y*fl*sin(tl*pi);
                #local P[j*nI+2]=pl+x*wl*(pow(sin(tl*pi),lpow)*(2+(pow(sin(tl*pi*ndents),2))))/3;
        #local j=j+1;        
#end   

//-----------------------------------------
// leaf triangles
//-----------------------------------------
#local q=0;
#while (q<(nI*(nJ-1)-1))
#local i=mod(q,nI);#local j=(q-i)/nI;
#if (i <nI-1)        
        #if (nlobes>1)
        
                #local j=0;
                #local a=-alobes/2;
                #while (j<nlobes)    
                        #local sc=qlobes+(1-qlobes)*sin(j*pi/(nlobes-1));
                        triangle{vaxis_rotate(P[q]*sc,z,a),vaxis_rotate(P[q+1]*sc,z,a),vaxis_rotate(P[q+nI+1]*sc,z,a)}
                        triangle{vaxis_rotate(P[q]*sc,z,a),vaxis_rotate(P[q+nI]*sc,z,a),vaxis_rotate(P[q+nI+1]*sc,z,a)}
                        #if (dofile=true)
                                #local tt1=vaxis_rotate(P[q]*sc,z,a);
                                #local tt2=vaxis_rotate(P[q+1]*sc,z,a);
                                #local tt3=vaxis_rotate(P[q+nI+1]*sc,z,a);
                                #local tt4=vaxis_rotate(P[q+nI]*sc,z,a);
                                #write(filehandle,"triangle{<",tt1.x,",",tt1.y,",",tt1.z,">,<",tt2.x,",",tt2.y,",",tt2.z,">,<",tt3.x,",",tt3.y,",",tt3.z,">}\n")
                                #write(filehandle,"triangle{<",tt1.x,",",tt1.y,",",tt1.z,">,<",tt4.x,",",tt4.y,",",tt4.z,">,<",tt3.x,",",tt3.y,",",tt3.z,">}\n")
                        #end
                        #local j=j+1;
                        #local a=a+alobes/(nlobes-1);
                #end
         #else

                triangle{P[q],P[q+1],P[q+nI+1]} triangle{P[q],P[q+nI],P[q+nI+1]}
                #if (dofile=true)
                        #write(filehandle,"triangle{<",P[q].x,",",P[q].y,",",P[q].z,">,<",P[q+1].x,",",P[q+1].y,",",P[q+1].z,">,<",P[q+nI+1].x,",",P[q+nI+1].y,",",P[q+nI+1].z,">}\n")
                        #write(filehandle,"triangle{<",P[q].x,",",P[q].y,",",P[q].z,">,<",P[q+nI].x,",",P[q+nI].y,",",P[q+nI].z,">,<",P[q+nI+1].x,",",P[q+nI+1].y,",",P[q+nI+1].z,">}\n")
                #end

         #end
         
#end
#local q=q+1;
#end        

//-----------------------------------------
// stalk
//-----------------------------------------
#if (ls>0) // only if a stalk length is >0
//-----------------------------------------
// stalk points
//-----------------------------------------
#local lseg=ls/seg;

        #local nI=2; #local nJ=seg+1;
        #local nP=nI*nJ; #local P=array[nP]
        #local pl=<0,0,0>; 
        #local j=0; 
        #while (j<nJ)          
                #local tl=j/seg;
                #if (j>0) #local pl=pl+lseg*vaxis_rotate(-z,x,-pow(tl,apow)*as);#end
                #local P[j*nI]=pl-x*ws;
                #local P[j*nI+1]=pl+x*ws;
                #local j=j+1;        
        #end   

//-----------------------------------------
// stalk triangles
//-----------------------------------------
#local q=0;
#while (q<(nI*(nJ-1)-1))
#local i=mod(q,nI);#local j=(q-i)/nI;
        #if (i <nI-1) 
                triangle{P[q],P[q+1],P[q+nI+1]} triangle{P[q],P[q+nI],P[q+nI+1]} 
                #if (dofile=true)
                        #write(filehandle,"triangle{<",P[q].x,",",P[q].y,",",P[q].z,">,<",P[q+1].x,",",P[q+1].y,",",P[q+1].z,">,<",P[q+nI+1].x,",",P[q+nI+1].y,",",P[q+nI+1].z,">}\n")
                        #write(filehandle,"triangle{<",P[q].x,",",P[q].y,",",P[q].z,">,<",P[q+nI].x,",",P[q+nI].y,",",P[q+nI].z,">,<",P[q+nI+1].x,",",P[q+nI+1].y,",",P[q+nI+1].z,">}\n")
                #end
                
        #end
#local q=q+1;
#end        
translate <0,-P[nP-1].y,-P[nP-1].z> // centers the leaf at origin
#if (dofile=true)
        #write (filehandle,"translate <0,",-P[nP-1].y,",",-P[nP-1].z,">\n")
#end
#end
//-----------------------------------------
// end of stalk
//-----------------------------------------
scale lsize  // resizes the mesh
rotate y*180 // required so that the leaf points at -z

#if(dofile = true)
        #write (filehandle,"scale ",lsize," rotate y*180}\n")
        #fclose filehandle
        #warning concat(ffname," individual leaf file created\n")

#end
}
#warning "end of leaf creation\n"
#end
//-----------------------------------------
// end of MakeLeaf macro
//=========================================


//=======================================
// macro vCone
//-----------------------------------------
// returns a normalized vector rotated on x (ax) and y (ay)
// along the P2-P1 vector
#macro vCone(P1,P2,ax,ay)
#local p = vaxis_rotate(vaxis_rotate(y,x,ax),y,ay);
#local yV1=vnormalize(P2-P1);
#local xV1=vnormalize(vcross(yV1,z));
#local zV1=vcross(xV1,yV1);
#local answer=vnormalize(<p.x*xV1.x + p.y*yV1.x + p.z*zV1.x,p.x*xV1.y + p.y*yV1.y + p.z*zV1.y,p.x*xV1.z + p.y*yV1.z + p.z*zV1.z>);
answer;
#end   
//-----------------------------------------
// end of vCone macro
//=========================================


//=========================================
// macro mAlign
//-----------------------------------------
// returns a matrix operation that aligns an object or texture along P2-P1
// the object is translated to P1
// translate to P2-P1 if you want the object to be on P2
#macro mAlign(P1,P2,dofile)
#local yV1=vnormalize(P2-P1);
#local xV1=vnormalize(vcross(yV1,z));
#local zV1=vcross(xV1,yV1);
        #if(dofile = true)
                #write(filehandle,"matrix <",xV1.x,",",xV1.y,",",xV1.z,",",yV1.x,",",yV1.y,",",yV1.z,",",zV1.x,",",zV1.y,",",zV1.z,",",P1.x,",",P1.y,",",P1.z,"> ")
        #else
                matrix <xV1.x,xV1.y,xV1.z,yV1.x,yV1.y,yV1.z,zV1.x,zV1.y,zV1.z,P1.x,P1.y,P1.z>

        #end
#end   
//-----------------------------------------
// end of mAlign macro
//=========================================

/*
//=========================================
// Tree macro Parameters
//-----------------------------------------
// These parameters must be declared before calling the tree macro
//-----------------------------------------
#declare dofile=false;    // true creates a tree file ; false otherwise
#declare dotexture=true;  // true creates a textured tree (with the texture following the branch); false otherwise
#declare ftname="gttree4.inc" // file name for tree
#declare fvname="gtleaves4.inc" // file name for tree foliage
#declare ffname="gtleaf4.inc" // file name for individual leaf (mesh)
#declare txtTree=texture{txtTree3} // Bark texture

//-----------------------------------------
// Random streams
// one stream for branches and another one for leaves
// so that the leafed tree has the same structure as the leafless one
//-----------------------------------------
#declare rsd=213;      // random seed
#declare rd=seed(rsd);  // random stream for branches
#declare rdl=seed(rsd); // separate random stream for leaves

//-----------------------------------------
// Tree structure parameters        
//-----------------------------------------
// test with low level0 and neseg0 
// High (>=6) recursion levels [level0] gives more complex trees
// High (>=6) segment numbers [nseg0] level gives smoother trees
//-----------------------------------------
#declare level0=4;      // recursion level
#declare nseg0=3;       // initial number of branch segments (decreases of one at each level)
#declare nb=3;          // max number of branches per level
#declare dotop=false;   // if true, generates an extra branch on top of trunk (sometimes necessary to give more verticality)

#declare lb0=20;        // initial branch length
#declare rb0=1;         // initial branch radius
#declare ab0=25;        // initial branch angle (x angle between the trunk and the first branch)
#declare qlb=0.7;       // branch length decrease ratio (1=constant length)
#declare qrb=0.6;       // branch radius decrease ratio (1=constant radius)
#declare qab=1;       // branch angle decrease ratio (1=constant angle)
#declare stdax=10;      // stdev of x angle (angle x = ax+(0.5-rand)*stdax)
#declare stday=10;      // stdev of y angle (angle y = ay+(0.5-rand)*stday)

#declare branchproba=0.9; // probability of branch apparition 
#declare jb=0.4;        // secondary branches start after this ratio of branch length

#declare fgnarl=0.3;    // gnarledness factor - keep it low <0.8
#declare stdlseg=0.5;     // stddev of branch segment length (0...1) (adds more randomness to branch length)

#declare twigproba=0; // probability to have a twig on a trunk segment

#declare v0=<0,1,0>;    // initial direction - change to give an initial orientation
#declare pos0=<0,0,0>;  // initial trunk position (no need to change this one)

//-----------------------------------------
// constraints parameters
//-----------------------------------------
#declare vpush=<0,-0.2,0>;// direction of push (wind, gravity...) <0,-1,0> = gravity ; <0,1,0> = antigravity
#declare fpush=0.7;       // force of push
#declare aboveground=4; // constrains the branches above this level 
#declare belowsky=1000;  // constrains the branches below this level 

//-----------------------------------------
// root parameters
//-----------------------------------------
#declare rootproba=1;   // probability of root 0=no root ; 1=all [nroot] roots
#declare nroot=5;      // number of main roots;
#declare vroot=<1,-0.4,0>; // initial direction of root 
#declare yroot=<0,0.5,0>;   // initial position of root above ground

//-----------------------------------------
// leaf position parameters
//-----------------------------------------
#declare leafproba=1;   // probability of leaf 0=no leaf ; 1=leaf on each segment
#declare leaflevel=4;   // level where the leaves start to appear
#declare alz0=100;       // max z angle for leaf
#declare alx0=-10;      // start x angle for leaf
#declare stdalx=20;     // stddev x angle for leaf
#declare stdlsize=0.1;  // stddev of leaf size 0=constant size; size = leafsize*(1+stdlsize*rand)

//-----------------------------------------
// leaf structure parameters
//-----------------------------------------
#declare txtLeaf=texture{txtLeaf1} // Leaf texture
#declare lsize=0.3;     // leaf size
#declare seg=10;        // nb of leaf segments and stalk segments : increase to smooth
#declare ll=5;          // leaf length
#declare wl=1;          // leaf width 
#declare fl=0.5;        // depth of leaf fold
#declare lpow=1;        // modifies the leaf shape : lpow=3 makes heart-shaped leaf
#declare al=100;        // leaf bending angle : the leaf bends backward until this angle
#declare apow=1;        // curve power, how slow the leaf bends
#declare ndents=0;      // dents in the leaf (8 = "oak" leaf). May require a high seg >20
#declare nlobes=1;      // number of lobes (individual leaves)
#declare alobes=0;    // angle made by all the lobes
#declare qlobes=1;    // size of the farest lobe (0.9 = size will be 0.9*leaf length)                               
#declare ls=3;          // stalk length (0=no stalk)
#declare ws=0.1;        // width of stalk
#declare as=10;         // stalk bending angle : the stalk bends forward
//-----------------------------------------
// end of parameters
//=========================================
*/
//MakeTree()
//background{color rgb <0.8,0.9,1>}
//plane{y,0 pigment{rgb<1,0.7,0.3>}}
