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

FUNCTION:  REGRIDLI   (derived from INTGRID)

This function is needed for finding isotherm depths, and in
similar cases where local linear interpolation is needed.  A
vector specified on an "old" grid is interpolated to a "new"
grid. Points on the new grid that are not bracketed by old
grid points are left unchanged.

Unlike the other members of the family (INTGRID and REGRID),
REGRIDLI can handle increasing or decreasing grids.  Like the
others, it does not check for monotonicity or gaps.  It
gives a warning if the function to be interpolated contains
bad values, but not if the grid vectors contain bad values.

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.

   Wed  10-05-1988
   Removed a bug that prevented regridding from occurring
   with a grid that decreases with index.  Clarified some of
   the code.

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

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

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

#if PROTOTYPE_ALLOWED
void regridli(float zo[], float fo[], float zn[], float fn[], int m, int n,
  int *n0_ptr, int *nout_ptr)
#else
void regridli(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 good point in fn   */
          *nout_ptr;     /* number of good points in fn       */
                         /*  (should be n)                    */
#endif
{
     int  imin,          /* first new gridpoint index for fn  */
          imax;          /* last new gridpoint index for fn   */
     int  i;
     int  jhigh;
     int  s;             /* 1 if z increases, -1 if it decreases */
     double zomin, zomax;
/****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 regridli;
     data on the old grid must match or overhang the part of the
     new grid for which interpolation will be done.
*********************************************************************/

     /* Find the extrema of the old grid. */
     zomin = BADFLOAT;
     zomax = -BADFLOAT;
     for (i=0; i<m; i++)
     {
          zomin = min_val(zomin, zo[i]);
          zomax = max_val(zomax, zo[i]);
     }

     /* Test whether the new grid increases or decreases with
     index.
     */

     /* 1999/07/13 Pierre Jaccard
        If there is one grid point in the new grid, zn[1] is undefined. The
        test below must therefore be modified
        
        Old code:
        if (zn[1]>zn[0])  s = 1;     
        
        New code:
     */
     if ((n==1) || (zn[1]>zn[0]))  s = 1;  /* z increases with index */

     else              s = -1;

     /* Find the index range of the new grid that is
        bracketed by the old grid.
     */
     imin = 0;
     imax = n-1;

     if (s==1)
     {
          while( (zn[imin] < zomin) && (imin < n-1) ) imin++;
          while( (zn[imax] > zomax) && (imax > 0) ) imax--;
          if ( (zn[imin] < zomin) || (zn[imax] > zomax) )
          {
               *nout_ptr = 0;
               return;
          }
     }
     else
     {
          while( (zn[imin] > zomax) && (imin < n-1) ) imin++;
          while( (zn[imax] < zomin) && (imax > 0) ) imax--;
          if ( (zn[imin] > zomax) || (zn[imax] < zomin) )
          {
               *nout_ptr = 0;
               return;
          }
     }

     *n0_ptr = imin;
     *nout_ptr = imax - imin +1;

/* If enabled, the following statements would set any uninterpolated
   ranges to BADFLOAT.  Instead, we will leave this to the calling
   routine.  In some cases it may be better to have the uninterpolated
   range set to zero, or left unchanged.  This is also consistent with
   REGRID.C.
*/
#if 0
     for (i=0; i<imin; i++) fn[i] = BADFLOAT;
     for (i=imax+1; i<n; i++) fn[i] = BADFLOAT;
#endif
     for ( jhigh=1, i=imin; i<=imax; i++)
     {
          if (s==1)
               while ( ((zo[jhigh] < zn[i]) 
                         || (zo[jhigh-1] > zn[i]))
                         && (jhigh < m-1) ) jhigh++;
          else
               while ( ((zo[jhigh] > zn[i])
                         || (zo[jhigh-1] < zn[i]))
                         && (jhigh < m-1) ) jhigh++;

          if (fo[jhigh-1]>=ADJ_BADFLOAT || fo[jhigh]>=ADJ_BADFLOAT)
          {
               fn[i] = BADFLOAT;
               printf("\nBADFLOAT found in input to REGRIDLI, i= %d\n", i);
          }
          else
          {
               fn[i] = interpolate(zo+jhigh-1, fo+jhigh-1, zn[i]);
          }
     }
     return;
}
/**************** end of regridli()  ************************************/