Inner Billiards

This notebook contains the basic Outer Billiards code along with some inner billiards code from the Billiards.nb notebook by Eric Weisstein available at  http://mathworld.wolfram.com/Billiards.html. (See Inner Billiards.pdf.) This code is fairly  primitive and only covers ellipses and regular polygons. We have extracted the portions that generate orbits. The Outer Billiards routines can be used to simulate the dual billiard map in the ‘ smooth’ case of a circle or ellipse.

Outer Billiards Routines

npoints=200;dplaces=30;HalfN=Floor[npoints/2];

Mc=Reverse[Table[{ 2*Cos[t],Sin[t]},{t,0,2 Pi,2 Pi/npoints}]];

Mc=Append[Mc,Mc[[1]]];

TangentMapRoutines=Module[{},Axy=Range[npoints];Mxy=Axy;Bxy=Axy;
Edges=Table[{Mc[[i]],Mc[[i+1]]},{i,1,npoints}];

Dxy=Table[{Edges[[i,2,1]]-Edges[[i,1,1]],Edges[[i,2,2]]-Edges[[i,1,2]]},{i,npoints}];(*Deltax and Deltay for each edge*)

Do[If[Abs[Dxy[[i,2]]]≤Abs[Dxy[[i,1]]],Axy[[i]]="True",Axy[[i]]="False"],{i,1,npoints}];(*we will swap x and y when slope>1*)

Do[If[Axy[[i]]=="True",Mxy[[i]]=Dxy[[i,2]]/Dxy[[i,1]];
If[Abs[Mxy[[i]]]<10^(-10),Mxy[[i]]=0],Mxy[[i]]=Dxy[[i,1]]/Dxy[[i,2]]];
If[Abs[Mxy[[i]]]<10^(-10),Mxy[[i]]=0],{i,1,npoints}]; (*the slopes of the edges*)

Do[If[Axy[[i]]=="True",Bxy[[i]]=Mc[[i,2]]-Mxy[[i]]*Mc[[i,1]],Bxy[[i]]=Mc[[i,1]]-Mxy[[i]]*Mc[[i,2]]],{i,1,npoints}];
(*intercepts*)
Do[f[i_,x_]:=N[Mxy[[i]],dplaces]*x+N[Bxy[[i]],dplaces],{i,1,npoints}];(*one function for each edge*)

Fc[p_]:=Module[{d=Range[npoints]},Do[d[[i]]=(p[[1]]-Mc[[i,1]])^2+(p[[2]]-Mc[[i,2]])^2,{i,1,npoints}];
Extract[Flatten[Position[d,Min[d]]],1]];(*Use nearest vertex for first corner*)

Cros[p_,i_]:=Dxy[[i,2]]*(Mc[[i,1]]-p[[1]])-Dxy[[i,1]]*(Mc[[i,2]]-p[[2]]);(*cross procuct to test for valid tangent*)

Nc[p_]:=Module[{},While[Cros[p,cc]>0,cc=Mod[cc+1,npoints,1]];
If[Cros[p,cc]<-10^(-14),Return[cc],Return[0]]];(*Use the'Right-hand' rule to determine if the tangent line is valid*)

Np[p_,i_]:=If[i>0,Flatten[{{-p+2*Mc[[i]]}}],p];(*next point routine:pnew=-p+2c*)

NextPoint[p_]:=N[Np[p,Nc[p]],dplaces];(*apply NP to corner to get next point in orbit*)

V[p_,k_]:=Module[{},cc=Fc[p];
NestList[NextPoint,p,k]]; (*End TangentMapRoutines*)];

Tau[p_]:= Last[V[p,1]]; (*The Tangent Map. For orbits use NestList with Tau but V above is faster*)
Tau2[p_]:= Last[V[p,2]];  (*Tau^2 is called the 'Return Map'*)

StepSequenceRoutines=Module[{},Ncp[p_]:=Module[{},For[i=cc,i≤npoints+1,i++,i=Mod[i,npoints,1];(*redo Nc to track skip sequences*)If[Cros[p,i]<-10^(-14),SD[[index]]=Mod[i-cc,npoints];
index=index+1;cc=i;Return[cc](*otherwise increment i*)]]];

Gp[p_]:=Np[p,Ncp[p]];(*For Vp,use old Np and Fc with Ncp*)(*redo Gc to track skip sequences*)

Vp[p_,k_]:=Module[{K},index=1;SD=Range[k+1];cc=Fc[p];
K=NestList[Gp,p,k+1];Lastq=Last[K];SD=Drop[SD,1]];
];

Cr[L_,xmin_,xmax_,ymin_,ymax_]:=Module[{smax=Length[L],SD=Range[Length[L]]},For[s=1,s≤smax,s++,SD[[s]]=s;
If[L[[s,1]]>xmax||L[[s,2]]>ymax,SD[[s]]=0];
If[L[[s,1]]<xmin||L[[s,2]]<ymin,SD[[s]]=0]];
Return[SD]];  (*used in Crop below *)

Crop[L_,xmin_,xmax_,ymin_,ymax_]:=Module[{},SD=Cr[L,xmin,xmax,ymin,ymax];
Ix=DeleteCases[SD,0];Return[L[[Ix]]];]; (*crops points to whatever rectangular region you specify*)

CropBox[A_]:=Crop[A,left,right,bottom,top]; (* crops to current values of left,right, bottom, top*)

RotateCorner[w_,nc_,center_]:=Module[{q=Range[nc]},q[[1]]=w;For[s=2,s≤nc,s++,q[[s]]=RotationTransform[(1-s)*2*Pi/nc,center][w]];Return[N[q,dplaces]]];  (*Note the use of RotationTransform here which will not compile for versions before 6.0, but it is easy to rewrite this without using RotationTransform. Since most Mathematica code is optimized it is usually better to use their routines when possible.*)

Period[L_]:=Module[{n=Length[L],f=First[L]},For[i=2,i≤n,i++,If[Abs[L[[i,1]]-f[[1]]]<10^(-14)&&Abs[L[[i,2]]-f[[2]]]<10^(-14),Return[i-1]]]];
(* Period[K] will search the elements of K looking for a match for K[[1]]. You might occasionally need to adjust the match criteria.*)

Periodq[L_,q_]:=Module[{n=Length[L]},For[i=2,i≤n,i++,If[Abs[L[[i,1]]-q[[1]]]<10^(-13)&&Abs[L[[i,2]]-q[[2]]]<10^(-13),Return[i-1]]]];
(*this is used in FindPeriod below for extended searches, using fixed initial point q *)

FindPeriod[startpoint_,steps_,length_]:=Module[{q=startpoint,post=0},q1=q;Unprotect[K];
For[j=1,j≤steps,j++,K=V[q1,length];
post=Periodq[K,q];Print[j,"   ",post];q1=Last[K];
K= .;If[post>1,Abort[]]]];

Datumcrop[startpoint_,steps_,length_,name_]:=Module[{q=startpoint},Unprotect[K];w=Range[steps];Wcrop={};
For[j=1,j≤steps,j++,K=V[q,length];
w[[j]]=Crop[K,left,right,bottom,top];
Print[j,"  "Length[w[[j]]]];
Wcrop=Join[Wcrop,w[[j]]];q=Last[K];K= .];Print[q];
Print[Length[Wcrop]];Save[name<>ToString[j-1]<>"Wcrop",Wcrop]];

web[stepsize_,length_,depth__,offset_]:=Module[{jmax=2*(length-offset)/stepsize,m,z=offset/stepsize},Jxy= Range[npoints];For[s=1,s≤npoints,s++,Jxy[[s]]={};
If[s>1,m=s-1,m=npoints];h[x_]:=N[f[m,x]];
If[Axy[[m]]=="True",x=Mc[[s,1]];
If[Dxy[[m,1]]>0,dx=stepsize,dx=-stepsize];
x=x+z*dx;(*to get sign right*)For[j=1,j<jmax,j=j+2,x=x+dx;
Jxy[[s]]=Join[Jxy[[s]],V[{x,h[x]},depth]]],y=Mc[[s,2]];If[Dxy[[m,2]]>0,dy=stepsize,dy=-stepsize];
y=y+z*dy;
For[j=1,j<jmax,j=j=j+2,y=y+dy;
Jxy[[s]]=Join[Jxy[[s]],V[{h[y],y},depth]];
(*If*)]   (*For*)]  (*For*)]  (*web*)];

WebPoints[stepsize_,length_,depth__]:=Module[{},web[stepsize,length,depth,0];WP={};
For[i=1,i≤npoints,i++,WP=Union[WP,Jxy[[i]]]];
Return[WP]];

Gr[depth_]:=Show[Graphics[{AbsolutePointSize[1.0],poly[Mom],Blue,Point[WebPoints[.01,6,depth]]},PlotRange->{{left,right},{bottom,top}}]];  (*Set box[center,length] to your region of interest and change Gr to match your needs*)

(*GraphicsGrid[{{Gr[0],Gr[1],Gr[2],Gr[3]},{Gr[4],Gr[5],Gr[6],Gr[7]}},Frame->All]*) (*four by 4 grid to different depths*)

box[center_,r_]:=Module[{},left=center[[1]]-r;right=center[[1]]+r;
top=center[[2]]+r;bottom=center[[2]]-r];
thinbox[center_,r_]:=Module[{},left=center[[1]]-r;right=center[[1]]+r;
top=center[[2]]+r/3;bottom=center[[2]]-r/3];
Mbox[Cx_]:=Module[{},left=Cx[[1]][[1]][[1]];
right=Cx[[1]][[2]][[1]];top=Cx[[1]][[2]][[2]];
bottom=Cx[[1]][[1]][[2]];];
(* box[{0,0},3] will set left,right,bottom,top to match a square of side 6 centered at {0,0} so you can say:
Show[Grahics[......, PlotRange->{{left,right},{bottom,top}}]], Mbox is for pasting (Control-C) a Mathematica-generated box to get a generic box with left,right,botom top instead of two corners*)

RadiusFromHeight[height_,n_]:=N[height/Cos[Pi/n],dplaces];
SideFromRadius[radius_,n_]:=N[2*radius*Sin[Pi/n],dplaces];
RadiusFromSide[side_,n_]:=N[side/(2*Sin[Pi/n]),dplaces];
HeightFromSide[side_,n_]:=N[side/(2*Tan[Pi/n]),dplaces];
HeightFromRadius[radius_,n_]:=N[radius*Cos[Pi/n],dplaces];
s0=SideFromRadius[1,npoints];h0=HeightFromRadius[1,npoints];
rDad=RadiusFromSide[s0,2*npoints];hDad=HeightFromRadius[rDad,2*npoints];
EvenScale[n_]:=N[(1-(Cos[Pi/n])^2)/Cos[Pi/n]^2,dplaces];
OddScale[n_]:=N[(1-Cos[Pi/n])/Cos[Pi/n],dplaces];
GenerationScale[n_]:=If[OddQ[n],OddScale[n],EvenScale[n]];

App[W_]:=Append[W,W[[1]]];poly[M_]:=Line[App[M]]; (*Line[M] will not draw a closed polygon*)

Projections = Module[{},Tx = Range[HalfN];Wc=Range[HalfN];
Tx=Table[Table[Mod[j,npoints,1],{j,1,npoints*k,k}],{k,1,HalfN}];
     Wc=Table[Table[Mc[[Tx[[k]][[j]]]],{j,1,npoints}],{k,1,HalfN}];(*remappings of vertices of Mom*)

IND[q_,maxn_]:=Module[{Ind=Range[maxn]},V[q,1];Ind[[1]]=Mod[cc,npoints,1];
St=Vp[q,maxn];
For[s=2,s≤maxn,s++,Ind[[s]]=Mod[Ind[[s-1]]+St[[s-1]],npoints,1]];
Return[Ind]];
PIM[q_,n_,j_]:=Module[{Pim=Range[n]},Pim[[1]]=q;
For[s=2,s≤n,s++,Pim[[s]]=Pim[[s-1]]+2*(Sum[(-1)^k*Wc[[j]][[Ind[[k]]]],{k,2*s-3,2*s-2}])];
Return[Pim]];];

Inner Billiards Routines

lineFromCorners[{p1:{x1_,y1_},p2:{x2_,y2_}}]:=With[{aslope=Divide@@Reverse[p1-p2]},y==y1+aslope (x-x1)];

BilliardPolygon[sides_,start_,angle_,n_]:=Module[{impderiv,corners,cushions,slopes,pos,dir,traj,slope,path,solns,solnrules,hit,eps=10^(-10)},impderiv[{x1_,y1_}]:=slopes[[Position[(First[#]-Last[#])&/@cushions/.Thread[{x,y}->{x1,y1}],_?(Abs[#]<eps&)][[1,1]]]];
corners=Partition[{Cos[#],Sin[#]}&/@Range[0.1,2 π+0.1,2 π/sides],2,1];
slopes=(Divide@@Reverse[#[[1]]-#[[2]]])&/@corners;
cushions=lineFromCorners/@corners;
traj={pos=start};
slope=Tan[dir=angle];
path=(y-pos[[2]]==slope*(x-pos[[1]]));
solnrules=Solve[{path,Or@@cushions},{x,y}];
solns=Select[{x,y}/.solnrules,Plus@@(#^2)≤1.0&];
hit=Select[solns,And@@Thread[#>pos]&];
AppendTo[traj,hit[[1]]];
Do[slope=Tan[2 π-(Apply[ArcTan,hit[[1]]-pos]-ArcTan[impderiv[hit[[1]]]])+ArcTan[impderiv[hit[[1]]]]];
pos=hit[[1]];
path=(y-pos[[2]]==slope*(x-pos[[1]]));
solnrules=Solve[{path,Or@@cushions},{x,y}];
solns=Select[{x,y}/.solnrules,Plus@@(#^2)≤1.0&];
hit=Select[solns,(#-pos).(#-pos)>eps&];
AppendTo[traj,hit[[1]]],{n}];Return[traj]];

Ellipse[{a_,b_},t_]:=N[{a Cos[t],b Sin[t]}];
foci[a_,b_]:={{Sqrt[a^2-b^2],0},{-Sqrt[a^2-b^2],0}}
slope1[{a_,b_},t0_,t1_]:=Divide@@Reverse[Ellipse[{a,b},t1]-Ellipse[{a,b},t0]]/;t0+t1≠2 Pi;
slope2[{a_,b_},t1_]:=-b Cot[t1]/a/;t1≠0&&t1≠Pi;
slope3[{a_,b_},t0_,t1_]:=Tan[2 ArcTan[slope2[{a,b},t1]]-ArcTan[slope1[{a,b},t0,t1]]]/;t1≠0&&t1≠Pi&&t0+t1≠2 Pi;
slope3[{a_,b_},t0_,t1_]:=-slope1[{a,b},t0,t1]/;t1==0||t1==Pi;
slope3[{a_,b_},t0_,t1_]:=Tan[2 ArcTan[slope2[{a,b},t1]]-Pi/2]/;t0+t1==2 Pi;

newpoint[{t0_,t1_},{a_,b_}]:=Module[{x,y,eps=1.*^-10,t2},t2=Mod[ArcTan[b x,a y]/.NSolve[{(x/a)^2+(y/b)^2==1,y-Ellipse[{a,b},t1][[2]]==(x-Ellipse[{a,b},t1][[1]]) slope3[{a,b},t0,t1]},{x,y}],2*Pi];
Select[t2,Abs[#-t1]>eps&][[1]]]

Example: Outer Billiards on a discrete ellipse called Mc with 200 vertices. The initial point is {2,.5}. Change the parameters then Shift-Enter to run again.

q1= {2,.5};K = V[q1,4];Graphics[{poly[Mc],Magenta,AbsolutePointSize[1.0],Arrow[K]}, Axes->True]

Vp[q1,4] (*the step sequence shows how many vertices wre skipped on each iteration*)

IND[q1,4] (*the vertices visited*)

(* Now Inner Billiards inside a regular pentagon*)

Orbit[n_]:=BilliardPolygon[5,{0,0},73Degree,n];Orb[n_]:=Graphics[{{Thickness[0.005],Blue,Line[{Cos[#],Sin[#]}&/@Range[0.1,2 π+0.1,2 π/5]]},Line[Orbit[n]] },AspectRatio->Automatic];
GraphicsGrid[{{Orb[0],Orb[3],Orb[6],Orb[9]},{Orb[12],Orb[15],Orb[20],Orb[100]}},Frame->All]

Graphics[Orb[3]]