function zp = zkpol(fctid,fctalp,xmin,xmax,aprmin,aprmax,ndeg,npts)
%
% fctid         = ID of interpolated function
% fctalp        = parameter alpha in interpolated function
% xmin,xmax     = interpolation interval (where we construct polynomials)
% aprmin,aprmax = evaluation interval (if > xmin,xmax, we axtrapolate)
% ndeg(1:3)     = degrees of interpolation polynomials
% npts          = number of points for plotting
%
%
%
% Properties of the polynomial interpolation
%
% Procedure POLINT from Numerical Recipes
%
% Ndeg1,ndeg2 or ndeg3 equidistant points between xmin and xmax are generated
%
% Y given by some function f = f(x;alpha)
%
% Interpolated and extrapolated values from aprmin to aprmax
% are calculated and compared to exact values

% fctid=1; fctalp=4.0;  xmin=0.0;  xmax=6.0;  npts = 30;  ndeg = [2 4 10];

[s1 s2] = size(ndeg);
if ((s1~=1)||(s2~=3))
    disp(' ndeg (6th argument) has to be a 1x3 array !!!');
    return;
end
    

 ypol = zeros(3,npts);
 %data = zeros(5,npts);
 xinp = zeros(3,max(ndeg(:))+1);
 yinp = xinp;

 % aprmin = 2*xmin - xmax;
 % aprmax = 2*xmax - xmin;
 % points where to evaluate (for output)
 dx = (aprmax - aprmin)/(npts-1);
 x  = (aprmin:dx:aprmax);
    
 % polynomial approximation
 for k=1:3
     dx = (xmax-xmin)/ndeg(k);
     xval = (xmin:dx:xmax);
     % evaluate the function for each value of vector xval
     yval = funval( xval , fctid , fctalp );

     % save the points
     xinp(k,1:ndeg(k)+1) = xval;
     yinp(k,1:ndeg(k)+1) = yval;

     % calculate the interpolation
     for i=1:npts
         [ypol(k,i),dy] = polint(xval,yval,ndeg(k)+1,x(i),0);
     end;
 end;

 % exact value (for output)
 fx = funval(x , fctid, fctalp);
 
 %% output to file and screen
 
 fileID = fopen('zkpol.dat','w');

 for i=1:npts
     fprintf(fileID,' %15.8e  %15.8e  %15.8e  %15.8e  %15.8e \n',...
                  x(i), fx(i), ypol(1:3,i));
     % fprintf(' %15.8e  %15.8e  %15.8e  %15.8e  %15.8e \n', x(i), fx(i), ypol(1:3,i));
     % data(i,1:5) = [x(i) fx(i) ypol(i,1:3)];
 end;
     
 fclose(fileID);

 %% read data from file 

 % fileID = fopen('zkpol.dat');
 % data = fscanf(fileID,'%g %g %g %g %g',[5 inf]);
 % fclose(fileID);
 % [bf,npts]=size(data);
 %
 % x  = data(:,1)
 % fx = data(:,2);
 % ypol(:,1)=data(:,3); ypol(:,2)=data(:,4); ypol(:,3)=data(:,5);
 %
 % xmin = x(1) +   (x(npts)-x(1))/3;
 % xmax = x(1) + 2*(x(npts)-x(1))/3;
 
 %% plot the graphs

 fig=figure(1);
 clf(fig,'reset')
 
 ymin = min(fx(:));  ymax = max(fx(:));
 bdy  = (ymax-ymin)/10;
 axis([x(1) x(npts) ymin-bdy ymax+bdy])
 hold on;

 % plot the function
 plot(x,fx,'k-');

 % plot the approximations
 plot(x,ypol(1,:),'mo');
 plot(x,ypol(2,:),'b+');
 plot(x,ypol(3,:),'rx');
 
 % plot lines at xmin and xmax
 plot([xmin xmin],[ymin-bdy ymax+bdy],'k--');
 plot([xmax xmax],[ymin-bdy ymax+bdy],'k--');
 
 leg = legend({'exact',['N = ' num2str(ndeg(1))], ...
                       ['N = ' num2str(ndeg(2))], ... 
                       ['N = ' num2str(ndeg(3))]},'Location','NorthWest');
 title(['xmin='  num2str(xmin) ', xmax= ' num2str(xmax) ', approx in ' ...
        num2str(npts) ' points']);

 hold off;

 %% plot the graphs separately long with input points
 fig=figure(2);
 clf(fig,'reset')
 
 col = ['m' 'b' 'r'];
 sty = ['o' '+' 'x'];

 
 for k=1:3
   subplot(3,1,k)
   axis([x(1) x(npts) ymin-bdy ymax+bdy])
   hold on;

   % plot the function
   plx = plot(x,fx,'k-');
   
   % plot the generating points
   plg = plot(xinp(k,1:ndeg(k)+1),yinp(k,1:ndeg(k)+1),[col(k) 'o'],'markerfacecolor',col(k));
   
   % plot the approximation
   pli = plot(x,ypol(k,:),[col(k) sty(k)]);
   
   % plot lines at xmin and xmax
   plot([xmin xmin],[ymin-bdy ymax+bdy],'k--');
   plot([xmax xmax],[ymin-bdy ymax+bdy],'k--');
   
   leg = legend([plx,plg,pli],{'exact','input','interp'},'Location','NorthWest');
   title(['N = ' num2str(ndeg(k))]);
   
   hold off;
 end
 
 
 return
    
%%%%%%%%%%%%

function fv = funval(x,id,alp);
%
% The function to be approximated. 
%
% Works for scalar as well as vector arguments
%
    if (id==1)
        fv = exp(x/alp);

    elseif (id==2)
        fv = sqrt( max(x,0.0) / alp );
    
    elseif (id==3)
        fv = 1 ./ (1+alp*x.^2) ;
    
    end;
    
    return;
    
%%%%%%%%%%%%

