% illustrates use of piecewise cubic Hermite interpolants 
% in the interpolation of a sine function
% NOTE: THE PROGRAM ASSUMES THAT THERE ARE AN ODD NUMBER OF NODES
clear
x=[0:0.02:3];
l = length(x);
xi=[0 0.2 0.5 1.2 1.6 2 2.4 2.8 3]; % nodes
li = length(xi);

fprintf('HERE ARE THE NODES: \n')
for i=1:li
  fprintf('%3.2f ',xi(i))
end
fprintf('\n')
for j = 1:li
  if j == 1  
    for i = 1:l
      if x(i) >=xi(j) & x(i) < xi(j+1)
	phi(j,i) = (x(i) - xi(j+1)).^2*(2*(x(i) - xi(j))+ ...
					   (xi(j+1)-xi(j)))/ ...
					   ( xi(j+1)-xi(j)).^3;
	psi(j,i) = (x(i) - xi(j+1)).^2*(x(i) - xi(j))/ ...
	    ( xi(j+1)-xi(j)).^2;	
      else
	phi(j,i) = 0;
      end
    end
  elseif j == li
    for i = 1:l
      if xi(j-1) <= x(i) & x(i) <= xi(j)
	phi(j,i) = (x(i) - xi(j-1)).^2*(2*(xi(j) - x(i))+ ...
					   (xi(j)-xi(j-1)))/ ...
					   ( xi(j)-xi(j-1)).^3;
	psi(j,i) = (x(i) - xi(j-1)).^2*(x(i) - xi(j))/ ...
	    ( xi(j)-xi(j-1)).^2;
	
      else
	phi(j,i) = 0;
      end
    end
  else
    for i = 1:l      
      if xi(j-1) <= x(i) & x(i) <= xi(j)
	phi(j,i) = (x(i) - xi(j-1)).^2*(2*(xi(j) - x(i))+ (xi(j)-xi(j-1)))/ ...
					   ( xi(j)-xi(j-1))^3;
	psi(j,i) = (x(i) - xi(j-1)).^2*(x(i) - xi(j))/ ...
					   ( xi(j)-xi(j-1)).^2;
      elseif x(i) >=xi(j) & x(i) < xi(j+1)
	phi(j,i) = (x(i) - xi(j+1)).^2*(2*(x(i) - xi(j))+ (xi(j+1)-xi(j)))/ ...
					   ( xi(j+1)-xi(j))^3;
	psi(j,i) = (x(i) - xi(j+1)).^2*(x(i) - xi(j))/ ...
	    ( xi(j+1)-xi(j)).^2;	
      else
	phi(j,i) = 0;
      end
    end
  end      
end
cc=sin(pi*xi);
dd=pi*cos(pi*xi);
for i = 1:l
  v(i) = 0;
  for j = 1:li
    v(i) = v(i) + cc(1,j)*phi(j,i) + dd(1,j)*psi(j,i);
  end
end
zz=zeros(1,l);
oo=ones(1,l);
hold on
for i=1:li
plot(x,phi(i,:),x,psi(i,:))
end
xlabel('X')
title('cubic Hermite interpolants')
plot(x,oo)
pause
hold
fprintf('Hit any key to continue')
plot(x,v,'-',x,sin(pi*x),x,zz)
xlabel('X')
title('piecewise cubic interpolation')
pause
close