//
//-----------------------------------------------------------------------
/// @copyright
///            (c) Copyright 2008 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.
///
//-----------------------------------------------------------------------
///
/// @file	Interpolate.cpp
///
/// @author	John Burton	<j.c.burton@gats-inc.com>
///
/// @date  	Thu May 29 14:31:06 2008
///
//-----------------------------------------------------------------------
//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//
#include "Interpolate.h"

//
//-----------------------------------------------------------------------
// Defines and Macros:
//-----------------------------------------------------------------------
//



//
//-----------------------------------------------------------------------
// Global Variables:
//-----------------------------------------------------------------------
//



//
//-----------------------------------------------------------------------
// Class Methods:
//-----------------------------------------------------------------------
//


static double linterp(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;
}

static double Interpolate(EventVar<int> &loc, EventVar<double> &yin,
			  EventVar<double> &xin, double xout)
{
  int ndx1, ndx2, len;
  double  yout;

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

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

  ndx2 = ndx1 + 1;

  if(loc[ndx1] == loc[ndx2])
    yout = linterp(xout,xin[ndx1],xin[ndx2],yin[ndx1],yin[ndx2]);
  else
    yout = yin[ndx1];

  return yout;
}

EventVar<double> VInterpolate(EventVar<int> &loc, EventVar<double> &yin, 
			      EventVar<double> &xin, EventVar<double> &xout)
{
  int nin  = xin.size();
  int nout = xout.size();
  EventVar<double> yout(nout);

  for(int i=0; i<nout; i++)
  {
    yout[i] = Interpolate(loc,yin,xin,xout[i]);
  }
  return yout;
}