/* Program for configuring Bush Robot symmetric pose */

#define CALCDEPTH 10		/* number of tree levels in internal calculation */
#define VRMLDEPTH 6		/* number of tree levels in VRML display */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "cmacs.h"

#define Title	"Bush Robot"

int Base;				    /* Base (branching) factor "B" */
double Scale;				/* Branch Scale factor: to_child/from_parent */
int Depth;					/* Number of levels above root */
double BSize;				/* Actual size of first level branches */

double logB;				/* log(B<ranching>) */
int Leaves, Twigs;

double radixr = 0.0, radixi = 1.4142135;  /* complex cover parameters */
double Radius = 1.414213562;
int FixNodes = 0;			/* Relax whole tree to root */


float *xp, *yp, *zp;		 /* endpoint co-ordinates */
/* A branch is defined by a vector V = [c0,c1,c2 ... cn-1] of
   subtree choices, with 0 <= ci < B.  The endpoint of V
   is located in the [xp, yp, zp] arrays at position 
   P(V) = (B^(n+1)-1)/(B-1) + sum(ci*B^(n-1-i),j=0..n-1)
   Its root is at P([c0,c1,c2 ... cn-2]).
   Given P,     n = log(P*(B-1)+1)/log(B)
   c[n-1] = (P - (B^n-1)/(B-1)) % B;
   */
float *xt[4], *yt[4], *zt[4];   /* vertices of branch bases of triangle skin */
int nodes, Leaves;  /* internal nodes and all nodes, including leaves */
int *bperm[3];  /* subbranch permutation for branch base tesselation */

FILE *ivskinf;   /* Output file pointer */


int NC=0, C[1000] = {0};   /* subtree path */

double saz[100], caz[100], sel[100], cel[100];

char Viewpoint[1000] = "";
   
double frand(void);  /* random double from 0 to 1 */

int Psize(int n);  /* number of points in depth n tree */
int Pend(int n);  /* endpoint index from n digits of subtree path */
int Nbranch(int P);   /* branching level of endpoint index */
int CNmake(int P);  /* set up subtree path from endpoint index */
double brlen(int a); /* length of big branch leading up to a */
int Parent(int P);  /* parent endpoint from endpoint */
int Child(int P);  /* first child of parent P; others are +1, +2 .. */

double xcoord(int depth);
double ycoord(int depth);
void setradix (double rad, double rot);
void radixcoord(int depth,  double *x, double *y);
float *ck(float *ma, char *mess);

void make_skin(void);


/* random double from 0 to 1 */
double frand(void) { return (rand()&0xfff)/((double) 0xfff); }

int Psize(int n)  /* number of vertices in depth n tree */
{  return((pow(Base,n+1)-1)/(Base-1)); }

int Pend(int n)  /* endpoint index from n digits of subtree path */
{	int i, PP;
PP = 0;  for (i=0; i<n; i++) PP = PP*Base + C[i];
return(PP + Psize(n-1));  }

int Nbranch(int P)  /* branching level of endpoint index */
{	int t;
t = 5e-8 + log(P*(Base-1)+1)/logB;
return(t); }

int CNmake(int P)  /* set up subtree path from endpoint index */
{	int n, t;	
NC = Nbranch(P);
t = P - Psize(NC-1);
for (n = NC-1; n>=0; n--)
  {
    C[n] = t % Base;
    t = t / Base;
  }	
}

double brlen(int a)  /* length of big branch leading up to a */
{ return(BSize*pow(Scale,Nbranch(a)-1)); }

int Child(int P)  /* node P's first child. add 0, 1, 2 ... to get all children */
{ int n, a;
  n = Nbranch(P)+1;
  a = (P - Psize(n-2))*Base + Psize(n-1);
  return(a);
  }


int Parent(int Ch)  /* node Ch's parent branch */
{ int n, n1, a, b; 
  n = Nbranch(Ch);  
  a = (Ch - Psize(n-1))/Base + Psize(n-2);
/* n1 = Nbranch(a);
  if (n-n1 != 1) printf ("n-1 = %d  n1 = %d \n",n-1,n1);
  b = Child(a);
  if (Ch-b<0 || Ch-b>=Base)  printf("D %d P %d -> %d -> %d\n",Ch-b,Ch,a,b); */
  return(a);
  }



double xcoord(int depth)
{
  double v;  int i;
  v=0.0;
  for (i = NC-1; i>=depth; i--)
    if ((i&1)==1)  { v = (C[i]+v)/Base; }
  return(v);
}

double ycoord(int depth)
{
  double v;  int i;
  v=0.0;
  for (i = NC-1; i>=depth; i--)
    if ((i&1)==0)  { v = (C[i]+v)/Base; }
  return(v);
}

void setradix (double rad, double rot)
{double mag; /* radix of complex number system for converting
		a number to a position on XY plane, in
		exponential notation  radix = rad*exp(%i*rot).
		In Base 2, rad should be sqrt(2) if twigs are
		to cover evenly.  rot = %pi/2 gives a simple
		rectangular grid.  rot = 3*%pi/4 gives an
		interesting "double dragon" fractal set.
		(Knuth v2 p190).
		In Base 3, rad should be sqrt(3), and rot=%pi/2
		gives a unique fractal covering
		In Base 4 rad should be sqrt(4), and rot=0
		gives a square grid */
Radius = rad;
radixr = cos(rot)/rad;  radixi = sin(rot)/rad;
}


void radixcoord(int depth,  double *x, double *y)
{double vr, vi, tr, ti, digit;  int i;
vr = vi = 0.0;
for (i=NC-1; i>=depth; i--)
  {
    digit = (C[i]*TWOPI)/Base; 
    vr += (Radius-1.0)*sin(digit);
    vi += (Radius-1.0)*cos(digit);
    tr = vr*radixr+vi*radixi;
    ti = vi*radixr-vr*radixi;
    vr = tr;  vi = ti;}
*x = tr;  *y = ti;
return;
}


float *ck(float *ma, char *mess)
{if (ma==NULL) {printf("%s malloc failed!\n",mess); exit(0);}
return(ma);
}


int main(argc, argv)
int argc; char *argv[];
{	char dummy[100];
int i, j, scene, dad;  
double scl, BSizeBest, ErrBest, Err, p;
double dx, dy, dz;   double xc, yc, tx, ty;
double xmin, xmax, ymin, ymax, zmin, zmax;

scene = 5;

printf("Scene %d, Calc depth %d, VRML depth %d\n\n",scene,CALCDEPTH,VRMLDEPTH);

Base = 3;
Scale = 1.0/sqrt(Base);
Depth = CALCDEPTH;
BSize = 1.5;  /* physical scale */

logB = log(Base);
Leaves = Psize(Depth);   /* node count of whole tree */
Twigs = Psize(Depth-1);  /* tree nodes not including final leaves */

xp = ck(malloc(Leaves*sizeof(float)),"x ");  /* vertex coords */
yp = ck(malloc(Leaves*sizeof(float)),"y ");
zp = ck(malloc(Leaves*sizeof(float)),"z ");

/* set up the heights of all the levels of the tree */

for (i=0; i<Leaves; i++)
  {
    p =  Nbranch(i);
    zp[i] = BSize*2.75*(pow(Scale,p) - Scale)/(1.0 - Scale);
  };


/* set up the x/y position of the leaves */

for (i=Twigs; i<Leaves; ++i)
  {
    CNmake(i);  
    setradix(sqrt(3.0),0.5*PI);
    radixcoord(0,&xc,&yc);
    xp[i] = BSize*4*xc;
    yp[i] = BSize*4*yc;
  }

FixNodes=1;

/* Interpolate x,y positions from the leaves towards the root */

for (i=Leaves-Base; i>0; i-=Base)
  {
    dad = Parent(i);
    if (dad>=FixNodes)
      {
	dx = dy = 0;
	for (j = 0; j < Base; j++)
	  { dx += xp[i+j];  dy += yp[i+j]; }
    
	xp[dad] = dx/Base;  yp[dad] = dy/Base;
      };
  }

xmax = xmin = xp[0];
ymax = ymin = yp[0];
zmax = zmin = zp[0];

for (i=0; i<Leaves; i++)
  {
    p =  Nbranch(i);
    if (xp[i]<xmin) xmin=xp[i];  if (xp[i]>xmax) xmax=xp[i];
    if (yp[i]<ymin) ymin=yp[i];  if (yp[i]>ymax) ymax=yp[i];
    if (zp[i]<zmin) zmin=zp[i];  if (zp[i]>zmax) zmax=zp[i];

  };

printf("Range (%.3g %.3g %.3g) ((%.2g %.2g) (%.2g %.2g) (%.2g %.2g))\n",
xmax-xmin, ymax-ymin, zmax-zmin,
xmin, xmax, ymin, ymax, zmin, zmax);

dx = xp[Leaves-1] - xp[Leaves-2];
dy = yp[Leaves-1] - yp[Leaves-2];
dz = zp[Leaves-1] - zp[Leaves-2];
dz = sqrt(dx*dx+dy*dy+dz*dz);

printf("Leaf separation %.3g = %.3g of height\n",dz,dz/(zmax-zmin));

printf("Making skin ... ");
make_skin();
printf("Done\n");

}


void MakeSiSkin() /* make simple inventor skin for bush */
{
char fname[20]="siskin00.iv";
int i,j,l,l1;

fname[7]='0'+(Depth%10);
fname[6]='0'+((Depth/10)%10);

ivskinf = fopen(fname,"w");
fprintf(ivskinf,"#Inventor V2.0 ascii\n\nBaseColor{ rgb .9 .7 .8 }\n\n");

/* make joins for skin, implicit in resulting vertex ordering */
for (i=0; i < nodes; i++)
  {
    for (l=0; l<3; l++)
      {
	l1 = l+1;  while (l1>=3) l1 -= 3;
	j = bperm[l][i];
	fprintf(ivskinf,"Coordinate3 {\n point [ \n");
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[3][i],yt[3][i],zt[3][i]);
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[1][j],yt[1][j],zt[1][j]);
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[l1][i],yt[l1][i],zt[l1][i]);
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[0][j],yt[0][j],zt[0][j]);
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[l][i],yt[l][i],zt[l][i]);
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[2][j],yt[2][j],zt[2][j]);
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[3][i],yt[3][i],zt[3][i]);
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[1][j],yt[1][j],zt[1][j]);
	fprintf(ivskinf," ]\n }\nTriangleStripSet{}\n\n");
      };
  };

/* Top endcap */
fprintf(ivskinf,"Coordinate3 {\n point [ \n");
fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[0][0],yt[0][0],zt[0][0]);
fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[2][0],yt[2][0],zt[2][0]);
fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[1][0],yt[1][0],zt[1][0]);
fprintf(ivskinf," ]\n }\nTriangleStripSet{}\n\n");

/* Bottom endcaps */
for (i=nodes; i < Leaves; i++)
  {
    fprintf(ivskinf,"Coordinate3 {\n point [ \n");
    fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[0][i],yt[0][i],zt[0][i]);
    fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[1][i],yt[1][i],zt[1][i]);
    fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[2][i],yt[2][i],zt[2][i]);
    fprintf(ivskinf," ]\n }\nTriangleStripSet{}\n\n");
  };


fclose(ivskinf);
}

void ivtri1(int t1,  int v1, int t2,  int v2, int t3,  int v3 )  
     /* write out iv triangle with coords {xyz}t[t?,v?] */
{
	fprintf(ivskinf,"Coordinate3 {\n point [ \n");
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[t1][v1],yt[t1][v1],zt[t1][v1]);
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[t2][v2],yt[t2][v2],zt[t2][v2]);
	fprintf (ivskinf,"  %.3g %.3g %.3g , \n", xt[t3][v3],yt[t3][v3],zt[t3][v3]);
	fprintf(ivskinf," ]\n }\nTriangleStripSet{}\n\n");
}

void ivtri(int t1,  int v1, int t2,  int v2, int t3,  int v3 )  
     /* write out iv triangle with coords {xyz}t[t?,v?] */
{	double ux, uy, uz, cx, cy, cz, ax, ay, az, a;

	fprintf(ivskinf,"Coordinate3 { point [ \n");
	fprintf (ivskinf,"  %.3g %.3g %.3g ,\n", xt[t1][v1],yt[t1][v1],zt[t1][v1]);
	fprintf (ivskinf,"  %.3g %.3g %.3g ,\n", xt[t2][v2],yt[t2][v2],zt[t2][v2]);
	fprintf (ivskinf,"  %.3g %.3g %.3g ", xt[t3][v3],yt[t3][v3],zt[t3][v3]);
	fprintf (ivskinf,"] }\n");

	ux = xt[t2][v2]-xt[t1][v1];  	
	uy = yt[t2][v2]-yt[t1][v1];	
	uz = zt[t2][v2]-zt[t1][v1];
	cx = xt[t3][v3]-xt[t2][v2];  	
	cy = yt[t3][v3]-yt[t2][v2];	
	cz = zt[t3][v3]-zt[t2][v2];

        /* a = ux x cx  gives orientation of triangle face  */
	ax = uy*cz-uz*cy;   ay = uz*cx-ux*cz;   az = ux*cy-uy*cx;
	a = 1.0/sqrt(ax*ax+ay*ay+az*az);  ax *= a;  ay *= a;  az *= a;

	fprintf(ivskinf,"Normal { vector [ \n");
	fprintf (ivskinf,"  %.3g %.3g %.3g ,\n", ax, ay, az);
	fprintf (ivskinf,"  %.3g %.3g %.3g ,\n", ax, ay, az);
	fprintf (ivskinf,"  %.3g %.3g %.3g ", ax, ay, az);
	fprintf (ivskinf,"] }\n");
	fprintf (ivskinf,"IndexedFaceSet { coordIndex [ 0, 1, 2, -1 ] }\n\n");

}

void MakeTrSkin() /* make separate triangle inventor skin for bush */
{
char fname[20]="trskin00.iv";
int i,j,l,l1;

fname[7]='0'+(Depth%10);
fname[6]='0'+((Depth/10)%10);

ivskinf = fopen(fname,"w");
fprintf(ivskinf,"#Inventor V2.0 ascii\n\nBaseColor{ rgb .9 .7 .8 }\n\n");

/* make joins for skin, implicit in resulting vertex ordering */
for (i=0; i < nodes; i++)
  {
    for (l=0; l<3; l++)
      {
	l1 = l+1;  while (l1>=3) l1 -= 3;
	j = bperm[l][i];
	ivtri(3,i, 1,j, l1,i);
	ivtri(l1,i, 1,j, 0,j);
	ivtri(l1,i, 0,j, l,i);
	ivtri(l,i, 0,j, 2,j);
	ivtri(l,i, 2,j, 3,i);
	ivtri(3,i, 2,j, 1,j);
      };
  };

/* Top endcap */
ivtri(0,0, 2,0, 1,0);

/* Bottom endcaps */
for (i=nodes; i < Leaves; i++)   ivtri(0,i, 1,i, 2,i);


fclose(ivskinf);
}

void MakeTrSkinRoute() /* make one path down separate triangle inventor skin */
{
char fname[20]="trskro00.iv";
int i,j,l,l1,n,f;

fname[7]='0'+(Depth%10);
fname[6]='0'+((Depth/10)%10);

ivskinf = fopen(fname,"w");
fprintf(ivskinf,"#Inventor V2.0 ascii\n\nBaseColor{ rgb .9 .7 .8 }\n\n");

/* make joins for skin, implicit in resulting vertex ordering */
for (i=0; i < nodes; i++)
  {
    CNmake(i); n=Nbranch(i);
    f = 1;  for (j=0; j<n; j++) f &= C[j] == 0;
    if (f && n>CALCDEPTH-4)
    for (l=0; l<3; l++)
      {
	l1 = l+1;  while (l1>=3) l1 -= 3;
	j = bperm[l][i];
	ivtri(3,i, 1,j, l1,i);
	ivtri(l1,i, 1,j, 0,j);
	ivtri(l1,i, 0,j, l,i);
	ivtri(l,i, 0,j, 2,j);
	ivtri(l,i, 2,j, 3,i);
	ivtri(3,i, 2,j, 1,j);
      };
  };

fclose(ivskinf);
}


void STLtri(int t1,  int v1, int t2,  int v2, int t3,  int v3 )  
     /* write out STL triangle with coords {xyz}t[t?,v?] */
{	double ux, uy, uz, cx, cy, cz, ax, ay, az, a;

	ux = xt[t2][v2]-xt[t1][v1];  	
	uy = yt[t2][v2]-yt[t1][v1];	
	uz = zt[t2][v2]-zt[t1][v1];
	cx = xt[t3][v3]-xt[t2][v2];  	
	cy = yt[t3][v3]-yt[t2][v2];	
	cz = zt[t3][v3]-zt[t2][v2];

        /* a = ux x cx  gives orientation of triangle face  */
	ax = uy*cz-uz*cy;   ay = uz*cx-ux*cz;   az = ux*cy-uy*cx;
	a = 1.0/sqrt(ax*ax+ay*ay+az*az);  ax *= a;  ay *= a;  az *= a;

	fprintf(ivskinf,"  facet normal  %.3f %.3f %.3f \n", ax, ay, az);
	fprintf(ivskinf,"    outer loop\n");
	fprintf (ivskinf,"      vertex %.5f %.5f %.5f\n", xt[t1][v1],yt[t1][v1],zt[t1][v1]);
	fprintf (ivskinf,"      vertex %.5f %.5f %.5f\n", xt[t2][v2],yt[t2][v2],zt[t2][v2]);
	fprintf (ivskinf,"      vertex %.5f %.5f %.5f\n", xt[t3][v3],yt[t3][v3],zt[t3][v3]);
	fprintf (ivskinf,"    endloop\n");
	fprintf (ivskinf,"  endfacet\n");

}

void MakeStlSkin() /* make STL skin for bush sculpture */
{
char fname[20]="bush00.stl";
int i,j,l,l1;

fname[5]='0'+(Depth%10);
fname[4]='0'+((Depth/10)%10);

ivskinf = fopen(fname,"w");
fprintf(ivskinf,"solid Bush%c%c\n",fname[4],fname[5]);

/* make joins for skin, implicit in resulting vertex ordering */
for (i=0; i < nodes; i++)
  {
    for (l=0; l<3; l++)
      {
	l1 = l+1;  while (l1>=3) l1 -= 3;
	j = bperm[l][i];
	STLtri(3,i, 1,j, l1,i);
	STLtri(l1,i, 1,j, 0,j);
	STLtri(l1,i, 0,j, l,i);
	STLtri(l,i, 0,j, 2,j);
	STLtri(l,i, 2,j, 3,i);
	STLtri(3,i, 2,j, 1,j);
      };
  };

/* Top endcap */
STLtri(0,0, 2,0, 1,0);

/* Bottom endcaps */
for (i=nodes; i < Leaves; i++)   STLtri(0,i, 1,i, 2,i);

fprintf(ivskinf,"endsolid\n");

fclose(ivskinf);
}


/* build a triangulated skin for a B=3  bush robot,
   for special case of regular posture */

void make_skin(void)
{
int i, j, k, l, l1, jl, jj, m;  
int kmax, lmax, mmax;
double apex, a, d, u, v, c, dmax, s30, c30, r;
double ux, uy, uz, vx, vy, vz, cx, cy, cz, ax, ay, az;
double Thik, da, dab, dav, err, ca, sa;
double XL,XH,YL,YH,ZL,ZH, Yx, Yy, Yz, Xx, Xy, Xz, Cx, Cy, Cz;
int bflag[3]; /* flag for branches that have been assigned roots */

Thik = 0.5;
apex = 0.4;   /* height of center apex as fraction of branch height */    

s30 = 1.0/2.0;   c30 = sqrt(3.0)/2.0;

Leaves = Psize(Depth);

if (Base != 3) {printf("make_skin called with B = %d\n",Base); exit(0); };

for (i=0; i<4; i++)
  {
    xt[i] = ck(malloc(Leaves*sizeof(float)),"xt ");  /* triangle vertex coords */
    yt[i] = ck(malloc(Leaves*sizeof(float)),"yt ");
    zt[i] = ck(malloc(Leaves*sizeof(float)),"zt ");
  };

for (i=0; i<3; i++)
  if (!(bperm[i] = malloc(Leaves*sizeof(int))))
    {printf("malloc(bperm) bombed\n"); exit(0);};

/* set up vertices xt,yt,zt[0..2][0] for base triangle */

printf("\n\n");
for (i=1; i<=3; i++)
  {
   printf("l%d  %.3g %.3g %.3g\n",
	  i,xp[i]-xp[0],yp[i]-yp[0],zp[i]-zp[0]);
  };

/* normal vector for base triangle plane */
ux = (xp[1]+xp[2]+xp[3])/3.0 - xp[0];
uy = (yp[1]+yp[2]+yp[3])/3.0 - yp[0];
uz = (zp[1]+zp[2]+zp[3])/3.0 - zp[0];
u = 1.0/sqrt(ux*ux+uy*uy+uz*uz);  ux *= u;  uy *= u;  uz *= u;
printf("u %.3g %.3g %.3g\n",ux,uy,uz);

/* one branch of the three */
vx = xp[1]-xp[0]; vy = yp[1]-yp[0]; vz = zp[1]-zp[0];
v = 1.0/sqrt(vx*vx+vy*vy+vz*vz);  vx *= v;  vy *= v;  vz *= v;
printf("v %.3g %.3g %.3g\n",vx,vy,vz);

/* cross product giving base triangle base direction */
cx = uy*vz-uz*vy;   cy = uz*vx-ux*vz;   cz = ux*vy-uy*vx;
c = 1.0/sqrt(cx*cx+cy*cy+cz*cz);  cx *= c;  cy *= c;  cz *= c;
printf("c %.3g %.3g %.3g\n",cx,cy,cz);

/* a = c x u  gives direction of base triangle altitude direction */
ax = uy*cz-uz*cy;   ay = uz*cx-ux*cz;   az = ux*cy-uy*cx;
a = 1.0/sqrt(ax*ax+ay*ay+az*az);  ax *= a;  ay *= a;  az *= a;

printf("a %.3g %.3g %.3g\n",ax,ay,az);
printf("\n");

/* now set up apex vertex, and two base vertices */
r = Thik*brlen(0);    /* r = radius of base triangle */
xt[0][0] =  -r*ax + xp[0];   yt[0][0] =  -r*ay + yp[0];   zt[0][0] =  -r*az + zp[0];
xt[1][0] =  r*(c30*cx+s30*ax) + xp[0];
yt[1][0] =  r*(c30*cy+s30*ay) + yp[0];
zt[1][0] =  r*(c30*cz+s30*az) + zp[0];
xt[2][0] = r*(-c30*cx+s30*ax) + xp[0];
yt[2][0] = r*(-c30*cy+s30*ay) + yp[0];
zt[2][0] = r*(-c30*cz+s30*az) + zp[0];

printf("\n");
for (i=0; i<3; i++)
  {
   printf("t%d  %.3g %.3g %.3g\n",
	  i,xt[i][0]-xp[0],yt[i][0]-yp[0],zt[i][0]-zp[0]);
  };

nodes = Psize(Depth-1);       printf(" %d nodes\n",nodes);

mmax = -1;

for (i=0; i < nodes; i++)
  {
    j = Child(i);
    /* set up tesselation of base triangle i, including apex point
       xt,yt,zt[3][i], branch permutation bperm[0..2], the 
       vertices for three branch triangles xt,yt,jt[0..2][j..j+2] */

    /* set up apex point for vertex i */
    xt[3][i] = (1.0-apex)*xp[i] + apex*(xp[j]+xp[j+1]+xp[j+2])/3.0;
    yt[3][i] = (1.0-apex)*yp[i] + apex*(yp[j]+yp[j+1]+yp[j+2])/3.0;
    zt[3][i] = (1.0-apex)*zp[i] + apex*(zp[j]+zp[j+1]+zp[j+2])/3.0;


    m = Nbranch(i);
    l = 
      m==0?0:
      m==1?1:
      m==2?0:
      m==3?1:
      m==4?1:
      m==5?1:
      m==6?0:
      m==7?1:
      m==8?1:
      m==9?1:
      m==10?0:  /* order 11 self-intersects */
      1;
        
    if (m>mmax) printf("Level %d, offset %d\n",mmax=m,l);
    bperm[0][i] = ((l+i+1)%3) +j;  
    bperm[1][i] = ((l+i+2)%3) +j;   
    bperm[2][i] = ((l+i+0)%3) +j;

    /* choose vertices for three equilateral branch triangles xt,yt,jt[0..2][j..j+2] */

    for (l = 0; l < 3; l++)  /* for each branch */
      { /* this branch goes to xp,yp,zp[bperm[l][i]] */

	jl = bperm[l][i];

	/* calculate cross product of branch direction with base triangle altitude */

	/* branch direction */
        if (jl<nodes)
	  {
	    jj = Child(jl);
	    ux = (xp[jj]+xp[jj+1]+xp[jj+2])/3.0 - xp[jl];
	    uy = (yp[jj]+yp[jj+1]+yp[jj+2])/3.0 - yp[jl];
	    uz = (zp[jj]+zp[jj+1]+zp[jj+2])/3.0 - zp[jl];
	  }
	else
	  {
	    /* ux = xp[jl]-xt[3][i];
	    uy = yp[jl]-yt[3][i];
	    uz = zp[jl]-zt[3][i];  */
	    ux =0; uy=0; uz = -1;
	  };
        u = 1.0/sqrt(ux*ux+uy*uy+uz*uz);  ux *= u;  uy *= u;  uz *= u;

	l1 = l+1;  while (l1>=3) l1 -= 3;

	/* base triangle altitude */
	vx = (xt[l][i]+xt[l1][i])/2.0 - xt[3][i];
	vy = (yt[l][i]+yt[l1][i])/2.0 - yt[3][i];
	vz = (zt[l][i]+zt[l1][i])/2.0 - zt[3][i];
/*
   vx = xt[l1][i] - xt[3][i];
   vy = yt[l1][i] - yt[3][i];
   vz = zt[l1][i] - zt[3][i]; 
*/
	v = 1.0/sqrt(vx*vx+vy*vy+vz*vz);  vx *= v;  vy *= v;  vz *= v;


	/* cross product in horizontal direction of new triangle */
	cx = uy*vz-uz*vy;   cy = uz*vx-ux*vz;   cz = ux*vy-uy*vx;
	c = 1.0/sqrt(cx*cx+cy*cy+cz*cz);  cx *= c;  cy *= c;  cz *= c;

	/* a = c x u  gives vertical direction of new triangle */
	ax = uy*cz-uz*cy;   ay = uz*cx-ux*cz;   az = ux*cy-uy*cx;
	a = 1.0/sqrt(ax*ax+ay*ay+az*az);  ax *= a;  ay *= a;  az *= a;

	/* now set up apex vertex, and two base vertices  of new triangle */

	r = Thik*brlen(jl);    /* r = radius of new triangle */
	xt[0][jl] =  -r*ax + xp[jl];
	yt[0][jl] =  -r*ay + yp[jl];
	zt[0][jl] =  -r*az + zp[jl];
	xt[1][jl] =  r*(c30*cx+s30*ax)  + xp[jl];
	yt[1][jl] =  r*(c30*cy+s30*ay)  + yp[jl];
	zt[1][jl] =  r*(c30*cz+s30*az)  + zp[jl];
	xt[2][jl] =  r*(-c30*cx+s30*ax) + xp[jl];
	yt[2][jl] =  r*(-c30*cy+s30*ay) + yp[jl];
	zt[2][jl] =  r*(-c30*cz+s30*az) + zp[jl];

      };

 };

/* Relieve twist in base joints */

for (i=0; i < nodes; i++)
  {
    /* get framework of base triangle for node i */
    /* centroid */
    Cx = (xt[0][i]+xt[1][i]+xt[2][i])/3.0;
    Cy = (yt[0][i]+yt[1][i]+yt[2][i])/3.0;
    Cz = (zt[0][i]+zt[1][i]+zt[2][i])/3.0;
    /* orthogonal vectors Y and X */
    Yx = xt[0][i] - Cx;
    Yy = yt[0][i] - Cy;
    Yz = zt[0][i] - Cz;
    Xx = (xt[2][i]-xt[1][i])/sqrt(3.0);
    Xy = (yt[2][i]-yt[1][i])/sqrt(3.0);
    Xz = (zt[2][i]-zt[1][i])/sqrt(3.0);    

    dab = 0.0;
    {
      dav = 1e9;
      for (da = -PI/2; da < PI/2;  da += .01)
      {
	err = 0.0;
	for (l=0; l<3; l++)
	{
	  a = (150.0 + l*120.0)*PI/180.0 + da;
	  /* l1 = l+1;  while (l1>=3) l1 -= 3; */
	  j = bperm[l][i];
	  ca = cos(a);  sa = sin(a);
	  /* dot product of triangle side with branch direction */
	  err += (Xx*ca + Yx*sa)*(xp[j]-xp[i]) +
	    (Xy*ca + Yy*sa)*(yp[j]-yp[i]) +
	    (Xz*ca + Yz*sa)*(zp[j]-zp[i]);
	};
	if (err < dav) {dav = err; dab = da;};
      };
    };

    /* printf("da %.3g ",dab*180.0/PI); */  dab = 45*PI/180;

    for (l=0; l<3; l++)
      {
	a = (90.0 + l*120.0)*PI/180.0 + dab;
	ca = cos(a);  sa = sin(a);
	xt[l][i] = Xx*ca + Yx*sa + Cx;
	yt[l][i] = Xy*ca + Yy*sa + Cy;
	zt[l][i] = Xz*ca + Yz*sa + Cz;
      };
  };

MakeTrSkinRoute();   /* make separate triangle inventor skin for bush */
if (0) MakeSiSkin();  /* make simple inventor skin for bush */
if (1) MakeTrSkin();   /* make separate triangle inventor skin for bush */
if (1) MakeStlSkin();   /* make STL skin for bush sculpture */


XL = 1e9;  XH = -1e9;   YL = 1e9;  YH = -1e9;   ZL = 1e9;  ZH = -1e9;
for (i=0; i<Leaves; i++)  for (l=0; l<3; l++)
  {
    if (xt[l][i]<XL) XL = xt[l][i];
    if (xt[l][i]>XH) XH = xt[l][i];

    if (yt[l][i]<YL) YL = yt[l][i];
    if (yt[l][i]>YH) YH = yt[l][i];

    if (zt[l][i]<ZL) ZL = zt[l][i];
    if (zt[l][i]>ZH) ZH = zt[l][i];
  };

}

/*
Inventor V2.0 ascii

BaseColor{ rgb .3 .7 .5 }

DEF SPRIG-0b
  Separator {
    Coordinate3 {
        point [   0.075    0.0    0.0,     # 0
                  0.0375   1.0    0.065,   # 4
                 -0.0375   0.0    0.065,   # 1
                 -0.075    1.0    0.0,     # 3
                 -0.0375   0.0   -0.065,   # 2
                  0.0375   1.0   -0.065,   # 5 
                  0.075    0.0    0.0,     # 0
                  0.0375   1.0    0.065,   # 4
               ]
    }
    TriangleStripSet {}
  }

  */

/*
     Coordinate3 {
         point [
             -0.8615 -0.3333 -0.4714,
             0 1 0,
             0.8165 -0.3333 -0.4714,
             0 -0.3333 0.9428,
         ]
     }
     Normal {
         vector [
             -0.8615 -0.3333 -0.4714,
             0 1 0,
             0.8165 -0.3333 -0.4714,
             0 -0.3333 0.9428,
         ]
     }
     IndexedFaceSet {
         coordIndex [
             0, 1, 2, -1,
             0, 1, 3, -1,
             2, 1, 3, -1,
             0, 3, 2, -1,
         ]
     }
*/

/*

solid name

     . . .

  facet normal  0.19 -0.92 -0.34
    outer loop
      vertex 8.24950 3.98440 0.25900
      vertex 7.94800 3.93060 0.23940
      vertex 7.94800 4.03320 -0.04250
    endloop
  endfacet

     . . .

endsolid

*/