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

static void     Sblur_position( int x );
void            Sdodyn3( void );
void            Sdodyn3_single( int index );
void            dodyn3_single( int index );

void MTSdodyn3( int x);

void constrain_volume(void);

void MTdodyn_kernal ( void);

static void
blur_velocity( int x )
{
	int             q;
	BASEHAIR        tmp;
	float           hit[15];
	int             y;

	if( hair[x].restlength > 0.0001 )
		for( y = 0; y < 15; y++ )
		{
			tmp.hv[y].x = 0;
			tmp.hv[y].y = 0;
			tmp.hv[y].z = 0;
			hit[y] = 0;
		}

	if( hair[x].restlength > 0.0001 )
		for( y = 0; y < 15; y++ )
		{
			for( q = y - 1; q <= y + 1; q++ )
			{
				int             q1;

				q1 = q;
				if( q1 < 0 )
					q1 = 0;
				if( q1 > 15 - 1 )
					q1 = 15 - 1;
				tmp.hv[q1].x += hair[x].velocity[y].x;
				tmp.hv[q1].y += hair[x].velocity[y].y;
				tmp.hv[q1].z += hair[x].velocity[y].z;
				hit[q1] += 1.0f;
			}
		}
	if( hair[x].restlength > 0.0001 )
		for( y = 0; y < 15; y++ )
		{
			if( hit[y] != 0 )
			{
				tmp.hv[y].x /= hit[y];
				tmp.hv[y].y /= hit[y];
				tmp.hv[y].z /= hit[y];
				hair[x].velocity[y] = tmp.hv[y];
			}
		}
}

static void
blur_position( int x )
{
	int             q;
	BASEHAIR        tmp;
	float           hit[15];
	int             y;

	if( hair[x].restlength > 0.0001 )
		for( y = 0; y < 15; y++ )
		{
			tmp.hv[y].x = 0;
			tmp.hv[y].y = 0;
			tmp.hv[y].z = 0;
			hit[y] = 0;
		}

	if( hair[x].restlength > 0.0001 )
		for( y = 0; y < 15; y++ )
		{
			for( q = y - 1; q <= y + 1; q++ )
			{
				int             q1;

				q1 = q;
				if( q1 < 0 )
					q1 = 0;
				if( q1 > 15 - 1 )
					q1 = 15 - 1;
				if( ( q1 != 0 ) && ( q1 != 14 ) && ( hair[x].pfID[y] == -1 ) )
				{
					tmp.hv[q1].x += hair[x].hv[y].x;
					tmp.hv[q1].y += hair[x].hv[y].y;
					tmp.hv[q1].z += hair[x].hv[y].z;
					hit[q1] += 1.0f;
				}
				else
				{
					tmp.hv[y].x += hair[x].hv[y].x;
					tmp.hv[y].y += hair[x].hv[y].y;
					tmp.hv[y].z += hair[x].hv[y].z;
					hit[y] += 1.0f;
				}
			}
		}
	if( hair[x].restlength > 0.0001 )
		for( y = 0; y < 15; y++ )
		{
			if( hit[y] != 0 )
			{
				tmp.hv[y].x /= hit[y];
				tmp.hv[y].y /= hit[y];
				tmp.hv[y].z /= hit[y];
				hair[x].hv[y] = tmp.hv[y];
			}
		}



}


static void
Sblur_position( int x )
{
	int             q;
	BASEHAIR        tmp;
	float           hit[15];
	int             y;

	if( Shair[x].restlength > 0.0001 )
		for( y = 0; y < 15; y++ )
		{
			tmp.hv[y].x = 0;
			tmp.hv[y].y = 0;
			tmp.hv[y].z = 0;
			hit[y] = 0;
		}

	if( Shair[x].restlength > 0.0001 )
		for( y = 0; y < 15; y++ )
		{
			for( q = y - 1; q <= y + 1; q++ )
			{
				int             q1;

				q1 = q;
				if( q1 < 0 )
					q1 = 0;
				if( q1 > 15 - 1 )
					q1 = 15 - 1;
				if( ( q1 != 0 ) && ( q1 != 14 ) && ( Shair[x].pfID[y] == -1 ) )
				{
					tmp.hv[q1].x += Shair[x].hv[y].x;
					tmp.hv[q1].y += Shair[x].hv[y].y;
					tmp.hv[q1].z += Shair[x].hv[y].z;
					hit[q1] += 1.0f;
				}
				else
				{
					tmp.hv[y].x += Shair[x].hv[y].x;
					tmp.hv[y].y += Shair[x].hv[y].y;
					tmp.hv[y].z += Shair[x].hv[y].z;
					hit[y] += 1.0f;
				}
			}
		}
	if( Shair[x].restlength > 0.0001 )
		for( y = 0; y < 15; y++ )
		{
			if( hit[y] != 0 )
			{
				tmp.hv[y].x /= hit[y];
				tmp.hv[y].y /= hit[y];
				tmp.hv[y].z /= hit[y];
				Shair[x].hv[y] = tmp.hv[y];
			}
		}



}
void MTdodyn(int x)
	{
	int nn;

//				if (x<totalguides)
				{
	if (total_splines>0) nn=total_splines;
	if (total_splines==0) nn=totalguides;
	if (total_splines==0)
	dodyn3_single( x );
//fprintf (stdout,"MTdodyn %d\n",x);fflush(stdout);

	if (total_splines>0)
	Sdodyn3_single( x );
if (Gslg==0)
		if( hair[x].restlength > 0.0001 )
		{
			blur_velocity( x );
			blur_position( x );
		}

		if (total_splines>0)
		if( Shair[x].restlength > 0.0001 )
		{
			Sblur_position( x );

		}

				}

	}


static void     dodyn3( void );
int             last_collider = -90;
void
dodyn( void )
{
	int             t;
	int             xx, x, y;
MTdodyn_kernal ();

}

void
dodyn3( void )
{

	int             x, y, qq;
	float           ytable[15], ytable3[15], iytable4[15];

//  BASEHAIR lh;
	float           stiff = ( float ) .1;
	float           maxrest[5];
	float           minrest[5];

	Jsrand( 0 );
	for( x = 0; x < 15; x++ )
		ytable[x] = ( float ) x / 15.0f;
	for( x = 0; x < 15; x++ )
		ytable3[x] = ytable[x] * ytable[x] * ytable[x];
	for( x = 0; x < 15; x++ )
		iytable4[x] = ( 1.0f - ytable[x] ) * ( 1.0f - ytable[x] ) * ( 1.0f - ytable[x] ) * ( 1.0f - ytable[x] );
	for( x = 0; x < totalverts; x++ )
		if( hair[x].restlength > 0.0001 )
			hair[x].restlength /= ( float ) restBOUNDLENGTH;
	for( x = 0; x < totalverts; x++ )
		if( hair[x].restlength > 0.0001 )
						{
							BASEHAIR       *hhh;

							hhh = &hair[x];
							for( qq = 0; qq < 15; qq++ )
							{
								hhh->hv[qq].x /= ( float ) restBOUNDLENGTH;
								hhh->hv[qq].y /= ( float ) restBOUNDLENGTH;
								hhh->hv[qq].z /= ( float ) restBOUNDLENGTH;
// need these for surflock !
								hhh->resthv[qq].x /= ( float ) restBOUNDLENGTH;
								hhh->resthv[qq].y /= ( float ) restBOUNDLENGTH;
								hhh->resthv[qq].z /= ( float ) restBOUNDLENGTH;

//      thishair[x].hv[qq].x/=(float)restBOUNDLENGTH;
//      thishair[x].hv[qq].y/=(float)restBOUNDLENGTH;
//      thishhh->hv[qq].z/=(float)restBOUNDLENGTH;
//                              hhh->noisev[qq].
//                                  x /= ( float )restBOUNDLENGTH;
//                              hhh->noisev[qq].
//                                  y /= ( float )restBOUNDLENGTH;
//                              hhh->noisev[qq].
//                                  z /= ( float )restBOUNDLENGTH;



								hhh->lasthv[qq].x /= ( float ) restBOUNDLENGTH;
								hhh->lasthv[qq].y /= ( float ) restBOUNDLENGTH;
								hhh->lasthv[qq].z /= ( float ) restBOUNDLENGTH;
//                              hhh->resthv[qq].
//                                  x /= ( float )restBOUNDLENGTH;
//                              hhh->resthv[qq].
//                                  y /= ( float )restBOUNDLENGTH;
//                              hhh->resthv[qq].
//                                  z /= ( float )restBOUNDLENGTH;
//  hhh->restlength/=(float)BOUNDLENGTH; 
							}
						}
#ifdef LIB
//	make_normals(  );
#endif
	maxrest[0] = 0;
	maxrest[1] = 0;
	maxrest[2] = 0;
	maxrest[3] = 0;
	maxrest[4] = 0;
	minrest[0] = 1110;
	minrest[1] = 1110;
	minrest[2] = 1110;
	minrest[3] = 1110;
	minrest[4] = 1110;
	for( x = 0; x < totalverts; x++ )
	{
		BASEHAIR       *hhh;

		hhh = &hair[x];
		if( hhh->mtl == head )
			if( hhh->restlength > maxrest[0] )
				maxrest[0] = hhh->restlength;

		if( hhh->mtl == splines )
			if( hhh->restlength > maxrest[4] )
				maxrest[4] = hhh->restlength;
		if( hhh->mtl == head )
			if( hhh->restlength < minrest[0] )
				minrest[0] = hhh->restlength;
		if( hhh->mtl == splines )
			if( hhh->restlength < minrest[4] )
				minrest[4] = hhh->restlength;
	}

//fprintf (stdout,"stiff slider %f  stiff vmap %f\n",sliders[8][0].value,hair[0].slider[8]);fflush(stdout);

	for( x = 0; x < totalverts; x++ )
	{
		BASEHAIR       *hhh;


		hhh = &hair[x];
		if( ( hhh->mtl == head ) || ( hhh->mtl == beard ) || ( hhh->mtl == eyelash ) || ( hhh->mtl == eyebrow ) || ( hhh->mtl == splines ) )
			if( hhh->restlength > 0.0001 )
						{
							int             slg = 0;
							BASEHAIR        rr;
							BASEHAIR        tmp;
							float           tmpstiff;
							float           stiff3 = 1.0f;
							float           root_stiffness, dampening_force;
							float           idampening_force;

							root_stiffness = sliders[21][0].value * hhh->slider[21];
							dampening_force = sliders[40][0].value * hhh->slider[40];
							idampening_force = 1.0f - dampening_force;
							if( head >= 0 )
								if( hhh->mtl == head )
								{
									slg = 0;
									stiff3 = sliders[8][0].value * hhh->slider[8];
								}
							if( splines >= 0 )
								if( hhh->mtl == splines )
								{
									slg = 4;
									stiff3 = sliders[8][4].value * hhh->slider[8];
								}
							tmpstiff = stiff3;
							if( hhh->restlength > 0.0001 )
								if( sliders[22][slg].value > 0 )
									hair[x] = displace_wiggle( hair[x], 82.0 * sliders[22][slg].value, sliders[23][slg].value );

							if( maxrest[slg] > 0 )
							{
								float           stf;

//              if ((minrest[slg]/maxrest[slg])<.7)
//  stf=(1.0f-stiff3)*(((hhh->restlength-minrest[slg])/(maxrest[slg]-minrest[slg])));
//  stiff3+=(1.0f-stiff3)*((hhh->restlength-0)/(maxrest[slg]-0));
							}
							if( stiff3 > 1.0f )
								stiff3 = 1.0f;
							if( stiff3 < .999 )
							{
								int             yy;
								float           dxx, dyy, dzz;
								float           ys;	// this is y from 0 to 1;
								float           iys;	// this is y from 1 to 0;

								stiff = stiff3;

// randomize stiffness								stiff = stiff * ( drand98(  ) * .6 - .3 ) + stiff;
								if( stiff < 0 )
									stiff = 0;
								if( stiff > 1 )
									stiff = 1;
								memcpy( &rr, &hair[x], sizeof( BASEHAIR ) );
								hhh->velocity[0].x = 0;
								hhh->velocity[0].y = 0;
								hhh->velocity[0].z = 0;
								hhh->velocity[1].x = 0;
								hhh->velocity[1].y = 0;
								hhh->velocity[1].z = 0;
//   tmp.hv[0]=hhh->hv[0];
//								for( yy = 0; yy < 15; yy++ )
//								{
//									hhh->velocity[yy].x *= idampening_force;
//									hhh->velocity[yy].y *= idampening_force;
//									hhh->velocity[yy].z *= idampening_force;
//								}


//lh=hair[x];     

//   tmp.hv[0]=hhh->hv[1];
								dxx = hhh->hv[0].x - hhh->lasthv[0].x;
								dyy = hhh->hv[0].y - hhh->lasthv[0].y;
								dzz = hhh->hv[0].z - hhh->lasthv[0].z;



								tmp = hair[x];	// back up wannabe position
								{
									float           stf2;

									stf2 = stiff * .8;
									for( y = 1; y < 15; y++ )
									{
//          VERT qqq;
										float           pw;
										float           ipw;

										ys = y/14.0f;
										iys = 1.0f - ys;
										pw = ytable3[y];
										ipw = iytable4[y];


										// this is the normalizing strength
										{
											float           damp;
											float           damp2;
damp=ipw;

											damp = iys*stiff3;	// root
											damp+= stiff3*stiff3;
											if (damp>1.0f) damp=1.0f;

rr.hv[y]=hhh->lasthv[y];
											rr.hv[y].x += ( hhh->hv[y].x - hhh->lasthv[y].x ) * damp;
											rr.hv[y].y += ( hhh->hv[y].y - hhh->lasthv[y].y ) * damp;
											rr.hv[y].z += ( hhh->hv[y].z - hhh->lasthv[y].z ) * damp;
											
											damp = ipw*ipw*root_stiffness;	// root
											rr.hv[y].x += ( hhh->hv[y].x - rr.hv[y].x ) * damp;
											rr.hv[y].y += ( hhh->hv[y].y - rr.hv[y].y ) * damp;
											rr.hv[y].z += ( hhh->hv[y].z - rr.hv[y].z ) * damp;

										}


										{


										hhh->hv[y]=rr.hv[y];
										}

    hhh->hv[y].x+=hhh->velocity[y].x*(1.0-stiff3); // move with velocity from last frame
    hhh->hv[y].y+=hhh->velocity[y].y*(1.0-stiff3);
    hhh->hv[y].z+=hhh->velocity[y].z*(1.0-stiff3);
									hhh->velocity[y].x+=(rr.hv[y].x-hhh->lasthv[y].x)*(1.0-stiff3); // accumulate velocity
									hhh->velocity[y].y+=(rr.hv[y].y-hhh->lasthv[y].y)*(1.0-stiff3); // from this one
									hhh->velocity[y].z+=(rr.hv[y].z-hhh->lasthv[y].z)*(1.0-stiff3);
									}	//y
								}
//#ifdef RENDERLW
//for (y=0;y<15;y++) hhh->hv[y]=hhh->lasthv[y];

#ifndef NOLIB
//       apply_wind(x);

#endif
								memcpy( &tmp, &hair[x], sizeof( BASEHAIR ) );	// tmp is wannabe + velocity
//                              recalc_hair(x);
								dyn_recalc_hair( x );	// gottabe
								if( Gcollision_hit )
								{
									Gcollision_hit = 0;
									dyn_recalc_hair( x );	// gottabe
								}
								if( Gcollision_hit )
								{
									Gcollision_hit = 0;
									dyn_recalc_hair( x );	// gottabe
								}
								for( y = 1; y < 15; y++ )	// sum in secondary motion
								{
//         float stiff2;
//         float stiff5;
									float           damp;

//used?                                 float           pw, ipw, pww;

//         int ii;
									ys = ( float ) ytable[y];
									iys = 1.0f - ys;
//          slow                        pw = pow( ys, 2.0f );
//      slow                            ipw = pow( ys, 2.0f );
///         used?                       pw=ys*ys;
////        used?                           ipw=pw;
									////    used?                           ipw = 1.0f - ipw;
									//// used?                          pww = 1.0f - pow( 1.0f - ys, 2.5 );
#ifndef NOLIB
#ifdef SOFTIMAGE
									{
										VERT tmpvel;
										float bent;
										bent=pow(ytable[y],3.5);
										hhh->velocity[y].x *= restBOUNDLENGTH;
										hhh->velocity[y].y *= restBOUNDLENGTH;
										hhh->velocity[y].z *= -restBOUNDLENGTH;
										tmpvel.x=0;
										tmpvel.y=0;
										tmpvel.z=0;
										SHAVEcoord_convertFROMSHAVE( &hhh->velocity[y] );

										if( hhh->restlength > 0.0001 )
											MAYAapply_cached_forces( x, y, &tmpvel );
//										if (hhh->vid==34)
//										{
//											printf ("tmpvel = %f %f %f\n",tmpvel.x,tmpvel.y,tmpvel.z);
//											fflush(stdout);
//										}
										tmpvel.x*=bent;
										tmpvel.y*=bent;
										tmpvel.z*=bent;
										hhh->velocity[y].x+=tmpvel.x;
										hhh->velocity[y].y+=tmpvel.y;
										hhh->velocity[y].z+=tmpvel.z;

										SHAVEcoord_convertTOSHAVE( &hhh->velocity[y] );
										hhh->velocity[y].x /= restBOUNDLENGTH;
										hhh->velocity[y].y /= restBOUNDLENGTH;
										hhh->velocity[y].z /= restBOUNDLENGTH;
										hhh->velocity[y].z *= -1;
									}
#endif
#endif
#ifndef NOLIB
									hhh->velocity[y].y -= gravity * .00451872253 * ys* ( 1.0f - stiff3 );	//gravity
#endif //secondary motion



									hhh->velocity[y].x += (  hhh->hv[y].x-tmp.hv[y].x ) * ( 1.0f - stiff3 ) * .9*(idampening_force);	//.15
									hhh->velocity[y].y += (  hhh->hv[y].y-tmp.hv[y].y ) * ( 1.0f - stiff3 ) * .9*(idampening_force);
									hhh->velocity[y].z += (  hhh->hv[y].z-tmp.hv[y].z ) * ( 1.0f - stiff3 ) * .9*(idampening_force);
//     damp=.75+(1.0f-stiff)*.252;
////       damp=.7+(1.0f-stiff)*.252;
//									damp = ( float ) .992;
									//  damp-=(ipw);
//									if( damp < 0 )
//										damp = 0;
//     damp=pww*(1.0f-stiff);
									// damp=1.0f;
//									hhh->velocity[y].x *= damp;	//*(.2/BOUNDLENGTH);// dampening
//									hhh->velocity[y].y *= damp;	//*(.2/BOUNDLENGTH);
//									hhh->velocity[y].z *= damp;	//*(.2/BOUNDLENGTH);

									hhh->velocity[y].x *= idampening_force*.99;	//*(.2/BOUNDLENGTH);// dampening
									hhh->velocity[y].y *= idampening_force*.99;	//*(.2/BOUNDLENGTH);
									hhh->velocity[y].z *= idampening_force*.99;	//*(.2/BOUNDLENGTH);

									
//										hhh->hv[y].x+=hhh->velocity[y].x;
//										hhh->hv[y].y+=hhh->velocity[y].y;
//										hhh->hv[y].z+=hhh->velocity[y].z;


//      for (ii=0;ii<4;ii++)
//      {
//      pww=pow((1.0f-ys),6.0f);
//      if (pww>1.0f) pww=1.0f;
//      damp=pww;
									//         damp=ipw*ipw*ipw; // normalize
//         if (damp>1) damp=1;
//  
//     hhh->hv[y].x=rr.hv[y].x*damp+hhh->hv[y].x*(1.0f-damp);
//     hhh->hv[y].y=rr.hv[y].y*damp+hhh->hv[y].y*(1.0f-damp);
//     hhh->hv[y].z=rr.hv[y].z*damp+hhh->hv[y].z*(1.0f-damp);
//      }




								}	//y


							}	// stiff < .99
						}		//x
	}
	//Sdodyn3(  );



	for( x = 0; x < totalverts; x++ )
	{
		BASEHAIR       *hhh;

		hhh = &hair[x];


					if( hhh->restlength > 0.0001 )


						for( qq = 0; qq < 15; qq++ )
						{
							hhh->hv[qq].x *= ( float ) restBOUNDLENGTH;
							hhh->hv[qq].y *= ( float ) restBOUNDLENGTH;
							hhh->hv[qq].z *= ( float ) restBOUNDLENGTH;
							hhh->resthv[qq].x *= ( float ) restBOUNDLENGTH;
							hhh->resthv[qq].y *= ( float ) restBOUNDLENGTH;
							hhh->resthv[qq].z *= ( float ) restBOUNDLENGTH;

							hhh->lasthv[qq].x *= ( float ) restBOUNDLENGTH;
							hhh->lasthv[qq].y *= ( float ) restBOUNDLENGTH;
							hhh->lasthv[qq].z *= ( float ) restBOUNDLENGTH;
//                          hhh->resthv[qq].x *= ( float )restBOUNDLENGTH;
//                          hhh->resthv[qq].y *= ( float )restBOUNDLENGTH;
//                          hhh->resthv[qq].z *= ( float )restBOUNDLENGTH;
//                          hhh->noisev[qq].x *= ( float )restBOUNDLENGTH;
//                          hhh->noisev[qq].y *= ( float )restBOUNDLENGTH;
//                          hhh->noisev[qq].z *= ( float )restBOUNDLENGTH;

//      thishhh->hv[qq].x*=(float)restBOUNDLENGTH;
//      thishhh->hv[qq].y*=(float)restBOUNDLENGTH;
//      thishhh->hv[qq].z*=(float)restBOUNDLENGTH;

//      lasthhh->hv[qq].x*=(float)restBOUNDLENGTH;
//      lasthhh->hv[qq].y*=(float)restBOUNDLENGTH;
//      lasthhh->hv[qq].z*=(float)restBOUNDLENGTH;

//      resthhh->hv[qq].x*=(float)restBOUNDLENGTH;
//      resthhh->hv[qq].y*=(float)restBOUNDLENGTH;
//      resthhh->hv[qq].z*=(float)restBOUNDLENGTH;
						}
	}
	for( x = 0; x < totalverts; x++ )
		if( hair[x].restlength > 0.0001 )
			hair[x].restlength *= restBOUNDLENGTH;
//#ifdef NOLIB


}								// end








void
Sdodyn3( void )
{

	int             x, y, qq;
	BASEHAIR        lh;
	float           stiff = ( float ) .1;
	float           maxrest[5];
	float           minrest[5];

	for( x = 0; x < total_splines; x++ )
		if( Shair[x].restlength > 0.0001 )
		{
			Shair[x].restlength = Distance( Shair[x].hv[1], Shair[x].hv[2] );
			Sresthair[x].restlength = Shair[x].restlength;
			Sthishair[x].restlength = Shair[x].restlength;
			Slasthair[x].restlength = Shair[x].restlength;

			Shair[x].restlength /= ( float ) restBOUNDLENGTH;
		}
	for( x = 0; x < total_splines; x++ )
		if( Shair[x].restlength > 0.0001 )
				{

					for( qq = 0; qq < 15; qq++ )
					{
						Shair[x].hv[qq].x /= ( float ) restBOUNDLENGTH;
						Shair[x].hv[qq].y /= ( float ) restBOUNDLENGTH;
						Shair[x].hv[qq].z /= ( float ) restBOUNDLENGTH;
						Sthishair[x].hv[qq].x /= ( float ) restBOUNDLENGTH;
						Sthishair[x].hv[qq].y /= ( float ) restBOUNDLENGTH;
						Sthishair[x].hv[qq].z /= ( float ) restBOUNDLENGTH;

						Slasthair[x].hv[qq].x /= ( float ) restBOUNDLENGTH;
						Slasthair[x].hv[qq].y /= ( float ) restBOUNDLENGTH;
						Slasthair[x].hv[qq].z /= ( float ) restBOUNDLENGTH;
						Sresthair[x].hv[qq].x /= ( float ) restBOUNDLENGTH;
						Sresthair[x].hv[qq].y /= ( float ) restBOUNDLENGTH;
						Sresthair[x].hv[qq].z /= ( float ) restBOUNDLENGTH;
//  hair[x].restlength/(float)BOUNDLENGTH;  
					}
				}

	maxrest[0] = 0;
	maxrest[1] = 0;
	maxrest[2] = 0;
	maxrest[3] = 0;
	maxrest[4] = 0;
	minrest[0] = 1110;
	minrest[1] = 1110;
	minrest[2] = 1110;
	minrest[3] = 1110;
	minrest[4] = 1110;

	for( x = 0; x < total_splines; x++ )
	{
		if( Shair[x].restlength > maxrest[4] )
			maxrest[4] = Shair[x].restlength;
		if( Shair[x].restlength < minrest[4] )
			minrest[4] = Shair[x].restlength;
	}
	for( x = 0; x < total_splines; x++ )
		if( Shair[x].restlength > 0.0001 )
					{
						int             slg;
						BASEHAIR        rr;
						BASEHAIR        tmp;
						float           tmpstiff;
						float           stiff3 = 1.0f;
						float           root_stiffness, dampening_force;

						if( splines >= 0 )
							if( Shair[x].mtl == splines )
							{
								slg = 4;
								stiff3 = sliders[8][4].value * Shair[x].slider[8];
							}
						tmpstiff = stiff3;
						if( Shair[x].restlength > 0.0001 )
							Shair[x] = displace_wiggle( Shair[x], 82.0 * sliders[22][slg].value, sliders[23][slg].value );

						root_stiffness = sliders[21][4].value * Shair[x].slider[21];
						dampening_force = sliders[40][4].value * Shair[x].slider[40];


						if( maxrest[slg] > 0 )
						{
							float           stf;

//              if ((minrest[slg]/maxrest[slg])<.7)
//  stf=(1.0f-stiff3)*(((hair[x].restlength-minrest[slg])/(maxrest[slg]-minrest[slg])));
//                          stf =
//                              ( 1.0f -
//                                stiff3 );/* ( ( ( Shair[x].restlength - 0 ) /
							//       ( maxrest[slg] - 0 ) ) ); */
//                          stf = 1.0f - stf;
//                          stf *= stf;
//                          stf *= stf;
//                          stiff3 = stf;
//  stiff3+=(1.0f-stiff3)*((hair[x].restlength-0)/(maxrest[slg]-0));
						}
						if( stiff3 > 1.0f )
							stiff3 = 1.0f;
						if( tmpstiff < .99 )
						{
							float           dxx, dyy, dzz;
							float           ys;	// this is y from 0 to 1;
							float           iys;	// this is y from 1 to 0;

							stiff = stiff3-.05;
// randomise stiffness							stiff = stiff * ( drand98(  ) * .6 - .3 ) + stiff;
							if( stiff < 0 )
								stiff = 0;
							if( stiff > 1 )
								stiff = 1;
							rr = Shair[x];
							Shairvelocity[x].hv[0].x = 0;
							Shairvelocity[x].hv[0].y = 0;
							Shairvelocity[x].hv[0].z = 0;
							Shairvelocity[x].hv[1].x = 0;
							Shairvelocity[x].hv[1].y = 0;
							Shairvelocity[x].hv[1].z = 0;
//   tmp.hv[0]=hair[x].hv[0];
							{

								int             yy;

								for( yy = 0; yy < 15; yy++ )
								{
									Shairvelocity[x].hv[yy].x *= ( 1.0f - dampening_force );
									Shairvelocity[x].hv[yy].y *= ( 1.0f - dampening_force );
									Shairvelocity[x].hv[yy].z *= ( 1.0f - dampening_force );
								}
							}

							lh = Slasthair[x];

//   tmp.hv[0]=hair[x].hv[1];
							dxx = Shair[x].hv[0].x - Slasthair[x].hv[0].x;
							dyy = Shair[x].hv[0].y - Slasthair[x].hv[0].y;
							dzz = Shair[x].hv[0].z - Slasthair[x].hv[0].z;
							tmp = Shair[x];	// back up wannabe position

							for( y = 1; y < 15; y++ )
							{
//          VERT qqq;
								float           pw;
								float           ipw;
								float           stf2;

//      stf2=stiff*.4+.25; root
								stf2 = stiff;	// root
								if( stf2 > 1 )
									stf2 = 1;
//stf2=1;
								ys = ( float ) ( y ) / 15.0f;
								iys = 1.0f - ys;
								pw = pow( ys, 2.7 );
								ipw = pow( iys, 4.0f );
								lh.hv[y] = Slasthair[x].hv[y];
								lh.hv[y].x += dxx * stf2 * .8;
								lh.hv[y].y += dyy * stf2 * .8;
								lh.hv[y].z += dzz * stf2 * .8;
								// this is the normalizing strength
								{
									float           damp;

//         float stiff2;
									ipw *= ( 1.0f - stiff );
									ipw += stiff;
									if( ipw > 1 )
										ipw = 1;
									damp = ipw * 0.937;	// primary motion
////           damp=(1.0f-stiff)*0.237; // primary motion
									damp += ( float ) .02;
									damp *= damp;
//         damp=ipw*0.137; // primary motion
									if( damp > 1 )
										damp = 1;
									damp *= root_stiffness;	// root


									Shairvelocity[x].hv[y].x += ( rr.hv[y].x - lh.hv[y].x ) * damp;
									Shairvelocity[x].hv[y].y += ( rr.hv[y].y - lh.hv[y].y ) * damp;
									Shairvelocity[x].hv[y].z += ( rr.hv[y].z - lh.hv[y].z ) * damp;

								}
								Shair[x].hv[y].x = lh.hv[y].x + Shairvelocity[x].hv[y].x;
								Shair[x].hv[y].y = lh.hv[y].y + Shairvelocity[x].hv[y].y;
								Shair[x].hv[y].z = lh.hv[y].z + Shairvelocity[x].hv[y].z;

////       Shair[x].hv[y].x=lh.hv[y].x+Shairvelocity[x].hv[y].x;
////       Shair[x].hv[y].y=lh.hv[y].y+Shairvelocity[x].hv[y].y;
////       Shair[x].hv[y].z=lh.hv[y].z+Shairvelocity[x].hv[y].z;
///    hair[x].hv[y].x=hair[x].hv[y].x+hairvelocity[x].hv[y].x;
///    hair[x].hv[y].y=hair[x].hv[y].y+hairvelocity[x].hv[y].y;
///    hair[x].hv[y].z=hair[x].hv[y].z+hairvelocity[x].hv[y].z;
//         hair[x].hv[y].x-=hairvelocity[x].hv[y].x;
//     hair[x].hv[y].y-=hairvelocity[x].hv[y].y;
//     hair[x].hv[y].z-=hairvelocity[x].hv[y].z;
								{
									int             ii;
									float           damp, pww;
									float           stf;

									stf = stiff + .15;
									if( stf > 1.0f )
										stf = 1.0f;
									for( ii = 0; ii < 1; ii++ )
									{
//      pww=pow((1.0f-(y/15.0f)),4.5);
										pww = ( 1.0f - ( y / 15.0f ) );
										pww += stiff;
										if( pww > 1.0f )
											pww = 1.0f;
										damp = pww * .2;	//+(stiff/18.0f)*pww;
										//         damp=ipw*ipw*ipw; // normalize
										if( damp > 1 )
											damp = 1;
										damp *= root_stiffness;	// root

										Shair[x].hv[y].x = rr.hv[y].x * damp + Shair[x].hv[y].x * ( 1.0f - damp );
										Shair[x].hv[y].y = rr.hv[y].y * damp + Shair[x].hv[y].y * ( 1.0f - damp );
										Shair[x].hv[y].z = rr.hv[y].z * damp + Shair[x].hv[y].z * ( 1.0f - damp );


									}


									damp = stiff;	//*stiff;  // normalize  again
									damp = damp + .05;	// normalize  again
									if( damp > 1 )
										damp = 1;
									Shair[x].hv[y].x = rr.hv[y].x * damp + Shair[x].hv[y].x * ( 1.0f - damp );
									Shair[x].hv[y].y = rr.hv[y].y * damp + Shair[x].hv[y].y * ( 1.0f - damp );
									Shair[x].hv[y].z = rr.hv[y].z * damp + Shair[x].hv[y].z * ( 1.0f - damp );
//     hair[x].hv[y].x= rr.hv[y].x*damp + hair[x].hv[y].x*(1.0f-damp);
//     hair[x].hv[y].y= rr.hv[y].y*damp + hair[x].hv[y].y*(1.0f-damp);
//     hair[x].hv[y].z= rr.hv[y].z*damp + hair[x].hv[y].z*(1.0f-damp);


								}


							}	//y

#ifndef NOLIB
//       Sapply_wind(x);

#endif
							memcpy( &tmp, &Shair[x], sizeof( BASEHAIR ) );	// tmp is wannabe + velocity
							Sdyn_recalc_hair( x );	// gottabe
							Sdyn_recalc_hair( x );	// gottabe
							Sdyn_recalc_hair( x );	// gottabe
							Sdyn_recalc_hair( x );	// gottabe

							for( y = 1; y < 15; y++ )	// sum in secondary motion
							{
//         float stiff2;
								float           stiff5;
								float           damp;
								float           pw, ipw, pww;

//         int ii;
								ys = ( float ) ( y ) / 15.0f;
								iys = 1.0f - ys;
								pw = pow( ys, 2.0f );
								ipw = pow( ys, 2.0f );
								ipw = 1.0f - ipw;
								pww = 1.0f - pow( 1.0f - ys, 2.5 );

#ifdef SOFTIMAGE
								{
									Shairvelocity[x].hv[y].x *= restBOUNDLENGTH;
									Shairvelocity[x].hv[y].y *= restBOUNDLENGTH;
									Shairvelocity[x].hv[y].z *= -restBOUNDLENGTH;
									SHAVEcoord_convertFROMSHAVE( &Shairvelocity[x].hv[y] );
									if( Shair[x].restlength > 0.0001 )
										MAYAapply_cached_forces( x, y, &Shairvelocity[x].hv[y] );

									SHAVEcoord_convertTOSHAVE( &Shairvelocity[x].hv[y] );
									Shairvelocity[x].hv[y].x /= restBOUNDLENGTH;
									Shairvelocity[x].hv[y].y /= restBOUNDLENGTH;
									Shairvelocity[x].hv[y].z /= restBOUNDLENGTH;
									Shairvelocity[x].hv[y].z *= -1.0f;
								}
#endif

#ifndef NOLIB
								Shairvelocity[x].hv[y].y -= gravity * .00251872253;	//gravity
#endif
								//secondary motion   
								Shairvelocity[x].hv[y].x -= ( tmp.hv[y].x - Shair[x].hv[y].x ) * ( 1.0f - stiff ) * .7;	//.15
								Shairvelocity[x].hv[y].y -= ( tmp.hv[y].y - Shair[x].hv[y].y ) * ( 1.0f - stiff ) * .7;
								Shairvelocity[x].hv[y].z -= ( tmp.hv[y].z - Shair[x].hv[y].z ) * ( 1.0f - stiff ) * .7;
//     damp=.75+(1.0f-stiff)*.252;
////       damp=.7+(1.0f-stiff)*.252;
								damp = ( float ) .95;
								//  damp-=(ipw);
								if( damp < 0 )
									damp = 0;
//     damp=pww*(1.0f-stiff);
								// damp=1.0f;

								Shairvelocity[x].hv[y].x *= damp;	//*(.2/BOUNDLENGTH);// dampening
								Shairvelocity[x].hv[y].y *= damp;	//*(.2/BOUNDLENGTH);
								Shairvelocity[x].hv[y].z *= damp;	//*(.2/BOUNDLENGTH);
//      for (ii=0;ii<4;ii++)
//      {
//      pww=pow((1.0f-ys),6.0f);
//      if (pww>1.0f) pww=1.0f;
//      damp=pww;
								//         damp=ipw*ipw*ipw; // normalize
//         if (damp>1) damp=1;
//  
//     hair[x].hv[y].x=rr.hv[y].x*damp+hair[x].hv[y].x*(1.0f-damp);
//     hair[x].hv[y].y=rr.hv[y].y*damp+hair[x].hv[y].y*(1.0f-damp);
//     hair[x].hv[y].z=rr.hv[y].z*damp+hair[x].hv[y].z*(1.0f-damp);
//      }




							}	//y

////            dyn_recalc_hair(x); // gottabe
							for( y = 1; y < 15; y++ )
							{
								int             ii;
								float           pww;
								float           damp;

								if( 0 == 1 )
									for( ii = 0; ii < 1; ii++ )
									{
										ys = ( float ) ( y ) / 15.0f;
										iys = 1.0f - ys;
										pww = pow( ( 1.0f - ys ), 6.0f );

										if( pww > 1.0f )
											pww = 1.0f;
										damp = pww * .008;
										//         damp=ipw*ipw*ipw; // normalize again
										if( damp > 1 )
											damp = 1;

										Shair[x].hv[y].x = rr.hv[y].x * damp + Shair[x].hv[y].x * ( 1.0f - damp );
										Shair[x].hv[y].y = rr.hv[y].y * damp + Shair[x].hv[y].y * ( 1.0f - damp );
										Shair[x].hv[y].z = rr.hv[y].z * damp + Shair[x].hv[y].z * ( 1.0f - damp );
									}
							}
							for( y = 1; y < 15; y++ )
							{
								int             ii;

								for( ii = 0; ii < 4; ii++ )
								{
									float           pww;
									float           damp;

//      pww=pow((1.0f-(y/15.0f)),(1.0f-stiff)*4.0+1);
									pww = ( float ) y / 15.0f;
//      pww*=1.8;
									pww *= ( float ) 3.14;
									if( pww > 3.14 )
										pww = ( float ) 3.14;
									pww = cos( pww ) * ( float ) .5 + ( float ) .5;
//      pww=1.0f-pww;
									if( pww > 1.0f )
										pww = 1.0f;
//      damp=pww*.98;//+(stiff/18.0f)*pww;
									damp = pww;
									{
										float           stf;

										stf = stiff * 1.5;
										if( stf > 1 )
											stf = 1;
										damp *= stf;
									}
									//         damp=ipw*ipw*ipw; // normalize spline
									if( damp > 1 )
										damp = 1;
									{
										float           stf;

										stf = stiff + .5;
										if( stf > 1 )
											stf = 1;
//     hair[x].hv[y].x=rr.hv[y].x*damp+hair[x].hv[y].x*(1.0f-damp);
//     hair[x].hv[y].y=rr.hv[y].y*damp+hair[x].hv[y].y*(1.0f-damp);
//     hair[x].hv[y].z=rr.hv[y].z*damp+hair[x].hv[y].z*(1.0f-damp);
										Shairvelocity[x].hv[y].x += ( rr.hv[y].x - Shair[x].hv[y].x ) * damp * stf;
										Shairvelocity[x].hv[y].y += ( rr.hv[y].y - Shair[x].hv[y].y ) * damp * stf;
										Shairvelocity[x].hv[y].z += ( rr.hv[y].z - Shair[x].hv[y].z ) * damp * stf;
									}
								}


								for( ii = 0; ii < 4; ii++ )
								{
									float           pww;
									float           damp;

//      pww=pow((1.0f-(y/15.0f)),(1.0f-stiff)*4.0+1);
									pww = ( float ) y / 15.0f;
									pww *= ( float ) 1.8;
									pww *= ( float ) 3.14;
									if( pww > ( float ) 3.14 )
										pww = ( float ) 3.14;
									pww = cos( pww ) * ( float ) .5 + ( float ) .5;
//      pww=1.0f-pww;
									if( pww > ( float ) 1.0f )
										pww = ( float ) 1.0f;
//      damp=pww*.98;//+(stiff/18.0f)*pww;
									damp = pww;
									damp *= stiff;
									//         damp=ipw*ipw*ipw; // normalize
									if( damp > 1 )
										damp = 1;
									{
										Shair[x].hv[y].x = rr.hv[y].x * damp + Shair[x].hv[y].x * ( 1.0f - damp );
										Shair[x].hv[y].y = rr.hv[y].y * damp + Shair[x].hv[y].y * ( 1.0f - damp );
										Shair[x].hv[y].z = rr.hv[y].z * damp + Shair[x].hv[y].z * ( 1.0f - damp );
//     hairvelocity[x].hv[y].x+=(rr.hv[y].x-hair[x].hv[y].x)*damp*stf;
//     hairvelocity[x].hv[y].y+=(rr.hv[y].y-hair[x].hv[y].y)*damp*stf;
//     hairvelocity[x].hv[y].z+=(rr.hv[y].z-hair[x].hv[y].z)*damp*stf;
									}
								}
							}

//#ifdef LIB
							//  dyn_recalc_hair(x); // gottabe
//#endif
//      blur_velocity(x);
							//  blur_postion(x);
						}		// stiff < .99

					}			//x

	for( x = 0; x < total_splines; x++ )
					if( Shair[x].restlength > 0.0001 )

						for( qq = 0; qq < 15; qq++ )
						{
							Shair[x].hv[qq].x *= ( float ) restBOUNDLENGTH;
							Shair[x].hv[qq].y *= ( float ) restBOUNDLENGTH;
							Shair[x].hv[qq].z *= ( float ) restBOUNDLENGTH;
							Sthishair[x].hv[qq].x *= ( float ) restBOUNDLENGTH;
							Sthishair[x].hv[qq].y *= ( float ) restBOUNDLENGTH;
							Sthishair[x].hv[qq].z *= ( float ) restBOUNDLENGTH;

							Slasthair[x].hv[qq].x *= ( float ) restBOUNDLENGTH;
							Slasthair[x].hv[qq].y *= ( float ) restBOUNDLENGTH;
							Slasthair[x].hv[qq].z *= ( float ) restBOUNDLENGTH;

							Sresthair[x].hv[qq].x *= ( float ) restBOUNDLENGTH;
							Sresthair[x].hv[qq].y *= ( float ) restBOUNDLENGTH;
							Sresthair[x].hv[qq].z *= ( float ) restBOUNDLENGTH;
						}
	for( x = 0; x < total_splines; x++ )
		if( Shair[x].restlength > 0.0001 )
			Shair[x].restlength *= restBOUNDLENGTH;
//#ifdef NOLIB
	for( x = 0; x < total_splines; x++ )
					if( Shair[x].restlength > 0.0001 )
						Slasthair[x] = Shair[x];
//#endif

}								// end







void constrain_volume(void)
{
	int f,fv,x;
for (f=0;f<totalfaces;f++)
if (face_end[f]-face_start[f]>=3)
	 for (fv=face_start[f];fv<face_end[f];fv++)
	 {
	 int vv;
	 int g,g1;
	 float restd;
	 float dx,dy,dz;
	 float d;
	 int ffv;
		float stiff;
		int y;

	 ffv=fv-1;
	 if (ffv<face_start[f]) ffv=face_end[f]-1;
	 g=facelist[ffv];
	 g1=facelist[fv];
	 stiff = sliders[8][0].value * hair[g1].slider[8];
//	 stiff/=2.0f;
	 stiff=.5;
for (y=1;y<15;y++)
{
	VERT vvv;
	float drd;
	 // measure rest length of seg level 
	   dx=hair[g1].resthv[y].x-hair[g].resthv[y].x;
	   dy=hair[g1].resthv[y].y-hair[g].resthv[y].y;
	   dz=hair[g1].resthv[y].z-hair[g].resthv[y].z;
		restd=sqrt(dx*dx+dy*dy+dz*dz);
	   dx=hair[g1].hv[y].x-hair[g].hv[y].x;
	   dy=hair[g1].hv[y].y-hair[g].hv[y].y;
	   dz=hair[g1].hv[y].z-hair[g].hv[y].z;
vvv.x=dx;
vvv.y=dy;
vvv.z=dz;
vvv=Vnorm(vvv);
		d=sqrt(dx*dx+dy*dy+dz*dz);
		drd=d/restd;
//		if (d<restd)
		{	
			VERT p;
		p.x=hair[g].hv[y].x+vvv.x*restd;	
		p.y=hair[g].hv[y].y+vvv.y*restd;	
		p.z=hair[g].hv[y].z+vvv.z*restd;
		hair[g1].hv[y].x=p.x*stiff+hair[g1].hv[y].x*(1.0-stiff);
		hair[g1].hv[y].y=p.y*stiff+hair[g1].hv[y].y*(1.0-stiff);
		hair[g1].hv[y].z=p.z*stiff+hair[g1].hv[y].z*(1.0-stiff);
		}
}
dyn_recalc_hair(g1);

	 }





}











void
dodyn3_single( int index )
{

	int             x, y, qq;
	float           ytable[15], ytable3[15], iytable4[15];

//  BASEHAIR lh;
	float           stiff = ( float ) .1;
	float           maxrest[5];
	float           minrest[5];

	Jsrand( 0 );
	for( x = 0; x < 15; x++ )
		ytable[x] = ( float ) x / 15.0f;
	for( x = 0; x < 15; x++ )
		ytable3[x] = ytable[x] * ytable[x] * ytable[x];
	for( x = 0; x < 15; x++ )
		iytable4[x] = ( 1.0f - ytable[x] ) * ( 1.0f - ytable[x] ) * ( 1.0f - ytable[x] ) * ( 1.0f - ytable[x] );
//	for( x = 0; x < totalverts; x++ )
		if( hair[index].restlength > 0.0001 )
			hair[index].restlength /= ( float ) restBOUNDLENGTH;

x=index;

//	for( x = 0; x < totalverts; x++ )
		if( hair[x].restlength > 0.0001 )
						{
							BASEHAIR       *hhh;

							hhh = &hair[x];
							for( qq = 0; qq < 15; qq++ )
							{
								hhh->hv[qq].x /= ( float ) restBOUNDLENGTH;
								hhh->hv[qq].y /= ( float ) restBOUNDLENGTH;
								hhh->hv[qq].z /= ( float ) restBOUNDLENGTH;
// need these for surflock !
								hhh->resthv[qq].x /= ( float ) restBOUNDLENGTH;
								hhh->resthv[qq].y /= ( float ) restBOUNDLENGTH;
								hhh->resthv[qq].z /= ( float ) restBOUNDLENGTH;

//      thishair[x].hv[qq].x/=(float)restBOUNDLENGTH;
//      thishair[x].hv[qq].y/=(float)restBOUNDLENGTH;
//      thishhh->hv[qq].z/=(float)restBOUNDLENGTH;
//                              hhh->noisev[qq].
//                                  x /= ( float )restBOUNDLENGTH;
//                              hhh->noisev[qq].
//                                  y /= ( float )restBOUNDLENGTH;
//                              hhh->noisev[qq].
//                                  z /= ( float )restBOUNDLENGTH;



								hhh->lasthv[qq].x /= ( float ) restBOUNDLENGTH;
								hhh->lasthv[qq].y /= ( float ) restBOUNDLENGTH;
								hhh->lasthv[qq].z /= ( float ) restBOUNDLENGTH;
//                              hhh->resthv[qq].
//                                  x /= ( float )restBOUNDLENGTH;
//                              hhh->resthv[qq].
//                                  y /= ( float )restBOUNDLENGTH;
//                              hhh->resthv[qq].
//                                  z /= ( float )restBOUNDLENGTH;
//  hhh->restlength/=(float)BOUNDLENGTH; 
							}
						}


/* pre-requsites ------
	maxrest[0] = 0;
	maxrest[1] = 0;
	maxrest[2] = 0;
	maxrest[3] = 0;
	maxrest[4] = 0;
	minrest[0] = 1110;
	minrest[1] = 1110;
	minrest[2] = 1110;
	minrest[3] = 1110;
	minrest[4] = 1110;
	for( x = 0; x < totalverts; x++ )
	{
		BASEHAIR       *hhh;

		hhh = &hair[x];
		if( hhh->mtl == head )
			if( hhh->restlength > maxrest[0] )
				maxrest[0] = hhh->restlength;

		if( hhh->mtl == splines )
			if( hhh->restlength > maxrest[4] )
				maxrest[4] = hhh->restlength;
		if( hhh->mtl == head )
			if( hhh->restlength < minrest[0] )
				minrest[0] = hhh->restlength;
		if( hhh->mtl == splines )
			if( hhh->restlength < minrest[4] )
				minrest[4] = hhh->restlength;
	}
*/
//fprintf (stdout,"stiff slider %f  stiff vmap %f\n",sliders[8][0].value,hair[0].slider[8]);fflush(stdout);

		x=index;
//	for( x = 0; x < totalverts; x++ )
	{
		BASEHAIR       *hhh;


		hhh = &hair[x];
		if( ( hhh->mtl == head ) || ( hhh->mtl == beard ) || ( hhh->mtl == eyelash ) || ( hhh->mtl == eyebrow ) || ( hhh->mtl == splines ) )
			if( hhh->restlength > 0.0001 )
						{
							int             slg = 0;
							BASEHAIR        rr;
							BASEHAIR        tmp;
							float           tmpstiff;
							float           stiff3 = 1.0f;
							float           root_stiffness, dampening_force;
							float           idampening_force;

							root_stiffness = sliders[21][0].value * hhh->slider[21];
							dampening_force = sliders[40][0].value * hhh->slider[40];
							idampening_force = 1.0f - dampening_force;
							if( head >= 0 )
								if( hhh->mtl == head )
								{
									slg = 0;
									stiff3 = sliders[8][0].value * hhh->slider[8];
								}
							if( splines >= 0 )
								if( hhh->mtl == splines )
								{
									slg = 4;
									stiff3 = sliders[8][4].value * hhh->slider[8];
								}
							tmpstiff = stiff3;
							if( hhh->restlength > 0.0001 )
								if( sliders[22][slg].value > 0 )
									hair[x] = displace_wiggle( hair[x], 82.0 * sliders[22][slg].value, sliders[23][slg].value );

							if( maxrest[slg] > 0 )
							{
								float           stf;

//              if ((minrest[slg]/maxrest[slg])<.7)
//  stf=(1.0f-stiff3)*(((hhh->restlength-minrest[slg])/(maxrest[slg]-minrest[slg])));
//  stiff3+=(1.0f-stiff3)*((hhh->restlength-0)/(maxrest[slg]-0));
							}
							if( stiff3 > 1.0f )
								stiff3 = 1.0f;
							if( stiff3 < .999 )
							{
								int             yy;
								float           dxx, dyy, dzz;
								float           ys;	// this is y from 0 to 1;
								float           iys;	// this is y from 1 to 0;

								stiff = stiff3;

// randomize stiffness								stiff = stiff * ( drand98(  ) * .6 - .3 ) + stiff;
								if( stiff < 0 )
									stiff = 0;
								if( stiff > 1 )
									stiff = 1;
								memcpy( &rr, &hair[x], sizeof( BASEHAIR ) );
								hhh->velocity[0].x = 0;
								hhh->velocity[0].y = 0;
								hhh->velocity[0].z = 0;
								hhh->velocity[1].x = 0;
								hhh->velocity[1].y = 0;
								hhh->velocity[1].z = 0;
//   tmp.hv[0]=hhh->hv[0];
//								for( yy = 0; yy < 15; yy++ )
//								{
//									hhh->velocity[yy].x *= idampening_force;
//									hhh->velocity[yy].y *= idampening_force;
//									hhh->velocity[yy].z *= idampening_force;
//								}


//lh=hair[x];     

//   tmp.hv[0]=hhh->hv[1];
								dxx = hhh->hv[0].x - hhh->lasthv[0].x;
								dyy = hhh->hv[0].y - hhh->lasthv[0].y;
								dzz = hhh->hv[0].z - hhh->lasthv[0].z;



								tmp = hair[x];	// back up wannabe position
								{
									float           stf2;

									stf2 = stiff * .8;
									for( y = 1; y < 15; y++ )
									{
//          VERT qqq;
										float           pw;
										float           ipw;

										ys = y/14.0f;
										iys = 1.0f - ys;
										pw = ytable3[y];
										ipw = iytable4[y];


										// this is the normalizing strength
										{
											float           damp;
											float           damp2;
damp=ipw;

											damp = iys*stiff3;	// root
											damp+= stiff3*stiff3;
											if (damp>1.0f) damp=1.0f;

rr.hv[y]=hhh->lasthv[y];
											rr.hv[y].x += ( hhh->hv[y].x - hhh->lasthv[y].x ) * damp;
											rr.hv[y].y += ( hhh->hv[y].y - hhh->lasthv[y].y ) * damp;
											rr.hv[y].z += ( hhh->hv[y].z - hhh->lasthv[y].z ) * damp;
											
											damp = ipw*ipw*root_stiffness;	// root
											rr.hv[y].x += ( hhh->hv[y].x - rr.hv[y].x ) * damp;
											rr.hv[y].y += ( hhh->hv[y].y - rr.hv[y].y ) * damp;
											rr.hv[y].z += ( hhh->hv[y].z - rr.hv[y].z ) * damp;

										}


										{


										hhh->hv[y]=rr.hv[y];
										}

    hhh->hv[y].x+=hhh->velocity[y].x*(1.0-stiff3); // move with velocity from last frame
    hhh->hv[y].y+=hhh->velocity[y].y*(1.0-stiff3);
    hhh->hv[y].z+=hhh->velocity[y].z*(1.0-stiff3);
									hhh->velocity[y].x+=(rr.hv[y].x-hhh->lasthv[y].x)*(1.0-stiff3); // accumulate velocity
									hhh->velocity[y].y+=(rr.hv[y].y-hhh->lasthv[y].y)*(1.0-stiff3); // from this one
									hhh->velocity[y].z+=(rr.hv[y].z-hhh->lasthv[y].z)*(1.0-stiff3);
									}	//y
								}
//#ifdef RENDERLW
//for (y=0;y<15;y++) hhh->hv[y]=hhh->lasthv[y];

#ifndef NOLIB
//       apply_wind(x);

#endif
								memcpy( &tmp, &hair[x], sizeof( BASEHAIR ) );	// tmp is wannabe + velocity
//                              recalc_hair(x);
								dyn_recalc_hair( x );	// gottabe
								if( Gcollision_hit )
								{
									Gcollision_hit = 0;
									dyn_recalc_hair( x );	// gottabe
								}
								if( Gcollision_hit )
								{
									Gcollision_hit = 0;
									dyn_recalc_hair( x );	// gottabe
								}
								for( y = 1; y < 15; y++ )	// sum in secondary motion
								{
//         float stiff2;
//         float stiff5;
									float           damp;

//used?                                 float           pw, ipw, pww;

//         int ii;
									ys = ( float ) ytable[y];
									iys = 1.0f - ys;
//          slow                        pw = pow( ys, 2.0f );
//      slow                            ipw = pow( ys, 2.0f );
///         used?                       pw=ys*ys;
////        used?                           ipw=pw;
									////    used?                           ipw = 1.0f - ipw;
									//// used?                          pww = 1.0f - pow( 1.0f - ys, 2.5 );
#ifndef NOLIB
#ifdef SOFTIMAGE
									{
										VERT tmpvel;
										float bent;
										bent=pow(ytable[y],3.5);
										hhh->velocity[y].x *= restBOUNDLENGTH;
										hhh->velocity[y].y *= restBOUNDLENGTH;
										hhh->velocity[y].z *= -restBOUNDLENGTH;
										tmpvel.x=0;
										tmpvel.y=0;
										tmpvel.z=0;
										SHAVEcoord_convertFROMSHAVE( &hhh->velocity[y] );

										if( hhh->restlength > 0.0001 )
											MAYAapply_cached_forces( index, y, &tmpvel );

//										if (hhh->vid==34)
//										{
//											printf ("tmpvel = %f %f %f\n",tmpvel.x,tmpvel.y,tmpvel.z);
//											fflush(stdout);
//										}
										tmpvel.x*=bent;
										tmpvel.y*=bent;
										tmpvel.z*=bent;
										hhh->velocity[y].x+=tmpvel.x;
										hhh->velocity[y].y+=tmpvel.y;
										hhh->velocity[y].z+=tmpvel.z;

										SHAVEcoord_convertTOSHAVE( &hhh->velocity[y] );
										hhh->velocity[y].x /= restBOUNDLENGTH;
										hhh->velocity[y].y /= restBOUNDLENGTH;
										hhh->velocity[y].z /= restBOUNDLENGTH;
										hhh->velocity[y].z *= -1;
									}
#endif
#endif
#ifndef NOLIB
									hhh->velocity[y].y -= gravity * .00451872253 * ys* ( 1.0f - stiff3 );	//gravity
#endif //secondary motion



									hhh->velocity[y].x += (  hhh->hv[y].x-tmp.hv[y].x ) * ( 1.0f - stiff3 ) * .9*(idampening_force);	//.15
									hhh->velocity[y].y += (  hhh->hv[y].y-tmp.hv[y].y ) * ( 1.0f - stiff3 ) * .9*(idampening_force);
									hhh->velocity[y].z += (  hhh->hv[y].z-tmp.hv[y].z ) * ( 1.0f - stiff3 ) * .9*(idampening_force);
//     damp=.75+(1.0f-stiff)*.252;
////       damp=.7+(1.0f-stiff)*.252;
//									damp = ( float ) .992;
									//  damp-=(ipw);
//									if( damp < 0 )
//										damp = 0;
//     damp=pww*(1.0f-stiff);
									// damp=1.0f;
//									hhh->velocity[y].x *= damp;	//*(.2/BOUNDLENGTH);// dampening
//									hhh->velocity[y].y *= damp;	//*(.2/BOUNDLENGTH);
//									hhh->velocity[y].z *= damp;	//*(.2/BOUNDLENGTH);

									hhh->velocity[y].x *= idampening_force*.99;	//*(.2/BOUNDLENGTH);// dampening
									hhh->velocity[y].y *= idampening_force*.99;	//*(.2/BOUNDLENGTH);
									hhh->velocity[y].z *= idampening_force*.99;	//*(.2/BOUNDLENGTH);

									
//										hhh->hv[y].x+=hhh->velocity[y].x;
//										hhh->hv[y].y+=hhh->velocity[y].y;
//										hhh->hv[y].z+=hhh->velocity[y].z;


//      for (ii=0;ii<4;ii++)
//      {
//      pww=pow((1.0f-ys),6.0f);
//      if (pww>1.0f) pww=1.0f;
//      damp=pww;
									//         damp=ipw*ipw*ipw; // normalize
//         if (damp>1) damp=1;
//  
//     hhh->hv[y].x=rr.hv[y].x*damp+hhh->hv[y].x*(1.0f-damp);
//     hhh->hv[y].y=rr.hv[y].y*damp+hhh->hv[y].y*(1.0f-damp);
//     hhh->hv[y].z=rr.hv[y].z*damp+hhh->hv[y].z*(1.0f-damp);
//      }




								}	//y


							}	// stiff < .99
						}		//x
	}
//	Sdodyn3_single(index);



//	for( x = 0; x < totalverts; x++ )
	x=index;
	{
		BASEHAIR       *hhh;

		hhh = &hair[x];


					if( hhh->restlength > 0.0001 )


						for( qq = 0; qq < 15; qq++ )
						{
							hhh->hv[qq].x *= ( float ) restBOUNDLENGTH;
							hhh->hv[qq].y *= ( float ) restBOUNDLENGTH;
							hhh->hv[qq].z *= ( float ) restBOUNDLENGTH;
							hhh->resthv[qq].x *= ( float ) restBOUNDLENGTH;
							hhh->resthv[qq].y *= ( float ) restBOUNDLENGTH;
							hhh->resthv[qq].z *= ( float ) restBOUNDLENGTH;

							hhh->lasthv[qq].x *= ( float ) restBOUNDLENGTH;
							hhh->lasthv[qq].y *= ( float ) restBOUNDLENGTH;
							hhh->lasthv[qq].z *= ( float ) restBOUNDLENGTH;
//                          hhh->resthv[qq].x *= ( float )restBOUNDLENGTH;
//                          hhh->resthv[qq].y *= ( float )restBOUNDLENGTH;
//                          hhh->resthv[qq].z *= ( float )restBOUNDLENGTH;
//                          hhh->noisev[qq].x *= ( float )restBOUNDLENGTH;
//                          hhh->noisev[qq].y *= ( float )restBOUNDLENGTH;
//                          hhh->noisev[qq].z *= ( float )restBOUNDLENGTH;

//      thishhh->hv[qq].x*=(float)restBOUNDLENGTH;
//      thishhh->hv[qq].y*=(float)restBOUNDLENGTH;
//      thishhh->hv[qq].z*=(float)restBOUNDLENGTH;

//      lasthhh->hv[qq].x*=(float)restBOUNDLENGTH;
//      lasthhh->hv[qq].y*=(float)restBOUNDLENGTH;
//      lasthhh->hv[qq].z*=(float)restBOUNDLENGTH;

//      resthhh->hv[qq].x*=(float)restBOUNDLENGTH;
//      resthhh->hv[qq].y*=(float)restBOUNDLENGTH;
//      resthhh->hv[qq].z*=(float)restBOUNDLENGTH;
						}
	}
//	for( x = 0; x < totalverts; x++ )
		if( hair[x].restlength > 0.0001 )
			hair[x].restlength *= restBOUNDLENGTH;
//#ifdef NOLIB


}								// end







void
Sdodyn3_single( int index )
{

	int             x, y, qq;
	BASEHAIR        lh;
	float           stiff = ( float ) .1;
	float           maxrest[5];
	float           minrest[5];
if (index<total_splines)
{
//	for( x = 0; x < total_splines; x++ )
x=index;
//fprintf (stdout,"Sdodyn %d\n",x);fflush(stdout);
		if( Shair[x].restlength > 0.0001 )
		{
			Shair[x].restlength = Distance( Shair[x].hv[1], Shair[x].hv[2] );
			Sresthair[x].restlength = Shair[x].restlength;
			Sthishair[x].restlength = Shair[x].restlength;
			Slasthair[x].restlength = Shair[x].restlength;

			Shair[x].restlength /= ( float ) restBOUNDLENGTH;
		}
//	for( x = 0; x < total_splines; x++ )
x=index;
		if( Shair[x].restlength > 0.0001 )
				{

					for( qq = 0; qq < 15; qq++ )
					{
						Shair[x].hv[qq].x /= ( float ) restBOUNDLENGTH;
						Shair[x].hv[qq].y /= ( float ) restBOUNDLENGTH;
						Shair[x].hv[qq].z /= ( float ) restBOUNDLENGTH;
						Sthishair[x].hv[qq].x /= ( float ) restBOUNDLENGTH;
						Sthishair[x].hv[qq].y /= ( float ) restBOUNDLENGTH;
						Sthishair[x].hv[qq].z /= ( float ) restBOUNDLENGTH;

						Slasthair[x].hv[qq].x /= ( float ) restBOUNDLENGTH;
						Slasthair[x].hv[qq].y /= ( float ) restBOUNDLENGTH;
						Slasthair[x].hv[qq].z /= ( float ) restBOUNDLENGTH;
						Sresthair[x].hv[qq].x /= ( float ) restBOUNDLENGTH;
						Sresthair[x].hv[qq].y /= ( float ) restBOUNDLENGTH;
						Sresthair[x].hv[qq].z /= ( float ) restBOUNDLENGTH;
//  hair[x].restlength/(float)BOUNDLENGTH;  
					}
				}

	maxrest[0] = 0;
	maxrest[1] = 0;
	maxrest[2] = 0;
	maxrest[3] = 0;
	maxrest[4] = 0;
	minrest[0] = 1110;
	minrest[1] = 1110;
	minrest[2] = 1110;
	minrest[3] = 1110;
	minrest[4] = 1110;

//	for( x = 0; x < total_splines; x++ )
	{
		if( Shair[x].restlength > maxrest[4] )
			maxrest[4] = Shair[x].restlength;
		if( Shair[x].restlength < minrest[4] )
			minrest[4] = Shair[x].restlength;
	}
//	for( x = 0; x < total_splines; x++ )

x=index;

		if( Shair[x].restlength > 0.0001 )
					{
						int             slg;
						BASEHAIR        rr;
						BASEHAIR        tmp;
						float           tmpstiff;
						float           stiff3 = 1.0f;
						float           root_stiffness, dampening_force;

						if( splines >= 0 )
							if( Shair[x].mtl == splines )
							{
								slg = 4;
								stiff3 = sliders[8][4].value * Shair[x].slider[8];
							}
						tmpstiff = stiff3;
						if( Shair[x].restlength > 0.0001 )
							Shair[x] = displace_wiggle( Shair[x], 82.0 * sliders[22][slg].value, sliders[23][slg].value );

						root_stiffness = sliders[21][4].value * Shair[x].slider[21];
						dampening_force = sliders[40][4].value * Shair[x].slider[40];


						if( maxrest[slg] > 0 )
						{
							float           stf;

//              if ((minrest[slg]/maxrest[slg])<.7)
//  stf=(1.0f-stiff3)*(((hair[x].restlength-minrest[slg])/(maxrest[slg]-minrest[slg])));
//                          stf =
//                              ( 1.0f -
//                                stiff3 );/* ( ( ( Shair[x].restlength - 0 ) /
							//       ( maxrest[slg] - 0 ) ) ); */
//                          stf = 1.0f - stf;
//                          stf *= stf;
//                          stf *= stf;
//                          stiff3 = stf;
//  stiff3+=(1.0f-stiff3)*((hair[x].restlength-0)/(maxrest[slg]-0));
						}
						if( stiff3 > 1.0f )
							stiff3 = 1.0f;
						if( tmpstiff < .99 )
						{
							float           dxx, dyy, dzz;
							float           ys;	// this is y from 0 to 1;
							float           iys;	// this is y from 1 to 0;

							stiff = stiff3-.05;
// randomise stiffness							stiff = stiff * ( drand98(  ) * .6 - .3 ) + stiff;
							if( stiff < 0 )
								stiff = 0;
							if( stiff > 1 )
								stiff = 1;
							rr = Shair[x];
							Shairvelocity[x].hv[0].x = 0;
							Shairvelocity[x].hv[0].y = 0;
							Shairvelocity[x].hv[0].z = 0;
							Shairvelocity[x].hv[1].x = 0;
							Shairvelocity[x].hv[1].y = 0;
							Shairvelocity[x].hv[1].z = 0;
//   tmp.hv[0]=hair[x].hv[0];
							{

								int             yy;

								for( yy = 0; yy < 15; yy++ )
								{
									Shairvelocity[x].hv[yy].x *= ( 1.0f - dampening_force );
									Shairvelocity[x].hv[yy].y *= ( 1.0f - dampening_force );
									Shairvelocity[x].hv[yy].z *= ( 1.0f - dampening_force );
								}
							}

							lh = Slasthair[x];

//   tmp.hv[0]=hair[x].hv[1];
							dxx = Shair[x].hv[0].x - Slasthair[x].hv[0].x;
							dyy = Shair[x].hv[0].y - Slasthair[x].hv[0].y;
							dzz = Shair[x].hv[0].z - Slasthair[x].hv[0].z;
							tmp = Shair[x];	// back up wannabe position

							for( y = 1; y < 15; y++ )
							{
//          VERT qqq;
								float           pw;
								float           ipw;
								float           stf2;

//      stf2=stiff*.4+.25; root
								stf2 = stiff;	// root
								if( stf2 > 1 )
									stf2 = 1;
//stf2=1;
								ys = ( float ) ( y ) / 15.0f;
								iys = 1.0f - ys;
								pw = pow( ys, 2.7 );
								ipw = pow( iys, 4.0f );
								lh.hv[y] = Slasthair[x].hv[y];
								lh.hv[y].x += dxx * stf2 * .8;
								lh.hv[y].y += dyy * stf2 * .8;
								lh.hv[y].z += dzz * stf2 * .8;
								// this is the normalizing strength
								{
									float           damp;

//         float stiff2;
									ipw *= ( 1.0f - stiff );
									ipw += stiff;
									if( ipw > 1 )
										ipw = 1;
									damp = ipw * 0.937;	// primary motion
////           damp=(1.0f-stiff)*0.237; // primary motion
									damp += ( float ) .02;
									damp *= damp;
//         damp=ipw*0.137; // primary motion
									if( damp > 1 )
										damp = 1;
									damp *= root_stiffness;	// root


									Shairvelocity[x].hv[y].x += ( rr.hv[y].x - lh.hv[y].x ) * damp;
									Shairvelocity[x].hv[y].y += ( rr.hv[y].y - lh.hv[y].y ) * damp;
									Shairvelocity[x].hv[y].z += ( rr.hv[y].z - lh.hv[y].z ) * damp;

								}
								Shair[x].hv[y].x = lh.hv[y].x + Shairvelocity[x].hv[y].x;
								Shair[x].hv[y].y = lh.hv[y].y + Shairvelocity[x].hv[y].y;
								Shair[x].hv[y].z = lh.hv[y].z + Shairvelocity[x].hv[y].z;

////       Shair[x].hv[y].x=lh.hv[y].x+Shairvelocity[x].hv[y].x;
////       Shair[x].hv[y].y=lh.hv[y].y+Shairvelocity[x].hv[y].y;
////       Shair[x].hv[y].z=lh.hv[y].z+Shairvelocity[x].hv[y].z;
///    hair[x].hv[y].x=hair[x].hv[y].x+hairvelocity[x].hv[y].x;
///    hair[x].hv[y].y=hair[x].hv[y].y+hairvelocity[x].hv[y].y;
///    hair[x].hv[y].z=hair[x].hv[y].z+hairvelocity[x].hv[y].z;
//         hair[x].hv[y].x-=hairvelocity[x].hv[y].x;
//     hair[x].hv[y].y-=hairvelocity[x].hv[y].y;
//     hair[x].hv[y].z-=hairvelocity[x].hv[y].z;
								{
									int             ii;
									float           damp, pww;
									float           stf;

									stf = stiff + .15;
									if( stf > 1.0f )
										stf = 1.0f;
									for( ii = 0; ii < 1; ii++ )
									{
//      pww=pow((1.0f-(y/15.0f)),4.5);
										pww = ( 1.0f - ( y / 15.0f ) );
										pww += stiff;
										if( pww > 1.0f )
											pww = 1.0f;
										damp = pww * .2;	//+(stiff/18.0f)*pww;
										//         damp=ipw*ipw*ipw; // normalize
										if( damp > 1 )
											damp = 1;
										damp *= root_stiffness;	// root

										Shair[x].hv[y].x = rr.hv[y].x * damp + Shair[x].hv[y].x * ( 1.0f - damp );
										Shair[x].hv[y].y = rr.hv[y].y * damp + Shair[x].hv[y].y * ( 1.0f - damp );
										Shair[x].hv[y].z = rr.hv[y].z * damp + Shair[x].hv[y].z * ( 1.0f - damp );


									}


									damp = stiff;	//*stiff;  // normalize  again
									damp = damp + .05;	// normalize  again
									if( damp > 1 )
										damp = 1;
									Shair[x].hv[y].x = rr.hv[y].x * damp + Shair[x].hv[y].x * ( 1.0f - damp );
									Shair[x].hv[y].y = rr.hv[y].y * damp + Shair[x].hv[y].y * ( 1.0f - damp );
									Shair[x].hv[y].z = rr.hv[y].z * damp + Shair[x].hv[y].z * ( 1.0f - damp );
//     hair[x].hv[y].x= rr.hv[y].x*damp + hair[x].hv[y].x*(1.0f-damp);
//     hair[x].hv[y].y= rr.hv[y].y*damp + hair[x].hv[y].y*(1.0f-damp);
//     hair[x].hv[y].z= rr.hv[y].z*damp + hair[x].hv[y].z*(1.0f-damp);


								}


							}	//y

#ifndef NOLIB
//       Sapply_wind(x);

#endif
							memcpy( &tmp, &Shair[x], sizeof( BASEHAIR ) );	// tmp is wannabe + velocity
							Sdyn_recalc_hair( x );	// gottabe
							Sdyn_recalc_hair( x );	// gottabe
							Sdyn_recalc_hair( x );	// gottabe
							Sdyn_recalc_hair( x );	// gottabe

							for( y = 1; y < 15; y++ )	// sum in secondary motion
							{
//         float stiff2;
								float           stiff5;
								float           damp;
								float           pw, ipw, pww;

//         int ii;
								ys = ( float ) ( y ) / 15.0f;
								iys = 1.0f - ys;
								pw = pow( ys, 2.0f );
								ipw = pow( ys, 2.0f );
								ipw = 1.0f - ipw;
								pww = 1.0f - pow( 1.0f - ys, 2.5 );

#ifdef SOFTIMAGE
								{
									Shairvelocity[x].hv[y].x *= restBOUNDLENGTH;
									Shairvelocity[x].hv[y].y *= restBOUNDLENGTH;
									Shairvelocity[x].hv[y].z *= -restBOUNDLENGTH;
									SHAVEcoord_convertFROMSHAVE( &Shairvelocity[x].hv[y] );
									if( Shair[x].restlength > 0.0001 )
										MAYAapply_cached_forces( x, y, &Shairvelocity[x].hv[y] );

									SHAVEcoord_convertTOSHAVE( &Shairvelocity[x].hv[y] );
									Shairvelocity[x].hv[y].x /= restBOUNDLENGTH;
									Shairvelocity[x].hv[y].y /= restBOUNDLENGTH;
									Shairvelocity[x].hv[y].z /= restBOUNDLENGTH;
									Shairvelocity[x].hv[y].z *= -1.0f;
								}
#endif

#ifndef NOLIB
								Shairvelocity[x].hv[y].y -= gravity * .00251872253;	//gravity
#endif
								//secondary motion   
								Shairvelocity[x].hv[y].x -= ( tmp.hv[y].x - Shair[x].hv[y].x ) * ( 1.0f - stiff ) * .7;	//.15
								Shairvelocity[x].hv[y].y -= ( tmp.hv[y].y - Shair[x].hv[y].y ) * ( 1.0f - stiff ) * .7;
								Shairvelocity[x].hv[y].z -= ( tmp.hv[y].z - Shair[x].hv[y].z ) * ( 1.0f - stiff ) * .7;
//     damp=.75+(1.0f-stiff)*.252;
////       damp=.7+(1.0f-stiff)*.252;
								damp = ( float ) .95;
								//  damp-=(ipw);
								if( damp < 0 )
									damp = 0;
//     damp=pww*(1.0f-stiff);
								// damp=1.0f;

								Shairvelocity[x].hv[y].x *= damp;	//*(.2/BOUNDLENGTH);// dampening
								Shairvelocity[x].hv[y].y *= damp;	//*(.2/BOUNDLENGTH);
								Shairvelocity[x].hv[y].z *= damp;	//*(.2/BOUNDLENGTH);
//      for (ii=0;ii<4;ii++)
//      {
//      pww=pow((1.0f-ys),6.0f);
//      if (pww>1.0f) pww=1.0f;
//      damp=pww;
								//         damp=ipw*ipw*ipw; // normalize
//         if (damp>1) damp=1;
//  
//     hair[x].hv[y].x=rr.hv[y].x*damp+hair[x].hv[y].x*(1.0f-damp);
//     hair[x].hv[y].y=rr.hv[y].y*damp+hair[x].hv[y].y*(1.0f-damp);
//     hair[x].hv[y].z=rr.hv[y].z*damp+hair[x].hv[y].z*(1.0f-damp);
//      }




							}	//y

////            dyn_recalc_hair(x); // gottabe
							for( y = 1; y < 15; y++ )
							{
								int             ii;
								float           pww;
								float           damp;

								if( 0 == 1 )
									for( ii = 0; ii < 1; ii++ )
									{
										ys = ( float ) ( y ) / 15.0f;
										iys = 1.0f - ys;
										pww = pow( ( 1.0f - ys ), 6.0f );

										if( pww > 1.0f )
											pww = 1.0f;
										damp = pww * .008;
										//         damp=ipw*ipw*ipw; // normalize again
										if( damp > 1 )
											damp = 1;

										Shair[x].hv[y].x = rr.hv[y].x * damp + Shair[x].hv[y].x * ( 1.0f - damp );
										Shair[x].hv[y].y = rr.hv[y].y * damp + Shair[x].hv[y].y * ( 1.0f - damp );
										Shair[x].hv[y].z = rr.hv[y].z * damp + Shair[x].hv[y].z * ( 1.0f - damp );
									}
							}
							for( y = 1; y < 15; y++ )
							{
								int             ii;

								for( ii = 0; ii < 4; ii++ )
								{
									float           pww;
									float           damp;

//      pww=pow((1.0f-(y/15.0f)),(1.0f-stiff)*4.0+1);
									pww = ( float ) y / 15.0f;
//      pww*=1.8;
									pww *= ( float ) 3.14;
									if( pww > 3.14 )
										pww = ( float ) 3.14;
									pww = cos( pww ) * ( float ) .5 + ( float ) .5;
//      pww=1.0f-pww;
									if( pww > 1.0f )
										pww = 1.0f;
//      damp=pww*.98;//+(stiff/18.0f)*pww;
									damp = pww;
									{
										float           stf;

										stf = stiff * 1.5;
										if( stf > 1 )
											stf = 1;
										damp *= stf;
									}
									//         damp=ipw*ipw*ipw; // normalize spline
									if( damp > 1 )
										damp = 1;
									{
										float           stf;

										stf = stiff + .5;
										if( stf > 1 )
											stf = 1;
//     hair[x].hv[y].x=rr.hv[y].x*damp+hair[x].hv[y].x*(1.0f-damp);
//     hair[x].hv[y].y=rr.hv[y].y*damp+hair[x].hv[y].y*(1.0f-damp);
//     hair[x].hv[y].z=rr.hv[y].z*damp+hair[x].hv[y].z*(1.0f-damp);
										Shairvelocity[x].hv[y].x += ( rr.hv[y].x - Shair[x].hv[y].x ) * damp * stf;
										Shairvelocity[x].hv[y].y += ( rr.hv[y].y - Shair[x].hv[y].y ) * damp * stf;
										Shairvelocity[x].hv[y].z += ( rr.hv[y].z - Shair[x].hv[y].z ) * damp * stf;
									}
								}


								for( ii = 0; ii < 4; ii++ )
								{
									float           pww;
									float           damp;

//      pww=pow((1.0f-(y/15.0f)),(1.0f-stiff)*4.0+1);
									pww = ( float ) y / 15.0f;
									pww *= ( float ) 1.8;
									pww *= ( float ) 3.14;
									if( pww > ( float ) 3.14 )
										pww = ( float ) 3.14;
									pww = cos( pww ) * ( float ) .5 + ( float ) .5;
//      pww=1.0f-pww;
									if( pww > ( float ) 1.0f )
										pww = ( float ) 1.0f;
//      damp=pww*.98;//+(stiff/18.0f)*pww;
									damp = pww;
									damp *= stiff;
									//         damp=ipw*ipw*ipw; // normalize
									if( damp > 1 )
										damp = 1;
									{







										Shair[x].hv[y].x = rr.hv[y].x * damp + Shair[x].hv[y].x * ( 1.0f - damp );
										Shair[x].hv[y].y = rr.hv[y].y * damp + Shair[x].hv[y].y * ( 1.0f - damp );
										Shair[x].hv[y].z = rr.hv[y].z * damp + Shair[x].hv[y].z * ( 1.0f - damp );
//     hairvelocity[x].hv[y].x+=(rr.hv[y].x-hair[x].hv[y].x)*damp*stf;
//     hairvelocity[x].hv[y].y+=(rr.hv[y].y-hair[x].hv[y].y)*damp*stf;
//     hairvelocity[x].hv[y].z+=(rr.hv[y].z-hair[x].hv[y].z)*damp*stf;
									}
								}
							}

//#ifdef LIB
							//  dyn_recalc_hair(x); // gottabe
//#endif
//      blur_velocity(x);
							//  blur_postion(x);
						}		// stiff < .99

					}			//x

//	for( x = 0; x < total_splines; x++ )

x=index;
					if( Shair[x].restlength > 0.0001 )

						for( qq = 0; qq < 15; qq++ )
						{
							Shair[x].hv[qq].x *= ( float ) restBOUNDLENGTH;
							Shair[x].hv[qq].y *= ( float ) restBOUNDLENGTH;
							Shair[x].hv[qq].z *= ( float ) restBOUNDLENGTH;
							Sthishair[x].hv[qq].x *= ( float ) restBOUNDLENGTH;
							Sthishair[x].hv[qq].y *= ( float ) restBOUNDLENGTH;
							Sthishair[x].hv[qq].z *= ( float ) restBOUNDLENGTH;

							Slasthair[x].hv[qq].x *= ( float ) restBOUNDLENGTH;
							Slasthair[x].hv[qq].y *= ( float ) restBOUNDLENGTH;
							Slasthair[x].hv[qq].z *= ( float ) restBOUNDLENGTH;

							Sresthair[x].hv[qq].x *= ( float ) restBOUNDLENGTH;
							Sresthair[x].hv[qq].y *= ( float ) restBOUNDLENGTH;
							Sresthair[x].hv[qq].z *= ( float ) restBOUNDLENGTH;
						}
//	for( x = 0; x < total_splines; x++ )
						x=index;
		if( Shair[x].restlength > 0.0001 )
			Shair[x].restlength *= restBOUNDLENGTH;
//#ifdef NOLIB
//	for( x = 0; x < total_splines; x++ )
					if( Shair[x].restlength > 0.0001 )
						Slasthair[x] = Shair[x];
//#endif
}
}								// end