/*
	Surface fitting program,  creates .z file from data points
*/

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define false (0)
#define true (!false)
#define gprint printf
int load_data(char *s);
int getstr(char *s);
#ifdef unix
#define VAXC 1
#endif
#ifdef VAXC
#define farfree free
#define farmalloc malloc
#else
#include <alloc.h>
#endif
typedef float real;
int idbvip_(long *md, long *ncp, long *ndp, real *xd, real *yd, real *zd
	,long *nip, real *xi, real *yi, real *zi, long *iwk, real *wk);
int setminmax(double v, double *min, double *max);
int sort_data(int npnts,float *xd,float *yd,float *zd);
int getrange(char *s, double  *min, double  *max, double  *step);
int xyz_alloc;
double getf(void);
float *pntxyz;
FILE *fp;
long npnts;
main()
{
	long md,ncp,ndp;
	/* real xd[8],yd[8],zd[8]; */
	char bfile[80],dfile[80],xrange[80];
	real *xd, *yd, *zd;
	double xmin,xmax,xstep,ymin,ymax,ystep,x,y;
	real *rx, *ry, *rz;
	long nx,ny;
	long nrp;
	/* real xi[888],yi[888],zi[888]; */
	long *iwk;
	real *wk;
	long i,j,k;
	
	ymin = xmin = 10e10;
	xmax = ymax = -xmin;

	printf("Data file conaining x y z data ? "); getstr(dfile);
	load_data(dfile);
	printf("Read %ld points from file {%s} \n",npnts/3,dfile);

	strcpy(bfile,dfile);
	if (strchr(bfile,'.')!=NULL) *strchr(bfile,'.') = 0;
	strcat(bfile,".z");
	printf("Will write data to {%s} \n",bfile);
	xd = malloc(sizeof(float) * (npnts/3+1));	
	yd = malloc(sizeof(float) * (npnts/3+1));	
	zd = malloc(sizeof(float) * (npnts/3+1));	
	printf("Number of points to use for contouring each point ? [3] "); ncp = getf();
	if (ncp==0) ncp = 3;
	
	for (j=0,i=0; i<npnts; i+=3) {
		xd[j] = pntxyz[i];
		yd[j] = pntxyz[i+1];
		zd[j] = pntxyz[i+2];
		setminmax(xd[j],&xmin,&xmax);
		setminmax(yd[j],&ymin,&ymax);
		j++;
	}
	farfree(pntxyz);
	npnts = npnts/3;
	sort_data(npnts,xd,yd,zd);

	for (i=0; i<npnts; i++) {
/*		printf("%g %g %g \n",xd[i],yd[i],zd[i]); */
		if (xd[i]==xd[i+1] && yd[i]==yd[i+1]) {
			printf("Duplicate points %g %g %g \n"
				,xd[i],yd[i],zd[i]);
		}
	}

	xstep = (xmax-xmin)/15.0;
	ystep = (ymax-ymin)/15.0;
	printf("Range of output x values [%g,%g,%g] ? ",xmin,xmax,xstep); getstr(xrange);
	getrange(xrange,&xmin,&xmax,&xstep);
	printf("Range of output y values [%g,%g,%g] ? ",ymin,ymax,ystep); getstr(xrange);
	getrange(xrange,&ymin,&ymax,&ystep);

	xyz_alloc = 0;
	nx = (xmax-xmin)/xstep + 1;
	ny = (ymax-ymin)/ystep + 1;
	rx = malloc(ny * nx * sizeof(float));
	ry = malloc(nx * ny * sizeof(float));
	rz = malloc(nx * ny * sizeof(float));
	
	k = 0;
	for (j=0, y=ymin; j<ny; y+= ystep, j++) {
		for (i=0, x=xmin; i<nx; x+= xstep, i++) {
			rx[k] = x;
			ry[k++] = y;
		}
	}	
	ndp = npnts;
	md = 1;
	nrp = nx*ny;
	k = 0;

	i = 27+ncp;
	if (i<31) i = 31;
	i = (i * ndp + nrp) * sizeof(*iwk) ;
	j = 8*ndp*sizeof(*wk);
	iwk = malloc(i);
	wk = malloc(j);
	printf("nx %ld, ny %ld, Work space iwk=%ld bytes wk=%ld bytes \n",nx,ny,i,j);
	if (iwk==NULL || wk == NULL) {
		printf("Unable to allocate enough workspace iwk=%ld wk=%ld  \n",i,j);
		exit(1);
	}

	idbvip_(&md,&ncp,&ndp,xd,yd,zd,&nrp,rx,ry,rz,iwk,wk);

	fp = fopen(bfile,"w");
	if (fp==NULL) {
		printf("Unable to Z file {%s} \n",bfile);
		exit(1);
	}

	fprintf(fp,"! nx %ld ny %ld xmin %g xmax %g ymin %g ymax %g\n",nx,ny,xmin,xmax,ymin,ymax);
	printf("Y = ");
	k = 0;
	for (j=0, y=ymin; j<ny; y+= ystep, j++) {
	  	printf("%g ",y);
		for (i=0, x=xmin; i<nx; x+= xstep, i++) {
			fprintf(fp,"%g ",rz[k++]);
		}
		fprintf(fp,"\n");
	}
	fclose(fp);
	printf("\n");
}
pnt_alloc(int size)
{
	void *d;
	if (size+10<xyz_alloc) return;
	size = size*2;
	d = farmalloc(size*sizeof(float));
	if (d==NULL) {
		gprint("Unable to allocate storage for data\n");
		exit(1);
	}
	if (xyz_alloc>0) memcpy(d,pntxyz,xyz_alloc*sizeof(float));
	xyz_alloc = size;
	pntxyz = d;
}
static char buff[2001];
load_data(char *fname)
{
	double v;
	char *s;
	FILE *df;
	int i,nd,nc;
	int done_err;
	done_err = false;

	pnt_alloc(30);

	df = fopen(fname,"r");
	if (df==NULL) {
		printf("Unable to open data file %s \n",fname);
		exit(1);
	}
	nd = 0;
	for (;!feof(df);) {
	  if (fgets(buff,2000,df)!=NULL) {
		if (strchr(buff,'!')!=NULL) *strchr(buff,'!') = 0;
		nc = 0;
		s = strtok(buff," \t\n,");
		for (;s!=NULL;) {
			v = atof(s);
			pnt_alloc(nd);
			if (isdigit(*s) || *s=='-' || *s=='+' || *s=='.') {
				pntxyz[nd++] = v; nc++;
			} else gprint("Not a number {%s} \n",s);
			s = strtok(NULL," \t\n,");
		}
		if (nc>0 && nc!=3) {
			if (!done_err) gprint("Expecting 3 columns in data file, found %d \n {%s}\n (FATAL ERROR)\n",nc,buff);
			exit(1);
		}
	  }
	}
	fclose(df);
	npnts = nd;
}
double getf()
{
	char buff[80];
	getstr(buff);
	return atof(buff);
}
getstr(char *s)
{
	char *ss;
	ss = gets(s);
	if (ss==NULL) {
		printf("\nERROR, End of input file \n");
		exit(1);
	}
	ss = strchr(s,'\n');
	if (ss!=NULL) *ss = 0;
	ss = strchr(s,'!');
	if (ss!=NULL) *ss = 0;
}
getrange(char *s, double  *min, double  *max, double  *step)
{
	s = strtok(s,", :\n\t");
	if (s!=NULL) *min = atof(s);
	s = strtok(NULL,", :\n\t");
	if (s!=NULL) *max = atof(s);
	*step = (*max - *min) / 15.0;
	s = strtok(NULL,", :\n\t");
	if (s!=NULL) *step = atof(s);
}
setminmax(double v, double *min, double *max)
{
	if (v< *min) *min = v;
	if (v> *max) *max = v;
}
float *xxx,*yyy,*zzz;
myswap(int i, int j)
{
	static float a;
	a = xxx[i]; xxx[i] = xxx[j]; xxx[j] = a;
	a = yyy[i]; yyy[i] = yyy[j]; yyy[j] = a;
	a = zzz[i]; zzz[i] = zzz[j]; zzz[j] = a;
}
mycmp(int i, float x1, float y1)
{
	static float a;
	if (xxx[i]<x1) return -1;
	if (xxx[i]>x1) return 1;
	if (yyy[i]<y1) return -1;
	if (yyy[i]>y1) return 1;
	return 0;
}
int quick_sort(int nitems, int (*fswap) (int i, int j), int (*fcmp) (int i, float x1, float x2));
sort_data(int npnts,float *xd,float *yd,float *zd)
{
	xxx = xd; yyy = yd; zzz = zd;
	quick_sort(npnts,myswap,mycmp);		
}

int (*ffswap) (int i, int j);
int (*ffcmp) (int i, float x1, float x2);
int qquick_sort(int left, int right);
quick_sort(int nitems, int (*fswap) (int i, int j), int (*fcmp) (int i,
	float x1, float x2))
{
	ffswap = fswap;
	ffcmp = fcmp;
	qquick_sort(0,nitems-1);
}
qquick_sort(int left, int right)
{
	int i,j,xx,yy;
	float x1,x2;
	i = left; j = right;
	xx = (left+right)/2;
	x1 = xxx[xx]; x2 = yyy[xx];
	do {
		while((*ffcmp)(i,x1,x2)<0 && i<right) i++;
		while((*ffcmp)(j,x1,x2)>0 && j>left) j--;
		if (i<=j) { (*ffswap)(i,j); i++; j--; }
	} while (i<=j);
	if (left<j) qquick_sort(left,j);
	if (i<right) qquick_sort(i,right);
}