/*bezier.c*/

#define MAXPTS 2000		/* sets max # points to read */
#define MAXBEZ 2000		/* sets max # points to write */

/*#include <m68881.h>*/		/* include for math coprocessor */
#include <math.h>
#include <stdio.h>
#include <string.h>

unsigned short int cnt=0,numpnts=0;
static double outarr[MAXBEZ][2];


main(argc,argv)
char *argv[];
{

double array[MAXPTS][2];
unsigned short int q,pstat=0,wb=1;

FILE *infile=stdin,*outfile=stdout,*fopen();

	if(argc==1)
		{
		fprintf(stderr," Use: %s [-nod] <# BezPoints> [infile][outfile]\n",argv[0]);
		exit(20);
		}
	for(q=1;q<argc;q++)
		{
		if(strncmp(argv[q],"-nod",4)==0)	/* this option allows files without DATA statement */
			{
			wb=0;
			continue;
			}
		else if(pstat==0)
			{
			if(sscanf(argv[q],"%hd",&numpnts)!=1)	/* number of points to generate */
				{
				fprintf(stderr," Strange number!\n");
				exit(20);
				}
			else ++pstat;
			}
		else if(pstat==1)
			{
			if((infile=fopen(argv[q],"r"))==NULL)
				{
				fprintf(stderr," %s cant open file %s for read\n",argv[0],argv[q]);
				exit(20);
				}
			else ++pstat;
			}
		else if(pstat==2)
			{
			if((outfile=fopen(argv[q],"w"))==NULL)
				{
				fprintf(stderr," %s cant open file %s for write\n",argv[0],argv[q]);
				exit(20);
				}
			else ++pstat;
			}
		}


	ReadData(infile,array,&cnt,wb);
	if(cnt==0)
		{
		fprintf(stderr," %s problem reading data\n",argv[0]);
		exit(20);
		}
	else fprintf(stderr," %s read %hd data pairs\n",argv[0],cnt);

	if(infile)fclose(infile);

	CalcBez(array,outarr,cnt,&numpnts);

	WriteData(outfile,outarr,numpnts,wb);

	fprintf(stderr,"\n %s wrote %hd data pairs\n",argv[0],numpnts);

	if(outfile)fclose(outfile);
	exit(0);
}



ReadData(ifp,data,count,byers)	/* this routine reads xy data pairs into array */

unsigned short int *count,byers;
double data[MAXPTS][2];
FILE *ifp;

{

register unsigned short int ctr;
char cm,key[50];

	while(byers)				/* if true search for keyword DATA */
		{					/* Byers is the twisted person who requires this format */

		if(fscanf(ifp,"%s",key)!=1)
			{
			fprintf(stderr," keyword DATA not found!\n");
			exit(20);
			}
		if(strcmp("DATA",key)==0)break;
		}

	for(ctr=0;ctr<MAXPTS;ctr++)
		if(fscanf(ifp,"%lf%c%lf",&data[ctr][0],&cm,&data[ctr][1])!=3)break;
	if(ctr==MAXPTS)fprintf(stderr," MAXPTS reached exiting ReadData\n");
	*count=ctr;
	return(0);
}


WriteData(ofp,bezdat,numb,byers)	/* this function writes array to file da */

unsigned short int numb,byers;
double bezdat[MAXBEZ][2];
FILE *ofp;

{
register unsigned short int q;

	if(byers)fprintf(ofp," DATA\n");
	for(q=0;q<numb;q++)fprintf(ofp,"%-14le%15le\n",bezdat[q][0],bezdat[q][1]);
	return(0);
}




CalcBez(orgarr,bezarr,orgcnt,bezcnt)	/* this function calculates an array of */
								/* bezier points from the orginal array */
unsigned short int *bezcnt,orgcnt;
double orgarr[MAXPTS][2], bezarr[MAXBEZ][2];

{
register unsigned short int i,j,k,n;

double t,step,mr[MAXPTS];
double coef[MAXPTS];


	if(*bezcnt>MAXBEZ)
		{
			*bezcnt=MAXBEZ;
			fprintf(stderr," bezpnts set to MAXBEZ\n");
		}

	CalcCoef(orgcnt,coef);

	n=orgcnt-1;
	k=1;

	step=1.0/(*bezcnt-1);

	for(j=0;j<2;j++)bezarr[0][j]=orgarr[0][j];

	for(t=step;t<=1.0;t+=step)		/* The big loop */
		{
		fprintf(stderr,"  t == %.3lf\r",t);

		for(i=0;i<=n;i++)
			{
			mr[i]=coef[i]*pow((1.0-t),(double)(n-i))*pow(t,(double)i);
			for(j=0;j<2;j++)bezarr[k][j]+=mr[i]*orgarr[i][j];
			}
		k++;
		}

	for(j=0;j<2;j++)bezarr[k][j]=orgarr[n][j];k++;

	return(0);
}



CalcCoef(numpts,coarr)	/* this function calcs coeficients from Pascals triangle */

unsigned short int numpts;
double coarr[MAXPTS];

{
register unsigned short int y,z;
static double temparr[MAXPTS];

	coarr[0]=coarr[1]=temparr[0]=1.0;

	for(y=1;y<=numpts;y++)
		{
		for(z=1;z<y;z++)temparr[z]=coarr[z-1]+coarr[z];
		for(z=1;z<=y;z++)coarr[z]=temparr[z];
		}
	return(0);
}