function integral = simpson_integration (fhandle, a, b, num_segments)

%SIMPSON_INTEGRATION  Computes a definite integral
%   INTEGRAL = SIMPSON_INTEGRATION ( FHANDLE, A, B, NUM_SEGMENTS )
%   determines the integral of the function defined by FHANDLE between the
%   bounds of A and B. The integration interval [A,B] is divided into
%   NUM_SEGMENTS segments of equal width and the multiple application
%   Simpson 1/3 integration rule is applied.  If the number of segments is
%   odd, then the last three segments are integrated by use of the Simpson
%   3/8 rule.
%
%   Inputs:
%      FHANDLE:   a Matlab function handle for the function to be
%                 integrated.
%      A:         Lower bound of the integration interval
%      B:         Upper bound of the integration interval
%      NUM_SEGMENTS: the number of segments into which the interval is
%                 divided to approximate the integral
%
%   Outputs:
%      INTEGRAL:  The result of the integration.
%
%   Assumptions:
%      - The function is continuous and smooth.
%      - the function pointed to by FHANDLE supports vector operations
%
%% ME 2016 Section
%% Fall 2007
%%
%% AUTHOR:  Chris Paredis
%%

% error checking:
% make sure the number of segments is a whole number larger than 1
num_segments = round(num_segments);
if num_segments<2
    error('the number of segments must be larger than or equal to 2');
end

% compute the function values.
% Note: the number of function values is num_segments+1
x_values = linspace(a,b,num_segments+1);
f_values = fhandle(x_values);
h = (b-a)/num_segments; % segment width

% compute the integral
if rem(num_segments,2) == 0
    % an even number of segments --> Simpson 1/3 can be applied
    integral = h/3 * (f_values(1) + 4*sum(f_values(2:2:end-1)) + ...
                      2*sum(f_values(3:2:end-2)) + f_values(end));
else
    % Simpson 1/3 applied to all but the last three segments
    integral13 = h/3 * (f_values(1) + 4*sum(f_values(2:2:end-4)) + ...
                        2*sum(f_values(3:2:end-5)) + f_values(end-3));
    % Simpson 3/8 for the last three segments
    integral38 = h*3/8 * (f_values(end-3) + ...
                          3*(f_values(end-2) + f_values(end-1)) + ...
                          f_values(end));
    % The total integral is the sum of the the two parts
    integral = integral13 + integral38;
end