(*
	file:   Chaos.m
	author: Hinke Osinga, The Geometry Center
	data:   April 15, 1998
*)


Print[""]
Print["     ---  Copyright by Hinke Osinga, 1998.  ---"]
Print[""]
Print["Chaos.m is a package designed for:"]
Print["      Math5337: One-Dimensional Dynamical Systems"]
Print["Course material can be found at:"]
Print["        http://www.geom.umn.edu/~math5337/ds/."]
Print["This package produces pictures of orbits, either by plotting"]
Print["iterates versus time, or by graphical iteration."]
Print["It can also produce bifurcation diagrams."]
Print[""]


<< Logistic.m

BeginPackage["Chaos`"]

Orbit::usage =
"Orbit computes the orbit of x0 up to the n-th iterate, where n = its,
and returns a list that contains the iterates."

PlotIterates::usage =
"PlotIterates first computes and then plots the orbit versus time."

PlotGraphical::usage =
"PlotGraphical first computes and then shows the graphical iteration
of the orbit. The last few iterates are shown in red. The graph of the
function (in blue) and the diagonal line y = x are drawn as well." 

HigherIts::usage =
"HigherIts[n] computes and shwos the graphical iteration of an orbit
for the n-th iterate of the map. Again, the last few iterates are
shown in red. The graph of the function (in blue) and the diagonal
line y = x are drawn as well. Add graphs of other iterates in
different colors with AddGraph."

ShowGraph::usage =
"ShowGraph[n, color] draws the graph of the n-th iterate in the
specified color together with the diagonal line, but without
iterates. Choose from colors: red, blue, green, yellow, purple, and
pink, or specify RGBColor[a, b, c] with a, b, c in [0,1]."

AddGraph::usage =
"AddGraph[picture, n, color] draws the graph of the n-th iterate in
the specified color inside the given picture. Choose from colors: red,
blue, green, yellow, purple, and pink, or specify RGBColor[a, b, c]
with a, b, c in [0,1]." 

OrbitDiagram::usage =
"OrbitDiagram[transient, actual] draws an orbit diagram by computing
the orbit of an arbitrary point up to 'actual' iterates, but only
starts displaying after 'transient' iterates."


Begin["`Private`"]

Orbit :=
 Module[{q, i, orb},
	q = Global`x0;
	orb = {q};
	For[i = 1, i <= Global`its, i++,
		q = Global`DynSys[q];
		(* Print[i, " : ", q]; *)
		orb = Append[orb, q]
	];
	Return[orb]
 ];

CritOrbit[p_, n_] :=
 Module[{q, i, orb},
	q = p;
	orb = {q};
	For[i = 1, i <= n, i++,
		q = Global`DynSys[q];
		(* Print[i, " : ", q]; *)
		orb = Append[orb, q]
	];
	Return[orb]
 ];

HigherOrb[n_] :=
 Module[{Fn, q, i, orb},
	Fn = Nest[Global`DynSys, #, n] &;
	q = Global`x0;
	orb = {q};
	For[i = 1, i <= Global`its, i++,
		q = Fn[q];
		(* Print[i, " : ", q]; *)
		orb = Append[orb, q]
	];
	Return[orb]
 ];

PlotIterates :=
 Module[{orb, pts},
	orb = Orbit;
	pts = Table[{i-1, orb[[i]]}, {i, Length[orb]}];
	Show[{ListPlot[pts, Prolog -> AbsolutePointSize[3], 
		            DisplayFunction -> Identity],
	      ListPlot[pts, PlotJoined -> True,
		            DisplayFunction -> Identity]},
	     PlotRange -> {{0, Global`its}, Global`xrange},
	     DisplayFunction -> $DisplayFunction]
 ];

GraphFunction[n_, color_] :=
	Plot[Evaluate[Nest[Global`DynSys, x, n], 
		      {x, Global`xrange[[1]], Global`xrange[[2]]}],
	     PlotStyle -> {AbsoluteThickness[2], color},
	     DisplayFunction -> Identity];

ShowGraph[n_, color_] :=
 Module[{xmin, xmax},
	xmin = Min[Global`xrange[[1]], Global`xrange[[2]]];
	xmax = Max[Global`xrange[[1]], Global`xrange[[2]]];
	Show[{Graphics[Line[{{xmin, xmin}, {xmax, xmax}}]],
	      GraphFunction[n, color]},
	     PlotRange -> {Global`xrange, Global`xrange},
	     DisplayFunction -> $DisplayFunction,
	     AspectRatio -> 1,
	     Frame -> True]
 ];

AddGraph[fig_, n_, color_] :=
	Show[{fig, GraphFunction[n, color]}, 
	     Frame -> True,
	     DisplayFunction -> $DisplayFunction]
	
AddGraph[n_, color_] := AddGraph[%, n, color];

PlotGraphical :=
 Module[{orb, pts, end, xmin, xmax},
	xmin = Min[Global`xrange[[1]], Global`xrange[[2]]];
	xmax = Max[Global`xrange[[1]], Global`xrange[[2]]];
	orb = Orbit;
	pts = Table[{{orb[[i]], orb[[i]]}, {orb[[i]], orb[[i+1]]}},
		    {i, Length[orb]-1}];
	pts = Flatten[pts, 1];
	end = Floor[0.6 Length[pts]];
	Show[{Graphics[Line[{{xmin, xmin}, {xmax, xmax}}]],
	      Graphics[{AbsolutePointSize[3], Point[First[pts]]}],
	      ListPlot[Take[pts, end], PlotJoined -> True,
	      	       PlotStyle -> GrayLevel[0.5], 
		       DisplayFunction -> Identity],
	      GraphFunction[1, RGBColor[0, 0, 1]],
	      ListPlot[Take[pts, {end, Length[pts]}], PlotJoined ->True,
	      	       PlotStyle -> RGBColor[1, 0, 0], 
		       DisplayFunction -> Identity],
	      ListPlot[Select[Take[pts, {end, Length[pts]}], #[[1]]==#[[2]]&], 
	      	       PlotStyle -> {RGBColor[1, 0, 0], AbsolutePointSize[3]}, 
		       DisplayFunction -> Identity]},
	     PlotRange -> {Global`xrange, Global`xrange},
	     DisplayFunction -> $DisplayFunction,
	     AspectRatio -> 1,
	     Frame -> True]
 ];

HigherIts[n_] :=
 Module[{orb, pts, end, xmin, xmax},
	xmin = Min[Global`xrange[[1]], Global`xrange[[2]]];
	xmax = Max[Global`xrange[[1]], Global`xrange[[2]]];
	orb = HigherOrb[n];
	pts = Table[{{orb[[i]], orb[[i]]}, {orb[[i]], orb[[i+1]]}},
		    {i, Length[orb]-1}];
	pts = Flatten[pts, 1];
	end = Floor[0.6 Length[pts]];
	Show[{Graphics[Line[{{xmin, xmin}, {xmax, xmax}}]],
	      Graphics[{AbsolutePointSize[3], Point[First[pts]]}],
	      ListPlot[Take[pts, end], PlotJoined -> True,
	      	       PlotStyle -> GrayLevel[0.5], 
		       DisplayFunction -> Identity],
	      GraphFunction[n, RGBColor[0, 0, 1]],
	      ListPlot[Take[pts, {end, Length[pts]}], PlotJoined ->True,
	      	       PlotStyle -> RGBColor[1, 0, 0], 
		       DisplayFunction -> Identity],
	      ListPlot[Select[Take[pts, {end, Length[pts]}], #[[1]]==#[[2]]&], 
	      	       PlotStyle -> {RGBColor[1, 0, 0], AbsolutePointSize[3]}, 
		       DisplayFunction -> Identity]},
	     PlotRange -> {Global`xrange, Global`xrange},
	     DisplayFunction -> $DisplayFunction,
	     AspectRatio -> 1,
	     Frame -> True]
 ];

OrbitDiagram[transient_, actual_, {lr_, xr_}] :=
 Module[{orb, pts, ptlst, param, step, stop, i, storelambda},
	storelambda = Global`lambda;
	step = N[(lr[[2]] - lr[[1]])/200];
	stop = Abs[step] / 2;
	ptlst = {};
	param = lr[[1]];
	i = 0;
	While[Abs[lr[[2]] - param] > stop, 
		If[Mod[i, 10] == 0, Print[200 - i, " orbits to go"]];
		Global`lambda = param;
		orb = CritOrbit[Global`critical, actual];
		pts = Map[{AbsolutePointSize[1], Point[{param, #}]} &, 
			        Take[orb, {transient, Length[orb]}]];
		ptlst = Append[ptlst, pts];
		param += step;
		i++
	];
	ptlst = Flatten[ptlst, 1];
	Global`lambda = storelambda;
	Show[Graphics[ptlst],
	     PlotRange -> {lr, xr},
	     Frame -> True,
	     AspectRatio -> 1]
 ];

OrbitDiagram[transient_, actual_] := 
	OrbitDiagram[transient, actual, {Global`lrange, All}];

End[ ]


EndPackage[ ]


red = RGBColor[1, 0, 0];
blue = RGBColor[0, 0, 1];
green = RGBColor[0, 1, 0];
pink = RGBColor[1, 0, 1];
purple = RGBColor[0.2, 0, 0.6];
yellow = RGBColor[1, 1, 0];