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

void make_basematrix(BASEHAIR *h,BASEMATRIX *bm)
{
int qv;
								VERT xv,yv,zv;
								VERT tangent[15];

					for( qv = 1; qv < 14; qv++ )
							{
								BASEHAIR       *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;
								yv=Vnorm(yv);
								yv2=Vnorm(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;
								yv=Vnorm(yv);
								tangent[qv]=yv;
							}
							tangent[0]=tangent[1];
							tangent[14]=tangent[13];



	
					for( qv = 0; qv < 15; qv++ )
							{
								BASEHAIR       *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;

								{
									xv = Vnorm( xv );
									yv = Vnorm( yv );
									zv = Vcross( xv, yv );
									zv = Vnorm( zv );
									xv = Vcross( yv, zv );
									xv = Vnorm( xv );

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



							}

}








void twist_basehair (BASEHAIR *clumpCenter, BASEHAIR *h,GLOBS *gs)
{ // rotate the matrix using clumprot
BASEMATRIX bm;
float rot_dir,offset_dir;
int qv;
//unsigned long lastSeed;
float randomize;
//lastSeed=gs->lastSeed;
MTJsrand(clumpCenter->index*5000,gs);
rot_dir=(MTdrand98( gs )-0.5f);
rot_dir=(float)SGN(rot_dir);

//rot_dir= .5;

//offset_dir=(MTdrand98( gs )-0.5f);
//offset_dir=(float)SGN(offset_dir);
offset_dir=1.0f;
randomize=clumpCenter->randomize*(MTdrand98( gs ))+(1.0f-clumpCenter->randomize);


//offset_dir = .5;
//fprintf (stdout,"clump = %d\n",clumpCenter->index);fflush(stdout);
//rot_dir=0.1f;
make_basematrix(clumpCenter,&bm);
for (qv=0;qv<15;qv++)
		{
			float qq;
			Matrix          rot;
			VERT xv,yv,zv,trans;
			float uu1;
			float uuu;
			float flat;
			uuu=(float)((float)qv/14.0f);
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]=vxm3x3(bm.world2zero[qv],h->hv[qv]);
			h->hv[qv].x+=offset_dir*h->clump_rot_offset*uuu;	// back here
			qq = rot_dir*h->clump_rot_strength*uuu;
			qq *= randomize;
			flat=1.0-h->clump_flatness;
			h->hv[qv].z*= flat*uu1+(1.0f-uu1); // flatness
roty(rot,qq);
					h->hv[qv] = vxm3x3( rot, h->hv[qv] ); // spin it

			h->hv[qv]=vxm3x3(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;

		}
//MTJsrand(lastSeed,gs);

}









static void RSxforminstMT(WFTYPE *ret,int ind, GLOBS *gs )
{
	int             x = 0;
	Matrix          pm;
	BASEMATRIX bm;
//BASEHAIR h;
	float           scal = 0.0f;
	int             l1 = 0, l2 = 0;
	VERT           *col, *vel;
	VERT            offs;
	float           sc = 0.0f;
	int             cc = 0;
	BASEHAIR h;
	WFTYPE tmp1,tmp2;
	CURVEINFO cinfo;
init_geomWF(ret);
init_geomWF(&tmp1);
init_geomWF(&tmp2);
ret->totalverts=gs->rs.freeze.totalverts;
ret->totalfaces=gs->rs.freeze.totalfaces;
ret->totalfverts=gs->rs.freeze.totalfverts;
alloc_geomWF(ret);
{
	VERT zero;
	zero.x=0.0f;
	zero.y=0.0f;
	zero.z=0.0f;
MTdraw_lotsWFRS(gs->Gslg,ind,zero,0,&h,&cinfo,gs);
}

 make_basematrix(&h,&bm);
init_geomWF(&tmp1);
init_geomWF(&tmp2);
 copy_geomWF(&tmp1,ret);
 copy_geomWF(&tmp2,ret);
// transform tmp1
 for (x=0;x<tmp1.totalverts;x++)
 {
	float u;
	int iu;
	float remainder,iremainder;
	u=tmp1.v[x].y*15.0f;
	iu = (int) floor((double)u);
	tmp1.v[x].y-=(float)iu;
	tmp1.v[x]=vxm(bm.zero2world[iu],tmp1.v[x]);
	tmp1.v[x].y+=(float)iu;
	remainder=u-(float)iu;
	iremainder=1.0f-remainder;
	iu = (int) ceil((double)u);
	tmp2.v[x].y-=(float)iu;
	tmp2.v[x]=vxm(bm.zero2world[iu],tmp2.v[x]);
	tmp2.v[x].y+=(float)iu;
	ret->v[x].x=tmp2.v[x].x*remainder+tmp1.v[x].x*iremainder;
	ret->v[x].y=tmp2.v[x].y*remainder+tmp1.v[x].y*iremainder;
	ret->v[x].z=tmp2.v[x].z*remainder+tmp1.v[x].z*iremainder;
 }
 free_geomWF(&tmp1);
 free_geomWF(&tmp2);
}



static void xforminstMT(WFTYPE *ret,int ind, GLOBS *gs )
{
	int             x = 0;
	Matrix          pm;
	BASEMATRIX bm;
	float ybound;
//	FILE *fpp;

//BASEHAIR h;
	float           scal = 0.0f;
	int             l1 = 0, l2 = 0;
	VERT           *col, *vel;
	VERT            offs;
	float           sc = 0.0f;
	int             cc = 0;
	BASEHAIR h;
	WFTYPE tmp1,tmp2;
	CURVEINFO cinfo;
//	fpp=fopen("c:\\tmp.txt","a");
//	fprintf(fpp,"making an instance\n");fflush(fpp);
init_geomWF(ret);
init_geomWF(&tmp1);
init_geomWF(&tmp2);
ret->totalverts=freeze.totalverts+100;
ret->totalfaces=freeze.totalfaces+100;
ret->totalfverts=freeze.totalfverts+100;
alloc_geomWF(ret);
ret->totalverts=freeze.totalverts;
ret->totalfaces=freeze.totalfaces;
ret->totalfverts=freeze.totalfverts;

for (x=0;x<freeze.totalverts;x++)
{
ret->v[x]=freeze.v[x];
}
for (x=0;x<freeze.totalfverts;x++)
{
ret->facelist[x]=freeze.facelist[x];
}
for (x=0;x<freeze.totalfaces;x++)
{
ret->face_start[x]=freeze.face_start[x];
ret->face_end[x]=freeze.face_end[x];
}

{
	VERT zero;
	zero.x=0.0f;
	zero.y=0.0f;
	zero.z=0.0f;
//	fprintf(fpp,"making an instance2\n");fflush(fpp);
MTdraw_lotsWF(gs->Gslg,ind,zero,0,&h,&cinfo,gs);

//	fprintf(fpp,"making an instance3\n");fflush(fpp);
}
for (x=0;x<freeze.totalverts;x++)
{
ret->v[x]=freeze.v[x];
//ret->color[x]=freeze.colorlock[x];
}

for (x=0;x<freeze.totalfverts;x++)
ret->facelist[x]=freeze.facelist[x];

for (x=0;x<freeze.totalfaces;x++)
{
ret->face_start[x]=freeze.face_start[x];
ret->face_end[x]=freeze.face_end[x];
ret->index[x]=h.index;
}
 make_basematrix(&h,&bm);
init_geomWF(&tmp1);
init_geomWF(&tmp2);
 copy_geomWF(&tmp1,ret);
 copy_geomWF(&tmp2,ret);
// transform tmp1
//	fprintf(fpp,"making an instance4\n");fflush(fpp);

 ybound= -1.0f;
for (x=0;x<freeze.totalverts;x++)
{
if (freeze.v[x].y>ybound) ybound=freeze.v[x].y;
}

if (ybound>0.0f)
for (x=0;x<freeze.totalverts;x++)
{
	float scal;
	float w;
//	scal = ( h.restlength * 14.0f ) / ybound;
	scal = ( 1.0f ) / ybound;
	w = scal*h.thickness*(0.3f);

tmp1.v[x].x*=w;
tmp2.v[x].x*=w;

tmp1.v[x].y*=scal;
tmp2.v[x].y*=scal;

tmp1.v[x].z*=w;
tmp2.v[x].z*=w;


//freeze.v[x].y*=14.0f;
if (tmp1.v[x].y<0.0f) tmp1.v[x].y=0.0f;
if (tmp2.v[x].y<0.0f) tmp2.v[x].y=0.0f;
}
if (ybound>0.0f)
 for (x=0;x<tmp1.totalverts;x++)
 {
	float u;
	int iu;
	int i2;
	float remainder,iremainder;
	u=(tmp1.v[x].y);
//	if (u>14.0f) u=14.0f;
	iu=(int)floor( u*14.0f);
	i2=iu+1;
	if (i2>14) i2=14;
	remainder=u*14.0f-(float)iu;
//	remainder/=14.0f;
	if (remainder>1.0f) remainder=1.0f;
	iremainder=1.0f-remainder;
//	if (iu>14) iu=14;

//	tmp1.v[x].x*=h.restlength*h.thickness*.3*(1.0/ybound);
//	tmp1.v[x].z*=h.restlength*h.thickness*.3*(1.0/ybound);
	tmp1.v[x].y-=(float)u;
	tmp1.v[x]=vxm(bm.zero2world[iu],tmp1.v[x]);
//	tmp1.v[x].y+=(float)iu;
	tmp1.v[x].x+=(float)h.hv[iu].x*iremainder+h.hv[i2].x*remainder;
	tmp1.v[x].y+=(float)h.hv[iu].y*iremainder+h.hv[i2].y*remainder;
	tmp1.v[x].z+=(float)h.hv[iu].z*iremainder+h.hv[i2].z*remainder;
	tmp1.velocity[x].x+=(float)h.velocity[iu].x*iremainder+h.velocity[i2].x*remainder;
	tmp1.velocity[x].y+=(float)h.velocity[iu].y*iremainder+h.velocity[i2].y*remainder;
	tmp1.velocity[x].z+=(float)h.velocity[iu].z*iremainder+h.velocity[i2].z*remainder;

	tmp1.color[x]=h.color[0];
//	iu = (int) ceil((double)u);
//	tmp2.v[x].y-=(float)iu;
//	iu=(int)(  ((float)freeze.envelope[x].link2/52.0f) * 14.0f );
//	iu++;
//	if (iu>14) iu=14;
//	tmp2.v[x].x*=h.restlength*h.thickness;
//	tmp2.v[x].z*=h.restlength*h.thickness;
//	tmp2.v[x].x*=h.restlength*h.thickness*.3f*(1.0f/ybound);
//	tmp2.v[x].z*=h.restlength*h.thickness*.3f*(1.0f/ybound);

	tmp2.v[x].y-=(float)u;
	tmp2.v[x]=vxm(bm.zero2world[i2],tmp2.v[x]);
	tmp2.v[x].x+=(float)h.hv[iu].x*iremainder+h.hv[i2].x*remainder;
	tmp2.v[x].y+=(float)h.hv[iu].y*iremainder+h.hv[i2].y*remainder;
	tmp2.v[x].z+=(float)h.hv[iu].z*iremainder+h.hv[i2].z*remainder;
	tmp2.velocity[x].x+=(float)h.velocity[iu].x*iremainder+h.velocity[i2].x*remainder;
	tmp2.velocity[x].y+=(float)h.velocity[iu].y*iremainder+h.velocity[i2].y*remainder;
	tmp2.velocity[x].z+=(float)h.velocity[iu].z*iremainder+h.velocity[i2].z*remainder;
	tmp2.color[x]=h.color[0];
	ret->v[x].x=tmp2.v[x].x*remainder+tmp1.v[x].x*iremainder;
	ret->v[x].y=tmp2.v[x].y*remainder+tmp1.v[x].y*iremainder;
	ret->v[x].z=tmp2.v[x].z*remainder+tmp1.v[x].z*iremainder;
	ret->velocity[x].x=tmp2.velocity[x].x*remainder+tmp1.velocity[x].x*iremainder;
	ret->velocity[x].y=tmp2.velocity[x].y*remainder+tmp1.velocity[x].y*iremainder;
	ret->velocity[x].z=tmp2.velocity[x].z*remainder+tmp1.velocity[x].z*iremainder;
	ret->color[x]=h.color[0];
	ret->v[x].z*= -1.0f;
	ret->velocity[x].z*= -1;

//	fprintf (fpp,"ret->v[%d] = %f %f %f\n",x,ret->v[x].x,ret->v[x].y,ret->v[x].z);fflush(fpp);
 }
 free_geomWF(&tmp1);
 free_geomWF(&tmp2);
//	fprintf(fpp,"making an instance5\n");fflush(fpp);
//	fclose(fpp);
}