/* -------------------------------------------------------------------
 *	$Id: mgd77sniffer.c,v 1.55 2006/05/08 07:49:08 pwessel Exp $
 *
 *	File:	mgd77sniffer.c
 *
 *	mgd77sniffer scans MGD77 files for errors in three ways: one, point-
 *	by-point scanning to determine if values are within reasonable limits;
 *	two, along-track scanning to find excessive changes in values; three,
 *	comparing ship gathered gravity and topography data with global
 *	reference grids for errors. Errors are reported to standard output
 *	by default with optional output of structured "E77" errata tables.
 *
 *	Authors:
 *		Michael Chandler and Paul Wessel
 *		School of Ocean and Earth Science and Technology
 *		University of Hawaii
 * 
 *	Date:	March 2006
 * 
 * ------------------------------------------------------------------*/

#include "mgd77.h"
#include "mgd77sniffer.h"

/*
#define HISTOGRAM_MODE 0
#define FIX 0
#define OUTPUT_TEST 0
*/

#if OUTPUT_TEST == 1
#define OR_TRUE || 1
#define AND_FALSE && 0
#else
#define OR_TRUE
#define AND_FALSE
#endif

int main (int argc, char **argv) {

	/* THE FOLLOWING VARIABLES DO NOT VARY FOR EACH CRUISE */
	int argno, n_alloc = GMT_CHUNK, n_cruises = 0, pos = 0, n_grids = 0, n_out_columns;
	unsigned int MGD77_this_bit[32], n_types[N_ERROR_CLASSES], n_bad_sections = 0;

	double time_factor = 1.0, distance_factor = 1.0, maxTime, west=0.0, east=0.0, north=0.0, south=0.0, adjustDC[32]; 
	double test_slope[5] = {0.1, 10.0, 0.5337, MGD77_METERS_PER_FATHOM, MGD77_FATHOMS_PER_METER}, adjustScale[32];
	double max_speed, MGD77_NaN, maxSlope[MGD77_N_NUMBER_FIELDS], maxGap, threshold = 1.0;
	float *f[8];
	time_t clock;

	char c, tmp_min[16], tmp_max[16], tmp_maxSlope[16], tmp_area[16], *derivative;
	char *custom_limit_file = NULL, custom_limit_line[BUFSIZ], arguments[BUFSIZ];
	char field_abbrev[8], *speed_units = "m/s", *distance_units = "km";
	char file[BUFSIZ], *display = NULL;

	BOOLEAN error = FALSE, nautical = FALSE, custom_max_speed = FALSE, simulate = FALSE, bilinear = FALSE, bad_sections = FALSE;
	BOOLEAN custom_warn = FALSE, warn[MGD77_N_WARN_TYPES], custom_maxGap = FALSE, decimate = TRUE, adjustData = FALSE;

	FILE *custom_fp, *fpout = NULL;

	struct MGD77_SNIFFER_DEFAULTS mgd77snifferdefs[MGD77_N_DATA_FIELDS] = {
#include "mgd77snifferdefaults.h"
	};

	struct BAD_SECTION BadSection[MAX_BAD_SECTIONS];
	
	/* THESE VARIABLES VARY FOR EACH CRUISE AND REQUIRE EXTRA CARE (RESET FOR EACH CRUISE) */
	int i, j, k, curr, nwords, nout, nvalues, distanceErrorCount, duplicates[MGD77_N_NUMBER_FIELDS], *iMaxDiff = NULL, n_nan;
	int noTimeCount, noTimeStart, timeErrorCount, timeErrorStart, distanceErrorStart, overLandStart, overLandCount, last_day;
	int **bin2d, ship_bin, grid_bin, n, npts, *offsetStart, rec, type, field, bccCode, col;
	unsigned int lowPrecision, lowPrecision5;

	double gradient, dvalue, dt, ds, **out, thisArea, speed, **G = NULL, min, *distance, date;
	double *offsetArea, stat[8], stat2[8], *ship_val, *grid_val, max, tcrit, se, range, range2, n_days;
	double thisLon, thisLat, lastLon, lastLat, *MaxDiff = NULL, **diff = NULL, *decimated_ship;
	double *offsetLength, *decimated_grid, S_xx, recommended_scale, *new_anom, *old_anom, IGRF[8];

	char timeStr[32], placeStr[64], errorStr[128], outfile[32], abbrev[8];

	BOOLEAN gotTime, landcruise, *offsetSign, newScale, mtf1;
	BOOLEAN *prevOffsetSign, prevFlag, prevType, decimated;
#ifdef FIX
	BOOLEAN deleteRecord = FALSE;
#endif

	/* INITIALIZE MEMORY FOR MGD77 DATA STRUCTURES */
	struct MGD77_DATA_RECORD *D;
	struct MGD77_HEADER H;
	struct MGD77_CONTROL M, Out;
	struct MGD77_GRID_INFO this_grid[8];
	struct MGD77_ERROR *E = NULL;
	struct tm *systemTime;
	struct MGD77_CARTER C;
	struct GMT_gcal cal;
	
	/* INITIALIZE GMT & MGD77 MACHINERY */
	argc = GMT_begin (argc, argv);
	gmtdefs.time_system = 4;
	clock = time (NULL);
	maxTime = GMT_dt_from_usert ((double) clock);
	systemTime = gmtime (&clock);
	MGD77_Init (&M, TRUE);
	MGD77_Init (&Out, TRUE);
	MGD77_carter_init (&C);
	Out.fp = GMT_stdout;
	GMT_make_dnan (MGD77_NaN);
	
	/* INITIALIZE E77 */
	n_types[E77_NAV] = N_NAV_TYPES;
	n_types[E77_VALUE] = N_DEFAULT_TYPES;
	n_types[E77_SLOPE] = N_DEFAULT_TYPES;
	n_types[E77_GRID] = N_DEFAULT_TYPES;

	/* TURN ON MGD77SNIFFER ERROR MESSAGES */
	for (i = 0; i<MGD77_N_WARN_TYPES; i++) warn[i] = TRUE;
	
	/* SET PROGRAM DEFAULTS */
	arguments[0] = 0;
	for (i = 0; i < MGD77_N_DATA_FIELDS; i++) MGD77_this_bit[i] = 1 << i;
	maxGap = (double) MGD77_MAX_DS;	/* 5 km by default */
	n_out_columns = 0;			/* No formatted output */
	M.verbose_level = 3;	/* 1 = warnings, 2 = errors, 3 = both */
	M.verbose_dest = 1;		/* 1 = stdout, 2 = stderr */
	max_speed = MGD77_MAX_SPEED;
	derivative = "SPACE";
	for (i = 0; i<MGD77_N_NUMBER_FIELDS; i++) {
		maxSlope[i] = mgd77snifferdefs[i].maxSpaceGrad; /* Use spatial gradients by default */
		adjustScale[i] = MGD77_NaN;
		adjustDC[i] = MGD77_NaN;
	}
	memset ((void *)this_grid, 0, 8 * sizeof (struct MGD77_GRID_INFO));

	/* READ COMMAND LINE ARGUMENTS */
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			sprintf(arguments,"%s %s",arguments,argv[i]);
			switch (argv[i][1]) {
				case 'b':
				case 'V':
				case 'R':
				case '\0':
					error += GMT_parse_common_options (argv[i], &west, &east, &south, &north);
					break;
				case 'A':	/* adjust slope and intercept */
					if (!error && sscanf (&argv[i][2], "%[^,]", abbrev) != 1) {
						fprintf (stderr, "%s: SYNTAX ERROR -A option: Give field abbreviation, slope and intercept\n", GMT_program);
						error = TRUE;
					}
					/* Find what column number this field corresponds to (i.e. depth == 11) */
					col = 0; 
					while (strcmp (abbrev, mgd77defs[col].abbrev) && col < MGD77_N_NUMBER_FIELDS)
						col++;
					if (col == MGD77_N_NUMBER_FIELDS) {
						fprintf (stderr, "%s: SYNTAX ERROR -A option: invalid field abbreviation\n", GMT_program);
						error = TRUE;
					}
					if (!error && sscanf (&argv[i][2], "%[^,],%lf,%lf", abbrev, &adjustScale[col], &adjustDC[col]) != 3) {
						fprintf (stderr, "%s: SYNTAX ERROR -A option: Give field abbreviation,slope,intercept\n", GMT_program);
						error = TRUE;
					}
					adjustData = TRUE;
					break;
				case 'C':	/* set max speed */
					max_speed = atof (&argv[i][2]);
					custom_max_speed = TRUE;
					break;
				case 'D':
					if (argv[i][2] == 'd') { /* cruise - grid differences */
						display = "DIFFS";
						n_out_columns = 6;
					}
					else if (argv[i][2] == 'e') { /* E77 error output */
						display = "E77";
						n_out_columns = 6;
					}					
					else if (argv[i][2] == 'f') { /* deltaZ and deltaS for each field */
						display = "DFDS";
						n_out_columns = 20;
					}					
					else if (argv[i][2] == 'l') { /* Sniffer limits */
						display = "LIMITS";
						n_out_columns = 4;
					}
					else if (argv[i][2] == 'm') { /* MGD77 output */
						display = "MGD77";
						n_out_columns = 27;
					}
					else if (argv[i][2] == 's') { /* gradients */
						display = "SLOPES";
						n_out_columns = 11;
					}
					else if (argv[i][2] == 'v') { /* values */
						display = "VALS";
						n_out_columns = 12;
					}
					else {
						fprintf (stderr, "%s: SYNTAX ERROR:  Unrecognized option -%c%c\n", GMT_program,\
						argv[i][1], argv[i][2]);
						error = TRUE;
					}
					/* Silence all warning messages for data dumps */
					for (j = 0; j<MGD77_N_WARN_TYPES; j++) warn[j] = FALSE;
					M.verbose_dest = 2;		/* 1 = stdout, 2 = stderr */
					break;
				case 'F':	/* fake mode (specify field and constant z value in -G - no grid reading */
					simulate = TRUE;
					break;
				case 'g':	/* Get grid filename and geophysical field name to compare with grid */
					this_grid[n_grids].format = 1; 	/* Mercator grid */
					if (sscanf (&argv[i][2], "%[^,],%[^,],%lf,%d,%lf", this_grid[n_grids].abbrev, this_grid[n_grids].fname, &this_grid[n_grids].scale, &this_grid[n_grids].mode, &this_grid[n_grids].max_lat) < 4) {
						fprintf (stderr, "%s: SYNTAX ERROR -g option: Give field abbreviation, grdfile, scale, mode [, and optionally max lat]\n", GMT_program);
						error = TRUE;
					}
				case 'G':	/* Get grid filename and geophysical field name to compare with grid */
					if (!error && this_grid[n_grids].format == 0 && sscanf (&argv[i][2], "%[^,],%s", this_grid[n_grids].abbrev, this_grid[n_grids].fname) != 2) {
						fprintf (stderr, "%s: SYNTAX ERROR -G option: Give field abbreviation and grdfile\n", GMT_program);
						error = TRUE;
					}
					else {
						/* Find what column number this field corresponds to (i.e. depth == 11) */
						this_grid[n_grids].col = 0; 
						while (strcmp (this_grid[n_grids].abbrev, mgd77defs[this_grid[n_grids].col].abbrev) &&\
						this_grid[n_grids].col < MGD77_N_NUMBER_FIELDS)
							this_grid[n_grids].col++;
						if (this_grid[n_grids].col == MGD77_N_NUMBER_FIELDS) {
							fprintf (stderr, "%s: SYNTAX ERROR -G option: invalid field abbreviation\n", GMT_program);
							error = TRUE;
						}
						if (!strcmp (this_grid[n_grids].abbrev,"depth")) this_grid[n_grids].sign = -1;
						else this_grid[n_grids].sign = 1;
						n_grids++;
					}
					break;
				case 'H':	/* Avoid decimation during grid comparison */
					decimate = FALSE;
					break;
				case 'I':	/* Pass ranges of data records to ignore for output to E77 */
					if (!error && sscanf (&argv[i][2], "%[^,],%d,%d", BadSection[n_bad_sections].abbrev, &BadSection[n_bad_sections].start, &BadSection[n_bad_sections].stop) != 3) {
						fprintf (stderr, "%s: SYNTAX ERROR -A option: Give field abbreviation,slope,intercept\n", GMT_program);
						error = TRUE;
					}
					/* Find what column number this field corresponds to (i.e. depth == 11) */
					col = 0; 
					while (strcmp (BadSection[n_bad_sections].abbrev, mgd77defs[col].abbrev) && col < MGD77_N_NUMBER_FIELDS)
						col++;
					if (col == MGD77_N_NUMBER_FIELDS) {
						fprintf (stderr, "%s: SYNTAX ERROR -I option: invalid field abbreviation\n", GMT_program);
						error = TRUE;
					}
					bad_sections = TRUE;
					BadSection[n_bad_sections].col = col;
					n_bad_sections++;
					if (n_bad_sections == MAX_BAD_SECTIONS) {
						fprintf (stderr, "%s: SYNTAX ERROR -I option: Max number of sections (%d) reached\n", GMT_program, MAX_BAD_SECTIONS);
						error = TRUE;
					}
					break;
				case 'L':	/* Overwrite default sniffer limits */
					custom_limit_file = &argv[i][2];
					break;
				case 'N':	/* Change to nautical units instead of metric */
					nautical = TRUE;
					speed_units = "knots";
					distance_units = "nm";
					break;
				case 'Q':	/* bilinear interpolation parameters */
					bilinear = TRUE;
					threshold = (argv[i][2]) ? atof (&argv[i][2]) : 1.0;
					break;
				case 'S':	/* Specify spatial, time, or simple difference gradients */
					if (argv[i][2] == 'd') {
						derivative = "DIFF";
						for (j = 0; j<MGD77_N_NUMBER_FIELDS; j++) maxSlope[j] = mgd77snifferdefs[j].maxTimeGrad;
					}
					else if (argv[i][2] == 's') {
						derivative = "SPACE";
						for (j = 0; j<MGD77_N_NUMBER_FIELDS; j++) maxSlope[j] = mgd77snifferdefs[j].maxSpaceGrad;
					}
					else if (argv[i][2] == 't') {
						derivative = "GMT_TIME";
						for (j = 0; j<MGD77_N_NUMBER_FIELDS; j++) maxSlope[j] = mgd77snifferdefs[j].maxTimeGrad;
					}
					else {
						fprintf (stderr, "%s: SYNTAX ERROR:  Unrecognized option -%c%c\n", GMT_program, argv[i][1],\
						argv[i][2]);
						error = TRUE;
					}	
					break;
				case 'T':	/* Specify maximum gap between records */
					custom_maxGap = TRUE;
					maxGap = atof (&argv[i][2]);
					if (maxGap < 0) {
						fprintf (stderr, "%s: SYNTAX ERROR -M option: max gap cannot be negative\n", GMT_program);
						error = TRUE;
					}
					break;					
				case 'W':	/* Choose which warning types to go to stdout (default - all) */
					for (j = 0; j<MGD77_N_WARN_TYPES; j++) warn[j] = FALSE;
					while (GMT_strtok (&argv[i][2], ",", &pos, &c)) {
						if (c == 'v')
							warn[VALUE_WARN] = TRUE;
						else if (c == 'g')
							warn[SLOPE_WARN] = TRUE;
						else if (c == 'o')
							warn[GRID_WARN] = TRUE;
						else if (c == 't')
							warn[TIME_WARN] = TRUE;
						else if (c == 's')
							warn[SPEED_WARN] = TRUE;
						else if (c == 'c')
							warn[TYPE_WARN] = TRUE;
						else if (c == 'x')
							warn[SUMMARY_WARN] = TRUE;
						else {
							fprintf (stderr, "%s: SYNTAX ERROR:  Unrecognized option -%c%c\n", GMT_program,\
								 argv[i][1], argv[i][2]);
							error = TRUE;
						}
					}
					custom_warn = TRUE;
					break;					
				default:
					fprintf (stderr, "%s: SYNTAX ERROR:  Unrecognized option -%c\n", GMT_program, argv[i][1]);
					error = TRUE;
					break;
			}
		}
		else
			n_cruises++;
	}

	/* DISPLAY ERROR PAGE */
	if (GMT_give_synopsis_and_exit || argc == 1) {	/* Display usage */
		fprintf(stderr,"mgd77sniffer %s - Along-track quality control of MGD77 cruises\n\n", MGD77_VERSION);
		fprintf (stderr, "usage:  mgd77sniffer <cruises> [-Afieldabbrev,scale,offset] [-Cmaxspd] [-Dd|e|f|l|m|s|v]\n");
		fprintf (stderr, "\t[-gfieldabbrev,imggrid,scale,mode[,latmax]] [-Gfieldabbrev,grid] [-H] [-Ifieldabbrev,rec1,recN] [-Lcustom_limits_file]\n");
		fprintf (stderr, "\t[-N] [-Q[<value>]] [%s] [-Sd|s|t] [-Tgap] [-Wc|g|o|s|t|v|x] [-V] [%s]\n\n", GMT_Rgeo_OPT, GMT_bo_OPT);
		if (GMT_give_synopsis_and_exit) exit (EXIT_FAILURE);
		fprintf (stderr, "\tScan MGD77 files for errors using point-by-point sanity checking,\n");
		fprintf (stderr, "\t\talong-track detection of excessive slopes and comparison of cruise\n");
		fprintf (stderr, "\t\tdata with global bathymetry and gravity grids.");
		fprintf (stderr, "\twhere <cruises> is one or more MGD77 legnames, e.g. 08010001 etc.\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-A Apply scale factor and DC adjustment to specified data field. Allows adjustment of\n");
		fprintf (stderr, "\t   cruise data prior to along-track analysis. CAUTION: data must be thoroughly examined\n");
		fprintf (stderr, "\t   before applying these global data adjustments. May not be used for multiple cruises.\n");		
		fprintf (stderr, "\t-C Set maximum ship speed (10 m/s by default, use -N to indicate knots)\n");
		fprintf (stderr, "\t-D Dump cruise data such as sniffer limits, values, gradients and mgd77 records.\n");
		fprintf (stderr, "\t  -Dd print out cruise-grid differences (requires -G option)\n");
		fprintf (stderr, "\t  -De output formatted error summary for each record. See E77 ERROR FORMAT below.\n");	
		fprintf (stderr, "\t  -Df for each field, output value change and distance since last observation.\n");	
		fprintf (stderr, "\t  -Dl print out mgd77sniffer default limits (requires no additional arguments)\n");
		fprintf (stderr, "\t  -Dm print out MGD77 format\n\t  -Ds print out gradients\n\t  -Dv print out values\n");
		fprintf (stderr, "\t-g Compare cruise data to the specified Sandwell/Smith Mercator grid. Requires valid MGD77\n");
		fprintf (stderr, "\t   field abbreviation followed by a comma, the path (if not in current directory)\n");
		fprintf (stderr, "\t   and grid filename, scale (0.1 or 1), and mode (see mgd77manage for details)\n");
		fprintf (stderr, "\t   Optionally, append max latitude in the IMG file [72.0059773539]\n");
		fprintf (stderr, "\t-G Compare cruise data to the specified GMT geographic grid. Requires valid MGD77 field abbreviation\n");
		fprintf (stderr, "\t   followed by a comma, then the path (if not in current directory) and grid filename.\n");
		fprintf (stderr, "\t   Excessive offsets are flagged according to maxArea threshold (use -L option to\n");
		fprintf (stderr, "\t   adjust maxArea). Useful for comparing faa or depth to global grids though any MGD77\n");
		fprintf (stderr, "\t   field can be compared to any GMT or IMG compatible grid. Multiple grid comparison is\n");
		fprintf (stderr, "\t   supported by  using  separate -G or -g calls for each grid. See GRID FILE INFO below.\n");
		fprintf (stderr, "\t-H (with -G|g only) disable data decimation during RLS regression..\n");
		fprintf (stderr, "\t-I Give one or more times to specify ranges of data record that should be flagged as bad\n");
		fprintf (stderr, "\t   prior to along-track analysis.  The flag information will be echoed out to E77 files.\n");		
		fprintf (stderr, "\t   May not be used for multiple cruises.\n");		
		fprintf (stderr, "\t-L Override mgd77sniffer default error detection limits. Supply path and filename of\n");
		fprintf (stderr, "\t   the custom limits file. Rows not beginning with a valid MGD77 field abbreviation are\n");
		fprintf (stderr, "\t   ignored. Field abbreviations are listed below in exact form under MGD77 FIELD INFO.\n");
		fprintf (stderr, "\t   Multiple field limits may be modified using one default file, one field per line.\n");
		fprintf (stderr, "\t   Field min, max, maxGradient and maxArea may be changed for each field. maxGradient\n");
		fprintf (stderr, "\t   pertains to the gradient type selected using the -S option. maxArea is used by the\n");
		fprintf (stderr, "\t   -G option as the threshold for flagging excessive offsets. Dump defaults (-Dd) to\n");
		fprintf (stderr, "\t   view syntax or to quickly create an editable custom limits file.\n");
		fprintf (stderr, "\t   Example custom default file contents (see below for field units):\n");
		fprintf (stderr, "\t\tdepth	0	11000	1000	4500\n");
		fprintf (stderr, "\t\tmag	-800	800	-	-\n");
		fprintf (stderr, "\t\tfaa	-250	250	100	2500\n");
		fprintf (stderr, "\t   Use a dash '-' to retain a default limit.\n");
		fprintf (stderr, "\t   Hint: to test your custom limits, try: mgd77sniffer -Dl -L<yourlimitsfile>\n");
		fprintf (stderr, "\t-N Use nautical units.\n");
		fprintf (stderr, "\t-Q Quick mode, use bilinear rather than bicubic interpolation (requires -g or -G).\n");
		fprintf (stderr, "\t   Optionally, append <value> in the 0 <= value <= 1 range.\n");
		fprintf (stderr, "\t   [Default = 1 requires all 4 nearest nodes to be non-NaN.], <value> = 0.5\n");
		fprintf (stderr, "\t   will interpolate about 1/2 way from a non-NaN to a NaN node, while\n");
		fprintf (stderr, "\t   0.1 will go about 90%% of the way, etc.\n");
		fprintf (stderr, "\t   -Q0 means only consider the nearest node (no interpolation).\n");
		fprintf (stderr, "\t-S Specify gradient type for along-track excessive slope  checking.\n");
		fprintf (stderr, "\t  -Sd Calculate change in z values along track (dz)\n");
		fprintf (stderr, "\t  -Ss Calculate spatial gradients (dz/ds) [default]\n");
		fprintf (stderr, "\t  -St Calculate time gradients (dz/dt)\n");
		fprintf (stderr, "\t-T Set maximum acceptable distance gap between records (km) [5].\n");
		fprintf (stderr, "\t   Set to zero to deactivate gap checking .\n");
		fprintf (stderr, "\t-W Print out only certain warning types. Comma delimit any combination of c|g|o|s|t|v|x:\n");
		fprintf (stderr, "\t   where (c) type code warnings, (g)radient out of range, (o)ffsets from grid (requires -G),\n");
		fprintf (stderr, "\t   (s)peed out of range, (t)ime warnings, (v)alue out of range, (x) warning summaries.\n");
		fprintf (stderr, "\t   By default ALL warning messages are printed. Not allowed with -D option.\n");
		fprintf (stderr, "\t-V runs in verbose mode.\n\n");
		fprintf (stderr, "\t-b output binary data for -D option.  Append d for double and s for single precision [double].\n\n");		
		fprintf (stderr, "\tMGD77 FIELD INFO:\n");
		fprintf (stderr, "\tField\t\t\tAbbreviation\t\tUnits\n");
		fprintf (stderr, "\tTwo-way Travel Time\ttwt\t\t\tsec\n");
		fprintf (stderr, "\tCorrected Depth \tdepth\t\t\tm\n");
		fprintf (stderr, "\tMag Total Field1\tmtf1\t\t\tnT\n");
		fprintf (stderr, "\tMag Total Field2\tmtf2\t\t\tnT\n");
		fprintf (stderr, "\tResidual Magnetic\tmag\t\t\tnT\n");
		fprintf (stderr, "\tDiurnal Correction\tdiur\t\t\tnT\n");
		fprintf (stderr, "\tMag Sensor Depth/Alt\tmsd\t\t\tm\n");		
		fprintf (stderr, "\tObserved Gravity\tgobs\t\t\tmGal\n");
		fprintf (stderr, "\tEotvos Correction\teot\t\t\tmGal\n");		
		fprintf (stderr, "\tfree-air Anomaly\tfaa\t\t\tmGal\n\n");
		fprintf (stderr, "\tGRID FILE INFO:\n");
		fprintf (stderr, "\t-g: Img files must be of Sandwell/Smith signed two-byte integer (i2) type with no header.\n");
		fprintf (stderr, "\t-G: Grid files can be any type of GMT grid file (native or netCDF) with header\n");
		fprintf (stderr, "\tA correctly formatted grid file can be generated as follows:\n");
		fprintf (stderr, "\t   e.g. gmtset GRIDFILE_SHORTHAND TRUE\n");
		fprintf (stderr, "\t\tCreate/edit .gmt_io file to include the following rows:\n");
		fprintf (stderr, "\t\t\t# GMT I/O shorthand file\n"); 
		fprintf (stderr, "\t\t\t# suffix   format_id scale offset       NaN\n");
		fprintf (stderr, "\t\t\tgrd             0       -       -       -\n");
		fprintf (stderr, "\t\t\ti2              2       -       -       32767\n");
		fprintf (stderr, "\t\tgrdraster 1 -R0/359:55/-90/90 -Getopo5_hdr.i2\n\n");
		fprintf (stderr, "E77 ERROR OUTPUT\n");
		fprintf (stderr, "\tError output is divided into (1) a header containing information globally\n");
		fprintf (stderr, "\tapplicable to the cruise and (2) individual error records summarizing all\n");
		fprintf (stderr, "\tall  errors  encountered in each cruise record. One error record is output\n");
		fprintf (stderr, "\tfor each  readable  cruise  record.\n");
		fprintf (stderr, "\tError Record Format: <time/distance>  <record  number>  <error code string> <description>\n");
		fprintf (stderr, "Example:\n\t# Cruise 06050010 ID L476WG MGD77 FILE VERSION: 19870415 N_RECS: 27268\n");
		fprintf (stderr, "\t# Examined: Tue Feb 21 16:33:34 2006\n");
		fprintf (stderr, "\t# Examiner: mtchandl\n");
		fprintf (stderr, "\t# Arguments:  -De -Gfaa,/data/GRIDS/sandwell_smith11_hdr.i2\n");
		fprintf (stderr, "\t# Errata: Header\n");
		fprintf (stderr, "\tN-faa-E-01: Regression scale 0.108643 different from 1. Recommended: [10]\n");
		fprintf (stderr, "\tN-faa-E-02: Regression offset -1.319411 different from 0\n");
		fprintf (stderr, "\tI-twt-W-03: More recent bathymetry correction table available\n");
		fprintf (stderr, "\tI-depth-W-04: Integer precision in depth\n");
		fprintf (stderr, "\t# Errata: Data\n");
		fprintf (stderr, "\t1976-06-26T08:15:00	2	C-0-0-0	 nav: excessive speed\n");
		fprintf (stderr, "\t1976-06-26T08:16:00	3	C-0-0-0	 nav: excessive speed\n");
		fprintf (stderr, "\n\tError Class Descriptions\n");
		fprintf (stderr, "\tNAV (navigation):\t0 --> fine\n");
		fprintf (stderr, "\t\tA --> time out of range\n");
		fprintf (stderr, "\t\tB --> time decreasing\n");
		fprintf (stderr, "\t\tC --> excessive speed\n");
		fprintf (stderr, "\t\tD --> above sea level\n");
		fprintf (stderr, "\t\tE --> lat undefined\n");
		fprintf (stderr, "\t\tF --> lon undefined\n\n");
		fprintf (stderr, "\tVAL (value)\t0 --> fine\n");
		fprintf (stderr, "\t\tK --> twt invalid\n");
		fprintf (stderr, "\t\tL --> depth invalid\n");
		fprintf (stderr, "\t\tO --> mtf1 invalid\n");
		fprintf (stderr, "\t\tetc.\n\n");
		fprintf (stderr, "\tGRAD (gradient):\t0 --> fine\n");
		fprintf (stderr, "\t\tK --> d[twt] excessive\n");
		fprintf (stderr, "\t\tL --> d[depth] excessive\n");
		fprintf (stderr, "\t\tO --> d[mtf1] excessive\n");
		fprintf (stderr, "\t\tetc.\n\n");
		fprintf (stderr, "\tGRID (grid comparison):\t0 --> fine\n");
		fprintf (stderr, "\t\tL --> depth offset\n");
		fprintf (stderr, "\t\tW --> faa offset\n");
		fprintf (stderr, "\t\tetc.\n");
		fprintf (stderr, "\nEXAMPLES:\n\tAlong-track excessive value and gradient checking:\n\t\tmgd77sniffer 08010001\n");
		fprintf (stderr, "\tDump cruise gradients:\n\t\tmgd77sniffer 08010001 -Ds\n"); 
		fprintf (stderr, "\tTo compare cruise depth with ETOPO5 bathymetry and gravity with Sandwell/Smith 2 min gravity version 11, try\n");
		fprintf (stderr, "\t\tmgd77sniffer 08010001 -Gdepth,/data/GRIDS/etopo5_hdr.i2 -gfaa,/data/GRIDS/grav.11.2.img,0.1,1\n\n"); 
		exit (EXIT_FAILURE);
	}
	/* ENSURE VALID USE OF OPTIONS */
	if (n_cruises == 0 && display != "LIMITS") {
		fprintf(stderr, "%s: ERROR: No cruises given\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	else if (GMT_io.binary[GMT_OUT] && !display) {
		fprintf(stderr, "%s: ERROR: -b option requires -D.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	else if (custom_warn && display) {
		fprintf(stderr, "%s: ERROR: Incompatible options -D and -W.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	else if (display == "DIFFS" && n_grids == 0) {
		fprintf(stderr, "%s: ERROR: -Dd option requires -G.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (east < west || south > north) {
		fprintf(stderr, "%s: ERROR: Region set incorrectly\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (adjustData && n_cruises > 1) {
		fprintf(stderr, "%s: ERROR: -A adjustments valid for only one cruise.\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	/* NAUTICAL CONVERSION FACTORS */
	if (nautical) {
		distance_factor = 1.0 / MGD77_METERS_PER_NM; /* meters to nm */
		time_factor = 1.0 / (3600); /* seconds to hours */
		/* Adjust max speed for unit change and user specified max speed */
		if (!custom_max_speed) max_speed = MGD77_MAX_SPEED * distance_factor / time_factor;
	}

	/* READ AND APPLY CUSTOM LIMITS FILE */
	mgd77snifferdefs[MGD77_YEAR].maxValue = (double) systemTime->tm_year + 1900;
	if (custom_limit_file) {
		if ((custom_fp = fopen (custom_limit_file, "r")) == NULL) {
			fprintf (stderr, "%s: Could not open custom limit file %s\n", GMT_program, custom_limit_file);
			exit (EXIT_FAILURE);
	 	}
		else {
			while (fgets (custom_limit_line, BUFSIZ, custom_fp)) {
				if (sscanf (custom_limit_line,"%s %s %s %s %s", field_abbrev, tmp_min, tmp_max, tmp_maxSlope, tmp_area) == 5) {
					i = 0;
					while (strcmp (mgd77snifferdefs[i].abbrev, field_abbrev) && i <= MGD77_N_NUMBER_FIELDS) i++;
					if (i <= MGD77_N_NUMBER_FIELDS) {
						if (strcmp(tmp_min, "-"))  mgd77snifferdefs[i].minValue = atof (tmp_min);
						if (strcmp(tmp_max, "-"))  mgd77snifferdefs[i].maxValue = atof (tmp_max);
						if (strcmp(tmp_maxSlope, "-")) maxSlope[i] = atof (tmp_maxSlope);
						if (strcmp(tmp_area, "-"))  mgd77snifferdefs[i].maxArea = atof (tmp_area);
					}
				}
				else {
					fprintf (stderr, "%s: Error in custom limits file [%s]\n", GMT_program, custom_limit_line);
					exit (EXIT_FAILURE);
				}
			}
		}
	}

	/* Use GMT time formatting */
	GMT_io.out_col_type[MGD77_TIME] = GMT_IS_ABSTIME;

	/* Adjust data dump for number of grids */
	if (display == "DIFFS")
		n_out_columns = 3+3*n_grids;

	/* PREPARE DUMP OUTPUT FORMATTING */
	if (display == "VALS" || display == "DIFFS") {
		for (i = MGD77_LATITUDE; i < (n_out_columns + MGD77_LATITUDE); i++) {
			/* Special lat formatting */
			if (i == MGD77_LATITUDE)
				GMT_io.out_col_type[i-MGD77_LATITUDE] = GMT_IS_LAT;

			/* Special lon formatting */
			else if (i == MGD77_LONGITUDE)
				GMT_io.out_col_type[i-MGD77_LATITUDE] = GMT_IS_LON;

			/* Everything else is float */
			else
				GMT_io.out_col_type[i-MGD77_LATITUDE] = GMT_IS_FLOAT;
		}
	}
	else if (display) {
		for (i = 0; i < n_out_columns; i++) {	
			/* All columns are floats */
			GMT_io.out_col_type[i] = GMT_IS_FLOAT;
		}
	}

	/* PRINT OUT DEFAULT GEOPHYSICAL LIMITS */
	if (display == "LIMITS") {
		fprintf (GMT_stdout, "#abbrev\tmin\tmax\tmaxSlope\tmaxArea\n");
		for (i = 0; i < MGD77_N_NUMBER_FIELDS; i++) {
			if ((1 << i) & (MGD77_GEOPHYSICAL_BITS + MGD77_CORRECTION_BITS)) {
				fprintf (GMT_stdout, "%s%s",mgd77defs[i].abbrev,gmtdefs.field_delimiter);
				GMT_ascii_output_one (GMT_stdout, mgd77snifferdefs[i].minValue, 0);
				fprintf (GMT_stdout, "%s",gmtdefs.field_delimiter);
				GMT_ascii_output_one (GMT_stdout, mgd77snifferdefs[i].maxValue, 1);
				fprintf (GMT_stdout, "%s",gmtdefs.field_delimiter);
				GMT_ascii_output_one (GMT_stdout, maxSlope[i], 2);
				fprintf (GMT_stdout, "%s",gmtdefs.field_delimiter);
				GMT_ascii_output_one (GMT_stdout, mgd77snifferdefs[i].maxArea, 3);
				fprintf (GMT_stdout, "\n");
			}
		}
		exit(EXIT_FAILURE);
	}

	/* PRINT HEADER FOR SNIFFER DATA DUMPS */
	if (display && !GMT_io.binary[GMT_OUT] && GMT_io.io_header[0]) {
		if (display == "SLOPES") {
			fprintf (GMT_stdout, "#speed(%s)%s",speed_units,gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[twt]%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[depth]%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[mtf1]%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[mtf2]%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[mag]%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[diur]%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[msd]%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[gobs]%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[eot]%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[faa]\n");
		}
		else if (display == "DFDS") {
			fprintf (GMT_stdout, "#d[twt]%sds%s",gmtdefs.field_delimiter,gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[depth]%sds%s",gmtdefs.field_delimiter,gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[mtf1]%sds%s",gmtdefs.field_delimiter,gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[mtf2]%sds%s",gmtdefs.field_delimiter,gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[mag]%sds%s",gmtdefs.field_delimiter,gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[diur]%sds%s",gmtdefs.field_delimiter,gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[msd]%sds%s",gmtdefs.field_delimiter,gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[gobs]%sds%s",gmtdefs.field_delimiter,gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[eot]%sds%s",gmtdefs.field_delimiter,gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "d[faa]%sds\n",gmtdefs.field_delimiter);
		}
		else if (display == "VALS") {
			fprintf (GMT_stdout, "#lat%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "lon%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "twt%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "depth%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "mtf1%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "mtf2%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "mag%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "diur%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "msd%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "gobs%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "eot%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "faa\n");
		}
		else if (display == "DIFFS" && n_grids > 0) {
			fprintf (GMT_stdout, "#lat%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "lon%s",gmtdefs.field_delimiter);
			fprintf (GMT_stdout, "dist%s",gmtdefs.field_delimiter);
			for (i = 0; i < n_grids; i++) {
				fprintf (GMT_stdout, "crs_%s%s",this_grid[i].abbrev,gmtdefs.field_delimiter);
				fprintf (GMT_stdout, "grd_%s%s",this_grid[i].abbrev,gmtdefs.field_delimiter);
				fprintf (GMT_stdout, "diff%s",gmtdefs.field_delimiter);
			}
			fprintf (GMT_stdout, "\n");
		}
	}

	/* Open grid files */
	for (i = 0; i < n_grids; i++) {
		if (!simulate)
			/* Open and store grid file */
			read_grid (&this_grid[i], &f[i], west, east, south, north, bilinear, threshold);
	}

	if (n_grids) {
		MaxDiff = (double *) GMT_memory (VNULL, (size_t)n_grids, sizeof (double), "mgd77sniffer");
		iMaxDiff = (int *) GMT_memory (VNULL, (size_t)n_grids, sizeof (int), "mgd77sniffer");
	}

	MGD77_Ignore_Format (MGD77_FORMAT_ANY);	/* Reset to all formats OK, then ... */
	MGD77_Ignore_Format (MGD77_FORMAT_CDF);	/* disallow netCDF MGD77+ files */
	MGD77_Ignore_Format (MGD77_FORMAT_TBL);	/* and plain ASCI tables */

	/* PROCESS CRUISES */
	for (argno = 1; argno < argc; argno++) {

		if (argv[argno][0] == '-') 
			continue;

		/* Verify valid cruise file */
  		if (MGD77_Get_Path (file, argv[argno], &M)) {
   			fprintf (stderr, "%s : Cannot find leg %s\n", GMT_program, argv[argno]);
			continue;
  		}

		/* Open cruise file */
		if (MGD77_Open_File (argv[argno], &M, 0)) 
			continue;

		if (gmtdefs.verbose) 
			fprintf (stderr, "%s: Processing file %s\n", GMT_program, argv[argno]);

		if (display == "E77") {
			sprintf (outfile,"%s.e77",M.NGDC_id);
			if ((fpout = fopen (outfile, "w")) == NULL) {
				fprintf (stderr, "%s: Could not open E77 output file %s\n", GMT_program, outfile);
				exit (EXIT_FAILURE);
			}
	 	}

		/* Read MGD77 header */
		if (MGD77_Read_Header_Record_asc (argv[argno], &M, &H))
			fprintf (stderr, "%s: Cruise %s has no header.\n", GMT_program, argv[argno]);

		/* Allocate memory for data records */
		D = (struct MGD77_DATA_RECORD *) GMT_memory (VNULL, (size_t)n_alloc, \
		sizeof (struct MGD77_DATA_RECORD), "mgd77sniffer");

		/* READ DATA RECORDS */
		n_alloc = GMT_CHUNK;
		gotTime = FALSE;
		nvalues = M.bit_pattern[0] = 0;
		lowPrecision = lowPrecision5 = 0;
		while (!MGD77_Read_Data_Record_m77 (&M, &D[nvalues])) {
			/* Increase memory allocation if necessary */
			if (nvalues == n_alloc - 1) {
				n_alloc += GMT_CHUNK;
				D = (struct MGD77_DATA_RECORD *) GMT_memory ((void *)D, (size_t)n_alloc, \
				sizeof (struct MGD77_DATA_RECORD), "mgd77sniffer");
			}
			if (!gotTime && !GMT_is_dnan(D[nvalues].time))
				gotTime = TRUE;
			M.bit_pattern[0] |= D[nvalues].bit_pattern;
			nvalues++;
		}
		
		/* Scale and DC adjust if selected */
		if (adjustData) {
			for (i=0;i<MGD77_N_NUMBER_FIELDS;i++) {
				if (!GMT_is_dnan(adjustScale[i]) || !GMT_is_dnan(adjustDC[i])) {
					for (j=0; j<nvalues;j++)
						D[j].number[i] = D[j].number[i] * adjustScale[i] + adjustDC[i];
					fprintf (stderr,"%s (%s) Scaled by %f and %f added\n", GMT_program, mgd77defs[i].abbrev,\
					adjustScale[i],adjustDC[i]);
				}
			}
		}

		/* Set user-specified flagged observations to NaN before analysis */
		if (bad_sections) {
			for (j=0;j<(int)n_bad_sections;j++) {	/* For each bad section */
				for (i=BadSection[j].start-1;i<BadSection[j].stop;i++) {	/* Loop over the flagged records (adjust -1 for C index) */
					D[i].number[BadSection[j].col] = MGD77_NaN;	/* and set them to NaN */
				}
			}
		}

		/* Output beginning of E77 header */
		if (display == "E77") {
			fprintf (fpout, "# Cruise %s ID %s MGD77 FILE VERSION: %.4d%2.2d%2.2d N_RECS: %d\n",argv[argno],D[0].word[0],\
			atoi(H.mgd77->File_Creation_Year),atoi(H.mgd77->File_Creation_Month),atoi(H.mgd77->File_Creation_Day),nvalues);
			fprintf (fpout,"# Examined: %s",ctime(&clock));
			fprintf (fpout,"# Examiner: %s\n",M.user);
			fprintf (fpout,"# Arguments: %s\n",arguments);
			fprintf (fpout,"# Errata: Header\n");
		}

		/* Scan MGD77 header */
		if (display == "E77" || warn[SUMMARY_WARN]) {
			MGD77_Verify_Prep_m77 (&M, &H.meta, D, nvalues);	/* First get key meta-data derived form data records */
			MGD77_Verify_Header (&M, &H, fpout);				/* Then verify the header information */
		}

		/* Re-set variables for this cruise */
		landcruise = FALSE;
		overLandCount = overLandStart = 0;
		timeErrorStart = noTimeStart = distanceErrorStart = -1;
		noTimeCount = timeErrorCount = distanceErrorCount = n_nan = bccCode = 0;
		offsetArea = (double *) GMT_memory (VNULL, (size_t)n_grids, sizeof (double), "mgd77sniffer");
		offsetStart = (int *)GMT_memory (VNULL, (size_t)n_grids, sizeof (int), "mgd77sniffer");
		offsetLength = (double *) GMT_memory (VNULL, (size_t)n_grids, sizeof (double), "mgd77sniffer");
		offsetSign = (BOOLEAN *)GMT_memory (VNULL, (size_t)n_grids, sizeof (BOOLEAN), "mgd77sniffer");
		prevOffsetSign = (BOOLEAN *)GMT_memory (VNULL, (size_t)n_grids, sizeof (BOOLEAN), "mgd77sniffer");
		range = range2 = S_xx = date = n_days = 0.0;
		prevFlag = FALSE;
		mtf1 = TRUE;
		for (i = 0; i<n_grids; i++) {
			offsetArea[i] = 0.0;
			offsetStart[i] = 0;
			offsetSign[i] = prevOffsetSign[i] = FALSE;
		}
		for (i = MGD77_LATITUDE; i<MGD77_N_NUMBER_FIELDS; i++) {
			if (i == MGD77_MSD) continue;
			/* Turn on low precision bits (will be turned off later if ok) */
			if ((M.bit_pattern[0] & (1 << i) & MGD77_FLOAT_BITS) && i != MGD77_MSD) lowPrecision |= (1 << i);
			if (M.bit_pattern[0] & (1 << i) & MGD77_FLOAT_BITS) lowPrecision5 |= (1 << i);
			duplicates[i] = 0;
		}
		/* Adjust along-track gradient type for time */
		if (derivative == "GMT_TIME" && !gotTime) {
			/*derivative = "SPACE";*/
			if (warn[TIME_WARN])
				fprintf (GMT_stdout, "%s Warning: cruise contains no time - time gradients invalid.\n",GMT_program);
		}

		/* Allocate memory for distance array */
		distance = (double *) GMT_memory (VNULL, (size_t)nvalues, \
		sizeof (double), "mgd77sniffer");
		distance[0] = 0;
		
		/* Setup output array */
		out = (double **)GMT_memory (VNULL, nvalues, sizeof (double *), GMT_program);
		
		/* Allocate memory for error array */
		E = (struct MGD77_ERROR *) GMT_memory ((void *)E, (size_t)nvalues, \
		sizeof (struct MGD77_ERROR), "mgd77sniffer");

		/* PROCESS GRID FILES */
		if (n_grids > 0) {

			/* Allocate memory for 2D arrays */
			G = (double **)GMT_memory (VNULL, n_grids, sizeof (double *), GMT_program);	/* grid z values */
			diff = (double **)GMT_memory (VNULL, n_grids, sizeof (double *), GMT_program);	/* cruise-grid differences */

			for (i = 0; i < n_grids; i++) {

				/* Initialize variables */
				for (k=0; k<6; k++) stat[k] = stat2[k] = 0.0;
				tcrit = se = 0;
				newScale = FALSE;
				MaxDiff[i] = 0.0;

				this_grid[i].g_pts = 0;
				/* Skip if cruise lacks data field */
				if (!(M.bit_pattern[0] & (1 << this_grid[i].col))) {
					fprintf (stderr, "%s: Warning:  %s field not present in MGD77 file\n", GMT_program, this_grid[i].abbrev);
					if (this_grid[i].col != MGD77_DEPTH)
						continue;
				}

				/* Allocate memory for current array */
				G[i] = (double *) GMT_memory (VNULL, (size_t)nvalues, sizeof (double), "mgd77sniffer");
				diff[i] = (double *) GMT_memory (VNULL, (size_t)nvalues, sizeof (double), "mgd77sniffer");

				/* Sample grid at each ship location */
				if (simulate) {
					this_grid[i].g_pts = nvalues;
					for (j = 0; j < nvalues; j++)
						G[i][j] = 0.0;
				}
				else
					this_grid[i].g_pts = sample_grid (&this_grid[i], D, G, f[i], i, nvalues);
					
				if (this_grid[i].g_pts < 2) {
						fprintf (stderr, "%s: Insufficient grid samples for %s comparison\n", GMT_program, this_grid[i].abbrev);
						continue;
				}

				/* Reverse grid sign if depth */
				if (this_grid[i].sign == -1) {
					for (j = 0; j < nvalues; j++)
						G[i][j] *= this_grid[i].sign;
				}

				for (j = 0; j < nvalues; j++) {
					/* Compute cruise - grid differences */
					diff[i][j] = D[j].number[this_grid[i].col] - G[i][j];				

					/* Count NaNs */
					if (GMT_is_dnan(D[j].number[this_grid[i].col]) || GMT_is_dnan(G[i][j]))
						n_nan++;
				}

				/* SHIP VS GRID RLS REGRESSION */
				/* Allocate memory for NaN-free arrays */
				ship_val = (double *) GMT_memory (VNULL, (size_t)nvalues-n_nan, sizeof (double), "mgd77sniffer");
				grid_val = (double *) GMT_memory (VNULL, (size_t)nvalues-n_nan, sizeof (double), "mgd77sniffer");

				/* Store grid/cruise pairs in NaN-free arrays */
				for (j = k = 0; j < nvalues-n_nan; j++) {
					if (!GMT_is_dnan(D[j].number[this_grid[i].col]) && !GMT_is_dnan(G[i][j])) {
						ship_val[k] = D[j].number[this_grid[i].col];
						grid_val[k] = G[i][j];
						k++;
					}
				}
					
				/* Get max and min z values for this field */
				min = mgd77snifferdefs[this_grid[i].col].minValue;
				max = mgd77snifferdefs[this_grid[i].col].maxValue;
				if (this_grid[i].col == MGD77_FAA) {
					min *= 2;
					max *= 2;
				}
					
				/* Create a 2-D bin table */
				n = irint ((max - min)/mgd77snifferdefs[this_grid[i].col].delta) + 1;
				bin2d = (int **)GMT_memory (VNULL, n, sizeof (int *), GMT_program);
				for (j = 0; j < n; j++) 
					bin2d[j] = (int *)GMT_memory (VNULL, n,  sizeof (int), GMT_program);

				/* Then loop over all the ship,cruise pairs */
				for (j = 0; j < nvalues-n_nan; j++) {
					/* Need to skip ship values that are outside of acceptable range */
					if (ship_val[j] >= min && ship_val[j] <= max) {
						ship_bin = irint ((ship_val[j] - min)/mgd77snifferdefs[this_grid[i].col].delta);
						grid_bin = irint ((grid_val[j] - min)/mgd77snifferdefs[this_grid[i].col].delta);
						bin2d[ship_bin][grid_bin]++;    /* Add up # of pairs in this bin */
					}
				}
					
				/* Then find how many binned pairs we got */
				for (ship_bin = npts = 0; ship_bin < n; ship_bin++) {
				    for (grid_bin = 0; grid_bin < n; grid_bin++) {
       					if (bin2d[ship_bin][grid_bin] > 0)
							npts++;
					}
				}

				/* When all ship data are outside grid range skip regression */
				if (npts > 2) {

					/* Then create arrays for passing to RLS */
					decimated_ship = (double *)GMT_memory (VNULL, npts, sizeof (double), GMT_program);
					decimated_grid = (double *)GMT_memory (VNULL, npts, sizeof (double), GMT_program);

					for (ship_bin  = k = 0; ship_bin < n; ship_bin ++) {
						for (grid_bin = 0; grid_bin < n; grid_bin ++) {
							if (bin2d[ship_bin][grid_bin]) {
								decimated_ship[k] = min + ship_bin * mgd77snifferdefs[this_grid[i].col].delta;
								decimated_grid[k] = min + grid_bin * mgd77snifferdefs[this_grid[i].col].delta;
								k++;	/* Count number of non-empty bins */
							}
						}
					}

					/* Decimation benefits marine gravity due to amplitude differences */
					/* between ship and satellite data and also broadens confidence  */
					/* intervals for any ship grid comparisons by reducing excessive */
					/* number of degrees of freedom */
					if (gmtdefs.verbose) 
						fprintf (stderr, "%s: Begin RLS Regression\n", GMT_program);
					if (!decimate) {
						regress_rls (grid_val, ship_val, nvalues-n_nan, stat, this_grid[i].col, S_xx);
						decimated = FALSE;
						tcrit = GMT_tcrit (0.975, (double)nvalues-n_nan - 2.0);
						npts=nvalues-n_nan;
					}
					else {
						regress_rls (decimated_grid, decimated_ship, npts, stat, this_grid[i].col, S_xx);
						decimated = TRUE;
						regress_rls (grid_val, ship_val, nvalues-n_nan, stat2, this_grid[i].col, S_xx);
						if (stat[4] < stat2[4] && stat2[5] == 1.0) {
							if (gmtdefs.verbose)
								fprintf (stderr, "%s: Regression on undecimated data due to better correlation\n",GMT_program);
							for (k=0; k<6; k++) stat[k] = stat2[k];
							npts=nvalues-n_nan;
							decimated = FALSE;
						}
						tcrit = GMT_tcrit (0.975, (double)npts - 2.0);
						if (stat[5] != 1.0) {
							if (gmtdefs.verbose)
								fprintf (stderr, "%s: RLS regression insignificant\n",GMT_program);
						}
					}

					/* Analyze regression if significant */
					/* stat[0] is rls slope, stat[1] is rls icept, stat[2] is standard deviation, stat[3] */
					/* is sum of squares, stat[4] is r, and stat[5] == 1.0 if this is a significant correlation. */
					if (stat[5] == 1.0) {
						/* Check if ship data are incorrectly scaled */
						range = (tcrit * stat[2]) / sqrt(stat[3]);	/* Draper 1.4.8 */
						if (1.0 <= (stat[0]-range) || 1.0 >= (stat[0]+range)) {
							for (j = 0; j < 5; j++) {
								if (test_slope[j] >= (stat[0]-range) &&\
									test_slope[j] <= (stat[0]+range)) {
									if (display == "E77") {
										if (!GMT_is_dnan(adjustScale[this_grid[i].col]) && fabs(adjustScale[this_grid[i].col]-1.0)>0.0)
											fprintf (fpout, "Y-%s-E-%.02d: Regression scale %.3f"\
											" different from 1. Recommended: [%g]\n",this_grid[i].abbrev,\
											E77_HDR_SCALE,stat[0],adjustScale[this_grid[i].col]);
										else {
											fprintf (fpout, "N-%s-E-%.02d: Regression scale %.3f"\
											" different from 1. Recommended: [%g]\n",this_grid[i].abbrev,\
											E77_HDR_SCALE,stat[0],1.0/test_slope[j]);
										}
									}
									else if (warn[SUMMARY_WARN])
										fprintf (GMT_stdout, "%s (%s) Slope %.3f different from 1."\
										" Recommended: [%g]\n",argv[argno],this_grid[i].abbrev,stat[0],1/test_slope[j]);
#ifdef FIX
									/* Apply RLS slope to current field if new scale. */
									for (k = 0; k < nvalues; k++) {
										D[k].number[this_grid[i].col] /= test_slope[j];
										diff[i][k] = D[k].number[this_grid[i].col] - G[i][k]; /* Re-compute cruise - grid differences */
									}
									if (gmtdefs.verbose)
										fprintf (stderr, "%s (%s) Warning: Scaled by %g for internal along-track analysis\n",\
										argv[argno],this_grid[i].abbrev, test_slope[j]);
#endif
									newScale = TRUE;
									/* If not depth comparison skip fathom check */
									if (j == 1 && this_grid[i].col != MGD77_DEPTH) break;
								}
							}
							if (!newScale) {
								if (display == "E77") {
									if (!GMT_is_dnan(adjustScale[this_grid[i].col]) && fabs(adjustScale[this_grid[i].col]-1.0)>0.0)
										fprintf (fpout, "Y-%s-E-%.02d: Regression scale %.3f"\
										" different from 1. Recommended: [%g]\n",this_grid[i].abbrev,\
										E77_HDR_SCALE,stat[0],adjustScale[this_grid[i].col]);
									else {
										recommended_scale = (!strcmp (this_grid[i].abbrev, "faa")) ? pow (10.0, rint (log10 (1.0/stat[0]))) : ((stat[0] > 1.5) ? MGD77_METERS_PER_FATHOM : 1.0);
										fprintf (fpout, "N-%s-E-%.02d: Regression scale %.3f"\
										" different from 1. Recommended: [%g]\n",this_grid[i].abbrev,E77_HDR_SCALE,stat[0],recommended_scale);
									}
								}
								else if (warn[SUMMARY_WARN])
									fprintf (GMT_stdout, "%s (%s) Slope %.3f is statistically different from 1\n",\
									argv[argno],this_grid[i].abbrev,stat[0]);
							}
						}

						/* Check for significant DC shift */
						range = tcrit * stat[2] * sqrt(S_xx/(npts*stat[3]));	/* Draper 1.4.11 */
						if (this_grid[i].col != MGD77_DEPTH && (0.0 <= (stat[1]-range) || 0.0 >= (stat[1]+range))) {
							if (display == "E77") {
								if (!GMT_is_dnan(adjustDC[this_grid[i].col])) {
									if (adjustDC[this_grid[i].col] == 0)
										fprintf (fpout, "N");
									else
										fprintf (fpout, "Y");
									fprintf (fpout, "-%s-E-%.02d: Regression offset %.2f different from 0. Recommended: [%g]\n",\
									this_grid[i].abbrev,E77_HDR_OFFSET,stat[1],adjustDC[this_grid[i].col]);
								}
								else
									fprintf (fpout, "N-%s-E-%.02d: Regression offset %.2f"\
									" different from 0. Recommended: [%g]\n",this_grid[i].abbrev,E77_HDR_OFFSET,stat[1],\
									-0.1*rint (10.0*stat[1]));
							}
							else if (warn[GRID_WARN])
								fprintf (GMT_stdout, "%s (%s) Offset different than 0 (%.2f)\n",argv[argno],this_grid[i].abbrev,stat[1]);
#ifdef FIX
							/* Apply DC shift to current field */
							for (k = 0; k < nvalues; k++) {
								D[k].number[this_grid[i].col] -= stat[1];
								diff[i][k] = D[k].number[this_grid[i].col] - G[i][k]; /* Re-compute cruise - grid differences */
							}
							if (gmtdefs.verbose)
								fprintf (stderr, "%s (%s) Warning: Offset corrected by %g for internal along-track analysis\n",\
								argv[argno],this_grid[i].abbrev, stat[1]);
#endif
						}
					}
					if (warn[SUMMARY_WARN])
						fprintf (GMT_stdout, "%s (%s) RLS scale: %.3f offset: %.2f r: %.2f significant: %d decimation: %d\n",
						argv[argno],this_grid[i].abbrev,stat[0],stat[1],stat[4],(int)stat[5],(int)decimated);

					GMT_free ((void *)decimated_ship);
					GMT_free ((void *)decimated_grid);
				}
				else {
					/* Turn off this empty field */
					M.bit_pattern[0] = M.bit_pattern[0] & (0 << this_grid[i].col);

					fprintf (stderr, "%s (%s) Warning: Cruise contains insufficient observations (%d found)\n",\
					argv[argno],this_grid[i].abbrev, npts);
				}
				/* Free up regression array memory */
				GMT_free ((void *)ship_val);
				GMT_free ((void *)grid_val);
			}
		}

		if (display == "E77") {	/* Echo out the ranges to be flagged as bad */
			for (i = 0; i < (int)n_bad_sections; i++) {
				fprintf (fpout, "Y-%s-E-%.02d: Record range with invalid data: [%d-%d]\n", BadSection[i].abbrev,E77_HDR_FLAGRANGE,BadSection[i].start,BadSection[i].stop);
			}
		}

		/* CHECK IF CORRECT FAA REFERENCE MODEL */
		if (M.bit_pattern[0] & (1 << MGD77_GOBS) &&  M.bit_pattern[0] & (1 << MGD77_FAA)) {
			for (i = n = 0; i < nvalues; i++)
				if (!GMT_is_dnan(D[i].number[MGD77_GOBS]) && !GMT_is_dnan(D[i].number[MGD77_FAA])) n++;
			new_anom = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), "mgd77sniffer");
			old_anom = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), "mgd77sniffer");			
			for (i = j = 0; i < nvalues; i++) {
				if (!GMT_is_dnan(D[i].number[MGD77_GOBS]) && !GMT_is_dnan(D[i].number[MGD77_FAA])) {
					new_anom[j] = D[i].number[MGD77_GOBS] - MGD77_Theoretical_Gravity (D[i].number[MGD77_LONGITUDE], D[i].number[MGD77_LATITUDE], 4);					
					if (!GMT_is_dnan(D[i].number[MGD77_EOT]))
						new_anom[j] += D[i].number[MGD77_EOT];
					old_anom[j] = D[i].number[MGD77_FAA];
					j++;
				}
			}
			regress_rls (new_anom, old_anom, n, stat, MGD77_FAA, S_xx);
			tcrit = GMT_tcrit (0.975, (double)n - 2.0);
			/* Check for significant scale */
			range = (tcrit * stat[2]) / sqrt(stat[3]);	/* Draper 1.4.8 */
			range2 = tcrit * stat[2] * sqrt(S_xx/(n*stat[3]));	/* Draper 1.4.11 */
			if (stat[5] == 1.0 && ((1.0 <= (stat[0]-range) || 1.0 >= (stat[0]+range)) || (0.0 <= (stat[1]-range) || 0.0 >= (stat[1]+range)))) {
				if (display == "E77")
					fprintf (fpout, "I-faa-W-%.2d: anomaly differs from expected (m: %.3f b: %.3f)\n",E77_HDR_ANOM_FAA,stat[0],stat[1]);
				else if (warn[SUMMARY_WARN])
					fprintf (GMT_stdout, "%s (faa) anomaly differs from expected (m: %.3f b: %.3f)\n",argv[argno],stat[0],stat[1]);
			}
			GMT_free (new_anom);
			GMT_free(old_anom);
		}
		/* CHECK CORRECT MAG REFERENCE MODEL */
		if ((M.bit_pattern[0] & (1 << MGD77_MTF1) || M.bit_pattern[0] & (1 << MGD77_MTF2)) &&  M.bit_pattern[0] & (1 << MGD77_MAG)) {
			if (M.bit_pattern[0] & (1 << MGD77_MTF2))
				mtf1 = FALSE;
			for (i = n = 0; i < nvalues; i++)
				if (!GMT_is_dnan(D[i].number[MGD77_MTF2-(int)mtf1]) && !GMT_is_dnan(D[i].number[MGD77_MAG])) n++;
			new_anom = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), "mgd77sniffer");
			old_anom = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), "mgd77sniffer");			
			for (i = j = 0; i < nvalues; i++) {
				if (!GMT_is_dnan(D[i].number[MGD77_MTF1-(int)mtf1]) && !GMT_is_dnan(D[i].number[MGD77_MAG])) {
					GMT_gcal_from_dt (D[i].time, &cal);	/* No adjust for TZ; this is GMT UTC time */
					n_days = (GMT_is_gleap (cal.year)) ? 366.0 : 365.0;	/* Number of days in this year */
					/* Get date as decimal year */
					date = cal.year + cal.day_y / n_days + (cal.hour * GMT_HR2SEC_I + cal.min * GMT_MIN2SEC_I + cal.sec) * GMT_SEC2DAY;
					MGD77_igrf10syn (0, date, 1, 0.0, D[i].number[MGD77_LONGITUDE], D[i].number[MGD77_LATITUDE], IGRF);
					new_anom[j] = D[i].number[MGD77_MTF2-(int)mtf1] - IGRF[MGD77_IGRF_F];			
					if (!GMT_is_dnan(D[i].number[MGD77_DIUR]))
						new_anom[j] += D[i].number[MGD77_DIUR];
					old_anom[j] = D[i].number[MGD77_MAG];
					j++;
				}
			}
			regress_rls (new_anom, old_anom, n, stat, MGD77_MAG, S_xx);
			tcrit = GMT_tcrit (0.975, (double)n - 2.0);
			/* Check for significant scale */
			range = (tcrit * stat[2]) / sqrt(stat[3]);	/* Draper 1.4.8 */
			range2 = tcrit * stat[2] * sqrt(S_xx/(n*stat[3]));	/* Draper 1.4.11 */
			if (stat[5] == 1.0 && ((1.0 <= (stat[0]-range) || 1.0 >= (stat[0]+range)) || (0.0 <= (stat[1]-range) || 0.0 >= (stat[1]+range)))) {
				if (display == "E77")
					fprintf (fpout, "I-mag-W-%.2d: anomaly differs from expected (m: %.3f b: %.3f)\n",E77_HDR_ANOM_MAG,stat[0],stat[1]);
				else if (warn[SUMMARY_WARN])
					fprintf (GMT_stdout, "%s (mag) anomaly differs from expected (m: %.3f b: %.3f)\n",argv[argno],stat[0],stat[1]);
			}
			GMT_free (new_anom);
			GMT_free(old_anom);
		}		

		/* CHECK SANITY ALONG-TRACK */
		lastLat = D[0].number[MGD77_LATITUDE];
		lastLon = D[0].number[MGD77_LONGITUDE];
		for (curr = 0; curr < nvalues; curr++) {

			thisLat = D[curr].number[MGD77_LATITUDE];
			thisLon = D[curr].number[MGD77_LONGITUDE];

			/* Compute distance (keep units in km) */
			/* Use GMT great circle distance function to calculate arc length between */
			/* current and previous record (in degrees) then convert to km and accumulate */
			if (curr > 0) {
				lastLat = D[curr-1].number[MGD77_LATITUDE];
				lastLon = D[curr-1].number[MGD77_LONGITUDE];
			}
			ds = GMT_great_circle_dist (lastLon, lastLat, thisLon, thisLat) * MGD77_DEG_TO_KM;
			distance[curr] = ds + distance[curr-1];

			/* Create the current time string formatted according to gmtdefaults */
			if (gotTime)
				GMT_ascii_format_one (timeStr, D[curr].time, GMT_io.out_col_type[MGD77_TIME]);
			else
				GMT_ascii_format_one (timeStr, distance[curr], GMT_io.out_col_type[GMT_IS_FLOAT]);
				
			/* Create the location portion of the verbose data warning string (not for E77) */
			sprintf (placeStr,"%s %s %d",argv[argno],timeStr,curr+1);

			/* Check for time out of range */
			if (D[curr].time > maxTime || D[curr].time < \
			GMT_rdc2dt(GMT_rd_from_gymd(irint(mgd77snifferdefs[MGD77_YEAR].minValue), 1, 1),0.0)) {
				E[curr].flags[E77_NAV] |= NAV_TIME_OOR;
				if (warn[TIME_WARN])
					fprintf (GMT_stdout, "%s - Time out of range\n",placeStr);
			}

			nwords = nout = 0;
			/* Initialize output array for this record */
			out[curr] = (double *)GMT_memory (VNULL, n_out_columns, sizeof (double), GMT_program);
			for (i = 0; i < n_out_columns; i++)
				out[curr][i] = MGD77_NaN;


			/* Store latitude and longitude in the output array */
			if (display == "VALS") {
				out[curr][nout] = D[curr].number[MGD77_LATITUDE];
				nout++;
				out[curr][nout] = D[curr].number[MGD77_LONGITUDE];
				nout++;
			}

			/* SCAN FOR VALUES OUT OF RANGE (POINT-BY-POINT Error Checking) */
			/* Check the 24 numeric fields (start with latitude)*/
   			for (i = MGD77_RECTYPE; i < MGD77_N_NUMBER_FIELDS; i++) {

				/* Store cruise values in the output array */
				if ((MGD77_this_bit[i] & (MGD77_GEOPHYSICAL_BITS | MGD77_CORRECTION_BITS)) && display == "VALS") {
					out[curr][nout] = D[curr].number[i];
					nout++;
				}

				/* Only scan fields present in this cruise */
				if (M.bit_pattern[0] & (1 << i)) {
					switch (i) {
						case (MGD77_RECTYPE):
							if ((int) D[curr].number[i] != 3 && (int) D[curr].number[i] != 5) {
								E[curr].flags[E77_VALUE] |= (1 << i);
								if (warn[TYPE_WARN])
									fprintf (GMT_stdout, "%s - Invalid code %s [%d]\n",\
									placeStr,mgd77defs[i].abbrev, (int) D[curr].number[i]);
							}
							break;
						case (MGD77_PTC):
						case (MGD77_BTC):
							if (((int) D[curr].number[i] != 1 && (int) D[curr].number[i] != 3 &&\
							(int) D[curr].number[i] != 9)) {
								E[curr].flags[E77_VALUE] |= (1 << i);
								if (warn[TYPE_WARN])
									fprintf (GMT_stdout, "%s - Invalid code %s [%d]\n",\
									placeStr,mgd77defs[i].abbrev, (int) D[curr].number[i]);
							}
							break;
						case (MGD77_MSENS):
							if ((int)D[curr].number[i] != 1 && (int)D[curr].number[i] != 2 &&\
							(int) D[curr].number[i] != 9) {
								E[curr].flags[E77_VALUE] |= (1 << i);
								if (warn[TYPE_WARN])
									fprintf (GMT_stdout, "%s - Invalid code %s [%d]\n", \
									placeStr,mgd77defs[i].abbrev, (int) D[curr].number[i]);
							}
							break;
						case (MGD77_NQC):
							if ((int) D[curr].number[i] != 5 && (int) D[curr].number[i] != 6 &&\
							(int) D[curr].number[i] != 9) {
								E[curr].flags[E77_VALUE] |= (1 << i);
								if (warn[TYPE_WARN])
									fprintf (GMT_stdout, "%s - Invalid code %s [%d]\n",placeStr,\
									mgd77defs[i].abbrev, (int) D[curr].number[i]);
							}
							break;
						case (MGD77_BCC):
							bccCode = irint (D[curr].number[i]);
							if (M.bit_pattern[0] & MGD77_TWT_BIT || M.bit_pattern[0] & MGD77_DEPTH_BIT) {
								if (bccCode < 1 || bccCode > 55) {
									switch ((int) bccCode) {
										case 59:	/* Matthews', no zone */
										case 60:	/* S. Kuwahara Formula */
										case 61:	/* Wilson Formula */
										case 62:	/* Del Grosso Formula */
										case 63:	/* Carter's Tables */
										case 88:	/* Other */
										case 98:	/* Unknown */
										case 99:	/* Unspecified */
											break;
										default:
											E[curr].flags[E77_VALUE] |= (1 << i);
											if (warn[TYPE_WARN])
												fprintf (GMT_stdout, "%s - Invalid code %s [%d]\n",placeStr,\
												mgd77defs[i].abbrev, (int) bccCode);
											break;
									}
								}
							}
							break;
						case (MGD77_DAY):	/* Separate case since # of days in a month varies */
							if (!GMT_is_dnan (D[curr].number[i])) {
								if ((E[curr].flags[E77_VALUE] & (1 << MGD77_YEAR)) || (E[curr].flags[E77_VALUE] & (1 << MGD77_MONTH)))
									last_day = irint (mgd77snifferdefs[i].maxValue);	/* Year or month has error so we use 31 as last day in this month */
								else
									last_day = GMT_gmonth_length (irint(D[curr].number[MGD77_YEAR]), irint(D[curr].number[MGD77_MONTH]));			/* Number of day in the specified month */

								if (GMT_is_dnan (D[curr].number[i]) && (D[curr].number[i] < mgd77snifferdefs[i].minValue || D[curr].number[i] > last_day)) {
									E[curr].flags[E77_VALUE] |= (1 << i);
									if (warn[VALUE_WARN])
										fprintf (GMT_stdout, "%s - %s out of range [%f]\n", placeStr, mgd77defs[i].abbrev, D[curr].number[i]);
								}
							}
							break;
						case (MGD77_LONGITUDE): /* Handle longitudes in 180-360 range */
							if (!GMT_is_dnan (D[curr].number[i]) && D[curr].number[i] > mgd77snifferdefs[i].maxValue && D[curr].number[i] <= 360.0) {
								if (warn[VALUE_WARN])
									fprintf (GMT_stdout, "%s - %s adjusted %f to +/- 180\n", placeStr, mgd77defs[i].abbrev, D[curr].number[i]);
								D[curr].number[i] -= 360.0;
							}
						default:
							/* Verify that measurements are within range */
							if (!GMT_is_dnan (D[curr].number[i]) && (D[curr].number[i] < mgd77snifferdefs[i].minValue ||\
							D[curr].number[i] > mgd77snifferdefs[i].maxValue)) {
								E[curr].flags[E77_VALUE] |= (1 << i);
								if (warn[VALUE_WARN])
									fprintf (GMT_stdout, "%s - %s out of range [%f]\n", placeStr, mgd77defs[i].abbrev, D[curr].number[i]);
							}
							if ((i == MGD77_LATITUDE || i == MGD77_LONGITUDE) && GMT_is_dnan(D[curr].number[i])) {
								E[curr].flags[E77_NAV] |= (1 << i);
								if (warn[VALUE_WARN])
									fprintf (GMT_stdout, "%s - %s cannot be nine-filled\n", placeStr, mgd77defs[i].abbrev);
							}
							break;
					}
				}
			}

			/* Along-Track Excessive Slope Error Checking */
			if (curr > 0) {

				/* Check dt */
				if (!GMT_is_dnan(D[curr].time) && !GMT_is_dnan(D[curr-1].time))
					dt = D[curr].time - D[curr-1].time;
				else {	/* Time not specified */
					dt = MGD77_NaN;
					if (noTimeStart == -1)
						noTimeStart = curr;
					noTimeCount++;
				}
				if (!GMT_is_dnan(dt) && dt <= 0) { /* Non-increasing time */
					if (dt == 0) {
						if (warn[TIME_WARN])
							fprintf (GMT_stdout, "%s - Time not monotonically increasing (%f sec.)\n",placeStr, dt);
						dt = MGD77_NaN;
					} else {
						E[curr].flags[E77_NAV] |= NAV_TIME_DECR;
						if (warn[TIME_WARN])
							fprintf (GMT_stdout, "%s - Time decreasing (%f sec.)\n",placeStr, dt);
					}
					if (timeErrorStart == -1)
						timeErrorStart = curr;
					timeErrorCount++;
				}

				/* Calculate speed */
				if (dt != 0)
					speed = (ds*1000*distance_factor)/(dt*time_factor);
				else
					speed = MGD77_NaN;

#ifdef HISTOGRAM_MODE
				/* Use this to print excessive slopes to stderr for
				tracking winning cruises when running histogram scripts */
				if (!GMT_is_dnan(speed) && speed > max_speed) {
					E[curr].flags[E77_NAV] |= NAV_HISPD;
					fprintf (stderr, "%s - Excessive speed %f %s\n",placeStr, speed, speed_units);
				}
#endif

				/* Check that ship speed is reasonable */
				if (!GMT_is_dnan(speed)) {
					if (speed > max_speed) {
						E[curr].flags[E77_NAV] |= NAV_HISPD;
						if (warn[SPEED_WARN]) 
							fprintf (GMT_stdout, "%s - Excessive speed %f %s\n",placeStr, speed, speed_units);
					}
				}

				/* Store speed in the output array */
				nout = 0;
				if (display == "SLOPES") {
					if (!GMT_is_dnan(speed))
						out[curr][nout] = speed;
					else
						out[curr][nout] = MGD77_NaN;
					nout++;
				}

				/* Check slope values for non-time geophysical measurements */
				for (i = MGD77_LATITUDE; i < MGD77_N_NUMBER_FIELDS; i++) {

					/* Only scan floating point fields present in this cruise */						
					if (M.bit_pattern[0] & (1 << i) & MGD77_FLOAT_BITS) {

						/* Compute the difference between current and previous value */
						if (GMT_is_dnan(D[curr].number[i])) {
							gradient = MGD77_NaN;
							dvalue = MGD77_NaN;
						}
						else {
							/* Search backward to find a non-empty record for the same field. */
							/* Hope it doesn't have to search too far */
							for (j = 1; GMT_is_dnan(D[curr-j].number[i]) && curr-j > 0; j++)
								if (j > MGD77_MAX_SEARCH) break;
							dvalue = D[curr].number[i]-D[curr-j].number[i];
							lastLat = D[curr-j].number[MGD77_LATITUDE];
							lastLon = D[curr-j].number[MGD77_LONGITUDE];
							/* Calculate ds & dt between current and last valid record */
							/* Note: ds may be different for each field */
							ds = GMT_great_circle_dist (lastLon, lastLat, thisLon, thisLat) * MGD77_DEG_TO_KM;
							dt = D[curr].time - D[curr-j].time;

							/* Set to nan if a gap is detected (unless gap skipping is turned off) */
							if (ds > maxGap && maxGap != 0)
								dvalue = MGD77_NaN;

							/* Compute gradient */
							if (derivative == "SPACE") {
								if (ds >= MGD77_NAV_PRECISION_KM)
									gradient = dvalue / ds;
								else
									gradient = MGD77_NaN;
							}
							else if (derivative == "DIFF")
								gradient = dvalue;
							else { /* Gradient with respect to time */
								if (dt > 0)
									gradient = dvalue / dt;
								else
									gradient = MGD77_NaN;
							}

							/* First Derivative Sanity Check */
							if (fabs(gradient) > maxSlope[i]) {
								E[curr].flags[E77_SLOPE] |= (1 << i);
								if (warn[SLOPE_WARN])
									fprintf (GMT_stdout, "%s - excessive %s gradient %f\n", placeStr, mgd77defs[i].abbrev,\
									gradient);
							}


#ifdef HISTOGRAM_MODE
							/* Use this to print excessive slopes to stderr for tracking winning
							   cruises when running histogram scripts */
							if (fabs(gradient) > maxSlope[i]) {
								E[curr].flags[E77_SLOPE] |= (1 << i);
								fprintf (stderr, "%s - excessive %s gradient %f\n", placeStr, mgd77defs[i].abbrev, gradient);
							}
#endif
						}
					}
					else
						gradient = dvalue = ds = MGD77_NaN;

					if ((1 << i) & (MGD77_GEOPHYSICAL_BITS | MGD77_CORRECTION_BITS)) {
						if (display == "SLOPES")
							/* Store slopes in output array */
							out[curr][nout++] = gradient;
						if (display == "DFDS") {
							out[curr][nout++] = fabs(dvalue);
							if (GMT_is_dnan(dvalue))
								out[curr][nout++] = MGD77_NaN;
							else
								out[curr][nout++] = ds;
						}
					}
				}
			}

			/* CRUISE - GRID COMPARISON */
			if (n_grids > 0) {
				for (i = 0; i < n_grids; i++) {

					if (this_grid[i].g_pts < 2)
						continue;

					/* Fill output array with cruise/grid differences */
					if (display == "DIFFS") {
						out[curr][0] = D[curr].number[MGD77_LATITUDE];
						out[curr][1] = D[curr].number[MGD77_LONGITUDE];
						out[curr][2] = distance[curr]*distance_factor;
						out[curr][3+i*3] = D[curr].number[this_grid[i].col];
						out[curr][4+i*3] = G[i][curr];
						out[curr][5+i*3] = D[curr].number[this_grid[i].col]-G[i][curr];
					}

					/* Check if the cruise went over land */
					if (!strcmp(this_grid[i].abbrev,"depth") && G[i][curr]*this_grid[i].sign > 0) {
						E[curr].flags[E77_NAV] = NAV_ON_LAND;
						if (!landcruise)
							overLandStart = curr;
						landcruise = TRUE;
						overLandCount++;
					}

					if (M.bit_pattern[0] & (1 << this_grid[i].col)) {
						/* Track min/max absolute difference between grid and cruise */
						if (fabs(diff[i][curr]) > fabs(MaxDiff[i])) {
							MaxDiff[i] = diff[i][curr];
							iMaxDiff[i] = curr;
						}

						/* Compare cruise and grid data to find offsets */
						if (curr > 0) {

							/* Areas are estimated for each offset by summing discrete rectangles while the offset sign
							   stays the same */
							/* Search backward to find a non-empty record for the same field. */
							/* Hope it doesn't have to search too far */
							for (j = 1; GMT_is_dnan(D[curr-j].number[this_grid[i].col]) && curr-j > 0; j++)
								if (j > MGD77_MAX_SEARCH) break;
							lastLat = D[curr-j].number[MGD77_LATITUDE];
							lastLon = D[curr-j].number[MGD77_LONGITUDE];							
							/* Calculate ds & dt between current and last valid record */
							/* Note: ds may be different for each field */
							ds = GMT_great_circle_dist (lastLon, lastLat, thisLon, thisLat) * MGD77_DEG_TO_KM;
							dt = D[curr].time - D[curr-j].time;

							/* Calculate area of this offset */
							if (!GMT_is_dnan(diff[i][curr])) {
								thisArea = 0.5 * (diff[i][curr] + diff[i][curr-j]) * ds;
								prevOffsetSign[i] = offsetSign[i];
								offsetSign[i] = thisArea > 0;
							}
							else
								thisArea = MGD77_NaN;

							/* Accumulate area of total offset (Attempt to skip big gaps in data) */
							/* Think about different types of gaps to avoid */
							/* 1. Lose GPS - time stays consistent but ds gets huge - max ds test */
							/* 2. Ship stops logging and moves to another region then logs again - max dt test */
							/* 3. No recorded time - use max ds test */
							if (!GMT_is_dnan(thisArea) && offsetSign[i] == prevOffsetSign[i] && \
							   (ds < maxGap || maxGap == 0)) {
								offsetArea[i] += thisArea;
								offsetLength[i] += ds;
							}

							/* Analyze offset after it ends (sign switch) */
							if (offsetSign[i] != prevOffsetSign[i] ||  curr == nvalues - 1) {
								/* Flag the offset if significant */
								if (fabs(offsetArea[i]) > mgd77snifferdefs[this_grid[i].col].maxArea) {
									if (warn[GRID_WARN])
										fprintf (GMT_stdout, "%s - excessive offset from grid: Area/Length/Height"\
										" %f\t%f\t%f\t%d\n", placeStr, fabs(offsetArea[i]),offsetLength[i],\
										fabs(offsetArea[i]/offsetLength[i]),curr);
									/* Go back and flag all records in this offset */
									for (k = offsetStart[i]; k <= curr; k++)
										E[k].flags[E77_GRID] |= (1 << this_grid[i].col);
								}
								offsetArea[i] = 0;
								offsetLength[i] = 0;
								offsetStart[i] = curr;
							}
						}
					}
				}
			}

			/* Check for lower precision values and duplicate records */
			for (i = MGD77_LATITUDE; i < MGD77_N_NUMBER_FIELDS; i++) {

				/* Only check floating point fields present in this cruise */						
				if (M.bit_pattern[0] & (1 << i) & MGD77_FLOAT_BITS) {

					/* Compute the difference between current and previous value */
					if (!GMT_is_dnan(D[curr].number[i])) {

						/* Check for lower precision */
						for (j = 1; GMT_is_dnan(D[curr-j].number[i]) && curr-j > 0; j++)
							if (j > MGD77_MAX_SEARCH) break;
						dvalue = D[curr].number[i]-D[curr-j].number[i];

						if (!GMT_is_dnan(dvalue) && dvalue != 0) {
							if (fmod(dvalue,1.0) != 0.0 && (lowPrecision & (1 << i)))
								lowPrecision -= (1 << i);
							if (fmod(dvalue,5.0) != 0.0 && (lowPrecision5 & (1 << i)))
								lowPrecision5 -= (1 << i);
						}

						/* Check for duplicate records */
						if (i == MGD77_MSD) continue;
						if (dvalue == 0)
							duplicates[i]++;
						else {
							if (duplicates[i] > MGD77_MAX_DUPLICATES) {
								if (warn[VALUE_WARN])
									fprintf (GMT_stdout, "%s - %d duplicate %s records\n",placeStr,\
									duplicates[i], mgd77defs[i].abbrev);
								for (k = curr-1; k >= curr-duplicates[i]; k--)
									E[k].flags[E77_VALUE] |= (1 << i);
							}
							duplicates[i] = 0;
						}
					}
				}
			}
#ifndef FIX
			/* Dump data row-by-row if requested */
			if (display != "E77" && display != NULL) {
				if (display == "MGD77")
					MGD77_Write_Data_Record_m77 (&Out, &D[curr]);
				else {
					if (curr > 0 || display == "VALS") {
						if (GMT_io.binary[GMT_OUT])
							/* Use GMT output machinery which can handle binary output */
							GMT_output (GMT_stdout, n_out_columns, out);
						else {
							for (i = 0; i < (n_out_columns); i++) {
								GMT_ascii_output_one (GMT_stdout, out[curr][i], i);
								if ((i+1) < n_out_columns) fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
							}
							fprintf (GMT_stdout, "\n");
						}
					}
				}
			}
#endif
		}

		if (warn[TYPE_WARN]) {
			if (bccCode != 63 && bccCode != 88 && (M.bit_pattern[0] & MGD77_TWT_BIT))
				fprintf (GMT_stdout, "%s - new bathymetry correction table available\n",placeStr);
			else if (bccCode > 0 && bccCode < 56 && !(M.bit_pattern[0] & MGD77_TWT_BIT))
				fprintf (GMT_stdout, "%s - twt may be extracted from depth for Carter correction\n",placeStr);
		}

#ifdef FIX
		/* Turn off fields if errors found */		
		for (rec = 0; rec < curr; rec++) {
			deleteRecord = FALSE;
			for (type = 0; type < N_ERROR_CLASSES; type++) {
				if (E[rec].flags[type]) { /*Error in this category */
					thisLon = D[rec].number[MGD77_LONGITUDE];
					thisLat = D[rec].number[MGD77_LATITUDE];
					switch (type) {
						case E77_NAV:
							/* 9-fill records with nav errors */
							for (i=MGD77_PTC; i<MGD77_N_NUMBER_FIELDS; i++) D[rec].number[i] = MGD77_NaN;
							deleteRecord = TRUE;
							break;
						default:
							for (field = MGD77_PTC; field < n_types[type]; field++) {
								if (E[rec].flags[type] & (1 << field))
									D[rec].number[field] = MGD77_NaN;
							}
							break;
					}
					D[rec].number[MGD77_LONGITUDE] = thisLon;
					D[rec].number[MGD77_LATITUDE] = thisLat;
				}
			}
			
			/* Dump "fixed" data row-by-row if requested */
			if (display != "E77" && display != NULL && !deleteRecord) {
				if (display == "MGD77")
					MGD77_Write_Data_Record_m77 (&Out, &D[rec]);
				else {
					if (rec > 0 || display == "VALS") {
						if (GMT_io.binary[GMT_OUT])
							/* Use GMT output machinery which can handle binary output */
							GMT_output (GMT_stdout, n_out_columns, out);
						else {
							for (i = 0; i < (n_out_columns); i++) {
								GMT_ascii_output_one (GMT_stdout, out[rec][i], i);
								if ((i+1) < n_out_columns) fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
							}
							fprintf (GMT_stdout, "\n");
						}
					}
				}
			}
		}
#endif

		/* OUTPUT E77 ERROR FORMAT */
		if (display == "E77") {
			/* E77 HEADER MESSAGES */
			if ((bccCode != 63 && bccCode != 88 && (M.bit_pattern[0] & MGD77_TWT_BIT)) OR_TRUE)
				fprintf (fpout, "I-twt-W-%.2d: More recent bathymetry correction table available\n",E77_HDR_BCC);
			if ((bccCode > 0 && bccCode < 56 && !(M.bit_pattern[0] & MGD77_TWT_BIT)) OR_TRUE)
				fprintf (fpout, "I-twt-W-%.2d: twt may be extracted from depth for Carter correction\n",E77_HDR_BCC);

			/* Output data precision warnings */
			if (lowPrecision || lowPrecision5 OR_TRUE) {
				for (k = MGD77_LATITUDE; k < MGD77_N_NUMBER_FIELDS; k++) {
					if (lowPrecision5 & (1 << k) OR_TRUE)
						fprintf(fpout, "I-%s-W-%.02d: Integer multiple of 5 precision in %s\n",mgd77defs[k].abbrev,\
						E77_HDR_PRECISION,mgd77defs[k].abbrev);
					if (lowPrecision & (1 << k) OR_TRUE)
						fprintf(fpout, "I-%s-W-%.02d: Integer precision in %s\n",mgd77defs[k].abbrev,E77_HDR_PRECISION,
						mgd77defs[k].abbrev);
				}
			}

			/* E77 ERROR RECORDS */
			fprintf (fpout,"# Errata: Data\n");
			for (rec = 0; rec < curr; rec++) {
				if (gotTime)
					GMT_ascii_format_one (timeStr, D[rec].time, GMT_io.out_col_type[MGD77_TIME]);
				else
					GMT_ascii_format_one (timeStr, distance[rec], GMT_io.out_col_type[GMT_IS_FLOAT]);
				sprintf (placeStr, "%s%s%d%s",timeStr,gmtdefs.field_delimiter,rec+1,gmtdefs.field_delimiter);
				errorStr[0]='\0';
				for (type = 0; type < N_ERROR_CLASSES; type++) {
					if (E[rec].flags[type] OR_TRUE) { /*Error in this category */
						for (field = 0; field < (int)n_types[type]; field++) {
							if (E[rec].flags[type] & (1 << field) OR_TRUE)
								sprintf (errorStr, "%s%c",errorStr,'A'+field);
						}
					}
					else
						sprintf (errorStr, "%s0",errorStr);
					if (type < N_ERROR_CLASSES-1)
						sprintf (errorStr, "%s-",errorStr);
				}
				if (!strcmp(errorStr,"0-0-0-0")) continue;
				fprintf (fpout, "%s%s%s",placeStr,errorStr,gmtdefs.field_delimiter);
				prevType = FALSE;
				for (type = 0; type < N_ERROR_CLASSES; type++) {
					if (E[rec].flags[type] OR_TRUE) { /*Error in this category */
						fprintf (fpout, " ");
						if (prevType && (E[rec].flags[type] OR_TRUE))
							fprintf(fpout, "- ");
						if (type == E77_NAV)
							fprintf (fpout, "NAV: ");
						if (type == E77_VALUE)
							fprintf (fpout, "VAL: ");
						if (type == E77_SLOPE)
							fprintf (fpout, "GRAD: ");
						if (type == E77_GRID)
							fprintf (fpout, "GRID: ");
						for (field=0; field < (int)n_types[type]; field++) {
							if (E[rec].flags[type] & (1 << field) OR_TRUE) {
								if (prevFlag)
									fprintf (fpout, ", ");
								switch (type) {
									case E77_NAV:
										switch (1 << field) {
											case NAV_TIME_OOR:
												fprintf (fpout, "time out of range");
												break;
											case NAV_TIME_DECR:
												fprintf (fpout, "decreasing time");
												break;
											case NAV_HISPD:
												fprintf (fpout, "excessive speed");
												break;
											case NAV_ON_LAND:
												fprintf (fpout, "on land");
												break;
											case NAV_LAT_UNDEF:
												fprintf (fpout, "latitude undefined");
												break;
											case NAV_LON_UNDEF:
												fprintf (fpout, "longitude undefined");
												break;
											default:
												fprintf (fpout, "undefined nav error!");
												break;
										}
										break;
									default:
										fprintf (fpout, "%s",mgd77defs[field].abbrev);
										break;
								}
								prevFlag = TRUE;
							}
						}
						prevFlag = FALSE;
						switch (type) {
							case E77_VALUE:
								fprintf (fpout, " invalid");
								break;
							case E77_SLOPE:
								fprintf (fpout, " excessive");
								break;
							case E77_GRID:
								fprintf (fpout, " offset");
								break;
							default:
								break;
						}
						prevType = TRUE;					
					}
				}
				fprintf (fpout, "\n");
			}
		}

		/* Print Error Summary */
		if (warn[SUMMARY_WARN]) {
			if (noTimeCount == curr - 1)			/* no valid time in data */
				fprintf (GMT_stdout, "%s (time) Cruise contains no time record\n", argv[argno]);

			if (noTimeCount >0 && noTimeCount < curr)	/* print the first and number of nine-filled time records */
				fprintf (GMT_stdout, "%s (time) %d records contain no time starting at record #%d\n",argv[argno],\
				noTimeCount, noTimeStart);

			if (timeErrorCount > 0 && timeErrorCount < curr)/* print the first and number of time errors */
				fprintf (GMT_stdout, "%s (time) %d records had time errors starting at record #%d\n",argv[argno],\
				timeErrorCount, timeErrorStart);

			if (distanceErrorCount > 0) /* print the first and number of distance errors */
				fprintf (GMT_stdout, "%s (dist) %d records had distance errors starting at record #%d\n",argv[argno],\
				distanceErrorCount,distanceErrorStart);

			for (i = MGD77_LATITUDE; i < MGD77_N_NUMBER_FIELDS; i++) {
				if ((lowPrecision & (1 << i)) && !(lowPrecision5 & (1 << i)))
					fprintf (GMT_stdout, "%s (%s) Data Precision Warning: only integer values found\n",argv[argno],mgd77defs[i].abbrev);
				if (lowPrecision5 & (1 << i))	/* low precision data */
					fprintf (GMT_stdout, "%s (%s) Data Precision Warning: only integer multiples of 5 found\n",argv[argno],\
					mgd77defs[i].abbrev);
			}
			
			if (n_grids) {
				for (i = 0; i < n_grids; i++)
					fprintf (GMT_stdout, "%s (%s) Max ship-grid difference [%.1f] at record %d\n",argv[argno],\
					mgd77defs[this_grid[i].col].abbrev,MaxDiff[i],iMaxDiff[i]);
			}
			
			if (landcruise)	/* vessel went over land (GPS error or other gross navigation) */
				fprintf (GMT_stdout, "%s (nav) Navigation Warning: %d records went over land starting at record %d\n",\
				argv[argno], overLandCount, overLandStart);
		}
		/* Clean-up after finishing this cruise */
		MGD77_Close_File (&M);
		GMT_free ((void *)D);
		GMT_free ((void *)distance);
		GMT_free ((void *)H.mgd77);
		if (n_grids > 0) {
			for (i = 0; i < n_grids; i++) {
				GMT_free ((void *)G[i]);
				GMT_free ((void *)diff[i]);
			}
			GMT_free ((void *)G);
			GMT_free ((void *)diff);
		}
		if (display == "E77")
			fclose (fpout);
	}
	/* De-allocate grid memory */
	if (n_grids > 0) {
		GMT_free ((void *)MaxDiff);
		GMT_free ((void *)iMaxDiff);
		for (i = 0; i < n_grids; i++)
			GMT_free ((void *)f[i]);
	}
	exit (EXIT_SUCCESS);
}

void regress_rls (double *x, double *y, int nvalues, double *stat, int col, double S_xx)
{
	int i, n;
	double y_hat, threshold, s_0, res, *xx, *yy;
	
	regress_lms (x, y, nvalues, stat, col);
	/* Get LMS scale and use 2.5 of it to detect regression outliers */
	s_0 = 1.4826 * (1.0 + 5.0 / nvalues) * sqrt (stat[2]);
	threshold = 2.5 * s_0;
	
	xx = (double *)GMT_memory (VNULL, nvalues, sizeof (double), "mgd77sniffer");
	yy = (double *)GMT_memory (VNULL, nvalues, sizeof (double), "mgd77sniffer");
	for (i = n = 0; i < nvalues; i++) {
		y_hat = stat[0] * x[i] + stat[1];
		res = y[i] - y_hat;
		if (fabs (res) > threshold) continue;	/* Skip outliers */
		xx[n] = x[i];
		yy[n] = y[i];
		n++;
	}
	/* Now do LS regression of the 'good' points */
	regress_ls (xx, yy, n, stat, col, S_xx);
	/*stat[4] = GMT_corrcoeff (xx, yy, n, 0);*/
	if (n > 2) {	/* Determine if correlation is significant at 95% */
		double t, tcrit;
		t = stat[4] * sqrt (n - 2.0) / sqrt (1.0 - stat[4] * stat[4]);
		tcrit = GMT_tcrit (0.95, (double)n - 2.0);
		stat[5] = (double)(t > tcrit);	/* 1.0 if significant, 0.0 otherwise */
	}
	else
		stat[5] = GMT_d_NaN;

	GMT_free ((void *)xx);
	GMT_free ((void *)yy);
}

void regress_ls (double *x, double *y, int n, double *stat, int col, double S_xx)
{
	int i;
	double S, sum_x, sum_y, sum_x2, sum_y2, sum_xy;
	double mean_x, mean_y, S_xy, S_yy, y_discrepancy = 0.0;
	
	sum_x = sum_y = sum_x2 = sum_y2 = sum_xy = 0.0;
	mean_x = mean_y = S_xx = S_xy = S_yy = 0.0;
	S = (double)n;

	for (i = 0; i < n; i++) {
		sum_x += x[i];
		sum_y += y[i];
		sum_x2 += x[i] * x[i];
		sum_y2 += y[i] * y[i];
		sum_xy += x[i] * y[i];
	}
	
	mean_x = sum_x / S;
	mean_y = sum_y / S;

	S_xy = sum_xy - S * mean_x * mean_y;
	S_xx = sum_x2 - S * mean_x * mean_x;
	S_yy = sum_y2 - S * mean_y * mean_y;
	if (col != MGD77_DEPTH) { /* Use LMS m & b for depth (since offset forced to 0) */
		stat[0] = S_xy / S_xx;						/* Slope */
		stat[1] = mean_y - stat[0] * mean_x;	/* Intercept */
	}
	
	for (i = 0; i < n; i++)
		y_discrepancy += pow (y[i] - stat[0] * x[i] - stat[1], 2);
	stat[2] = sqrt (y_discrepancy / (S-1));			/* Standard deviation */
	stat[3] = S_xx;									/* Sum of squares */
	stat[4] = sqrt(S_xy * S_xy / (S_xx * S_yy));	/* Correlation (r) */
}

void regress_lms (double *x, double *y, int nvalues, double *stat, int col)
{

	double d_angle, limit, a, old_error, d_error, angle_0, angle_1;
	int n_angle;
	void regresslms_sub (double *x, double *y, double angle0, double angle1, int nvalues, int n_angle, double *stat, int col);

	d_angle = 1.0;
	limit = 0.1;
	n_angle = irint ((180.0 - 2 * d_angle) / d_angle) + 1;
	regresslms_sub (x, y, -90.0 + d_angle, 90.0 - d_angle, nvalues, n_angle, stat, col);
	old_error = stat[2];
	d_error = stat[2];

	while (d_error > limit) {
		d_angle = 0.1 * d_angle;
		a = atan (stat[0]) * 180 / M_PI;
		angle_0 = floor (a / d_angle) * d_angle - d_angle;
		angle_1 = angle_0 + 2.0 * d_angle;
		regresslms_sub (x, y, angle_0, angle_1, nvalues, 21, stat, col);
		d_error = fabs (stat[2] - old_error);
		old_error = stat[2];
	}
	if (gmtdefs.verbose)
		fprintf (stderr, "\n");
}

void regresslms_sub (double *x, double *y, double angle0, double angle1, int nvalues, int n_angle, double *stat, int col)
{
	double da, *slp, *icept, *z, *sq_misfit, *angle, *e, emin = DBL_MAX;
	int i, j = 0;

	slp = (double *) GMT_memory (VNULL, (size_t) n_angle, sizeof (double), "mgd77sniffer");
	icept = (double *) GMT_memory (VNULL, (size_t) n_angle, sizeof (double), "mgd77sniffer");
	angle = (double *) GMT_memory (VNULL, (size_t) n_angle, sizeof (double), "mgd77sniffer");
	e = (double *) GMT_memory (VNULL, (size_t) n_angle, sizeof (double), "mgd77sniffer");
	z = (double *) GMT_memory (VNULL, (size_t) nvalues, sizeof (double), "mgd77sniffer");
	sq_misfit = (double *) GMT_memory (VNULL, (size_t) nvalues, sizeof (double), "mgd77sniffer");

	for (i=0; i < 4; i++)
		stat[i] = 0;
	for (i=0; i < n_angle; i++)
		slp[i] = icept[i] = angle[i] = e[i] = 0;
	da = (angle1 - angle0) / (n_angle - 1);
	
	for (i = 0; i < n_angle; i++) {
		angle[i] = angle0 + i * da;
		if (gmtdefs.verbose) 
			fprintf (stderr, "%s: Trying angle...%+03.0f\r", GMT_program, angle[i]);
		slp[i] = tan (angle[i] * M_PI / 180.0);
		for (j = 0; j < nvalues; j++)
			z[j] = y[j] - slp[i] * x[j];
		if (col == MGD77_DEPTH)
			icept[i] = 0.0;
		else
			icept[i] = lms (z, nvalues);
		for (j = 0; j < nvalues; j++)
			sq_misfit[j] = pow(z[j]-icept[i], 2);
		e[i] = median (sq_misfit, nvalues);
	}
	for (i = 0; i < n_angle; i++) {
		if (e[i] < emin || i == 0) {
			emin = e[i];
			j = i;
		}
	}
	stat[0] = slp[j];
	stat[1] = icept[j];
	stat[2] = e[j];

	GMT_free ((void *)slp);
	GMT_free ((void *)icept);
	GMT_free ((void *)angle);
	GMT_free ((void *)e);
	GMT_free ((void *)z);
	GMT_free ((void *)sq_misfit);
}

double lms (double *x, int n)
{
	double mode;
	int GMT_n_multiples = 0;

	GMT_mode (x, n, n/2, 1, 0, &GMT_n_multiples, &mode);
	return mode;
}

double median (double *x, int n)
{
	double *sorted, med;
	int i;
	
	sorted = (double *) GMT_memory (VNULL, (size_t) n, sizeof (double), "mgd77sniffer");
	for (i=0; i < n; i++)
		sorted[i] = x[i];
	qsort ((void *) sorted, n, sizeof(double), GMT_comp_double_asc);
	med = (i%2) ? sorted[i/2] : 0.5*(sorted[(i-1)/2]+sorted[i/2]);
	GMT_free ((void *)sorted);
	return med;
}

/* Read Grid Header (from Smith & Wessel grdtrack.c) */
void read_grid (struct MGD77_GRID_INFO *info, float **grid, double w, double e, double s, double n, BOOLEAN bilinear, double threshold) {

	if (strlen(info->fname) == 0) return;	/* No name */
	
	GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 2;  /* GMT_pad is declared automatically in gmt.h */
	
	if (info->format == 0) {	/* GMT geographic grid with header */
		if (GMT_read_grd_info (info->fname, &info->grdhdr)) {
				fprintf (stderr, "%s: Error opening file %s\n", GMT_program, info->fname);
				exit (EXIT_FAILURE);
		}

		/* Get grid dimensions */
		info->one_or_zero = (info->grdhdr.node_offset) ? 0 : 1;
		info->nx = irint ( (info->grdhdr.x_max - info->grdhdr.x_min) / info->grdhdr.x_inc) + info->one_or_zero;
		info->ny = irint ( (info->grdhdr.y_max - info->grdhdr.y_min) / info->grdhdr.y_inc) + info->one_or_zero;

		/* Allocate grid memory */
		*grid = (float *) GMT_memory (VNULL, (size_t)((info->nx + 4) * (info->ny + 4)), sizeof (float), GMT_program);
		
		/* Read into memory with 2 rows/cols as padding */
		if (GMT_read_grd (info->fname, &info->grdhdr, *grid, w, e, s, n, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error reading file %s\n", GMT_program, info->fname);
			exit (EXIT_FAILURE);
		}
	}
	else {	/* Read a Mercator grid Sandwell/Smith style */
		GMT_read_img (info->fname, &info->grdhdr, grid, w, e, s, n, info->scale, info->mode, info->max_lat, TRUE);
		if (fabs (info->grdhdr.x_max - info->grdhdr.x_min - 360.0) < GMT_CONV_LIMIT) GMT_boundcond_parse (&info->edgeinfo, "g");
	}

	info->mx = info->grdhdr.nx + 4;
	info->interpolate = (threshold > 0.0);
	if (fabs (info->grdhdr.x_max - info->grdhdr.x_min - 360.0) < GMT_CONV_LIMIT) GMT_boundcond_parse (&info->edgeinfo, "g");
	
	GMT_boundcond_param_prep (&info->grdhdr, &info->edgeinfo);

	/* Initialize bcr structure with 2 row/col boundaries:  */

	GMT_bcr_init (&info->grdhdr, GMT_pad, bilinear, threshold, &info->bcr);
	
	/* Set boundary conditions  */
	GMT_boundcond_set (&info->grdhdr, &info->edgeinfo, GMT_pad, *grid);
}

/* Sample Grid at Cruise Locations (from Smith & Wessel grdtrack.c) */
int sample_grid (struct MGD77_GRID_INFO *info, struct MGD77_DATA_RECORD *D, double **g, float *grid, int n_grid, int n) {

	int rec, pts = 0, ii, jj;
	double MGD77_NaN, x, y;
	GMT_make_dnan (MGD77_NaN);

	/* Get grid values at cruise locations */
	for (rec = 0; rec < n; rec++) {

		if (info->format == 1)	/* Mercator IMG grid - get Mercator coordinates x,y */
			GMT_geo_to_xy (D[rec].number[MGD77_LONGITUDE], D[rec].number[MGD77_LATITUDE], &x, &y);
		else {		/* Regular geographic grd, just copy lon,lat to x,y */
			x = D[rec].number[MGD77_LONGITUDE];
			y = D[rec].number[MGD77_LATITUDE];
			/* Adjust cruise longitude if necessary; We know cruise is between +/-180 */ 
			if (info->grdhdr.x_min >= 0.0 && D[rec].number[MGD77_LONGITUDE] < 0.0)
				GMT_lon_range_adjust(0, &x);	/* Adjust to 0-360 range */
		}
		g[n_grid][rec] = MGD77_NaN;	/* Default value if we are outside the grid domain */
		if (y < info->grdhdr.y_min || y > info->grdhdr.y_max) continue;	/* Outside latitude range */

		/* If point is outside grd area, shift it using periodicity or skip if not periodic. */
		while ( (x < info->grdhdr.x_min) && (info->edgeinfo.nxp > 0))
			x += (info->grdhdr.x_inc * info->edgeinfo.nxp);
		if (x < info->grdhdr.x_min)
			continue;  /* West of our area */

		while ((x > info->grdhdr.x_max) && (info->edgeinfo.nxp > 0)) 
			x -= (info->grdhdr.x_inc * info->edgeinfo.nxp);
		if (x > info->grdhdr.x_max)
			continue;  /* East of our area */

		/* Get the value from the grid - it could be NaN */
		if (info->interpolate) {	/* IMG has been corrected, and GRD is good to go */
			g[n_grid][rec] = GMT_get_bcr_z (&info->grdhdr, x, y, grid, &info->edgeinfo, &info->bcr);
		}
		else {	/* Take IMG nearest node and do special stuff (values already set during read) */
			ii = GMT_x_to_i (x, info->grdhdr.x_min, info->grdhdr.x_inc, info->grdhdr.xy_off, info->grdhdr.nx);
			jj = GMT_y_to_j (y, info->grdhdr.y_min, info->grdhdr.y_inc, info->grdhdr.xy_off, info->grdhdr.ny);
			g[n_grid][rec] = grid[(jj+GMT_pad[3])*info->mx+ii+GMT_pad[0]];
		}
		pts++;
	}
	return pts;
}
