/* imgbm_etopo_diff.c <etopo5.img> [<threshhold>] < binary input > ascii output
 *
 * Read 8-byte binary records consisting of two 4-byte ints;
 * the first is the index to the img pixel, and the second is the data value.
 * Look up difference between these values and values in etopo5.img, and write
 * difference to stdout. 
 * If threshhold given, only write difference if greater than threshhold.
 * Otherwise, use a default threshhold of 350m, which should find travel time
 * errors of 1/2 second.
 * 
 * Later could add option to dump percentage differences, in order to find
 * problems in shallow water.
 *
 * For now, I have not done any clever buffering of the array, so you should
 * run it on seasat.
 *
 * W H F Smith,  16 July 1996.
 *
 * 25 July 1996 REVISED by WHFS to dump value when either:  1. difference
 * in depth is greater than threshold, or 2.  percentage difference is
 * greater than 10%.  Also, it now dumps xy coordinates of Mercator map,
 * instead of lat lon.  Plot with -R0/360/0/211.2 -Jx
 *
 *
 */

#include "../../../include/gmt.h"

#define NX 10800
#define NY 6336

struct IMGPAIR {
	int	index;
	int	depth;
};

main(argc, argv)
int	argc;
char	**argv;
{
	
	int	i, j, diff, thresh;
	struct IMGPAIR imgpair;	
	FILE	*fp;
	short	*etopo5;
	double	radius, y, radlat, lon, lat, mean, pct;
	int	bad_pct;

	if (argc < 2 || (fp = fopen(argv[1], "r")) == NULL) {
		fprintf(stderr,"usage:  imgbm_etopo_diff <etopo5.img> [<threshhold>] < <binary input> > <ascii output>\n\n");
		exit(-1);
	}
	if (argc == 3 && (sscanf(argv[2], "%ld", &thresh)) == 1) {
		thresh = abs(thresh);
	}
	else {
		thresh = 350;
	}
	radius = (NX/2)/M_PI;
	

	etopo5 = (short *)memory(CNULL, NX*NY, 2, "imgbm_etopo_diff");

	if ( (fread((char *)etopo5, 2, NX*NY, fp)) != NX*NY) {
		fprintf(stderr,"imgbm_etopo_diff:  etopo5 read failure.\n");
		exit(-1);
	}
	fclose(fp);

	while ( (fread((char *)&imgpair, sizeof(struct IMGPAIR), 1, stdin)) == 1) {
		diff = imgpair.depth - etopo5[imgpair.index];
		mean = 0.5 * (imgpair.depth + etopo5[imgpair.index]);
		pct = (mean == 0) ? 100.0 : 100.0 * diff/mean;
		bad_pct =  (fabs(pct) > 20.0 || (fabs(pct) > 10.0 && fabs(mean) > 200.0) );
		if (abs(diff) > thresh || bad_pct) {
			i = imgpair.index % NX;
			j = imgpair.index / NX;
		  /*	y = ((NY/2) - j - 0.5)/radius;
			radlat = 2.0 * atan(exp(y))-M_PI_2;
			lon = (i + 0.5)/30.0;
			lat = 180.0 * radlat / M_PI;   */
			y = (NY - (j + 0.5)) / 30.0;
			printf("%.8lg\t%.8lg\t%ld\t%.8lg\n", (i+0.5)/30.0, y, diff, pct);
		}
	}
	free((char *)etopo5);
	exit(0);
}