F:=FreeGroup("w","x", "y", "z");;
x:=F.2;; y:=F.3;; z:=F.4;;w:=F.1;;
#Gives the orbifold with volume v1/6
#for this one index of link group should be 24
#G:=F/[x^3,y^2,z^2,(y*z)^4,(z*x)^2,(x*y)^4, w];;

#Gives an orbifold with volume v1/2 that covers the orbifold of volume v1/4
#for this one index of link group should be 8
#G:=F/[x^4,y^2,z^2,(y*z^(-1))^4,(z*x^(-1))^2,(x*y^(-1))^4, w];;

#Gives the orbifold with volume v1/2 that covers the orbifold of volume v1/2
#for this one index of link group should be 16
# the relation should either be w*x*w*z*y^-1 or w*x*w*y*z^(-1)
G:=F/[x^4,y^2,z^2,w^2,(y*z^(-1))^4,(z*x^(-1))^2,(x*y^(-1))^4, w*x*w*z*y^(-1),w*y*w*x*z^(-1),w*z*w*z ];;

#the index of the link group in G
IndexLinkInG :=16;

#the index of the Orbifold With Torus Cusp in G 
IndexOWTCInG :=IndexLinkInG/4;


SetInfoLevel( InfoFpGroup, 2 );
lis:=LowIndexSubgroupsFpGroup(G, IndexLinkInG);;

possibleLinks := [];;
possibleLinksNicePresentation := [];;

SetInfoLevel( InfoFpGroup, 0 );
for counter in [1..Size(lis)] do
#indexForPossibleLinks = 0;;
#Print("counter is ", counter, ".\n");
        if Index(G,lis[counter])=IndexLinkInG then 
                #Print("Index is ",Index(G,lis[counter]), ".\n");
                hom:=MaximalAbelianQuotient(lis[counter]);;
                if AbelianInvariants(Image(hom))=[0,0] then
                #Print("Counter is ",counter, "\n");
                #passes the result to a finitely presented group
                fp:=Image(IsomorphismFpGroup(lis[counter]));;
                #passes result to a presentation group in Tietze prese
                p:=PresentationFpGroup(fp);
                #passes back to a finitely presented group
                fp:=FpGroupPresentation(p);
                Add(possibleLinks,lis[counter]);;
                Add(possibleLinksNicePresentation,fp);;

                fi;
        fi;
od;

Print("There are ", Size(possibleLinks), " possible links. \n");
Print(possibleLinks, "\n");
#Print(IsGroup(possibleLinks[1]));

size:=Size(possibleLinks);

for counter in [1..(size)] do

        Print("Group ", counter, " is ", possibleLinks[counter], "\n");
        Print(RelatorsOfFpGroup(possibleLinksNicePresentation[counter]), "\n");

od;

for counter1 in [1..(size)] do
        Print("Group ", counter1, " is ", possibleLinks[counter], "\n");
        SetInfoLevel( InfoFpGroup, 2);
        #gives the subgroups of index IndexOWTCInG that contain our link groups
        lis2:=LowIndexSubgroupsFpGroup(G, possibleLinks[counter1], IndexOWTCInG);
        
        for counter2 in [1..(Size(lis2))] do
        if Index(G,lis2[counter2])=IndexOWTCInG then
        Print("Group is ", lis2[counter2],"\n");
        hom:=MaximalAbelianQuotient(lis2[counter2]);
        Print("Abelization is ",  AbelianInvariants(Image(hom)), "\n");
        fi;
        od;
od;





#for counter1 in [1..(size-1)] do

 #       for counter2 in [(counter1+1)..size] do
 #               Print("Group ", counter1, " and ", counter2, " have isomorphism ", IsomorphismGroups(possibleLinks[counter1], possibleLinks[counter2]));
                #Print("counter 1 ", counter1);
                #Print("counter 2 ", counter2, "\n");       
#        od;

#od;



## next need to find out which subgroups are isomorphic to the Whitehead Sister
## then compute all subgroups of index 6 that are supergroups of
######### the whitehead sister 

quit;