(* :Package Version: 2.0, 12.30.04 *)

(* :Copyright: Copyright 2004, Selwyn Hollis *)

(* :Name: Graphics`DEGraphics` *)
 
(* :Author: Selwyn Hollis *)

(*:Summary:

This package defines the functions NDPlot, PhasePlot, DEPlot, TimeStatePlot,
PoincareTimeSection, ViewProjections, and PlotImplicit.

This is the last version of DEGraphics.m that will
work properly with Mathematica versions 4.x and 5.0.

*)

(* Shadowing prevention *)

If[Not[ValueQ[$degraphicspackagenames]],
$degraphicspackagenames ={"NDPlot", "DEPlot", "PhasePlot", "TimeStatePlot",
	"PoincareTimeSection", "ViewProjections", "PlotImplicit","Arrows","ScaleArrows",
	"ArrowColor","InitialPoints","Gap","MaxCurves","StepSize","ShowStartPoints",
	"ShowEquilibriumPoints","PrintEquilibriumPoints","UseVectorField",
	"ShowNullclines","Isoclines","ShadowColors","ShadowPositions","OuterFrame","OuterFrameStyle",
	"IsoclinePoints","IsoclineStyle","NullclineStyle"};
$trashvaluelesssymbols[name$list:{_String ..}]:=
    Select[ToExpression/@
       Intersection[name$list,ToString/@ToExpression/@name$list],({UpValues[#],DownValues[#]}==={{},{}}&)];
Remove@@($trashvaluelesssymbols[$degraphicspackagenames])
];
Remove[$trashvaluelesssymbols, name$list];


BeginPackage["DiffEqs`DEGraphics`", "Graphics`PlotField`", "Graphics`Arrow`",
"Graphics`Graphics3D`", "Graphics`Graphics`","Utilities`FilterOptions`", "Graphics`Common`GraphicsCommon`"]

Unprotect["DiffEqs`DEGraphics`*"];
ClearAll["DiffEqs`DEGraphics`*"];

$DiffEqsVersion="2.0, 12.30.04";

NDPlot::usage = "NDPlot[DEs, fns, {t,tmin,tmax}, options].
Uses NDSolve to compute a numerical solution of a system of first-order differential equations
and, if there are 1, 2, or 3 equations, plots the solution. 
Returns {solution, -Graphics-}. Argument syntax is the same as for NDSolve.
Uses the same options as NDSolve, Plot, ParametricPlot, and ParametricPlot3D.";

PhasePlot::usage = "PhasePlot[{f,g}, {t,tmin,tmax}, {x,xmin,xmax}, {y,ymin,ymax}, options].
By default, four solution curves are drawn. The InitialPoints option can be
used to specify a list of starting points for curves drawn by PhasePlot. Each initial
point may give rise to multiple curves as specified by the MaxCurves option.";

DEPlot::usage = "DEPlot[f, {t,tmin,tmax}, {y,ymin,ymax}, options].
By default, four solution curves are drawn. The InitialPoints option can be
used to specify a list of starting points for curves drawn by DEPlot. Each initial
point may give rise to multiple curves as specified by the MaxCurves option.";

InterpolatingFunctionPlot::usage = " "

TimeStatePlot::usage = "TimeStatePlot[{x[t],y[t]}, {t,tmin,tmax}, options]. Creates a 3D plot
of the time-state trajectory (t,x,y) with projections onto the tx, ty, and xy planes.";

PoincareTimeSection::usage = "PoincareTimeSection[{f,g}, {t,tmin,tmax,dt}, {x,x0}, {y,y0}, options].
Creates a Poincare time section plot for the system x'=f, x(t0)=x0, y'=g, y(t0)=y0, with time step dt.";

ViewProjections::usage = "ViewProjections[{f[t],g[t],h[t]}, {t,tmin,tmax}, {x,y,z}, options].
Creates a GraphicsArray of 2-D projections of the curve {f[t],g[t],h[t]} onto the xy, xz, and yz coordinate planes.";

PlotImplicit::usage = "PlotImplicit[F[t,y], {t,tmin,tmax}, {y,ymin,ymax}, options].
Creates a contour plot of the function F."

Arrows::usage = "An option of PhasePlot and DEPlot. Specifies the number of arrows to draw in each
direction. The default is {11,9}.";

ScaleArrows::usage = "An option of PhasePlot and DEPlot. A scale factor for the length of the
arrows. The default is 0.25.";

ArrowColor::usage = "An option of PhasePlot and DEPlot. Specifies the color of the arrows. The
default is RGBColor[0,.6,.6].";

InitialPoints::usage = "An option of PhasePlot and DEPlot. A list of one or more initial points through
which solution curves are drawn. (See MaxCurves.)" 

Gap::usage = "An option of PhasePlot and DEPlot. Controls the spacing between initial points on 
automatically generated curves. The default is 0.1. The actual spacing is Gap*diag, where
diag=Sqrt[(xmax-xmin)^2+(ymax-ymin)^2]";

MaxCurves::usage = "An option of PhasePlot and DEPlot. Specifies the number of automatically drawn
curves arising from each seed point. The default is {1}.";

StepSize::usage = "An option of PhasePlot and DEPlot. Specifies the stepsize used in the numerical
method or, in case Method->Automatic, the length of the intervals over which NDSolve is used to 
compute the solution. The default is 0.1.";

ShowStartPoints::usage = "An option of PhasePlot and DEPlot. If set to True, then the starting
point for each curve is plotted as a red dot. The default is False.";

ShowEquilibriumPoints::usage = "An option of PhasePlot and DEPlot. If set to True, then equilibrium
points are plotted as blue dots. The default is False. Equilibrium points are computed
with Solve, which will fail in many cases";

PrintEquilibriumPoints::usage = "An option of PhasePlot. The default is False.";

ReturnFunctions::usage = "An option of PhasePlot. The default is False.";

UseVectorField::usage = "An option of PhasePlot and DEPlot. If set to True, then the vector field
is drawn rather than the direction field. The default is False."

ShowDensityPlot::usage = "An option of PhasePlot and DEPlot. The default is False."

ShowNullclines::usage = "An option of PhasePlot. The default is False.";

Isoclines::usage = "An option of DEPlot. The default is False. If set to True, then 10 isoclines
are drawn. If Isoclines->{m1,m2,...,mn}, then isoclines f[t,y]=mi, i=1,...,n, are drawn.";

ShadowColors::usage = "An option of TimeStatePlot. Specifies the colors of the projections onto 
the tx, ty, and xy planes. The default is {RGBColor[1,0,0], RGBColor[0,1,0], RGBColor[0,0,1]}.";

ShadowPositions::usage = "An option of TimeStatePlot. Specifies the positions of the planes containing 
the planar projections. Contains values passed on to Shadow (in Graphics`Graphics3D`) as XShadowPosition,
YShadowPosition, ZShadowPosition. The default is {-.67,.8,-.67}.";

OuterFrame::usage = "An option of ViewProjections. The default is False.";
OuterFrameStyle::usage = "An option of ViewProjections. The default is None.";

Options[PhasePlot] := {ArrowColor-> RGBColor[0, 0.6, 0.6], Arrows-> {11, 9}, AspectRatio-> 1/GoldenRatio,
	Axes-> True, AxesLabel-> None, AxesOrigin-> Automatic, AxesStyle-> Automatic,
	Background-> Automatic, ColorOutput-> Automatic, DefaultColor-> Automatic, Epilog-> {}, Frame-> False, 
	FrameLabel-> None, FrameStyle-> Automatic, FrameTicks-> Automatic, Gap-> 0.1, GridLines-> None, 
	ImageSize-> Automatic, InitialPoints-> Automatic, MaxCurves-> {1}, 
	Method-> Automatic, MaxStepSize->1, MaxSteps->10^6, MaxBend->1, NullclineStyle->GrayLevel[.6], 
	NullclinePoints->67, PlotJoined-> True, PlotLabel-> None, PlotRange-> All, PlotRegion-> Automatic, 
	PlotStyle-> Thickness[0.005], Prolog-> {}, ReturnFunctions->False, RotateLabel-> True, ScaleArrows-> 1, 
	ShowEquilibriumPoints-> False, ShowField-> True, ShowNullclines->False, ReturnFunctions->False, 
	ShowStartPoints-> False, ShowLastPoints-> False, Ticks-> Automatic, 
	UseVectorField-> False, ShowDensityPlot->False, ColorFunction->GrayLevel, ContourShading->False, 
	DefaultFont :> $DefaultFont, DisplayFunction :> $DisplayFunction, 
	FormatType :> $FormatType, TextStyle :> $TextStyle};
  
Options[DEPlot] := {ArrowColor-> RGBColor[0, 0.6, 0.6], Arrows-> {11, 9}, AspectRatio-> 1/GoldenRatio,
	Axes-> True, AxesLabel-> None, AxesOrigin-> Automatic, AxesStyle-> Automatic,
	Background-> Automatic, ColorOutput-> Automatic, DefaultColor-> Automatic, Epilog-> {}, Frame-> False, 
	FrameLabel-> None, FrameStyle-> Automatic, FrameTicks-> Automatic, Gap-> 0.1, GridLines-> None, 
	ImageSize-> Automatic, InitialPoints-> Automatic, Isoclines-> False, IsoclineStyle->GrayLevel[.6],
	IsoclinePoints->67, MaxCurves-> {1}, Method-> Automatic, ContourShading->False, ColorFunction->GrayLevel,
	PlotJoined-> True, PlotLabel-> None, PlotRange-> All, PlotRegion-> Automatic, ShowDensityPlot->False,
	PlotStyle-> Thickness[0.006], Prolog-> {}, RotateLabel-> True, ReturnFunctions->False, ScaleArrows-> 1, 
	ShowEquilibriumPoints-> False, ShowField-> True, ReturnFunctions->False, ShowStartPoints-> False,
	Ticks-> Automatic, UseVectorField-> False, DefaultFont :> $DefaultFont,
	DisplayFunction :> $DisplayFunction, FormatType :> $FormatType, TextStyle :> $TextStyle};
  
Options[TimeStatePlot] :={AspectRatio-> Automatic, Axes-> True, AxesEdge-> Automatic, AxesLabel-> None,
	AxesStyle-> Automatic, Background-> Automatic, Boxed-> True, BoxRatios-> {2, 1, 1},
	BoxStyle-> Automatic, ColorOutput-> Automatic, Compiled-> True, DefaultColor-> Automatic, Epilog-> {},
	FaceGrids-> None, ImageSize-> Automatic, Lighting-> True, PlotLabel-> None, PlotPoints-> Automatic,
	PlotRange-> All, PlotRegion-> {{0.025, 0.975}, {0.025, 0.975}}, PlotStyle-> Automatic, Prolog-> {},
	RenderAll-> True, ShadowColors-> {RGBColor[1, 0, 0], RGBColor[0, 1, 0], RGBColor[0, 0, 1]},
	ShadowPositions-> {-0.67, 0.8, -0.67}, Ticks-> Automatic, ViewCenter-> Automatic, ViewPoint-> {2, -3, 1},
	DefaultFont :> $DefaultFont, DisplayFunction :> $DisplayFunction, FormatType :> $FormatType,
	TextStyle :> $TextStyle};
 
Options[PlotImplicit] := 
	{AspectRatio-> 1/GoldenRatio, Axes-> False, AxesLabel-> None, AxesOrigin-> Automatic,
	AxesStyle-> Automatic, Background-> Automatic, ColorFunction->(Hue[ArcTan[#]]&), ColorFunctionScaling-> True,
	ColorOutput-> Automatic, Compiled-> True, ContourLines-> True, Contours-> 10, ContourShading-> False,
	ContourSmoothing-> True,  ContourStyle-> Automatic, DefaultColor-> Automatic, Epilog-> {},
	Frame-> True, FrameLabel-> None, FrameStyle-> Automatic, FrameTicks-> Automatic, ImageSize-> Automatic,
	Mesh-> False, MeshStyle-> Automatic, PlotLabel-> None, PlotPoints-> 15, PlotRange-> Automatic, 
	PlotRegion-> Automatic, Prolog-> {}, RotateLabel-> True, ShowDensity-> False, Ticks-> Automatic,
	DefaultFont:> $DefaultFont, DisplayFunction:> $DisplayFunction, FormatType:> $FormatType,
	TextStyle:> $TextStyle};

Options[ViewProjections] := {AspectRatio-> 1/GoldenRatio, Axes-> Automatic, 
	AxesOrigin-> Automatic, AxesStyle-> Automatic, Background-> Automatic, ColorOutput-> Automatic,
	Compiled-> True, DefaultColor-> Automatic, Epilog-> {}, Frame-> False, FrameLabel-> None,
	FrameStyle-> Automatic, FrameTicks-> Automatic, GraphicsSpacing-> 0.1, GridLines-> None,
	ImageSize->{72*6, Automatic}, MaxBend-> 10., OuterFrame-> False, OuterFrameStyle-> {}, PlotDivision-> 30.,
	PlotLabel-> None, PlotPoints-> 25, PlotRange-> Automatic, PlotRegion-> Automatic, PlotStyle-> Automatic,
	Prolog-> {}, Ticks-> Automatic, DefaultFont:> $DefaultFont, DisplayFunction:> $DisplayFunction,
	FormatType:> $FormatType, TextStyle:> $TextStyle};

Options[NDPlot]= Union[Options[NDSolve], Options[Plot], Options[ParametricPlot], Options[ParametricPlot3D]];
	
Options[PoincareTimeSection] = Union[Options[ListPlot], Options[NDSolve]];
 
IsoclinePoints::usage = "An option of DEPlot.";
NullclinePoints::usage = "An option of PhasePlot.";
IsoclineStyle::usage = "An option of DEPlot.";
NullclineStyle::usage = "An option of PhasePlot.";

Begin["`Private`"]


Off[Graphics::"gptn"];


(*** Definition of NDPlot ***)
         
NDPlot[de:{_Equal..}, y:(_[t_]|{_[t_]..}|{{_[t_]..}..}), {t_Symbol, t1_?NumericQ, t2_?NumericQ}, opts___?OptionQ] :=
    Module[{soln, ndsolveout, dim = Last[Dimensions[y]]}, 
         If[Head[ndsolveout = NDSolve[de, Flatten[y], {t, t1, t2}, FilterOptions[NDSolve, opts], 
                FilterOptions[NDSolve, Options[NDPlot]]]] === NDSolve,
           Return[$Failed], 
           soln = y /. Flatten[ndsolveout];
           If[t2 != Flatten[{soln}][[1,0,1,1,2]], Return[$Failed], 
              soln = If[Head[soln] =!= List, {soln}, soln];
              {soln, Which[
                     dim == 1,
                     Plot[Evaluate[soln], {t,t1,t2},
                          Evaluate[FilterOptions[Plot, opts], FilterOptions[Plot,Options[NDPlot]]]],
                     dim == 2, 
                     ParametricPlot[Evaluate[soln], {t,t1,t2},
                          Evaluate[FilterOptions[ParametricPlot, opts], FilterOptions[ParametricPlot, Options[NDPlot]]]],
                     dim == 3,
                     ParametricPlot3D[Evaluate[Append[soln, PlotStyle/.{opts}/.{PlotStyle->{}}]], {t, t1, t2}, 
                          Evaluate[FilterOptions[ParametricPlot3D, opts], AspectRatio -> Automatic, 
                          FilterOptions[ParametricPlot3D, Options[NDPlot]]]]]}]
            ]
           ] /; MemberQ[{1,2,3}, Last[Dimensions[y]]] && IntegerQ[Length[de]/Length[Flatten[y]]]       
 

 

(*** Definition of PhasePlot ***)


PhasePlot::"tint"="Time interval must include zero.";
PhasePlot::"num"="Vector field does not provide numerical values.";
PhasePlot::"method"="Method must be Euler, Heun, RungeKutta, DormandPrince, or Automatic.";
PhasePlot::"badseeds"="InitialPoints must be Automatic, None, or a list of one or more points. Reverting to default.";


PhasePlot[{ff_, gg_}, {t_Symbol, t1_?NumericQ, t2_?NumericQ},
		  {x_Symbol, x1_?NumericQ, x2_?NumericQ}, {y_Symbol, y1_?NumericQ, y2_?NumericQ},
			opts___Rule] :=
Module[{arrows, autonomous, drawfield, endpoints, eqpoints, eqpts={}, f, field, plot1,  
	fuzz, g, gap, h, hh, i, isoclines, isopoints, isos, isostyle, lastpoints, k, m, 
	maxcurves, method, n, ndsopts, nonautonomouspair, nullstyle, nullpoints, numcurves, perp, plotrange, 
	plotregion, printeqpts, retfuncs, scalearrows, scalefunction, seeds, showeqpts, showdensity, plotstyle,
	showfield, showisoclines, showlast, shownullclines, slopes, start, startpoints, step, 
	stepsize, tmax, tmin, u, usevectorfield, v, w, xmax, xmaxminus, xmaxplus, xmin, 
	xminminus, xminplus, xnulls, ymax, ymaxminus, ymaxplus, ymin, yminminus, yminplus, 
	ynulls, a, b, c, d, fc, gc, uc, DPstep, allpoints, returnpoints},

If[t1<=0<=t2, {}, Message[PhasePlot::"tint"];Return[$Failed]];
If[And@@NumericQ/@({ff,gg}/.{t->0, x->N[x1], y->N[y1]}), {}, Message[PhasePlot::"num"];Return[$Failed]];

f = ff/.{x[t]-> x, y[t]-> y}; g = gg/.{x[t]-> x, y[t]-> y};

fc=Compile[{tt,xx,yy},Evaluate[f//.{t-> tt, x-> xx, y-> yy}]];
gc=Compile[{tt,xx,yy},Evaluate[g//.{t-> tt, x-> xx, y-> yy}]];

showfield = ShowField /. {opts} /. Options[PhasePlot];
usevectorfield = UseVectorField /. {opts} /. Options[PhasePlot];
usevectorfield = If[usevectorfield=!=True && usevectorfield=!=False, False, usevectorfield];
arrows = Arrows /. {opts} /. Options[PhasePlot];
arrowcolor = ArrowColor /. {opts} /. Options[PhasePlot];
scalearrows = .25 ScaleArrows /. {opts} /. Options[PhasePlot];
scalefunction = If[usevectorfield, None, (1.&)];

showstartpoints = ShowStartPoints /. {opts} /. If[t===x, Options[DEPlot], Options[PhasePlot]];
showeqpts = ShowEquilibriumPoints /. {opts} /. Options[PhasePlot];
showlast = ShowLastPoints /. {opts} /. Options[PhasePlot];
showlast = If[showlast=!=True, False, True];
shownullclines = ShowNullclines /. {opts} /. Options[PhasePlot];
showdensity = ShowDensityPlot /. {opts} /. Options[PhasePlot];
retfuncs = ReturnFunctions /. {opts} /. Options[PhasePlot];

isoclines = Isoclines /. {opts};
isoclines = If[isoclines===Isoclines, False, isoclines];
showisoclines = If[isoclines===False || isoclines===None, False, True];
showisoclines = If[x=!=t, False, showisoclines];
slopes = If[isoclines===True || isoclines===Automatic, 10, isoclines];

isostyle = IsoclineStyle /. {opts} /. Options[DEPlot];
nullstyle = NullclineStyle /. {opts} /. Options[PhasePlot];
nullstyle = If[Head[nullstyle]=!=List || VectorQ[nullstyle], {nullstyle, nullstyle}, nullstyle];
isopoints = IsoclinePoints /. {opts} /. Options[DEPlot];
nullpoints = NullclinePoints /. {opts} /. Options[PhasePlot];

(* secret options *) 

	printeqpts = TrueQ[PrintEquilibriumPoints /. {opts} /. {PrintEquilibriumPoints->False}];

maxcurves = MaxCurves /. {opts} /. Options[PhasePlot];
maxcurves= If[Head[maxcurves]===Integer, {maxcurves}, maxcurves];
gap = Gap /. {opts} /. Options[PhasePlot];

seeds = InitialPoints /. {opts} /. Options[PhasePlot];
seeds=If[maxcurves===None || maxcurves===0 || maxcurves==={0} || seeds==={} || seeds==={{}}, None, seeds];
seeds=If[seeds===Automatic,
N[{{(3x1+x2)/4,(2y1+y2)/3},{(3x1+x2)/4,(y1+2y2)/3},{(x1+3x2)/4,(2y1+y2)/3},{(x1+3x2)/4,(y1+2y2)/3}}], seeds];
seeds=If[seeds=!=None && Not[MatrixQ[seeds]], Message[PhasePlot::"badseeds"];
	N[{{(3x1+x2)/4,(2y1+y2)/3},{(3x1+x2)/4,(y1+2y2)/3},{(x1+3x2)/4,(2y1+y2)/3},{(x1+3x2)/4,(y1+2y2)/3}}],
	seeds];

stepsize = Stepsize /. {opts};
stepsize=If[NumberQ[N[stepsize]], stepsize, StepSize/.{opts}/.Options[PhasePlot]];


(* Process PlotRange. *)

{tmin, tmax, xmin, xmax, ymin, ymax} = N[{t1, t2, x1, x2, y1, y2}];

plotrange= PlotRange/.{opts}/.{PlotRange->N[{{x1, x2}, {y1, y2}}]};

{{xmin, xmax}, {ymin, ymax}}= N[
	Which[	MatrixQ[plotrange],
				plotrange,
			Head[plotrange]===List && plotrange[[2]]===Automatic,
				{plotrange[[1]],{y1,y2}},
			Head[plotrange]===List && plotrange[[1]]===Automatic,
				{{x1,x2},plotrange[[2]]},
			VectorQ[plotrange],
				{{x1,x2},plotrange},
			1==1,
				{{x1, x2}, {y1, y2}}
		]];

{{xminplus, xmaxminus},{yminplus, ymaxminus}}= {{xmin,xmax},{ymin,ymax}}.{{0.975,0.025},{0.025,0.975}};
{{xminminus,xmaxplus},{yminminus,ymaxplus}}= {{xmin,xmax},{ymin,ymax}}.{{1.025,-0.025},{-0.025,1.025}};

plotrange= PlotRange/.{opts}/.{PlotRange->{{xminminus,xmaxplus},{yminminus,ymaxplus}}};
plotrange= 
	Which[	Head[plotrange]===List && plotrange[[2]]===Automatic,
				{plotrange[[1]],{yminminus,ymaxplus}},
			Head[plotrange]===List && plotrange[[1]]===Automatic,
				{{xminminus,xmaxplus},plotrange[[2]]},
			VectorQ[plotrange],
				{{xminminus,xmaxplus},plotrange},
			plotrange===Automatic,
				{{xminminus,xmaxplus},{yminminus,ymaxplus}},
			1==1,
				plotrange
		];
		
{a,b,c,d}={Max[x1,xmin],Min[x2,xmax],Max[y1,ymin],Min[y2,ymax]};

(* Finished processing options. Now make the field plot. *)

nonautonomouspair = If[t===x, False, (D[f,t]=!= 0 || D[g,t]=!= 0)];
drawfield=If[nonautonomouspair||arrows===0||arrows==={0,0}||arrows==={0}, False, True];
showeqpts=If[nonautonomouspair, False, showeqpts];
shownullclines=If[nonautonomouspair, False, shownullclines];

showfield=If[showfield=!=True, False, True];    

field= If[drawfield && showfield,
	PlotVectorField[{f, g}, {x, xminplus, xmaxminus}, {y, yminplus, ymaxminus},
		PlotPoints-> arrows, ScaleFactor-> scalearrows,
		ScaleFunction-> scalefunction, ColorFunction-> (arrowcolor &),
		DisplayFunction-> Identity, Graphics`Arrow`HeadCenter->.67],
	{}];

(* Find equilibria if requested	*)
	
If[TrueQ[showeqpts],
    
    If[t=!=x,
      eqpts=N@Chop[{x,y}/.NSolve[{f==0,g==0},{x,y}]];
      eqpts=
        If[eqpts==={x,y},{},
          Select[eqpts,(Element[N[#[[1]]],Reals] && Element[N[#[[2]]], Reals])&]];
      If[printeqpts&&eqpts=!={}, Print["Equilibrium points are ", eqpts]],
      
      eqpts=N@Chop[{0.,y}/.NSolve[{g==0},y]];
      If[eqpts==={0.,y}, eqpts={},
      				eqpts=Select[eqpts, Element[#[[2]], Reals]&];
          			If[printeqpts && eqpts =!= {0.,y}, Print["Equilibria are ", eqpts] ]
          	]
      ];
    eqpts=If[MatchQ[eqpts,{{_Real,_Real}..}],
    			DeleteCases[eqpts,_?((#[[1]]<xmin||#[[1]]>xmax||#[[2]]<ymin||#[[2]]>ymax)&)],
    			{}]
    ]; 

eqpoints=If[TrueQ[showeqpts] && eqpts=!={},
				Graphics[Join[{PointSize[0.02], RGBColor[0,0,1]}, Point/@eqpts]], {} ];
	
(*
	Beginning of Big If:
*)

If[seeds===None, plot1={},

(*
   Otherwise generate starting points based on seeds, maxcurves, and gap.
*)

SetAttributes[LessEqual,Listable];

w = gap*Sqrt[(xmax - xmin)^2 + (ymax - ymin)^2];

numcurves =
If[Length[maxcurves] == 1, Table[maxcurves[[1]], {Length[seeds]}], maxcurves];

perp = {-g,f}/Sqrt[f^2 + g^2 + 10^(-6)] /. t-> 0;

start = {};
Do[
	start = Join[{seeds[[i]]}, start, {seeds[[i]]}];
	m = Length[start];

	While[
		Length[start] - m + 1 <= numcurves[[i]]/2 &&
			First[Union[{a,c} <= Last[start] <= {b,d}]],
		v = perp /. {x-> Last[start][[1]], y-> Last[start][[2]]};
		start = Append[start, Last[start] + w*v]];
	
	m = Length[start];

	While[
		Length[start] - m + 1 <= numcurves[[i]]/2 &&
			First[Union[{a,c} <= First[start] <= {b,d}]],
		v = perp /. {x-> start[[1,1]], y-> start[[1,2]]};
		start = Prepend[start, start[[1]] - w*v]],
	{i, 1, Length[seeds]}];
	
start= DeleteCases[start, _?((#[[1]]<a || #[[1]]>b || #[[2]]<c || #[[2]]>d)&)];

start = Union[start, If[t===x, eqpts, {}]]; 

(*
	Build plots.
*)

n = Length[start];
h = N[stepsize];

ndsolns = Join[
If[t2==0, {},
   If[x=!=t,
	{x[t],y[t]}/.First@NDSolve[
		{x'[t]==(ff/.{x->x[t],y->y[t]}), y'[t]==(gg/.{x->x[t],y->y[t]}), Thread[{x[0],y[0]}==#]},
		{x[t],y[t]}, {t,0,t2}, StoppingTest->Not[x1<x[t]<x2 && y1<y[t]<y2],
                               FilterOptions[NDSolve,opts]] & /@start,
	y[t]/.First@NDSolve[
		{y'[t]==(gg/.{y->y[t]}), Thread[y[0]==#]},
		y[t], {t,0,t2}, StoppingTest->Not[x1<t<x2 && y1<y[t]<y2],
				FilterOptions[NDSolve,opts,Options[PhasePlot]]] & /@ Last[Transpose[start]]	
	]],
If[t1==0, {},
  If[x=!=t,
	{x[t],y[t]}/.First@NDSolve[
		{x'[t]==(ff/.{x->x[t],y->y[t]}), y'[t]==(gg/.{x->x[t],y->y[t]}), Thread[{x[0],y[0]}==#]},
		{x[t],y[t]}, {t,t1,0}, StoppingTest->Not[x1<x[t]<x2 && y1<y[t]<y2],
                                FilterOptions[NDSolve,opts]] & /@start,
	y[t]/.First@NDSolve[
		{y'[t]==(gg/.{y->y[t]}), Thread[y[0]==#]},
		y[t], {t,t1,0}, StoppingTest->Not[x1<t<x2 && y1<y[t]<y2],
				FilterOptions[NDSolve,opts,Options[PhasePlot]]] & /@ Last[Transpose[start]]	
	]]
	];

Block[{$DisplayFunction=Identity},
	plot1 = InterpolatingFunctionPlot[ndsolns, PlotStyle->(PlotStyle/.Join[{opts}, Options[PhasePlot]]), opts]/.
		pts:{{_Real, _Real} ..} :> (clip[#,{{x1,x2},{y1,y2}}]&[pts])]; (*clip is defined below*)

lastpoints = DeleteCases[First[plot1] /. {___,{Line[pts : {{_, _} ..}]}} :> Last[pts],
             _?((#[[1]] <= a || #[[1]] >= b || #[[2]] <= c || #[[2]] >= d) &)];


ClearAttributes[LessEqual,Listable];	

]; (* end of big If *)

	
(*
	Create extra stuff for the plot as requested.
*)
	
endpoints=If[showlast, Graphics[Join[{PointSize[0.015], RGBColor[.65,0,.65]}, Point /@ lastpoints]], {}];

startpoints = If[showstartpoints, Graphics[Join[{PointSize[0.015], RGBColor[0,.7,0]}, Point /@ start]], {}];	

xnulls:= If[shownullclines,
		ContourPlot[Evaluate[g], {x,xmin,xmax}, {y,ymin,ymax},
			Frame->False, Contours->{0}, ContourShading->False,
			ContourStyle->First@nullstyle, PlotPoints->nullpoints, DisplayFunction->Identity],
		{}];
		
ynulls:= If[shownullclines,
		ContourPlot[Evaluate[f], {x,xmin,xmax}, {y,ymin,ymax},
			Frame->False, Contours->{0}, ContourShading->False,
			ContourStyle->nullstyle[[2]], PlotPoints->nullpoints, DisplayFunction->Identity],
		{}];

isos := If[showisoclines,
		ContourPlot[Evaluate[g], {x,xmin,xmax}, {y,ymin,ymax},
			Frame->False, Contours->slopes, ContourStyle->isostyle, PlotPoints->isopoints,
			DisplayFunction->Identity, Evaluate@FilterOptions[ContourPlot, opts, Options[DEPlot]] ],
		{}]; 
		
densplot := If[showdensity,
		DensityPlot[Evaluate[If[t===x, g^2, Abs[g]/(Abs[f]+Abs[g]+10^(-4))]], {x,xmin,xmax}, {y,ymin,ymax},
			Mesh->False, PlotPoints->200, PlotRange->Automatic, DisplayFunction->Identity,
			Evaluate@FilterOptions[DensityPlot, opts],
			Evaluate@FilterOptions[DensityPlot, Options[PhasePlot]] ],
		{}];

(*
	Put it all together in the Big Show. 
*)

If[x===t,

If[TrueQ@retfuncs,
	
	{ndsolns, Show[{densplot,isos, field, plot1, endpoints, startpoints, eqpoints},
					FilterOptions[Graphics, opts],
					Sequence@@If[showdensity||TrueQ[ContourShading/.{opts}/.If[t===x,Options[DEPlot],Options[PhasePlot]]], {Axes->False, Frame->True},{}],
					FilterOptions[Graphics, Options[DEPlot]]]},
	
	Show[{densplot,isos, field, plot1, endpoints, startpoints, eqpoints},
		FilterOptions[Graphics, opts], 
		Sequence@@If[showdensity||TrueQ[ContourShading/.{opts}/.If[t===x,Options[DEPlot],Options[PhasePlot]]], {Axes->False, Frame->True},{}],
		FilterOptions[Graphics, Options[DEPlot]]]
	],
	
If[TrueQ@retfuncs,
	
	{ndsolns, Show[{densplot, xnulls, ynulls, field, plot1, endpoints, startpoints, eqpoints, isos},
					FilterOptions[Graphics, opts],
					Sequence@@If[showdensity||TrueQ[ContourShading/.{opts}/.If[t===x,Options[DEPlot],Options[PhasePlot]]], {Axes->False, Frame->True},{}],
					FilterOptions[Graphics, Options[PhasePlot]]]},
	
	Show[{densplot, xnulls, ynulls, field, plot1, endpoints, startpoints, eqpoints},	
		FilterOptions[Graphics, opts],
		Sequence@@If[showdensity||TrueQ[ContourShading/.{opts}/.If[t===x,Options[DEPlot],Options[PhasePlot]]], {Axes->False, Frame->True},{}],
		FilterOptions[Graphics, Options[PhasePlot]] ]
] ] 
	
];  (* end of module *)



InterpolatingFunctionPlot[f:{(_)..} /; MatchQ[f, {InterpolatingFunction[__][_]..}], (opts___)?OptionQ] := 
   DisplayTogether[(Plot[Evaluate[#1], Evaluate[#1 /. _[{{t1_, t2_}}, __][t_] :> {t, t1, t2}],
   Evaluate@FilterOptions[Plot,opts]] & ) /@ f, FilterOptions[DisplayTogether,opts]]; 
   
InterpolatingFunctionPlot[f_ /; MatchQ[f, InterpolatingFunction[__][_]], (opts___)?OptionQ] := 
   InterpolatingFunctionPlot[{f}, opts]; 
   
InterpolatingFunctionPlot[f:{{_, _}..} /; MatchQ[f, {{InterpolatingFunction[__][_]..}..}], (opts___)?OptionQ] := 
  DisplayTogether[
    (ParametricPlot[Evaluate[#1], Evaluate[First[#1 /. _[{{t1_, t2_}}, __][t_] :> {t, t1, t2}]],
    Evaluate@FilterOptions[ParametricPlot,opts]] & ) /@ f, FilterOptions[DisplayTogether,opts]]; 
    
InterpolatingFunctionPlot[f:{{_, _, _}..} /; MatchQ[f, {{InterpolatingFunction[__][_]..}..}], (opts___)?OptionQ] := 
  DisplayTogether[
    (ParametricPlot3D[Evaluate[#1], Evaluate[First[#1 /. _[{{t1_, t2_}}, __][t_] :> {t, t1, t2}]],
    Evaluate@FilterOptions[ParametricPlot3D,opts]] & ) /@ f, FilterOptions[DisplayTogether,opts]]; 


(** line clipping function used in PhasePlot **) 
	
	Attributes[pull]={HoldFirst};

	pull[z_, end_, xory_, k_, bounds_]:= With[
		{r=(bounds[[xory,k]]-z[[2end,xory]])/(z[[end,xory]]-z[[2end,xory]])},
		z[[end]]=r*z[[end]]+(1-r)*z[[2end]]];
											 (* Danke schoen, Herr Hartmut Wolf! *)
	clip[zIn_, bounds_]:=
	Module[{z=zIn},
	Scan[Function[end,
			If[z[[end,1]]<bounds[[1,1]], pull[z,end,1,1,bounds],
			If[z[[end,1]]>bounds[[1,-1]], pull[z,end,1,-1,bounds]]];
			If[z[[end,-1]]<bounds[[-1,1]], pull[z,end,-1,1,bounds],
			If[z[[end,-1]]>bounds[[-1,-1]], pull[z,end,-1,-1,bounds]]];],
			{1,-1}]; z];

(*** definition of DEPlot ***)

DEPlot[g_, {t_Symbol, t1_?NumericQ, t2_?NumericQ},
			{y_Symbol, y1_?NumericQ, y2_?NumericQ}, opts___Rule] :=
	Module[{ta, tb, x1, x2, seeds, plotrange},
	DEPlot::"tint"="Time interval must include zero.";
		If[t1<=0<=t2, {},
		Message[DEPlot::"tint"];Return[$Failed]];
		
		DEPlot::"num"="Slope function does not provide numerical values.";
		If[Not[NumericQ[g/.{t->0., y->N[y1]}]], Message[DEPlot::"num"]; Return[$Failed]];

		maxcurves = MaxCurves /. {opts} /. Options[DEPlot];
		seeds = InitialPoints /. {opts} /. Options[DEPlot];
		plotrange = PlotRange /. {opts} /. Options[DEPlot];
		seeds=If[maxcurves===None || maxcurves===0 || maxcurves==={0} || seeds==={} || seeds==={{}}, None, seeds];
		seeds=Which[
      		seeds===Automatic && t1<=0<=t2,
      			Transpose[{{0,0,0,0},{y1,y2}.{{1,2,3,4},{4,3,2,1}}/5.}],
      		seeds===Automatic,
      			Transpose[{(t1+t2){1,1,1,1}/2.,{y1,y2}.{{1,2,3,4},{4,3,2,1}}/5.}],
     		 1==1,
      			seeds];
		ta=If[seeds===None, 0, Min[First@(seeds//Transpose)]];
		tb=If[seeds===None, 0, Max[First@(seeds//Transpose)]];
		{x1,x2}=
			Which[MatrixQ[plotrange],
					plotrange[[1]], 
				Head[plotrange]===List && plotrange[[2]]===Automatic,
					plotrange[[1]],
				1==1,
					{t1,t2}
				];
		PhasePlot[{1, g}, {t, t1-tb, t2-ta}, {t,x1,x2}, {y,y1,y2}, InitialPoints->seeds, opts]
	];

   
    
(* definition of TimeStatePlot *)

TimeStatePlot[{xsoln_, ysoln_}, {t_Symbol, tmin_?NumericQ, tmax_?NumericQ}, opts___?OptionQ] := 
   Module[{curve, curvex, curvey, curvez, xshadow, yshadow, zshadow, colors,
   			shadows, style, trip, tplotrange, t1, t2, tt},
   		
		TimeStatePlot::"num"="Functions do not provide numerical values.";
		If[And@@NumericQ/@({xsoln,ysoln}/.{t->0.}),{}, Message[TimeStatePlot::"num"];Return[$Failed]];   		

		colors = ShadowColors /. {opts} /. Options[TimeStatePlot]; 
		shadows = ShadowPositions /. {opts} /. Options[TimeStatePlot]; 
		style = PlotStyle /. {opts} /. Options[TimeStatePlot]; 
		style=If[style===Automatic, {}, style];
		
		tplotrange= PlotRange /. {opts} /.{PlotRange->{tmin,tmax}};
		tplotrange= If[Length[tplotrange]==3, tplotrange[[1]], tplotrange];
		tplotrange= If[tplotrange===Automatic,{tmin,tmax}, tplotrange];
		{t1,t2}=tplotrange;
		tt= (t2(t-tmin)-t1(t-tmax))/(tmax-tmin);
		
	trip := {t, xsoln, ysoln};

	curve := ParametricPlot3D[Evaluate[{{tt, 0, 0}, Append[trip, Flatten[style]]}], 
       {t, tmin, tmax}, DisplayFunction-> Identity,
       Evaluate[FilterOptions[ParametricPlot3D,opts]],
       Evaluate[FilterOptions[ParametricPlot3D,Options[TimeStatePlot]]]]; 
     
	curvex := ParametricPlot3D[Evaluate[Append[trip, Flatten[{style, colors[[1]]}]]],
		{t, tmin, tmax}, DisplayFunction-> Identity,
		Evaluate[FilterOptions[ParametricPlot3D,opts]],
		Evaluate[FilterOptions[ParametricPlot3D,Options[TimeStatePlot]]]]; 
	curvey := ParametricPlot3D[Evaluate[{{tt, 0, 0}, Append[trip, Flatten[{style, colors[[2]]}]]}],
		{t, tmin, tmax}, DisplayFunction-> Identity,
		Evaluate[FilterOptions[ParametricPlot3D,opts]],
		Evaluate[FilterOptions[ParametricPlot3D,Options[TimeStatePlot]]]]; 
	curvez := ParametricPlot3D[Evaluate[{{tt, 0, 0}, Append[trip, Flatten[{style, colors[[3]]}]]}],
		{t, tmin, tmax}, DisplayFunction-> Identity,
		Evaluate[FilterOptions[ParametricPlot3D,opts]],
		Evaluate[FilterOptions[ParametricPlot3D,Options[TimeStatePlot]]]]; 
     
	xshadow := Shadow[curvex, XShadowPosition-> shadows[[1]], ZShadow-> False, YShadow-> False,
		DisplayFunction-> Identity]; 
	yshadow := Shadow[curvey, YShadowPosition-> shadows[[2]], ZShadow-> False, XShadow-> False,
		DisplayFunction-> Identity]; 
	zshadow := Shadow[curvez, ZShadowPosition-> shadows[[3]], XShadow-> False, YShadow-> False,
		 DisplayFunction-> Identity]; 

	Show[xshadow, yshadow, zshadow, curve,
		 FilterOptions[Graphics3D,opts],
		 FilterOptions[Graphics3D,Options[TimeStatePlot]]
		 ]

]; 

(*** Definition of PoincareTimeSection ***)

PoincareTimeSection[{f_, g_}, {t_Symbol, t0_?NumericQ, tend_?NumericQ, dt_?NumericQ},
				{x_Symbol, x0_?NumericQ}, {y_Symbol, y0_?NumericQ}, opts___Rule] := 
  Module[{xsoln, ysoln, maxsteps}, 
  
  PoincareTimeSection::"num"="Vector field does not provide numerical values.";
  If[And@@NumericQ/@({f,g}/.{t->0.,x->N[x0],y->N[y0]}),{}, Message[PoincareTimeSection::"num"];Return[$Failed]];
  
  maxsteps= MaxSteps/.{opts}/.{MaxSteps->100(tend-t0)};
  {xsoln, ysoln} = {x, y} /.
    NDSolve[{x'[t] == (f /. {x-> x[t], y-> y[t]}), y'[t] == (g /. {x-> x[t], y-> y[t]}),
    x[0] == x0, y[0] == y0}, {x, y}, {t, t0, tend}, MaxSteps->maxsteps,
    FilterOptions[NDSolve,opts], FilterOptions[NDSolve, Options[PoincareTimeSection]]][[1]];
    ListPlot[Table[{xsoln[t], ysoln[t]}, {t, t0, tend, dt}],
    FilterOptions[ListPlot,opts],
    FilterOptions[ListPlot, Options[PoincareTimeSection]]]
  ];
    

(*** Definition of ViewProjections ***)
    
ViewProjections[{f_, g_, h_}, {t_Symbol, t0_?NumericQ, tmax_?NumericQ}, {x_, y_, z_}, opts___Rule] := 
  Module[
  {xyProj, xzProj, yzProj, xx, yy, zz, outerframe, outerframestyle, optpp}, 
  
    ViewProjections::"num"="Functions not provide numerical values.";
    If[And@@NumericQ/@({f,g,h}/.{t->N[t0]}),{}, Message[ViewProjections::"num"];Return[$Failed]];

	xx = StyleForm[x, FontSlant-> "Italic"];
	yy = StyleForm[y, FontSlant-> "Italic"]; 
	zz = StyleForm[z, FontSlant-> "Italic"];
	
	optpp = FilterOptions[ParametricPlot, opts];
   
	outerframe = OuterFrame /.{opts} /.{OuterFrame-> False}; 
	outerframe= If[outerframe=!=True, False, True];
	outerframestyle = OuterFrameStyle /.{opts} /.{OuterFrameStyle-> None};
	
			
	xyProj = ParametricPlot[Evaluate[{f,g}], {t, t0, tmax}, PlotLabel-> None, 
		AxesLabel-> {xx, yy}, DisplayFunction-> Identity, Evaluate[optpp]]; 
	xzProj = ParametricPlot[Evaluate[{f,h}], {t, t0, tmax}, PlotLabel-> None, 
		AxesLabel-> {xx, zz}, DisplayFunction-> Identity, Evaluate[optpp]]; 
	yzProj = ParametricPlot[Evaluate[{g,h}], {t, t0, tmax}, PlotLabel-> None, 
		AxesLabel-> {yy, zz}, DisplayFunction-> Identity, Evaluate[optpp]];
      
    Show[GraphicsArray[{xyProj, xzProj, yzProj}],
    	Frame->outerframe, FrameTicks->None,
    	FrameStyle->outerframestyle, Axes->False, AspectRatio->Automatic,
    	FilterOptions[GraphicsArray, opts],
    	FilterOptions[GraphicsArray, Options[ViewProjections]]]
  ];
    
    
(*** Definition of PlotImplicit ***)

PlotImplicit[f_, {x_Symbol, xmin_?NumericQ, xmax_?NumericQ},
			{y_Symbol, ymin_?NumericQ, ymax_?NumericQ}, opts___Rule] := 
  Module[{curves, density, plotpoints, showdensity},

	plotpoints= PlotPoints/. {opts}/. Options[PlotImplicit];
    showdensity = ShowDensity /. {opts} /. Options[PlotImplicit]; 
    showdensity = If[showdensity =!= True && showdensity =!= False, False, showdensity];
	
   	curves = ContourPlot[Evaluate[f], {x, xmin, xmax}, {y, ymin, ymax},
   		DisplayFunction-> Identity,
   		Evaluate[FilterOptions[ContourPlot,opts]],
     	Evaluate[FilterOptions[ContourPlot,Options[PlotImplicit]]]];
	density =
		If[showdensity, 
			DensityPlot[Evaluate[f], {x, xmin, xmax}, {y, ymin, ymax}, 
      			DisplayFunction-> Identity,
				Evaluate[FilterOptions[DensityPlot,opts]],
				PlotPoints-> plotpoints+25,
				Evaluate[FilterOptions[DensityPlot,Options[PlotImplicit]]]],
		{}];
	 
    Show[density, curves,
    	FilterOptions[Graphics, opts],
    	FilterOptions[Graphics, Options[PlotImplicit]]]
]


End[]

EndPackage[]

Protect @@ ToExpression /@ $degraphicspackagenames ;
Unprotect[StepSize];