/********************************************************************************
% format: to run in matlab 
%   [top_lat, top_lon, topog] = get_topog(latmin, latmax, lonmin, lonmax];
%   output:
%          top_lat: m x 1
%	  top_lon: n x 1
%	  topog:   m x n
**********************************************************************************/

#include <stdio.h>
#include <malloc.h>
#include <mex.h>
#include "/home/noio/programs/codas3/include/matfile.h"

#define LONMIN            -180
#define LONMAX             180
#define LATMIN             -78
#define LATMAX              90
#define GRID               5
#define NX              4321
#define NY              2017
#define MASK            65535
/* #define INPUTFILE       "/home/noio/programs/bathy/dbdb5.fin" */

/***** header structure, although not used  ********/
struct header{
  float min_max[4];
  float grid;
  int nx, ny;
  float zeros[4314];
};

#ifdef USE_FUNCTION_PROTOTYPES
void mexFunction(int nlhs,Matrix *plhs[], int nrhs,Matrix *prhs[])
#else
mexFunction(nlhs, plhs, nrhs, prhs)
     int nlhs, nrhs;
     Matrix *plhs[], *prhs[];
#endif
{
  FILE *infile;
  char filename[1000];
  char errMsg[200];
  unsigned short int *buffer;     /******* store all the data read from the input file ***/
  struct header h;
  float latmin, latmax, lonmin, lonmax;
  float del_lat, del_lon;     /***** del_lat, del_lon *****/
  int i,j, x1, x2, y1, y2, x_range, y_range, n, m, k; 
  double *pr0, *pr1, *pr2;      /******* output matrix pointer ********/
  char *m2s();

/***** 4 input arguments(latmin, latmax, lonmin, lonmax ****/

  if(mxGetString(prhs[0], filename, 1000))
    mexErrMsgTxt("Error in the filename");
  latmin = (float)*mxGetPr(prhs[1]);
  latmax = (float)*mxGetPr(prhs[2]);
  lonmin = (float)*mxGetPr(prhs[3]);
  lonmax = (float)*mxGetPr(prhs[4]);

  switch(nrhs){
  case 7:
    n = (int)*mxGetPr(prhs[6]);
    m = (int)*mxGetPr(prhs[5]);
    break;
  case 6:
    n = (int)*mxGetPr(prhs[5]);
    m = (int)*mxGetPr(prhs[5]);
    break;
  default:
    n = 1;
    m = 1;
  }

/*  if(latmin > latmax || lonmin > lonmax){
    mexErrMsgTxt(" Note: lonmin must be less than lonmax(similarly latmin < latmax is required)\n	If you must cross the date line (e.g. 178W to 178E) use \n       lonmin = -182, lonmax = -178 or lonmin=178, lonmax=182");
  }
*/
  x1=( (lonmin - LONMIN) * 60 / GRID) + 0.5;
  x2=( (lonmax - LONMIN) * 60 / GRID) + 0.5;
  y1=( (latmin - LATMIN) * 60 / GRID) + 0.5;
  y2=( (latmax - LATMIN) * 60 / GRID) + 0.5;

  del_lat = (float)GRID / 60 * m;
  del_lon = (float)GRID / 60 * n;

  y_range=(y2-y1) / m + 1;
  x_range=(x2-x1) / n + 1;

  mexPrintf("  get_topog is extracting %d X %d (lat X Lon) matrix\n", y_range, x_range);

/****** initialize output matrix ******/

  plhs[0] = mxCreateFull(y_range, 1, REAL);
  plhs[1] = mxCreateFull(x_range, 1, REAL);
  plhs[2] = mxCreateFull(y_range, x_range, REAL); 

  pr0 = mxGetPr(plhs[0]);
  pr1 = mxGetPr(plhs[1]);
  pr2 = mxGetPr(plhs[2]);

  if(!(infile = fopen(filename,"r"))){
    sprintf(errMsg,"%s can't be open to read\n", filename); 
    mexErrMsgTxt(errMsg);
  }

/***** read from the input file the whole data to buffer *******/

  if(!(buffer=(unsigned short int *)malloc((y2-y1+1) * NX * sizeof(unsigned short int)))){
    mexErrMsgTxt("Not enough memory\n");
  }
  fread((char *)&h, sizeof(struct header),1,infile);
  if(fseek(infile, (y1+1)*NX*sizeof(unsigned short int), 0) == -1){
    mexErrMsgTxt("fseek error\n");
  }
  fread((char *)buffer, (y2-y1+1) * NX * sizeof(unsigned short int), 1, infile);
/*  mexPrintf("%d %d\n", buffer[0], buffer[1]);*/
  fclose(infile);

/***** write the specific range to output matrix ******/

  k=0;
/*  for(j=0;j<x_range;j++)
    for(i=0;i<y_range;i++){
      if((*(buffer + i * m * NX + x1+j*n)) >= MASK)
	pr2[k] = -1.0;
      else
	pr2[k] = buffer[i * m * NX + x1+j*n];
      k++;
    }
  */

  for(j=0;j<x_range;j++)
    for(i=0;i<y_range;i++){
      tonan_d(&pr2[k], &buffer[i*m*NX+x1+j*n]);
      k++;
    }

  for(i=0; i < y_range; i++)
    *(pr0+i)=latmin + i * del_lat;
  for(i=0; i < x_range; i++)
    *(pr1+i)=lonmin + i * del_lon;
  
  free(buffer); 
}


/* 
**   m2s() - function to return string from MEX Matrix argument
*/

char   *m2s(mat)
   Matrix	*mat;
{

   double  *pr;
   char    *str,*p;
   int     len;
   
   len = mxGetN(mat);
   str = (char *)calloc(len+1,sizeof(char));
   p = str;
   pr = mxGetPr(mat);
   while (len--)
      *p++ = (char)*pr++;
   *p = '\0';
   return(str);

} /* end of function m2s */

void tonan_d(dest, source)
double *dest;
unsigned short *source;
{
   int i;
   union { double d; long l[2]; } NaN;

   NaN.l[0] = NaN_d_0;
   NaN.l[1] = NaN_d_1;

/*mexPrintf("here %d \n", *source); */
   *dest =  (*source < MASK) ? *source : NaN.d; 

}