% function a = fit(x, ain, fnvals)
%
%     x      abscissa values of x to do the fit at
%   ain      initial coefficient values
%            polyfit gives a useful set
%   fnvals   function to be approximated values at the abscissa locations
%
%   returns the coefficients of the equiripple polynomial

function a = fit(x, ain, fnvals);

% minimax polynomial fitter
%
% based on Hildebrand's Introduction to
% Mathematical Analysis text .. Dover
%
% 20Feb00 K.Metzger..quick and dirty..needs cleaning

a = ain;              % copy initial polynomial coefficients
N=length(x);          % determine number of x values
Np = length(ain)-1;   % determine polynomial order
fig=1;
for iter = [1:3]      % iterate the procedure this many times
   y = polyval(a,x);
   diff = y-fnvals;
   %figure(fig); fig=fig+1; clf
   delta = sign(diff(1:N-1) - diff(2:N));
   %plot(diff);
   test = (delta == 0);
   delta = delta+test;
   index=[];
   j = 1; 
   for i=1:N-2          % use changes of sign in slope to find maxima
      if delta(i)*delta(i+1) < 0
          index(j)=i;   % save indices of maxima locations
          j=j+1;
      end
   end
   if length(index)==(Np)    % need to add both ends
       index = [1 index N];
   elseif length(index)==(Np+1)  % need to bigest end
       if abs(diff(1))>abs(diff(N))
           index = [1 index];
       else
           index = [index N];
       end
   end
   xn = x([index]).';
   xm = zeros(Np+2, Np+2);
   xpwr = ones(Np+2, 1);
   signx = 1;
   for i = 1:Np+1     % form matrix of maxima location powers
      xm(:, Np+2-i) = xpwr;
      xm(i, Np+2) = signx;  % also need to include E term
      signx = -signx;       % E term has alternating signs
      xpwr = xpwr.*xn;      % advance power
   end
   xm(Np+2, Np+2) = signx;
   fm = fnvals([index]).';  % function values at chosen maxima
   a_next =  inv(xm)*fm;    % solve for updated coefficients
   a = a_next(1:Np+1);      % copy into working array
   E = a_next(Np+2);        % also solved for attained error
end
 E % print out the attained error value