function [out] = interp_cubic(xin,valin,xout,der)
% [out] = interp_cubic(x,val,xout,der);
% integrate quantity given on 4-points grid using cubic fit
%  x: input grid
%  val: value at input grid
%  xout: output grid
%  der: derivative order: 0: value, 1.2 first/second derivative, -1, integral

assert(size(valin,1)==4);
assert(numel(xin)==4);
assert(iscolumn(xout));

nout = numel(xout);
% four points on grid given, fit a cubic spline through them (no loss
% of accuracy for 3rd order splines.

% recursively fill matrices: A,C=[ones,x,x^2,x^3]
%        D  = [0,1,2*x,3*x^2]
%        Ci = x,1/2*x^2,1/3*x^3,1/4*x^4]
xin2 = xin.^2; xin3 = xin2.*xin;
xout2 = xout.^2; xout3 = xout2.*xout; xout4=xout2.*xout2;

A = [ones(4,1),xin,xin2,xin3];

switch der
  case -1
    I = [xout,xout2/2,xout3/3,xout4/4];
    M = bsxfun(@minus,I,I(1,:)); % subtract first output grid point to enforce zero integration constant
  case 0
    M = [ones(nout,1),xout,xout2,xout3];
  case 1
    M = [zeros(nout,1),ones(nout,1),2*xout,3*xout2];
  otherwise
    error('undefined for this der')
end

% solve linear problem for coefficients and give integral output
P = (A\valin); % coefficieents
out  = M*P;

end