***********MoSS - Mophogenetic Surface Structures***********
*        Markus Kristian Kangas -  mkkangas@mit.edu        *
*       Design under Supervision of Una-May O'Reilly       *
*              and Associate Prof. Peter Testa             *
*               (c)Emergent Design Group, MIT              *
*        THE FULLY DOCUMENTED VERSION  June 12, 1999       *
************************************************************


/*********************** THE CODE *************************/
//NOTE: this program contains 3 seperate functions:
//1. L-system creation using tilt/pitch/roll frame
//2. L-system creation using x/y/z frame
//3. Fractal Tree Generator

//currently the x/y/z frame generator has been disabled, but all the code
//for its functionality remains intact below.


#include <AlCoordinateSystem.h>
#include <AlUniverse.h>
#include <AlPolysetNode.h>
#include <AlPolyset.h>
#include <AlPolygon.h>
#include <AlPolysetVertex.h>
#include <AlIterator.h>
#include <AlShader.h>

#include <AlTM.h>

#include <AlLiveData.h>
#include <AlFunction.h>
#include <AlFunctionHandle.h>

#include <math.h>

#include <AlVector.h>

#include <string.h>



extern char * makeAltPath(const char *dirName, const char *suffix);


#define INIT_OPTIONBOX_FILE	"tree_init.scm"
#define OPTIONBOX_FILE		"tree.scm"


float threshhold=1.5;
int shared=3;
int links=0;
float color[15000];

int FILTER=0; //toggle hexagonal filter 0/1 
int FILTER2=0; //toggle universal filter
int FILTER3=0; //filter3
int MESH=0; //toggle mesh on or off as final process 0/1
int TREE=0; //toggle fractal tree mode on or off to replace L-systems

#define TOP 10 //take TOP closest vertices to make mesh with... the
               //higher this is, the more smoothing takes place

#define MAX 15000 //maximum number of vertices allowed

#define BLACK 12
#define GREEN 10 //archaic color definition of a shader
#define RED 11 //same here

#define CHILD_MAX 20 //maximum children allowed per L-system (i.e. branches
                     //of L-system at the same generation depth)
int CHILDREN=3; //number of children for tree mode to use
#define RMAX 10 //maximum number of repellers/attractors allowed


//BOUNDING BOX VARS. (x1,y1,z2) is far bottom corner, (x2,y2,z2) is near top corner
float X1=-15;
float Y1=-15;
float Z1=-15;
float X2=15;
float Y2=15;
float Z2=15;


//STARING POINT VARS
float SX[100];
float SY[100];
float SZ[100];
int SNUM=1;
#define RESOURCES 8000 //limits max number of surfaces that can be generated.

//for the CURRENT PROGRAM these are valid:
float ALPHA=29.99; //rotation angle for +/-
float BETA=29.99; //pitch angle for &/^
float GAMMA=29.99; //yaw angle for </>

//for the VARIANT of this program based on pruzniwicz work, these are valid
//representing the rotations about x,y, and z axis rather than pitch, roll,yaw
//code for this functionality still exists but is not used as described later.
#define ANGLE 90 //+ and -
#define ROTATE_ANGLE 90 //& and ^
#define ROLL_ANGLE 90 //< and >

#define BASE 0 //thickness of branch segments i.e. calls to F in the L-system
int GENERATION_LIMIT=13; //maximum number of generations allowed
#define TRUNK_RATIO 0.75 //thickness ratio applied to sucessive branches for
                         //fractal tree mode
#define LENGTH_RATIO 1 //length ratio applied for fractal tree mode.
#define SIDES 1 //sides the branches have. i.e. 4 is a cubic trunk geometry,
                // 3 is a triangular trunk shape etc.
#define SHRINK_FACTOR 1 //same as length_ratio... but used for creating
                        //non-tree fractal images. i.e. the archaic 
                        //VARIANT mode...

#define GRAMMAR_MAX 10 //maximum number of different grammar sets allowed


char* START="X"; //start string or character of L-system generation
                 //to which grammars are applied


char* INPUT="F[+F][-F]F"; //no longer used

int FACET_SIZE=6; //number of edges per surface facet. in hexagonal geometry
                  //this is 6. this number could be easily determined 
                  //by a quick parse of the surface string in the grammar. 

int resources=RESOURCES; //resource counter (counts down each time a surface
                         //is made)
int surface=0; //toggles 0/1 whether a surface is currently being created
               //or not. important for attractor repeller effects


//no longer important...
const float  COS=cos(ANGLE*M_PI/180);
const float  SIN=sin(ANGLE*M_PI/180);

//used for quick creation of segment trunk geometries.
const float  RCOS=cos(2*M_PI/SIDES);
const float  RSIN=sin(2*M_PI/SIDES);


//NOTE: structs are used often in this code rather than classes since I had
//      some problems getting alias to accept the compiled plugin.

//repeller/attractor struct
struct repeller {
  float x[RMAX];
  float y[RMAX];
  float z[RMAX];
  float power[RMAX];
  int size;
};

//struct to store vertex information for quick filtering/post-processing
struct Vert {
  float x,y,z;
  float ox,oy,oz; //the offset caused by attractors or repellers
  int vertex;
  int polygon;
  int generation;
  int fixed;
};

struct Link {
  int v[10];
    int vnum;
};
struct Link link[5000];

int Vnum=0; //keeps track of number of vertices created
int pnum=0; //keeps track of number of polygons created
struct Vert V[MAX]; //array for storing vertex info

//struct for storing an individual grammar
struct grammar {
  int lines;
  char key[GRAMMAR_MAX];
  char value[GRAMMAR_MAX][80];
};

//struct for storing segment information while growing an L-sytem or
//fractal tree
struct segment {
  int children;
  float bx,by,bz,tx,ty,tz,length,bot_width,top_width;
  int last_angle, last_rotate_angle, last_roll_angle;
  segment* child[CHILD_MAX];
};

//the main tree class. adapted from Alias binary tree plugin.
class Tree
{
public:

  //archaic variables used to store input for the plugin...
	int seed;
	int max_branch;
	double bend;
	double fork_angle;
	double branch_frequency;
	double gravity;
	double branch_length;
	double trunk_taper;
	double min_branch_size;
	double twist;

public:
	Tree();

	statusCode plant( AlPolyset *polyset ) const;
	// routine to fetch the parameters
	statusCode getParams();

};

static inline double DegreesToRads( double degrees )
{
	return degrees * M_PI / 180.0;
}

Tree::Tree()
{
  //dont worry about these vars, no longer important...
	seed			= 0;
	max_branch		= 15;	// this should be stopped by the taper ..
							// but just in case ...
	bend			= 0.2;
	fork_angle	= DegreesToRads( 30 );
	branch_frequency	= 0.5;
	gravity			= -0.1;
	branch_length	= 5;
	trunk_taper		= 0.95;
	min_branch_size	= 0.025;
	twist 			= DegreesToRads( 45 );
}

statusCode Tree::getParams()
{
	if( sSuccess == AlGetDouble( "gravity",				gravity ) &&
		sSuccess == AlGetDouble( "twist",				twist ) &&
		sSuccess == AlGetDouble( "bend",				bend ) &&
		sSuccess == AlGetDouble( "trunk_taper",			trunk_taper ) &&
		sSuccess == AlGetDouble( "branch_length",		branch_length ) &&
		sSuccess == AlGetDouble( "branch_frequency",	branch_frequency ) &&
		sSuccess == AlGetDouble( "fork_angle",		fork_angle ) &&
		sSuccess == AlGetDouble( "min_branch_size",	min_branch_size ))
	{
		if( sSuccess == AlGetInteger( "seed", seed ) &&
			sSuccess == AlGetInteger( "max_branch", max_branch ))
		{
			// convert to radians
			twist = DegreesToRads( twist );
			fork_angle = DegreesToRads( fork_angle );
			branch_frequency /= 100.0;
			trunk_taper /= 100.0;
			return sSuccess;
		}
	}
	else
		AlPrintf( kPrompt, "Invalid size parameter!");
	return sFailure;
}



// builds a quadrateral
static int addQuad( AlPolyset *polyset, int v1, int v2, int v3, int v4, int p )
{
	int polygon_index = polyset->newPolygon();

	if( polygon_index != -1 )
	{
		AlPolygon *polygon = polyset->polygon( polygon_index );
		if( polygon )
		{
			polygon->addVertex( v1 );
			polygon->addVertex( v2 );
			polygon->addVertex( v3 );
			polygon->addVertex( v4 );
			if(p<0) p=0;
			if(p>12) p=12;
			polygon->setShaderIndex(p);
			delete polygon;
		}
	}
	else
		AlPrintf( kPrompt, "Error: cant create polygon");
	return polygon_index;
}


//builds a triangle
static int addTri( AlPolyset *polyset, int v1, int v2, int v3, int p)
{
	int polygon_index = polyset->newPolygon();

	if( polygon_index != -1 )
	{
		AlPolygon *polygon = polyset->polygon( polygon_index );
		if( polygon )
		{
			polygon->addVertex( v1 );
			polygon->addVertex( v2 );
			polygon->addVertex( v3 );
			if(p<0) p=0; 
			if(p>12) p=12;
			polygon->setShaderIndex(p);
		   
			delete polygon;
		}
	}
	else
		AlPrintf( kPrompt, "Error: cant create polygon");
	return polygon_index;
}

//



//normalizes a vector and returns the magnitude before normalization
float normalize(float *a, float *b, float *c) {
  float mag=sqrt(*a * *a + *b * *b + *c * *c);
  *a=*a/mag;
  *b=*b/mag;
  *c=*c/mag;
  return mag;
}

//no longer used
void fix(float *a) {
  *a=0;
  //  if(*a<-10) *a=-10;
  //if(*a>10) *a=10;
  //if(*a==0) *a=1;
}

//creates a blank attractor or repeller.
repeller repeller_new(void) {
  repeller a;
  for(int t=0;t<RMAX;t++) {
    a.x[t]=0;
    a.y[t]=0;
    a.z[t]=0;
    a.power[t]=0;
  }
  a.size=0;
  return a;
}

//returns a grammar defined as below
//commented out lines are remannts of past tests

grammar grammar_new(void) {
  grammar a;
  a.lines=5;

  char t[80]="f&&&f>>>>>>>f[K]>>[fK]>>[K]>>[fK]>>[K]>>[fK]";
  strcpy(a.value[0],t);
   a.key[0]='X';
  //  a.value[0]="f&&&f>>>>>>fK<<<<f<<K";
  
  //  a.value[0]="f&&&fK>>>>K>>>>K>>fK>>>>>>f>>fK>>>>>>f>>fK";
   //a.value[0]="f&&&f>>>>>>>f [K]>>[fK]>>[K]>>[fK]>>[K]>>[fK]";

  char u[80]="Y>>f<<fK";
  a.key[3]='K';
  //  a.value[3]="Y[>>f<<K]>>>>>>Y";
  strcpy(a.value[3],u);
  //a.value[3]="Y>>f<<fK";

  char v[80]="S<<f>>fY";
  strcpy(a.value[1],v);
  a.key[1]='Y';
  //  a.value[1]="S<<f>>fY";

  //  a.value[1]="STfY";
  //a.value[1]="Sf<<S>>Y";

  char w[80]="<<f>>f>>f>>f>>f>>f";
  strcpy(a.value[2],w);
  a.key[2]='S';
  //  a.value[2]="f>>>>>>f>>>>f";

  //  a.value[2]="<<f>>f>>f>>f>>f>>f";

  a.key[4]='T';
  char x[80]="f<<<<<<f<<<<f";
  strcpy(a.value[4],x);
  //  a.value[4]="f<<<<<<f<<<<f";

  //  a.value[2]="f>>>>>>f>>f>>f>>f>>f";

  //STRUCTURE
  //a.value[0]="f&&&[fnS]";
  //BLOSSOM 
  //a.value[0]="f&&&[fn+S]++[fn+S]++[fn+S]++[fn+S]++[fn+S]++[fn+S]";

  //[fff&&&fK>>>>>>K]";
  //a.value[3]="Y[>>ff<<K]>>>>>>Y";
  //a.value[1]="S[>>Sf<<S][<<S]fY";

  //a.value[2]="fS>>>>fS>>>>f";

  
  //  a.value[1]="f&&&&&&+++f++&&&&+f[ff^fS]++++++^^^^^^S";
  //STRUCTURE
  //a.value[1]="fS+-&&&&&&+++fS++&&&&+f+++^^^^S";
  //BLOSSOM
  //a.value[1]="fS&&&&&&+++f&&&&+++f+S";
  
  //  a.value[2]="++f++^f++^f++^f";
  //  a.value[1]="^F^F^F^F^F^F^F^F^F^F^F^F&F&F&F&F&F&F&F&F&F&F&F&F";
  //  a.value[1]="F&-F+F+F";
  
  return a;
}


//returns a new segment with the inputted values
//bx,by,bz refer to the bottom coords of the segment
//tx,ty,tz refer to the top coords
//top_width,bot_width refer to top and bottom widths for segments with thickness
//last_angle,last_rotate_angle,last_roll_angle are archaic and were used for
//the VARIANT of the program using x,y,z frame rotations....

segment segment_new(float bx, float by, float bz, float tx, float ty, float tz, float top_width, float bot_width,  int last_angle, int last_rotate_angle, int last_roll_angle) {
  segment a;
  //  bx=0;by=0;bz=0;
  a.bx=bx;a.by=by;a.bz=bz;
  a.tx=tx;a.ty=ty;a.tz=tz;

  a.length=sqrt((tx-bx)*(tx-bx) + (ty-by)*(ty-by) + (tz-bz)*(tz-bz));
  a.top_width=top_width;
  a.bot_width=bot_width;
  a.last_angle=last_angle;
  a.last_rotate_angle=last_rotate_angle;
  a.last_roll_angle=last_roll_angle;
  a.children=0;
  for(int temp=0;temp<CHILD_MAX;temp++) 
    a.child[temp]=NULL;
  //a.child[temp]=(segment*) malloc(sizeof(segment));
  return a;
}


//copies a segment 
segment segment_copy(segment *a) {
  segment b;
  b.bx=a->bx;
  b.by=a->by;
  b.bz=a->bz;
  b.tx=a->tx;
  b.ty=a->ty;
  b.tz=a->tz;
  b.length=a->length;
  b.top_width=a->top_width;
  b.bot_width=a->bot_width;
  b.last_angle=a->last_angle;
  b.last_rotate_angle=a->last_rotate_angle;
  b.last_roll_angle=a->last_roll_angle;
  b.children=0;
  for(int temp=0;temp<CHILD_MAX;temp++) 
    b.child[temp]=NULL;
  return b;
}


//USED exclusively for TREE generations... spanws CHILDREN number of
//segments from a parent until generation max is reached..
statusCode segment_spawn(segment *a, int num) {
  float bx,by,bz;
  float matrix_bot[9];
  float fork_COS=cos(M_PI*2/num);
  float fork_SIN=sin(M_PI*2/num);

  float ax=a->tx-a->bx;
  float ay=a->ty-a->by;
  float az=a->tz-a->bz;
  float len=normalize(&ax,&ay,&az);

  matrix_bot[0]=ax*ax*(1-fork_COS) + fork_COS;
  matrix_bot[1]=ay*ax*(1-fork_COS) - az*fork_SIN;
  matrix_bot[2]=az*ax*(1-fork_COS) + ay*fork_SIN;
  matrix_bot[3]=ax*ay*(1-fork_COS) + az*fork_SIN;
  matrix_bot[4]=ay*ay*(1-fork_COS) + fork_COS;
  matrix_bot[5]=az*ay*(1-fork_COS) - ax*fork_SIN;
  matrix_bot[6]=ax*az*(1-fork_COS) - ay*fork_SIN;
  matrix_bot[7]=ay*az*(1-fork_COS) + ax*fork_SIN;
  matrix_bot[8]=az*az*(1-fork_COS) + fork_COS;


  float tx=a->tx-a->bx;
  float ty=a->ty-a->by;
  float tz=a->tz-a->bz;



  float td=sqrt(tx*tx + ty*ty);
  if(td!=0) {
    tx=tx*(COS*td - SIN * tz)/(td);
    ty=ty*(COS*td - SIN * tz)/(td);
    tz=SIN*td + COS * tz;      
  } else {
    tx=SIN;
    ty=0;
    tz=COS;
  }

  normalize(&tx,&ty,&tz);

  tx*=len*LENGTH_RATIO;
  ty*=len*LENGTH_RATIO;
  tz*=len*LENGTH_RATIO;

  for(int temp=0;temp<num;temp++) {
    a->child[temp]=(segment*) malloc(sizeof(segment)); //returns new segment pointer
    bx=a->tx;
    by=a->ty;
    bz=a->tz;
  
    *a->child[temp]=segment_new(bx,by,bz, tx+a->tx, ty+a->ty, tz+a->tz, a->top_width*TRUNK_RATIO,a->top_width,  temp*2*M_PI/num+a->last_angle, 0,0); 


    float temp_x=matrix_bot[0]*tx + matrix_bot[1]*ty + matrix_bot[2]*tz;
    float temp_y=matrix_bot[3]*tx + matrix_bot[4]*ty + matrix_bot[5]*tz;
    float temp_z=matrix_bot[6]*tx + matrix_bot[7]*ty + matrix_bot[8]*tz;
    tx=temp_x;
    ty=temp_y;
    tz=temp_z;

    a->children++;
  }
  return(sSuccess);
}



//USED EXCLUSIVELY for the VARIANT version. this procedure positions a
//new segment remembering the last rotate, turn, and roll angles
//create toggles whether an actual segment should be created 'F' or
//whether only a movement in that direction should be done 'f'
statusCode segment_position(segment *a, float bx,float by,float bz, float tx, float ty, float tz, int angle, int rotate_angle, int roll_angle, float *fx, float *fy, float *fz, int create) {

  float *matrix_rotate=(float*)malloc(9*sizeof(float));
  float *matrix_turn=(float*)malloc(9*sizeof(float));
  float *matrix_roll=(float*)malloc(9*sizeof(float));
  float ax=1;
  float ay=0;
  float az=0;
  //angle=0;
  //rotate_angle=0;
  roll_angle=0;
  float tcos=cos(angle*ANGLE*M_PI/180);
  float tsin=sin(angle*ANGLE*M_PI/180);
  matrix_turn[0]=ax*ax*(1-tcos) + tcos;
  matrix_turn[1]=ay*ax*(1-tcos) - az*tsin;
  matrix_turn[2]=az*ax*(1-tcos) + ay*tsin;
  matrix_turn[3]=ax*ay*(1-tcos) + az*tsin;
  matrix_turn[4]=ay*ay*(1-tcos) + tcos;
  matrix_turn[5]=az*ay*(1-tcos) - ax*tsin;
  matrix_turn[6]=ax*az*(1-tcos) - ay*tsin;
  matrix_turn[7]=ay*az*(1-tcos) + ax*tsin;
  matrix_turn[8]=az*az*(1-tcos) + tcos;    

  ax=0;
  ay=1;
  az=0;
  float acos=cos(rotate_angle*ROTATE_ANGLE*M_PI/180);
  float asin=sin(rotate_angle*ROTATE_ANGLE*M_PI/180);
  matrix_rotate[0]=ax*ax*(1-acos) + acos;
  matrix_rotate[1]=ay*ax*(1-acos) - az*asin;
  matrix_rotate[2]=az*ax*(1-acos) + ay*asin;
  matrix_rotate[3]=ax*ay*(1-acos) + az*asin;
  matrix_rotate[4]=ay*ay*(1-acos) + acos;
  matrix_rotate[5]=az*ay*(1-acos) - ax*asin;
  matrix_rotate[6]=ax*az*(1-acos) - ay*asin;
  matrix_rotate[7]=ay*az*(1-acos) + ax*asin;
  matrix_rotate[8]=az*az*(1-acos) + acos;    

  ax=0;
  ay=0;
  az=1;
  float rcos=cos(roll_angle*ROLL_ANGLE*M_PI/180);
  float rsin=sin(roll_angle*ROLL_ANGLE*M_PI/180);
  matrix_roll[0]=ax*ax*(1-rcos) + rcos;
  matrix_roll[1]=ay*ax*(1-rcos) - az*rsin;
  matrix_roll[2]=az*ax*(1-rcos) + ay*rsin;
  matrix_roll[3]=ax*ay*(1-rcos) + az*rsin;
  matrix_roll[4]=ay*ay*(1-rcos) + rcos;
  matrix_roll[5]=az*ay*(1-rcos) - ax*rsin;
  matrix_roll[6]=ax*az*(1-rcos) - ay*rsin;
  matrix_roll[7]=ay*az*(1-rcos) + ax*rsin;
  matrix_roll[8]=az*az*(1-rcos) + rcos;    

  float length=LENGTH_RATIO*SHRINK_FACTOR*a->length;
  a->child[a->children]=(segment*) malloc(sizeof(segment)); //returns new segment pointer

  float gx=(tx-bx);
  float gy=(ty-by);
  float gz=(tz-bz);
  normalize(&gx,&gy,&gz);

  ax=matrix_turn[0]*gx + matrix_turn[1]*gy + matrix_turn[2]*gz;
  ay=matrix_turn[3]*gx + matrix_turn[4]*gy + matrix_turn[5]*gz;
  az=matrix_turn[6]*gx + matrix_turn[7]*gy + matrix_turn[8]*gz;

  gx=matrix_rotate[0]*ax + matrix_rotate[1]*ay + matrix_rotate[2]*az;
  gy=matrix_rotate[3]*ax + matrix_rotate[4]*ay + matrix_rotate[5]*az;
  gz=matrix_rotate[6]*ax + matrix_rotate[7]*ay + matrix_rotate[8]*az;

  ax=matrix_roll[0]*gx + matrix_roll[1]*gy + matrix_roll[2]*gz;
  ay=matrix_roll[3]*gx + matrix_roll[4]*gy + matrix_roll[5]*gz;
  az=matrix_roll[6]*gx + matrix_roll[7]*gy + matrix_roll[8]*gz;  


  ax*=length;
  ay*=length;
  az*=length;
  

  *fx=ax+tx;
  *fy=ay+ty;
  *fz=az+tz;
  
  // bx=0;
  // by=0;
  //bz=0;

  float ttx=ax+tx;
  float tty=ay+ty;
  float ttz=az+tz;

  free(matrix_rotate);
  free(matrix_turn);
  free(matrix_roll);

  //bx,by,bz    fx,fy,fz
  if(create) {
    *a->child[a->children]=segment_new(tx,ty,tz, ttx, tty, ttz, a->top_width,a->top_width, angle, rotate_angle, roll_angle);        
    a->children++;
  }
  return(sSuccess);
}



float square(float a) {
  return(a*a);
}


//a proxy call to segment_position for the VARIANT mode only
statusCode segment_empty_place(segment *a, float bx,float by,float bz, float tx, float ty, float tz, int angle, int rotate_angle, int roll_angle, float *fx, float *fy, float *fz) {

  segment_position(a,bx,by,bz,tx,ty,tz, angle,rotate_angle,roll_angle,fx,fy,fz,0);
  return(sSuccess);
}

//a proxy call to segemnt_position for the VARIANT mode only
statusCode segment_place(segment *a, float bx,float by,float bz, float tx, float ty, float tz, int angle, int rotate_angle, int roll_angle) {
  float fx,fy,fz;
  segment_position(a,bx,by,bz,tx,ty,tz,angle,rotate_angle,roll_angle,&fx,&fy,&fz,1);
  
  return(sSuccess);
}


//places a segment according the ALPHA, BETA, and GAMMA angle inputs
//attract toggles whether attractors and repellers should have an effect
//on the segment being placed. 0 means there is an affect, 1 means
//ignore attractors since a surface is being created which would be 
//unwantingly mutated. 
int segment_put(segment *a, repeller R, float* bx, float* by, float* bz, float* tx, float* ty, float* tz, int alpha, int beta, int gamma, int attract) {
  float ax=*tx-*bx;
  float ay=*ty-*by;
  float az=*tz-*bz;
  float A=(alpha*ALPHA)*M_PI/180;
  float B=beta*BETA*M_PI/180;
  float G=gamma*GAMMA*M_PI/180;
  
  if(ax>0) {
    if(ay>0) {
      //      if(az>=0) B=-B;
      
    } else {
      //  if(az>=0) B=-B;
      //else B=-B;
    }
    //    B=-B;
    if(az<-0.005||az>0.005) B=-B;
  } else {
    //B=-B;
    if(az<-0.005||az>0.005) B=-B;
    if(ay>0) {
      if(az<0) B=B;
    } else {
      if(az<0) {B=B;}
      else {};
    }
  }


  float len=normalize(&ax,&ay,&az);
  //len=1;
  float gx=ax;
  float gy=ay;
  float gz=az;

  float acos=cos(A);
  float asin=sin(A);
  float* matrix_alpha=(float *)malloc(9*sizeof(float));
  matrix_alpha[0]=ax*ax*(1-acos) + acos;
  matrix_alpha[1]=ay*ax*(1-acos) - az*asin;
  matrix_alpha[2]=az*ax*(1-acos) + ay*asin;
  matrix_alpha[3]=ax*ay*(1-acos) + az*asin;
  matrix_alpha[4]=ay*ay*(1-acos) + acos;
  matrix_alpha[5]=az*ay*(1-acos) - ax*asin;
  matrix_alpha[6]=ax*az*(1-acos) - ay*asin;
  matrix_alpha[7]=ay*az*(1-acos) + ax*asin;
  matrix_alpha[8]=az*az*(1-acos) + acos;    


  float bcos=cos(B);
  float bsin=sin(B);

  float td=sqrt(gx*gx + gy*gy);
  if(td!=0) {
    gx=gx*(bcos*td - bsin * gz)/(td);
    gy=gy*(bcos*td - bsin * gz)/(td);
    gz=bsin*td + bcos*gz;      
  } else {
    gx=0;
    gy=bsin;
    gz=bcos;
  }

  if(td!=0) {
    ax=-ax*az/td;
    ay=-ay*az/td;
    az=td;
  } else {
    ax=0;
    ay=1;
    az=0;
  }

  float gcos=cos(G);
  float gsin=sin(G);
  float* matrix_gamma=(float *)malloc(9*sizeof(float));
  matrix_gamma[0]=ax*ax*(1-gcos) + gcos;
  matrix_gamma[1]=ay*ax*(1-gcos) - az*gsin;
  matrix_gamma[2]=az*ax*(1-gcos) + ay*gsin;
  matrix_gamma[3]=ax*ay*(1-gcos) + az*gsin;
  matrix_gamma[4]=ay*ay*(1-gcos) + gcos;
  matrix_gamma[5]=az*ay*(1-gcos) - ax*gsin;
  matrix_gamma[6]=ax*az*(1-gcos) - ay*gsin;
  matrix_gamma[7]=ay*az*(1-gcos) + ax*gsin;
  matrix_gamma[8]=az*az*(1-gcos) + gcos;    


  normalize(&gx,&gy,&gz);

  float temp_x=matrix_alpha[0]*gx + matrix_alpha[1]*gy + matrix_alpha[2]*gz;
  float temp_y=matrix_alpha[3]*gx + matrix_alpha[4]*gy + matrix_alpha[5]*gz;
  float temp_z=matrix_alpha[6]*gx + matrix_alpha[7]*gy + matrix_alpha[8]*gz;

  gx=matrix_gamma[0]*temp_x + matrix_gamma[1]*temp_y + matrix_gamma[2]*temp_z;
  gy=matrix_gamma[3]*temp_x + matrix_gamma[4]*temp_y + matrix_gamma[5]*temp_z;
  gz=matrix_gamma[6]*temp_x + matrix_gamma[7]*temp_y + matrix_gamma[8]*temp_z;
  

  float fx,fy,fz;
  //****EFFECT OF ATTRACTOR REPELLER****
  if(R.size>0&&attract==0&&FILTER3==0) {
    float net_x=0;
    float net_y=0;
    float net_z=0;
    for(int t=0;t<R.size;t++) {
      float dist=((*bx-R.x[t])*(*bx-R.x[t]) + (*by-R.y[t])*(*by-R.y[t]) + (*bz-R.z[t])*(*bz-R.z[t]));
      net_x+=R.power[t]*(R.x[t]-*bx)/dist;
      net_y+=R.power[t]*(R.y[t]-*by)/dist;
      net_z+=R.power[t]*(R.z[t]-*bz)/dist;
    }
    //    float mag=normalize(&net_x,&net_y,&net_z);

  

    gx=gx + 0.4*net_x;
    gy=gy + 0.4*net_y;
    gz=gz + 0.4*net_z;
    normalize(&gx,&gy,&gz);
  }

  gx=len*gx;
  gy=len*gy;
  gz=len*gz;

  
  a->child[a->children]=(segment*) malloc(sizeof(segment)); //returns new segment pointer

  gx=gx + *tx;
  gy=gy + *ty;
  gz=gz + *tz;


  //tx,ty,tz is the new bottom 
  //gx,gy,gz is the new top

  //*****************************************
  //apply boundary conditions as last step...
  //*****************************************

  int limit=0;
  float k=0;
  
  if(gx<X1) {k=(X1-gx)/(*tx-gx);limit=1;}
  if(gx>X2) {k=(X2-gx)/(*tx-gx);limit=1;}
  if(gy<Y1) {k=(Y1-gy)/(*ty-gy);limit=1;}
  if(gy>Y2) {k=(Y2-gy)/(*ty-gy);limit=1;}
  if(gz>Z2) {k=(Z2-gz)/(*tz-gz);limit=1;}
  if(gz<Z1) {k=(Z1-gz)/(*tz-gz);limit=1;}
  
  
  //if(k<0) k=0;
  //if(k>1) k=1;
  //k=0;
  //gx-= k* *tx;
  //gy-= k* *ty;
  //gz-= k* *tz; 
  
  
  float x=gx-*tx;
  float y=gy-*ty;
  float z=gz-*tz;

  /*  if(gz<Z1) {
    float a1=x*x + y*y;
    float b1=2*x* *tx + 2*y* *ty;
    float c1=*tx**tx + *ty**ty - len + square(Z1-*tz);
   k=(-b1+sqrt(b1*b1-4*a1*c1))/(2*a1);
    gz=Z1;
    gx=*tx+k*x;
    gy=*ty+k*y;
  }*/
  //  if(gz<Z1) {len-square(Z1-tz)}


  //  if(create) {
  *a->child[a->children]=segment_new(*tx,*ty,*tz, gx, gy, gz, a->top_width,a->top_width, 0, 0, 0);        
  a->children++;
  //}
  *bx=*tx;*by=*ty;*bz=*tz;
  *tx=gx;*ty=gy;*tz=gz;
  return(limit);
}


//cleansup a segment and all its children from memory
statusCode segment_cleanup(segment *a) {
  for(int temp=0;temp<a->children;temp++) {
    if(a->child[temp]!=NULL) segment_cleanup(a->child[temp]);
    a->child[temp]=NULL;
  }
  free(a);
  return(sSuccess);
}

//defines a BLACK Alshader
statusCode make_BLACK(AlPolyset *polyset) {
  AlShader shader;
  shader.create();
  shader.setParameter( kFLD_SHADING_COMMON_COLOR_R, 255);
  shader.setParameter( kFLD_SHADING_COMMON_COLOR_G, 255 );
  shader.setParameter( kFLD_SHADING_COMMON_COLOR_B, 255 );
  polyset->assignShader(12,&shader);

  return(sSuccess);
}

statusCode make_RED(AlPolyset *polyset) {
  AlShader shader;
  shader.create();
  shader.setParameter( kFLD_SHADING_COMMON_COLOR_R, 255);
  shader.setParameter( kFLD_SHADING_COMMON_COLOR_G, 0 );
  shader.setParameter( kFLD_SHADING_COMMON_COLOR_B, 255);
  polyset->assignShader(RED,&shader);

  return(sSuccess);
}


//defines a set of Alshaders from red to blue.
statusCode make_COLORS(AlPolyset *polyset) {

  int NUM=11;
  for(int temp=0;temp<NUM;temp++) {
    AlShader shader;
    shader.create();
    if(temp!=0) {
      //    shader.setParameter( kFLD_SHADING_COMMON_COLOR_R, 255-temp*(255/NUM));
      shader.setParameter( kFLD_SHADING_COMMON_COLOR_R, 255);
      shader.setParameter( kFLD_SHADING_COMMON_COLOR_G, temp*(255/NUM));
      //    shader.setParameter( kFLD_SHADING_COMMON_COLOR_B, temp*(255/NUM));
      shader.setParameter( kFLD_SHADING_COMMON_COLOR_B, 0);
    } else {
      shader.setParameter( kFLD_SHADING_COMMON_COLOR_R, 0);
      shader.setParameter( kFLD_SHADING_COMMON_COLOR_G, 0);
      shader.setParameter( kFLD_SHADING_COMMON_COLOR_B, 0);
    }
    polyset->assignShader(temp,&shader);
  }
  return(sSuccess);
}


//draws a segment with thickness (not surfaces)
statusCode segment_draw(segment *a, AlPolyset *polyset, int level) {
  if(level>GENERATION_LIMIT||TREE>=0) {
    printf("*");
    float ax=(a->tx - a->bx);
    float ay=(a->ty - a->by);
    float az=(a->tz - a->bz);
    normalize(&ax,&ay,&az);
    
    // base vectors
    float bx=az;
    float by=0;
    float bz=-ax;
    if(ax==00&&az==0) {
      bx=0;
      by=-ax;
      bz=ay;
    }
    normalize(&bx,&by,&bz);
    
    
    //rotation_matrices
    float *matrix_top=(float*)malloc(9*sizeof(float));
    
    
    matrix_top[0]=ax*ax*(1-RCOS) + RCOS;
    matrix_top[1]=ay*ax*(1-RCOS) - az*RSIN;
    matrix_top[2]=az*ax*(1-RCOS) + ay*RSIN;
    matrix_top[3]=ax*ay*(1-RCOS) + az*RSIN;
    matrix_top[4]=ay*ay*(1-RCOS) + RCOS;
    matrix_top[5]=az*ay*(1-RCOS) - ax*RSIN;
    matrix_top[6]=ax*az*(1-RCOS) - ay*RSIN;
    matrix_top[7]=ay*az*(1-RCOS) + ax*RSIN;
    matrix_top[8]=az*az*(1-RCOS) + RCOS;
    
      
    int *base=(int*)malloc(4*sizeof(int));
    
    for(int cnt=0;cnt<SIDES;cnt++) {     
      printf("1");
      //      printf("\n<%f %f %f %f %f %f %f %f>\n",a->tx,a->ty,a->tz,bx,by,bz,a->top_width,a->bot_width);
      
      float gx=a->tx+bx*a->top_width;
      float gy=a->ty+by*a->top_width;
      float gz=a->tz+bz*a->top_width;
      printf("<%f %f %f | %f %f %f>\n",a->tx,a->ty,a->tz,gx,gy,gz);
      base[0] = polyset->newVertex(gx, gy, gz);
      
      printf(".");
      base[1] = polyset->newVertex(  a->bx+bx*a->bot_width, a->by+by*a->bot_width, a->bz+bz*a->bot_width);
      
      printf("2");
      //rotate (360/SIDES) degrees
      
      //top
      float tx=matrix_top[0]*bx + matrix_top[1]*by + matrix_top[2]*bz;
      float ty=matrix_top[3]*bx + matrix_top[4]*by + matrix_top[5]*bz;
      float tz=matrix_top[6]*bx + matrix_top[7]*by + matrix_top[8]*bz;
      bx=tx;
      by=ty;
      bz=tz;
      
      printf("3");
      base[2] = polyset->newVertex(  a->bx+bx*a->bot_width, a->by+by*a->bot_width, a->bz+bz*a->bot_width);
      printf(".");
      base[3] = polyset->newVertex(  a->tx+bx*a->top_width, a->ty+by*a->top_width, a->tz+bz*a->top_width);
      
      printf("4");
      addQuad( polyset, base[0], base[1], base[2], base[3], 5 );  
      
      printf("5");
    }
    free(base);
    //    if(level<2) AlUniverse::redrawScreen();
  }
  return(sSuccess); 
}



//draw a box at a point

void draw_box(AlPolyset *polyset) {
  double x4,y4,z4;
  AlPolysetVertex* vertex;
  int base[8];
  int p=5;
  float offset=0.2;
  int loops=polyset->numberOfVertices();
  //int loops=8;
  
  for(int w=0;w<loops;w=w+6) {
    vertex=polyset->vertex(w);
    vertex->worldPosition(x4,y4,z4);
    base[0]=polyset->newVertex(x4-offset,y4-offset,z4-offset);
    base[1]=polyset->newVertex(x4+offset,y4-offset,z4-offset);
    base[2]=polyset->newVertex(x4+offset,y4+offset,z4-offset);
    base[3]=polyset->newVertex(x4-offset,y4+offset,z4-offset);
    
    base[4]=polyset->newVertex(x4-offset,y4-offset,z4+offset);
    base[5]=polyset->newVertex(x4+offset,y4-offset,z4+offset);
    base[6]=polyset->newVertex(x4+offset,y4+offset,z4+offset);
    base[7]=polyset->newVertex(x4-offset,y4+offset,z4+offset);

    addQuad( polyset, base[0], base[1], base[2], base[3], p);
    addQuad( polyset, base[4], base[5], base[6], base[7], p);
    addQuad( polyset, base[0], base[1], base[4], base[5], p);
    addQuad( polyset, base[2], base[3], base[7], base[6], p);
    addQuad( polyset, base[3], base[0], base[4], base[7], p);
    addQuad( polyset, base[1], base[2], base[6], base[5], p);
  }
}



//draw the bounding box

void draw_bounding(AlPolyset *polyset) {
  int base[48];
  int p=0;
  float offset=((sqrt((((X1-X2)*(X1-X2))+((Y1-Y2)*(Y1-Y2)))+((Z1-Z2)*(Z1-Z2))))/80);
  //boxpoints
  base[0]=polyset->newVertex(X1,Y1,Z1);
  base[1]=polyset->newVertex(X2,Y1,Z1);
  base[2]=polyset->newVertex(X2,Y2,Z1);
  base[3]=polyset->newVertex(X1,Y2,Z1);
  
  base[4]=polyset->newVertex(X1,Y1,Z2);
  base[5]=polyset->newVertex(X2,Y1,Z2);
  base[6]=polyset->newVertex(X2,Y2,Z2);
  base[7]=polyset->newVertex(X1,Y2,Z2);

  //lowerface
  base[8]=polyset->newVertex(X1+offset,Y1+offset,Z1);
  base[9]=polyset->newVertex(X2-offset,Y1+offset,Z1);
  base[10]=polyset->newVertex(X2-offset,Y2-offset,Z1);
  base[11]=polyset->newVertex(X1+offset,Y2-offset,Z1);
  addQuad( polyset, base[0], base[1], base[9], base[8], p);
  addQuad( polyset, base[1], base[2], base[10], base[9], p);
  addQuad( polyset, base[2], base[3], base[11], base[10], p);
  addQuad( polyset, base[3], base[0], base[8], base[11], p);

  //rightface
  base[12]=polyset->newVertex(X1+offset,Y1,Z1+offset);
  base[13]=polyset->newVertex(X2-offset,Y1,Z1+offset);
  base[14]=polyset->newVertex(X2-offset,Y1,Z2-offset);
  base[15]=polyset->newVertex(X1+offset,Y1,Z2-offset);
  addQuad( polyset, base[0], base[1], base[13], base[12], p);
  addQuad( polyset, base[1], base[5], base[14], base[13], p);
  addQuad( polyset, base[5], base[4], base[15], base[14], p);
  addQuad( polyset, base[4], base[0], base[12], base[15], p);

  //leftface
  base[16]=polyset->newVertex(X1+offset,Y2,Z1+offset);
  base[17]=polyset->newVertex(X2-offset,Y2,Z1+offset);
  base[18]=polyset->newVertex(X2-offset,Y2,Z2-offset);
  base[19]=polyset->newVertex(X1+offset,Y2,Z2-offset);
  addQuad( polyset, base[3], base[2], base[17], base[16], p);
  addQuad( polyset, base[2], base[6], base[18], base[17], p);
  addQuad( polyset, base[6], base[7], base[19], base[18], p);
  addQuad( polyset, base[7], base[3], base[16], base[19], p);

  //topface
  base[20]=polyset->newVertex(X1+offset,Y1+offset,Z2);
  base[21]=polyset->newVertex(X2-offset,Y1+offset,Z2);
  base[22]=polyset->newVertex(X2-offset,Y2-offset,Z2);
  base[23]=polyset->newVertex(X1+offset,Y2-offset,Z2);
  addQuad( polyset, base[4], base[5], base[21], base[20], p);
  addQuad( polyset, base[5], base[6], base[22], base[21], p);
  addQuad( polyset, base[6], base[7], base[23], base[22], p);
  addQuad( polyset, base[7], base[4], base[20], base[23], p);

  //frontface
  base[24]=polyset->newVertex(X1,Y1+offset,Z1+offset);
  base[25]=polyset->newVertex(X1,Y1+offset,Z2-offset);
  base[26]=polyset->newVertex(X1,Y2-offset,Z2-offset);
  base[27]=polyset->newVertex(X1,Y2-offset,Z1+offset);
  addQuad( polyset, base[0], base[4], base[25], base[24], p);
  addQuad( polyset, base[4], base[7], base[26], base[25], p);
  addQuad( polyset, base[7], base[3], base[27], base[26], p);
  addQuad( polyset, base[3], base[0], base[24], base[27], p);

  //backface
  base[28]=polyset->newVertex(X2,Y1+offset,Z1+offset);
  base[29]=polyset->newVertex(X2,Y1+offset,Z2-offset);
  base[30]=polyset->newVertex(X2,Y2-offset,Z2-offset);
  base[31]=polyset->newVertex(X2,Y2-offset,Z1+offset);
  addQuad( polyset, base[1], base[5], base[29], base[28], p);
  addQuad( polyset, base[5], base[6], base[30], base[29], p);
  addQuad( polyset, base[6], base[2], base[31], base[30], p);
  addQuad( polyset, base[2], base[1], base[28], base[31], p);

}




//filters the drawn L-system by averaging nearby vertices and repositioning
//them...
statusCode filter2(AlPolyset *polyset) {

  //***************************
  //find nearest vertex of differnet polygon with same generation number
  //find nearest vertex of different polygon with 1+/- generation number
  //fix them in that order


  AlPolysetVertex* vertex;
  AlPolysetVertex* vertex2;
  AlPolygon* polygon;
  double x,y,z;

  int numPolygons=polyset->numberOfPolygons();
  int numVertices=polyset->numberOfVertices();
 
  
  int v=0;
  for(int index=0;index<numPolygons;index++) {
    polygon=polyset->polygon(index);
    int vnum=polygon->numberOfVertices();
    for(int v2=0;v2<vnum;v2++) {
      vertex=polygon->vertex(v2);
      vertex->worldPosition(x,y,z);
      V[v].x=x;
      V[v].y=y;
      V[v].z=z;
      V[v].fixed=0;
      V[v].vertex=v2;
      V[v].polygon=index;
      v++;
    }
  }

  int frac=numPolygons/20;
  for(index=0;index<numPolygons;index++) {
    if(index % frac==0) 
      AlPrintf( kPrompt, "Filter2ing %d percent %d %d complete", (100*index)/numPolygons, numPolygons, numVertices);
    polygon=polyset->polygon(index);
    int vnum=polygon->numberOfVertices();    
    int base[100];
    float color=0;
    for(v=0;v<vnum;v++) {      
      vertex=polygon->vertex(v);
      vertex->worldPosition(x,y,z);

      int vert[100];
      int numlist=0;
      float disavg=0;
      for(int v2=0;v2<numVertices;v2++) {
        double x2,y2,z2;
        x2=V[v2].x;
        y2=V[v2].y;
        z2=V[v2].z;
        //      vertex2=polyset->vertex(v2);
        //vertex2->worldPosition(x2,y2,z2);
        float dis=(x-x2)*(x-x2) + (y-y2)*(y-y2);
        if(dis<0.2&&((z-z2)*(z-z2)<1)) {
           //&&((z-z2)*(z-z2)<0.5)) {
          disavg+=dis;
          vert[numlist]=v2;
          numlist++;
        }
      }
      printf("%d\n",numlist);
      disavg=disavg/numlist;
      x=0;
      y=0;
      z=0;
      color+=(disavg*2*10)/1;
      for(int i=0;i<numlist;i++) {
        double x2,y2,z2;
        vertex2=polyset->vertex(vert[i]);
        vertex2->worldPosition(x2,y2,z2);
        x+=x2/numlist;
        y+=y2/numlist;
        z+=z2/numlist;
      }
      base[v]=polyset->newVertex(x,y,z);
    }
    
    
    //MAKE POLYGON FROM new vertices
    AlPolygon *polygon=polyset->polygon(polyset->newPolygon());
    for(v=0;v<vnum;v++) {
      polygon->addVertex(base[v]);
    }
    if(color>9) color=9;
    polygon->setShaderIndex((int)color);
    delete polygon;
  }

  AlPrintf( kPrompt, "Removing old polygons and vertices...");
  //now we delete the old polygons (0..numPolygons-1)  
  //note: new polygons are (numPolygons...numPolygons*2-1)
  for(index=0;index<numPolygons;index++) {
    polygon=polyset->polygon(index);
    polyset->deletePolygon(polygon->index());
  }


  draw_bounding(polyset);

  //now we delete the old vertices since they are obsolete...
  for(index=0;index<numVertices;index++) {
    vertex=polyset->vertex(index);
    polyset->deleteVertex(vertex->index());
  }

  AlPrintf( kPrompt, "Filter2 process complete");
  return sSuccess;
}

/*** NEW FILTER IDEA - EXPERIMENTAL ***/
statusCode filter(AlPolyset *polyset) {

  //***************************
  //find nearest vertex of differnet polygon with same generation number
  //find nearest vertex of different polygon with 1+/- generation number
  //fix them in that order


  AlPolysetVertex* vertex;
  AlPolygon* polygon;
  double x,y,z;

  int numPolygons=polyset->numberOfPolygons();
  int numVertices=polyset->numberOfVertices();
   
  int v=0;
  for(int index=0;index<numPolygons;index++) {
    polygon=polyset->polygon(index);
    int vnum=polygon->numberOfVertices();
    for(int v2=0;v2<vnum;v2++) {
      vertex=polygon->vertex(v2);
      vertex->worldPosition(x,y,z);
      V[v].x=x;
      V[v].y=y;
      V[v].z=z;
      V[v].fixed=0;
      V[v].vertex=v2;
      V[v].polygon=index;
      v++;
    }
  }

  float color[15000];
  float best=1000;
  int frac=numVertices/20;
  for(v=0;v<numVertices;v++) {
    if(v % frac==0) 
      AlPrintf( kPrompt, "Experimenting %d percent %d %d complete", (100*v)/numVertices, v, numVertices);
    if(!V[v].fixed) {
      int i=-1;
      float ceffect=0;
      best=1000;
      for(int j=0;j<numVertices;j++) if(V[j].polygon!=V[v].polygon && V[j].fixed==0) {
	float val=square(V[v].x-V[j].x)+square(V[v].y-V[j].y);
	if(val<best&&square(V[v].z-V[j].z)<0.5) {
	  best=val;
	  i=j;index=V[j].polygon;
	}
      }
      ceffect+=best;
      if(best>threshhold) i=-1; 
      best=1000;
      int index2=0;
      int i2=-1;
      if(i>0) for(j=0;j<numVertices;j++) if(V[j].polygon!=V[v].polygon && V[j].polygon!=index && V[j].fixed==0)
	{
	  float val=square(V[v].x-V[j].x)+square(V[v].y-V[j].y);
	  if(val<best&&square(V[v].z-V[j].z)<0.5) {
	    best=val;
	    i2=j;index2=V[j].polygon;
	  }
	}
      if(shared<=2) i2=-1;
      if(best>threshhold) i2=-1;
      ceffect+=best;
      best=1000;
      int i3=-1;
      if(i2>0) for(j=0;j<numVertices;j++) if(V[j].polygon!=V[v].polygon && V[j].polygon!=index && V[j].polygon!=index2 && V[j].fixed==0)
	{
	  float val=square(V[v].x-V[j].x)+square(V[v].y-V[j].y);
	  if(val<best&&square(V[v].z-V[j].z)<0.5) {
	    best=val;
	    i3=j;
	  }
	}
      if(shared<=3) i3=-1;
      if(best>threshhold) i3=-1;
      ceffect+=best;

      color[v]=ceffect;
      V[v].fixed=1;      
      int num=1;
      if(i>0) {
	num++;
	color[i]=ceffect;
	V[v].x+=V[i].x;V[v].y+=V[i].y;V[v].z+=V[i].z;
      }
      if(i2>0) {
	num++;
	color[i2]=ceffect;
	V[v].x+=V[i2].x;V[v].y+=V[i2].y;V[v].z+=V[i2].z;
      }
      if(i3>0) {
	num++;
	color[i3]=ceffect;
	V[v].x+=V[i3].x;V[v].y+=V[i3].y;V[v].z+=V[i3].z;
      }
      V[v].x/=num;
      V[v].y/=num;
      V[v].z/=num;
      
      if(i>0) {
	V[i].x=V[v].x;
	V[i].y=V[v].y;
	V[i].z=V[v].z;
	V[i].fixed=1;
      }
      if(i2>0) {
	V[i2].x=V[v].x;      
	V[i2].y=V[v].y;
	V[i2].z=V[v].z;
	V[i2].fixed=1;
      }
      if(i3>0) {
	V[i3].x=V[v].x;      
	V[i3].y=V[v].y;
	V[i3].z=V[v].z;
	V[i3].fixed=1;
      }
    }
    
  }
  
  int base[50];
  int temp=0;
  
  for(v=0;v<numVertices;v++) {
    base[temp]=polyset->newVertex(V[v].x,V[v].y,V[v].z);
    temp++;
    int paint=5;
    if(temp % FACET_SIZE==0) {  
      temp=0;
      AlPolygon *polygon=polyset->polygon(polyset->newPolygon());
      for(int v2=0;v2<FACET_SIZE;v2++) {
	polygon->addVertex(base[v2]);
      }
      paint=(int)(9-color[v]*10);
      if(paint>11) paint=11;
      if(paint<1) paint=1;
      polygon->setShaderIndex(paint);
      delete polygon;    
    }
  }
  
  
  AlPrintf( kPrompt, "Removing old polygons and vertices...");
  //now we delete the old polygons (0..numPolygons-1)  
  //note: new polygons are (numPolygons...numPolygons*2-1)
  for(index=0;index<numPolygons;index++) {
    polygon=polyset->polygon(index);
    polyset->deletePolygon(polygon->index());
  }
  draw_bounding(polyset);
  //  draw_box(polyset);

    //now we delete the old vertices since they are obsolete...
  for(index=0;index<numVertices;index++) {
    vertex=polyset->vertex(index);
    polyset->deleteVertex(vertex->index());
  }
  
  AlPrintf( kPrompt, "Experiment complete");
  return sSuccess;
}



/*** OFFSET AFTERWARDS FILTER PROPOSED BY AKILIAN ***/
statusCode filter3(AlPolyset *polyset) {

  //***************************
  //find nearest vertex of differnet polygon with same generation number
  //find nearest vertex of different polygon with 1+/- generation number
  //fix them in that order

  AlPolysetVertex* vertex;
  AlPolygon* polygon;
  double x,y,z;

  int numPolygons=polyset->numberOfPolygons();
  int numVertices=polyset->numberOfVertices();
  links=0;
  int index=0;
  float best=1000;
  int frac=numVertices/20;
  for(int v=0;v<Vnum;v++) {
    if(v % frac==0) 
      AlPrintf( kPrompt, "links=%d Linking %d percent %d %d complete", links,(100*v)/numVertices, v, numVertices);
    if(!V[v].fixed) {
      int i=-1;
      float ceffect=0;
      best=1000;
      for(int j=0;j<numVertices;j++) if(V[j].polygon!=V[v].polygon && V[j].fixed==0) {
	float val=square(V[v].x-V[j].x)+square(V[v].y-V[j].y)+square(V[v].z-V[j].z);
	if(val<best) {
	  best=val;
	  i=j;index=V[j].polygon;
	}
      }
      ceffect+=best;
      if(best>threshhold) i=-1; 
      best=1000;
      int index2=0;
      int i2=-1;
      if(i>0) for(j=0;j<numVertices;j++) if(V[j].polygon!=V[v].polygon && V[j].polygon!=index && V[j].fixed==0)
	{
	  float val=square(V[v].x-V[j].x)+square(V[v].y-V[j].y)+square(V[v].z-V[j].z);
	  if(val<best) {
	    best=val;
	    i2=j;index2=V[j].polygon;
	  }
	}
      if(shared<=2) i2=-1;
      if(best>threshhold) i2=-1;
      ceffect+=best;
      best=1000;
      int i3=-1;
      if(i2>0) for(j=0;j<numVertices;j++) if(V[j].polygon!=V[v].polygon && V[j].polygon!=index && V[j].polygon!=index2 && V[j].fixed==0)
	{
	  float val=square(V[v].x-V[j].x)+square(V[v].y-V[j].y)+square(V[v].z-V[j].z);
	  if(val<best) {
	    best=val;
	    i3=j;
	  }
	}
      if(shared<=3) i3=-1;
      if(best>threshhold) i3=-1;
      ceffect+=best;

      color[v]=ceffect;
      if(i>0) {
	color[i]=ceffect;
      }
      if(i2>0) {
	color[i2]=ceffect;
      }
      if(i3>0) {
	color[i3]=ceffect;
      }

      link[links].v[0]=v;
      link[links].vnum=1;
      V[v].fixed=1;      
      if(i>0) {link[links].v[1]=i;link[links].vnum++;V[i].fixed=1;}
      if(i2>0) {link[links].v[2]=i2;link[links].vnum++;V[i2].fixed=1;}
      if(i3>0) {link[links].v[3]=i3;link[links].vnum++;V[i3].fixed=1;}
      //printf("#%d - %d\n",links,link[links].vnum);
      links++; //incremenet number of links...
    }
    
  }
  AlDebug::cleanUpUniverse();

  return sSuccess;
}

statusCode filter3b(AlPolyset *polyset) {
  int base[50];
  int temp=0;
  int numVertices=polyset->numberOfVertices();
  int numPolygons=polyset->numberOfPolygons();


  for(int c=0;c<links;c++) {
    float tx=0;
    float ty=0;
    float tz=0;
    for(int t=0;t<link[c].vnum;t++) {
      tx+=V[link[c].v[t]].x;
      ty+=V[link[c].v[t]].y;
      tz+=V[link[c].v[t]].z;
    }
    tx/=link[c].vnum;
    ty/=link[c].vnum;
    tz/=link[c].vnum;
    for(t=0;t<link[c].vnum;t++) {
      V[link[c].v[t]].x=tx;
      V[link[c].v[t]].y=ty;
      V[link[c].v[t]].z=tz;
    }
    
  }

  AlPrintf( kPrompt, "links=%d Vnum=%d..",links,Vnum);
  for(int v=0;v<Vnum;v++) {
    base[temp]=polyset->newVertex(V[v].x,V[v].y,V[v].z);
    temp++;
    int paint=5;
    if(temp % FACET_SIZE==0) {  
      temp=0;
      AlPolygon *polygon=polyset->polygon(polyset->newPolygon());
      for(int v2=0;v2<FACET_SIZE;v2++) {
	polygon->addVertex(base[v2]);
      }
      paint=(int)(9-color[v]*10);
      if(paint>11) paint=11;
      if(paint<1) paint=1;
      polygon->setShaderIndex(paint);
      delete polygon;    
    }
  }
  
  AlPolygon *polygon;
  AlPolysetVertex *vertex;
  //  AlPrintf( kPrompt, "Removing old polygons and vertices...");
  //now we delete the old polygons (0..numPolygons-1)  
  //note: new polygons are (numPolygons...numPolygons*2-1)
  for(int index=0;index<numPolygons;index++) {
    polygon=polyset->polygon(index);
    polyset->deletePolygon(polygon->index());
  }
  draw_bounding(polyset);
  //  draw_box(polyset);

    //now we delete the old vertices since they are obsolete...
  for(index=0;index<numVertices;index++) {
    vertex=polyset->vertex(index);
    polyset->deleteVertex(vertex->index());
  }
  
  AlPrintf( kPrompt, "Filter-3 complete");
  return sSuccess;
}

//creates an overlying mesh over a drawn L-system structure... an x/y plane
//projection...
statusCode mesh(AlPolyset *polyset) {
  struct Vert v[MAX]; 
  int num=0;

  AlPolysetVertex* vertex;
  AlPolygon* polygon;

  int numPolygons=polyset->numberOfPolygons();
  int numVertices=polyset->numberOfVertices();

  for(int index=0;index<numPolygons;index++) {
    polygon=polyset->polygon(index);
    int vnum=polygon->numberOfVertices();
    /*    v[num].x=0;
    v[num].y=0;
    v[num].z=0;*/
    double x,y,z;
    for(int i=0;i<vnum;i++) {
      vertex=polygon->vertex(i);
      vertex->worldPosition(x,y,z);
      v[num].x=x;
      v[num].y=y;
      v[num].z=z;
      num++;
    }
    //    num++;
  }

  struct Vert z[200][200];
  int xpos=0;
  int ypos=0;
  int size=0;
  int frac=50;

  for(float x=-6;x<6;x=x+0.2) {
    int temp=100*x;
    if(temp % frac==0) 
      AlPrintf( kPrompt, "Meshing %d percent complete", (100*xpos)/30);
    for(float y=-6;y<6;y=y+0.2) {
      z[xpos][ypos].x=x;
      z[xpos][ypos].y=y;
      int win[TOP];
      float best[TOP];
      for(int t=0;t<TOP;t++) {
	win[t]=-100;best[t]=1000;
      }
      for(index=0;index<num;index++) {
	float dis=(x-v[index].x)*(x-v[index].x) + (y-v[index].y)*(y-v[index].y);
	//best[TOP-1] contains the closest, best[0] the farthest that is acceptible
	//if(dis>1) dis=10000;
	if(dis<best[0]) {
	  win[0]=index;
	  best[0]=dis;
	  int pos=TOP-1;
	  for(int t=1;t<TOP;t++) {
	    if(dis>best[t]) {pos=t-1;t=TOP;};
	  }
	  for(int u=0;u<pos;u++)
	    {best[u]=best[u+1];win[u]=win[u+1];}
	  best[pos]=dis;
	  win[pos]=index;
	}
      }
      z[xpos][ypos].z=0;
      for(t=0;t<TOP;t++) {
	if(best[t]<2) z[xpos][ypos].z+=v[win[t]].z/TOP;
	else z[xpos][ypos].z=-5;
      }
      ypos++;
    }
    xpos++;
    ypos=0;
  }
  size=xpos;
  
  int base[4];
  
  for(int X=0;X<size-1;X++)
    for(int Y=0;Y<size-1;Y++) {
      base[0]=polyset->newVertex(z[X][Y].x,z[X][Y].y,z[X][Y].z);
      base[1]=polyset->newVertex(z[X+1][Y].x,z[X+1][Y].y,z[X+1][Y].z);
      base[2]=polyset->newVertex(z[X+1][Y+1].x,z[X+1][Y+1].y,z[X+1][Y+1].z);
      base[3]=polyset->newVertex(z[X][Y+1].x,z[X][Y+1].y,z[X][Y+1].z);
      
      //MAKE POLYGON FROM new vertices
      AlPolygon *polygon=polyset->polygon(polyset->newPolygon());
      for(int i=0;i<4;i++) {
	polygon->addVertex(base[i]);
      }
      int color=9;
      polygon->setShaderIndex(color);
      delete polygon;
    }
  
  AlPrintf( kPrompt, "Removing old polygons and vertices...");
  //now we delete the old polygons (0..numPolygons-1)  
  for(index=0;index<numPolygons;index++) {
    polygon=polyset->polygon(index);
    polyset->deletePolygon(polygon->index());
  }
  //now we delete the old vertices...
  for(index=0;index<numVertices;index++) {
    vertex=polyset->vertex(index);
    polyset->deleteVertex(vertex->index());
  }
  AlPrintf( kPrompt, "Meshing complete");
  return sSuccess;
}



//draw repellers as octahedrons of various sizes and colors
//the stronger it is, the bigger the octahedron
//attractors are blue, repellers are red
void draw_repellers(repeller R, AlPolyset *polyset) {
  int base[24];
  int p=0;
  float rx,ry,rz,h,w;
  
  for(int t=0;t<R.size;t++) {
    if(R.power[t]>0) p=2; else p=9;
    h=R.power[t]/2;
    if(h<0) h=-h;
    w=h/4;
    rx = R.x[t];
    ry = R.y[t];
    rz = R.z[t];
            
    if(R.power[t]>0) {

      base[0]=polyset->newVertex(rx-w,ry-w,rz-h);
      base[1]=polyset->newVertex(rx+w,ry-w,rz-h);
      base[2]=polyset->newVertex(rx+w,ry+w,rz-h);
      base[3]=polyset->newVertex(rx-w,ry+w,rz-h);
      
      base[4]=polyset->newVertex(rx-w,ry-w,rz+h);
      base[5]=polyset->newVertex(rx+w,ry-w,rz+h);
      base[6]=polyset->newVertex(rx+w,ry+w,rz+h);
      base[7]=polyset->newVertex(rx-w,ry+w,rz+h);
      
      base[8]=polyset->newVertex(rx-w,ry-h,rz-w);
      base[9]=polyset->newVertex(rx+w,ry-h,rz-w);
      base[10]=polyset->newVertex(rx+w,ry-h,rz+w);
      base[11]=polyset->newVertex(rx-w,ry-h,rz+w);
      
      base[12]=polyset->newVertex(rx-w,ry+h,rz-w);
      base[13]=polyset->newVertex(rx+w,ry+h,rz-w);
      base[14]=polyset->newVertex(rx+w,ry+h,rz+w);
      base[15]=polyset->newVertex(rx-w,ry+h,rz+w);
      
      base[16]=polyset->newVertex(rx-h,ry+w,rz-w);
      base[17]=polyset->newVertex(rx-h,ry-w,rz-w);
      base[18]=polyset->newVertex(rx-h,ry-w,rz+w);
      base[19]=polyset->newVertex(rx-h,ry+w,rz+w);
      
      base[20]=polyset->newVertex(rx+h,ry+w,rz-w);
      base[21]=polyset->newVertex(rx+h,ry-w,rz-w);
      base[22]=polyset->newVertex(rx+h,ry-w,rz+w);
      base[23]=polyset->newVertex(rx+h,ry+w,rz+w);
    
      addQuad( polyset, base[0], base[1], base[2], base[3], p);
      addQuad( polyset, base[4], base[5], base[6], base[7], p);
      addQuad( polyset, base[0], base[1], base[5], base[4], p);
      addQuad( polyset, base[1], base[2], base[6], base[5], p);
      addQuad( polyset, base[2], base[3], base[7], base[6], p);
      addQuad( polyset, base[3], base[0], base[4], base[7], p);

      addQuad( polyset, base[8], base[9], base[10], base[11], p);
      addQuad( polyset, base[12], base[12], base[13], base[14], p);
      addQuad( polyset, base[8], base[9], base[13], base[12], p);
      addQuad( polyset, base[10], base[9], base[13], base[14], p);
      addQuad( polyset, base[11], base[10], base[14], base[15], p);
      addQuad( polyset, base[8], base[11], base[15], base[12], p);

      addQuad( polyset, base[16], base[17], base[18], base[19], p);
      addQuad( polyset, base[20], base[21], base[22], base[23], p);
      addQuad( polyset, base[16], base[17], base[21], base[20], p);
      addQuad( polyset, base[18], base[17], base[21], base[22], p);
      addQuad( polyset, base[19], base[18], base[22], base[23], p);
      addQuad( polyset, base[19], base[16], base[20], base[23], p);
    }
    else {
      rz=rz+(2*w);
      rx=rx+(2*w);
      base[8]=polyset->newVertex(rx-w,ry-h,rz-w);
      base[9]=polyset->newVertex(rx+w,ry-h,rz-w);
      base[10]=polyset->newVertex(rx+w,ry-h,rz+w);
      base[11]=polyset->newVertex(rx-w,ry-h,rz+w);
      
      base[12]=polyset->newVertex(rx-w,ry+h,rz-w);
      base[13]=polyset->newVertex(rx+w,ry+h,rz-w);
      base[14]=polyset->newVertex(rx+w,ry+h,rz+w);
      base[15]=polyset->newVertex(rx-w,ry+h,rz+w);

      rz=rz-(4*w);
      rx=rx-(4*w);
      
      base[16]=polyset->newVertex(rx-h,ry+w,rz-w);
      base[17]=polyset->newVertex(rx-h,ry-w,rz-w);
      base[18]=polyset->newVertex(rx-h,ry-w,rz+w);
      base[19]=polyset->newVertex(rx-h,ry+w,rz+w);
      
      base[20]=polyset->newVertex(rx+h,ry+w,rz-w);
      base[21]=polyset->newVertex(rx+h,ry-w,rz-w);
      base[22]=polyset->newVertex(rx+h,ry-w,rz+w);
      base[23]=polyset->newVertex(rx+h,ry+w,rz+w);
      
      addQuad( polyset, base[8], base[9], base[10], base[11], p);
      addQuad( polyset, base[12], base[12], base[13], base[14], p);
      addQuad( polyset, base[8], base[9], base[13], base[12], p);
      addQuad( polyset, base[10], base[9], base[13], base[14], p);
      addQuad( polyset, base[11], base[10], base[14], base[15], p);
      addQuad( polyset, base[8], base[11], base[15], base[12], p);

      addQuad( polyset, base[16], base[17], base[18], base[19], p);
      addQuad( polyset, base[20], base[21], base[22], base[23], p);
      addQuad( polyset, base[16], base[17], base[21], base[20], p);
      addQuad( polyset, base[18], base[17], base[21], base[22], p);
      addQuad( polyset, base[19], base[18], base[22], base[23], p);
      addQuad( polyset, base[19], base[16], base[20], base[23], p);
    }
  }
}



//***********************************************************
//****** SEGMENT_GROW ********************************
//this is the primary recursive procedure which interprets L-systems,
//fractal trees, and formerly the VARIANT mode.
statusCode segment_grow(segment *a, repeller R, grammar G, AlPolyset *polyset, int level, char input[80], float obx, float oby, float obz, float otx, float oty, float otz, int old_ang, int old_ang2, int old_ang3) 
{

  
  if(TREE) segment_draw(a,polyset,level);
  float tx=otx;  //BX since it begins at fractal start not end
  float ty=oty;
  float tz=otz;
  float bx=obx;
  //	a->bx-(a->tx-a->bx)*SHRINK_FACTOR*LENGTH_RATIO;
  float by=oby;
  //	a->by-(a->ty-a->by)*SHRINK_FACTOR*LENGTH_RATIO;
  float bz=obz;
  //	a->bz-(a->tz-a->bz)*SHRINK_FACTOR*LENGTH_RATIO;
  if(level<=GENERATION_LIMIT) {

    //  AlUniverse::redrawScreen();
    if(TREE<1) {

      
      //----FOLLOWING ORDERS-----
      
      //F] is next, spawn child with [full string] (default 0)
      //F[ is next, spawn child, (enclosed in [])
      //[ is next,  spawn child, ([enclosed in []])
      //+ is next, rotate ang
      //- is next, rotate ang
      //] is next, pop ang
      
      int ang=old_ang;
      int ang2=old_ang2;
      int ang3=old_ang3;
      //ang=0;
      //ang2=0;
      int vnum=0;
      int base[20];
      int was=0;

      printf("\nL:%d <%s>",level,input);
      float oldx=a->bx;
      float oldy=a->by;
      float oldz=a->bz;
      int pos=0;
      while(pos<strlen(input)) {
	for(int line=0;line<G.lines;line++) {
	  if(input[pos]==G.key[line]) {
	    a->bx=tx;
	    a->by=ty;
	    a->bz=tz;
	    float old_length=a->length;
	    float oldx=a->bx;
	    float oldy=a->by;
	    float oldz=a->bz;
	    a->length*=SHRINK_FACTOR*LENGTH_RATIO;
	    //printf("G<%c>",G.key[line]);
	    a->child[a->children]=(segment*)malloc(sizeof(segment));
	    //printf("#");
	    *a->child[a->children]=segment_copy(a);
	    if(surface==1) was=1; else was=0;
	    if((input[pos]=='S'||input[pos]=='T')&&resources>0) {surface=1;resources--;}
	    if(resources>0) segment_grow(a->child[a->children],R,G,polyset,level+1, G.value[line] ,bx,by,bz,tx,ty,tz,ang,ang2,ang3);
	    if(!was) surface=0;
	    a->children++;
	    if(input[pos]!='F') pos++;
	 //	    printf("C<%d %s>",level,input);
	    //	    if(level<=0) {printf("\nDONE!!!\n");return sSuccess;}
	    a->bx=oldx;
	    a->by=oldy;
	    a->bz=oldz;
	    a->length=old_length;
	  }
	}
	switch(input[pos]) {
	case '+': pos++;ang++;break;  //turn
	case '-': pos++;ang--;break;
	case '&': pos++;ang2--;break; //pitch
	case '^': pos++;ang2++;break;
	case '<': pos++;ang3--;break; //roll
	case '>': pos++;ang3++;break;
	case 'n': pos++;ang=0;break;
	case '[': 
	  char* substring=(char*)malloc(strlen(input));
	  substring="*";
	  int num=0;
	  //	  input=&input[1];
	  pos++;
	  while(input[pos]!=']') {
	    substring[num]=input[pos];
	    num++;
	    pos++;
	    //	    input=&input[1];
	  }
	  substring[num]='*';
	  a->bx=tx;
	  a->by=ty;
	  a->bz=tz;
	  //	  printf("g");
	  segment_grow(a,R,G,polyset,level,substring,bx,by,bz,tx,ty,tz,ang,ang2,ang3);
	  //printf("c");
	  a->bx=oldx;
	  a->by=oldy;
	  a->bz=oldz;
	  //	  input=&input[1];
	  pos++;
	  break;
	case ']':
	  pos++;
	  //	  input=&input[1];

	  break;
	case 'F':
	case 'f':
	  //printf("f");
	  //segment_place(a,bx,by,bz,tx,ty,tz,ang,ang2,ang3);
	  //tx=a->tx;
	  //ty=a->ty;
	  //tz=a->tz;

	  //NOTE: IF surface==1, then no attractor/repeller effects will take place
	  int limit=segment_put(a,R,&bx,&by,&bz, &tx,&ty,&tz,ang,ang2,ang3,surface);
	  if(limit) {level=100;}

	  ang2=0;
	  ang3=0;
	  ang=old_ang;

	  
	  if(surface>0&&vnum<FACET_SIZE) {
	    base[vnum] = polyset->newVertex(tx,ty,tz);
	    V[Vnum].x=tx;
	    V[Vnum].y=ty;
	    V[Vnum].z=tz;
	    V[Vnum].generation=level;
	    V[Vnum].vertex=vnum;
	    V[Vnum].polygon=pnum;
	    V[Vnum].fixed=0;
	    Vnum++;
	    vnum++;
	  }
	  if(vnum==FACET_SIZE) {
	    pnum++;
	    if(FILTER3==0) {
	      int polygon_index = polyset->newPolygon();
	      AlPolygon *polygon = polyset->polygon( polygon_index );
	      if( polygon )
		{
		  for(int t=0;t<FACET_SIZE;t++) {
		    polygon->addVertex(base[t]);
		    
		  }
		  polygon->setShaderIndex(1);
		  delete polygon;
		}
	    }
	    vnum=0;
	  }
	  if(input[pos]=='F') segment_draw(a->child[a->children-1],polyset,level+1);
	  //	  input=&input[1];
	  pos++;
	  break;
	default:printf("\n{%s}\n",input);pos++;break;
	}  
      }
      //      printf("D");
    }
    if(TREE) segment_spawn(a,CHILDREN);
    if(TREE) 
      for(int temp=0;temp<a->children;temp++) segment_grow(a->child[temp],R,G,polyset,level+1,INPUT,bx,by,bz,tx,ty,tz,0,0,0);
  } //END GENERATION CHECK
  
  return(sSuccess);
}



statusCode read_input_file(grammar *G, repeller *R) {
  FILE *in;
  in=fopen("/mit/edgsrc/moss/devynmoss/data.txt","rt");
  char c;
  char C[80]={};
  int len=0;
  int data=0;
  int size=0;
  char key='!';
  int current=0;
  G->lines=0;
  for(int q=0;q<GRAMMAR_MAX;q++)
    for(int u=0;u<80;u++)
      G->value[q][u]='\0';
  while(!feof(in)) {
    c=getc(in);
    C[len]=c;
    len++;
    if(len==80) len=0;
    if(c=='=') {
      for(int t=0;t<len-1;t++) {
	//	printf("%c",C[t]);  //i.e. GENERATION
      }
      if(C[2]=='=') {data=16;key=C[1];printf("");}
      if(C[3]=='-'&&C[4]=='X') {data=18;size=C[2]-48;}
      if(C[3]=='-'&&C[4]=='Y') {data=19;size=C[2]-48;}
      if(C[3]=='-'&&C[4]=='Z') {data=20;size=C[2]-48;}
      if(C[3]=='-'&&C[4]=='P') {data=21;size=C[2]-48;}
      if(C[1]=='G'&&C[2]=='E'&&C[3]=='N') {printf("=");data=1;}
      if(C[1]=='F'&&C[2]=='I'&&C[3]=='L') {printf("=");data=2;}
      if(C[1]=='M'&&C[2]=='E'&&C[3]=='S') {printf("=");data=3;}
      if(C[1]=='T'&&C[2]=='R'&&C[3]=='E') {printf("=");data=4;}
      if(C[1]=='C'&&C[2]=='H'&&C[3]=='I') {printf("=");data=5;}
      if(C[1]=='A'&&C[2]=='L'&&C[3]=='P') {printf("=");data=6;}
      if(C[1]=='B'&&C[2]=='E'&&C[3]=='T') {printf("=");data=7;}
      if(C[1]=='G'&&C[2]=='A'&&C[3]=='M') {printf("=");data=8;}
      if(C[1]=='X'&&C[2]=='1') {printf("=");data=9;}
      if(C[1]=='Y'&&C[2]=='1') {printf("=");data=10;}
      if(C[1]=='Z'&&C[2]=='1') {printf("=");data=11;}
      if(C[1]=='X'&&C[2]=='2') {printf("=");data=12;}
      if(C[1]=='Y'&&C[2]=='2') {printf("=");data=13;}
      if(C[1]=='Z'&&C[2]=='2') {printf("=");data=14;}
      if(C[1]=='F'&&C[2]=='I'&&C[7]=='2') {printf("=");data=15;}
      if(C[1]=='L'&&C[2]=='I'&&C[3]=='N') {printf("=");data=17;}
      if(C[1]=='F'&&C[2]=='A'&&C[3]=='C') {printf("=");data=22;}
      if(C[1]=='S'&&C[2]=='H'&&C[3]=='A') {printf("=");data=23;}
      if(C[1]=='T'&&C[2]=='H'&&C[3]=='R') {printf("=");data=24;}
      if(C[1]=='S'&&C[2]=='X') {SNUM=C[3]-48;printf("=");data=25;}
      if(C[1]=='S'&&C[2]=='Y') {printf("=");data=26;}
      if(C[1]=='S'&&C[2]=='Z') {printf("=");data=27;}
      if(C[1]=='F'&&C[2]=='I'&&C[7]=='3') {printf("=");data=28;}
      len=0;
      for(t=0;t<80;t++) C[t]='\0';
    }
    if(c=='.') {
      int sum=0;int factor=1;
      if(data>0) {
	if(data!=16) for(int t=len-2;t>=0;t--) {
	  if(C[t]=='-') sum=-sum; else 
	    sum=sum+(C[t]-48)*factor;
	  factor=factor*10;
	}
	printf("\n\n");
	if(data==16) {
	  C[len-1]='\0';
	  char temp[80]="";
	  strcpy(G->value[current],temp);
	  //	  for(int t=len-1;t<40;t++)
	  //  C[t]='\0';
	  strcat(G->value[current],(const char*)C);
	  G->key[current]=key;
	  printf("%c=>%s\n",G->key[current],G->value[current]);
	  current++;
	}
	for(int t=0;t<80;t++) C[t]='\0';
	switch(data) {
	case 1: GENERATION_LIMIT=sum;break;
	case 2: FILTER=sum;break;
	case 3: MESH=sum;break;
	case 4: TREE=sum;break;
	case 5: CHILDREN=sum;break;
	case 6: ALPHA=((float)sum)/10;break;
	case 7: BETA=((float)sum)/10;break;
	case 8: GAMMA=((float)sum)/10;break;
	case 9: X1=sum/10;break;
	case 10: Y1=sum/10;break;
	case 11: Z1=sum/10;break;
	case 12: X2=sum/10;break;
	case 13: Y2=sum/10;break;
	case 14: Z2=sum/10;break;
	case 15: FILTER2=sum;break;
	case 17: G->lines=sum;break;
	case 18: R->x[size]=sum/10;break;
	case 19: R->y[size]=sum/10;break;
	case 20: R->z[size]=sum/10;break;
	case 21: R->power[size]=sum/10;break;
	case 22: FACET_SIZE=sum;break;
	case 23: shared=sum;break;
	case 24: threshhold=((float)sum)/10;break;
	case 25: SX[SNUM]=((float)sum)/10;break;
	case 26: SY[SNUM]=((float)sum)/10;break;
	case 27: SZ[SNUM]=((float)sum)/10;break;
	case 28: FILTER3=sum;break;
	default: break;
	}
	if(data!=16) printf("%d",sum);
      }
      data=0;len=0;
    }
  }
  G->lines=current;
  R->size=size+1;
  SNUM++;
  printf("\nlines=%d\n",G->lines);
  printf("\n\n\n****************************\n\n");
  for(int t=0;t<G->lines;t++) {
    //printf("%c=%s\n",G->key[t],G->value[t]);

  }
  fclose(in);
  return sSuccess;
}


//initalizes the system and initally calls segment_grow...
statusCode Tree::plant( AlPolyset *polyset ) const
{
  srandom( seed );

  struct grammar G;
  //grammar_new();
  struct repeller R=repeller_new();
  R.x[0]=2;
  R.y[0]=0;
  R.z[0]=4;
  R.power[0]=-2;
  R.x[1]=-2;
  R.y[1]=0;
  R.z[1]=3.5;
  R.power[1]=-2;
  R.size=2;

  SNUM=1;
  SX[0]=0;
  SY[0]=0;
  SZ[0]=0;

  read_input_file(&G, &R);

  struct segment root=segment_new(0.0,0.0,0.0, 0, 0 ,1, TRUNK_RATIO*BASE, BASE, 0,0,0);



  // start off the random number generator
  
  // start off the main area
  make_COLORS(polyset);
  // make_BLACK(polyset);
  //make_RED(polyset);

  //**************STARTING POINT GROWTH**********
  Vnum=0;
  for(int t=0;t<SNUM;t++) 
    segment_grow(&root, R, G, polyset,0,START,  SX[t],SY[t],SZ[t]-0.5,  SX[t]+0.00000010,SY[t]+0.000000001,SZ[t],  0,0,0);
  //***APPLY FILTERS / POST-PROCESSES****
  if(FILTER) filter(polyset);
  if(FILTER2) filter2(polyset);
  if(FILTER3) {
    filter3(polyset);
    FILTER3=0;
    Vnum=0;
    for(t=0;t<SNUM;t++) 
      segment_grow(&root, R, G, polyset,0,START,  SX[t],SY[t],SZ[t]-0.5,  SX[t]+0.00000010,SY[t]+0.000000001,SZ[t],  0,0,0);
    filter3b(polyset);
    FILTER3=1;
  }
  if(MESH) mesh(polyset);
  
  draw_repellers(R,polyset);
  printf("\nDrawing...\n");
  AlUniverse::redrawScreen();
  
  return sSuccess;
  //  return segment_cleanup(root);
}

//***********************************************************
//PLUGIN STUFF.... NOT TOO IMPORTANT...

static void do_tree()
{
	Tree* tree=(Tree*) malloc(sizeof(tree));
	tree = new Tree;
	AlPrintf( kPrompt, "MOSS L-system generating...");
	if( tree->getParams() == sSuccess )
	{
		AlPolyset * polyset=new AlPolyset;

		if( polyset->create() != sSuccess )
		{
			AlPrintf( kPrompt, "Couldn't create polyset");
		}
		else
		{
		  tree->plant( polyset );
		}

	}
}

static AlFunctionHandle h;
static AlMomentaryFunction hFunc; 

extern "C"
int plugin_init( const char *dirName )
{
	char *dirScm;
		
	AlUniverse::initialize( kZUp );

	dirScm = makeAltPath(dirName,"scm");

	AlInvokeSchemeFile( INIT_OPTIONBOX_FILE, dirScm );

	hFunc.create( "pl_FractalTree", do_tree );
	h.create( "Poly Fractal Tree", &hFunc ); 
	h.setAttributeString( "tree" );

	if ( sSuccess != h.setOptionBox( OPTIONBOX_FILE, "tree.options", dirScm ) ) {
		AlPrintf( kPrompt, "Tree plug-in unable to find .scm file for option box" );
		return 1;
	}
	h.setIconPath( makeAltPath( dirName, NULL ) );
	h.appendToMenu( "mp_objtools" );

	AlPrintf( kPrompt, "Fractal Tree installed on Palette 'Object Edit'");

	return 0;
}

extern "C"
int plugin_exit( void )
{
	h.removeFromMenu();
	h.deleteObject();
	hFunc.deleteObject();

	// do nothing
	return 0;
}