BeginPackage["Caustic`Geodesic`"]

(*  File name: Geodesic.m                                  *)

(*  This package Geodesic.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 and not for release.                    *)

(*  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                        *)

Figure9              ::usage="";

Begin["`Private`"] 

HalfSphere:=ParametricPlot3D[{r Cos[t],r Sin[t],Sqrt[1-r^2]},{r,0,1},{t,0,2Pi},
  Axes->False,AspectRatio->1, Boxed->False, PlotPoints->40];

(* Parametric Plot of geodescis  *)
c={0,0,0,0,1};  P=N[Pi]; p=5;  h=-0.5;
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]}];
Tisch[t_]:={Re[TR[t+P]]/p,Im[TR[t+P]]/p-1,h}; 

P1[alpha_,s_]:={Cos[alpha+P/2] Sin[(P/2)*s],
                Sin[alpha+P/2] Sin[(P/2)*s],
                               Cos[(P/2)*s]}; 
P2[alpha_,s_]:= (Cos[(P/2)*s]^2)*P1[alpha,1] +
                (Sin[(P/2)*s]^2)*Tisch[alpha]; 
P3[alpha_,s_]:= (1-s) Tisch[alpha] + s Tisch[alpha+P]; 
P4[alpha_,s_]:= (Sin[(P/2)*s]^2)*P5[alpha,0] +
               (Cos[(P/2)*s]^2)*Tisch[alpha+P]; 
P5[alpha_,s_]:={Cos[alpha+3 P/2] Sin[(P/2)*(1-s)],
                Sin[alpha+3 P/2] Sin[(P/2)*(1-s)], 
                                 Cos[(P/2)*(1-s)]};           

Geodesic[alpha_,s_]:=If[0  <=Mod[s,1]<1/5,  P1[alpha,Mod[5 s      ,1]],
                     If[1/5<=Mod[s,1]<2/5,  P2[alpha,Mod[5 (s-1/5),1]],
                     If[2/5<=Mod[s,1]<3/5,  P3[alpha,Mod[5 (s-2/5),1]],
                     If[3/5<=Mod[s,1]<4/5,  P4[alpha,Mod[5 (s-3/5),1]],
                     If[4/5<=Mod[s,1]<1,    P5[alpha,Mod[5 (s-4/5),1]]]]]]];  

(* For smoother curves, we average with a Bernoulli distribution  *)
Geodesic1[alpha_,s_]:=(6 Geodesic[alpha,Mod[s,     1]]
                      +4*Geodesic[alpha,Mod[s-0.03,1]]
                      +4*Geodesic[alpha,Mod[s+0.03,1]]
                        +Geodesic[alpha,Mod[s-0.06,1]]
                        +Geodesic[alpha,Mod[s+0.06,1]])/16;

Geo[alpha_]:=ParametricPlot3D[Evaluate[Geodesic[alpha,s]],{s,0,1},
  Boxed->False,
  AspectRatio->1,
  Axes->False,
  DisplayFunction->Identity,
  PlotPoints->200]; 
MyGeo[alpha_]:=Graphics3D[{Thickness[0.0001],
  Line[Table[Geodesic1[alpha,k/300],{k,0,300}]]}, 
  Boxed->False,AspectRatio->1,
  Axes->False,
  ViewPoint->{5,1,1}];
Family[n_]:=Show[Table[Geo[P k/n],{k,0,n-1}],
  DisplayFunction->$DisplayFunction]
MyFamily[n_]:=Show[Table[MyGeo[P k/n],{k,0,n-1}],
  DisplayFunction->$DisplayFunction, 
  ViewPoint->{5,1,1}]

Surface:=ParametricPlot3D[Evaluate[Geodesic[alpha,s]],{s,0,1},{alpha,0,P/2},
  Boxed->False,
  AspectRatio->1,
  Axes->False,
  Shading->False,
  PlotPoints->20,
  ViewPoint->{5,5,-3}];

Figure9:=MyFamily[40];             (* newfigure9.ps (  thinner lines) *)   

End[]             
EndPackage[]

