// Shave and a Haircut
// (c) 2019 Epic Games
// US Patent 6720962

void put_staticguideRH(MINIKERNAL *m, STATICGUIDE *g, int hairnum);
void fetch_staticguideRH(MINIKERNAL *m, STATICGUIDE *g, int hairnum);
float randcl(long seed);
void color_a_hairK(MINIKERNAL *m,STATICROOTPT *pt,STATICGUIDE *h,int seed);


#ifndef GPU_COMPILE
void get_rootptK(STATICROOTPT *pt,MINIKERNAL *m,int hairID)
{
int x;
for (x=0;x<SHAVE_NUM_PARAMS;x++) 
pt->baked_tex[x]=1.0f; // empty channel = 1.0
for (x=0;x<m->hair_roots[hairID].total_tex;x++) // unroll the linked list into the proper channel directly from the ROOTPT version
  pt->baked_tex[m->hair_roots[hairID].tex_link[x]]=m->hair_roots[hairID].baked_tex[x];
for (x=0;x<SHAVE_NUM_PARAMS;x++)
{
	if (x!=CUTLENGTH) // cutlength doesn't have a slider
   pt->baked_tex[x]*=m->slider_val[x];
}

pt->baked_tex[TIPRAD]*= 4.0f * m->RESTBOUND / 1500.0f;
pt->baked_tex[BASERAD]*= 4.0f * m->RESTBOUND / 1500.0f;
pt->baked_tex[ROOT_COLOR_R]/=255.0f;
pt->baked_tex[ROOT_COLOR_G]/=255.0f;
pt->baked_tex[ROOT_COLOR_B]/=255.0f;
pt->baked_tex[TIP_COLOR_R]/=255.0f;
pt->baked_tex[TIP_COLOR_G]/=255.0f;
pt->baked_tex[TIP_COLOR_B]/=255.0f;


pt->bary=m->hair_roots[hairID].bary;
pt->closest_clump=m->hair_roots[hairID].closest_clump;
pt->distance=m->hair_roots[hairID].distance;
pt->flip=m->hair_roots[hairID].flip;
pt->id=m->hair_roots[hairID].id;
pt->killme=m->hair_roots[hairID].killme;
pt->pntid1=m->hair_roots[hairID].pntid1;
pt->pntid2=m->hair_roots[hairID].pntid2;
pt->pntid3=m->hair_roots[hairID].pntid3;
pt->resthv=m->hair_roots[hairID].resthv;
pt->restlength=m->hair_roots[hairID].restlength;
}

void fetch_staticguide(MINIKERNAL *m, STATICGUIDE *g, int guideID)
{
int x;
g->handle=m->guides[guideID].handle;
g->normal=m->guides[guideID].normal;
g->rest_handle=m->guides[guideID].rest_handle;
g->rest_normal=m->guides[guideID].rest_normal;
 g->uu=m->guides[guideID].uu;
 g->vv=m->guides[guideID].vv;
 g->smoothing_group=m->guides[guideID].smoothing_group;
 g->splitmerge=m->guides[guideID].splitmerge;
for (x=0;x<m->guide_segs;x++)
 {
 g->resthv[x]=m->guides[guideID].resthv[x];
 g->hv[x]=m->guides[guideID].hv[x];
 g->velocity[x]=m->guides[guideID].velocity[x];
 g->noisev[x]=m->guides[guideID].noisev[x];
 }
}

void fetch_staticguideRH(MINIKERNAL *m, STATICGUIDE *g, int guideID)
{
int x;
g->rootcolor=m->return_hair[guideID].rootcolor;
g->tipcolor=m->return_hair[guideID].tipcolor;
g->handle=m->return_hair[guideID].handle;
g->normal=m->return_hair[guideID].normal;
g->rest_handle=m->return_hair[guideID].rest_handle;
g->rest_normal=m->return_hair[guideID].rest_normal;
 g->uu=m->return_hair[guideID].uu;
 g->vv=m->return_hair[guideID].vv;
 g->smoothing_group=m->return_hair[guideID].smoothing_group;
 g->splitmerge=m->return_hair[guideID].splitmerge;
for (x=0;x<m->guide_segs;x++)
 {
 g->resthv[x]=m->return_hair[guideID].resthv[x];
 g->hv[x]=m->return_hair[guideID].hv[x];
 g->velocity[x]=m->return_hair[guideID].velocity[x];
 g->noisev[x]=m->return_hair[guideID].noisev[x];
 }
}
void put_staticguideRH(MINIKERNAL *m, STATICGUIDE *g, int guideID)
{
int x;
m->return_hair[guideID].rootcolor=g->rootcolor;
m->return_hair[guideID].tipcolor=g->tipcolor;
m->return_hair[guideID].normal=g->normal;
m->return_hair[guideID].rest_handle=g->rest_handle;
m->return_hair[guideID].rest_normal=g->rest_normal;
 m->return_hair[guideID].uu=g->uu;
 m->return_hair[guideID].vv=g->vv;
 m->return_hair[guideID].smoothing_group=g->smoothing_group;
for (x=0;x<m->guide_segs;x++)
 {
 m->return_hair[guideID].resthv[x]=g->resthv[x];
 m->return_hair[guideID].hv[x]=g->hv[x];
 m->return_hair[guideID].velocity[x]=g->velocity[x];
 m->return_hair[guideID].noisev[x]=g->noisev[x];
 }
}


#endif

float bary_interpf(float va, float vb, float vc, float a, float b, float c);

static VERT VcrossK( VERT a, VERT b );
static double   det4x4K( Matrix m );
static double   det2x2K( double a, double b, double c, double d );
static double   det3x3K( double, double, double, double, double, double, double, double, double );
static void     adjointK( Matrix in, Matrix out );
static void     inverseK( Matrix in, Matrix out );
static int FroundK( float f );
float VdistanceK(VERT a, VERT b);
void VnormalizeK( VERT *ret );
float getnoiseK( MINIKERNAL *m, float u, float v, float z );
void gen_hairK (MINIKERNAL *m, STATICROOTPT *inpt, STATICGUIDE *ret);
void matinitK( Matrix mat);

void displace_scaleK( STATICGUIDE * h, float rs, MINIKERNAL *m);

VERT vxm3x3K( Matrix m, VERT v	);
void rotxK( Matrix m,float th);
void rotyK( Matrix m,float th);
void rotzK( Matrix m,float th);
void mMatrix2K( Matrix m, VERT xx, VERT yy, VERT zz, VERT b );
void displace_kinkyK( STATICGUIDE * h, MINIKERNAL *m,float *texture ,float restlength, long seed);
float getnoise_octavesK(MINIKERNAL *m,float x, float y, float z, int octaves);


void make_basematrixK(STATICGUIDE *h,BASEMATRIXK *bm,MINIKERNAL *m);

void twist_hairK (STATICGUIDE *clumpCenter, STATICGUIDE *h,MINIKERNAL *m, float *texture, long seed);
void displace_randscaleK( STATICGUIDE * h, MINIKERNAL *m,int seed,float rss );

static void    displace_frizzK( STATICGUIDE * h, MINIKERNAL * m, STATICROOTPT *pp, float restlength, int seed);

static VERT vxmK( Matrix m, VERT v);
float randcl(long seed);

void bary_interpv(VERT *ret,VERT *va, VERT *vb, VERT *vc, float a, float b, float c);

static VERT vxmK( Matrix m,					/* input matrix (in) */
	 VERT v						/* input vector (in) */
	 )
{
	VERT            r;

	r.x = v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0];
	r.y = v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1];
	r.z = v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2];
	return ( r );
}


// these funcs are used by both cl and c since they're just ansii
// and use no old style types

static VERT VcrossK( VERT a, VERT b )
{
	VERT            c;

	c.x = ( a.y * b.z ) - ( a.z * b.y );
	c.y = ( a.z * b.x ) - ( a.x * b.z );
	c.z = ( a.x * b.y ) - ( a.y * b.x );
	return ( c );
}
void matinitK( Matrix mat				/* matrix to initialize (in-out) */
	 )
{
	register int    i, j;

	for( i = 0; i < 4; i++ )
		for( j = 0; j < 4; j++ )
			mat[i][j] = ( i == j ? 1.0 : 0.0 );
}


/*
Matrix Inversion
by Richard Carling
from "Graphics Gems", Academic Press, 1990
*/
#define SMALL_NUMBER	1.e-8
/* 
 *   inverse( original_matrix, inverse_matrix )
 * 
 *    calculate the inverse of a 4x4 matrix
 *
 *     -1     
 *     A  = ___1__ adjoint A
 *         det A
 */
//#include "joemath.h" 
static void inverseK( Matrix in, Matrix out )
{
	int             i, j;
	double          det;

	/* calculate the adjoint matrix */
	adjointK( in, out );
	/*  calculate the 4x4 determinant
	 *  if the determinant is zero, 
	 *  then the inverse matrix is not unique.
	 */
	det = det4x4K( in );
	if( fabs( det ) < SMALL_NUMBER )
	{
//        printf("Non-singular matrix, no inverse!\n");
		/*exit(1); */
		matinitK( out );
	}
	/* scale the adjoint matrix to get the inverse */
	for( i = 0; i < 4; i++ )
		for( j = 0; j < 4; j++ )
			out[i][j] = out[i][j] / det;
}

/* 
 *   adjoint( original_matrix, inverse_matrix )
 * 
 *     calculate the adjoint of a 4x4 matrix
 *
 *      Let  a   denote the minor determinant of matrix A obtained by
 *           ij
 *
 *      deleting the ith row and jth column from A.
 *
 *                    i+j
 *     Let  b   = (-1)    a
 *          ij            ji
 *
 *    The matrix B = (b  ) is the adjoint of A
 *                     ij
 */
void
adjointK( Matrix in, Matrix out )
//Matrix in; Matrix out;
{
	double          a1, a2, a3, a4, b1, b2, b3, b4;
	double          c1, c2, c3, c4, d1, d2, d3, d4;

	/* assign to individual variable names to aid  */
	/* selecting correct values  */
	a1 = in[0][0];
	b1 = in[0][1];
	c1 = in[0][2];
	d1 = in[0][3];
	a2 = in[1][0];
	b2 = in[1][1];
	c2 = in[1][2];
	d2 = in[1][3];
	a3 = in[2][0];
	b3 = in[2][1];
	c3 = in[2][2];
	d3 = in[2][3];
	a4 = in[3][0];
	b4 = in[3][1];
	c4 = in[3][2];
	d4 = in[3][3];

	/* row column labeling reversed since we transpose rows & columns */
	out[0][0] = det3x3K( b2, b3, b4, c2, c3, c4, d2, d3, d4 );
	out[1][0] = -det3x3K( a2, a3, a4, c2, c3, c4, d2, d3, d4 );
	out[2][0] = det3x3K( a2, a3, a4, b2, b3, b4, d2, d3, d4 );
	out[3][0] = -det3x3K( a2, a3, a4, b2, b3, b4, c2, c3, c4 );

	out[0][1] = -det3x3K( b1, b3, b4, c1, c3, c4, d1, d3, d4 );
	out[1][1] = det3x3K( a1, a3, a4, c1, c3, c4, d1, d3, d4 );
	out[2][1] = -det3x3K( a1, a3, a4, b1, b3, b4, d1, d3, d4 );
	out[3][1] = det3x3K( a1, a3, a4, b1, b3, b4, c1, c3, c4 );

	out[0][2] = det3x3K( b1, b2, b4, c1, c2, c4, d1, d2, d4 );
	out[1][2] = -det3x3K( a1, a2, a4, c1, c2, c4, d1, d2, d4 );
	out[2][2] = det3x3K( a1, a2, a4, b1, b2, b4, d1, d2, d4 );
	out[3][2] = -det3x3K( a1, a2, a4, b1, b2, b4, c1, c2, c4 );

	out[0][3] = -det3x3K( b1, b2, b3, c1, c2, c3, d1, d2, d3 );
	out[1][3] = det3x3K( a1, a2, a3, c1, c2, c3, d1, d2, d3 );
	out[2][3] = -det3x3K( a1, a2, a3, b1, b2, b3, d1, d2, d3 );
	out[3][3] = det3x3K( a1, a2, a3, b1, b2, b3, c1, c2, c3 );
}

/*
 * double = det4x4( Matrix )
 * 
 * calculate the determinant of a 4x4 matrix.
 */
double
det4x4K( Matrix m )
{
	double          ans;
	double          a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4;

	/* assign to individual variable names to aid selecting */
	/*  correct elements */
	a1 = m[0][0];
	b1 = m[0][1];
	c1 = m[0][2];
	d1 = m[0][3];
	a2 = m[1][0];
	b2 = m[1][1];
	c2 = m[1][2];
	d2 = m[1][3];
	a3 = m[2][0];
	b3 = m[2][1];
	c3 = m[2][2];
	d3 = m[2][3];
	a4 = m[3][0];
	b4 = m[3][1];
	c4 = m[3][2];
	d4 = m[3][3];
	ans = a1 * det3x3K( b2, b3, b4, c2, c3, c4, d2, d3, d4 ) - b1 * det3x3K( a2, a3, a4, c2, c3, c4, d2, d3, d4 ) + c1 * det3x3K( a2, a3, a4, b2, b3, b4, d2, d3, d4 ) - d1 * det3x3K( a2, a3, a4, b2, b3, b4, c2, c3, c4 );
	return ans;
}

/*
 * double = det3x3(  a1, a2, a3, b1, b2, b3, c1, c2, c3 )
 * 
 * calculate the determinant of a 3x3 matrix
 * in the form
 *
 *     | a1,  b1,  c1 |
 *     | a2,  b2,  c2 |
 *     | a3,  b3,  c3 |
 */
double
det3x3K( double a1, double a2, double a3, double b1, double b2, double b3, double c1, double c2, double c3 )
{
	double          ans;

//    double det2x2();
	ans = a1 * det2x2K( b2, b3, c2, c3 ) - b1 * det2x2K( a2, a3, c2, c3 ) + c1 * det2x2K( a2, a3, b2, b3 );
	return ans;
}

/*
 * double = det2x2( double a, double b, double c, double d )
 * 
 * calculate the determinant of a 2x2 matrix.
 */
double
det2x2K( double a, double b, double c, double d )
{
	double          ans;

	ans = a * d - b * c;
	return ans;
}






static int FroundK( float f )
{
	int             ret;

	if( f > 0 )
		ret = ( int ) ( f + 0.5 );	// Posotive number
	else
		ret = ( int ) ( f - 0.5 );	// Negative number
	return ( ret );
}





float VdistanceK(VERT a, VERT b) {
	VERT d;
	d.x = a.x - b.x;
	d.y = a.y - b.y;
	d.z = a.z - b.z;
	return sqrtf(d.x*d.x + d.y*d.y + d.z*d.z);
}
void VnormalizeK( VERT *v )
{
	VERT            ret;
	float           len;

	len = ( float ) sqrt( ( double ) ( v->x * v->x + v->y * v->y + v->z * v->z ) );
	if( len != 0.0f )
	{
		v->x /= len;
		v->y /= len;
		v->z /= len;
	}
	if( len == 0.0 )
	{
		v->x = .0000001f;
		v->y = 0;
		v->x = 0;
	}
}

float getnoiseK( MINIKERNAL *m, float u, float v, float z )
{
	int             ui, vi, zi, ui2, vi2, zi2;
	float           b1, a, b, c, a1, b2, a2, ur, vr, zr, iur, ivr, izr;

	if( u < 0 )
		u = -u;
	if( v < 0 )
		v = -v;
	if( z < 0 )
		z = -z;

	ui = ( int ) floor( u );
	vi = ( int ) floor( v );
	zi = ( int ) floor( z );
	zr = z - zi;
	ur = u - ui;
	vr = v - vi;
	iur = 1.0f-ur;
	ivr = 1.0f-vr;
	izr = 1.0f-zr;

	ui = ui % 17;
	vi = vi % 17;
	zi = zi % 17;
	ui2 = ui + 1;
	vi2 = vi + 1;
	zi2 = zi + 1;
	ui2 = ui2 % 17;
	vi2 = vi2 % 17;
	zi2 = zi2 % 17;
	a = m->noise[ui][vi][zi];
	b = m->noise[ui2][vi][zi];
	a1 = a * iur + b * ur;
	a = m->noise[ui][vi2][zi];
	b = m->noise[ui2][vi2][zi];
	a2 = a * iur + b * ur;
	b1 = a1 * ivr + a2 * vr;

	a = m->noise[ui][vi][zi2];
	b = m->noise[ui2][vi][zi2];
	a1 = a * iur + b * ur;
	a = m->noise[ui][vi2][zi2];
	b = m->noise[ui2][vi2][zi2];
	a2 = a * iur + b * ur;
	b2 = a1 * ivr + a2 * vr;
	c = b1 * izr + b2 * zr;
	if( c > 1 )
		c = 1;
	if( c < 0 )
		c = 0;
	return ( c );
}




float getnoise_octavesK(MINIKERNAL *m,float x, float y, float z, int octaves)
{
int xx;
float accum=0.0f;
float amp=1.0f;
	x/=2.0f; y/=2.0f ; z/=2.0f;
for (xx=0;xx<octaves;xx++)
	{
	accum+=(getnoiseK(m,x,y,z)-.5)*amp;
	x*=2.0f; y*=2.0f ; z*=2.0f;
	amp/=2.0f;
	}

return(accum);

}

void displace_kinkyK( STATICGUIDE * h, MINIKERNAL *m,float *texture,float restlength,long seed )
{
	int             a;
	int             xx;
	float           ss1, ss2;
	float           freqqx, freqqy, freqqz;
	freqqx = texture[KINKFREQX] * .01;
	freqqy = texture[KINKFREQY] * .01;
	freqqz = texture[KINKFREQZ] * .01;

restlength*=m->guide_segs; 
restlength/=15.0f;

	ss1 = texture[KINK] * 3.0f;
	ss2 = texture[KINKROOT] * 3.0f;

	if( m->interpolation_type == SPLINE_GROWTH )
		freqqx *= 56.0f;
	if( m->interpolation_type == SPLINE_GROWTH )
		freqqy *= 56.0f;
	if( m->interpolation_type == SPLINE_GROWTH )
		freqqz *= 56.0f;


	if( ( ss1 > 0.0f ) || ( ss2 > 0.0f ) )
	{
		for( xx = 1; xx < m->guide_segs; xx ++ )
			if( xx != 0 )
			{
				float           sc1;
				float           sc2;
				float           amp = 1.0f;
				float           freqx, freqy, freqz;
				float           hr;
				VERT            dv;

				hr = ( float ) 9999999.0f;
				dv.x = 0;
				dv.y = 0;
				dv.z = 0;

				hr = restlength;	//yikes!!
				sc1 = ( hr ) * ss1;
				sc1 *= 2.0f;
				sc1 *= ( xx / (float)m->guide_segs );
				sc1 /= 15.0f; // this is really not right , but that's how it was
//sc1/=9.0f;

				sc2 = ( hr ) * ss2;	// root
				sc2 *= 2.0f;
				sc2 *= ( ( (float)m->guide_segs - ( float ) xx ) / (float)m->guide_segs );
//sc2/=9.0f;

				freqx = ( 1.0 / restlength );
				freqy = ( 1.0 / restlength );
				freqz = ( 1.0 / restlength );

				freqx *= freqqx;
				freqy *= freqqy;
				freqz *= freqqz;
				if( ( sc1 != 0.0f ) || ( sc2 != 0.0f ) )
				{

					dv.x += ( float ) ( getnoiseK(m, h->noisev[xx].x * freqx, h->noisev[xx].y * freqy, h->noisev[xx].z * freqz ) - ( float ) .5f ) * amp;
					dv.y += ( float ) ( getnoiseK(m, h->noisev[xx].x * freqx + ( float ) 72.3, h->noisev[xx].y * freqy + ( float ) 16.08, h->noisev[xx].z * freqz + ( float ) 3.0f ) - ( float ) .5f ) * amp;
					dv.z += ( float ) ( getnoiseK(m, h->noisev[xx].x * freqx + ( float ) 15.0f, h->noisev[xx].y * freqy + ( float ) 19.99f, h->noisev[xx].z * freqz + ( float ) 12.0f ) - ( float ) .5f ) * amp;
				}
				amp /= ( float ) 4.0f;
				freqx *= ( float ) 4.0f;
				freqy *= ( float ) 4.0f;
				freqz *= ( float ) 4.0f;
				if( ( sc1 != 0.0f ) || ( sc2 != 0.0f ) )
					{
						dv.x += ( float ) ( getnoiseK(m, h->noisev[xx].x * freqx + ( float ) 72.3, h->noisev[xx].y * freqy + ( float ) 16.08, h->noisev[xx].z * freqz + ( float ) 3.0f ) - ( float ) .5f ) * amp;
						dv.y += ( float ) ( getnoiseK(m, h->noisev[xx].x * freqx + ( float ) 72.3, h->noisev[xx].y * freqy + ( float ) 16.08, h->noisev[xx].z * freqz + ( float ) 3.0f ) - ( float ) .5f ) * amp;
						dv.z += ( float ) ( getnoiseK(m, h->noisev[xx].x * freqx + ( float ) 15.0f, h->noisev[xx].y * freqy + ( float ) 19.99, h->noisev[xx].z * freqz + ( float ) 12.0f ) - ( float ) .5f ) * amp;
						amp /= 2;
					}
				if( sc1 != 0.0f )
				{
					h->hv[xx].x += dv.x * sc1;
					h->hv[xx].y += dv.y * sc1;
					h->hv[xx].z += dv.z * sc1;
				}
				if( sc2 != 0.0f )
				{
					h->hv[xx].x += dv.x * sc2;
					h->hv[xx].y += dv.y * sc2;
					h->hv[xx].z += dv.z * sc2;
				}

			}
	}
}
void mMatrix2K( Matrix m, VERT xx, VERT yy, VERT zz, VERT b )
{
	int x,y;
	for (x=0;x<4;x++)
		for (y=0;y<4;y++)
			m[x][y]=0.0f;

	m[0][0] = xx.x;
	m[0][1] = xx.y;
	m[0][2] = xx.z;
	m[1][0] = yy.x;
	m[1][1] = yy.y;
	m[1][2] = yy.z;
	m[2][0] = zz.x;
	m[2][1] = zz.y;
	m[2][2] = zz.z;
	m[0][3] = 1.0f;
	m[1][3] = 1.0f;
	m[2][3] = 1.0f;
	m[3][3] = 1.0f;
	m[3][0] = 0.0f;
	m[3][1] = 0.0f;
	m[3][2] = 0.0f;

}
void rotzK( Matrix m,					/* resulting rotation matrix (out) */
	  float th					/* rotation vector (in) */
	 )
{
	m[0][0] = cos( th );
	m[1][0] = -sin( th );
	m[2][0] = 0.0f;
	m[3][0] = 0.0f;
	m[0][1] = sin( th );
	m[1][1] = cos( th );
	m[2][1] = 0.0f;
	m[3][1] = 0.0f;
	m[0][2] = 0.0f;
	m[1][2] = 0.0f;
	m[2][2] = 1.0f;
	m[3][2] = 0.0f;
	m[0][3] = 0.0f;
	m[1][3] = 0.0f;
	m[2][3] = 0.0f;
	m[3][3] = 1.0f;

}
void rotyK( Matrix m,					/* resulting rotation matrix (out) */
	  float th					/* rotation vector (in) */
	 )
{
	m[0][0] = cos( th );
	m[1][0] = 0.0f;
	m[2][0] = sin( th );
	m[3][0] = 0.0f;
	m[0][1] = 0.0f;
	m[1][1] = 1.0f;
	m[2][1] = 0.0f;
	m[3][1] = 0.0f;
	m[0][2] = -sin( th );
	m[1][2] = 0;
	m[2][2] = cos( th );
	m[3][2] = 0.0f;
	m[0][3] = 0.0f;
	m[1][3] = 0.0f;
	m[2][3] = 0.0f;
	m[3][3] = 1.0f;

}
void rotxK( Matrix m,					/* resulting rotation matrix (out) */
	  float th					/* rotation vector (in) */
	 )
{
	m[0][0] = 1.0f;
	m[1][0] = 0.0f;
	m[2][0] = 0.0f;
	m[3][0] = 0.0f;
	m[0][1] = 0.0f;
	m[1][1] = cos( th );
	m[2][1] = -sin( th );
	m[3][1] = 0.0f;
	m[0][2] = 0.0f;
	m[1][2] = sin( th );
	m[2][2] = cos( th );
	m[3][2] = 0.0f;
	m[0][3] = 0.0f;
	m[1][3] = 0.0f;
	m[2][3] = 0.0f;
	m[3][3] = 1.0f;

}


VERT vxm3x3K( Matrix m,					/* input matrix (in) */
	 VERT v						/* input vector (in) */
	 )
{
	VERT            r;

	r.x = v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] ;
	r.y = v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] ;
	r.z = v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] ;
	return ( r );
}


void displace_scaleK( STATICGUIDE * h, float rs, MINIKERNAL *m)
{
	int             xx;

	if( rs != 1.0f )
		for( xx = 1; xx < m->guide_segs; xx++ )
		{
			h->hv[xx].x -= h->hv[0].x;
			h->hv[xx].y -= h->hv[0].y;
			h->hv[xx].z -= h->hv[0].z;
			h->hv[xx].x *= rs;
			h->hv[xx].y *= rs;
			h->hv[xx].z *= rs;
			h->hv[xx].x += h->hv[0].x;
			h->hv[xx].y += h->hv[0].y;
			h->hv[xx].z += h->hv[0].z;
		}

}



void make_basematrixK(STATICGUIDE *h,BASEMATRIXK *bm,MINIKERNAL *m)
{
int qv;
								VERT xv,yv,zv;
								VERT tangent[15];

					for( qv = 1; qv < m->guide_segs; qv++ )
							{
								STATICGUIDE       *bh, *bh1;
								VERT zero;
								int             q1, q2, q0;
								VERT            yv2,b;

								zero.x=0.0f;
								zero.y=0.0f;
								zero.z=0.0f;
								q0=qv-1;
								q1 = qv;
								q2 = qv + 1;
								
								yv.x = h->hv[q2].x - h->hv[q1].x;
								yv.y = h->hv[q2].y - h->hv[q1].y;
								yv.z = h->hv[q2].z - h->hv[q1].z;

								yv2.x = h->hv[q1].x - h->hv[q0].x;
								yv2.y = h->hv[q1].y - h->hv[q0].y;
								yv2.z = h->hv[q1].z - h->hv[q0].z;
								VnormalizeK(&yv);
								VnormalizeK(&yv2);
								yv.x = ( yv.x + yv2.x ) / 2.0f;
								yv.y = ( yv.y + yv2.y ) / 2.0f;
								yv.z = ( yv.z + yv2.z ) / 2.0f;
								VnormalizeK(&yv);
								tangent[qv]=yv;
							}
							tangent[0]=h->normal;
							tangent[m->guide_segs-1]=tangent[m->guide_segs-2];



	
					for( qv = 0; qv < m->guide_segs; qv++ )
							{
								STATICGUIDE       *bh, *bh1;
								VERT zero;
								int             q1, q2, qb, q3, q4;
								VERT            yv2,b;

								zero.x=0.0f;
								zero.y=0.0f;
								zero.z=0.0f;
								bm->pos[qv]=h->hv[qv];


yv=tangent[qv];
if (qv==0)
										xv = h->handle;

								{
									 VnormalizeK( &xv );
									VnormalizeK( &yv );
									zv = VcrossK( xv, yv );
									 VnormalizeK( &zv );
									xv = VcrossK( yv, zv );
									VnormalizeK( &xv );

									b = h->hv[qv];
									mMatrix2K( bm->zero2world[qv], xv, yv, zv, zero );
									inverseK(bm->zero2world[qv],bm->world2zero[qv]);
									
								}



							}

}





void twist_hairK (STATICGUIDE *clumpCenter, STATICGUIDE *h, MINIKERNAL *m, float *texture,long seed)
{ // rotate the matrix using clumprot
BASEMATRIXK bm;
float rot_dir,offset_dir;
int qv;
float randomize;
seed=seed*5000;
rot_dir=(randcl( seed++ )-0.5f);
rot_dir=(float)SGN(rot_dir);

offset_dir=1.0f;
randomize=texture[RANDOMIZE]*(randcl(seed++))+(1.0f-texture[CLUMP_RANDOMIZE]);


make_basematrixK(clumpCenter,&bm,m);
for (qv=0;qv<m->guide_segs;qv++)
		{
			float qq;
			Matrix          rot;
			VERT xv,yv,zv,trans;
			float uu1;
			float uuu;
			float flat;
			uuu=(float)((float)qv/(float)m->guide_segs);
uu1=uuu;
			uuu*=uuu;
			trans.x=0.0f;
			trans.y=0.0f;
			trans.z=0.0f;
			h->hv[qv].x-=clumpCenter->hv[qv].x;
			h->hv[qv].y-=clumpCenter->hv[qv].y;
			h->hv[qv].z-=clumpCenter->hv[qv].z;
			h->hv[qv]=vxm3x3K(
				bm.world2zero[qv],
				h->hv[qv]);
			h->hv[qv].x+=offset_dir*texture[CLUMP_ROT_OFFSET]*uuu;	// back here
			qq = rot_dir*texture[CLUMP_ROT_STRENGTH]*uuu;
			qq *= randomize;
			flat=1.0-texture[CLUMP_FLATNESS];
			h->hv[qv].z*= flat*uu1+(1.0f-uu1); // flatness
rotyK(rot,qq);
			h->hv[qv] = vxm3x3K( rot, h->hv[qv] ); // spin it

			h->hv[qv]=vxm3x3K(bm.zero2world[qv],h->hv[qv]);// transformn it back
			h->hv[qv].x+=clumpCenter->hv[qv].x;
			h->hv[qv].y+=clumpCenter->hv[qv].y;
			h->hv[qv].z+=clumpCenter->hv[qv].z;

		}

}

void displace_randscaleK( STATICGUIDE * h, MINIKERNAL *m,int seed,float rss )
{
	int             xx;
	float           rs = 1.0f;
//	float rss;
//	rss=texture[RANDSCALE];
	rss*=0.7f;

	rs = randcl( seed++ );
	rs = ( rss ) * rs + ( 1.0 - rss ) * 1.0f;

	if( rs != 1.0f )
		for( xx = 1; xx < m->guide_segs; xx++ )
		{
			h->hv[xx].x -= h->hv[0].x;
			h->hv[xx].y -= h->hv[0].y;
			h->hv[xx].z -= h->hv[0].z;
			h->hv[xx].x *= rs;
			h->hv[xx].y *= rs;
			h->hv[xx].z *= rs;
			h->hv[xx].x += h->hv[0].x;
			h->hv[xx].y += h->hv[0].y;
			h->hv[xx].z += h->hv[0].z;
// possible look difference - we didn't scale the rest verts before
//			h->resthv[xx].x -= h->resthv[0].x;
//			h->resthv[xx].y -= h->resthv[0].y;
//			h->resthv[xx].z -= h->resthv[0].z;
//			h->resthv[xx].x *= rs;
//			h->resthv[xx].y *= rs;
//			h->resthv[xx].z *= rs;
//			h->resthv[xx].x += h->resthv[0].x;
//			h->resthv[xx].y += h->resthv[0].y;
//			h->resthv[xx].z += h->resthv[0].z;

		}

}



static void    displace_frizzK( STATICGUIDE * h, MINIKERNAL * m, STATICROOTPT *pt, float restlength,int seed)
{

	int             xx;
	float           freqqx, freqqy, freqqz;
	VERT            dv;
	VERT            dv2;
	Matrix          mr1, mr2, mr3;
	float           freqx, freqy, freqz;
	int             oc;
	int pp=0;
	float           amp = 1.0f;
	float wt;
	int fly=0;
	float tfrizz;
	float *texture;
texture= pt->baked_tex;
tfrizz=0.8f*texture[TIPFRIZZ]; // concatinated alread


// the old restlength was 1/15th the length of the hair, now it's realative to 1/guide_segs.

restlength*=m->guide_segs; 
restlength/=15.0f;

if (texture[FLYAWAY]>0.0f)
	{
		float wt;
		for (xx=0;xx<20*randcl(seed++)+10;xx++) // necessary - optimise this later
		wt=randcl(seed++);
		if (wt<texture[FLYAWAY]) 
								{
								fly=1;
								tfrizz/=2.0f;
								}
				else fly=0;
	}



	freqqx = texture[FRIZZFREQX] * 0.5f;
	freqqy = texture[FRIZZFREQY] * 0.5f;
	freqqz = texture[FRIZZFREQZ] * 0.5f;


	if (texture[FRIZZANIM]>0.0f||(texture[ROOTFRIZZ]>0.0f)||texture[TIPFRIZZ]>0.0f||texture[FLYAWAY]>0.0f||texture[MESS]>0.0f)
	{
		VERT animdir;
		VERT            offs;
		dv.x = 0;
		dv.y = 0;
		dv.z = 0;
		{
			amp = 1.0f;
			freqx = freqqx / ( restlength * 111.1f );
			freqy = freqqy / ( restlength * 111.1f );
			freqz = freqqz / ( restlength * 111.1f );

			if( m->interpolation_type == SPLINE_GROWTH )
				freqx *= 56.0f;
			if( m->interpolation_type == SPLINE_GROWTH )
				freqy *= 56.0f;
			if( m->interpolation_type == SPLINE_GROWTH )
				freqz *= 56.0f;

#ifdef EXTERNAL_COLLISION
// backwards compatibility
			if( m->interpolation_type == SPLINE_GROWTH )
			{
				freqx *= 5.0f;
				freqy *= 5.0f;
				freqz *= 5.0f;	// unlock
			}

#endif

			animdir.x = 0.0f; // need to follow up on this
			animdir.y = 1.0f;
			animdir.z = 0.0f;
			offs.x = 0;
			offs.y = 0;
			offs.z = 0;

			dv.x = (getnoise_octavesK(m,h->noisev[pp].x*freqx+offs.x+.4,
										h->noisev[pp].y*freqy+offs.y,
										h->noisev[pp].z*freqz+offs.z-0.5f,3));

			dv.y = (getnoise_octavesK(m,h->noisev[pp].x*freqx+offs.x+0.1f,
										h->noisev[pp].y*freqy+offs.y+0.5f,
										h->noisev[pp].z*freqz+offs.z-0.2f,3));

			dv.z = (getnoise_octavesK(m,h->noisev[pp].x*freqx+offs.x,
										h->noisev[pp].y*freqy+offs.y-0.1f,
										h->noisev[pp].z*freqz+offs.z+0.5f,3));


			dv.x *= ( 1.0 - texture[FRIZZANIM] );
			dv.y *= ( 1.0 - texture[FRIZZANIM] );
			dv.z *= ( 1.0 - texture[FRIZZANIM] );

			animdir.x = 0.0f;
			animdir.y = 1.0f;
			animdir.z = 0.0f;

			offs.x = m->time * freqx * texture[FRIZZANIMSPEED] * animdir.x;
			offs.y = m->time * freqy * texture[FRIZZANIMSPEED] * animdir.y;
			offs.z = m->time * freqz * texture[FRIZZANIMSPEED] * animdir.z;

			dv2.x = 0;
			dv2.y = 0;
			dv2.z = 0;
			if( texture[FRIZZANIM] != 0.0f )
			{
			dv2.x = (getnoise_octavesK(m,h->noisev[pp].x*freqx+offs.x+.9,
										h->noisev[pp].y*freqy+offs.y-.2,
										h->noisev[pp].z*freqz+offs.z-.1,2));

			dv2.y = (getnoise_octavesK(m,h->noisev[pp].x*freqx+offs.x+0.11f,
										h->noisev[pp].y*freqy+offs.y+0.15f,
										h->noisev[pp].z*freqz+offs.z-0.62f,2));

			dv2.z = (getnoise_octavesK(m,h->noisev[pp].x*freqx+offs.x,
										h->noisev[pp].y*freqy+offs.y-0.71f,
										h->noisev[pp].z*freqz+offs.z+0.25f,2));

				dv2.x *= texture[FRIZZANIM];
				dv2.y *= texture[FRIZZANIM];
				dv2.z *= texture[FRIZZANIM];

			}

		}


			dv.x += dv2.x;
			dv.y += dv2.y;
			dv.z += dv2.z;

			dv.x = dv.x + ( dv.x * .05 ) * randcl( seed++ );
			dv.y = dv.y + ( dv.y * .05 ) * randcl( seed++ );
			dv.z = dv.z + ( dv.z * .05 ) * randcl( seed++ );


			dv.x = -dv.x;
			dv.y = -dv.y;
			dv.z = -dv.z;



			if ((texture[ROOTFRIZZ]>0.0f)||texture[TIPFRIZZ]||texture[FLYAWAY]>0.0f||texture[MESS]>0.0f)
			{
				float st;
				VERT dv3;
				float           rff, tff;
// tip frizz and root frizz are reversed
				rff = texture[ROOTFRIZZ] / 6.0f;
				tff = tfrizz / 6.0f;
				rff*=7.5f;//*15.0f;
				tff*=7.5f;//*15.0f;
				if( m->interpolation_type == SPLINE_GROWTH ) // SPLINE_GROWTH
				{
					rff /= 4.0f;
					tff /= 4.0f;
				}
				for( xx = 1; xx < m->guide_segs; xx++ )
				{
					float           qq;

					qq = ( xx / ( float ) m->guide_segs );
					qq = ( ( float ) 1.0 - qq ) * rff + qq * tff;
					qq *= ( float ) 3.14 / ( float ) 180.0f;
					qq *= ( float ) .5;
//if (h->its_flyaway)
					if (fly)
					qq += texture[FLYAWAY]*10.0f*( xx / ( float ) m->guide_segs );

					{
						rotxK( mr1, dv.x*qq ); // load up rotation matricies
						rotzK( mr3, dv.z*qq );
						rotyK( mr2, dv.y*qq );
					}
					h->hv[xx].x -= h->hv[0].x; // bring it back to origin
					h->hv[xx].y -= h->hv[0].y;
					h->hv[xx].z -= h->hv[0].z;

					h->hv[xx] = vxmK( mr1, h->hv[xx] );
					h->hv[xx] = vxmK( mr3, h->hv[xx] );
					h->hv[xx] = vxmK( mr2, h->hv[xx] );

					if (fly)
					{
					float hm;
					float ddx,ddy,ddz;
					hm=texture[MESS]*restlength;

					h->hv[xx].x+=hm*(randcl(seed++)-0.5f); // mess
					h->hv[xx].y+=hm*(randcl(seed++)-0.5f);
					h->hv[xx].z+=hm*(randcl(seed++)-0.5f);
						{
						h->hv[xx].x*=(1.0f+texture[FLYAWAY]*1.2f); // flyaway - multipliers
						h->hv[xx].y*=(1.0f+texture[FLYAWAY]*1.2f);
						h->hv[xx].z*=(1.0f+texture[FLYAWAY]*1.2f);
						}
					}

					h->hv[xx].x += h->hv[0].x;
					h->hv[xx].y += h->hv[0].y;
					h->hv[xx].z += h->hv[0].z;
				}
			}

		}
}






float randcl(long seed)
{
float ret;
// this algorithm makes use of some overflow, hopefully it works same on gpu
 static long a ; // could be made the seed value
 a=seed *500;
a = (a * 32719 + 3) % 32749;
return (   fabs((float)(a % RAND_MAX) + 1)/(float)RAND_MAX  );
}

void bary_interpv(VERT *ret,VERT *va, VERT *vb, VERT *vc, float a, float b, float c)
{
if (c>=0.0f) // it's 3 guide interp C== -1 if it's 2
       {
       ret->x=va->x*a+vb->x*b+vc->x*c;
       ret->y=va->y*a+vb->y*b+vc->y*c;
       ret->z=va->z*a+vb->z*b+vc->z*c;
	}
else
       {
       ret->x=va->x*a+vb->x*b;
       ret->y=va->y*a+vb->y*b;
       ret->z=va->z*a+vb->z*b;
       }


}

float bary_interpf(float va, float vb, float vc, float a, float b, float c)
{
	float ret;
if (c>=0.0f) // it's 3 guide interp C== -1 if it's 2
       ret=va*a+
	   vb*b+
	   vc*c;
else
       ret=va*a+vb*b;
return(ret);
}





void gen_hairK (MINIKERNAL *m, STATICROOTPT *inpt, STATICGUIDE *ret)
{
STATICGUIDE g1,g2,g3;
float a,b,c;
int x;
int breakit;
int qq;
float flip;
long seed;
a=inpt->bary.x;
b=inpt->bary.y;
c=inpt->bary.z;



fetch_staticguide(m, &g1, inpt->pntid1);
fetch_staticguide(m, &g2, inpt->pntid2);
if (c>0)
fetch_staticguide(m, &g3, inpt->pntid3);

breakit=1;
if (m->interpolation_type==MESH_GROWTH)
{
if ( (g1.smoothing_group==g2.smoothing_group)&&
 (g2.smoothing_group==g3.smoothing_group) )
     breakit=0;
if (g1.splitmerge+g2.splitmerge+g3.splitmerge > 0)
breakit = 0;
}
else
{
if ( (g1.smoothing_group==g2.smoothing_group))
     breakit=0;
if (g1.splitmerge+g2.splitmerge > 0)
breakit = 0;
}

qq=m->guide_segs;
if (breakit) qq=1;

seed=inpt->id; // need to add the offset

flip=randcl(seed);

// interpolate guide data using bary's from roots
for (x=0;x<qq;x++)
    {
// optimisation note, we can make a mask to only interp the values
// we're going to need, ie rest hairs only
	bary_interpv(&ret->hv[x],&g1.hv[x],&g2.hv[x],&g3.hv[x],a,b,c);
	bary_interpv(&ret->resthv[x],&g1.resthv[x],&g2.resthv[x],&g3.resthv[x],a,b,c);
	bary_interpv(&ret->velocity[x],&g1.velocity[x],&g2.velocity[x],&g3.velocity[x],a,b,c);
	bary_interpv(&ret->noisev[x],&g1.noisev[x],&g2.noisev[x],&g3.noisev[x],a,b,c);
    }



if ((breakit)&&(m->interpolation_type==SPLINE_GROWTH))
{
	VERT tmp;
for (x=1;x<m->guide_segs;x++)
 {
if (flip<.5)
	{
	STATICGUIDE *gg;
	gg=&g1;
	tmp=gg->hv[x];
	tmp.x-=gg->hv[0].x;
	tmp.y-=gg->hv[0].y;
	tmp.z-=gg->hv[0].z;
	ret->hv[x].x=ret->hv[0].x+tmp.x;
	ret->hv[x].y=ret->hv[0].y+tmp.y;
	ret->hv[x].z=ret->hv[0].z+tmp.z;

	tmp=gg->resthv[x];
	tmp.x-=gg->resthv[0].x;
	tmp.y-=gg->resthv[0].y;
	tmp.z-=gg->resthv[0].z;
	ret->resthv[x].x=ret->resthv[0].x+tmp.x;
	ret->resthv[x].y=ret->resthv[0].y+tmp.y;
	ret->resthv[x].z=ret->resthv[0].z+tmp.z;

	tmp=gg->velocity[x];
	tmp.x-=gg->velocity[0].x;
	tmp.y-=gg->velocity[0].y;
	tmp.z-=gg->velocity[0].z;
	ret->velocity[x].x=ret->velocity[0].x+tmp.x;
	ret->velocity[x].y=ret->velocity[0].y+tmp.y;
	ret->velocity[x].z=ret->velocity[0].z+tmp.z;

	tmp=gg->noisev[x];
	tmp.x-=gg->noisev[0].x;
	tmp.y-=gg->noisev[0].y;
	tmp.z-=gg->noisev[0].z;
	ret->noisev[x].x=ret->noisev[0].x+tmp.x;
	ret->noisev[x].y=ret->noisev[0].y+tmp.y;
	ret->noisev[x].z=ret->noisev[0].z+tmp.z;

	}
if (flip>=.5)
	{
	STATICGUIDE *gg;
	gg=&g2;
	tmp=gg->hv[x];
	tmp.x-=gg->hv[0].x;
	tmp.y-=gg->hv[0].y;
	tmp.z-=gg->hv[0].z;
	ret->hv[x].x=ret->hv[0].x+tmp.x;
	ret->hv[x].y=ret->hv[0].y+tmp.y;
	ret->hv[x].z=ret->hv[0].z+tmp.z;

	tmp=gg->resthv[x];
	tmp.x-=gg->resthv[0].x;
	tmp.y-=gg->resthv[0].y;
	tmp.z-=gg->resthv[0].z;
	ret->resthv[x].x=ret->resthv[0].x+tmp.x;
	ret->resthv[x].y=ret->resthv[0].y+tmp.y;
	ret->resthv[x].z=ret->resthv[0].z+tmp.z;

	tmp=gg->velocity[x];
	tmp.x-=gg->velocity[0].x;
	tmp.y-=gg->velocity[0].y;
	tmp.z-=gg->velocity[0].z;
	ret->velocity[x].x=ret->velocity[0].x+tmp.x;
	ret->velocity[x].y=ret->velocity[0].y+tmp.y;
	ret->velocity[x].z=ret->velocity[0].z+tmp.z;

	tmp=gg->noisev[x];
	tmp.x-=gg->noisev[0].x;
	tmp.y-=gg->noisev[0].y;
	tmp.z-=gg->noisev[0].z;
	ret->noisev[x].x=ret->noisev[0].x+tmp.x;
	ret->noisev[x].y=ret->noisev[0].y+tmp.y;
	ret->noisev[x].z=ret->noisev[0].z+tmp.z;
	}
 }
}




if ((breakit)&&(m->interpolation_type==MESH_GROWTH))
{
	VERT tmp;
for (x=1;x<m->guide_segs;x++)
 {
if (flip<.3)
	{
	STATICGUIDE *gg;
	gg=&g1;
	tmp=gg->hv[x];
	tmp.x-=gg->hv[0].x;
	tmp.y-=gg->hv[0].y;
	tmp.z-=gg->hv[0].z;
	ret->hv[x].x=ret->hv[0].x+tmp.x;
	ret->hv[x].y=ret->hv[0].y+tmp.y;
	ret->hv[x].z=ret->hv[0].z+tmp.z;

	tmp=gg->resthv[x];
	tmp.x-=gg->resthv[0].x;
	tmp.y-=gg->resthv[0].y;
	tmp.z-=gg->resthv[0].z;
	ret->resthv[x].x=ret->resthv[0].x+tmp.x;
	ret->resthv[x].y=ret->resthv[0].y+tmp.y;
	ret->resthv[x].z=ret->resthv[0].z+tmp.z;

	tmp=gg->velocity[x];
	tmp.x-=gg->velocity[0].x;
	tmp.y-=gg->velocity[0].y;
	tmp.z-=gg->velocity[0].z;
	ret->velocity[x].x=ret->velocity[0].x+tmp.x;
	ret->velocity[x].y=ret->velocity[0].y+tmp.y;
	ret->velocity[x].z=ret->velocity[0].z+tmp.z;

	tmp=gg->noisev[x];
	tmp.x-=gg->noisev[0].x;
	tmp.y-=gg->noisev[0].y;
	tmp.z-=gg->noisev[0].z;
	ret->noisev[x].x=ret->noisev[0].x+tmp.x;
	ret->noisev[x].y=ret->noisev[0].y+tmp.y;
	ret->noisev[x].z=ret->noisev[0].z+tmp.z;

	}
if ((flip>=.3)&&(flip<.6))
	{
	STATICGUIDE *gg;
	gg=&g2;
	tmp=gg->hv[x];
	tmp.x-=gg->hv[0].x;
	tmp.y-=gg->hv[0].y;
	tmp.z-=gg->hv[0].z;
	ret->hv[x].x=ret->hv[0].x+tmp.x;
	ret->hv[x].y=ret->hv[0].y+tmp.y;
	ret->hv[x].z=ret->hv[0].z+tmp.z;

	tmp=gg->resthv[x];
	tmp.x-=gg->resthv[0].x;
	tmp.y-=gg->resthv[0].y;
	tmp.z-=gg->resthv[0].z;
	ret->resthv[x].x=ret->resthv[0].x+tmp.x;
	ret->resthv[x].y=ret->resthv[0].y+tmp.y;
	ret->resthv[x].z=ret->resthv[0].z+tmp.z;

	tmp=gg->velocity[x];
	tmp.x-=gg->velocity[0].x;
	tmp.y-=gg->velocity[0].y;
	tmp.z-=gg->velocity[0].z;
	ret->velocity[x].x=ret->velocity[0].x+tmp.x;
	ret->velocity[x].y=ret->velocity[0].y+tmp.y;
	ret->velocity[x].z=ret->velocity[0].z+tmp.z;

	tmp=gg->noisev[x];
	tmp.x-=gg->noisev[0].x;
	tmp.y-=gg->noisev[0].y;
	tmp.z-=gg->noisev[0].z;
	ret->noisev[x].x=ret->noisev[0].x+tmp.x;
	ret->noisev[x].y=ret->noisev[0].y+tmp.y;
	ret->noisev[x].z=ret->noisev[0].z+tmp.z;
	}
if (flip>=0.6f)
	{
	STATICGUIDE *gg;
	gg=&g3;
	tmp=gg->hv[x];
	tmp.x-=gg->hv[0].x;
	tmp.y-=gg->hv[0].y;
	tmp.z-=gg->hv[0].z;
	ret->hv[x].x=ret->hv[0].x+tmp.x;
	ret->hv[x].y=ret->hv[0].y+tmp.y;
	ret->hv[x].z=ret->hv[0].z+tmp.z;

	tmp=gg->resthv[x];
	tmp.x-=gg->resthv[0].x;
	tmp.y-=gg->resthv[0].y;
	tmp.z-=gg->resthv[0].z;
	ret->resthv[x].x=ret->resthv[0].x+tmp.x;
	ret->resthv[x].y=ret->resthv[0].y+tmp.y;
	ret->resthv[x].z=ret->resthv[0].z+tmp.z;

	tmp=gg->velocity[x];
	tmp.x-=gg->velocity[0].x;
	tmp.y-=gg->velocity[0].y;
	tmp.z-=gg->velocity[0].z;
	ret->velocity[x].x=ret->velocity[0].x+tmp.x;
	ret->velocity[x].y=ret->velocity[0].y+tmp.y;
	ret->velocity[x].z=ret->velocity[0].z+tmp.z;

	tmp=gg->noisev[x];
	tmp.x-=gg->noisev[0].x;
	tmp.y-=gg->noisev[0].y;
	tmp.z-=gg->noisev[0].z;
	ret->noisev[x].x=ret->noisev[0].x+tmp.x;
	ret->noisev[x].y=ret->noisev[0].y+tmp.y;
	ret->noisev[x].z=ret->noisev[0].z+tmp.z;
	}
 }
}


ret->uu=bary_interpf(g1.uu,g2.uu,g3.uu,a,b,c);
ret->vv=bary_interpf(g1.vv,g2.vv,g3.vv,a,b,c);
bary_interpv(&ret->normal,&g1.normal,&g2.normal,&g3.normal,a,b,c);
bary_interpv(&ret->rest_normal,&g1.rest_normal,&g2.rest_normal,&g3.rest_normal,a,b,c);
bary_interpv(&ret->handle,&g1.handle,&g2.handle,&g3.handle,a,b,c);
VnormalizeK(&ret->normal);
VnormalizeK(&ret->rest_normal);
VnormalizeK(&ret->handle);

// need to add dice rolls for splitgroups

}

extern void the_kernal (MINIKERNAL *m,STATICGUIDE *ret,STATICROOTPT *pt,int hairnum)
{
STATICGUIDE clumphair;
STATICGUIDE scratch;
STATICGUIDE h;
STATICGUIDE ttt;
STATICROOTPT clump_pt;
VERT t1[GUIDE_MAX_SEGS];
VERT t2[GUIDE_MAX_SEGS];
BASEMATRIXK bm;
int x;
float restlength;
int killme=0;
	int y;

long seed;

if (hairnum<m->total_roots)
{
get_rootptK(pt,m, hairnum); // pull in rootpt from global data
seed=(long) pt->id;
h.killme=pt->killme;
if (pt->killme==0)
{
gen_hairK(m,pt, &h); // tmphair is output

restlength=VdistanceK(h.resthv[0],h.resthv[1]);
if (pt->baked_tex[CUTLENGTH]<1.0f/(float)m->guide_segs) killme=1;
if (restlength<0.001f) killme=1;
if (killme==1)
pt->killme=1;
x=0;

if (!killme)
{
make_basematrixK(&h,&bm,m);

displace_scaleK(&h,pt->baked_tex[SCALE],m);			
seed=(long) pt->id;
displace_randscaleK( &h, m,seed,pt->baked_tex[RANDSCALE] );

if (pt->closest_clump>=0)
if (m->total_clumps>0)
{
float rs;
get_rootptK(&clump_pt, m,pt->closest_clump); 
gen_hairK(m,&clump_pt, &clumphair);
displace_scaleK(&clumphair,pt->baked_tex[SCALE],m);			
displace_randscaleK( &clumphair, m,seed,pt->baked_tex[RANDSCALE] );

													for (x=1;x<m->guide_segs;x++) // apply clumps (test)
													{
														float ii,aa,bb;
														float flat;
														float ccl;
														ii=((float)x/(float)m->guide_segs);
														aa=(1.0-ii);
														bb=(ii);
														flat=1.0f-pt->baked_tex[CLUMP_FLATNESS];		
														ccl=pt->baked_tex[CLUMP_STRENGTH];
														h.hv[x].x=(h.hv[x].x*aa+clumphair.hv[x].x*bb)*ccl+h.hv[x].x*(1.0-ccl);
														h.hv[x].y=(h.hv[x].y*aa+clumphair.hv[x].y*bb)*ccl+h.hv[x].y*(1.0-ccl);
														h.hv[x].z=(h.hv[x].z*aa+clumphair.hv[x].z*bb)*ccl+h.hv[x].z*(1.0-ccl);
//														h.hv[qv]=clumphair.hv[qv]; // lets look at the center
													}
			seed=(long) pt->id;
													twist_hairK (&clumphair,&h,m,pt->baked_tex,seed);
									
}


for (x=0;x<m->guide_segs;x++)
	t2[x]=h.hv[x];

for( x = 0; x < m->guide_segs; x++ )
   {										   t1[x].x = 0;
   t1[x].y = ( ( float ) x / (float)m->guide_segs )*restlength* (float)m->guide_segs;
   t1[x].z = 0;
   ttt.hv[x]=t1[x];
   ttt.noisev[x]=h.noisev[x];
   ttt.resthv[x]=h.resthv[x];
   }

seed=pt->id;


seed=pt->id;
if ((pt->baked_tex[KINKROOT]!=0.0f)||(pt->baked_tex[KINK]!=0.0f))
displace_kinkyK( &ttt, m,pt->baked_tex, restlength,seed);
if ((pt->baked_tex[ROOTFRIZZ]!=0.0f)||(pt->baked_tex[TIPFRIZZ]!=0.0f))
displace_frizzK( &ttt, m,pt,restlength,seed );
for( x = 0; x < m->guide_segs; x++ )
{
h.hv[x].x=ttt.hv[x].x-t1[x].x;
h.hv[x].y=ttt.hv[x].y-t1[x].y;
h.hv[x].z=ttt.hv[x].z-t1[x].z;



h.hv[x] = vxmK( bm.zero2world[x], h.hv[x]);					
h.hv[x].x += t2[x].x;
h.hv[x].y += t2[x].y;
h.hv[x].z += t2[x].z;
}



// put it into returnhair
displace_randscaleK( &h, m,seed,pt->baked_tex[RANDSCALE] );


memcpy(&scratch,&h,sizeof(STATICGUIDE));
	scratch.handle.z*= -1.0f;
	scratch.normal.z*= -1.0f;
	scratch.rest_handle.z*= -1.0f;
for (x=0;x<m->guide_segs;x++)
{
	scratch.hv[x].z*= -1.0f;
	scratch.velocity[x].z*= -1.0f;
	scratch.noisev[x].z*= -1.0f;
	scratch.resthv[x].z*= -1.0f;
}
	/* poll external forces

if (m->ext_forces!=0.0f)
for (x=0;x<guide_segs;x++)
{
float uu;
VERT vect;
uu=(float)x/(float)m->guide_segs;
vect.x=0.0f;vect.y=0.0f;vect.z=0.0f;
h.hv[x].x*= -1.0f;
h.hv[x].y*= -1.0f;
h.hv[x].z*= -1.0f;
MAYAexternal_forces( &h.hv[x], &vect, x );
h.hv[x].x+=vect.x*uu*m->ext_forces;
h.hv[x].y+=vect.y*uu*m->ext_forces;
h.hv[x].z+=vect.z*uu*m->ext_forces;
h.hv[x].x*= -1.0f;
h.hv[x].y*= -1.0f;
h.hv[x].z*= -1.0f;
}
*/





if (pt->baked_tex[CUTLENGTH]<.01) pt->killme=1; // optimisation note - put this at the top

if (pt->baked_tex[CUTLENGTH]<1.0f)
	{
// TODO	resize_static_guide_hair_to_newguide(&h,&m->return_hair[x],m->segs,pt->baked_tex[CUTLENGTH]);
	}

if ((pt->killme)||(killme))
for (x=0;x<m->guide_segs;x++)
	{
scratch.hv[x]=scratch.hv[0];
scratch.velocity[x]=scratch.velocity[0];

	}


if (m->multiplier<=1)
{
	color_a_hairK(m,pt,&scratch,pt->id);
	put_staticguideRH(m, &scratch, hairnum);
}

if (m->multiplier>1)
{
seed=pt->id;
make_basematrixK(&h,&bm,m);

for (x=0;x<m->guide_segs;x++)
	t2[x]=h.hv[x];
for (y=0;y< pt->baked_tex[DREADCOUNT];y++)
{
	int hn;
for( x = 0; x < m->guide_segs; x++ )
   {
   float u;
   t1[x].x = 0;
   t1[x].y = 0;
   t1[x].z = 0;
   u=( ( float ) x / (float)m->guide_segs )  ;
   ttt.hv[x]=t1[x];// + displacement
  ttt.noisev[x]=h.noisev[x];
   ttt.resthv[x]=h.resthv[x];
 	}



													{
														float ss;
														ss=(float) randcl(seed);
														seed++;
														if (ss<.5) ss= -1.0f;
														if (ss>=.5) ss= 1.0f;

														{
															float ang, mag;
															Matrix          pm;
															VERT            foffs;
															float           sc = 0.0f;
															int             ll = 0;

															{

																ang= (float) (randcl( seed )*6.284-3.14);
																														seed++;
																mag= (float) randcl( seed );
																														seed++;

																foffs.x=cos(ang)*mag;
																foffs.z=sin(ang)*mag;
																foffs.y=0.0f;
															}



															if( m->interpolation_type==MESH_GROWTH )
																sc = pt->restlength * 6.0f;
															if( m->interpolation_type==SPLINE_GROWTH )
																sc = pt->restlength / 5.0f;	// * 6.0f;
//			                                                 sc=5.0f;
															foffs.x *= sc;
															foffs.y *= sc;
															foffs.z *= sc;
															ang= (float) (randcl( seed )*6.284-3.14);
															seed++;
															mag= (float) randcl( seed );
															seed++;
															for( ll = 0; ll < m->guide_segs; ll++ )
															{
																float           vs;
																float           ivs;

																vs = ( float ) ll / (float)m->guide_segs;
																if( vs > 1.0 )
																	vs = 1.0f;
																ivs = ( float ) ll / (float)m->guide_segs;
																if( ivs > 1.0 )
																	ivs = 1.0f;
																ivs = 1.0f - ivs;
																vs = ( ( float ) vs ) * pt->baked_tex[DREADTIP] + ( ivs ) * pt->baked_tex[DREADROOT];
																if( vs < 0 )
																	vs = 0;
																{
																	float ang2;
																	float offx;
																	float uu;
																	uu=(float)ll/(float)m->guide_segs;
																	uu=pow(uu,1.25);
																	ang2=uu*1.28f*pt->baked_tex[MULTASP]*ss;
																	foffs.x=cos(ang)*mag*vs*sc;
																	foffs.z=sin(ang)*mag*vs*sc;
																	foffs.z*=pt->baked_tex[ASPECT];
																	foffs.x+=(pt->baked_tex[OFFSET]*sc*uu);

																	foffs.y=0;
																	{
																		VERT off2;
																		off2=foffs;
																		foffs.x=off2.x*cos(ang2)-off2.z*sin(ang2);
																		foffs.z=off2.x*sin(ang2)+off2.z*cos(ang2);
																	}

																	foffs = vxmK( bm.zero2world[ll], foffs );
																}

																ttt.hv[ll].x = foffs.x;// * vs;
																ttt.hv[ll].y = foffs.y;// * vs;
																ttt.hv[ll].z = foffs.z;// * vs;
															}
														}
													}









for( x = 0; x < m->guide_segs; x++ )
	{
	scratch.hv[x].x=ttt.hv[x].x;//-t1[x].x;
	scratch.hv[x].y=ttt.hv[x].y;//-t1[x].y;
	scratch.hv[x].z=ttt.hv[x].z;//-t1[x].z;
	scratch.hv[x] = vxmK( bm.zero2world[x], scratch.hv[x]);					
	scratch.hv[x].x += t2[x].x;
	scratch.hv[x].y += t2[x].y;
	scratch.hv[x].z += t2[x].z;
	scratch.hv[x].z*= -1.0f;
	scratch.velocity[x].z*= -1.0f;

	}
	hn=m->total_roots*y+hairnum;
	color_a_hairK(m,pt,&scratch,pt->id*y);
	put_staticguideRH(m, &scratch, hn);



}


}
}
}// killme=0
} // hairnum<total_roots
if (hairnum>=m->total_roots) pt->killme=1;
}


void color_a_hairK(MINIKERNAL *m,STATICROOTPT *pt,STATICGUIDE *h,int seed)
{
	float d;
	float rv1;
	VERT rc;
	h->rootcolor.x=(float)pt->baked_tex[ROOT_COLOR_R];//255.0f;   
	h->rootcolor.y=(float)pt->baked_tex[ROOT_COLOR_G];//255.0f;   
	h->rootcolor.z=(float)pt->baked_tex[ROOT_COLOR_B];//255.0f;   
	h->tipcolor.x=(float)pt->baked_tex[TIP_COLOR_R];//255.0f;   
	h->tipcolor.y=(float)pt->baked_tex[TIP_COLOR_G];//255.0f;   
	h->tipcolor.z=(float)pt->baked_tex[TIP_COLOR_B];//255.0f;

	d=randcl(seed++);
if (m->interpolation_type!=SPLINE_GROWTH) // no density on spline growths
	if (pt->baked_tex[DENSITY]*pt->baked_tex[DENSITY]<d) h->killme=1;
	else h->killme=0;

	rv1 = randcl(seed++);
//	rv1 = smoothstep(rv1);
//	rv1 = smoothstep(rv1);

	rc.x = randcl(seed++);
	rc.y = randcl(seed++);
	rc.z = randcl(seed++);

//if (0==1)
	{
		float           h1, s1, l1;
		float           rv;
		float           rr, gg, bb;
		float           hue1;

		hue1 = pt->baked_tex[HUE_VARIATION] / 100.0f;
		hue1 *= .15;
		rr = rc.x;
		gg = rc.y;
		bb = rc.z;
		// this is a rand hue with same val/sat
		// blend it into the original color
		h->rootcolor.x = ( h->rootcolor.x ) * ( 1.0 - ( hue1 ) ) + ( rr * ( hue1 ) );
		h->rootcolor.y = ( h->rootcolor.y ) * ( 1.0 - ( hue1 ) ) + ( gg * ( hue1 ) );
		h->rootcolor.z = ( h->rootcolor.z ) * ( 1.0 - ( hue1 ) ) + ( bb * ( hue1 ) );
		rv = ( pt->baked_tex[HUE_VARIATION] / 100.0f ) ;
		h->rootcolor.x = h->rootcolor.x*(1.0-rv) + h->rootcolor.x * rv1 * rv;
		h->rootcolor.y = h->rootcolor.y*(1.0-rv) + h->rootcolor.y * rv1 * rv;
		h->rootcolor.z = h->rootcolor.z*(1.0-rv) + h->rootcolor.z * rv1 * rv;

		h->tipcolor.x = ( h->tipcolor.x ) * ( 1.0 - ( hue1 ) ) + ( rr * ( hue1 ) );
		h->tipcolor.y = ( h->tipcolor.y ) * ( 1.0 - ( hue1 ) ) + ( gg * ( hue1 ) );
		h->tipcolor.z = ( h->tipcolor.z ) * ( 1.0 - ( hue1 ) ) + ( bb * ( hue1 ) );
		rv = ( pt->baked_tex[VALUE_VARIATION] / 100.0f ) ;
		h->tipcolor.x = h->tipcolor.x*(1.0-rv) + h->tipcolor.x * rv1 * rv;
		h->tipcolor.y = h->tipcolor.y*(1.0-rv) + h->tipcolor.y * rv1 * rv;
		h->tipcolor.z = h->tipcolor.z*(1.0-rv) + h->tipcolor.z * rv1 * rv;
		if (h->tipcolor.x>1.0f) 1.0f;
		if (h->tipcolor.y>1.0f) 1.0f;
		if (h->tipcolor.z>1.0f) 1.0f;
		if (h->rootcolor.x>1.0f) 1.0f;
		if (h->rootcolor.y>1.0f) 1.0f;
		if (h->rootcolor.z>1.0f) 1.0f;

	}


	
}