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

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

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

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

index9FP:=Group(A, Zg*A*Zg^(-1), Zg^(-1)*A*Zg);

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

index9FP:=Group(A, Zg*A*Zg^-1, Zg*B*Zg^-1*B^-1, B*A*B^-1);

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

index9FCA:=FactorCosetAction(GammaBarFP,index9FP);

# List of finite order elements

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

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

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

# Check indexes

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

Index(GammaBarFP,index9FP) = 9;

# Normality

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

# Abelianizations

AbelianInvariants(GammaBarFP) = [3,3];

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

AbelianInvariants(index9FP) = [7];

# 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^index9FCA);

# Check that the index~9 subgroup is torsion free

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

# Check inclusions

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

# Check intersections

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

# Function to calculate fundamental group

FundGp:=function(G,FOList)
  local e1,e2,e3,e4,FCA,GFO,G0;
  e1:=GeneratorsOfGroup(GammaBarFP);
  e2:=Concatenation(e1,List(e1,Inverse));
  Add(e2,One(G));
  e3:=ListX(e2,e2,\*);
  FCA:=FactorCosetAction(GammaBarFP,G);
  GFO:=Concatenation(
    List(FOList,fo->
      List(Difference([1..Index(GammaBarFP,G)],MovedPoints(fo^FCA)),
        j -> fo^First(e3,e->j^(e^FCA) = 1)
      )
    )
  );
  Print("Size(GFO)=",Size(GFO),"\n");
  e1:=GeneratorsOfGroup(G);
  e2:=Concatenation(e1,List(e1,Inverse));
  Add(e2,One(G));
  e3:=ListX(e2,e2,\*);
  e4:=ListX(e2,e3,\*);
  G0:=Group(ListX(e4,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)) = "C7";
StructureDescription(FundGp(index3bFP,FOList)) = "1";
StructureDescription(FundGp(index3cFP,FOList)) = "1";
StructureDescription(FundGp(index3dFP,FOList)) = "1";

# 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(index9FP)) = "C3 x C3";