/*--------------------------------------------------------------------
 *	$Id: mgd77list.c,v 1.3 2004/08/19 21:29:22 pwessel Exp $
 *
 *    Copyright (c) 2004 by P. Wessel
 *    See README file for copying and redistribution conditions.
 *--------------------------------------------------------------------*/
/*
 * mgd77list produces ASCII listings of <legid>.mgd77 files. The *.mgd77
 * files distributed from NGDC contain along-track geophysical observations
 * such as bathymetry, gravity, and magnetics, and the user may extract
 * any combination of these parameters as well as four generated quantities
 * such as distance (in km), heading, velocity (m/s), and weight by using
 * the -F option.  The order of the choices given to the -F option is used
 * to determine the sequence in which the parameters will be printed out.
 * If -F is not specified, the default will output all the data.  E.g. to
 * create an input file for surface, use -Flon,lat,depth (for bathymetry).
 *
 * To select a sub-section of the track, specify the start/endpoints by:
 *	1) Start-time (yyyy-mm-ddT[hh:mm:ss]) OR start-distance (km)
 *	2) Stop-time  (yyyy-mm-ddT[hh:mm:ss]) OR stop-distance (km)
 * To select data inside an area, use the -R option.
 * To start output with a header string, use -H.
 * To separate each data set with a multisegmnt header string, use -M.
 *
 * Author:	Paul Wessel
 * Date:	19-JUN-2004
 * Version:	1.0 Based somewhat on the old gmtlist.c
 *
 *
 */
 
#include "mgd77.h"

main (int argc, char **argv)
{
	int i, rec, n_out = 0, argno, n_cruises = 0;
	
	unsigned int this_pattern;
	
	BOOLEAN error = FALSE, exact = FALSE, nautical = FALSE, string_output = FALSE;
	
	char *start_date = CNULL, *stop_date = CNULL, file[BUFSIZ];

	double west, east, south, north, weight, start_time, stop_time, start_dist;
	double stop_dist, dist_scale, vel_scale, this_dist, this_time, this_lon, this_speed;
	double this_lat, last_time, last_lon, last_lat, dx, dy, dlon, ds, dt, this_heading, *out;
	
	struct MGD77_CONTROL M;
	struct MGD77_HEADER_RECORD H;
	struct MGD77_DATA_RECORD D;
	
	FILE *fp;
	
	argc = GMT_begin (argc, argv);		/* Initialize GMT Machinery */
	
	/* Initialize MGD77 output order and other parameters*/
	
	MGD77_Init (&M, TRUE);			/* Initialize MGD77 Machinery */
	west = 0.0;	east = 360.0;	south = -90.0;	north = 90.0;
	start_time = start_dist = 0.0;
	stop_time = stop_dist = DBL_MAX;
	weight = 1.0;

	for (i =1; !error && i < argc; i++) {	/* Process input options */
		if (argv[i][0] == '-') {
			switch(argv[i][1]) {
			
				case 'H':
				case 'R':
				case 'V':
				case 'b':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;
										
				case 'C':	/* Crossover correction table */
					/* NOT IMPLEMENTED YET */
					break;
					
				case 'D':		/* Assign start/stop times for sub-section */
					if (argv[i][2] == 'a') {		/* Start date */
						start_date = &argv[i][3];
					}
					else if (argv[i][2] == 'b')	 {	/* Stop date */
						stop_date = &argv[i][3];
					}
					else
						error = TRUE;
					break;

				case 'N':	/* Nautical units (knots, nautical miles) */
					nautical = TRUE;
					break;
					
				case 'E':	/* Exact parameter match */
					exact = TRUE;
					break;
					
				case 'F':	/* Selected output fields */
					MGD77_Select_Columns (&argv[i][2], &M);
					break;
						
				case 'M':               /* Multiple line segments */
					GMT_multisegment (&argv[i][2]);
					break;

				case 'S':		/* Assign start/stop position for sub-section */
					if (argv[i][2] == 'a')		/* Start position */
						start_dist = atof (&argv[i][3]);
					else if (argv[i][2] == 'b')	/* Stop position */
						stop_dist = atof (&argv[i][3]);
					else
						error = TRUE;
					break;
					
				case 'W':		/* Assign a weight to these data */
					weight = atof (&argv[i][2]);
					break;
					
				default:		/* Options not recognized */
					error = TRUE;
					break;
			}
		}
		else
			n_cruises++;
	}
	
	/* Check that the options selected are mutually consistent */
	
	string_output = ((M.bit_pattern & MGD77_STRING_BITS) > 0);	/* TRUE if we request one of more strings */
	if (string_output && GMT_io.binary[1]) {
		fprintf(stderr, "%s: ERROR: Cannot specify binary output AND the text items ID, sln, or sspn\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (start_date && start_dist > 0.0) {
		fprintf(stderr, "%s: ERROR: Cannot specify both start time AND start distance\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (stop_date && stop_dist < DBL_MAX) {
		fprintf(stderr, "%s: ERROR: Cannot specify both stop time AND stop distance\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 (weight <= 0.0) {
		fprintf(stderr, "%s: ERROR: -W weight must be nonzero\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (start_date) GMT_atoft (start_date, &start_time);
	if (stop_date)  GMT_atoft (stop_date,  &stop_time);
	if (start_time > stop_time) {
		fprintf(stderr, "%s: ERROR: Start time exceeds stop time!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (start_dist > stop_dist) {
		fprintf(stderr, "%s: ERROR: Start distance exceeds stop distance!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	if (GMT_quick || argc == 1) {	/* Display usage */
		fprintf(stderr,"usage: mgd77list <cruise(s)> [-Da<startdate>] [-Db<stopdate>] [-E] [-F<dataflags>] [-H] [-M[<flag>]] [-N]\n");
		fprintf(stderr,"	[-R<west>/<east>/<south>/<north>] [-Sa<startdist>] [-Sb<stopdist>] [-V] [-W<Weight>] [-bo]\n\n");
         
		if (GMT_quick) exit (EXIT_FAILURE);
              
		fprintf(stderr,"	<cruises> is one or more MGD77 legnames, e.g., 01010083 etc.\n");
		fprintf(stderr,"	OPTIONS:\n\n");
/*		fprintf(stderr,"	-C<file> applies crossover corrections to the data. If no file name\n");
		fprintf(stderr,"	   is given, the default correction file is assumed\n"); */
		fprintf(stderr,"	-Da<date> lists from date (given as yyyy-mm-ddT[hh:mm:ss]) [Start of cruise]\n");
		fprintf(stderr,"	-Db<date> lists up to date (given as yyyy-mm-ddT[hh:mm:ss]) [End of cruise]\n");
		fprintf(stderr,"	-E Only records that exactly matches the requested geophysical information in -F will be used.\n");
		fprintf(stderr,"	   [Default will output all record that matches at least one column]\n");
		fprintf(stderr,"	-F Dataflags is a string made up of one or more of these abbreviations:\n");
		fprintf(stderr,"	   >Track information:\n");
		fprintf(stderr,"	     time:    Absolute time (formatted according to OUTPUT_DATE_FORMAT, OUTPUT_CLOCK_FORMAT)\n");
		fprintf(stderr,"	     lon:     Longitude (formatted according to OUTPUT_DEGREE_FORMAT)\n");
		fprintf(stderr,"	     lat:     Latitude (formatted according to OUTPUT_DEGREE_FORMAT)\n");
		fprintf(stderr,"	     ID:      Cruise ID [TEXTSTRING]\n");
		fprintf(stderr,"	     dist:    Along-track distance (km)\n");
		fprintf(stderr,"	     azim:    Track azimuth (Degrees east from north)\n");
		fprintf(stderr,"	     vel:     Ship velocity (m/s)\n");
		fprintf(stderr,"	   >Geophysical Observations:\n");
		fprintf(stderr,"	     twt:    Two-way traveltime (s)\n");
		fprintf(stderr,"	     depth:  Corrected bathymetry (m)\n");
		fprintf(stderr,"	     mtf1:   Magnetic Total Field Sensor 1 (nTesla)\n");
		fprintf(stderr,"	     mtf2:   Magnetic Total Field Sensor 2 (nTesla)\n");
		fprintf(stderr,"	     mag:    Magnetic residual anomaly (nTesla)\n");
		fprintf(stderr,"	     gobs:   Observed gravity (mGal)\n");
		fprintf(stderr,"	     faa:    Free-air gravity anomaly (mGal)\n");
		fprintf(stderr,"	   >Codes, Corrections, and Information:\n");
		fprintf(stderr,"	     ptc:    Position type code\n");
		fprintf(stderr,"	     bcc:    Bathymetric correction code\n");
		fprintf(stderr,"	     btc:    Bathymetric type code\n");
		fprintf(stderr,"	     msens:  Magnetic sensor for residual field\n");
		fprintf(stderr,"	     msd:    Magnetic sensor depth/altitude (m)\n");
		fprintf(stderr,"	     diur:   Magnetic diurnal correction (nTesla)\n");
		fprintf(stderr,"	     eot:    Eotvos correction (mGal)\n");
		fprintf(stderr,"	     sln:    Seismic line number string [TEXTSTRING]\n");
		fprintf(stderr,"	     sspn:   Seismic shot point number string [TEXTSTRING]\n");
		fprintf(stderr,"	     weight: Give weight specified in -W\n");
		fprintf(stderr,"	     nqc:    Navigation quality code\n");
		fprintf(stderr,"	  The data are written in the order specified in <dataflags>\n");
		fprintf(stderr,"	  [Default is all data types and all records]\n");
		fprintf(stderr,"	  Use UPPER CASE <dataflag> to only request records where this col is present.\n");
		fprintf(stderr,"	  (Note that -E is a shorthand to set all geophysical observations to upper case).\n");
		fprintf(stderr,"	  Optionally, append comma-separated tests that columns must pass to be output.\n");
		fprintf(stderr,"	  Format is <flag><OP><value>, where flag is any of the dataflags above, and <OP> is\n");
		fprintf(stderr,"	  one of the operators <, <=, =, >=, >, and !=.  <value> is the limit you are testing,\n");
		fprintf(stderr,"	  including NaN (with = and != only).  If <flag> is UPPERCASE the test MUST be passed;\n");
		fprintf(stderr,"	  else at least one of the tests must pass for output to take place.\n");
		fprintf(stderr,"	-H write one header record\n");
		fprintf(stderr,"	-M write multisegment header records for each cruise\n");
		fprintf(stderr,"	-N Use Nautical units (knots, nautical miles [Default is m/s and km]\n");
		fprintf(stderr,"	-R only return data inside the specified region [0/360/-90/90]\n");
		fprintf(stderr,"	-Sa<dist> lists from dist (in km or nm (-N)) [Start of the cruise]\n");
		fprintf(stderr,"	-Sb<dist> lists up to dist (in km or nm (-N)) [End of the cruise]\n");
		fprintf(stderr,"	-V verbose, report progress\n");
		fprintf(stderr,"	-W sets weight for these data [1]\n");
		GMT_explain_option ('o');
		exit(EXIT_FAILURE);
	}

	if (n_cruises == 0) {
		fprintf(stderr, "%s: ERROR: No cruises given\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	dist_scale = (nautical) ? 1000.0 * DEG2KM / MGD77_METERS_PER_NM : DEG2KM;
	vel_scale = (nautical) ? 3600.0 : 1000.0;	/* Since ds is in nm or km and dt is in sec */
	if (!string_output) out = (double *) GMT_memory (VNULL, (size_t)M.n_out_columns, sizeof(double), GMT_program);

	for (i = 0; i < M.n_out_columns; i++) {	/* Prepare GMT output formatting machinery */
		switch (M.order[i]) {
			case MGD77_TIME:	/* Special time formatting */
				GMT_io.out_col_type[i] = GMT_IS_ABSTIME;
				break;
			case MGD77_LONGITUDE:	/* Special lon formatting */
				GMT_io.out_col_type[i] = GMT_IS_LON;
				break;
			case MGD77_LATITUDE:	/* Special lat formatting */
				GMT_io.out_col_type[i] = GMT_IS_LAT;
				break;
			default:		/* Everything else is float (not true for the 3 strings though) */
				GMT_io.out_col_type[i] = GMT_IS_FLOAT;
				break;
		}
	}
		
	if (!GMT_io.binary[GMT_OUT] && gmtdefs.io_header[GMT_OUT]) {	/* Write out header record */
		fprintf (GMT_stdout, "# ");
		for (i = 0; i < M.n_out_columns; i++) {
			switch (M.order[i]) {
				case MGD77_TIME:	/* Print out time header */
					fprintf (GMT_stdout, "time");
					break;
				case MGD77_DISTANCE:	/* Print out distance */
					(nautical) ? fprintf (GMT_stdout, "dist(nm)") : fprintf (GMT_stdout, "dist(km)");
					break;
				case MGD77_HEADING:	/* Print out azimuth */
					printf ("azim");
					break;
				case MGD77_SPEED:	/* Print out velocity */
					(nautical) ? fprintf (GMT_stdout, "speed(kts)") : fprintf (GMT_stdout, "speed(m/s)");
					break;
				case MGD77_WEIGHT:	/* Print out weights */
					fprintf (GMT_stdout, "weight");
					break;
				default:	/* Print out named items */
					fprintf (GMT_stdout, "%s", mgd77defs[M.order[i]].abbrev);
					break;
			}
			if ((i+1) < M.n_out_columns) fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
		}
		fprintf (GMT_stdout, "\n");
	}

	for (argno = 1; argno < argc; argno++) {	/* Loop over all the files */
	
		if (argv[argno][0] == '-') continue;
  		if (MGD77_Get_Path (file, argv[argno], &M)) {
   			fprintf (stderr, "%s : Cannot find leg %s\n", argv[argno], GMT_program);
     			continue;
  		}
		if ((fp = GMT_fopen (file, "r")) == NULL) {
			fprintf (stderr, "%s: Could not open %s\n", file, GMT_program);
			continue;
		}

		if (gmtdefs.verbose) fprintf (stderr, "%s: Now processing cruise %s\n", GMT_program, argv[argno]);
		
		if (!MGD77_Read_Header_Record (fp, &H)) {
			fprintf (stderr, "%s: Error reading header sequence for cruise %s\n", GMT_program, argv[argno]);
			exit (EXIT_FAILURE);
		}
		if (GMT_io.multi_segments) {	/* Write multisegment header between each cruise */
			sprintf (GMT_io.segment_header, "%c %s\n", GMT_io.EOF_flag, argv[argno]);
			GMT_write_segmentheader (GMT_stdout, M.n_out_columns);
		}
		this_dist = this_lon = this_lat = this_time = 0.0;
		rec = 0;
	
		/* Start reading data from file */
	
		while (MGD77_Read_Data_Record (fp, &D)) {		/* While able to read a data record */
		
			/* Compute accumulated distance along track (Flat Earth) */
		
			last_lon  = this_lon;
			last_lat  = this_lat;
			last_time = this_time;
			this_lon  = D.number[MGD77_LONGITUDE];
			this_lat  = D.number[MGD77_LATITUDE];
			this_time = D.time;
			if (rec == 0) {	/* Initialize various counters */
				ds = dt = 0.0;
				this_heading = this_speed = GMT_d_NaN;
			}
			else {		/* Need a previous point to calculate distance, speed, and heading */
				dlon = this_lon - last_lon;
				if (fabs (dlon) > 180.0) dlon = copysign ((360.0 - fabs (dlon)), dlon);
				dx = dlon * cosd (0.5 * (this_lat + last_lat));
				dy = this_lat - last_lat;
				ds = dist_scale * hypot (dx, dy);
				dt = this_time - last_time;
				if (M.use_column[MGD77_HEADING]) {
					this_heading = (dx == 0.0 && dy == 0.0) ? GMT_d_NaN : 90.0 - R2D * atan2 (dy, dx);
					if (this_heading < 0.0) this_heading += 360.0;
				}
				if (M.use_column[MGD77_SPEED]) this_speed = (dt == 0.0) ? GMT_d_NaN : vel_scale * ds / dt;
				this_dist += ds;
			}
			rec++;
			
			/* Check if record has the required fields */
			
			if (exact) {	/* Want to match the geophysical fields exactly */
				this_pattern = M.bit_pattern & MGD77_GEOPHYSICAL_BITS;		/* Knock off the bits for distance, time etc since not part of the data */
				if ((D.bit_pattern & this_pattern) != this_pattern) continue;
			}
			
			/* Check if time or this_dist falls outside specified range */
		
			if (this_dist < start_dist || this_dist > stop_dist) continue;
			if (this_time < start_time || this_time > stop_time) continue;
		
			/* Check is lat/lon is outside specified area */
		
			if (this_lat < south || this_lat > north) continue;
			while (this_lon > east) this_lon -= 360.0;
			while (this_lon < west) this_lon += 360.0;
			if (this_lon > east) continue;
		
			/* Check if it passes any given column data constraints */
			
			if (!MGD77_pass_record (&D, &M)) continue;	/* Failed the test */

			/* This record will now be printed out */
		
			if (string_output) {	/* Must do it col by col and deal with the requested string(s) */
				for (i = 0; i < M.n_out_columns; i++) {
					switch (M.order[i]) {
						case MGD77_TIME:	/* Print out time header */
							GMT_ascii_output_one (GMT_stdout, this_time, i);
							break;
						case MGD77_DISTANCE:	/* Print out distance */
							GMT_ascii_output_one (GMT_stdout, this_dist, i);
							break;
						case MGD77_HEADING:	/* Print out azimuth */
							GMT_ascii_output_one (GMT_stdout, this_heading, i);
							break;
						case MGD77_SPEED:	/* Print out velocity */
							GMT_ascii_output_one (GMT_stdout, this_speed, i);
							break;
						case MGD77_WEIGHT:	/* Print out weights */
							GMT_ascii_output_one (GMT_stdout, weight, i);
							break;
						case MGD77_ID:
						case MGD77_SLN:
						case MGD77_SSPN:
							fprintf (GMT_stdout, "%s", D.word[M.order[i]-MGD77_ID]);
							break;
						default:	/* Print out named floating point items */
							GMT_ascii_output_one (GMT_stdout, D.number[M.order[i]], i);
							break;
					}
					if ((i+1) < M.n_out_columns) fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
				}
				fprintf (GMT_stdout, "\n");
			}
			else {	/* Use GMT output machinery which can handle binary output, if requested */
				for (i = 0; i < M.n_out_columns; i++) {
					switch (M.order[i]) {
						case MGD77_TIME:	/* Print out time header */
							out[i] = this_time;
							break;
						case MGD77_DISTANCE:	/* Print out distance */
							out[i] = this_dist;
							break;
						case MGD77_HEADING:	/* Print out azimuth */
							out[i] = this_heading;
							break;
						case MGD77_SPEED:	/* Print out velocity */
							out[i] = this_speed;
							break;
						case MGD77_WEIGHT:	/* Print out weights */
							out[i] = weight;
							break;
						default:	/* Print out named floating point items */
							out[i] = D.number[M.order[i]];
							break;
					}
				}
				GMT_output (GMT_stdout, M.n_out_columns, out);
			}
			n_out++;
		}
		GMT_fclose (fp);
	}
	
	if (!string_output) GMT_free ((void *)out);
	
	if (gmtdefs.verbose) fprintf (stderr, "%s: Returned %d output records\n", GMT_program, n_out);
	
	exit (EXIT_SUCCESS);
}
