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

% Last Modified by GUIDE v2.5 20-Sep-2004 17:56:20

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 = 4;
    handles.const_double = 9;
    handles.quad_data = [handles.const_n_init,handles.const_double]; 
    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

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

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

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

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);
Integrate_GUI('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))


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

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

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);
Integrate_GUI('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

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

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);
Integrate_GUI('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

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

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

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);
Integrate_GUI('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

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

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

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

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

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

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];
guidata(h,handles);
Integrate_GUI('ChooseIntegrand_Callback',handles.ChooseIntegrand,[],guidata(handles.ChooseIntegrand))

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

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

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

% --------------------------------------------------------------------
% --- Executes on button press in ErrorEstimates.
function ErrorEstimates_Callback(hObject, eventdata, handles)
% hObject    handle to ErrorEstimates (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


ListErrorData(handles);

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

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','integrate_gui_graph') 
close(h_invisible);
disp('The graph is stored under the name "integrate_gui_graph.eps"')

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

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

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

switch handles.method
case 'Trapezoidal'
    fid = fopen('output_Trapezoidal','w');
    disp('The table is saved under the name "output_Trapezoidal"')
case 'Midpoint'
    fid = fopen('output_Midpoint','w');
    disp('The table is saved under the name "output_Midpoint"')
case 'Simpson'
    fid = fopen('output_Simpson','w');
    disp('The table is saved under the name "output_Simpson"')
case 'Boole'
    fid = fopen('output_Boole','w');
    disp('The table is saved under the name "output_Boole"')
end

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

Integrate_GUI_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 varargout = TrapezoidalPushbutton_Callback(h, eventdata, handles, varargin)

handles.method = 'Trapezoidal';
guidata(h,handles);
Calculate_Integral(h,handles);

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

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

handles.method = 'Midpoint';
guidata(h,handles);

Calculate_Integral(h,handles);

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

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

handles.method = 'Simpson';
guidata(h,handles);

Calculate_Integral(h,handles);

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

% --- Executes on button press in BoolePushbutton.
function BoolePushbutton_Callback(h, eventdata, handles, varargin)

handles.method = 'Boole';
guidata(h,handles);

Calculate_Integral(h,handles);

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

function Calculate_Integral(h,handles)

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

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

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

switch handles.method
case 'Trapezoidal'
    [integral,difference,ratio_diff]=trapezoidal(f_data,quad_data);
    celldata(1) = cellstr('                          Trapezoidal method output');
case 'Midpoint'
    [integral,difference,ratio_diff]=midpoint(f_data,quad_data);
    celldata(1) = cellstr('                          Midpoint method output');
case 'Simpson'
    [integral,difference,ratio_diff]=simpson(f_data,quad_data);
    celldata(1) = cellstr('                          Simpson''s method output');
case 'Boole'
    [integral,difference,ratio_diff]=boole(f_data,quad_data);
    celldata(1) = cellstr('                          Boole''s method output');
end

% 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_init*2^j;
    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);

% Do asymptotic error results
error_rich = zeros(doubles+1,1);
error_asym = zeros(doubles+1,1);
celldata_error = cellstr(num2str(zeros(doubles+4,1)));

n = n_init*(2.^[0:doubles])';
stepsize = (b-a)./n;

switch handles.method
case 'Trapezoidal'
    error_rich = (1/3)*difference;
    constant = -(dfcn(b,1,f_data)-dfcn(a,1,f_data))/12;
    error_asymp = constant*stepsize.^2;
    celldata_error(1) = cellstr('                                    Trapezoidal method output');
case 'Midpoint'
    error_rich = (1/3)*difference;
    constant = (dfcn(b,1,f_data)-dfcn(a,1,f_data))/24;
    error_asymp = constant*stepsize.^2;
    celldata_error(1) = cellstr('                                    Midpoint method output');
case 'Simpson'
    error_rich = (1/15)*difference;
    constant = -(dfcn(b,3,f_data)-dfcn(a,3,f_data))/180;
    error_asymp = constant*stepsize.^4;
    celldata_error(1) = cellstr('                                    Simpson''s method output');
case 'Boole'
    error_rich = (1/63)*difference;
    constant = -(2/945)*(dfcn(b,5,f_data)-dfcn(a,5,f_data));
    error_asymp = constant*stepsize.^6;
    celldata_error(1) = cellstr('                                    Boole''s method output');
end

celldata_error(2) = cellstr('    n     True       Richardson       Richardson         R. Extrap.   Asymptotic       Asymptotic         A. Extrap.');
celldata_error(3) = cellstr('          Error        Error         Extrapolation         Error        Error         Extrapolation         Error');
extrapolation_rich = integral + error_rich;
extrapolation_asymp = integral + error_asymp;
error_R = true_integral - extrapolation_rich;
error_A = true_integral - extrapolation_asymp;

for j=0:doubles
    celldata_error(j+4) = cellstr(sprintf('%5i%13.3e%13.3e%23.14e%13.3e%13.3e%23.14e%13.3e',...
       n(j+1),true_error(j+1),error_rich(j+1),extrapolation_rich(j+1),error_R(j+1),...
              error_asymp(j+1),extrapolation_asymp(j+1),error_A(j+1)));
end

StringData_Error = char(celldata_error);
handles.StringData_Error = StringData_Error;
guidata(h,handles);

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

function [integral,difference,ratio]=trapezoidal(f_data,quad_data)
%
%  function [integral,difference,ratio]=trapezoidal(f_data,quad_data)
%
%  This uses the trapezoidal rule with n subdivisions to  
%  integrate the function f over the interval [a,b].  The 
%  values of n used are
%        n = n0,2*n0,4*n0,...,(2^doubles)*n0
%  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.  The parameter index_f allows the user 
%  to do calculations with multiple integrands.

% Set problem parameters
indexfcn = f_data(1);
a = f_data(2);
b = f_data(3);
c1 = f_data(4);
c2 = f_data(5);

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

% Initialize for trapezoidal rule.
sumend = (fcn(a,f_data) +fcn(b,f_data))/2;
h = (b-a)/n0;

%  Initialize for case of n0 > 1.
if n0 > 1
    x = a + h*[1:n0-1];
    accum_sum = sum(fcn(x,f_data));
else
    accum_sum = 0;
end

integral(1) = h*(sumend + accum_sum);

%  Calculate the numerical integrals, doing each
%  by appropriately modifying the preceding case.
for i=1:doubles
    n = n0*(2^i);
    h = (b-a)/n;
    x = a + h*[1:2:n-1];
    accum_sum = accum_sum + sum(fcn(x,f_data)); 
    integral(i+1) = h*(sumend + accum_sum);
end

%  Calculate the differences of the successive 
%  trapezoidal rule 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 [integral,difference,ratio]=midpoint(f_data,quad_data)
%
%  function [integral,difference,ratio]=midpoint(f_data,quad_data)
%
%  This uses the midpoint rule with n subdivisions to  
%  integrate the function f over the interval [a,b].  The 
%  values of n used are
%        n = n0,2*n0,4*n0,...,(2^doubles)*n0
%  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.  The parameter index_f allows the user 
%  to do calculations with multiple integrands.

% Set problem parameters
indexfcn = f_data(1);
a = f_data(2);
b = f_data(3);
c1 = f_data(4);
c2 = f_data(5);

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

%  Calculate the numerical integrals, doing each
%  by appropriately modifying the preceding case.
for i=1:doubles+1
    n = n0*(2^(i-1));
    h = (b-a)/n;
    x = a:h:b-h;
    x = x + (h/2);
    integral(i) = h*sum(fcn(x,f_data));
end

%  Calculate the differences of the successive 
%  trapezoidal rule 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 [integral,difference,ratio]=simpson(f_data,quad_data)
%
%  function [integral,difference,ratio]=simpson(f_data,quad_data)
%
%  This uses Simpson's rule with n subdivisions to  
%  integrate the function f over the interval [a,b].  The 
%  values of n used are
%        n = n0,2*n0,4*n0,...,(2^doubles)*n0
%  The initial value of n0 must be a positive even integer.
%  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.  The parameter index_f allows the user 
%  to do calculations with multiple integrands.

% Set problem parameters

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

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

if n0 ~= 2*floor(n0/2)
    disp('Simpson''s rule error: n_init must be even and >= 2.')
    return
end

% Initialize for Simpson's rule.
sum_endpts = fcn(a,f_data) +fcn(b,f_data);
h = (b-a)/n0;

if n0 > 2
    x = a + h*[2:2:n0-2];
    sum_even = sum(fcn(x,f_data));
else
    sum_even = 0;
end

%  Calculate the numerical integrals, doing each
%  by appropriately modifying the preceding case.
for i=1:doubles+1
    n = n0*(2^(i-1));
    h = (b-a)/n;
    x = a + h*[1:2:n-1];
    sum_odd = sum(fcn(x,f_data)); 
    integral(i) = (h/3)*(sum_endpts + 2*sum_even +4*sum_odd);
    sum_even = sum_even + sum_odd;
end

%  Calculate the differences of the successive 
%  trapezoidal rule 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 [integral,difference,ratio]=boole(f_data,quad_data)
%
%  function [integral,difference,ratio]=boole(f_data,quad_data)
%
%  This uses Boole's rule with n subdivisions to  
%  integrate the function f over the interval [a,b].  The 
%  values of n used are
%        n = n0,2*n0,4*n0,...,(2^doubles)*n0
%  The initial value n0 must be a positive integral multiple of 4.
%  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.  The parameter index_f allows the user 
%  to do calculations with multiple integrands.

% Set problem parameters

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

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

if n0 ~= 4*floor(n0/4)
    disp('Boole''s rule error: n_init must be a multiple of 4.')
    return
end

% Initialize for Boole's rule.
sum_endpts = fcn(a,f_data) +fcn(b,f_data);
h = (b-a)/n0;

if n0 == 4
    sum_2 = fcn(a+2*h,f_data);
    sum_4 = 0;
elseif n0 > 4
    x = a + h*[2:4:n0-2];
    sum_2 = sum(fcn(x,f_data));
    x = a + h*[4:4:n0-4];
    sum_4 = sum(fcn(x,f_data));
end

%  Calculate the numerical integrals, doing each
%  by appropriately modifying the preceding case.
for i=1:doubles+1
    n = n0*(2^(i-1));
    h = (b-a)/n;
    x = a + h*[1:2:n-1];
    sum_odd = sum(fcn(x,f_data)); 
    integral(i) = (2*h/45)*(7*sum_endpts + 12*sum_2 + 14*sum_4 +32*sum_odd);
    sum_4 = sum_4 + sum_2;
    sum_2 = sum_odd;
end

%  Calculate the differences of the successive 
%  trapezoidal rule 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


% ========================================================
% ========================================================
% ========================================================
%
% 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 dfcn_val = dfcn(x,order,f_data)

% Calculate the derivative of order "order".
% It is assumed that "x" is a simple variable, not an array.

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)
    switch order
    case 1
        dfcn_val = c1*cos(c1*x);
    case 3
        dfcn_val = -(c1^3)*cos(c1*x);
    case 5
        dfcn_val = (c1^5)*cos(c1*x);
    end
case 2  % f(x) = cos(c1*x)
    switch order
    case 1
        dfcn_val = -c1*sin(c1*x);
    case 3
        dfcn_val = (c1^3)*sin(c1*x);
    case 5
        dfcn_val = -(c1^5)*sin(c1*x);
    end
case 3  % f(x) = exp(c1*x)
    switch order
    case 1
        dfcn_val = c1*exp(c1*x);
    case 3
        dfcn_val = (c1^3)*exp(c1*x);
    case 5
        dfcn_val = (c1^5)*exp(c1*x);
    end
case 4  % f(x) = abs(x)^c1 : Need c1 > -1 unless a >= 0.
    abs_x = abs(x);
    sign_x = sign(x);
    switch order
    case 1
        if c1 < 1 & x == 0
            dfcn_val = inf;
        else
            dfcn_val = c1*sign_x*abs_x^(c1-1);
        end
    case 3
        if c1 < 3 & x == 0
            dfcn_val = inf;
        else
            dfcn_val = c1*(c1-1)*(c1-2)*sign_x*abs_x^(c1-3);
        end
    case 5
        if c1 < 5 & x == 0
            dfcn_val = inf;
        else
            dfcn_val = c1*(c1-1)*(c1-2)*(c1-3)*(c1-4)*sign_x*abs_x^(c1-5);
        end
    end
case 5  % f(x) = 1/(c1+sin(x))
    cs = cos(x); sn = sin(x);
    d = c1 + sn;
    switch order
    case 1
        dfcn_val = -cs/d^2;;
    case 3
        dfcn_val = (cs*c1^2 - 6*cs^3 - 5*cs*sn^2 - 4*cs*sn*c1)/d^4;
    case 5
        term1 = -cs/d^2;
        term2 = 30*cs*sn/d^3;
        term3 = (60*cs^3 - 90*cs*sn^2)/d^4;
        term4 = -240*(cs^3)*sn/d^5;
        term5 = -120*(cs^5)/d^6;
        dfcn_val = term1+term2+term3+term4+term5;
    end
case 6  % f(x) = exp(c1*x)*cos(c2*x)
    cs = cos(c2*x); sn = sin(c2*x); ex = exp(c1*x);
    switch order
    case 1
        dfcn_val = ex*(c1*cs-c2*sn);
    case 3
        d1 = c1^3 - 3*c1*c2^2; 
        d2 = 3*c2*c1^2 - c2^3;
        dfcn_val = ex*(d1*cs-d2*sn);
    case 5
        d1 = c1^3 - 3*c1*c2^2; 
        d2 = 3*c2*c1^2 - c2^3;
        e1 = d1*c1^2 - 2*c1*c2*d2 - d1*c2^2;
        e2 = d2*c1^2 + 2*c1*c2*d1 -d2*c2^2;
        dfcn_val = ex*(e1*cs-e2*sn);
    end
case 7  % f(x) = 1/(c1 + (x-c2)^2)
    t = x-c2;
    d = c1 + (x-c2)^2;
    switch order
    case 1
        dfcn_val = -2*t/d^2;
    case 3
        dfcn_val = 24*t*(c1 - t^2)/d^4;
    case 5
        dfcn_val = 240*t*(10*c1*t^2 - 3*c1^2 - 3*t^4)/d^6;
    end
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