(*warning messages*)

riemannSum3D::badarg="Something is wrong with your input. Check it carefully
	and try again.";

riemannSum3D::badHeight="Something is wrong with the Height specification. Check
it carefully and try again.";

(*options*)

Options[riemannSum3D]={
Height->Midpoint,
XSubintervals->10,
YSubintervals->10,
Printout->True,
ZMinimum->0,
TopsOnly->False
};

(*definition of auxiliary functions*)

(* make expr a list if it isn't already *)
makeList[expr_]:=Switch[expr,_List,expr,_,{expr}]

(* Add the x-coordinate "x" to each of the y or{y,z} coordinates in "coords" list *)
addX[x_,coords_List]:=Map[Prepend[makeList[#],x]&,coords]

(* Add the y-coordinate "y" to each of the x or{x,z} coordinates in "coords" list *)
addY[y_,coords_List]:=Map[Insert[makeList[#],y,2]&,coords]

(* Add the z-coordinate "z" to each of the {x,y} coordinates in "coords" list *)
addZ[z_,coords_List]:=Map[Append[makeList[#],z]&,coords]

(*definition of riemannSum3D*)

Clear[riemannSum3D];riemannSum3D[f_,{x_,xmin_?NumericQ,xmax_?NumericQ,delx_:Automatic},{y_,ymin_?NumericQ,ymax_?NumericQ,dely_:Automatic},opts___Rule]:=Module[{dx=delx,dy=dely,xvals,yvals,zpts,xn,yn,a,b,xtops,ytops,xycoords,xyzcoords,zmin,grid,displayItems},

{xsubint,ysubint,ht,pntout,zmin,to}={XSubintervals,YSubintervals,Height,Printout,ZMinimum,TopsOnly}/.{opts}/.Options[riemannSum3D];

(* xvals are the subinterval endpoints on the x axis*)
xvals=Which[NumericQ[dx]&&xsubint==XSubintervals/.Options[riemannSum3D],Union[N@Range[xmin,xmax,dx],{N@xmax}],
NumericQ[dx]&&IntegerQ[xsubint]&&xsubint>=1,
Union[N[xmin]+Table[i*dx,{i,0,Min[xsubint-1,(xmax-xmin)/dx]}],{N[xmax]}],
dx==Automatic&&IntegerQ[xsubint]&&xsubint>=1,Range[xmin,xmax,(xmax-xmin)/N[xsubint]],
Head[xsubint]===List,Union[{N[xmin]},N[xsubint],{N[xmax]}],
True,Return[Message[riemannSum3D::badarg]]];

(* yvals are the subinterval endpoints on the y axis*)
yvals=Which[NumericQ[dy]&&ysubint==YSubintervals/.Options[riemannSum3D],Union[N@Range[ymin,ymax,dy],{N@ymax}],
NumericQ[dy]&&IntegerQ[ysubint]&&ysubint>=1,
Union[N[ymin]+Table[i*dy,{i,0,Min[ysubint-1,(ymax-ymin)/dy]}],{N[ymax]}],
dy==Automatic&&IntegerQ[ysubint]&&ysubint>=1,Range[ymin,ymax,(ymax-ymin)/N[ysubint]],
Head[ysubint]===List,Union[{N[ymin]},N[ysubint],{N[ymax]}],
True,Return[Message[riemannSum3D::badarg2]]];

xn=Length[xvals];
yn=Length[yvals];

(* zpts is the single list of values of f that are the heights of the rectangular solids;
The length of zpts is (Length[xvals]-1)(Length[yvals]-1) and the order of the values
illustrated as follows: zpts = {{xmin,ymin},{xmin,ymin+dx},...,{xmin,ymax},{xmin+dx,ymin},...}  *)
{a,b}=Which[ht===Midpoint||ht===MidPoint,{.5,.5},
Head[ht]==List&&Length[ht]==2&&0<=ht[[1]]<=1&&0<=ht[[2]]<=1,ht,
True,Return[Message[riemannSum3D::badHeight]]
];
zpts=f/.Thread[{x,y}->#]&/@Flatten[Outer[List,(1-a) Drop[xvals,-1]+a Drop[xvals,1],(1-b) Drop[yvals,-1]+b Drop[yvals,1]],1];

(* If Printout set to True or Long, print stuff out *)
Which[pntout===True,
Print["Riemann Sum = ",Flatten[Outer[Times,(Drop[xvals,1]-Drop[xvals,-1]),(Drop[yvals,1]-Drop[yvals,-1])]].zpts],pntout===Long,Print["x partition = ",xvals];
Print["y partition = ",yvals];
Print["Each sub-rectangle is sampled at the point ",100 a,"% from left and ",100 b,"% from bottom."];
Print["Rectangle heights = ",Transpose[Reverse/@Partition[zpts,yn-1]]//TableForm];
Print["Riemann Sum = ",Flatten[Outer[Times,(Drop[xvals,1]-Drop[xvals,-1]),(Drop[yvals,1]-Drop[yvals,-1])]].zpts]
];

(* the polygons for the tops of the rectangular solids *)
xtops=(Flatten[Join[#,Reverse[#]]]&)/@Partition[xvals,2,1];
ytops=Apply[Join,Partition[Transpose[{yvals,yvals}],2,1],1];
xycoords=Flatten[Table[Transpose[{xtops[[i]],ytops[[j]]}],{i,xn-1},{j,yn-1}],1];
xyzcoords=MapThread[addZ,{zpts,xycoords},1];
If[to,
(* grid if TopsOnly *)
grid=Join[{Hue[1]},
Line/@Table[addX[xvals[[i]],{{ymin,zmin},{ymax,zmin}}],{i,xn}],
Line/@Table[addY[yvals[[i]],{{xmin,zmin},{xmax,zmin}}],{i,yn}]
];
displayItems={Polygon/@xyzcoords,grid},

(* else: display rectangular solids *)
displayItems=Cuboid[#[[1]],ReplacePart[#[[3]],zmin,3]]&/@xyzcoords
];

(* show it *)
Show[Graphics3D[displayItems],FilterRules[{opts},Options[Graphics3D]],Axes->True,PlotLabel->Style[Row[{StringJoin["Riemann sum = ",ToString[Chop@PaddedForm[Flatten[Outer[Times,(Drop[xvals,1]-Drop[xvals,-1]),(Drop[yvals,1]-Drop[yvals,-1])]].zpts,{6,4},ExponentFunction->(Null&)]]],StringJoin["double integral = ",ToString[PaddedForm[NIntegrate[f,{x,-2,2},{y,-2,2}],{6,4},ExponentFunction->(Null&)]]]},Style["  |  ",GrayLevel[.5]]],"Label",10],ImageSize->340,SphericalRegion->True]
]

Manipulate[
riemannSum3D[f,{x,-2,2},{y,-2,2},XSubintervals->xVal,YSubintervals->yVal,Axes->True,PlotRange->Switch[f,(x^2)-(y^2),{{-2,2},{-2,2},{-4,4}},(x^2)+(y^2),{{-2,2},{-2,2},{0,8}},Cos[(x^2)-(y^2)],{{-2,2},{-2,2},{-1,1}},Sin[(x^2)+(y^2)],{{-2,2},{-2,2},{-1,1}}],BoxRatios->{1,1,1},Printout->False,Height->pt],
{{f,(x^2)+(y^2),"f(x,y) ="},{(x^2)-(y^2)->"(x^2) - (y^2)",(x^2)+(y^2)->"(x^2)+ (y^2)",Cos[(x^2)-(y^2)]->"cos((x^2) - (y^2))",Sin[(x^2)+(y^2)]->"sin((x^2) + (y^2))"},ControlType->PopupMenu,ControlPlacement->Left},{{xVal,10,"subdivision of x range"},{5,10,20,40},ControlPlacement->Left},{{yVal,10,"subdivision of y range"},{5,10,20,40},ControlPlacement->Left},{{pt,{.5,.5},"point in rectangle used for evaluation"},{0,0},{1,1},ImageSize->Tiny},SaveDefinitions->True,AutorunSequencing->{1,2}]

ContourPlot3D[z==x^2+y^2, {x,-2,2}, {y,-2,2},{z,0,10},MeshStyle->None, AxesLabel->Automatic]
ContourPlot3D[z==x^2-y^2, {x,-2,2}, {y,-2,2},{z,-4,4},MeshStyle->None, AxesLabel->Automatic]
ContourPlot3D[z==Cos[x^2-y^2], {x,-2,2}, {y,-2,2},{z,-1,1},MeshStyle->None, AxesLabel->Automatic]
ContourPlot3D[z==Sin[x^2+y^2], {x,-2,2}, {y,-2,2},{z,-1,1},MeshStyle->None, AxesLabel->Automatic]