% CHEBLOOP	Used by CHEBDEM, Chebyshev demonstration.
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
%                                          % latest change: February 13, '90
echo on
% For the complete levelling, we use the  Remez algorithm. This means that we 
% construct new tau's as the extrema of our current approximation to  c  
% and try again.  We find the new interior  tau  as the zeros of  c' .

     pause;                      % touch any key to continue

cp = fnder(c);

% We take the zeros of the corresponding control polygon as our first
% guess. For this, we must take apart the spline  cp .

[knots,coefs,np,kp] = spbrk(cp);

% The control polygon has the vertices  (tstar(i),coefs(i)) , with
%  tstar  the knot averages for the spline as supplied by

tstar=aveknt(knots,kp);

% Here are the zeros of the resulting control polygon of  cp :
npp=[1:np-1];
guess=tstar(npp) - coefs(npp).*(diff(tstar)./diff(coefs));
plot(guess,zeros(1,np-1),'o');pause

% This provides already a very good first guess for the actual zeros

     pause;                      % touch any key to continue


% Now we compute the first derivative at both these points
points=ones(4,1)*tau(2:n-1);
points(1,:)=guess;
values=0*points;
values(1,:)= fnval(cp,points(1,:));
values(2,:)= fnval(cp,points(2,:));

% use two steps of the secant method, but guard against division by zero:
for j=2:3
   rows=[j,j-1];cpd=diff(values(rows,:));
   index=find(cpd==0);cpd(index)=ones(1,length(index));
   points(j+1,:)=points(j,:)-values(j,:).*(diff(points(rows,:))./cpd);
   values(j+1,:)=fnval(cp,points(j+1,:));
end
   max(abs(values'))

     pause;                      % touch any key to continue

% This is a good improvement. Now we take these  points  as our  tau :
tau = [tau(1) points(4,:) tau(n)]
% .. and check the extreme values of our current approximation there
extremes=abs(fnval(c,tau));
% The difference
max(extremes)-min(extremes)
% is an estimate of how far we are from total levelling. 
echo off       % for a quick check
if ~exist('answer'),
% format long
% ans
% format short
zcheck(ans-0.69047027491413);end;echo on
% Here is the picture of the new spline corresponding to that
c=spapi(t,tau,b);
points=sort([tau [0:100]*(t(n+1)-t(k))/100]);
values=fnval(c,points);plot(points,values);pause

% If this is not close enough, simply try again, starting from these new  tau .
%
echo off
  fprintf(' If you want to continue, type  1 , then RETURN.')
  answer = input(' Otherwise, just hit RETURN:  ');
if ~isempty(answer)&(answer==1), chebloop
else,hold off
end
