/**
      @file ValInterp.cpp
      @author Brian Magill
      @datecreated 4/07/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
*/
//----------------------------------------------------------------------
//

#include "ValInterp.h"

using namespace std;

int Vlocate(valarray<double> xin, double 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;
    }
  


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

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

    return y;
    }


double Interpol(valarray<double> const &yin, valarray<double> const &xin, double xout)
    {
    int ndx1, ndx2, len;                                                      
    double 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;
    }

valarray<double> VInterpol(valarray<double> const &yin, 
                            valarray<double> const &xin, valarray<double> const &xout)
    {
    int i;                                  
    int cnt = xout.size();
    valarray<double> yout(cnt);

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

    return yout;
    }