#include <stdio.h>
#include <stdlib.h>		/* calloc(), free() */
#include <string.h>		/* memset() */
#include <math.h>		/* fabs(), rint() */

#ifndef BATH_UTIL
#include "bath-util.h"
#endif

/*
 * Peter Daly
 * Ocean Acoustics Group
 * Massachusetts Institute of Technology
 * Room 5-435
 * Cambridge, Massachusetts, USA
 *
 * This program extracts bathymetric information between the FFP array and
 * any specified point.
 *
 * Copyright (C) 1998, Massachusetts Institute of Technology
 *                     All Rights Reserved
 */
int
main (int argc, char **argv)
{
  int arg_counter;
  char *bathymetric_filename;
  Bathymetry *b;
  FILE *ofp = stdout, *bfp;

  UTM utms, utmd, step;
  LL lls, lld;

  double sample_interval, bearing;
  double r, range, depth;

  unsigned char range_independent_flag = 0;

  /* Initialize coordinates */

  memset(&(utms), 0, sizeof(UTM));
  memset(&(utmd), 0, sizeof(UTM));
  memset(&(lls), 0, sizeof(LL));
  memset(&(lld), 0, sizeof(LL));

  /* Parse the command line */
  
  if (argc == 1) {
    fprintf (stderr, "usage: %s bathymetry-filename [-ri] [ -d sample_interval ] [ -su utmx utmy ] [ -sl lat long ] [ -du umtmx utmy ] [ -dl lat log]\n",
	     argv[0]);
    return(-1);
  }

  sample_interval = 500.0;	/* meters, default */
  arg_counter = 0;

  if (strcmp(argv[++arg_counter], "-ri") == 0) 
      range_independent_flag = 1;
  else
    arg_counter--;

  bathymetric_filename = argv[++arg_counter];

  while (++arg_counter < argc) {

    if (strcmp(argv[arg_counter], "-d") == 0) {
      if (arg_counter < argc) {
	if (sscanf(argv[++arg_counter], "%lf", &sample_interval) != 1) {
	  fprintf (stderr, "sample interval should be a floating point number, in meters\n");
	  fprintf (stderr, "not [%s]\n",
		   argv[arg_counter]);
	  return(-1);
	}
      } else {
	fprintf (stderr, "must provide floating point number with sample interval (-d) flag\n");
	return(-1);
      }
    } /* sample interval */

    if (strcmp(argv[arg_counter], "-su") == 0) {
      if (arg_counter < argc) {
	if (sscanf(argv[++arg_counter], "%lf", &(utms.x)) != 1) {
	  fprintf (stderr, "source UTM X coordinate must be a floating point number, in meters; not [%s]\n",
		   argv[arg_counter]);
	  return(-1);
	}
      } else {
	fprintf (stderr, "must provide floating point number for source UTM X coordinate\n");
	return(-1);
      }

      if (arg_counter < argc) {
	if (sscanf(argv[++arg_counter], "%lf", &(utms.y)) != 1) {
	  fprintf (stderr, "source UTM Y coordinate must be a floating point number, in meters; not [%s]\n",
		   argv[arg_counter]);
	  return(-1);
	}
      } else {
	fprintf (stderr, "must provide floating point number for source UTM Y coordinate\n");
	return(-1);
      }
    } /* utm source coordinates */

    if (strcmp(argv[arg_counter], "-du") == 0) {
      if (arg_counter < argc) {
	if (sscanf(argv[++arg_counter], "%lf", &(utmd.x)) != 1) {
	  fprintf (stderr, "destination UTM X coordinate must be a floating point number, in meters; not [%s]\n",
		   argv[arg_counter]);
	  return(-1);
	}
      } else {
	fprintf (stderr, "must provide floating point number for destination UTM X coordinate\n");
	return(-1);
      }

      if (arg_counter < argc) {
	if (sscanf(argv[++arg_counter], "%lf", &(utmd.y)) != 1) {
	  fprintf (stderr, "destination UTM Y coordinate must be a floating point number, in meters; not [%s]\n",
		   argv[arg_counter]);
	  return(-1);
	}
      } else {
	fprintf (stderr, "must provide floating point number for destination UTM Y coordinate\n");
	return(-1);
      }
    } /* utm dest coordinates */

    if (strcmp(argv[arg_counter], "-sl") == 0) {
      if (arg_counter < argc) {
	if (sscanf(argv[++arg_counter], "%lf", &(lls.latitude)) != 1) {
	  fprintf (stderr, "source latitude must be a floating point number, in meters; not [%s]\n",
		   argv[arg_counter]);
	  return(-1);
	}
      } else {
	fprintf (stderr, "must provide floating point number for latitude\n");
	return(-1);
      }

      if (arg_counter < argc) {
	if (sscanf(argv[++arg_counter], "%lf", &(lls.longitude)) != 1) {
	  fprintf (stderr, "source longitude must be a floating point number, in meters; not [%s]\n",
		   argv[arg_counter]);
	  return(-1);
	}
      } else {
	fprintf (stderr, "must provide a floating point number for source longitude\n");
	return(-1);
      }

      /* Now convert the lat/long to UTM */

      ll2utm (&(lls), &(utms), WGS_84_DATUM);

    } /* lat/long source coordinates */

    if (strcmp(argv[arg_counter], "-dl") == 0) {
      if (arg_counter < argc) {
	if (sscanf(argv[++arg_counter], "%lf", &(lld.latitude)) != 1) {
	  fprintf (stderr, "destination latitude must be a floating point number, in meters; not [%s]\n",
		   argv[arg_counter]);
	  return(-1);
	}
      } else {
	fprintf (stderr, "must provide a floating point number for destination latitude\n");
	return(-1);
      }

      if (arg_counter < argc) {
	if (sscanf(argv[++arg_counter], "%lf", &(lld.longitude)) != 1) {
	  fprintf (stderr, "destination longitude must be a floating point number, in meters; not [%s]\n",
		   argv[arg_counter]);
	  return(-1);
	}
      } else {
	fprintf (stderr, "must provide a floating point number for destination longitude\n");
	return(-1);
      }

      /* Now convert the lat/long to UTM */

      ll2utm (&(lld), &(utmd), WGS_84_DATUM);

    } /* lat/long dest coordinates */

  } /* while */

  /* Read the bathymetric database */

  bfp = fopen (bathymetric_filename, "r");
  if (bfp == NULL) {
    perror ("fopen bfp");
    fprintf (stderr, "error opening bathymetric database [%s] for reading, sorry!\n",
	     bathymetric_filename);
    return(-1);
  }
  b = load_bathymetric_database (bfp);
  if (b == NULL) {
    fprintf (stderr, "error loading bathymetric database, from [%s], sorry!\n",
	     bathymetric_filename);
    return(-1);
  }
  fclose (bfp);

  /* Calculate the range between the array and the desired point */

  range = sqrt((utmd.x - utms.x) * (utmd.x - utms.x) + (utmd.y - utms.y) * (utmd.y - utms.y));

  /* Calculate the bearing between the desired point and the array */

  bearing = atan2(utmd.y - utms.y, utmd.x - utms.x);

  /* Output the start/end UTM coordinates */

  fprintf (ofp, "%% Start: %f %f (UTM WGS84)\n",
	   utms.x, utms.y);
  fprintf (ofp, "%% End:   %f %f (UTM WGS84)\n",
	   utmd.x, utmd.y);
  
  /* Output the range */

  fprintf (ofp, "%% Range: %f meters Bearing: %f radians %f degrees\n", 
	   range, bearing,
	   bearing * 180.0 / M_PI);

  /* Now step through range points */

  for (r = 0; r < range; r+=sample_interval) {
    if (range_independent_flag) {
      step.x = utms.x;
      step.y = utms.y;
    } else {
      step.x = utms.x + r * cos(bearing);
      step.y = utms.y + r * sin(bearing);
    }

    depth = find_bathymetry (step, b);

    fprintf (ofp, "%f %f %f %f\n",
	     step.x, step.y,
	     r, depth);
  }

  return(0);
}