(*:Mathematica Version: 3.0 *)

(*:Package Version: 2.0 *)

(*:Name: DDS1` *)

(*:Title: Exploring Discrete Dynamical Systems of One Variable *) 

(*:Authors: Silvia Heubach and Jen Chen*)

(*:Keywords:  Iterative model equation, equilibrium, cobb-web graph *)


(*:Summary:
This package can be used in conjunction with the palette
         DDS1P. It works for discrete dynamical systems of the form
	x(n+1) = f(x(n)). Functions compute system values from this
	equation, and the model values can be displayed as regular
	graph or as cobb-web graph (statically and dynamically)
*)



BeginPackage["DDS1`","Graphics`MultipleListPlot`"] 

UnProtect[IteratedValSeq,ListGraph,Equil,SolveIt,
			MapIt,LiveMap]


IteratedValueSeq::usage = "IteratedValueSeq[f,var,init,nmax,s] gives the values of x(n) for n = 0 to n = nmax, in steps of size s, for the DDS x(n+1)=f(x(n)) where init is the initial value for the DDS and the function f is expressed in terms of the variable var. The step size s is optional; if no value is given, the step size will be 1. To just get the value of x(nmax), use stepsize s = nmax."

ListGraph::usage = "The function ListGraph[list,joint] is used to graph a list
of paired values. The variable joint is optional and has default value False. If joint = True, then the points will be connected."

Equil::usage = "The function Equil[c,d] is used to calculate the equilibrium of the linear DDS x(n+1)=c*x(n)+d."
 
SolveIt::usage = "The function SolveIt[lhs,rhs,var,approx]is used to solve the equation lhs = rhs for the variable var using the built-in functions Solve and FindRoot (with starting value approx). This value is optional; if no value is given, approx = 0."

MapIt::usage = "The function MapIt[{f, x}, xinit, n, {xmin, xmax}, focus, percent] displays the first n iterations of the DDS x(n+1)=f(x(n)), where f is expressed in terms of x. The iteration starts at x = xinit and is displayed for x values  from xmin to xmax. When zooming in, the graph is centered at x = focus. The value of percent determines the amount of zooming: For 0 < percent < 1, we zoom in, for percent > 1, we zoom out."

LiveMap::usage = "The function LiveMap[{f, x}, xinit, n, {xmin, xmax}, {ymin, ymax}] creates the graphics for an animation of n iterations of the DDS x(n+1)=f(x(n)), where f is expressed in terms of x. The iteration starts at x = xinit and is displayed for input values from xmin to xmax and output values from ymin to ymax. " 

(* Focus::usage="The function 'Focus[{func_,var_},{xmin_,xmax_},{ymin_,ymax_},
percent_:2]' is used to zoom into the center of a picture; where the 4
corners of the picture are generated by the values of {xmin,xmax}, and
{ymin,ymax}. The value of 'percent' must be greater than 1, with percent=2
means normal size. If 1 < percent <= 2, then we will zoom into the center
of the picture. If the value of percent > 2, then we will zoom out of the
picture. Note that 'func' is the user's input function; 'var' is the variable
of the function; 'xmin' and 'xmax' are the values on the two ends of the
a-axis; 'ymin' and 'ymax' are the values on the two ends of the y-axis; and
'percent' is the zoom value mentioned before." *)



Begin["`Private`"] 

    Off[RuleDelayed::rhs];
    Off[SetDelayed::write];


 IteratedValueSeq[f_,var_,init_,n_,step_:1]:=
 		Module[{func,val,seq},
				func[x_]:= f /.{var ->x};
				vals = NestList[func,init,n];
				seq = Table[{i-1,vals[[i]]},{i,1,n+1,step}];
				Prepend[seq,{"n","x(n)"}]//TableForm]
		(* End of Module *)	
 
		

 ListGraph[name_,joint_:False]:=Module[
			{ls, xlabel=" ",ylabel=" "},
	If[ListQ[name]==False,
			ls = name[[1]];
			If[NumberQ[N[ls[[1,1]]]]==False,
	           xlabel = ls[[1,1]];
	           ylabel = ls[[1,2]];
	           ls = Rest[ls]];
	            ListPlot[ls, PlotStyle->{RGBColor[1,0,0],PointSize[0.02]},
				PlotJoined->joint,	AxesLabel ->{xlabel,ylabel}],
	(*else*)
				ls = name;
Print[ls];
			   If[NumberQ[N[ls[[1,1]]]]==False,
	              xlabel = ls[[1,1]];
	               ylabel = ls[[1,2]];
	              ls = Rest[ls]];
			   ListPlot[ls,PlotStyle->{RGBColor[1,0,0],PointSize[0.02]},
				PlotJoined->joint]]
	]     (* End of Module  *)
 

Equil[x_,y_]:=
    Module[{equilibrium},
          equilibrium = y/(1 - x)
          ]
          
          
SolveIt[equa1_,equa2_:0,var_,approx_:0]:=
	Module[{equation},
     Off[Solve::tdep];
     Off[Out::intm];
		Off[FindRoot::jsing];
         
If[ListQ[Solve[equa1-equa2==0,var]]==True,
Solve[equa1-equa2==0,var],
{
Print["The solution of your equation can not be found by
using the built-in function Solve[]. We will try to solve it
by using the function FindRoot[]. This function will display
one solution (if it found one) at a time."];

      FindRoot[equa1==equa2,{var,approx}]
  }
   ]   (* End of IF statement     *)

]   (* End of MODULE  *)	


MapIt[{funct_,var_},xinit_,numofiteration_:10,{xmin_,xmax_},focus_:0,
    percent_:1]:=
   Module[{easyline,Newfunc,xdistance,newxmin,newxmax},
		
xdistance=Abs[xmax-xmin];
		If[focus==0,
			 newxmin=xmin;
			  newxmax=xmax   ,
		newxmin=focus-xdistance*percent/2 ;
		  newxmax=focus+xdistance*percent/2]  ;     (* End of IF *)
		
      Newfunc[var_]:=funct;
	  easyline[w_]:={Line[{{w,w},{w,Newfunc[w]}}],
		 Line[{{w,Newfunc[w]},{Newfunc[w],Newfunc[w]}}]} ;
				
	Plot[{funct,var},{var,newxmin,newxmax},
		 Epilog->{
			{PointSize[0.02],Point[{xinit,xinit}]},
			 {RGBColor[1,0,0],	
				Map[easyline,NestList[Newfunc,xinit,numofiteration-1]]}
				}    (*  End of EPILOG  *)
	]            (*  End of PLOT *)
	]            (* End of Module *)



LiveMap[{func_,var_},xinit_,
    numofiteration_,{xmin_,xmax_},{ymin_,ymax_}]:=
      Module[{allines,Newfunc},

        Newfunc[var_]:= func;
		allines[w_]:={Line[{{w,w},{w,Newfunc[w]}}],
		      Line[{{w,Newfunc[w]},{Newfunc[w],Newfunc[w]}}]} ;
		
		 Plot[{func,var},{var,xmin,xmax},
				         PlotRange->{ymin,ymax},
				             Epilog->{
					             {PointSize[0.02],Point[{xinit,xinit}]}
				}
			];
		
		 Plot[{func,var},{var,xmin,xmax},
		        PlotRange->{ymin,ymax},
		            Epilog->{
		             {PointSize[0.02],Point[{xinit,xinit}]},
			        	{
						RGBColor[0,0,1],
          Take[Flatten[Map[allines,NestList[Newfunc,xinit,0]]],{1}]
					}
					}
	]   ;    (* End of PLOT  *)
		
		For[k=0,k < numofiteration,k++,
		     Plot[{func,var},{var,xmin,xmax},
		         PlotRange->{ymin,ymax},
		             Epilog->{
		             {PointSize[0.02],Point[{xinit,xinit}]},
			        	{RGBColor[0,0,1],
		                 Map[allines,NestList[Newfunc,xinit,k]]}
						}
	]       (* End of PLOT  *)
		]     (* End of FOR *)
		]       (* End of MODULE *)

	
	
 Focus[{func_,var_},{xmin_,xmax_},{ymin_,ymax_},percent_:2]:=
	Module[{Bigger,xzoom,yzoom,newxmin,newxmax,
			newymin,newymax},
	If[percent<=1,Print["Please enter a value greater than 1. A 
	value of 2 is equivalent to the default value for display"],
		  xzoom= (Abs[xmax-xmin]*percent)/2;
	        yzoom= (Abs[ymax-ymin]*percent)/2;
      newxmin=xmin+xzoom;
		newxmax=xmax-xzoom;
		  newymin=ymin+yzoom;
		     newymax=ymax-yzoom;
		Bigger=Plot[func,{var,newxmin,newxmax},PlotRange->
					{newymin,newymax}]
		]     (* End of IF  *)
	]
		
	
End[]
Protect[IteratedValue,ValueSeq,ListGraph,Equil,SolveIt,
			MapIt,LiveMap]

EndPackage[]