function output = csder(x,y,yp,xx)

% values = csder(x,y,yp[,xx])
%
% returns the values at xx of the cubic piecewise polynomial that interpolates
% the given points (x,y) knowing the derivatives yp.
%
% There can be many columns in y to handle multiple y values. y can be complex.
% The alternative call, pp = csder(x,y,yp), returns the pp-form instead.
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% argument consistency check 
nx = length(x); [n,ny] = size(y);
if n ~= nx
 if ny == nx
  y = y.'; yp = yp.'; [n,ny] = size(y);
 else
  error('Abscissa and ordinate data should be of the same length.')
 end
end
if n < 2, error('There should be at least two data points.'), end
if any([n,ny] ~= size(yp)), error('Inconsistent derivative data size.'), end

[x,k]=sort(x(:)); y=y(k,:);
if any(diff(x)  ==  0), error('The data abscissae should be distinct.'), end
dx = diff(x); dy = diff(y) ./ dx(:,ones(ny,1)).^2;

d = (yp(2:n,:) + yp(1:n-1,:)) ./ dx(:,ones(ny,1)) - 2*dy;
c = dy - yp(1:n-1,:) ./ dx(:,ones(ny,1)) - d;
d = d ./ dx(:,ones(ny,1));
pp = ppmak(x.', reshape([d.' c.' yp(1:n-1,:).' y(1:n-1,:).'],(n-1)*ny,4), ny);

if nargin == 3
 output = pp;
else
 output = ppual(pp, xx);
end
