Double Riemann Sums

Code

graphics[xSubdivisions_,ySubdivisions_,xPoint_,yPoint_,xLeft_,xRight_,yLeft_,yRight_,showSurface_,showPoints_]:=(
deltaX=(xRight-xLeft)/xSubdivisions;
deltaY=(yRight-yLeft)/ySubdivisions;
points=Join@@Table[{x,y,f[x,y]},{x,xLeft+deltaX*xPoint,xRight,deltaX},{y,yLeft+deltaY*yPoint,yRight,deltaY}];
riemannSum=Total[#[[3]]&/@points]*deltaX*deltaY;
definiteIntegral=NIntegrate[f[x,y],{x,xLeft,xRight},{y,yLeft,yRight}];

riemannPlot=DiscretePlot3D[f[x,y],{x,xLeft+deltaX*xPoint,xRight,deltaX},{y,yLeft+deltaY*yPoint,yRight,deltaY},ExtentSize->{{deltaX*xPoint,deltaX*(1-xPoint)},{deltaY*yPoint,deltaY*(1-yPoint)}}];
surface=If[showSurface,Plot3D[f[x,y],{x,xLeft,xRight},{y,yLeft,yRight},Mesh->None,PlotStyle->Opacity[0.75]],{}];
pointsPlot=If[showPoints,Graphics3D[{Green,PointSize[0.02],Point[points]}],{}];
Show[riemannPlot,surface,pointsPlot,Axes->Automatic,AxesLabel->{x,y,z},BoxRatios->{2,2,3},ViewPoint->{20,4,8},PlotRange->{{xLeft,xRight},{yLeft,yRight},Automatic},PlotLabel->Style[Row[{"Riemann Sum: ",riemannSum,"\nDefinite Integral: ",definiteIntegral}],FontSize->16]])

Graphics

Modify the function and options below to change the Riemann sum.
xPoint specifies where in each x interval to take the sample point; a left sum corresponds to xPoint=0, a right sum to xPoint=1, and a midpoint to xPoint=0.5.  (And similarly for yPoint.)
showPoints specifies whether to draw the sample points on the surface.

f[x_,y_]:=18-x^2-y^2;
xSubdivisions=6;
ySubdivisions=8;
xPoint=0.5;
yPoint=0.5;
xLeft=0;
xRight=3;
yLeft=-1;
yRight=3;
showSurface=True;
showPoints=True;
graphics[xSubdivisions,ySubdivisions,xPoint,yPoint,xLeft,xRight,yLeft,yRight,showSurface,showPoints]

Animation

This animation shows the Riemann sum as the number of subdivisions in each direction increases from startingSubdivions to endingSubdivisions.
All options except xSubdivisions and ySubdivions are taken from the Graphics section above.
It may take some time to run this command, as all frames of the animation are precomputed so that it can run smoothly.

startingSubdivisions=1;
endingSubdivisions=24;

frames=Table[Image[graphics[t,t,xPoint,yPoint,xLeft,xRight,yLeft,yRight,showSurface,showPoints]],{t,startingSubdivisions,endingSubdivisions,1}];
ListAnimate[frames]