BeginPackage["Caustic`Caustic`"]

(*  File name: Caustic.m                                   *)

(*  This package Caustic.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                        *)


Figure2a             ::usage=" ";  
Figure2b             ::usage=" ";  
Figure8a             ::usage=" ";  
Figure8b             ::usage=" ";  
Figure4a             ::usage=" ";  
Figure4b             ::usage=" ";  
Figure4c             ::usage=" ";  
Figure5aold          ::usage=" ";
Figure5bold          ::usage=" ";
Figure5a             ::usage=" ";
Figure5b             ::usage=" ";
Figure6a             ::usage=" ";
Figure6b             ::usage=" ";
Figure7              ::usage=" ";
Example0             ::usage=" ";
Example1             ::usage=" ";
Example2             ::usage=" ";
Example3             ::usage=" ";
Example4             ::usage=" ";
Example5             ::usage=" ";
Example6             ::usage=" ";
FC1                  ::usage=" ";
FC2                  ::usage=" ";
CausticMovie         ::usage=" "; 

Begin["`Private`"] 
Caustic            ::usage="Caustic[p,c] draws the caustic of
                    a table of equal width given by the 
                    Fourier expansion p + sum_k 2 c_k cos(k alpha) 
                    of the radius of curvature rho. For example, 
                    if  rho[t_]=p + 2 Cos[3t], then c={0,0,1} ";

T                  ::usage="The image of the function T[t] is a table of 
                            equal width in the complex plane. 
                            The table is given by the Fourier expansion"; 
TablePlot          ::usage="Plots a table of equal width";
rho                ::usage="Radius of curvature of the table given by 
                            a Fourier expansion";
WholeRuledSurface  ::usage="WholeRuledSurface[p0_,M_,c0_] plots the ruled
                            surface over the table of equal thickness. M is
                            the number of lines used to draw the surface";
HalfRuledSurface   ::usage="HalfRuledSurface[p0_,M_,c0_] plots the ruled
                            surface over the table of equal thickness but only 
                            up to the line of striction.  M is
                            the number of lines used to draw the surface";
OldWholeRuledSurface ::usage=" ";
OldHalfRuledSurface  ::usage=" ";  
WaveFront            ::usage=" "; 
WaveFrontMovie       ::usage=" "; 
WeierStrassCaustic   ::usage=" "; 
ZoomIt               ::usage=" "; 
FindRoughCaustic     ::usage=" ";
Cheby                ::usage=" ";
FindTheta            ::usage=" ";
Tisch                ::usage=" "; 
Caustic              ::usage=" "; 
TableAndCaustic      ::usage=" "; 
RuledSurface         ::usage=" "; 
Tisch3D              ::usage=" "; 
U3D                  ::usage=" "; 
CC3D                 ::usage=" "; 
Caustic3D            ::usage=" "; 
Striction            ::usage=" "; 
Proj3D               ::usage=" "; 

Caustic[p0_,c_List]:=Module[{w}, p=p0; 
 w[t_]:= N[I*Sum[(k/(k+1))*c[[k]]   (Exp[I( k+1)t]-1),{k,2,Length[c]}]+
      I*Sum[(-k/(-k+1))*Conjugate[c[[k]]] (Exp[I(-k+1)t]-1),{k,2,Length[c]}]+
      I*(p+2*Sum[c[[k]],{k,2,Length[c]}])];
  ParametricPlot[{Re[w[t]],Im[w[t]]},{t,0,N[2Pi]},
  Axes->False,
  AspectRatio->1,
  PlotPoints->300,
  Background->GrayLevel[1.0],
  PlotRange->All,
  DefaultColor->RGBColor[0.0,1.0,0.0],
  DisplayFunction->Identity,
  PlotStyle->AbsoluteThickness[1]]];

T[t_]:=-I*p*(Exp[I t]-1) -
   I*Sum[          c[[k]] *(1/( k+1))*(Exp[I( k+1)*t]-1),{k,2,Length[c]}] -
   I*Sum[Conjugate[c[[k]]]*(1/(-k+1))*(Exp[I(-k+1)*t]-1),{k,2,Length[c]}];

TablePlot[p0_,c0_List]:=Module[{}, p=p0; c=c0; 
  ParametricPlot[{Re[T[t]],Im[T[t]]},{t,0,N[2Pi]},
  Axes->False,
  AspectRatio->1,
  Background->GrayLevel[1.0],
  DisplayFunction->Identity,
  PlotPoints->300]];

rho[t_]:=p+Sum[ c[[k]] Exp[I k t] + 
         Conjugate[c[[k]]] Exp[-I k t],{k,2,Length[c]}];

OldWholeRuledSurface[p0_,c0_List]:=Module[{},  p=p0; c=c0; 
  ParametricPlot3D[{Re[T[t]+I s*(rho[0]+rho[Pi])*Exp[I t]],
                    Im[T[t]+I s*(rho[0]+rho[Pi])*Exp[I t]],
                    s (rho[0]+rho[Pi])},{t,0,2Pi},{s,0,1},
                    PlotPoints->50,
                    Boxed->False,
                    Axes->False]];
OldHalfRuledSurface[p0_,c0_List]:=Module[{}, p=p0; c=c0;  
  ParametricPlot3D[{Re[T[t]+I s*rho[t]*Exp[I t]],
                    Im[T[t]+I s*rho[t]*Exp[I t]],
                    s*rho[t]},{t,0,2Pi},{s,0,1},
                    PlotPoints->40,
                    Boxed->False,
                    Axes->False]]; 

WholeRuledSurface[p0_,M_,c0_]:=Module[{},p=p0;P=N[Pi]; c=c0;  
  DDiam=rho[0]+rho[P];
  x[t_,s_]:=N[Re[T[t]+I*s*DDiam*Exp[I t]]];
  y[t_,s_]:=N[Im[T[t]+I*s*DDiam*Exp[I t]]];
  z[t_,s_]:=N[s*DDiam];         
  u1=Table[Line[{{x[2P k/M,0],y[2P k/M,0],z[2P k/M,0]},
                 {x[2P k/M,1],y[2P k/M,1],z[2P k/M,1]}}], {k,M}];
  u2=Table[Line[{{x[2P k/M,0],y[2P k/M,0],z[2P k/M,0]},
                 {x[2P (k+1)/M,0],y[2P (k+1)/M,0],z[2P (k+1)/M,0]}}], {k,M}];
  u3=Table[Line[{{x[2P k/M,1],y[2P k/M,1],z[2P k/M,1]},
                 {x[2P (k+1)/M,1],y[2P (k+1)/M,1],z[2P (k+1)/M,1]}}], {k,M}];
  Show[Graphics3D[Join[{u1,u2,u3}]],
  Boxed->False,
  Axes->False ]];

HalfRuledSurface[p0_,M_,c0_]:=Module[{},p=p0;P=N[Pi]; c=c0;  
  DDiam=rho[0]+rho[P]; 
  x[t_,s_]:=N[Re[T[t]+I*s*rho[t]*Exp[I t]]]; 
  y[t_,s_]:=N[Im[T[t]+I*s*rho[t]*Exp[I t]]]; 
  z[t_,s_]:=N[s*rho[t]];  
  u1=Table[Line[{{x[2P k/M,0],y[2P k/M,0],z[2P k/M,0]},
                 {x[2P k/M,1],y[2P k/M,1],z[2P k/M,1]}}], {k,M}];
  u2=Table[Line[{{x[2P k/M,0],y[2P k/M,0],z[2P k/M,0]},
                 {x[2P (k+1)/M,0],y[2P (k+1)/M,0],z[2P (k+1)/M,0]}}], {k,M}];
  u3=Table[Line[{{x[2P k/M,1],y[2P k/M,1],z[2P k/M,1]},
                 {x[2P (k+1)/M,1],y[2P (k+1)/M,1],z[2P (k+1)/M,1]}}], {k,M}];
  Show[Graphics3D[Join[{u1,u2,u3}]],
  Boxed->False,
  Axes->False ]];

WaveFront[p0_,M_,c0_,s0_]:=Module[{},p=p0;P=N[Pi]; c=c0; 
  DDiam=rho[0]+rho[P];MM=2M;
  x[t_,s_]:=N[Re[T[t]+I*s*DDiam*Exp[I t]]];
  y[t_,s_]:=N[Im[T[t]+I*s*DDiam*Exp[I t]]];
  u1=Table[Line[{{x[2P k/MM,s0],    y[2P k/MM,s0]},
                 {x[2P (k+1)/MM,s0],y[2P (k+1)/MM,s0]}}], {k,MM}];
  u2=Table[Line[{{x[2P k/M,s0],y[2P k/M,s0]},
                 {x[2P k/M,s0+0.02],y[2P k/M,s0+0.02]}}],{k,M}];
  u3=Table[Line[{{x[2P k/MM,0],y[2P k/MM,0]},
                 {x[2P (k+1)/MM,0],y[2P (k+1)/MM,0]}}], {k,MM}];
  Show[Graphics[Join[{u1,u2,u3}]],
  Frame->False,Axes->False,AspectRatio->1,PlotRange->{{-p,p},{0,2p}}]];

WaveFrontMovie:=Module[{},
  WF[s1_]:=WaveFront[4,30,{0,0,1},s1]; 
  Show[GraphicsArray[Partition[Table[WF[k/17],{k,1,16}],4]]]]; 

WeierStrassCaustic:=Module[{},
  roughness=0.75;      tt=roughness;
  bb={2^(-tt),0,2^(-tt 2),0,0,0,2^(-tt 3),0,0,0,0,0,0,0,2^(-tt 4),
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2^(-tt 5),0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2^(-tt 6)};
  Show[Caustic[10.0,bb],DisplayFunction->$DisplayFunction]]; 

ZoomIt:=Module[{},
  WC=WeierStrassCaustic; 
  WC0=Show[WC,Frame->True,FrameTicks->None,DisplayFunction->Identity];
  WC1=Show[WC,PlotRange->{{0  ,1.0},{9.0,10.0}},
      Frame->True,FrameTicks->None,DisplayFunction->Identity];
  WC2=Show[WC,PlotRange->{{0.2,0.6},{9.6,10. }},
      Frame->True,FrameTicks->None,DisplayFunction->Identity];
  WC3=Show[WC,PlotRange->{{0.3,0.4},{9.8,9.9 }},
      Frame->True,FrameTicks->None,DisplayFunction->Identity]; 
  Dummy:=Graphics[{},Frame->True,FrameTicks->None]; 
  WCR0=Show[{WC0,Graphics[Rectangle[{0,9.0},{1.0,10.0},Dummy]]},
     DisplayFunction->Identity]; 
  WCR1=Show[{WC1,Graphics[Rectangle[{0.2,9.6},{0.6,10.0},Dummy]]},
     DisplayFunction->Identity]; 
  WCR2=Show[{WC2,Graphics[Rectangle[{0.3,9.8},{0.4,9.9},Dummy]]},
     DisplayFunction->Identity]; 
  WCR3=Show[{WC3,Graphics[Rectangle[{0.3,9.8},{0.4,9.9},Graphics[{}]]]},
     DisplayFunction->Identity]; 
  Show[GraphicsArray[{{WCR0,WCR1},{WCR2,WCR3}}],
    DisplayFunction->$DisplayFunction]           ]; 

FindRoughCaustic[n_]:=Module[{c={0}},
   RR:=2 Random[]-1; 
   Do[c=Append[c,0];c=Append[c,RR/k^2],{k,n}];
   Print[c];Caustic[10.0,c]]; 

Cheby[theta_]:= Tan[n theta] - n Tan[theta];
u:=DSolve[ {x'[t] == (1+p Cos[n t]) Cos[t],x[0]==0}, x[t],t]; (* For formulas *)
v:=DSolve[ {y'[t] == (1+p Cos[n t]) Sin[t],y[0]==0}, y[t],t];
x[t_]:=(-2*Sin[t] + 2*n^2*Sin[t] - p*Sin[t - n*t] - n*p*Sin[t - n*t] -
        p*Sin[t + n*t] + n*p*Sin[t + n*t])/(2*(-1 + n)*(1 + n));
y[t_]:=(-1 + n^2 - p)/(-1 + n^2) +
      (2*Cos[t] - 2*n^2*Cos[t] + p*Cos[t - n*t] + n*p*Cos[t - n*t] +
        p*Cos[t + n*t] - n*p*Cos[t + n*t])/(2*(-1 + n)*(1 + n));
FindTheta[d_]:=theta /. FindRoot[Cheby[theta]==0,{theta,d}];
Tisch:=ParametricPlot[{x[t],y[t]},{t,0,N[2Pi]},
  Axes->False,AspectRatio->1, PlotRange->{{-1.1,1.1},{-0.1,2.1}},
  AspectRatio->1, 
  DisplayFunction->Identity];
b[alpha_]:=(1+p Cos[n alpha]) Sin[theta0];
U[alpha_]:={Cos[alpha+theta0],Sin[alpha+theta0]};
Caustic[t_]:={x[t],y[t]}+ b[t] U[t];
Caus[pp_,nn_]:=Module[{},n=nn;p=pp;
  ParametricPlot[Caustic[t],{t,0,N[2Pi]}, Axes->False,
    AspectRatio->1, 
    PlotRange->{{-1.1,1.1},{-0.1,2.1}},
    DisplayFunction->Identity ]];
TableAndCaustic[p0_,n0_,d_]:=Module[{},n=n0; p=p0; 
  theta0=FindTheta[d];
  Show[{Tisch,Caus[p,n]},  
  DisplayFunction->$DisplayFunction,  
  Background->GrayLevel[0.0], 
  DefaultColor->RGBColor[1.0,0.0,0.0]]];  

RuledSurface[pp_,nn_,d_]:=Module[{}, n=nn; p=pp; theta0=FindTheta[d];
  t1=ParametricPlot3D[{x[t],y[t],0},{t,0,N[2Pi]},
  Axes->False,AspectRatio->1, PlotRange->{{-1.1,1.1},{-0.1,2.1},{-0.1,1.1}},
  DisplayFunction->Identity];
  t2=ParametricPlot3D[{x[t],y[t],1},{t,0,N[2Pi]},
  Axes->False,AspectRatio->1, PlotRange->{{-1.1,1.1},{-0.1,2.1},{-0.1,1.1}},
  DisplayFunction->Identity];
  li[t_]:= {Hue[t/N[2Pi]],
           Line[{{x[t],y[t],0},{x[t+2theta0],y[t+2theta0],1}}]};
  lis=Table[li[N[2Pi]*k/100],{k,1,100}];
  Show[{t1,t2,Graphics3D[lis]}, DisplayFunction->$DisplayFunction,
   Background->GrayLevel[1.0],
   Boxed->False,
   DefaultColor->Hue ]];


l[t_]:=b[t]+b[t+2 theta0]; 
u1:=ParametricPlot3D[{x[t]+s*l[t]*Cos[t+theta0],
  y[t]+s*l[t]*Sin[t+theta0],(s*l[t])/l[t]},
  {t,0,N[2Pi]},{s,0,1},
  PlotPoints->40, 
  Shading->True, 
  Background->GrayLevel[1.0],
  Boxed->False,
  Axes->False]; 
u2:=ParametricPlot3D[{x[t]+b[t]*Cos[t+theta0],y[t]+b[t]*Sin[t+theta0],b[t]/l[t]},
   {t,0,N[2Pi]},
   PlotPoints->60, 
   Background->GrayLevel[1.0],
   Boxed->False, 
   Axes->False]; 

Tisch3D[aa_]:=ParametricPlot3D[{x[t],y[t],aa},{t,0,N[2Pi]},
  Axes->False,AspectRatio->1, PlotRange->{{-1.1,1.1},{-0.1,2.1},{-0.01,2.01}},
  AspectRatio->1,Boxed->False, DisplayFunction->Identity ];
U3D[alpha_]:={Cos[alpha+theta0],Sin[alpha+theta0],0};
Caustic3D[t_]:={x[t],y[t],0}+ b[t] U3D[t];
Caustic3D[pp_,nn_]:=Module[{},n=nn;p=pp;
  ParametricPlot3D[Caustic3D[t],{t,0,N[2Pi]}, Axes->False,
    AspectRatio->1, Boxed->False,
    DisplayFunction->Identity,
    PlotRange->{{-1.1,1.1},{-0.1,2.1},{-0.01,2.01}} ]];
CC3D[pp_,nn_,d_]:=Module[{},
  n=nn; p=pp; theta0=FindTheta[d];
  Show[{Tisch3D[0.0],Tisch3D[1.0],Tisch3D[2.0],Caustic3D[p,n]},
  DisplayFunction->Identity]];
Striction:=
  ParametricPlot3D[{x[t]+b[t]*Cos[t+theta0],y[t]+b[t]*Sin[t+theta0],1+b[t]/l[t]},
   {t,0,N[2Pi]},PlotPoints->60, Background->GrayLevel[1.0],
   PlotRange->{{-1.1,1.1},{-0.1,2.1},{-0.01,2.01}},
   DisplayFunction->Identity,
   Boxed->False, Axes->False];
Proj3D:=CC3D[0.8,5,0.91];                 

(* Examples          *)

Figure2a:=Module[{},p=4; 
   C103:=Caustic[4.0,{0,0,1}];                   
   C105:=Caustic[4.0,{0,0,0,0,1}];              
   C107:=Caustic[4.0,{0,0,0,0,0,0,1}];           
   C109:=Caustic[4.0,{0,0,0,0,0,0,0,0,1}];       
   E103:=Show[{C103,TablePlot[4,{0,0,1}]}];
   E105:=Show[{C105,TablePlot[4,{0,0,0,0,1}]}];
   E107:=Show[{C107,TablePlot[4,{0,0,0,0,0,0,1}]}];
   E109:=Show[{C109,TablePlot[4,{0,0,0,0,0,0,0,0,1}]}]; 
   Show[GraphicsArray[{{E103,E105},{E107,E109}}]]];

Figure2b:=Module[{}, p=10;
   C103:=Caustic[10.0,{0,0,1}];                    
   C105:=Caustic[10.0,{0,0,0,0,1}];                
   C107:=Caustic[10.0,{0,0,0,0,0,0,1}];            
   C109:=Caustic[10.0,{0,0,0,0,0,0,0,0,1}];    
   D103:=Show[{C103,TablePlot[10,{0,0,1}]}];
   D105:=Show[{C105,TablePlot[10,{0,0,0,0,1}]}];
   D107:=Show[{C107,TablePlot[10,{0,0,0,0,0,0,1}]}];
   D109:=Show[{C109,TablePlot[10,{0,0,0,0,0,0,0,0,1}]}]; 
   Show[GraphicsArray[{{D103,D105},{D107,D109}}]]];

Figure8a:=Module[{},
   C3a[g_]:=Caustic[10.0,{0,0,0,0,1-g,0,g}];
   Show[GraphicsArray[Partition[Table[C3a[k/7],{k,1,6}],3]]]];

Figure8b:=Module[{},
    C3b[g_]:=Caustic[10.0,{0,0,0,0,0,0,1-g,0,0,0,0,0,0,0,0,0,g}];
    Show[GraphicsArray[Partition[Table[C3b[k/7],{k,1,6}],3]]]]; 
CausticMovie:=Module[{},
    C3b[g_]:=Caustic[10.0,{0,0,0,0,0,0,1-g,0,0,0,0,0,0,0,0,0,g}];
    Do[Show[C3b[k/11],DisplayFunction->$DisplayFunction],{k,10}]];

Figure4a:=Module[{},                               
  n=5; p=0.8; theta0=FindTheta[0.91]; u1]; 
Figure4b:=TableAndCaustic[0.8,5,0.91];                            
Figure4c:=Show[{Proj3D,Striction},
  ViewPoint->{1.3,-2.4,1.0},
  DisplayFunction->$DisplayFunction];   

Figure5aold:=OldWholeRuledSurface[4,{0,0,1}];                    
Figure5bold:=OldHalfRuledSurface[4,{0,0,1}];                    

Figure5a:=WholeRuledSurface[4,80,{0,0,1}];                    
Figure5b:=HalfRuledSurface[4,80,{0,0,1}];                    

Figure6a:=WeierStrassCaustic;
Figure6b:=ZoomIt; 
Figure7:=WaveFrontMovie; 

Example0:= TableAndCaustic[1,3,0.91];
Example1:= TableAndCaustic[0.3,5,0.91]; 
Example2:= TableAndCaustic[1,5,0.91]; 
Example3:= TableAndCaustic[1,17,0.91]; 
Example4:= TableAndCaustic[1,7,1.1119];
Example5:= TableAndCaustic[0.1,7,0.8];
Example6:=Do[p=Random[];n=Random[Integer,{3,10}]; d=0.4+Random[];
          Print[p,n,d]; TableAndCaustic[p,n,d],{20}];

FC1:=RuledSurface[0.3,5,0.91]; 
FC2:=RuledSurface[1  ,5,0.91]; 

End[] 
EndPackage[]

