function varargout = Spline_GUI(varargin)
% SPLINE_GUI Application M-file for Spline_GUI.fig
%    FIG = SPLINE_GUI launch Spline_GUI GUI.
%    SPLINE_GUI('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.5 21-Sep-2004 14:42:24

if nargin == 0  % LAUNCH GUI

	fig = openfig(mfilename,'reuse');

	% Generate a structure of handles to pass to callbacks, and store it. 
	handles = guihandles(fig);
    handles.indexfcn=1;
    handles.indexbc=1;
    handles.subdivisions=4;
    handles.const_a=0;
    handles.const_b=1;
	handles.problem_data = [handles.indexfcn, handles.subdivisions,...
        handles.const_a, handles.const_b,handles.indexbc];
    guidata(fig, handles)

	if nargout > 0
		varargout{1} = fig;
	end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch
		disp(lasterr);
	end

end


%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and 
%| sets objects' callback properties to call them through the FEVAL 
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the 
%| callback type separated by '_', e.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(hObject, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.



% --------------------------------------------------------------------
function varargout = popup_fcn_menu_Callback(hObject, eventdata, handles, varargin)

indexfcn = get(hObject,'Value');
handles.indexfcn = indexfcn;
	handles.problem_data = [handles.indexfcn, handles.subdivisions,...
        handles.const_a, handles.const_b,handles.indexbc];
guidata(hObject,handles);     % Save the handles structure after adding data

% --------------------------------------------------------------------
function varargout = SelectSubdiv_Callback(hObject, eventdata, handles, varargin)

subdiv = eval(get(handles.SelectSubdiv,'string'));
handles.subdivisions = subdiv;
	handles.problem_data = [handles.indexfcn, handles.subdivisions,...
        handles.const_a, handles.const_b,handles.indexbc];
guidata(hObject,handles);     % Save the handles structure after adding data

% --------------------------------------------------------------------
function varargout = LeftHand_Callback(hObject, eventdata, handles, varargin)

const_a = eval(get(handles.LeftHand,'string'));
handles.const_a = const_a;
	handles.problem_data = [handles.indexfcn, handles.subdivisions,...
        handles.const_a, handles.const_b,handles.indexbc];
guidata(hObject,handles);     % Save the handles structure after adding data

% --------------------------------------------------------------------
function varargout = RightHand_Callback(hObject, eventdata, handles, varargin)

const_b = eval(get(handles.RightHand,'string'));
handles.const_b = const_b;
	handles.problem_data = [handles.indexfcn, handles.subdivisions,...
        handles.const_a, handles.const_b,handles.indexbc];
guidata(hObject,handles);     % Save the handles structure after adding data

% --------------------------------------------------------------------
function varargout = ErrorOutput_Callback(hObject, eventdata, handles, varargin)


% --------------------------------------------------------------------
function varargout = CalculatePushbutton_Callback(hObject, eventdata, handles, varargin)

prob_data = handles.problem_data;
indexfcn = prob_data(1);
subdiv=prob_data(2);
a = prob_data(3);
b = prob_data(4);
indexbc = prob_data(5);

stepsize = (b-a)/subdiv;
nodes = a:stepsize:b;
x_eval = a:(b-a)/(20*subdiv):b;

% Calculate interpolating spline function.  Evaluate it
% and the true function at the points given in x_eval.
[y_eval,true] = fcn(x_eval,prob_data);


max_y = max([max(y_eval),max(true)]); min_y = min([min(y_eval),min(true)]);
delta_y = (max_y - min_y)/10;
delta_x = (b-a)/10;

% Graph f(x) vs. P(x)
axes(handles.axes1)
cla
axis([a-delta_x,b+delta_x,min_y-delta_y,max_y+delta_y])
hold on
plot(x_eval,true,x_eval,y_eval,'r:')
h_step=(b-a)/subdiv;
x=a:h_step:b;
y=f(x,indexfcn);
plot(x,y,'.','MarkerSize',10)
legend('Original fcn','Interpolant',0)
hold off

error = true-y_eval;
max_y = max(error); min_y = min(error);
delta_y = (max_y - min_y)/10;

% Graph error in P(x)
axes(handles.axes2)
cla
axis([a-delta_x,b+delta_x,min_y-delta_y,max_y+delta_y])
hold on
plot(x_eval,error)
plot([a-delta_x b+delta_x/2],[0 0],'k')
text(b+.6*delta_x,0,'\itx')
hold off
max_error = max(abs(true-y_eval));
StringData=sprintf('Maximum Absolute Error = %12.3e',max_error);
set(handles.ErrorOutput,'string',StringData)
guidata(hObject,handles);



% --------------------------------------------------------------------
function varargout = Menu_file_Callback(hObject, eventdata, handles, varargin)

if isempty(get(handles.axes1,'Children'))
    set(handles.Submenu_save_axes1,'Enable','off')
    set(handles.Submenu_save_axes2,'Enable','off')
else
    set(handles.Submenu_save_axes1,'Enable','on')
    set(handles.Submenu_save_axes2,'Enable','on')
end

% --------------------------------------------------------------------
function varargout = Submenu_save_axes1_Callback(hObject, eventdata, handles, varargin)

h_invisible = figure('visible','off'); 
hax = copyobj(handles.axes1,h_invisible);
set(hax,'units',get(0,'defaultaxesunits'),'position',get(0,'defaultaxesposition'))
print(h_invisible,'-depsc','spline') 
close(h_invisible);
disp('The graph is stored under the name "spline.eps"')

% --------------------------------------------------------------------
function varargout = Submenu_save_axes2_Callback(hObject, eventdata, handles, varargin)

h_invisible = figure('visible','off'); 
hax = copyobj(handles.axes2,h_invisible);
set(hax,'units',get(0,'defaultaxesunits'),'position',get(0,'defaultaxesposition'))
print(h_invisible,'-depsc','spline_error') 
close(h_invisible);
disp('The graph is stored under the name "spline_error.eps"')

% --------------------------------------------------------------------
function varargout = Submenu_close_Callback(hObject, eventdata, handles, varargin)

delete(handles.figure1)

% --------------------------------------------------------------------

% --- Executes on selection change in BoundaryCondition.
function BoundaryCondition_Callback(hObject, eventdata, handles)
% hObject    handle to BoundaryCondition (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns BoundaryCondition contents as cell array
%        contents{get(hObject,'Value')} returns selected item from BoundaryCondition

indexbc = get(hObject,'Value');
handles.indexbc = indexbc;
	handles.problem_data = [handles.indexfcn, handles.subdivisions,...
        handles.const_a, handles.const_b,handles.indexbc];
guidata(hObject,handles);     % Save the handles structure after adding data

% --------------------------------------------------------------------

% --- Executes during object creation, after setting all properties.
function BoundaryCondition_CreateFcn(hObject, eventdata, handles)
% hObject    handle to BoundaryCondition (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --------------------------------------------------------------------

function [y_eval,true] = fcn(x_eval,prob_data)

indexfcn = prob_data(1);
subdiv=prob_data(2);
a = prob_data(3);
b = prob_data(4);
indexbc = prob_data(5);

true = f(x_eval,indexfcn);

% Construct data for interpolating spline.
stepsize=(b-a)/subdiv;
x=a:stepsize:b;
y=f(x,indexfcn);

switch indexbc
    case 1 % natural boundary conditions
        y_eval = spline_natural(x,y,x_eval,a,b,stepsize,subdiv);
    case 2 % clamped boundary conditions
        y_eval = spline_clamped(x,y,x_eval,a,b,stepsize,subdiv,indexfcn);
    case 3 % not-a-knot boundary conditions
        y_eval = spline(x,y,x_eval);
    case 4 % not-a-knot boundary conditions, with first and last 
           % subintervals halved.
           x = [a,a+stepsize/2,x(2:subdiv),b-stepsize/2,b];
           y = [y(1),f(a+stepsize/2,indexfcn),y(2:subdiv),f(b-stepsize/2,indexfcn),y(subdiv+1)];
           y_eval = spline(x,y,x_eval);
end

% --------------------------------------------------------------------

function y_eval = spline_natural(x,y,x_eval,a,b,stepsize,subdiv)

% x and y have length n = subdiv+1.
n = subdiv+1;
% Set up the tridiagonal linear system for obtaining the second derivative
% values associated with the natural cubic interpolating spline for the 
% data given in x and y. Then solve using tridiag.
avec = zeros(n,1);
avec(2:n-1) = (stepsize/6)*ones(n-2,1);
cvec = avec;
bvec = ones(n,1);
bvec(2:n-1) = (2*stepsize/3)*bvec(2:n-1);
rhs = zeros(n,1);
rhs(2:n-1) = (y(3:n)-2*y(2:n-1)+y(1:n-2))/stepsize;
[mvec, alpha, beta, ier] = tridiag(avec,bvec,cvec,rhs,n,0);

y_eval = spline_evaluate(x,y,a,b,subdiv,stepsize,x_eval,mvec);

% --------------------------------------------------------------------

function y_eval = spline_clamped(x,y,x_eval,a,b,stepsize,subdiv,indexfcn)

% x and y have length n = subdiv+1.
n = subdiv+1;
% Obtain derivative values at endpoints.
dy = df([a,b],indexfcn);
% Set up the tridiagonal linear system for obtaining the second derivative
% values associated with the natural cubic interpolating spline for the 
% data given in x and y. Then solve using tridiag.
avec = zeros(n,1);
avec(2:n) = (stepsize/6)*ones(n-1,1);
cvec = zeros(n,1);
cvec(1:n-1) = (stepsize/6)*ones(n-1,1);
bvec = ones(n,1);
bvec(2:n-1) = (2*stepsize/3)*bvec(2:n-1);
bvec(1) = stepsize/3; bvec(n) = stepsize/3;
rhs = zeros(n,1);
rhs(2:n-1) = (y(3:n)-2*y(2:n-1)+y(1:n-2))/stepsize;
rhs(1) = (y(2)-y(1))/stepsize - dy(1);
rhs(n) = dy(2) - (y(n)-y(n-1))/stepsize;
[mvec, alpha, beta, ier] = tridiag(avec,bvec,cvec,rhs,n,0);

y_eval = spline_evaluate(x,y,a,b,subdiv,stepsize,x_eval,mvec);

% --------------------------------------------------------------------

function y_eval = spline_evaluate(x,y,a,b,subdiv,stepsize,x_eval,mvec);
nx = length(x_eval);
for i=1:nx
    t = x_eval(i);
    j = 1 + floor((t-a)/stepsize + 2*eps);
    if j == subdiv+1
        j = subdiv;
    end
    term1 = (mvec(j)*(x(j+1)-t).^3+mvec(j+1)*(t-x(j)).^3)/(6*stepsize);
    term2 = (y(j)*(x(j+1)-t)+y(j+1)*(t-x(j)))/stepsize;
    term3 = (stepsize/6)*(mvec(j)*(x(j+1)-t)+mvec(j+1)*(t-x(j)));
    y_eval(i) = term1 + term2 - term3;
end

% --------------------------------------------------------------------

function f_value=f(x,indexfcn)

switch indexfcn
case 1  % f(x) = exp(x)
    f_value = exp(x);
case 2  % f(x) = cos x
    f_value = cos(x);
case 3  % f(x) = arctan x
    f_value = atan(x);
case 4  % f(x) = 1/(1+x^2)
    f_value = 1./(1+x.^2);
case 5  % f(x) = log x
    f_value = log(x);
case 6  % f(x) = tan x
    f_value = tan(x);
case 7  % f(x) = exp(-x^2)
    f_value = exp(-x.^2);
case 8  % f(x) = 1/(2+cos(x))
    f_value = 1 ./(2+cos(x));
case 9  % f(x) = sqrt(abs(x))
    f_value = sqrt(abs(x));
end

% --------------------------------------------------------------------

function df_value=df(x,indexfcn)

switch indexfcn
case 1  % f(x) = exp(x)
    df_value = exp(x);
case 2  % f(x) = cos x
    df_value = -sin(x);
case 3  % f(x) = arctan x
    df_value = 1./(1+x.^2);
case 4  % f(x) = 1/(1+x^2)
    df_value = -2*x./(1+x.^2).^2;
case 5  % f(x) = log x
    df_value = 1./x;
case 6  % f(x) = tan x
    sec_x = 1./cos(x);
    df_value = sec_x.^2;
case 7  % f(x) = exp(-x^2)
    df_value = -2*x.*exp(-x.^2);
case 8  % f(x) = 1/(2+cos(x))
    df_value = sin(x)./(2+cos(x)).^2;
case 9  % f(x) = sqrt(abs(x))
    df_value = sign(x)./sqrt(abs(x));
end

% --------------------------------------------------------------------

function [x, alpha, beta, ier] = tridiag(a,b,c,f,n,iflag)

% function [x, alpha, beta, ier] = tridiag(a,b,c,f,n,iflag)
%
% Solve a tridiagonal linear system M*x=f
%
% INPUT:
% The order of the linear system is given in n.
% The subdiagonal, diagonal, and superdiagonal of M are given
% by the arrays a,b,c, respectively.  More precisely,
%     M(i,i-1) = a(i), i=2,...,n
%     M(i,i)   = b(i), i=1,...,n
%     M(i,i+1) = c(i), i=1,...,n-1
% iflag=0 means that the original matrix M is given as specified above.
% iflag=1 means that the LU factorization of M is already known and is
% stored in a,b,c.  This will have been accomplished by a previous call
% to this routine.  In that case, the vectors alpha and beta should 
% have been substituted for a and b in the calling sequence.
%
% OUTPUT:
% Upon exit, the LU factorization of M is already known and is stored
% in alpha,beta,c.  The solution x is given as well.
% ier=0 means the program was completed satisfactorily.
% ier=1 means that a zero pivot element was encountered and the 
% solution process was abandoned.

a(1) = 0;
if iflag == 0
   % Compute LU factorization of matrix M.
   for j=2:n
      if b(j-1) == 0
         ier = 1;
         return
      end
      a(j) = a(j)/b(j-1);
      b(j) = b(j) - a(j)*c(j-1);
   end
   if b(n) == 0 
      ier = 1;
      return
   end
end

% Compute solution x to M*x = f.
% Do forward substitution to solve lower triangular system.
for j=2:n
   f(j) = f(j) - a(j)*f(j-1);
end

% Do backward substitution to solve upper triangular system.
f(n) = f(n)/b(n);
for j=n-1:-1:1
   f(j) = (f(j) - c(j)*f(j+1))/b(j);
end

% Set output variables.
ier = 0;
x = f;
alpha = a; beta = b;

% --------------------------------------------------------------------

function Help_GUI_Callback(hObject, eventdata, handles)
% hObject    handle to Help_GUI (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

Spline_GUI_Help;