
(* --------------------Billiards in R^2----------------------- *)
(*		Ambient data: table = {sides,glue}             *)
(* glue[[i]] =  j:  for 1-form *)
(* glue[[i]] = -j:  for quad diff'l *)
(* glue[[i]] =  i:  for reflector *)
(*		Operations					*)
(* DiTable[] - draws table					*)
(* Shots[s] - draws trajectories at slope s 			*)
(* Moduli[s] - gives ratios of moduli of annuli			*)
(* DiGp[] - draws GL(X,w) assuming reflection (x,-y) is present *)
(* DiL[d] - draws GL(X,w) for symmetric L-table 1:(1+sqrt[d])/2 *)
(* DiS, DiX --- similar for genus 3,4 *)

(* Bounces[slope,i,p,table,N] 
	returns the billiard path of given slope, starting at vertex i (when p=0),
	inside polygon table, of length N.
	When p>0 the path starts internal to the edge with initial vertex i,
	at positive p in [0,1] from that vertex. 
*)

(*   Shot[slope,[i,p,N]] displays same *)
(*   Shots[slope[,N]]     displays paths from all vertices *)

(*---------------------Verbosity ----------------------*) 
(*  2 or more draws pictures; 1 prints to terminal *)
Verbosity=2

(*---------------------Numerical accuracy ----------------------*) 
Digits = 200
Eps := 10^(-12)
Off[N::meprec]
Off[$MaxExtraPrecision::meprec]
RationalQ[x_]:=TrueQ[(Head[x]==Integer || Head[x]==Rational)]

(*---------------------Graphics directives ----------------------*) 
(* Thickness of lines for table and for shot *)

ShotThickness=Thickness[1/600] 
TableThickness=Thick

(*---------------------Lines and Points----------------------*) 
(* An oriented line in R^2 is described by {p,v} 
   p lies on L, v points along L *)

(* Unit vector of given slope *)
UnitVec[s_]:=If[s==Infinity,{0,1},{1,s}/Sqrt[1+s^2]]

(* Rotate by 90 degrees *)
RotV[{x_,y_}]:={-y,x}

(* Slope of a vector *)
Slope[{0,y_}]:=Infinity
Slope[{x_,y_}]:=y/x

(* Reflect point p through a vector v *)
RefVec[p_,v_]:= 2*(v.p)*v/(v.v) - p

(* Compare using first coordinates *)
SortComp[x_,y_]:=x[[1]] < y[[1]]

(* Sort points along a line *)
SortPs[ps_,{p_,v_}]:=Transpose[
	Sort[Transpose[{ps.v,ps}],SortComp]
			][[2]]

(* Intersection of 2 lines *)
(* Must have positive intersection number! *)
MeetL[{p_,v_},{q_,w_}]:=Block[{m,r},
		m={RotV[v],RotV[w]};
		If[Det[m] <= Eps,Return[{}]];
		r = Inverse[m].{m[[1]].p,m[[2]].q};
		r = FullSimplify[r];
		Return[r]
	]

(* Determine if p lies between a and b *)
(* given that p is on the line from a to b *)
Betwixt[p_,{a_,b_}]:= -Eps < (p-a).(b-a) < (b-a).(b-a)+Eps


(* Intersection of a line and an edge *)
MeetE[{p_,v_},{a_,b_}]:=Block[{q},
		q=MeetL[{p,v},{a,b-a}];
		If[q=={},Return[{}]];
		If[!Betwixt[q,{a,b}],Return[{}]];
		Return[q]
]

Vertices[{sides_,glue_}]:=Transpose[sides,{2,1}][[1]];

(* Check if point is at a vertex *)
AtVertex[p_,table_]:= (Min[Map[(p-#).(p-#)&,Vertices[table]]] < Eps)

(* Check if 2 (pts, vectors) coincide *)
AtSame[{p_,v_},{q_,w_}]:= Max[(p-q).(p-q),(v-w).(v-w)] < Eps

(* Return side(s) on which p lies for given vector *)
AtSide[p_,table_]:= Select[Range[Length[table[[1]]]],
			AtSide[p,table,#]&]

AtSide[p_,table_,i_]:=Block[{a,b}, 
			{a,b} = table[[1,i]];
			If[Abs[Det[{a-b,p-b}]] > Eps, Return[False]];
			Return[Betwixt[p,{a,b}]]
]

(* Check if vector points inward from a vertex *)
Inward[v_,i_,table_]:=Block[{det,inner,j,n,v0,v1},
	n = Length[table[[1]]];
	j = i-1;
	If[j < 1, j=n];
	v0 = table[[1,j,2]]-table[[1,j,1]];
	v1 = table[[1,i,2]]-table[[1,i,1]];
	det = Det[{v0,v1}];
	inner = v0.v1;
	If[det == 0, 
		Return[(inner<0) || (Det[{v0,v}] > 0 && Det[{v1,v}] > 0)]];
	Return[
	If[Det[{v0,v1}] > 0, 
		Det[{v0,v}] > 0 && Det[{v1,v}] > 0,
		Det[{v0,v}] > 0 || Det[{v1,v}] > 0]
	]
]

(*---------------------Billiards----------------------*) 
(* Next intersection with polygon *)
(* Return {pt,iside} or {} *)
NextPI[{p_,v_},table_]:=Block[{ds,es,i,is,qs,sides},
		sides=table[[1]];
		is=Range[Length[sides]];
		qs=Map[MeetE[{p,v},#]&,sides];
		qs=Transpose[{qs,is}];
		qs=Select[qs,#[[1]] != {}&];
		If[qs == {}, Return[{}]];
		ds=Map[v.(#[[1]]-p)&,qs];
		ds=FullSimplify[ds];
		es=Select[ds,#>Eps&];
		If[es=={},Return[{}]];
		i=Position[ds,Min[es],1][[1,1]]; 
		Return[qs[[i]]]
]

(* Using gluing data, get new starting point and new vector *)
(* Return {p,w} *)
Glue[{q_,v_},i_,{sides_,glue_}]:= Block[{a,b,c,d,j,p,w},
		{p,w}={q,v};
		j = glue[[i]];
		{a,b}=sides[[i]];
		{c,d}=sides[[Abs[j]]];
		If[i==j, w = RefVec[w,a-b]];
		If[i!=j>0, p = p + (d-a)];
		If[i!=j<0, w = -w; p = d-(p-a)];
		Return[{p,w}]
]

(* Find interior points of billiard path and list of corresponding sides *)
(* Return {{p,i}} *)
Hits[vin_,vertex_,table_,num_]:= Block[{glue,hits,i,p,q,v},
		hits={};
		v = vin;
		p = table[[1,vertex,1]];
		glue = table[[2]];
		If[!Inward[v,vertex,table],Return[hits]];
		Do[	q= NextPI[{p,v},table];
			If[q=={},Return[hits]];
			{q,i} = Simplify[q];
			If[AtVertex[q,table],Return[hits]];
			hits=Append[hits,{q,i}];
			{p,v}=Glue[{q,v},i,table];
			{p,v}=Simplify[{p,v}];
			hits=Append[hits,{p,glue[[i]]}],
		{num}];
		Return[hits]
]

(* Determine point and vector {p,v} to initiate billiard path, given vertex, p and slope *)
(* p in [0,1] gives location along an edge; p=0 starts at vertex *)
(* v=0 if there is no trajectory at the given slope *)
Lift[slope_,vertex_,edgeloc_,table_]:=Block[{e,p,p0,p1,v},
		v = UnitVec[slope];
		p0 = table[[1,vertex,1]];
		p1 = table[[1,vertex,2]];
		e = p1-p0;
		p = p0 + edgeloc*e;
		If[edgeloc != 0 && Det[{e,v}] < 0,v=-v];
		If[edgeloc==0 && !Inward[v,vertex,table],v=-v];
		If[edgeloc==0 && !Inward[v,vertex,table],v=0];
		Return[{p,v}]
]

(* --------------------Generate billiard path ------------------------*)
Bounces[slope_,vertex_,table_,num_]:= Bounces[slope,vertex,0,table,num]

Bounces[slope_,vertex_,edgeloc_,table_,num_]:= Block[{i,p,p0,v0,q,segs,v},
		segs={};
		{p,v} = {p0,v0} = 
			Lift[slope,vertex,edgeloc,table];
		If[v==0,Return[segs]];
		Do[	q= NextPI[{p,v},table];
			If[q=={},Return[segs]];
			{q,i} = Simplify[q];
			segs=Append[segs,{p,q}];
			If[AtVertex[q,table],Return[segs]];
			{p,v}=Glue[{q,v},i,table];
			{p,v}=Simplify[{p,v}];
			If[AtSame[{p,v},{p0,v0}],Return[segs]],
		{num}];
		Return[segs]
]

(*---------------------Graphics----------------------*) 
SetOptions[Graphics,AspectRatio->Automatic]

GraphTable[]:=Join[{TableThickness},Map[Line[#]&,table[[1]]]];

GraphShot[ps_]:=Join[{ShotThickness},Map[Line[#]&,ps]];

DiChain[ps_]:= Block[{g},
		g=Graphics[Join[GraphTable[],GraphShot[ps]]];
		Export["bild.pdf",g];
		If[Verbosity>1,Show[g]]
]

DiTable[]:=DiChain[{}]

DiTable[t_]:=(table=t; DiTable[]);

(*---------------------Trajectories----------------------*) 
NTable[]:=Block[{glue,sides},
	{sides,glue} = table;
	sides = Chop[N[sides,Digits]];
	table={sides,glue};
]

NTable[t_]:=(table=t; NTable[])

(* Shot[slope,length,starting vertex or {x,y}] *)
MaxBounce = 100
SegLength[{a_,b_}]:=Norm[a-b]
SegsLength[ss_]:=Plus @@ Map[SegLength,ss]

Shot[] := Shot[1,1,0,MaxBounce]
Shot[slope_]:=Shot[slope,1,0,MaxBounce]
Shot[slope_,vertex_]:=Shot[slope,vertex,0,MaxBounce]
Shot[slope_,vertex_,edgeloc_] := Shot[slope,vertex,edgeloc,MaxBounce]
Shot[slope_,vertex_,edgeloc_,num_] := Block[{len,p,segs,v},
	NTable[];
	segs = Bounces[slope,vertex,edgeloc,table,num];
	len  = SegsLength[segs]//N;
	If[Verbosity > 0,
	Print[" Slope: ",N[slope]," : ",
		Length[segs]," segments, Length ",
		len];
	If[Length[segs]==num,
		Print["Stopped after reaching ",num," bounces"]];
	];
	DiChain[Join[table[[1]],segs]];
	Return[len]
]

(* Draw shot in both directions when it starts from the middle of an edge *)
(* Handles vertex connections *)
BiShot[slope_,vertex_,edgeloc_] := BiShot[slope,vertex,edgeloc,MaxBounce]
BiShot[slope_,vertex_,edgeloc_,num_] := Block[{a,b,flip,len,p,segs1,seg2, segs,v},
	NTable[];
	{a,b} = table[[1,vertex]];
	flip = Slope[RefVec[UnitVec[slope],a-b]];
	segs1 = Bounces[slope,vertex,edgeloc,table,num];
	segs2 = Bounces[flip,vertex,edgeloc,table,num];
	segs = Join[segs1,segs2];
	len  = SegsLength[segs]//N;
	If[Verbosity > 0,
	Print[" Slope: ",N[slope]," : ",
		Length[segs]," segments, Length ",
		len];
	If[Length[segs]==num,
		Print["Stopped after reaching ",num," bounces"]];
	];
	DiChain[Join[table[[1]],segs]];
	Return[len]
]

Shots[slope_]:=Shots[slope,MaxBounce,0]

Shots[slope_,n_]:=Shots[slope,n,0]

MidShots[slope_]:=Shots[slope,MaxBounce,1/2]

Shots[slope_,num_,edgeloc_]:=Block[{bs,ls,segs},
	NTable[];
	Print["Trajectories at slope ",slope];
	bs = Table[Bounces[slope,i,edgeloc,table,num],
			{i,1,Length[table[[1]]]}];
	segs = Join @@ bs;
	DiChain[Join[table[[1]],segs]];
	ls = Map[Length,bs];
	If[Max[ls]==num,
		Print["MaxBounce = ",num," attained"]];
	Return[ls]
]

(*----------------Saturation------------------*)
Chain[xs_]:=Drop[Transpose[{xs,RotateLeft[xs]}],-1]

(* Add vertices to a given side *)
AddVs[ps_,table_,i_]:=Block[{a,b,qs},
		{a,b} = table[[1,i]];
		qs = Sort[ps,#1.(b-a) < #2.(b-a) &];
		qs = Join[{a},qs,{b}];
		Return[Chain[qs]];
]

(* Add vertices (p,i) to a table *)
AddVs[ps_,table_]:=Block[{i,is,g,glue,map,n,ns,qs,sides},
	(* Compute sides *)
	n = Length[table[[1]]];
	qs = Table[Map[#[[1]]&,Select[ps,#[[2]]==i&]],{i,1,n}];
	ns = Map[Length,qs];
	sides = Table[AddVs[qs[[i]],table,i],{i,1,n}];
	sides = Join @@ sides;
	(* Compute glue *)
	map = FoldList[Plus,1,1+Drop[ns,-1]];
	glue = table[[2]];
	If[Or @@ (ns != ns[[glue]]), 
		Print["Gluing problem!"];
		Return[table]];
	glue = Table[Table[map[[glue[[i]]]]+j,{i,1,n},{j,ns[[i]],0,-1}]];
	glue = Join @@ glue;
	Return[{sides,glue}]
]


(* Saturate table for a given slope *)
Saturate[s_,table_]:= Block[{i,n,ntable,ps,v},
	v = UnitVec[s];
	n = Length[table[[1]]];
	ps = Table[Hits[v,i,table,100000],{i,1,n}];
	ps = Join @@ ps;
	Return[AddVs[ps,table]]
]

(*----------------Cylinders------------------*)

(* Jump to next side of saturated table *)
(* {jump-vector,newvertex} *)
Jumps[s_,table_]:= 
	Transpose[Table[Jump[s,table,i],{i,1,Length[table[[1]]]}]]

Jump[s_,table_,i_]:= Block[{a,b,glue,j,p,q,sides,v},
	v = UnitVec[s];
	{sides,glue} = table;
	{a,b} = sides[[i]];
	If[Det[{b-a,v}] <= 0,Return[{{0,0}, glue[[i]]}]];
	p = MidPt[table,i];
	{q,j} = NextPI[{p,v},table];
	Return[{q-p,j}]
]

MidPt[table_,i_]:=(table[[1,i,1]] + table[[1,i,2]])/2

(* Return moduli of annuli *)
Moduli[s_]:=Moduli[s,table]

Moduli[s_,table_]:=Block[{hs,ws},
	{hs,ws} = Transpose[Annuli[s,table]];
	mods = hs/ws;
	mods = mods/mods[[1]];
	mods = FullSimplify[mods];
	Return[Rationalize[mods]]
]

(* Returns list of {ht,width} with one entry for each annulus *)
Annuli[s_]:=Annuli[s,table]

Annuli[s_,table_]:=Block[{js,cs,hts,list,perm,sides,stable,u,widths},
	stable = Saturate[s,table];
	Print["Saturated table has ",Length[stable[[1]]]," sides"];
	sides=stable[[1]];
	{js,perm} = Jumps[s,stable];
	u = UnitVec[s];
	js = js.u;
	cs = ToCycles[Abs[perm]];
	widths = Map[Plus @@ js[[#]]&,cs];
	hts = Map[stable[[1,#[[1]]]]&,cs];
	hts = Map[(#[[2]] - #[[1]]).RotV[u]&,hts];
	hts = Abs[hts];
	list = Transpose[{hts,widths}];
	list = Select[list,#[[2]]>Eps&];
	Return[list]
]

(*---------------------Look for periodic slopes-------------------*)
(* Parabolic transformation z -> z+t, rotated to fix s *)
Parabolic[s_,t_,z_]:=Block[{u,v,w},
	{v,u} = UnitVec[s];
	w = If[z==Infinity,-u/v,(u*z + v)/(-v*z + u)];
	w = w + t;
	w = (u*w - v)/(v*w + u);
	Return[Simplify[w]]
]

(* Numerical LCM *)
(* Returns NULL if could not rationalize *)
NLCM[xs_]:=Block[{qs},
	qs = Rationalize[xs/xs[[1]]];
	If[Or @@ Map[!RationalQ[#]&,qs], 
		Print["No LCM for ",qs];
		Return[Null]];
	Return[xs[[1]] * LCM @@ qs]
]

(* Returns Null if next slope cannot be found *)
NextSlope[{s1_,s2_}]:=Block[{h,m,s3,t,w},
	{h,w} = Transpose[Annuli[s2]];
	m = h/w;
	t = NLCM[1/m];
	If[t==Null,Return[Null]];
	s3 = Parabolic[s2,t/2,s1];
	Return[s3]
]

(* Basic slope of table > 1 *)
Slope[]:= 2*NextSlope[{0,Infinity}]

(* Returns {ss,flag} where flag=True if new slope was found *)
Slopes[ss_]:=Block[{snew},
	snew = NextSlope[Take[ss,-2]];
	snew = Chop[snew];
	If[snew==Null,
		Print["Group is not a lattice"];
		Return[{ss,False}]];
	Return[{Append[ss,snew],True}]
]

(*---------------------Produce input to limit set program-------------------*)
LimFile[mirrors_]:=Block[{out,r},
	r = Abs[Flatten[mirrors]];
	r = Select[r,# != Infinity&];
	r = Max[r];
	r = 1.2*r;
	out = OpenWrite["lim.run"]; 
	WriteString[out,"#\n"];
	WriteString[out,"#  Lim input produced from billiards\n"];
	WriteString[out,"#\n"];
	WriteString[out,"lim -w -0.1 0 ",
		ToString[r]," ",
		ToString[r*3/4],
		" -e 0.001 -d 30 <<eof > lim.ps\n"];
	Map[LimLine[#,out]&,mirrors];
	WriteString[out,"eof\n"];
	WriteString[out,"ghostview lim.ps &\n"];
	Close[out];
	Run["chmod +x lim.run"]
]

LimLine[{a_,b_},out_]:=Block[{},
	If[a==Infinity,
		WriteString[out,"c ",ToString[N[b,20]]," 0 -0.5\n"];
		WriteString[out,"r ",ToString[N[b,20]]," 0 -0.5\n"];
		Return[]];
	If[b==Infinity,
		WriteString[out,"c ",ToString[N[a,20]]," 0 -0.5\n"];
		WriteString[out,"r ",ToString[N[a,20]]," 0 -0.5\n"];
		Return[]];
	WriteString[out,"c ",ToString[N[(a+b)/2,20]]," 0 ",
		ToString[N[Abs[a-b]/2,20]],"\n"];
	WriteString[out,"r ",ToString[N[(a+b)/2,20]]," 0 ",
		ToString[N[Abs[a-b]/2,20]],"\n"];
]

(* Calculate the Veech group up to n reflections *)
(* Return a list of mirrors (x1,x2) *)
(* SLX uses {0,Infinity} and {Infinity,0} as initial reflectors *)
SLX[n_]:=Block[{f1,f2,mirrors,ss1,ss2},
	ss1 = {0,Infinity};
	ss2 = {Infinity,0};
	{ss1,f1} = Slopes[ss1];
	{ss2,f2} = Slopes[ss2];
	ss2 = Drop[ss2,1];
	While[(ss1[[-1]] > Eps-ss2[[-1]]) &&
		(Length[ss1]+Length[ss2] < n+2) &&
		(f1 && f2), 
		Print["Slopes:  ",N[ss1[[-1]]]," > ",N[-ss2[[-1]]]];
		If[Length[ss1]==Length[ss2],
			{ss1,f1} = Slopes[ss1],
			{ss2,f2} = Slopes[ss2]
		];
	];
	Print["Slopes:  ",N[ss1[[-1]]]," : ",N[-ss2[[-1]]]];
	mirrors = Join[Chain[ss1],Chain[ss2]];
	If[Length[mirrors] == n,
		Print["MaxMirrors achieved"]];
	Return[mirrors]
]

(* SLX uses {0,Infinity} and {-1,1} as initial reflectors *)
SLY[n_]:=Block[{f1,f2,mirrors,ss1,ss2},
	ss1 = {0,Infinity};
	ss2 = {-1,1};
	{ss1,f1} = Slopes[ss1];
	{ss2,f2} = Slopes[ss2];
	While[(ss1[[-1]] > Eps+1/ss2[[-1]]) &&
		(Length[ss1]+Length[ss2] < n+2) &&
		(f1 && f2),
		Print["Slopes:  ",N[ss1[[-1]]]," > ",N[1/ss2[[-1]]]];
		If[Length[ss1]==Length[ss2],
			{ss1,f1} = Slopes[ss1],
			{ss2,f2} = Slopes[ss2]
		];
	];
	Print["Slopes:  ",N[ss1[[-1]]]," : ",N[1/ss2[[-1]]]];
	mirrors = Join[Chain[ss1],Chain[ss2]];
	If[Length[mirrors] == n,
		Print["MaxMirrors achieved"]];
	Return[mirrors]
]

MaxMirrors=7

DiGp[t_]:=(table=t; DiGp[])

DiGp[]:=Block[{mirrors},
	NTable[];
	DiM[SLX[MaxMirrors]];
]

DiL[d_]:=Block[{mirrors},
	table=LTable[d];
	NTable[];
	DiTable[];
	mirrors = SLY[MaxMirrors];
	DiM[mirrors];
]

DiS[d_]:=Block[{mirrors},
	table=STable[d];
	NTable[];
	DiTable[];
	mirrors = SLX[MaxMirrors];
	DiM[mirrors];
]

DiX[d_]:=Block[{mirrors},
	table=XTable[d];
	NTable[];
	DiTable[];
	mirrors = SLY[MaxMirrors];
	DiM[mirrors];
]

DiM[mirrors_]:=Block[{},
	LimFile[mirrors];
	Run["lim.run"]
]

(*---------------------Tables-------------------------------*)
(* table = {Sides,Glue}				            *) 
(* The sides are given as (pi,qi) enclosing table counterclockwise *)
(* BTable[] turns glued table into billiards table *)
BTable[]:=(table[[2]] = Range[Length[table[[2]]]];)

(* Make n vertices from n vectors *)
Vertices[deltas_]:= FoldList[Plus,{0,0},Drop[deltas,-1]]

(* Make n sides from n vertices *)
Sides[vertices_]:=Simplify[Transpose[{vertices,RotateLeft[vertices]}]]

(*---------------Transform table--------------------------*)
GLTable[m_,table_]:=Block[{glue,sides},
                {sides,glue}=table;
                sides = Map[m.#&,sides,{2}];
                Return[{sides,glue}]
]

(*---------------------L-shaped tables----------------------*) 
(*
   |-1-|
   |   |
   a   |----|
   |        1
   |---b----|
*)
ClassicLTable[d_]:=ClassicLTable[(1+Sqrt[d])/2,(1+Sqrt[d])/2]
ClassicLTable[a_,b_]:= 
	{Sides[{{0,0}, {1,0}, {a,0}, {a,1}, {1,1}, {1,b}, {0,b}, {0,1}}],
	{6,4,8,2,7,1,5,3}}

(* Billiards in Unit square + attached rectangle *)
ModernLTable[a_,b_]:=
	{Sides[{{0,0}, {1,0}, {1+a,0}, {1+a,b}, {1,b}, {1,1}, {0,1}, {0,b}}],
	Range[8]}

(* Table from prototype *)
ProtoLTable[b_,e_]:= Block[{glue,lambda,sides},
	lambda = (e + Sqrt[e^2 + 4b])/2;
	sides = Sides[{{0,0},{lambda,0},{b,0},{b,1},{lambda,1},
			{lambda,1+lambda},{0,1+lambda},{0,1}}];
	glue = {6,4,8,2,7,1,5,3};
	FieldDef[sides];
	Return[{sides,glue}]
]

(*---------------------Barrier tables----------------------*) 
(*
   |---2---|
   |       | 
   1   |   1
   |   h   |
   |-1-|-1-|
*)
BTable[h_]:= 
	{Sides[{{0,0}, {1,0}, {1,h}, {1,0}, {2,0}, {2,1}, {1,1}, {0,1}}],
	{1,2,3,4,5,6,7,8}}
(*	{7,2,3,6,8,4,1,5}} *)

(*---------------------Z-shaped tables----------------------*) 
(* 
   |-a---1-|
   b       b
   |-a-|   |-c-|
       1       1
       |-1---c-|
*)
ZTable[d_]:=ZTable[ d-Sqrt[d], Sqrt[d], (1+d)(Sqrt[d]-1)]
ZTable[a_,b_,c_]:= 
		{Sides[{{0,0},{1,0},{1+c,0},{1+c,1},{1,1},{1,1+b},
			{0,1+b},{-a,1+b},{-a,1},{0,1}}],
		{6,4,10,2,8,1,9,5,7,3}}

(*---------------------N-gon tables----------------------*) 
(* Vertices at roots of unity *)
NgonTable[n_]:= 
		{Sides[
		Table[{Cos[t],Sin[t]},{t,0,2Pi-2Pi/n,2Pi/n}]],
		Range[n]}

(* Resting on [0,1] *)
NgonRestingTable[n_]:=
		{Sides[
			Table[Pt[(Exp[2Pi I k/n]-1)/(Exp[2Pi I/n]-1)],{k,0,n-1}]],
		Range[n]}

(* Opposite sides glued.  N should be even! *)
NgonSurface[n_]:=Switch[Mod[n,2],
			0,NgonSurfaceEven[n],
			1,NgonSurfaceOdd[n]
]

NgonSurfaceOdd[n_]:=Block[{bary,gap,glue,sides0,sides1,sides},
		sides0 = NgonRestingTable[n][[1]];
		gap = Max[Transpose[Flatten[sides0,1]][[1]]]-1/2;
		gap = 2*gap + 0.5;
		bary = (Plus @@ Transpose[sides0][[1]])/n;
		sides1 = Map[#+{gap,0}+2*bary&,-sides0,{2}];
		sides = Join[sides0,sides1];
		glue = Join[n+Range[n],Range[n]];
		Return[{sides,glue}]
]

NgonSurfaceEven[n_]:={NgonRestingTable[n][[1]], 
		Join[n/2+Range[n/2],Range[n/2]]}

(*---------------------Quaternionic / Eierlengende Table----------------------*) 
QTable[]:=
		{Sides[{{0,0},{1,0},{2,0},{3,0},{4,0},{4,1},{5,1},{6,1},{7,1},
		        {7,2},{6,2},{5,2},{4,2},{3,2},{3,1},{2,1},{1,1},{0,1}}],
		{12,13,10,11,18,15,16,17,14,3,4,1,2,9,6,7,8,5}}


(*---------------------Unfolded N-gon tables----------------------*) 
UnfoldedNgonTable[n_]:= Block[{glue,j,k,ps,qs,sides},
		j = (1-n)/2;
		k = (n-1)/2;
		ps = Table[{Cos[2Pi t/n]-Cos[2Pi j/n],Sin[2Pi t/n]},{t,j,k}];
		qs = Map[{-1,1}*#&,ps];
		qs = Drop[Drop[qs,-1],1];
		qs = Reverse[qs];
		ps = Join[ps,qs];
		sides = Sides[ps];
		glue = Join[n-1+Range[n-1],Range[n-1]];
		Return[{sides,glue}]
]

(*---------------------Triangular tables----------------------*) 
(* Interior Angles (a,b,Pi-a-b) *)
TriTable[a_,b_]:= 
		{Sides[
		{{0,0},{1,0},{(1+Tan[a]/Tan[b])^-1,
				(Cot[a]+Cot[b])^-1}}],
		Range[3]};

(*---------------------Rectangular tables----------------------*) 
(* a x b *)
RecTable[a_,b_]:= 
		{Sides[
		{{0,0},{a,0},{a,b},{0,b}}],
		Range[4]};

(*---------------------Torus tables----------------------*) 
(* a x b *)
TorusTable[a_,b_]:= 
		{Sides[
		{{0,0},{a,0},{a,b},{0,b}}],
		{3,4,1,2}}


(*---------------------General polygons----------------------*) 
Polar[{r_,theta_}]:=r*{Cos[theta],Sin[theta]}

(* Internal angles proportional to {a_i} and first n-2 sides given *)
(* First edge at angle GenFirstAngle *)
GenlFirstAngle = 0;
GenlTable[as_,rsin_]:=Block[{deltas,n,rs=rsin,table,ts,vs,x1,x2},
		n = Length[as];
		denom = (Plus @@ as)/(n-2);
		ts = FoldList[Plus,GenlFirstAngle,Pi (1-Drop[as,1]/denom)];
		rs = Join[rs,{x1,x2}];
		deltas = Map[Polar,Transpose[{rs,ts}]];	
		deltas = deltas/.Solve[0==Plus @@ deltas,{x1,x2}][[1]];
		table = {Sides[Vertices[deltas]],Range[n]};
		Return[table]
]

(* Triangles have no moduli *)
TriangleTable[as_]:=GenlTable[as,{1}]

(*---------------------(1,1,1,9) Quad table ----------------------*)
(* a,b meet at 90 degrees *)
Q1119Table[a_,b_]:=GenlTable[{1,9,1,1},{a,b}]

(*---------------------(1,1,1,9) table unfolded ----------------------*)
U1119Table[a_,b_]:=Block[{abs,ang,dang,glue,vs},
        abs = Flatten[Table[{a,b,b,a},{6}]];
        dang = Flatten[Table[{-Pi/2,2Pi/3},{12}]];
        ang = FoldList[Plus,0,Drop[dang,-1]];
        deltas = Map[Polar,Transpose[{abs,ang}]];
        vs = Vertices[deltas];
        glue = {20,7,22,9,24,11,2,13,4,15,6,17,8,19,10,21,12,23,14,1,16,3,18,5};
        table = {Sides[vs],glue};
        Return[table]
]

(*---------------------(1,1,2,8) Quad table ----------------------*)
(* a,b meet at 120 degrees *)
Q1128Table[a_,b_]:= GenlTable[{1,8,2,1},{a,b}]


(*---------------------(1,1,2,8) table unfolded ----------------------*)
(* b = 120 degree edge; a = 60 degree edge *)
U1128Table[a_,b_]:=Block[{abs,ang,dang,glue,vs},
        abs = Flatten[Table[{a,b,b,a},{6}]];
        dang = Flatten[Table[{-Pi/3,Pi/3,-Pi/3,2Pi/3},{6}]];
        ang = FoldList[Plus,0,Drop[dang,-1]];
        deltas = Map[Polar,Transpose[{abs,ang}]];
        vs = Vertices[deltas];
        glue = {20,11,18,9,24,15,22,13,4,19,2,17,8,23,6,21,12,3,10,1,16,7,14,5};
        table = {Sides[vs],glue};
        Return[table]
]

(*---------------------(1,2,1,8) Quad table ----------------------*)
(* a, b incident to obtuse vertex *)
Q1218Table[a_,b_]:= GenlTable[{1,8,1,2},{a,b}]

(*---------------------(1,2,1,8) table unfolded ----------------------*)
U1218Table[a_,b_]:=Block[{abs,ang,dang,glue,vs1,vs2},
        abs = Flatten[Table[{a,b,b,a},{3}]];
        dang = Flatten[Table[{-Pi/3,2Pi/3},{6}]];
        ang = FoldList[Plus,0,Drop[dang,-1]];
        deltas = Map[Polar,Transpose[{abs,ang}]];
        vs1 = Vertices[deltas];
        abs = Flatten[Table[{b,a,a,b},{3}]];
        deltas = Map[Polar,Transpose[{abs,ang}]];
        vs2 = Map[vs1[[7]]+#&,Vertices[deltas]];
        glue = Flatten[Transpose[{{10,12,2,4,6,8},{5,7,9,11,1,3}}]];
	glue = Join[12+glue,glue];
        table = {Join[Sides[vs1],Sides[vs2]],glue};
        Return[table]
]

(*---------------------Bouw-Moeller (2,2,n,3n-4) table --------------------*)
(* Called T[4,n,inf]; n >= 5, odd *)
BMT4[n_]:=GenlTable[{3n-4,2,2,n},{1/2, Cos[Pi/n] + Cos[Pi/4] }]

(*---------------------Bouw-Moeller (1,2,2,4n-5) table --------------------*)
(* Called T[5,n,inf]; n >= 7, odd, not = 0(5) *)
BMT5[n_]:=GenlTable[{4n-5,1,2,2},{Cos[Pi/(2n)], Cos[Pi/n] + Cos[Pi/5] }]

(*---------------------Tables from Coxeter diagrams-------------------*)
(*---------------------Square seed ----------------------*) 
(*		--a--		*)
(*		b   b		*)
(*		--a--		*)

UniSeed[a_,b_]:= {Sides[{{0,0}, {a,0},{a,b},{0,b}}], 
		{3,4,1,2}, {2}}

(*---------------------Three sides seed ----------------------*) 
(*		--b-c-d--		*)
(*		a       a		*)
(*		--b-c-d--		*)

TriSeed[a_,b_,c_,d_]:= {Sides[{{0,0}, {b,0}, {b+c,0}, {b+c+d,0},
			{b+c+d,a},{b+c,a},{b,a},{0,a}}], 
		{7,6,5,8,3,2,1,4},
		{7,2,5}}

(*---------------------Seed + Sprouts ----------------------*) 

(*---------------------Add a new rectangle ----------------------*) 
Grow[{sides_,glue_,leaves_}] :=
	Grow[{sides,glue,leaves},1,1,1]
Grow[{sides_,glue_,leaves_},ht_] :=
	Grow[{sides,glue,leaves},ht,1,1]
Grow[{sides_,glue_,leaves_},ht_,l_] :=
	Grow[{sides,glue,leaves},ht,l,1]

Grow[{sides_,glue_,leaves_},ht_,l_,parity_]:= 
	Block[{a,b,c,d,i,map,nsides,nglue,nleaves,ntree,v},
	i = leaves[[l]];
	{a,d} = sides[[i]];
	v = Sign[Chop[N[RotV[a-d]]]]*ht;
	b = a+v;
	c = d+v;
	nsides=Take[sides,i-1];
	nsides=Join[nsides,{{a,b},{b,c},{c,d}}];
	nsides=Join[nsides,Drop[sides,i]];
	map= Join[Range[i],2+Range[i+1,Length[sides]]];
	nglue = map[[glue]];
	j = nglue[[i]];
	nglue = Join[Take[nglue,i-1],{i+2,j,i},Drop[nglue,i]];
	nglue[[j]] = i+1;
	nleaves = map[[leaves]];
	nleaves[[l]] = If[v[[1+Mod[parity,2]]]==0,i,i+2];
	ntree = {nsides,nglue,nleaves};	
	Return[ntree]
]

(*-----------Find Heights Using Intersections------------*)
(* Convert Coxeter matrix to intersection matrix *)
CoxToInt[q_]:=Map[If[#==-1,1,0]&,q,{2}]

(* Heights from intersection or Coxeter matrix, and moduli *)
(* 	intermat_ij = i(C_i,C_j) *)
(*			    	 *)
Hts[q_]:=Hts[q,Table[1,{Length[q]}]]

Hts[q_,ms_]:=Block[{a,es,hs,i,i1,i2,j,ns,n1,n2,x1,x2},
	a = If[q[[1,1]] == 2, CoxToInt[q],q];
        ns = TwoColor[a];
        a = DiagonalMatrix[ms] . a;
        es=Eigensystem[a];
        i=Position[Re[es[[1]]],Max[Re[es[[1]]]],1][[1,1]];
        hs=es[[2,i]];
        If[ns=={},
                Print["Hts:  No two coloring available"];
                Return[hs]];
        {n1,n2} = ns;
	x1 = Chop[N[hs[[n1]]]];
	x2 = Chop[N[hs[[n2]]]];
	i1=Position[x1,Max[x1]][[1,1]];
	i2=Position[x2,Max[x2]][[1,1]];
        hs[[n1]] = hs[[n1]]/hs[[n1[[i1]]]];
        hs[[n2]] = hs[[n2]]/hs[[n2[[i2]]]];
        hs = FullSimplify[hs];
	FieldDef[hs];
        Return[hs]
]

(*---------------------A_n trees --------------*)
TA[n_Integer]:=TA[n,Table[1,{n}]]

TA[n_,ms_]:=Block[{hts,i,table,tree},
	hts = Hts[QA[n],ms];
	tree= UniSeed[hts[[1]],hts[[2]]];
	Do[tree = Grow[tree,hts[[i]],1,0],{i,3,n}];
	table=Take[tree,2];
	Return[table]
]

TAE[n_Integer]:=TAE[n,Table[1,{n}]]
TAE[n_,ms_]:=Block[{hts,i,glue,table,tree},
	If[Mod[n,2] != 0,
		Print["TAE only works for even n"];
		Return[{}]];
	hts = Hts[QAE[n],ms];
	tree= UniSeed[hts[[1]],hts[[2]]];
	Do[tree = Grow[tree,hts[[i]],1,0],{i,3,n}];
	tree = Grow[tree,hts[[1]],1,0];
	glue = tree[[2]];
	glue[[{1,n,n+2,2n+1}]] = {n+2,2n+1,1,n};
	tree[[2]] = glue;
	table=Take[tree,2];
	Return[table]
]

TAH[n_Integer]:=TAH[n,Table[1,{n}]]
TAH[n_,ms_]:=Block[{hts,i,glue,sides,tree},
	If[Mod[n,2] != 1,
		Print["TAE only works for odd n"];
		Return[{}]];
	hts = Hts[QAH[n],ms];
	tree= TriSeed[hts[[n-1]],hts[[n]],hts[[n-2]],hts[[1]]];
	Do[tree = Grow[tree,hts[[i]],3,0],{i,2,n-3}];
	tree = Grow[tree,hts[[n-3]],2,0];
	sides = tree[[1]];
	sides = Simplify[sides];
	glue = tree[[2]];
	glue[[{2,4,n+2,n+4}]] = {n+2,n+4,2,4};
	Return[{sides,glue}]
]

(*---------------------Y_{a,b,c} trees --------------*)
TY[{a_,b_,c_}]:=TY[{a,b,c},Table[1,{a+b+c-2}]]

TY[{a_,b_,c_},ms_]:= Block[{ha,hb,hc,hr,hts,i,tree},
	hts = Hts[QY[{a,b,c}],ms];
	hr = hts[[1]];
	ha = Take[hts,{2,a}];
	hb = Take[hts,{a+1,a+b-1}];
	hc = Take[hts,{a+b,a+b+c-2}];
	tree= TriSeed[hr,ha[[1]],hb[[1]],hc[[1]]];
	Do[tree = Grow[tree,ha[[i]],1],{i,2,a-1}];
	Do[tree = Grow[tree,hb[[i]],2],{i,2,b-1}];
	Do[tree = Grow[tree,hc[[i]],3],{i,2,c-1}];
	table=Take[tree,2];
	Return[table]
]
(*---------------------Exception diagrams E_n --------------*)
TE[n_]:=TY[{2,3,n-3}]
TE[n_,ms_]:=TY[{2,3,n-3},ms]

(*-------------- Genus 2:  L Tables -------------------*)
(* w:  1 zero, Coxeter diagram A_4                     *)
(* q:  4 zeros, Coxeter diagram a=b=c-d                *)
(*                                   -e                *)
LTable[n_Integer]:=TA[4,{n,1,1,n}]

LTable[ms_]:=TA[4,ms]

LTableQ[n_Integer]:=LTableQ[{n,1,1,2*n,2*n}]

LTableQ[ms_]:=Block[{a,b,c,d,e,glue,q,sides},
	q = {	{0,2,0,0,0},
		{2,0,2,0,0},
		{0,2,0,1,1},
		{0,0,1,0,0},
		{0,0,1,0,0}
	};
	{a,b,c,d,e} = Hts[q,ms];
	a = {a,0}; c = {c,0}; 
	b = {0,b}; d = {0,d}; e = {0,e};
	sides = Sides[Vertices[
		{c,a,-e,c,e,a,b,-a,-c,-a,d,-c,-d,-b}]];
	glue = {12,-6,5,9,3,-2,14,-10,4,-8,13,1,11,7};
	Return[{sides,glue}]
]

(*---------------Genus 3:  S Tables----------------------*)
(* w: 1 zero; Coxeter diagram E_6  *)
(* q: 1 zero; Coxeter diagram Ah_5 *)

(* SSetup[n_Integer]:= SSetup[{1/2,n,1,1,1,1}] *)
(* SSetup[n_Integer]:= SSetup[{1,n,1,1,1,1}]  *)
SSetup[n_Integer]:= SSetup[{1,n,1,n,1,n}]  

SSetup[ms_]:=Block[{a,b,c,d,e,f},
	{a,b,c,d,e,f} = Hts[QY[{2,3,3}],ms];
	a = {a,0}; d = {d,0}; f = {f,0};
	b = {0,b}; c = {0,c}; e = {0,e};
	Return[{a,b,c,d,e,f}]
]

STable[n_]:=Block[{a,b,c,d,e,f,glue,sides},
	{a,b,c,d,e,f} = SSetup[n];
	sides = Sides[Vertices[	
		{d,a,e,b,f,c,-f,-a,-c,-b,-d,-e}]];
	glue = {11, 8, 12, 10, 7, 9, 5, 2, 6, 4, 1, 3};
	Return[{sides,glue}]
]

(* Genus three q, with one zero *)
STableQ[n_]:=Block[{a,b,c,d,e,f,glue,sides},
	{a,b,c,d,e} = SSetup[n];
	sides = Sides[Vertices[	
		{d,a,e,b,f,c,-f,-a,-c,-b,-d,-e}]];
	glue = {-5,8,12,10,-1,9,-11,2,6,4,-7,3};
	Return[{sides,glue}]
]

AltSTable[n_]:=TY[{2,3,3},{1/2,n,1,1,1,1}]

(*-------------- Genus 4:  X Tables-------------------     *)
(* w:  1 zero, Coxeter diagram pinwheel on 8 vertices      *)
(* q:  1 zero, Coxeter diagram K_3,3 with one edge removed *)

AX[]:= {	{0,1,0,1,1,0,0,0},
		{1,0,1,0,0,1,0,0},
		{0,1,0,1,0,0,1,0},
		{1,0,1,0,0,0,0,1},
		{1,0,0,0,0,0,0,0},
		{0,1,0,0,0,0,0,0},
		{0,0,1,0,0,0,0,0},
		{0,0,0,1,0,0,0,0}
	};

XSetup[n_Integer]:=XSetup[{1,1,1,1,n,n,n,n}]

XSetup[ms_]:=Block[{a,b,c,d,e,f,g,h},
	{a,b,c,d,e,f,g,h} = Hts[AX[],ms];
	b = {b,0}; d = {d,0}; e = {e,0}; g = {g,0};
	a = {0,a}; c = {0,c}; f = {0,f}; h = {0,h};
	Return[{a,b,c,d,e,f,g,h}]
]

XTable[n_]:=Block[{a,b,c,d,e,f,g,h,glue,sides},
	{a,b,c,d,e,f,g,h} = XSetup[n];
	sides = Sides[Vertices[
		{b,f,d,g,c,-g,a,h,-d,-h,-b,-e,-a,e,-c,-f}]];
	glue = {11,16,9,6,15,4,13,10,3,8,1,14,7,12,5,2};
	Return[{sides,glue}]
]

XTableQ[n_]:=Block[{a,b,c,d,e,f,g,h,glue,sides},
	{a,b,c,d,e,f,g,h} = XSetup[n];
	sides = Sides[Vertices[
		{b,f,d,g,c,-g,a,h,-d,-h,-b,-e,-a,e,-c,-f}]];
	glue = {11,-8,9,-14,15,-12,13,-2,3,-16,1,-6,7,-4,5,-10};
	Return[{sides,glue}]
]

(*-------------------- Table for (1/p,1/q,--) triangle -----------------*)
TTri[p_,q_]:=Block[{a,b,c,delta,glue,k,n,sides,vs},
	a = 1/p;
	b = 1/q;
	n = LCM[p,q];
	If[!IntegerQ[n/2] || p!=n,Print["TTri[2n,div of 2n] required"];
			Return[]
	];
	c = 1 - a - b;
	vs = Table[
		{Sin[Pi c] * {Cos[Pi k a],Sin[Pi k a]},
		 Sin[Pi b] * {Cos[Pi (k+1) a],Sin[Pi (k+1) a]}},
		{k,0,2n-2,2} ];
	vs = Flatten[vs,1];
	sides=Sides[vs];
	delta = 2*n*b+1;
	glue=Table[1+Mod[{k-delta,k+1+delta},2*n],{k,0,2*n-2,2}];
	glue=Flatten[glue]; 
	Return[{sides,glue}]
]


(* ---------------- Cathedral tables --------------------*)
Pt[z_]:={Re[z],Im[z]}
Oct[k_]:=I^k * (1-I)/2
CathedralTable[a_,b_]:={Sides[Vertices[Map[Pt,
       {I,Oct[2],-2a,Oct[3],-a,b,Oct[2],Oct[3],-b,-a,-I,a,-b,Oct[4],Oct[1],b,a,Oct[4],2a,Oct[1]}]]],
               {11,14,19,15,17,9,18,20,6,12,1,10,16,2,4,13,5,7,3,8}}

(*---------------------Coxeter routines -------------------------------*)

(* -----------------Matrices---------------- *)
Id[v_]:=IdentityMatrix[Length[v]]

(* Turn a row of matrices into a single matrix *)
UnTensor1[r_]:=Transpose[Join @@ Map[Transpose,r]]

(* Turn a matrix of matrices into a single matrix *)
UnTensor[m_]:=Join @@ Map[UnTensor1,m]

(* p x q zero matrix *)
Zero[p_,q_]:=Table[0,{p},{q}] 

(* Direct sum of 2 square matrices *)
DirSum[m_,q_]:=UnTensor[ {{m,Zero[Length[m],Length[q]]},
			{Zero[Length[q],Length[m]],q}}]

Diag[q_]:=Map[q[[#,#]]&,Range[Length[q]]]

(*------------Two coloring of intermatrix---------------------*) 
TwoColor[q_]:=Block[{M,n,ns},
        n=Length[q];
        M=q.q;
        M=Sign[MatrixPower[M,n]];
        ns={Select[Range[n],M[[1,#]] != 0&],
            Select[Range[n],M[[1,#]] == 0&]};
        If[ns[[1]]=={} || ns[[2]]=={},Return[{}]];
        Return[ns]
]

(* -----------------Diagram operations---------------- *)
(* Inner product for wt. n in Coxeter diagram *)
Wt[n_]:=-2 Cos[Pi/n]

QDrop[q_,j_]:=Drop[Transpose[Drop[q,{j}]],{j}]

(* Glue node i to j in two Coxeter diagrams *)
(* Default is to glue FIRST nodes together; result remains first *)
QJoin[u_,v_]:=QJoin[u,1,v,1]
QJoin[u_,i_,v_,j_]:=Block[{q},
		q=DirSum[u,QDrop[v,j]];
		q[[Length[u]+Range[Length[v]-1],i]] =
		q[[i,Length[u]+Range[Length[v]-1]]] =
			Drop[v[[j]],{j}];
		Return[q]
]

(* Multiple joins - glue 1st nodes together *)
QJoins[{u_,v_}]:=QJoin[u,v]
QJoins[qs_]:=QJoin[qs[[1]],QJoins[Drop[qs,1]]]

(* Linear diagram *)
(*   a   b   c     z   *)
(* O - O - O - ... - O *)
QL[{}]:={{2}}
QL[{n_}]:={{2,Wt[n]},{Wt[n],2}}
QL[ns_]:=QJoin[QL[Take[ns,1]],2,QL[Drop[ns,1]],1]

(* Join ends of diagram *)
QClose[u_]:=Block[{q=u},
		q[[Length[q],1]]=q[[1,Length[q]]]=-1;
		Return[q]
]

QEdge={{2,-1},{-1,2}}
AddNode[q_,i_]:=QJoin[q,i,QEdge,1]
AddNode[q_]:=AddNode[q,Length[q]]

QDouble[q_]:=QJoin[q,Length[q],q,Length[q]]

(* -----------------Spherical Diagrams---------------- *)
(* A_n O--O--O--O--...--O *)
QA[1]:={{2}}
QA[n_]:=AddNode[QA[n-1]]

(* B_n O==O--O--O--O--...--O *)
QB[2]:={{1,-1},{-1,2}}
QB[n_]:=QJoin[QB[2],2,QA[n-1],1]

(* C_n O==O--O--O--O--...--O *)
QC[2]:={{4,-2},{-2,2}}
QC[n_]:=QJoin[QC[2],2,QA[n-1],1]

(* D_n O--O--O...O--O *)
(*        |           *)
(*        O           *)
QD[3]:={{ 2, 0,-1},
	{ 0, 2,-1},
	{-1,-1, 2}}
QD[n_]:=QJoin[QD[3],3,QA[n-2],1]

(* E_n: O--O--O--O--...--O    *)
(*            |               *)
(*            O               *)
QE[6]:=QY[{3,3,2}]
QE[7]:=QY[{2,4,3}]
QE[8]:=QY[{2,3,5}]
QE[n_]:=QY[{2,3,n-3}]
QY[ns_]:=QJoins[Map[QA,ns]]

(* F_4 O--O=>=O--O *)
QF[3]:={{ 4,-2, 0},
	{-2, 4,-2},
	{ 0,-2, 2}}
QF[4]:= AddNode[QF[3]]

(* G_2:  O-6-O *)
QG[2]:={{6,-3},{-3,2}}

(* H_3, H_4: O-5-O--O--O *)
QH[3]:=QL[{5,3}]
QH[4]:=QL[{5,3,3}]

(* -----------------Euclidean Diagrams---------------- *)
QAE[n_]:=QClose[QA[n]]
QBE[n_]:=QJoin[QB[n-2],n-2,QD[3],3]
QCE[n_]:=QJoin[QC[n-1],n-1,QC[2],2]
QDE[n_]:=QJoin[QD[n-2],n-2,QD[3],3]
QEE[n_]:=AddNode[QE[n-1]]
QFE[5]:=AddNode[QF[4]]
QGE[3]:=AddNode[QG[2]]

(* -----------------Hyperbolic Diagrams---------------- *)
QAH[n_]:=AddNode[QAE[n-1]]
QBH[n_]:=AddNode[QBE[n-1]]
QDH[n_]:=AddNode[QDE[n-1]]
QEH[n_]:=AddNode[QEE[n-1]]
QK[ns_]:=AddNode[QL[ns],Length[ns]]

(* Q_n:  Loop with 1 edge O==O *)
QQ[n_]:=QClose[QL[Join[{4},Table[3,{n-2}]]]]

(* X_5,6 *) 
QX[5]:=QJoins[{QA[2],QA[2],QA[2],QL[{4}]}]
QX[6]:=QY[{2,2,2,2,2}]

(*----------- Unimodular Lattices --------------*)
(* I_{n,1} *)
QSL[2,1]=QL[{3,Infinity}]

(* II_{n,1} *)
QII[9,1]:= QE[10]
QII[17,1] := QDouble[QE[10]]

(* Complete bipartite graphs *)
(* The first m nodes join the last n nodes *)
QCB[m_,n_]:=Table[If[Min[i,j]<=m && Max[i,j]>m,-1,0],
			{i,m+n},{j,m+n}] +
		2*IdentityMatrix[m+n]

(* The complete graph *)
QCG[n_]:=Table[If[i != j,-1,0],{i,n},{j,n}] +
		2*IdentityMatrix[m+n]

	
(*---------------ToCycles---------------------- *)
ToCycles[p_] :=
        Module[{k, j, first, np = p, q = Table[0, {Length[p]}]},
               DeleteCases[
                   Table[If[np[[i]] == 0,
                            {},
                            j = 1; first = np[[i]]; np[[i]] = 0;
                            k = q[[j++]] = first;
                            While[np[[k]] != 0, q[[j++]] = np[[k]]; np[[k
]] = 0; k = q[[j-1]]];
                            Take[q, j-1]
                         ],
                         {i, Length[p]}
                   ],
                   _?(#==={}&)
               ]
        ]

(*---------------Guess field of definition---------------------- *)
FieldDef[xs_]:=Block[{ps},
		ps = Position[xs,Root];
		ps = Map[Drop[#,-1]&,ps];
		If[ps != {}, Print["Field: Q(",
			Extract[xs,ps[[1]]],")"];
			Return[]];
		ps = Position[xs,Sqrt[x_]];
		If[ps != {}, Print["Field: Q(",
			Extract[xs,ps[[1]]],")"];
			Return[]];
		Print["Field: Q"]
]

(* Characteristic Polynomial *)
CharPoly[q_] := CharPoly[q,Table[1,{Length[q]}]]

(* Discriminant of quadratic polynomial *)
QuadDisc[p_]:=Simplify[
	Coefficient[p,x,1]^2
	-4*Coefficient[p,x,2]*Coefficient[p,x,0]]


(* Perron-Frobenius matrix for weighted Thurston system *)
(* Its leading eigenvalue is mu^2 *)
(* Dehn twists are {{1,mu},{0,1}},{{1,0},{-mu,1}} *)
PFMat[q_,ms_]:= Block[{a,ns,p},
	a = If[q[[1,1]] == 2, CoxToInt[q],q];
	ns = TwoColor[a];
	a = DiagonalMatrix[ms].a; 
	a = a.a;
	a = a[[ns[[1]],ns[[1]]]];
	Return[a]
]

PFPoly[q_,ms_]:= Block[{a,ns,p},
	a = PFMat[q,ms];
	p = Det[x IdentityMatrix[Length[a]] - a];
	Return[Factor[p]]
]

(* lambda + 1/lambda = tau = mu^2-2; return polynomial for tau *)
PFTrace[q_,ms_]:= Factor[PFPoly[q,ms]/.(x->x+2)]

PFSalem[q_,ms_]:=Factor[Expand[x^Length[PFMat[q,ms]] * (PFTrace[q,ms]/.(x->(x+1/x)))]]

PFNo[q_,ms_]:=Max[Abs[x/.Solve[PFSalem[q,ms]==0]]]//N

PFWts[q_]:= Table[1,{Length[q]}]
PFMat[q_]:= PFMat[q,PFWts[q]]
PFPoly[q_]:= PFPoly[q,PFWts[q]] 
PFSalem[q_]:= PFSalem[q,PFWts[q]] 
PFTrace[q_]:= PFTrace[q,PFWts[q]] 
PFNo[q_]:= PFNo[q,PFWts[q]]

(*-------------Billiards through point inside of table--------------*) 
(* Given slope and point, find {i,x} such that a shot started at *)
(* {i,x} passes through the point.  *) 

(* Barycenter of table *)
TableBary[]:=TableBary[table]
TableBary[table_]:=(Plus @@ Flatten[table[[1]],1])/(2*Length[table[[1]]])

EdgeLength[i_]:=Norm[table[[1,i,2]] - table[[1,i,1]]]

(* Turn {p,i} into {i,x} *)
ShotIn[{p_,i_}]:={i,Norm[p-table[[1,i,1]]]/EdgeLength[i]} // Simplify

(* Given slope and point, find ShotIn *)

ShotIn[s_,p_]:=ShotIn[NextPI[{p,UnitVec[s]},table]]

ShotCenter[s_]:=Block[{i,x},
		NTable[];
		{i,x} = ShotIn[s,TableBary[]];
		Shot[s,i,x]
]
