/*
 * Program for creatin histogram data from dump files from BIT1
 * 
 * By J.K. 27.2.2008 
 */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

/*#define DEBUG 0 */
void getMinMaxValues( FILE *fp, int col, float *x1, float *x2, float *v1, float *v2 );
void dispHelp();

/*
float x1 = 100;
float x2 = 700;
#define  NX  40
float yy1 = 150;
float yy2 = 500;
#define NY 30
*/

int main( int argc, char *argv[] ) {
	float x1, x2, v1, v2, xmin, xmax, deltax, dx, stepx;
	float x, v[3], deltav;
	float* stripXpos[2];
	int i, j, nStripX, nStripV, loop, ypos, col, seq;
	int *histoData;
	char file[80], vel[3];
	FILE *fp;
	
	/* init values */
	x1 = x2 = v1 = v2 = 0;
	
	if( argc < 8 ) {
		/* not all arguments are given; display usage */
		dispHelp();
	}
	
	/* read parameters from command line */
	strcpy(file, argv[1]);
	fp = fopen(file, "r");
	
	sscanf( argv[2], "%s", vel );			/* velocity */
	sscanf( argv[3], "%d", &nStripV );		/* number of cells in V[x,z,y] direction */
	sscanf( argv[4], "%f", &xmin );			/* lower limit of area of intrest */
	sscanf( argv[5], "%f", &xmax );			/* upper limit of area of intrest */
	sscanf( argv[6], "%d", &nStripX );		/* number of strips in <dx> direction */
	sscanf( argv[7], "%f", &deltax );		/* width of strip */
	
	if( strncmp( vel, "vx", 2 ) == 0 ) {
		/* we calcualte for vx */
		col = 3;
	} else if( strncmp( vel, "vy", 2 ) == 0 ) {
		/* we calcualte for vy */
		col = 4;
	} else if( strncmp( vel, "vz", 2 ) == 0 ) {
		/* we calcualte for vx */
		col = 5;
	} else {
		/* unknown velocity was given, display warning and use vx */
		printf( "!!== Unknown velocity given, we will use vx\n\n" );
		strcpy( vel, "vx" );
		/* vx is in column 3 */
		col = 3;
	}
	
	/* display input parameters */
	printf("Parameters:\n");
	printf("  file: %s\n", file );
	printf("  velocity: %s (column: %d)\n", vel, col );
	printf("  numX: %d; numVx: %d\n", nStripX, nStripV);
	printf("  x: %f -> %f\n", xmin, xmax );
	
	/* first get min and max in x and vx */
	getMinMaxValues( fp, col, &x1, &x2, &v1, &v2 );
	
	/* display limits for which histogram will be made */
	printf("\nData from file:\n  x: %.3f -> %.3f; vx: %.3f -> %.3f\n", x1, x2, v1, v2 );
	
	stepx = (xmax-xmin) / nStripX;
	
	stripXpos[0] = (float *) malloc( nStripX * sizeof(float) );
	stripXpos[1] = (float *) malloc( nStripX * sizeof(float) );
	
	for(i=0 ; i < nStripX ; i++ ){
		/* calculate begin and end position of every strip */
		stripXpos[0][i] = xmin + (i+0.5)*stepx-deltax/2;
		stripXpos[1][i] = xmin + (i+0.5)*stepx+deltax/2;
	}
	
	/* allocate space for histogram data and fill it with zeros */
	histoData = (int *) malloc( nStripX * nStripV * sizeof(int) );
	memset( histoData, 0, nStripX * nStripV * sizeof(int) );
	
	while (!feof(fp)) {
		/* OK, lets roll.... */
		loop = 1;
		fscanf(fp, "%d%f%f%f%f", &seq, &x, &v[0], &v[1], &v[2]);
	
#ifdef DEBUG
		printf( "x: %f, %f, %f, %f : ", x, v[0], v[1], v[2] );
#endif

		for( i=0 ; i < nStripX && loop ; i++ ){
		    /* find in which strip (if in any) current point falls */
#ifdef DEBUG
			printf( ".%d.", i );
#endif

		    if( x >= stripXpos[0][i] && x <= stripXpos[1][i] ) {
				/* in current strip */
				/* find in which cell it falls */
				ypos = (v[col-3]-v1)/(v2-v1)*(nStripV-1);
				/* add one count to right cell */
				histoData[i*nStripV+ypos] ++;

#ifdef DEBUG
				printf("--> %f < %f < %f -- %d", stripXpos[0][i], x, stripXpos[1][i], ypos );
#endif

				/* stop loop on strips */
				loop = 0;
		    }
		}

#ifdef DEBUG
		printf("\n");
#endif
	}
	
	fclose(fp);
	
	printf( "\nResults:\n");
	deltav = (v2 - v1) / nStripV;
	
	for (j = 0; j < nStripV; j++) {
		/* display the result */
		printf( "%g", v1 + (j+0.5)*deltav  );
	
		for (i = 0; i < nStripX; i++) {
		    printf("%6d ", *(histoData+i*nStripV+j) );
		}
	
		printf("\n");
	}
}


/*
  read values from file and define "bounding" box = min and max values for x and vx
*/
void getMinMaxValues( FILE *fp, int col, float *x1, float *x2, float *v1, float *v2 ){
	float x, v[3];
	int seq;
	
	if( col < 3 && col > 5 ) {
		/* check for input column number; not right, return column 3 (vx) */
		col = 3;
	}
	
	/* go back to begining */
	rewind( fp );
	
	if( !feof(fp) ){
		/* initial values */
		fscanf(fp, "%d%f%f%f%f", &seq, &x, &v[0], &v[1], &v[2]);
		*x1 = *x2 = x;
		*v1 = *v2 = v[col-3];
	}
	
	while( !feof(fp) ){
		/* pase whole file */
		fscanf(fp, "%d%f%f%f%f", &seq, &x, &v[0], &v[1], &v[2]);
		
		if( x < *x1 ) *x1 = x;
		if( x > *x2 ) *x2 = x;
		if( v[col-3] < *v1 ) *v1 = v[col-3];
		if( v[col-3] > *v2 ) *v2 = v[col-3];
	}
	
	/* go back to begining */
	rewind( fp );
}


/*
  display help about using program
*/
void dispHelp(){
	printf( "=== Parameter(s) missing !! ===\n\nParameters: <input file> <col> <ny> <xmin> <xmax> <nx> <dx>\n");
	
	printf( "where:\n");
	printf( "  <input file> = file from which to read input\n");
	printf( "  <col> = which velocity to use (vx, vy, xz; other values give vx)\n");
	printf( "  <ny> = number of cells in Vx direction\n");
	printf( "  <xmin>, <xmax> = limits of area of interest\n");
	printf( "  <nx> = number of strips in area of interes, each is <dx> wide\n");
	printf( "  <dx> = width of strip\n\n");
	
	exit(1);
}
