function demo_hermite_interpolation
% plots the Hermite interpolation of a function
% and computes the error

% Initialize data
global AMB_Data;
AMB_Data=[];
AMB_Data.npts = 301;  
AMB_Data.Fun = @(x) x.*sin(x);  
AMB_Data.Funstring = '@(x) x.*sin(x)';  
AMB_Data.DFun = @(x) sin(x)+x.*cos(x);  
AMB_Data.DFunstring = '@(x) sin(x)+x.*cos(x)';  
AMB_Data.interval = [0 10];
AMB_Data.x=[0 10];

% Initialize GUI
clf reset;
set(gcf,'WindowButtonDownFcn',@myDownfct);

% Text for title
uicontrol('Units','normalized', ...
  	'BackgroundColor',get(gcf,'Color'), ...
	'HorizontalAlignment','center', ...
        'FontWeight','bold',...
        'FontSize',20,...
	'Position',[.05,.92,.9,.07],...
	'String','Hermite Interpolation',...
	'Style','text');

% Text for function
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'FontWeight','bold',...
        'HorizontalAlignment','left', ...
	'Position',[.05,.85,.3,.05],...
	'String','Function:',...
	'Style','text');


% Input for function
AMB_Data.uicontrols(1)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'FontWeight','bold',...
        'HorizontalAlignment','left', ...
	'Position',[.2,.85,.4,.05],...
	'String',sprintf('%s',AMB_Data.Funstring),...
	'Style','edit',...
	'Callback',@function_update);

% Text for derivative
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'FontWeight','bold',...
        'HorizontalAlignment','left', ...
	'Position',[.05,.8,.3,.05],...
	'String','Derivative:',...
	'Style','text');


% Input for derivative
AMB_Data.uicontrols(2)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'FontWeight','bold',...
        'HorizontalAlignment','left', ...
	'Position',[.2,.8,.4,.05],...
	'String',sprintf('%s',AMB_Data.DFunstring),...
	'Style','edit',...
	'Callback',@function_update);

% Text for lower interval bound
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.65,.85,.05,.05],...
	'String','a:',...
	'Style','text');

% Input for lower interval bound
AMB_Data.uicontrols(3)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.7,.85,.1,.05],...
	'String',sprintf('%5.2f',AMB_Data.interval(1)),...
	'Style','edit',...
	'Callback',@interval_update);

% Text for upper interval bound
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.65,.8,.05,.05],...
	'String','b:',...
	'Style','text');

% Input for upper interval bound
AMB_Data.uicontrols(4)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.7,.8,.1,.05],...
	'String',sprintf('%5.2f',AMB_Data.interval(2)),...
	'Style','edit',...
	'Callback',@interval_update);

% Restart button
uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'String','Restart', ...
	'Position',[.85,.85,.1,.05], ...
	'Callback',@restart,...
	'Style','pushbutton');

% Text for error
AMB_Data.uicontrols(5)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.05,0,.3,.05],...
	'String',sprintf('Error: %5.3f',100),...
	'Style','text');

% Text for click
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'ForegroundColor',[0,0,0], ...
        'FontWeight','bold',...
        'HorizontalAlignment','right', ...
	'Position',[.5,0,.45,.05],...
	'String','Mouse: left = insert point ',...
	'Style','text');


% Axis for function display
AMB_Data.axis(1)=subplot('position',[.05,.1,.9,.65]);

% Display first approximation
myplot;


%%%  Callback Functions %%%

% new function input
function function_update(varargin)
global AMB_Data;
AMB_Data.Funstring = get(AMB_Data.uicontrols(1),'String');
AMB_Data.Fun= eval(AMB_Data.Funstring);
AMB_Data.DFunstring = get(AMB_Data.uicontrols(2),'String');
AMB_Data.DFun= eval(AMB_Data.DFunstring);
restart;

% new errorbound input
function interval_update(varargin)
global AMB_Data;
AMB_Data.interval(1) = str2double(get(AMB_Data.uicontrols(3),'String'));
AMB_Data.interval(2) = str2double(get(AMB_Data.uicontrols(4),'String'));
restart;

% button restart pressed
function restart(varargin)
global AMB_Data;
AMB_Data.x=AMB_Data.interval;
myplot;

% mouse button pressed 
function myDownfct(varargin)
global AMB_Data;
if strcmp(get(gcf,'SelectionType'),'normal') % left mouse: new point    
  p=get(AMB_Data.axis(1),'CurrentPoint');
  if p(1) >= AMB_Data.interval(1) & p(1)  <= AMB_Data.interval(2)
    np=round(100*p(1))/100;
    AMB_Data.x= unique([AMB_Data.x,np]);
    myplot;
  end    
end

%%% Internal Functions %%%

function myplot;
global AMB_Data;
subplot(AMB_Data.axis(1))
cla;hold on

% compute Hermite interpolation
fx=AMB_Data.Fun(AMB_Data.x);
dfx=AMB_Data.DFun(AMB_Data.x);
[t,p,data]= hermite(AMB_Data.x,fx,dfx,AMB_Data.npts);

% plot function and visualize prescribed tolerance
plot(t,AMB_Data.Fun(t),'-','linewidth',2,'color',[0,0,.6]) 
ymin=min(AMB_Data.Fun(t));
ymax=max(AMB_Data.Fun(t));
ydelta=(ymax-ymin)*1.5;
axis([t(1),t(end),(ymin+ymax-ydelta)/2,(ymin+ymax+ydelta)/2]);
axis manual

% plot Hermite approximation 
h = diff(AMB_Data.x); 
xn = ones(4,1)*AMB_Data.x(1:end-1) + ([0:3]'/3)*h; 
plot(t,p,'k-','linewidth',2);
plot(xn,data,'--k');
set(gca,'xtick',AMB_Data.x);
a=axis;
plot(ones(2,1)*AMB_Data.x,[fx;a(3)*ones(size(AMB_Data.x))],':k');
plot(AMB_Data.x,AMB_Data.Fun(AMB_Data.x),'ok');

% display error and check for bound
e=norm(AMB_Data.Fun(t)-p,inf);
set(AMB_Data.uicontrols(5),'String',sprintf('Error: %5.3f',e));

%%% Functions from AMB_Programs (included for standalone use) %%% 

function [t,p,c] = hermite(x,f,d,m) 
% HERMITE: evaluation of a Hermite spline
% p = hermite_spline(x,f,d,m) 
% m+1 values p(j) per interval at equally spaced points t(j) and Bernstein coefficients c(1:4,:)
% of a Hermite spline with data f(k), d(k) at x(k)

% convert to column vectors (if necessary)
x = x(:); f = f(:); d = d(:);

% evaluation points
h = diff(x); 
t = ones(m+1,1)*x(1:end-1)' + ([0:m]'/m)*h'; 
t = t(:);

% Bernstein coefficients
c = [f(1:end-1), f(1:end-1)+d(1:end-1).*h/3, ...
   f(2:end)-d(2:end).*h/3, f(2:end)]';

% evaluation 
b = bernstein(3,[0:m]'/m);
p = b*c; 
p = p(:);
 

function b = bernstein(n,t)
% BERNSTEIN: Bernstein polynomials
% b = bernstein(n,t)
% values b(j,k) of b^n_{k-1} at t(j)

% change to column vector (if necessary)
t=t(:);

% initialization for degree 0 
b = zeros(length(t),2*n+1); 
b(:,n+1) = 1;

% recursion for Bernstein polynomials
for k=1:n
   tk = repmat(t,1,2*n+1-k);
   b = (1-tk).*b(:,2:end) + tk.*b(:,1:end-1);
end

