/****************************************************************************
 *                  clothray.cpp
 *
 *  Author: Christophe Bouffartigue (tofbouf@oreka.com)
 *
 * from Persistence of Vision(tm) Ray Tracer version 3.6.
 * Copyright 1991-2003 Persistence of Vision Team
 * Copyright 2003-2004 Persistence of Vision Raytracer Pty. Ltd.
 *---------------------------------------------------------------------------
 * NOTICE: This source code file is provided so that users may experiment
 * with enhancements to POV-Ray and to port the software to platforms other
 * than those supported by the POV-Ray developers. There are strict rules
 * regarding how you are permitted to use this file. These rules are contained
 * in the distribution and derivative versions licenses which should have been
 * provided with this file.
 *
 * These licences may be found online, linked from the end-user license
 * agreement that is located at http://www.povray.org/povlegal.html
 *---------------------------------------------------------------------------
 * This program is based on the popular DKB raytracer version 2.12.
 * DKBTrace was originally written by David K. Buck.
 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
 *---------------------------------------------------------------------------
 *
 *===========================================================================
 * This file is part of MegaPOV, a modified and unofficial version of POV-Ray
 * For more information on MegaPOV visit our website:
 * http://megapov.inetart.net/
 *===========================================================================
 *
 * $RCSfile: clothray.cpp,v $
 * $Revision: 1.19 $
 * $Author: chris $
 *
 *****************************************************************************/

#include "frame.h"

#ifdef CLOTHRAY_PATCH

#define MAX_STR 1500



#include <algorithm>
#include "parse.h"
#include "colour.h"
#include "express.h"
#include "parsestr.h"
#include "vector.h"
#include "objects.h"
#include "pigment.h"
#include "tokenize.h"
#include "parstxtr.h"
#include "ray.h"
#include "pov_util.h"
#include "povray.h"
#include "povmsend.h"
#include "clothray.h"

BEGIN_POV_NAMESPACE



static int iCounter, n1Counter;
static CLOTH_SYS *current_sys;


/*****************************************************************************

* Local preprocessor defines

******************************************************************************/

#define AccessPoint(i,j) (cloth_sys->cloth->points[((i)*cloth_sys->cloth->nb2)+(j)])
#define AccessVel(i,j)   (cloth_sys->cloth->vel[((i)*cloth_sys->cloth->nb2)+(j)])
#define AccessConstraint(i,j)   (cloth_sys->cloth->constraint[((i)*cloth_sys->cloth->nb2)+(j)])
#define AccessResul(i,j) (resul[((i)*cloth_sys->cloth->nb2)+(j)])


/*****************************************************************************
* Local typedefs
******************************************************************************/


/*****************************************************************************
* Local typedefs
******************************************************************************/

static const DBL length_mul[3][3] = {
	{0.0, 1.0        , 2.0        },  /* 0,       1,       2 */
	{1.0, 1.414213562, 2.236067977},  /* 1, sqrt(2), sqrt(5) */
	{2.0, 2.236067977, 2.828427125},  /* 2, sqrt(5), sqrt(8) */
};


static const char format_triangle[] = "triangle { <%.9f, %.9f, %.9f>, <%.9f, %.9f, %.9f>, <%.9f, %.9f, %.9f>";
static const char format_smooth_triangle[] = "smooth_triangle { <%.9f, %.9f, %.9f>, <%.9f, %.9f, %.9f>,\n<%.9f, %.9f, %.9f>, <%.9f, %.9f, %.9f>,\n<%.9f, %.9f, %.9f>, <%.9f, %.9f, %.9f>";
static const char format_uv_vectors[] = "\nuv_vectors <%.9f, %.9f>, <%.9f, %.9f>, <%.9f, %.9f>";
static const int  La1[8][2] = {{-1,-1}, {0,-1}, {1,-1}, {1,0}, {1,1}, {0,1}, {-1,1}, {-1,0}};


/*****************************************************************************
* Static functions
******************************************************************************/

static void Iterate_Cloth_Simulation(CLOTH_SYS *cloth_sys, VECTOR *resul);
static void Detect_Collision(OBJECT * Object, VECTOR Vel, DBL Disp, VECTOR Start, VECTOR Result, DBL fric_obj);
static void GetNormal(CLOTH_SYS *cloth_sys, int i, int j, VECTOR Result);
static bool ReadClothFile(CLOTH_SYS *cloth_sys);
static bool WriteClothFile(CLOTH_SYS *cloth_sys);
static bool WriteClothMesh(CLOTH_SYS *cloth_sys);

static int Make_Cloth_Simulation(CLOTH_SYS *cloth_sys);

/*****************************************************************************
* Local variables
******************************************************************************/

static OStream *ClothMeshFile=NULL;
static  OStream *ClothFile=NULL;

/*****************************************************************************
* FUNCTION   void Terminate_Cloth_Sys();
* 
* If we do a user abort, the file is left open and a second render will not work
******************************************************************************/

void Terminate_Cloth_Sys()
{
	if ( ClothMeshFile)
	{
		delete ClothMeshFile;
		ClothMeshFile=NULL;
	}
	if ( ClothFile )
	{
		delete ClothFile;
		ClothFile=NULL;
	}
}

/*****************************************************************************
* FUNCTION   Create_Cloth_Sys()
******************************************************************************/

static CLOTH_SYS * Create_Cloth_Sys(void)

{
	CLOTH_SYS * New = (CLOTH_SYS *)POV_MALLOC(sizeof(CLOTH_SYS), "cloth system");
	New->cloth = NULL;
	Make_Vector(New->gravity, 0, 0, 0);
	New->wind = NULL;
	New->object = NULL;
	New->neighbors = 1;
	New->smooth_mesh = 0;
	New->internal_collision = 0;
	New->friction = 1.0;
	New->viscosity = 0.0;
	New->damping = 0.90;
	New->dt = 0.05;
	New->nbit = 1;
	New->input = NULL;
	New->output = NULL;
	New->mesh_output = NULL;
    return New;
}

static void Destroy_Cloth_Sys(CLOTH_SYS *cloth_sys)
{
	if (cloth_sys)
	{
		if (cloth_sys->cloth)
		{
			if (cloth_sys->cloth->points) 
				POV_FREE(cloth_sys->cloth->points);
			if (cloth_sys->cloth->vel) 
				POV_FREE(cloth_sys->cloth->vel);
			if (cloth_sys->cloth->constraint) 
				POV_FREE(cloth_sys->cloth->constraint);
			POV_FREE(cloth_sys->cloth);
		}

		if (cloth_sys->wind) 
			Destroy_Pigment(cloth_sys->wind);

		if (cloth_sys->object) 
			POV_FREE(cloth_sys->object);

		if (cloth_sys->input) 
			POV_FREE(cloth_sys->input);
		if (cloth_sys->output)
			POV_FREE(cloth_sys->output);
		if (cloth_sys->mesh_output) 
			POV_FREE(cloth_sys->mesh_output);
		
		POV_FREE(cloth_sys);
	}
}


/*****************************************************************************
* FUNCTION   Parse_Cloth_Sys()
******************************************************************************/
void Parse_Cloth_Sys(void)
{
	CLOTH_SYS * cloth_sys = Create_Cloth_Sys();

	Parse_Begin();

	EXPECT
		/* obstacle */
		CASE (ENVIRONMENT_TOKEN)
			cloth_sys->object = Parse_Object();
			if(cloth_sys->object == NULL)
				Expectation_Error ("object or directive");
		END_CASE

		/* gravity */
		CASE (GRAVITY_TOKEN)
			Parse_Vector(cloth_sys->gravity);
		END_CASE

		/* wind */
		CASE (WIND_TOKEN)
			cloth_sys->wind = Create_Pigment();
			Parse_Begin();
			Parse_Pigment(&(cloth_sys->wind));
			Parse_End();
			Post_Pigment(cloth_sys->wind);
		END_CASE

		/* type */
		CASE2 (NEIGHBORS_TOKEN,NEIGHBOURS_TOKEN)
			cloth_sys->neighbors = (int)Allow_Float(1.0);
		END_CASE

		/* internal_collision */
		CASE (INTERNAL_COLLISION_TOKEN)
			cloth_sys->internal_collision = Allow_Float(1.0) > 0.0;
		END_CASE

		/* friction */
		CASE (FRICTION_TOKEN)
			cloth_sys->friction = Allow_Float(1.0);
		END_CASE

		/* viscosity */
		CASE (VISCOSITY_TOKEN)
			cloth_sys->viscosity = Allow_Float(1.0);
		END_CASE

		/* damping */
		CASE (DAMPING_TOKEN)
			cloth_sys->damping = Allow_Float(1.0);
		END_CASE

		/* intervals */
		CASE (INTERVALS_TOKEN)
			cloth_sys->dt = Parse_Float();
			if (cloth_sys->dt < 0)
			{
				cloth_sys->dt = -cloth_sys->dt;
				Warning(0,"negative value for intervals; intervals negated\n");
			}
			else if (cloth_sys->dt == 0)
			{
				cloth_sys->dt = 0.05;
				Warning(0,"0.0 value for intervals; using default value (0.05)\n");
			}
		END_CASE

		/* iterations */
        CASE (ITERATIONS_TOKEN)
            cloth_sys->nbit = (long )Parse_Float();
			if (cloth_sys->nbit <= 0)
			{
				cloth_sys->dt = 1;
				Warning(0,"negative or null value for iterations; using default value (1)\n");
			}
        END_CASE

		/* input */
        CASE (INPUT_TOKEN)
			cloth_sys->input =Parse_C_String(true);
		END_CASE

		/* output */
        CASE (OUTPUT_TOKEN)
			cloth_sys->output = Parse_C_String(true);
		END_CASE

		/* mesh_output */
        CASE (MESH_OUTPUT_TOKEN)
			cloth_sys->mesh_output = Parse_C_String(true);
		END_CASE

		/* smooth_mesh */
        CASE (SMOOTH_MESH_TOKEN)
			cloth_sys->smooth_mesh = Allow_Float(1.0) > 0.0;
		END_CASE

		/* uv_mesh */
        CASE (UV_MESH_TOKEN)
			cloth_sys->uv_mesh = Allow_Float(1.0) > 0.0;
		END_CASE

 
		OTHERWISE
			UNGET
			EXIT
		END_CASE
	END_EXPECT

	Parse_End();

	/* doing simulation */
	if (cloth_sys->input != NULL)
	{
		if (Make_Cloth_Simulation(cloth_sys))
		{
			Debug_Info("Cloth simulation OK.\n");
		}
		else
		{
			Debug_Info("Cloth simulation failed !!!!\n");
		}
	}
	else
	{
		Error("Cloth simulation: no input file !!!!");
	}
    //back to parsing
  	 Send_Progress("Parsing", PROGRESS_PARSING);

	Destroy_Cloth_Sys(cloth_sys);
}




/*****************************************************************************
* FUNCTION   Detect_Collision()
******************************************************************************/
static void Detect_Collision(OBJECT * Object, VECTOR Vel, DBL dt, VECTOR Start, VECTOR Result, DBL fric_obj)
{ 
 INTERSECTION Intersect;
 RAY Ray;
 VECTOR Norm;
 VECTOR Dir;
 VECTOR IPt;
 VECTOR posDisp;
 DBL dist;
 DBL speed;
 DBL n;

	Initialize_Ray_Containers(&Ray);

	VNormalize(Dir, Vel);
	Assign_Vector(Ray.Direction, Dir);
	Assign_Vector(Ray.Initial, Start);

	if(Intersection(&Intersect, Object, &Ray))
	{
		Assign_Vector(IPt, Intersect.IPoint);
		VLength(speed, Vel);
		VDist(dist, IPt, Start);
		if (dist < (dt*speed))
		{
			Normal(Norm, Intersect.Object, &Intersect);
			VNormalizeEq(Norm);

			/*Set new velocity*/
			VDot(n, Vel, Norm);
			n = -n;
			VAddScaledEq(Vel, n, Norm);
			VScaleEq(Vel, fric_obj);
			VAddScaledEq(Vel, n*0.1, Norm);

			/*Calculate new position*/
			VAddScaled(Result, IPt, (dist*.05), Norm);
			VInverseScale(posDisp, Vel, speed); VScaleEq(posDisp, ((dt*speed)-dist));
			VAddEq(Result, posDisp);

			if (Inside_Object(Result, Object))
			{
				VAddScaled(Result, IPt, (dist*.05), Norm);
			}
		}
		else
		{
			VAddScaled(Result, Start, dt, Vel);
			if (Inside_Object(Result, Object))
			{
				Assign_Vector(Result, Start);
			}
		}
	}
	else
	{
		VAddScaled(Result, Start, dt, Vel);
		if (Inside_Object(Result, Object))
		{
			Assign_Vector(Result, Start);
		}
	}
}

/*****************************************************************************
* FUNCTION   Iterate_Cloth_Simulation()
******************************************************************************/
static void Iterate_Cloth_Simulation(CLOTH_SYS *cloth_sys, VECTOR *resul)
{
 int i, j;
 int k, l;
 int ecart, ecart_col, factor;
 VECTOR String;
 VECTOR Force;
 VECTOR Move;
 VECTOR ApWind;
 COLOUR Wind;
 VECTOR resultmp;
 DBL windforce;
 DBL length,  nlength;
 DBL test_length, test_length2;
 DBL ForceScaler;


	ecart = 1 + cloth_sys->neighbors;
	/* don't calculate internal collision for points joined by springs and their neighbors */
	ecart_col = ecart + 1;
	
	for (i=0; i<cloth_sys->cloth->nb1; i++)
	{
		/* Calculation of the force exerced by the 24 nearest neighbours */
		for (j=0; j<cloth_sys->cloth->nb2; j++)
		{
			/* Gravity */
			Assign_Vector(Force, cloth_sys->gravity);


			for(k=-ecart; k<=ecart; k++)
			{
				for(l=-ecart; l<=ecart; l++)
				{
					if ( ((k != 0) || (l != 0)) &&
					     ((i+k) < cloth_sys->cloth->nb1) &&
						 ((i+k) >= 0) &&
						 ((j+l) < cloth_sys->cloth->nb2) &&
						 ((j+l) >= 0) )
					{
						VSub(String, AccessPoint(i+k, j+l), AccessPoint(i, j));
						VLength(length, String);
						if (length == 0)
						{
							Make_Vector(String,0,0,0);
						}
						else
						{
							nlength = length_mul[abs(k)][abs(l)]*(cloth_sys->cloth->nlng);
							ForceScaler = (cloth_sys->cloth->ks)*(length - nlength)/length;
							VScaleEq(String, ForceScaler);
						}
						VAddEq(Force, String);
					}
				}
			}

			if (cloth_sys->viscosity)
			{
				/* Calculate Force exerced by the wind and/or air friction */
				GetNormal(cloth_sys, i, j, String);
				VNormalizeEq(String);
				factor = 0;
				if (i>0) factor++;
				if (j>0) factor++;
				if ((i+1)<cloth_sys->cloth->nb1) factor ++;
				if ((j+1)<cloth_sys->cloth->nb2) factor ++;

				if (cloth_sys->wind != NULL)
				{
					Compute_Pigment(Wind, cloth_sys->wind, AccessPoint(i,j), NULL);
					Make_Vector(ApWind, Wind[0], Wind[1], Wind[2]);
					VSubEq(ApWind, AccessVel(i,j));
				}
				else
				{
					VScale(ApWind, AccessVel(i,j), -1);
				}
				VDot(windforce, ApWind, String);
				windforce = cloth_sys->viscosity*factor*windforce/4.0;
				VAddScaledEq(Force, windforce, String);
			}

			if (cloth_sys->internal_collision)
			{
				test_length = .75*.75*cloth_sys->cloth->nlng*cloth_sys->cloth->nlng;
				test_length2 = .75*cloth_sys->cloth->nlng;

				/* Calculate internal collision repulsive force */
				for (k=0; k<cloth_sys->cloth->nb1; k++)
				{
					for (l=0; l<cloth_sys->cloth->nb2; l++)
					{
						if ((abs(i-k)>ecart_col) || (abs(j-l)>ecart_col))
						{
							/* kind of bounding boxes test */
							if ((fabs(AccessPoint(i, j)[X] - AccessPoint(k, l)[X]) < test_length2) &&
							    (fabs(AccessPoint(i, j)[Y] - AccessPoint(k, l)[Y]) < test_length2) &&
							    (fabs(AccessPoint(i, j)[Z] - AccessPoint(k, l)[Z]) < test_length2)  )
							{
								VSub(String, AccessPoint(i, j), AccessPoint(k, l));
								VDot(length, String, String);
								if ((length < test_length) && (length != 0))
								{
									length = length / test_length;
									length = max(length, 0.15);
									ForceScaler = (cloth_sys->cloth->ks)*(1.0-length)/(length*length);
									VScaleEq(String, ForceScaler);
									VAddEq(Force, String);
								}
							}
						}
					}
				}
			}

			/* Update velocity */
			VScaleEq(Force, cloth_sys->dt);
			VScaleEq(AccessVel(i, j), cloth_sys->damping);
			VAddEq(AccessVel(i, j), Force);
			/* apply constraint */
			VEvaluateEq(AccessVel(i, j), AccessConstraint(i, j));
			VScale(Move, AccessVel(i, j), cloth_sys->dt);

			/* Collision dectection */
			if (cloth_sys->object != NULL)
			{
				VLength(length, Move);
				Detect_Collision(cloth_sys->object, AccessVel(i, j), cloth_sys->dt, AccessPoint(i, j), resultmp, cloth_sys->friction);
			}
			else
			{
				VAdd(resultmp, AccessPoint(i, j), Move);
			}
			Assign_Vector(AccessResul(i,j), resultmp);
		}
	}
}


/*****************************************************************************
* FUNCTION   ReadClothFile()
******************************************************************************/
static bool ReadClothFile(CLOTH_SYS *cloth_sys)
{
  IStream *fd;
  char line[MAX_STR];
	 SNGL x, y, z, vx, vy, vz;
	 int i,j;

	fd = New_IStream(cloth_sys->input, POV_File_Unknown);
	if ( fd == NULL)
	{
		Error("Cannot open input cloth file...");
		return false;
	}

    if ( fd->getline (line, MAX_STR).eof () )
    {
      Error( "Error reading cloth file header...(unexpected file end)");
		delete fd;
      return false;
	}
	
	if (sscanf(line, "%i, %i, %f, %f,", &(cloth_sys->cloth->nb1), &(cloth_sys->cloth->nb2),
	                                   &(cloth_sys->cloth->nlng), &(cloth_sys->cloth->ks)) != 4)
	{
		Error("Error reading cloth file header...");
		delete fd;
		return false;
	}

	if (cloth_sys->cloth->nb1 < 2) 
		cloth_sys->cloth->nb1 = 2;
	if (cloth_sys->cloth->nb2 < 2) 
		cloth_sys->cloth->nb2 = 2;
	if (cloth_sys->cloth->nlng <= 0) 
		cloth_sys->cloth->nlng = 0.0001;
	if (cloth_sys->cloth->ks < 0.1) 
		cloth_sys->cloth->ks = 0.1;

/*	Status_Info("n1 = %i, n2 = %i, nlng = %f, ks = %f\n", cloth->nb1, cloth->nb2,
	                                   cloth->nlng, cloth->ks);  */

	cloth_sys->cloth->points = (VECTOR *)POV_MALLOC(sizeof(VECTOR)*(cloth_sys->cloth->nb1)*(cloth_sys->cloth->nb2), "cloth points");
	cloth_sys->cloth->vel = (VECTOR *)POV_MALLOC(sizeof(VECTOR)*(cloth_sys->cloth->nb1)*(cloth_sys->cloth->nb2), "cloth points velocity");
	cloth_sys->cloth->constraint = (VECTOR *)POV_MALLOC(sizeof(VECTOR)*(cloth_sys->cloth->nb1)*(cloth_sys->cloth->nb2), "cloth constraints");

	for (i=0; i<cloth_sys->cloth->nb1; i++)
	{
		for (j=0; j<cloth_sys->cloth->nb2; j++)
		{
			*line=0;
		   fd->getline (line, MAX_STR);
			if (sscanf(line, "%f, %f, %f, %f, %f, %f,", &x, &y, &z, &vx, &vy, &vz) == 6)
			{
				AccessPoint(i, j)[X] = x;
				AccessPoint(i, j)[Y] = y;
				AccessPoint(i, j)[Z] = z;
				AccessVel(i, j)[X] = vx;
				AccessVel(i, j)[Y] = vy;
				AccessVel(i, j)[Z] = vz;
			}
			else
			{
		      Warning(0, "Error reading cloth points...");
				delete fd;
				return false;
			}
		}
	}

	/* constraints init. */
	for (i=0; i<(cloth_sys->cloth->nb1*cloth_sys->cloth->nb2); i++)
	{
		cloth_sys->cloth->constraint[i][X] = 1.0;
		cloth_sys->cloth->constraint[i][Y] = 1.0;
		cloth_sys->cloth->constraint[i][Z] = 1.0;
	}

/*	Status_Info("\nReading input cloth file constraints...\n"); */
	/* read constraints */
	j=0;
    while ( ! fd->eof() && fd->getline (line, MAX_STR) )
    {
		if (sscanf(line, "%i, %f, %f, %f,", &i, &vx, &vy, &vz) == 4)
		{
			if ( (i < 0) || (i > cloth_sys->cloth->nb1*cloth_sys->cloth->nb2))
				Error("Constrained index too high.");
			if (i<((cloth_sys->cloth->nb1)*(cloth_sys->cloth->nb2)))
			{
				cloth_sys->cloth->constraint[i][X] = vx;
				cloth_sys->cloth->constraint[i][Y] = vy;
				cloth_sys->cloth->constraint[i][Z] = vz;
				j++;
			}
		}
	}
	Debug_Info("%i constraints read\n", j);
	
	delete fd;
	return true;
}

/*****************************************************************************
* FUNCTION   WriteClothFile()
******************************************************************************/
static bool WriteClothFile( CLOTH_SYS *cloth_sys)
{
   DBL x, y, z, vx, vy, vz;
	int i,j;
	ClothFile = New_OStream(cloth_sys->output, POV_File_Unknown,false);
	if ( ClothFile == NULL)
	{
		Warning(0,"Cannot open output cloth file...");
		return false;
	}
	Debug_Info("Writing output cloth file header\n"); 
	ClothFile->printf( "%i, %i, %f, %f,\n", cloth_sys->cloth->nb1, cloth_sys->cloth->nb2,
	                                   cloth_sys->cloth->nlng, cloth_sys->cloth->ks);

   Debug_Info("Writing output cloth file points\n"); 
	for (i=0; i<cloth_sys->cloth->nb1; i++)
	{
		Do_Cooperate(1);
		for (j=0; j<cloth_sys->cloth->nb2; j++)
		{
			x = AccessPoint(i, j)[X];
			y = AccessPoint(i, j)[Y];
			z = AccessPoint(i, j)[Z];
			vx = AccessVel(i, j)[X];
			vy = AccessVel(i, j)[Y];
			vz = AccessVel(i, j)[Z];
			ClothFile->printf( "%.12f, %.12f, %.12f, %.12f, %.12f, %.12f,\n", x, y, z, vx, vy, vz);
		}
	}

	/* writing constraints */
	for (i=0; i<((cloth_sys->cloth->nb1)*(cloth_sys->cloth->nb2)); i++)
	{
		vx = cloth_sys->cloth->constraint[i][X];
		vy = cloth_sys->cloth->constraint[i][Y];
		vz = cloth_sys->cloth->constraint[i][Z];
		if ( (vx != 1.0 ) ||   (vy != 1.0 ) ||     (vz != 1.0 ) )
			ClothFile->printf( "%i, %.12f, %.12f, %.12f,\n", i, vx, vy, vz);
	}
    delete ClothFile;
    ClothFile=NULL;
	
	return true;
}



/*****************************************************************************
* FUNCTION   WriteClothMesh()
******************************************************************************/
static bool WriteClothMesh(CLOTH_SYS *cloth_sys)
{
 DBL x1, y1, z1, x2, y2, z2, x3, y3, z3;
 int j,  n2;
 VECTOR *normals=NULL;

	ClothMeshFile = New_OStream(cloth_sys->mesh_output, POV_File_Unknown,false);
	if ( ClothMeshFile == NULL)
	{
		Error("Cannot open output mesh file...");
		return false;
	}


	n1Counter = cloth_sys->cloth->nb1;
	n2 = cloth_sys->cloth->nb2;

	if (cloth_sys->smooth_mesh)
	{
		normals = (VECTOR *)POV_MALLOC(sizeof(VECTOR)*n1Counter*n2, "normal vectors");
		
		for (iCounter=0; iCounter<n1Counter; iCounter++)
		{
			for (j=0; j<n2; j++)
			{
				GetNormal(cloth_sys, iCounter, j, normals[iCounter*n2+j]);
			}
		}
	}

  Send_Progress("Writing output mesh triangles", PROGRESS_WRITING_CLOTHRAY_MESH);
	for (iCounter=0; iCounter<n1Counter; iCounter++)
	{
		Send_ProgressUpdate(PROGRESS_WRITING_CLOTHRAY_MESH);
          Do_Cooperate(0);
		
		for (j=0; j<n2; j++)
		{
          Do_Cooperate(2);
			x1 = AccessPoint(iCounter, j)[X];
			y1 = AccessPoint(iCounter, j)[Y];
			z1 = AccessPoint(iCounter, j)[Z];

			if (((j+1) < n2) && ((iCounter+1) < n1Counter))
			{
				x2 = AccessPoint(iCounter, j+1)[X];
				y2 = AccessPoint(iCounter, j+1)[Y];
				z2 = AccessPoint(iCounter, j+1)[Z];
				x3 = AccessPoint(iCounter+1, j)[X];
				y3 = AccessPoint(iCounter+1, j)[Y];
				z3 = AccessPoint(iCounter+1, j)[Z];

				if (cloth_sys->smooth_mesh)
				{
					DBL nx1, ny1, nz1, nx2, ny2, nz2, nx3, ny3, nz3;
					nx1 = normals[iCounter*n2+j][X];
					ny1 = normals[iCounter*n2+j][Y];
					nz1 = normals[iCounter*n2+j][Z];
					nx2 = normals[iCounter*n2+j+1][X];
					ny2 = normals[iCounter*n2+j+1][Y];
					nz2 = normals[iCounter*n2+j+1][Z];
					nx3 = normals[(iCounter+1)*n2+j][X];
					ny3 = normals[(iCounter+1)*n2+j][Y];
					nz3 = normals[(iCounter+1)*n2+j][Z];

					ClothMeshFile->printf((char*) format_smooth_triangle, x1, y1, z1, nx1, ny1, nz1,
					                                              x2, y2, z2, nx2, ny2, nz2,
					                                              x3, y3, z3, nx3, ny3, nz3);
				}
				else
				{
					ClothMeshFile->printf( (char*)format_triangle, x1, y1, z1, x2, y2, z2, x3, y3, z3);
				}
				if (cloth_sys->uv_mesh)
				{
					 ClothMeshFile->printf( (char*)format_uv_vectors, (SNGL)iCounter/(n1Counter-1), (SNGL)j/(n2-1), (SNGL)iCounter/(n1Counter-1), ((SNGL)j+1)/(n2-1), ((SNGL)iCounter+1)/(n1Counter-1), (SNGL)j/(n2-1));
				}
				 ClothMeshFile->printf( " }\n");

			}

			if ((j > 0) && (iCounter > 0))
			{
				x2 = AccessPoint(iCounter, j-1)[X];
				y2 = AccessPoint(iCounter, j-1)[Y];
				z2 = AccessPoint(iCounter, j-1)[Z];
				x3 = AccessPoint(iCounter-1, j)[X];
				y3 = AccessPoint(iCounter-1, j)[Y];
				z3 = AccessPoint(iCounter-1, j)[Z];

				if (cloth_sys->smooth_mesh)
				{
					DBL nx1, ny1, nz1, nx2, ny2, nz2, nx3, ny3, nz3;
					nx1 = normals[iCounter*n2+j][X];
					ny1 = normals[iCounter*n2+j][Y];
					nz1 = normals[iCounter*n2+j][Z];
					nx2 = normals[iCounter*n2+j-1][X];
					ny2 = normals[iCounter*n2+j-1][Y];
					nz2 = normals[iCounter*n2+j-1][Z];
					nx3 = normals[(iCounter-1)*n2+j][X];
					ny3 = normals[(iCounter-1)*n2+j][Y];
					nz3 = normals[(iCounter-1)*n2+j][Z];

					ClothMeshFile->printf( (char*)format_smooth_triangle, x1, y1, z1, nx1, ny1, nz1,
					                                              x2, y2, z2, nx2, ny2, nz2,
					                                              x3, y3, z3, nx3, ny3, nz3);
				}
				else
				{
					ClothMeshFile->printf( (char*)format_triangle, x1, y1, z1, x2, y2, z2, x3, y3, z3);
				}
				if (cloth_sys->uv_mesh)
				{
					ClothMeshFile->printf((char*) format_uv_vectors, (SNGL)iCounter/(n1Counter-1), (SNGL)j/(n2-1), (SNGL)iCounter/(n1Counter-1), ((SNGL)j-1)/(n2-1), ((SNGL)iCounter-1)/(n1Counter-1), (SNGL)j/(n2-1));
				}
				ClothMeshFile->printf( " }\n");
			}
		}
	}
	
	Send_ProgressUpdate(PROGRESS_WRITING_CLOTHRAY_MESH);
    Do_Cooperate(0);

    delete ClothMeshFile;
    ClothMeshFile=NULL;

	return true;
}

/*****************************************************************************
* FUNCTION   GetNormal()
******************************************************************************/
static void GetNormal(CLOTH_SYS *cloth_sys, int i, int j, VECTOR Result)
{
 int k, n1, n2, indi1, indi2, indj1, indj2;
 VECTOR temp1, temp2, temp3;

	Make_Vector(Result, 0, 0, 0);
	n1 = cloth_sys->cloth->nb1;
	n2 = cloth_sys->cloth->nb2;

	for (k=0; k<8; k++)	
	{
		indi1 = i+La1[k][0];
		indj1 = j+La1[k][1];
		indi2 = i+La1[(k+1)%8][0];
		indj2 = j+La1[(k+1)%8][1];
		if ( (indi1<n1) &&
		     (indi1>=0) &&
		     (indj1<n2) &&
		     (indj1>=0) &&
		     (indi2<n1) &&
		     (indi2>=0) &&
		     (indj2<n2) &&
		     (indj2>=0) )
		{
			VSub(temp1, AccessPoint(i,j), AccessPoint(indi1, indj1));
			VSub(temp2, AccessPoint(i,j), AccessPoint(indi2, indj2));
			VCross(temp3, temp1, temp2);
			VAddEq(Result, temp3);
		}
	}
	VNormalizeEq(Result);
}


/*****************************************************************************
*
* FUNCTION   Make_Cloth_Simulation()
*
*   
*
* INPUT    
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*
* DESCRIPTION
*
*
* CHANGES
*
*
******************************************************************************/
static int Make_Cloth_Simulation(CLOTH_SYS *cloth_sys)
{
 	VECTOR *resul, *temp;

	/* ================== */
	/* Reading input file */
	/* ================== */
	cloth_sys->cloth = (CLOTH *)POV_MALLOC(sizeof(CLOTH), "cloth structure");
	cloth_sys->cloth->points = NULL;
	cloth_sys->cloth->vel = NULL;
	cloth_sys->cloth->constraint = NULL;

	if (!ReadClothFile(cloth_sys))
		return false;


	resul = (VECTOR *)POV_MALLOC(sizeof(VECTOR)*(cloth_sys->cloth->nb1)*(cloth_sys->cloth->nb2), "result cloth points");


	/* ====================== */
	/* running the simulation */
	/* ====================== */
 	current_sys = cloth_sys;
 	if ( cloth_sys->nbit)
        {
          Send_Progress("Cloth Simulation", PROGRESS_RUNNING_CLOTH_SIMULATION);
          START_TIME

          for(n1Counter=0; n1Counter<cloth_sys->nbit; n1Counter++)
          {
            Send_ProgressUpdate(PROGRESS_RUNNING_CLOTH_SIMULATION);
            Do_Cooperate(0);

            /* update positions */
            Iterate_Cloth_Simulation(cloth_sys, resul);

            /* swap buffers */
            temp = cloth_sys->cloth->points;
            cloth_sys->cloth->points = resul;
            resul = temp;
          }
 	
          //make sure 100% is displayed as well
          Send_ProgressUpdate(PROGRESS_RUNNING_CLOTH_SIMULATION, 0);
          Do_Cooperate(0);

          STOP_TIME
          tcloth = TIME_ELAPSED;
          tcloth_frame = tcloth;
          tcloth_total += tcloth;
	}


	POV_FREE(resul);

	/* =================== */
	/* Writing output file */
	/* =================== */

/*	Status_Info("\nOpening output cloth file...\n"); */
	if (cloth_sys->output != NULL)
	{
		if (  !WriteClothFile(cloth_sys))
			return false;
	}

	if (cloth_sys->mesh_output != NULL)
	{
		return WriteClothMesh(cloth_sys);
	}

	return true;
}
int GetClothrayStat(POVMSType a)
{
	switch(a)
	{
		case kPOVAttrib_ClothRayMeshWrite:
			return (int)(iCounter*100/n1Counter);
		case kPOVAttrib_ClothRaySimulation:
			return (int)(n1Counter*100/current_sys->nbit);
	}
	return 0;
}

END_POV_NAMESPACE
#endif
