% illustrates use of piecewise quadratic lagrange 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+2)
	phi(j,i) = (xi(j+1) - x(i))*(xi(j+2) - x(i))/ ...
	    (( xi(j+1)-xi(j))*( xi(j+2)-xi(j)));
      else
	phi(j,i) = 0;
      end
    end
  elseif j == 2
    for i = 1:l
      if xi(j-1) <= x(i) & x(i) <= xi(j+1)
	phi(j,i) = (x(i) - xi(j-1))*(xi(j+1) - x(i))/ ...
	    (( xi(j)-xi(j-1))*( xi(j+1)-xi(j)));
      else
	phi(j,i) = 0;
      end
    end
  elseif j == li-1
 for i = 1:l
      if mod(j,2)==1 % j odd
	if xi(j-2) <= x(i) & x(i) <= xi(j)
	  phi(j,i) = (x(i) - xi(j-1))*(x(i) - xi(j-2))/ ...
	      (( xi(j)-xi(j-1))*( xi(j)-xi(j-2)));
	elseif x(i) >=xi(j)
	  phi(j,i) = (xi(j+1) - x(i))*(xi(j+2) - x(i))/ ...
	      (( xi(j+1)-xi(j))*( xi(j+2)-xi(j)));
	else
	  phi(j,i) = 0;
	end
      else % j even
	if xi(j-1) <= x(i) & x(i) <= xi(j+1)
	  phi(j,i) = (x(i) - xi(j-1))*(xi(j+1) - x(i))/ ...
	      (( xi(j)-xi(j-1))*( xi(j+1)-xi(j)));
	else
	  phi(j,i) = 0;
	end
      end  % j
    end    
  elseif j == li
     for i = 1:l
      if mod(j,2)==1 % j odd
	if xi(j-2) <= x(i) & x(i) <= xi(j)
	  phi(j,i) = (x(i) - xi(j-1))*(x(i) - xi(j-2))/ ...
	      (( xi(j)-xi(j-1))*( xi(j)-xi(j-2)));
	else
	  phi(j,i) = 0;
	end
      else % j even
	if xi(j-1) <= x(i)
	  phi(j,i) = (x(i) - xi(j-1))*(xi(j+1) - x(i))/ ...
	      (( xi(j)-xi(j-1))*( xi(j+1)-xi(j)));
	else
	  phi(j,i) = 0;
	end
      end  % j
    end
  else
    for i = 1:l
      if mod(j,2)==1 % j odd
	if xi(j-2) <= x(i) & x(i) <= xi(j)
	  phi(j,i) = (x(i) - xi(j-1))*(x(i) - xi(j-2))/ ...
	      (( xi(j)-xi(j-1))*( xi(j)-xi(j-2)));
	elseif x(i) >=xi(j) & x(i) < xi(j+2)
	  phi(j,i) = (xi(j+1) - x(i))*(xi(j+2) - x(i))/ ...
	      (( xi(j+1)-xi(j))*( xi(j+2)-xi(j)));
	else
	  phi(j,i) = 0;
	end
      else % j even
	if xi(j-1) <= x(i) & x(i) <= xi(j+1)
	  phi(j,i) = (x(i) - xi(j-1))*(xi(j+1) - x(i))/ ...
	      (( xi(j)-xi(j-1))*( xi(j+1)-xi(j)));
	else
	  phi(j,i) = 0;
	end
      end  % j
    end
  end      
end
cc=sin(pi*xi);
for i = 1:l
  v(i) = 0;
  for j = 1:li
    v(i) = v(i) + cc(1,j)*phi(j,i);
  end
end
zz=zeros(1,l);
oo=ones(1,l);
hold on
for i=1:li
plot(x,phi(i,:))
end
xlabel('X')
title('piecewise quadratic 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 quadratic interpolation')
pause
close