/*
 *  mexeps.c   
 *      MEX-file which provides MATLAB interface to netCDF and EPIC 
 *      System data files.
 *
 *      The time axis returned is in the form of YYMMDDhhmmss.fff.   
 *      This function  provides the access to EPIC system files 
 *      (netCDF, Classic EPIC file and ASCII files), and many other 
 *      netCDF file conventions, such as those produced by Ferret, 
 *      WHOI, etc.
 *
 *      To compile:
 *
 *        % cmex mexeps.c -I/opt/matlab/extern/include \
 *           -L/opt/matlab/extern/lib/sun4 -lmat -leps -lnetcdf -lm
 *
 *      (Makefile is available.)
 *
 *      To use this function, you must set MATLABPATH to include the
 *      directory  that  contains mexeps.m and mexeps.mexARCH, where
 *      ARCH is `4' for SunOS or `sol2' for Solaris.  At  PMEL,  for
 *      SunOS,  this  directory  is /usr/local/matlab_tools, and for
 *      Solaris, it is /opt/src/matlab_tools.
 *
 *      For a complete description of MEXEPS usage, type
 *
 *          help mexeps
 *
 *      within the MATLAB.
 *
 *--------------------------------------------------------------------
 * Copyright 1994 
 * Willa Zhu                           Phone:    206-526-6208
 * NOAA/PMEL/OCRD                      FAX:      206-526-6744
 * 7600 Sand Point Way NE              OMNET:    TAO.PMEL
 * Seattle, WA 98115                   Internet: willa@pmel.noaa.gov
 *--------------------------------------------------------------------
 *      All Rights Reserved.  Version 1.0, December 1994
 *--------------------------------------------------------------------
 */

#include <stdio.h>
#include <math.h>
#include "mex.h"
#include "epsystem.h"

Dbase *epf;
int
free_array(cptr)
     char **cptr;
{
  if(*cptr == (char *)NULL){
    *cptr = (char *)NULL;
    return 0;
  } else
    free(*cptr);

  *cptr = (char *)NULL;
  return 1;
}

static char *
mat2str(mat)
mxArray  *mat;
{
    mxArray **in;
    
    char *fcn;
    int   buflen=mxGetN(mat)+1;
    int   status, i;

    fcn=(char *)mxCalloc(buflen,sizeof(char));
    status=mxGetString(mat,fcn,buflen);
    return(fcn);
    
/*
   double *pr;
   char   *p;
   char   *str; 
   int    len, i;

   fcn=(char *)mxCalloc(buflen,sizeof(char));
   status=mxGetString(prhs[0],fcn,buflen);

   len =  mxGetM(mat)*mxGetN(mat);
   printf("len = 5d\n",len);
   str = (char *) calloc(len+1, sizeof(char));
 
   p = str;
   pr = mxGetPr(mat);
   for (i = 0; i < len; i++)      {
      *p++ = (char) *pr++;
	printf("*p = %c\n",*p);
       }
   *p = '\0';
   return(str);
*/
}

static double *
mat2double(mat)
mxArray   *mat;
{
  double *pr;
  double *out_double;
  double *p;
  int	len;
  int	i;

   pr = mxGetPr(mat);
   len =  mxGetM(mat)*mxGetN(mat);
	
   out_double = (double *) calloc(len, sizeof(double));
   p = out_double;
	
   for (i = 0; i < len; i++) *p++ = (double) *pr++;
	
   return (out_double);
}

static float *
mat2float(mat)
mxArray   *mat;
{
  double *pr;
  float *out_float;
  float *p;
  int   len;
  int   i;

   pr = mxGetPr(mat);
   len =  mxGetM(mat)*mxGetN(mat);
        
   out_float = (float *) calloc(len, sizeof(float));
   p = out_float;
        
   for (i = 0; i < len; i++) *p++ = (float) *pr++;
        
   return (out_float);
}
 
static int *
mat2int(mat)
mxArray   *mat;
{
  double *pr;
  int	 *out_int;
  int	 *p;
  int	len;
  int	i;

   pr = mxGetPr(mat);
   len =  mxGetM(mat)*mxGetN(mat);
	
   out_int = (int *) calloc(len, sizeof(int));
   p = out_int;
	
   for (i = 0; i < len; i++) *p++ = (int) *pr++;
	
   return (out_int);
}

void
time2eps(in_t,ept) 
double *in_t;
long *ept;
{  
  int year, mon, day, hour, min;
  float sec;
  double t = *in_t;

  year = (int) (t/1e+10);
  t = t-year*(1e+10);
  mon = (int) (t/1e+8);
  t = t-mon*(1e+8);
  day = (int) (t/1e+6);
  t = t-day*(1e+6);
  hour = (int) (t/1e+4);
  t = t-hour*(1e+4);
  min = (int) (t/1e+2);
  sec = (float) t- (min*(1e+2));

  if(year < 1900) year += 1900;
  mdyhms_to_ep_time(mon,day,year,hour,min,sec,ept);
}

void
check_n_input(eps_func,nrhs,lower,upper)
char *eps_func;
int nrhs, lower, upper;
{
   char msg[256];

   if (!(nrhs == lower || nrhs == upper)) {
      if(lower == upper) 
	 sprintf(msg,"mexeps: %s requires %d input arguments.",eps_func,lower);
      else
         sprintf(msg,"mexeps: %s requires %d or %d input arguments.",
		eps_func,lower,upper);
      mexErrMsgTxt(msg);
      } 
}

void
check_err(eps_f,op)
Dbase *eps_f;
char *op;
{
   char *err_msg;
   char message[256];
   char msg2[50];

   sprintf(msg2,"\nType \'help %s\' for more information.",op);
   if(ep_inq_error(eps_f) > 0) {
  	ep_message(eps_f,&err_msg);
	strcat(strcpy(message,err_msg),msg2);
	mexErrMsgTxt(message);
        }
}

void
mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
	static int file_opened = 0;
   	void    *data;
	int 	i, j, m, n, len;
	char 	*filename, *dir, *dirname;
	double 	eps_f;
	char 	*eps_func;
	int   	varlist[100], nvar;
	int 	*varcode, *varcode1, *axid, *axcode, *range;
        int     *err_opts;
	double  *dvarlist, *dshape, *daxis;
	int 	uci[4], lci[4], dims[4], origin[4];
	double 	*ddate, *ddata;
	double  *double_data;
	float   *fdata;
 	short   *sdata;
	char    *bdata;
	int	*idata;
	char  	*date;
        char 	snam[10],lnam[50],gnam[40],frm[20],unit[20];
   	char  	*sname, *lname, *gname, *frmt,*units, *type, *attname;
   	char  	*axname, *axfrmt, *axunit, *axtype;
	char    *tstr;
	int  	dtype, xid;
   	float 	spcing;
	float 	*ax_val;
   	long  	*itime, *eptime;
	int     natt, *iatt, attlen, attype;
	char	*attchar;
	float   *attfloat;
	double  *attdouble;
	int     *attint;
	int     file_opt;
	char    timestr[20];
	double  basetime;
	long    tbase[2];
	double  *todouble;
	char	*tochar;
	int     *row_axis, *col_axis;
	char *newfile;
	char    tmp_str[20];

        ep_err_opts(EPCONTW+EPPRINTERR);
	if(nrhs < 1) {
           mexErrMsgTxt("Error using ==> mexeps\nNot enough input arguments.");
	   }
	eps_func = (char *) mat2str(prhs[0]);

	/* remove following line after the testing.  - Willa */
	printf("first argument = |%s|\n",eps_func);
	if(!(
	   strcasecmp(eps_func,"OPENR") == 0 || 
	   strcasecmp(eps_func,"OPENW") == 0 ||  
	   strcasecmp(eps_func,"OPENE") == 0 ||  
	   strcasecmp(eps_func,"TOTARRAY") == 0 ||  
	   strcasecmp(eps_func,"TOTARRAY19") == 0 ||  
	   strcasecmp(eps_func,"ARRAY19TOT") == 0 ||  
	   strcasecmp(eps_func,"TOTSTR") == 0 ||  
	   strcasecmp(eps_func,"INQCODE") == 0) && !file_opened) 
	   mexErrMsgTxt("mexeps: You must open the file first with 'openr'");

	if(strcasecmp(eps_func,"OPENW") == 0) {
	   if(file_opened) ep_close(epf);
	   check_n_input(eps_func,nrhs,2,2);
/*
	   if(strcasecmp((char *) mat2str(prhs[1]),"netcdf") == 0
*/
	   file_opt = EPCREATE + EPCDF;
	   tochar = (char *)  mat2str(prhs[1]);
	   ep_open(&epf, tochar, file_opt);	
	   free_array(&tochar);
	   check_err(epf,eps_func);
	   file_opened = 1;
	   tochar = (char *) mat2str(prhs[1]);
 	   plhs[0] = mxCreateString(tochar);
	   free_array(&tochar);
	   }
	else if(strcasecmp(eps_func,"OPENR") == 0) {
	   if(file_opened) ep_close(epf);
	   check_n_input(eps_func,nrhs,2,2);
	   tochar = (char *)  mat2str(prhs[1]);
           ep_err_opts(EPCONTW+EPPRINTERR);
	   ep_open(&epf, tochar, EPREAD);	
	   free_array(&tochar);
	   check_err(epf,eps_func);
	   file_opened = 1;
     	   filename = (char *)ep_set_next(epf, (char *) NULL);
	   check_err(epf,eps_func);
	   if(nlhs == 1) 
 	      plhs[0] = mxCreateString(filename);
	   else if(nlhs == 2) {
	      ep_inq_dir(epf, &dir);
	      dirname = (char *) strdup(dir);
 	      plhs[0] = mxCreateString(dirname);
 	      plhs[1] = mxCreateString(filename);
	      }
	   }
        else if(strcasecmp(eps_func,"OPENE") == 0) {
	   if(file_opened) ep_close(epf);
	   check_n_input(eps_func,nrhs,2,2);
	   tochar = (char *)  mat2str(prhs[1]);
	   ep_open(&epf, tochar, EPEDIT);	
	   free_array(&tochar);
	   check_err(epf,eps_func);
	   file_opened = 1;
     	   filename = (char *)ep_set_next(epf, (char *) NULL);
	   check_err(epf,eps_func);
	   if(nlhs == 1) 
 	      plhs[0] = mxCreateString(filename);
	   else if(nlhs == 2) {
	      ep_inq_dir(epf, &dir);
	      dirname = (char *) strdup(dir);
 	      plhs[0] = mxCreateString(dirname);
 	      plhs[1] = mxCreateString(filename);
	      }
           }
	else if(strcasecmp(eps_func,"SETNEXTW") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,2,2);
	   tochar = (char *) mat2str(prhs[1]);
     	   (char *)ep_set_next(epf, tochar);
	   free_array(&tochar);
	   check_err(epf,eps_func);
	   }
	else if(strcasecmp(eps_func,"SETNEXTR") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,1,1);
	   if(!ep_last_datafile(epf)) { 
     	     filename = (char *)  ep_set_next(epf, (char *) NULL);
	     check_err(epf,eps_func);
	     if(nlhs == 1) 
 	       plhs[0] = mxCreateString(filename);
	     else if(nlhs == 2) {
	       ep_inq_dir(epf, &dir);
	       dirname = (char *) strdup(dir);
 	       plhs[0] = mxCreateString(dirname);
 	       plhs[1] = mxCreateString(filename);
	       }
	     }
	   else { 
	     printf("Warning: no more data files\n");
 	     plhs[0] = mxCreateString(" ");
	     if(nlhs == 2) 
 	       plhs[1] = mxCreateString(" ");
	     }
	   }
	else if(strcasecmp(eps_func,"VARLIST") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,1,1);
           ep_get_var_list(epf,varlist,100,&nvar);
	   check_err(epf,eps_func);
 	   plhs[0] = mxCreateFull(1, nvar, mxREAL);
           dvarlist = mxGetPr(plhs[0]);
	   for(i= 0; i<nvar; i++) {
	      sprintf(tmp_str,"%d",varlist[i]);  
	      sscanf(tmp_str,"%lf",&dvarlist[i]);
/*
	      dvarlist[i] = (double) varlist[i];
*/
	      }
	   }
	else if(strcasecmp(eps_func,"VARSHAPE") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,2,2);
	   varcode = mat2int(prhs[1]);
	   ep_get_var_shape(epf, *varcode, lci, uci);
	   check_err(epf,eps_func);
 	   plhs[0] = mxCreateFull(1, 4, mxREAL);
           dshape = mxGetPr(plhs[0]);
	   for(i= 0; i<4; i++) {
	      sprintf(tmp_str,"%d",uci[i]);
	      sscanf(tmp_str,"%lf",&dshape[i]);
/*
	      dshape[i] = (double) uci[i];
*/	      }
	   free_array((char **) &varcode);
	   }
        else if(strcasecmp(eps_func,"INQAXIS") == 0) {
           check_n_input(eps_func,nrhs,3,3);
           varcode = mat2int(prhs[1]);
           axid = mat2int(prhs[2]);
           ep_inq_axis(epf, *varcode, *axid, &xid,&sname,&frmt,&units,&type,&spcing);
           check_err(epf,eps_func);
           if(*axid > 3 || *axid < 1)
             mexErrMsgTxt("mexeps: INQAXIS 3rd input argument must be 1, 2 or 3");
           if(nlhs == 0) {
             printf("Axis Code %d: %s : %s : %s : %s : %f\n",xid,sname,frmt,units,type,spcing);
             }
           else if(nlhs >= 1 && nlhs <=6) {
             plhs[0] = mxCreateFull(1, 1, mxREAL);
             ddata = mxGetPr(plhs[0]);
	     sprintf(tmp_str,"%d",xid);
	     sscanf(tmp_str,"%lf",&ddata[0]);
/*
             ddata[0] = (double) xid;
*/
             if(nlhs >= 2) plhs[1] = mxCreateString(sname);
             if(nlhs >= 3) plhs[2] = mxCreateString(frmt);
             if(nlhs >= 4) plhs[3] = mxCreateString(units);
             if(nlhs >= 5) plhs[4] = mxCreateString(type);
             if(nlhs >= 6) {
               plhs[5] = mxCreateFull(1, 1, mxREAL);
               daxis = mxGetPr(plhs[5]);
	       sprintf(tmp_str,"%f",spcing);
	       sscanf(tmp_str,"%lf",&daxis[0]);
/*
               daxis[0] = (double) spcing;
*/
               }
             }
           else
             mexErrMsgTxt("mexeps: INQAXIS requires 0 to 6 output arguments");
	   free_array((char **) &varcode);
	   free_array((char **) &axid);
          }
	else if(strcasecmp(eps_func,"INQVAR") == 0) { 
	   check_n_input(eps_func,nrhs,2,2);
	   varcode = mat2int(prhs[1]);
	   ep_inq_var(epf, *varcode, &sname,&lname,&frmt,&units,&dtype);
	   check_err(epf,eps_func);
	   if(nlhs == 0) { 
             printf("Variable %d : %s : %s : %s :%s \n",
			*varcode,sname,lname,frmt,units);
	     }
	   else if(nlhs == 4) {
 	     plhs[0] = mxCreateString(sname);
 	     plhs[1] = mxCreateString(lname);
 	     plhs[2] = mxCreateString(frmt);
 	     plhs[3] = mxCreateString(units);
	     }
	   else 
             mexErrMsgTxt("mexeps: INQVAR requires 0 or 4 output arguments");
	   free_array((char **) &varcode);
	   }
	else if(strcasecmp(eps_func,"INQCODE") == 0) { 
	   check_n_input(eps_func,nrhs,2,2);
	   varcode = mat2int(prhs[1]);
           ep_inq_code(*varcode,snam,lnam,gnam,frm,unit);
	   if(nlhs == 0) { 
             printf("Variable %d : %s : %s : %s : %s : %s\n",
			*varcode,snam,lnam,gnam,frm,unit);
	     }
	   else if(nlhs == 5) {
 	     plhs[0] = mxCreateString(snam);
 	     plhs[1] = mxCreateString(lnam);
 	     plhs[2] = mxCreateString(gnam);
 	     plhs[3] = mxCreateString(frm);
 	     plhs[4] = mxCreateString(unit);
	     }
	   else 
             mexErrMsgTxt("mexeps: INQCODE requires 0 or 5 output arguments");
	   free_array((char **) &varcode);
	   }
	else if(strcasecmp(eps_func,"GETAXIS") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,3,4);
	   varcode = mat2int(prhs[1]);
	   axid = mat2int(prhs[2]);
           ep_get_var_shape(epf, *varcode, lci, uci);
	   check_err(epf,eps_func);
 	   if(nrhs == 4) {
 	      range = mat2int(prhs[3]);
	      if(range[0] >= 1 && range[1] <= uci[*axid-1] 
		  && range[0] <= range[1]) {
	        lci[*axid-1] = range[0];
	        uci[*axid-1] = range[1];
	        }
	      }
	   n =  uci[*axid-1] - lci[*axid-1] + 1;
	   if(*axid < 4) {  /* spacial axis */
              ax_val = (float *) malloc(n*sizeof(float));
              ep_get_axis(epf, *varcode, *axid, lci, uci, ax_val, n);
	      check_err(epf,eps_func);
 	      plhs[0] = mxCreateFull(n, 1, mxREAL);
              daxis = mxGetPr(plhs[0]);
	      for(i=0; i<n; i++) { 
		 sprintf(tmp_str,"%f",ax_val[i]);
	         sscanf(tmp_str,"%lf",&daxis[i]);
/*
		  daxis[i] = (double) ax_val[i];
*/
		 }
	      free(ax_val);
	      }
	   else 
             mexErrMsgTxt("mexeps: GETAXIS 3rd input argument must be 1, 2 or 3.");
	   free_array((char **) &varcode);
	   }
        else if(strcasecmp(eps_func,"COPYGRID") == 0 && file_opened) {
           check_n_input(eps_func,nrhs,3,3);
           varcode = mat2int(prhs[1]);
           varcode1 = mat2int(prhs[2]);
           ep_copy_grid(epf, *varcode, epf, *varcode1);
	   check_err(epf,eps_func);
	   }
        else if(strcasecmp(eps_func,"GETTAXIS") == 0 && file_opened) {
           check_n_input(eps_func,nrhs,2,3);
           varcode = mat2int(prhs[1]);
           ep_get_var_shape(epf, *varcode, lci, uci);
           check_err(epf,eps_func);
           if(nrhs == 3) {
              range = mat2int(prhs[3]);
              if(range[0] >= 1 && range[1] <= uci[*axid-1]
                  && range[0] <= range[1]) {
                lci[*axid-1] = range[0];
                uci[*axid-1] = range[1];
                }
	      free_array((char **) &range);
              }  
           n =  uci[3] - lci[3] + 1;
           itime = (long *) malloc(n*2*sizeof(long));
           ep_get_t_axis(epf, *varcode, lci, uci, itime, n);
           check_err(epf,eps_func);
           plhs[0] = mxCreateFull(n, 1, mxREAL);
           ddate = mxGetPr(plhs[0]);
           for (i=0; i<n; i++) {
             ep_time_to_str(itime+i*2, "YYMMDDhhmmss.fff", &date);
             ddate[i] = (double) strtod(date, (char **) NULL);
	     }
	   free_array((char **) &varcode);
	   }
	else if(strcasecmp(eps_func,"GETVAR") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,2,3);
	   varcode = mat2int(prhs[1]);
	   ep_inq_var(epf, *varcode, &sname,&lname,&frmt,&units,&dtype);
           ep_get_var_shape(epf, *varcode, lci, uci);
	   check_err(epf,eps_func);
	   if(nrhs == 3) {
 	       range = mat2int(prhs[2]);
	       for(i=0;i<4;i++) {
		   if(*(range+i) >= 1 && *(range+i+4) <= uci[i]) {
		     lci[i] = *(range+i);
		     uci[i] = *(range+i+4);
		     }
		}
	       free_array((char **) &range);
	       }
	   for(i=0;i<4;i++)
	      dims[i] =  uci[i]-lci[i]+1;
	   m = dims[0]*dims[1]*dims[2]*dims[3];
           if(dtype == EPDOUBLE) {
             data = (double *) malloc(m*sizeof(double));
             double_data = data;
             }
           else if(dtype == EPREAL) {
             data = (float *) malloc(m*sizeof(float));
             fdata = data;
             }
           else if(dtype == EPINT) {
             data = (int *) malloc(m*sizeof(int));
             idata = data;
             }
           else if(dtype == EPSHORT) {
             data = (short *) malloc(m*sizeof(short));
             sdata = data;
             }
           else if(dtype == EPCHAR) {
             data = (char *) malloc(m*sizeof(char));
             bdata = data;
             }

           ep_get_var(epf, *varcode, lci, uci, (void *) data, dims);
	   check_err(epf,eps_func);
 	   plhs[0] = mxCreateFull(m, 1, mxREAL);
           ddata = mxGetPr(plhs[0]);
           if(dtype == EPDOUBLE) 
	     for(i=0; i<m; i++) ddata[i] = (double) double_data[i];
           else if(dtype == EPREAL) 
	     for(i=0; i<m; i++) {
	        sprintf(tmp_str,"%f",fdata[i]);
		sscanf(tmp_str,"%lf",&ddata[i]);
/*
		ddata[i] = (double) fdata[i];
*/
		}
           else if(dtype == EPINT) 
	     for(i=0; i<m; i++) {
	        sprintf(tmp_str,"%d",idata[i]);
		sscanf(tmp_str,"%lf",&ddata[i]);
/*
		ddata[i] = (double) idata[i];
*/
		}
           else if(dtype == EPSHORT) 
	     for(i=0; i<m; i++) {
	        sprintf(tmp_str,"%d",sdata[i]);
		sscanf(tmp_str,"%lf",&ddata[i]);
/*
		ddata[i] = (double) sdata[i];
*/
	        }
           else if(dtype == EPCHAR) 
	     for(i=0; i<m; i++) {
	        sprintf(tmp_str,"%d",bdata[i]);
		sscanf(tmp_str,"%lf",&ddata[i]);
/*
		ddata[i] = (double) bdata[i];
*/
	        }
	   free(data);
	   free_array((char **) &varcode);
	   }
        else if(strcasecmp(eps_func,"GETSLICE") == 0 && file_opened) {
           check_n_input(eps_func,nrhs,5,5);
           varcode = mat2int(prhs[1]);
           ep_get_var_shape(epf, *varcode, lci, uci);
           check_err(epf,eps_func);
           range = mat2int(prhs[2]);
           for(i=0;i<4;i++) {
             if(*(range+i) >= 1 && *(range+i+4) <= uci[i]) {
                lci[i] = *(range+i);
                uci[i] = *(range+i+4);
                }
             }
           free_array((char **) &range);
           row_axis = mat2int(prhs[3]);
           col_axis = mat2int(prhs[4]);
	   if(*row_axis == *col_axis) 
             mexErrMsgTxt("mexeps: GETSLICE input argument 4 and 5 must not be equal.");
           for(i=0;i<4;i++) {
	     if(!((i+1) == *row_axis || (i+1) == *col_axis)) {
	       if(uci[i] != lci[i]) {
		 uci[i] = lci[i];
		 dims[i] = 1;
		 printf("warning: GETSLICE: range of axis %d reset to [%d %d]\n",i+1,lci[i],uci[i]);
		 }
	       }
             dims[i] =  uci[i]-lci[i]+1;
	     }
           m = dims[*row_axis-1];
           n = dims[*col_axis-1];
           fdata = (float *) malloc(m*n*sizeof(float));
           ep_get_var(epf, *varcode, lci, uci, (void *) fdata, dims);
           check_err(epf,eps_func);
           plhs[0] = mxCreateFull(m, n, mxREAL);
           ddata = mxGetPr(plhs[0]);
	   if(*row_axis < *col_axis) {
             for(i=0; i<n; i++)
               for(j=0; j<m; j++) {
		 sprintf(tmp_str,"%f",fdata[i*m+j]);
		 sscanf(tmp_str,"%lf",&ddata[i*m+j]);
/*
                 ddata[i*m+j] = (double) fdata[i*m+j];
*/
		 }
	     }
	   else {
             for(i=0; i<n; i++)
               for(j=0; j<m; j++) {
		 sprintf(tmp_str,"%f",fdata[j*n+i]);
		 sscanf(tmp_str,"%lf",&ddata[i*m+j]);
/*
                 ddata[i*m+j] = (double) fdata[j*n+i];
*/
	   	 }
	     }
	     
           free(fdata);
           free_array((char **) &varcode);
           }

	else if(strcasecmp(eps_func,"GETATTNUM") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,1,1);
           ep_get_att_num(epf, &natt);
	   check_err(epf,eps_func);
 	   plhs[0] = mxCreateFull(1, 1, mxREAL);
           ddata = mxGetPr(plhs[0]);
	   sprintf(tmp_str,"%d",natt);
	   sscanf(tmp_str,"%lf",&ddata[0]);
/*
	   ddata[0] = (double) natt;
*/
	   }
	else if(strcasecmp(eps_func,"ATTNAME") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,2,2);
	   iatt = mat2int(prhs[1]);
           ep_att_name(epf, *iatt, &attname);
	   check_err(epf,eps_func);
 	   plhs[0] = mxCreateString(attname);
	   free_array((char **) &iatt);
	   }
	else if(strcasecmp(eps_func,"DELATT") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,2,2);
	   ep_del_att(epf,(char *) mat2str(prhs[1]));
	   check_err(epf,eps_func);
	  }
	else if(strcasecmp(eps_func,"PUTATTVAL") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,3,3);
	   if(mxIsString(prhs[2])) {
	      ep_put_att_value(epf,(char *) mat2str(prhs[1]), EPCHAR,
			strlen((char *) mat2str(prhs[2])),
			(char *) mat2str(prhs[2]));
	      }
	   else {
	      todouble = mat2double(prhs[2]);
	      ep_put_att_value(epf,(char *) mat2str(prhs[1]), EPDOUBLE,
   			mxGetM(prhs[2])*mxGetN(prhs[2]),
			todouble);
	      free_array((char **) &todouble);
	      }
	   check_err(epf,eps_func);
	  }
	else if(strcasecmp(eps_func,"PUTVARATTVAL") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,4,4);
 	   varcode = mat2int(prhs[1]);
	   if(mxIsString(prhs[3])) {
	      ep_put_varatt_value(epf, *varcode, (char *) mat2str(prhs[2]), EPCHAR,
			strlen((char *) mat2str(prhs[3])),
			(char *) mat2str(prhs[3]));
	      }
	   else {
	      todouble = mat2double(prhs[3]);
	      ep_put_varatt_value(epf, *varcode, (char *) mat2str(prhs[2]), EPDOUBLE,
   			mxGetM(prhs[3])*mxGetN(prhs[3]),
			todouble);
	      free_array((char **) &todouble);
	      }
	   check_err(epf,eps_func);
           free_array((char **) &varcode);
	  }
	else if(strcasecmp(eps_func,"PUTAXATTVAL") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,4,4);
	   axname = (char *) mat2str(prhs[1]);
	   if(mxIsString(prhs[3])) {
	      ep_put_axatt_value(epf, axname, (char *) mat2str(prhs[2]), EPCHAR,
			strlen((char *) mat2str(prhs[3])),
			(char *) mat2str(prhs[3]));
	      }
	   else {
	      todouble = mat2double(prhs[3]);
	      ep_put_axatt_value(epf, axname, (char *) mat2str(prhs[2]), EPDOUBLE,
   			mxGetM(prhs[3])*mxGetN(prhs[3]),
			todouble);
	      free_array((char **) &todouble);
	      }
	   check_err(epf,eps_func);
	   free_array(&axname);
	  }
        else if(strcasecmp(eps_func,"SETVAR") == 0 && file_opened) {
           check_n_input(eps_func,nrhs,2,7);
           varcode = mat2int(prhs[1]);
	   if(nrhs == 7) {
             sname = (char *) mat2str(prhs[2]);
             lname = (char *) mat2str(prhs[3]);
             gname = (char *) mat2str(prhs[4]);
             frmt = (char *) mat2str(prhs[5]);
             units = (char *) mat2str(prhs[6]);
             ep_set_var(epf,*varcode,sname,lname,gname,frmt,units,EPREAL);
	     free_array(&sname);
	     free_array(&lname);
	     free_array(&gname);
	     free_array(&frmt);
	     free_array(&units);
	     }
	   else   
             ep_set_var(epf,*varcode," "," "," "," "," ",EPREAL);
           check_err(epf,eps_func);
	   free_array((char **) &varcode);
	   }
        else if(strcasecmp(eps_func,"PUTVAR") == 0 && file_opened) {
           check_n_input(eps_func,nrhs,4,5);
           varcode = mat2int(prhs[1]);
           fdata = mat2float(prhs[2]);
           range = mat2int(prhs[3]);
           for(i=0;i<4;i++) {
	      lci[i]=origin[i]=1;
	      dims[i] = *(range+i);
	      uci[i] = dims[i];
	      }
	   free_array((char **) &range);
           if(nrhs == 5) {
               range = mat2int(prhs[4]);
               for(i=0;i<4;i++) {
                   if(*(range+i) >= 1 && *(range+i+4) <= uci[i]) {
                     lci[i] = *(range+i);
                     uci[i] = *(range+i+4);
                     }
                }
	       free_array((char **) &range);
               } 
           ep_put_var(epf, *varcode, origin, lci, uci, fdata, dims);
           check_err(epf,eps_func);
	   free_array((char **) &varcode);
	   free_array((char **) &fdata);
	   } 
        else if(strcasecmp(eps_func,"REWRITEVAR") == 0 && file_opened) {
           check_n_input(eps_func,nrhs,3,3);
           varcode = mat2int(prhs[1]);
	   ep_inq_var(epf, *varcode, &sname,&lname,&frmt,&units,&dtype);
           ep_get_var_shape(epf, *varcode, lci, uci);
	   strcpy(epf->unlimited, "time");
           if(dtype == EPDOUBLE) {
             todouble = mat2double(prhs[2]);
             ep_overwrite_var(epf, *varcode, lci, uci, todouble);
	     free_array((char **) &todouble);
             }
           else if(dtype == EPREAL) {
             fdata = mat2float(prhs[2]);
             ep_overwrite_var(epf, *varcode, lci, uci, fdata);
	     free_array((char **) &fdata);
             }
           else if(dtype == EPINT) {
             idata = mat2int(prhs[2]);
             ep_overwrite_var(epf, *varcode, lci, uci, idata);
	     free_array((char **) &idata);
            }
           check_err(epf,eps_func);
	   epf->unlimited[0] = (char) NULL;
	   } 
	else if(strcasecmp(eps_func,"PUTAXIS") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,9,10);
	   varcode = mat2int(prhs[1]);
	   axid = mat2int(prhs[2]);
	   axcode = mat2int(prhs[3]);
	   axname = (char *) mat2str(prhs[4]);
	   axfrmt = (char *) mat2str(prhs[5]);
	   axunit = (char *) mat2str(prhs[6]);
	   axtype = (char *) mat2str(prhs[7]);
	   ep_set_axis(epf,*varcode,*axid,*axcode,axname,axfrmt,axunit,axtype);
	   check_err(epf,eps_func);
	   origin[0]=origin[1]=origin[2]=origin[3]=1;
	   lci[0]=lci[1]=lci[2]=lci[3]=1;
	   uci[0]=uci[1]=uci[3]=uci[3]=1;
	   fdata = mat2float(prhs[8]);
	   uci[*axid-1] = mxGetM(prhs[8])*mxGetN(prhs[8]);
	   if(nrhs == 10) {
              range = mat2int(prhs[9]);
              if(range[0] >= 1 && range[1] <= uci[*axid-1]
                  && range[0] <= range[1]) {
                lci[*axid-1] = range[0];
                uci[*axid-1] = range[1];
                }
	      free_array((char **) &range);
              }  
	   ep_put_axis(epf, *varcode, *axid, origin, lci, uci, fdata); 
	   free_array(&fdata);
	   check_err(epf,eps_func);
	   free_array((char **) &varcode);
	   free_array((char **) &axid);
	   free_array((char **) &axcode);
	   free_array(&axname);
	   free_array(&axfrmt);
	   free_array(&axunit);
	   free_array(&axtype);
	  }
        else if(strcasecmp(eps_func,"PUTTAXIS") == 0 && file_opened) {
           check_n_input(eps_func,nrhs,7,8);
           varcode = mat2int(prhs[1]);
           axname = (char *) mat2str(prhs[2]);
           axfrmt = (char *) mat2str(prhs[3]);
           axunit = (char *) mat2str(prhs[4]);
           axtype = (char *) mat2str(prhs[5]);
           ep_set_t_axis(epf,*varcode,axname,axfrmt,axunit,axtype);
           check_err(epf,eps_func);
           origin[0]=origin[1]=origin[2]=origin[3]=1;
           lci[0]=lci[1]=lci[2]=lci[3]=1;
           uci[0]=uci[1]=uci[2]=1;
           todouble = mat2double(prhs[6]);
           uci[3] = mxGetM(prhs[6])*mxGetN(prhs[6]);

           eptime = (long *) malloc(uci[3]*2*sizeof(long));
	   for(i=0;i<uci[3];i++) {
/*
	     sprintf(timestr,"%16.3f",todouble[i]);
	     str_to_ep_time("YYMMDDhhmmss.fff",timestr,(eptime+i*2));
*/
	     time2eps(todouble+i,(eptime+i*2));
	     }
	   free_array((char **)&todouble);
           if(nrhs == 8) {
              range = mat2int(prhs[9]);
              if(range[0] >= 1 && range[1] <= uci[*axid-1]
                  && range[0] <= range[1]) {
                lci[3] = range[0];
                lci[3] = range[1];
                }
	      free_array((char **) &range);
              }  
           ep_put_t_axis(epf, *varcode, origin, lci, uci, eptime);
           check_err(epf,eps_func);
	   free(eptime);
	   free_array((char **) &varcode);
	   free_array(&axname);
	   free_array(&axfrmt);
	   free_array(&axunit);
	   free_array(&axtype);
	  }
	else if(strcasecmp(eps_func,"GETATTVAL") == 0 && file_opened) { 
	   check_n_input(eps_func,nrhs,2,2);
	   ep_inq_att(epf, (char *) mat2str(prhs[1]), &attlen, &attype);	
	   check_err(epf,eps_func);
	   if(attype == EPCHAR) {
	     attchar = (char *) malloc(sizeof(char)*attlen);
	     ep_get_att_value(epf, (char *) mat2str(prhs[1]), attchar);	
	     check_err(epf,eps_func);
 	     plhs[0] = mxCreateString(attchar);
	     free(attchar);
	     }
	   else if(attype == EPINT) {
	     attint = (int *) malloc(sizeof(int)*attlen);
	     ep_get_att_value(epf, (char *) mat2str(prhs[1]), attint);	
	     check_err(epf,eps_func);
	     plhs[0] = mxCreateFull(1, attlen, mxREAL);
             ddata = mxGetPr(plhs[0]);
	     for(i=0; i<attlen; i++) { 
	       sprintf(tmp_str,"%d",attint[i]);
	       sscanf(tmp_str,"%lf",&ddata[i]);
/*
	       ddata[i] = (double) attint[i];
*/
	       }
	     free(attint);
	   }
	   else if(attype == EPREAL) {
	     attfloat = (float *) malloc(sizeof(float)*attlen);
	     ep_get_att_value(epf, (char *) mat2str(prhs[1]), attfloat);	
	     check_err(epf,eps_func);
	     plhs[0] = mxCreateFull(1, attlen, mxREAL);
             ddata = mxGetPr(plhs[0]);
	     for(i=0; i<attlen; i++) { 
	       sprintf(tmp_str,"%f",attfloat[i]);
	       sscanf(tmp_str,"%lf",&ddata[i]);
/*
	       ddata[i] = (double) attfloat[i];
*/
	       }
	     free(attfloat);
	   }
	   else if(attype == EPDOUBLE) {
	     attdouble = (double *) malloc(sizeof(double)*attlen);
	     ep_get_att_value(epf, (char *) mat2str(prhs[1]), attfloat);	
	     check_err(epf,eps_func);
	     plhs[0] = mxCreateFull(1, attlen, mxREAL);
             ddata = mxGetPr(plhs[0]);
	     for(i=0; i<attlen; i++) 
	       ddata[i] = (double) attdouble[i];
	     free(attdouble);
	   }
         }
  	else if(strcasecmp(eps_func,"TOTSTR") == 0) {
	   check_n_input(eps_func,nrhs,3,3);
           todouble = mat2double(prhs[1]);      /* input value */
	   tochar = (char *) mat2str(prhs[2]);  /* output format */
	   sprintf(timestr,"19%16.3f",todouble[0]);
           str_to_ep_time("YYYYMMDDhhmmss.fff",timestr,tbase);
	   ep_time_to_str(tbase,tochar,&tstr);
 	   plhs[0] = mxCreateString(tstr);
	   free_array((char **)&todouble);
	   free_array(&tochar);
	   }
  	else if(strcasecmp(eps_func,"TOTARRAY") == 0) {
	   check_n_input(eps_func,nrhs,3,4);
           todouble = mat2double(prhs[1]);
	   if(nrhs >= 4) {
		ddate = mat2double(prhs[3]);
		basetime = ddate[0];
		free(ddate);
		}
	   else  basetime = todouble[0];
	   len = mxGetM(prhs[1])*mxGetN(prhs[1]); 

           units = (char *) mat2str(prhs[2]);
           eptime = (long *) malloc(len*2*sizeof(long));
           for(i=0;i<len;i++) {
	     time2eps(todouble+i,(eptime+i*2));
             }
	   free_array((char **)&todouble);
	   time2eps(&basetime,tbase);
	   fdata = (float *) malloc(len*sizeof(float));
           eptime_to_t_array(eptime,len,EPREAL,units,tbase,fdata);
           free(eptime);
	   free_array(&units);

           plhs[0] = mxCreateFull(mxGetM(prhs[1]), mxGetN(prhs[1]), mxREAL);
           ddata = mxGetPr(plhs[0]);
           for(i=0; i<len; i++) {
	      sprintf(tmp_str,"%f",fdata[i]);
	      sscanf(tmp_str,"%lf",&ddata[i]);
/*
              ddata[i] = (double) fdata[i];
*/
	      }
	   free(fdata);
	   }
        else if(strcasecmp(eps_func,"ARRAY19TOT") == 0) {
           check_n_input(eps_func,nrhs,4,4);
           todouble = mat2double(prhs[1]);
           units = (char *) mat2str(prhs[2]);
           ddate = mat2double(prhs[3]);
           basetime = ddate[0];    
           free(ddate);
	   time2eps(&basetime,tbase);

           len = mxGetM(prhs[1])*mxGetN(prhs[1]);
           fdata = (float *) malloc(len*sizeof(float));
           for(i=0;i<len;i++) {
	     fdata[i] = (float) todouble[i];
	     } 
           free_array((char **)&todouble);

           eptime = (long *) malloc(len*2*sizeof(long));
           t_array_to_eptime(fdata,len,EPREAL,units,tbase,eptime);
           free(fdata);

           plhs[0] = mxCreateFull(mxGetM(prhs[1]), mxGetN(prhs[1]), mxREAL);
           ddate = mxGetPr(plhs[0]);
           for(i=0;i<len;i++) {
             ep_time_to_str(eptime+i*2, "YYYYMMDDhhmmss.fff", &date);
             ddate[i] = (double) strtod(date, (char **) NULL);
             }
           free(eptime);
           free_array(&units);
           }
  	else if(strcasecmp(eps_func,"TOTARRAY19") == 0) {
	   check_n_input(eps_func,nrhs,3,4);
           todouble = mat2double(prhs[1]);
	   if(nrhs >= 4) {
		ddate = mat2double(prhs[3]);
		basetime = ddate[0];
		free(ddate);
		}
	   else  basetime = todouble[0];
	   len = mxGetM(prhs[1])*mxGetN(prhs[1]); 

           units = (char *) mat2str(prhs[2]);
           eptime = (long *) malloc(len*2*sizeof(long));
           for(i=0;i<len;i++) {
	     time2eps(todouble+i,(eptime+i*2));
/*
	     sprintf(timestr,"%18.3f",todouble[i]);
             str_to_ep_time("YYYYMMDDhhmmss.fff",timestr,(eptime+i*2));
*/
             }
	   free_array((char **)&todouble);
/*
	   sprintf(timestr,"%18.3f",basetime);
           str_to_ep_time("YYYYMMDDhhmmss.fff",timestr,tbase);
*/
	   time2eps(&basetime,tbase);
	   fdata = (float *) malloc(len*sizeof(float));
           eptime_to_t_array(eptime,len,EPREAL,units,tbase,fdata);
           free(eptime);
	   free_array(&units);

           plhs[0] = mxCreateFull(mxGetM(prhs[1]), mxGetN(prhs[1]), mxREAL);
           ddata = mxGetPr(plhs[0]);
           for(i=0; i<len; i++) {
	      sprintf(tmp_str,"%f",fdata[i]);
	      sscanf(tmp_str,"%lf",&ddata[i]);
	      }
	   free(fdata);
	   }

	else if(strcasecmp(eps_func,"CLOSE") == 0 && file_opened) { 
           ep_close(epf);
	   check_err(epf,"eps_close");
	   file_opened = 0;
	   }
	else {
           mexErrMsgTxt("Error using ==> mexeps\nInvalid input arguments.");
	   }

	free_array(&eps_func);

	/* Assign pointers to the various parameters */

	return;
}