(* :Title: ImplicitPlot *)

(* :Authors: Jerry B. Keiper, Wolfram Research, Inc.,
		contour plot method: Theo Gray, Jerry Glynn, Dan Grayson *)

(* :Summary:
The built-in function Plot requires one to specify an explicit function.
Many simple graphs (e.g., circles, ellipses, etc.) are not functions.
ImplicitPlot allows one to plot figures defined by equations.
*)

(* :Context: Graphics`ImplicitPlot` *)

(* :Package Version: 2.1 *)

(* :History:
	V2.0 by Jerry B. Keiper, April 1991.
	V2.1 Modifications by John M. Novak, April 1992.
*)

(* :Keywords: solution set, graphics *)

(* :Sources: The contour plot alternate method is from:
	Gray, Theodore and Glynn, Jerry, Exploring Mathematics
		with Mathematica, (Addison-Wesley, 1991) *)

(* :Warning: *)

(* :Mathematica Version: 2.0 *)
(* :Copyright: Copyright 1991-1999, Wolfram Research, Inc.*)

(* :Limitation:
	ImplicitPlot relies on Solve for much of the work.  If Solve
	fails, no plot can be made.

	Subscripted variables (e.g., x[1], x[2]) cannot be used.
*)

BeginPackage["Graphics`ImplicitPlot`","Utilities`FilterOptions`"]

ImplicitPlot::usage =
"ImplicitPlot[eqn, {x, a, b}] draws a graph of the set of points
that satisfy the equation eqn.  The variable x is associated with
the horizontal axis and ranges from a to b.  The remaining
variable in the equation is associated with the vertical axis.
ImplicitPlot[eqn, {x, a, x1, x2, ..., b}] allows the user to specify
values of x where special care must be exercised.
ImplicitPlot[{eqn1, eqn2, ...}, {x, a, b}] allows more than one equation
to be plotted, with PlotStyles set as in the Plot function.
ImplicitPlot[eqn, {x, a, b}, {y, a, b}] uses a contour plot method of
generating the plot.  This form does not allow specification
of intermediate points."

Options[ImplicitPlot] =
	{AspectRatio -> Automatic, Axes -> Automatic, AxesLabel -> None,
	AxesOrigin -> Automatic, AxesStyle -> Automatic,
	Background -> Automatic, ColorOutput -> Automatic,
	DefaultColor -> Automatic, Epilog -> {}, Frame -> False,
	FrameLabel -> None, FrameStyle -> Automatic, FrameTicks -> Automatic,
	GridLines -> None, PlotLabel -> None, PlotPoints -> 25,
	PlotRange -> Automatic, PlotRegion -> Automatic,
	PlotStyle -> Automatic, Prolog -> {}, RotateLabel -> True,
	Ticks -> Automatic, DefaultFont :> $DefaultFont,
	DisplayFunction :> $DisplayFunction}

Begin["`Private`"]
Unprotect[ImplicitPlot];

ImplicitPlot[eqns:{__Equal}, xr:{_,_,_},yr:{_,_,_},opts___] :=
	Module[{ps,df},
	{ps,df} = {PlotStyle,DisplayFunction}/.{opts}/.Options[ImplicitPlot];
	ps = cyclestyles[ps,Length[eqns]];
	gr = MapThread[ImplicitPlot[#1,xr,yr,
		ContourStyle->#2,DisplayFunction->Identity,opts]&, {eqns,ps}];
	gr = Select[gr,Head[#] === ContourGraphics &];
	Show[gr,FilterOptions[ContourGraphics,opts,
		Sequence @@ Options[ImplicitPlot]],DisplayFunction->df]/;
			gr =!= {}
	]
		
ImplicitPlot[eqns:{__Equal}, {x_,a_,m___,b_}, opts___] :=
    Module[{ps, df, gr, ln},
	{ps,df} = {PlotStyle,DisplayFunction}/.{opts}/.Options[ImplicitPlot];
	ps = cyclestyles[ps,Length[eqns]];
	gr = MapThread[makegr[#1, {x,a,m,b},
		PlotStyle->#2,DisplayFunction->Identity,opts]&, {eqns,ps}];
	gr = Select[gr, (# =!= $Failed)&];
	Show[Graphics[gr], FilterOptions[Graphics, opts,
		Sequence @@ Options[ImplicitPlot]],DisplayFunction->df]/;
			gr=!={}
	]

ImplicitPlot[lhs_ == rhs_, xr:{_,_,_},yr:{_,_,_},opts___] :=
	With[{ps = PlotStyle/.{opts}/.Options[ImplicitPlot],
			copts = FilterOptions[ContourPlot,opts,
				Sequence @@ Options[ImplicitPlot]]},
	ContourPlot[lhs - rhs,xr,yr,
		copts,
		ContourStyle->ps,
		Contours->{0},
		ContourLines->True,
		ContourShading->False,
		ContourSmoothing->True]
	]

ImplicitPlot[eqn_Equal, {x_,a_,m___,b_}, opts___] :=
    Module[{ps, df, gr},
	{ps,df} = {PlotStyle,DisplayFunction}/.{opts}/.Options[ImplicitPlot];
	gr = makegr[eqn, {x,a,m,b}, PlotStyle->ps, opts];
	Show[Graphics[gr], FilterOptions[Graphics, opts,
		Sequence @@ Options[ImplicitPlot]], DisplayFunction->df]/;
			gr=!=$Failed
	]

cyclestyles[ps_,ln_] :=
	Module[{style = ps},
		If[Head[ps] =!= List,
			style = {ps},
			If[Length[ps] == 0, 
				style = {{}}]
		];
		While[Length[style] < ln, style = Join[style,style]];
		Take[style,ln]
	]

ImplicitPlot::var =
"Equation `1` does not have a single variable other than `2`."

findy[f_, x_] :=
	Module[{nf},
	nf = Select[Union[Cases[f,
			(_Symbol | _[(_?NumberQ)...]),
				Infinity]],
		(!(NumberQ[N[#]] || #===x))&];
	If[Length[nf] == 1,
		nf[[1]],
		(* else *)
		Message[ImplicitPlot::var, f, x];
		$Failed
		]
	]

ImplicitPlot::epfail = "Equation `1` could not be solved for points to plot."

makegr[eqn_Equal, {x_, a_, m___, b_}, opts___] :=
    Module[{f = eqn[[1]] - eqn[[2]], ranges, plots, ar, y},
	If[(y = findy[eqn, x]) === $Failed, Return[$Failed]];
	ranges = Solve[f == 0 && D[f, y] == 0, {x, y}];
    If[ListQ[ranges] && Length[ranges] > 0, ranges = N[x /. ranges]];
	If[!VectorQ[ranges, NumberQ],
	    Message[ImplicitPlot::epfail, eqn];
	    Return[$Failed]];
	ranges = Select[Chop[ranges], FreeQ[#, Complex]&];
	ranges = Sort[Select[ranges, (a < # < b)&]];
	ranges = Union[Sort[Join[ranges, N[{a, m, b}]]]];
	ar = N[b-a]/10^8;
	ranges = Transpose[{Drop[ranges+ar, -1], Drop[ranges-ar, 1]}];
	(* ranges is now a (sorted) list of disjoint intervals with small
	    gaps between them where singularities probably exist. *)
	plots = Map[rangeplot[f, x, y, #, opts]&, ranges]
	];

distx[{x_, y_List}] :=
    Module[{yy = Sort[Select[Chop[y], FreeQ[#, Complex]&]]},
		Transpose[{Table[x, {Length[yy]}], yy}]
	]

evenup[{}] = {};

evenup[xys_] :=
    Module[{ll = Length /@ xys},
		If[Max[ll] == Min[ll], Return[xys]];
		If[Length[ll] < 3, {}, Drop[Drop[xys, 1], -1]]
	]

rangeplot[f_, x_, y_, {a_, b_}, opts___] :=
    Module[{pp, ps, j, multipoints, mdpt, len},
	{pp, ps} = {PlotPoints - 1, PlotStyle} /. {opts} /.
						Options[ImplicitPlot];
	If[ps === Automatic,ps = {}];
	mdpt = (a+b)/2;
	len = (b-a)/2;
	multipoints = N[{#, y /. Solve[f==0 /. x -> #, y]}& /@
			Table[N[mdpt + len Cos[j Pi/pp]], {j, pp, 0, -1}]];
	multipoints = evenup[Select[distx /@ multipoints, (Length[#] > 0)&]];
	If[Length[Dimensions[multipoints]] =!= 3,
		multipoints = {}];
	(* connect the dots to form the various curves *)
	If[Length[multipoints] > 0,
		Flatten[{ps,Line[#]}]& /@ Transpose[multipoints, {2,1,3}],
	    (* else *)
		{}]
	];

Protect[ImplicitPlot];

End[] (* "`Private`" *)

EndPackage[]  (* "Graphics`ImplicitPlot`" *)

(* :Tests: *)
(* :Examples:

ImplicitPlot[x^2 + 2 y^2 == 3, {x, -2, 2}] (* ellipse *)
ImplicitPlot[(x^2 + y^2)^2 == (x^2 - y^2), {x, -2, 2}] (*lemniscate *)
ImplicitPlot[(x^2 + y^2)^2 == 2 x y, {x, -2, 2}] (* lemniscate *)
ImplicitPlot[x^3 + y^3 == 3 x y, {x, -3, 3}] (* folium of Descarte *)
ImplicitPlot[x^2 + y^2 == x y + 3, {x, -3, 3}] (* ellipse *)
ImplicitPlot[x^2 + y^2 == 3 x y + 3, {x, -10, 10},
	PlotRange -> {{-10,10},{-10,10}}] (* hyperpola *)
ImplicitPlot[{(x^2 + y^2)^2 == (x^2 - y^2),
	(x^2 + y^2)^2 == 2 x y}, {x,-2,2},
	PlotStyle->{GrayLevel[0],Hue[0]}] (* combined plots *)
ImplicitPlot[{(x^2 + y^2)^2 == (x^2 - y^2),
	(x^2 + z^2)^2 == 2 x z}, {x,-2,2},
	PlotStyle->{GrayLevel[0],Dashing[{.01}]}] (* combined plots *)
ImplicitPlot[{a == b, x^2 + 2 y^2 == 3}, {x, -1, 1}] (* one bad plot *)
ImplicitPlot[x^2 + y^2 == Pi, {x, -2, 2}] (* OK eqn with 3 symbols *)
ImplicitPlot[Sin[x] == Cos[y], {x, 1.5, Pi/2, 1.7}]
(* contour method *)
ImplicitPlot[Sin[2 x] + Cos[3 y] == 1,{x,-2 Pi,2 Pi},{y,-2 Pi,2 Pi}]
ImplicitPlot[x^2 + x y + y^2 == 1,{x,-2Pi,2Pi},{y,-2Pi,2Pi}]
ImplicitPlot[x^3 + x y + y^2 == 1,{x,-2Pi,2Pi},{y,-2Pi,2Pi}]
ImplicitPlot[x^3 - x^2 == y^2 - y,{x,-1,2},{y,-1,2}]
(* failure cases *)
ImplicitPlot[a == b, {x, -1, 1}] (* bad plot *)
ImplicitPlot[x^y == y^x, {x, -1, 1}] (* bad plot *)
ImplicitPlot[{a == b, c == d}, {x, -1, 1}] (* bad plots *)
ImplicitPlot[x^2 + y^2 == z, {x, -2, 2}] (* bad eqn with 3 vars *)
ImplicitPlot[Sin[x] == y, {x, -3, 3}] (* Solve fails... *)
ImplicitPlot[Sin[x] == Cos[y], {x, -5, 5}]
ImplicitPlot[x^y == y^x, {x, -3, 3}]
*)
