%=============================================================================
%(0) Vorbemerkungen
%Alle B-Splines k-ter Ordnung (Grad k-1) werden zu vorgegebener Knotenfolge
%Knoten (ein Vektor beliebiger L\"ange, aufsteigend geordnet, mehrfache Knoten 
%erlaubt), berechnet. Wir vergr\"o\ss ern (oder verkleinern) ggf. die
%Knotenfolge an beiden Enden so, da\ss\ die beiden Endknoten die gleichen 
%Vielfachheiten k=Ordnung haben. Beginnt man also nur mit zwei Knoten, 
%so erh\"alt man die Bernstein-Polynome. Ist i die Anzahl der gegebenen, 
%inneren Knoten (mehrfache mehrfach gez\"ahlt), so ist die Gesamtanzahl 
%der verwendeten Knoten lenge:=i+2k. Jeder Spline braucht k+1 Knoten. 
%Das Ergebnis ist eine Matrix BB(j,s). Mit j werden die ausgerechneten 
%B-Splines der Ordnung k numeriert, s steht f\"ur die Anzahl der angegebenen 
%Funktionswerte x. Das Programm hat die Form
%
%             function BB=bsplinesfct(Ordnung,Knoten,x);
%
%Die Größe x enthält die Argumente, für die die B-Splines ausgerechnet 
%werden sollen. Sind die Knoten nicht geordnet, werden sie sortiert, 
%ist die Anzahl der Knoten zu klein, wird die Rechnung abgebrochen.
%Die individuellen Splines sind (mindestens) rechtsseitig stetig, der letzte 
%Spline (mit Nummer lenge-Ordnung) ist am rechten Rand auch linksseitig stetig.
%Der erste, bei a:=min Knoten beginnende Spline B(1,:) mit Tr\"ager 
%[a,Knoten(k+1)] verschwindet nicht, da a vom Programm her die Vielfachheit
%k bekommt und a<Knoten(k+1). Das kann aber ausgeschaltet werden. 
%Man vergleiche weiter unten unter Punkt (2).  
%Programm: Gerhard Opfer, Juli 2000.
%Quelle: Carl de Boor: B(asic)-Spline Basics, ftp-server Madison, 1996.

function BB=bsplinesfct(Ordnung,Knoten,x);
%=============================================================================
%(1) Vorbesetzungen
t=Knoten;
if length(t)<2, error('Zu wenig Knoten'); end; %Schlu\ss 
talt=t;
t=sort(t);
if talt  ~= t
  disp('Knoten nicht aufsteigend geordnet, Umordnung vorgenommen!');
end; %if
talt=[];
%=============================================================================
%(2) Randknoten auf jeder Seite auf genau Ordnung St\"uck erg\"anzen oder 
%reduzieren. Hinweis: Dieser Teil kann in gewissen Fällen unerwünscht sein. 
%Man kann ihn problemlos weglassen (auskommentieren).
li=t(1);re=t(length(t));
links=find(t==li);
rechts=find(t==re);
t([links,rechts])=[];
rnd=ones(1,Ordnung);
t=[li*rnd,t,re*rnd];
%Randknoten angepa\ss t.
%=============================================================================
%(3) Vorbesetzungen
lenge=length(t); %lenge = Anzahl der Knoten
l=length(x);     %l = Anzahl der Funktionswerte
tt=x;
B=zeros(lenge-1,l); 
n=zeros(1,lenge);
n(1)=1;
%Vorbesetzungen zu Ende
%=============================================================================
%tic;  %Die Stoppuhr wird gestartet.
%=============================================================================
%(4) Zuerst st\"uckweise konstante Splines, das ist der komplizierteste Teil.
for j=1:lenge-1
  n(j+1)=n(j);
  if t(j+1)>t(j)
    if j<lenge-Ordnung
      ax=find(x>=t(j) & x<t(j+1));  %Anzahl der x-Werte in [t_j,t_{j+1}[
    elseif j==lenge-Ordnung
      ax=find(x>=t(j) & x<=t(j+1)); %Anzahl der x-Werte in
                          %[t_{lenge-Ordnung},t_{lenge-Ordnung+1}]
                          %linkssseitige Stetigkeit am rechten Intervallpunkt
    end; % if
    n(j+1)=n(j)+length(ax);
    J=n(j):n(j+1)-1;
    if isempty(ax)==0  %nicht leer
      B(j,J)=ones(1,length(J));
    end; %if is
  end; %if t(...
end; %for j
%Stückweise konstante Splines B_1,...,B_{lenge-1} berechnet.
%=============================================================================
%(5) Globalschleife f\"ur alle Ordnungen
for k=2:Ordnung
  %===========================================================================
  %(5a) Bestimmung der Omegas
  for j=1:lenge-k+1
    if t(j+k-1)>t(j)
      omega(j,:)=(tt-t(j))/(t(j+k-1)-t(j));
      else
      omega(j,:)=zeros(1,l);
    end; %if t(
  end; %for j
    %===========================================================================
  %(5b) Bestimmen der Splines: Rekursionsformel
  for j=1:lenge-k
    B(j,:)=omega(j,:).*B(j,:)+(1-omega(j+1,:)).*B(j+1,:);
  end; %for j
%=============================================================================
end; %for k
%Globalschleife f\"ur alle Ordnungen zu Ende
%=============================================================================
%toc;  %Die Stoppuhr wird gestoppt.
%=============================================================================
%(6) Weglassen überflüssiger Teile
%Da zuerst lenge - 1  stückweise konstante Splines hergestellt werden, werden 
%insgesamt lenge-1 B-Splines berechnet. Tatsächlich werden am Schluß aber nur 
%lenge-Ordnung B-Splines der gegebenen Ordnung berechnet. Die auch berechneten 
%B-Splines der Ordnungen 1,2,...,Ordnung - 1 werden jetzt weggelassen.
B(lenge-Ordnung+1:lenge-1,:)=[];
%=============================================================================
%(7) Schlus\ss
BB=B;
return
%=============================================================================
