Contour Initialization

Needs["Graphics`Colors`"];

colorQ[a_]:=(Head[a]===RGBColor);
grayQ[a_]:=(Head[a]===GrayLevel);
thickQ[a_]:=(Head[a]===Thickness);
absthickQ[a_]:=(Head[a]===AbsoluteThickness);
notruleQ[u_] := !(Head[u]===Rule);
ruleQ[u_] := (Head[u]===Rule);

Norm::usage = "Norm[expr] returns the Euclidean length of the
		vector, 'expr'.";

Norm[a_?VectorQ] := Sqrt[a.a];

Angle[a_?VectorQ,b_?VectorQ] := ArcCos[Simplify[(a.b)/(Norm[a] Norm[b])]];

Proj[a_?VectorQ,b_?VectorQ] := ((a.b)/(b.b))b;

tol = .5*10^-6;

OrthUp[a_?VectorQ] :=
	Module[{b,c,d,f},
		d = N[a/Norm[a]];
		If[Chop[d[[3]]^2 - 1,tol] == 0, b={1,0,0} ,
		b = {0,0,1} - d[[3]] d; b = b/Norm[b];];
	Return[b]];
	
OrthSet[a_?VectorQ] :=
	Module[{b,c,d},
		d = N[a/Norm[a]];
		If[Chop[d[[1]],tol]==0,b={1,0,0},
			 If[Chop[d[[2]],tol]==0,b={0,1,0},
			  b={-d[[2]],d[[1]],0}]];
			b=b/Norm[b];c = Cross[d,b];
	Return[{b,c/Norm[c],d}]];

				
Rotate[a_?VectorQ,theta_] :=
	Module[{u,v,r},
		u = OrthSet[a];
		v = Transpose[u];
		r = N[{{Cos[theta],-Sin[theta],0},{Sin[theta],Cos[theta],0},{0,0,1}}];
	Return[v.r.u]];

Vector2D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics[Point[{0,0}]]];
	,						
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = 0.075 (ntip-ntail) ;
	up = 0.15 {-Punit[[2]],Punit[[1]]};
	arrow = Line[{ntail,ntip,ntip-Punit+up,ntip,ntip-Punit-up}];
	]];
Return[Graphics[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

Vector3D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,u2,u3,u4,rotmat,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics3D[Point[{0,0,0}]]];					
	,		
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics3D[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = (ntip-ntail);
	rotmat = Rotate[Punit,N[Pi/2]];
	up = OrthUp[Punit];
	Punit = .075 Punit;
	up = .15 Norm[Punit] up;
	u2 = rotmat.up;
	u3 = rotmat.u2;
	u4 = rotmat.u3;
	arrow = Line[{ntail,ntip,
				ntip - Punit + up,ntip,
				ntip - Punit + u2,ntip,
				ntip - Punit + u3,ntip,
				ntip - Punit + u4}];
	]];
Return[Graphics3D[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

aRCtAN[y_,x_] := Module[{v},
	If[ Chop[N[x],tol] == 0 ,If[ N[y] < 0, v = -N[Pi/2], v = N[Pi/2]],
		v = N[ArcTan[y/x]];
		If[ x < 0 && y < 0, v = N[Pi -  v],v]];
	Return[v]];

Arc2D[X_,Y_,opts___] := Graphics[{opts,
		Disk[{0,0},1,{Min[
							aRCtAN[
						N[X[[2]]],N[X[[1]]]
								]
							,
							aRCtAN[
							N[Y[[2]]],N[Y[[1]]]
								]
						]
						,
					Max[
							aRCtAN[
						N[X[[2]]],N[X[[1]]]
								]
							,
							aRCtAN[
							N[Y[[2]]],N[Y[[1]]]
								]
						]}
				]}
				];
Arc3D[X_,Y_] := 
	Module[{t,U,V},
		U[V_]:= N[(1/Sqrt[V.V]) V];
		Return[Graphics3D[
		Polygon[Append[Table[U[t X + (1-t) Y],{t,0,1,0.1}],{0,0,0}]]
				]
		]
	];

(*Zero Arrow causes trouble. Put in Vector?*)

Vector[tail_,tip_,opts___]:=
	If[Length[N[tail]]==2,Return[Vector2D[tail,tip,opts]],
		Return[Vector3D[tail,tip,opts]]];
		
Arrow[vec_,opts___]:= 
	Module[{a,b},
				If[Length[N[vec]]==2,
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0}};
				Return[Vector2D[a,a+vec,opts]],
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0,0}};
				Return[Vector3D[a,a+vec,opts]]]
				];

Frame2D[delta_] := 
      Graphics[{
     Green,Line[{{-delta/4,0},{delta,0}}],
                    Text["x",{1.01 delta,0}],Line[{{0,-delta/4},{0,delta}}],Text["y",{0, 1.01 delta}]}];

Frame3D[delta_] := 
      Graphics3D[{Green,Line[{{-delta/4,0,0},{delta,0,0}}],Text["x",{1.01 delta,0,0}],Line[{{0,-delta/4,0},{0,delta,0}}],Text["y",{0, 1.01 delta,0}],Line[{{0,0,-delta/4},{0,0,delta}}],Text["z",{0,0,1.01 delta}]}];

View2D[stuff_] := Show[stuff,Frame2D[1]];
View3D[stuff_] := Show[stuff,Frame3D[1]];

SetOptions[Plot3D,ViewPoint->{5,3,1},
	BoxRatios ->Automatic,PlotPoints->30,
	AxesLabel->{"x","y","z"}];
SetOptions[ContourPlot,PlotPoints->30,
	AxesLabel->{"x","y"}];

spaceCurve[ fun:{_, _, _}, ul:{_, u0_, u1_, du_}, opts___] :=
	Show[ Graphics3D[Line[Table[N[fun], ul]]], opts,ViewPoint->{5,3,1} ] /;
			NumberQ[N[u0]] && NumberQ[N[u1]] && NumberQ[N[du]]

spaceCurve[fun_, {u_, u0_, u1_}, opts___] :=
Module[{plotpoints},
    plotpoints = PlotPoints /. {opts} /. Options[Plot3D];
	spaceCurve[ fun, {u, u0, u1, (u1-u0)/(plotpoints-1)}, opts ]
    ];

\(\(sliceSurface[z_, xLo_, xHi_, yLo_, yHi_, zLo_, zHi_, opts___] := 
        Module[{web, box, 
            slicePlane}, 
box = {{xLo -  .01, 
                xHi +  .01}, {yLo -  .01, yHi +  .01}, {zLo -  .01, 
                zHi +  .01}}; 
web = {}; 

          Do[AppendTo[web, 
              spaceCurve[{x, y, z}, {y, yLo, yHi}, PlotRange -> box, 
                DisplayFunction -> 
                  Identity]]; 
slicePlane = 
              Graphics3D[
                Polygon[{{x, yLo, zLo}, {x, yHi, zLo}, {x, yHi, zHi}, {x, 
                      yLo, zHi}, {x, yLo, zLo}}]]; 
            Show[slicePlane, web, opts, ViewPoint -> {5, 3, 1}, 
              PlotRange -> box, 
              DisplayFunction -> $DisplayFunction];
	, {x, xLo, 
              xHi, \((xHi - xLo)\)/5}]; 

          Do[AppendTo[web, 
              spaceCurve[{x, y, z}, {x, xLo, xHi}, PlotRange -> box, 
                DisplayFunction -> 
                  Identity]]; 
slicePlane = 
              Graphics3D[
                Polygon[{{xLo, y, zLo}, {xHi, y, zLo}, {xHi, y, zHi}, {xLo, 
                      y, zHi}, {xLo, y, zLo}}]]; 
            Show[web, slicePlane, opts, ViewPoint -> {5, 3, 1}, 
              PlotRange -> box, 
              DisplayFunction -> $DisplayFunction];
	, {y, yLo, 
              yHi, \((yHi - yLo)\)/5}]; 
Show[web, opts, 
            ViewPoint -> {5, 3, 1}, PlotRange -> box, 
            DisplayFunction -> $DisplayFunction];];\)\)

plot3DContour[f_,xCnt_,yCnt_,scale_,opts___]:= 
Module[{zSs,surface,contour,cont,baseContour,combined,crvs3Dplus,fPlus,crvs3Dminus,fMinus},

zSs = Table[N[f[x,y]],{x,xCnt-scale,xCnt+scale,scale/30},
			{y,yCnt-scale,yCnt+scale,scale/30}];
zMIN = Min[zSs];

contour = ContourPlot[ f[x,y] ,{x,xCnt-scale,xCnt+scale},
			{y,yCnt-scale,yCnt+scale},PlotPoints->30,ContourStyle->Blue,
			DisplayFunction -> Identity];
topo=Show[contour,ContourShading->False,
      ContourStyle->AbsoluteThickness[1],DisplayFunction->Identity];
      			
surface = Plot3D[f[x,y] ,{x,xCnt-scale,xCnt+scale},{y,yCnt-scale,yCnt+scale},opts,ViewPoint->{5,3,1},Mesh->False,BoxRatios->Automatic,DisplayFunction->Identity];

fPlus[x_,y_]=f[x,y]+scale/100;
crvs3Dplus = Show[Graphics3D[First[Graphics[topo]]/.Line[pts_]:>(val=Apply[fPlus,First[pts]];
Line[Map[Append[#,val]&,pts]]
)],DisplayFunction->Identity];

fMinus[x_,y_]=f[x,y]-scale/100;
crvs3Dminus = Show[Graphics3D[First[Graphics[topo]]/.Line[pts_]:>(val=Apply[fMinus,First[pts]];
Line[Map[Append[#,val]&,pts]]
)],DisplayFunction->Identity];

cont = Show[contour,ContourShading->False,
				DisplayFunction -> Identity];

baseContour = First[Graphics[cont]];
					
baseContour = Graphics3D[baseContour /. Line[x_]:>
			Line[Map[Append[#,zMIN]&,x]]];

combined = Show[surface,baseContour,crvs3Dplus,crvs3Dminus,Graphics3D[{Red,AbsolutePointSize[4],Point[{xCnt,yCnt,f[xCnt,yCnt]}]}],DisplayFunction->$DisplayFunction];
Return[combined]
		];

sliceAt[f_,xCnt_,yCnt_,scale_,z_,zMIN_,optns___] :=
	Module[{surface,plane,contour,pLz,contLz,w},

	surface = Plot3D[ f[x,y] , {x,xCnt-scale,xCnt+scale},
						{y,yCnt-scale,yCnt+scale},
				PlotPoints->30,BoxRatios->Automatic,Mesh->False,optns,
				AxesLabel->{"x","y","z"},
				ViewPoint->{5,3,1},
				DisplayFunction->Identity];


	plane[w_] := Graphics3D[Polygon[{{xCnt-scale,yCnt-scale,w},{xCnt+scale,yCnt-scale,w},{xCnt+scale,yCnt+scale,w},{xCnt-scale,yCnt+scale,w},{xCnt-scale,yCnt-scale,w}}]];
						
	contour[w_] := ContourPlot[ f[x,y] ,
		{x,xCnt-scale,xCnt+scale},{y,yCnt-scale,yCnt+scale},
		PlotPoints->30, Contours->{w},ContourStyle->Blue,DisplayFunction->Identity];

				
	pLz = plane[z];

	contLz = Show[Evaluate[contour[z]],
			ContourShading->False,
			DisplayFunction -> Identity];

	contLz = First[Graphics[contLz]];
					
	contLz = Graphics3D[contLz /. Line[x_]:>
			Line[Map[Append[#,zMIN]&,x]]];

Show[surface,pLz,contLz,
		DisplayFunction->$DisplayFunction];]

A more complicated level curve

Let us consider the  relation

This equation defines some curve as a level curve of the function

Here is a countour plot of the function near the origin (0,0) with the one contour line (or level curve) drawn in red.

The intersections of the blue and green lines are the points where the gradient is zero.  Our curve does not contain any of the four points where the gradient of f is zero.

Clear[x, y, F, a, b, c, d, k];
f[x_, y_] =  
        (x^3) -  6 x +  5 y - 3 (y^3) + (y^2);
fX[x_, y_] := Evaluate[(∂_x) f[x, y]];
fY[x_, y_] := Evaluate[(∂_y) f[x, y]];
gradf[x_, y_]  :=  {fX[x, y],  
          fY[x, y]};
a = 4;
plot1  = 
        ContourPlot[f[x, y], {x, -a, a}, {y, -a, a}, 
          Contours ->  50,  
          DisplayFunction  ->  Identity];
plot2  =  
        ContourPlot[
          f[x, y], {x, -a, a}, {y, -a, 
            a}, 
Contours ->  { 4},  
          ContourShading ->  False,  ContourStyle -> Red,  
          DisplayFunction  ->  Identity];
plot3  =  
        ContourPlot[
          fX[x, y], {x, -a, a}, {y, -a, 
            a}, 
Contours ->  { 0},  
          ContourShading ->  False,  ContourStyle -> Green,  
          DisplayFunction  ->  Identity];
plot4  =  
        ContourPlot[
          fY[x, y], {x, -a, a}, {y, -a, 
            a}, 
Contours ->  { 0},  
          ContourShading ->  False,  ContourStyle -> Blue,  
          DisplayFunction  ->  Identity];
Show[{plot1,  plot2,  plot3,  plot4},  
      DisplayFunction  ->  $DisplayFunction]

We don't know how to explicitly solve for y in terms of x or to solve for x in terms of y, but we can see that the red  curve can be split up into pieces such that each piece is the graph of a function -- either y is a function of x, or x is a function of y.    

Moreover, the red curve is smooth, according the the implicit function theorem.