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

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

# Finitely presented version of $\bar\Gamma$

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

# Important subgroups of $\bar\Gamma$

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

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

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

# List of finite order elements

# These are Donald's
#T1:=C*A;
#T2:=B^(-1)*A*A*C^(-1)*B^(-1);
#T3:=A*C^(-1)*A*B*B;

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

FOList:=[T1,T2,T3];

# Check the indexes of the subgroups

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

# Check normality of the subgroups

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

# Check that normalizers are as expected

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

# Calculate abelianizations of the subgroups

AbelianInvariants(GammaBarFP) = [2,3,3];

AbelianInvariants(index3aFP) = [2,3,13];
AbelianInvariants(index3bFP) = [2,3];
AbelianInvariants(index3cFP) = [2,3];
AbelianInvariants(index3dFP) = [2,3,3];
AbelianInvariants(index3xFP) = [2,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^index3xFCA);

# Check that four of the index 3 groups are torsion free

ForAll(FOList,k->MovedPoints(k^index3aFCA) = [1..3]);
ForAll(FOList,k->MovedPoints(k^index3bFCA) = [1..3]);
ForAll(FOList,k->MovedPoints(k^index3cFCA) = [1..3]);
ForAll(FOList,k->MovedPoints(k^index3dFCA) = [1..3]);
ForAll(FOList,k->MovedPoints(k^index3xFCA) = []);

# 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)) = "C6";

# 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)) = "1";
StructureDescription(AutGp(index3xFP)) = "C3";