%======================Vorspann und Kommentar=================================
%Zu gegebener Knotenverteilung Knoten und f\"ur gegebene Werte X werden 
%alle B-Splines der Ordnung k rekursiv ausgerechnet. Das Ergebnis ist 
%eine Matrix B mit length(Knoten)-k Zeilen und mit length(X) Spalten. 
%Die Spalten entsprechen den Werten von X, die Zeilen enthalten die
%verschiedenen B-Splines.

%Wird zum Zeichnen verwendet in: bsplines.m
%Wird zum Interpolieren und Zeichnen verwendet in: bsplinesinterpol.m

%Aufruf:                [B]=l06_deBoor(k,Knoten,X)

%Dies ist i.w. der Algorithmus von de Boor, S. 134,
%geschneidert nach der bekannten Rekursionsformel f\"ur B-Splines.
%Das Ergebnis ist f\"ur jeden reellen Wert x ein k-Vektor b, 
%in dem die Werte aller k B-Splines an der Stelle x, immer zugeordnet  
%zum Knoten Knoten(i), gespeichert werden. Also
%b(1)=B_{i-k+1,k}(x), b(2)=B_{i-k+2,k}(x), ..., b(k)=B_{i,k}(x).
%Um die B-Splines zu zeichnen, w\"ahle man geeignet viele x.
%Will man Interpolationsprobleme (i.e. Gleichungssysteme) l\"osen, so mu\ss\ der
%Vektor Knoten die L\"ange n+k und X die L\"ange n haben. Die resultierende 
%Matrix B ist dann quadratisch mit Ordnung n.

%Man vgl. auch den Kommentar in l06_indexlauf.doc
%Programm: Gerhard Opfer, 13. Mai 1999

%===================================Kopf=====================================
function [B]=deBoor(k,Knoten,X)

%=================================Vorspann===================================
n=length(Knoten);dis=length(X);
B=zeros(n+k-2,dis); %In B werden die B-Splines gesammelt.
%Test 1, ob gen\"ugend Knoten vorhanden sind.
if n<=k 
  error(['Ordnung = ',int2str(k),', Anzahl der Knoten = ',int2str(n),...
         ', zu wenig  Knoten!']);
end; %if
%Test 2, ob die Vielfachheit der Knoten nicht \"uberschritten wird.
for j=1:n-k
  if Knoten(j)>=Knoten(j+k)  
    error(['Knoten ',int2str(j),' hat zu hohe Vielfachheit']);
  end; %if
end; %for j

%==============================Vorspann Ende=================================

%================================B-Splines===================================
sp=0;
for x=X;
  sp=sp+1; %wird unten gebraucht.
  b=zeros(k,1); %Platz f\"ur die B-Splines f\"ur ein einziges x
  deltaR=b;deltaL=b;
  %x (ungleich letzter Knoten) definiert eindeutig ein Intervall 
  %x in   [Knoten(i),Knoten(i+1)[.
  %Ist x der letzte Knoten, so ist i=n-Vielfachheit des letzten Knotens.
  ih=find(x>=Knoten);
  i=max(ih);
  if i==n, %wenn x der letzte Knoten ist
    ih=find(Knoten==x);
    i=n-length(ih);
  end; %if
  b(1)=1;
  for j=1:k-1
    if i+j>n
      deltaR(j)=-x; else %Sonderma\ss nahme f\"ur nicht vorhandene Knoten
      deltaR(j)=Knoten(i+j)-x;
    end; %if
    if i+1-j<=0
      deltaL(j)=x;  else %Sonderma\ss nahme f\"ur nicht vorhandene Knoten
      deltaL(j)=x-Knoten(i+1-j);
    end; %if
    jp1=j+1;
    saved=0;
    for r=1:j
      term=b(r)/(deltaR(r)+deltaL(jp1-r));
      b(r)=saved+deltaR(r)*term;
      saved=deltaL(jp1-r)*term;
    end; %for r
    b(jp1)=saved;
  end; %for j
  B(i:i+k-1,sp)=b;   
end; %for x
B([1:k-1,n:n+k-2],:)=[]; %\"Uberschu\ss\ wird gestrichen.
%Die Matrix B hat schlie\ss lich n-k Zeilen und dis Spalten.
 
%=======================Ende B-Splines erzeugen==============================
