function r = demo_spline_interpolation(nn,tol)

% demo_spline_interpolation berechnet die Spline-Interpolation f"ur
% Splines 0., 1. und 3. Ordnung sowie das Newton-Interpolationspolynom
% mit Hilfe von dividierten Differenzen.
% Aufruf:
%   r = demo_spline_interpolation(nn,tol)
% nn - Anzahl der St"utzstellen
% tol - Schwankung der St"utzwerte um 0.5 (+- tol/2)
% x : St"utzstellen
% f : St"utzwerte
% Die Matrix der dividierten Differenzen wird im Feld d gespeichert.
% f"ur n>50 sehen die Bilder des Interpolationspolynoms nicht mehr gut aus

clf

for i=1:nn+1
    x(i) = (i-1)/nn; 
end
f = 0.5+tol*(rand(1,nn+1)-0.5)

if abs(tol) < 1e-4
     tol = 1e-4
end

% berechne Dimension der Vektoren
[n,m] = size(x);
[nf,mf] = size(f);

% konvertiere, falls n"otig, die Vektoren zu Spaltenvektoren
if n==1
  x = x';
  n = m;
end

if nf==1
  f = f';
  nf = mf;
end

% kontrolliere, ob gleiche Anzahl von Werten vorhanden ist
if n ~= nf
  error('Die Dimensionen von x und f passen nicht')
end

% erniedrige die Dimension, damit die Formeln die gleiche Gestalt
% wie in der Vorlesung haben (Vorlesung: Indizes von 0 bis n;
% hier, Indizes von 1 bis n)
n = n-1;

% berechne Dimensionen f"ur das Bild
x_min = 1e10;
x_max = -1e10;
f_min = 1e10;
f_max = -1e10;

for i=1:n+1;
  if (x(i) < x_min)
    x_min = x(i);
  end
  if (x(i) > x_max)
    x_max = x(i);
  end
  if (f(i) < f_min)
    f_min = f(i);
  end
  if (f(i) > f_max)
    f_max = f(i);
  end
end

%%%%%%%%%%%%%%%%%%%%%% Newton-Interpolation %%%%%%%%%%%%%%%%%%
% erzeuge die Matrix f"ur die dividierten Differenzen
close all
d = zeros(n+1,n+1);

% belege erste Spalte mit der rechten Seite
d(:,1) = f;

% berechne die dividierten Differnzen
for k=2:n+1
  for i=k:n+1
    d(i,k) = (d(i,k-1)-d(i-1,k-1))/(x(i)-x(i-k+1));
  end
end

% belege den Vektor a mit der Hauptdiagonalen von d
a = diag(d);

% zeichne Bild des Interpolationspolynoms
 close all

% berechne Dimensionen f"ur das Bild
x_min = 1e10;
x_max = -1e10;
f_min = 1e10;
f_max = -1e10;

for i=1:n+1;
  if (x(i) < x_min)
    x_min = x(i);
  end
  if (x(i) > x_max)
    x_max = x(i);
  end
  if (f(i) < f_min)
    f_min = f(i);
  end
  if (f(i) > f_max)
    f_max = f(i);
  end
end

plot(x,f,'ok');
pause()

hold off
clf

% berechne Werte f"ur das Interpolationspolynom
steps = 1000;
h = (x_max-x_min+2)/steps;

for i=1:steps+1;
  % Abszisse
  xp(i) = x_min-1 + (i-1)*h;
  % Verwende Formel (7) aus der Vorlesung
  % Initialisiere Funktionswert und Knotenpolynom
  fp(i) = 0;
  omega = 1;
  % berechne Funktionswert
  for j = 1:n+1
    fp(i) = fp(i) + a(j)*omega;
    omega = omega * (xp(i)-x(j));
  end
end


plot(xp(335:667),fp(335:667),'-k',x,f,'ok');
%plot(xp,fp,'-k',x,f,'ok');
pause();

hold off;
clf;


plot(xp,fp,'-k',x,f,'ok');
axis([x_min-0.1,x_max+0.1,f_min-0.1,f_max+0.1])
pause()
hold on;

%%%%%%%%%%%%%%%%%%% lineare Splines %%%%%%%%%%%%%%%%%%%%%%
%close all

% zeichne die St"utzstellen
plot(x,f,'o');
axis([x_min-0.1,x_max+0.1,f_min-0.1,f_max+0.1])
hold on;

% zeichne den st"uckweise konstanten Spline

for i=1:n
  vx = [x(i),x(i+1)];
  vy = [f(i),f(i+1)];
  poly0  = line(vx,vy);
  set(poly0,'linestyle','-')
  hold on;
end

pause();
clf
% zeichne die St"utzstellen
plot(x,f,'o');
axis([x_min-0.1,x_max+0.1,f_min-tol/2,f_max+tol/2])
hold on;

% zeichne den st"uckweise konstanten Spline

for i=1:n
  vx = [x(i),x(i+1)];
  vy = [f(i),f(i+1)];
  poly0  = line(vx,vy);
  set(poly0,'linestyle','-')
  hold on;
end

pause();

%%%%%%%%%%%%%%%%%%% kubische Splines %%%%%%%%%%%%%%%%%%%%%%
%close all

% zeichne die St"utzstellen
plot(x,f,'o');
axis([x_min-0.1,x_max+0.1,f_min-tol/2,f_max+tol/2])
hold on;

% definiere Felder
A=zeros(4*n,4*n);
rhs=zeros(4*n,1);

% fill matrix with the conditions for the 
% cubic splines

% linker Rand
% Funktionswert
A(1,1) = 1;
rhs(1) = f(1);
% erste und zweite Ableitung gleich Null
% A(2,2) = 1;
A(2,2) = 2;

% innere St"utzstellen
for j=1:n-1
  A(3*j,4*j-3) = 1;
  A(3*j,4*j-2) = x(j+1) - x(j);
  A(3*j,4*j-1) = (x(j+1) - x(j))^2;
  A(3*j,4*j) = (x(j+1) - x(j))^3;
  A(3*j,4*j+1) = -1;
  A(3*j+1,4*j-2) = 1;
  A(3*j+1,4*j-1) = 2*(x(j+1) - x(j));
  A(3*j+1,4*j) = 3*(x(j+1) - x(j))^2;
  A(3*j+1,4*j+2) = -1;
  A(3*j+2,4*j-1) = 2;
  A(3*j+2,4*j) = 6*(x(j+1) - x(j));
  A(3*j+2,4*j+3) = -2;  
end

% rechter Rand
% Funktionswert
A(3*n,4*n-3) = 1;
A(3*n,4*n-2) = x(n+1)-x(n);
A(3*n,4*n-1) = (x(n+1)-x(n))^2;
A(3*n,4*n) = (x(n+1)-x(n))^3;
rhs(3*n) = f(n+1);
% erste und zweite Ableitung gleich Null 
%A(3*n+2,4*n-2) = 1;
%A(3*n+2,4*n-1) = 2*(x(n+1)-x(n));
%A(3*n+2,4*n) = 3*(x(n+1)-x(n))^2;
A(3*n+1,4*n-1) = 2;
A(3*n+1,4*n) = 6*(x(n+1)-x(n));

% Gleichungen f"ur die Werte in den inneren St"utzstellen
for i=2:n
  A(3*n+i,4*i-3) = 1;
  rhs(3*n+i) = f(i);
end
    
% l"ose das lineare System 
sol = A\rhs;

% plot computed reconstruction
step = 10;

for j=1:n
  h = (x(j+1) - x(j))/step;
  for i=1:step+1
    xx(i) =  x(j)+(i-1)*h;
    fx(i) =  sol(4*j-3) +  sol(4*j-2)* (xx(i) - x(j)) + ...
             sol(4*j-1)* (xx(i) - x(j))^2 + sol(4*j)* (xx(i) - x(j))^3;
  end 
            
  % plot all functions  
  plot(xx(:),fx(:),'r-');
  hold on
end


