(* GraphicalAnalysis.m Package,  By Thomas LoFaro, tlofaro@gustavus.edu *)
(* Copyright 2001, Thomas LoFaro. Permission is hereby granted to modify*)
(* and/or make copies of this file for any purpose other than direct    *)
(* profit, or as part of a commercial product, provided this copyright  *)
(* notice is left intact.  Sale, other than for the cost of media, is   *)
(* prohibited.  Permission is hereby granted to reproduce part or all   *)
(* of this file, provided the source is acknowledged.                   *)

BeginPackage["GraphicalAnalysis`"]

ListIterates::usage = "ListIterates[expr, {x,x0}, opts] iterates the 
expression expr with initial condition x0 to generate a list of
iterates.  ListIterates[expr, {x,x0,x1,...,xn}, opts] generates lists
of iterates with initial conditions x0,x1,...,xn.  Options are 
Transients and Iterations."

Options[ListIterates] = {Transients -> 0, Iterations -> 10};
 
Iterate::usage = "Iterate[f, start, color, opts]
iterates the function f with initial conditions start to
generate a stairstep diagram.  The output is a Graphics object 
of color c.  If argument c is omitted then the default color
is GrayLevel[.3].  Options are Transients and Iterations."

Options[Iterate] = {Transients -> 0, Iterations -> 10};

PlotIterates::usage = "PlotIterates[expr, {x,x0,x1}, ics, opts]
plots the graph of expr on the interval {x0,x1} and the orbits
with initial conditions given in the list ics.  Options are
Transients, Iterations, and Composition."

Options[PlotIterates] = {Transients -> 0, Iterations -> 10,
	Composition -> 1};

BifIterate::usage = "BifIterate[expr, ics, pars, opts]
iterates the function of one variable starting at at the
initial condition given in the list ics {x, x0} with
parameter values given in the list pars {p, p0}.
Options include Transients and Iterations."

Options[BifIterate] = {Transients -> 20, Iterations -> 40,
	Tolerance -> 10};
	
Bifurcation::usage = "Bifurcation[expr, ics, pars, opts]
computes a bifurcation diagram for the function expr.  
ics_ is a list of initial conditions of the form {x,x0,...,xn}
and p is a list of parameters of the form {p,p_start,p_end}.
Options include Transients, Iterations and Slices."

Options[Bifurcation] = 
	{Transients -> 20, Iterations -> 40, Slices -> 100,
	 Color -> True, Tolerance -> 10};

Begin["`Private`"]

(* ListIterates generates lists of iterates of expr with  *)
(* initial conditions x0.                                *)
ListIterates[expr_, {x_,x0__}, opts___] := 
  Module[{func,res,trans,its}, 
     trans = Transients /. {opts} /. Options[ListIterates];
     its = Iterations /. {opts} /. Options[ListIterates];
		
		func = Function[x, expr];
     res = Map[Nest[func,#,trans]&,{x0}];
		res = Map[NestList[func, #, its]&, res]
   ]
   
(* define a function to iterate a function f with initial *)
(* condition st                                           *)
Iterate[func_, st_, color_:GrayLevel[0.3], opts___] := 
  Module[{res,trans,its}, 
    	trans = Transients /. {opts} /. Options[Iterate];
    	its = Iterations /. {opts} /. Options[Iterate];
    	res = NestList[func, Nest[func,st,trans], its]; 
    	res = Rest[Flatten[Transpose[{res, res}]]];
	res = Partition[res, 2, 1]; 
    	Graphics[{color, Line[res]}]
	];
	
PlotIterates[expr_, limits_List, ics_List, opts___] :=
	Module[{pf, pics, colors,trans,its,comp,newexpr},
		trans = Transients /. {opts} /. Options[PlotIterates];
    	its = Iterations /. {opts} /. Options[PlotIterates];
    	comp = Composition /. {opts} /. Options[PlotIterates];
    	
    	(* compute the list of colors *)
    	colors = Table[Hue[1 -  i/Length[ics] ],{i,Length[ics]} ];
    	
    	(* generate function to be plotted based on        *)
    	(* Composition option                              *)
    	newexpr = Nest[
    		Evaluate[expr /.limits[[1]] -> #] &,
    			limits[[1]],
    			comp
    	];
    	 
    	(* graph function being iterated *)
    	pf = {Plot[{limits[[1]], newexpr},  limits,
    		AspectRatio -> Automatic,
    		DisplayFunction -> Identity]};
    		
    	(* generate list of iteration pictures *)
     	pics = Inner[
     		Iterate[
     			Evaluate[newexpr /.limits[[1]] -> #] &, #1, #2, opts
     		] &,
            ics, colors, List
        ];
    	(* show it all together *)
		Show[Prepend[pics,pf],
			DisplayFunction -> $DisplayFunction,
			AspectRatio -> Automatic]
	];


(* define function which iterates.  The arguments are the *)
(* expression to iterate and lists of the form            *)
(* {x, x0},{p,p0}.  This function returns a list of       *)
(* the form {{p,x[n]},{p,x[n+1]}, ... {p, x[N]}}.         *)
(*                                                   	  *)
BifIterate[expr_, {x_,x0_}, {p_, p0_}, opts___] :=
Module[{trns, numits, func, xstart, tol},
	trns = Transients /. {opts} /. Options[Bifurcation]; 
	numits = Iterations /. {opts} /. Options[Bifurcation];
	tol - Tolerance /. {opts} /. Options[Bifurcation];
	
	func = Function[x, Evaluate[expr /. p->p0]];
	
		(* get rid of transients *)
        xstart = N[Nest[func, N[x0], trns]];
        (* check to see if tol exceeded *)
        If[Abs[xstart] < 10, 
    		
    		Flatten[
				Outer[List, {p0},
        			NestList[func, xstart, numits - trns]
				],
			1],
		
		(* Plot zeros if tol exceeded *)
			Flatten[
				Outer[List, {p0},
					NestList[0&,0,numits - trns]
				],
			1]
		]
	]

(* This function computes a bifurcation diagram for the    *)
(* function expr.  ics_ is a list of initial conditions    *)
(* of the form {x,x0,...,xn} and p is a list of parameters *)
(* of the form {p,p_start,p_end}.                          *)

Bifurcation[expr_, ics_List,{p_, pstart_, pend_}, opts___] := 
  Module[{trns,numits,slcs,pars,color,CorG,start,pics,tol},
    trns = Transients /. {opts} /. Options[Bifurcation];
    numits = Iterations /. {opts} /. Options[Bifurcation]; 
    slcs = Slices /. {opts} /. Options[Bifurcation];
    CorG = Color /. {opts} /. Options[Bifurcation];
    tol = Tolerance /. {opts} /. Options[Bifurcation];

    (*                                         *)
    (* determine the list of parameter values  *)
    (*                                         *)
    pars = Table[N[pstart + i (pend - pstart)/slcs], {i,0,slcs}];
    (*                                         *)
    (* determine the list of colors indexed by *)
    (* initial conditions.                     *)
    (*                                         *)
    If[CorG == True,
	Table[
		color[i] = Hue[i/Length[ics]],
			{i,Length[ics]-1}
			],
	Table[
		color[i] = GrayLevel[.9 (i-1)/(Length[ics])],
			{i,Length[ics]-1}
			]
	];
    (*                                                *)
    (* Construct list of initial conditions           *)
    (* and associated parameter values of             *)
    (* the form                                       *)
    (* {{{f1[p[[1]]],p[[1]]},...,{f1[p[[n]]],p[[n]]}},*)
    (* ...,{{fk[p[[1]]],p[[1]]}...}}                  *)
    (*                                                *)
    start =Table[
		Inner[List,
	      	      Map[Function[ p, N[ics[[k]]] ], pars],
	      	      pars,
	      	      List
	     	     ],
		{k,2,Length[ics]}
		];
    (*                                         *)
    (* compute bifurcation diagram             *)
    (*                                         *)
    pics = Table[
		ListPlot[ Flatten[
		    Map[ 
		    BifIterate[expr,{ics[[1]],#[[1]]},{p,#[[2]]},opts]&,
		    start[[k]] ],
		    1],
		PlotStyle -> {color[k], AbsolutePointSize[1.5]},
		DisplayFunction -> Identity],
	   	{k, 1, Length[ics]-1} ];
    Show[pics, DisplayFunction -> $DisplayFunction]
    ]
End[ ]

EndPackage[ ]
Protect[Iterate, PlotIterates, BifIterate, Bifurcation];
(*


^*)