function ppstruct = simple_cubic_setup(xdat,ydat)
% This function sets up a structure that can be used for simple
% cubic interpolation of the data. The inputs xdat and ydat are
% just the data x and y values.
%
% To perform interpolations for an array of x values, use
%     y = ppval(ppstruct, x)
%
% The main task of simple_cubic_setup() is to work out coefficients
% for the piecwise cubic polynomials. The matrix of coefficients is 
% then fed to mkpp() to create an object that ppval() will understand.
%
% For examples showing the structure of the coefficient matrix,
% check out the Matlab help page for mkpp:
%  http://www.mathworks.com/access/helpdesk/help/techdoc/ref/mkpp.html

N = numel(xdat); % Find out how many data points there are.

% Calculate derivatives using the 3-point derivative formula.
dydx = ( ydat(3:N) - ydat(1:N-2) ) ./ ( xdat(3:N) - xdat(1:N-2) );
% Derivatives at ends...
dydx = [ ( ydat(2) - ydat(1) ) / ( xdat(2) - xdat(1) ), dydx, ...
           ( ydat(N) - ydat(N-1) ) ./ ( xdat(N) - xdat(N-1) ) ];

h = xdat(2:N) - xdat(1:N-1); % data spacings

% Coefficients for each Hermite basis function, from highest
% power in (x - x(i))/h(i).
phi1 = [-2, 3, 0, 0];   % controls value at end of interval.
phi0 = [2, -3, 0, 1];   % controls value at start.
psi1 = [1, -1, 0, 0];   % controls slope at end.
psi0 = [-1, 2, -1, 0];  % controls slope at start.

% Fill out the coefficients matrix for each interval in xdat. 
% Each row is an interval, and the four columns are powers of 
% x - x(i) for the cubic polynomial, from 3 down to 0.
for int = 1:N-1  % loop through all intervals
   hp = h(int).^[3,2,1,0]; % powers of h(m) needed with phi & psi
   coeffs(int,1:4) =  ydat(int)   * phi0 ./ hp ...
                    + ydat(int+1) * phi1 ./ hp ...
                    - dydx(int)   * h(int) * psi0 ./ hp ...
                    + dydx(int+1) * h(int) * psi1 ./ hp;
end

ppstruct = mkpp(xdat, coeffs);


return