%Alle Splines k. Ordnung (Grad k-1) sollen berechnet und gezeichnet werden. 
%Wegen der Zeichnungen ist das Programm im Moment so gemacht, daß nur ganz- 
%zahlige Knoten zugelassen sind und alle Funktionswerte simultan ausgerechnet 
%werden. Will man nur für einen einzigen Funktionswert
%die B-Splines ausrechnen, könnte das Programm vereinfacht werden.
%Programm: Gerhard Opfer, Juni 2000.
%Quelle: Carl de Boor: B(asic)-Spline Basics, ftp-server Madison, 1996.
%=============================================================================
%Vorbesetzungen
clear
Ordnung=6; %vorgegebene Ordnung
t=[0,2,5,10,14,14,16]; %gegebene Knotenfolge
if length(t)<2, error('Zu wenig Knoten'); end; %error
talt=t;
t=sort(t);
if talt  ~= t
  disp('Knoten nicht aufsteigend geordnet, Umordnung vorgenommen!');
end; %if
talt=[];
if t ~= round(t)
  error('Knoten nicht ganzzahlig'); 
end; %if 
%Wir vergrößern (oder verkleinern) ggf. die Knotenfolge an beiden Enden 
%so, daß die beiden Endknoten die gleichen Vielfachheiten k=Ordnung haben.
%Beginnt man also nur mit zwei Knoten, so erhält man die Bernstein-Polynome.
%Ist i die Anzahl der inneren Knoten (mehrfache mehrfach gezählt), 
%so ist die Gesamtzahl i+2k. Jeder Spline braucht k+1 Knoten
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];
lenge=length(t);
h=100;
l=h*(t(length(t))-t(1))+1;
tt=linspace(t(1),t(length(t)),l);
B=zeros(lenge-1,l); 
n(1)=1;
%=============================================================================
tic;  %Die Stoppuhr wird gestartet.
%=============================================================================
%Zuerst stückweise konstante Splines
for j=1:lenge-1
  n(j+1)=n(j)+(t(j+1)-t(j))*h;
  if t(j+1)>t(j)
    if j==lenge-Ordnung %am rechten Rand linksseitig stetig!
      J=n(j):n(j+1);
      else
      J=n(j):n(j+1)-1;
    end; %if
    B(j,J)=ones(1,length(J));
  end; %if t(...
end; %for j
%=============================================================================
%Globalschleife f\"ur alle Ordnungen
for k=2:Ordnung
  %===========================================================================
  %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
%=============================================================================
%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
%=============================================================================
toc;  %Die Stoppuhr wird gestoppt.
%=============================================================================
%Globalschleife f\"ur alle Ordnungen zu Ende 
%Alle B-Slines sind berechnet, sie sind in einer Matrix B(j,l) zusammengefaßt.
%j=1,2,...,lenge-ordnung; lenge=Anzahl der Knoten, inkl. der zusätzlichen
%Randknoten, l=100*(b-a)+1, dabei werden die Randknoten a, b und alle 
%anderen Knoten als ganzzahlig vorausgesetzt.
%=============================================================================
%Zeichnen der Splines B(j,:), j=1,2,...,lenge-Ordnung 
%=============================================================================
%Zuerst Knoten zeichnen
figure(1)
for j=1:lenge
  plot([t(j),t(j)],[0,0.5],'r');  %Trennstriche zwischen den Knoten
                                  %Die Zahl 0.5 muß evtl angepaßt werden
  if j==1,hold on, end; %if
end; %for j
plot(t,zeros(size(t)),'ro');
%=============================================================================
%Alle (oder einige) Splines der Ordnung k in eine Zeichnung
for j=6%:lenge-Ordnung 
  plot(tt,B(j,:),'k'); %'k' heißt "schwarz"
end; %for j
%=============================================================================
%Ende der Zeichnung
%=============================================================================
%Will man die Polynomstücke zwischen benachbarten Knoten wiederfinden, so ist 
%das am einfachsten, wenn man ein Interpolationsproblem mit polyfit löst.
%Beispiel für den Fall der Aufgabe 17.
Stellen=[1,10,20,40,80,160]; %6 beliebige Werte zwischen Null und zwei
%Diese Stellen entsprechen den Argumenten [1,10,20,40,80,160]/100.
u=tt(Stellen);v=B(6,Stellen);
p=polyfit(u,v,5);
z=linspace(0,6.2);
w=polyval(p,z);
plot(z,w,'--m');   %gestrichelt mit Farbe "magenta"
text(6,0.48,'#1'); %das muß ausprobiert werden
hold off;
