/************************************************************************
* readgrd routine to read a grd file in and                             *
* interpolate to a prescribed latitude and longitude                    *
*************************************************************************/
/************************************************************************
* Creator: David T. Sandwell    Scripps Institution of Oceanography     *
* Date   : 10/02/06                                                     *
************************************************************************/
/************************************************************************
* Modification history:                                                 *
************************************************************************/

#include "../../include/gmt.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
# include <string.h>
char    *argsav[2000]; /* a hack to make gips stuff link */

int interp_topo(filein,rln,rlt,topo)

char  *filein;		/* filename of input file */
double *rln;		/* longitude -180 to 180 */
double *rlt;		/* latitude -90 to 90  */
double *topo;		/* topography(m) */

{
   static char fnameold[120]=" ";
   static float *rdat;            /* real array for output */
   static int nx;                 /* number of x points */
   static int ny;                 /* number of y points */
   static double rlt0;            /* starting latitude */
   static double rln0;            /* starting longitude */
   static double dlt;             /* latitude spacing */
   static double dln;             /* longitude spacing */
   static double rdum=9999.;      /* dummy value */
   static char  title[80];        /* title */
   int i,j,mx,my,i0,j0,indx;
   int argc2 = 1;
   char *argv2[2]= {"dummy",0};
   struct GRD_HEADER grd;

/* if this is the first call to this routine then malloc the memory 
   and read the grid */

   if(strcmp(filein,fnameold) != 0) {

/* execute GMT_begin */

	argc2 = GMT_begin (argc2, argv2);

/*  read the header */

	if (GMT_read_grd_info (filein, &grd)) {
	   fprintf (stderr, "Error opening file %s\n", filein);
	   exit (EXIT_FAILURE);
	}

/*  malloc the memory for the grid */

	nx = grd.nx;
	ny = grd.ny;
	mx = nx + 4;
	my = ny + 4;
	rdat = (float *) GMT_memory (VNULL, (size_t)(mx * my), sizeof (float), GMT_program);

/*  read the grid */

	if(GMT_read_grd(filein, &grd, rdat, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE )) {
	   fprintf (stderr, "Error reading file %s\n", filein);
	   exit (EXIT_FAILURE);
	}

	strcpy(fnameold,filein);
	dlt  = -grd.y_inc;
	dln  =  grd.x_inc;
	if(grd.node_offset == 1) {
	   rlt0 = grd.y_max;
	   rln0 = grd.x_min; 
	}
	else {
	   rlt0 = grd.y_max + 0.5*grd.y_inc;
	   rln0 = grd.x_min - 0.5*grd.x_inc; 
	}
	strncpy(title,grd.title,80); 

/*  Reset NaN to dummy values. */

	for (i = 0; i < nx * ny; i++) {
	   if(isnan (rdat[i])) rdat[i] = rdum; 
	}
   }

/* interpolate to the nearest cell */

   i0=floor((*rlt-rlt0)/dlt)+1;
   j0=floor((*rln-rln0)/dln);
   indx=i0*nx+j0;
   *topo = rdum;
   if(indx >=0 && indx < (nx*ny)) *topo = rdat[indx];
}