/* img2min_gradstat.c
 *
 * read a 2 minute img file and compute as if averaging into
 * 4 minute sampling, but compute the gradient in the N to S
 * direction for illumination.
 * Does not include a test for 32767, because gradient should
 * be computed on world_grav.img.7.2 before landmasking.
 * Assumes file has tracks encoded.
 *
 * LATER MODIFIED to comment out track encoding fix.  9 July 96
 *
 * After one run to check the min (-710) and max (743) I set
 * it up to make a CDF in 5 unit (1/2 mgal/pixel) steps from
 * -750 to 750.
 *
 * LATER MODIFIED to min (-983) and max (908) of downward-
 * continued, bandpassed version.  9 July 96
 *
 *
 * 30 January 1996
 */

#include "/usr/local/gmt/src/gmt.h"

#define NXIN 10800		/* 360 x 30  */
#define NTOTIN 21600		/* NXIN * 2  */
#define NYOUT 3168
#define NXOUT 5400		/* 360 x 15  */
#define NBOXES 400
#define BOXSIZE 5
#define BOXSTART 1000

short int data[NTOTIN];
int cum[NBOXES];

main(argc, argv)
int     argc;
char    **argv;
{
	int	i,j, ii, k, east, south, deriv;
	short	a, b, c, d, mingrad, maxgrad;
	FILE	*fp;

	if (argc != 2) {
		fprintf(stderr,"usage:  img2min_gradstat.c <imgfile>\n");
		exit(-1);
	}

	if ( (fp = fopen(argv[1], "r") ) == NULL) {
		fprintf(stderr,"ERROR opening file %s.\n", argv[1]);
		exit(-1);
	}

	mingrad = maxgrad = 0;
	for (k = 0; k < NBOXES; k++) cum[k] = 0;

	for (j = 0, k = 0; j < NYOUT; j++) {
		if ( (fread((char *)data, 2, NTOTIN, fp)) != NTOTIN) {
			fprintf(stderr,"ERROR reading input.\n");
			exit(-1);
		}
		/* Check data and demux the tracks:  */
		for (i = 0; i < NTOTIN; i++) {
			if (data[i] == 32767) {
				fprintf(stderr,"ERROR values have 32767.\n");
				exit(-1);
			}
/*			if (data[i]%2) data[i]--;     */
		}
		for (i = 0; i < NXOUT; i++, k++) {
			ii = 2*i;
			a = data[ii + 1];
			b = data[ii];
			c = data[ii+NXIN];
			d = data[ii+NXIN+1];
			east  = ((a+d)-(b+c))/2;
			south = ((c+d)-(a+b))/2;
			if (abs(east)>abs(south)) {
				deriv = east;
			}
			else {
				deriv = south;
			}
			if (mingrad > deriv) mingrad = deriv;
			if (maxgrad < deriv) maxgrad = deriv;
			ii = ceil((double)(deriv + BOXSTART) / (double)BOXSIZE);
			for (k = ii; k < NBOXES; k++) cum[k]++;
		}
	}
	fclose (fp);

	fprintf(stderr,"min and max grad:  %hd %hd\n", mingrad, maxgrad);

	for (k = 0; k < NBOXES; k++) printf("%ld\t%.8lg\n", k*BOXSIZE - BOXSTART, (double)cum[k]/(double)(NXOUT*NYOUT));

	exit (0);
}