function [y,dy] = polint(xa,ya,n,x,verbose);
%
% Calculates value Y of the interpolation polynomial of order N-1
% at point X
%    
% dY is the error estimate of Y
%
% Given function values are xa(n), ya(n)
%
    
c = zeros(n);
d = zeros(n);
    
ns  = 1;
dif = abs(x-xa(1));
if(verbose==1)
    fprintf('Step  1 \n');
end;
    
for i = 1:n
    dift = abs(x-xa(i));
    if (dift<dif)
        ns  = i;
        dif = dift;
    end;
    c(i) = ya(i);
    d(i) = ya(i);
    if(verbose==1)
        fprintf(' c(%2d) = %12.7e,  d(%2d) = %12.7e \n',i,c(i),i,d(i));
    end;
end;
y = ya(ns);
if(verbose==1)
    fprintf(' ns = %2d.\n',ns);
    fprintf(' Value y in step  1 is %12.7e \n',y);
end;
ns = ns-1;
for m=1:n-1
    if(verbose==1)
        pause;
        fprintf('Step %2d \n',m+1)
    end;
    for i=1:n-m
        ho  = xa(i)  -x;
        hp  = xa(i+m)-x;
        w   = c(i+1)-d(i);
        den = ho-hp;
        if (den==0.0)
            fprintf('pause in routine POLINT');
            pause;
        end;
        den  = w/den;
        d(i) = hp*den;
        c(i) = ho*den;
        if(verbose==1)
            fprintf(' c(%2d) = %12.7e,  d(%2d) = %12.7e \n',i,c(i),i,d(i));
        end;
    end;

    if (2*ns < n-m) 
        dy = c(ns+1);
        if(verbose==1)
            fprintf(' dy = c(%2d) = %12.7e \n',ns+1,dy);
        end;
    else
        dy = d(ns);
        if(verbose==1)
            fprintf(' dy = d(%2d) = %12.7e \n',ns,dy);
        end;
        ns = ns-1;
    end;
    y = y+dy;
    if(verbose==1)
        fprintf(' Value y in step %2d is %12.7e \n',m+1,y);
    end;
end;
return;
    
