BeginPackage["Caustic`ThreeDim`"]

(*  File name: ThreeDim.m                                  *)

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

Figure14a              ::usage="";
Figure14b              ::usage="";
Figure14c              ::usage="";
Figure14d              ::usage="";
Figure14e              ::usage="";
Figure14f              ::usage="";
Figure14g              ::usage="";
Figure14h              ::usage="";

Begin["`Private`"] 

P=N[Pi]; RR:=P*Random[];
p=4; c={0,0,1};
TR[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]}];
 
f[alpha_]:=4 Sin[alpha] +(1/2) Sin[2alpha] +(1/4) Sin[4alpha];
g[alpha_]:= 3.75 - 4 Cos[alpha] +(1/2)  Cos[2alpha] - (1/4) Cos[4alpha];
 
T[alpha_,phi_]:={Cos[phi] f[alpha], -Sin[phi] f[alpha], g[alpha]}; (* table   *)
n[gamma_,delta_]:=(T[gamma+P,delta]-T[gamma,delta]);               (* normal  *)

Le[vv_] := Sqrt[Sum[vv[[i]]^2, {i,3}]];
r[alpha_,phi_]:=Module[{},
  Point1=T[alpha,phi];
  Sqrt[Point1[[1]]^2+Point1[[2]]^2]];  
Diam=N[Le[T[0.2,0]-T[0.2+Pi,0]]]; 
CrossProd[{v1_,v2_,v3_},{w1_,w2_,w3_}]:={v2*w3-v3*w2,
                                         v3*w1-v1*w3, 
                                         v1*w2-v2*w1};
Crossed[{v1_,v2_,v3_},{w1_,w2_,w3_}]:=(v2*w3-v3*w2)^2 +
                                    (v3*w1-v1*w3)^2 + 
                                    (v1*w2-v2*w1)^2;
Clairot[{alpha1_,phi1_,alpha2_,phi2_}]:=Module[{}, 
  zAxes={0,0,1};
  FirstTwo[{a_,b_,c_}]:={a,b,0}; 
  Impact1=T[alpha1,phi1];   Impact2=T[alpha2,phi2]; 
  Incoming=Impact2-Impact1;
  UnitIncoming=Incoming/Le[Incoming]; 
  ToAxes=Impact1-{0,0,Impact1[[3]]};  UnitToAxes=ToAxes/Le[ToAxes]; 
  CosAngle=(FirstTwo[UnitIncoming].UnitToAxes)/Le[FirstTwo[UnitIncoming]]; 
  r[alpha1,phi1]*Sqrt[1-Min[CosAngle^2,1]]]*Le[FirstTwo[UnitIncoming]];
 
ReflectAt[v_,w_]:=(2/Diam^2)*N[(v.w)]*v - w;
 
S[{alpha1_,phi1_,alpha2_,phi2_}]:= Module[{},
  Incoming = N[T[alpha2,phi2]-T[alpha1,phi1]]; 
          UnitIncoming=Incoming/Le[Incoming]; 
  NormalV  = N[n[alpha2,phi2]];              
  Outgoing = N[ReflectAt[NormalV,-Incoming]]; 
          UnitOutgoing=Outgoing/Le[Outgoing]; 
        (* ss={}; ss=Append[ss,Line[{T[alpha1,phi1],T[alpha2,phi2]}]];
          ss=Append[ss,
             {Thickness[0.01],Line[{T[alpha2,phi2],T[alpha2+P,phi2]}]}];
          ss=Append[ss,
             Line[{T[alpha2,phi2], T[alpha2,phi2]+Outgoing}]];
          Show[Graphics3D[ss],Boxed->True];          *)
  {alpha3,phi3}={x,y}/.
  FindMinimum[Crossed[T[alpha2,phi2]-T[x,y],Outgoing],
             {x,2alpha2-alpha1},{y,2phi2-phi1}][[2]];
  While[Le[T[alpha3,phi3]-T[alpha2,phi2]] < 10^(-5), (* Print["."]; *)
  {alpha3,phi3}={x,y} /.
  FindMinimum[Crossed[T[alpha2,phi2]-T[x,y],Outgoing],{x,RR},{y,RR}][[2]]];
  (* Print["test ",(T[alpha3,phi3]-T[alpha2,phi2])/UnitOutgoing];  *)
  {Mod[alpha2,2P],Mod[phi2,2P],Mod[alpha3,2P],Mod[phi3,2P]}];

Orbit[{a_,b_,c_,d_},nn_]:=Module[{u={a,b,c,d}}, s={}; 
  Do[u=S[u]; (* Print[j]; *) s=Append[s,u],{j,nn}];s];
ShowTraj[{a_, b_, c_, d_}, nn_] :=
  Module[{u = {a, b, c, d}}, s = {};
    Do[u = S[u]; (* Print[j, " ", u]; *) 
      s = Append[s, T[u[[1]], u[[2]]]], 
    {j, nn}];
    Show[Graphics3D[{Thickness[0.0001],Line[s]}], 
      Boxed -> False, 
      Axes  -> False, 
      AspectRatio -> 1] ];
 
ShowNormals[{a_,b_,c_,d_},nn_]:=Module[{u={a,b,c,d}},s={};
  Do[u=S[u]; (* Print[j," ",u]; *)
    s=Append[s,Line[{T[u[[1]],u[[2]]],T[u[[1]],u[[2]]]+n[u[[1]],u[[2]]]}]],
  {j,nn}];
  Show[Graphics3D[s], 
    Boxed->False, 
    Axes->False, 
    AspectRatio->1]];
 
ShowTheSurface:=ParametricPlot3D[T[a,b],{a,0,2P},{b,0,P},PlotPoints->20,
   Shading->False,
   RenderAll->True,
   Boxed->False]; 

ClairotTest:=Module[{uuu={RR,RR,RR,RR}},Clairot[S[uuu]]/Clairot[uuu]]; 
Period2Test:=Module[{u1 = RR, u2 = RR},Sin[S[S[{u1,u2,u1+P,u2}]]-{u1,u2,u1,u2}]];

ShowTraj[nn_]:=ShowTraj[{RR,RR,RR,RR},nn]
ShowNormals[nn_]:=ShowNormals[{RR,RR,RR,RR},nn]

Figure14a:=ShowTraj[{1.73513, 2.17327, 1.96578, 3.08649},600];
Figure14b:=ShowTraj[{1.73513, 2.17327, 1.66578, 3.08649},600];
Figure14c:=ShowTraj[{1.44518, 0.488606, 1.54747, 4.29294},200];
Figure14d:=ShowTraj[{0.4,0.4,0.4+P,0.45},200];
Figure14e:=ShowTraj[{0.02,0.1,P+0.09,0.2},600];
Figure14f:=ShowTraj[{P/2+0.02,0.1,P/2+P+0.06,0.3},600];
Figure14g:=ShowTraj[{P/2+0.02,0.1,P/2+P+0.09,0.3},600];
Figure14h:=ShowTraj[{P/2+0.02,0.1,P/2+P+0.12,0.3},600];

End[]             
EndPackage[]

