(* ::Package:: *)

(* ::Title:: *)
(*DrawFunDoms*)


(* ::Subsubtitle:: *)
(*Version 1.1 (June 2022)*)


(* ::Text:: *)
(*This package was written in 2020/2021 by Paul Kainberger as part of his master's thesis under supervision of Univ.-Prof. Dr. Peter Paule (RISC). The thesis, which describes the program and provides mathematical theory behind it, is available to download for free at https://www3.risc.jku.at/research/combinat/software/DrawFunDoms/.*)


(* ::Subsection::Closed:: *)
(*License*)


(* ::Text:: *)
(*Copyright (C) 2021 Paul Kainberger*)
(**)
(*This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, version 3.*)
(**)
(*This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.*)
(**)
(*You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/. *)


(* ::Subsection:: *)
(*Version History*)


(* ::Subsubsection:: *)
(*Version 1.1 (June 2022)*)


(* ::Text:: *)
(*Bug fix for cusp width: The width of the cusps for congruence subgroups which do not contain the negative identity matrix is now calculated correctly. Thank you Univ.-Prof. Dr. Peter Paule (RISC) for pointing that out.*)


(* ::Subsubsection:: *)
(*Version 1.0 (August 2021)*)


(* ::Text:: *)
(*First release.*)


(* ::Chapter::Closed:: *)
(*Begin Package*)


BeginPackage["FundamentalDomains`"]
Unprotect@@Names["FundamentalDomains`*"];
ClearAll@@Names["FundamentalDomains`*"];


(* ::Chapter::Closed:: *)
(*Usage*)


$Id::usage="{{1,0},{0,1}}";
$S::usage="{{0, -1},{1, 0}}";
$T::usage="{{1, 1},{0, 1}}";
$SInv::usage="{{0,1},{-1,0}}";
$TInv::usage="{{1,-1},{0,1}}";
$listOfSubgroups::usage="List of all subgroups in form of strings.";


(* ::Section::Closed:: *)
(*Subgroups*)


SubgroupNotation::usage=StringJoin[
"SubgroupNotation[group] returns the notation of ",
"the subgroup in for of a string."
];


ArgumentsNotation::usage=StringJoin[
"ArgumentsNotation[group] returns the notation of ",
"the subgroup's arguments in for of a string."
];


SubgroupMemberQ::usage=StringJoin[
"SubgroupMemberQ[group,{args},A] checks if the matrix A ",
"is a member of the congruence subgroup group(args).\n",
"Returns True or False.\n",
"Example: SubgroupMemberQ[\"Gamma\",{2},{{1,0},{0,1}}]."
];


EquivalentMatricesQ::usage=StringJoin[
"EquivalentMatricesQ[group,{args},M1,M2] checks ",
"whether the matrices M1 \[Element] SL2(\!\(\*TemplateBox[{},\n\"Integers\"]\)) and M2 \[Element] SL2(\!\(\*TemplateBox[{},\n\"Integers\"]\)) ",
"are equivalent in the ",
"congruence subgroup group(args).\n",
"Returns True or False.\n",
"Example: EquivalentMatricesQ[\"Gamma\",{2},{{1,0},{0,1}},{{1,1},{0,1}}]."
];


FindMatrix::usage=StringJoin[
"FindMatrix[z] for ",
"z \[Element] upper half complex plane ",
"returns a matrix M \[Element] SL2(\!\(\*TemplateBox[{},\n\"Integers\"]\)) such that ",
"Inverse[M] z \[Element] standard fundamental domain."
];


RandomSL2Z::usage=StringJoin[
"RandomSL2Z[n] for n \[Element] \!\(\*TemplateBox[{},\n\"PositiveIntegers\"]\) ",
"returns a random matrix M \[Element] SL2(\!\(\*TemplateBox[{},\n\"Integers\"]\)) ",
"such that MaxNorm[M] \[LessEqual] n."
];


GenerateMatrices::usage=StringJoin[
"GenerateMatrices[n] gives a list of all matrices in ",
"\!\(\*SubscriptBox[\(SL\), \(2\)]\)(\!\(\*TemplateBox[{},\n\"Integers\"]\)), ",
"with MaxNorm \[LessEqual] n"];


STExponents::usage=StringJoin[
"Input: M \[Element] SL2(\!\(\*TemplateBox[{},\n\"Integers\"]\)).\n",
"Example: STExponents[{{1,0},{-2,1}}]={{S,-1},{T,2},{S,1}}, ",
"because {{1,0},{-2,1}}=Inverse[S].T.T.S."
];


STExponentsToMatrix::usage=
"STExponentsToMatrix is the inverse function of STExponents.";


STWord::usage=StringJoin[
"STWord[M] for M \[Element] SL2(\!\(\*TemplateBox[{},\n\"Integers\"]\)) returns a string of ",
"concatenated powers of matrices S and T such that ",
"the product of these powers of matrices equals M.\n",
"Example: STWord[{{1,0},{-2,1}}]=\!\(\*SuperscriptBox[\(S\), \(-1\)]\)\!\(\*SuperscriptBox[\(T\), \(2\)]\)S."
];


SubgroupIndex::usage=StringJoin[
"SubgroupIndex[group,{args}] returns the index for the ",
"congruence subgroup group(args) for group \[Element] ",
"{\[CapitalGamma](N),\!\(\*SubscriptBox[\(\[CapitalGamma]\), \(0\)]\)(N), 
\!\(\*SubscriptBox[\(\[CapitalGamma]\), \(1\)]\)(N),
\!\(\*SuperscriptBox[\(\[CapitalGamma]\), \(0\)]\)(N),
\!\(\*SuperscriptBox[\(\[CapitalGamma]\), \(1\)]\)(N),
\!\(\*SuperscriptBox[SubscriptBox[\(\[CapitalGamma]\), \(0\)], \(0\)]\)(N),
\!\(\*SubscriptBox[\(\[CapitalGamma]\), \(\[CapitalTheta]\)]\)}. ",
"For other subgroups, use SubgroupIndex[FunDom[group,{args}]]."
];


(* ::Section::Closed:: *)
(*Calculations*)


GroupAction::usage=StringJoin[
"GroupAction[M,z] for ",
"M={{a,b},{c,d}} \[Element] SL2(\!\(\*TemplateBox[{},\n\"Integers\"]\)), ",
"z \[Element] \!\(\*TemplateBox[{},\n\"Complexes\"]\) ",
"returns (az+b)/(cz+d)."
];


FunDom::usage=StringJoin[
"FunDom[group,{args}] returns a set of coset repersentatives ",
"for the congruence subgroup group(args) together ",
"with a list of adjacencies for each matrix. ",
"Use this function's output as input for the functions ",
"ExchangeMatrix, CosetRepresentatives, CosetRepresentativesWithoutNegatives, ",
"SubgroupIndex, Cusps, NumberOfEllipticPointsOrder2, ",
"NumberOfEllipticPointsOrder3 and Genus.\n",
"Example: FunDom[\"Gamma\",{2}]"
];


ExchangeMatrix::usage=StringJoin[
"ExchangeMatrix[group, {args}, FunDom[group, {args}], M] ",
"exchanges the matrix M \[Element] SL2(\!\(\*TemplateBox[{},\n\"Integers\"]\)) with the equivalent matrix ",
"from FunDom[group, {args}] with respect to the ",
"congruence subgroup group(args).\n",
"Use this function's output as input for the functions ",
"ExchangeMatrix, CosetRepresentatives, CosetRepresentativesWithoutNegatives, ",
"SubgroupIndex, Cusps, NumberOfEllipticPointsOrder2, ",
"NumberOfEllipticPointsOrder3 and Genus."
];


CosetRepresentatives::usage=StringJoin[
"CosetRepresentatives[FunDom[group,{args}]] returns ",
"a list of coset representatives for the congruence ",
"subgroup group(args)."
];


CosetRepresentativesWithoutNegatives::usage=StringJoin[
"CosetRepresentativesWithoutNegatives[FunDom[group,{args}]] returns ",
"a list of coset representatives for the congruence ",
"subgroup group(args), where for each coset representative M, ",
"-M does not appear in that list."
];


Cusps::usage=StringJoin[
"Cusps[FunDom[group,{args}]] returns the cusps ",
"for the congruence subgroup group(args) ",
"together with their widths."
];


NumberOfEllipticPointsOrder2::usage=StringJoin[
"NumberOfEllipticPointsOrder2[FunDom[group,{args}]] ",
"returns the number of elliptic points of order 2 ",
"for the congruence subgroup group(args)."
];


NumberOfEllipticPointsOrder3::usage=StringJoin[
"NumberOfEllipticPointsOrder3[FunDom[group,{args}]] ",
"returns the number of elliptic points of order 3 ",
"for the congruence subgroup group(args)."
];


Genus::usage=StringJoin[
"Genus[FunDom[group,{args}]] ",
"returns the genus ",
"for the congruence subgroup group(args)."
];


(* ::Section::Closed:: *)
(*Drawing*)


Tessellation::usage=StringJoin[
"Tessellation[d,r] with d,r \[Element] \!\(\*TemplateBox[{},\n\"NonNegativeIntegers\"]\) returns ",
"a set of Circles and Lines which represent a modular ",
"tessellation with depth d and range r."
];


DrawMatrixImage::usage=StringJoin[
"DrawMatrixImage[M,n] returns a Polygon with 3n vertices ",
"which represents the triangular image of the standard ",
"fundamental domain under the matrix M."
];


DrawFunDom::usage=StringJoin[
"DrawFunDom[FunDom[group,{args}]] returns objects of type ",
"Polygon which together represent a fundamental domain for ",
"the congruence subgroup group(args)."
];


(* ::Section::Closed:: *)
(*DrawFundamentalDomains*)


DrawFundamentalDomains::usage=StringJoin[
"DrawFundamentalDomains[] opens a graphical user interface ",
"for drawing fundamental domains for several congruence ",
"subgroups."
];


(* ::Chapter:: *)
(*Private*)


Begin["`Private`"]


Id=({{1, 0},{0, 1}});
S=({{0, -1},{1, 0}});
T=({{1, 1},{0, 1}});
SInv=Inverse[S];
TInv=Inverse[T];
TPower[n_]:={{1,n},{0,1}};


$Id=Id;
$S=S;
$T=T;
$SInv=SInv;
$TInv=TInv;


(* ::Section::Closed:: *)
(*SL2(Z) and Congruence Subgroups*)


(* ::Subsection::Closed:: *)
(*Add Congruence Subgroups*)


(* ::Subsubsection::Closed:: *)
(*Step 1 ($listOfSubgroups)*)


$listOfSubgroups={
	"SL2Z",
	"Gamma",
	"Gamma0",
	"Gamma1",
	"GammaUpper0",
	"GammaUpper1",
	"Gamma0Upper0",
	"GammaNMP",
	"GammaTheta"
	(*,"Gamma0CapGamma1"*)
};


(* ::Subsubsection::Closed:: *)
(*Step 2 (SubgroupMemberQ)*)


(* Defining the group Subscript[SL, 2](\[DoubleStruckCapitalZ]) *)
SubgroupMemberQ["SL2Z",{},{{a_,b_},{c_,d_}}]:=
	AllTrue[{a,b,c,d},IntegerQ]&&
	a d-b c==1;

(* Defining the subgroups *)
(* \[CapitalGamma](N) (Principal Congruence Subgroup) *)
SubgroupMemberQ["Gamma",{n_},{{a_,b_},{c_,d_}}]:=
	SubgroupMemberQ["SL2Z",{},{{a,b},{c,d}}]&&
	Mod[{a,b,c,d},n]==Mod[{1,0,0,1},n];

(* Subscript[\[CapitalGamma], 0](N) (Hecke Congruence Subgroup) *)
SubgroupMemberQ["Gamma0",{n_},{{a_,b_},{c_,d_}}]:=
	SubgroupMemberQ["SL2Z",{},{{a,b},{c,d}}]&&
	Mod[c,n]==0;

(* Subscript[\[CapitalGamma], 1](N) *)
SubgroupMemberQ["Gamma1",{n_},{{a_,b_},{c_,d_}}]:=
	SubgroupMemberQ["SL2Z",{},{{a,b},{c,d}}]&&
	Mod[{a,c},n]=={Mod[1,n],0};

(* \[CapitalGamma]^0(N) *)
SubgroupMemberQ["GammaUpper0",{n_},{{a_,b_},{c_,d_}}]:=
	SubgroupMemberQ["SL2Z",{},{{a,b},{c,d}}]&&
	Mod[b,n]==0;

(* \[CapitalGamma]^1(N) *)
SubgroupMemberQ["GammaUpper1",{n_},{{a_,b_},{c_,d_}}]:=
	SubgroupMemberQ["SL2Z",{},{{a,b},{c,d}}]&&
	Mod[{b,d},n]=={0,Mod[1,n]};

(* Subscript[\[CapitalGamma], 0]^0(N) *)
SubgroupMemberQ["Gamma0Upper0",{n_},{{a_,b_},{c_,d_}}]:=
	SubgroupMemberQ["SL2Z",{},{{a,b},{c,d}}]&&
	Mod[b,n]==0&&
	Mod[c,n]==0;

(* \[CapitalGamma](N,M,P) *)
SubgroupMemberQ["GammaNMP",{n_,m_,p_},{{a_,b_},{c_,d_}}]:=
	SubgroupMemberQ["SL2Z",{},{{a,b},{c,d}}]&&
	Mod[a,m]==Mod[1,m]&&
	Mod[b,p]==0&&
	Mod[c,n]==0&&
	Mod[d,m]==Mod[1,m]&&
	Mod[n p,m]==0;

(* Subscript[\[CapitalGamma], \[Theta]] (Modular Group Lambda) *)
SubgroupMemberQ["GammaTheta",{},{{a_,b_},{c_,d_}}]:=
	SubgroupMemberQ["SL2Z",{},{{a,b},{c,d}}]&&
	Mod[a c,2]==0&&
	Mod[b d,2]==0;

(*SubgroupMemberQ["Gamma0CapGamma1",{N1_,N2_},{{a_,b_},{c_,d_}}]:=
	AllTrue[{a,b,c,d},IntegerQ]&&
	a d-b c==1&&
	Mod[c,N1]\[Equal]0&&
	Mod[c,N2]\[Equal]0&&
	Mod[a,N2]\[Equal]Mod[1,N2];*)


(* ::Subsubsection::Closed:: *)
(*Step 3 (SubgroupEmptyQ)*)


SubgroupEmptyQ["GammaNMP",{n_,m_,p_}]/;Mod[n p,m]!=0:=True;
(*SubgroupEmptyQ["Gamma0CapGamma1",_]:=False;*)

SubgroupEmptyQ[_,_]:=False;


(* ::Subsubsection::Closed:: *)
(*Step 4 (SubgroupNotation)*)


SubgroupNotation["SL2Z"]:="\!\(\*SubscriptBox[\(SL\), \(2\)]\)(\!\(\*TemplateBox[{},\n\"Integers\"]\))";
SubgroupNotation["Gamma"]:="\[CapitalGamma]";
SubgroupNotation["Gamma0"]:="\!\(\*SubscriptBox[\(\[CapitalGamma]\), \(0\)]\)";
SubgroupNotation["Gamma1"]:="\!\(\*SubscriptBox[\(\[CapitalGamma]\), \(1\)]\)";
SubgroupNotation["GammaUpper0"]:="\!\(\*SuperscriptBox[\(\[CapitalGamma]\), \(0\)]\)";
SubgroupNotation["GammaUpper1"]:="\!\(\*SuperscriptBox[\(\[CapitalGamma]\), \(1\)]\)";
SubgroupNotation["Gamma0Upper0"]:="\!\(\*SuperscriptBox[SubscriptBox[\(\[CapitalGamma]\), \(0\)], \(0\)]\)";
SubgroupNotation["GammaNMP"]:="\[CapitalGamma]";
SubgroupNotation["GammaTheta"]:="\!\(\*SubscriptBox[\(\[CapitalGamma]\), \(\[Theta]\)]\)";
(*SubgroupNotation["Gamma0CapGamma1"]:="\!\(\*SubscriptBox[\(\[CapitalGamma]\), \(0\)]\) \[Intersection] \!\(\*SubscriptBox[\(\[CapitalGamma]\), \(1\)]\)";*)

SubgroupNotation[group_]:=group;


(* ::Subsubsection::Closed:: *)
(*Step 5 (ArgumentsNotation)*)


ArgumentsNotation["SL2Z"]:="";
ArgumentsNotation["Gamma"]:="N";
ArgumentsNotation["Gamma0"]:="N";
ArgumentsNotation["Gamma1"]:="N";
ArgumentsNotation["GammaUpper0"]:="N";
ArgumentsNotation["GammaUpper1"]:="N";
ArgumentsNotation["Gamma0Upper0"]:="N";
ArgumentsNotation["GammaNMP"]:="N,M,P";
ArgumentsNotation["GammaTheta"]:="";
(*ArgumentsNotation["Gamma0CapGamma1"]:="N1,N2";*)

ArgumentsNotation[_]:="";


(* ::Subsubsection::Closed:: *)
(*Step 6 (ArgumentsNumber)*)


ArgumentsNumber["SL2Z"]:=0;
ArgumentsNumber["Gamma"]:=1;
ArgumentsNumber["Gamma0"]:=1;
ArgumentsNumber["Gamma1"]:=1;
ArgumentsNumber["GammaUpper0"]:=1;
ArgumentsNumber["GammaUpper1"]:=1;
ArgumentsNumber["Gamma0Upper0"]:=1;
ArgumentsNumber["GammaNMP"]:=3;
ArgumentsNumber["GammaTheta"]:=0;
(*ArgumentsNumber["Gamma0CapGamma1"]:=2;*)

ArgumentsNumber[_]:=0;


(* ::Subsubsection::Closed:: *)
(*Step 7 (ArgumentsMaxValues)*)


ArgumentsMaxValues["Gamma"]:={6};
ArgumentsMaxValues["Gamma0"]:={107};
ArgumentsMaxValues["Gamma1"]:={16};
ArgumentsMaxValues["GammaUpper0"]:={107};
ArgumentsMaxValues["GammaUpper1"]:={16};
ArgumentsMaxValues["Gamma0Upper0"]:={11};
ArgumentsMaxValues["GammaNMP"]:={5,5,5};
(*ArgumentsMaxValues["Gamma0CapGamma1"]:={10,5};*)

ArgumentsMaxValues[group_]:=Table[5,ArgumentsNumber[group]];


(* ::Subsection::Closed:: *)
(*Group Action*)


GroupAction[{{a_,b_},{c_,d_}},t_]:=(a t+b)/(c t+d);


GroupActionReal[{{a_,b_},{c_,d_}},{x_,y_}]:=
	{a d x+b(d+c x)+a c(x^2+y^2),y}/((c x+d)^2+c^2 y^2);


(* ::Subsection::Closed:: *)
(*Equivalent Matrices*)


EquivalentMatricesQ[group_,n_,A_?MatrixQ,B_?MatrixQ]:=SubgroupMemberQ[group,n,A.Inverse[B]];
EquivalentMatricesQ[group_,n_,{A_,{_,_,_},_},{B_,{_,_,_},_}]:=EquivalentMatricesQ[group,n,A,B];
EquivalentMatricesQ[group_,n_,R___,A_?MatrixQ]:=AnyTrue[R,EquivalentMatricesQ[group,n,#,A]&];


(* ::Subsection::Closed:: *)
(*FindMatrix, FindMatrixSTWord, FindMatrixReal*)


FindMatrix[z_]/;Im[z]<=0:=
	Nothing;
FindMatrix[z_]:=
	Module[
		{b,tPower,tPowerInverse,matrix,z1=z},
		(* apply T as often as possible *)
		b=Floor[Re[z]+0.5]; (* Round[0.5] = 0 != 1 *)
		tPower=TPower[b];
		tPowerInverse=TPower[-b];
		z1=GroupAction[tPowerInverse,z1];
		matrix=tPower;

		While[Re[z1]<-0.5||0.5<=Re[z1]||Im[z1]<Sqrt[1-Re[z1]^2],
			(* apply S once *)
			z1=GroupAction[S,z1];
			matrix=matrix.S;

			(* apply T as often as possible *)
			b=Floor[Re[z1]+0.5];
			tPower=TPower[b];
			tPowerInverse=TPower[-b];
			z1=GroupAction[tPowerInverse,z1];
			matrix=matrix.tPower;
		];
		matrix
	];


FindMatrixSTWord[z_]:=
	Module[
		{tpower,matrix=Id,word={},z1=z},
		(* apply T as often as possible *)
		tpower=TPower[-Floor[Re[z]+0.5]];
		z1=GroupAction[tpower,z1];
		AppendTo[word,Superscript["T",tpower[[1,2]]]];
		matrix=tpower.matrix;

		While[Re[z1]<-0.5||0.5<=Re[z1]||Im[z1]<Sqrt[1-Re[z1]^2],
			(* apply S once *)
			z1=GroupAction[S,z1];
			AppendTo[word,"S"];
			matrix=S.matrix;

			(* apply T as often as possible *)
			tpower=TPower[-Floor[Re[z1]+0.5]];
			z1=GroupAction[tpower,z1];
			AppendTo[word,Superscript["T",tpower[[1,2]]]];
			matrix=tpower.matrix;
		];
		{matrix,word}
	];


FindMatrixReal[{x_,y_}]:=
	FindMatrix[x+y I];


(* ::Subsection::Closed:: *)
(*Random Matrix*)


RandomSL2Z[n_:1000]/;n>0:=
	Module[{a,b,c,d,g},
		a=RandomInteger[{-n,n}];
		b=RandomInteger[{-n,n}];
		While[\[Not]CoprimeQ[a,b]\[And]\[Not]((a==0\[And]b==1)\[Or](a==1\[And]b==1)),
			a=RandomInteger[{-n,n}];
			b=RandomInteger[{-n,n}];
		];
		{g,{d,c}}=ExtendedGCD[a,b];
		{{a,b},{-c,d}}];


(* ::Subsection::Closed:: *)
(*Generate Matrices*)


Coprimes[a_]:=Coprimes[a,a];
Coprimes[a_,n_]:=Module[{c=Select[Range[Abs[n]],CoprimeQ[a,#]&]},
Union[c,-c]];


GenerateMatrices[n_]/;n>0:=
	Flatten[Table[GenerateMatrices[i,n],{i,0,n}],1];

GenerateMatrices[0,n_]:=
	Join[{{0,-1},{1,#}}&/@Range[1,n],{{0,1},{-1,#}}&/@Range[0,n]];

GenerateMatrices[a_,n_]:=
	Module[{b,c,d,coprimes,output={}},
		coprimes=Append[Coprimes[a,n],0];
		For[i=1,i<=Length[coprimes],i++,
			For[j=1,j<=Length[coprimes],j++,
				d=(1+coprimes[[i]] coprimes[[j]])/a;
				If[Abs[d]<=n&&IntegerQ[d],
				output=Append[output,{{a,coprimes[[i]]},{coprimes[[j]],d}}]
				]
			]
		];
	output];


(* ::Subsection::Closed:: *)
(*STWord, STExponents, STExponentsToMatrix*)


STExponents[m_]/;SubgroupMemberQ["SL2Z",{},m]:=STExponents[m,{}];

(* Case 1 *)
STExponents[{{1,0},{0,1}},st_List]:=st;
STExponents[{{-1,0},{0,-1}},{}]:={{"S",2}};
STExponents[{{1,b_},{0,1}},st_List]:=Append[st,{"T",b}];
STExponents[{{-1,b_},{0,-1}},{}]:={{"S",2},{"T",-b}};
STExponents[{{-1,b_},{0,-1}},st_List]:=
	With[{last=Last[st]},
		If[last[[1]]=="S",
			STExponents[{{1,-b},{0,1}},Append[Most[st],{"S",last[[2]]+2}]],
			STExponents[{{1,-b},{0,1}},Append[st,{"S",2}]]
		]
	];

(* Case 2a *)
STExponents[{{a_,b_},{c_,d_}},st_List]/;Abs[a]<Abs[c]:=
	STExponents[{{-c,-d},{a,b}},Append[st,{"S",-1}]];

(* Case 2b *)
STExponents[{{a_,b_},{c_,d_}},st_List]:=
	Module[{qr=QuotientRemainder[a,c],q,r},
		q=qr[[1]];
		r=qr[[2]];
		STExponents[{{r,b-q d},{c,d}},Append[st,{"T",q}]]
	];


STExponentsToMatrix[list__]:=
	Module[{listMatrices,matrixPowers},
		listMatrices=Replace[list,{"S"->S,"T"->T},2];
		matrixPowers=MatrixPower[#[[1]],#[[2]]]&/@listMatrices;
		Fold[Dot,Id,matrixPowers]
	];


STWord::badargs="Bad argument. A 2x2 matrix in the group \!\(\*SubscriptBox[\(SL\), \(2\)]\)(\!\(\*TemplateBox[{},\n\"Integers\"]\)) is expected.";
STWord[{front___,{st_String,pot_Integer},back___}]:=
	Row[If[#[[2]]==1,#[[1]],Superscript[#[[1]],#[[2]]]]&/@{front,{st,pot},back}];
STWord[{{a_,b_},{c_,d_}}]/;SubgroupMemberQ["SL2Z",{},{{a,b},{c,d}}]:=
	STWord[STExponents[{{a,b},{c,d}},{}]];
STWord[{}]:=1;
STWord[_]:="nothing"/;Message[STWord::badargs]


(* ::Subsection::Closed:: *)
(*MatrixPowerOfTQ*)


MatrixPowerOfTQ[{{1,_},{0,1}}]:=True;
MatrixPowerOfTQ[{{-1,_},{0,-1}}]:=True;
MatrixPowerOfTQ[{{_,_},{_,_}}]:=False;


(* ::Section:: *)
(*Fundamental Domains*)


(* ::Subsection::Closed:: *)
(*FunDom*)


(* ::Text:: *)
(*The main idea of this algorithm was done by Helena A. Verrill.*)


(* Information stored in the FunDom list *)
(*
   group ... the group name as String
n ... the group's arguments, e.g. {5} for \[CapitalGamma](5)
{M_1,M_2,...,M_m} ... the list of found coset representatives together with their information
d ... the distance of the matrices we are currently looking at. Throughout the algorithm, this number increases
True/False ... this boolean tells us whether -Id is a group element.
                If that is not the case, everytime a new coset rep is found, its negative form is also added to the list of coset reps.
*)

(* Information for each M_i *)
(*
{
{{a,b},{c,d}},
{
{checked in dir T?, for which j do we have equivalence?, equivalence => equality?},
{checked in dir T^{-1}?, for which j do we have equivalence?, equivalence => equality?},
{checked in dir S?}, for which j do we have equivalence?, equivalence => equality?}
}
distance to I
}
e.g. { {{1,1},{1,0}}, {{True,2,False},{True,2,True},{True,3,False}}, 2 }
*)

(* If there are no matrices in the group, the fundamental domain spans the whole upper half complex plane *)
FunDom[group_,args_]/;SubgroupEmptyQ[group,args]:="SL2Z";

(* STEP 1 *)
FunDom[group_,n_]/;SubgroupMemberQ[group,n,-Id]:=
	FunDom[group,n,{{Id,{{False,0,False},{False,0,False},{False,0,False}},0}},0,True];
FunDom[group_,n_]:=
	FunDom[group,n,{{Id,{{False,0,False},{False,0,False},{False,0,False}},0},
	{-Id,{{True,0,False},{True,0,False},{True,0,False}},0}},0,False];

(* STEP i+1 *)
(* Debug stop at d = 30 *)
FunDom[group_,n_,R___,30,isNegIdInGroup_]:=-1;
(* Check all matrices having the current distance d in direction T *)
FunDom[group_,n_,{front___,{M_,{{False,0,False},dirTInv_,dirS_},d_},back___},d_,isNegIdInGroup_]:=
	Module[{MT=M.T,reps,equivMatrixPos},
	reps=Map[#[[1]]&,{front,{M,{{False,0,False},dirTInv,dirS},d},back}];
	equivMatrixPos=FirstPosition[reps,x_/;EquivalentMatricesQ[group,n,MT,x],{-1}][[1]];
	If[\[Not]isNegIdInGroup&&EvenQ[equivMatrixPos],equivMatrixPos--]; (* Problem with cusps at e.g. FunDom["Gamma1",{4}]: STTS.T ~ -STTS *)
	If[Or[
		MemberQ[reps,MT],(* MT is already a found rep *)
		equivMatrixPos!=-1],(* MT is equivalent to at least one rep *)
		FunDom[group,n,
		{front,{M,{{True,equivMatrixPos,Or[MT==reps[[equivMatrixPos]],-MT==reps[[equivMatrixPos]]]},dirTInv,dirS},d},back},
		d,isNegIdInGroup], (* No new coset rep *)
		FunDom[group, n,
		{front,{M,{{True,Length[{front}]+2+Length[{back}],True},dirTInv,dirS},d},back,
			{MT,{{False,0,False},{True,Length[{front}]+1,True},{False,0,False}},d+1},
			If[\[Not]isNegIdInGroup,
				{-MT,{{True,0,False},{True,0,False},{True,0,False}},d+1},
				Nothing]},
		d,isNegIdInGroup]
	](* new coset rep found *)
];

(* Check all matrices having the current distance d in direction T^{-1} *)
FunDom[group_,n_,{front___,{M_,{dirT_,{False,0,False},dirS_},d_},back___},d_,isNegIdInGroup_]:=
	Module[{MTInv=M.TInv,reps,equivMatrixPos},
	reps=Map[#[[1]]&,{front,{M,{dirT,{False,0,False},dirS},d},back}];
	equivMatrixPos=FirstPosition[reps,x_/;EquivalentMatricesQ[group,n,MTInv,x],{-1}][[1]];
	If[\[Not]isNegIdInGroup&&EvenQ[equivMatrixPos],equivMatrixPos--];
	If[Or[
		MemberQ[reps,MTInv],(* MTInv is already a found rep *)
		equivMatrixPos!=-1],
		FunDom[group,n,
			{front,{M,{dirT,{True,equivMatrixPos,Or[MTInv==reps[[equivMatrixPos]],-MTInv==reps[[equivMatrixPos]]]},dirS},d},back},
			d,isNegIdInGroup], (* No new coset rep *)
		FunDom[group, n,
			{front,{M,{dirT,{True,Length[{front}]+2+Length[{back}],True},dirS},d},back,
			{MTInv,{{True,Length[{front}]+1,True},{False,0,False},{False,0,False}},d+1},
			If[\[Not]isNegIdInGroup,
				{-MTInv,{{True,0,False},{True,0,False},{True,0,False}},d+1},Nothing]},
				d,isNegIdInGroup]
	](* new coset rep found *)
];

(* Check all matrices having the current distance d in direction S *)
FunDom[group_,n_,{front___,{M_,{dirT_,dirTInv_,{False,0,False}},d_},back___},d_,isNegIdInGroup_]:=
	Module[{MS=M.S,reps,equivMatrixPos},
	reps=Map[#[[1]]&,{front,{M,{dirT,dirTInv,{False,0,False}},d},back}];
	equivMatrixPos=FirstPosition[reps,x_/;EquivalentMatricesQ[group,n,MS,x],{-1}][[1]];
	If[\[Not]isNegIdInGroup&&EvenQ[equivMatrixPos],equivMatrixPos--];
	If[Or[
		IntersectingQ[reps,{MS,-MS}],(* \pm MS is already a found rep *)
		equivMatrixPos!=-1],
	FunDom[group,n,
		{front,{M,{dirT,dirTInv,{True,equivMatrixPos,Or[MS==reps[[equivMatrixPos]],-MS==reps[[equivMatrixPos]]]}},d},back},
		d,isNegIdInGroup], (* No new coset rep *)
	FunDom[group, n,
		{front,{M,{dirT,dirTInv,{True,Length[{front}]+2+Length[{back}],True}},d},back,
		{MS,{{False,0,False},{False,0,False},{True,Length[{front}]+1,True}},d+1},
		If[\[Not]isNegIdInGroup,{-MS,{{True,0,False},{True,0,False},{True,0,False}},d+1},Nothing]},
		d,isNegIdInGroup]
	](* new coset rep found *)
];

(* All matrices having the current distance d are checked in all directions. We increment the current distance by 1 *)
FunDom[group_,n_,R__,d_,isNegIdInGroup_]:=FunDom[group,n,R,d+1,isNegIdInGroup];

(* All matrices are checked in all directions. We have a set of representatives *)
FunDom[group_,n_,
	{front___,{M_,
	{
		{True,equivMatrixTPos_,equivMatrixT_},
		{True,equivMatrixTInvPos_,equivMatrixTInv_},
		{True,equivMatrixSPos_,equivMatrixS_}
	},
	d_}},d_,isNegIdInGroup_]:=
		Delete[#,{{2,1,1},{2,2,1},{2,3,1}}]&/@
		Join[Most/@{front},
		{{M,
		{
			{True,equivMatrixTPos,equivMatrixT},
			{True,equivMatrixTInvPos,equivMatrixTInv},
			{True,equivMatrixSPos,equivMatrixS}
		}}}];


(* ::Subsection::Closed:: *)
(*Exchange Matrix for Equivalent One*)


ExchangeMatrix[group_,n_,funDom_,m_]/;SubgroupMemberQ["SL2Z",{},m]:=
	Module[
	{
		reps=#[[1]]&/@funDom,
		negIdInGroup=SubgroupMemberQ[group,n,-Id],
		equivMatrixPos,equivMatrixDirTPos,equivMatrixDirTInvPos,equivMatrixDirSPos,
		posT,posTInv,posS,
		mT=m.T,mTInv=m.TInv,mS=m.S,
		result=funDom
	},

	(* Check if m is already an existing matrix in the set of representatives *)
	If[MemberQ[reps,m],Return[funDom]];

	(* Find equivalent matrix and exchange *)
	equivMatrixPos=FirstPosition[reps,x_/;EquivalentMatricesQ[group,n,m,x],{-1}][[1]];
	If[\[Not]negIdInGroup,
		If[EvenQ[equivMatrixPos],
			equivMatrixPos--];
		result[[equivMatrixPos+1,1]]=-m];
	result[[equivMatrixPos,1]]=m;

	(* Update list of adjacencies for new matrix *)
	(* Direction T *)
	equivMatrixDirTPos=FirstPosition[reps,x_/;EquivalentMatricesQ[group,n,mT,x],{-1}][[1]];
	If[\[Not]negIdInGroup,
		If[EvenQ[equivMatrixDirTPos],
			equivMatrixDirTPos--]];
	result[[equivMatrixPos,2,1,1]]=equivMatrixDirTPos;
	If[result[[equivMatrixDirTPos,1]]==mT\[Or]-result[[equivMatrixDirTPos,1]]==mT,
		result[[equivMatrixPos,2,1,2]]=True,
		result[[equivMatrixPos,2,1,2]]=False];

	(* Direction TInv *)
	equivMatrixDirTInvPos=FirstPosition[reps,x_/;EquivalentMatricesQ[group,n,mTInv,x],{-1}][[1]];
	If[\[Not]negIdInGroup,
		If[EvenQ[equivMatrixDirTInvPos],
			equivMatrixDirTInvPos--]];
	result[[equivMatrixPos,2,2,1]]=equivMatrixDirTInvPos;
	If[result[[equivMatrixDirTInvPos,1]]==mTInv\[Or]-result[[equivMatrixDirTInvPos,1]]==mTInv,
		result[[equivMatrixPos,2,2,2]]=True,
		result[[equivMatrixPos,2,2,2]]=False];

	(* Direction S *)
	equivMatrixDirSPos=FirstPosition[reps,x_/;EquivalentMatricesQ[group,n,mS,x],{-1}][[1]];
	If[\[Not]negIdInGroup,
		If[EvenQ[equivMatrixDirSPos],
			equivMatrixDirSPos--]];
	result[[equivMatrixPos,2,3,1]]=equivMatrixDirSPos;
	If[result[[equivMatrixDirSPos,1]]==mS\[Or]-result[[equivMatrixDirSPos,1]]==mS,
		result[[equivMatrixPos,2,3,2]]=True,
		result[[equivMatrixPos,2,3,2]]=False];

	(* Update lists of adjacencies for other matrices *)
	(* Matrix having T-neighbour equivalent to m *) 
	posT=FirstPosition[#[[2,1,1]]&/@funDom,x_/;x==equivMatrixPos][[1]];
	If[posT!=equivMatrixPos,
		If[funDom[[posT,1]].T==m\[Or]-funDom[[posT,1]].T==m,
			result[[posT,2,1,2]]=True,
			result[[posT,2,1,2]]=False
		]
	];

	(* Matrix having TInv-neighbour equivalent to m *) 
	posTInv=FirstPosition[#[[2,2,1]]&/@funDom,x_/;x==equivMatrixPos][[1]];
	If[posTInv!=equivMatrixPos,
		If[funDom[[posTInv,1]].TInv==m\[Or]-funDom[[posTInv,1]].TInv==m,
			result[[posTInv,2,2,2]]=True,
			result[[posTInv,2,2,2]]=False
		]
	];

	(* Matrix having S-neighbour equivalent to m *) 
	posS=FirstPosition[#[[2,3,1]]&/@funDom,x_/;x==equivMatrixPos][[1]];
	If[posS!=equivMatrixPos,
		If[funDom[[posS,1]].S==m\[Or]-funDom[[posS,1]].S==m,
			result[[posS,2,3,2]]=True,
			result[[posS,2,3,2]]=False
		]
	];

	result
	];


(* ::Subsection::Closed:: *)
(*Coset Representatives*)


(* ::Subsubsection::Closed:: *)
(*With negative Matrices if -Id is not in the Subgroup*)


CosetRepresentatives[group_,n_]:=#[[1]]&/@FunDom[group,n];
(* If there are no matrices in the group, the set of coset representatives equals the whole group SL2(Z) *)
CosetRepresentatives[group_,n_]/;SubgroupEmptyQ[group,n]:="SL2Z";
CosetRepresentatives[list__]:=#[[1]]&/@list;


(* ::Subsubsection::Closed:: *)
(*Without negative Matrices*)


CosetRepresentativesWithoutNegatives[group_,n_]/;SubgroupMemberQ[group,n,-Id]:=#[[1]]&/@FunDom[group,n];
CosetRepresentativesWithoutNegatives[group_,n_]:=#[[1]]&/@Cases[FunDom[group,n],x_/;x[[2,1,1]]>0];

(* In case we already have the result of FunDom[group,n] *)
(* The boolean tells if -Id is in the group *)
CosetRepresentativesWithoutNegatives[list_,True]:=#[[1]]&/@list;
CosetRepresentativesWithoutNegatives[list_,False]:=#[[1]]&/@Cases[list,x_/;x[[2,1,1]]>0];
CosetRepresentativesWithoutNegatives[list_]:=CosetRepresentativesWithoutNegatives[list,False];
(* If there are no matrices in the group, the set of coset representatives equals the whole group SL2(Z) *)
CosetRepresentativesWithoutNegatives["SL2Z",_]:="SL2Z";


(* The matrices of the Form T^n will be in front of the list *)
(* We just use the coset representatives without negative matrices, since the image of those would just be duplicated *)
CosetRepresentativesSplit[group_,n_]:=
	With[{list=GatherBy[CosetRepresentativesWithoutNegatives[group,n],MatrixPowerOfTQ]},
		If[MatrixPowerOfTQ[list[[1,1]]],
			list,
			Reverse[list]
		]
	];


(* ::Subsection::Closed:: *)
(*Subgroup Index*)


SubgroupIndex["Gamma",{n_}]:=n^3 Product[1-(1/p^2),{p,Select[Divisors[n],PrimeQ]}];

SubgroupIndex["Gamma0",{n_}]:=n Product[1+1/p,{p,Select[Divisors[n],PrimeQ]}];

SubgroupIndex["Gamma1",{n_}]:=n^2 Product[1-(1/p^2),{p,Select[Divisors[n],PrimeQ]}];

SubgroupIndex["GammaUpper0",{n_}]:=n Product[1+1/p,{p,Select[Divisors[n],PrimeQ]}]; (* = SubgroupIndex["Gamma0",{n}] *)

SubgroupIndex["GammaUpper1",{n_}]:=n^2 Product[1-(1/p^2),{p,Select[Divisors[n],PrimeQ]}]; (* = SubgroupIndex["Gamma1",{n}]; *)

SubgroupIndex["Gamma0Upper0",{n_}]:=n^2 Product[1+1/p,{p,Select[Divisors[n],PrimeQ]}];

SubgroupIndex["GammaTheta",{}]:=3;

SubgroupIndex[group_,args_]:=Length[FunDom[group,args]];

SubgroupIndex[list__]:=Length[list];


(* ::Subsection:: *)
(*Cusps*)


Cusps["SL2Z"]:={{Rationals,\[Infinity]},{ComplexInfinity,\[Infinity]}};

Cusps[list__]:=
	Module[{unused,isNegIdInGroup,result,startIndex,startMatrix,nextIndex,nextMatrix,cusp,width},

	(* this drops all double counted matrices (we either consider the matrix M or -M, not both) *)
	(* if -Id \[NotElement] subgroup, then we later double the widths *)
	unused=Flatten@Position[#[[2,1,1]]&/@list,x_/;x>0];
	isNegIdInGroup=Length[list]==Length[unused];
	result={};

	While[Length[unused]>0, (* one run for each cusp, calculating the cusp and its width *)
		startIndex=First[unused];
		unused=Rest[unused];
		startMatrix=list[[startIndex]];
		If[startMatrix[[1,2,1]]==0,
			cusp=ComplexInfinity,
			cusp=startMatrix[[1,1,1]]/startMatrix[[1,2,1]]];
		width=1;

		nextIndex=startMatrix[[2,1,1]];
		While[nextIndex!=startIndex, (* one run for each matrix in that cusp *)
			unused=DeleteCases[unused,x_/;x==nextIndex];
			width++;
			nextMatrix=list[[nextIndex]];
			nextIndex=nextMatrix[[2,1,1]];
		];

		result=Append[result,{cusp,width}];
	];

	SortBy[result,First]];


(* ::Subsection::Closed:: *)
(*Number of Elliptic Points*)


(* ::Subsubsection::Closed:: *)
(*Order 2*)


NumberOfEllipticPointsOrder2[list_]:=
	Module[{listWithoutNegatives,isNegIdInGroup,e2=0},
		listWithoutNegatives=Cases[list,x_/;x[[2,1,1]]>0];
		isNegIdInGroup=Length[listWithoutNegatives]==Length[list];
		For[i=1,i<=Length[listWithoutNegatives],i++,
			If[list[[i,2,3,1]]==i,e2++]];
		e2
	];


(* ::Subsubsection::Closed:: *)
(*Order 3*)


NumberOfEllipticPointsOrder3[list_]:=
	Module[{listWithoutNegatives,isNegIdInGroup,e3=0,setOfMatrices,nextMatrix},
		listWithoutNegatives=Cases[list,x_/;x[[2,1,1]]>0];
		isNegIdInGroup=Length[listWithoutNegatives]==Length[list];
		For[i=1,i<=Length[list],i++,
			nextMatrix=list[[i,2,1,1]]; (* .T *)
			setOfMatrices={i,nextMatrix};
			nextMatrix=list[[nextMatrix,2,3,1]]; (* .S *)
			setOfMatrices=Append[setOfMatrices,nextMatrix];
			nextMatrix=list[[nextMatrix,2,1,1]]; (* .T *)
			setOfMatrices=Append[setOfMatrices,nextMatrix];
			nextMatrix=list[[nextMatrix,2,3,1]]; (* .S *)
			setOfMatrices=Append[setOfMatrices,nextMatrix];
			nextMatrix=list[[nextMatrix,2,1,1]]; (* .T *)
			setOfMatrices=Append[setOfMatrices,nextMatrix];
			If[Length[DeleteDuplicates[setOfMatrices]]<=2,e3++];
			If[\[Not]isNegIdInGroup,i++]
		];
		e3
	];


(* ::Subsection::Closed:: *)
(*Genus*)


Genus[list_]:=
	1+(Length[Cases[list,x_/;x[[2,1,1]]>0]]- (* projective index *)
		6 Length[Cusps[list]]- (* number of projective cusps *)
		3 NumberOfEllipticPointsOrder2[list]-
		4 NumberOfEllipticPointsOrder3[list]
	)/12;


(* ::Subsection::Closed:: *)
(*Links*)


(* ::Subsubsection::Closed:: *)
(*Open/Closed Boundaries*)


ClosedLineDirT[matrix_,precision_,0]:=
	ClosedLineDirT[matrix,precision,1];
ClosedLineDirT[matrix_,precision_,t_]:=
	If[MatrixPowerOfTQ[matrix],
		Line[{{-0.5+matrix[[1,2]],Sqrt[3]/2},{-0.5+matrix[[1,2]],100}}],
		Line[GroupActionReal[matrix.SInv,#]&/@
			(ReIm/@Table[Exp[m Pi I/3]+1,{m,Table[3-(-m+1)^t,{m,0,1,1/precision}]}])]
	];


ClosedLineDirS[matrix_,precision_]:=
	ClosedLineDirS[matrix,precision,1];
ClosedLineDirS[matrix_,precision_]:=
	Line[GroupActionReal[matrix.SInv,#]&/@
		(ReIm/@Table[Exp[m Pi I/3],{m,Table[m+1,{m,0,0.5,1/precision}]}])];


(* ::Subsubsection::Closed:: *)
(*Arrows*)


ArrowDirTInv[matrix_,0]:=
	Line[ReIm/@(GroupAction[matrix.SInv,#]&/@{Exp[4.98\[Pi] I/6]+1.031,Exp[4.9\[Pi] I/6]+1,Exp[4.92 \[Pi] I/6]+0.965})];
ArrowDirTInv[matrix_,1]:=
	Line[ReIm/@(GroupAction[matrix.SInv,#]&/@{Exp[5.07\[Pi] I/6]+1.06,Exp[4.9\[Pi] I/6]+1,Exp[4.955 \[Pi] I/6]+0.935})];
ArrowDirTInv[matrix_,2]:=
	Line[ReIm/@(GroupAction[matrix.SInv,#]&/@{Exp[5.15\[Pi] I/6]+1.08,Exp[4.9\[Pi] I/6]+1,Exp[5 \[Pi] I/6]+0.91})];
ArrowDirTInv[matrix_,___]:=
	ArrowDirTInv[matrix,1];

ArrowDirS[matrix_,0]:=
	Line[ReIm/@(GroupAction[matrix.SInv,#]&/@{Exp[4.87/4 \[Pi] I/3]+0.055I,Exp[4.75/4 \[Pi] I/3],Exp[5/4 \[Pi] I/3]-0.05I})];
ArrowDirS[matrix_,1]:=
	Line[ReIm/@(GroupAction[matrix.SInv,#]&/@{Exp[5/4 \[Pi] I/3]+0.095I,Exp[4.75/4 \[Pi] I/3],Exp[5.15/4 \[Pi] I/3]-0.09I})];
ArrowDirS[matrix_,2]:=
	Line[ReIm/@(GroupAction[matrix.SInv,#]&/@{Exp[5.2/4 \[Pi] I/3]+0.15I,Exp[4.75/4 \[Pi] I/3],Exp[5.3/4 \[Pi] I/3]-0.13I})];
ArrowDirS[matrix_,___]:=
	ArrowDirS[matrix,1];


(* ::Section::Closed:: *)
(*Drawing*)


(* ::Subsection::Closed:: *)
(*Modular Tessellation*)


(* ::Text:: *)
(*The algorithm for the following code was created by Jerzy Kocik (https://arxiv.org/pdf/1912.05768.pdf).*)


KocikCircle[k_,n_,t_:0]:=
	Nothing;
KocikCircle[0,2,t_:0]:=
	Line[Scaled[#,{(t+0.5),0}]&/@{{0,0},{0,1}}];
KocikCircle[k_,n_,t_]/;And[EvenQ[n],Mod[n,8]==0]:=
	With[{x=(k^2-1)/n},
		If[OddQ[x],
			Circle[{k/n+t,0},1/n,{0,Pi}]
		]
	];
KocikCircle[k_,n_,t_:0]/;OddQ[n]:=
	If[Mod[k^2-1,n]==0,Circle[{k/n+t,0},1/n,{0,Pi}]];


KocikK[n_,t_:0]:=
	DeleteCases[Table[KocikCircle[k,n,t],{k,0,n-1}],Null];


Tessellation[n_,t_:0]:=
	Flatten[Table[
		Flatten[Table[KocikK[m,translation],{m,0,n}],1],
		{translation,-(t+1),t}
	],1];


(* ::Subsection::Closed:: *)
(*Generate Vertices for Polygons*)


PointsTop[n_,0]:=
	PointsTop[n,1];
PointsTop[n_,t_:1]/;t>0:=
	Table[Exp[m Pi I/3],{m,Table[m^t+1,{m,0,1,1/n}]}];


PointsLeft[n_,0]:=
	PointsLeft[n,1];
PointsLeft[n_,t_:1]/;t>0:=
	Table[Exp[m Pi I/3]-1,{m,Table[m^t,{m,0,1,1/n}]}];


(* The mirrored function of x^t is 1-(-x+1)^t *)
PointsRight[n_,0]:=
	PointsRight[n,1];
PointsRight[n_,t_:1]/;t>0:=
	Table[Exp[m Pi I/3]+1,{m,Table[3-(-m+1)^t,{m,0,1,1/n}]}];


Points[n_,0]:=
	Points[n,1];
Points[n_,t_:1]:=
	Join[Most[PointsTop[n,1]],Most[Reverse[PointsLeft[n,t/2]]],Most[Reverse[PointsRight[n,t/2]]]];

(* We also allow different options for each edge depending on whether the edge is connected to another matrix image. *)
Points[{precisionGood_,precisionBad_},{c_,d_},{dirT_,dirTInv_,dirS_}]:=
	Join[
		Most[If[dirS,PointsTop[precisionBad,1],PointsTop[precisionGood,1]]],
		Most[Reverse[If[dirT,PointsLeft[precisionBad,Abs[d]],PointsLeft[precisionGood,Abs[d]]]]],
		Most[Reverse[If[dirTInv,PointsRight[precisionBad,Abs[d-c]],PointsRight[precisionGood,Abs[d]]]]]
	];


(* ::Subsection::Closed:: *)
(*Draw Polygons*)


DrawInf[M_,precision_Integer,drawToInf_:True]:=
	Module[
		{infPol,infRec,t=M[[1,1]] M[[1,2]]},
		If[drawToInf,
		(
			(* Draw rectangle with an imaginary part of up to infinity *)
			infPol=Polygon[Join[{t,0}+#&/@ReIm[PointsTop[precision]],{{-0.5+t,1.1},{0.5+t,1.1}}]];
			infRec=Rectangle[Scaled[{0,0},{-1/2+t,1}],Scaled[{0,1},{1/2+t,1}]];
			Return[{infPol,infRec}]
		)];
		(* Draw rectangle with a max imaginary part of 100 *)
		Polygon[Join[{{0.5+t,100}},{t,0}+#&/@ReIm[PointsTop[precision]],{{-0.5+t,100}}]]
];

(* For a faster program, we also allow a list of real points as input *)
DrawInf[M_,points___]:=
	Module[
		{t=M[[1,1]] M[[1,2]]},
		Polygon[Join[{{0.5+t,100}},{t,0}+#&/@points,{{-0.5+t,100}}]]
	];

(* dirS is a boolean telling whether there is a directly connecting neighbour in S direction,
meaning that the precision may be very low at this edge *)
DrawInf[{{{a_,b_},{c_,d_}},{{_,_},{_,_},{_,dirS_}}},{precisionGood_,precisionBad_}]:=
	Module[
		{t=a b},
		Polygon[Join[{{0.5+t,100}},{t,0}+#&/@ReIm[PointsTop[If[dirS,precisionBad,precisionGood],1]],{{-0.5+t,100}}]]
	];


DrawNotInf[{{a_,b_},{c_,d_}},precision_Integer]:=
	With[
		{M={{a,b},{c,d}}.SInv},
		Polygon[GroupActionReal[M,#]&/@ReIm[Points[{precision,precision},{Abs[d],Abs[d]},{True,False,False}]]]
	];

(* For a faster program, we also allow a list of real points as input *)
DrawNotInf[{{a_,b_},{c_,d_}},points___]:=
	With[
		{M={{a,b},{c,d}}.SInv},
		Polygon[GroupActionReal[M,#]&/@points]
	];

(* dirX is a boolean telling whether there is a directly connecting neighbour in X direction,
meaning that the precision may be very low at this edge *)
DrawNotInf[{{{a_,b_},{c_,d_}},{{_,dirT_},{_,dirTInv_},{_,dirS_}}},{precisionGood_,precisionBad_}]:=
	With[
		{M={{a,b},{c,d}}.SInv},
		Polygon[GroupActionReal[M,#]&/@ReIm[Points[{precisionGood,precisionBad},{c,d},{dirT,dirTInv,dirS}]]]
	];


DrawMatrixImage[{{-1,b_},{0,-1}},precision_:10,drawToInf_:True]/;precision>0:=
	DrawInf[{{1,-b},{0,1}},precision,drawToInf];
DrawMatrixImage[{{1,b_},{0,1}},precision_:10,drawToInf_:True]/;precision>0:=
	DrawInf[{{1,b},{0,1}},precision,drawToInf];
DrawMatrixImage[M_,precision_:10,drawToInf_:True]:=
	DrawNotInf[M,precision];


(* ::Subsection::Closed:: *)
(*Draw Fundamental Domain*)


DrawFunDom[list_]:=
Module[{listWithoutNegatives,listSplit,precisionGood,precisionBad},
	precisionGood=20;
	precisionBad=40;
	listWithoutNegatives=Cases[list,x_/;x[[2,1,1]]>0];
	listSplit=GatherBy[listWithoutNegatives,MatrixPowerOfTQ[#[[1]]]&];
	If[Length[listSplit[[1]]]==0,
		If[Length[listSplit[[2]]]==0,
			{}, (* there are no matrices in the set of reps *)
			{DrawNotInf[#,{precisionGood,precisionBad}]&/@listSplit[[2]]} (* there are no powers of T in the set listSplit[[2]] *)
		],
		If[Length[listSplit[[2]]]==0,
			{DrawInf[#,{precisionGood,precisionBad}]&/@listSplit[[1]]}, (* there are solely powers of T in the set listSplit[[1]] *)
			{
				DrawInf[#,{precisionGood,precisionBad}]&/@listSplit[[1]],
				DrawNotInf[#,{precisionGood,precisionBad}]&/@listSplit[[2]]
			}
		]
	]
]


(* ::Section::Closed:: *)
(*DrawFundamentalDomains (Main Program)*)


DrawFundamentalDomains[]:=

DynamicModule[
{
	debug,fontSize,fontSizePopupMenu,fontSizeMMatrix,buttonSize,buttonSizeMMatrix,popupMenuSize,
	precisionBad,precisionGood,args,
	isNegIdInGroup,list,listWithoutNegatives,listSplit,repsWithoutNegatives,reps,repsSplit,
	cusps,cuspPoints,cuspPointsWithInf,drawCuspInf,genus,
	matricesOpenInDirT,matricesOpenInDirS,matricesOpenInDirTInv,
	linksPtT={-0.5,1.86603},linksPtS={-0.258819,0.965926},linksMatricesT,linksMatricesS,linksPointsT,linksPointsS,links,
	closedLinesDirT,closedLinesDirS,calculateClosedLinesVertices,closedLinesDirTVertices,closedLinesDirSVertices,arrows,
	points,pointsTop,inf,notInf,funDom,
	currentGroup=-1,currentArgs=-1,matrix=Id,
	(* Tabs *)
	selectedTab=pSelectPoint,currentTab=pSelectPoint,pointIsGood=True,
	(* Tab select point *)
	pMatrix=Id,pMatrixInverse,pEquivMatrix,pTransformationMatrix,pInverseImage,pInverseTransformationMatrix,pMsgPointNotGood,
	pCurrentPoint={0,1},
	btnPromotePMatrixPressed=False ,
	(* Tab select matrix *)
	btnPromoteMMatrixPressed=False,mMatrix=Id,mCurrentPoint={0,1}
},

(* Global adjustments *)
debug=False;
fontSize=10;
fontSizePopupMenu=10;
fontSizeMMatrix=9;
buttonSize={50,25};
buttonSizeMMatrix={55,27};
popupMenuSize={50,25};
precisionGood=30;
precisionBad=15;
drawCuspInf=False;
calculateClosedLinesVertices=True;
pMsgPointNotGood="Imaginary part of selected point is too small.\nPlease select a point with imaginary part >= 0.001";

Manipulate[

	(* Store the control variables in the variable args, which is the input argument for the subgroup *)
	(* Reset input values, if bigger than MaxInputValue[group] *)
	Switch[ArgumentsNumber[group],
		0,args={},
		1,If[ArgumentsMaxValues[group][[1]]<n,
			n=1;m=1;p=1];
			args={n},
		2,If[ArgumentsMaxValues[group][[1]]<n\[Or]
			ArgumentsMaxValues[group][[2]]<m,
			n=1;m=1;p=1];
			args={n,m},
		3,If[ArgumentsMaxValues[group][[1]]<n\[Or]
			ArgumentsMaxValues[group][[2]]<m\[Or]
			ArgumentsMaxValues[group][[3]]<p,
			n=1;m=1;p=1];
			args={n,m,p}
	];

	(*****************************************************************)
	(**************************  F U N D O M  **************************)
	(*****************************************************************)

	(* These calculations are redone, whenever the subgroup or the arguments change,
	 or if the button "Promote selected Matrix to Coset Representative" gets pressed *)
	If[currentGroup=!=group\[Or]currentArgs=!=args\[Or]
		btnPromotePMatrixPressed\[Or]btnPromoteMMatrixPressed,

		currentGroup=group;
		currentArgs=args;

		If[SubgroupEmptyQ[group,args],
		(
			If[debug,Print["___________/nGroup is empty"]];
			list=FunDom[group,args];
			listWithoutNegatives={};
			listSplit={};
			repsWithoutNegatives={};
			reps=list;
			repsSplit={{},{}};
			cusps=Cusps[list];
			cuspPointsWithInf={};
			cuspPoints={};
			links={};
			funDom={Rectangle[{-100,0},{100,10}]};
		),
		(
			isNegIdInGroup=SubgroupMemberQ[group,args,-Id];
			If[debug,Print[{"___________/n",group,args,isNegIdInGroup}]];

			(* Calculate representants *)
			If[debug,Print["calculate reps"]];
			(* Initialise list *)
			If[btnPromotePMatrixPressed\[Or]btnPromoteMMatrixPressed,
			(
				(* Button "Rep." gets pressed *)
				If[debug,Print["Button pressed"]];
				If[btnPromotePMatrixPressed,btnPromotePMatrixPressed=False;matrix=pMatrix,
				btnPromoteMMatrixPressed=False;matrix=mMatrix];
				If[\[Not]SubgroupEmptyQ[group,args]\[And]\[Not]MemberQ[reps,matrix],
					(* Update representants *)
					list=ExchangeMatrix[group,args,list,matrix];
					If[debug,Print[{"Button finished"}]];
				]
			),
				(* Button not pressed *)
				list=FunDom[group,args];
			];
			listWithoutNegatives=Cases[list,x_/;x[[2,1,1]]>0];
			listSplit=GatherBy[listWithoutNegatives,MatrixPowerOfTQ[#[[1]]]&];
			repsWithoutNegatives=CosetRepresentativesWithoutNegatives[list,isNegIdInGroup];
			reps=If[isNegIdInGroup,repsWithoutNegatives,Union[repsWithoutNegatives,-repsWithoutNegatives]];
			repsSplit=GatherBy[repsWithoutNegatives,MatrixPowerOfTQ];
			If[Length[repsSplit]==1,(* there are either solely matrices of a power of T or none of these *)
				listSplit=Append[listSplit,{}];
				repsSplit=Append[repsSplit,{}]];
			If[MatrixPowerOfTQ[repsSplit[[1,1]]],Nothing,
				listSplit=Reverse[listSplit];
				repsSplit=Reverse[repsSplit]]; (* e.g. reps = { {Id,T,T.T,TInv.TInv,T.T.T.T.T,...}, {S,S.T,T.S,S.TInv.S.T,...} } *)
			If[debug,Print[{list,listWithoutNegatives,listSplit,repsWithoutNegatives,reps,repsSplit}]];
			If[debug,Print["reps finished\n"]];

			(* Calculate cusps and genus *)
			If[debug,Print["calculate cusps"]];
			cusps=Cusps[list];
			cuspPointsWithInf=cusps/.{x_,_Integer}->{x,0};
			cuspPoints=cuspPointsWithInf/.{ComplexInfinity,_}->If[drawCuspInf,{0,3},Nothing];
			genus=Genus[list];
			If[debug,Print[{cusps,cuspPointsWithInf,cuspPoints,genus}]];
			If[debug,Print["cusps and genus finished\n"]];

			(* Calculate links *)
			If[debug,Print["calculate links"]];
			(* matrices having no neighbour in direction T/S *)
			matricesOpenInDirT=Cases[listWithoutNegatives,x_/;x[[2,1,2]]==False];
			matricesOpenInDirTInv=Cases[listWithoutNegatives,x_/;x[[2,2,2]]==False];
			matricesOpenInDirS=Cases[listWithoutNegatives,x_/;x[[2,3,2]]==False];
			(* linksMatricesT = {first matrix in direction T not in FunDom, equivalent matrix in FunDom} *)
			linksMatricesT={#[[1]].T,list[[#[[2,1,1]]]][[1]]}&/@matricesOpenInDirT; 
			linksMatricesS={#[[1]].S,list[[#[[2,3,1]]]][[1]]}&/@matricesOpenInDirS;
			linksPointsT={GroupActionReal[#[[1]],linksPtT],GroupActionReal[#[[2]],linksPtT]}&/@linksMatricesT;
			linksPointsS={GroupActionReal[#[[1]],linksPtS],GroupActionReal[#[[2]],linksPtS]}&/@linksMatricesS;
			links=BezierCurve[
			{
				#[[1]],
				{x424242,y424242}/.SortBy[NSolve[{x424242,y424242}\[Element]Circle[#[[1]],EuclideanDistance[#[[1]],#[[2]]]/1.5]&&
				{x424242,y424242}\[Element]Circle[#[[2]],EuclideanDistance[#[[1]],#[[2]]]/1.5],{x424242,y424242}],#[[2]]&][[2]],
				#[[2]]
			}]&/@Join[linksPointsT,linksPointsS];
			(* Calculate closed boundaries *)
			closedLinesDirT=ClosedLineDirT[#,precisionGood,Abs@#[[2,2]]]&/@(#[[1]]&/@matricesOpenInDirTInv);
			closedLinesDirS=ClosedLineDirS[#,precisionGood]&/@(#[[1]]&/@matricesOpenInDirS);
			closedLinesDirTVertices={};
			closedLinesDirSVertices={};
			If[calculateClosedLinesVertices,
				closedLinesDirTVertices=Flatten[{#[[1,1]],#[[1,-1]]}&/@closedLinesDirT,1];
				closedLinesDirSVertices=Flatten[{#[[1,1]],#[[1,-1]]}&/@closedLinesDirS,1]];
			(* Calculate arrows *)
			arrows=Dynamic[Join[ArrowDirTInv[#,arrowsSize]&/@Flatten[linksMatricesT,1],ArrowDirS[#,arrowsSize]&/@Flatten[linksMatricesS,1]]];
			If[debug,Print[{linksPointsT,linksPointsS,links}]];
			If[debug,Print[{matricesOpenInDirTInv,matricesOpenInDirS}]];
			If[debug,Print["links finished\n"]];

			(* Draw Polygons *)
			If[debug,Print["draw polygons"]];
			(* Initializing a variable inside DynamicModule with Nothing causes an error, see *)
			(* https://mathematica.stackexchange.com/questions/239447/dynamicmodule-error-when-assigning-nothing-to-variable *)
			If[Length[listSplit[[1]]]==0,
				If[Length[listSplit[[2]]]==0,
					funDom={}, (* there are no matrices in the set of reps *)
					funDom={DrawNotInf[#,{precisionGood,precisionBad}]&/@listSplit[[2]]} (* there are no powers of T in the set listSplit[[2]] *)
				],
				If[Length[listSplit[[2]]]==0,
					funDom={DrawInf[#,{precisionGood,precisionBad}]&/@listSplit[[1]]}, (* there are solely powers of T in the set listSplit[[1]] *)
					funDom={
						DrawInf[#,{precisionGood,precisionBad}]&/@listSplit[[1]],
						DrawNotInf[#,{precisionGood,precisionBad}]&/@listSplit[[2]]}
				]
			];
			If[debug,Print[funDom]];
			If[debug,Print["drawing finished\n"]];
		)
		]
	];

	(*****************************************************************)
	(**************************  BTN SET REP  **************************)
	(*****************************************************************)

	(* Button "Promote selected Matrix to Coset Representative" gets pressed *)
	If[btnPromotePMatrixPressed||btnPromoteMMatrixPressed,
		If[debug,Print[{"Button pressed"}]];
		If[btnPromotePMatrixPressed,btnPromotePMatrixPressed=False;
			matrix=pMatrix,
			btnPromoteMMatrixPressed=False;
			matrix=mMatrix];
		If[\[Not]SubgroupEmptyQ[group,args]&&\[Not]MemberQ[reps,matrix],
			(* Update representants *)
			list=ExchangeMatrix[group,args,list,matrix];
			If[debug,Print[{"Button finished"}]];
		]
	];

	(*****************************************************************)
	(****************************  T A B S  ***************************)
	(*****************************************************************)

	pointIsGood=(0<pt[[2]]);

	(* Switching between tabs *)
	If[currentTab=!=selectedTab,
		currentTab=selectedTab;
		If[debug,Print["Switching tabs"]];
		If[debug,Print[{pCurrentPoint,mCurrentPoint,pt}]];
		If[selectedTab==pSelectPoint,
			mCurrentPoint=pt;
			pt=pCurrentPoint;
			pointIsGood=(0<pt[[2]])
		];
		If[selectedTab==mSelectMatrix,
			pCurrentPoint=pt;
			pt=mCurrentPoint;
			pointIsGood=(0<pt[[2]])
		]
	];

	(* Select point calculations *)
	If[selectedTab===pSelectPoint&&pointIsGood,
		pMatrix=FindMatrix[Complex@@pt];
		pMatrixInverse=Inverse[pMatrix];
		If[SubgroupEmptyQ[group,args],
		(
			pEquivMatrix=pMatrix;
			pTransformationMatrix=Id;
			pInverseImage=pt;
			pInverseTransformationMatrix=Id;
		),
		(
			pEquivMatrix=SelectFirst[reps,EquivalentMatricesQ[group,args,pMatrix,#]&];
			pTransformationMatrix=pEquivMatrix.pMatrixInverse;
			pInverseImage=GroupActionReal[pTransformationMatrix,pt];
			pInverseTransformationMatrix=Inverse[pTransformationMatrix];
		)
		]
	];

	(* Select matrix calculations *)
	If[selectedTab===mSelectMatrix&&pointIsGood&&mCurrentPoint!=pt,
		mCurrentPoint=pt;
		mMatrix=FindMatrix[Complex@@pt]];

	(*****************************************************************)
	(************************  G R A P H I C S  ************************)
	(*****************************************************************)

	(* Graphics *)
	Graphics[
	{
	(* Fundamental Domain *)
		If[drawFunDom,Dynamic[{colorFunDom,funDom}]],
	(* Modular Tessellation *)
		Dynamic[Tessellation[t1,t2]],
	(* Tab 1: Matrix *)
		If[selectedTab===pSelectPoint&&pFillMatrix,
			{pColorMatrixImage,DrawMatrixImage[pMatrix,ControlActive[5,precisionGood],False]}],
	(* Tab 1: Equiv. Matrix *)
		If[selectedTab===pSelectPoint&&pFillEquivMatrix,
			{pColorEquivMatrixImage,DrawMatrixImage[pEquivMatrix,ControlActive[5,precisionGood],False]}],
	(* Tab 2: Matrix *)
		If[selectedTab===mSelectMatrix&&mFillMatrix,
			{mColorMatrix,DrawMatrixImage[mMatrix,precisionGood,False]}],
	(* Boundaries *)
		If[drawBoundaries,
			{colorBoundaries,Thickness[boundariesThickness],closedLinesDirT,closedLinesDirS,
			PointSize[0.009],Point[closedLinesDirTVertices],Point[closedLinesDirSVertices]}],
	(* Links *)
		If[drawLinks,
			{colorLinks,Dashing[0.005],links}],
	(* Arrows *)
		If[drawArrows,{colorArrows,arrows}],
	(* Cusps *)
		If[drawCusps,
			If[SubgroupEmptyQ[group,args],
				{colorCusps,Thickness[pPointSize],Line[{{-100,0},{100,0}}],If[drawCuspInf,PointSize[pPointSize],Point[{0,3}]]},
				{colorCusps,PointSize[pPointSize],Point[cuspPoints]}
			]
		],
	(* Tab 1: Inverse Point *)
		If[selectedTab===pSelectPoint&&pDrawInverseImage,
			{pColorPointInverseImage,PointSize[pPointSize],Point[pInverseImage]}],
	(* Tab 1: Point *)
		If[selectedTab===pSelectPoint&&pDrawSelectedPoint,
			{pColorPoint,PointSize[pPointSize],Point[pt]}]
	},
		PlotRange->Dynamic[{{-15+x*zoom^(0.1 zoom),15+x*zoom^(0.1 zoom)}/zoom^(0.1 zoom),
			{-.1-i*zoom^(0.1 zoom),9-i*zoom^(0.1 zoom)}/zoom^(0.1 zoom)}],
		Axes->True,
		PlotLabel->Style[Framed[SubgroupNotation[group]~~
			If[ArgumentsNumber[group]>0,"("~~
			ToString[Row[args,","]]~~
			")",""]~~
			"       Genus:"~~
			ToString[genus]],
			11,Background->LightOrange],
		ImageSize->Dynamic[imagesize]
	],


	(*****************************************************************)
	(*****************************************************************)
	(************************  C O N T R O L S  ************************)
	(*****************************************************************)
	(*****************************************************************)

	(* Top *)
	Dynamic[Control[{{pt,{0,1}},Locator,Appearance->None}]],
	Framed[Grid[
	{
		{
			(* Subgroups Panel *)
			Column[
			{
				Style["\!\(\*SubscriptBox[\(SL\), \(2\)]\)(\!\(\*TemplateBox[{},\n\"Integers\"]\)) and its Congruence Subgroups",16,Bold],
				Grid[
				{
					{
						Tooltip[Style["Subgroup",FontSize->fontSize],"Congruence Subgroup",TooltipDelay->0.5],
						Tooltip[Button[Style["Copy",fontSize],CopyToClipboard[group];Alignment->Center,ImageSize->buttonSize],
							"Copy the congruence subgroup to clipboard",TooltipDelay->0.5],
						Tooltip[Control[
							{
								{group,$listOfSubgroups[[1]], ""},
								Table[$listOfSubgroups[[i]]->SubgroupNotation[$listOfSubgroups[[i]]]~~
									If[ArgumentsNumber[$listOfSubgroups[[i]]]>0,"("~~ArgumentsNotation[$listOfSubgroups[[i]]]~~")",""],
									{i,1,Length[$listOfSubgroups]}],
								ControlType->PopupMenu,
								ImageSize->popupMenuSize+{2 popupMenuSize[[1]],0},
								ContentPadding->False,
								MenuStyle->{FontSize->fontSize},
								FieldSize->9
							}],
							"Select a congruence subgroup",TooltipDelay->0.5],
						SpanFromLeft
					},
					{
						Dynamic[Style[ArgumentsNotation[group],FontSize->fontSize]],
						Tooltip[Button[Style["Copy",fontSize],CopyToClipboard[args];Alignment->Center,ImageSize->buttonSize],
							"Copy the subgroup's arguments to clipboard",TooltipDelay->0.5],
						Dynamic[Grid[
						{
							{
								Dynamic[If[ArgumentsNumber[group]>=1,Identity,Invisible]@Control[
								{
									{n,1,""},
									Table[i,{i,1,If[ArgumentsNumber[group]>=1,ArgumentsMaxValues[group][[1]],1]}],
									ControlType->PopupMenu,
									FieldSize->4,
									ImageSize->If[ArgumentsNumber[group]==3,
										popupMenuSize,
										If[ArgumentsNumber[group]==2,
											popupMenuSize+{popupMenuSize[[1]]/2,0},
											popupMenuSize+{2 popupMenuSize[[1]],0}]],
									ContentPadding->False,
									MenuStyle->{FontSize->fontSizePopupMenu},
									FrameMargins->{{-1,-3},{0,0}}
								}]],
								Dynamic[If[ArgumentsNumber[group]>=2,Identity,Invisible]@Control[
								{
									{m,1,""},
									Table[i,{i,1,If[ArgumentsNumber[group]>=2,ArgumentsMaxValues[group][[2]],1]}],
									ControlType->PopupMenu,
									FieldSize->4,
									ImageSize->If[ArgumentsNumber[group]==3,
										popupMenuSize,
										If[ArgumentsNumber[group]==2,
											popupMenuSize+{popupMenuSize[[1]]/2,0},
											{0,0}]],
									ContentPadding->False,
									MenuStyle->{FontSize->fontSizePopupMenu}
								}]],
								Dynamic[If[ArgumentsNumber[group]>=3,Identity,Invisible]@Control[
									{{p,1,""},Table[i,{i,1,If[ArgumentsNumber[group]>=3,ArgumentsMaxValues[group][[3]],1]}],
									ControlType->PopupMenu,
									FieldSize->4,
									ImageSize->If[ArgumentsNumber[group]==3,
										popupMenuSize,
										{0,0}],
									ContentPadding->False,
									MenuStyle->{FontSize->fontSizePopupMenu}
								}]]
							}
						},Spacings->-1
						]],
						SpanFromLeft
					},
					{
						Tooltip[Style["FunDom",FontSize->fontSize],"Fundamental Domain",TooltipDelay->0.5],
						Tooltip[Button[Style["Copy",fontSize],CopyToClipboard[reps];Alignment->Center,ImageSize->buttonSize],
							"Copy the list of coset representatives to clipoard",TooltipDelay->0.5],
						Tooltip[Control[{{colorFunDom,Orange,""},ColorSetter,ImageSize->buttonSize}],
							"Select a colour for the fundamental domain",TooltipDelay->0.5],
						Tooltip[Control[{{drawFunDom,True,Style["Draw",FontSize->fontSize]},{True,False}}],
							"Display the fundamental domain?",TooltipDelay->0.5]
					},
					{
						Style["Cusps",FontSize->fontSize],
						Tooltip[Button[Style["Copy",fontSize],CopyToClipboard[cusps];Alignment->Center,ImageSize->buttonSize],
							"Copy the cusps to clipoard",TooltipDelay->0.5],
						Tooltip[Control[{{colorCusps,RGBColor[0.7,0,0.9],""},ColorSetter,ImageSize->buttonSize}],
							"Select a colour for the cusps",TooltipDelay->0.5],
						Tooltip[Control[{{drawCusps,False,Style["Draw",FontSize->fontSize]},{True,False}}],
							"Display the cusps?",TooltipDelay->0.5]
					},
					{
						Tooltip[Style["Links",FontSize->fontSize],
							"Links connect equivalent boundaries",TooltipDelay->0.5],
						Tooltip[Button[Style["Copy",fontSize],
								CopyToClipboard[Flatten[#.{{1},{I}}]&/@Join[linksPointsT,linksPointsS]];
								Alignment->Center,
								ImageSize->buttonSize],
							"Copy the links to clipoard",TooltipDelay->0.5],
						Tooltip[Control[{{colorLinks,Black,""},ColorSetter,ImageSize->buttonSize}],
							"Select a colour for the links",TooltipDelay->0.5],
						Tooltip[Control[{{drawLinks,False,Style["Draw",FontSize->fontSize]},{True,False}}],
							"Display the links?",TooltipDelay->0.5]
					},
					{
						Style["Boundaries",FontSize->fontSize],
						Tooltip[Control[{{boundariesThickness,0.005,""},
								{0.0025->"S",0.005->"M",0.0085->"L"},
								PopupMenu,MenuStyle->{FontSize->fontSize},
								ImageSize->popupMenuSize,
								ContentPadding->False,
								FrameMargins->{{0,-5},{0,0}}}],
							"Select a thickness of the closed boundaries",TooltipDelay->0.5],
						Tooltip[Control[{{colorBoundaries,Black,""},ColorSetter,ImageSize->buttonSize}],
							"Select a colour for the closed boundaries",TooltipDelay->0.5],
						Tooltip[Control[{{drawBoundaries,False,Style["Draw",FontSize->fontSize]},{True,False}}],
							"Display the boundaries?",TooltipDelay->0.5]
					},
					{
						Style["Arrows",FontSize->fontSize],
						Tooltip[Control[{{arrowsSize,1,""},
								{0->"S",1->"M",2->"L"},
								PopupMenu,
								MenuStyle->{FontSize->fontSize},
								FieldSize->4,
								ImageSize->popupMenuSize,
								ContentPadding->False,
								FrameMargins->{{0,-5},{0,0}}}],
							"Select a size of the arrows",TooltipDelay->0.5],
						Tooltip[Control[{{colorArrows,Black,""},ColorSetter,ImageSize->buttonSize}],
							"Select a colour for arrows",TooltipDelay->0.5],
						Tooltip[Control[{{drawArrows,False,Style["Draw",FontSize->fontSize]},{True,False}}],
							"Display the arrows?",TooltipDelay->0.5]
					}
				},
				Alignment->{{Left,Center,Center,Center},Center},ItemSize->{{8,8,11,11},Table[3.5,7]}
				]
			},
				ItemSize->{Automatic,{4,3.5*7+4}},
				Alignment->{Center,Center}
			],
			(* Tab view *)
			TabView[
			{
				(* Select Point Tab *)
				{
					pSelectPoint,Tooltip[Style["Select Point",16,Bold],"Select a point in the graphic.",TooltipDelay->0.5]->Grid[
					{
						{
							Tooltip[Style["Selected Point",FontSize->fontSize],
								"Selected point in the graphic",TooltipDelay->0.5],
							Text[If[pointIsGood,Style[Dynamic[Complex@@pt],FontSize->fontSize],""]],
							Tooltip[Button[Style["Copy",fontSize],
									If[pointIsGood,CopyToClipboard[Complex@@pt],MessageDialog[pMsgPointNotGood]];
									Alignment->Center,
									ImageSize->buttonSize],
								"Copy the selected point to the clipboard",TooltipDelay->0.5],
							Tooltip[Control[{{pColorPoint,Red,""},ColorSetter,ImageSize->buttonSize}],
								"Select a colour for the selected point",TooltipDelay->0.5],
							Tooltip[Control[{{pDrawSelectedPoint,True,Style["Draw",FontSize->fontSize]},{True,False}}],
								"Display the selected point?",TooltipDelay->0.5],
							Tooltip[Control[{{pPointSize,0.015,""},
									{0.01->"S",0.015->"M",0.025->"L"},
									PopupMenu,
									MenuStyle->{FontSize->fontSize},
									FieldSize->2,
									ImageSize->popupMenuSize,
									ContentPadding->False,
									FrameMargins->{{0,-5},{0,0}}}],
								"Select the size of the selected point, its inverse image and the cusps",TooltipDelay->0.5]
						},
						{
							Tooltip[Style["Inverse Image",FontSize->fontSize],
								"The selected point's equivalent point inside the standard fundamental domain",TooltipDelay->0.5],
							Text[If[pointIsGood,Style[Dynamic[NumberForm[Complex@@pInverseImage,3]],FontSize->fontSize],"bla"]],
							Tooltip[Button[Style["Copy",fontSize],
									If[pointIsGood,CopyToClipboard[Complex@@pInverseImage],MessageDialog[pMsgPointNotGood]];
									Alignment->Center,
									ImageSize->buttonSize],
								"Copy the inverse image point to the clipboard",TooltipDelay->0.5],
							Tooltip[Control[{{pColorPointInverseImage,Blue,""},ColorSetter,ImageSize->buttonSize}],
								"Select a colour for the inverse image point",TooltipDelay->0.5],
							Tooltip[Control[{{pDrawInverseImage,True,Style["Draw",FontSize->fontSize]},{True,False}}],
								"Display the inverse image point?",TooltipDelay->0.5],
							SpanFromAbove
						},
						{
							Tooltip[Style["Transformation",FontSize->fontSize],
								"The transformation needed to map the inverse image point to the selected point",TooltipDelay->0.5],
							Text[If[pointIsGood,
								Style[Dynamic[
								RightTeeArrow["z",(pInverseTransformationMatrix[[1,1]] "z"+pInverseTransformationMatrix[[1,2]])/(pInverseTransformationMatrix[[2,1]] "z"+pInverseTransformationMatrix[[2,2]])]],
								FontSize->fontSize],""]],
							Tooltip[Button[Style["Copy",fontSize],
									If[pointIsGood,
										CopyToClipboard[With[{z=Unique["z"]},(pInverseTransformationMatrix[[1,1]] z+pInverseTransformationMatrix[[1,2]])/(pInverseTransformationMatrix[[2,1]] z+pInverseTransformationMatrix[[2,2]])]],
										MessageDialog[pMsgPointNotGood]];
									Alignment->Center,
									ImageSize->buttonSize],
								"Copy the transformation to the clipboard",TooltipDelay->0.5]
						},
						{
							Tooltip[Style["Matrix",FontSize->fontSize],
								"The matrix needed to map the inverse image point to the selected point",TooltipDelay->0.5],
							Text[If[pointIsGood,Style[Dynamic[MatrixForm[pMatrix]],FontSize->fontSize-2],""]],
							Tooltip[Button[Style["Copy",fontSize],
									If[pointIsGood,CopyToClipboard[pMatrix],MessageDialog[pMsgPointNotGood]];
									Alignment->Center,
									ImageSize->buttonSize],
								"Copy the matrix to the clipboard",TooltipDelay->0.5],
							Tooltip[Control[{{pColorMatrixImage,RGBColor[0,0.75,0],""},ColorSetter,ImageSize->buttonSize}],
								"Select a colour for the matrix's triangular area",TooltipDelay->0.5],
							Tooltip[Control[{{pFillMatrix,False,Style["Draw",FontSize->fontSize]},{True,False}}],
								"Display the matrix's triangular area?",TooltipDelay->0.5],
							Tooltip[Button[Style["Rep.",fontSize],
									If[pointIsGood,btnPromotePMatrixPressed=True,MessageDialog[pMsgPointNotGood]];
									Alignment->Center,
									ImageSize->buttonSize,
									FrameMargins->{{-2,-5},{0,0}}],
								"Promote the matrix to a coset representative",TooltipDelay->0.5]
						},
						{
							Tooltip[Style["Equiv. matrix",FontSize->fontSize],
								"The matrix's equivalent matrix in the fundamental domain",TooltipDelay->0.5],
							Text[If[pointIsGood,Style[Dynamic[MatrixForm[pEquivMatrix]],FontSize->fontSize-2],""]],
							Tooltip[Button[Style["Copy",fontSize],
									If[pointIsGood,CopyToClipboard[pEquivMatrix],MessageDialog[pMsgPointNotGood]];
									Alignment->Center,
									ImageSize->buttonSize],
								"Copy the equivalent matrix to the clipboard",TooltipDelay->0.5],
							Tooltip[Control[{{pColorEquivMatrixImage,RGBColor[0,0.75,0],""},ColorSetter,ImageSize->buttonSize}],
								"Select a colour for the equivalent matrix's triangular area",TooltipDelay->0.5],
							Tooltip[Control[{{pFillEquivMatrix,False,Style["Draw",FontSize->fontSize]},{True,False}}],
								"Display the equivalent matrix's triangular area?",TooltipDelay->0.5]
						}
					},
						Alignment->{{Left,Center,Center,Center,Center},{Center}},
						ItemSize->{{13,15,8,8,10,8},Table[3.5,5]}
					]
				},
				(* Select Matrix Tab *)
				{
					mSelectMatrix,Tooltip[Style["Select Matrix",16,Bold],"Select a triangular area in the graphic.",TooltipDelay->0.5]->Grid[
					{
						{
							Tooltip[Button[Style["Id",fontSize],mMatrix=Id,Alignment->Center,ImageSize->buttonSizeMMatrix],
								"Reset the matrix to the identity matrix",TooltipDelay->0.5],
							Tooltip[Button[Style["Copy",fontSize],CopyToClipboard[mMatrix],Alignment->Center,ImageSize->buttonSizeMMatrix],
								"Copy the matrix to the clipboard",TooltipDelay->0.5],
							Tooltip[Control[{{mColorMatrix,RGBColor[0,0.75,0],""},ColorSetter,ImageSize->buttonSize}],
								"Select a colour for the matrix's triangular area",TooltipDelay->0.5],
							Tooltip[Control[{{mFillMatrix,True,Style["Draw",FontSize->fontSize]},{True,False}}],
								"Display the matrix's triangular area?",TooltipDelay->0.5],
							Tooltip[Button[Style["Rep.",fontSize],btnPromoteMMatrixPressed=True,Alignment->Center,ImageSize->buttonSizeMMatrix,FrameMargins->{{-2,-5},{0,0}}],
								"Promote the matrix to a coset representative",TooltipDelay->0.5]
						},
						{
							Tooltip[Button[Style["S .",fontSizeMMatrix],mMatrix=S.mMatrix,Alignment->Right,ImageSize->buttonSizeMMatrix,FrameMargins->{{0,10},{-5,-5}}],
								"Left side matrix multiplication by the matrix S",TooltipDelay->0.5],
							Tooltip[InputField[Dynamic[mMatrix,
								(
									mMatrix=
									Which[MatchQ[#,{{_Integer,_Integer},{_Integer,_Integer}}]&&SubgroupMemberQ["SL2Z",{},#],#,True,Id]
								)&],FieldHint->"{{a,b},{c,d}}",ImageSize->{175,25},BaseStyle->{FontSize->fontSize+2.5}],
								"Enter a matrix from \!\(\*SubscriptBox[\(SL\), \(2\)]\)(\!\(\*TemplateBox[{},\n\"Integers\"]\)) here",TooltipDelay->0.5],
							SpanFromLeft,
							SpanFromLeft,
							Tooltip[Button[Style[". S",fontSizeMMatrix],mMatrix=mMatrix.S,Alignment->Left,ImageSize->buttonSizeMMatrix,FrameMargins->{{10,0},{-5,-5}}],
								"Right side matrix multiplication by the matrix S",TooltipDelay->0.5]
						},
						{
							Tooltip[Button[Style["\!\(\*SuperscriptBox[\(S\), \(-1\)]\) .",fontSizeMMatrix],mMatrix=SInv.mMatrix,Alignment->Right,ImageSize->buttonSizeMMatrix,FrameMargins->{{0,10},{-5,-5}}],
								"Left side matrix multiplication by the matrix \!\(\*SuperscriptBox[\(S\), \(-1\)]\)",TooltipDelay->0.5],
							Style[Dynamic[MatrixForm[mMatrix]],FontSize->fontSize+8],
							SpanFromLeft,
							SpanFromLeft,
							Tooltip[Button[Style[". \!\(\*SuperscriptBox[\(S\), \(-1\)]\)",fontSizeMMatrix],mMatrix=mMatrix.SInv,Alignment->Left,ImageSize->buttonSizeMMatrix,FrameMargins->{{10,0},{-5,-5}}],
								"Right side matrix multiplication by the matrix \!\(\*SuperscriptBox[\(S\), \(-1\)]\)",TooltipDelay->0.5]
						},
						{
							Tooltip[Button[Style["T .",fontSizeMMatrix],mMatrix=T.mMatrix,Alignment->Right,ImageSize->buttonSizeMMatrix,FrameMargins->{{0,10},{-5,-5}}],
								"Left side matrix multiplication by the matrix T",TooltipDelay->0.5],
							SpanFromAbove,
							SpanFromBoth,
							SpanFromBoth,
							Tooltip[Button[Style[". T",fontSizeMMatrix],mMatrix=mMatrix.T,Alignment->Left,ImageSize->buttonSizeMMatrix,FrameMargins->{{10,0},{-5,-5}}],
								"Right side matrix multiplication by the matrix T",TooltipDelay->0.5]
						},
						{
							Tooltip[Button[Style["\!\(\*SuperscriptBox[\(T\), \(\(-1\)\(\\\ \)\)]\).",fontSizeMMatrix],mMatrix=TInv.mMatrix,Alignment->Right,ImageSize->buttonSizeMMatrix,FrameMargins->{{0,10},{-5,-5}}],
								"Left side matrix multiplication by the matrix \!\(\*SuperscriptBox[\(T\), \(-1\)]\)",TooltipDelay->0.5],
							Tooltip[Style[Dynamic[STWord[mMatrix]],FontSize->fontSize],
								"ST-Word of the matrix",TooltipDelay->0.5],
							SpanFromLeft,
							SpanFromLeft,
							Tooltip[Button[Style[". \!\(\*SuperscriptBox[\(T\), \(-1\)]\)",fontSizeMMatrix],mMatrix=mMatrix.TInv,Alignment->Left,ImageSize->buttonSizeMMatrix,FrameMargins->{{10,0},{-5,-5}}],
								"Right side matrix multiplication by the matrix \!\(\*SuperscriptBox[\(T\), \(-1\)]\)",TooltipDelay->0.5]
						}
					},
						ItemSize->{{10,10,7,10,10},{5,5,3.5,5,3.5}},Dividers->{{},{False,True}},Alignment->{{Right,Center,Center,Center,Left},Center}]
				}
			},
				Dynamic[selectedTab],ImageSize->Automatic(*{470,280}*),Alignment->{Center,Center}]
		}
	},
		ItemSize->Automatic(*{{50,75},{36,36}}*),
		Alignment->{Center,Center},
		Dividers->{{False,True,False},False}
	],
		RoundingRadius->10,ContentPadding->False,FrameMargins->10
	],

	(* Left *)
	Grid[
	{{
		Control[{{zoom,7,""},1,20,ControlType->VerticalSlider,ImageSize->{Automatic,imagesize/2.5}}],
		Control[{{i,0,""},0,-3,ControlType->VerticalSlider,ImageSize->{Automatic,imagesize/2.5}}]},
		{Style["Zoom",FontSize->fontSize],Style["i",FontSize->fontSize]
	}}
	],

	(* Bottom *)
	{{x,0,Style["x",FontSize->fontSize]},-10,10,ImageSize->imagesize},
	{{imagesize,630,Style["Image size",FontSize->fontSize]},100,1158,ImageSize->imagesize},

	(* Right *)
	Grid[
	{
		{
			Tooltip[Control[{{t1,10,""},Range[0,150],ControlType->VerticalSlider,ImageSize->{Automatic,imagesize/2.5}}],
				"Depth",TooltipDelay->0.5],
			Tooltip[Control[{{t2,5,""},Range[0,15],ControlType->VerticalSlider,ImageSize->{Automatic,imagesize/2.5}}],
				"Range",TooltipDelay->0.5]
		},
		{Tooltip[Style["Tessellation",FontSize->fontSize],"Modular Tessellation",TooltipDelay->0.5],SpanFromLeft}
	}],

	ControlPlacement->Join[
		Table[Top,30],Table[Left,3],Table[Bottom,2],Table[Right,3]],

	ContentSize->Automatic
	]
];


(* ::Chapter:: *)
(*End Package*)


End[]

Protect@@Names["FundamentalDomains`*"];
EndPackage[]
