%
% Exercise 3-5.32 
% Use Algorithm 3.5
%
n       = 8;
x       = zeros(n+1,1);
fun_val = zeros(n+1,1); 

x(1)    = 1;
x(2)    = 2;
x(3)    = 5;
x(4)    = 6;
x(5)    = 7;
x(6)    = 8;
x(7)    = 10;
x(8)    = 13;
x(9)    = 17;
% x(10)   = 20;
% x(11)   = 23;
% x(12)   = 24;
% x(13)   = 25;
% x(14)   = 27;
% x(15)   = 27.7;
% x(16)   = 28;
% x(17)   = 29;
% x(18)   = 30;

fun_val(1 ) = 3.0;
fun_val(2 ) = 3.7;
fun_val(3 ) = 3.9;
fun_val(4 ) = 4.2;
fun_val(5 ) = 5.7;
fun_val(6 ) = 6.6;
fun_val(7 ) = 7.1;
fun_val(8 ) = 6.7;
fun_val(9 ) = 4.5;
% fun_val(10) = 7.0;
% fun_val(11) = 6.1;
% fun_val(12) = 5.6;
% fun_val(13) = 5.8;
% fun_val(14) = 5.2;
% fun_val(15) = 4.1;
% fun_val(16) = 4.6;
% fun_val(17) = 4.1;
% fun_val(18) = 3.0;
df_val      = [ 1  -2/3 ];

[b, c, d] = cubic_spline_clamp(x, fun_val, df_val);

% for ii = 1:n
%     h(ii) = x(ii+1) - x(ii);
% end
% 
% A(1,1) = 2 * h(1);
% A(1,2) = h(1);
% rhs(1) = 3 * ( fun_val(2) - fun_val(1) ) / h(1) - 3 * df_val(1);
% for ii  = 1:n-1
%     A(ii+1,ii+1) = ( h(ii,1) + h(ii+1,1) ) * 2;
%     A(ii+1,ii  ) = h(ii,1);
%     A(ii+1,ii+2) = h(ii+1,1);
%     rhs(ii+1,1)  = 3 * ( fun_val(ii+2,1) - fun_val(ii+1,1) ) / h(ii+1,1);
%     rhs(ii+1,1)  = rhs(ii+1,1) - 3 * ( fun_val(ii+1,1) - fun_val(ii,1) ) / h(ii,1);
% end
% A(n+1,n  ) = h(n);
% A(n+1,n+1) = 2 * h(n);
% rhs(n+1)   = 3 * df_val(2) - 3 * ( fun_val(n+1) - fun_val(n) ) / h(n);
% 
% c = A \ rhs;
% 
% b = zeros(n,1);
% d = zeros(n,1);
% for ii = 1:n
%     b(ii,1) = ( fun_val(ii+1,1) - fun_val(ii,1) ) / h(ii,1);
%     b(ii,1) = b(ii,1) - h(ii,1) * ( c(ii+1,1) + 2 * c(ii,1) ) / 3;
%     d(ii,1) = ( c(ii+1,1) - c(ii,1) ) / ( 3 * h(ii,1) );
% end

%
% Plot figure of 
%      S(x) = fun_val(j)+b(j)(x-x(j))+c(j)(x-x(j))^2+d(j)(x-x(j))^3
% on [ x(j), x(j+1) ] for j = 1, 2, ..., n
%
m       = 10;
sub_x   = zeros(n*m,1);
y       = zeros(n*m,1); 
for ii = 1:n
    sub_x((ii-1)*m+1:ii*m,1) = linspace(x(ii,1), x(ii+1,1), m);
    y((ii-1)*m+1:ii*m,1) = fun_val(ii) * ones(m,1) + ...
                 b(ii) * ( sub_x((ii-1)*m+1:ii*m,1) - x(ii) * ones(m,1) );
    y((ii-1)*m+1:ii*m,1) = y((ii-1)*m+1:ii*m,1) + ...
        c(ii) * (sub_x((ii-1)*m+1:ii*m,1)-x(ii)*ones(m,1)).^2;
    y((ii-1)*m+1:ii*m,1) = y((ii-1)*m+1:ii*m,1) + ...
        d(ii) * (sub_x((ii-1)*m+1:ii*m,1)-x(ii)*ones(m,1)).^3; 
end
figure(1)
plot(sub_x,y,'LineWidth',2) 
hold on

n       = 6;
x       = zeros(n+1,1);
fun_val = zeros(n+1,1); 

x(1)    = 17;
x(2)   = 20;
x(3)   = 23;
x(4)   = 24;
x(5)   = 25;
x(6)   = 27;
x(7)   = 27.7;
% x(16)   = 28;
% x(17)   = 29;
% x(18)   = 30;

fun_val(1 ) = 4.5;
fun_val(2) = 7.0;
fun_val(3) = 6.1;
fun_val(4) = 5.6;
fun_val(5) = 5.8;
fun_val(6) = 5.2;
fun_val(7) = 4.1;
% fun_val(16) = 4.6;
% fun_val(17) = 4.1;
% fun_val(18) = 3.0;
df_val      = [ 3  -4 ];

[b, c, d] = cubic_spline_clamp(x, fun_val, df_val);

sub_x   = zeros(n*m,1);
y       = zeros(n*m,1); 
for ii = 1:n
    sub_x((ii-1)*m+1:ii*m,1) = linspace(x(ii,1), x(ii+1,1), m);
    y((ii-1)*m+1:ii*m,1) = fun_val(ii) * ones(m,1) + ...
                 b(ii) * ( sub_x((ii-1)*m+1:ii*m,1) - x(ii) * ones(m,1) );
    y((ii-1)*m+1:ii*m,1) = y((ii-1)*m+1:ii*m,1) + ...
        c(ii) * (sub_x((ii-1)*m+1:ii*m,1)-x(ii)*ones(m,1)).^2;
    y((ii-1)*m+1:ii*m,1) = y((ii-1)*m+1:ii*m,1) + ...
        d(ii) * (sub_x((ii-1)*m+1:ii*m,1)-x(ii)*ones(m,1)).^3; 
end
figure(1)
plot(sub_x,y,'LineWidth',2)  

n       = 3;
x       = zeros(n+1,1);
fun_val = zeros(n+1,1); 

x(1)   = 27.7;
x(2)   = 28;
x(3)   = 29;
x(4)   = 30;

fun_val(1) = 4.1;
fun_val(2) = 4.6;
fun_val(3) = 4.1;
fun_val(4) = 3.0;
df_val     = [ 1/3  -3/2 ];

[b, c, d] = cubic_spline_clamp(x, fun_val, df_val);

sub_x   = zeros(n*m,1);
y       = zeros(n*m,1); 
for ii = 1:n
    sub_x((ii-1)*m+1:ii*m,1) = linspace(x(ii,1), x(ii+1,1), m);
    y((ii-1)*m+1:ii*m,1) = fun_val(ii) * ones(m,1) + ...
                 b(ii) * ( sub_x((ii-1)*m+1:ii*m,1) - x(ii) * ones(m,1) );
    y((ii-1)*m+1:ii*m,1) = y((ii-1)*m+1:ii*m,1) + ...
        c(ii) * (sub_x((ii-1)*m+1:ii*m,1)-x(ii)*ones(m,1)).^2;
    y((ii-1)*m+1:ii*m,1) = y((ii-1)*m+1:ii*m,1) + ...
        d(ii) * (sub_x((ii-1)*m+1:ii*m,1)-x(ii)*ones(m,1)).^3; 
end
figure(1)
plot(sub_x,y,'LineWidth',2)  
axis([0 30 0 30])