#include <stdio.h>
#include <stdlib.h>		/* calloc(), qsort() */
#include <math.h>		/* fabs(), sqrt() */
#include <string.h>		/* memset() */
#include "bath-util.h"

Bathymetry *
load_bathymetric_database (FILE *fp)
{
  int l, w;
  Bathymetry *b;
  int i;

  /* Allocate the Bathymetry structure */

  b = (Bathymetry *) calloc (1, sizeof(Bathymetry));
  if (b == NULL) {
    perror ("calloc bathymetry in load_bathymetry_database()");
    return(NULL);
  }

  /* Read in the number of points (length x width) */

  if (fscanf (fp, "%d %d", &l, &w) != 2) {
    fprintf (stderr, "error getting database dimensions from bathymetry file, sorry!\n");
    return(NULL);
  }

  b->npts = l * w;

  /* Allocate space for the bathymetric points */

  b->p = (BathymetricPoint *) calloc (b->npts, sizeof(BathymetricPoint));
  if (b->p == NULL) {
    perror ("calloc b->p");
    return(NULL);
  }

  /* Read in the bathymetric points */

  for (i = 0; i < b->npts; i++) {
    if (fscanf(fp, "%lf %lf %lf %s",
	       &(b->p[i].x),
	       &(b->p[i].y),
	       &(b->p[i].z),
	       b->grid_zone) != 4) {
      fprintf (stderr, "error reading in bathymetric point %d from file\n",
	       i);
    }
  }

  /* Sort the grid points */

  qsort (b->p, b->npts, sizeof(BathymetricPoint), bathymetric_compare);

  /* All set! */

  return(b);
}

int
bathymetric_compare (BathymetricPoint *p1, BathymetricPoint *p2)
{
  if (p1->x < p2->x)
    return(-1);
  if (p1->x > p2->x)
    return(1);
  if (p1->y < p2->y)
    return(-1);
  if (p1->y > p2->y)
    return(1);

  return(0);
}

double
find_bathymetry (UTM utm, Bathymetry *b)
{
  int start_el, center_el, end_el;
  double distance, min_distance;
  int distance_index;
  int i;

  /* Find which grid point is closest to the current x,y coordinate */

  start_el = 0;
  end_el = b->npts - 1;
  
  do {
    center_el = (start_el + end_el) * 0.5;
    if (b->p[center_el].x > utm.x)
      end_el = center_el;
    else
      start_el = center_el;
  } while (end_el - start_el > 1);

  i = center_el;
  distance_index = i;
  min_distance = (utm.x - b->p[i].x) * (utm.x - b->p[i].x) + 
    (utm.y - b->p[i].y) * (utm.y - b->p[i].y);

  while (++i < b->npts && fabs(b->p[i].x - utm.x) < GRID_TOLERANCE) {
    distance = (utm.x - b->p[i].x) * (utm.x - b->p[i].x) + 
      (utm.y - b->p[i].y) * (utm.y - b->p[i].y);
    if (distance < min_distance) {
      min_distance = distance;
      distance_index = i;
    }
  } 
 
  i = center_el;
  while (--i >= 0  && fabs(b->p[i].x - utm.x) < GRID_TOLERANCE) {
    distance = (utm.x - b->p[i].x) * (utm.x - b->p[i].x) + 
      (utm.y - b->p[i].y) * (utm.y - b->p[i].y);
    if (distance < min_distance) {
      min_distance = distance;
      distance_index = i;
    }
  } 

  /* Get the distance to the closest point */

  distance = sqrt(min_distance);

  return(b->p[distance_index].z);
}


Bathymetric *
obtain_bathymetric_data (double utm_sx, double utm_sy,
			 double utm_dx, double utm_dy,
			 unsigned int range_independent_flag)
{
  /*
   * Obtains the bathymetric data and returns it in
   * a series of range, depth pairs
   */
  int cmd_length;
  char *cmd_line, *output_line;
  Bathymetric *b;
  FILE *pp;

  /* Allocate the bathymetric structure */

  b = (Bathymetric *) calloc (1, sizeof(Bathymetric));
  if (b == NULL) {
    perror ("calloc b in obtain_bathymetric_data()");
    return(NULL);
  }
  b->bp = (BathymetricPoint *) calloc (b->nr + 1, sizeof(BathymetricPoint));
  if (b->bp == NULL) {
    perror ("calloc b->bp in obtain_bathymetric_data()");
    return(NULL);
  }

  /* Form the command line */

  cmd_length = strlen(BATHYMETRIC_PROG) + strlen(BATHYMETRIC_FILE) + BUFFER_LENGTH;
  
  cmd_line = (char *) calloc (cmd_length, sizeof(char));
  if (cmd_line == NULL) {
    perror ("calloc cmd_line in obtain_bathymetric_data()");
    return(NULL);
  }

  if (range_independent_flag)

    sprintf (cmd_line, "%s -ri %s -d %f -su %f %f -du %f %f",
	     BATHYMETRIC_PROG, 
	     BATHYMETRIC_FILE,
	     BATHYMETRIC_SAMPLE_INTERVAL,
	     utm_sx,
	     utm_sy,
	     utm_dx,
	     utm_dy);
  else
    sprintf (cmd_line, "%s %s -d %f -su %f %f -du %f %f",
	     BATHYMETRIC_PROG, 
	     BATHYMETRIC_FILE,
	     BATHYMETRIC_SAMPLE_INTERVAL,
	     utm_sx,
	     utm_sy,
	     utm_dx,
	     utm_dy);

  /* Call the bathymetry extraction program */

  output_line = (char *) calloc (BUFFER_LENGTH, sizeof(char));
  if (output_line == NULL) {
    perror ("calloc output_line in obtain_bathymetric_data()");
    return(NULL);
  }

  fprintf (stderr, "running [%s]\n", cmd_line);

  pp = popen (cmd_line, "r");
  if (pp == NULL) {
    perror ("pp in obtain_bathymetric_data()");
    return(NULL);
  }

  while (!feof(pp)) {
    memset (output_line, 0, BUFFER_LENGTH);
    if (fgets (output_line, BUFFER_LENGTH, pp) != NULL) {
      /*
       * Check to make sure the line is not a comment
       */
      if (output_line[0] != '%') {
	if (sscanf(output_line, "%lf %lf %lf %lf",
		   &(b->bp[b->nr].x),
		   &(b->bp[b->nr].y),
		   &(b->bp[b->nr].r),
		   &(b->bp[b->nr].z)) != 4) {
	  fprintf (stderr, "malformed response from bathymetric program: ");
	  fprintf (stderr, "[%s]\n",
		   output_line);
	} else {
	  (b->nr)++;

	  b->bp = (BathymetricPoint *)
	    realloc (b->bp, (b->nr + 1) * sizeof(BathymetricPoint));
	  if (b->bp == NULL) {
	    perror ("calloc b->bp");
	    return(NULL);
	  }

	  /* Check to see if the bathymetry goes above zero */

	  if (b->bp[b->nr - 1].z < 0.0)
	    b->bp[b->nr - 1].z = 0.0;

	} /* good bathymetric point */
      } /* not a comment line */
    } /* good input line */
  } /* while !feof(pp) */

  /* Close the program */

  pclose(pp);

  fprintf (stderr, "done\n");

  /* Free the cmd and output lines */

  free (cmd_line);
  free (output_line);

  /* All done; return bathymetric information */

  return(b);
}

Bathymetric *
reverse_bathymetric_data (char *reverse_bathymetry_file,
			  double sx, double sy,
			  double bx, double by)
{
  /*
   * This routine reads in a file which contains a bathymetric slope.
   * It assumes the beginning of the file (first coordinate) is the bottom
   * of the FFP array (bx, by).  It then re-arranges the slope so the FFP array
   * is the _end_ of the slope.  This is done so MFP runs will have consistent bathymetry
   * between different trial source locations
   */

  /*
   * Obtains the bathymetric data and returns it in
   * a series of range, depth pairs
   */
  Bathymetric *tb, *b;
  FILE *fp;
  char *output_line;
  double range;
  int i;

  /* Allocate bathymetric structures (both returned and temporary) */

  b = (Bathymetric *) calloc (1, sizeof(Bathymetric));
  if (b == NULL) {
    perror ("calloc b in reverse_bathymetric_data()");
    return(NULL);
  }
  b->bp = (BathymetricPoint *) calloc (b->nr + 1, sizeof(BathymetricPoint));
  if (b->bp == NULL) {
    perror ("calloc b->bp in reverse_bathymetric_data()");
    return(NULL);
  }

  tb = (Bathymetric *) calloc (1, sizeof(Bathymetric));
  if (tb == NULL) {
    perror ("calloc tb in reverse_bathymetric_data()");
    return(NULL);
  }
  tb->bp = (BathymetricPoint *) calloc (tb->nr + 1, sizeof(BathymetricPoint));
  if (tb->bp == NULL) {
    perror ("calloc tb->bp in reverse_bathymetric_data()");
    return(NULL);
  }

  /* Allocate output line */

  output_line = (char *) calloc (BUFFER_LENGTH, sizeof(char));
  if (output_line == NULL) {
    perror ("calloc output_line in reverse_bathymetric_data()");
    return(NULL);
  }

  fp = fopen (reverse_bathymetry_file, "r");
  if (fp == NULL) {
    perror ("fopen in reverse_bathymetric_data()");
    fprintf (stderr, "unable to open reverse bathymetry file: [%s]\n",
	     reverse_bathymetry_file);
    return(NULL);
  }

  while (!feof(fp)) {
    memset (output_line, 0, BUFFER_LENGTH);
    if (fgets (output_line, BUFFER_LENGTH, fp) != NULL) {
      /*
       * Check to make sure the line is not a comment
       */
      if (output_line[0] != '%') {
	if (sscanf(output_line, "%lf %lf %lf %lf",
		   &(tb->bp[tb->nr].x),
		   &(tb->bp[tb->nr].y),
		   &(tb->bp[tb->nr].r),
		   &(tb->bp[tb->nr].z)) != 4) {
	  fprintf (stderr, "malformed response from bathymetric file: ");
	  fprintf (stderr, "[%s]\n",
		   output_line);
	} else {

	  (tb->nr)++;

	  tb->bp = (BathymetricPoint *)
	    realloc (tb->bp, (tb->nr + 1) * sizeof(BathymetricPoint));
	  if (tb->bp == NULL) {
	    perror ("calloc tb->bp in reverse_bathymetric_data()");
	    return(NULL);
	  }

	  /* Check to see if the bathymetry goes above zero */

	  if (tb->bp[tb->nr - 1].z < 0.0)
	    tb->bp[tb->nr - 1].z = 0.0;

	} /* good bathymetric point */
      } /* not a comment line */
    } /* good input line */
  } /* while !feof(fp) */

  /* Close the file */

  fclose(fp);

  /* Free the cmd and output lines */

  free (output_line);

  /* Sanity check the file */

  if (fabs (bx - tb->bp[0].x) > 1.0 ||
      fabs (by - tb->bp[0].y) > 1.0) {
    fprintf (stderr, "WARNING: reverse bathymetry from file [%s] does not match!\n",
	     reverse_bathymetry_file);
    fprintf (stderr, "file start:    %f/%f\n",
	     tb->bp[0].x, tb->bp[0].y);
    fprintf (stderr, "desired start: %f/%f\n",
	     bx, by);
    return(NULL);
  }

  /* Reverse the file */

  range = sqrt((sx - bx) * (sx - bx) + (sy - by) * (sy - by));

  for (i = tb->nr - 1; i >= 0; i--) {
    if ((range - tb->bp[i].r) > 0.0) {
      b->bp[b->nr].x = tb->bp[i].x;
      b->bp[b->nr].y = tb->bp[i].y;
      b->bp[b->nr].z = tb->bp[i].z;
      b->bp[b->nr].r = range - tb->bp[i].r;

      /* Reallocate */

      (b->nr)++;
      b->bp = (BathymetricPoint *) realloc (b->bp, ((b->nr + 1) * sizeof(BathymetricPoint)));
      if (b->bp == NULL) {
	perror ("realloc b->bp in reverse_bathymetric_file()");
	return(NULL);
      }
    } /* range check loop */
  } /* i loop */

  /* Free temporary bathymetric structure */

  free (tb->bp);
  free (tb);
    
  /* All done; return bathymetric information */

  return(b);
}

Bathymetric *
read_bathymetric_data (char *bathymetry_file,
		       double sx, double sy)

{
  /*
   * This routine reads in a file which contains a bathymetric slope.
   * It assumes the beginning of the file (first coordinate) is the bottom
   * of the FFP array (sx, sy).  This is done to minimize the number
   * of times the (slow) obtain_bathymetry program needs to be run.
   */

  /*
   * Obtains the bathymetric data and returns it in
   * a series of range, depth pairs
   */
  Bathymetric *b;
  FILE *fp;
  char *output_line;

  /* Allocate bathymetric structure */

  b = (Bathymetric *) calloc (1, sizeof(Bathymetric));
  if (b == NULL) {
    perror ("calloc b in read_bathymetric_data()");
    return(NULL);
  }
  b->bp = (BathymetricPoint *) calloc (b->nr + 1, sizeof(BathymetricPoint));
  if (b->bp == NULL) {
    perror ("calloc b->bp in read_bathymetric_data()");
    return(NULL);
  }

  /* Allocate output line */

  output_line = (char *) calloc (BUFFER_LENGTH, sizeof(char));
  if (output_line == NULL) {
    perror ("calloc output_line in read_bathymetric_data()");
    return(NULL);
  }

  fp = fopen (bathymetry_file, "r");
  if (fp == NULL) {
    perror ("fopen in read_bathymetric_data()");
    fprintf (stderr, "unable to open bathymetry file: [%s]\n",
	     bathymetry_file);
    return(NULL);
  }

  while (!feof(fp)) {
    memset (output_line, 0, BUFFER_LENGTH);
    if (fgets (output_line, BUFFER_LENGTH, fp) != NULL) {
      /*
       * Check to make sure the line is not a comment
       */
      if (output_line[0] != '%') {
	if (sscanf(output_line, "%lf %lf %lf %lf",
		   &(b->bp[b->nr].x),
		   &(b->bp[b->nr].y),
		   &(b->bp[b->nr].r),
		   &(b->bp[b->nr].z)) != 4) {
	  fprintf (stderr, "malformed response from bathymetric file: ");
	  fprintf (stderr, "[%s]\n",
		   output_line);
	} else {

	  (b->nr)++;

	  b->bp = (BathymetricPoint *)
	    realloc (b->bp, (b->nr + 1) * sizeof(BathymetricPoint));
	  if (b->bp == NULL) {
	    perror ("calloc b->bp in reverse_bathymetric_data()");
	    return(NULL);
	  }

	  /* Check to see if the bathymetry goes above zero */

	  if (b->bp[b->nr - 1].z < 0.0)
	    b->bp[b->nr - 1].z = 0.0;

	} /* good bathymetric point */
      } /* not a comment line */
    } /* good input line */
  } /* while !feof(fp) */

  /* Close the file */

  fclose(fp);

  /* Free the cmd and output lines */

  free (output_line);

  /* Sanity check the file */

  if (fabs (sx - b->bp[0].x) > 1.0 ||
      fabs (sy - b->bp[0].y) > 1.0) {
    fprintf (stderr, "WARNING: bathymetry from file [%s] does not match!\n",
	     bathymetry_file);
    fprintf (stderr, "file start:    %f/%f\n",
	     b->bp[0].x, b->bp[0].y);
    fprintf (stderr, "desired start: %f/%f\n",
	     sx, sy);
    return(NULL);
  }
    
  /* All done; return bathymetric information */

  return(b);
}