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

% Last Modified by GUIDE v2.5 03-Oct-2004 22:18:00

if nargin == 0  % LAUNCH GUI

	fig = openfig(mfilename,'reuse');

	% Generate a structure of handles to pass to callbacks, and store it. 
	handles = guihandles(fig);
    % Initialization of fcn_data = [indexfcn,a,b,c1,c2]
    handles.indexfcn = 1;
    handles.const_a = 0;
    handles.const_b = pi;
    handles.const_c1 = 1;
    handles.const_c2 = 1;
    handles.fcn_data = [handles.indexfcn,handles.const_a,handles.const_b,handles.const_c1,handles.const_c2]; 
    % Initialization of quad_data = [n_init,doubles]
    handles.const_n_init = 2;
    handles.const_double = 4;
    handles.increment = 1;
    handles.quad_data = [handles.const_n_init,handles.const_double,handles.increment];
    handles;
	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(h, 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.

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

% --- Executes on selection change in ChooseIntegrand.
function ChooseIntegrand_Callback(h, eventdata, handles)
% h          handle to ChooseIntegrand (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(h,'String') returns ChooseIntegrand contents as cell array
%        contents{get(h,'Value')} returns selected item from ChooseIntegrand

% disp('ChooseIntegrand_Callback - Begin')
indexfcn = get(h,'Value');
a = handles.const_a;
b = handles.const_b;
c1 = handles.const_c1;
c2 = handles.const_c2;
handles.fcn_data = [indexfcn,a,b,c1,c2];
guidata(h,handles); % Save the handles structure after adding data
% disp('ChooseIntegrand_Callback - End')

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

% --- Executes during object creation, after setting all properties.
function ChooseIntegrand_CreateFcn(h, eventdata, handles)
% h          handle to ChooseIntegrand (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(h,'BackgroundColor','white');
else
    set(h,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

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

function varargout = LeftEndEditText_Callback(h, eventdata, handles, varargin)

% disp('LeftEndEditText_Callback - Begin')

const_a = eval(get(handles.LeftEndEditText,'string'));
handles.const_a = const_a;
handles.fcn_data = [handles.indexfcn,handles.const_a, handles.const_b, handles.const_c1, handles.const_c2];
guidata(h,handles);
Gaussian_Quad('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

% disp('LeftEndEditText_Callback - End')

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

function varargout = RightEndEditText_Callback(h, eventdata, handles, varargin)

% disp('RightEndEditText_Callback - Begin')

const_b = eval(get(handles.RightEndEditText,'string'));
handles.const_b = const_b;
handles.fcn_data = [handles.indexfcn,handles.const_a, handles.const_b, handles.const_c1, handles.const_c2];
guidata(h,handles);
Gaussian_Quad('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

% disp('RightEndEditText_Callback - End')

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

function varargout = Constant_c1_EditText_Callback(h, eventdata, handles, varargin)

const_c1 = eval(get(handles.Constant_c1_EditText,'string'));
handles.const_c1 = const_c1;
handles.fcn_data = [handles.indexfcn,handles.const_a, handles.const_b, handles.const_c1, handles.const_c2];
guidata(h,handles);
Gaussian_Quad('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

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

function varargout = Constant_c2_EditText_Callback(h, eventdata, handles, varargin)

% disp('Constant_c2_EditText_Callback - Begin')

const_c2 = eval(get(handles.Constant_c2_EditText,'string'));
handles.const_c2 = const_c2;
handles.fcn_data = [handles.indexfcn,handles.const_a, handles.const_b, handles.const_c1, handles.const_c2];
guidata(h,handles);
Gaussian_Quad('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

% disp('Constant_c2_EditText_Callback - End')
% --------------------------------------------------------------------

function varargout = PlotPushbutton_Callback(h, eventdata, handles, varargin)

% disp('PlotPushbutton_Callback - Begin')

f_data = handles.fcn_data;
indexfcn = f_data(1);
a = f_data(2);
b = f_data(3);
c1 = f_data(4);
c2 = f_data(5);
x = a:(b-a)/100:b ;
y = fcn(x,f_data);
axes(handles.axes1)
plot(x,y)
grid on

% disp('PlotPushbutton_Callback - End')

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

function varargout = Choose_Ninit_Callback(h, eventdata, handles, varargin)

% disp('Choose_Ninit_Callback - Begin')

const_n_init = eval(get(handles.Choose_Ninit,'string'));
handles.const_n_init = const_n_init;
handles.quad_data = [handles.const_n_init,handles.const_double,handles.increment];
guidata(h,handles);
Gaussian_Quad('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

% disp('Choose_Ninit_Callback - End')

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

function varargout = ChooseDoubles_Callback(h, eventdata, handles, varargin)

% disp('ChooseDoubles_Callback - Begin')

const_double = eval(get(handles.ChooseDoubles,'string'));
handles.const_double = const_double;
handles.quad_data = [handles.const_n_init,handles.const_double,handles.increment];
guidata(h,handles);
Gaussian_Quad('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

% disp('ChooseDoubles_Callback - End')

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

function varargout = IntegrationOutput_Callback(h, eventdata, handles, varargin)

% --------------------------------------------------------------------
function varargout = file_menu_Callback(h, eventdata, handles, varargin)
if isempty(get(handles.axes1,'Children'))
    set(handles.print_submenu,'Enable','off')
else
    set(handles.print_submenu,'Enable','on')
end

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

function varargout = print_submenu_Callback(h, 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','gaussian_quad_graph') 
close(h_invisible);
disp('The graph is stored under the name "gaussian_quad_graph.eps"')

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

function varargout = table_submenu_Callback(h, eventdata, handles, varargin)

% Load data from handles structure
StringData = handles.StringData;

fid = fopen('output_Gauss_Quad','w');
disp('The table is saved under the name "output_Gauss_Quad"')

[num_rows,num_cols] = size(StringData);

for i=1:num_rows
    fprintf(fid,'%c',StringData(i,1:num_cols));
    fprintf(fid,'\n');
end

fclose(fid);

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

function varargout = close_submenu_Callback(h, eventdata, handles, varargin)
delete(handles.figure1)

% --------------------------------------------------------------------
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)

Gaussian_Quad_Help

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

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

% Hint: edit 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 arithmetic_increment_Callback(h, eventdata, handles)
% h          handle to arithmetic_increment (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(h,'String') returns contents of arithmetic_increment as text
%        str2double(get(h,'String')) returns contents of arithmetic_increment as a double

% disp('arithmetic_increment_Callback - Begin')

increment = eval(get(handles.arithmetic_increment,'string'));
handles.increment = increment;
handles.quad_data = [handles.const_n_init,handles.const_double,handles.increment];
guidata(h,handles);
Gaussian_Quad('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

% disp('arithmetic_increment_Callback - End')

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

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

% Hint: edit 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 varargout = GQ_Geometric_Callback(h, eventdata, handles, varargin)

% disp('GQ_Geometric_Callback - Begin')

handles.method = 'geometric';

quad_data = handles.quad_data;
n_init = quad_data(1);
doubles = quad_data(2);
increment = quad_data(3);

n_vec = n_init*2.^[0:doubles];
handles.n_vec = n_vec;    
guidata(h,handles);

Calculate_Integral(h,handles);
Gaussian_Quad('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

% disp('GQ_Geometric_Callback - End')

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

function varargout = GQ_Arithmetic_Callback(h, eventdata, handles, varargin)

% disp('GQ_Arithmetic_Callback - Begin')

handles.method = 'arithmetic';

quad_data = handles.quad_data;
n_init = quad_data(1);
doubles = quad_data(2);
increment = quad_data(3);

n_vec = n_init + increment*[0:doubles];
handles.n_vec = n_vec;    
guidata(h,handles);

Calculate_Integral(h,handles);
Gaussian_Quad('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

% disp('GQ_Arithmetic_Callback - End')

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

function Calculate_Integral(h,handles)

% disp('Calculate_Integral - Begin')

f_data = handles.fcn_data;

quad_data = handles.quad_data;
doubles = quad_data(2);

n_vec = handles.n_vec;

celldata = cellstr(num2str(zeros(doubles+3,1)));

[integral,difference,ratio_diff]=gauss_quad(f_data,quad_data,n_vec);
celldata(1) = cellstr('                          Gaussian quadrature output');

% Calculate true errors and their ratios.
true_integral = true_value(f_data);
true_error = true_integral - integral;
error_ratio = zeros(size(integral));
if doubles > 0
    for i=2:doubles+1
        if true_error(i) ~= 0
            error_ratio(i) = true_error(i-1)./true_error(i);
        else
            error_ratio(i) = inf;
        end
    end
end

celldata(2) = cellstr('    n          I_n             difference     ratio      true error     ratio');
for j=0:doubles
    n = n_vec(j+1);
    celldata(j+3) = cellstr(sprintf('%5i%23.14e%13.3e%13.3e%13.3e%13.3e',...
       n,integral(j+1),difference(j+1),ratio_diff(j+1),true_error(j+1),error_ratio(j+1)));
end
handles.StringData = char(celldata);
guidata(h,handles);

set(handles.IntegrationOutput,'string',handles.StringData);

% disp('Calculate_Integral - End')

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

function [integral,difference,ratio]=gauss_quad(f_data,quad_data,n_vec)
%
%  function [integral,difference,ratio]=gauss_quad(f_data,quad_data,n_vec)
%
%  This uses Gaussian quadrature with n node points, to  
%  integrate the function f over the interval [a,b].  The 
%  values of n used are taken from the input vector n_vec.
%  The corresponding numerical integrals are returned in the 
%  vector integral.  The differences of successive numerical 
%  integrals are returned in the vector difference:
%        difference(i) = integral(i)-integral(i-1),  i=2,...,doubles
%  The entries in ratio give the rate of decrease in these
%  differences.
%
%  In using this program, define the integrand using the 
%  function given below.  

% Set problem parameters
n0 = quad_data(1);
doubles = quad_data(2);

%  Initialize output vectors.
integral = zeros(doubles+1,1);
difference = zeros(doubles+1,1);
ratio = zeros(doubles+1,1);

% Do the Gaussian quadrature.
for j=0:doubles
    n = n_vec(j+1);
    integral(j+1) = gaussint(n,f_data);
end

%  Calculate the differences of the successive 
%  Gaussian quadrature integrals and the ratio
%  of decrease in these differences.
if doubles > 0
    difference(2:doubles+1) = integral(2:doubles+1)-integral(1:doubles);
end
if doubles > 1
    for i=3:doubles+1
        if difference(i) ~= 0 
            ratio(i) = difference(i-1)/difference(i);
        else
            ratio(i) = inf;
        end
    end
end

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

function val=gaussint(n,f_data)

% [val,bp,wf]=gaussint(fun,a,b,n) integrates
% a function from a to b using an n-point
% Gauss rule which is exact for a polynomial
% of degree 2*n-1. Concepts on page 93 of
% 'Methods of Numerical Integration' by
% Philip Davis and Philip Rabinowitz yield
% the base points and weight factors.
%
% a,b   -  integration limits
% n     -  order of formula 
% val   -  value of the integral
% bp,wf -  Gauss base points and weight factors on [-1,1]

a = f_data(2);
b = f_data(3);

u = (1:n-1)./sqrt((2*(1:n-1)).^2-1);
[vc,bp] = eig(diag(u,-1)+diag(u,1));
[bp,k] = sort(diag(bp));
wf = 2*vc(1,k)'.^2; 
x = (a+b)/2+((b-a)/2)*bp;
f = fcn(x,f_data)*(b-a)/2; 
val = wf(:)'*f(:);

% ========================================================
% ========================================================
% ========================================================
%
% This is the place to add new functions to the GUI.  You must
% also make a change in the "ChooseIntegrand" portion of the fig file, 
% adding the name of the new function to the list shown when 
% the table is activated.
% --------------------------------------------------------------------

function fcn_val = fcn(x,f_data)

indexfcn = f_data(1);
a = f_data(2);
b = f_data(3);
c1 = f_data(4);
c2 = f_data(5);

switch indexfcn
case 1  % f(x) = sin(c1*x)
    fcn_val = sin(c1*x);
case 2  % f(x) = cos(c1*x)
    fcn_val = cos(c1*x);
case 3  % f(x) = exp(c1*x)
    fcn_val = exp(c1*x);
case 4  % f(x) = abs(x)^c1 : Need c1 > -1 unless a >= 0.
    if c1 <= -1 & ((a <= 0 & b >= 0))
        disp('Error in f(x) = abs(x)^c1 : Need c1 > -1')
        fcn_val = zeros(size(x));
        return
    end
    L = length(x);
    fcn_val = zeros(size(x));
    for i=1:L
        if c1 < 0 & x(i) == 0
            fcn_val(i) = 0;
        else
            fcn_val(i) = abs(x(i))^c1;
        end
    end
case 5  % f(x) = 1/(c1+sin(x))
    if c1 <= 1
        disp('Error in f(x) = 1/(c1+sin(x)) : Need c1 > 1')
        fcn_val = zeros(size(x));
        return
    end
    fcn_val = 1./(c1+sin(x));
case 6  % f(x) = exp(c1*x)*cos(c2*x)
    fcn_val = exp(c1*x).*cos(c2*x);
case 7  % f(x) = 1/(c1 + (x-c2)^2)
    if c1 <= 0
        disp('Error in f(x) = 1/(c1 + (x-c2)^2) : Need c1 > 0')
        fcn_val = zeros(size(x));
        return
    end
    fcn_val = 1./(c1 + (x-c2).^2);
end

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

function ans = true_value(f_data)

indexfcn = f_data(1);
a = f_data(2);
b = f_data(3);
c1 = f_data(4);
c2 = f_data(5);

switch indexfcn
case 1  % f(x) = sin(c1*x)
    ans = (cos(c1*a)-cos(c1*b))/c1;
case 2  % f(x) = cos(x)
    ans = (sin(c1*b)-sin(c1*a))/c1;
case 3  % f(x) = exp(c1*x)
    ans = (exp(c1*b)-exp(c1*a))/c1;
case 4  % f(x) = abs(x)^c1 : Use only on [a,b] with a >= 0.
    if b < 0
        ans = ((-a)^(c1+1) - (-b)^(c1+1))/(c1+1);
    elseif a >= 0
        ans = (b^(c1+1) - a^(c1+1))/(c1+1);
    else
        ans = (b^(c1+1)+(-a)^(c1+1))/(c1+1);
    end
case 5  % f(x) = 1/(c1+sin(x))
    % Check whether b-a is longer than 2*pi.  If so, reduce the length of
    % the interval to have the length be less than 2*pi.  Later add in an
    % appropriate amount to compensate the integral.
    k = floor((b-a)/(2*pi) + 10*eps);
    b_mod = b - 2*pi*k;
    if -pi-10*eps <= a & b_mod <= pi+10*eps
        a_mod = a;
        neg = 0;
    elseif b_mod > pi+10*eps
        a_mod = b_mod - 2*pi;
        b_mod = a;
        neg = 1;
    elseif a < -pi-10*eps
        a_mod = b_mod;
        b_mod = a + 2*pi;
        neg = 1;
    end
     temp1 = GFcn(b_mod,c1) - GFcn(a_mod,c1);
     if neg == 1 | k > 0
        temp2 = 2*pi/sqrt(c1*c1 - 1);
    end
    if neg == 1
        temp1 = temp2 - temp1;
    end
    if k > 0
        ans = temp1 + k*temp2;
    else
        ans = temp1;
    end
case 6  % f(x) = exp(c1*x)*cos(c2*x)
    temp1 = (c1*cos(b*c2)+c2*sin(b*c2))*exp(b*c1);
    temp2 = (c1*cos(a*c2)+c2*sin(a*c2))*exp(a*c1);
    ans = (temp1-temp2)/(c1^2+c2^2);
case 7  % f(x) = 1/(c1 + (x-c2)^2)
    con = sqrt(c1);
    temp1 = atan((b-c2)/con);
    temp2 = atan((a-c2)/con);
    ans = (1/con)*(temp1 - temp2);
end

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

function ans = GFcn(x,c)
%
% This evaluates the anti-derivative of the integrand 1/(c1+sin(x)).  
% It is restricted to the case that -pi <= x <= pi. Also, the
% parameter c is required to satisfy c > 1.
%
temp1 = sqrt(c*c-1);
temp2 = atan(1/temp1);
if -pi+10*eps < x & x < pi-10*eps
    temp3 = c*tan(x/2) + 1;
    temp4 = atan(temp3/temp1);
    ans   = (2/temp1)*(temp4 - temp2);
elseif abs(x+pi) < 20*eps
    ans = (-2/temp1)*(pi/2 + temp2);
elseif abs(x-pi) < 20*eps
    ans = (2/temp1)*(pi/2 - temp2);
else
    disp('GFcn: x is outside of [-pi,pi]')
end