function g=SplineInterpolation(x,y,ys,z);
% SPLINEINTERPOLATION generic interpolation with cubic splines
%   g=SplineInterpolation(x,y,ys,z); interpolates at the scalar
%   location z the data (x,y,ys) with a cubic spline function.
%   Here ys is the desired derivative at x.

n=length(x);
a=1; b=n; i=a;
while a+1~=b,                      
  i=floor((a+b)/2);
  if x(i)<z, a=i; else b=i; end
end
i=a; h=(x(i+1)-x(i));
t=(z-x(i))/h;
a0=y(i); a1=y(i+1)-a0; a2=a1-h*ys(i);
a3= h*ys(i+1)-a1; a3=a3-a2;
g=a0+(a1+(a2+a3*t)*(t-1))*t;