function area = riemann( f, a, b, n )
%RIEMANN Integrate f from a to b as a Riemann sum of n plotted rectangles.
% (Uses left rectangles, but that doesn't really matter.)


% Set up x and f(x) values. Don't assume f vectorizes.

x = linspace( a, b, n+1 );

y = [];
for i = 1 : n
    y(i) = f( x(i) );
end


% Calculate the area under the curve

deltaX = ( b - a ) / n;
area = sum( deltaX * y );


% Plot the rectangles by building a list of (x,y) coordinates for each
% corner of each rectangle, then plotting a line connecting those coordinates

rectX = [];
rectY = [];

for i = 1 : n
    rectX = [ rectX x(i) x(i) x(i+1) x(i+1) ];
    rectY = [ rectY 0 y(i) y(i) 0 ];
end

clf;
plot( rectX, rectY, 'r' );
hold on;


% Plot the function using more densely spaced samples. Again, don't
% assume the function vectorizes.

nSamples = 200;

funcX = linspace( a, b, nSamples );
funcY = [];

for i = 1 : nSamples
    funcY(i) = f( funcX(i) );
end

plot( funcX, funcY, 'b' );

end

