%
% -----------------------------------------------------
%
%   Approximation  im  Wintersemester 2013/14
%   H.J. Oberle
%
%   Aufgabe 21:
%   Interpolation durch C^1-Splines nach Ellis, McLain                  
%
% -----------------------------------------------------
%
%
function aufg21
clear all
close all
%
for iteil=1:4
%
if iteil == 1 
% 
   figure
   x  =  [-5 -4 -3 -2 -1 0 1 2 3 4 5];
   y  =  [0 0 0 0 0 1 0 0 0 0 0];
   n  =  length(x);
   n1 =  n-1;
%
   plot(x,y,'ko')
   set(gca,'FontSize',12)
   axis([-5 5 -0.2 1.1])
   hold
%
   [b,c,d] = ellis(x,y);
%
%  Kubischer C^2-Spline (not a knot Randbedingungen)
%
   x1  =[-5:0.01:5];
   y1  =  spline(x,y,x1);
   plot(x1,y1,'b--')
   set(gca,'FontSize',12)
   xlabel('Aufgabe 21 a)','FontSize',14)
%
   [x2,y2] =  eval(x,y,b,c,d,50);
%
   plot(x2,y2,'b')
   hold off
%
   print -depsc aufg21a
end
%   --------------------------------------------------
if iteil == 2 
%
   figure
   x = [0 0.1 0.2 0.3 4 5 6 7 8 9 10];
   y = [-1 1 -1 1 -1 1 -1 1 -1 1 -1];
   n  =  length(x);
   n1 =  n-1;
%
   plot(x,y,'ko')
   set(gca,'FontSize',12)
   axis([0 10 -2 15])
   hold
%
   [b,c,d] = ellis(x,y);
%
%  Kubischer C^2-Spline
%
   x1  =[0:0.01:10];
   y1  =  spline(x,y,x1);
   plot(x1,y1,'b--')
   set(gca,'FontSize',12)
   xlabel('Aufgabe 21 b)','FontSize',14)
%
   [x2,y2] =  eval(x,y,b,c,d,50);
%
   plot(x2,y2,'b')
   hold off
%
   print -depsc aufg21b
end
%   --------------------------------------------------
if iteil == 3 
%
   figure
   x = [-1 -0.8 -0.55 -0.5 -0.25 0 0.25 0.5 0.55 0.8 1];
   y = [0 0 0 1 1 1 1 1 0 0 0];
   n  =  length(x);
   n1 =  n-1;
%
   plot(x,y,'ko')
   set(gca,'FontSize',12)
   axis([-1 1 -1 2])
   hold
%
   [b,c,d] = ellis(x,y);
%
%  Kubischer C^2-Spline
%
   x1  =[-1:0.01:1];
   y1  =  spline(x,y,x1);
   plot(x1,y1,'r--')
   set(gca,'FontSize',12)
   xlabel('Aufgabe 21 c)','FontSize',14)
%
   [x2,y2] =  eval(x,y,b,c,d,50);
%
   plot(x2,y2,'b')
   hold off
%
   print -depsc aufg21c
end
%   --------------------------------------------------
if iteil == 4 
%
   figure
   x = [-5:1:5];
   y = [0 0 0 0 0 1 1 1 1 1 1];
   n  =  length(x);
   n1 =  n-1;
%
   plot(x,y,'ko')
   set(gca,'FontSize',12)
   axis([-5 5 -0.2 1.2])
   hold
%
   [b,c,d] = ellis(x,y);
%
%  Kubischer C^2-Spline mit den berechneten Randsteigungen
%
   x1  =[-5:0.05:5];
   y1  =  spline(x,[b(1) y b(n)],x1);
   plot(x1,y1,'r--')
   set(gca,'FontSize',12)
   xlabel('Aufgabe 21 d)','FontSize',14)
%
   [x2,y2] =  eval(x,y,b,c,d,50);
%
   plot(x2,y2,'b')
   hold off
%
   print -depsc aufg21d
end
%
end
%
% -------------------------------------------------------------------
%
function [b,c,d] = ellis(x,y);
%
   n    =  length(x);
   n1   =  n - 1;
   n2   =  n - 2;
   n3   =  n - 3;
   h    =  zeros(n,1);
   b    =  zeros(n,1);
   c    =  zeros(n,1);
   d    =  zeros(n,1);
%
   for i=1:n1,
      h(i) = x(i+1)-x(i);
      d(i) = (y(i+1)-y(i))/h(i);
   end
%
   y1    = y(1);
   y12   = d(1);
   y123  = (d(2)-d(1))/(x(3)-x(1));
   y234  = (d(3)-d(2))/(x(4)-x(2));
   y1234 = (y234 - y123)/(x(4)-x(1));
%
   yn    =  y(n);
   yn1n  =  d(n1);
   yn2n  =  (d(n1)-d(n2))/(x(n)-x(n2));
   yn3n1 =  (d(n2)-d(n2-1))/(x(n1)-x(n2-1));  
   yn3n  =  (yn2n-yn3n1)/(x(n)-x(n2-1));
%
   b(1)  =  y12 - y123*h(1) + y1234*h(1)*(h(1)+h(2));
   b(2)  =  y12 + y123*h(1) - y1234*h(1)*h(2);
   b(n)  =  yn1n + yn2n*h(n1) + yn3n*h(n1)*(h(n2)+h(n1));
   b(n1) =  yn1n - yn2n*h(n1) - yn3n*h(n1)*h(n2); 
%
%  Innere Steigungen nach Ellis, McLain
%
   for i=3:n2
      i1 = i - 1;
      ai = (h(i)*d(i1)+h(i1)*d(i))/(h(i1)+h(i));
      bi = (d(i)-d(i1))/(h(i1)+h(i));
      f1 = y(i-2) - y(i) + (x(i)-x(i-2))*(ai-bi*(x(i)-x(i-2))); 
      f2 = y(i+2) - y(i) - (x(i+2)-x(i))*(ai+bi*(x(i+2)-x(i)));
      w1 = h(i-2)*(h(i1)+h(i-2))*(x(i-2)-x(i+1));  
      w2 = (x(i+2)-x(i1))*(h(i)+h(i+1))*h(i+1);  
      ci = (w1*f1+w2*f2)/(w1*w1+w2*w2);
      b(i) = ai - ci*h(i1)*h(i); 
end
%
%   Berechnung der restlichen Spline-Koeffizienten nach (5.5)
%
   for i=1:n1
      i1   = i + 1;
      c(i) = (3*d(i) - 2*b(i) - b(i1))/h(i);
      d(i) = (b(i) + b(i1) - 2*d(i))/(h(i)*h(i));
   end
%
   return;
% -------------------------------------------------------------------
%
%
%
   function [t,s] = eval(x,y,b,c,d,m);
%
   t    =  [ ];
   s    =  [ ];
   n    =  length(x);
   n1   =  n - 1;
   for i=1:n1,
       h   =  (x(i+1) - x(i))/m;
       ti  =  [x(i):h:x(i+1)];
       t1  =  ti - x(i);
       si  =  y(i) + t1.*(b(i)+t1.*(c(i)+t1*d(i)));
       t   =  [t ti];
       s   =  [s si];
   end
%
return;
% -------------------------------------------------------------------





