function yy=rat_interpolation(x,y,xx)
% Function for interpolation using
% rational functions.
% x and y are input data points
% yy are the function values of the interpolant
% at points xx.
% From Fausett's book

[k1,k]=size(x); [kk1,kk]=size(xx);
test=0;
for h=1:kk
    for i=1:k
        dd=(xx(h)-x(i));
        if dd==0
        yy(h)=y(i);
        test=1;
        end
    end
    if test==0
        R=zeros(k);
        R(:,1)=y(:);
        for i=1:k-1
            D=R(i+1,1)-R(i,1);
            rr=(xx(h)-x(i))/(xx(h)-x(i+1));
            denom=rr*(1-D/R(i+1,1))-1;
            R(i,2)=R(i+1,1)+D/denom;
        end
        for j=3:k
            for i=1:k-j+1
                D=R(i+1,j-1)-R(i,j-1);
                rr=(xx(h)-x(i))/(xx(h)-x(i+j-1));
                if D==0
                    R(i,j)=R(i+1,j-1);
                else
                    DD=R(i+1,j-1)-R(i+1,j-2);
                    denom=rr*(DD-D)-DD;
                    R(i,j)=R(i+1,j-1)+D*DD/denom;
                end
            end
        end
        yy(h)=R(1,k);
    end
    test=0;
end