csinterp( x, y, xi )

function [yi,yxi,yxxi] = csinterp( x, y, xi )

% CSINTERP 1-D piecewise cubic spline interpolation of tabulated   
%    data and calculation of first and second derivatives. 
%    [YI,YXI,YXXI] = CSINTERP(X,Y,XI) interpolates function values, 
%    first and second derivatives using the rows/columns (X,Y) at 
%    the nodes specified in XI. X and Y should be row or column 
%    vectors, each containing at least three values, otherwise the 
%    function will return an empty output. 
%    The vector X specifies the points at which the data Y is given.  
%    Output data YI, YXI and YXXI are same size as XI.
%    X and XI can be non-uniformly spaced, but repeated values in X 
%    will be removed and the resulting vector will be sorted. 
%    This function also extrapolates the input data (no NaNs are 
%    returned); use extrapolated values at your own risk. 
%
%    Requires bindex.m (available at Matlab FEX). 
% 
% Examples: 
%          x = [0.1:0.1:0.3 1:3]; y = sin( x );
%          xi = linspace(0,3,200); 
%          [yi,yxi,yxxi] = csinterp(x,y,xi);
%	   figure(1)
%	   plot(x,y,'o',xi,yi,'r-'), box on, grid on
%          xlabel('x')
%          ylabel('y(x)')
%          title('Continuity test')
%          legend('exact','csinterp')
%         
%          n = 101; m = n; 
%          x = linspace(0,2*pi,n); 
%          y = sin( x ); yx = cos( x ); yxx = -y; 
%          xi = 2*pi*rand(1,m); 
%          [yi,yxi,yxxi] = csinterp(x,y,xi);
%	   figure(2)
%          subplot(311)
%	   plot(x,y,xi,yi,'o'), box on, grid on
%          title('y(x)')
%          legend('exact','csinterp')
%	   subplot(312)
%	   plot(x,yx,xi,yxi,'o'), box on, grid on
%	   title('dy/dx')
%          legend('exact','csinterp')
%          subplot(313)
%	   plot(x,yxx,xi,yxxi,'o'), box on, grid on
%	   title('d^2y/dx^2')
%          legend('exact','csinterp')
% 
%          n = 2001; m = 5001; 
%          x = linspace(0,2*pi,n); 
%          y = sin( x ); yx = cos( x ); yxx = -y; 
%          xi = sort( 2*pi*rand(1,m) );
%          [yi,yxi,yxxi] = csinterp(x,y,xi);
%	   figure(3)
%          subplot(311)
%	   plot(x,y,'o',xi,yi,'r-'), box on, grid on
%          ylabel('y(x)')
%          title('Large dataset test')
%          legend('exact','csinterp')
%	   subplot(312)
%	   plot(x,yx,'o',xi,yxi,'r-'), box on, grid on
%	   ylabel('dy/dx')
%          legend('exact','csinterp')
%          subplot(313)
%	   plot(x,yxx,'o',xi,yxxi,'r-'), box on, grid on
%	   ylabel('d^2y/dx^2')
%          legend('exact','csinterp') 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First  version: 20/11/2007
% 
% Contact: [e-mail]
% 
% Any suggestions to improve the performance of this 
% code will be greatly appreciated. 
% 
% Reference: 
% C. Moler, Numerical Computing with Matlab (NCM)
% http://www.mathworks.com/moler
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Let's go: 

% Declare an empty output just in case calculations won't be 
% performed: 

  yi = [];
 yxi = [];
yxxi = [];

% Error checking: 

sx  = size(x);
sy  = size(y);
sxi = size(xi); 

nx = length(x); 
ny = length(y); 

if ( nargin ~= 3 )
error('This function requires three input arguments...');
end

if nx ~= ny 
error('The length of x must match the length of y...')
end 

if min(sx) ~= 1 
error('Input x should be a row or column vector...')
end

if min(sy) ~= 1 
error('Input y should be a row or column vector...')
end 

% Convert x and y to column vectors to avoid problems 
% with mixed row/column input data: 

x = x(:); 
y = y(:); 

% Remove repeated elements: 

[x,ui] = unique(x);

y = y(ui);

nx = length(x); 

% Proceed with calculations only when the # of function points 
% is greater than 2:

if nx > 2 
    
       yi = zeros( sxi );
      yxi = zeros( sxi );
     yxxi = zeros( sxi );
     
       dy = diff( y );
       
        h = diff( x ); hsq = h.*h;
      
       h1 = h(1);
       h2 = h(2);
       h3 = h(nx-2);
       h4 = h(nx-1);
       hb = h(1:nx-2);
       hf = h(2:nx-1); 

% Bracketting: 

       ix = bindex(xi,x,1); 
      
% Calculate the local variable: 
      
      if sxi(1) == 1 
      s = xi(:) - x(ix);
      else
      s = xi - x( ix );
      end
      
% Calculate the slopes:
        
      slopes = dy./h;
      
          s1 = slopes(1);
	  s2 = slopes(2);
          s3 = slopes(nx-2);
	  s4 = slopes(nx-1);	  
          sb = slopes(1:nx-2);
          sf = slopes(2:nx-1); 

% Calculate the vectors filling the matrix of the system: 
      
      lodiag(1:nx-2) = hf;
      lodiag( nx-1 ) = h3 + h4; l4 = lodiag( nx-1 ); 

      cendiag(   1  ) = h2;
      cendiag(2:nx-1) = 2*( hb + hf );
      cendiag(  nx  ) = h3;

       updiag(  1   ) = h1 + h2; u1 = updiag(1);
       updiag(2:nx-1) = hb;

% Right side of the system (not-a-knot end conditions): 

      r1 = ( ( h1 + 2*u1 )*h2*s1 + h1^2*s2 )/u1;
      rn = ( ( h4 + 2*l4 )*h3*s4 + h4^2*s3 )/l4;
      r  = 3*( hf.*sb + hb.*sf );
      r  = [r1;r;rn];
      
% Fill the matrix of the system (it is tridiagonal, so let us use a sparse matrix): 

      M = sparse( diag( cendiag ) + diag( lodiag, -1 ) + diag( updiag, 1 ) ); 

%  Solve the linear system:

      d = M\r;
      
     db = d(1:nx-1);
     df = d(2: nx );
     
% Chapter 3, page 14: 

% Calculate the remaining spline coefficients 
            
      c = ( 3*slopes - 2*db     - df )./h  ;
      b = (   db     - 2*slopes + df )./hsq;
    
% Interpolate the function and its derivatives 
      
      s2 =  s.*s; 
      s3 = s2.*s;
      
      yi = y( ix ) + s.*d( ix ) +   s2.*c( ix ) +   s3.*b( ix );
     yxi =              d( ix ) +  2*s.*c( ix ) + 3*s2.*b( ix );
    yxxi =                            2*c( ix ) +  6*s.*b( ix );
          
      if sxi(1) == 1 
         yi =   yi(:)';
        yxi =  yxi(:)';
       yxxi = yxxi(:)';
      end

end