#include	<math.h>
#include	"misc.h"
#define		BINS	102


/* BINSTRM -- binned statistics of input stream */

binstrm (input, hflag, aflag, cflag, tflag, sflag, rflag, m,
    xmin, xmax, xbins, ymin, ymax, ybins)

char	input[];
int	*hflag, *aflag, *cflag, *tflag, *sflag, *rflag, *m;
float	*xmin, *xmax, *ymin, *ymax;
int	*xbins, *ybins;

{
    FILE    *fd;
    float   temp;
    float   xbin[BINS], ybin[BINS];
    float   count[BINS][BINS], sum[BINS][BINS], sum2[BINS][BINS];
    float   zmin[BINS][BINS], zmax[BINS][BINS];
    int     nxbins, nybins, i, j, k, l;
    float   x, y, z, dx, dy;

    if (strcmp (input, "STDIN") == 0)
	fd = getstdin;
    else if ((fd = fopen (input, "r")) == NULL) {
	printf ("BINSTRM: Can't read file %s\n", input);
	return;
    }

    nxbins = *xbins;
    nybins = *ybins;
    if (nxbins >= BINS - 1)
	nxbins = BINS - 2;
    if (nybins >= BINS - 1)
	nybins = BINS - 2;
    if (*m > 0) {
	dx = (*xmax - *xmin) / nxbins;
	for (i = 0; i <= nxbins+1; i++)
	    xbin[i] = (i - 1) * dx + *xmin;
    }
    if (*m > 1) {
	dy = (*ymax - *ymin) / nybins;
	for (i = 0; i <= nybins+1; i++)
	    ybin[i] = (i - 1) * dy + *ymin;
    }

    for (i = 0; i <= nxbins; i++) {
	for (j = 0; j <= nybins; j++) {
	    count[i][j] = 0.;
	    sum[i][j] = 0.;
	    sum2[i][j] = 0.;
	    zmin[i][j] = 1e30;
	    zmax[i][j] = -1e30;
	}
    }

    for (i = 1, j = 1;;) {
	switch (*m) {
	    case 1: 
		scanf ("%hf", &x);
		if (feof (fd))
		    break;
		if ((x < xbin[0]) || (x >= (xbin[nxbins+1])))
		    x = xbin[0];
		for (i = 0; (x >= xbin[i + 1]) && (i < nxbins); i++);
		break;
	    case 2: 
		scanf ("%hf", &x);
		if (feof (fd))
		    break;
		if ((x < xbin[0]) || (x >= (xbin[nxbins+1])))
		    x = xbin[0];
		for (i = 0; (x >= xbin[i + 1]) && (i < nxbins); i++);
		scanf ("%hf", &y);
		if (feof (fd))
		    break;
		if ((y < ybin[0]) || (y >= (ybin[nybins+1])))
		    y = ybin[0];
		for (j = 0; (y >= ybin[j + 1]) && (j < nybins); j++);
		break;
	}
	scanf ("%hf", &z);
	if (feof (fd))
	    break;
	count[i][j] += 1;
	sum[i][j] += z;
	sum2[i][j] += z * z;
	if (z < zmin[i][j])
	    zmin[i][j] = z;
	if (z > zmax[i][j])
	    zmax[i][j] = z;
    }

    fclose (fd);

    /* print header */
    if (*hflag) {
	if (*m > 0)
	    printf ("           X          ");
	if (*m > 1)
	    printf ("           Y          ");
	if (*cflag)
	    printf ("      Count");
	if (*tflag)
	    printf ("        Sum");
	if (*aflag)
	    printf ("    Average");
	if (*sflag)
	    printf ("    Std Dev");
	if (*rflag)
	    printf ("    Minimum");
	if (*rflag)
	    printf ("    Maximum");
	printf ("\n");
    }

    /* print statistics */
    for (i = 1; i <= nxbins; i++) {
	for (j = 1; j <= nybins; j++) {
	    if ((count[i][j] == 0) && (*cflag == 0))
		continue;
	    if (*m > 0)
		printf (" %10g %10g", xbin[i], xbin[i+1]);
	    if (*m > 1)
		printf (" %10g %10g", ybin[j], ybin[j+1]);
	    if (*cflag)
		printf (" %10.0f", count[i][j]);
	    if (*tflag)
		printf (" %10g", sum[i][j]);
	    if (*aflag)
		if (count[i][j] > 0)
		    printf (" %10g", sum[i][j] / count[i][j]);
		else
		    printf (" %10g", 0.);
	    if (*sflag)
		if (count[i][j] > 1) {
		    temp = (count[i][j] * sum2[i][j] - sum[i][j] * sum[i][j]) /
			(count[i][j] * (count[i][j] - 1));
		    printf (" %10g", temp>0.?sqrt(temp):0.);
		} else
		    printf ("          0");
	    if (*rflag)
		printf (" %10g %10g", zmin[i][j], zmax[i][j]);
	    printf ("\n");
	}
    }
}