function [newbrk,distfn]=newknt(pp,newl)
% NEWKNT	Improve knot sequence distribution.
%
%        [newbrk,distfn]=newknt(pp,newl)
%
% returns the breakpoint sequence  newbrk  which cuts the basic interval
% for  pp  into  newl  pieces in such a way that a certain p.linear monotone
% function (whose pp-form is in  distfn ) related to the high derivative of  
%  pp  is equidistributed.
%
%   The intent is to choose a break sequence suitable to the fine 
% approximation of a function  f  whose rough approximation in  pp  is assumed
% to contain enough information about  f  to make this feasible.

% C. de Boor / latest change: May 24, 1989
% C. de Boor / latest change: Nov.14, 1990  (add missing final shift)
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.

[breaks,coefs,l,k,d] = ppbrk(pp);
if d>1, fprintf('newknt only functions for functions.\n'),end
%  it would be feasible to have it function for curves as well by choosing
%   distfn  to have as its derivative the maximum over all component
%  functions.

if l==1,   % if there is only one piece, return uniform break sequence
   newbrk=breaks(1) + [0:newl]*(breaks(l+1)-breaks(1))/newl;
   return
end

% The distribution function  distfn  is constructed as the integral of the
% k-th root of the absolute value of the k-th derivative of  pp .
%  Since  pp  is of degree < k , this requires some approximation. The ap-
% proximation to the absolute value of the k-th derivative is found as the 
% piecewise constant on the same break sequence whose value on  
% breaks(i)..breaks(i+1)  is the derivative at  
%   (break(i-1)+3break(i)+3break(i+1)+break(i+2))/8  
% of the parabola which agrees with the variation of  the (k-1)st derivative 
% of  pp  at the three points breaks(i-1/2), breaks(i+1/2), breaks(i+3/2) .

abs(diff(coefs(:,1))'./(breaks(3:l+1)-breaks(1:l-1)));
(ans([1 1:l-1])+ans([1:l-1 l-1])).^(1/k);
distfn = fnint(ppmak(breaks,ans(:)'));

% The total variation of  distfn  is its value at  break(l+1) , i.e,
var = fnval(distfn,breaks(l+1));
if var==0 ,
   newbrk=breaks(1) + [0:newl]*(breaks(l+1)-breaks(1))/newl;
   return
end

% The break sequence is to be chosen so that  
%        distfn(newbrk) = 
steps = [0:newl]*(var/newl);
newbrk=zeros(1,newl+1);newbrk([1 newl+1])=breaks([1 l+1]);

% For this, do the inverse interpolation:
[junk,coefs]=ppbrk(distfn);
flats=find(coefs(:,1)==0);
if (~isempty(flats)),coefs(flats,1)=coefs(flats,1)-1;end
breaks=breaks-breaks(1);
invbrk=[coefs(:,2)' var];
coefs=[ones(l,1)./coefs(:,1),breaks(1:l)'];
if (~isempty(flats)),coefs(flats,1)=(breaks(flats)+breaks(flats+1))'/2;end

newbrk(2:newl)=newbrk(1) + fnval(ppmak(invbrk,coefs,1),steps(2:newl));
