time   calls  line
                  1 function output = spline(x,y,xx)
                  2 %SPLINE Cubic spline data interpolation.
                  3 %   PP = SPLINE(X,Y) provides the piecewise polynomial form of the 
                  4 %   cubic spline interpolant to the data values Y at the data sites X,
                  5 %   for use with the evaluator PPVAL and the spline utility UNMKPP.
                  6 %   X must be a vector.
                  7 %   If Y is a vector, then Y(j) is taken as the value to be matched at X(j), 
                  8 %   hence Y must be of the same length as X  -- see below for an exception
                  9 %   to this.
                 10 %   If Y is a matrix or ND array, then Y(:,...,:,j) is taken as the value to
                 11 %   be matched at X(j),  hence the last dimension of Y must equal length(X) --
                 12 %   see below for an exception to this.
                 13 %
                 14 %   YY = SPLINE(X,Y,XX) is the same as  YY = PPVAL(SPLINE(X,Y),XX), thus
                 15 %   providing, in YY, the values of the interpolant at XX.  For information
                 16 %   regarding the size of YY see PPVAL.
                 17 %
                 18 %   Ordinarily, the not-a-knot end conditions are used. However, if Y contains
                 19 %   two more values than X has entries, then the first and last value in Y are
                 20 %   used as the endslopes for the cubic spline.  If Y is a vector, this
                 21 %   means:
                 22 %       f(X) = Y(2:end-1),  Df(min(X))=Y(1),    Df(max(X))=Y(end).
                 23 %   If Y is a matrix or N-D array with SIZE(Y,N) equal to LENGTH(X)+2, then
                 24 %   f(X(j)) matches the value Y(:,...,:,j+1) for j=1:LENGTH(X), then
                 25 %   Df(min(X)) matches Y(:,:,...:,1) and Df(max(X)) matches Y(:,:,...:,end).
                 26 %
                 27 %   Example:
                 28 %   This generates a sine-like spline curve and samples it over a finer mesh:
                 29 %       x = 0:10;  y = sin(x);
                 30 %       xx = 0:.25:10;
                 31 %       yy = spline(x,y,xx);
                 32 %       plot(x,y,'o',xx,yy)
                 33 %
                 34 %   Example:
                 35 %   This illustrates the use of clamped or complete spline interpolation where
                 36 %   end slopes are prescribed. In this example, zero slopes at the ends of an 
                 37 %   interpolant to the values of a certain distribution are enforced:
                 38 %      x = -4:4; y = [0 .15 1.12 2.36 2.36 1.46 .49 .06 0];
                 39 %      cs = spline(x,[0 y 0]);
                 40 %      xx = linspace(-4,4,101);
                 41 %      plot(x,y,'o',xx,ppval(cs,xx),'-');
                 42 %
                 43 %   Class support for inputs x, y, xx:
                 44 %      float: double, single
                 45 %
                 46 %   See also INTERP1, PCHIP, PPVAL, MKPP, UNMKPP.
                 47 
                 48 %   Carl de Boor 7-2-86
                 49 %   Copyright 1984-2010 The MathWorks, Inc.
                 50 %   $Revision: 5.18.4.6 $  $Date: 2010/09/02 13:36:29 $
                 51 
                 52 
                 53 % Check that data are acceptable and, if not, try to adjust them appropriately
  0.02       7   54 [x,y,sizey,endslopes] = chckxy(x,y); 
             7   55 n = length(x); yd = prod(sizey); 
                 56 
                 57 % Generate the cubic spline interpolant in ppform
                 58 
             7   59 dd = ones(yd,1); dx = diff(x); divdif = diff(y,[],2)./dx(dd,:);  
             7   60 if n==2 
                 61    if isempty(endslopes) % the interpolant is a straight line
                 62       pp=mkpp(x,[divdif y(:,1)],sizey);
                 63    else         % the interpolant is the cubic Hermite polynomial
                 64       pp = pwch(x,y,endslopes,dx,divdif); pp.dim = sizey;
                 65    end
  0.01       7   66 elseif n==3&&isempty(endslopes) % the interpolant is a parabola 
                 67    y(:,2:3)=divdif;
                 68    y(:,3)=diff(divdif')'/(x(3)-x(1));
                 69    y(:,2)=y(:,2)-y(:,3)*dx(1); 
                 70    pp = mkpp(x([1,3]),y(:,[3 2 1]),sizey);
             7   71 else % set up the sparse, tridiagonal, linear system b = ?*c for the slopes 
             7   72    b=zeros(yd,n); 
             7   73    b(:,2:n-1)=3*(dx(dd,2:n-1).*divdif(:,1:n-2)+dx(dd,1:n-2).*divdif(:,2:n-1)); 
             7   74    if isempty(endslopes) 
             7   75       x31=x(3)-x(1);xn=x(n)-x(n-2); 
             7   76       b(:,1)=((dx(1)+2*x31)*dx(2)*divdif(:,1)+dx(1)^2*divdif(:,2))/x31; 
             7   77       b(:,n)=... 
                 78       (dx(n-1)^2*divdif(:,n-2)+(2*xn+dx(n-1))*dx(n-2)*divdif(:,n-1))/xn;
                 79    else
                 80       x31 = 0; xn = 0; b(:,[1 n]) = dx(dd,[2 n-2]).*endslopes;
                 81    end
             7   82    dxt = dx(:); 
  0.03       7   83    c = spdiags([ [x31;dxt(1:n-2);0] ... 
                 84         [dxt(2);2*(dxt(2:n-1)+dxt(1:n-2));dxt(n-2)] ...
                 85         [0;dxt(2:n-1);xn] ],[-1 0 1],n,n);
                 86 
                 87    % sparse linear equation solution for the slopes
  0.01       7   88    mmdflag = spparms('autommd'); 
  0.01       7   89    spparms('autommd',0); 
  0.01       7   90    s=b/c; 
             7   91    spparms('autommd',mmdflag); 
                 92 
                 93    % construct piecewise cubic Hermite interpolant
                 94    % to values and computed slopes
  0.01       7   95    pp = pwch(x,y,s,dx,divdif); pp.dim = sizey; 
                 96 
             7   97 end 
                 98 
  0.81       7   99 if nargin==2, output = pp; else output = ppval(pp,xx); end