time   calls  line
                  1 function v=ppval(pp,xx)
                  2 %PPVAL  Evaluate piecewise polynomial.
                  3 %   V = PPVAL(PP,XX) returns the value, at the entries of XX, of the 
                  4 %   piecewise polynomial f contained in PP, as constructed by PCHIP, SPLINE,
                  5 %   INTERP1, or the spline utility MKPP.
                  6 %
                  7 %   V is obtained by replacing each entry of XX by the value of f there.
                  8 %   If f is scalar-valued, then V is of the same size as XX. XX may be ND.
                  9 %
                 10 %   If PP was constructed by PCHIP, SPLINE or MKPP using the orientation of
                 11 %   non-scalar function values specified for those functions, then:
                 12 %
                 13 %   If f is [D1,..,Dr]-valued, and XX is a vector of length N, then V has
                 14 %   size [D1,...,Dr, N], with V(:,...,:,J) the value of f at XX(J).
                 15 %   If f is [D1,..,Dr]-valued, and XX has size [N1,...,Ns], then V has size
                 16 %   [D1,...,Dr, N1,...,Ns], with V(:,...,:, J1,...,Js) the value of f at
                 17 %   XX(J1,...,Js).
                 18 %
                 19 %   If PP was constructed by INTERP1 using the orientation of non-scalar
                 20 %   function values specified for that function, then:
                 21 %
                 22 %   If f is [D1,..,Dr]-valued, and XX is a scalar, then V has size
                 23 %   [D1,...,Dr], with V the value of f at XX.
                 24 %   If f is [D1,..,Dr]-valued, and XX is a vector of length N, then V has
                 25 %   size [N,D1,...,Dr], with V(J,:,...,:) the value of f at XX(J).
                 26 %   If f is [D1,..,Dr]-valued, and XX has size [N1,...,Ns], then V has size
                 27 %   [N1,...,Ns,D1,...,Dr], with V(J1,...,Js,:,...,:) the value of f at
                 28 %   XX(J1,...,Js).
                 29 %
                 30 %   Example:
                 31 %   Compare the results of integrating the function cos and its spline
                 32 %   interpolant:
                 33 %
                 34 %     a = 0; b = 10;
                 35 %     int1 = quad(@cos,a,b);
                 36 %     x = a:b; y = cos(x); pp = spline(x,y); 
                 37 %     int2 = quad(@(x)ppval(pp,x),a,b);
                 38 %
                 39 %   int1 provides the integral of the cosine function over the interval [a,b]
                 40 %   while int2 provides the integral, over the same interval, of the piecewise
                 41 %   polynomial pp that approximates the cosine function by interpolating the
                 42 %   computed x,y values.
                 43 %
                 44 %   Class support for input X and the fields of PP:
                 45 %      float: double, single
                 46 %
                 47 %   See also SPLINE, PCHIP, INTERP1, MKPP, UNMKPP.
                 48 
                 49 %   Carl de Boor
                 50 %   Copyright 1984-2010 The MathWorks, Inc.
                 51 %   $Revision: 5.16.4.12 $  $Date: 2010/07/02 16:14:37 $
                 52 
            25   53 if isstruct(xx) % we assume that ppval(xx,pp) was used 
                 54    temp = xx; xx = pp; pp = temp;
                 55 end
                 56 
                 57 %  obtain the row vector xs equivalent to XX
            25   58 sizexx = size(xx); lx = numel(xx); xs = reshape(xx,1,lx); 
                 59 %  if XX is row vector, suppress its first dimension
  0.01      25   60 if length(sizexx)==2&&sizexx(1)==1, sizexx(1) = []; end 
                 61 
                 62 % take apart PP
            25   63 [b,c,l,k,dd]=unmkpp(pp); 
                 64 
                 65 % for each evaluation site, compute its breakpoint interval
                 66 % (mindful of the possibility that xx might be empty)
  0.57      25   67 if lx, [~,index] = histc(xs,[-inf,b(2:l),inf]); 
                 68 else index = ones(1,lx);
                 69 end
                 70 
                 71 % adjust for troubles, like evaluation sites that are NaN or +-inf
  0.09      25   72 infxs = find(xs==inf); if ~isempty(infxs), index(infxs) = l; end 
  0.05      25   73 nogoodxs = find(index==0); 
            25   74 if ~isempty(nogoodxs), xs(nogoodxs) = NaN; index(nogoodxs) = 1; end 
                 75 
                 76 % now go to local coordinates ...
  0.16      25   77 xs = xs-b(index); 
                 78 
  0.01      25   79 d = prod(dd); 
            25   80 if d>1 % ... replicate xs and index in case PP is vector-valued ... 
                 81    xs = reshape(xs(ones(d,1),:),1,d*lx);
                 82    index = d*index; temp = (-d:-1).';
                 83    index = reshape(1+index(ones(d,1),:)+temp(:,ones(1,lx)), d*lx, 1 );
            25   84 else 
            25   85    if length(sizexx)>1, dd = []; else dd = 1; end 
            25   86 end 
                 87 
                 88 % ... and apply nested multiplication:
  0.11      25   89 v = c(index,1); 
            25   90 for i=2:k 
  0.74      75   91    v = xs(:).*v + c(index,i); 
            75   92 end 
                 93 
                 94 % If evaluating a piecewise constant with more than one piece at NaN, return
                 95 % NaN.  With one piece return the constant.
            25   96 if ~isempty(nogoodxs) && k==1 && l>1 
                 97    v = reshape(v,d,lx); v(:,nogoodxs) = NaN;
                 98 end
            25   99 v = reshape(v,[dd,sizexx]); 
                100 
            25  101 if isfield(pp,'orient') && strcmp(pp.orient,'first') 
                102     % spline orientation returns    size(yi) == [d1 ... dk m1 ... mj]
                103     % but the interp1 usage prefers size(yi) == [m1 ... mj d1 ... dk]
                104     if ~(isempty(dd) || (isscalar(dd) && dd == 1))
                105         % The function is non-scalar valued
                106         if isvector(xx)&&~isscalar(xx)
                107             permVec = [ndims(v) 1:(ndims(v)-1)];
                108         else
                109             ndimsxx = ndims(xx);
                110             permVec = [(ndims(v)-ndimsxx+1) : ndims(v) 1:(ndims(v)-ndimsxx)];
                111         end
                112         v = permute(v,permVec);
                113     end
                114 end

Other subfunctions in this file are not included in this listing.