/*
Name:    trimake - Performs Delaunay triangulation of a set of points.
                   Coincident points ate rejected.

Date of the version:   Pisa, 10-JAN-91

Authors: G. Macedonio and M.T. Pareschi

Modified by R. Ambroziak, USGS/OEMG, Reston, VA 703-648-6168

Usage:
       trimake ptsfile trifile

where:
       ptsfile     -  input file vertices
       trifile     -  output file triangle
Note:
       MAXPTS = Maximum number of points
       MAXTRI = Maximum number of triangles
                set MAXTRI = 2*MAXPTS
       MAXOPT = Maximum number of recursive calls
                of the routine that performs deluaney
                optimization (optim).  This numbe is
                always less than MAXTRI/3.  You can set
                any value.  Greater numbers assure
                that the deluaney triangulation is well
                performed but the program may require
                too much RAM memory.  Usually the number
                of recursive calls is less than ten.

*/
#include <stdio.h>
#include <math.h>
#include <string.h>

/* 
	MAXTRI = 2*MAXPTS therefore 1 MAXPTS = 34 bytes of memory */
#define MAXPTS 12001
#define MAXTRI 24002
#define MAXOPT    50
#define DBERR 1.e-15

extern void exit();    			/* exits main program */
static int jnd[3]={1,2,0};		/* Auxiliary vector */
static int jrd[3]={2,0,1};		/* Auxiliary vector */
int huge vt[3][MAXTRI];				/* list of verticies */
int huge nt[3][MAXTRI];				/* list of neighbor triangles */
int dcd[MAXPTS];					/* list of ordered points */
int npts; 							/* Number of vertices */
float px[MAXPTS],py[MAXPTS];	/* Coordiates of vertices */
int ntri;							/* number of triangles -1 */
int vtm[3],ntm[3],ntr[3];		/* Auxiliary pointers */
int loptim=0;						/* Optimization level */
float Zscale=1.0;					/* Vertical exageration */
float Thickness=0;
int K[3]={0,2,1};

int comment_wvf(FILE *fp,char name[],float scale,
		float xmin,float xmid,float xmax,
		float ymin,float ymid,float ymax,
		float zmin,float zmid,float zmax);

main(int argc,char *argv[])
{
	int nptf,nptfs,nptfi;			/* Hull verticies (tot,up,low) */
	unsigned i,k,l,nl;				/* Auxiliary variables */
	int it;								/* Pointer to triangle */
	int nrej1;							/* Number of rejected points */
	int nrej2;							/* Number of coincident points */
	double xmin,xmax,ymin,ymax;	/* Min and max coordinates */
	double strip;						/* External frame fraction */
	double dxstrip,dystrip;			/* Dimensions of external frame */
	double px1,px2,py1,py2,aux;	/* Auxiliary variables */
	double m,q,x0;						/* Params straight line */
	int swapv();						/* Swaps two memory locations */
	int ckturn();						/* Check counterclockwise turn */
	int trifind();						/* Finds external triangle */
	int cklop();						/* Check L.O.P. */
	int optim();						/* Optimizes two triangles */
	int swaptr();						/* Swaps two diagonals */
	void qsort();						/* C-Language sort, if available */
	void fsort();						/* Quick sort */
	int xcmp();							/* Compares two abscissa */
	int bcmp();							/* Compares two bin numbers */
	double pow2();						/* squares a number */
	float xx,yy,zz,xbar=0.0,ybar=0.0,zbar=0.0,fnum=0.0,xscale,yscale,scale;
	float zmin,zmax,xmid,ymid,zmid;
	char string[100],name[100],inname[50];
	FILE *ptstrm;						/* Stream file list of vertices */
	FILE *trstrm;						/* Strean file list of triangles */
	FILE *fptxt;						/* Test output from thin */
	FILE *fpwvf;						/* wave front file */
	float zscale,thick;
	int mm;

	strip=0.05;
	if(argc<3)
	{
		printf("trimake:  Invalid number of arguments\n\n");
		printf("trimake [ptsfile] [trifile] <vertexag> </d> <thickness>\n\n");
		printf("where:\n");
		printf(
"   ptsfile   --  input file vertices (file must be of type '.v)'\n");
		printf(
"   trifile   --  output file triangle (files will be of type '.w' & '.wvf')\n");
		printf(
"\n   vertexag  --  vertical exaggeration (optional -- default = 1)\n");
		printf(
"         /d  --  indicator for a solid object with double sides (optional)\n");
		printf(
"  thickness  --  thickness of object\n");
		exit(1);
	}
	if(argc>3)
	{
		if(sscanf(argv[3],"%f",&zscale)==1)
			Zscale=zscale;
		for(i=3;i<argc-1;i++)
			if(strcmpi(argv[i],"/d")==0)
				if(sscanf(argv[i+1],"%f",&thick)==1)
				{
					Thickness=thick;
					printf(
				"Double thickness object requested -- Thickness = %0.3f\n",
						Thickness);
				}
	}
	strcpy(inname,argv[1]);
	for(k=0;k<strlen(inname);k++)
		if(inname[k]=='.')
			inname[k]='\0';
	sprintf(string,"%s.v",inname);
	if((ptstrm=fopen(string,"r"))==NULL)
	{
		printf("Cannot open '%s' to read vertices.\n\n",argv[1]);
		exit(1);
	}
	strcpy(name,argv[2]);
	for(k=0;k<strlen(name);k++)
		if(name[k]=='.')
			name[k]='\0';
	sprintf(string,"%s.wvf",name);
	if((fpwvf=fopen(string,"wt"))==NULL)
	{
		printf("Cannot open '%s' to write wavefront file.\n\n",string);
		exit(1);
	}
	sprintf(string,"%s.txt",inname);
	fptxt=fopen(string,"rt");
	if(fptxt)
		while(fgets(string,100,fptxt))
			fprintf(fpwvf,"%s",string);
	fprintf(fpwvf,"#\n");
	fprintf(fpwvf,"# USGS/OEMG PC 3-D Visualization System\n");
	fprintf(fpwvf,"# Phone 703-648-6168 of 703-648-6481");
	fprintf(fpwvf,"#\n");
	fprintf(fpwvf,"# file %s.wvf created by program 'trimake.c'\n",name);
	sprintf(string,"%s.w",inname);
	if((trstrm=fopen(string,"wt"))==NULL)
	{
		printf("Cannot open '%s' to write triangles.\n\n",argv[2]);
		exit(1);
	}
	/* read 1st input point */
	printf("Reading data ...\n");
	npts=0;
	if((i=fscanf(ptstrm,"%f%f%f",&px[npts],&py[npts],&zz))!=EOF)
	{
		fnum+=1.0;
		xbar+=px[npts];
		ybar+=py[npts];
		zbar+=zz;
		xmax=xmin=px[npts];
		ymax=ymin=py[npts];
		zmax=zmin=zz;
		npts+=1;
	}
	else
		exit(1);
	/* read rest of points */
	while((i=fscanf(ptstrm,"%f%f%f",&px[npts],&py[npts],&zz))!=EOF)
	{
		fnum+=1.0;
		xbar+=px[npts];
		ybar+=py[npts];
		zbar+=zz;
		xmin=(xmin<px[npts])?xmin:px[npts];
		ymin=(ymin<py[npts])?ymin:py[npts];
		zmin=(zmin<zz)?zmin:zz;
  		xmax=(xmax>px[npts])?xmax:px[npts];
		ymax=(ymax>py[npts])?ymax:py[npts];
		zmax=(zmax>zz)?zmax:zz;
		npts+=1;
		if(npts==MAXPTS)
		{
			printf("trimake:  Number of points are more than %d\n",(int)MAXPTS);
			exit(1);
		}
	}
	printf("Number of input points = %d\n",npts);
	if(npts<3)
		exit(1);
	printf("Printing vertices ... \n");
	rewind(ptstrm);
	if(fnum<=0.0)
	{
		printf("fnum invalid...\n\n");
		exit(0);
	}
	xbar/=fnum;
	ybar/=fnum;
	zbar/=fnum;
	xmid=(xmax+xmin)/2.0;
	ymid=(ymax+ymin)/2.0;
	zmid=(zmax+zmin)/2.0;
	xscale=xmax-xmin;
	yscale=ymax-ymin;
	scale=(xscale>yscale)?xscale:yscale;
	scale/=9.0;
	printf("scale = %f\n",scale);
	printf("xmid = %f ymid = %f zmid= %f\n",xmid,ymid,zmid*1000.0);
	printf("\n Do you want to change them?  (y or n)\n\n");
	if(getch()=='y')
	{
		do
		{
			printf("Give new scale:  ");
			scanf("%f",&scale);
			printf("Give new xmid:  ");
			scanf("%f",&xmid);
			printf("Give new ymid:  ");
			scanf("%f",&ymid);
			printf("Give new zmid:  ");
			scanf("%f",&zmid);
			zmid/=1000.0;
			printf("scale = %f\n",scale);
			printf("xmid = %f ymid = %f zmid= %f\n",xmid,ymid,zmid*1000.0);
			printf("\n Do you want to change them again?  (y or n)\n\n");
		}while(getch()=='y');
	}
	comment_wvf(fpwvf,argv[1],scale,xmin,xmid,xmax,ymin,ymid,ymax,
				zmin,zmid,zmax);
	while((i=fscanf(ptstrm,"%f%f%f",&xx,&yy,&zz))!=EOF)
	{
		fprintf(fpwvf,"v %f %f %f\n",
			(xx-xmid)/scale,(yy-ymid)/scale,(zz-zmid)/scale*Zscale*1000.0);
		if(Thickness>0.0)
				fprintf(fpwvf,"v %f %f %f\n",
					(xx-xmid)/scale,(yy-ymid)/scale,
					(zz-zmid)/scale*Zscale*1000.0-Thickness);
	}

	/* Compute external frame */
	dxstrip=strip*(xmax-xmin);
	dystrip=strip*(ymax-ymin);
	px1=xmin-dxstrip;
	py1=ymin-dystrip;
	aux=npts;
	nl=1+pow(aux,(double)0.25);
	px2=(xmax-xmin+2*dxstrip)/nl;
	py2=(ymax-ymin+2*dystrip)/nl;
	for(i=0;i<npts;i++)
	{
		int nx,ny,k1,k2;

		nx=1+(px[i]-px1)/px2;
		ny=1+(py[i]-py1)/py2;
		k1=(ny%2==0)?1:-1;
		k2=(1-k1)/2;
		vt[0][i]=nl*ny-k1*nx-k2*(nl+1)+1;
		vt[1][i]=i;
		vt[2][i]=i;
	}
	/* sort bin number by increasing abscissa */
	printf("Sorting ...1\n");
	qsort((char *)(vt[2]),npts,sizeof(vt[2][0]),xcmp);
	printf("Sorting ...2\n");
	qsort((char *)(vt[1]),npts,sizeof(vt[1][0]),bcmp);
	/* Creates the vortex hull */
	printf("Creating convex hull ...\n");
	x0=px[vt[2][0]];
	q=py[vt[2][0]];
	m=(py[vt[2][npts-1]]-q)/(px[vt[2][npts-1]]-x0);
	nptfi=1;
	nptfs=npts;
	nt[0][0]=vt[2][0];
	for(i=1;i<npts;i++)
	{
		int ind;

		ind=vt[2][i];
		if((py[ind]-q)>=m*(px[ind]-x0))
		{
			nptfs--;
			nt[0][nptfs]=ind;
		}
		else
		{
			nt[0][nptfi]=ind;
			nptfi++;
		}
	}
	nptf=1;
	for(i=2;i<npts;i++)
	{
		while((ckturn(nt[0],nptf,i)>0)&&(nptf>=1))
			nptf--;
		nptf++;
		nt[0][nptf]=nt[0][i];
	}
	while((ckturn(nt[0],nptf,0)>0)&&(nptf>=1))
		nptf--;
	/* builds ordered list of vertices */
	for(i=0;i<=nptf;i++)
	{
		int ind;

		ind=nt[0][i];
		dcd[i]=ind;
		vt[0][ind]=0;
	}
	k=nptf;
	for(i=0;i<npts;i++)
	{
		int  ind;

		ind=vt[1][i];
		if(vt[0][ind]!=0)
		{
			k++;
			dcd[k]=ind;
		}
	}
	/* Beginning triangulation procedure */
		/* Builds first triangle */
	ntri=nptf-2;
	for(i=0;i<=ntri;i++)
	{
		vt[0][i]=dcd[0];
		vt[1][i]=dcd[i+1];
		vt[2][i]=dcd[i+2];
		nt[0][i]=i+1;
		nt[1][i]=i-1;
		nt[2][i]= -1;
		optim(i,2);
	}
	nt[0][ntri]= -1;
	/* Beginning of loop */
	printf("Beginning triangulation ...\n");
	it=0;
	nrej2=nrej1=0;
	for(i=nptf+1;i<npts;i++)
	{
		int it1;

		nl=dcd[i];
		it1=trifind(nl,it);
		if(it1>=0)
		{
			it=it1;
			ntr[0]=it;
			ntr[1]=ntri+1;
			ntr[2]=ntri+2;
			for(k=0;k<3;k++)
			{
				vtm[k]=vt[k][it];
				ntm[k]=nt[k][it];
			}
			for(k=0;k<3;k++)
			{
				int ntrk,knd,krd;

				ntrk=ntr[k];
				knd=jnd[k];
				krd=jrd[k];
				vt[0][ntrk]=nl;
				vt[1][ntrk]=vtm[k];
				vt[2][ntrk]=vtm[knd];
				nt[0][ntrk]=ntr[knd];
				nt[1][ntrk]=ntr[krd];
				nt[2][ntrk]=ntm[knd];
			}
			for(k=1;k<3;k++)
			{
				int ind;

				ind=nt[2][ntr[k]];
				if(ind>=0)
				{
					for(l=0;l<3;l++)
					{
						if(nt[l][ind]==ntr[0])
						{
							nt[l][ind]=ntr[k];
						}
					}
				}
			}
			ntri+=2;
			/* Local optimization Procedure (LOP) */
			optim(ntr[0],0);
			optim(ntr[1],0);
			optim(ntr[2],0);
		}
		else
		{
			nrej1++;
			if(it1==-2)
				nrej2++;
		}
	}
	printf("Number of coincident points = %d\n",nrej2);
	printf("  Number of rejected points = %d\n",nrej1);
	if(Thickness==0.0)
		printf("        Number of triangles = %d\n",ntri+1);
	else
		printf("        Number of triangles = %ld\n",(long)(ntri+1)*(long)2);
/*
	PRINTING OF TRIANGLES LIST
*/
	printf("Writing data ...\n");
	for(i=0;i<=ntri;i++)
	{
		fprintf(fpwvf,"f  ");
		if(Thickness==0.0)
		{
			for(k=0;k<3;k++)
			{
				fprintf(trstrm," %d",vt[k][i]);
				fprintf(fpwvf," %d",vt[k][i]+1);
			}
			for(k=0;k<3;k++)
				fprintf(trstrm," %d",nt[k][i]);
			fputc('\n',trstrm);
			fputc('\n',fpwvf);
		}
		else
		{
			for(k=0;k<3;k++)
			{
				fprintf(trstrm," %ld",(long)vt[k][i]*(long)2);
				fprintf(fpwvf," %ld",(long)vt[k][i]*(long)2+(long)1);
			}
			for(k=0;k<3;k++)
			{
				if(nt[k][i]>=0)
					fprintf(trstrm," %ld",(long)nt[k][i]*(long)2);
				else
					fprintf(trstrm," -1");
			}
			fputc('\n',trstrm);
			fputc('\n',fpwvf);
			fprintf(fpwvf,"f  ");
			for(k=0;k<3;k++)
			{
				fprintf(trstrm," %ld",(long)vt[K[k]][i]*(long)2+(long)1);
				fprintf(fpwvf," %ld",(long)vt[K[k]][i]*(long)2+(long)2);
			}
			fputc('\n',fpwvf);
			for(k=0;k<3;k++)
			{
				if(nt[K[k]][i]>=0)
					fprintf(trstrm," %ld",(long)nt[K[k]][i]*(long)2+(long)1);
				else
					fprintf(trstrm," -1");
			}
			fputc('\n',trstrm);
			for(k=0;k<3;k++)
			{
				if(nt[k][i]<0)		/* edge triangle */
				{
					mm=k-1;
					if(mm<0)
						mm=2;
					fprintf(fpwvf,"f    %ld %ld %ld\n",
						(long)vt[k][i]*(long)2+(long)1,
						(long)vt[mm][i]*(long)2+(long)1,
						(long)vt[mm][i]*(long)2+(long)2);
					fprintf(fpwvf,"f    %ld %ld %ld\n",
						(long)vt[k][i]*(long)2+(long)1,
						(long)vt[mm][i]*(long)2+(long)2,
						(long)vt[k][i]*(long)2+(long)2);
				}
			}
		}
	}
	exit(0);
}

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

int comment_wvf(FILE *fp,char name[],float scale,
		float xmin,float xmid,float xmax,float ymin,float ymid,float ymax,
		float zmin,float zmid,float zmax)

{
	fprintf(fp,"# input file is '%s'\n",name);
	fprintf(fp,"# Scale = %f km/in  (xyz below in inches)\n",scale);
	fprintf(fp,"# Midpoints:\n");
	fprintf(fp,"# X0 = %11.5f km\n",xmid);
	fprintf(fp,"# Y0 = %11.5f km\n",ymid);
	fprintf(fp,"# Z0 = %11.5f m\n",zmid*1000.0);
	fprintf(fp,"# Before adjustment and scaling (real data values):\n");
	fprintf(fp,"# X min = %10.4f km  X max = %10.4f km  X range = %10.4f km\n",
		xmin,xmax,xmax-xmin);
	fprintf(fp,"# Y min = %10.4f km  Y max = %10.4f km  Y range = %10.4f km\n",
		ymin,ymax,ymax-ymin);
	fprintf(fp,"# Z min = %10.4f m   Z max = %10.4f m   Z range = %10.4f m\n",
		zmin*1000.0,zmax*1000.0,(zmax-zmin)*1000.0);
	fprintf(fp,"# After adjustment and scaling ('.wvf' file values):\n");
	fprintf(fp,"# x min = %f  x max = %f  x range = %f\n",
		(xmin-xmid)/scale,(xmax-xmid)/scale,(xmax-xmin)/scale);
	fprintf(fp,"# y min = %f  y max = %f  y range = %f\n",
		(ymin-ymid)/scale,(ymax-ymid)/scale,(ymax-ymin)/scale);
	if(Zscale>0.0)
		fprintf(fp,"# z min = %f  z max = %f  z range = %f\n",
			(zmin-zmid)/scale*Zscale,(zmax-zmid)/scale*Zscale,
			(zmax-zmin)/scale*Zscale);
	else
		fprintf(fp,"# z min = %f  z max = %f  z range = %f\n",
			(zmax-zmid)/scale*Zscale*1000.0,(zmin-zmid)/scale*Zscale*1000.0,
			(zmax-zmin)/scale*Zscale*1000.0);
	fprintf(fp,"# Vertical Exaggeration = %f\n",Zscale);
	fprintf(fp,"# x = (X - X0)/Scale\n");
	fprintf(fp,"# y = (Y - Y0)/Scale\n");
	fprintf(fp,"# z = (Z - Z0)/Scale*Vert.Exagg./1000.0\n");
	if(Thickness>0.0)
		fprintf(fp,"# Double thickness object -- Thickness = %0.3f\n",
			Thickness);
	fprintf(fp,"#\n");
	fprintf(fp,
"# use programs 'imvis.c' and 'latlon.c' to make '.raw' files\n");
	fprintf(fp,
"# to create input for this program use 'thinraw.c' or 'thinpts.c'.\n");
	fprintf(fp,"# use programs 'waveprep.c' and 'wave.c' to view data.\n");
	fprintf(fp,"#\n");
}





/**********************************************************************
**
**		swaps 2 memory integer locations
**
**********************************************************************/

int swapv(int pv[],int i1,int i2)

{
	int pt;

	pt=pv[i1];
	pv[i1]=pv[i2];
	pv[i2]=pt;
	return(1);
}

/**********************************************************************
**
**	checks if 3 vertices are:
**
**		alligned               [ 0]
**    counterclockwise turn  [-1]
**		clockwise turn         [ 1]
**		coincident points      [ 2]
**
**********************************************************************/

int ckturn(int pv[],int i2,int i3)

{
	double aux;
	int i1;

	i1=i2-1;
	aux=(py[pv[i2]]-py[pv[i1]])*(px[pv[i3]]-px[pv[i1]])-
			(py[pv[i3]]-py[pv[i1]])*(px[pv[i2]]-px[pv[i1]]);
	if(aux>DBERR)
		return(1);
	if(aux<-DBERR)
		return(-1);
	aux=(px[pv[i3]]-px[pv[i2]])*(px[pv[i2]]-px[pv[i1]])+
			(py[pv[i3]]-py[pv[i2]])*(py[pv[i2]]-py[pv[i1]]);
	if(aux<=-DBERR)
		return(1);
	if(aux>=DBERR)
		return(0);
	return(2);
}


/**********************************************************************
**
**	Finds a triangle containing the new vertice pointed by  ind
**
**		tr = starting triangle
**		trifind returns the number of the triangle or:
**		-1 if the triangle is not fond
**		-2 if the vertice coincides with an existing one
**
**********************************************************************/

int trifind(int ind,int tr)

{
	int i,aux,newtr;

	newtr=tr;
	do
	{
		aux=0;
		for(i=0;i<=2;i++)
		{
			double area,x21,x31,y21,y31;
			int i1,ind1,ind2;

			i1=jnd[i];
			ind1=vt[i][newtr];
			ind2=vt[i1][newtr];
			x31=(px[ind]-px[ind1]);
			x21=(px[ind2]-px[ind1]);
			y31=(py[ind]-py[ind1]);
			y21=(py[ind2]-py[ind1]);
			area=y31*x21-y21*x31;
			if(area<0.0)
			{
				newtr=nt[i1][newtr];
				if(newtr==-1)
					return(newtr);
				aux=1;
				break;
			}
			if(area==0.0)
				if(x31==0.0&&y31==0.0)
					return(-2);
		}
	}while(aux==1);
	return(newtr);
}

/**********************************************************************
**
**	checks if 4 points satisfy the local optimization
**	cklop = 1 if local optimization is satisfied
**	cklop = 0 if local optimization is not satisfied
**
**********************************************************************/

int cklop(int i1,int i2,int i3,int i4)

{
	double x2sq,y2sq,a,b,c,d,e,f,xc,yc,rsq;

	x2sq=pow2(px[i2]);
	y2sq=pow2(py[i2]);
	a=2*(px[i2]-px[i1]);
	b=2*(py[i2]-py[i1]);
	c=y2sq-pow2(py[i1])+x2sq-pow2(px[i1]);
	d=2*(px[i3]-px[i2]);
	e=2*(py[i3]-py[i2]);
	f=pow2(py[i3])-y2sq+pow2(px[i3])-x2sq;
	x2sq=d*b-e*a;
	rsq=pow2(x2sq);
	if(rsq>0.0)
	{
		xc=(f*b-c*e)/x2sq;
		yc=(d*c-f*a)/x2sq;
		rsq=pow2(xc-px[i1])+pow2(yc-py[i1]);
		if((pow2(py[i4]-yc)+pow2(px[i4]-xc))>rsq)
			return(1);
	}
	return(0);
}


/**********************************************************************
**
**	optimizes 2 triangles according to L.O.P.
**
**********************************************************************/

int optim(int it,int k1)

{
	int xnear,i,i1,k2,k3;

	if(loptim==MAXOPT)
		return(1);
	k2=jnd[k1];
	k3=jrd[k1];
	xnear=nt[k3][it];
	if(xnear==-1)
		return(2);
	loptim++;
	for(i=0;i<3;i++)
		if(nt[i][xnear]==it)
			break;
	if(i==3)
	{
		printf("ERROR IN OPTIM\n");
		exit(0);
	}
	i1=jnd[i];
	if(cklop(vt[k3][it],vt[k1][it],vt[k2][it],vt[i1][xnear])==0)
	{
		swaptr(it,xnear,k1,i);
		optim(it,k1);
		optim(xnear,i);
	}
	loptim--;
	return(3);
}


/**********************************************************************
**
**	Swaps the diagonals of a quadrilateral
**
**********************************************************************/

int swaptr(int t1,int t2,int k1,int k2)

{
	int save,i,ind1,ind2;

	vt[jrd[k1]][t1]=vt[jnd[k2]][t2];
	vt[k2][t2]=vt[k1][t1];
	save=nt[k1][t1];
	nt[k1][t1]=t2;
	nt[jrd[k1]][t1]=nt[jnd[k2]][t2];
	nt[jnd[k2]][t2]=t1;
	nt[k2][t2]=save;
	ind1=nt[jrd[k1]][t1];
	ind2=nt[k2][t2];
	for(i=0;i<3;i++)
	{
		if(ind1!=-1&&nt[i][ind1]==t2)
			nt[i][ind1]=t1;
		if(ind2!=-1&&nt[i][ind2]==t1)
			nt[i][ind2]=t2;
	}
	return(0);
}

/**********************************************************************
**
**	returns the square of a number
**
**********************************************************************/

double pow2(double x)

{
	return(x*x);
}

/**********************************************************************
**
**		routine called by qsort
**		compares two abscissa
**		if they are equal it compares the ordinates
**
**********************************************************************/

int xcmp(int *i1,int *i2)

{
	double aux;

	aux=px[*i1]-px[*i2];
	if(aux>0.0)
		return(1);
	if(aux<0.0)
		return(-1);
	aux=py[*i1]-py[*i2];
	if(aux>0.0)
		return(1);
	if(aux<0.0)
		return(-1);
	return(0);
}

/**********************************************************************
**
**		routine called by qsort
**		compares two bin numbers
**
**********************************************************************/

int bcmp(int *i1,int *i2)

{
	return(vt[0][*i1]-vt[0][*i2]);
}