function ds = demo_spline_atan()

    
%% Aproximovana funkce
%% a jeji prvni a druha derivace analyticky (jen pro vykreslovani)
f      = inline('atan(x/alpha)','x','alpha');
fder   = inline('alpha./(x.^2+alpha^2)','x','alpha');
fder2  = inline('-2*alpha.*x./(x.^2+alpha^2).^2','x','alpha');

%% Parametr funkce
alpha = 5;

%% Interval pro aproximaci
x_min = -5.0;
x_max =  5.0;

%% Okrajove podminky: zadavame hodnoty derivace nebo prirozeny spline
%
%% Pripad 1: spravne derivace
yder_left  = 0.1; 
yder_right = 0.1; 
natural_left  = false;
natural_right = false;
%
%% Pripad 2: nespravne derivace - generuje oscilace u okraje
%yder_left  = 0.0; 
%yder_right = 0.0; 
%natural_left  = false;
%natural_right = false;
%
%% Pripad 3: prirozeny spline
%yder_left  = 0.0; % hodnota se nepouziva, jen plnime kvuli volani
%yder_right = 0.0; % hodnota se nepouziva, jen plnime kvuli volani
%natural_left  = true;
%natural_right = true;

%% Pocty bodu
bodu_input  = 20;   % Pocet "zadanych" bodu (vygeneruje se z f(x))
bodu_approx = 100;  % V kolika bodech chceme aproximovat
bodu_exact  = 200;  % Pocet bodu pro vykresleni analyt. funkci f, f' a f''


%% 1)
%% Aproximace funkce z BODU_INPUT bodu na BODU_APPROX bodech pomoci splinu.

% vygenerujeme diskretni hodnoty
dx = (x_max-x_min) / (bodu_input-1);
x  = (x_min:dx:x_max);
y  = f(x,alpha);

% napocitame hodnoty druhych derivaci pro konstrukci splinu
yder2 = NUMREC_spline( x , y , natural_left , yder_left, natural_right, yder_right );

% zvolime body a napocitame v nich aproximovane hodnoty
dx_approx = (x_max-x_min) / (bodu_approx-1);
x_approx  = (x_min:dx_approx:x_max);
y_approx  = zeros(size(x_approx));
for i=1:bodu_approx
    y_approx(i) = NUMREC_splint(x, y, yder2, x_approx(i));
end;

% napocitame hustou sit s presnym (analytickym) resenim
dx_exact = (x_max-x_min) / (bodu_exact-1);
x_exact  = (x_min:dx_exact:x_max);
y_exact  = f(x_exact,alpha);


% stanovime rozmery obrazku dle velikosti obrazovky
scr = get( 0, 'Screensize' );
fig_wid = floor(scr(3)/3.5);
fig_off = floor(scr(3)/10);
fig_hei = floor(fig_wid/1.5);


% nakreslime
hf = figure;
set(hf,'Position',[fig_off fig_hei fig_wid fig_hei]);

ex = plot(x_exact,y_exact,'b-');
hold on;
in = plot(x,y,'kx','linewidth',3,'markersize',10);
ap = plot(x_approx,y_approx,'ro','linewidth',2);
hold off;
pl_max = max(y_exact(:));
pl_min = min(y_exact(:));
dpl = (pl_max-pl_min)/20;
axis([x_min x_max pl_min-dpl pl_max+dpl]);
title('f(x)','fontsize',16);
legend([ex, in, ap],{'exact','input','approx'},'Location','southeast','fontsize',12)


%% 2)
%% Aproximace prvni derivace na BODU_INPUT bodech pomoci derivace splinu.

% napocitame analyticke hodnoty prvni derivace
yder_exact  = fder(x_exact,alpha);

% v bodech X_APPROX napocitame derivace splinu
yder_approx = zeros(size(x_approx));
for i=1:bodu_approx
    yder_approx(i) = NUMREC_splder(x, y, yder2, x_approx(i));
end;

% nakreslime
hf = figure;
set(hf,'Position',[fig_off+fig_wid fig_hei fig_wid fig_hei]);

ex = plot(x_exact,yder_exact,'b-');
hold on;
ap = plot(x_approx,yder_approx,'ro','linewidth',2);
hold off;
pl_max = max(yder_exact(:));
pl_min = min(yder_exact(:));
dpl = (pl_max-pl_min)/20;
axis([x_min x_max pl_min-dpl pl_max+dpl]);
title('f''(x)','fontsize',16);
legend([ex, ap],{'exact','approx'},'Location','northeast','fontsize',12)


%% 3)
%% Vykresleni vypoctenych druhych derivaci pouzitych ke konstrukci splinu.

% napocitame analyticke hodnoty druhe derivace
yder2_exact = fder2(x_exact,alpha);

% nakreslime je spolu s hodnotami druhe derivace spoctenymi pri
% konstrukci splinu
hf = figure;
set(hf,'Position',[fig_off+2*fig_wid fig_hei fig_wid fig_hei]);

ex = plot(x_exact,yder2_exact,'b-');
hold on;
ap = plot(x,yder2,'r+','linewidth',2);
hold off;
pl_max = max(yder2_exact(:));
pl_min = min(yder2_exact(:));
dpl = (pl_max-pl_min)/20;
axis([x_min x_max pl_min-dpl pl_max+dpl]);
title('f"(x)','fontsize',16);
legend([ex, ap],{'exact','approx'},'fontsize',12)


end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Subrutiny pro konstrukci a vypocet splinu %
%         podle Numerical Recipes           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Subrutina pro vypocet druhych derivaci
function [yder2] = NUMREC_spline(x, y, l_natur, l_yder, r_natur, r_yder)
   % Calculate second derivatives for the interpolation by cubic spline
   %
   % Written according to routine SPLINE from Numerical Recipes
   %
   % INPUT:
   % - Arrays X and Y of length NPT containing a tabulated function, 
   %   i.e., y_i = f(xi), with x_1 < x_2 < .. . < x_NPT
   % - On each boundary either flag to use natural spline
   %   or value of first derivative      
   %    
   % OUTPUT: 
   % - array YDER2 of length npt which contains the second derivatives
   %   of the interpolating function at the tabulated points X_i. 
       
   npt   = length(x);
   u     = zeros(size(x));
   yder2 = zeros(size(x));
   %
   if (l_natur)  % Left BC natural
       yder2(1) = 0.0;
       u(1)     = 0.0;
   else          % Left BC specified by first derivative l_yder
       yder2(1)  = -0.5;
       u(1)      = (3.0/(x(2)-x(1))) * ((y(2)-y(1)) / (x(2)-x(1)) - l_yder);
   end;
   % Decomposition loop of the tridiagonal algorithm. 
   % yder2 and u are used for temporary storage of the decomposed factors.
   for i = 2:npt-1
       sig      = (x(i)-x(i-1)) / (x(i+1)-x(i-1));
       p        = sig * yder2(i-1) + 2.0;
       yder2(i) = (sig-1.0) / p;
       u(i)     = ( 6. * (  (y(i+1)-y(i)) / (x(i+1)-x(i)) ...
                           -(y(i)-y(i-1)) / (x(i)-x(i-1)) ) ...
                       / (x(i+1)-x(i-1)) - sig*u(i-1) ...
                  ) / p;
   end;
   %
   if (r_natur)  % Right BC natural
       qn = 0.0;
       un = 0.0;
   else          % Right BC specified by first derivative r_yder
       qn = 0.5;
       un = (3.0/(x(npt)-x(npt-1))) ...
            * (r_yder - (y(npt)-y(npt-1))/(x(npt)-x(npt-1)));
   end;
   % Backsubstitution loop of the tridiagonal algorithm.
   yder2(npt) = (un-qn*u(npt-1)) / (qn*yder2(npt-1) + 1.0);
   for k=npt-1:-1:1
      yder2(k) = yder2(k) * yder2(k+1) + u(k);
   end;
end


%%% Subrutina pro vypocet interpolovane hodnoty
function y = NUMREC_splint(xa, ya, yader2, x)
   % Calculate the interpolated value at point x by cubic spline
   %
   % Written according to routine SPLINT from Numerical Recipes
   %
   % INPUT:
   % - Arrays XA and YA of length NPT containing a tabulated function, 
   %   i.e., y_i = f(xi), with x_1 < x_2 < .. . < x_NPT
   % - Arrays XADER2 containing the second derivatives
   % - Point X
   %    
   % OUTPUT: 
   % - Value of Y(X) approximated by cubic spline

   npt   = length(xa);
   %    
   % We will find the right place in the table by means of bisection.
   % This is optimal if sequential calls to this routine are at random
   % values of x. If sequential calls are in order, and closely spaced,
   % one would do better to store previous values of klo and khi
   % and test if they remain appropriate on the next call.
   klo = 1;
   khi = npt;
   while (khi-klo > 1)
       k = floor ((khi+klo)/2);
       if (xa(k) > x) 
           khi = k; 
       else
           klo = k;
       end;
   end;
   % klo and khi now bracket the input value of x.
   h = xa(khi)-xa(klo);
   if (h==0.0)
       disp(['Error in routine NUMREC_splint. Points of XA must be distinct.']);
       return;
   end; 
   % Cubic spline polynomial is now evaluated.
   a = (xa(khi)-x) / h; 
   b = (x-xa(klo)) / h;
   y = a*ya(klo) + b*ya(khi) ...
       + ( (a^3-a)*yader2(klo) + (b^3-b)*yader2(khi) ) * h^2 / 6.0;
end


%%% Subrutina pro vypocet hodnoty derivace splinu
function yder = NUMREC_splder(xa, ya, yader2, x)
   % Calculate the derivative of cubic spline at point x
   %
   % Written according to routine SPLDER from Numerical Recipes
   %
   % INPUT:
   % - Arrays XA and YA of length NPT containing a tabulated function, 
   %   i.e., y_i = f(xi), with x_1 < x_2 < .. . < x_NPT
   % - Arrays XADER2 containing the second derivatives
   % - Point X
   %    
   % OUTPUT: 
   % - Value of Y'(X) approximated by cubic spline

   npt   = length(xa);
   %    
   % We will find the right place in the table by means of bisection.
   % This is optimal if sequential calls to this routine are at random
   % values of x. If sequential calls are in order, and closely spaced,
   % one would do better to store previous values of klo and khi
   % and test if they remain appropriate on the next call.
   klo = 1;
   khi = npt;
   while (khi-klo > 1)
       k = floor ((khi+klo)/2);
       if (xa(k) > x) 
           khi = k; 
       else
           klo = k;
       end;
   end;
   % klo and khi now bracket the input value of x.
   h = xa(khi)-xa(klo);
   if (h==0.0)
       disp(['Error in routine NUMREC_splder. Points of XA must be distinct.']);
       return;
   end; 
   % Cubic spline polynomial is now evaluated.
   a = (xa(khi)-x) / h; 
   b = (x-xa(klo)) / h;
   yder = ( ya(khi)-ya(klo) ) / h ...
          + ( -(3*a^2-1)*yader2(klo) + (3*b^2-1)*yader2(khi) ) * h / 6.0;
end




