MATHEMATICA IN ACTION for issue 9-1 of Mathematica in Education and Research

Initialization, Setting of Various Options

Get["ContoursAndFlows`"];

SetOptions[#, AspectRatio->Automatic, Frame->True, Axes->None]& /@ 
   {Plot, Graphics, ParametricPlot};
   
SetOptions[#, ColorOutput->GrayLevel] & /@ 
{Plot3D, ParametricPlot3D, Graphics3D, ContourAndFlowPlot};

SetOptions[#, Boxed -> False] & /@ 
{Plot3D, ParametricPlot3D, Graphics3D, ContourAndFlowPlot};

SetOptions[ContourPlot, ContourShading->False];

$DefaultFont = {"Courier", 8};
Off[Graphics3D::nlist3,ContourGraphics::ctsm, ContourGraphics::ctpnt]

Adaptive Surface Plotting,
Continued

How to use contours and flows to define customized coordinate systems.

by  Ron Goetz and Stan Wagon

Four years ago [Goetz and Wagon 1996], we showed how, in some simple cases, the surface defined by z = f(x, y) could be displayed using a special-purpose coordinate system based on the contours of the surface and the curves (orthogonal trajectories, or flows) perpendicular to them. This idea has its roots in the use of polar coordinates to view certain surfaces, since often the radial lines are in fact contours, which makes the resulting image of the surface very smooth.

The examples we succeeded with back in 1996 had the property that the contour curves could be obtained in a closed form. But we know that Mathematica can generate contour plots and that the data of those plots can be used to set up and numerically solve differential equations for the flows. Thus we have the long-term goal of getting the contour-and-flow idea to work on any equation z = f(x, y). This dream still seems far away, in part because of the varying number and types of critical points a surface can have. But we have succeeded in developing a package that allows one to attack a particular surface, perhaps using some trial and error to get the best choice of contours, and generate what we call a contour-and-flow surface. We think the results are stunning: This is surely the right way to view many surfaces, and we hope that our package and description will inspire readers to generate contour-and-flow renderings of their favorite surfaces.

. Euler’s Mistake

Euler, seeking to generalize the Fermat equation, conjectured that three fourth powers could not sum to a fourth power. He was wrong, as proved by Noam Elkies in 1988. Here is a counterexample to Euler’s conjecture (due to Frye; slightly smaller than that of Elkies).

(((95800^4) + (217519^4) + (414560^4))^(1/4))

Such examples can be viewed as being nontrivial rational points on the surface (x^4) + (y^4) + (z^4) = 1; this surface will be our first example of a contour-and-flow plot, where we will use the functional form \(TraditionalForm
      \`f(x, y) = \(±\((1 - x\^4 - y\^4)\)\^\(1/4\)\)\). Plot3D fails badly because of the steepness at the border. We define f in the usual way, and also on lists, for convenience.

f[x_, y_] := (((1 - (x^4) - (y^4)))^(1/4))
f[v_List]  :=  f  @@  v

Show[
        Block[{$DisplayFunction = 
              Identity}, Plot3D[#*f[x, y], {x, -1, 1}, {y, -1, 1}, 
                PlotPoints -> 
                  31]  &  /@  {1, -1}]];

Code for figure

\(\(Show[
        Block[{$DisplayFunction = 
              Identity}, {\(Plot3D[#\ f[x, y], {x, \(-1\), 1}, {y, \(-1\), 
                    1}, PlotPoints -> 31]\  &\)\  /@ \ {1, \(-1\)}}], 
        AxesEdge -> {{\(-1\), \(-1\)}, {1, \(-1\)}, {1, 1}}, 
	
        BoxRatios -> {1, 1, 1}, 
        Ticks -> {{\(-1\), 0, 1}, {\(-1\), 0.5}, {\(-0.5\), 
              1}}];\)\)

The problem is that Plot3D’s standard x-y coordinate system cannot handle the steep areas. Our approach is to set up a coordinate system using contours and their orthogonal trajectories. First we define the heights that will be used. In general finding good heights will require some thought and trial-and-error; in the present case we use a fourth power and stop slightly short of 1 (otherwise the differential equation solver will run into problems because of the flatness at 1).

hts = 1 - Range[1. ,  1/15, -1/15]^4

Figure  shows a contour plot with the 0.9 contour highlighted. Note that we asked for 15 contours, but got only 12; this is because the steepness at the edge defeats the resolution of ContourPlot. However, a side benefit of our flow method is that we can get an accurate contour plot even in the steep areas. We use contours = Cases[…] to store the contour data for later use.

\(\(SetOptions[ContourPlot, ContourShading -> False];\)
\)
\(\(contours = 
        Cases[Show[
            Block[{$DisplayFunction = 
                  Identity}, 
			\(ContourPlot[
                    f[x, y], {x, \(-1\), 1}, {y, \(-1\), 1}, 
                    Contours -> #[\([1]\)], 
                    ContourStyle -> #[\([2]\)], 
                    PlotPoints -> 
                      60] &\)\  /@ \ {{hts, \ {}}, \ {{0.9}, \ 
                    Thickness[
                      0.013]}}]], 
		_Line, ∞];\)\)

Code for figure

\(\(SetOptions[ContourPlot, ContourShading -> False];\)
\)
\(\(contours = 
        Cases[Show[
            Block[{$DisplayFunction = 
                  Identity}, 
			\(ContourPlot[
                    f[x, y], {x, \(-1\), 1}, {y, \(-1\), 1}, 
                    Contours -> #[\([1]\)], 
                    ContourStyle -> #[\([2]\)], 
                    PlotPoints -> 
                      60] &\)\  /@ \ {{hts, \ {}}, \ {{0.9}, \ 
                    Thickness[0.013]}}], 
			FrameTicks -> {{\(-1\), 
                  0, 1}, {\(-1\), 0, 1}, None, None}, 
            PlotRange -> {{\(-1.001\), 1}, {\(-1.001\), 
                  1}}], 
		_Line, ∞];\)\)

Now, to get the initial conditions for the differential equations that will define the flow, we need a nice set of points on the 0.9 contour: they should be somewhat evenly spaced and there should not be too many of them. Our strategy is to first remove duplicates from the points generated by ContourPlot, and then thin the set down to every third point. The package utilities EliminateDuplicates and Thin allow us to do this. We first get the contour plot; then we convert it to a standard Graphics object, use Cases to pick out the Line objects and strip off the Line header, and eliminate the duplicates and thin the list by a factor of three.

Needs["ContoursAndFlows`"]

\(\(points = 
        Thin[EliminateDuplicates[
            First[
	
              Cases[Graphics[
                  ContourPlot[f[x, y], {x, \(-1.1\), 1}, {y, \(-1.1\), 1}, 
                    Contours -> {0.9}, PlotPoints -> 25, 
                    DisplayFunction -> Identity]], 
	
                Line[x_] :> x, ∞]], 0.01], 
          3];\)\)

Show[Graphics[{PointSize[0.02], Point/@ points}]];

Code for Figure

Show[Graphics[{PointSize[0.02], Point /@  points}], 
        PlotRange -> {{-1.1,  1.1},  {-1.1,  1.1}},  
        FrameTicks -> {{-0.9, 0, 0.9},  {-0.9, 0, 0.9},  
            None,  None}];

Now the idea is to take these points as starting values and solve the differential equation that defines the flow that follows the gradient; this will yield a curve that is orthogonal to all the contours. A key point is that we use the height z as the independent variable in the differential equation. We therefore turn our attention to getting the correct differential equation.

Theoretical Digression

Assume that we have a surface S in 3-space given locally as the graph of a smooth function \(TraditionalForm
      \`f : ℝ\^2
            ⟶ℝ\), with z = f(x, y). We further assume that we have some solution curve \(TraditionalForm
      \`α : \([t\_0, t\_1]\) -> ℝ\^2\) of the initial-value problem α((Subscript[t, 0])) = ((Subscript[x, 0]), (Subscript[y, 0])), \(TraditionalForm
      \`\(α\^′\)(t) = ∇\(f(α(t))\)\). Then the curve α lifts to a curve \(TraditionalForm
      \`\(α\& ∼ \) : \([t\_0, t\_1]\) -> S\), \(TraditionalForm
      \`\(α\& ∼ \)(t) = \((x(t), y(t), f[x(t), y(t)])\)\), where α(t) = (x(t), y(t)). We wish to reparametrize this lifted curve by the z-coordinates of the lifted point.

We begin by noticing that the third coordinate of (Overscript[α,  ∼])  is the function \(TraditionalForm
      \`f∘\(α : \([t\_0, t\_1]\)\) -> 
        ℝ\), and the derivative of this map is

({{(d/d)t 
                  (f∘α) (t) = 
                (Subscript[f, x])(x(t), y(t)) (x^′)(t) + 
                  (Subscript[f, y])(x(t), y(t)) (y^′)(t)
                      }, {\(\( = 
                  ∇f·α\^′\) = 
                \(∇f·∇f = 
                  \( || ∇\(f(x(t), y(t))\)\( || \^2\)\)\)\)}})

From this it follows that f∘α either maps [(Subscript[t, 0]), (Subscript[t, 1])] to the point f((Subscript[x, 0]), (Subscript[y, 0])) (which occurs only when ((Subscript[x, 0]), (Subscript[y, 0])) is a critical point of f), or it is an orientation-preserving diffeomorphism, \(TraditionalForm
      \`f∘\(α : \([t\_0, t\_1]\)\) -> 
        \([z\_0, z\_1]\)\). In the latter case, f∘α has a smooth inverse, which we will denote by g : [(Subscript[z, 0]), (Subscript[z, 1])] -> [(Subscript[t, 0]), (Subscript[t, 1])] for convenience. We use this diffeomorphism to obtain the desired reparametrization of (Overscript[α,  ∼]) , namely \(TraditionalForm
      \`β = 
        \(α\& ∼ \)∘\(g : \([z\_0, z\_1]\)\) 
          -> S\), where

({{β(z) = 
                ( x(g(z)), y(g(z)), 
                  (f∘α) (g(z)) )}, { = 
                ( (x∘g) (z), 
                    (y∘g) (z), z ) = 
                  ( X(z), Y(z), z)}})

It remains to find a differential equation satisfied by the functions X = x∘g and Y = y∘g. But

({{(d/d)z 
                    (x∘g) (z) = 
                  (x^′)(g(z))·+
                       (g^′)(z)}, { =  
                  ((Subscript[f, x])(x(g(z)), y(g(z))))·(1
                    /((((f∘α))^′) 
                      (g(z))))}})
                        = \(\(\(f\_x\)(X(z), Y(z))\)\/\( || 
                ∇\(f(x(g(z)), y(g(z)))\)\( || \^2\)\) = \ 
            \(\(f\_x\)(X(z), Y(z))
                \)\/\(\( || ∇\(f(X(z), Y(z))\)\( || \^2\)\)\ \  . 
                \)\)

Similarly, \(TraditionalForm
      \`\(d\/dz\) \(Y(z)\) = 
        \(\(f\_y\)(X(z), Y(z))\)\/\( || 
            ∇\(f(X(z), Y(z))\)\( || \^2\)\)\).

Returning to the implementation, is is easy to set up the differential equation for X and Y and use NDSolve to find the family of curves.

(Subscript[f, x])[x_, y_] := Evaluate[(∂_x) f[x, y]]
(Subscript[f, y])[x_, y_] := Evaluate[(∂_y) f[x, y]]
∇f[a_, b_] := Evaluate[{(Subscript[f, x])[a, b], (Subscript[f, y])[a, b]}]
\(normSq[v_] := \(\(v . v\)\(
\)\)\)
{(Subscript[x, 0]), (Subscript[y, 0])} = 
        points[[1]];
(Subscript[z, 0]) = f[(Subscript[x, 0]), (Subscript[y, 0])];
{zmin, zmax} = {0, 0.99998};

Some small imaginaries show up in the interpolating functions; we use Re to sweep them away.

flow =Re[ {x[z],  y[z] }]/. First[NDSolve[
					Thread[{x'[z],y'[z],x[(Subscript[z, 0])],y[(Subscript[z, 0])]}==Join[∇f[x[z], y[z]]/normSq[∇f[x[z], y[z]]], {(Subscript[x, 0]), (Subscript[y, 0])}]],{x[z], y[z]},{z, zmin, zmax}]]

A big advantage of using z as the independent variable is that the flows can be very easily evaluated at different z-values. Figure  shows the flow emanating from the specified point in both directions, together with the points on it corresponding to the values in hts. This data will make it easy to generate very accurate contour plots, including the contours that were missing when ContourPlot was used.

\(\(ParametricPlot[
        Evaluate[flow], {z, hts[[1]], 
          hts[[\(-1\)]]}, 
	\ PlotRange -> {{\(-0.9\), 0.9}, {\(-1.1\), 
              0.85}}, 
\ \ \ Epilog -> {PointSize[0.02], 
            Point /@ points, 
            Line[Append[points, 
                points[[1]]]], 
            PointSize[0.018], 
            Point /@ Transpose[flow /. z -> hts]}];\)\)

Code for Figure

\(\(ParametricPlot[
        Evaluate[flow], {z, hts[[1]], 
          hts[[\(-1\)]]}, 
	
        PlotPoints -> 100, MaxBend -> 3, 
	
        PlotRange -> {{\(-0.9\), 0.9}, {\(-1.1\), 0.85}}, 
        FrameTicks -> {{0, 1}, {\(-1\), 0, 1}, None, 
            None}, 
\ \ \ Epilog -> {PointSize[0.02], Point /@ points, 
            Line[Append[points, 
                points[[1]]]], 
            PointSize[0.018], 
            Point /@ Transpose[flow /. z -> hts]}];\)\)

Now we put it all together by finding flows starting from each of the chosen points.

flowCurves=Re[{x[z], y[z]}]/.
		(First[NDSolve[Thread[
				{x'[z],y'[z],x[f[#]],y[f[#]]}==
							   Join[∇f[x[z], y[z]]/normSq[∇f[x[z], y[z]]], #]],
					{x[z], y[z]},{z, 0, 0.99998}]]&) /@ 
			points;

\(ParametricPlot[Evaluate[flowCurves], 
	\ \ {z, zmin, zmax}, 
        Epilog -> {
			{GrayLevel[0.5], contours}, \ 
            PointSize[0.02], Point\  /@ \ points}]\);

You can see from this figure that NDSolve had no problem obtaining a solution down to z = 0. Those points can be used to get the 0-contour, and other missing contours in figure . Indeed, we can dispense with the contour data entirely and just evaluate the flows at the values in hts to get new, and more accurate, contour data. NDSolve will have difficulty in the flat region and that is why we stopped at z = 0.99998. However, we can simply move the nearest points to the critical point to the critical point itself. It took some effort, but we have been able to automate the whole process. The following code shows how one get the surface in one step, provided the heights (hts) are provided and the options understood.

The first two arguments to ContourAndFlowPlot are obvious: the function being graphed and the variables. The third argument is a pair consisting of (Subscript[z, 0]), the z-value that will be used to generate the initial values for the flows, and the x-y plotrange for the contour plot for (Subscript[z, 0]). The fourth argument is the pair of minimum and maximum z-values and the fifth argument is the set of z-values that should be used.

Now for the options:

• ClosedCurve indicates that the contours in the plotting region are closed.
• ContourPlotPoints controls the resolution of the contour plot for (Subscript[z, 0]).
• SymbolicCriticalPoints asks that symbolic methods (Solve) be used to locate the critical points. For many examples one would have to use numerical methods.
• ShowGridBox asks that the data regarding the critical points be displayed.
• EliminateDuplicatesToleranceFactor sets the tolerance for the elimination of duplicates from the (Subscript[z, 0]) contour data
• ThinDataStep asks that the contour data be thinned to every third point.
• ShowStartPoints, when True, places points at the starting values of the flows.
• ContourStyle and FlowStyle set the style of the contour and flow lines. A natural choice is brown for contours and blue for flows (because of flowing water).
• RefineHeightsFactor asks that extra values be used to create a smoother-looking surface. The final surface will be a mix of polygons with edges and polygons without edges. This option allows us to smooth the flow lines without adding additional contours.
• MergingTolerance sets the tolerance for replacing points close to a critical point with the critical point itself.

\(\(tophalf = 
        ContourAndFlowPlot[
		f[x, y], {x, 
            y}, {0.9, {{\(-1\), 1}, {\(-1\), 1}}}, 
          hts[[{1, \(-1\)}]], 
          hts, 
		ClosedCurve -> True, 
          ContourPlotPoints -> 25, 
		SymbolicCriticalPoints -> 
            True, 
		EliminateDuplicatesToleranceFactor -> 
            0.005, 
		ThinDataStep -> 3, 
          ShowStartPoints -> False, 
		MergingTolerance -> 
            0.1, 
		ContourStyle -> Brown, 
          FlowStyle -> Blue, 
		RefineHeightsFactor -> 3, 
          ShowGridBox -> 
            True, 
		ViewPoint -> {1.5, \(-2.5\), 1.5}];\)\)

Code for Figure

\(\(tophalf = 
        ContourAndFlowPlot[
		f[x, y], {x, 
            y}, {0.9, {{\(-1\), 1}, {\(-1\), 1}}}, 
          hts[[{1, \(-1\)}]], 
          hts, 
		ClosedCurve -> True, 
          ContourPlotPoints -> 25, 
		SymbolicCriticalPoints -> 
            True, 
		EliminateDuplicatesToleranceFactor -> 
            0.005, 
		ThinDataStep -> 3, 
          ShowStartPoints -> False, 
		MergingTolerance -> 
            0.1, 
		Ticks\  -> \ {{\(- .5\), 0,  .5, 1}, {\(- .5\), 
                0, 1}, {0,  .5,  .9}}, 
		ContourStyle -> 
            GrayLevel[0.35], 
          FlowStyle -> GrayLevel[0], 
		RefineHeightsFactor -> 
            3, ShowGridBox -> 
            True, 
		ViewPoint -> {1.5, \(-2.5\), 1.5}];\)\)

The result is quite good. The top of the surface has been merged into the maximum, and the resolution at the steep areas is good.

When working on a particular example it can be useful to see only the planar view of the contours and flows. That is accomplished with the ShowPlanarViewOnly option. To illustrate the importance of the merging tolerance, we change it to a too-small value of 0.01 for the planar view in figure .

\(\(\(ContourAndFlowPlot[
        f[x, y], {x, y}, {0.9, {{\(-1\), 1}, {\(-1\), 1}}}, 
        hts[[{1, \(-1\)}]], 
        hts, 
		ClosedCurve -> True, 
        ContourPlotPoints -> 25, 
		SymbolicCriticalPoints -> 
          True, 
		EliminateDuplicatesToleranceFactor -> 
          0.005, 
		ThinDataStep -> 3, 
        MergingTolerance -> 0.05, 
		RefineHeightsFactor -> 1, 
        Axes -> False, 
\ \ \ PlanarViewOnly -> 
          True];\)\(
\)\)\)

Code for Figure

Off[ContourAndFlowPlot::"elimdp"];

\(ContourAndFlowPlot[
        f[x, y], {x, y}, {0.9, {{\(-1\), 1}, {\(-1\), 1}}}, 
        hts[[{1, \(-1\)}]], 
        hts, 
		ClosedCurve -> True, 
        ContourPlotPoints -> 25, 
		SymbolicCriticalPoints -> 
          True, 
		EliminateDuplicatesToleranceFactor -> 
          0.005, 
		ThinDataStep -> 3, 
		MergingTolerance -> 
          0.05, 
		RefineHeightsFactor -> 1, 
        ShowGridBox -> False, 
        Axes -> False, 
\ \ \ PlanarViewOnly -> True]\);

On[ContourAndFlowPlot::"elimdp"]

A consequence of our implementation is that we can get an accurate contour plot by suppressing the flow lines. That is done by using the PlanarViewOnly option as just illustrated but adding the option FlowLines -> False. We increase the contour plot points and eliminate the data thinning option for a smoother image.

\(\(\(ContourAndFlowPlot[
        f[x, y], {x, y}, {0.9, {{\(-1\), 1}, {\(-1\), 1}}}, 
        hts[[{1, \(-1\)}]], 
        hts, 
		ClosedCurve -> 
          True, 
		SymbolicCriticalPoints -> 
          True, 
		EliminateDuplicatesToleranceFactor -> 
          0.005, 
\ \ \ \ \ FlowLines -> False, 
        Axes -> False, 
\ \ \ \ PlanarViewOnly -> 
          True];\)\(
\)\)\)

Code for Figure

Off[ContourAndFlowPlot::"elimdp"];
\(\(ContourAndFlowPlot[
        f[x, y], {x, y}, {0.9, {{\(-1\), 1}, {\(-1\), 1}}}, 
        hts[[{1, \(-1\)}]], 
        hts, 
		ClosedCurve -> True, 
        ContourPlotPoints -> 50, 
		SymbolicCriticalPoints -> 
          True, 
		EliminateDuplicatesToleranceFactor -> 
          0.005, 
		ThinDataStep -> 3, 
        FlowLines -> False, 
\ \ \ \ ShowGridBox -> False, 
        Axes -> False, 
\ \ \ \ PlanarViewOnly -> True];\)\)
On[ContourAndFlowPlot::"elimdp"]

And it is easy to get the entire surface.

Show[tophalf, 
      tophalf /. {x_?NumericQ, y_, z_} :> {x, y, -z}];

Code for figure

Show[tophalf, 
        tophalf /. {x_?NumericQ, y_, z_} :> {x, y, -z}, 
        PlotRange -> {-1.1, 1.1},  
        BoxRatios -> {1, 1, 1}];

The complexity is equivalent to a 31×31 grid, which means that the comparison with Plot3D’s effort in figure  is a fair one.

(Sqrt[N])[Count[tophalf, _Polygon, ∞]]

For any single example there are often several additional steps one could take to improve the image. In the present case, one can get a closed-form polar expression for the 0.9-contour

Solve[9/10 == f[r Cos[θ], r Sin[θ]], r]
        [[4]]

Then one could use a customized function of θ to generate contour data that bunched up at the ridges. We do not pursue this here.

. The Plot Thickens

Our next example shows what must be done to show a surface with maxima, minima, and saddle points. This example is in fact handled quite well by Plot3D. But getting a contour-and-flow rendering illustrates several important points about our method, the most important of which is that one must manually decompose the domain into regions containing a single critical point. A quick contour plot shows that there is a minimum at the origin, maxima above and below it, and saddles to its left and right.

Clear[f]; 
    f[x_, y_] := ((x^2) + 3 (y^2)) (E^(1 - (x^2) - (y^2)))

ContourPlot[f[x, y], {x, -2, 2}, {y, -2.5, 2.5}, 
        PlotPoints -> 60];

The package function CriticalPointData uses numerical techniques to find and classify the critical points. It performs well on this example

CriticalPointData[
      f[x, y],  {x, -2, 2},  {y, -2.5, 2.5}]

The results are stored in SaddlePoints, Maxima, Minima, and ExtremalPoints for later use.

{SaddlePoints,  Maxima,  Minima,  ExtremalPoints}  // 
      ColumnForm

The first step is to capture the contour data for the height of the saddles. The package function SingleContour does that.

data = 
        SingleContour[f[x, y], {x, -2, 2}, {y, -2, 2}, 1, 
          PlotPoints -> 50];

Show[Graphics[Map[Point, data, {2}]]];

Now we decompose the contours at the saddles. We use the package utility DecomposeContours; if the BreakPoints option is set to Automatic (the default), then it uses the values in SaddlePoints to determine where to break the contour data. In the present case that yields six curves.

pieces = 
        DecomposeContours[f[x, y], {x, -2, 2}, {y, -2.5, 2.5}, 
          data];

Length  /@  pieces
\(\(Show[
        Graphics[
          Transpose[{
				Thickness /@ {0.01, \ 0.0001, \ 0.02, \ 
                  0.03, \ 0.02, \ 0.0001}, 
              Line /@ \ pieces}]]];\)\)

Code for Figure

\(\(\(
\)\(Show[
        Graphics[
          Transpose[{
				Thickness /@ {0.01, \ 0.0001, \ 0.02, \ 
                  0.03, \ 0.02, \ 0.0001}, Line /@ \ pieces}]], 

        AspectRatio -> 1/GoldenRatio, \ 
        FrameTicks\  -> \ {{\(-1\), 0, 1}, {\(-1\), 0, 1}, \ None, \ None}, \ 
        PlotRange -> {{\(-1.5\), \ 1.5}, \ {\(-2.2\), 2.2}}];\)\)\)

We now go through and define various data sets. First we get the border of a region containing the upper maximum; the data is improved by deleting the points near the saddles and replacing them with approximations such as {-0.99,  0.01}; when using such approximations one must be careful to move away from the saddle so that the flow will go in the proper direction. And then symmetry gives us the other region with a maximum. The exterior region, regionext, is obtained by pasting together the right pieces in the right order, again using a saddle approximation.

\(\(regionmaxpos = 
        EliminateDuplicates[
          Join[
			Rest[
              Drop[Thin[
                  Rest[pieces[[4]]], 
                  3], \(-1\)]], 
			{{\(-0.99\), 0.01}}, 
			Thin[
              Drop[Drop[
                  Reverse[pieces[[1]]], 
                  3], \(-3\)], 3], 
			{{0.99,  .01}}], \ 
          0.05];\)
\)
\(\(regionmaxneg = 
        regionmaxpos /. 
		\ \ {x_?NumericQ, 
              y_} :> {x, \(-y\)};\)
\)
\(\(regionext = 
        Join[Take[regionmaxpos, 
            26], 
		\ \ {{\(-1.01\), \ 0.01}}, 
		\ \ Reverse[
            Take[regionmaxneg, 26]], 
		\ {{1.01, 0.01}}];\)\)

Next we get the region containing the central minimum, using symmetry and approximations to the saddles.

\(\(regionmin = \((quad2 = 
            Sort[Prepend[
                Select[
						regionmaxpos, #[[2]] < 
                        1/2 && \(-1\) ≤ #[[1]] < 0 &], 
					\ \ \ \ {\(-0.99\), 
                  0.01}]]; 
		\ \ Join[quad2, 
            Reverse[quad2] /. {x_, y_} :> {\(-x\), 
                  y}])\);\)
\)
\(\(regionmin\  = 
        EliminateDuplicates[
          Join[
			regionmin, \ 
            Reverse[regionmin]\  /. \ 
				\ \ {x_, \ 
                  y_}\  :> \ {x, \(-y\)}], \ 0.05];\)\)

A look at these regions shows that we have a good approximation to the regions we want.

\(\(Show[
        Graphics[{
	GrayLevel[0.8], 
            Polygon[regionmaxneg], 
\ \ GrayLevel[0.3], 
            Polygon[regionmaxpos], 
\ \ GrayLevel[0], 
            Polygon[regionmin], 
\ \ PointSize[0.02], 
            Point /@ regionmaxpos, \ Point\  /@ \ regionmaxneg, 
	
            GrayLevel[1], 	Point\  /@ \ ExtremalPoints}]];\)\)

Code for Figure

\(\(Show[
        Graphics[{
			{GrayLevel[0.8], 
              Polygon[regionmaxneg]}, 
			{GrayLevel[0.3], 
              Polygon[regionmaxpos]}, 
			{GrayLevel[0], 
              Polygon[regionmin]}, 
			PointSize[0.02], 
            Point /@ regionmaxpos, \ 
            Point\  /@ \ regionmaxneg, 
		GrayLevel[1], 	
            Point\  /@ \ ExtremalPoints}], Frame -> True\ , \ 
        AspectRatio -> 1/GoldenRatio, \ 
        FrameTicks\  -> \ {{\(-1\), 0, 1}, \ {\(-1\), 0, 1}, None, \ None}, \ 
        PlotRange -> {{\(-1.4\), \ 1.4}, \ {\(-2.1\), \ 2.1}}];\)\)

The idea now is to run the flow down from the exterior points and down from the boundary of the minimum and up from the points surrounding the two maxima. Recall from section 2 that the third argument to FlowLineSurface can be a contour and plot range, which ContourPlot will use. But it can also be a set of contour data that the user has obtained manually, as we have done in this example.

\(\(SetOptions[ContourAndFlowPlot, 
	
        EliminateContourDuplicates -> False, 
	
        Saddles -> SaddlePoints, 
	Extrema -> ExtremalPoints, 
	
        ContourStyle -> GrayLevel[0.35], 
	
        FlowStyle -> GrayLevel[0], 
	
        DisplayFunction -> Identity, 
	ClosedCurve -> True, \ 
	
        MergingTolerance\  -> \ 0.07];\)\)

hts  =  Join[Range[0. , 3, 1/ 9],  {2.99, 0.005}];

\(\(\(Block[{$DisplayFunction = Identity}, 
	\ \ exterior = 
          ContourAndFlowPlot[f[x, y], {x, y}, regionext, {0.5, 1.0}, 
            hts\ ]; 
	\ \ maxneg = 
          ContourAndFlowPlot[f[x, y], {x, y}, regionmaxneg, {1, 2.99}, 
            hts]; 
	\ \ maxpos = 
          ContourAndFlowPlot[f[x, y], {x, y}, regionmaxpos, {1, 2.99}, 
            hts]; 
	\ \ min = 
          ContourAndFlowPlot[f[x, y], {x, y}, regionmin, {0.005, 1.01}, 
            hts]];\)\(
\)\(	\)\)\)

We now show all parts, enhancing the plot by the use of LiftToViewPoint (described in my 8:2 column) to add the 1-contour as a thick line moved a little bit in the direction of the viewpoint.

\(pr = \ {{\(-2.2\), 2.2}, {\(-2.6\), 2.6}, {\(-0.2\), 3.2}}; 

    br = \ {1, 1, 0.5}; 
vp = {1.3, \(- .3\), 0.5}; 


    Show[min, maxneg, maxpos, exterior, 
	
      Graphics3D[{GrayLevel[0.2], AbsoluteThickness[2], 
			
          Line[Append[regionext, First[regionext]] /. 

              						\ {x_, y_}\  :> 
                LiftToViewPoint[{x, y, f[x, y]}, vp, br, pr, 0.05]], 
			
          Line[Append[regionmin, First[regionmin]] /. 

              						\ {x_, y_}\  :> 
                LiftToViewPoint[{x, y, f[x, y]}, vp, br, pr, 0.05]]}], 
\ \ 
      BoxRatios -> br, ViewPoint -> vp, 
	PlotRange -> pr, \ 
      DisplayFunction -> $DisplayFunction]; \)

Code for Figure

pr =  {{-2.2, 2.2}, {-2.6, 2.6}, {-0.2, 3.2}};
br =  {1, 1, 0.5};
vp = {1.3, - .3, 0.5};

Show[min,maxneg,maxpos,exterior,
	\(Graphics3D[{GrayLevel[0.2], 
              AbsoluteThickness[2], 
			Line[
                Append[regionext, 
                    First[regionext]] /. 
						\ {x_, 
                      y_}\  :> 
                    LiftToViewPoint[{x, y, f[x, y]}, vp, br, pr, 
                      0.05]], 
			Line[
                Append[regionmin, 
                    First[regionmin]] /. 
						\ {x_, 
                      y_}\  :> 
                    LiftToViewPoint[{x, y, f[x, y]}, vp, br, pr, 0.05]]}]\),
  BoxRatios -> br,ViewPoint -> vp,
	Axes -> None,
	AxesLabel -> {"x", "y", "z"},
	AxesEdge  ->  {{-1, -1}, {1, -1}, {1, -1}},
	PlotRange -> pr, DisplayFunction  ->  $DisplayFunction];

. A Singular Example

Here is a function studied in detail in [Hern et al 1998], where it appears as an example related to the failure of Fubini’s theorem on multiple integration. Its graph is quite attractive and Plot3D has no big problems. Still, there are inaccuracies near the origin, and a much more attractive plot can be obtained by using the contour-and-flow method. We reload the package to reset the options to their defaults.

Get["ContoursAndFlows`"];

\(

\(\(SetOptions[#, 
	
            PlotRange -> {{\(-1.2\), 1.2}, {\(-1.2\), 1.2}, {\(-9.05\), 
                  9.05}}, 
	BoxRatios -> {1, 1, 1.2}, 
            Boxed -> False, 
	
            ViewPoint\  -> {0.5, \ \(-2\), \ 1}\ ] &\)\  /@ \ {Plot3D, 
          Graphics3D, ContourAndFlowPlot};\)
\)
\(\(SetOptions[ContourAndFlowPlot, 
        DisplayFunction -> Identity];\)
\)
f[x_, y_] := (((x^2) - (y^2))/((((x^2) + (y^2)))^2))

Plot3D[f[x, y], {x, -1.2, 1.2}, {y, -1.2, 1.2}, 
      PlotPoints -> 60, ClipFill -> None];

Code for figure

Plot3D[f[x, y], {x, -1.2, 1.2}, {y, -1.2, 1.2}, 
        PlotPoints -> 60, ClipFill -> None, 
        Ticks -> {{-1, 0, 1}, {0,  .5}, {-6, 0, 9}}, 
        AxesEdge -> {{-1, -1}, {-1, -1}, {-1, -1}}];

To use the contour-and-flow method we will exploit the symmetry of the surface to the maximum extent possible. Thus the crux of the matter will be obtaining some contour data and getting an image in the lower half of the second quadrant. We could do this by using a closed form for the z = (1/2) contour, but we will leave that avenue to the reader. Or we could just ask ContourAndFlowPlot to get the contour data as in earlier examples. But we will get more control by using ContourPlot directly to generate the data. We want the contour data to include a point close to ((Sqrt[2]), 0); adding this exact point would lead to problems in NDSolve because the flow from that point would follow the ridge exactly and never reach 0. We get around the problem by adding the nearby flow-start point (-(Sqrt[2]), -0.001). Then, by including ((Sqrt[2]), 0) in the singularity list, we guarantee that the added initial point is merged with (-(Sqrt[2]), 0) at the final rendering. Of course, the origin is added to the singularity list so that points near the z-axis (i.e., above the origin) are moved to the z-axis. We also force the data point nearest the origin to be (-0.02,  0.02).

\(\(data = 
        Join[{{\(-0.02\), 0.02}}, \ 
		\ \ \ \ Rest[
            Drop[Thin[
                EliminateDuplicates\ [
				First[
                    SingleContour[
                      f[x, y], 
						{x, \(-1.5\), 0}, {y, 0.001, 
                        0.6}, 1/2, PlotPoints -> 50]], 0.01], 
                6], \ \(-1\)]], \ {{\(-\@2\), \  .001}}];\)\)

Show[Graphics[{PointSize[0.02],  Point/@data}]];

Code for Figure

\(\(Show[Graphics[{PointSize[0.02], \ Point /@ data}], 
	
        PlotRange -> {{\(-1.5\), \  .1}, \ {\(-0.1\), \ 0.6}}, 
	
        FrameTicks -> {{\(-\@2\), \(-1\), 0}, {0, 0.5}, \ None, \ 
            None}];\)\)

Here are some heights appropriate for this surface.

hts = 
      Sort[Select[
          Range[0. , (9^(1/3)), ((2  (9^(1/3)))/40)]^3, # == 
                0  ∨  Abs[#] ≥ 0.12  &]]

Now we look at what ContourAndFlowPlot can do with the data. We set Saddles and Extrema to the empty set to avoid the Automatic setting, which would entail a search for critical points. The PlanarViewOnly option gives us the 2-dimensional view.

\(\(eighthOfSurface = 
        DeleteCases[
		ContourAndFlowPlot[f[x, y], {x, y}, data, {0, 10}, 
            hts, 
		Extrema -> {}, Saddles -> {}, 
            Singularities -> {{0, 0}, \ {\(-\@2\), 
                  0}}, \ 
				ShowStartPoints\  -> \ True, 
            RefineHeightsFactor -> 2, 
		PlanarViewOnly -> True, 
            Axes -> 
              None, 
		PlotRange -> {{\(-1.5\), \ 0.1}, {\(-0.1\), 
                  1.2}}, 
	\ MergingTolerance -> 
              0.08], _Point, ∞];\)\)

Code for Figure

\(\(eighthOfSurface = 
        DeleteCases[
		ContourAndFlowPlot[f[x, y], {x, y}, data, {0, 10}, 
            hts, 
		Extrema -> {}, Saddles -> {}, 
            Singularities -> {{0, 0}, \ {\(-\@2\), 
                  0}}, \ 
				ShowStartPoints\  -> \ True, 
            RefineHeightsFactor -> 2, 
		PlanarViewOnly -> 
              True, 
			Axes -> None, \ 
            FrameTicks -> {{\(-1.4\), \ \(-1\), \ \(- .5\), \ 
                  0}, {0, \  .5, \ 1}, None, \ 
                None}, 
		PlotRange -> {{\(-1.5\), \ 0.1}, {\(-0.1\), 
                  1.2}}, 
	\ MergingTolerance -> 
              0.08], _Point, ∞];\)\)

This looks good. Note how the origin is properly included, thanks to the slightly larger merging tolerance and the inclusion of the origin as a singularity. When this image is reflected to get the full surface, the number of polygons will be near 3600, as we can discover using the package’s CountPolygons function. This means that the PlotPoints setting of 60 in figure  yields a fair comparison.

8 CountPolygons[
        eighthOfSurface, {-1.2, 1.2}, {-1.2, 1.2}]

Now we can see the entire surface, which fits together perfectly around the z-axis.

\(\(quarterSurface = {eighthOfSurface, \ 
		\ eighthOfSurface\  /. \ {x_?NumericQ, y_, z_} :> {x, \(-y\), z}};\)\)

\(\(Show[quarterSurface, \ 
        quarterSurface /. \ {
			{{x_?NumericQ, y_, 
                  z_} :> {\(-x\), y, 
                  z}}, 
		\ \ {{x_?NumericQ, y_, z_} :> {y, 
                  x, \(-z\)}}, 
			{{x_?NumericQ, y_, 
                  z_} :> {y, \(-x\), \(-z\)}}}, 
	
        DisplayFunction\  -> \ $DisplayFunction\ ];\)\)

. One More Example

Each surface presents its own problems. In this section, we present only the final outcome of a surface-generation project, without the code (it is in the electronic version) showing how that outcome was obtained. The reader interested in using the our package should pick a function of his or her own and try to get a contour-and-flow image of it.

The example is z = (((x^2)(1 + y))^3) + (y^2), a surface that S.W. has been pondering for ten years. This polynomial, discovered by Calvert and Vamanamurthy [CV] is of interest because it has a local minimum at (0, 0) that is not an absolute minimum, yet it has no additional critical points. This sort of thing cannot happen for functions y = f(x). However, getting a good image of the surface takes a bit of work to get all the parameters right, and it must be done in pieces because of the topology of the surface. The contour-and-flow plot shown in figure  shows the minimum clearly and the rollover out of the local minimum is easy to comprehend thanks to the flow lines. Note how the flow lines descending from the right part of the image either drop off the ridge steeply down to negative infinity or drop off on the other side into a river that flows into the origin.

Anyone who generates many surfaces might well run into the following lighting issue: a polygon can show up in a color that seems totally inappropriate, given the lighting of its neighbors. This can happen because a quadrilateral in 3-space will typically not be planar, and so not have an obvious normal vector to use to define the lighting. Mathematica splits such quadrilaterals into two triangles and lights each triangle according to a normal direction that is defined by an averaging process, and so is the same for the two triangles. This can cause problems for some awkward quadrilaterals for which the averaging process is inappropriate. The fix is to turn all the quadrilaterals into triangles first. The following routine does this, and was used on the example of figure  to eliminate some bad coloring effects.

\(\(\(QuadsToTriangles[
        g_]\)\(\ \)\(:=\)\(\ \)\({Graphics3D[EdgeForm[]], 
	
        g\  /. \ Polygon[{p1_, \ p2_, \ p3_, \ 
                p4_}]\  :> 
	{Polygon[{p1, \ p2, \ p3}], \ 
              Polygon[{\ p3, \ p4, \ p1}]}}\)\(
\)\)\)

Code for Figure   (include in e-version)

Defining f, setting of options, and defining the heights

Get["ContoursAndFlows`"];

Clear[f];    f[x_, y_] := (x^2) (((1 + y))^3) + 7 (y^2)
vp = {-1.1, -0.4, 0.5};
br  =  {1, 1,  .4};
pr  =  {{-10, 10}, {-4, 3.5}, {-30, 81}};
\(\(SetOptions[ContourAndFlowPlot, ShowStartPoints -> True, 
	
        PlotRange -> pr, 
	
        StartPointStyle -> AbsolutePointSize[3], 
	Axes -> Automatic, 
        EliminateContourDuplicates -> False, 
	
        AxesLabel -> {"x", "y", "z"}, 
        BoxRatios -> br, 
	ViewPoint -> vp, 
        Boxed -> False, \ DisplayFunction -> Identity, 
	
        ContourStyle -> \((GrayLevel[0.35]; \ Brown)\), 
        FlowStyle -> \((GrayLevel[0]; Blue)\)];\)
\)
\(\(surfopts\  = \ 
        Sequence[	
          PlotRange -> {{\(-10\), 10}, {\(-4\), 3.5}, {\(-30\), 
                81}}, 
	Axes -> Automatic, 
          AxesLabel -> {"x", "y", "z"}, 
          BoxRatios -> {1, 1, 0.4}, 
	
          ViewPoint -> {\(-1.1\), \(-0.4\), 0.5}, 
          Boxed -> 
            False, 
		\ DisplayFunction -> $DisplayFunction];\)
\)
\(\(SetOptions[FindCriticalPoints, 
        SymbolicCriticalPoints -> True];\)
\)
\(\(SetOptions[Graphics3D, 
        FilterOptions[Graphics3D, \ 
          Options[ContourAndFlowPlot]]];\)
\)
\(\(zvals = {\(-30\), \(-22\), \(-15\), \(-7\), 0, 0.005, 0.055, 0.439, 
          1.48, 3.5, 5, 6, 7, 
		10, 12, 15, 21, 28, 37, 47, 58, 69, 
          80};\)\)

A special set of points of the form (x,  -1) with f(x,-1)=7

f[x, -1]

\(data7manual = \(({#, \(-1\)} &)\) /@ 
		Sort[
          Join[Range[\(-50\), \(-10\), 5], Range[10, 50, 5], 
            Range[\(-9\), 9, 1.2]]]\)

Note that some of the start points are OUTSIDE the plot range; they are needed to start some flows to give good coverage. The z-value of 6 is deleted from the heights for this patch. Now we generate the surface emanating from the straight line just defined.

\(\(surffrom7 = 
        ContourAndFlowPlot[f[x, y], {x, y}, 
		data7manual, {0.005, 80}, 
          Complement[zvals, {6}], LiftingFactor -> 0.02, 
          RefineHeightsFactor -> 3];\)\)

Show[surffrom7,  DisplayFunction -> $DisplayFunction,  
        Boxed -> True];

Generating the left half (positive y) part of the figure), from z = 80 down to z = 0

Here I fiddle with the data to get a nice red border flow at x  =  10 and x  =  -10.

yspec[x_,  c_,  y0_]  :=   y  /.  
        FindRoot[f[x, y]  ==  c,  {y,  y0}]

\(\(data\  = \ 
        SingleContour[
		f[x, y], \ {x, \(-18\), \ 18}, \ {\ 
            y, \ \(-1\), \ 3.5}, 80, \ ContourPlotPoints -> 100];\)\)
datanew = Thin[data[[1]], 12];

datanew[[{2, 3, 4, 
          5, -5, -4, -3, -2, -1}]]

datanew  =  
        ReplacePart[
          datanew,  {-9.94,  
            yspec[-9.94,  80,  datanew[[2, 2]]]},  2];
datanew  =  
        ReplacePart[
          datanew,  {10.06,  
            yspec[10.06,  80,  
              datanew[[-2, 2]]]},  -2];

\(\(surf80to0 = 
        ContourAndFlowPlot[f[x, y], {x, y}, datanew, 
			{0.005, 80}, 
          Complement[zvals, {6, 7, 10, 15}], 
\ LiftingFactor -> 0.005, 	
          RefineHeightsFactor -> 2];\)\)

Show[surf80to0,  DisplayFunction -> $DisplayFunction,  
        Boxed -> True];

Generating the part going to z = -∞

\(\(data\  = \ 
        SingleContour[
		f[x, y], \ {x, \(-18\), \ 0}, \ {\ 
            y, \ \(-4\), \ \(-1\)}, 7, \ 
          ContourPlotPoints -> 60];\)\)
datanew = 
        Thin[data[[1]], 
          6];

datanew[[{1, 2, 3, 4, 5, 6}]]

datanew  =  
        ReplacePart[
          datanew,  {-9.94,  
            yspec[-9.94,  7,  datanew[[6, 2]]]},  6];

\(\(surf7down1 = \ 
        ContourAndFlowPlot[f[x, y], {x, y}, 
          datanew, 
		\ \ \ {\(-40\), 80}\ , 
          Select[zvals, \ # ≤ 0\  || # ≥ 7\  &], 
          RefineHeightsFactor -> 2];\)\)

Show[surf7down1, 
        DisplayFunction -> $DisplayFunction];

And now we symmetrize to get the reflection.

surf7down2 = 
        surf7down1 /. {x_?NumericQ, y_?NumericQ, z_} :> {-x, 
              y, z};

A subtle point arises in mathematica's method of lighting of surfaces. A quadrilateral in a surface might not lie in a plane and so Mathematica divides the quadrilateral into two triangles for rendering, lighting each triangles according to a normal vector defined from the quadrilateral. This can be considered as an averaging process and generally improves the lighting. But there are examples, generally where the quadrilateral is seriously onnplanar, when this can cause a polygon to be lit in an inappropriate away. The user can spot this by, say, a black polygon in a region of gray polygons. One can fix this by turning all the quadrilaterals in a surface into pairs of triangles; this will cause a normal vector appropriate to each triangle to be used. The following function does this in a way that destroys the borders of the polygons, since that is the sort of image that ContourAndFlowPlot will produce (recall that the various curves, in order to be colored, are added to the surface atop the borderless polygons)...xxx...check to see if this is discussed earlier....perhaps do so where CFPlot is discussed.

QuadsToTriangles[g_] := {Graphics3D[EdgeForm[]], 
        g /. Polygon[{p1_, p2_, p3_, 
                p4_}] :> {Polygon[{p1, p2, p3}], 
              Polygon[{p3, p4, p1}]}}

A planar view to see how the curves fit together

PlanarView[surf7down1,  surf7down2, surffrom7, surf80to0,  
        PlotRange -> {{-10, 10},  {-4, 3.5}}];

Putting it all together

We delete the black points that were shown in the partial surfaces as an aid, we turn the quadrilaterals into triangles to improve the lighting (necessary for one of the parts only) and show everything together. The rendering of this is a little slow; make sure your memory is adequate! Applying QuadsToTriangles to all four pieces yields a slightly more pleasing coloring scheme, so we do that here.

\(\(\( (*\ if\ low\ on\ memory, \ 
      this\ is\ better\ *) \)\(
\)\(
\)\(Show[
        DeleteCases[{surf7down1, surf7down2, \ QuadsToTriangles[\ surf80to0], 
            surffrom7}, \ _Point, \ ∞], Axes -> None, \ 
        surfopts];\)\)\)

Show[
        QuadsToTriangles[
          DeleteCases[{surf7down1, surf7down2,   surf80to0, 
              surffrom7},  _Point,  ∞]], Axes -> None,  
        surfopts];

References

[Calvert and M. K. Vamanamurthy 1980]  B. Calvert and M. K. Vamanamurthy, Local and global extrema for functions of several variables, Journal of the Australian Mathematical Society (Series A) 29 (1980) 362–368.

[Goetz and Wagon, 1996] Ron Goetz and Stan Wagon, Adaptive surface plotting: a beginning, Mathematica in Education and Research 5:3 (1996) 74–83.

[Hern et al 1998] Thomas Hern, Cliff Long, and Andy Long, Looking at order of integration and a minimal surface, College Math. Journal 29 (1998) 128–133.

About the Authors

Ron Goetz is an enthusiastic proponent of Mathematica as a visualization tool in education and an alumnus of Wagon’s Mathematica course in the Colorodo mountains. He spends his spare time working on a dissertation in differential topology at Washington University in St. Louis. 

Ron Goetz
Department of Mathematics
St. Louis Community College at Meramec
11333 Big Bend Blvd, 
St. Louis, MO  63122
rgoetz@mcmail.stlcc.cc.mo.us

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition, and is especially appreciative of the new ways Mathematica allows us to look at mathematical objects. He is author or coauthor of several books, including Animating Calculus, VisualDSolve, Which Way Did the Bicycle Go?, and the just-released and much expanded second edition of Mathematica in Action. He teaches a high-altitude course, Rocky Mountain Mathematica, every summer in the mountains of Colorado.

Stan Wagon
Department of Mathematics and Computer Science
Macalester College
St. Paul, MN 55105
wagon@macalester.edu

Code for Electronic Version (* ed: possibly add code from the article too?*)

The ContourAndFlow package

(* remove FlowLinesOnly and ContoursOnly -- now commented out *)
(* add message that lifting factor requires specific PlotRange done I think*)

(* Add LIghting Fix from Rootcellar file!!! *)


(* All tolerances, etc MUST be made relative!! *)

Off[General::spell1]
BeginPackage["ContoursAndFlows`", 
   {"Graphics`Colors`", "Utilities`FilterOptions`"}]; 

(* BottomlessTetrahedron::usage = "BottomlessTetrahedron[pt,size] gives the three upright faces of a tetrahedron centered at pt. Its color is solid red and can be used to mark points on a surface."; *)

PlanarViewOnly::usage = "PlanarViewOnly is an option to ContourAndFlowPlot that causes the graphic shown to be the planar view of the surface, even though the output is still the surface object.";

ContourPlotPoints::usage = "ContourPlotPoints is an option to SingleContour and ContourAndFlowPlot and sets the PlotPoints in the call to ContourPlot used to generate the contour data.";

FlowStep::usage = "";
refine::usage = "";
HighResolutionFlows::usage = "";
FlowResolution::usage = "";

LiftToViewPoint::usage = "LiftToViewPoint[p, vp, br, pr, lf] moves the point p an amount lf toward the viewpoint vp, where the box ratios and plot range are given by br and vp, respectively.";

ShowFlowSurface::usage = "";
ContoursToShow::usage = "";
FlowLines::usage = "";
Polygons::usage = "";

LiftingFactor::usage = "";
FontBackground::usage = "";
BreakPoints::usage = "";

SymbolicCriticalPoints::usage = "";

ThinDataStep::usage = "ThinDataStep is an option to ContourAndFlowPlot that, when set to n, causes the contour data to be thinned to every nth point.";

RefineHeightsFactor::usage = "RefineHeightsFactor is an option to ContourAndFlowPlot that refines the heights for the purposes of the meshless polygons only";

ShowPoints::usage = "ShowPoints is an option to PlanarView that, when False, causes any points in the 3-dimensional object to be deleted, as opposed to projected into the plane.";

EliminateDuplicatesToleranceFactor::usage = "EliminateDuplicatesToleranceFactor gives a multipler that decides how close is close enough to be eliminated as a duplicate. The factor multiplies xmax - xmin where these are the x-bounds on the contour data.";

CountPolygons::usage = "CountPolygons[g, {a, b}, {c, d}] counts the number of polygons in g that have any of their x and y coordinates within the interval {a,b}, and {c, d}, resp.";

MergingTolerance::usage = "";
Singularities::usage = "";
FlowStyle::usage = "";
PlanarView::usage = "";
EliminateContourDuplicates::usage = "";

ShowStartPoints::usage = "ShowStartPoints is an option to ContourAndFlowPlot that adds the starting points from the contour data to the surface.";

StartPointStyle::usage = "StartPointStyle is an option to ContourAndFlowPlot that sets the style of the starting points if they are shown."; 

FlowLinesOnly::usage = "";
ContoursOnly::usage = "";

Thin::usage = "Thin[a, n] takes every nth element of the list a, starting with the first.";

Up::usage = "";
Down::usage = "";

ShowGridBox::usage = "ShowGridBox is an option to CriticalPointData that causes a gridbox display of the critical points and their type to be displayed.";

ConvergingFlowDirection::usage = "";

EliminateDuplicates::usage = "EliminateDuplicates[a,t] removes duplicates from the set a of points in the plane, using the tolerance t.";

SaveSolutionsName::usage = "";
ContourLabels::usage = "";

Maxima::usage = "Maxima contains the local maxima of the function under consideration once CriticalPointData has been run.";

Minima::usage = "Minima contains the local minima of the function under consideration once CriticalPointData has been run.";

ExtremalPoints::usage = "ExtremalPoints contains the local extrema of the function under consideration once CriticalPointData has been run.";

SaddlePoints::usage = "SaddlePoints contains the saddle points of the function under consideration once CriticalPointData has been run.";

Failures::usage = "";

LabeledContourPlot::usage = "";

SingleContour::usage = ""; 

ContourAndFlowPlot::usage = "ContourAndFlowPlot[f, {x, a,b},{y,c,d}, contourdata, {zmin, zmax}, zvalues] generates a surface starting from the points in contourdata (which are assumed to lie on a contour) and following the flows between z-values zmin and zmax. a,b,c,d are the range for the critical point search. If (x,y) values are found that are near an extremum or saddle, then such points are replaced by the actual extremum or saddle. If contourdata has the form {z,{x0,x1},{y0,y1}} then an attempt is made to find the contour data appropriate for height z, using xi and yi as the plotrange."; 

ImproveContours::usage = "ImproveContours[f, {x,y}, data, z, wp] takes data that comes close to being a contour set for f at height z and improve the points by solving a differential equation based on the gradient. The wp argument is working precision, which defaults to machine precision.";

ImproveContourPoints::usage = "";

normsq::usage = "normsq[p] gives p.p for a point p in the plane."; 

VerticalPlotPoints::usage = "I think this is obsolete. Delete from file where it occurs";

ClosedCurve::usage = "ClosedCurve is an option to ContourAndFlowPlot that indicates that the contour data is to be viewed as a closed curve.";

DecomposeContours::usage = "DecomposeContours[f, {x, a, b}, {y, c, d}, data]...";

CriticalPointType::usage = "";

Saddles::usage = "Saddles is an option to ContourAndFlowPlot.";
Extrema::usage = "Extrema is an option to ContourAndFlowPlot";

CriticalPointData::usage = "CriticalPointData[f, {x, a, b},
   {y, c, d}] returns the minima, maxima, saddles, and extreme points (min and max). It stores the results in Minima, Maxima, SaddlePoints, and ExtremePoints. It also displays a grid box showing the critica points.";

GradientFlow::usage = "GradientFlow[{f, g}, {zmin, zmax}, 
  points,";

FixingTolerance::usage = "FixingTolerance is an option for FixingRules...";
SeparationOfExtrema::usage = "SeparationOfExtrema is an option for FixingRules...";


FindCriticalPoints::usage="FindCriticalPoints[funcs,{x,xmin,xmax},{y,ymin,ymax}] finds all nontangential solutions to {f=0, g=0}.";

lift::usage = "";

CriticalPointData::symfl="The preceding messages indicate that a symbolic attempt to find the critical points failed. A numerical attempt is now being made.";

ContourAndFlowPlot::elimdp = "Points that are close have been eliminated from the contour data. `` points were discarded.";

ContourAndFlowPlot::elimdp1 = "Points that are close have been eliminated from the contour data. One point was discarded.";

ContourAndFlowPlot::bdcntr = "The contour data has `` branches. You will have to generate the data set you want manually (use SingleContour) and enter that data as an argument to ContourAndFlowPlot.";

ContourAndFlowPlot::badpr = "When the LiftingFactor is nonzero, a full plot range (three pairs) must be provided.";


GradientFlow::fail="NDSolve ran into problems at starting value ``.";

ShowFlowSurface::badbox="Automatic is not a valid setting for the BoxRatios option when used with ShowFlowSurface when LiftingFactor is nonzero.";

Begin["`Private`"]; 

Clear[SingleContour, fixrule, CriticalPointData, CriticalPointDataRaw];

		
Options[SingleContour] = {ContourPlotPoints -> 50, 
  ImproveContourPoints -> False}; 

SingleContour[f_, {x_, a_, b_}, {y_, c_, d_}, z_?NumericQ, opts___] := 
SingleContour[f, {x, a, b}, {y, c, d}, z, opts] = Module[
 {iQ,pp,out},
  {iQ, pp} = {ImproveContourPoints, ContourPlotPoints} /. {opts} /. Options[SingleContour];
Off[ContourGraphics::"ctpnt"];

 out = First /@ Cases[
   Graphics[ContourPlot[f, {x, a, b}, {y, c, d}, 
   PlotPoints -> pp, 
   Evaluate[FilterOptions[ContourPlot, opts]], 
   ContourShading -> False, ContourSmoothing -> False, 
   DisplayFunction -> Identity, Contours -> {z}]],
      _Line, Infinity];
 On[ContourGraphics::"ctpnt"]; 
 If[iQ, out /. (l_List)?(VectorQ[#1[[1]]] & ) :> 
    ImproveContours[f, {x, y}, l, z], 
  out]]
     
normsq = Compile[{{x, _Real, 1}}, x . x]; 
 


Options[GradientFlow] = {SaveSolutionsName->False};

GradientFlow[{f_, g_}, {zmin_, zmax_}, points_, opts___] := 
   Module[{x, y, z, sol},
   
   saveQ = SaveSolutionsName /. {opts} /. Options[GradientFlow]; 
   
   If[StringQ[saveQ], Clear[Evaluate[StringJoin["Global`",saveQ]]]];
  Catch[ Flatten[sol = (  
   Check[NDSolve[
     {(*x'[z] == f[x[z], y[z]]/(f[x[z], y[z]]^2 + g[x[z], y[z]]^2), 
   y'[z] == g[x[z], y[z]]/(f[x[z], y[z]]^2 + g[x[z], y[z]]^2),*)
      
   x'[z] == f[x[z], y[z]]/(temp=
       f[x[z], y[z]]^2 + g[x[z], y[z]]^2), 
   y'[z] == g[x[z], y[z]]/temp, 
   x[FFF @@ #] == First[#], y[FFF @@ #] == Last[#]}, 
     {x, y}, {z, zmin, zmax}, 
       FilterOptions[NDSolve, opts]],
       Message[GradientFlow::fail, #]; Throw[hf]])& /@ 
        points, 1];
 
If[saveQ =!= False, 
Set[Evaluate[ToExpression[StringJoin["Global`",ToString[saveQ]]]],
   Flatten[{x, y} /. sol, 1]]];
    Flatten[{x, y} /. sol, 1]]]; 
    
Options[ContourAndFlowPlot] = Join[{
 FlowLines->True,
 FlowStyle -> {},
 ContourLines->True,
 ContourStyle->{},
 FlowStep->1,
 HighResolutionFlows->False,
 ContoursToShow -> Automatic,
 Axes->True,
 EliminateDuplicatesToleranceFactor -> 0.01, 
 VerticalPlotPoints -> 20, ClosedCurve -> False, 
 Saddles->Automatic, SeparationOfExtrema -> 0.5,
 MergingTolerance -> 0.025,
 Extrema->Automatic, 
 Singularities -> {},
 BoxRatios->{1,1,.4}, 
 ShowStartPoints->False,
 PlanarViewOnly->False,
 FlowResolution -> 30,
 EliminateContourDuplicates->True,
 ThinDataStep->1,
 RefineHeightsFactor->1,
 LiftingFactor->0,
 StartPointStyle->{PointSize[0.022]},
 BoxRatios->{1,1,.4}}, Options[Graphics3D]]; 


ContourAndFlowPlot[expr_, {x_, y_}, contourdata_, {zmin_, zmax_}, 
  zvalues_List:{}, opts___Rule] := 
 Module[{fx, fy, n,ss}, (* FFF and hf NOT moduled!! *)
{clQ, zDivs, sad, ext, tol, condir, ptQ, ptsty,
edupQ,singularities, mergetol, eduptol,thinstep,pvq, rhfac, csty, fsty, conQ, 
flQ,lf,br, hrfq, flowppts, fstep} = 
 {ClosedCurve, VerticalPlotPoints, Saddles, Extrema,
 SeparationOfExtrema, ConvergingFlowDirection,
 ShowStartPoints, StartPointStyle, EliminateContourDuplicates,
 Singularities, MergingTolerance, EliminateDuplicatesToleranceFactor,
 ThinDataStep, PlanarViewOnly, RefineHeightsFactor, ContourStyle, FlowStyle,
 ContourLines, FlowLines, LiftingFactor,
 BoxRatios,
 HighResolutionFlows, FlowResolution, FlowStep} /. {opts}
     /. Options[ContourAndFlowPlot]; 


zvalues1 = refine[zvactual = Sort[Join[
   Cases[zvalues, _?NumericQ],
   temprules = (Map[Thread[{x,y}->#] &, 
      (zpts=DeleteCases[zvalues, _?NumericQ])]);
    If[temprules != {},expr /.temprules, {} ]]], rhfac];
    
    
contourdata1 =If[Union[Length /@ contourdata] == {2}, N[contourdata],
 (ss = SingleContour[expr,
      N@{x,contourdata[[2,1,1]], contourdata[[2,1,2]]},
	  N@{y,contourdata[[2,2,1]], contourdata[[2,2,2]]},
	  N[contourdata][[1]],
	  FilterOptions[SingleContour, opts]])[[1]]];



hf = HoldForm[ContourAndFlowPlot[expr, {x, y}, 
    "<< >>", {zmin, zmax}, zvalues, opts]];
    
vp = ViewPoint /. {opts} /. Options[Graphics3D];
pr = PlotRange /. {opts} /. Options[Graphics3D];
br = BoxRatios /. {opts} /. Options[ContourAndFlowPlot];

If[!MatchQ[pr, {{_?NumericQ,_},{_,_},{_,_}}]  && lf != 0, 
   Message[ContourAndFlowPlot::badpr]; Return[hf]];
    
If[ListQ[ss] && Length[ss] != 1, Message[ContourAndFlowPlot::bdcntr,
Length[ss]]; Return[hf]];

If[!ListQ[ptsty], ptsty = {ptsty}];
If[!ListQ[fsty], fsty = {fsty}];
If[!ListQ[csty], csty = {csty}];
flQ = flQ || hrfq;
If[hrfQ, lf = Max[0.01, lf]];

{{a,b}, {c,d}} = DataBounds[contourdata1];

cdata = If[edupQ, 
EliminateDuplicates[contourdata1, (b-a) * eduptol], contourdata1]; 
 
ll = Length[contourdata1]-Length[cdata];
If[ll>1, Message[ContourAndFlowPlot::elimdp,ll]]; 
If[ll==1, Message[ContourAndFlowPlot::elimdp1]];

cdata = Thin[cdata, thinstep];
n = Length[cdata];
If[ext===Automatic || sad === Automatic,
   {Minima, Maxima, SaddlePoints, ExtremalPoints, Failures} = 
  CriticalPointData[expr, {x,a,b}, {y,c,d},
   FilterOptions[FindCriticalPoints, opts],
   FilterOptions[CriticalPointData, opts],
   ShowGridBox->False]];
sad = sad /. Automatic -> SaddlePoints;
ext = ext /. Automatic -> ExtremalPoints;

  FFF = Evaluate[expr /. {x -> #1, y -> #2}] & ; 
  fx = Evaluate[Simplify[D[FFF[x, y], x]] /. {x -> #1, y -> #2} ] &; 
  fy = Evaluate[Simplify[D[FFF[x, y], y]] /. {x -> #1, y -> #2}] &;
  solutions = GradientFlow[{fx, fy}, {zmin, zmax}, cdata, 
      FilterOptions[NDSolve, opts],
      FilterOptions[GradientFlow,opts]];
      
If[Head[solutions]===HoldForm, Return[solutions]]; 
 Do[{xsol[i], ysol[i]} = solutions[[i]], {i, n}]; 

fd = If[hrfq, 
 Cases[ParametricPlot[Evaluate[Table[{
			xsol[i][z], ysol[i][z]}, {i, 1, Length[solutions], fstep}]],
	{z, zmin,zmax}, DisplayFunction->Identity,
	PlotPoints -> flowppts, Evaluate[FilterOptions[ParametricPlot, opts]]], 
	  _Line, Infinity] ,{}] /. {xx_?NumericQ, yy_} :> {xx, yy, expr /. {x->xx, y->yy}};

zVals = Union[Sort[Select[N[zvalues1],
  zmin <= # <= zmax & ] /. 
  {} -> Range[N[zmin], zmax, (zmax - zmin)/(zDivs - 1.)]]]; 
  
(* the following seems to require a lifting factor....so not sure...*)

(*fdd = Table[Line[{xsol[i][#], ysol[i][#]}& /@ zVals], 
   {i, 1, Length[solutions], fstep}]/. {xx_?NumericQ, yy_} :> {xx, yy, FFF[xx,yy]};
*)

pointData = Outer[{xsol[#2][#1], ysol[#2][#1], #1} &, zVals, 
      Range[Length[cdata]]]; 

If[clQ, pointData = (Append[#1, First[#1]] & ) /@ pointData]; 
pointDatanew = (Partition[#1, 2, 1] & ) /@ pointData; 
ef = If[{csty, fsty, conQ, flQ} == {{},{}, True, True},
     {},EdgeForm[]]; 

liftrule = If[lf == 0, {}, {xx_?NumericQ, yy_, zz_} :> 
   LiftToViewPoint[{xx, yy, zz}, vp, br,pr,lf]];


gr = Chop[{
  polys = MapThread[Polygon[Flatten[{#1, Reverse[#2]}, 1]]&, 
 {Drop[pointDatanew, -1], Drop[RotateLeft[pointDatanew], -1]}, 2];

( {DeleteCases[
  If[conQ === True, Join[csty, Union[Map[Sort,Flatten[polys /. 
    Polygon[{aa_, bb_, cc_, dd_}] :> 
       {Line[{aa, bb}], Line[{cc, dd}]},1],{2}]]], {}], 
    Line[{aa_, bb_}] /; !MemberQ[N[zvactual], aa[[-1]]] ,
         Infinity], 
  If[flQ === True,   
     Join[fsty, If[hrfq, fd, Union[Map[Sort,Flatten[polys /. 
        Polygon[{aa_, bb_, cc_, dd_}] :> 
     {Line[{bb, cc}], Line[{aa, dd}]},1],{2}]]]], {}]
         }) /. liftrule, ef, polys}, 10^-5];

ends1 = Chop[Table[{xsol[i][zmin], ysol[i][zmin]}, {i,n}], 10^-6];
ends2 = Chop[Table[{xsol[i][zmax], ysol[i][zmax]}, {i,n}], 10^-6];
ends = If[ringmeshmax[ends1] < ringmeshmax[ends2], ends1, ends2];


fixrule = Join[({xx_?NumericQ, yy_?NumericQ, _?NumericQ} /; (
  (normsq[{xx, yy} - #]) <= mergetol^2) :>
     (Append[#, expr /. {x->#[[1]], y->#[[2]]}])) & /@ sad,
  
   (({xx_?NumericQ, yy_?NumericQ, zz_?NumericQ} /; 
    normsq[{xx, yy} - #] <= mergetol^2) :> 
     Append[#, zz])& /@ singularities, 
  If[ext == {}, {},
    {{xx_?NumericQ, yy_?NumericQ, _?NumericQ} /;
        (MemberQ[ends, {xx,yy}] &&
     (First[Sort[(normsq[{xx, yy} - #1] & ) /@ ext]] < mergetol^2)) :> (mm = Sort[ext, normsq[#1 -{xx,yy}] < normsq[#2 - {xx, yy}]&][[1]];
    AppendTo[mm, expr /. {x->mm[[1]], y->mm[[2]]}])} ]]; 


If[pvq,
  PlanarView[Graphics3D[{Chop[gr, 10^-6],
   If[ptQ, 
     Join[ptsty, Point /@ (lift[Join[cdata, zpts], FFF,
        (zmax - zmin)/100])],{}]}] /. fixrule, 
FilterOptions[PlanarView, opts],
Sequence @@ Options[PlanarView],
FilterOptions[Graphics, opts],
AxesLabel->None,
FilterOptions[Graphics, Options[ContourAndFlowPlot]]];
  df = Identity, df = $DisplayFunction];

Block[{$DisplayFunction = df},
Show[Graphics3D[{Chop[gr, 10^-6],
If[ptQ, Join[ptsty, Point /@ (lift[Join[cdata, zpts], FFF,
        (zmax - zmin)/100])],{}]}] /. (fixrule),
FilterOptions[Graphics3D, Sequence@@ 
 DeleteCases[{opts},
PlotRange->{{_,_},{_,_}}]], 
FilterOptions[Graphics3D, Sequence @@ Options[ContourAndFlowPlot]]]]]; 


pos[d_, p_] := Position[d, 
   First[Select[d, 
      normsq[# - p] == Min[(normsq[# - p] & ) /@ d]&, 1]]][[1,1]];

Options[DecomposeContours] = {BreakPoints->Automatic};

DecomposeContours[expr_, {x_, a_, b_}, {y_, c_, d_}, data_?MatrixQ, opts___Rule] := DecomposeContours[expr, {x, a, b}, {y, c, d}, {data}, opts]

DecomposeContours[expr_, {x_, a_, b_}, {y_, c_, d_}, data_, opts___Rule] := Module[{},   
  sadd = BreakPoints /.{opts} /. Options[DecomposeContours];
  sadd = sadd /. Automatic :> 
  CriticalPointData[expr, {x,a,b}, {y,c,d}, ShowGridBox->False,
   FilterOptions[CriticalPointData, opts],
   FilterOptions[FindCriticalPoints, opts],
   FilterOptions[FindRoot, opts]][[3]];
If[Depth[sadd]==2, sadd = {sadd}];
data1 = If[Depth[data]==3, {data}, data];
new = (0&) /@ data1; 
Do[len = Length[data1[[i]]]; 

  ppp = Sort[(pos[data1[[i]], #] &) /@ sadd]; 
  new[[i]] = 
    (Take[data1[[i]], #] &) /@ (zzz=
( Cases[ Partition[Append[ppp, First[ppp] + len], 2, 1]/.
  {xx_, yy_}  /; yy>len :>{ {xx,len},{1, yy-len-1}}, 
  {_Integer, _Integer}, Infinity]);
  zzz), 
     {i, Length[data1]}]; 
  Flatten[new, 1]]; 
     
CriticalPointType[f_, {x_, y_}, {a_, b_}] := Module[{dxx, disc}, 
  dxx = D[f, x, x] /. {x -> a, y -> b}; 
  disc = dxx D[f, y, y] - D[f, x, y]^2 /. {x -> a, y -> b}; 
   Which[
     dxx < 0 && disc > 0, "Maximum",
     dxx > 0 && disc > 0, "Minimum",
     disc < 0, "SaddlePoint", 
     disc == 0, NCriticalPointType[f, {x,y}, {a,b}]]]
     
NCriticalPointType[f_, {x_,y_}, cp_] := Module[{rad, d,s},
 	rad = N[0.05];
    d=Table[f /. Thread[{x,y} -> cp + rad {Cos[t], Sin[t]}], 
      {t, 0, 2 Pi, 2 Pi /15}];
    s=Union[Sign[	d - f /. Thread[{x,y} -> cp]]];
        Which[s=={-1}, "Maximum", 
    s == {1}, "Minimum",
    s == {-1,1}, "SaddlePoint",
    True, "Failure"]	]
      
lift[{a_?NumericQ, b_}, f_, d_:0] := {a, b, f[a, b]+d};

lift[a_List, f_, d_:0] := (lift[#, f, d] & ) /@ a

Options[FindCriticalPoints] = {PlotPoints -> 59,
SymbolicCriticalPoints->False}; 

FindCriticalPoints[funcs_, {x_, xmin_, xmax_}, 
   {y_, ymin_, ymax_}, opts___] := Module[
 {seeds, tt, 
 fy = Compile[{x, y}, Evaluate[funcs[[2]]]], 
 contourData, symsol = $Failed},

symtry = SymbolicCriticalPoints /. {opts} /. Options[FindCriticalPoints];
If[symtry, symsol = Check[Solve[{funcs[[1]]==0, funcs[[2]]==0}, {x,y}], $Failed]];

If[symsol =!= $Failed,
Return[  Select[Chop[N[Union[{x,y} /. symsol]]], 
    Im[#[[1]]] == 0 && Im[#[[2]]] == 0 &&xmin<= #[[1]]<= xmax &&
   xmin <= #[[2]] <= xmax &]],
 If[symtry,Message[CriticalPointData::symfl]]];

contourData = (
  Off[ContourGraphics::"ctpnt"];
  First /@ Cases[Graphics[ContourPlot[
  funcs[[1]], {x, xmin, xmax}, {y, ymin, ymax},
     Evaluate[FilterOptions[ContourPlot, opts]], 
    Contours -> {0}, ContourShading -> False, ContourSmoothing->False,PlotPoints -> 59, 
DisplayFunction -> Identity]], _Line, Infinity]);

 
On[ContourGraphics::"ctpnt"];
seeds = Flatten[(#1[[1 + Flatten[Position[Rest[tt = 
  Sign[Apply[fy, #1, 2]]]*Rest[RotateRight[tt]], -1]]]]&)/@
    contourData, 1];
 If[seeds == {}, {}, Chop[Select[Union[({x, y} /. 
   FindRoot[{funcs[[1]] == 0, funcs[[2]] == 0}, {x, #1[[1]]}, 
   {y, #1[[2]]}, Evaluate[FilterOptions[FindRoot, opts]]] & ) /@ 
   seeds], xmin < #1[[1]] < xmax && ymin < #1[[2]] < ymax &], 10^(-5)]  ]]


Options[CriticalPointData] = {ShowGridBox->True, SymbolicCriticalPoints->False};

CriticalPointDataRaw[f_, {x_, a_, b_}, {y_, c_, d_}, 
 opts___] := 
(CriticalPointDataRaw[f, {x, a, b}, {y, c, d},___] = (
 CriticalPoints = FindCriticalPoints[
  {D[f, x], D[f, y]}, 
   {x, a-Random[Real, {.9, 1}]10^-1 (b-a), b}, 
   {y, c-Random[Real, {.9, 1}]10^-1 (b-a), d},
   FilterOptions[FindCriticalPoints, opts],
   FilterOptions[FindRoot, opts]]; 
  
Minima = Select[CriticalPoints,
  CriticalPointType[f, {x, y}, #] == "Minimum" &]; 
Maxima = Select[CriticalPoints,
  CriticalPointType[f,{x, y},#] == "Maximum" &]; 
SaddlePoints = Union[Select[CriticalPoints,
  CriticalPointType[f,{x,y},#] == "SaddlePoint" &]];
ExtremalPoints = Join[Maxima,Minima];
Failures = Select[CriticalPoints,
   CriticalPointType[f, {x, y}, #] == "Failure" &];
  
{Minima, Maxima, SaddlePoints, ExtremalPoints, Failures}))   
   
CriticalPointData[f_, {x_, a_, b_}, 
   {y_, c_, d_}, opts___] := (
{Minima, Maxima, SaddlePoints, ExtremalPoints, Failures} = CriticalPointDataRaw[
   f, {x, a, b}, {y, c, d},      
     FilterOptions[FindCriticalPoints, opts],
      FilterOptions[FindRoot, opts]];
 
 grQ= ShowGridBox/.{opts}/.Options[CriticalPointData];
 
 If[grQ, Print, Identity][ DisplayForm[StyleBox[
  GridBox[Map[MakeBoxes[#]&, 
   Prepend[Transpose[{
    {"Minima", "Maxima", "SaddlePoints", 
       "ExtremalPoints", "Failures"}, 
   col2={Minima, Maxima, SaddlePoints, 
      ExtremalPoints, Failures} ,
   col2 /. {xx_?NumericQ,yy_} :> (f /. 
         Thread[{x,y}->{xx,yy}])}],
   {"Type", "Points", "Function Values"}] /.
 {} -> "", {2}],
    ColumnAlignments -> Left,
    GridFrame -> 1, 
    RowLines -> {2,1},
    ColumnLines -> 1],
    FontFamily -> Times, 
    Background -> GrayLevel[0.85]]]];
    
 {Minima, Maxima, SaddlePoints, ExtremalPoints, Failures})


ImproveContours[f_, {x_, y_}, data_List, ht_?NumericQ, wp_Integer:$MachinePrecision, 
   opts___Rule] := Module[
   {FF, fx, fy, z, data1, ht1}, 
   data1 = (SetPrecision[#1, wp] & ) /@ data; ht1 = SetPrecision[ht, wp]; 
    FF[xx_, yy_] := f /. {x -> xx, y -> yy}; 
    fx[xx_, yy_] := Evaluate[D[f, x] /. {x -> xx, y -> yy}]; 
    fy[xx_, yy_] := Evaluate[D[f, y] /. {x -> xx, y -> yy}]; 
    ((zz = FF @@ #1; {x[z], y[z]} /. 
           First[NDSolve[{x'[z] == 
               fx[x[z], y[z]]/(fx[x[z], y[z]]^2 + fy[x[z], y[z]]^2), 
              y'[z] == fy[x[z], y[z]]/(fx[x[z], y[z]]^2 + fy[x[z], y[z]]^2), 
              x[zz] == #1[[1]], y[zz] == #1[[2]]}, {x[z], y[z]}, {z, zz, ht1}, 
             FilterOptions[NDSolve, opts], WorkingPrecision -> wp]] /. z -> ht) & ) /@ data1]

nearest[l_, v_] := First[Sort[l, Abs[#1 - v] < Abs[#2 - v] &]]
 
median[l_] := l[[Round[Length[l]/2]]]

Options[LabeledContourPlot] = {Contours -> Range[-10, 10, 2], ContourLabels -> {}, FontColor->RGBColor[1,1,.6], FontBackground->RGBColor[0,0,1], 
FontWeight -> "Plain", FontSize->10}; 

LabeledContourPlot[f_, {x_, a_, b_}, {y_, c_, d_}, opts___] := 
   (zv = ContourLabels /. {opts} /. Options[LabeledContourPlot]; 

cp = ContourPlot[f, {x, a, b}, {y, c, d}, Evaluate[FilterOptions[ContourPlot, opts]], DisplayFunction -> Identity]; 

{zvals, fc, tbg,fs,fw} = {Contours, FontColor, FontBackground,FontSize, FontWeight} /. {opts} /. Options[LabeledContourPlot]; ld = Cases[Graphics[cp], _Line, Infinity]; 

Show[cp, Graphics[{((xy = median[#1[[1]]]; 
   near = nearest[zvals, f /. Thread[{x, y} -> xy]]; 
     If[MemberQ[Join[zv, N[zv]], near], 
         Text[near, xy, Background -> tbg, 
         TextStyle -> 
           {FontSize -> fs, FontWeight -> fw, 
           FontColor -> fc}], {}]) & ) /@ ld}], 
      DisplayFunction -> $DisplayFunction, 
      FilterOptions[Graphics, opts]])
      
EliminateDuplicates[l_, tol_] := Module[{ll},
 ll = l //.( {a___, b_?VectorQ, c_, d___}/; normsq[b-c]<tol^2 ):> {a,b,d};
 If[normsq[l[[1]]-l[[-1]]] < tol^2, Rest[ll],ll]]
  
Thin[l_, s_] := l[[Range[1, Length[l], s]]]

ringmeshmax[l_] := Max[Map[
  normsq[#[[1]]-#[[2]]]&,Partition[l,2,1]]]
     
Options[PlanarView] = {FlowStyle -> {RGBColor[1,0,0]},
ContourStyle->{}, Frame->True, ShowPoints->True,
ContourLines -> True,FlowLines -> True,
DisplayFunction->$DisplayFunction};

PlanarView[x_Graphics, opts___] := Show[x, opts]

PlanarView[surf_List, opts___] := Module[{a,b,c,d},
{fsty,csty, ptsQ, conQ, flowQ} =
 {FlowStyle, ContourStyle, ShowPoints, ContourLines, FlowLines} /. {opts} /. Options[PlanarView];
 
If[!ListQ[fsty], fsty={fsty}];
If[!ListQ[csty], csty={csty}];
 
linerule = Polygon[{a_, b_, c_, d_}] :> 
  Evaluate[If[Count[surf, _Line, Infinity] > 10, 
  {}, 
  {If[conQ, Join[csty, {Line[{a, b}], Line[{c, d}]}],{}],
   If[flowQ, Join[fsty, {Line[{b, c}], Line[{a, d}]}],{}]}]];
 
Show[Graphics[
 DeleteCases[surf /. Graphics3D[xx__] :> {xx} /. 
   {s1_?NumericQ, s2_, _} :> {s1, s2} /. linerule,
    _Rule | _RuleDelayed | If[!ptsQ, _Point], Infinity]], 
  Sequence@@({FilterOptions[Graphics, opts, AxesLabel->None,
      Axes->None,
        Sequence @@ Options[PlanarView]]} /. 
    (PlotRange->{xr_,yr_, z_}) :> (PlotRange->{xr,yr}))]] /;
       Complement[Head /@ {opts}, {Rule, RuleDelayed}]=={} 
    
PlanarView[surf__, opts___] := (
PlanarView[
Join[{surf} /. xx_Graphics3D :> {xx}], opts]) /; Complement[
  Head /@ {surf}, {Graphics3D, List}]=={} && Complement[Head /@ {opts}, {Rule, RuleDelayed}]=={}

    
CountPolygons[g_, {a_, b_}, {c_, d_}] := 
  Length[Cases[g, Polygon[{p_, q_, r_, s_}] /; 
     a <= p[[1]] <= b && c <= p[[2]] <= d || 
      a <= q[[1]] <= b && c <= q[[2]] <= d || 
      a <= r[[1]] <= b && c <= r[[2]] <= d || 
      a <= s[[1]] <= b && c <= s[[2]] <= d, Infinity]]

bd[x_] := Sign[x] 2. Ceiling[Abs[x]] /; x != 0
bd[0] = 1;
 
bdl[0] = -1;
bdl[x_] := Sign[x] (Floor[Abs[x]]-1)/2. /; x != 0

DataBounds[cd_] := (xmin = Min[First /@ cd];
xmax=Max[First /@ cd];
ymin=Min[Last /@ cd];
ymax=Max[Last /@ cd];
{Which[xmin<=0<=xmax, {bd[xmin], bd[xmax]},
    0 <= xmin, {bdl[xmin], bd[xmax]},
    xmax <= 0, {bd[xmin], bdl[xmax]}],
Which[ymin<=0<=ymax, {bd[ymin], bd[ymax]},
    0 <= ymin, {bdl[ymin], bd[ymax]},
    ymax <= 0, {bd[ymin], bdl[ymax]}]})
   
   
LiftToViewPoint[pt_, vp_, br_, pr_, lf_] := (
 vec = (Apply[Plus, pr, {1}]/2 + 
    (-vp*Apply[Subtract, pr, {1}])Max[br]/br - pt);
   pt + lf* vec/Sqrt[normsq[vec]])  


Options[ShowFlowSurface] = 
   {ContourStyle -> {RGBColor[0,0,1]}, 
   FlowStyle -> {RGBColor[1,0,0]}, 
   Polygons->True,
   Contours -> True, FlowLines -> True, LiftingFactor -> 0};

ShowFlowSurface[gg__, opts___Rule] := Module[{hf=HoldForm[ShowFlowSurface[gg,opts]]},
    ggg = Flatten[{gg}]; 
  {cQ, fQ, fsty, csty,lf,polQ} = 
   {Contours, FlowLines, FlowStyle, ContourStyle, 
   LiftingFactor, Polygons} /. 
         {opts} /. Options[ShowFlowSurface]; 
pr = PlotRange /. {opts};
vp = ViewPoint /. {opts} /. FullOptions[ggg[[1]]]; 
br = BoxRatios /. {opts} /. FullOptions[ggg[[1]]]; 
prs = Map[Flatten, PlotRange /@ 
         Cases[ggg, _Graphics3D, Infinity], 1]; 
      xmin = Min[First /@ prs]; 
      xmax = Max[(#1[[2]] & ) /@ prs]; 
      ymin = Min[(#1[[3]] & ) /@ prs]; 
      ymax = Max[(#1[[4]] & ) /@ prs]; 
      zmin = Min[(#1[[5]] & ) /@ prs];
      zmax = Max[Last /@ prs];

If[pr === All && lf != 0, 
   Message[ShowFlowSurface::badpr]; Return[hf]];

If[pr === PlotRange,
  pr = {{xmin,xmax},{ymin,ymax},{zmin,zmax}},
  If[Length[pr]==3,      
  {xmin, xmax, ymin, ymax, zmin, zmax} = Flatten[pr],
        
  pr ={{xmin,xmax},{ymin,ymax},pr}] ];

If[br === Automatic && lf != 0, 
   Message[ShowFlowSurface::badbox]; Return[hf]];



liftrule = If[lf==0, {}, {x_?NumericQ, y_, z_} :> 
   LiftToViewPoint[{x, y, z}, vp, br,pr,lf]];
 
If[!ListQ[fsty], fsty = {fsty}]; 
If[ !ListQ[csty], csty = {csty}]; 
allpolys = Cases[ggg, Polygon[{a_, b_, c_, d_}] /; 
    (xmin <= a[[1]] <= xmax && ymin <= a[[2]] <= ymax) || 
    (xmin <= b[[1]] <= xmax && ymin <= b[[2]] <= ymax) || 
    (xmin <= c[[1]] <= xmax && ymin <= c[[2]] <= ymax) || 
    (xmin <= d[[1]] <= xmax && ymin <= d[[2]] <= ymax), 
        Infinity];
Show[Graphics3D[{
      {If[cQ === True, 
            Join[csty, Union[Map[Sort,Flatten[allpolys /. 
              Polygon[{a_, b_, c_, d_}] :> 
       {Line[{a, b}], Line[{c, d}]},1],{2}]]], {}], 
           If[fQ === True, 
            Join[fsty, Union[Map[Sort,Flatten[allpolys /. 
              Polygon[{a_, b_, c_, d_}] :> 
               {Line[{b, c}], Line[{a, d}]},1],{2}]]], {}]
            } /. 
          liftrule,
          {EdgeForm[{}], If[polQ===True, allpolys, {}]}
}], Evaluate[FilterOptions[Graphics3D, opts]], 
       Evaluate[Sequence @@ 
         (#1 -> (#1 /. FullOptions[ggg[[1]]]) & ) /@ 
          {Boxed, BoxRatios, PlotRange, Lighting, 
           ViewPoint}]]]/;
       Complement[Head /@ {gg}, {List, Graphics3D}] == {}
       
  
 refine[l_List, n_Integer] := 
  Union[N[Append[Flatten[Partition[l, 2, 1] /. 
      {a_?NumericQ, b_} :> Range[a, b, (b - a)/n], 1], 
    Last[l]]]]
    
End[]; EndPackage[];

Code for Figure   (include in e-supplement)

Defining f, setting of options, and defining the heights

Get["ContoursAndFlows`"];

Clear[f];    f[x_, y_] := (x^2) (((1 + y))^3) + 7 (y^2)
vp = {-1.1, -0.4, 0.5};
br  =  {1, 1,  .4};
pr  =  {{-10, 10}, {-4, 3.5}, {-30, 81}};
\(\(SetOptions[ContourAndFlowPlot, ShowStartPoints -> True, 
	
        PlotRange -> pr, 
	
        StartPointStyle -> AbsolutePointSize[3], 
	Axes -> Automatic, 
        EliminateContourDuplicates -> False, 
	
        AxesLabel -> {"x", "y", "z"}, 
        BoxRatios -> br, 
	ViewPoint -> vp, 
        Boxed -> False, \ DisplayFunction -> Identity, 
	
        ContourStyle -> \((GrayLevel[0.35]; \ Brown)\), 
        FlowStyle -> \((GrayLevel[0]; Blue)\)];\)
\)
\(\(surfopts\  = \ 
        Sequence[	
          PlotRange -> {{\(-10\), 10}, {\(-4\), 3.5}, {\(-30\), 
                81}}, 
	Axes -> Automatic, 
          AxesLabel -> {"x", "y", "z"}, 
          BoxRatios -> {1, 1, 0.4}, 
	
          ViewPoint -> {\(-1.1\), \(-0.4\), 0.5}, 
          Boxed -> 
            False, 
		\ DisplayFunction -> $DisplayFunction];\)
\)
\(\(SetOptions[FindCriticalPoints, 
        SymbolicCriticalPoints -> True];\)
\)
\(\(SetOptions[Graphics3D, 
        FilterOptions[Graphics3D, \ 
          Options[ContourAndFlowPlot]]];\)
\)
\(\(zvals = {\(-30\), \(-22\), \(-15\), \(-7\), 0, 0.005, 0.055, 0.439, 
          1.48, 3.5, 5, 6, 7, 
		10, 12, 15, 21, 28, 37, 47, 58, 69, 
          80};\)\)

A special set of points of the form (x,  -1) with f(x,-1)=7

f[x, -1]

\(data7manual = \(({#, \(-1\)} &)\) /@ 
		Sort[
          Join[Range[\(-50\), \(-10\), 5], Range[10, 50, 5], 
            Range[\(-9\), 9, 1.2]]]\)

Note that some of the start points are OUTSIDE the plot range; they are needed to start some flows to give good coverage. The z-value of 6 is deleted from the heights for this patch. Now we generate the surface emanating from the straight line just defined.

\(\(surffrom7 = 
        ContourAndFlowPlot[f[x, y], {x, y}, 
		data7manual, {0.005, 80}, 
          Complement[zvals, {6}], LiftingFactor -> 0.02, 
          RefineHeightsFactor -> 3];\)\)

Show[surffrom7,  DisplayFunction -> $DisplayFunction,  
        Boxed -> True];

Generating the left half (positive y) part of the figure), from z = 80 down to z = 0

Here I fiddle with the data to get a nice red border flow at x  =  10 and x  =  -10.

yspec[x_,  c_,  y0_]  :=   y  /.  
        FindRoot[f[x, y]  ==  c,  {y,  y0}]

\(\(data\  = \ 
        SingleContour[
		f[x, y], \ {x, \(-18\), \ 18}, \ {\ 
            y, \ \(-1\), \ 3.5}, 80, \ ContourPlotPoints -> 100];\)\)
datanew = Thin[data[[1]], 12];

datanew[[{2, 3, 4, 
          5, -5, -4, -3, -2, -1}]]

datanew  =  
        ReplacePart[
          datanew,  {-9.94,  
            yspec[-9.94,  80,  datanew[[2, 2]]]},  2];
datanew  =  
        ReplacePart[
          datanew,  {10.06,  
            yspec[10.06,  80,  
              datanew[[-2, 2]]]},  -2];

\(\(surf80to0 = 
        ContourAndFlowPlot[f[x, y], {x, y}, datanew, 
			{0.005, 80}, 
          Complement[zvals, {6, 7, 10, 15}], 
\ LiftingFactor -> 0.005, 	
          RefineHeightsFactor -> 2];\)\)

Show[surf80to0,  DisplayFunction -> $DisplayFunction,  
        Boxed -> True];

Generating the part going to z = -∞

\(\(data\  = \ 
        SingleContour[
		f[x, y], \ {x, \(-18\), \ 0}, \ {\ 
            y, \ \(-4\), \ \(-1\)}, 7, \ 
          ContourPlotPoints -> 60];\)\)
datanew = 
        Thin[data[[1]], 
          6];

datanew[[{1, 2, 3, 4, 5, 6}]]

datanew  =  
        ReplacePart[
          datanew,  {-9.94,  
            yspec[-9.94,  7,  datanew[[6, 2]]]},  6];

\(\(surf7down1 = \ 
        ContourAndFlowPlot[f[x, y], {x, y}, 
          datanew, 
		\ \ \ {\(-40\), 80}\ , 
          Select[zvals, \ # ≤ 0\  || # ≥ 7\  &], 
          RefineHeightsFactor -> 2];\)\)

Show[surf7down1, 
        DisplayFunction -> $DisplayFunction];

And now we symmetrize to get the reflection.

surf7down2 = 
        surf7down1 /. {x_?NumericQ, y_?NumericQ, z_} :> {-x, 
              y, z};

A subtle point arises in mathematica's method of lighting of surfaces. A quadrilateral in a surface might not lie in a plane and so Mathematica divides the quadrilateral into two triangles for rendering, lighting each triangles according to a normal vector defined from the quadrilateral. This can be considered as an averaging process and generally improves the lighting. But there are examples, generally where the quadrilateral is seriously onnplanar, when this can cause a polygon to be lit in an inappropriate away. The user can spot this by, say, a black polygon in a region of gray polygons. One can fix this by turning all the quadrilaterals in a surface into pairs of triangles; this will cause a normal vector appropriate to each triangle to be used. The following function does this in a way that destroys the borders of the polygons, since that is the sort of image that ContourAndFlowPlot will produce (recall that the various curves, in order to be colored, are added to the surface atop the borderless polygons).

QuadsToTriangles[g_] := {Graphics3D[EdgeForm[]], 
        g /. Polygon[{p1_, p2_, p3_, 
                p4_}] :> {Polygon[{p1, p2, p3}], 
              Polygon[{p3, p4, p1}]}}

A planar view to see how the curves fit together

PlanarView[surf7down1,  surf7down2, surffrom7, surf80to0,  
        PlotRange -> {{-10, 10},  {-4, 3.5}}];

Putting it all together

We delete the black points that were shown in the partial surfaces as an aid, we turn the quadrilaterals into triangles to improve the lighting (necessary for one of the parts only) and show everything together. The rendering of this is a little slow; make sure your memory is adequate! Applying QuadsToTriangles to all four pieces yields a slightly more pleasing coloring scheme, so we do that here.

\(\(\( (*\ if\ low\ on\ memory, \ 
      this\ is\ better\ *) \)\(
\)\(
\)\(Show[
        DeleteCases[{surf7down1, surf7down2, \ QuadsToTriangles[\ surf80to0], 
            surffrom7}, \ _Point, \ ∞], Axes -> None, \ 
        surfopts];\)\)\)

Show[
        QuadsToTriangles[
          DeleteCases[{surf7down1, surf7down2,   surf80to0, 
              surffrom7},  _Point,  ∞], Axes -> None,  
          surfopts]];