/*
 * xyzcarter takes an  lon, lat, uncorrected depth file and
 * calculates Carter corrected depth.
 *
 * The methodology used to calculate the corrected depth is that
 * used in mgd77togmt from the original gmt database program.  This
 * method was found to be correct.
 *
 */
#include <time.h>
#include "/usr/local/gmt/src/gmt.h"

#define N_BINS 64800
#define N_CARTER_ZONES 85
#define N_CARTER_OFFSETS 86
#define N_CARTER_CORRECTIONS 5812
#define MAX_VALID_TWT   16000   /* Two-way time (msec) to 12000 nominal m  */

int     carter_setup();
int     get_bin();
int     get_carter_zone();
int     get_depth_from_twt();
int     get_twt_from_depth();

short int       carter_zone[N_BINS];
short int       carter_offset[N_CARTER_OFFSETS];
short int       carter_correction[N_CARTER_CORRECTIONS];

int     carter_not_initialized = 1;
int error_count = 0;
int bad_lon = 0;
int bad_lat = 0;
int bin_out_error = 0;
int zone_out_error = 0;
int bad_twt = 0;
int big_twt = 0;
int bad_depth = 0;
int big_depth = 0;
 
main (argc, argv)
int argc;
char **argv; {

	int lon, lat, bin, zone;
	short int twt;
	short int itcor;
	int n_in, n_out, n_edit; 
        float rlon, rlat, depth0, depthf;
	float veloc = 1500.;

        if (argc < 1) {
                fprintf(stderr,"usage:  xyzcarter  < input_xyz > output_xyz \n\n");
                exit(-1);
        }

	n_in = 0;
	n_out = 0;
	n_edit = 0;

	while ( (scanf("%f %f %f ", &rlon, &rlat, &depth0) != EOF)) {
                n_in++;
		if(rlon < 0.) rlon = rlon + 360.;
		if(rlon > 360.) rlon = rlon - 360.;
		lon = (int)(1000000.*rlon);
		lat = (int)(1000000.*rlat);

/* convert to two-way travel time in milliseconds */
		twt = (short int)(-2000. * depth0/veloc);
/*		fprintf(stderr,"%d %d %d \n",lon, lat, twt);*/

/* Convert the twt to Carter-corrected depth  */
	   if ((get_bin_6(lat, lon, &bin))||(get_carter_zone(bin, &zone))||(get_depth_from_twt(zone, twt, &itcor)) ) {
/*		fprintf(stderr,"xyzcarter:  ERROR in Carter correction system.\n");*/
		n_edit++;
		depthf = depth0;
	   } else {
		depthf = -itcor;
	   }
		if(depthf < -1) {
			n_out++;
			fprintf(stdout,"%f %f %f \n",rlon, rlat, depthf);
		}
	}
	fprintf(stderr,"%s %d %d %d \n","#in #out #no_carter",n_in,n_out,n_edit);
}

/* The following carter depth calculation code is
   substantially the same as from the original mgd77togmt.c code
   by Walter Smith et al.
*/

int	get_bin_6(lat,lon, bin)
int	lat;
int	lon;
int	*bin;
{
	/* Given signed long ints in the 1.0e06 times decimal degree range, 
		-90000000 <= lat < 90000000, 0 <= lon < 360000000, set bin
		number.
		Returns 0 if OK, -1 if error.  */

	int	latdeg, londeg;

	if (lat < -90000000 || lat > 90000000) {
		/*fprintf(stderr, "getbin_6:  Latitude domain error.\n");*/
        error_count++;
        bad_lat++;
		return(-1);
	}
	if (lon < 0 || lon > 360000000) {
		/*fprintf(stderr, "getbin_6:  Longitude domain error.\n");*/
        error_count++;
        bad_lon++;
		return(-1);
	}
	latdeg = (lat + 90000000)/1000000;
	if (latdeg == 180) latdeg = 179;	/* Map north pole to previous row  */

	londeg = lon/1000000;
	*bin = 360 * latdeg + londeg;

	return(0);
}

int	carter_setup()
{
	/* This routine must be called once before using carter table stuff.
	It reads the carter_corr.b, carter_offset.b, carter_zone.b files and
	loads the appropriate arrays.  It sets carter_not_initialized = FALSE
	upon successful completion and returns 0.  If failure occurrs, it
	returns -1.  */

	FILE	*fp = NULL;
	char	buffer [512];
	int	i;

	/* Read the correction table:  */

	sprintf(buffer,"%s/carter_corr.b", "/usr/local/gmt/src");
	if ( (fp = fopen(buffer, "r")) == NULL) {
                fprintf(stderr,"carter_setup:  Cannot open r %s\n", buffer);
		carter_not_initialized = TRUE;
                return(-1);
        }
        for (i = 0; i < N_CARTER_CORRECTIONS; i++) {
                if ((fread((char *)&carter_correction[i], sizeof(short int), 1, fp)) != 1) {
                        fprintf(stderr, "carter_setup:  Cannot read record %ld from %s\n", i, buffer);
		carter_not_initialized = TRUE;
                return(-1);
                }
        }
        fclose(fp);


	/* Read the offset table:  */

	sprintf(buffer,"%s/carter_offset.b", "/usr/local/gmt/src");
	if ( (fp = fopen(buffer, "r")) == NULL) {
                fprintf(stderr,"carter_setup:  Cannot open r %s\n", buffer);
		carter_not_initialized = TRUE;
                return(-1);
        }
        for (i = 0; i < N_CARTER_OFFSETS; i++) {
                if ((fread((char *)&carter_offset[i], sizeof(short int), 1, fp)) != 1) {
                        fprintf(stderr, "carter_setup:  Cannot read record %ld from %s\n", i, buffer);
		carter_not_initialized = TRUE;
                return(-1);
                }
        }
        fclose(fp);


	/* Read the zone table:  */

	sprintf(buffer,"%s/carter_zone.b", "/usr/local/gmt/src");
	if ( (fp = fopen(buffer, "r")) == NULL) {
                fprintf(stderr,"carter_setup:  Cannot open r %s\n", buffer);
		carter_not_initialized = TRUE;
                return(-1);
        }
        for (i = 0; i < N_BINS; i++) {
                if ((fread((char *)&carter_zone[i], sizeof(short int), 1, fp)) != 1) {
                        fprintf(stderr, "carter_setup:  Cannot read record %ld from %s\n", i, buffer);
		carter_not_initialized = TRUE;
                return(-1);
                }
        }
        fclose(fp);

	/* Get here when all is well.  */

	carter_not_initialized = FALSE;
	return(0);
}


int	get_carter_zone(bin, zone)
int	bin;
int	*zone;
{
	/* Sets value pointed to by zone to the Carter zone corresponding to
		the bin "bin".  Returns 0 if successful, -1 if bin out of
		range.  */

	if (carter_not_initialized && carter_setup() ) {
		fprintf(stderr,"get_carter_zone:  Initialization failure.\n");
		return(-1);
	}

	if (bin < 0 || bin >= N_BINS) {
		/*fprintf(stderr,"get_carter_zone:  bin out of range.\n");*/
        error_count++;
        bin_out_error++;
		return(-1);
	}
	*zone = carter_zone[bin];
	return(0);
}


int	get_depth_from_twt(zone, twt_in_msec, depth_in_corr_m)
int	zone;
short int	twt_in_msec;
short int *	depth_in_corr_m;
{
	/* Given two-way travel time of echosounder in milliseconds, and
		Carter Zone number, finds depth in Carter corrected meters.
		Returns (0) if OK, -1 if error condition.  */

	int	i, nominal_z1500, low_hundred, part_in_100;

	if (carter_not_initialized && carter_setup() ) {
		fprintf(stderr,"get_depth_from_twt:  Initialization failure.\n");
		return(-1);
	}
	if (zone < 1 || zone > N_CARTER_ZONES) {
		/*fprintf(stderr,"get_depth_from_twt:  Zone out of range.\n");*/
        error_count++;
        zone_out_error++;
		return(-1);
	}
	if (twt_in_msec < 0) {
		/*fprintf(stderr,"get_depth_from_twt:  Negative twt.\n");*/
        error_count++;
        bad_twt++;
		return(-1);
	}

	nominal_z1500 = irint(0.75 * twt_in_msec);

	if (nominal_z1500 <= 100) {
		/* There is no correction in water this shallow.  */
		*depth_in_corr_m = nominal_z1500;
		return(0);
	}

	low_hundred = nominal_z1500 / 100;
	i = carter_offset[zone - 1] + low_hundred - 1;	/* -1 'cause .f indices */
	
	if (i >= (carter_offset[zone] - 1) ) {
		/*fprintf(stderr, "get_depth_from_twt:  twt too big.\n");*/
        error_count++;
        big_twt++;
		return(-1);
	}

	part_in_100 = nominal_z1500%100;

	if (part_in_100) {	/* We have to interpolate the table  */

		if ( i == (carter_offset[zone] - 2) ) {
			/*fprintf(stderr, "get_depth_from_twt:  twt too big.\n");*/
            error_count++;
            big_twt++;
			return(-1);
		}

		*depth_in_corr_m = (short) irint (carter_correction[i] +
			0.01 * part_in_100 * (carter_correction[i+1]
				- carter_correction[i]) );
		return(0);
	}
	else {
		*depth_in_corr_m = carter_correction[i];
		return(0);
	}
}


int	get_twt_from_depth(zone, depth_in_corr_m, twt_in_msec)
int	zone;
short int	depth_in_corr_m;
short int *	twt_in_msec;
{
	/*  Given Carter zone and depth in Carter corrected meters,
	finds the two-way travel time of the echosounder in milliseconds.
	Returns -1 upon error, 0 upon success.  */

	int	min, max, guess;
	double	fraction;

	if (carter_not_initialized && carter_setup() ) {
		fprintf(stderr,"get_twt_from_depth:  Initialization failure.\n");
		return(-1);
	}
	if (zone < 1 || zone > N_CARTER_ZONES) {
		/*fprintf(stderr,"get_twt_from_depth:  Zone out of range.\n");*/
        error_count++;
        zone_out_error++;
		return(-1);
	}
	if (depth_in_corr_m < 0) {
		/*fprintf(stderr,"get_twt_from_depth:  Negative depth.\n");*/
        error_count++;
        bad_depth++;
		return(-1);
	}

	if (depth_in_corr_m <= 100) {	/* No correction applies.  */
		*twt_in_msec = irint(1.33333 * depth_in_corr_m);
		return(0);
	}

	max = carter_offset[zone] - 2;
	min = carter_offset[zone - 1] - 1;

	if (depth_in_corr_m > carter_correction[max]) {
		/*fprintf(stderr,"get_twt_from_depth:  Depth too big.\n");*/
        error_count++;
        big_depth++;
		return(-1);
	}

	if (depth_in_corr_m == carter_correction[max]) {
		/* Hit last entry in table exactly  */
		*twt_in_msec = irint(133.333 * (max - min) );
		return(0);
	}


	guess = (depth_in_corr_m / 100) + min;
	if (guess > max) guess == max;
	while (guess < max && carter_correction[guess] < depth_in_corr_m) guess++;
	while (guess > min && carter_correction[guess] > depth_in_corr_m) guess--;

	if (depth_in_corr_m == carter_correction[guess]) {
		/* Hit a table value exactly  */
		*twt_in_msec = irint(133.333 * (guess - min) );
		return(0);
	}
	fraction = ( (double)(depth_in_corr_m - carter_correction[guess])
			/ (double)(carter_correction[guess + 1] - carter_correction[guess]) );
	*twt_in_msec = irint(133.333 * (guess - min + fraction) );
	return(0);
}