/********************************************************************************
% 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>

#define LONMIN            -180
#define LONMAX             180
#define LATMIN             -78
#define LATMAX              90
#define GRID               5
#define NX              4321
#define NY              2017
#define MASK            100000

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

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  FILE *infile;
  char *filename;
  char errMsg[200];
  int filename_len, status;
  float *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; 
  double *pr0, *pr1, *pr2;      /******* output matrix pointer ********/

  filename_len = 1 + (mxGetM(prhs[0]) * mxGetN(prhs[0]));
  filename = mxCalloc(filename_len, sizeof(char));
  status = mxGetString(prhs[0], filename, filename_len);
  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] = mxCreateDoubleMatrix(y_range, 1, mxREAL);
  plhs[1] = mxCreateDoubleMatrix(x_range, 1, mxREAL);
  plhs[2] = mxCreateDoubleMatrix(y_range, x_range, mxREAL); 

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

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

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

  if(!(buffer=(float *)mxCalloc((y2-y1+1) * NX,  sizeof(float)))){
    mexErrMsgTxt("Not enough memory\n");
  }
  fread((char *)&h, sizeof(struct header),1,infile);
  if(fseek(infile, y1*NX*sizeof(float), 0) == -1){
    mexErrMsgTxt("fseek error\n");
  }
  fread((char *)buffer, (y2-y1+1) * NX * sizeof(float), 1, infile);
  
  fclose(infile);

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


  for(j=0;j<x_range;j++)
    for(i=0;i<y_range;i++){
      if((*(buffer + i * m * NX + x1+j*n)) >= MASK)
	*pr2 = -1.0;
      else
	*pr2 = *(buffer + i * m * NX + x1+j*n);
      pr2++;
    }
  
  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); 
}