// Interpolate.cc
// Thu Dec 20 06:00:03 PST 2007 Kevin Karplus

#include <assert.h>
#include "Interpolate.h"

// Do linear interpolation based on a table
// The table consists of two vectors, x_tab and y_tab, of the same
//	length, sorted in increasing order of x.
// If x is outside the range of x_tab, then extrapolate.
// Clipping can be done instead of extrapolation by duplicating the
//	first or last point of the table
double interpolate(double x, 
	const vector<double> &x_tab,
	const vector<double> &y_tab)
{
    assert(x_tab.size()==y_tab.size());
    
    int use_point;	// interpolate on line through points
    			// (x_tab[use_point], y_tab[use_point])
    			// (x_tab[use_point+1], y_tab[use_point+1])
    
    switch(x_tab.size())
    { case 0:
    	assert(x_tab.size()>0);
	// BUG: should throw an error rather than just failing the assertion.
	return 0;
      case 1:
      	// with only one data point in table, we can't interpolate or extrapolate,
	//	so just return that value (constant function)
	return y_tab[0];
      case 2:
        // with only two data points, do interpolation/extrapolation 
	// on line through them
    	use_point=0;
	break;
      case 3:
      	use_point =  (x>=x_tab[1])? 1: 0;
	break;
      default:
      	// more than 3 points, first check if x out of range:
	if (x<=x_tab[0])
	{   use_point=0;
	    break;
	}
	if (x>=x_tab[x_tab.size()-1])
	{   use_point=x_tab.size()-2;
	    break;
	}
	
	// search for use_point such that
	//	x_tab[use_point] <= x < x_tab[use_point+1]
	int left=0, right=x_tab.size()-1;
	while (right>left +1)
	{   int mid = (left+right)/2;
	    if (x< x_tab[mid])  {right=mid;}
	    else {left=mid;}
	}
	use_point=left;
    }
    
    double x_diff = x_tab[use_point+1] - x_tab[use_point];
    if (x_diff==0) return y_tab[use_point];
    double frac = (x-x_tab[use_point])/x_diff;
    return y_tab[use_point] + (y_tab[use_point+1] - y_tab[use_point])*frac;
    
}
	


// Do linear interpolation based on a table
// The table consists of two arrays, x_tab and y_tab, of the same
//	length, sorted in increasing order of x.
// If x is outside the range of x_tab, then extrapolate.
// Clipping can be done instead of extrapolation by duplicating the
//	first or last point of the table
double interpolate(double x, 
	unsigned int size,
	const double *x_tab,
	const double *y_tab)
{
    int use_point;	// interpolate on line through points
    			// (x_tab[use_point], y_tab[use_point])
    			// (x_tab[use_point+1], y_tab[use_point+1])
    
    switch(size)
    { case 0:
    	assert(size>0);
	// BUG: should throw an error rather than just failing the assertion.
	return 0;
      case 1:
      	// with only one data point in table, we can't interpolate or extrapolate,
	//	so just return that value (constant function)
	return y_tab[0];
      case 2:
        // with only two data points, do interpolation/extrapolation 
	// on line through them
    	use_point=0;
	break;
      case 3:
      	use_point =  (x>=x_tab[1])? 1: 0;
	break;
      default:
      	// more than 3 points, first check if x out of range:
	if (x<=x_tab[0])
	{   use_point=0;
	    break;
	}
	if (x>=x_tab[size-1])
	{   use_point=size-2;
	    break;
	}
	
	// search for use_point such that
	//	x_tab[use_point] <= x < x_tab[use_point+1]
	int left=0, right=size-1;
	while (right>left +1)
	{   int mid = (left+right)/2;
	    if (x< x_tab[mid])  {right=mid;}
	    else {left=mid;}
	}
	use_point=left;
    }
    
    double x_diff = x_tab[use_point+1] - x_tab[use_point];
    if (x_diff==0) return y_tab[use_point];
    double frac = (x-x_tab[use_point])/x_diff;
    return y_tab[use_point] + (y_tab[use_point+1] - y_tab[use_point])*frac;
    
}


// CHANGE LOG:
// Thu Dec 20 13:20:13 PST 2007 Kevin Karplus
//	Added array version of interpolate() and
//	added case 3: as special case to speed things up a little.
// Sat May 31 20:01:14 PDT 2008 Kevin Karplus
//	Added assert.h