%-----------------------------------------------------------------------
% Find cubic spline fit for given data
%-----------------------------------------------------------------------
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear all

% Data to be interpolated with cubic spline ----------------------------
      x=[0.0,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.8,  1.0 ]';
      y=[0.00, 0.95, 0.90, 0.95, 0.10, 0.05, 0.05, 0.20, 1.00]';
      n=size(x,1);     % Number of points in the given data vector

% Call Matlab's SPLINE function, which does not return the coefficients directly
      pp=spline(x,y);    %piecewise polynomial in Matlab's internal code

% Matlab's UNMKPP function returns the following: ----------------------
%   xnode ... the values of x
%   coeff ... (n-1)x4 matrix containing the cubic polynomial coefficients
%             n-1 pieces,
%             4 coefficients for each piece, starting with the highest order,
%             expanded around xnode

      [xnode, coeff]=unmkpp(pp)   %break into coefficients of constituent piecewise polynomials

% Plot with the help of Matlab's PPVAL function ------------------------
      xx=[0: 0.01: 1]';
      yy=zeros(size(xx));
      plot(xx, ppval(pp,xx), '-', x, y, 'o')

% To confirm the meaning of the coefficients returned by UNMKPP,
% we plot the old-fashioned, hard way ----------------------------------
      for i=1:size(xx,1)
        for j=1:n-1
          if xx(i) <= xnode(j+1)
            yy(i)=coeff(j,4);
            for iorder=1:3
              yy(i) = yy(i) + coeff(j,iorder)*(xx(i)-xnode(j))^(4-iorder);
            end
            break
          end
        end
      end
      plot(xx, yy, '-', x, y, 'o')
      title('Cubic Spline')