Mathematica/3DContours

Mathematica supports 3D contours on surfaces using the option MeshFunctions. MeshFunctions draws ISO lines of an arbitary function of the x,y,z values, and for some functions additional parameters. For a conventional set of altitude contours we just use the identity function of the z value (the third parameter) using MeshFunctions->{#3&} eg Plot3D[Sin[x y], {x, 0, 3}, {y, 0, 3}, MeshFunctions -> {#3 &}] In addition the region between the contours can be colored using the option MeshShading. The following example divides the image into 10 altitude contours and colors them in gray scale. Plot3D[Sin[x y], {x, 0, 3}, {y, 0, 3}, Mesh -> 10, MeshFunctions -> {#3 &}, MeshShading -> GrayLevel /@ Range[0, 1, 0.1], Lighting -> "Neutral"]