# Andriy Zhugayevych, azh@ukr.net
# Point Groups package
# created 21.09.2004


PointGroups:=module() option package;

export ModuleLoad, Setup,
       angles2xyz, EulerParam, pAngle, RepDL, ChDL, # Rotation group
       PG_c, PG_e, PG_u2h, PG_u2d, PG_u2K, PG_u3K, PG_sh, PG_sv, PG_svK, PG_u5Y, # Symmetry elements
       GF_C, GF_D, GF_Dh, GF_TO, GF_Ch, GF_T, GF_Th, GF_O, GF_Oh, GF_Y, GF_Yh, # Generating functions of groups
       GFC_C, GFC_D, GFC_T, GFC_O, GFC_Y, NC_C, NC_D,  # Generating functions of conjugacy classes
       genR_C, genR_D, genR_T, genR_O, genR_Y, # Generators of irreducible representations
       tbPG; # Table of point groups (Group Records)

local  version, SimRad, AntiSym;


ModuleLoad:=proc()
local v,u,w;
version:=20191201;
# FiniteGroups package is required
if not(member('FiniteGroups',packages())) then try with(FiniteGroups) catch: error "Cannot find the required FiniteGroups package: %1",lastexception end end;
if (FiniteGroups[Setup]('version')<20191117) then WARNING("Update FiniteGroups package to version from 17.11.2019 at least") end;
# Default settings
FiniteGroups[Setup]('setting'="Matrix");
NULL
end:


Setup:=proc({printout::boolean:=false})
local knowngroups,sq;
knowngroups:=[indices(tbPG,nolist)];
sq:=ProcessSetupArgs([_rest],[],['version','knowngroups']);
sq
end:


##########################################################
### Auxiliary procedures

SimRad:=e->signum(e)*sqrt(expand(rationalize((radnormal(e)^2)))):

AntiSym:=proc(M::Matrix)
local n,T,i;
n:=[op(1,M)][1];
T:=Matrix(n,(i,j)->`if`(i>j,0,1)):
for i from 1 to n-1 do T[i+1,i]:=-i end:
for i from 1 to n do T[1..-1,i]:=LinearAlgebra[Normalize](T[1..-1,i],2) end:
FiniteGroups[FGS](LinearAlgebra[Transpose](T).M.T)[1..n-1,1..n-1]
end:



##########################################################
### Rotation group

angles2xyz:=psi->simplify(factor(subs(r^2=x^2+y^2+z^2,simplify(algsubs(sin(theta)^2=(x^2+y^2)/r^2,subs(cos(theta)=z/r,subs(sin(phi)=y/sin(theta)/r,cos(phi)=x/sin(theta)/r,expand(evalc(psi),trig)))))))):

EulerParam:=proc(g)
local alpha1,n,theta,phi,p,a,b,c;
if type(g,string) then
 if type(op(sscanf(g,"c%d")),posint) then
  n:=op(sscanf(g,"c%d"));
  if (g=cat("c",n)) then return 1,0,0,2*Pi/n end
 elif type(op(sscanf(g,"ci%d")),posint) then
  n:=op(sscanf(g,"ci%d"));
  if (g=cat("ci",n)) then return -1,0,0,2*Pi/n end
 elif (g="u2h") then return 1,Pi/2,Pi,-Pi/2
 elif (g="u2K") then return 1,Pi,Pi/2,0
 elif (g="u3K") then return 1,Pi/2,Pi/2,0
 elif (g="ui3K") then return -1,Pi/2,Pi/2,0
 elif (g="u5Y") then return 1,arctan((sqrt(5)+1)/2),2*Pi/5,-arctan((sqrt(5)+1)/2)
 elif (g="ui5Y") then return -1,arctan((sqrt(5)+1)/2),2*Pi/5,-arctan((sqrt(5)+1)/2)
 elif (g="i") then return -1,0,0,0
 elif (g="sh") then return -1,0,0,Pi
 elif (g="sv") then return -1,Pi/2,Pi,-Pi/2
 elif (g="svK") then return -1,Pi,Pi/2,0  end
elif type(g,Matrix) then
 p,alpha1,n:=RotationParam(g);
 theta:=simplify(arccos(n[3]/sqrt(n[1]^2+n[2]^2+n[3]^2)));
 phi:=simplify(arctan(n[2],n[1]));
 return p,Rot2EulerY(alpha1,theta,phi)
else
 error "Unrecognized argument",g  end;
error "Cannot determine Euler parameters"
end:


pAngle:=proc(g)
local p,a,n,theta,phi,e;
if type(g,string) then
 if type(op(sscanf(g,"c%d")),posint) then
  e:=sscanf(g,"c%d^%d");
  if (nops(e)=1) then return 1,2*Pi/e[1] else return 1,2*Pi*e[2]/e[1] end
 elif type(op(sscanf(g,"ci%d")),posint) then
  e:=sscanf(g,"ci%d^%d");
  if (nops(e)=1) then return -1,2*Pi/e[1] else return (-1)^e[2],2*Pi*e[2]/e[1] end
 elif type(op(sscanf(g,"ic%d")),posint) then
  e:=sscanf(g,"ic%d^%d");
  if (nops(e)=1) then return -1,2*Pi/e[1] else return -1,2*Pi*e[2]/e[1] end
 elif (g[1..2]="u2") then return 1,Pi
 elif (g[1..2]="u3") then return 1,2*Pi/3
 elif (g[1..3]="ui3") then return -1,2*Pi/3
 elif (g="e") then return 1,0
 elif (g="i") then return -1,0
 elif (g[1]="s") then return -1,Pi  end
elif type(g,Matrix) then
 return op([RotationParam(g)][1..2])
else
 error "Unrecognized argument",g  end;
error "Cannot determine rotation angle"
end:


RepDL:=proc(g,L::nonnegint)
local det,a,b,c,p;
if (nargs=2) then
 p:=(-1)^L
else
 p:=args[3];
 if (p<>1 and p<>-1) then error "Incorrect p",p  end  end;
det,a,b,c:=EulerParam(g);
p:=`if`(p=-1 and det=-1,-1,1);
Matrix(2*L+1,(i,j)->simplify(convert(p*WignerD(L,L-i+1,L-j+1,a,b,c),radical)))
end:


ChDL:=proc(g,L::nonnegint)
local det,a,p;
if (nargs=2) then
 p:=(-1)^L
else
 p:=args[3];
 if (p<>1 and p<>-1) then error "Incorrect p",p  end  end;
det,a:=pAngle(g);
p:=`if`(p=-1 and det=-1,-1,1);
p*`if`(type(a/2/Pi,integer), 2*L+1, simplify(convert(sin((L+1/2)*a)/sin(a/2),radical)) )
end:



##########################################################
### Point groups

# Symmetry elements
PG_c:=n->map(SimRad,map(convert,RotationM(2*Pi/n,<0,0,1>),radical)):
PG_e:=PG_c(1):
PG_u2h:=RotationM(Pi,<1,0,0>):
PG_u2d:=RotationM(Pi,<1,1,0>):
PG_u2K:=RotationM(Pi,<1,0,1>):
PG_u3K:=RotationM(2*Pi/3,<1,1,1>):
PG_sh:=ReflectionM(<0,0,1>):
PG_sv:=ReflectionM(<1,0,0>):
PG_svK:=ReflectionM(<1,0,1>):
PG_u5Y:=map(SimRad,map(convert,RotationM(2*Pi/5,<sqrt(5)+1,2,0>),radical)):

# Universal generating functions
GF_C:=n->proc(gen::list) FiniteGroups[FG_cyclic](gen[1],n) end:
GF_D:=n->proc(gen::list) FiniteGroups[FG_prod]( FiniteGroups[FG_cyclic](gen[1],n), FiniteGroups[FG_cyclic](gen[2],2) ) end:
GF_Dh:=n->proc(gen::list) FiniteGroups[FG_prod](FiniteGroups[FG_prod]( FiniteGroups[FG_cyclic](gen[1],n), FiniteGroups[FG_cyclic](gen[2],2) ), FiniteGroups[FG_cyclic](gen[3],2) ) end:
GF_TO:=(n,m)->proc(gen::list) FiniteGroups[FG_prod](FiniteGroups[FG_prod]( FiniteGroups[FG_cyclic](gen[1],n), FiniteGroups[FG_cyclic](gen[2],2) ), FiniteGroups[FG_cyclic](foldl(FiniteGroups[FG_mul],gen[1],gen[2],gen[1],gen[1]),m) ) end:

# Specific generating functions of nonisomorphic groups
GF_Ch:=GF_D:
GF_T :=GF_TO(3,2):
GF_Th:=GF_TO(6,2):
GF_O :=GF_TO(4,3):
GF_Oh:=GF_TO(4,6):
GF_Y :=proc(gen::list)  local a,b,C5,S,g;
 a,b:=op(gen);
 C5:=FiniteGroups[FG_cyclic](a,5);
 S:=[C5[1],seq(FiniteGroups[FGS](g.b),g=C5)];
 FiniteGroups[FG_prod](S, FiniteGroups[FG_prod]( C5, FiniteGroups[FG_cyclic](FiniteGroups[FGS](a^2 .b.a^3 .b.a^2 .b),2) )) end:
GF_Yh:=proc(gen::list)  local Y,i;
 i:=FiniteGroups[FGS](gen[1]^5);
 Y:=GF_Y([i.gen[1],gen[2]]);
 [op(Y),seq(i.g,g=Y)] end:

# Generating functions of conjugacy classes
GFC_C:=n->proc(gen::list)  local a;  a:=gen[1];
 if type(n,odd) then
  map(FiniteGroups[FGS],[seq(a^p,p=0),seq(op([a^p,a^(-p)]),p=1..(n-1)/2)])
 else
  map(FiniteGroups[FGS],[seq(a^p,p=[0,n/2]),seq(op([a^p,a^(-p)]),p=1..n/2-1)]) end end:
GFC_D:=n->proc(gen::list)  local a,b;  a,b:=op(gen[1..2]);
 if type(n,odd) then
  map(FiniteGroups[FGS],[seq(a^p,p=0),b,seq(a^p,p=1..(n-1)/2)])
 else
  map(FiniteGroups[FGS],[seq(a^p,p=[0,n/2]),b,a.b,seq(a^p,p=1..n/2-1)]) end end:
GFC_T:=proc(gen::list)
  map(FiniteGroups[FGS],[seq(gen[1]^p,p=0..2),gen[2]]) end:
GFC_O:=proc(gen::list)
  map(FiniteGroups[FGS],[seq((gen[1].gen[2])^p,p=0..1),gen[1]^2,gen[1],gen[2]]) end:
GFC_Y:=proc(gen::list)
  map(FiniteGroups[FGS],[seq(gen[1]^p,p=0..2),gen[2],gen[1].gen[2]]) end:

# Number of elements in conjugacy classes
NC_C:=n->[1$n]:
NC_D:=n->`if`(type(n,odd),[1,n,2$((n-1)/2)],[1,1,n/2,n/2,2$(n/2-1)]):

# Basic generators of irreducible representations
genR_C:=n->`if`(type(n,odd),
 [[<<1>>],
   seq(op([[<<exp(I*2*Pi*p/n)>>],[<<exp(-I*2*Pi*p/n)>>]]),p=1..(n-1)/2)],
 [[<<1>>],[<<-1>>],
   seq(op([[<<exp(I*2*Pi*p/n)>>],[<<exp(-I*2*Pi*p/n)>>]]),p=1..n/2-1)]):
genR_D:=n->`if`(type(n,odd),
 [[<<1>>,<<1>>],[<<1>>,<<-1>>],
   seq([<<exp(I*2*Pi*p/n),0>|<0,exp(-I*2*Pi*p/n)>>,
   <<0,1>|<1,0>>],p=1..(n-1)/2)],
 [[<<1>>,<<1>>],[<<1>>,<<-1>>],[<<-1>>,<<1>>],[<<-1>>,<<-1>>],
   seq([<<exp(I*2*Pi*p/n),0>|<0,exp(-I*2*Pi*p/n)>>,
   <<0,1>|<1,0>>],p=1..n/2-1)]):
genR_T:=[
 [<<1>>, <<1>>],
 [<<exp( 2*Pi*I/3)>>, <<1>>],
 [<<exp(-2*Pi*I/3)>>, <<1>>],
 [ PG_u3K, PG_c(2)]]:
genR_O:=[
 [<<1>>, <<1>>], [<<-1>>, <<-1>>],
 [<<0|exp(2*Pi*I/3)>,<exp(-2*Pi*I/3)|0>>, <<0|1>,<1|0>>],
 [ PG_c(4), PG_u2K],
 [-PG_c(4), PG_svK]]:
genR_Y:=[
 [<<1>>,<<1>>],
 [PG_u5Y,PG_u2h],
 [FiniteGroups[FGS]((PG_u2h.PG_u5Y^2)^2),PG_u2h],
 [AntiSym(Matrix(5,{(1,5)=1,(5,4)=1,(4,3)=1,(3,2)=1,(2,1)=1})),
  AntiSym(Matrix(5,{(1,1)=1,(2,3)=1,(3,2)=1,(4,5)=1,(5,4)=1}))],
 [AntiSym(Matrix(6,{(1,1)=1,(2,6)=1,(6,5)=1,(5,4)=1,(4,3)=1,(3,2)=1})),
  AntiSym(Matrix(6,{(1,2)=1,(2,1)=1,(3,6)=1,(6,3)=1,(4,4)=1,(5,5)=1}))]]:

##########################################################
# Table of point groups (Group Records)

tbPG:=table():

tbPG["C"]:=n->Record(
 'name'=cat("C",n),
 'namegen'=[cat("c",n)],
 'gen'=[PG_c(n)],
 'GF'=GF_C(n),
 'nameC'=`if`(type(n,odd),
   ["e",seq(op([cat("c",n,"^",p),cat("c",n,"^",-p)]),p=1..(n-1)/2)],
   ["e","c2",seq(op([cat("c",n,"^",p),cat("c",n,"^",-p)]),p=1..n/2-1)]),
 'GFC'=GFC_C(n),
 'NC'=NC_C(n),
 'nameR'=`if`(type(n,odd),
   ["A",seq(cat("E",p)$2,p=1..(n-1)/2)],
   ["A","B",seq(cat("E",p)$2,p=1..n/2-1)]),
 'genR'=genR_C(n)):

tbPG["Ci"]:=n->Record(
 'name'=cat("Ci",n),
 'namegen'=[cat("ci",n)],
 'gen'=[-PG_c(n)],
 'GF'=GF_C(`if`(type(n,even),n,2*n)),
 'nameC'=`if`(type(n,even),
   ["e","c2",seq(op([cat("ci",n,"^",p),cat("ci",n,"^",-p)]),p=1..n/2-1)],
   ["e","i",seq(op([cat("ci",n,"^",p),cat("ci",n,"^",-p)]),p=1..n)]),
 'GFC'=GFC_C(`if`(type(n,even),n,2*n)),
 'NC'=NC_C(`if`(type(n,even),n,2*n)),
 'nameR'=`if`(type(n,even),
   ["A","B",seq(cat("E",p)$2,p=1..n/2-1)],
   ["Ag","Au",seq(op([cat("E",p,"u"),cat("E",p,"u"),cat("E",p,"g"),cat("E",p,"g")]),p=1..(n-1)/2)]),
 'genR'=genR_C(`if`(type(n,even),n,2*n))):

tbPG["Ch"]:=n->`if`(type(n,odd),sprintf("See Ci[%d]",2*n),Record(
 'name'=cat("Ch",n),
 'namegen'=[cat("c",n),"i"],
 'gen'=[PG_c(n),-PG_e],   ### strictly speaking [PG_c(n),PG_sh]
 'GF'=GF_Ch(n),
 'nameC'=
   ["e","c2",seq(op([cat("c" ,n,"^",p),cat("c" ,n,"^",-p)]),p=1..n/2-1),
    "i","sh",seq(op([cat("ic",n,"^",p),cat("ic",n,"^",-p)]),p=1..n/2-1)],
 'GFC'=proc(gen::list) local ls;
   ls:=GFC_C(n)(gen);
   [op(ls),seq(gen[2].M,M=ls)] end,
 'NC'=[op(NC_C(n)),op(NC_C(n))],
 'nameR'=[
   "Ag","Bg",seq(cat("E",p,"g")$2,p=1..n/2-1),
   "Au","Bu",seq(cat("E",p,"u")$2,p=1..n/2-1)],
 'genR'=[seq([e[1],seq( e[1]^p,p=0)],e=genR_C(n)),
         seq([e[1],seq(-e[1]^p,p=0)],e=genR_C(n))])):

tbPG["D"]:=n->Record(
 'name'=cat("D",n),
 'namegen'=[cat("c",n),"u2h"],
 'gen'=[PG_c(n),PG_u2h],
 'GF'=GF_D(n),
 'nameC'=`if`(type(n,odd),
   ["e","u2",seq(cat("c",n,"^",p),p=1..(n-1)/2)],
   ["e","c2","u2","u2'",seq(cat("c",n,"^",p),p=1..(n-1)/2)]),
 'GFC'=GFC_D(n),
 'NC'=NC_D(n),
 'nameR'=`if`(type(n,odd),
   ["A1","A2",seq(cat("E",p),p=1..(n-1)/2)],
   ["A1","A2","B1","B2",seq(cat("E",p),p=1..n/2-1)]),
 'genR'=genR_D(n)):

tbPG["Di"]:=n->Record(
 'name'=cat("Di",n),
 'namegen'=[cat("ci",n),"u2h"],
 'gen'=[-PG_c(n),PG_u2h],
 'GF'=GF_D(`if`(type(n,even),n,2*n)),
 'nameC'=`if`(type(n,even),`if`(type(n/2,even),
   ["e","c2","u2","u2'",seq(cat("ci",n,"^",p),p=1..(n-1)/2)],
   ["e","sh","u2","sv'",seq(cat("ci",n,"^",p),p=1..(n-1)/2)]),
   ["e","i","u2","sv'",seq(cat("ci",n,"^",p),p=1..n)]),
 'GFC'=GFC_D(`if`(type(n,even),n,2*n)),
 'NC'=NC_D(`if`(type(n,even),n,2*n)),
 'nameR'=`if`(type(n,even),`if`(type(n/2,even),
   ["A1","A2","B1","B2",seq(cat("E",p),p=1..n/2-1)],
   ["A1'","A2'","A1''","A2''",seq(cat("E",p),p=1..n/2-1)]),
   ["A1g","A2g","A1u","A2u",seq(op([cat("E",p,"u"),cat("E",p,"g")]),p=1..(n-1)/2)]),
 'genR'=genR_D(`if`(type(n,even),n,2*n))):

tbPG["Dh"]:=n->`if`(type(n,odd),sprintf("See Di[%d]",2*n),Record(
 'name'=cat("Dh",n),
 'namegen'=[cat("c",n),"u2h","i"],
 'gen'=[PG_c(n),PG_u2h,-PG_e],    ### strictly speaking [PG_c(n),PG_u2h,PG_sh]
 'GF'=GF_Dh(n),
 'nameC'=
   ["e","c2","u2","u2'",seq(cat("c" ,n,"^",p),p=1..(n-1)/2),
    "i","sh","sv","sv'",seq(cat("ic",n,"^",p),p=1..(n-1)/2)],
 'GFC'=proc(gen::list) local ls;
   ls:=GFC_D(n)(gen);
   [op(ls),seq(gen[3].M,M=ls)] end,
 'NC'=[op(NC_D(n)),op(NC_D(n))],
 'nameR'=[
   "A1g","A2g","B1g","B2g",seq(cat("E",p,"g"),p=1..n/2-1),
   "A1u","A2u","B1u","B2u",seq(cat("E",p,"u"),p=1..n/2-1)],
 'genR'=[seq([e[1],e[2],seq( e[1]^p,p=0)],e=genR_D(n)),
         seq([e[1],e[2],seq(-e[1]^p,p=0)],e=genR_D(n))])):

tbPG["Cv"]:=n->Record(
 'name'=cat("Cv",n),
 'namegen'=[cat("c",n),"sv"],
 'gen'=[PG_c(n),PG_sv],
 'GF'=GF_D(n),
 'nameC'=`if`(type(n,odd),
   ["e","sv",seq(cat("c",n,"^",p),p=1..(n-1)/2)],
   ["e","c2","sv","sv'",seq(cat("c",n,"^",p),p=1..(n-1)/2)]),
 'GFC'=GFC_D(n),
 'NC'=NC_D(n),
 'nameR'=`if`(type(n,odd),
   ["A1","A2",seq(cat("E",p),p=1..(n-1)/2)],
   ["A1","A2","B1","B2",seq(cat("E",p),p=1..n/2-1)]),
 'genR'=genR_D(n)):

tbPG["T"]:=Record(
 'name'="T",
 'namegen'=["u3K","c2"],
 'gen'=[PG_u3K,PG_c(2)],
 'GF'=GF_T,
 'nameC'=["e","c3","c3^2","u2"],
 'GFC'=GFC_T,
 'NC'=[1,4,4,3],
 'nameR'=["A","E","E","F"],
 'genR'=genR_T):

tbPG["Th"]:=Record(
 'name'="Th",
 'namegen'=["ui3K","c2"],
 'gen'=[-PG_u3K,PG_c(2)],
 'GF'=GF_Th,
 'nameC'=["e","c3","c3^2","u2","i","ci3","ci3^5","sv"],
 'GFC'=proc(gen::list)  local i,ls;
   i:=FiniteGroups[FGS](gen[1]^3);
   ls:=GFC_T([i.gen[1],gen[2]]);
   [op(ls),seq(i.M,M=ls)] end,
 'NC'=[1,4,4,3,1,4,4,3],
 'nameR'=["Ag","Eg","Eg","Fg","Au","Eu","Eu","Fu"],
 'genR'=[op(genR_T),seq([-e[1],e[2]],e=genR_T)]):

tbPG["Td"]:=Record(
 'name'="Td",
 'namegen'=["ci4","svK"],
 'gen'=[-PG_c(4), PG_svK],
 'GF'=GF_O,
 'nameC'=["e","c3","u2","ci4","sv"],
 'GFC'=GFC_O,
 'NC'=[1,8,3,6,6],
 'nameR'=["A1","A2","E","F1","F2"],
 'genR'=genR_O):

tbPG["O"]:=Record(
 'name'="O",
 'namegen'=["c4","u2K"],
 'gen'=[PG_c(4),PG_u2K],
 'GF'=GF_O,
 'nameC'=["e","c3","c2","c4","u2"],
 'GFC'=GFC_O,
 'NC'=[1,8,3,6,6],
 'nameR'=["A1","A2","E","F1","F2"],
 'genR'=genR_O):

tbPG["Oh"]:=Record(
 'name'="Oh",
 'namegen'=["ci4","u2K"],
 'gen'=[-PG_c(4),PG_u2K],
 'GF'=GF_Oh,
 'nameC'=["e","c3","c2","c4","u2","i","ci3","sh","ci4","sv"],
 'GFC'=proc(gen::list)  local i,ls;
   i:=FiniteGroups[FGS]((gen[1].gen[2])^3);
   ls:=GFC_O([i.gen[1],gen[2]]);
   [op(ls),seq(i.M,M=ls)] end,
 'NC'=[1,8,3,6,6,1,8,3,6,6],
 'nameR'=["A1g","A2g","Eg","F1g","F2g","A1u","A2u","Eu","F1u","F2u"],
 'genR'=[op(genR_O),seq([-e[1],e[2]],e=genR_O)]):

tbPG["Y"]:=Record(
 'name'="Y",
 'namegen'=["u5Y","u2h"],
 'gen'=[PG_u5Y,PG_u2h],
 'GF'=GF_Y,
 'nameC'=["e","c5","c5^2","u2","u3"],
 'GFC'=GFC_Y,
 'NC'=[1,12,12,15,20],
 'nameR'=["A","F1","F2","G","H"],
 'genR'=genR_Y):

tbPG["Yh"]:=Record(
 'name'="Yh",
 'namegen'=["ui5Y","u2h"],
 'gen'=[-PG_u5Y,PG_u2h],
 'GF'=GF_Yh,
 'nameC'=["e","c5","c5^2","u2","u3","i","ci5","ci5^3","sv","ui3"],
 'GFC'=proc(gen::list)  local i,ls;
   i:=FiniteGroups[FGS]((gen[1].gen[2])^3);
   ls:=GFC_Y([i.gen[1],gen[2]]);
   [op(ls),seq(i.M,M=ls)] end,
 'NC'=[1,12,12,15,20,1,12,12,15,20],
 'nameR'=["Ag","F1g","F2g","Gg","Hg","Au","F1u","F2u","Gu","Hu"],
 'genR'=[op(genR_Y),seq([-e[1],e[2]],e=genR_Y)]):


ModuleLoad()

end module: