/*
	Contouring program, creates gle files.
*/

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define false (0)
#define true (!false)

char *clab[100];
float cval[100];
char outlab[80];
char outdat[80];
char outgle[80];
char outkey[80];
char outfile[80];
int smooth=false,smoothsub;
FILE *fp,*flab,*fdat;

#ifdef unix 
#define VAXC 1
#endif
#ifdef VAXC
char *strdup(char *s);
#else
#include <alloc.h>
#endif
int gen_next(char *pat, double *v);
int gen_cval(char *cstr, float cval[], int *ncv) ;
FILE *myfopen(char *fname, char *mode);
int read_data(char *fname, float *zmin,float *zmax,int *nx,int *ny);
int docontour(float zz[], long nrz, long nx, long ny
	, float cval[], char *clab[], long ncv, float zmax);

int addvect(int m, double x, double y);
FILE *myfopen(char *fname, char *mode)
{
	FILE *f;
	f = fopen(fname,mode);
	if (f==NULL) {
		printf("Unable to open {%s} \n",fname); perror(""); exit(1);
	}
	return f;
}
int dflt(char *s, char *d);
double dfltv( char *d);
double dfltv( char *d)
{
	char buff[80];
	dflt(buff,d);
	return atof(buff);
}
dflt(char *s, char *d)
{
	printf(" [%s] ? ",d); gets(s);
	if (strcmp(s,"")==0) strcpy(s,d);
}
#ifdef VAXC
#define farmalloc malloc
#endif
static double dxmin,dymin,dxmax,dymax;
static float *zz;
static 	double nnx,nny;
main()
{
	float zmin = 1e10, zmax = -1e10;
	char infile[80],cdflt[80],cstr[80],lablet[80];
	int nx,ny,i,ncv,dokey;
	int zdim=51;
	char buff[90];
	
	printf("Data file containing Z values");  dflt(infile,"");
	read_data(infile,&zmin,&zmax,&nx,&ny);
	zdim = nx;
	if (dxmin==dxmax) {dxmax = nx; dxmin = 1;}
	if (dymin==dymax) {dymax = ny; dymin = 1;}

	sprintf(cdflt,"%g:%g:%g",zmin,zmax,(zmax-zmin)/10.0);
	printf("Specify values to contour at, e.g. 1 2 3 4 or 1:4:.5\n");
	dflt(cstr,cdflt);
	printf("Label contours with letters A...Z");  dflt(lablet,"YES");
	if (toupper(lablet[0]) == 'Y') dokey = true; else dokey = false;
	gen_cval(cstr, cval, &ncv);
	for (i=0; i<ncv; i++) {
	  if (dokey) {
		sprintf(buff,"%c",i+'A');
		clab[i] = strdup(buff);
		printf("{%s} ",clab[i]);		
	  } else {
		sprintf(buff,"%g",cval[i]);
		clab[i] = strdup(buff);
	  }
	}	 
	if (toupper(lablet[0]) == 'Y') printf("\n");
	printf("Contour values ");
	for (i=0; i<ncv; i++)	printf("%g ",cval[i]);
	printf("\n");
	printf("Warning: Smoothing is SLOW\n");
	printf("Smoothing factor (points per grid step 3-10, 0=no smoothing) "); dflt(lablet,"0");
	smoothsub = atoi(lablet);
	if (smoothsub>0) smooth = true;

	printf("Output file names (E.G. xxx.gle, xxx.dat, xxx.lab, xxx.key");
	dflt(outfile,"xxx");
	
	strcpy(outgle,outfile); strcat(outgle,".gle");
	strcpy(outdat,outfile); strcat(outdat,".dat");
	strcpy(outlab,outfile); strcat(outlab,".lab");
	strcpy(outkey,outfile); strcat(outkey,".key");

	fp = myfopen(outgle,"w");
	fprintf(fp,"size 24 18\n");
	fprintf(fp,"begin graph\n");
	fprintf(fp,"	size 24 18\n");
	fprintf(fp,"	xaxis min %g max %g  \n",dxmin,dxmax);
	fprintf(fp,"	yaxis min %g max %g  \n",dymin,dymax);
	fprintf(fp,"	d1 lstyle 1 \n",nx,ny);
	fprintf(fp,"	d1 bigfile %s \n",outdat);
	fprintf(fp,"end graph\n");
	fprintf(fp,"save theorigin\n");
	fprintf(fp,"set hei .3 just cc\n");
	fprintf(fp,"bigfile %s\n",outlab);
	fprintf(fp,"move theorigin\n");
	fprintf(fp,"set hei .3 just left \n");
	if (dokey) fprintf(fp,"bigfile %s\n",outkey);
	fclose(fp);

	flab = myfopen(outlab,"w");
	fdat = myfopen(outdat,"w");
	fprintf(flab,"set lwidth .0000001\n");
	nnx = nx; nny = ny;

	docontour(zz,zdim,nx,ny,cval,clab,ncv,zmax);

	fclose(flab);
	fclose(fdat);

	fp = myfopen(outkey,"w");
	fprintf(fp,"rmove .2 .2 \n");
	fprintf(fp,"set lwidth 0\n");
	fprintf(fp,"save keystart \n");
	fprintf(fp,"begin box fill white add .2 \n");
	for (i=ncv-1; i>=0; i--) {
		fprintf(fp,"text %s\n",clab[i]);
		fprintf(fp,"rmove .6 0 \n");
		fprintf(fp,"text %g\n",cval[i]);
		fprintf(fp,"rmove -.6 .5\n");
	}
	fprintf(fp,"end box \n");
	fprintf(fp,"move keystart \n");
	for (i=ncv-1; i>=0; i--) {
		fprintf(fp,"text %s\n",clab[i]);
		fprintf(fp,"rmove .6 0 \n");
		fprintf(fp,"text %g\n",cval[i]);
		fprintf(fp,"rmove -.6 .5\n");
	}
	fprintf(fp,"move keystart \n");
	fclose(fp);
}
docontour(float zz[], long nrz, long nx, long ny
	, float cval[], char *clab[], long ncv, float zmax)
{
	extern int draw_(float *x, float *y, long *iflag), plot_(float *, float *, int *);
	static long *work, i, j;
	extern  gcontr_(float *, long *, long *, long *, float *, long *, float *, long *, 
	   int draw_(float *x, float *y, long *iflag) );
/*     dimension of work is large enough to contain 	*/
/*     2*(dimension of c)*(total dimension of z) useful bits.  */
	work = farmalloc(( (long) sizeof(long) * 2 * ncv * nx * ny)/31 + 10);
	if (work==NULL) {
		printf("Unable to allocate storage for work array\n");
		exit(1);
	}
	zmax = zmax + 100; /* no clipping */
    	gcontr_(zz, &nrz, &nx, &ny, cval, &ncv, &zmax, work, draw_);
}

/*#define sx(v) (dxmin + (dxmax-dxmin)*(v)/nnx) */
#define sx(v) (dxmin + (dxmax-dxmin)*(v-1)/(nnx-1))
#define sy(v) (dymin + (dymax-dymin)*(v-1)/(nny-1))
int draw_(float *x, float *y, long *iflag)
{
	static float xcur,ycur;
	static int ih, il;

	ih = *iflag / 10;
	il = *iflag - ih * 10;
	switch (il) {
	  case 6:   /* Get current point */
		*x = xcur;
	    	*y = ycur;
	    	break;
	  case 2: /* start contour at bondary */
	  case 3: /* start contour not at bondary */
		if (smooth) addvect(1,sx(*x),sy(*y));
		else {
			fprintf(fdat,"* * \n",*x,*y);
			fprintf(fdat,"%g %g \n",sx(*x),sy(*y));
		}
		fprintf(flab,"amove xg(%g) yg(%g) \n",sx(*x),sy(*y));
		fprintf(flab,"begin box fill white add .07 \n");
		fprintf(flab," text %s \n",clab[ih-1]);
		fprintf(flab,"end box \n");
		fprintf(flab,"text %s \n",clab[ih-1]);
		break;
	  case 1: /* continue a contour */
		if (smooth) addvect(2,sx(*x),sy(*y));
		else fprintf(fdat,"%g %g \n",sx(*x),sy(*y));
		break;
	  case 4: /* finish contour at a boundary */
		if (smooth) addvect(4,sx(*x),sy(*y));
		else fprintf(fdat,"%g %g \n",sx(*x),sy(*y));
		break;
	  case 5: /* finish close contour */
		if (smooth) addvect(3,sx(*x),sy(*y));
		else fprintf(fdat,"%g %g \n",sx(*x),sy(*y));
		break;
	}
	xcur = *x;
	ycur = *y;
} 
gen_cval(char *cstr, float cval[], int *ncv)
{
	double v;
	int j;
	*ncv = 0;
	gen_next(cstr,&v);
	for (;gen_next(NULL,&v);) {
		cval[*ncv] = v;
		(*ncv)++;
	}
}
gen_next(char *pat, double *v)
{
	static char *s,p[200];
	static double v1,v2,v3;
	static int instep;
	char *c1,*c2;
	if (pat!=NULL) {
		strcpy(p,pat);
		s = strtok(p,", ");
		return true;
	}
	if (instep) {
	       v1 += v3;
	       *v = v1;
	       if (v1<=(v2+v3/10)) return true;
	       else instep = false;
	}
	if (s==NULL) return false;
	c2 = c1 = strchr(s,':');
	if (c1!=NULL) c2 = strchr(c1+1,':');

	if (c1==NULL) {
		*v = atof(s);
		s = strtok(NULL,", ");
		return true;
	}
	v1 = atof(s);
	v2 = atof(c1+1);
	if (c2==NULL) v3 = 1; else v3 = atof(c2+1);
	*v = v1;
	instep = true;
	s = strtok(NULL,", ");
	if (v1>v2) return false;
	return true;
}
static char buff[2032];
FILE *myfopen(char *fname, char *mode);
int alloc_zdata(int nx, int ny);
alloc_zdata(int nx, int ny)
{
	if (zz!=NULL) free(zz);
	zz = farmalloc((long) nx * ((long) ny+1) * sizeof(float));
	if (zz==NULL) {
		printf("Unable to allocate enough memory for datafile\n");
		return true;
	}
	return false;
}
double getkeyval(char *buff,char *k);
double getkeyval(char *buff,char *k)
{
	char *s;
	s = strstr(buff,k);
	if (s!=NULL) return atof(s+strlen(k));
	return 0.0;
}
FILE *df;
read_data(char *fname, float *zmin,float *zmax,int *nx,int *ny)
{	
	double v;
	long x,y;
	int c,b;
	char *s;
	int i;
	x = y = 0;
	*nx = 0; *ny = 0;

	df = myfopen(fname,"r");
	if (df==NULL) {*nx = 0; *ny = 0; return;}
	for (;!feof(df);) {
	  if (fgets(buff,2000,df)!=NULL) {
		if (*nx==0) {
			strupr(buff);
			*nx  = getkeyval(buff,"NX");
			*ny  = getkeyval(buff,"NY");
			dxmin = getkeyval(buff,"XMIN");
			dymin = getkeyval(buff,"YMIN");
			dxmax = getkeyval(buff,"XMAX");
			dymax = getkeyval(buff,"YMAX");
			if (*nx==0 || *ny==0) {
				printf("Expecting to find  ! nx 11 ny 11 xmin 0 xmax 10 ymin 0 ymax 10 \n");
				printf("on first line of data file,  Please enter this data \n");
				printf("Number of columns (nx) ");  *nx = dfltv("");
				printf("Number of rows (ny) ");  *ny = dfltv("100");
				printf("Min x ");  dxmin = dfltv("0");
				printf("Max x ");  dxmax = dfltv("1");
				printf("Min y ");  dymin = dfltv("0");
				printf("Max y ");  dymax = dfltv("1");
			} else {
				fgets(buff,2000,df);
			}
			if (alloc_zdata(*nx,*ny)) return;
		}
check_again:
		b = strlen(buff);
		c = buff[b-1];
		if (strchr(" \n\t",c)==NULL) { /* we're halfway thru a number */
			buff[b] = getc(df);
			buff[b+1] = 0;
			goto check_again;
		}
		s = strchr(buff,'!');
		if (s!=NULL) *s = 0;
		s = strtok(buff," \t\n,");
		for (;s!=NULL;) {
			v = atof(s);
			if (isdigit(*s) || *s=='-' || *s=='+' || *s=='.') {
				if (x>= *nx) {
					x = 0; y++;
				}
				if (y>= *ny) {
					printf("Too much data in data file %ld %d \n",y,*ny);
					return;
				}
				if (v < *zmin) *zmin = v;
				if (v > *zmax) *zmax = v;
				zz[x + y * *nx] = v;		
				x++;
			} else printf("Not a number {%s} \n",s);
			s = strtok(NULL," \t\n,");
		}	
	  }
	}
	fclose(df);
	*ny = y+1;
	return;
abort_data:
	fclose(df);
	*nx = 0; *ny = 0;
	return;

}

#ifndef __TURBOC__
char *strdup(char *s)
{
	char *ss;
	ss = malloc(strlen(s)+1);
	strcpy(ss,s);
	return ss;
}
strupr(char *s)
{
	for (;*s!=0;s++) *s = toupper(*s);
}
#endif
plot_(float *a, float *b, int *c){}


int nsm,nout_alloc;
float xin[400],yin[400],*xout,*yout;
addvect(int m, double x, double y)
{
	long mode,ninp,nsub,nout;
	int i;
	int snipends=false;
	if (m==1) nsm = 0;
	if (nsm>0) if (xin[nsm-1]==x && yin[nsm-1]==y) {
		return;
	} 
	xin[nsm] = x;
	yin[nsm++] = y;
	if (nsm>299) {
		printf("Too many points to smooth, don't use smooth option\n");
		exit(1);
	}
	/* Connect the ends by adding some points if contour is right sort */
	/* 4=boundary, 3=closed shape */
	if (m==3 || m==4) { 
		if (nsm<2) {
			for (i=0;i<nsm;i++) {
				fprintf(fdat,"%g %g \n",xin[i],yin[i]);
			}
			nsm = 0;
			return;
		}
		if (m==3) {
			snipends = true;
			for (i=nsm; i>0; i--) { 
				xin[i] = xin[i-1];
				yin[i] = yin[i-1]; 
			}
			xin[0] = xin[nsm-1];
			yin[0] = yin[nsm-1];
			xin[nsm+1] = xin[2];
			yin[nsm+1] = yin[2];
			nsm += 2;
			/*printf("\n");
			for (i=0; i<nsm; i++) { 
				printf("%g %g \n",xin[i],yin[i]);
			}*/
		}

		ninp = nsm;		/* Number of input points */
		mode = 2; 		/* multi valued */
		nsub = smoothsub;	/* user parameter */
		if (nsub<3) nsub = 3;
		nout = (ninp-1)*nsub+1;
		if (nout>nout_alloc) {
			xout = malloc(nout*sizeof(*xout));
			yout = malloc(nout*sizeof(*xout));
			nout_alloc = nout;
		}
		glefitcf_(&mode,xin,yin,&ninp,&nsub,xout,yout,&nout);
		nsm = 0;

		fprintf(fdat,"* * \n");
		if (snipends) {
		 for (i=nsub;i<nout-nsub;i++) {
			fprintf(fdat,"%g %g \n",xout[i],yout[i]);
		 }
		} else {
		 for (i=0;i<nout;i++) {
			fprintf(fdat,"%g %g \n",xout[i],yout[i]);
		 }
		}
	}
}