/********************************************************************
* This is a matlab interface to read NetCDF files 
* Senya Basin, 1994
********************************************************************/
#include <sys/types.h>
#include <malloc.h>
#include <string.h>
#include <math.h>

#include "/usr/local/include/netcdf.h"
#include "mex.h"

#define ODB_MAXNAME 64
#define SizeMN(x) (mxGetM(x)*mxGetN(x))

static void to_double(n, type, obj)
long    n;
nc_type type;
void *obj;
{
  char   *pch;
  short  *psh;
  long   *plo;
  float  *pfl;
  double *pdb = (double *)obj;
  register i;

  switch (type) {
  case NC_BYTE:
  case NC_CHAR:
    for (pch = (char *)obj,  i=n-1; i>=0; i--) pdb[i] = (double)pch[i];
    break;
  case NC_SHORT:
    for (psh = (short *)obj, i=n-1; i>=0; i--) pdb[i] = (double)psh[i];
    break;
  case NC_LONG:
    for (plo = (long *)obj,  i=n-1; i>=0; i--) pdb[i] = (double)plo[i];
    break;
  case NC_FLOAT:
    for (pfl = (float *)obj, i=n-1; i>=0; i--) pdb[i] = (double)pfl[i];
    break;
  case NC_DOUBLE:
  default:
    break;
  }
}

static unpack_sxy(nseg, ny, nx, pmask, pin, pout)
int nseg, ny, nx, *pmask;
double *pin, *pout;
{
  register i,k,m;
  double NN = sqrt(-1.);

  for (i = 0; i < nx*ny; i++) pout[i] = NN;

  for (m = i = 0; i < nseg; i+=2) 
    for (k = pmask[i]-1; k < pmask[i+1]; k++,m++) 
      pout[k] = pin[m];
}

static void var_sxy_read(idf, idv, nin, in, out)
int    idf, idv, nin;
Matrix *out[], *in[];
{
  register i;
  char str[120];
  int idm, ndim, nmask, pdim[10], dims[10];
  long size, *pmask, count[10], start[10];
  double *pntr, *data;
  nc_type type;

  ncvarinq(idf, idv, NULL, &type, &ndim, NULL, NULL); 

  if (ncattget(idf, idv, "mask_index", str) == -1)
    {printf("ODB: mask_index attribute missed in NetCDF file !!!\n");return;}
    
  ncattinq(idf, idv, "mask_index", NULL, &size);
  str[size] = '\0';
  ncvarinq(idf, (idm = ncvarid(idf, str)), NULL, NULL, &size, pdim, NULL); 
  ncdiminq(idf, pdim[0], NULL, &nmask);
  
  pmask = (long *)calloc((size_t)nmask, (size_t)sizeof(long));
  start[0] = 0L; count[0] = nmask;
  if (ncvarget(idf, idm, start, count, (void *)pmask) == -1) {
    free(pmask); return;
  }

  ncvarinq(idf, idv, NULL, NULL, NULL, pdim, NULL); 
  for (i = 0; i < ndim; i++) 
    ncdiminq(idf, pdim[i], NULL, &count[i]), start[i] = 0L;
    
  ncdiminq(idf, ncvarid(idf,"Y"), NULL, &dims[1]);
  ncdiminq(idf, ncvarid(idf,"X"), NULL, &dims[0]);
  if (dims[0]*dims[1] < count[ndim-1]) {
    printf("??? ODB:packed size %d > NX*NY=%d !!!\n",count[ndim-1],dims[0]*dims[1]);
    free(pmask); return;
  }
    
  pntr = (double *)calloc((size_t)count[ndim-1], (size_t)sizeof(double));

  if (ndim > 1) {
    if(nin < 4)
      {printf("ODB - index array missed !!!\n");return;} 
    else if ( SizeMN(in[3]) != ndim-1)
      {printf("??? index array of dimension %d expected\n", ndim-1);return;}
    else  
      for (data = mxGetPr(in[3]), i = 0; i < ndim-1; i++) {
	if ((int)data[i] > count[i] || (int)data[i] <= 0)
	  {printf("??? index: %d out of range !!!\n",(int)data[i]);return;}
	else
	  start[i] = (int)data[i]-1, count[i] = 1L;
      }
  }

  for (size = i = 1; i <= ndim; i++) size *= count[i-1];

  if (size != count[ndim-1])
    {printf("ODB: only one packed slice permitted at a time !!!\n");return;}

  if (ncvarget(idf, idv, start, count, (void *)pntr) != -1) 
    to_double (size, type, (void *)pntr);
  else 
    {free(pntr);free(pmask);return;}
    
  data = mxGetPr(out[0] = mxCreateFull(dims[0], dims[1], 0L));
    
  unpack_sxy (nmask, dims[0], dims[1], pmask, pntr, data);
	
  free(pntr); free(pmask);
}

static odb_help(key)
int key;
{
  printf("  ODB - a simplified matlab interface to read NetCDF files\n\n");
  
  if ( key == 0) {
    printf("    [object] = odb( [idf], 'COMMAND', [args...])\n\n");
    printf("    COMMAND: 'open', 'dim', 'var','attr','help'...,\n");
    printf("    may be arbitrary abbreviated. See odb('help') for details.\n");
  }
  else if (key == 1) { 
    printf("    idf = odb('open','name')\n");
    printf("          - opens file <name> for reading, returns a file ID\n");
    printf("            which may be used in further references\n");
    
    printf("    nx = odb(idf,'dim','X')\n");
    printf("          - returns a value of a dimension 'NX' from\n");
    printf("            a previously opened file with ID=idf\n");
    
    printf("    v1 = odb(idf,'var','v1')\n");
    printf("          - returns values of 1-D variable V1 as a vector v1\n");
    
    printf("    a12 = odb(idf,'var','A12')\n");
    printf("          - returns values of 2D variable <A12> as a matrix a12\n");
    printf("            (if A12 was A12(NY,NX) variable in Netcdf file then\n");
    printf("            matrix a12 has NX rows and NY columns)\n");
    
    printf("    v1 = odb(idf,'var','A12','X', ix)\n");
    printf("          - extracts a column from a 2D variable <A12>, which\n");
    printf("            corresponds to the index <ix> of the dimension <X>\n");
    
    printf("    v = odb(idf,'var','V', [ i1 i2 ... iN] )\n");
    printf("          - extracts values of N-dimensional variable <V> \n");
    printf("            as a 1-2D matrix in accordance with an index array.\n");
    printf("            If value of index in that array equal 0 the whole\n");
    printf("            dimension will be retrieved, otherwise only values\n");
    printf("            for specified index will be read.\n");
    
    printf("    a = odb(idf,'attr','V','A')\n");
  printf("          - reads the value(s) of attribute <a> of a variable <v>\n");
    
    printf("    odb('help');\n");
    printf("          - prints this help\n");

    printf("    odb('example');\n");
    printf("          - an example of using odb\n");
  }  
  else if (key == 2) {
    printf("\nExample: to plot temperature contours from the following NetCDF file:");
    printf("\nnetcdf sst {\n");
    printf("dimensions:\n");
    printf("        X = 90 ;\n");
    printf("        Y = 38 ;\n");
    printf("variables:\n");
    printf("        float X(X) ;\n");
    printf("                X:units = \"longitude\" ;\n");
    printf("        float Y(Y) ;\n");
    printf("                Y:units = \"latitude\" ;\n");
    printf("        float sst(Y, X) ;\n");
    printf("                bath:filename = \"sst\" ;\n");
    printf("                bath:description = \"Sea Surface Temperature\";\n");
    printf("                bath:missing_value = NaNf ;\n");
    printf("}\n");
    printf(" ...one may want to issue the commands' sequence:\n");
    printf(" id = odb('open','sst');\n");
    printf(" x  = odb(id,'var','X');\n");
    printf(" y  = odb(id,'var','Y');\n");
    printf(" sst = odb(id,'var','sst');\n");
    printf(" contour(x,y,sst');\n");
    printf(" title(odb(id,'attr','sst','description'));\n");
  } 
  printf("\n    Senya Basin, 1994\t\t(senya@ldeo.columbia.edu)\n");
}

mexFunction (nout, out, nin, in)
int      nin, nout;
Matrix  *in[], *out[];
{
  char str[ODB_MAXNAME], name[120];
  double *data, *pntr;
  int i, idf, idv, ndim, natt, pdim[MAX_VAR_DIMS];
  nc_type type;
  long size, start[MAX_VAR_DIMS], count[MAX_VAR_DIMS], dims[MAX_VAR_DIMS];

  if (!nin) {odb_help(0); return;}

  if ( !mxGetString(in[0], str, ODB_MAXNAME) ) {
    if (!strncasecmp(str, "open", strlen(str))) {
/**** Open *********/
      if (nin < 2 || mxGetString(in[1], name, 120) )
	mexErrMsgTxt("ODB: idf=modb('open','file') - file name expected !!!"); 
      
      data = mxGetPr(out[0] = mxCreateFull(1L, 1L, 0L));
      
      ncopts = NC_VERBOSE;
      *data = (double)ncopen(name, NC_NOWRITE);
      ncopts = NC_VERBOSE | NC_FATAL;
    }
    else if (!strncasecmp(str, "help", strlen(str))) { 
      odb_help(1); return;}
    else if (!strncasecmp(str, "example", strlen(str))) { 
      odb_help(2); return;}
    else {
      printf ("??? %s ", str);
      mexErrMsgTxt("ODB: unknown command (try odb('help')) !!!"); 
    }
  }
  else if( mxGetString(in[1], str, ODB_MAXNAME) ) 
    mexErrMsgTxt("ODB: command string expected"); 

  else if (!strncasecmp(str, "close", strlen(str))) {
/**** Close *********/
    ncopts = NC_VERBOSE;
    *data = (double)ncclose((int)mxGetScalar(in[0]));
    ncopts = NC_VERBOSE | NC_FATAL;
  }
  else if (!strncasecmp(str, "dim", strlen(str))) {
/**** Dimension *********/
    if (nin < 2 || mxGetString(in[2],name,120) )
    mexErrMsgTxt("ODB: modb(idf, 'dim','name') - dimension name expected !!!"); 

    ncopts = NC_VERBOSE;

    if (ncdiminq((int)mxGetScalar(in[0]), 
             ncdimid((int)mxGetScalar(in[0]), name), (char *)0, &size) != -1) {
      data = mxGetPr(out[0] = mxCreateFull(1L, 1L, 0L));
      *data = (double)size;
    }
    else if (nout) 
      data = mxGetPr(out[0] = mxCreateFull(1L, 1L, 0L));
    
    ncopts = NC_VERBOSE | NC_FATAL;
  }

  else if (!strncasecmp(str, "attr", strlen(str))) {
/**** Attribute *********/
    if (nin < 4 || mxGetString(in[2],name,64) || mxGetString(in[3],str,120))
    mexErrMsgTxt("ODB: modb(idf, 'attr','var','name') - variable & attribute names expected !!!"); 

    ncopts = NC_VERBOSE;
    if ( (idv = ncvarid(idf = (int)mxGetScalar(in[0]), name)) == -1 ||
	(ncattinq(idf, idv, str, &type, &size)) == -1 ){
      ncopts = NC_VERBOSE | NC_FATAL; return;} 
    
    if (type == NC_CHAR) 
      data = (double *) calloc((size_t)size, sizeof(double));
    else 
      data = mxGetPr(out[0] = mxCreateFull(1L, (long)size, 0L));
    
    if ( ncattget(idf, idv, str, data) == -1) {
      ncopts = NC_VERBOSE | NC_FATAL; return;} 
    else if (type == NC_CHAR) { 
      out[0] = mxCreateString((char *)data);
      free (data);
    }
    else
      to_double (size, type, (void *)data);

    ncopts = NC_VERBOSE | NC_FATAL; 
  }
  else if (!strncasecmp(str, "var", strlen("var"))) {
/**** Variable *********/
    if (nin < 2 || mxGetString(in[2],name,120) )
     mexErrMsgTxt("ODB: modb(idf, 'var','name') - variable name expected !!!"); 

    ncopts = NC_VERBOSE;
    if ( (idv = ncvarid(idf = (int)mxGetScalar(in[0]), name)) == -1) return; 

    ncvarinq(idf, idv, NULL, &type, &ndim, NULL, &natt); 

    ncopts = 0;
    if (natt && (ncattget(idf, idv, "packed", &size) != -1) && size) {
/******* packed data ***********/
      ncopts = NC_VERBOSE;
      if ( (ncattget(idf, idv, "compression", name) != -1) && 
	  !strncmp(name, "SXY", strlen("SXY")) ) { 

	var_sxy_read(idf, idv, nin, in, out); 
	ncopts = NC_VERBOSE | NC_FATAL;
      }
      else { 
	ncopts = NC_VERBOSE | NC_FATAL;
	mexErrMsgTxt("ODB: such a compression not implemented yet !!!");
      }
    }
    else {
/******* NOT-packed data ***********/
      ncopts = NC_VERBOSE;

      if      (ndim == 1) {
/******* 1D data ***********/
	ncvarinq(idf, idv, NULL, NULL, NULL, pdim, NULL); 
	ncdiminq(idf, pdim[0], NULL, &size);
	start[0] = 0L, count[0] = size;
	data = mxGetPr(out[0] = mxCreateFull(1L, size, 0L));
	if (ncvarget(idf, idv, start, count, (void *)data) != -1) 
	  to_double (size, type, (void *)data);
	else if (nout)
	  mxFreeMatrix(out[0]);
	ncopts = NC_VERBOSE | NC_FATAL;
      }
      else if (ndim == 2) {
/******* 2D data ***********/
	ncvarinq(idf, idv, NULL, NULL, NULL, pdim, NULL); 
	start[0] = start[1] = 0L;
	ncdiminq(idf, pdim[0], NULL, &count[0]);
	ncdiminq(idf, pdim[1], NULL, &count[1]);
	dims[0] = count[1]; dims[1] = count[0];

	if (nin > 3 && !mxGetString(in[3],name,120)) { 

	  ncdiminq(idf, pdim[0], str, NULL); 
	  if ( !strncasecmp(str, name, strlen(str)) ) {
	    if (nin < 5 || !mxIsNumeric(in[4])) 
	      mexErrMsgTxt("ODB:modb(idf,'var','varname','dimname',index) - dimension index expected !!"); 
	    start[0] = (long)mxGetScalar(in[4]) - 1;
	    dims[0] = count[0] = 1;
	    dims[1] = count[1];
	  }
	  else {
	    ncdiminq(idf, pdim[1], str, NULL); 
	    if ( !strncasecmp(str, name, strlen(str)) ) 
	      if (nin < 5 || !mxIsNumeric(in[4])) 
		mexErrMsgTxt("ODB:modb(idf,'var','varname','dimname',index) - dimension index expected !!"); 
	    start[1] = (long)mxGetScalar(in[4]) - 1;
	    dims[0] = count[1] = 1; 
	  }
	}
	    
	data = mxGetPr(out[0] = mxCreateFull(dims[0], dims[1], 0L));

	if (ncvarget(idf, idv, start, count, (void *)data) != -1) 
	  to_double (SizeMN(out[0]), type, (void *)data);
	else
	  mxFreeMatrix(out[0]);
	ncopts = NC_VERBOSE | NC_FATAL;
      }
      
      else if (ndim >= 3) {
/******* 3D & up data ***********/
	ncvarinq(idf, idv, NULL, NULL, NULL, pdim, NULL); 
	for (i = 0; i < ndim; i++) 
	  ncdiminq(idf, pdim[i], NULL, &size), start[i] = 0L, count[i] = size;
	
	if (nin != 4)
	  mexErrMsgTxt("ODB:modb(idf,'var','name',[indx]) - index array missed!!"); 
	else if (SizeMN(in[3]) != ndim ) {
	  printf("??? index array of dimension %d expected\n", ndim);
	  mexErrMsgTxt(NULL); }

	for (pntr = mxGetPr(in[3]),size = i = 0; i < ndim; i++) {
	  if  ( !(int)pntr[i] )
	    dims[size++] = count[i];
	  else if ((int)pntr[i] > count[i] || (int)pntr[i] < 0)
	    printf("??? index: %d out of range: %d.\n",(int)pntr[i],count[i]),
	    mexErrMsgTxt(NULL);
	  else
	    start[i] = (int)pntr[i]-1, count[i] = 1L;
	}
	
	if      (size == 0) dims[0] = dims[1] = 1;
	else if (size == 1) dims[1] = 1;
	else if (size > 2)
	  mexErrMsgTxt("ODB:modb(idf,'var',..) - output slab should have less then 3 dimensions !!!"); 

	data = mxGetPr(out[0] = mxCreateFull(dims[1],dims[0],0L));
	
	if (ncvarget(idf, idv, start, count, (void *)data) != -1) 
	  to_double (SizeMN(out[0]), type, (void *)data);
	else
	  mxFreeMatrix(out[0]);
	ncopts = NC_VERBOSE | NC_FATAL;
      }
    }
  }

  else {
    printf ("??? %s ", str);
    mexErrMsgTxt("ODB: unknown command (try odb('help'))!!!"); 
  }
}