function cubic_spline_demo
%
% This demonstrates the behavior of cubic spline interpolating 
% functions on design data (e.g., for part of a car door) that 
% suggest a discontinuous derivative at a certain point. This  
% demo illustrates the following options for interpolating 
% the data: 
%  (1) a single cubic spline, 
%  (2) two cubic splines joined at the derivative discontinuity, 
%  (3) two cubic splines joined at the derivative discontinuity 
%      and having specified derivative end conditions, 
%  (4) a piecewise cubic Hermite interpolating polynomial 
%      determined by MATLAB's pchip. 
%
% The data are from Example 6.16, p. 291, S. C. Conte and 
% C. de Boor, ``Elementary Numerical Analysis: An Algorithmic 
% Approach'' (3rd ed.), McGraw-Hill, 1980. 
% 

x = zeros(1,17);
y = x;

x(1:6) = 0:1:5;
x(7:10) = 6.1:.2:6.7;
x(11:13) = 7:1:9;
x(14:17) = 12:2:18; 

y(1) = 2.51; y(2) = 3.30; y(3) = 4.04; 
y(4) = 4.70; y(5) = 5.22; y(6) = 5.54; 
y(7) = 5.80; y(8) = 5.55; y(9) = 5.44; 
y(10) = 5.40; y(11) = 5.40; y(12) = 5.57; 
y(13) = 5.70; y(14) = 5.84; y(15) = 5.75; 
y(16) = 5.48; y(17) = 4.90;

figure(1); plot(x,y,'r.','LineWidth',0.5,'MarkerSize',6.0,'MarkerFaceColor','r');
title('Interpolation Data'); 
axis equal; 

disp(' '); 
disp('These are the data to be interpolated'); 
disp(' '); 

while 1 == 1

  ind = input(['\n Choose an option: \n   0 => quit ',...
           '\n   1 => single cubic spline \n   2 => two cubic splines',...
	   '\n   3 => two cubic splines with derivative end conditions',...
	   '\n   4 => piecewise cubic determined by pchip.\n']);

  if ind <= 0, break; end 

  if ind == 1 
    disp(' '); 
    disp('Fit the data with a single cubic spline:'); 
    disp(' '); 

    x = zeros(1,17);
    y = x;

    x(1:6) = 0:1:5;
    x(7:10) = 6.1:.2:6.7;
    x(11:13) = 7:1:9;
    x(14:17) = 12:2:18; 

    y(1) = 2.51; y(2) = 3.30; y(3) = 4.04; 
    y(4) = 4.70; y(5) = 5.22; y(6) = 5.54; 
    y(7) = 5.80; y(8) = 5.55; y(9) = 5.44; 
    y(10) = 5.40; y(11) = 5.40; y(12) = 5.57; 
    y(13) = 5.70; y(14) = 5.84; y(15) = 5.75; 
    y(16) = 5.48; y(17) = 4.90;

    xleft = x(1); xright = x(17);
    h = (xright-xleft)/100;
    xx = [xleft:h:xright];
    yy = spline(x,y,xx); 

    figure(1); plot(x,y,'r.',xx,yy,'-b');
    title('Single Cubic Spline Interpolation');
    axis equal; 
  end 

  if ind == 2
    disp(' '); 
    disp('Fit the data with two cubic splines.'); 
    disp(' '); 

    x = zeros(1,18);
    y = x;

    x(1:6) = 0:1:5; x(7) = 6.1; 

    x(8:11) = 6.1:.2:6.7;
    x(12:14) = 7:1:9;
    x(15:18) = 12:2:18; 

    y(1) = 2.51; y(2) = 3.30; y(3) = 4.04; 
    y(4) = 4.70; y(5) = 5.22; y(6) = 5.54; 
    y(7) = 5.80; 

    y(8) = 5.80; y(9) = 5.55; y(10) = 5.44; 
    y(11) = 5.40; y(12) = 5.40; y(13) = 5.57; 
    y(14) = 5.70; y(15) = 5.84; y(16) = 5.75; 
    y(17) = 5.48; y(18) = 4.90;

    xleft = x(1); xright = x(7);
    h = (xright-xleft)/50;
    xxleft = [xleft:h:xright];
    yyleft = spline(x(1:7),y(1:7),xxleft); 

    xleft = x(8); xright = x(18);
    h = (xright-xleft)/50;
    xxright = [xleft:h:xright];
    yyright = spline(x(8:18),y(8:18),xxright); 

    figure(1); plot(x,y,'r.',xxleft,yyleft,'b-',xxright,yyright,'b-');
    title('Interpolation with Two Cubic Splines.'); 
    axis equal; 
  end 

  if ind == 3
    disp(' '); 
    disp('Fit the data with two splines having specified end derivatives.'); 
    disp(' '); 

    x = zeros(1,18);
    y = x;

    x(1:6) = 0:1:5; x(7) = 6.1; 
    x(8:11) = 6.1:.2:6.7;
    x(12:14) = 7:1:9;
    x(15:18) = 12:2:18; 

    y(1) = 2.51; y(2) = 3.30; y(3) = 4.04; 
    y(4) = 4.70; y(5) = 5.22; y(6) = 5.54; 
    y(7) = 5.80; 

    y(8) = 5.80; y(9) = 5.55; y(10) = 5.44; 
    y(11) = 5.40; y(12) = 5.40; y(13) = 5.57; 
    y(14) = 5.70; y(15) = 5.84; y(16) = 5.75; 
    y(17) = 5.48; y(18) = 4.90;

    xleft = x(1); xright = x(7);
    h = (xright-xleft)/50;
    xxleft = [xleft:h:xright];
    yyleft = spline(x(1:7),[.8, y(1:7), .2],xxleft); 

    xleft = x(8); xright = x(18);
    h = (xright-xleft)/50;
    xxright = [xleft:h:xright];
    yyright = spline(x(8:18),[-1.8,y(8:18),-.35],xxright); 

    figure(1); plot(x,y,'r.',xxleft,yyleft,'b-',xxright,yyright,'b-');
    title('Interpolation with Two Cubic Splines Having Specified End Derivatives'); 
    axis equal; 
  end

  if ind == 4
    disp(' '); 
    disp('Fit the data using MATLABs pchip.'); 
    disp(' '); 
    
    x = zeros(1,17);
    y = x;

    x(1:6) = 0:1:5;
    x(7:10) = 6.1:.2:6.7;
    x(11:13) = 7:1:9;
    x(14:17) = 12:2:18; 

    y(1) = 2.51; y(2) = 3.30; y(3) = 4.04; 
    y(4) = 4.70; y(5) = 5.22; y(6) = 5.54; 
    y(7) = 5.80; y(8) = 5.55; y(9) = 5.44; 
    y(10) = 5.40; y(11) = 5.40; y(12) = 5.57; 
    y(13) = 5.70; y(14) = 5.84; y(15) = 5.75; 
    y(16) = 5.48; y(17) = 4.90;

    xleft = x(1); xright = x(17);
    h = (xright-xleft)/100;
    xx = [xleft:h:xright];
    yy = pchip(x,y,xx); 

    figure(1); plot(x,y,'r.',xx,yy,'b-');
    title('PCHIP Interpolation');
    axis equal; 
  end

end
