/****************************************************************
   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
2000/11/22 - EF - updated for Matlab 5/6 compatibility

***************************************************************/
#include "mex.h"
#include "matrix.h"


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


void regridli(double zo[], double fo[], double zn[], double fn[]);
void interpolate(double *z, double *f, double *fi, double zi);
void to_NaN(double *fi);

/* 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 double NAN_TYPE;

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


/* Specialized function to set all columns in a given row
   to NaN:
*/
void to_NaN(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;
   }
}

/* Specialized function to interpolate all columns of a
   given row:
*/
void interpolate(double *z, double *f, double *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 (mxIsNaN(f0) || mxIsNaN(f1))
         *(fi+off_new) = NaN;
      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(double zo[], double fo[], double zn[], double fn[])
/*        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()  ************************************/


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
   extern NAN_TYPE NaN;
   extern int nold, nnew, ncol;
   double *xold, *yoldr, *yoldi, *xnew, *ynewr, *ynewi;
   int nr, nc;

   NaN = mxGetNaN();


   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] = mxCreateDoubleMatrix(nnew, ncol, mxCOMPLEX);
      ynewi = mxGetPi(plhs[0]);
   }
   else
   {
      plhs[0] = mxCreateDoubleMatrix(nnew, ncol, mxREAL);
   }
   ynewr = mxGetPr(plhs[0]);

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