BeginPackage["Caustic`Dynamic`"]

(*  File name: Dynamic.m                                   *)

(*  This package Dynamic.m was used for the                *)
(*  illustrations for my paper  O. Knill:                  *)
(*  "On nonconvex caustics of convex billiards"            *)
(*  which appeared in "Elemente der Mathematik, 1998       *)
(*  The package is not tuned for efficiency and clarity    *)
(*  and not documented.                                    *)

(*  The first Version was developed at Caltech. The        *)
(*  second Version at the university of Arizona            *)

(*  Oliver Knill                                           *)
(*  Department of Mathematics                              *)
(*  University of Texas, Austin, TX                        *)
(*  e-mail:  knill@math.utexas.edu                         *)
(*  http://www.ma.utexas.edu/~knill                        *)



Figure1              ::usage="";
Figure3a             ::usage="";
Figure3b             ::usage="";
Figure3c             ::usage="";
Figure3d             ::usage="";
Figure3aT            ::usage="";
Figure3bT            ::usage="";
Figure3cT            ::usage="";
Figure3dT            ::usage="";
Figure13             ::usage="";

Begin["`Private`"] 

f                    ::usage="";
T                    ::usage=""; 
Coord                ::usage=""; 
Orbit                ::usage="";
DisplayOrb           ::usage="";
Orb                  ::usage="";
TR                   ::usage="";
Tisch                ::usage="";
Tisch3D              ::usage="";
ShowTisch            ::usage="";
SchowTisch3D         ::usage="";
Trajectory           ::usage="";
Trajectory3D         ::usage="";
ShowTraj             ::usage="";
ShowTraj3D           ::usage="";
TischAndTraj         ::usage="";
MultiTisch           ::usage="";
ThreeDTisch          ::usage="";
TischAndTraj         ::usage="";
TischAndTraj3D       ::usage="";

P=N[Pi]; p=6; 
f[t_,y_,x_]:=-Cos[3t- y]/4 - p*Cos[y] + Cos[3*t + y]/2 + 
   Cos[t - 4*x + y]/4 + p*Cos[t - x + y] - Cos[t + 2*x + y]/2; 
T[{t_,y_}]:=Module[{},
  t1=FindRoot[f[t,y,x]==0,{x,t+2y}][[1,2]];
  While[(t1-t)^2<10^(-10),t1=FindRoot[f[t,y,x]==0,{x,t+Random[]}][[1,2]]];
  {Mod[t1,2P],Mod[t1-t-y,P]}];
Coord[{a_,b_}]:={a,Cos[b]};  
Orbit[b_,n_]:=Module[{u={Random[],N[ArcCos[b]]}},
  Do[u=T[u]; (* Print[j]; *)
   s=Append[s,Coord[u]],{j,n}];s]; 
Orbit[a_,b_,n_]:=Module[{u={a,N[ArcCos[b]]}},Do[u=T[u];(* Print[j]; *)
  s=Append[s,Coord[u]],{j,n}];s];
DisplayOrb[a_,b_,n_]:=Module[{}, s={}; 
  OO=ListPlot[Orbit[a,b,n],DisplayFunction->Identity];
  Show[{Graphics[Line[{{0,0},{2P,0}}]],OO},
    Axes->False,
    Frame->True,
    FrameTicks->None,
    DisplayFunction->$DisplayFunction,
    PlotRange->{{0,2P},{-1,1}}]];
Orb[nn_,mm_]:=Module[{},s={};
  Do[Orbit[N[1-2*kk/mm],nn];
    (* Print["orbit number ",kk]  *),
  {kk,mm-1}];
  OO=Show[Graphics[{PointSize[0.00001],Map[Point,s]}],
    DisplayFunction->Identity];
  Show[{Graphics[Line[{{0,0},{2P,0}}]],OO},
    Axes->False,
    Frame->True,
    DisplayFunction->$DisplayFunction,
    PlotRange->{{0,2P},{-1,1}}]]; 

TR[t_]:= -I/2*(1 - E^(-2*I*t)) - I/4*(-1 + E^(4*I*t)) - I*(-1 + E^(I*t))*p;
Tisch[t_]:={Re[TR[t]],Im[TR[t]]};
Tisch3D[t_]:={Re[TR[t]],0,Im[TR[t]]};
Tisch3D[t_,a_]:={{ Cos[a],Sin[a],0},
                 {-Sin[a],Cos[a],0},
                 {      0,     0,1}}.{Re[TR[t]],0,Im[TR[t]]};
ShowTisch:=ParametricPlot[Tisch[t],{t,0,2P},
  Axes->False,
  AspectRatio->1,
  DisplayFunction->Identity]; 
ShowTisch3D:=ParametricPlot3D[Tisch3D[t],{t,0,2P},
  Axes->False,
  AspectRatio->1,
  DisplayFunction->Identity]; 
Trajectory[O_]:={Thickness[0.0001],
  Line[Map[Tisch,Table[O[[i,1]],{i,Length[O]}]]]}; 
Trajectory3D[O_]:={Thickness[0.0001],
  Line[Map[Tisch3D,Table[O[[i,1]],{i,Length[O]}]]]}; 
ShowTraj[b_,n_]:=Module[{},s={}; OOO=Orbit[b,n]; 
  Show[Graphics[Trajectory[OOO]],
  DisplayFunction->Identity]]; 
ShowTraj[a_,b_,n_]:=Module[{},s={}; OOO=Orbit[a,b,n];
  Show[Graphics[Trajectory[OOO]],
  DisplayFunction->Identity]];
ShowTraj3D[a_,b_,n_]:=Module[{},s={}; OOO=Orbit[a,b,n];
  Show[Graphics3D[Trajectory3D[OOO]],
  DisplayFunction->Identity]];
TischAndTraj[a_,n_]:=Show[{ShowTraj[a,n],ShowTisch},
  Axes->False,
  DisplayFunction->$DisplayFunction,
  AspectRatio->1]; 
TischAndTraj[a_,b_,n_]:=Show[{ShowTraj[a,b,n],ShowTisch},
  Axes->False,
  DisplayFunction->$DisplayFunction,
  AspectRatio->1];
MultiTisch[m_]:=ParametricPlot3D[Tisch3D[t,a],{t,Pi,2Pi},{a,0,Pi},
  PlotPoints->m,
  DisplayFunction->Identity];
ThreeDTisch[m_]:=ParametricPlot3D[Tisch3D[t,a],{t,Pi,2Pi},{a,0,2Pi},
  PlotPoints->m,
  Boxed->False];
TischAndTraj3D[a_,b_,n_]:=Show[{ShowTraj3D[a,b,n],ShowTisch3D,MultiTisch[20]},
  Boxed->False,Axes->False,DisplayFunction->$DisplayFunction,AspectRatio->1,
  Shading->False (* ,ViewPoint->{6,-6,7} *)];

Figure1:=Orb[300,40]; 

Figure3a:=TischAndTraj[1.0,0.470,200];     
Figure3b:=TischAndTraj[1.0,0.060,200];      
Figure3c:=TischAndTraj[1.0,0.939,200];      
Figure3d:=TischAndTraj[1.8,0.500,200];     

Figure3aT:=DisplayOrb[1.0,0.470,200];       
Figure3bT:=DisplayOrb[1.0,0.060,200];
Figure3cT:=DisplayOrb[1.0,0.939,200];
Figure3dT:=DisplayOrb[1.8,0.500,200];

Figure13:=TischAndTraj3D[1.8,0.500,200];   

(* Derivation and simplifications of the             *)
(* intersection Condition  (Appendix of Article)     *)

TR[t_]:= -I/2*(1 - E^(-2*I*t)) - I/4*(-1 + E^(4*I*t)) - I*(-1 + E^(I*t))*p;
Simplify[(TR[t] - TR[x])*Exp[-I (t+y)]/I];
f1[t_,y_,x_]:=Re[-E^(3*I*t - I*y)/4 - p*E^(-I*y) + E^(-I*(3*t + y))/2 +
   E^(-I*(t - 4*x + y))/4 + p*E^(-I*(t - x + y)) - E^(-I*(t + 2*x + y))/2];

End[]             
EndPackage[]

