/****************************************************************
                     87-30-87

FUNCTION:  REGRID
Interpolate an array f(zo) onto a new grid, f(zn).  Instead of
simple linear interpolation, this routine approximates the function
as a set of line segments connecting the input data points.
Now, consider the output grid, and a grid of midpoints between
the output grid points.  The function at an output gridpoint is
then approximated as the average of the function (line-segment
approx. on the old grid) from one midpoint to the next (new grid).

regrid does not check for monotonically increasing grids or for
bad input array values; these are left as the responsibility of
the calling routine.

regrid alters only those elements in the output array for which
new values are calculated.   Therefore it can be used repeatedly
with different sections of the grid.  For example, if there are bad
data values, a function might scan f(zo) and call regrid once for
each section of contiguous good values.  If the output array had
been first initialized to all bad values, then the result of this
process would be an output array with interpolated data wherever
there were no gaps in the input array, and bad value flags elsewhere.

MODIFIED: Tue  09-27-1988
   Removed a bug causing extrapolation when the function is
   called with no overlap between the new and old grids.
   Changed the location of initialization of imin, imax,
   for readability.

***************************************************************/

#include "common.h"
#include "dbext.h"  /* min_val */
#include "vector.h"

/************      regrid     *****************************/

#if PROTOTYPE_ALLOWED
void regrid(float zo[], float fo[], float zn[], float fn[], int m, int n,
  int *n0_ptr, int *nout_ptr)
#else
void regrid(zo, fo, zn, fn, m, n, n0_ptr, nout_ptr)
float     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)       */
int       m,             /* number of points in input arrays  */
          n;             /*    "   "    "     " output  "     */
                         /* additional output:                */
                         /* calling routine supplies the space*/
int       *n0_ptr,       /* index of first interp. pt. in fn  */
          *nout_ptr;     /* number of interp. points in fn    */
#endif
{
     int  imin,          /* first new gridpoint index for fn  */
          imax;          /* last new gridpoint index for fn   */
     int  i, j;          /* index counters; i:new; j:old      */
     int  jlow, jhigh;
     float     zlow,     /* midway between zn(i) and zn(i-1)  */
               zhigh,    /* midway between zn(i) and zn(i+1)  */
               flow,     /* f(zlow)  (interpolated)           */
               fhigh;    /* f(zhigh) (interpolated)           */
     float     integral;

/****************************************************************
                 zo(jlow)              zo(jhigh)
         ----------|--------------------|-------------       zo
        ----|--------|--------|--------|---------|------   zn
          zn(i-1)   zlow    zn(i)     zhigh    zn(i+1)
******************************************************************/

/****Find limits of gridpoint index for new grid. *****************

           |-----------------------------------------|        zn (case a)
        -----------------------------------------------       zo
     ---|---------------------------------------------|------ zn (case b)
       imin                                         imax

     There can be no extrapolation in this version of regrid;
     data on the old grid must match or overhang the part of the
     new grid for which interpolation will be done.
*********************************************************************/

     imin = 0;
     imax = n-1;

     while( (zn[imin] < zo[0]) && (imin < n-1) ) imin++;
     while( (zn[imax] > zo[m-1]) && (imax > 0) )  imax--;

     /*** If there is no overlap between the old and
          new grids, there can be no interpolation.
          Set number of points regridded to zero, and
          quit.  If this occurs, then either imin = n-1
          or imax = 0; the converse is NOT true.
     ***/
     if ( (zn[imin] < zo[0]) || (zn[imax] > zo[m-1]) )
     {
          *nout_ptr = 0;
          return;
     }

     /*** Otherwise, the first interpolated point will be
          at index imin, and the last will be at imax.
     ***/
     *n0_ptr = imin;
     *nout_ptr = imax - imin + 1;


     /*** Initialize ?high to start the loop.
          In the loop ?high is immediately replaced by ?low,
          and a new ?high is calculated.
     ***/
     if (imin==0) zhigh = zn[imin];
     else
     {
          zhigh = average( zn[imin], zn[imin+1] );
          zhigh = max_val( zhigh, zo[0] );
     }
     jhigh = 1;
     while( (zo[jhigh] < zhigh) && (jhigh < m-1) ) jhigh++;
     /* This is one past initial jlow.  */
     fhigh = interpolate( zo+jhigh-1, fo+jhigh-1, zhigh);

     for ( i=imin; i<=imax; i++)
     {          /* I think this could be rearranged for compactness */
          zlow = zhigh;
          if ( i == imax)
          {
               if (imax==n-1)  zhigh = zn[imax];
               else
               {
                    zhigh = average(zn[i], zn[i+1]);
                    zhigh = min_val( zhigh, zo[m-1] );
               }
          }
          else  zhigh = average( zn[i], zn[i+1] );

          jlow = jhigh-1;      /* new jlow is 1 to left of old jhigh */
          while ( (zo[jhigh] < zhigh) && (jhigh < m-1) ) jhigh++;
          flow = fhigh;
          fhigh = interpolate(zo+jhigh-1, fo+jhigh-1, zhigh);
          if ( jhigh == (jlow+1) )
               fn[i] = interpolate(zo+jlow, fo+jlow, zn[i]);
          else                               /* all other cases       */
          {
               /** end pieces are needed first **/
               integral = average(flow,fo[jlow+1]) * (zo[jlow+1]-zlow)
                        + average(fhigh,fo[jhigh-1]) * (zhigh-zo[jhigh-1]);
               /** now the middle sections  **/
               for (j=jlow+1; j<jhigh-1; j++)
                    integral += average(fo[j], fo[j+1]) * (zo[j+1]-zo[j]);
               fn[i] = integral / (zhigh - zlow);
          }
     }
     return;
}
/**************** end of regrid()  ************************************/