FreeGp:=FreeGroup("Z","A","B");
Zf:=FreeGp.1;
Af:=FreeGp.2;
Bf:=FreeGp.3;

Rels:=[
Zf^3,
(Bf*Zf)^3,
Af^2*Zf*Bf*Zf^-1*Af*Zf^-1*Af^3*Bf^-1*Af*Zf,
Af^-1*Zf^-1*Bf*Zf*Af^-1*Zf^-1*Bf*Zf*Af^-1*Bf,
Bf*Af^-1*Bf*Af^-2*Zf*Bf^-2*Zf*Af*Zf^-1*Af*Zf^-1,
Zf*Af*Zf*Af^2*Bf^-1*Af*Zf*Af^-1*Zf*Af^-1*Zf^-1*Bf,
Bf*Zf*Af^-1*Bf^-1*Zf*Af^-1*Zf*Bf^-1*Zf*Af^-1*Zf^-1*Bf,
Bf^-1*Af^-1*Bf*Af^-2*Bf*Zf^-1*Af*Zf^-1*Af*Zf^-1*Bf^-1,
Zf^-1*Af^-1*Bf*Af^-2*Zf^-1*Bf^-2*Zf*Af^-1*Bf^-1*Zf*Af*Bf^-1,
Bf*Af^-1*Zf^-1*Bf*Af*Zf^-1*Bf^-1*Zf*Af^-1*Bf*Zf^-1*Af*Zf^-1*Bf*Af^-1,
Bf^-1*Zf*Af^-1*Zf*Bf^-1*Zf*Af*Zf^-1*Bf^-2*Zf*Af^2*Zf^-1*Bf^-2*Zf*Af
];

# Finitely presented version of $\bar\Gamma$

GammaBarFP:=FreeGp/Rels;
Zg:=GammaBarFP.1;
A:=GammaBarFP.2;
B:=GammaBarFP.3;

# Important subgroups of $\bar\Gamma$

# These are Donald's
#index3aFP:=Group(A, Zg, B*A*B^-1, B^3, B*Zg*B^-1);
#index3bFP:=Group(A, B^-1*Zg^-1, Zg^-1*B^-1, B*A*B^-1);
#index3cFP:=Group(A, Zg*B^-1, Zg^-1*B);
#index3dFP:=Group(A, B, Zg*A*Zg^-1);
#index9aFP:=Group(A, B*A*B^-1, B^3, Zg*A*Zg^-1, Zg*B*Zg^-1*B^-1, B^-1*A*B);
#index9bFP:=Group(A, B*Zg*A^-1*B^-1, B*A^-1*B, B*Zg^-1*B, Zg*B*Zg^-1);

# These are Tim's
index3aFP:=Group(Zg, A, B*Zg*B^-1, B*A*B^-1);
index3bFP:=Group(A, B*Zg, Zg*B, Zg*A*Zg^-1);
index3cFP:=Group(A, Zg*B^-1, Zg^-1*B);
index3dFP:=Group(A, B, Zg*A*Zg^-1);
index9aFP:=Group(A, Zg*A*Zg^-1, Zg*B*Zg^-1*B^-1, B*A*B^-1, B^3, Zg^-1*A*Zg);
index9bFP:=Group(A, Zg*A*Zg, Zg*B*Zg^-1, B*Zg*A^-1*B^-1);

index3aFCA:=FactorCosetAction(GammaBarFP,index3aFP);
index3bFCA:=FactorCosetAction(GammaBarFP,index3bFP);
index3cFCA:=FactorCosetAction(GammaBarFP,index3cFP);
index3dFCA:=FactorCosetAction(GammaBarFP,index3dFP);

index9aFCA:=FactorCosetAction(GammaBarFP,index9aFP);
index9bFCA:=FactorCosetAction(GammaBarFP,index9bFP);

# List of finite order elements

# These are Donald's
#T1:=Zg;
#T2:=B*Zg;
#T3:=Zg*A*A*B^-1*A;
#T4:=B^-1*Zg*A^-1*Zg*B^-1*Zg;

# These are Tim's
T1:=Zg;
T2:=B*Zg;
T3:=Zg*A^2*B^-1*A;
T4:=B^-1*Zg*A^-1*Zg*B^-1*Zg;

FOList:=[T1,T2,T3,T4];

# Check the indexes of the subgroups

Index(GammaBarFP,index3aFP) = 3;
Index(GammaBarFP,index3bFP) = 3;
Index(GammaBarFP,index3cFP) = 3;
Index(GammaBarFP,index3dFP) = 3;

Index(GammaBarFP,index9aFP) = 9;
Index(GammaBarFP,index9bFP) = 9;

# Check normality of the subgroups

IsNormal(GammaBarFP,index3aFP);
IsNormal(GammaBarFP,index3bFP);
IsNormal(GammaBarFP,index3cFP);
IsNormal(GammaBarFP,index3dFP);

IsNormal(GammaBarFP,index9aFP);
not IsNormal(GammaBarFP,index9bFP);

# Check that normalizers are as expected

Normalizer(GammaBarFP,index3aFP) = GammaBarFP;
Normalizer(GammaBarFP,index3bFP) = GammaBarFP;
Normalizer(GammaBarFP,index3cFP) = GammaBarFP;
Normalizer(GammaBarFP,index3dFP) = GammaBarFP;

Normalizer(GammaBarFP,index9aFP) = GammaBarFP;
Normalizer(GammaBarFP,index9bFP) = index9bFP;

# Check inclusions

IsSubgroup(index3aFP,index9aFP);
IsSubgroup(index3bFP,index9aFP);
IsSubgroup(index3cFP,index9aFP);
IsSubgroup(index3dFP,index9aFP);

not IsSubgroup(index3aFP,index9bFP);
not IsSubgroup(index3bFP,index9bFP);
not IsSubgroup(index3cFP,index9bFP);
not IsSubgroup(index3dFP,index9bFP);

# Check intersections

Intersection(index3aFP,index3bFP) = index9aFP;
Intersection(index3aFP,index3cFP) = index9aFP;
Intersection(index3aFP,index3dFP) = index9aFP;
Intersection(index3bFP,index3cFP) = index9aFP;
Intersection(index3bFP,index3dFP) = index9aFP;
Intersection(index3cFP,index3dFP) = index9aFP;

# Calculate abelianizations of the subgroups

AbelianInvariants(GammaBarFP) = [3,3];

AbelianInvariants(index3aFP) = [3,13];
AbelianInvariants(index3bFP) = [3];
AbelianInvariants(index3cFP) = [3];
AbelianInvariants(index3dFP) = [2,2,3];

AbelianInvariants(index9aFP) = [2,2,13];
AbelianInvariants(index9bFP) = [2,3,3];

# Check which elements of finite order belong to which conjugates of
# which finite index subgroups 

List(FOList,k->k^index3aFCA);
List(FOList,k->k^index3bFCA);
List(FOList,k->k^index3cFCA);
List(FOList,k->k^index3dFCA);

List(FOList,k->k^index9aFCA);
List(FOList,k->k^index9bFCA);

# Check that the index 9 groups are torsion free

ForAll(FOList,k->MovedPoints(k^index9aFCA) = [1..9]);
ForAll(FOList,k->MovedPoints(k^index9bFCA) = [1..9]);

# Function to calculate fundamental group

FundGp:=function(G,FOList)
  local e1,e2,e3,GFO,G0,G0FCA;
  e1:=GeneratorsOfGroup(GammaBarFP);
  e2:=Concatenation(e1,List(e1,Inverse));
  Add(e2,One(G));
  e3:=ListX(e2,e2,\*);
  GFO:=Filtered(FOList,fo->fo in G);
  G0:=Group(ListX(e3,GFO,function(elt,fo) return fo^elt; end));
  Print(ForAll(GFO,fo->fo in G0),"\n");
  Print(IsNormal(G,G0),"\n");
  return FactorGroup(G,G0);
end;

StructureDescription(FundGp(GammaBarFP,FOList)) = "1";
StructureDescription(FundGp(index3aFP,FOList)) = "C13";
StructureDescription(FundGp(index3bFP,FOList)) = "1";
StructureDescription(FundGp(index3cFP,FOList)) = "1";
StructureDescription(FundGp(index3dFP,FOList)) = "Q8";
AbelianInvariants(FundGp(index3dFP,FOList)) = [2,2];

# This calculates the automorphism group of the corresponding surface
AutGp:=G->FactorGroup(Normalizer(GammaBarFP,G),G);

StructureDescription(AutGp(index3aFP)) = "C3";
StructureDescription(AutGp(index3bFP)) = "C3";
StructureDescription(AutGp(index3cFP)) = "C3";
StructureDescription(AutGp(index3dFP)) = "C3";
StructureDescription(AutGp(index9aFP)) = "C3 x C3";
StructureDescription(AutGp(index9bFP)) = "1";