function area = riemannsum(f, a, b, n, method, c, varargin)
%RIEMANNSUM computes and graphs the Riemann sum of an integrable function 
%   f, on an interval [a,b] using n equal rectangular polygons.
%
% Inputs:
%   f           a handle to an integrable function
%   a           left end point of interval of integration 
%   b           right end point of interval of integration
%   n           the number of rectangles
%   method      the method for computing the riemann sum is a single character
%               string chosen from 'l','r','m','s','t' for 
%               'left sum', 'right sum', 'middle sum', 'random', 'trapezoid'
%               option 's' chooses a pseudo-random value in each interval of
%               the partition of [a,b]. 
%   c           color must be a single character string chosen from the list 
%               'r','g','b', 'c','m','y','w','k', or 
%               an RGB row vector triple, [r g b], the rectangles are 
%               filled with the constant specified color.
%   varargin    the variable argument list modifies one or more Name,Value
%               argument pairs to set the patch properties of the filled
%               rectangles of the graph. For example
%               FaceColor - this is the fill color you set with C
%               FaceAlpha - alpha blending the faces for transparency [0,1]
%               EdgeColor - an RGB row vector triple line color specification 
%               LineStyle - for plotting edges
%
% Outputs:
%   area    the area between the curve and the x-axis
%
% Example: 
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%
% See also: polyval, isreal, isnumeric, zeros, error, sprintf
% 
% Author: Instructor
% Math 451 
% student id: 
% email: 
% Website: https://blackboard.wichita.edu/
% 14-Sept-2022; Last revision: 16-Sept-2022         
%------------- BEGIN CODE --------------

    if nargin < 4, n = 100; end
    if nargin < 5, method = 's'; end
    if nargin < 6, c = []; end

    switch true
        case ~isa(f, 'function_handle')
            error('riemannsum:InvalidInput', 'function handle invalid')
            
        case ~isa([a, b], 'double')
            error('riemannsum:InvalidInput', 'an end point is invalid')
            
        case ~isfinite([a, b]) 
            error('riemannsum:InvalidInput', 'end points must be finite')
            
        case ~isfinite(n) || ~(floor(n) == n) || ~(n > 0)
            error('riemannsum:InvalidInput', 'number of subintervals must be a positive integer')
            
        case ~isa(method, 'char') 
            error('riemannsum:InvalidInput', 'method option must be a single character')
            
        case ~isempty(c) && ~isa(c, 'char') && ~isa(c, 'double') && ~(length(c) == 3)
            error('riemannsum:InvalidInput', 'color option invalid; see function header')
            
        otherwise
            
            
            area = GraphRiemannSum(f, a, b, n, method, c, varargin{:});
    end
end

function area = GraphRiemannSum(f, a, b, n, method, c, varargin)
%GraphRiemannSum  computes and graphs the Riemann sum of an integrable function 
%   f, on an interval [a,b] using n equal rectangular polygons.
%
% Inputs:
%   f           a handle to an integrable function
%   a           left end point of interval of integration 
%   b           right end point of interval of integration
%   n           the number of rectangles
%   method      the method for computing the riemann sum is a single character
%               string chosen from 'l','r','m','s','t' for 
%               'left sum', 'right sum', 'middle sum', 'random', 'trapezoid'
%               option 's' chooses a pseudo-random value in each interval of
%               the partition of [a,b]. 
%   c           color must be a single character string chosen from the list 
%               'r','g','b', 'c','m','y','w','k', or 
%               an RGB row vector triple, [r g b], the rectangles are 
%               filled with the constant specified color.
%   varargin    the variable argument list modifies one or more Name,Value
%               argument pairs to set the patch properties of the filled
%               rectangles of the graph. For example
%               FaceColor - this is the fill color you set with C
%               FaceAlpha - alpha blending the faces for transparency [0,1]
%               EdgeColor - an RGB row vector triple line color specification 
%               LineStyle - for plotting edges
%
% Outputs:
%   area    the area between the curve and the x-axis
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%
% See also: polyval, isreal, isnumeric, zeros, error, sprintf
% 
% Author: Instructor
% Math 451 
% student id: 
% email: 
% Website: https://blackboard.wichita.edu/
% 14-Sept-2022; Last revision: 16-Sept-2022         
%------------- BEGIN CODE --------------

    [area, X, Y] = ComputeRiemannSum(f, a, b, n, method);

    % if we need to create a graph
    if ~isempty(c)
        t = linspace(a, b, n*41);
        y = f(t);
        
        axes = [a b min(y) max(y)];
        axis(axes);
        hold on

        for i=1:n
            fill(X(i,:), Y(i,:), c, varargin{:});
        end

        plot(t, y, 'k')
        hold off
    end
end
        
function [s, X, Y] = ComputeRiemannSum(f, a, b, n, method)
%ComputeRiemannSum  computes the Riemann sum of an integrable function 
%   f, on an interval [a,b] using n equal rectangular polygons.
%
% Inputs:
%   f           a handle to an integrable function
%   a           left end point of interval of integration 
%   b           right end point of interval of integration
%   n           the number of rectangles
%   method      the method for computing the riemann sum is a single character
%               string chosen from 'l','r','m','s','t' for 
%               'left sum', 'right sum', 'middle sum', 'random', 'trapezoid'
%               option 's' chooses a pseudo-random value in each interval of
%               the partition of [a,b]. 
%
% Outputs:
%   s           the area between the curve and the x-axis
%   X           patch vertices for graphing trapezoids
%   Y           patch vertices for graphing trapezoids
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%
% See also: polyval, isreal, isnumeric, zeros, error, sprintf
% 
% Author: Instructor
% Math 451 
% student id: 
% email: 
% Website: https://blackboard.wichita.edu/
% 14-Sept-2022; Last revision: 16-Sept-2022         
%------------- BEGIN CODE --------------

    % partition of interval [a,b]
    x = linspace(a, b, n+1);
    
    % Patch vertices
    X = [ x(1:n); x(2:n+1); x(2:n+1); x(1:n) ].';

    deltax = (b - a)/n;
    
    switch lower(method)
        case 'l' 
            % values at left end points
            x = x(1:end-1);

        case 'r'
            % values at right end points
            x = x(2:end);

        case 'm'
            % values at midpoints
            x = (x(1:end-1) + x(2:end))/2;

        case 's'
            % uniform random value in each interval [x(k), x(k+1)]
            x = x(1:end-1) + rand(1,n)*deltax;

        case 't'
            % trapezoid rule needs all n+1 values of y = f(x)

        otherwise
            % error unrecognized method
            error('riemannsum:InvalidInput', 'invalid method option %s', method)
    end

    z = zeros(1,n);
    y = f(x);
    
    % area
    if method == 't'
        % trapezoid rule
        s = deltax * y * [1; 2*ones(n-1,1); 1] / 2;
        
        % Patch vertices
        Y = [ z; z; y(2:n+1); y(1:n)].';
    else
        s = deltax * sum(y);

        % Patch vertices
        Y = [ z; z; y; y ].';
    end
end
