function output = csmod(x,y,ec,xx)

% values = csmod(x,y[,ec][,xx])
%
% returns the values at xx of the cubic spline interpolant to the
% given data (x,y), using the following edge condition :
%
% ec = 'a' : not-a-knot 3rd der. continuous at inner knots (default)
%      'n' : natural    2nd der. zero at outer knots
%      's' : symetric   y(x1+x) = y(x1-x)
%      'p' : periodic   y(x+xN) = y(x)
%
% There can be many columns in y to handle multiple y values. y can be complex.
% If xx is missing, returns the pp-form of the cubic spline interpolant instead.
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% argument consistency check
if nargin == 2
 ec = 'ap';
elseif nargin ==3
 if isstr(ec)
  ec = [ec 'p'];
 else
  xx = ec;
  ec = 'a';
 end
end
nx = length(x); [n,ny] = size(y);
if n ~= nx
 if ny == nx
  y = y.'; [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

[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));

if (n == 2), % the interpolant is a straight line
 pp = ppmak(x.', [dy.' y(1,:).'], ny);

elseif n == 3, % the interpolant is a parabola
 c = diff(dy) / (x(3) - x(1));
 b = dy(1,:) - c * dx(1);
 pp = ppmak([x(1) x(3)], [c.' b.' y(1,:).'], ny);

else % set up the linear system for solving for the slopes.
 A = diag([dx(2:n-1);0],-1) + 2*diag([0;dx(2:n-1)+dx(1:n-2);0]) + diag([0;dx(1:n-2)],1);
 B = zeros(n,ny);
 B(2:n-1,:) = 3*(dx(2:n-1,ones(ny,1)).*dy(1:n-2,:) + dx(1:n-2,ones(ny,1)).*dy(2:n-1,:));

 % edge condition
 if     ec(1) == 'a' % not-a-knot
  x31 = x(3) - x(1); xn = x(n) - x(n-2);
  A(1,1:2) = [dx(2) x31];
  A(n,n-1:n) = [xn dx(n-2)];
  B(1,:) = ((dx(1)   + 2*x31)*dx(2)  *dy(1,:)   + dx(1)^2  *dy(2,:))  /x31;
  B(n,:) = ((dx(n-1) + 2*xn )*dx(n-2)*dy(n-1,:) + dx(n-1)^2*dy(n-2,:))/xn;
 elseif ec(1) == 'n' % natural
  A(1,1:2) = dx(1)*[2 1]; A(n,n-1:n) = dx(n-1)*[1 2]; % could have been 1,
  B(1,:) = 3*dx(1)  *dy(1,:); 
  B(n,:) = 3*dx(n-1)*dy(n-1,:);                       % but better scaled this way  
 elseif ec(1) == 's' % symetric
  xn = x(n) - x(n-2);
  A(1,1) = dx(1);
  A(n,n-1:n) = [xn dx(n-2)];
  B(n,:) = ((dx(n-1) + 2*xn )*dx(n-2)*dy(n-1,:) + dx(n-1)^2*dy(n-2,:))/xn;
 elseif ec(1) == 'p' % periodic
  A(1,1:2) = [2*(dx(1)+dx(n-1)) dx(n-1)]; A(1,n-1) = dx(1);
  A(n,1) = -dx(1); A(n,n) = dx(1);
  B(1,:) = 3*(dx(1)*dy(n-1,:) + dx(n-1)*dy(1,:));
 end

 % convert to pp form
 b = A\B;
 d = (b(1:n-1,:) + b(2:n,:) - 2*dy(1:n-1,:))./dx(:,ones(ny,1));
 c = (dy(1:n-1,:) - b(1:n-1,:))./dx(:,ones(ny,1)) - d;
 pp = ppmak(x.', reshape([(d./dx(:,ones(ny,1))).' c.' b(1:n-1,:).' y(1:n-1,:).'],(n-1)*ny,4), ny);
end

if length(ec) == 2
 output = pp;
else
 output = ppual(pp,xx);
end
