BeginPackage["Caustic`LpBill`"]

(*  File name: LpBill.m                                    *)

(*  This package LpBill.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    *)
(*  It is 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                        *)

Figure0     ::usage=""; 

Begin["`Private`"] 

(* Illustration billiards in x^4+y^4=1                                    *)
(* unitball in l^p, Figure in introduction                                *)

<<"Graphics`ImplicitPlot`"

Graylevel=1.0;   (* put it 0.0 if a black background is needed   *)
ort[v_]:={-v[[2]],v[[1]]}; 
Normed[ww_]:=ww/Sqrt[ww[[1]]^2+ww[[2]]^2];
l[R_,S_]:=Sqrt[(R[[1]]-S[[1]])^2+(R[[2]]-S[[2]])^2];
f[P_,p_]:=Abs[P[[1]]]^p+Abs[P[[2]]]^p -1;
Grad[P_,p_] :=Sign[P]*(Abs[P]^(p-1));
T[{P_,v_},p_]:=Module[{R=P,H=P+v,M=P},
  While[f[H,p]<0,H=H+v];c=f[H,p];
  While[Abs[c]>10^(-15),M=(H+R)/2;c=f[M,p];If[c<0,R=M,H=M]];
  w=Normed[ort[Grad[M,p]]];{M,Normed[2 (v.w) w-v]}];
Coord[{P_,v_},p_]:=Module[{w=Normed[ort[Grad[P,p]]]},
  {Mod[N[Arg[P[[1]]+I P[[2]]]],N[2Pi]],w.v}];
PlotOrbit[orb_,p_]:=ListPlot[Table[Coord[orb[[i]],p],{i,Length[orb]}],
  Frame->True,
  Axes->False, 
  DisplayFunction->Identity,
  PlotRange->{{0,N[2Pi]},{-1,1}}];
Track[orb_]:=Show[Graphics[Line[Table[orb[[i,1]],{i,Length[orb]}]]],
  PlotRange->{{-1.3,1.3},{-1.3,1.3}},
  Axes->False,
  AspectRatio->1,
  Background->GrayLevel[Graylevel],
  AspectRatio->1, 
  DisplayFunction->Identity,
  DefaultColor->RGBColor[0,0,0]];
Curve4=x^4+y^4-1;
ShowCurve:=ImplicitPlot[Curve4==0,{x,-1,1},
  PlotRange->{{-1.3,1.3},{-1.3,1.3}},
  Axes->False,AspectRatio->1,
  AspectRatio->1,
  DisplayFunction->Identity,
  Background->GrayLevel[Graylevel], 
  DefaultColor->RGBColor[0,0,0]];  
FindPoint[{x_,y_},p_]:=Module[{P={0.,0.},
  v={Cos[x],Sin[x]}/2,H,R,M,c,w},
  H=P+v; R=P; M=P; While[f[H,p]<0,H=H+v];c=f[H,p];
  While[Abs[c]>10^(-15),M=(H+R)/2;c=f[M,p];If[c<0,R=M,H=M]];
  w=Normed[Grad[H,p]];{H,-{{Cos[y],Sin[y]},{-Sin[y],Cos[y]}}.w}];
Orbit[z_,p_,n_]:=Module[{Q,S={}}, 
  Q=FindPoint[z,p]; 
  Do[Q=T[Q,p]; S=Append[S,Q],{n}];S];
Coordmod[{P_,v_},p_]:=Module[{w=Normed[ort[Grad[P,p]]]},
  {Mod[N[Arg[P[[1]]+I P[[2]]]],N[2Pi/4]],Abs[w.v]}];
PlotOrbitmod[orb_,p_]:=ListPlot[Table[Coordmod[orb[[i]],p],{i,Length[orb]}],
  Frame->True,
  Axes->False,
  DisplayFunction->Identity,
  PlotRange->{{0,N[2Pi/4]},{0,1}}];
Curve[yy_]:=Table[{N[2 Pi] k/800,yy},{k,0,800}];
TT[{x_,y_},p_]:=Coord[T[FindPoint[{x,y},p],p],p];
ImageCurve[p_,a_,b_]:= Module[{}, TT[{x_,y_}]:=TT[{x,y},p];
  ListPlot[Join[Map[TT,Curve[a]],Map[TT,Curve[b]]]]];
Orb1:=Orbit[{0.4,0.3},4.0,200];
Figure0:=Show[{ShowCurve,Track[Orb1]},
  DisplayFunction->$DisplayFunction]; 

End[]             
EndPackage[]

