/****************************************************************
   TABLE1.C: source for Mex function to take the place of
   the table1 Matlab m-file.  The function should be very
   similar, but much faster.

   Difference from table1.m:  out of range values of x 
   generate error conditions in table1.m but NaNs in the output
   of table1.c.


   Eric Firing
   Sun  08-22-1993  17:01:25


95/11/10 - EF - Cleaned up, commented, and debugged.
98/09/11 - jules - added comments for linux compilation
     - still compliing with: "mex -V4 table1.c" 
     - linux is a PC_COMPATIBLE_HOST and handles prototypes 
***************************************************************/
#include "mex.h"

/* The following defines are largely taken from codas3/include/dbhost.h
   and matfile.h.
   They have been extracted and put here to make table1.c independent
   of codas.
*/
#define PC_COMPATIBLE_HOST   0 /* including linux on an x86 */
#define VAXD_COMPATIBLE_HOST 1
#define SUN3_COMPATIBLE_HOST 2 /* including Alliant, HPUX */

/* #define USE_FUNCTION_PROTOTYPES */ /* uncomment this for linux*/


/* #define HOST_ENVIRONMENT PC_COMPATIBLE_HOST */
#define HOST_ENVIRONMENT SUN3_COMPATIBLE_HOST

/* The following define is needed to avoid conflict with the
   declaration in mex.h as of Matlab 4.2; for compilers that
   don't handle prototypes, it declares 
      mexFunction();
   instead of
      void mexFunction();
*/
/*#define VOID void */ /* for compilers that handle prototypes */ 
#define VOID       /* for Sunos 4.x cc */

#if HOST_ENVIRONMENT == PC_COMPATIBLE_HOST
#define NaN_d_0 0L
#define NaN_d_1 0xFFF80000L

#elif HOST_ENVIRONMENT == SUN3_COMPATIBLE_HOST
#define NaN_d_0 0x7FFFFFFFL
#define NaN_d_1 0xFFFFFFFFL

#elif HOST_ENVIRONMENT == VAXD_COMPATIBLE_HOST
#define NaN_d_0 0x1D97FF5DL
#define NaN_d_1 0x687CD88DL
#endif


#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))


void regridli();
void interpolate();
void to_NaN();
int ISNAN();

/* In the following, "columns" refers to Matlab
   usage--numbers are stored consecutively in a column.
   The following variables are global in this module for
   efficiency and compactness.
*/
static int    nold,      /* number of rows in input columns       */
              nnew,      /*    "   "    "  "  output  "           */
              ncol;      /* number of columns to be interpolated  */

typedef union
{
  double d;
  long l[2];
} NAN_TYPE;

NAN_TYPE NaN; /* values are set in the mexFunction */

/* Testing for NaN can't be done directly with the double
   precision NaN itself, so we need an equivalent to the
   Matlab "isnan" function.
*/
int ISNAN(x)
double x;
{
   extern NAN_TYPE NaN;
   NAN_TYPE xn;

   xn.d=x;
   return ( xn.l[0] == NaN.l[0] && xn.l[1] == NaN.l[1]  );
}

/* Specialized function to set all columns in a given row
   to NaN:
*/
void to_NaN(fi)
double     *fi;
{
   extern NAN_TYPE NaN;
   extern int nnew, ncol;
   int i, offset;

   for (i=0; i<ncol; i++)
   {
      offset = i*nnew;
      *(fi+offset) = NaN.d;
   }
}

/* Specialized function to interpolate all columns of a
   given row:
*/
void interpolate(z, f, fi, zi)
double     *z, *f, *fi;
double    zi;
{
   extern int nold, nnew, ncol;
   double    dzratio, f0, f1;
   int i, off_old, off_new;

   dzratio =   (zi - *(z-1)) / (*z - *(z-1));
   for (i=0; i<ncol; i++)
   {
      off_new = i*nnew;
      off_old = i*nold;
      f0 = *(f+off_old-1);
      f1 = *(f+off_old);
      if (ISNAN(f0) || ISNAN(f1))
         *(fi+off_new) = NaN.d;
      else
         *(fi+off_new) = f0 + dzratio * (f1 - f0);
      /* Explicit handling of NaN is just in case some
      compiler does not meet IEEE standard.
      */
   }

}



/************      regridli     *****************************/

void regridli(zo, fo, zn, fn)
double    zo[],          /* input array of m depths: "z_old"  */
          fo[],          /* input array of m; f(z_old)        */
          zn[],          /* input array of n depths: "z_new". */
          fn[];          /* output array of n; f(z_new)       */
{
   extern int nold, nnew;
   int  i;
   int  jhigh;
   int  increasing;    /* 1 if zo increases, 0 if it decreases */
   int       jmax;
   double    zz;


   if (nold < 2)
   {
      mexErrMsgTxt("Need more than one row in the table.");
   }
   if (zo[1] > zo[0])
      increasing = 1;
   else
      increasing = 0;

   /* Make sure the old grid is monotonic: */
   if (increasing)
   {
      for (i = 1; i < nold; i++)
      {
         if (zo[i] <= zo[i-1])
         {
            mexErrMsgTxt("Table is not monotonic.");
         }
      }
   }
   else  /* decreasing */
   {
      for (i = 1; i < nold; i++)
      {
         if (zo[i] >= zo[i-1])
         {
            mexErrMsgTxt("Table is not monotonic.");
         }
      }
   }

   /* Interpolate each row: */
   jmax = nold - 1;
   if (increasing)
   {
      jhigh = 1;
      for (i=0; i<nnew; i++)
      {
         zz = zn[i];
         while (zo[jhigh] < zz && jhigh < jmax)
            jhigh++;
         while (zo[jhigh-1] > zz && jhigh > 1)
            jhigh--;
         if ( zo[jhigh] < zz || zo[jhigh-1] > zz )
         {
            to_NaN(fn+i);
         }
         else
         {
            interpolate(zo+jhigh, fo+jhigh, fn+i, zz);
         }
      }
   }
   else /* decreasing */
   {
      jhigh = jmax;
      for (i=0; i<nnew; i++)
      {
         zz = zn[i];
         while (zo[jhigh] > zz && jhigh < jmax)
            jhigh++;
         while (zo[jhigh-1] < zz && jhigh > 1)
            jhigh--;
         if ( zo[jhigh] > zz || zo[jhigh-1] < zz )
         {
            to_NaN(fn+i);
         }
         else
         {
            interpolate(zo+jhigh, fo+jhigh, fn+i, zz);
         }
      }
   }
   return;
}
/* The search method in regridli is sequential, starting
from the last position found.  This should be reasonably
efficient when the new grid is sorted, and very efficient
when the new grid is as fine as, or finer than, the old
grid.  If these conditions are not true, then a binary
search would be superior.
*/
/**************** end of regridli()  ************************************/


#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
{
   extern NAN_TYPE NaN;
   extern int nold, nnew, ncol;
   double *xold, *yoldr, *yoldi, *xnew, *ynewr, *ynewi;
   int nr, nc;

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


   if (nrhs != 2)
   {
      mexErrMsgTxt("Need 2 input arguments.");
   }
   if (nlhs > 1)
   {
      mexErrMsgTxt("Only 1 output argument permitted.");
   }

   /* Second argument: new x grid. */
   nr = mxGetM(prhs[1]);  nc = mxGetN(prhs[1]);
   if (!(nc == 1 | nr == 1))
   {
      mexErrMsgTxt("Need single row or column in the new x vector.");
   }
   nnew = MAX(nr, nc);
   xnew = mxGetPr(prhs[1]);

   /* First argument: table = [old_x old_y]. */
   nr = mxGetM(prhs[0]);  nc = mxGetN(prhs[0]);
   if (nc < 2)
   {
      mexErrMsgTxt("Need at least two columns in the table.");
   }
   nold = nr;
   ncol = nc-1;  /* First input column is x, not y. */
   xold = mxGetPr(prhs[0]);
   yoldr = xold + nold;
   yoldi = mxGetPi(prhs[0]);

   if (yoldi)
   {
      yoldi += nold;
      plhs[0] = mxCreateFull(nnew, ncol, COMPLEX);
      ynewi = mxGetPi(plhs[0]);
   }
   else
   {
      plhs[0] = mxCreateFull(nnew, ncol, REAL);
   }
   ynewr = mxGetPr(plhs[0]);

   regridli(xold, yoldr, xnew, ynewr);
   if (yoldi)
   {
      regridli(xold, yoldi, xnew, ynewi);
   }
}