#ifndef __VECTOR_LINEAR_INTERPOLATION_TEMPLATE_CLASS__
#define __VECTOR_LINEAR_INTERPOLATION_TEMPLATE_CLASS__ 
/**
      @file VectorLinInterp.hpp
      @author Brian Magill
      @datecreated 7/04/2006

	  $Date:$
      $Revision:$

      @copyright
      (©) Copyright 2006 by GATS Inc.
      11864 Canon Blvd., Suite 101, Newport News, VA 23606

      All Rights Reserved. No part of this software or publication may be
      reproduced, stored in a retrieval system, or transmitted, in any form
      or by any means, electronic, mechanical, photocopying, recording, or
      otherwise without the prior written permission of GATS Inc.
 
      @brief Class for linear interpolation
*/
//----------------------------------------------------------------------
//

template <typename T>
class VectorLinInterp {

public:

    VectorLinInterp() {};
    ~VectorLinInterp() {};
    
    int Vlocate(std::vector<T> const  & xin, T const  x)
        {
        bool asc;
        int hi, mid, lo;

        lo  = -1;
        hi  = xin.size();
        asc = (xin[hi-1] > xin[0] );

        while((hi - lo) > 1)
            {
            mid = (hi + lo) / 2;
            if ( (x > xin[mid]) == asc )
                lo = mid;
            else
                hi = mid;
            }
            
        return lo;
        };
  


    T linear_interpolate(T const x, T const xa, T const xb, T const ya, T const yb)
        {
        T dx, y;
        T a;

        dx = xb - xa;
        if(dx != 0)
        {
          a = (x - xa) / dx;
          y = ya + a * (yb - ya);
        }else
          y = (ya + yb) / 2.0;

        return y;
        };


    T Interpol(std::vector<T> const &yin, std::vector<T> const  &xin, T const &xout)
        {
        int ndx1, ndx2, len;                                                    
        T yout;

        len  = xin.size();
        ndx1 = Vlocate(xin, xout);

        if(ndx1 <= 0)
          ndx1 = 0;
        if(ndx1 >= len-1)
          ndx1 = len-2;

        ndx2 = ndx1 + 1;

        yout = linear_interpolate(xout, xin[ndx1], xin[ndx2], yin[ndx1], yin[ndx2]);
  
        return yout;
        };

    std::vector<T> Interpol(std::vector<T> const  &yin, std::vector<T> const  &xin, std::vector<T>  &xout)
        {                                         
        int i;                                
        int cnt = xout.size();
        std::vector<T> yout(cnt);

        for(i=0; i<cnt; i++)
          yout[i] = Interpol(yin, xin, xout[i]);

        return yout;
        };
};

#endif