(* :Title: Filled Plot *)

(* :Context: Graphics`FilledPlot` *)

(* :Author: John M. Novak *)

(* :Summary: 
This package allows one to fill the space between a plotted function
and the x-axis or between a pair of plotted functions with a
color.
*)

(* :Package Version: 2.0.2 *)

(* :Mathematica Version: 4.0 *)

(* :History:
	V1.0 by John M. Novak, May 1991.
	V2.0 by John M. Novak, April 1994.
	V2.0.1 by John M. Novak, February 1997, bug fix for DisplayString. 
    V2.0.2 by John M. Novak, January 1998, bug fix for Background->None.
*)

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

(* :Keywords:
	graphics, Plot, axes, area under the curve
*)

(* :Limitations: Uses rectangles to place the axes in front
	of the filled area;  this should work in all cases, but
	is not guaranteed.  (Handled by the AxesFront option.)
*)

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

FilledPlot::usage =
"FilledPlot[function,{var,varmin,varmax}] generates a plot
with the area between the curve and the var axis filled.
FilledPlot[{f1,f2,...},range] generates a plot with the
areas between curve f1 and f2, f2 and f3, etc. filled.
The shade of the fill and other specifications can be
given by the Fills option.";

FilledListPlot::usage =
"FilledListPlot[data] generates a plot with the area between the
x axis and the curve given by the data filled.
FilledListPlot[data1, data2, ...] generates a plot with the areas
between the curves given by data1 and data2, data2 and data3, etc.
filled.  The shade of the fill and other specifications can be given
by the Fills option.";

Fills::usage =
"Fills is an option to FilledPlot.  There are two forms for this
option;  a list of color primitives (i.e., Hue[0], etc.),
or a list of pairs, with the first element being another
pair that specifies two curves, and the other being a color
primitive (i.e., {{2,3},Hue[.3]}, etc.).  In this second form,
each curve is specified by an integer, indicating the
position of the curve in the list of curves given to
FilledPlot.  Also, the symbol Axis (meaning the horizontal or x axis) 
can be specified.  If just the list of colors is given,
this is equivalent to {{{1,2},color1},{{2,3},color2},etc.}.
If there are more colors than curves, extra colors are
ignored. If there are as many colors as curves, the last
color is the fill between the final curve and the Axis.
If there are fewer colors, then only the fills up through
the given colors will be plotted.";

Curves::usage =
"Curves is an option to FilledPlot.  There are three possible values for this
option: None, Front, or Back.  None specifies that the curves are not displayed.
Front specifies that all the curves are drawn in front of the fills.  Back
specifies that curves are drawn in the same layer as the corresponding fill.
Thus, specifying Back may cause some curves to be covered by later fills.";

Front::usage =
"Front is a value for the option Curves.";

Back::usage =
"Back is a value for the option Curves.";

Axis::usage =
"Axis is used in the value for the option Fills to specify that
one of the curves to be used in a fill is the horizontal (or x) axis.";

AxesFront::usage =
"AxesFront is an option to Graphics that allows axes to be drawn
in front of a graphic.  This option may be set to True or False.";

(* Backwards compatibility *)
ListFilledPlot::usage =
"ListFilledPlot is the old name for FilledListPlot, and is now obsolete."

ListFilledPlot = FilledListPlot;

Begin["`Private`"]

(* numberQ utility *)
numberQ[n_] := NumberQ[N[n]]

(* Curves option takes None, Front, or Back. *)

Options[FilledPlot] = Join[{
		Fills->Automatic,
		Curves->Back,
		AxesFront -> True},
	Options[Plot]];

FilledPlot::badfill =
"The Fills option has been given an incorrect form: `1`.";

FilledPlot::badcurv =
"The Curves option has been given bad value `1`; using
'Back' in its place.";

FilledPlot[funcs_List,{x_Symbol,xmin_,xmax_},opts___] :=
	Module[{pl,lines,fills,gopts,disp,curve,
			ln = Length[funcs],n,gr,part, front},
		{fills, curve, front, disp} = {Fills, Curves,
				AxesFront, DisplayFunction}/.
			Flatten[{opts,Options[FilledPlot]}];
		pl = Plot[funcs,{x,xmin,xmax},
			DisplayFunction->Identity,
			Evaluate[FilterOptions[Plot,##]& @@
			            Flatten[{opts, Options[FilledPlot]}]
			]
		     ];
		lines = Map[Last,pl[[1]]];
		gopts = List @@ Drop[pl,1];
		curve = checkcurve[curve, FilledPlot];
		fills = checkfills[fills, ln, FilledPlot];
		polys = Map[dopoly[#,lines,{xmin,xmax}]&,fills];
		gr = Switch[curve,
			Back,{pl[[1]], MapThread[{#1,
				Map[If[# =!= Axis && # =!= 0,
					pl[[1,#]],{}]&,#2[[1]]]}&,{polys,fills}] },
			None,polys,
			Front,{polys,pl[[1]]}];
		Show[Graphics[gr,Flatten[{DisplayFunction -> disp,
			gopts, AxesFront -> front}]]]
	]

FilledPlot[func_,rng_List,opts___] :=
	FilledPlot[{func},rng,opts]

(* FilledListPlot *)
FilledListPlot::badfill =
"The Fills option has been given an incorrect form: `1`.";

FilledListPlot::badcurv =
"The Curves option has been given bad value `1`; using
'Back' in its place.";

Options[FilledListPlot] = Join[{
	Curves -> Back,
	Fills -> Automatic,
	AxesFront -> True,
	PlotStyle -> Automatic},
	Options[Graphics]
];

SetOptions[FilledListPlot, Axes -> Automatic];

FilledListPlot[lists:(({_?numberQ..} | {{_?numberQ, _?numberQ}..})..),
				opts:((_Rule | _RuleDelayed)...)] :=
	Module[{lines, curve, fills, front, polys, gopts, gr, xmin, xmax,
				xcoords, linestyle},
		{curve, fills, front, linestyle} = {Curves, Fills,
					AxesFront, PlotStyle}/.
				Flatten[{opts,Options[FilledListPlot]}];
		lines = Map[If[MatchQ[#, {_?numberQ..}],
						Transpose[{Range[Length[#]], #}],
						#]&,
					{lists}
				];
		xcoords = First[Transpose[Join @@ lines]];
		lines = Map[splitatcomplex, Transpose[{Map[Line, lines]}]];
		{xmin, xmax} = {Min[xcoords], Max[xcoords]};
		If[Head[linestyle] =!= List || Length[linestyle] === 0,
			linestyle = {GrayLevel[0.]}];
		linestyle = CycleValues[linestyle, Length[lines]];
		curve = checkcurve[curve, FilledListPlot];
		fills = checkfills[fills, Length[lines], FilledListPlot];
		gopts = FilterOptions[Graphics, ##]& @@
		        Flatten[{opts, AxesFront -> front,
					Options[FilledListPlot]}];
		polys = Map[dopoly[#,lines,{xmin,xmax}]&,fills];
		lines = Map[Flatten,Transpose[{linestyle, lines}]];
		gr = Switch[curve,
			Back,{lines, MapThread[{#1,
				Map[If[# =!= Axis && # =!= 0,
					lines[[#]],{}]&,#2[[1]]]}&,{polys,fills}]},
			None,polys,
			Front,{polys,lines}];
		Show[Graphics[gr, gopts]]
	]

		
(* Auxiliary Functions *)

dopoly[{{first_,second_},style_},lines_,rng_] :=
	Flatten[{style,
		Outer[linestopolygon,getline[first,lines,rng],
		getline[second,lines,rng]]}]

getline[0,lines_,rng_] := getline[Axis,lines,rng]

getline[Axis,lines_,{min_,max_}] := {Line[N[{{min,0},{max,0}}]]}

getline[ln_,lines_,rng_] := lines[[ln]]/;
	ln <= Length[lines] && ln >= 1

getline[_,_,{min_,max_}] := {Line[N[{{min,0},{max,0}}]]}

(* delete to intersecting point, leading; then, delete
	to intersecting point, trailing; then, create polygon
	(no checking needed.) *)

linestopolygon[Line[line1_List],Line[line2_List]] :=
	{}/;First[Last[line1]] <= First[First[line2]] ||
		First[Last[line2]] <= First[First[line1]]

linestopolygon[Line[line1_List],Line[line2_List]] :=
	With[{f1 = First[First[line1]],f2 = First[First[line2]]},
		If[f1 > f2,
			linestopolygon[Line[line1],chopline[line2,f1,True]],
			linestopolygon[chopline[line1,f2,True],Line[line2]]]
	]/;First[First[line1]] != First[First[line2]]

linestopolygon[Line[line1_List],Line[line2_List]] :=
	With[{f1 = First[Last[line1]],f2 = First[Last[line2]]},
		If[f1 < f2,
			linestopolygon[Line[line1],chopline[line2,f1,False]],
			linestopolygon[chopline[line1,f2,False],Line[line2]]]
	]/; First[Last[line1]] != First[Last[line2]]

linestopolygon[Line[line1_List],Line[line2_List]] :=
	Polygon[Join[line1,Reverse[line2]]]

chopline[line_,pt_,front_:True] :=
	Module[{pl = Partition[line,2,1],x1,y1,x2,y2,pos},
		pos = Position[pl,
			x_?(#[[1,1]] < pt && #[[2,1]] >= pt &),
				{1},Heads->False][[1,1]];
		{{x1,y1},{x2,y2}} = pl[[pos]];
		If[TrueQ[front],
			Line[Prepend[Drop[line,pos],
				{pt,(y2 - y1)/(x2 - x1) (pt - x1) + y1}]],
			Line[Append[Drop[line,-(Length[line] - pos)],
				{pt,(y2 - y1)/(x2 - x1) (pt - x1) + y1}]]]
	]

(* check form of curve *)
checkcurve[curve_, callingfun_] :=
	If[!MatchQ[curve,(Front | Back | None)],
		Message[callingfun::badcurv,curve];Back,
		curve
	]

(* Checking the form of fills *)

checkfills[fills_, ln_, callingfun_] :=
	Module[{fillsout = fills},
		If[MatchQ[fillsout,_Hue | _GrayLevel | _RGBColor | _CMYKColor],
			fillsout = {fillsout}];
		If[!MatchQ[fillsout,({{{(_Integer | Axis),(_Integer | Axis)},
				(_Hue | _GrayLevel | _RGBColor | _CMYKColor)}..} |
				{(_Hue | _GrayLevel | _RGBColor | _CMYKColor)..} |
				Automatic)],
			Message[callingfun::badfill,fillsout]; fillsout = Automatic];
		If[MatchQ[fillsout,{(_Hue | _GrayLevel | _RGBColor | _CMYKColor)..}],
			If[Length[fillsout] > ln,
				fillsout = Take[fillsout,ln]];
			part = Partition[Range[1,ln],2,1];
			If[Length[fillsout] === ln,
				PrependTo[part,{1,Axis}]];
			fillsout = Transpose[{part,fillsout}]
		];
		If[fillsout === Automatic,
			If[ln === 1,
				fillsout = {{{1,Axis}, GrayLevel[.5]}},
				fillsout = Transpose[{
					Partition[Range[1,ln],2,1],
					Table[Hue[n/ln],{n,1,ln-1}]}]
			]
		];
		fillsout
	]

(* For FilledListPlot; split line segments containig complex values
   at the points with the complex numbers. Equivalent behavior to
   FilledPlot when it hits complex values (handled automatically by
   Plot. *)
splitatcomplex[{Line[{{_,_Complex}..}]}] := {}

splitatcomplex[{l:Line[{___, {_, _Complex}, ___}]}] :=
   Flatten[{l//.{
        Line[{{_,_Complex}, a___}] :> Line[{a}],
        Line[{a___, {_, _Complex}}] :> Line[{a}],
        Line[{a__, {_,_Complex}, b__}] :> {Line[{a}], Line[{b}]}
   }}]

splitatcomplex[any_] := any

(* The following is a useful internal utility function to be
used when you have a list of values that need to be cycled to
some length (as PlotStyle works in assigning styles to lines
in a plot).  The list is the list of values to be cycled, the
integer is the number of elements you want in the final list. *)

CycleValues[{},_] := {}

CycleValues[list_List, n_Integer] :=
	Module[{hold = list},
		While[Length[hold] < n,hold = Join[hold,hold]];
		Take[hold,n]
	]

CycleValues[item_,n_] := CycleValues[{item},n]


(* Adding the AxesFront option to Graphics *)

Unprotect[Graphics];

Options[Graphics] = Append[Options[Graphics],
	AxesFront -> False];

Protect[Graphics];

Unprotect[Display];

Display[out_, gr_, args___] :=
	Module[{tgr, opts, defaultaf, posn},
		defaultaf = AxesFront/.Options[Graphics];
		Unprotect[Graphics];
		Options[Graphics] = 
			Select[Options[Graphics], First[#] =!= AxesFront &];
		Display[out,
			MapAt[layeraxes[#,defaultaf]&,
				gr,
				Position[gr, _Graphics]],
			args
		];
		Unprotect[Graphics];
		Options[Graphics] =
			Append[Options[Graphics], AxesFront -> defaultaf];
		Protect[Graphics];
		gr
	]/;Count[gr,AxesFront,Infinity] != 0 ||
		TrueQ[AxesFront/.Options[Graphics]]

Display[out_, gr:(_Graphics | _GraphicsArray | _Sound), args___] :=
	Module[{result, defaultaf},
		defaultaf = AxesFront/.Options[Graphics];
		Unprotect[Graphics];
		Options[Graphics] =
			Select[Options[Graphics], First[#] =!= AxesFront &];
		result = Display[out,gr,args];
		Unprotect[Graphics];
		Options[Graphics] =
			Append[Options[Graphics], AxesFront -> defaultaf];
		Protect[Graphics];
		result
	]/;(AxesFront/.Options[Graphics]) =!= AxesFront

Protect[Display];

Unprotect[DisplayString];

DisplayString[gr_, args___] :=
  Module[{tgr, opts, defaultaf, posn, strn},
	defaultaf = AxesFront/.Options[Graphics];
	Unprotect[Graphics];
	Options[Graphics] = 
		Select[Options[Graphics], First[#] =!= AxesFront &];
	strn = DisplayString[
	  MapAt[layeraxes[#,defaultaf]&,
			gr,
			Position[gr, _Graphics]],
	  args
	];
	Unprotect[Graphics];
	Options[Graphics] =
		Append[Options[Graphics], AxesFront -> defaultaf];
	Protect[Graphics];
	strn
  ]/;Count[gr,AxesFront,Infinity] != 0 ||
	TrueQ[AxesFront/.Options[Graphics]]

DisplayString[gr:(_Graphics | _GraphicsArray | _Sound), args___] :=
  Module[{result, defaultaf},
	defaultaf = AxesFront/.Options[Graphics];
	Unprotect[Graphics];
	Options[Graphics] =
		Select[Options[Graphics], First[#] =!= AxesFront &];
	result = DisplayString[gr, args];
	Unprotect[Graphics];
	Options[Graphics] =
		Append[Options[Graphics], AxesFront -> defaultaf];
	Protect[Graphics];
	result
  ]/;(AxesFront/.Options[Graphics]) =!= AxesFront

Protect[DisplayString]

$V3Fix = True;

Unprotect[FullGraphics];

FullGraphics[gr_Graphics] :=
	Module[{oldopts, result},
		oldopts = Options[Graphics];
		defAF = AxesFront/.Options[Graphics];
		Unprotect[Graphics];
		Options[Graphics] =
			fgallGraphics[
				Select[Options[Graphics],
					First[#] =!= AxesFront &],
				defAF
			];
		result = FullGraphics[fgallGraphics[gr, defAF]];
		Unprotect[Graphics];
		Options[Graphics] = oldopts;
		Protect[Graphics];
		result
	]/;Count[gr,AxesFront,Infinity] != 0 ||
		(AxesFront/.Options[Graphics]) =!= AxesFront

Protect[FullGraphics];

fgallGraphics[gr_, def_] :=
	MapAt[fgoneGraphics[#,def]&, gr, Position[gr, _Graphics]]

fgoneGraphics[Graphics[gr_, opts___], def_] :=
	Module[{fopts, pr, ngr, agr},
		fopts = Select[Flatten[{opts}], First[#] =!= AxesFront &];
		pr = PlotRange[Graphics[gr, fopts]];
		ngr = FullGraphics[Graphics[gr, fopts]];
		agr = FullGraphics[Graphics[{},
			Flatten[{PlotRange -> pr, fopts}]]];
		Graphics[{ngr[[1]], agr[[1]]}, Flatten[List @@ Rest[ngr]]]
	]/;TrueQ[AxesFront/.Flatten[{opts, AxesFront -> def}]]

fgoneGraphics[Graphics[gr_, opts___], _] :=
	FullGraphics[Graphics[gr, 
		Select[Flatten[{opts}], First[#] =!= AxesFront &]
	]]

Unprotect[FullAxes];

FullAxes[gr_] :=
	Module[{oldopts, result},
		oldopts = Options[Graphics];
		Unprotect[Graphics];
		Options[Graphics] =
			removeGraphicsOption[
				Select[Options[Graphics],
					First[#] =!= AxesFront &],
				AxesFront
			];
		result = FullAxes[removeGraphicsOption[gr,AxesFront]];
		Unprotect[Graphics];
		Options[Graphics] = oldopts;
		Protect[Graphics];
		result
	]/;Count[gr,AxesFront,Infinity] != 0 ||
		(AxesFront/.Options[Graphics]) =!= AxesFront


Protect[FullAxes];

Unprotect[PlotRange];

PlotRange[gr_] :=
	Module[{oldopts, result},
		oldopts = Options[Graphics];
		Unprotect[Graphics];
		Options[Graphics] =
			removeGraphicsOption[
				Select[Options[Graphics],
					First[#] =!= AxesFront &],
				AxesFront
			];
		result = PlotRange[removeGraphicsOption[gr,AxesFront]];
		Unprotect[Graphics];
		Options[Graphics] = oldopts;
		Protect[Graphics];
		result
	]/;Count[gr,AxesFront,Infinity] != 0 ||
		(AxesFront/.Options[Graphics]) =!= AxesFront

Protect[PlotRange];

removeGraphicsOption[gr_, opt_] :=
	MapAt[removegoption[#, opt]&,
		gr,
		Position[gr, _Graphics]
	]

removegoption[Graphics[gr_, opts___], ropt_] :=
	Graphics[gr, Select[Flatten[{opts}], First[#] =!= ropt &]]

(* warning - I may be taking advantage of a bug in clearing
	 the epilog and prolog of the overlay. It appears that these are not
	 counted in plotrange generation, so we should be OK, and this eliminates
	 bugs of the epilog and prolog being rendered in the wrong coordinate
	 system a second time... *)
    
layeraxes[Graphics[pic_, opts___], def_] :=
	Module[{fopts, pr, asp, back, defc},
		fopts = Select[Flatten[{opts}], First[#] =!= AxesFront &];
		{pr, asp, back, defc} = FullOptions[Graphics[pic,fopts],
			{PlotRange,AspectRatio, Background, DefaultColor}];
		Graphics[{
			Rectangle[{0,0},{1,1},
				Graphics[pic,
					Join[{PlotRange -> pr}, fopts]
				]
			],
			Rectangle[{0,0},{1,1},
				Graphics[{PointSize[0]},
					Join[{PlotRange -> pr,
						Epilog -> {}, Prolog -> {},
						DefaultColor -> fixdefcolor[back, defc],
					    Background -> None}, fopts]
				]
			]
		}, Join[{ Axes -> False,
				Frame -> False,
				GridLines -> None,
				PlotLabel -> None,
				AspectRatio -> asp,
				Epilog -> {},
				Prolog -> {},
				PlotRange -> {{0,1},{0,1}} },
			fopts
			]
		]
	]/;TrueQ[AxesFront/.Flatten[{opts, AxesFront -> def}]]

layeraxes[Graphics[pic_, opts___], def_] :=
	Graphics[pic, Select[Flatten[{opts}],First[#] =!= AxesFront &]]

fixdefcolor[GrayLevel[l_?(#<.5 &)], Automatic] :=
   GrayLevel[1.]

fixdefcolor[_,def_] := def

End[]

EndPackage[]
