function func = plotpint ( ifcn, iknot, n );
% PlotPInt
%
% This m-file is part of the material of the lectures on
%   Numerical Mathematics
% at the Math Department of Halle University, Germany.
%
% Author :      Prof. Dr. M. Arnold, arnold@mathematik.uni-halle.de
% Version of :  Apr 5, 2005
%
% Plot of test results on Polynomial Interpolation
%
% Parameters:
%   ifcn   (input)  : control flag "function"
%                       =1 ... sin(x)
%                       =2 ... 1/(1+x^2)  [Runge]
%   iknot  (input)  : control flag "knots"
%                       =1 ... equidistant
%                       =2 ... Tschebyscheff
%   n      (input)  : determines the number n+1 of data points
%
% Examples:
%   plotpint ( 1, 1, 2 );   print -dpng ../png/plotpint1102
%   plotpint ( 1, 1, 6 );   print -dpng ../png/plotpint1106
%   plotpint ( 2, 1, 2 );   print -dpng ../png/plotpint2102
%   plotpint ( 2, 1, 6 );   print -dpng ../png/plotpint2106
%   plotpint ( 2, 1, 8 );   print -dpng ../png/plotpint2108
%   plotpint ( 2, 2, 6 );   print -dpng ../png/plotpint2206
%   plotpint ( 2, 2, 8 );   print -dpng ../png/plotpint2208

graphwidth  = 24.0;
graphheight = 12.0;
fac         =  0.85;
shift       =  0.6;
papershift  =  0.0;

ActFontSize   = 14;
ActLineWidth  =  2;
ActMarkerSize =  8;
RatioFlag     =  0;

border      = 1.0;

%
% -> Bestimme Groesse des Bildschirms und Skalierungsfaktor fuer die Ausgabe
root = 0;
set ( root, 'Units', 'centimeters' );
screensize = get ( root, 'ScreenSize' );
left   = screensize(1);
bottom = screensize(2);
width  = screensize(3);
height = screensize(4);
set ( root, 'Units', 'pixels' );
if RatioFlag==0,
  screenscale = 1.0;
elseif RatioFlag==1,
  screenscale = (width-2*border) / graphwidth;
end;

%
% -> Definiere den Bildschirmausschnitt
if length ( get ( root, 'Children' ) ) == 0,
  h = figure;
else
  h = gcf;
  clf;
end;
set ( h, 'Units',    'centimeters' );
set ( h, 'Position', ...
  [ left+width-screenscale*graphwidth-border height-2*border-screenscale*graphheight ...
    screenscale*graphwidth                   screenscale*graphheight ] );
set ( h, 'Units',    'pixels' );

%
% -> Skalierung der Ausgabe in PostScript
%
set ( gcf, 'PaperUnits',    'centimeters' );
set ( gcf, 'PaperPosition', [ papershift papershift graphwidth graphheight ] );

% -> control data 
if ifcn==1,
    xmin =  0.0;
    xmax =  1.0;
    ymin = -1.0;
    ymax =  2.0;
    fstr = 'sin(\pi{}x)';
elseif ifcn==2,
    xmin = -5.0;
    xmax =  5.0;
    ymin = -1.0;
    ymax =  2.0;
    fstr = '1/(1+x^2)';
end;

errmin = 1.0e-6;
errmax = 10.0;

% -> data grid
if iknot==1,
    x = linspace ( xmin, xmax, n+1 );
    kstr = sprintf ( 'Knoten: %2i (aequidistant)', n+1 );
elseif iknot==2,
    x = ( xmin + xmax + ( xmax - xmin ) * cos(pi*(2*(1:(n+1))-1)/(2*(n+1))) ) / 2;
    kstr = sprintf ( 'Knoten: %2i (Tschebyscheff)', n+1 );
end;

% -> data points
if ifcn==1,
    f = sin(pi*x);
elseif ifcn==2,
    f = 1 ./ ( 1 + x.^2 );
end;

% -> difference scheme
a = zeros ( n+1, 1 ); 
for i=n:-1:0, 
    a(i+1) = f(i+1); 
    for j=i+1:n, 
        a(j+1) = (a(j+1)-a(j)) / (x(j+1)-x(i+1));
    end;
end;

% -> output data
nxx = 1001;
xx  = linspace ( -5, 5, nxx );
if ifcn==1,
    ff = sin(pi*xx);
elseif ifcn==2,
    ff = 1 ./ ( 1 + xx.^2 );
end;

pp  = zeros ( size(xx) );

for ixx=1:nxx, 
    pp(ixx) = a(n+1); 
    for i=n-1:-1:0, 
        pp(ixx) = pp(ixx) * (xx(ixx)-x(i+1)) + a(i+1);
    end;
end;

% -> plot
subplot ( 1, 2, 1 );
plot ( xx, ff, 'g', xx, pp, 'r', x, f, '*b'  ...
     , 'LineWidth',  ActLineWidth            ...
     , 'MarkerSize', ActMarkerSize           ...
     )
hold on;
plot ( [ xmin xmax ], [ 0 0 ], ':k' );
plot ( [ 0 0 ], [ ymin ymax ], ':k' );
hold off;
axis ( [ xmin xmax ymin ymax ] );
xlabel ( 'x' );
ylabel ( 'f' );
title ( [ 'f(x) = ', fstr, ', ', kstr ] );
subplot ( 1, 2, 2 );
semilogy ( xx, abs(pp-ff), 'b' );
axis ( [ xmin xmax errmin errmax ] );
xlabel ( 'x' );
ylabel ( 'Interpolationsfehler' );
title ( [ 'f(x) = ', fstr, ', ', kstr ] );

%
% -> Korrektur der Groesse und Lage der Koordinatensysteme
allax = get ( gcf, 'Children' )';
for ax = allax,
  rect    = get ( ax, 'Position' );
  rect(2) = rect(2) + shift * (1-fac) * rect(4);
  rect(4) = fac * rect(4);
  set ( ax, 'Position', rect );
end;

%
% -> Aenderung der Schriftgroesse
for hf = allax,
  set ( hf, 'FontSize', [ActFontSize] );
  set ( get ( hf, 'title' ), 'FontSize', [ActFontSize] );
  set ( get ( hf, 'xlabel' ), 'FontSize', [ActFontSize] );
  set ( get ( hf, 'ylabel' ), 'FontSize', [ActFontSize] );
end;

