## FileName: SimpleLieAlgebra.mpl
## Contents: Maple program for the simple Lie algebras
## Creation: Hideo Kodama 2012/10/25
## 2020/4/23: plete version 1.00
## 2929/9/16; Open access version 2.00
## 2020/10/11: modification of proc:SimpleRootsbyH
## 2020/10/20: proc: mkembM is included.
## 2020/10/21: version 2.10
## 2020/12/14: proc:mkembM is revised: allowing the Dynkin type D3 as a subalgebra and type C2 as the mother algebra => version 2.11
## 2021/2/7: bug fix  in proc:SimpleRootsbyH.  The typo was corrected as 'D' => 'E' 
## 2021/2/8: =>  version 2.12
## Last update: 2021/2/8

######>Initialization
with(LinearAlgebra):
with(StringTools):
with(plots):

#####>Global variables
KM:='KM': # Killing metric (possitive)
CM:='CM': # Cartan Matrix
alpha:='alpha': # simple root array
####>Initialise the plotsetup
plotsetup(default,plotoutput=terminal):

######>Job controle
Job_all:=true:

Job_SU2:=true:
Job_SU3:=false:
 Job_SU3_SubG:=true:
   Job_SU3_SubG_SO3:=true:
   Job_SU3_SubG_SU2U1:=true:
Job_SU4:=false:
  Job_SU4_WeightPlot:=true:
  Job_SU4_SubG_SU3:=true:
  Job_SU4_SubG_SU2SU2U1:=true:
  Job_SU4_SubG_SU2SU2:=true:
  Job_SU4_SubG_Sp2:=true:
Job_SU5:=false:
  Job_SU5_SubG:=true:
    Job_SU5_SubG_SU4U1:=true:
    Job_SU5_SubG_SU3SU2U1:=true:
    Job_SU5_SubG_SO5:=true:
Job_SU6:=false:
  Job_SU6_SubG_SU5U1:=true:
  Job_SU6_SubG_SU4SU2U1:=true:
  Job_SU6_SubG_SU3SU3U1:=true:
  Job_SU6_SubG_SU4:=true:
  Job_SU6_SubG_SU3:=true:
  Job_SU6_SubG_SU3SU2:=true:
  Job_SU6_SubG_Sp3:=true:
Job_SO5:=false:
 Job_SO5_SubG_SU2U1:=true:
 Job_SO5_SubG_SU2SU2:=true:
 Job_SO5_SubG_SU2max:=true:
Job_SO7:=false:
  Job_SO7_SubG_SU4:=true:
  Job_SO7_SubG_SO5U1:=true:
  Job_SO7_SubG_SU2SU2SU2:=true:
  Job_SO7_SubG_G2:=true:
Job_Sp3:=false:
  Job_Sp3_SubG_SU3U1:=true:
  Job_Sp3_SubG_Sp2SU2:=true:
  Job_Sp3_SubG_SU2:=true:
  Job_Sp3_SubG_SU2SU2:=true:
Job_SO8:=false:
  Job_SO8_SubG_SU4U1:=true:
  Job_SO8_SubG_SU2SU2SU2SU2:=true:
  Job_SO8_SubG_SO7:=true:
  Job_SO8_SubG_SU3:=true:
  Job_SO8_SubG_Sp2SU2:=true:
Job_SO9:=false:
  Job_SO9_SubG_SO8:=true:
  Job_SO9_SubG_SO7U1:=true:
  Job_SO9_SubG_SU4SU2:=true:
  Job_SO9_SubG_Sp2SU2SU2:=true:
  Job_SO9_SubG_SU2:=true:
  Job_SO9_SubG_SU2SU2:=true:
Job_Sp4:=false:
  Job_Sp4_SubG_SU4U1:=true:
  Job_Sp4_SubG_Sp3SU2:=true:
  Job_Sp4_SubG_Sp2Sp2:=true:
  Job_Sp4_SubG_SU2:=true:
  Job_Sp4_SubG_SU2SU2SU2:=true:
Job_Sp5:=false:
   Job_Sp5_SubG_SU5U1:=true:
Job_SO10:=false:
   Job_SO10_SubG_SO8U1:=true:
   Job_SO10_SubG_SU5U1:=true:
   Job_SO10_SubG_SU4SU2SU2:=true:
   Job_SO10_SubG_SO9:=true:
   Job_SO10_SubG_SO7SU2:=true:
   Job_SO10_SubG_Sp2Sp2:=true:
   Job_SO10_SubG_Sp2:=true:
Job_E6:=false:
   Job_E6_SubG_SO10U1:=true:
   Job_E6_SubG_SU6SU2:=true:
   Job_E6_SubG_SU3SU3SU3:=true:
   Job_E6_SubG_SU3:=true:
   Job_E6_SubG_G2:=true:
   Job_E6_SubG_G2SU3:=true:
   Job_E6_SubG_Sp4:=true:
   Job_E6_SubG_F4:=true:
Job_F4:=false:
   Job_F4_SubG_SO9:=true:
   Job_F4_SubG_SU3SU3:=true:
   Job_F4_SubG_Sp3SU2:=true:
   Job_F4_SubG_G2SU2:=true:
Job_G2:=false:
   Job_G2_SubG_SU3:=true:

if Job_all then
   Job_SU3:=true:
  Job_SU4:=true:
  Job_SU5:=true:
  Job_SU6:=true:
  Job_SO5:=true:
  Job_SO7:=true:
  Job_SO8:=true:
  Job_SO9:=true:
  Job_SO10:=true;
  Job_Sp3:=true:
  Job_Sp4:=true:
  Job_Sp5:=true:
  Job_G2:=true:
  Job_F4:=true:
  Job_E6:=true:
  Job_F4:=true:
end if:

######>Basic procedures
#####>Dynkin type procs
####>proc:resolvedt0
#% resolvedt0(dt): dt=Xr => ["X", r]: X=A/B/C/D/E/F/G, r=0,1,...
resolvedt0:=proc()

local dt, X,t, r, n,len,i:

dt:=_passed[1]:
if is(dt, list) then
   if nops(dt)=2 and is(dt[1],string) and is(dt[2],integer) and dt[2]>=0 
      and member(dt[1],["A","B","C","D","E","F","G"]) then
      dt;
   else
     error("dt=Xr or dt::list=[X::string, r]: X=A/B/C/D/E/F/G, r=0,1,...");
   end if:
else
  dt:=convert(dt,string):
  X:=substring(dt,1..1):
  len:=length(dt):
  r:=0: i:='i': for i from 2 to len do 
     r:=10*r + Ord(dt[i])-Ord("0"):
  end do:
  if member(X,["A","B","C","D","E","F","G"]) and r>=0 then
    [X,r];
  else
     error("dt=Xr or dt::list=['X', r]: X=A/B/C/D/E/F/G, r=0,1,...");
  end if:
end if:
end proc:

####>proc:resolvedt
#% resolvedt(dt0): dt0=Xr[i] => ["X",r,i]
resolvedt:=proc()
local dt0,dt;

dt0:=_passed[1]:
if is(dt0,list) then
   if nops(dt0)=2 then 
      return(resolvedt0(dt0));
   elif nops(dt0)>=3 then
      return [op(resolvedt0(dt0[1..2])),dt0[3]];
   else
      error(dt0,"Invalid Dynkin type notation");
   end if:
elif is(op(0,dt0)=`symbol`) then
   return resolvedt0(dt0);
end if:
dt:=op(0,dt0):
return [op(resolvedt0(dt)),op(dt0)];
end proc:

####>proc:abbrevdt0
#% abbrevdt0(dt): dt=["X",r]=>Xr
abbrevdt0:=proc()
local dt,dt1,X,t,r,n,len,i:

dt:=_passed[1]:
if is(dt,symbol) then
   return ds;
elif is(dt,list) and nops(dt)=2 then 
  if member(dt[1],["A","B","C","D","E","F","G"]) and is(dt[2],integer) and dt[2]>0 then
    return convert(cat(dt[1],dt[2]),symbol);
  end if:
end if:
error("Invalid Dynkin type in the list representation");
end proc:

####>proc:abbrevdt
#% abbrevdt(dt0): dt0=["X",r,i] => Xr[i]
abbrevdt:=proc()
local dt0, dt;

dt0:=_passed[1]:
if not type(dt0,list) then
  return dt0;
end if:
if nops(dt0)<2 then 
  error("usage: abbrevdt0([\"X\",r,i])");
end if:
dt:=abbrevdt0([dt0[1],dt0[2]]):
if nops(dt0)=2 then
  return dt;
else
  return dt[dt0[3]];
end if:  
end proc:

#####>proc:SimpleRootsbyH
#% SimpleRootsbyH(dt)
#% dt=[t,n] => Matrix A[i,j]: A_{ij}=alpha_j(h_i):
#%  h_j: canonical basis of the Cartan subalgebra
#% for the algebra embedding  f: H' -> H, f(h_i')=B_i^j h_j,
#% f^{-1}(alpha_j)=alpha'_i M^i_j;  M=(A')^(-1) B A
#% dt[1]="A/B/C/D": 

SimpleRootsbyH:=proc()

  local dt, t,rn,i, AM:

  dt:=resolvedt(_passed[1]):
  t:=dt[1]: rn:=dt[2]:
  if not member(t,{"A","B","C","D","E","F","G"}) then
    error("Dynkin type should be in [A,B,C,D,E,F,G]");
  end if:
  
  if member(t,{"G","F","E"}) then
    AM:=IdentityMatrix(rn):
  else
     AM:=Matrix(rn):
    i:='i': for i from 1 to rn do AM[i,i]:=1: end do:
    i:='i': for i from 1 to rn-1 do AM[i,i+1]:=-1: end do:
  
    if t="A" then
      i:='i': for i from 1 to rn do AM[rn,i]:=AM[rn,i]+1: end do:
    elif t="C" then
      AM[rn,rn]:=2:
    elif t="D" then
      AM[rn,rn-1]:=1:
    end if:
  end if:
  return Transpose(AM);
end proc:

#####>proc:Kmetric
#% Kmetric(dt) => KM::Matrix  
# dt=[t,n]: t=type("A","B","C","D","E","F","G"), n=rank
#% -Killing metric in the root space KM[i,j]=gamma_{ij} wrt the simple root basis
#% A constant rescaling of KM[i,j] does not affect the final results.

Kmetric:=proc()

local dt,t,n,i, KM:

dt:=resolvedt(_passed[1]):
t:=dt[1]: n:=dt[2]:
KM:=Matrix(n,shape=symmetric):
if t="A" or (t="B" and n=1) or (t="C" and n=1) #or (t="D" and n<=3) 
  then
  for i from 1 to n-1 do
    KM[i,i]:=2:
    KM[i,i+1]:=-1:
  end do:
  KM[n,n]:=2:
elif t="B" then
  for i from 1 to n-1 do
    KM[i,i]:=2:
    KM[i,i+1]:=-1:
  end do:
  KM[n,n]:=1:
elif t="C" then
  for i from 1 to n-2 do
    KM[i,i]:=1:
    KM[i,i+1]:=-1/2:
  end do:
  KM[n-1,n-1]:=1:
  KM[n-1,n]:=-1:
  KM[n,n]:=2:
elif t="D" and n>=3 then
  for i from 1 to n-2 do
    KM[i,i]:=2:
    KM[i,i+1]:=-1:
  end do:
  KM[n-1,n-1]:=2:
  KM[n,n]:=2:
  KM[n-2,n]:=-1:
elif t="G" and n=2 then
  KM[1,1]:=3:
  KM[1,2]:=-3/2:
  KM[2,2]:=1:
elif t="F" and n=4 then
 KM[1,1]:=2:
 KM[1,2]:=-1:
 KM[2,2]:=2:
 KM[2,3]:=-1:
 KM[3,3]:=1:
 KM[3,4]:=-1/2:
 KM[4,4]:=1:
elif t="E" and n>=6 and n<=8 then
 for i from 1 to n-2 do
   KM[i,i]:=2:
   KM[i,i+1]:=-1:
 end do:
 KM[n-1,n-1]:=2:
 KM[n,n]:=2:
 KM[3,n]:=-1:
else
  error("No such simple group");
end if:
KM;
end proc:

#####>proc:Cmatrix
#% Cmatrix(dt) => C: Cartan matrix  
# dt=[t,n]]: t=type(A,B,C,D,E,F,G), n=rank
#% Cartan matrix  C[i,j] 

Cmatrix:=proc()

local dt,t,n,i,j, KM,CM:

dt:=resolvedt(_passed[1]):
t:=dt[1]: n:=dt[2]:
KM:=Kmetric(dt):
CM:=Matrix(n):
for i from 1 to n do
  j:='j':
  for j from 1 to n do
    CM[i,j]:=2*KM[i,j]/KM[j,j]:
  end do:
end do:
CM;
end proc:

#####>proc:Gmetric
#% Gmetric(dt) => G_{ij}: metric in the Dynkin basis
Gmetric:=proc()
local dt,t,n,i, KM:

dt:=resolvedt(_passed[1]):
return Cmatrix(dt)^(-1).Kmetric(dt).Transpose(Cmatrix(dt)^(-1));
end proc:

#####>proc:IPDB
#% IPDB( dl1, dl2, dt) => the inner product (dl1,dl2) wrt Dynkin basis
IPDB:=proc()
local dl1,dl2,dt;

dl1:=_passed[1]: dl2:=_passed[2]:
dt:=resolvedt(_passed[3]):
return DotProduct(Vector(dl1),Gmetric(dt).Vector(dl2));
end proc:

#####>proc:RIP
#% RIP(r1,r2,dt) # r1,r2::list = root vectors, dt=Dynkin type =[t,n]
#% Killing inner product for roots
RIP:=proc()

local r1,r2,dt,KM,v1,v2:

r1:=_passed[1]: r2:=_passed[2]:
dt:=resolvedt(_passed[3]):
KM:=Kmetric(dt):
v1:=Vector(r1): v2:=Vector(r2):
Transpose(v2).KM.v1;

end proc:

#####>proc:RCP
#% RCP(r1,r2,dt) # r1,r2::list = root vectors, dt=Dynkin type =[t,n]
#% Cartan product for roots: 2 (r1,r2)/(r2,r2)

RCP:=proc()

local r1,r2,dt:

r1:=_passed[1]: r2:=_passed[2]:
dt:=resolvedt(_passed[3]):

2*RIP(r1,r2,dt)/RIP(r2,r2,dt);

end proc:

#####>proc:HighestRoot
#% HighestRoot(dt) highest root::list (in the simple root basis)
# dt=(t,n): t=Lie algebra type, n=rank
HighestRoot:=proc()

local dt, t,n,i, hr,CM,CMI:

dt:=resolvedt(_passed[1]):
t:=dt[1]: n:=dt[2]:
CM:=Cmatrix(dt): 
CMI:=CM^(-1):

if t="A" or (t="B" and n=1) or (t="C" and n=1) #or (t="D" and n<=3) 
then
  hr:=[seq(CMI[1,i]+CMI[n,i],i=1..n)]:
elif t="B" then
  hr:=[seq(CMI[2,i],i=1..n)]:
  if n=2 then hr:=2*hr: end if:
elif t="C" then
  hr:=[seq(2*CMI[1,i],i=1..n)]:
elif t="D" and n>=3 then
  if n=3 then
    hr:=[seq(CMI[2,i]+CMI[3,i],i=1..n)]:
  else
    hr:=[seq(CMI[2,i],i=1..n)]:
  end if:
elif t="G" and n=2 then
  hr:=[seq(CMI[1,i],i=1..n)]:
elif t="F" and n=4 then
  hr:=[seq(CMI[1,i],i=1..n)]:
elif t="E" and n>=6 and n<=8 then
  if n=6 then
    hr:=[seq(CMI[6,i],i=1..n)]:
  elif n=7 then
    hr:=[seq(CMI[1,i],i=1..n)]:
  else
    hr:=[seq(CMI[7,i],i=1..n)]:
  end if: 
else
  error("No such simple group");
end if:
hr;
end proc:

#####>proc:Rlevel
#% Rlevel( r) # r::list=root vector in the SRB => level height
Rlevel:=proc()

local r,nr, i:

r:=_passed[1]:
nr:=nops(r):

add(r[i],i=1..nr);

end proc:

#####>proc:RootSystem
####>proc:RootSystem0
##% RootSystem0(dt[,o]) => Root system for the Dynkin type dt
#% dt=[t,nr]: t=Dynkin type, nr=rank, o=output style (1/0)
RootSystem0:=proc() 

local dt,t,nr, prj, i,j, hrv, hl, RootSys,beta, l,l1, rtv,rtv1,rtv2,p,q,pmq,dim:

dt:=resolvedt(_passed[1]):
t:=dt[1]: nr:=dt[2]:
if _npassed=2 then
  prj:=_passed[2]:
else
  prj:=0:
end if:

hrv:=HighestRoot([t,nr]):
i:='i':
hl:=add(hrv[i],i=1..nr):
RootSys[hl]:=[hrv]:
if prj>0 then
  print(hl=RootSys[hl]);
end if:
dim:=1:
l:='l':
for l from hl-1 to 1 by -1 do
  RootSys[l]:=[]:
  for rtv in RootSys[l+1] do
    i:='i':
    for i from 1 to nr do
       j:='j':
       beta:=[seq(0,j=1..nr)]: beta[i]:=1:
       pmq:=RCP(rtv,beta,[t,nr]):
       q:=0:
       rtv1:=rtv:
       rtv1[i]:=rtv1[i]+1:
       rtv2:=rtv:
       rtv2[i]:=rtv2[i]-1: 
       if not member(rtv2,RootSys[l]) then
         l1:='l1':
         for l1 from l+1 to hl do
           if member(rtv1,RootSys[l1]) then
             rtv1[i]:=rtv1[i]+1:
             q:=q+1:
           end if:  
         end do:
         p:=q+pmq:
         if p>0 then
           RootSys[l]:=[op(RootSys[l]),rtv2]:
         end if:
       end if:
     end do:
  end do:
  RootSys[l]:=sort(RootSys[l]):
  dim:=dim+nops(RootSys[l]):
  if prj>0 then
    print(l=RootSys[l]);
  end if:
end do:
if prj>0 then
  print('dimension'=2*dim+nr);
end if:
RootSys["type"]:=[t,nr]:
RootSys["dim"]:=2*dim+nr:
RootSys["hl"]:=hl:
RootSys:
end proc:

####>proc:SRBtoDB
#% SRBtoDB( srwv, dt)  
# srwv= weight vector in the simple root basis, 
# dt: Dynkin type=Ar, ... or ["A",r],..
SRBtoDB:=proc()

local wv,dt, i,j,CM,nr:

wv:=_passed[1]: dt:=resolvedt(_passed[2]):
nr:=dt[2]:
CM:=Cmatrix(dt):
[seq(add(wv[j]*CM[j,i],j=1..nr),i=1..nr)];
end proc:

####>proc:DBtoSRB
#% DBtoSRB( dwv, dt)  
# dwv= weight vector in the Dynkin basis,
# dt: Dynkin type=Ar, ... or ["A",r],..
DBtoSRB:=proc()

local wv,dt, i,j,CMI,nr:

wv:=_passed[1]: dt:=resolvedt(_passed[2]):
nr:=dt[2]:
CMI:=Cmatrix(dt)^(-1):
[seq(add(wv[j]*CMI[j,i],j=1..nr),i=1..nr)];
end proc:

####>proc:RootSystem
#% RootSystem(dt[,prn,basis]) => Root system in the Dynkin basis
#% dt=Dynkin type, prn=output control flag, baiss="S" or "D"
RootSystem:=proc() 

local dt,prn, BT,
      t,n,RS,RDS,hl, output,
      i,u;

dt:=resolvedt(_passed[1]):
t:=dt[1]: n:=dt[2]:
BT:="S":
if _npassed>=2 then
  BT:=_passed[2]:
end if:
prn:=0:
if _npassed>=3 then
  prn:=_passed[3]:
end if:


if BT="D" then 
  RS:=RootSystem0(dt,0):
  hl:=RS["hl"]:
  i:='i':
  for i from hl to 1 by -1 do
    RDS[i]:=map(u->SRBtoDB(u,dt),RS[i]):
    if prn>0 then
      print(i=RDS[i]);
    end if:
  end do:
  RDS["hl"]:=hl:
  RDS["dim"]:=RS["dim"]:
  RDS["type"]:=[t,n]:
  if prn>0 then
    print('dimension'=RDS["dim"]);
  end if:
  return RDS;
else
  return RootSystem0(dt,prn);
end if:
end proc:

#####>Str Constin for the Weyl Basis
####>proc:isposlist
#% isposlist(x) => true if x>0, false othewise
#% x: list of number 
#% x>0 <=> x_1>=0 ... x_n>=0  & x<>[0,...,0]
isposlist:=proc()
local x,i,zerolist:


if _npassed>=1 then
  x:=simplify(_passed[1]): 
  if not is(x,list) then 
    return false:
  end if:
  i:='i': zerolist:=[seq(0,i=1..nops(x))]:
  if x=zerolist then
    return false:
  else
    i:='i': for i from 1 to nops(x) do
     if is(x[i],numeric) and x[i]<0 then 
        return false:
     end if:
    end do:
  return true:
  end if
else
  error("usage: listgr(x::list, y::list)");
end if:
end proc:

####>proc:BianchiID
#% BianchiID(x,y,z,NC) => [[Ex,Ey],Ez]+[[Ey,Ez],Ex]+[[Ez,Ex],Ey] /E_x+y+z
#%
BianchiID:=proc()
local x, y, z, NN,NNtable,
      RootSet,w,zeroroot:
      
      x:=simplify(_passed[1]): 
      y:=simplify(_passed[2]):
      z:=simplify(_passed[3]):
      NNtable:=_passed[4]:
      NN:=NNtable["StrConst"]:
      RootSet:=NNable["RootSet"]:
      zeroroot:=simplify(0*RootSet[1]):

     if member(simplify(x+y+z),RootSet) then
         w:=0:
         print(x+y, y+z, z+x);
         if member(simplify(x+y),RootSet) then
           w:=w + NN[x,y]*NN[x+y,z]:
         elif simplify(x+y)=zeroroot then
           w:=w+DotProduct(Vector(NN[x,y]),Vector(z)):
         end if:
         if member(simplify(y+z), RootSet) then
           w:=w + NN[y,z]*NN[y+z,x]:
         elif simplify(y+z)=zeroroot then
           w:=w+DotProduct(Vector(NN[y,z]),Vector(x)):
         end if:
         if member(simplify(z+x),RootSet) then
           w:=w + NN[z,x]*NN[z+x,y]:
         elif simplify(z+x)=zeroroot then
           w:=w+DotProduct(Vector(NN[z,x]),Vector(y)):
         end if:
         if op(0,w)<>`+` then
            print("single term:", [x,y,z]);
         end if:
         return simplify(w);
      else
        return 0;
      end if: 
end proc:

####>proc:StrConstWB0
#% StrConstWB0([alpha,beta],dt[,RootSet]) => N2=NN(alpha,beta)^2
#% [E_alpha,E_beta]=NN(alpha,beta)E_(alpha+beta)
#%
StrConstWB0:=proc()
local dt, BR, RootSet,
      RootSys, hl, el, x, y, p,q, v:

#% Set parameters    
if _npassed>=2 then
  BR:=_passed[1]:
  if not type(BR,list) or nops(BR)<>2 then 
    error("usage: StrConstWB0([E_alpha,E_beta],DynkinType, RootSet)");
  end if:
  x:=BR[1]: y:=BR[2]:
  dt:=resolvedt(_passed[2]):
else
    error("usage: StrConstWB0([E_alpha,E_beta],DynkinType, RootSet)");
end if:

if _npassed>=3 then
  RootSet:=_passed[3]:
else 
  RootSys:=RootSystem(dt,"S",0):
  hl:=RootSys["hl"]:
  RootSet:=[]:
  el:='el': for el from 1 to hl do
    RootSet:=[op(RootSet),op(RootSys[el]),op(map(`*`,RootSys[el],-1))]:
  end do:
end if:


#%  Calculate NN
  q:='q': for q from 0 by 1 while member(simplify(y+(q+1)*x),RootSet) 
          do end do:
  p:='p': for p from 0 by 1 while member(simplify(y-(p+1)*x),RootSet)
           do end do:
  v:=Vector(x):
  simplify(1/2*q*(p+1)*DotProduct(v,Kmetric(dt).v)); #=N2=StrConst^2
end proc:
 
 
####>proc:StrConstWB
#% StrConstWB( dt[,sw]) => Weyl Basic CCR coeff.  NN=table(N[wv1,wv2])
#% 
StrConstWB:=proc()
local dt, sw,
      sgn,epsilon, RootSys,hl, RootSet, NN,sgnidx,
      eqs,sgnsol,sgnsol1,zeroroot,NNtable,
      i,j,k,el,x,y,z,w, p,q,xv,N2,n;

if _npassed>=1 then
  dt:=resolvedt(_passed[1]):
  sw:=0:
  if _npassed>=2 then
     sw:=_passed[2]:   #% message control switch
  end if:
else
  error("usage: StrConstWB( dt[, sw=0/1])");
end if:
 

###% construct root system of the Dynkin type dt

RootSys:=RootSystem(dt,"S",0):
hl:=RootSys["hl"]:
RootSet:=[]:
el:='el': for el from 1 to hl do
  RootSet:=[op(RootSet),op(RootSys[el]),op(map(`*`,RootSys[el],-1))]:
end do:

#print("RootSet"=RootSet);

###% construct structure constant table
sgn:=epsilon:
sgnidx:=0:
el:='el': for el from 1 to hl do
  x:='x': for x in RootSys[el] do
     y:='y': for y in RootSet do
        if y=simplify(-x) then
           xv:=Vector(x):
           NN[x,-x]:=convert( Kmetric(dt).xv, list):
           NN[-x,x]:=-1*NN[x,-x]:
       else  
        z:=simplify(x+y): 
        if not assigned(NN[x,y]) then 
           q:='q':
           for q from 0 by 1 while member(simplify(y+(q+1)*x),RootSet) 
           do end do:
           p:='p': 
           for p from 0 by 1 while member(simplify(y-(p+1)*x),RootSet)
           do end do:
           xv:=Vector(x):
           N2:=1/2*q*(p+1)*DotProduct(xv,Kmetric(dt).xv):
           if N2>0 then
              sgnidx:=sgnidx+1:
              w:=sqrt(N2)*sgn[sgnidx]:
           else
              w:=0:
           end if:
           NN[x,y]:=w:
           NN[y,x]:=-w:
           NN[-x,-y]:=-w:
           NN[-y,-x]:=w:
           if sw>0 and w<>0 then
              print("NN",[x,y]=NN[x,y]);
           end if:
        end if:
      end if:
     end do:
  end do:
end do:
x:='x': for x in RootSet do
  y:='y': for y in RootSet do
    if not assigned(NN[x,y]) then 
       NN[x,y]:=0:
    end if:
  end do:
end do:

###% determine signs epsilon[*] of the str. consts
n:=nops(RootSet):
zeroroot:=simplify(0*RootSet[1]):
#sgn[1]:=1:
eqs:={}:
i:='i': for i from 1 to n-2 do
  j:='j': for j from i+1 to n-1 do
    k:='k': for k from j+1 to n do
      x:=RootSet[i]: y:=RootSet[j]: z:=RootSet[k]:
      if member(simplify(x+y+z),RootSet) then
          w:=0:
         if member(simplify(x+y),RootSet) then
           w:=w + NN[x,y]*NN[x+y,z]:
         elif simplify(x+y)=zeroroot then
           w:=w+DotProduct(Vector(NN[x,y]),Vector(z)):
         end if:
         if member(simplify(y+z), RootSet) then
           w:=w + NN[y,z]*NN[y+z,x]:
         elif simplify(y+z)=zeroroot then
           w:=w+DotProduct(Vector(NN[y,z]),Vector(x)):
         end if:
         if member(simplify(z+x),RootSet) then
           w:=w + NN[z,x]*NN[z+x,y]:
         elif simplify(z+x)=zeroroot then
           w:=w+DotProduct(Vector(NN[z,x]),Vector(y)):
         end if:
         w:=simplify(w):
          if w<>0 then
           eqs:={op(eqs), w}:
#           
#           if op(0,w)<>`+` then
#             print("single term: [x,y,z]"=[x,y,z]);
#           end if:
         end if:
      end if:
    end do:
  end do:
end do:
# print("eqs"=eqs);

i:='i': sgnsol:=solve(eqs,[seq(sgn[i],i=1..sgnidx)]):
#print(sgnsol);
i:='i':
if is(sgnsol,list) and nops(sgnsol)>1 then
  sgnsol1:=sgnsol[1]:
else
  sgnsol1:=sgnsol:
end if:
#print("sgnsol1"=sgnsol1);
NNtable["StrConst"]:=eval(NN, sgnsol1):
NNtable["RootSystem"]:=RootSys:
NNtable["RootSet"]:=RootSet:

return NNtable;
end proc:

####>proc:printStrConst
#% printStrConst(NNtable) => show the list of non-vanishing SC
printStrConst:=proc()
local NNtable,
      NN, RootSet, zeroroot,E, H,CRlist,
      n,i,j,k,x,y,z,c:

NNtable:=_passed[1]:
NN:=NNtable["StrConst"]:
RootSet:=sort(NNtable["RootSet"],(x,y)->op(1,x)+op(2,x)>op(1,y)+op(2,y)):
zeroroot:=simplify(0*RootSet[1]):
n:=nops(RootSet):
CRlist:=[]:
i:='i': for i from 1 to n-1 do
  j:='j': for j from i+1 to n do
     x:=RootSet[i]: y:=RootSet[j]:
     c:=NN[x,y]:
     if c<>0 then 
         z:=simplify(x+y):
        if z=zeroroot then
          k:='k':
          CRlist:=[ op(CRlist), [E[op(x)],E[op(y)]]=add(c[k]*H[k],k=1..nops(zeroroot))];
        else
          CRlist:=[ op(CRlist), [E[op(x)],E[op(y)]]= NN[x,y]*E[op(z)]];
        end if;
     end if:
  end do:
end do:
CRlist;
end proc:


#####>proc:WeylTrf
#% WeylTrf(dt, rv[, wt]) => Weyl trf matrix wrt the root vector rv
#% dt=Dynkin type,  rv=root vector, wt="S" or "D"
WeylTrf:=proc()
local dt, rv, wt,
      rank, i,j,k,IDM,WT,alpha:

if _npassed<2 or not is(_passed[2],list) then
  error("Usage: WeylTrf(Dynkin type, root vector)");
end if:
dt:=resolvedt(_passed[1]):
rank:=dt[2]:
rv:=_passed[2]:
if nops(rv)<> rank then
  error("The lenght of the root vector should be the same as the rank");
end if:
wt:="S":
if _npassed>=3 then
  wt:=_passed[3]:
end if:
#if wt="D" then
#   rv:=DBtoSRB(rv):
#end if:

#% unit basis vectors alpha[i]
i:='i': for i from 1 to rank do
  j:='j': k:='k':
  alpha[i]:=[seq(0,j=1..(i-1)),1,seq(0,k=1..(rank-i))]:
end do:

i:='i':
IDM:=DiagonalMatrix([seq(1,i=1..rank)]):
i:='i':
WT:= IDM - Matrix(rank,1,rv).Matrix(1,rank,[seq(RCP(alpha[i],rv,dt),i=1..rank)]):
if wt="D" then 
  WT:=Transpose(Cmatrix(dt)).WT.Transpose(Cmatrix(dt)^(-1)):
end if:
return WT;
end proc:

#####>WeightSystem
####>proc:WeightSystem
#% WeightSystem(dt,hdw,basis, sw) => WeightSystem
#% hdw=hight weight in the Dynkin basis
#% dt=[t,nr], t:Dynkin type =A,B,C,D,E,F,G
#% basis="D"(Dynkin) or "S"(Simple Root)
#% sw= 0 (no message) 1 (with message) 2 (outpu only weight list)
WeightSystem:=proc()

local dt, hdw,t,sw,basis, 
      nr, CM,CMI, i,j, hwv, hl, WeightSys,RS, WL,
      wtv,wtv1,wtv2,p,q,pmq,dim,
      wtd,wtd1,DWS,WS,dw,k,ml,u,x,uwv,m,rv,
        l,l1, ghl, prt, beta, kmax:

dt:=resolvedt(_passed[1]):
t:=dt[1]: nr:=dt[2]:
hdw:=_passed[2]:
if nops(hdw)<>nr then
  error("The weight vector does not match the rank");
end if: 
basis:="D":
if _npassed>=3 then
  basis:=_passed[3]:
end if:
sw:=0:
if _npassed>=4 then
  sw:=_passed[4]:
end if:

CM:=Cmatrix(dt):
CMI:=CM^(-1):
i:='i': j:='j':
hwv:=[seq(add(hdw[i]*CMI[i,j],i=1..nr),j=1..nr)]:
#if dt[1]="B" and nr=2 and hdw=[0,1] then
#  hwv:=2*hwv:
#end if:

####% Construction of weight vector lists WeightSys & DWS 
i:='i':
hl:=add(hwv[i],i=1..nr):
WeightSys[hl]:=[hwv]: #% WeightSys =weight vectors in the simple root basis
DWS[hl]:=[hdw]:       #% DWS = weight vectors in the Dykin basis
for l from hl-1 to -hl by -1 do
  WeightSys[l]:=[]: DWS[l]:=[]:
  for wtv in WeightSys[l+1] do
    i:='i':
    for i from 1 to nr do
       j:='j':
       beta:=[seq(0,j=1..nr)]: beta[i]:=1: #% i-th simple root
       pmq:=RCP(wtv,beta,dt):
       q:=0:
       wtv1:=wtv:
       wtv1[i]:=wtv1[i]+1:
       wtv2:=wtv:
       wtv2[i]:=wtv2[i]-1:
       if not member(wtv2,WeightSys[l]) then
         l1:='l1':
         for l1 from l+2 to hl do  # changed!!
           if member(wtv1,WeightSys[l1]) then
              q:=q+1:
              wtv1[i]:=wtv1[i]+1:
           end if:  
         end do:
         p:=q+pmq:
         if p>0 then
           WeightSys[l]:=[op(WeightSys[l]),wtv2]:
           j:='j'; k:='k':
           dw:=[seq(add(wtv2[k]*CM[k,j],k=1..nr),j=1..nr)]:
           DWS[l]:=[op(DWS[l]),dw]:
         end if:
       end if:
     end do:
  end do:
#  print(l=DWS[l]);
end do:

####% simple weight list output
if sw=2 then
  WL:=[]: l:='l': for l from hl to -hl by -1 do
    if basis="D" then
       for x in DWS[l] do
         WL:=[op(WL),x]:
        end do:
    else
      for x in WeightSys[l] do
        WL:=[op(WL),x]:
      end do:
    end if:
  end do:
  return WL;
end if:

####% Calculation of weight multiplicity
RS:=RootSystem(dt):
#% uwv= 1/2 x the sum of all positive root vectors
ghl:= Rlevel(HighestRoot(dt)):
i:='i':
uwv:=[seq(0,i=1..nr)]:
l:='l': 
for l from ghl by -1 while l>0  do
 rv:='rv':
 for rv in RS[l] do
   uwv:=uwv+1/2*rv:
 end do:
end do:  
#% Highest level
DWS[hl][1]:=[1,DWS[hl][1]]:  # adding the multiplicity info
WS[hl][1]:=[1,WeightSys[hl][1]]:
WS[hl]:=convert(WS[hl],list):
DWS["hl"]:=hl:
WS["hl"]:=hl:
if sw=1 then
 if basis="D" then
    print(hl=DWS[hl]);
 else 
    print(hl=WS[hl]);
 end if:
end if:
#% Lower levels
dim:=1:
l:='l':
for l from hl-1 to -hl by -1 do
   i:='i': 
   for i from 1 to nops(WeightSys[l]) do
      wtv:=WeightSys[l][i]:
      x:=0:
      l1:='l1':
      for l1 from 1 to ghl do     #
         for prt in RS[l1] do     # for all positive root  prt
            k:='k':
            kmax:= trunc((hl-l)/l1):
            for k from 1 to kmax do  # for all k st wtv+k*prt in Delta^+
              m:='m':
              if member(eval(wtv+k*prt),WeightSys[l+k*l1],m) then
                 x:=x+2*DWS[l+k*l1][m][1]*RIP(wtv+k*prt,prt,dt):
              end if:
            end do:
         end do:
      end do:
      ml:=x/(RIP(hwv+uwv,hwv+uwv,dt)-RIP(wtv+uwv,wtv+uwv,dt)):
      DWS[l][i]:=[ml,DWS[l][i]]:
      WS[l][i]:=[ml,WeightSys[l][i]]:
      dim:=dim+ml:
   end do:
   WS[l]:=convert(WS[l],list):
  if sw=1 then
     if basis="D" then
       print(l=DWS[l]);
     else
       print(l=WS[l]);
     end if:  
  end if:
end do:     
DWS["dt"]:=abbrevdt(dt): WS["dt"]:=abbrevdt(dt):
DWS["hdw"]:=hdw: WS["hdw"]:=hdw:
DWS["dim"]:=dim:  WS["dim"]:=dim:
if sw>=1 then
  print('dim'=dim,  'hl'=hl);
end if:
if sw<=1 then
  if basis="D" then
    return DWS:
  else
    return WS:
  end if:
end if:
end proc:

####>proc:printWS
#% printWS(WS/DWS)
printWS:=proc()

local WS, 
      hl, i,j;
      
WS:=_passed[1]:
hl:=WS["hl"]:
i:='i': for i from hl to -hl by -1 do
  print(i=WS[i]);
end do:
print('dim'=WS["dim"], 'hl'=hl);
end proc:


####>proc:isDominantWeight
#% isDominantWeight(dw) => true/false
isDominantWeight:=proc()
local dw, ans, x:

dw:=_passed[1]:
ans:=true:
for x in dw while ans do
  if x<0 then
    ans:=false:
  end if:
end do:  
return ans:
end proc:

####>proc:mklattice
#% mklattice(rn::integer,[n1,n2]) => lattice point list
mklattice:=proc()
local rn,rng,
      sublattice,newlattice,x,i:

rn:=_passed[1]:
if _npassed<2 or not type(rn,integer) or rn<1 then
   error("usage: mklattice(rn::positive integer, [n1::integer,n2::integer])");
end if:
rng:=_passed[2]:
if not type(rng,list) or nops(rng)<2 or not type(rng[1],integer) or not type(rng[2],integer) then
   error("usage: mklattice(rn::positive integer, [n1::integer,n2::integer])");
end if:

if rn=1 then
  i:='i': newlattice:=[seq([i],i=rng[1]..rng[2])]:
else
  sublattice:=mklattice(rn-1,rng):
  newlattice:=[]:
  for x in sublattice do
    i:='i':newlattice:=[op(newlattice),seq([op(x),i],i=rng[1]..rng[2])]:
  end do:
end if:
return newlattice;
end proc:

####>proc:mklevellattice
#% mklevellattice(rn::integer>0, level::integer>=0)=> level lattice point set 
mklevellattice:=proc()
local rn, level,
      x,lattice,
      levellattice:
if _npassed<2 then
  error("usage: mklevellattice(rn::intger>0,level::integer>=0)");
end if:
rn:=_passed[1]:
level:=_passed[2]:
lattice:=mklattice(rn,[0,level]):
levellattice:=[]:
x:='x': for x in lattice do
  if Rlevel(x)=level then
    levellattice:=[op(levellattice),x]:
  end if:
end do:
return levellattice;
end proc:

####>proc:findHighestWeight0
#% findHighestWeight0(dw0,dt) => dominant weight list
findHighestWeight0:=proc()
local  dw0,dt,
       rn,RS,SimpleRoots,maxrlevel, found,rlevel,rootlist, 
       dw, dwlist,
      i,j,k, x, y;

dw0:=_passed[1]:
dt:=resolvedt(_passed[2]):
rn:=dt[2]:
RS:=RootSystem(dt,"D",0):
SimpleRoots:=RS[1]:

found:=false:
rlevel:='rlevel':
maxrlevel:=10:
dwlist:=[]:
for rlevel from 1 to maxrlevel while not found do
 rootlist:=mklevellattice(rn,rlevel):
 x:='x': for x in rootlist do
    i:='i': dw:=simplify(dw0+add(x[i]*SimpleRoots[i],i=1..rn)):
    if isDominantWeight(dw) then
       dwlist:=[op(dwlist),dw]:
       found:=true:
    end if:
  end do:
end do:

if found then
  return dwlist;
else
  print("Not found");
end if:
end proc:

####>proc:findHighestWeight
#% findHighestWeight(dw0,dt) => dominant weight list
findHighestWeight:=proc()
local  dw0,dt,
       rn,RS,w0,level0,SimpleRootDB,maxrlevel, found,rlevel,
       dw, w, q,qmax,bd,
      i,j,k, x, y;

dw0:=_passed[1]:
dt:=resolvedt(_passed[2]):
rn:=dt[2]:
w0:=DBtoSRB(dw0,dt):
level0:=`+`(op(w0)):
i:='i': j:='j':
SimpleRootDB:=[seq([seq(Cmatrix(dt)[i,j],i=1..rn)],j=1..rn)]:

maxrlevel:=50:
rlevel:=1:
found:=false:
dw:=dw0:
while not found and rlevel<maxrlevel do
  w:=DBtoSRB(dw,dt):
  qmax:=0: bd:=1:
  i:='i': for i from 1 to rn do
    j:='j': q:=-add(w[j]*Cmatrix(dt)[j,i],j=1..rn):
    if q> qmax then
      qmax:=q:  bd:=i:
    end if:
  end do:
  rlevel:=rlevel+qmax:
  dw:=dw+qmax*SimpleRootDB[bd]:
  if isDominantWeight(dw) then
     found:=true:
  end if:
end do:

return dw;

end proc:

#####>Rep. Product Reduction
####>proc:ProductOfRep
#% ProductOfRep(dt,hdw1,hdw2[, prn,basis])
#% dt=[t,rn] or e.g. A4
ProductOfRep:=proc() 

local hdw1,hdw2,dt,t,nr,prn,basis,
  hdw,DWS1,DWS2, CDWS,CSWS, hl1,hl2,chl,
  IRDWS,i,j,k,l,wd1,wd2,wd,TDWS,cw,cm,chl0,
  cdw,ml,tmp,ReprList,dim;

dt:=resolvedt(_passed[1]): 
t:=dt[1]: nr:=dt[2]:#% Dynkin type
hdw1:=_passed[2]: #% 1st highest Dynkin weight
hdw2:=_passed[3]: #% 2nd ...
if _npassed>=4 then
  prn:=_passed[4]:
  basis:="D":
  if _npassed>=5 then
    basis:=_passed[5];
  end if;
else
  prn:=0:
end if:

hdw:=hdw1+hdw2:
DWS1:=WeightSystem(dt,hdw1):
DWS2:=WeightSystem(dt, hdw2):
hl1:=DWS1["hl"]:
hl2:=DWS2["hl"]:
chl:=hl1+hl2:
if prn>=2 then
  if basis="S" then
     print("Product weight system (Simple root basis):",'dim'=DWS1["dim"]*DWS2["dim"]);
  else
     print("Product weight system (Dynkin basis):",'dim'=DWS1["dim"]*DWS2["dim"]);
  end if:
end if:

i='i':
for i from chl to -chl by -1 do
  CDWS[i]:=[]:  #% the list of the composite Dynkin weights with multipicity
  CSWS[i]:=[]:  #% CDWS[i] in the simple root basis
  TDWS[i]:=[]:  #% the simple list of the Dynkin weight sum
  j:='j':
  for j from hl1 to -hl1 by -1 do
     k:=i-j:
     if is(DWS1[j],list) and is(DWS2[k],list) then
        wd1:='wd1':
        for wd1 in DWS1[j] do
           wd2:='wd2':
           for wd2 in DWS2[k] do
               cw:=wd1[2]+wd2[2]: cm:=wd1[1]*wd2[1]:
               l:='l':
               if member(cw,TDWS[i],l) then
                  CDWS[i][l][1]:=CDWS[i][l][1]+cm:
                  CSWS[i][l][1]:=CSWS[i][l][1]+cm:
               else
                  TDWS[i]:=[op(TDWS[i]),cw]:
                  CDWS[i]:=[op(CDWS[i]),[cm,cw]]:
                  CSWS[i]:=[op(CSWS[i]),[cm,DBtoSRB(cw,dt)]]:
               end if:
           end do:
        end do:
     end if:
  end do:
  if prn=2 then
    if basis="S" then
      print(i=CSWS[i]);
    else
      print(i=CDWS[i]);
    end if:
  end if:
end do:

if false then
dim:=0:
l:='l':
for l from chl to -chl by -1 do
  wd:='wd':
  for wd in CDWS[l] do
    dim:=dim+wd[1]:
  end do:
end do:
print("Check ml sum"=dim);
end if:

if prn>=1 then
  print("Component representations:");
end if:
ReprList:=[]:
cdw:=hdw:
chl0:=chl:
while chl>=0 do
# print('chl'=chl);
# print(CDWS[chl]);
 cdw:='cdw':
 for cdw in CDWS[chl] do
 ml:=cdw[1]:
 if ml>0 then
  IRDWS:=WeightSystem(dt,cdw[2]):
  i:='i':
  for i from chl to -chl by -1 do
     wd:='wd':
     for wd in IRDWS[i] do 
       l:='l':
       if member(wd[2],TDWS[i],l) and CDWS[i][l][1]>0 then 
          CDWS[i][l][1]:=CDWS[i][l][1]-ml*wd[1]:
       end if:
     end do:
  end do:
  ReprList:=[op(ReprList),[cdw,'dim'=IRDWS["dim"],'ml'=ml]]:
  if prn>=1 then
    print(cdw[2], 'dim'=IRDWS["dim"],'ml'=ml);
#  print("CDWS",sort(CDWS));
  end if:
 end if:
 end do:
 chl:=chl-1:
end do:
#i:='i': 
#for i from chl0 to -chl0 by -1 do
#  print(CDWS[i]);
#end do;
if prn<=1 then
return ReprList:
end if:      
end proc:

#####> Subgroup reduction
####>proc: sbtDWSm 
#% sbtDWSm:=proc(dt, DWS, dw) => DWS2: reduced dw list
#% dt=[[t1,rn1],..,[tp,rnp]]
#% DWS=[ h=[[ml,dw],...], ...], dim, hl]
#% dw=[dwv1,..,dwvp]: highest Dynkin weight vector pair in DWS
#% output: [ml, dim, DWS2[Q]=DWS[Q] - DWS[Q][h][dwv]]
sbtDWSm:=proc()

local dt, DWS, dw, sgc,
   dim0,hl0, hls,hv, hvlist,dim1, dimlist,hl1, hl2, 
   DWS1,DWS2,dws2, dws3,ml1,DWSs,dwsv,hlssum, dwcount,
  x,y,i,j,h,dw1,ml,flag,found,tmp:

if _npassed<3 then
  error("usage: sbtDWSm(dt,DWS, dw)");
end if:

sgc:=nops(_passed[1]):
dt:=map(resolvedt,_passed[1]): 
DWS:=_passed[2]: dw:=_passed[3]; 

#print(DWS);

##% constructing DWS1 to subtract for dw
dim0:=DWS["dim"]: hl0:=DWS["hl"]:
i:='i': for i from 1 to sgc do
  DWSs[i]:=WeightSystem(dt[i],dw[i]):
  hls[i]:=DWSs[i]["hl"]:
end do:
x:='x': i:='i':
tmp:=[seq(x[i],i=1..sgc)]:
i:='i': for i from 1 to sgc do
  j:='j': y:='y':
  tmp:=seq(eval(tmp,x[i]=y),y in [seq(hls[i]-j,j=0..2*hls[i])]):
end do:
hvlist:=[tmp]:

#print(hvlist);

hv:='hv': for hv in hvlist do
   i:='i': h:=add(hv[i],i=1..sgc):
   if not assigned(DWS1[h]) then DWS1[h]:=[]: end if:
   x:='x';
   tmp:=[seq(x[i],i=1..sgc)]:
   i:='i': for i from 1 to sgc do
     j:='j': y:='y':
     tmp:=seq(eval(tmp,x[i]=y), y in DWSs[i][hv[i]]):
   end do:
   dwsv:=[tmp]:
   x:='x': for x in dwsv do
      i:='i': ml:=product(x[i][1],i=1..sgc):
      i:='i': dw1:=[seq(x[i][2],i=1..sgc)]:
            DWS1[h]:=[op(DWS1[h]),[ml,dw1]]:
   end do:
end do:

i:='i':
dimlist:=[seq(DWSs[i]["dim"],i=1..sgc)]:
i:='i':
dim1:=product(dimlist[i],i=1..sgc):
DWS1["dim"]:=dim1:

#print(DWS1);


##% making the top level list of the new DWS2
i:='i': hlssum:=add(hls[i],i=1..sgc):
i:='i':
if hl0<>hlssum or not is(map(nops, dw)=map2(op,2,dt)) then 
  error("DWS and dw should have the same height and rank");
end if:

if nops(DWS[hl0])=1 then
  ml1:=DWS[hl0][1][1]:
  hl2:=hl0-1:
else
  DWS2[hl0]:=[]: hl2:=hl0:
  x:='x': for x in DWS[hl0] do
    if x[2]<>dw then 
       DWS2[hl0]:=[op(DWS2[hl0]),x]:
    else
       ml1:=x[1]:
    end if:
  end do:
end if:


##% making the lower level lists of DWS2
h:='h':
for h from hl0 to -hl0 by -1 do
   if assigned(DWS[h]) then
       dwcount:=nops(DWS[h]):
       dws2:=Array(1..dwcount):
       i:='i': for i from 1 to dwcount do
         dws2[i]:=DWS[h][i]:
       end do:
        x:='x': for x in DWS1[h] do
        i:='i': found:=false:
        for i from 1 to dwcount while not found do
           if DWS[h][i][2]=x[2] then 
              found:=true:
              dws2[i][1]:= DWS[h][i][1]-ml1*x[1]:
           end if:
        end do:
        if not found then
          error("DWS does not contain the irr. rep dwv");
        end if:
    end do:
    x:='x': dws3:=[]:
    i:='i': for i from 1 to dwcount do
      x:=dws2[i]:
      if x[1]>0 then dws3:=[op(dws3),x]: end if:
    end do:
#    print(dws3);
    if nops(dws3)>0 then
       DWS2[h]:=dws3:
     elif hl2=h then 
        hl2:=h-1:
     end if:
  end if:
end do:

h:='h':
for h from hl0-1/2 to -hl0+1/2 by -1 do
   if assigned(DWS[h]) and nops(DWS[h])>0 then 
      DWS2[h]:=DWS[h]:
      if h> hl2 then
         hl2:=h:
      end if:
   end if:
end do:

DWS2["hl"]:=hl2: 
DWS2["dim"]:=dim0-ml1*dim1:

return [ml1, dimlist, DWS2];

end proc:




####>proc: SubGrdm
#% SubGrdm:=proc(dt0,dwv0,dts,embMDs[,sw,dws0]) => DWL: Dynkin weight list
#% Irr. decomposition of an irrep. of Lie algebra (L,H) 
#% wrt a subalgebra (L',H')=(L1+..+Lp+U1,H1+..+Hp+HU1)
#% dt0:=[t0,rn0]: Dynkin type of (L,H)
#% dwv=highest Dynkin weight of (L,H) of type t0
#% dts=[[t1,rn1],..,[tp,rnp](,U1_1,...,U1_m)] = Dynkin types of the subalgebras
#% (L1+..l+Lp,H1+...+Hp) 
#% specified by the embeding matrix 
#% embMDs=Matrix[embMD[1],...,embMD[p](,QvecD)]: H^* -> H1^*+...+Hp^*(+U1): merged to a single matrix
#% output: the list of highest Dynkin weights for 
#% (L1+...+Lp(+U1+...+U1),H1+...+Hp(+HU1+...+HU1))
#%   DWL::table=[DWL(Q), Q in Qlist],
#%   DWL(Q)::list=[[dwv1,...,dwvp], dim, ml],...,[...], 
#% "Representation"=[[t0,rn0],dwv], 
#%  "SubGroup"=[dts,embMD] ] 
SubGrdm:=proc()

local dt0, dwv0,dts,embMDs,sgc0, sgc,sw,dws0,
      print_DWS1, print_DWS01,
      t0,rn0,ts,rns,CM0,CM,DWS0, ml0, 
      dwvs,dwvec0,dwvec,wvecs,Q,embMD,  
      U1count, U1plist, QvecDs, sgname, Qlist,QHLlist,
       DWS1,DWS2, hl0,dim0,dim1,hls, dwlist, DWS01,
       DWL1,DWL,dws,dws1, dwv1,dw1,sbtx,tdim,
      p,q, dwv, ml, x,y, h,h1,hl1,hs,i, j,k, flag, U1sw,tmp:

###% set parameters & create DWS0 for (dt0, dwv0)
if _npassed<4 then
  error("usage: SubGsU1reduction(dt0,dwv0,dts,embMDs)");
elif type(embMDs,list) and nops(_passed[3])<>nops(_passed[4]) then
  error("Subgroup list and embMDs are not consistent.");
end if:

###% option parameters
print_DWS1:=false: print_DWS01:=false: sw:=0:
if _npassed>=5 then
   sw:=_passed[5]:
   if sw=2 then
      print_DWS1:=true:
   elif sw>=3 then
      print_DWS01:=true:
   end if:
   if sw>=4 and _npassed>=6 then
      dws0:=_passed[6]:
   else
      dws0:=NULL:
   end if:
end if:


###% U1 list
dt0:=resolvedt(_passed[1]): dwv0:=_passed[2]: 
t0:=dt0[1]: rn0:=dt0[2]:
embMDs:=_passed[4]:

sgc0:=nops(_passed[3]):
U1plist:=[]:
if member(U1,_passed[3]) then
  U1sw:=true:
  i:='i': for i from 1 to nops(_passed[3]) do
    if _passed[3][i]=U1 then
       U1plist:=[op(U1plist),i]:
    end if:
  end do:
  U1count:=nops(U1plist):
  dts:=map(resolvedt,remove(has,_passed[3],U1)):
else
  U1sw:=false:
  U1count:=0:
  dts:=map(resolvedt,_passed[3]):
end if:


####% embMD/QvecD list
q:=1:
j:=0: QvecDs:=[]: 
i:='i': for i from 1 to sgc0 do
  if member(i,U1plist)  then
    if type(embMDs,list) then
      QvecDs:=[op(QvecDs),convert(embMDs[i],list)]:
    else
      QvecDs:=[op(QvecDs),convert(embMDs[q,1..rn0],list)]:
      q:=q+1:
    end if:
  else
    j:=j+1:
    ts[j]:=dts[j][1]: rns[j]:=dts[j][2]:
    CM[j]:=Cmatrix(dts[j]):
    if type(embMDs,list) then
       embMD[j]:=embMDs[i]:
    else
       embMD[j]:=embMDs[q..(q+rns[j]-1),1..rn0]:
       q:=q+rns[j]:
    end if:
  end if:
end do:
sgc:=sgc0-U1count:
if U1count=0 then
  QvecDs:=[convert(Vector(rn0),list)]:
  U1count:=1:
end if:

#i:='i': print([seq(embMD[i],i=1..sgc),convert(QvecD,list)]);

#print("Input params"=[dt0, dwv0, dts, embMDs,sgc]):

###% construct the irr rep weight system DWS0 for [dt0, dwv0]
if nops(dwv0)<>rn0 then
  error("The weight vector does not match the rank:", "dw vector"=dwv0, "rank"=rn0);
end if:
DWS0:=WeightSystem(dt0,dwv0):
hl0:=DWS0["hl"]: dim0:=DWS0["dim"]:

CM0:=Cmatrix(dt0):
dwvec0:=Vector(dwv0):

if sw>0 then
  print("Original Rep."=[cat(t0,rn0),dwv0,'dim'=dim0]);
end if:
if sgc>0 then
  sgname:=cat(ts[1],rns[1]):
  i:='i': for i from 2 to sgc do
    sgname:=cat(sgname,cat("x",cat(ts[i],rns[i]))):
  end do:
else
  sgname:="":
end if:
if sw>0 then
i:='i':
if U1sw then
  if sgc>0 then
     sgname:=sprintf(cat(cat(sgname,"xU(1)"),"^%a"),U1count):
  else
     sgname:=sprintf(cat(cat(sgname,"U(1)"),"^%a"),U1count):
  end if:
  print('Subgroup'=[sgname, 'embedding'=[seq(embMD[i],i=1..sgc),QvecDs]]);
else
  print('Subgroup'=[sgname, 'embedding'=[seq(embMD[i],i=1..sgc)]]);
end if:
end if:

#print(DWS0);

###% Irr DWS0 for [dt0, dwv0] to DWS1 for 
#% dts=[[t1,rn1],,,[tp,rnp]](xU(1))
h:='h':
for h from hl0 to -hl0 by -1 do
   i:='i': 
  for i from 1 to nops(DWS0[h]) do 
    ml0:=DWS0[h][i][1]:
    dwv:=DWS0[h][i][2]:
    j:='j': Q:=[seq(DotProduct(Vector(QvecDs[j]),Vector(dwv)),j=1..U1count)]:
    j:='j': for j from 1 to sgc do
      dwvec[j]:=embMD[j].Vector(dwv):
      wvecs[j]:=Transpose(CM[j]^(-1)).dwvec[j]:
       x:='x': hs[j]:=add(x, x in wvecs[j]):
    end do:
    
    j:='j':
    dwvs:=[seq(convert(dwvec[j],list),j=1..sgc), Q];
       k:='k': h1:=add(hs[k],k=1..sgc):      
       if assigned(hl1[Q]) then
         if h1>hl1[Q] then hl1[Q]:=h1; end if;
         dim1[Q]:=dim1[Q]+ml0:
       else
         hl1[Q]:=h1:
         dim1[Q]:=ml0:
       end if:
     
    if assigned(DWS1[Q][h1]) then
       j:='j': flag:=true:
        for j from 1 to nops(eval(DWS1[Q][h1])) do
         if DWS1[Q][h1][j][2]=dwvs[1..sgc] then
            DWS1[Q][h1][j][1]:=DWS1[Q][h1][j][1]+ml0:
             DWS01[Q][h1][j][1]:=eval(DWS01[Q][h1][j][1])+ml0:
            y:=eval(DWS01[Q][h1][j][3]);
            DWS01[Q][h1][j][3]:=[op(y),dwv]:
            flag:=false:
         end if:
       end do:
       if flag then
          j:='j': y:=[ml0,[seq(dwvs[j],j=1..sgc)]]:
          DWS1[Q][h1]:=[op(DWS1[Q][h1]), y]:
          DWS01[Q][h1]:=[op(DWS01[Q][h1]), [op(y),[dwv]]]:
       end if:
    else
       j:='j': y:=[ml0,[seq(dwvs[j],j=1..sgc)]]:
       DWS1[Q][h1]:=[y]:
       DWS01[Q][h1]:=[[op(y),[dwv]]]:
   end if:
  end do:   
end do:

Qlist:={}:
x:='x': for x in op(op(DWS1)) do
  Qlist:={op(Qlist),op(1,x)}:
end do:
x:='x': y:='y':
Qlist:=sort(convert(Qlist,list),(x,y)->op(1,x)>op(1,y)):
# DWS1["Qlist"]:=Qlist:


#print(DWS1);
if print_DWS1 or print_DWS01 then

print("Qlist"=Qlist);
#print("DWS1"=DWS1);
print("Derived Dynkin weight system:");
Q:='Q': for Q in Qlist do
  print("Q"=Q);
  x:='x': y:='y':
  if print_DWS01 then
    tmp:=sort(op(eval(DWS01[Q])),(x,y)->op(1,x)>op(1,y)):
  else
    tmp:=sort(op(eval(DWS1[Q])),(x,y)->op(1,x)>op(1,y)):
  end if:
  x:='x': for x in tmp do
    if is(2*op(1,x),integer) then
       print("level"=op(1,x));
       y:='y': for y in sort(op(2,x)) do
         if sw<=3 or member(op(2,y),dws0) then
            print(y);
         end if:
       end do:
    end if:
  end do:
end do:
end if:

Q:='Q': for Q in Qlist do
  DWS1[Q]["hl"]:=hl1[Q]:
  DWS1[Q]["dim"]:=dim1[Q]:
end do:


###% Irr decomposition of DWS1 => DWL
if sw>0 then
  print(cat("Derived Irr. Reps. wrt ",sgname));
end if:
QHLlist:=[]:
tdim:=0:
Q:='Q':
for Q in Qlist do
  DWS2:=DWS1[Q]:
   h:='h':
  for h from hl1[Q] by -1/2 while h>=0 do
    while assigned(DWS2[h]) and nops(DWS2[h])>0 do
       dw1:=DWS2[h][1][2]: ml:=DWS2[h][1][1]:
       sbtx:=sbtDWSm(dts,DWS2,dw1):
       if sgc=0 and U1sw then
          sbtx[2]:=1:
       end if:
       if assigned(DWL1[Q][h]) then 
           DWL1[Q][h]:=[op(DWL1[Q][h]),[dw1, sbtx[2], ml]]:
       else
           DWL1[Q][h]:=[[dw1, sbtx[2], ml]]:
           QHLlist:=[op(QHLlist),[Q,h]]:
       end if:
       i:='i': tdim:=tdim +ml*product(sbtx[2][i],i=1..sgc):
       if sw>0 then
         print([Q,h]=[dw1,sbtx[2],ml]);
       end if:
       DWS2:=sbtx[3]:
     end do:
  end do:
end do:
if tdim<>dim0 then
  error("Unknown error!  The total dimension count is not consistent.");
end if:
DWL["Induced Reps"]:=eval(DWL1):
DWL["Original Rep."]:=[cat(t0,rn0),dwv0]:
DWL["SubGroup"]:=[sgname, embMDs]:
DWL["Qlist"]:=Qlist:
DWL["QHLlist"]:=QHLlist:
return DWL;

end proc:

######>Embedding/Projection Matrices
#####>Common part
####>variables
embM:='embM':
Qvec:='Qvec':
####>Cartan matrix
#% CM[X]=Transpose(Cmatrix(X)):
SL:='SL': A:='A': 
r:='r': for r from 1 to 10 do
  CM[cat(SL,r+1)]:=Transpose(Cmatrix(cat(A,r))):
end do:
SO:='SO': B:='B': 
r:='r': for r from 1 to 10 do
  CM[cat(SO,2*r+1)]:=Transpose(Cmatrix(cat(B,r))):
end do:
r:='r': for r from 3 to 10 do
  CM[cat(SO,2*r)]:=Transpose(Cmatrix(cat('D',r))):
end do:
Sp:='Sp': C:='C': 
r:='r': for r from 2 to 10 do
  CM[cat(Sp,r)]:=Transpose(Cmatrix(cat(C,r))):
end do:
CM[G2]:=Transpose(Cmatrix(G2)):
CM[F4]:=Transpose(Cmatrix(F4)):
CM[E6]:=Transpose(Cmatrix(E6)):
CM[E7]:=Transpose(Cmatrix(E7)):
CM[E8]:=Transpose(Cmatrix(E8)):
####>Simple roots by dual Cartan basis
#% [alpha_1,...,alpha_r]=[h^1,...,h^r] SR
SL:='SL': A:='A': 
r:='r': for r from 1 to 10 do
  SR[cat(SL,r+1)]:=SimpleRootsbyH(cat(A,r)):
end do:
SO:='SO': B:='B': 
r:='r': for r from 1 to 10 do
  SR[cat(SO,2*r+1)]:=SimpleRootsbyH(cat(B,r)):
end do:
r:='r': for r from 3 to 10 do
  SR[cat(SO,2*r)]:=SimpleRootsbyH(cat('D',r)):
end do:
Sp:='Sp': C:='C': 
r:='r': for r from 1 to 10 do
  SR[cat(Sp,r)]:=SimpleRootsbyH(cat(C,r)):
end do:
SR[G2]:=DiagonalMatrix([1,1]):
SR[E6]:=Kmetric(E6):

####>Maximal Sugalgebra List
###>proc:MaxSubGlist
#% MaxSubGlist(dt) => the list of maximal semisimple subalgebras
MaxSubGlist:=proc()
local dt:

if _npassed=0 then 
  error("Usage: MaxSubGlist(dt)");
end if:

dt:=resolvedt(_passed[1]):
if dt[2]=1 then 
     return [[U1]],   [];
elif dt[2]=2 then
   if dt[1]="A" then
      return [[A1,U1]],  [];
   elif dt[1]="B" or  dt[1]="C" then
      return [[A1[1],A1[2]], [A1,U1]],  [[A1]];
   elif dt[1]="G" then
      return [[A2],[A1[1],A1[2]]],  [[A1]];
   end if:
elif dt[2]=3 then
   if dt[1]="A" or dt[1]="D" then
      return [[A2,U1],[A1[1],A1[2],U1]],  [[A1[1],A1[2]],[C2]];
   elif dt[1]="B" then
      return [[A3],[A1[1],A1[2],A1[3]],[C2,U1]],  [[G2]];
   elif dt[1]="C" then
      return [[A2,U1],[C2,A1]],  [[A1],[A1[1],A1[2]]];
   end if:
elif dt[2]=4 then
   if dt[1]="A" then
      return [[A3,U1],[A2,A1,U1]], [[C2]];
   elif dt[1]="B" then
      return [[D4],[C2,A1[1],A1[2]],[A3,A1],[B3,U1]],  [[A1],[A1[1],A1[2]]];
   elif dt[1]="C" then
      return [[A3,U1],[C3,A1],[C2[1],C2[2]]], [[A1],[A1[1],A1[2],A1[3]]];
   elif dt[1]="D" then
      return [[A1[1],A1[2],A1[3],A1[4]],[A3,U1]],  [[A2],[B3],[C2,A1]];
   elif dt[1]="F" then
      return[[B4],[A2[1],A2[2]],[C3,A1]],  [[G2,A1],[A1]];
   end if:
elif dt[2]=5 then
   if dt[1]="A" then
      return [[A4,U1],[A3,A1,U1],[A2[1],A2[2],U1]],  [[C3],[A3],[A2],[A2,A1]];
   elif dt[1]="B" then
      return [[D5],[D4,A1],[A3,C2],[B3,A1[1],A1[2]],[B4,U1]],  [[A1]];
   elif dt[1]="C" then
      return [[A4,U1],[C4,A1],[C3,C2]],  [[A1],[C2,A1]];
   elif dt[1]="D" then
      return [[A4,U1],[A3,A1[1],A1[2]],[D4,U1]],  [[B4],[B3,A1],[C2[1],C2[2]],[C2]];
   end if:
elif dt[2]=6 then
   if dt[1]="A" then
      return [[A5,U1],[A4,A1,U1],[A3,A2,U1]],  [[B3]];
   elif dt[1]="B" then
      return [[D6],[D5,A1],[D4,C2],[A3,B3],[B4,A1[1],A1[2]],[B5,U1]], [[A1]];
   elif dt[1]="C" then
      return [[A5,U1],[C5,A1],[C4,C2],[C3[1],C3[2]]], [[C2,A1],[A3,A1],[A1]];
   elif dt[1]="D" then
      return [[D5,U1],[A5,U1],[D4,A1[1],A1[2]],[A3[1],A3[2]]], 
             [[B5],[B4,A1],[B3,C2],[C3,A1],[A1[1],A1[2],A1[3]]];
   elif dt[1]="E" then
      return [[D5,U1],[A5,A1],[A2[1],A2[2],A2[3]]],  
             [[A2],[G2],[C4],[F4],[G2,A2]];
   end if:
elif dt[2]=7 then
   if dt[1]="A" then
      return [[A6,U1],[A5,A1,U1],[A4,A2,U1],[A3[1],A3[2],U1]],
             [[D4],[C4],[A3,A1]];
   elif dt[1]="B" then
      return [[D7],[D6,A1],[D5,C2],[D4,B3],[B4,A3],[B5,A1[1],A1[2]],[B6,U1]],
             [[A3],[C2,A1],[A1]];
   elif dt[1]="C" then
      return [[A6,U1],[C6,A1],[C5,C2],[C4,C3]],  [[B3,A1],[A1]];
   elif dt[1]="D" then
      return [[A6,U1],[D5,A1[1],A1[2]],[D4,A3],[D6,U1]], 
             [[C3],[C2],[G2],[B6],[B5,A1],[B4,C3],[B3[1],B3[2]]];
   elif dt[1]="E" then
      return [[E6,U1],[A7],[D6,A1],[A5,A2]],
             [[A1],[A1],[A2],[A1[1],A1[2]],[G2,A1],[F4,A1],[C3,G2]];
   end if:
elif dt[2]=8 then
   if dt[1]="A" then
      return [[A7,U1],[A6,A1,U1],[A5,A2,U1],[A4,A3,U1]],  
             [[B4],[A2[1],A2[1]]];
   elif dt[1]="B" then
      return [[D8],[D7,A1],[D6,C2],[D5,B3],[B4,D4],[B5,A3],
             [B6,A1[1],A1[2]],[B7,U1]],  [[A1]];
   elif dt[1]="C" then
      return [[A7,U1],[C7,A1],[C6,C2],[C5,C3],[C4[1],C4[1]]],
             [[C3],[D4,A1],[A1]];
   elif dt[1]="D" then
      return [[A7,U1],[D6,A1[1],A1[2]],[D5,A3],[D4[1],D4[2]],[D7,U1]],
             [[B4],[C4,A1],[C2[1],C2[2]],[B7],[B6,A1],[B5,C2],[B4,B3]];
   elif dt[1]="E" then
      return [[D8],[A4[1],A4[2]],[E6,A2],[E7,A1],[A8]], 
             [[A1],[A1],[A1],[F4,G2],[A2,A1],[C2]];
   end if;
end if: 
return "No infomation";
end proc:

###>proc:MaxSubGlist0
#% MaxSubGlist0(dt) => flat list of maximal subalgebras
MaxSubGlist0:=proc()
local dt,subglist:

if _npassed=0 then 
  error("Usage: MaxSubGlist0(dt)");
end if:

dt:=resolvedt(_passed[1]):
subglist:=MaxSubGlist(dt):
return [op(subglist[1]),op(subglist[2])];
end proc:

####>Symmetry breaking pattern search
###>Subalgebras
##>proc:MaySubG
#% MaySubG(dts1,dts0) => maplist : [dts0[1]=[dts1[1],..], dts0[2]=[dts1[3],..]]
MaySubG:=proc()
local dts1,dts0,
      ranklist1,ranklist0,maplist,maplist1,newmap,
      x,y,ranklis,i,j,k,dt0:

##% set parameters
if _npassed>=2 then
  dts1:=_passed[1]: dts0:=_passed[2]:
else
  error("usage: MaySubG(dts1::list, dts0::list)");
end if:
if not type(dts0,list) or not type(dts1,list) then
  error("usage: MaySubG(dts1::list, dts0::list)");
end if:
x:='x': y:='y':
dts1:=map(resolvedt,dts1):
dts1:=sort(dts1,(x,y)->x[2]>y[2]):
dts0:=map(resolvedt,dts0):
dts0:=sort(dts0,(x,y)->x[2]>y[2]):

#print("Check pt 1");
#print("dts0"=dts0, "dts1"=dts1);

##% make maplist
x:='x':
ranklist0:=map(x->x[2],dts0):
ranklist1:=map(x->x[2],dts1):
maplist:=[table(["res. rank"=ranklist0])]:
i:='i': for i from 1 to nops(ranklist1) do
  x:=ranklist1[i]:
  maplist1:=[]:
  y:='y': for y in maplist do # y= a table repr. a map
    ranklis:=y["res. rank"]:
    j:='j': for j from 1 to nops(ranklis) do
      if ranklis[j]>=x then
        newmap:=copy(y):
        dt0:=abbrevdt(dts0[j]):
        if assigned(newmap[dt0]) then
           newmap[dt0]:=[op(newmap[dt0]),abbrevdt(dts1[i])]:
        else
           newmap[dt0]:=[abbrevdt(dts1[i])]:
        end if:
        newmap["res. rank"]:=ranklis:
        newmap["res. rank"][j]:=ranklis[j]-x:
        maplist1:=[op(maplist1),copy(newmap)]:
      end if:
    end do:
  end do:
  maplist:=maplist1:
end do:

maplist1:=maplist:
i:='i': for i from 1 to nops(maplist1) do
  x:=maplist1[i]:
  j:='j': for j from 1 to nops(ranklist0) do
    if x["res. rank"][j]=ranklist0[j] then
       maplist[i][abbrevdt(dts0[j])]:=[]:
    end if:
  end do:
end do:
  
##% output
return maplist;
end proc:


###>proc:mkSBchain
#% mkSBchain(SBClist0[,sw]) => complete SB chain list
#% SBClist0= [ SBchains]: SBchain=[[A_0,[B_1,...]],,..[[C_1,[D_1,..]],[C_2,[D_2,...],...], [[D_1=[X1,.],D_2=[X_2,..]..], [X_3,U_1,...]] ]
mkSBchain:=proc()
local SBClist0,sw,
      lastp,SBClist1,SBC, algmap0,algmap1,algmap,SBC1,
      unbrokenlist,unblis1, newSBlist,addlist,oldSBlist,oldSBC, finished,
      i,j,k,x,x1,y,X,dt0,dts0,dts1,dts1r,dts2,tbl,oldunb,oldmaps,oldSBs:

##% set parameters
SBClist0:=_passed[1]: #% SB chain list
if  not type(SBClist0,list) then
  error("Usage: mkSBchain (SBClist0)");
end if:
sw:=0:
if _npassed>1 then
  sw:=1:
end if:

##% make one-step advanced SB chain list
SBClist1:=[]:
SBC:='SBC': for SBC in SBClist0 do  ## SBC do-loop
    # SBC=[[[E6,[D5,U1]]],[D5,[A3,A1,U1]],[[A3=[A2,A1],..],[A1,U1,..]]]
# print("SBC"=SBC);
 lastp:=nops(SBC):
 algmap0:=SBC[lastp][1]:
       # algmap0=[A3=[A2,A1,..],...]
 algmap1:=algmap0:
 algmap0:=[]:
 X:='X': for X in algmap1 do
   if nops(op(2,X))>1 or abbrevdt(resolvedt(op(1,X))[1..2])<>op(op(2,X)) then
      algmap0:=[op(algmap0),X]:
   end if:
 end do:
 if lastp>1 then
    oldSBC:=SBC[1..(lastp-1)]:   # oldSBC=[[E6,[D5,U1]],[D5,[A3,A1,U1]]
 else 
    oldSBC:=[]:
 end if:
 unbrokenlist:=SBC[lastp][2]: # unbrokenlist=[A1,U1,...]
 SBC1:=[op(oldSBC),[algmap0,unbrokenlist]]:
 newSBlist:=[[[],[[],[]]]]: 
 X:='X': for X in algmap0 do ## Algfactor X do-loop
   dt0:=resolvedt(op(1,X)):
   dt0:=dt0[1..2]:
   dts0:=op(2,X):
   addlist:=[]:
   dts1:='dts1': for dts1 in MaxSubGlist0(dt0) do  ## subalg dts1 do-loop
                             # A4 -> [[A3,U1],[A2,A1,U1],..]
     dts1r:=remove(has,dts1,U1):
     algmap:= MaySubG(dts0,dts1r): # algmap=[table("res. rank"=[0,1],A3=[A1,A2],...),..]
     if nops(algmap)>0 then 
        tbl:='tbl': for tbl in algmap do   ## algmap do-loop
                    # tbl=table("res. rank"=[0,1],A3=[A1,A2],..)
          algmap1:=[]: unblis1:=[]: 
          x:='x': for x in dts1 do
            x1:=x:
            if x<>U1 and nops(resolvedt(x))=3 then
                x1:=op(0,x):
            end if:
            if assigned(tbl[x]) then
               if x1<>op(tbl[x]) then
                 algmap1:=[op(algmap1),x=tbl[x]]:
              end if:
            else
              unblis1:=[op(unblis1),x1]:
            end if:
          end do:
          dts2:=[]:
          x:='x': for x in dts1 do
            if x=U1 then
              dts2:=[op(dts2),U1]:
            else
              dts2:=[op(dts2),abbrevdt(resolvedt(x)[1..2])]:
            end if:
          end do:
          addlist:=[op(addlist),[abbrevdt(dt0)=dts2,algmap1,unblis1]]:
        end do:   ## algmap do-loop
     end if:
   end do:  ## subalg dts1 do-loop
#  print("addlist"=addlist);
   oldSBlist:=newSBlist:
   newSBlist:=[]:
   x:='x': for x in oldSBlist do
     y:='y': for y in addlist do
       if nops(x[1])=0 then
          oldSBs:=NULL:
       else 
          oldSBs:=op(x[1]):
       end if:
       if nops(x[2][1])=0 then
         oldmaps:=NULL:
       else
         oldmaps:=op(x[2][1]):
       end if:
       if nops(x[2][2])=0 then
         oldunb:=NULL:
       else
         oldunb:=op(x[2][2]):
       end if:
       newSBlist:=[op(newSBlist),[[oldSBs,y[1]],[[oldmaps,op(y[2])],[oldunb,op(y[3])]] ]]:
     end do:
   end do:
 end do:   ## algfactor X do-loop
 if nops(algmap0)=0 then
    SBClist1:=[op(SBClist1),SBC1]:
 else
    oldSBlist:=newSBlist:
    newSBlist:=[]:
    x:='x': for x in oldSBlist do
        newSBlist:=[op(newSBlist),
          [op(oldSBC),x[1],[x[2][1],[op(unbrokenlist),op(x[2][2])]]] ]:
    end do: 
    if nops(newSBlist)>0 then
       SBClist1:=[op(SBClist1),op(newSBlist)]: 
    end if:
  end if:
end do:  ## SBC choice do-loop
finished:=true:
SBC:='SBC': for SBC in SBClist1 while finished do
  lastp:=nops(SBC):
  if nops(SBC[lastp][1])<>0 then
     finished:=false:
  end if:
end do:
if finished or sw=1 then
   return SBClist1;
else 
   mkSBchain(SBClist1);
end if:
end proc:


###>proc:SBpattern
#% SBpattern(dt0,dts[,outsw=0/1]) => Symmetry breaking chains from dt0 to dts=[dt1,...] 
#% dt0: a single Dynkin type
#% dts: a list of Dynkin types
SBpattern:=proc()
local dt0,dts,outsw,
      U1count,SBClist0,SBClist1,SBClist2,SBCtable,SBCnum,SBlist, unbrkn,resrank,
      SBC,lastp,x,finished,completeSBC:

if _npassed<2 or not type(_passed[2],list) then
  error("Usage: SBpattern(dt0,[dt1,dt2,...])");
end if:
dt0:=resolvedt(_passed[1]):
dts:=remove(has,_passed[2],U1):
U1count:=nops(_passed[2])-nops(dts):
dts:=map(resolvedt,dts):
outsw:=0:
if _npassed>=3 then
  outsw:=1:
end if:

SBClist0:=[[[[dt0=dts],[]]]]:
finished:=false:
while not finished do
  SBClist1:=mkSBchain(SBClist0,1):
  completeSBC:=true:
  SBC:='SBC': for SBC in SBClist1 while completeSBC do
    lastp:=nops(SBC):
    if nops(SBC[lastp][1])<>0 then
       completeSBC:=false:
    end if:
  end do:
  if completeSBC then 
     finished:=true:
  else 
     SBClist0:=SBClist1:
  end if:
end do:


SBClist2:={}: SBCtable:=table(): SBCnum:=0:
SBC:='SBC': for SBC in SBClist1 do
    lastp:=nops(SBC):
    SBlist:=SBC[1..(lastp-1)]:
    unbrkn:=SBC[lastp][2]:
    resrank:=0: 
    x:='x': for x in unbrkn do
      if x=U1 then
        resrank:=resrank+1:
      else
        resrank:=resrank+resolvedt(x)[2]:
      end if:
    end do:
    if resrank>=U1count then
       if not member(SBC,SBClist2) then
         SBCnum:=SBCnum+1:
         SBCtable[SBCnum]:=SBC:
         SBClist2:={op(SBClist2),SBC}:
         if outsw=1 then
         printf(" (%a) %a => [%a,%a]\n",SBCnum,SBlist,map(abbrevdt,dts),unbrkn);
         end if:
       end if:
     end if:
end do:
return eval(SBCtable):
end proc:
     
           
#####>Type A
####>SU(2)
#% Begin: Job_SU2 
if Job_SU2 then
printf("Groups with the info on the maximal quasi-semisimple subgroups and projection/embedding matrices preloaded:\n");
printf("  SU(2),");
#%
###>U1 embedding
embM[H][U1_SU2][0]:=Matrix([[1]]):
embM[S][U1_SU2][0]:=embM[H][U1_SU2][0].SR[SL2]:
embM[D][U1_SU2][0]:=embM[S][U1_SU2][0].CM[SL2]^(-1):
#%
end if:
#% End: Job_SU2
####>SU(3)
#% Begin: Job_SU3
if Job_SU3 then
printf("SU(3),");
#%
#%
###> Basis of su(3)
E:='E': H:='H': A:='A': S:='S':
BasisDefSU3:=[
  H[1]=I*(E[1,1]-E[3,3]), H[2]=I*(E[2,2]-E[3,3]),
  A[1]=E[2,3]-E[3,2], A[2]=E[3,1]-E[1,3], A[3]=E[1,2]-E[2,1],
  S[1]=I*(E[2,3]+E[3,2]), S[2]=I*(E[3,1]+E[1,3]), S[3]=I*(E[1,2]+E[2,1])]:
#
H:='H': WB:='WB':
WeylBasisDef:=[
H[1]=1/(3*I)*(2*H[1]-H[2]), H[2]=1/(3*I)*(H[1]+H[2]),
WB[1]=A[3]-I*S[3], WB[-1]=A[3]+I*S[3],
WB[2]=A[1]-I*S[1], WB[-2]=A[1]+I*S[1],
WB[3]=A[2]+I*S[2], WB[-3]=A[2]-I*S[2]
]:
#%
#% Matrix representation
HM[1]:=Matrix([[I,0,0],[0,0,0],[0,0,-I]]):
HM[2]:=Matrix([[0,0,0],[0,I,0],[0,0,-I]]):
AM[1]:=Matrix([[0,0,0],[0,0,1],[0,-1,0]]):
AM[2]:=Matrix([[0,0,-1],[0,0,0],[1,0,0]]):
AM[3]:=Matrix([[0,1,0],[-1,0,0],[0,0,0]]):
SM[1]:=Matrix([[0,0,0],[0,0,I],[0,I,0]]):
SM[2]:=Matrix([[0,0,I],[0,0,0],[I,0,0]]):
SM[3]:=Matrix([[0,I,0],[I,0,0],[0,0,0]]):
###>Subgroup reduction
#% Begin: Job_SU3_SubG
if Job_SU3_SubG then
#print("SU(3) subgroup reduction");
#%
#%
##>U1xU1->SU3
Y:=U1U1_SU3: Z:=0:
embM[H][Y][Z]:=DiagonalMatrix([1,1]):
embM[S][Y][Z]:=embM[H][Y][Z].SR[SL3]:
embM[D][Y][Z]:=embM[S][Y][Z].CM[SL3]^(-1):
#%
#%
##>SO(3)->SU(3)
#% Begin: Job_SU3_SubG_SO3
if Job_SU3_SubG_SO3 then
#print("Job: SU(3) -> SO(3) reduction");
#%
EmbDef:=[Hs[1]=H[1]]:
embM[H][SO3SU3]:=Matrix(1,2):
i:='i': for i from 1 to 2 do
  embM[H][SO3SU3][1,i]:=coeff(eval(Hs[1],EmbDef),H[i]):
end do:
embM[S][SO3SU3]:=SR[SO3]^(-1).embM[H][SO3SU3].SR[SL3]:
embM[D][SO3SU3]:=CM[SO3].embM[S][SO3SU3].CM[SL3]^(-1):
#DWL:=SubGreduction(A2,[1,1],B1,embM[D][SO3SU3]);
Y:=SO3_SU3: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][SO3SU3]:
end do:
Y:=SO3_SU3: Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][SO3SU3]:
end do:#%
end if:
#% End: Job_SU3_SubG_SO3
##>SU(2)xU(1) -> SU(3)
#% Begin: Job_SU3_SubG_SU2U1
if Job_SU3_SubG_SU2U1 then
#print("Job: SU(3) -> SU(2)xU(1)");
#%
#%
#>normal embedding 
embM[H][SU2SU3][n1]:=Matrix([[1,-1]]):
embM[S][SU2SU3][n1]:=SR[SL2]^(-1).embM[H][SU2SU3][n1].SR[SL3]:
embM[D][SU2SU3][n1]:=CM[SL2].embM[S][SU2SU3][n1].CM[SL3]^(-1):
#%
Qvec[H][SU3][n1]:=Matrix([[1,1]]):
Qvec[S][SU3][n1]:=Qvec[H][SU3][n1].SR[SL3]:
Qvec[D][SU3][n1]:=Qvec[S][SU3][n1].CM[SL3]^(-1):
#%
X:='X':  for X in [H,S,D] do
 Y:=SU2U1_SU3:
 embM[X][Y][n]:=Matrix(2,2):
 embM[X][Y][n][1,1..2]:=embM[X][SU2SU3][n1]:
 embM[X][Y][n][2,1..2]:=Qvec[X][SU3][n1]:
end do:
#>Canonical embedding
embM[D][SU2SU3][c1]:=Matrix([[1,1]]):
embM[S][SU2SU3][c1]:=CM[SL2]^(-1).embM[D][SU2SU3][c1].CM[SL3]:
embM[H][SU2SU3][c1]:=SR[SL2].embM[S][SU2SU3][c1].SR[SL3]^(-1):
#%
Qvec[H][SU3][c1]:=Matrix([[-1,2]]):
Qvec[S][SU3][c1]:=Qvec[H][SU3][c1].SR[SL3]:
Qvec[D][SU3][c1]:=Qvec[S][SU3][c1].CM[SL3]^(-1):
#%
X:='X':  for X in [H,S,D] do
 Y:=SU2U1_SU3:
 embM[X][Y][c]:=Matrix(2,2):
 embM[X][Y][c][1,1..2]:=embM[X][SU2SU3][c1]:
 embM[X][Y][c][2,1..2]:=Qvec[X][SU3][c1]:
end do:
#%
Y:=SU2U1_SU3: Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][SU2U1_SU3][n]:
end do:
#%
#DWL:=SubGU1reduction(A2, [1,0], [A1,U1], [embM[D][SU2SU3][0],Qvec[D][SU3][0]);
#%
end if:
#% End:Job_SU3_SubG_SU2U1
end if:
#% End: Job_SU3_SubG
end if:
#% End: Job_SU3
#%
####>SU(4)
#% Begin: Job_SU4
if Job_SU4 then
printf("SU(4),");
#%
###>RootSystem
SRbyONB:=[ alpha[1]=[sqrt(2),0,0], 
           alpha[2]=[-1/sqrt(2),sqrt(3/2),0],
           alpha[3]=[0,-sqrt(2/3),2/sqrt(3)]
         ]:
i:='i': j:='j':  CMI[SU4]:=Cmatrix(A3)^(-1):
f:='f':
DBbySR:=[seq(f[i]=add(CMI[SU4][i,j]*alpha[j],j=1..3),i=1..3)]:
DBbyONB:=eval(DBbySR,SRbyONB):
RootSystemp:=[alpha[1],alpha[2],alpha[3],
              alpha[1]+alpha[2],alpha[2]+alpha[3],
              alpha[1]+alpha[2]+alpha[3]
             ]:
RootSystemm:=-1*RootSystemp:
RVSp:=eval(RootSystemp,SRbyONB):
RVSm:=eval(RootSystemm,SRbyONB):
#% Dynkin basis
FVS[1]:=eval([f[1],-f[1]+f[2],-f[2]+f[3],-f[3]],DBbyONB):
FVS[2]:=eval([f[2],f[1]-f[2]+f[3],f[1]-f[3],-f[1]+f[3],-f[1]+f[2]-f[3],-f[2]], DBbyONB):
FVS[3]:=eval([f[3], f[2]-f[3],f[1]-f[2], -f[1]],DBbyONB):
#% projection to the SU3 hyperplane
nv:=eval(alpha[1]+2*alpha[2]+3*alpha[3],SRbyONB):  
x:='x': y:='y':
SU3HPeq:= [-(nv[1]*x+nv[2]*y)/nv[3]]:
#% plot
##>Weight Plot
#% Begin: Job_SU4_WeightPlot
if Job_SU4_WeightPlot  then
#%
colorlist:=[green,blue, yellow]:
i:='i': for i from 1 to 3 do
  v:='v':
  FVSplot[i]:=[seq(arrow(v,color=colorlist[i]),v in FVS[i])]:
end do:
v:='v': RSplot:=[seq(arrow(v,color=red),v in RVSp), 
                 seq(arrow(v,color=black),v in RVSm) ]:
i:='i': SRplot:=[seq(arrow(RVSp[i],color=red),i=1..3), 
                 seq(arrow(RVSm[i],color=black),i=1..3) ]:
SU3HP:=plot3d( op(SU3HPeq), x=-1..1, y=-1..1):
i:='i':
#print(display([op(SRplot),seq(op(FVSplot[i]),i=1..3), SU3HP]));
end if:
#% End; Job_SU4_WeightPlot
###>Subgroups
##>SU3xU1->SU4
#% Begin: Job_SU4_SubG_SU3
if Job_SU4_SubG_SU3 then
#print("Job: embedding of SU(3)xU(1) to SU(4)");
#%
#>standard embedding (Slansky)
# HA2[1]=HA3[1]-HA3[3], HA3[2]=HA3[2]-HA3[3]
  embM[H][SU3SU4][n1]:=Matrix([[1,0,-1],[0,1,-1]]):
  embM[S][SU3SU4][n1]:=SR[SL3]^(-1).embM[H][SU3SU4][n1].SR[SL4]:
  embM[D][SU3SU4][n1]:=CM[SL3].embM[S][SU3SU4][n1].CM[SL4]^(-1):
#%
  Qvec[H][SU4][n1]:=Matrix([[1,1,1]]):
  Qvec[S][SU4][n1]:=Qvec[H][SU4][n1].SR[SL4]:
  Qvec[D][SU4][n1]:=Qvec[S][SU4][n1].CM[SL4]^(-1):
#%
Y:=SU3U1_SU4: Z:=s:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU3SU4][n1],-1/3*Qvec[X][SU4][n1]>>:
end do:
#%
#>canonical embedding
embM[D][SU3SU4][c1]:=Matrix([[1,1,0],[0,0,1]]):
embM[S][SU3SU4][c1]:=CM[SL3]^(-1).embM[D][SU3SU4][c1].CM[SL4]:
embM[H][SU3SU4][c1]:=SR[SL3].embM[S][SU3SU4][c1].SR[SL4]^(-1):
#%
Qvec[D][SU4][c1]:=Matrix([[1,-2,-1]]):
Qvec[S][SU4][c1]:=Qvec[D][SU4][c1].CM[SL4]:
Qvec[H][SU4][c1]:=Qvec[S][SU4][c1].SR[SL4]^(-1): #=[1,-3,1]
#%
Y:=SU3U1_SU4:
X:='X': for X in [H,S,D] do
 embM[X][Y][c]:=<<embM[X][SU3SU4][c1],Qvec[X][SU4][c1]>>:
end do:
#%
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# DWL:=SubGrd(A3,[1,0,0],[A2,U1],[embM[D][SU3SU4,Qvec[D][SU4][1]]);
#%
end if:
#% End: Job_SU4_SubG_SU3
#%
##>SU2xSU2xU1->SU4
#% Begin: Job_SU4_SubG_SU2SU2U1
if Job_SU4_SubG_SU2SU2U1 then
#print("Job: embedding of SU(2)xSU(2)xU(1) to SU(4)");
#%
#>standard embedding
  embM[H][SU2SU4][1]:=Matrix([[1,-1,0]]):
  embM[H][SU2SU4][2]:=Matrix([[0,0,1]]):
  j:='j': for j from 1 to 2 do
    embM[S][SU2SU4][j]:=SR[SL2]^(-1).embM[H][SU2SU4][j].SR[SL4]:
    embM[D][SU2SU4][j]:=CM[SL2].embM[S][SU2SU4][j].CM[SL4]^(-1):
  end do:
  Qvec[H][SU4][2]:=Matrix([[1,1,-1]]):
  Qvec[S][SU4][2]:=Qvec[H][SU4][2].SR[SL4]:
  Qvec[D][SU4][2]:=Qvec[S][SU4][2].CM[SL4]^(-1):
#%
Y:=SU2SU2U1_SU4:
X:='X': for X in [H,S,D] do
 embM[X][Y][n]:=<<embM[X][SU2SU4][1],embM[X][SU2SU4][2],Qvec[X][SU4][2]>>:
end do:  
#%
#>canonical embedding
embM[D][SU2SU4][3]:=Matrix([[1,1,0]]):
embM[D][SU2SU4][4]:=Matrix([[0,1,1]]):
  j:='j': for j from 3 to 4 do
    embM[S][SU2SU4][j]:=CM[SL2]^(-1).embM[D][SU2SU4][j].CM[SL4]:
    embM[H][SU2SU4][j]:=SR[SL2].embM[S][SU2SU4][j].SR[SL4]^(-1):
end do:
Qvec[D][SU4][3]:=Matrix([[1,0,1]]):
Qvec[S][SU4][3]:=Qvec[D][SU4][3].CM[SL4]:
Qvec[H][SU4][3]:=Qvec[S][SU4][3].SR[SL4]^(-1):
#%
Y:=SU2SU2U1_SU4:
X:='X': for X in [H,S,D] do
 embM[X][Y][c]:=<<embM[X][SU2SU4][3],embM[X][SU2SU4][4],Qvec[X][SU4][3]>>:
end do:  
#%
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# DWL:=SubGrd(A3,[1,0,0],[A1,A1,U1],[embM[D][SU2SU4][1],embM[D][SU2SU4][2],Qvec[D][SU4][2]]);
#%
end if:
#% End: Job_SU4_SubG_SU2SU2U1
#%
##>SU2xSU2->SU4
#% Begin: Job_SU4_SubG_SU2SU2
if Job_SU4_SubG_SU2SU2 then
#print("Job: embedding of SU2xSU2 to SU4");
#%
#>canonical embedding
#% this is the Weyl trf w[0,0,0,1]w[0,0,1,0] of the normal embedding
embM[H][SU2SU4][5]:=Matrix([[1 , 1 , -1]]):
embM[H][SU2SU4][6]:=Matrix([[1 , -1, 1]]):
j:='j': for j from 5 to 6 do
embM[S][SU2SU4][j]:=SR[SL2]^(-1).embM[H][SU2SU4][j].SR[SL4]:
embM[D][SU2SU4][j]:=CM[SL2].embM[S][SU2SU4][j].CM[SL4]^(-1):
end do:
#%
Y:=SU2SU2_SU4:
X:='X': for X in [H,S,D] do
 embM[X][Y][c]:=<<embM[X][SU2SU4][5],embM[X][SU2SU4][6]>>:
end do:
#%
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][c]:
end do:
#%
# DWL:=SubGrd(A3,[1,0,0],[A1,A1],[embM[D][SU2SU4][5],embM[D][SU2SU4][6]]);
#%
end if:
#% End: Job_SU4_SubG_SU2SU2
#%
##>Sp2->SU4
#% Begin: Job_SU4_SubG_Sp2
if Job_SU4_SubG_Sp2 then
#print("Job: embedding of Sp2 to SU4");
#%
#% SO5->SO6
embM[H][SO5SO6]:=Matrix([[1,0,0],[0,1,0]]):
SR1[SO6]:=<<0,1,-1>|<1,-1,0>|<0,1,1>>:
embM[S][SO5SO6]:=SR[SO5]^(-1).embM[H][SO5SO6].SR[SO6];
embM[D][SO5SO6]:=CM[SO5].embM[S][SO5SO6].CM[SO6]^(-1);
#%
#% Sp2->SU4
embM[H][Sp2SU4]:=Matrix([[1,1],[1,-1]]).embM[H][SO5SO6].Matrix([[1,1,0],[1,0,-1],[0,1,-1]])^(-1):
embM[S][Sp2SU4]:=SR[Sp2]^(-1).embM[H][Sp2SU4].SR[SL4];
embM[D][Sp2SU4]:=CM[Sp2].embM[S][Sp2SU4].CM[SL4]^(-1);
#%
Y:=Sp2_SU4:Z:=n:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=embM[X][Sp2SU4]:
end do:  
#%
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
end if:
#% End: Job_SU4_SubG_Sp2
#%
end if:
#% End: Job_SU4
####>SU(5)
#% Begin: Job_SU5
if Job_SU5 then
printf("SU(5),");
#
###>Weyl trf
rv:='rv': WTD[A4]:= rv -> WeylTrf(A4,rv,"D"):
###>Root System of SU(5)
Cartan[SL5]:=[
  HA4[1]=E[1,1]-E[5,5], HA4[2]=E[2,2]-E[5,5], HA4[3]=E[3,3]-E[5,5],   
  HA4[4]=E[4,4]-E[5,5]
  ]:
###>Subgroup
#% Begin: Job_SU5_SubG
if Job_SU5_SubG then
#print("Job: Subgroups of SU(5)");
#%
##>SU(4)xU(1)
#% Begin: Job_SU5_SubG_SU4U1
if Job_SU5_SubG_SU4U1 then
#print("Job: SU(5) => SU(4)xU(1)");
#%
#>normal embedding
embM[H][SU4SU5][n1]:=Matrix([[1,0,0,-1],[0,1,0,-1],[0,0,1,-1]]):
embM[S][SU4SU5][n1]:=SR[SL4]^(-1).embM[H][SU4SU5][n1].SR[SL5]:
embM[D][SU4SU5][n1]:=CM[SL4].embM[S][SU4SU5][n1].CM[SL5]^(-1):
Qvec[H][SU5][SU4U1][n1]:=Matrix([[1,1,1,1]]):
Qvec[S][SU5][SU4U1][n1]:=Qvec[H][SU5][SU4U1][n1].SR[SL5]:
Qvec[D][SU5][SU4U1][n1]:=Qvec[S][SU5][SU4U1][n1].CM[SL5]^(-1):
Y:=SU4U1_SU5: Z:=n:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU4SU5][n1],Qvec[X][SU5][SU4U1][n1]>>:
end do:
#%
#>canonical embedding
embM[D][SU4SU5][c1]:=<<1,0,0>|<0,1,0>|<0,1,0>|<0,0,1>>:
embM[S][SU4SU5][c1]:=CM[SL4]^(-1).embM[D][SU4SU5][c1].CM[SL5]:
embM[H][SU4SU5][c1]:=SR[SL4].embM[S][SU4SU5][c1].SR[SL5]^(-1):
Qvec[D][SU5][SU4U1][c1]:=Matrix([[1,2,-2,-1]]):
Qvec[S][SU5][SU4U1][c1]:=Qvec[D][SU5][SU4U1][c1].CM[SL5]:
Qvec[H][SU5][SU4U1][c1]:=Qvec[S][SU5][SU4U1][c1].SR[SL5]^(-1):
Y:=SU4U1_SU5: Z:=c:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU4SU5][c1],Qvec[X][SU5][SU4U1][c1]>>:
end do:
#%
#>Slansky embedding
WTDs:=WTD[A4]([0,0,0,1]).WTD[A4]([0,1,1,0]).WTD[A4]([0,0,1,0]):
Z:=s:
embM[D][Y][Z]:=embM[D][Y][c].WTDs:
embM[S][Y][Z]:=DiagonalMatrix([CM[SL4],1])^(-1).embM[D][Y][Z].CM[SL5]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL4],1]).embM[S][Y][Z].SR[SL5]^(-1):
#%
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(A4,[1,0,0,0], [A3,U1],mmbM[D][SU4U1_SU5][s]):
#%
end if:
#% End: Job_SU5_SubG_SU4U1
#%
##>SU(3)xSU(2)xU(1)
#% Begin: Job_SU5_SubG_SU3SU2U1 
if Job_SU5_SubG_SU3SU2U1  then
#print("SU(5) => SU(3)xSU(2)xU(1)");
#%
#> SU(3)xSU(2)xU(1) standard model
if false then
quark[UL]:=[2/3, 1/2]: quark[DL]:=[-1/3,-1/2]: 
quark[URc]:=[-2/3,0]: quark[DRc]:=[1/3,0]:
lepton[nuL]:=[0,1/2]: lepton[eL]:=[-1,-1/2]: lepton[eRc]:=[1,0]:
x:='x': y:='y':
Ydef:=x*Qem+y*I3: 
ps:='ps': for ps in {UL, DL, URc, DRc} do
  Y[ps]:=eval(Ydef,{Qem=quark[ps][1],I3=quark[ps][2]}):
end do:
ps:='ps': for ps in {nuL, eL, eRc} do
  Y[ps]:=eval(Ydef,{Qem=lepton[ps][1],I3=lepton[ps][2]}):
end do:
sol:=solve({Y[UL]=Y[DL],Y[nuL]=Y[eL]},{x,y});
end if:
#%
#>normal embedding
#%
##% embedding: SU(3)->SU(5)
#% HA3[1]=HA4[1]-HA4[3], HA3[2]=HA4[2]-HA4[3]
  embM[H][SU3SU5][n]:=Matrix([[1,0,-1,0],[0,1,-1,0]]):
  embM[S][SU3SU5][n]:=SR[SL3]^(-1).embM[H][SU3SU5][n].SR[SL5]:
  embM[D][SU3SU5][n]:=CM[SL3].embM[S][SU3SU5][n].CM[SL5]^(-1):
##% embedding: SU(2) ->SU(5)
embM[H][SU2SU5][n]:=Matrix([[0,0,0,1]]):
embM[S][SU2SU5][n]:=SR[SL2]^(-1).embM[H][SU2SU5][n].SR[SL5]:
embM[D][SU2SU5][n]:=CM[SL2].embM[S][SU2SU5][n].CM[SL5]^(-1):
##% embedding: U(1) ->SU(5)
q:=-1/3:
Qvec[H][SU5][SU3SU2U1][n]:=q*Matrix([[2,2,2,-3]]):
Qvec[S][SU5][SU3SU2U1][n]:=Qvec[H][SU5][SU3SU2U1][n].SR[SL5]:
Qvec[D][SU5][SU3SU2U1][n]:=Qvec[S][SU5][SU3SU2U1][n].CM[SL5]^(-1):
#Qvec[D][SU5][SU3SU2U1]:=convert(Qvec[D][SU5][SU3SU2U1],list):
Y:=SU3SU2U1_SU5: Z:=n:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU3SU5][Z],embM[X][SU2SU5][Z],Qvec[X][SU5][SU3SU2U1][Z]>>:
end do:
#%
#%
#> Canonical embedding
##% embedding: SU(3)->SU(5)
embM[D][SU3SU5][c]:=<<1,0>|<1,0>|<0,1>|<0,1>>:
embM[S][SU3SU5][c]:=CM[SL3]^(-1).embM[D][SU3SU5][c].CM[SL5]:
embM[H][SU3SU5][c]:=SR[SL3].embM[S][SU3SU5][c].SR[SL5]^(-1):

#%
##%embedding: SU(2) ->SU(5)
embM[D][SU2SU5][c]:= Matrix([[0,1,1,0]]):
embM[S][SU2SU5][c]:=CM[SL2]^(-1).embM[D][SU2SU5][c].CM[SL5]:
embM[H][SU2SU5][c]:=SR[SL2].embM[S][SU2SU5][c].SR[SL5]^(-1):
#%
##%  embedding: U(1) ->SU(5)
QemD[SU5]:=1/3*[-1,2,1,1]:
I3D:=1/2*[0,1,1,0]:
Qvec[D][SU5][SU3SU2U1][c]:=Matrix(simplify(2*(QemD[SU5]-I3D))):
Qvec[S][SU5][SU3SU2U1][c]:=Qvec[D][SU5][SU3SU2U1][c].CM[SL5]:
Qvec[H][SU5][SU3SU2U1][c]:=Qvec[S][SU5][SU3SU2U1][c].SR[SL5]^(-1):
#%
Y:=SU3SU2U1_SU5: Z:=c:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU3SU5][Z],embM[X][SU2SU5][Z],Qvec[X][SU5][SU3SU2U1][Z]>>:
end do:
#%
#>Slansky type
Y:=SU3SU2U1_SU5: Z:=s:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU3SU5][c],embM[X][SU2SU5][c],Qvec[X][SU5][SU3SU2U1][c]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
#>decomposition: SU(3)xSU(2)xU(1)
# DWL:=SubG2U1reduction(A4,[0,0,0,1],[A2,A1],embM[D][SU3SU2U1_SU5][n]):
##print(DWL);
#>SU3xU1em
Y:=SU3U1em_SU5: Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=DiagonalMatrix([1,1,Matrix([[1/2, 1/2]])]).embM[X][SU3SU2U1_SU5][Z]:
end do:
#%
end if:
#% End: Job_SU5_SubG_SU3SU2U1
#%
##>SO(5)
#% Begin: Job_SU5_SubG_SO5
if Job_SU5_SubG_SO5 then
#print("SU(5) => SO(5)");
#%
embM[H][SO5SU5]:=Matrix([[1, 0 ,0,0],[0,1,0, -1]]):
embM[S][SO5SU5]:=SR[SO5]^(-1).embM[H][SO5SU5].SR[SL5]:
embM[D][SO5SU5]:=CM[SO5].embM[S][SO5SU5].CM[SL5]^(-1):
#%
Y:=SO5_SU5: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][SO5SU5]:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
#%
# SubGrd(A4,[1,0,0,0],[B2],embM[D][SO5SU5]):
#
end if:
#% End: Job_SU5_SubG_SO5
end if:
#% End: Job_SU5_SubG
end if:
#% End: Job_SU5
#%
####>SU(6)
#% Begin: Job_SU6
if Job_SU6 then
printf("SU(6),");
#
###>Weyl Trf
rv:='rv': WTD[A5]:= rv -> WeylTrf(A5,rv,"D"):
###>Root System of SU(6)
i:='i':
Cartan[SL6]:=[seq(HA[i]=E[i,i]-E[6,6],i=1..5)]:
###>Subalgebras
##>SU5xU1
#% Begin: Job_SU6_SubG_SU5U1
if Job_SU6_SubG_SU5U1 then
#print("Subalgbra: su5+u1");
#%
#>normal embedding =Slansky
embM[D][SU5SU6]:=<<0,0,0,0>|<1,0,0,0>|<0,1,0,0>|<0,0,1,0>|<0,0,0,1>>:
c:='c': Qv:=[5,c[1],c[2],c[3],c[4]];
i:='i': eqs:=[seq(DotProduct(Vector(Qv),Kmetric(A5).Transpose(embM[D][SU5SU6][i,1..5])),i=1..4)]:
i:='i': sol:=solve(eqs,[seq(c[i],i=1..4)])[1]:
Qvec[D][SU5SU6]:=Matrix([eval(Qv,sol)]):
Y:=SU5U1_SU6: Z:=n:
embM[D][Y][Z]:=<<embM[D][SU5SU6],Qvec[D][SU5SU6]>>:
embM[S][Y][Z]:=DiagonalMatrix([CM[SL5],1])^(-1).embM[D][Y][Z].CM[SL6]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL5],1]).embM[S][Y][Z].SR[SL6]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
end if:
#% End: Job_SU6_SubG_SU5U1
#%
##>SU4xSU2xU1
#% Begin: Job_SU6_SubG_SU4SU2U1
if Job_SU6_SubG_SU4SU2U1 then
#print("Subalgbra: su4+su2");
#%
#>normal embedding 
embM[D][SU4SU2SU6]:=<<1,0,0,0>|<0,1,0,0>|<0,0,1,0>|<0,0,0,0>|<0,0,0,1>>:
c:='c': Qv:=[c[1],c[2],c[3],4,c[4]];
i:='i': eqs:=[seq(DotProduct(Vector(Qv),Kmetric(A5).Transpose(embM[D][SU4SU2SU6][i,1..5])),i=1..4)]:
i:='i': sol:=solve(eqs,[seq(c[i],i=1..4)])[1]:
Qvec[D][SU4SU2SU6]:=Matrix([eval(Qv,sol)]):
Y:=SU4SU2U1_SU6: Z:=n:
embM[D][Y][Z]:=<<embM[D][SU4SU2SU6],Qvec[D][SU4SU2SU6]>>:
embM[S][Y][Z]:=DiagonalMatrix([CM[SL4],CM[SL2],1])^(-1).embM[D][Y][Z].CM[SL6]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL4],SR[SL2],1]).embM[S][Y][Z].SR[SL6]^(-1):
#%
#>Slansky embedding
WTDs:=WTD[A5]([1,1,1,0,0]).WTD[A5]([0,0,1,1,1]).WTD[A5]([0,1,1,0,0]).WTD[A5]([0,0,1,1,0]):
Y1:=SU4SU2U1_SU6: Z:=s:
embM[D][Y1][Z]:=embM[D][Y][n].WTDs:
embM[S][Y1][Z]:=DiagonalMatrix([CM[SL4],CM[SL2],1])^(-1).embM[D][Y1][Z].CM[SL6]:
embM[H][Y1][Z]:=DiagonalMatrix([SR[SL4],SR[SL2],1]).embM[S][Y1][Z].SR[SL6]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
end if:
#% End: Job_SU6_SubG_SU4SU2U1
#%
##>SU3xSU3xU1
#% Begin: Job_SU6_SubG_SU3SU3U1
if Job_SU6_SubG_SU3SU3U1 then
#print("Subalgbra: su3+su3+u1");
#%
#>normal embedding 
embM[D][SU3SU3SU6]:=<<1,0,0,0>|<0,1,0,0>|<0,0,0,0>|<0,0,1,0>|<0,0,0,1>>:
c:='c': Qv:=[c[1],c[2],3,c[3],c[4]];
i:='i': eqs:=[seq(DotProduct(Vector(Qv),Kmetric(A5).Transpose(embM[D][SU3SU3SU6][i,1..5])),i=1..4)]:
i:='i': sol:=solve(eqs,[seq(c[i],i=1..4)])[1]:
Qvec[D][SU3SU3SU6]:=Matrix([eval(Qv,sol)]):
Y:=SU3SU3U1_SU6: Z:=n:
embM[D][Y][Z]:=<<embM[D][SU3SU3SU6],Qvec[D][SU3SU3SU6]>>:
embM[S][Y][Z]:=DiagonalMatrix([CM[SL3],CM[SL3],1])^(-1).embM[D][Y][Z].CM[SL6]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL3],SR[SL3],1]).embM[S][Y][Z].SR[SL6]^(-1):
#%
#>Slansky embedding
WTDs:=WTD[A5]([1,1,0,0,0]).WTD[A5]([0,1,1,0,0]).WTD[A5]([0,0,0,1,0]):
Z:=s:
embM[D][Y][Z]:=embM[D][Y][n].WTDs:
embM[S][Y][Z]:=DiagonalMatrix([CM[SL3],CM[SL3],1])^(-1).embM[D][Y][Z].CM[SL6]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL3],SR[SL3],1]).embM[S][Y][Z].SR[SL6]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
end if:
#% End: Job_SU6_SubG_SU3SU3U1
#%
##>SU4 (special)
#% Begin: Job_SU6_SubG_SU4
if Job_SU6_SubG_SU4 then
#print("Subalgbra: su4 (special)");
#%
#>normal embedding 
Y:=SU4_SU6: Z:=n:
embM[S][Y][Z]:=<<0,1,0>|<1,0,0>|<-1,0,1>|<1,0,0>|<0,1,0>>:
embM[H][Y][Z]:=SR[SL4].embM[S][Y][Z].SR[SL6]^(-1):
embM[D][Y][Z]:=CM[SL4].embM[S][Y][Z].CM[SL6]^(-1):
#%
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
end if:
#% End: Job_SU6_SubG_SU4
##>SU3 (special)
#% Begin: Job_SU6_SubG_SU3
if Job_SU6_SubG_SU3 then
#print("Subalgbra: su3 (special)");
#%
#>normal embedding 
if false then
#% 6-dim irrep of su6 
WSA5:=WeightSystem(A5,[1,0,0,0,0],"S"):
hl:=WSA5["hl"]:  
i:='i': for i from 1 to 5 do
  HV6[i]:=[]: 
end do:
i:='i': for i from hl to -hl by -1 do
  x:='x': for x in WSA5[i] do
    j:='j': for j from 1 to 5 do
       HV6[j]:=[op(HV6[j]),op(2,x)[j]]:
    end do:
  end do:
end do:
#%
#% 6-dim irrep of su3 
WSA2:=WeightSystem(A2,[2,0],"S"):
hl:=WSA2["hl"]:  
HV3[1]:=[]: HV3[2]:=[]:
i:='i': for i from hl to -hl by -1 do
  x:='x': for x in WSA2[i] do
    j:='j': for j from 1 to 2 do
       HV3[j]:=[op(HV3[j]),op(2,x)[j]]:
    end do:
  end do:
end do:
#%
#% HV3 by HV6
i:='i': for i from 1 to 2 do
  c:='c': j:='j': k:='k':
  eqs:= [seq(add(c[k]*HV6[k][j],k=1..5)=HV3[i][j],j=1..6)]:
  k:='k': sol[i]:=solve(eqs,[seq(c[k],k=1..5)]):
end do:
end if:
#%
Y:=SU3_SU6: Z:=n:
embM[S][Y][Z]:=<<1,0>|<1,0>|<-1,1>|<1,0>|<0,1>>:
embM[D][Y][Z]:=CM[SL3].embM[S][Y][Z].CM[SL6]^(-1):
embM[H][Y][Z]:=SR[SL3].embM[S][Y][Z].SR[SL6]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
end if:
## End: Job_SU6_SubG_SU3
##>SU3xSU2 (special)
#% Begin: Job_SU6_SubG_SU3SU2
if Job_SU6_SubG_SU3SU2 then
#print("Subalgbra: su3+su2 (special)");
#%
#>normal embedding 
if false then
#% 6-dim irrep of su6 
WSA5:=WeightSystem(A5,[1,0,0,0,0],"S"):
hl:=WSA5["hl"]:  
i:='i': for i from 1 to 5 do
  HV6[i]:=[]: 
end do:
i:='i': for i from hl to -hl by -1 do
  x:='x': for x in WSA5[i] do
    j:='j': for j from 1 to 5 do
       HV6[j]:=[op(HV6[j]),op(2,x)[j]]:
    end do:
  end do:
end do:
#% 6-dim irrep of su3+su2
HV32[1]:=[2/3,-1/3,2/3,-1/3,-1/3,-1/3]:
HV32[2]:=[1/3,1/3,1/3,-2/3,1/3,-2/3]:
HV32[3]:=[1/2,1/2,-1/2,1/2,-1/2,-1/2]:
#%
#% HV32 by HV6
i:='i': for i from 1 to 3 do
  c:='c': j:='j': k:='k':
  eqs:= [seq(add(c[k]*HV6[k][j],k=1..5)=HV32[i][j],j=1..6)]:
  k:='k': sol[i]:=solve(eqs,[seq(c[k],k=1..5)]):
end do:
end if:
#%
Y:=SU3SU2_SU6: Z:=n:
embM[S][Y][Z]:=<<1,0,0>|<-1,0,1>|<1,1,-1>|<0,-1,1>|<0,1,0>>:
embM[D][Y][Z]:=DiagonalMatrix([CM[SL3],CM[SL2]]).embM[S][Y][Z].CM[SL6]^(-1):
embM[H][Y][Z]:=DiagonalMatrix([SR[SL3],SR[SL2]]).embM[S][Y][Z].SR[SL6]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
end if:
#%End: Job_SU6_SubG_SU3SU2
#%
##>Sp3 (special)
#% Begin: Job_SU6_SubG_Sp3
if Job_SU6_SubG_Sp3 then
#print("Subalgbra: sp3 (special)");
#%
#>normal embedding 
Y:=Sp3_SU6: Z:=n:
embM[S][Y][Z]:=<<1,0,0>|<0,1,0>|<0,0,1>|<0,1,0>|<1,0,0>>:
embM[H][Y][Z]:=SR[Sp3].embM[S][Y][Z].SR[SL6]^(-1):
embM[D][Y][Z]:=CM[Sp3].embM[S][Y][Z].CM[SL6]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
end if:
#% End: Job_SU6_SubG_Sp3
end if:
#% End: Job_SU6
#####>Type B
####>SO(5)
#% Bein: Job_SO5
if Job_SO5 then
printf("SO(5), ");
#%
###>Algebra 
#% Basis
H:='H': E:='E':
CartanBasis[SO5]:=[
  H[1]=-I*(E[1,2]-E[2,1]),  H[2]=-I*(E[3,4]-E[4,3])
  ]:
###>Subgroups
##>SU2xU1
#% Begin: Job_SO5_SubG_SU2U1
if Job_SO5_SubG_SU2U1 then
#print("Job: Regular subgroup SU2xU1");
#%
#>normal Embedding 
embM[H][SU2SO5][n1]:=Matrix([[0,2]]):
embM[S][SU2SO5][n1]:=SR[SL2]^(-1).embM[H][SU2SO5][n1].SR[SO5]:
embM[D][SU2SO5][n1]:=CM[SL2].embM[S][SU2SO5][n1].CM[SO5]^(-1):
Qvec[H][SO5][n1]:=Matrix([[2,0]]):
Qvec[S][SO5][n1]:=Qvec[H][SO5][n1].SR[SO5]:
Qvec[D][SO5][n1]:=Qvec[S][SO5][n1].CM[SO5]^(-1):
#%
Y:=SU2U1_SO5: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU2SO5][n1],Qvec[X][SO5][n1]>>:
end do:
#%
#%
#>canonical embedding
embM[D][SU2SO5][c1]:=Matrix([[2,1]]):
embM[S][SU2SO5][c1]:=CM[SL2]^(-1).embM[D][SU2SO5][c1].CM[SO5]:
embM[H][SU2SO5][c1]:=SR[SL2].embM[S][SU2SO5][c1].SR[SO5]^(-1):
#5
Qvec[D][SO5][c1]:=Matrix([[0,1]]):
Qvec[S][SO5][c1]:=Qvec[D][SO5][c1].CM[SO5]:
Qvec[H][SO5][c1]:=Qvec[S][SO5][c1].SR[SO5]^(-1):
#%
Y:=SU2U1_SO5: Z:=c:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU2SO5][c1],Qvec[X][SO5][c1]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%#%
# SubGrd(B2,[1,0],[B1,U1],[embM[D][SU2SO5][1],Qvec[D][SO5]]);
#%
end if:
#% End: Job_SO5_SubG_SU2U1
#%
##>SU2xSU2
#% Begin: Job_SO5_SubG_SU2SU2
if Job_SO5_SubG_SU2SU2 then
#print("Job: Regular subgroup SU2xSU2=SO4");
#%
#% Embedding
embM[H][SU2SO5][2]:=Matrix([[1,-1]]):
embM[H][SU2SO5][3]:=Matrix([[1,1]]):
i:='i': for i from 2 to 3 do
  embM[S][SU2SO5][i]:=SR[SL2]^(-1).embM[H][SU2SO5][i].SR[SO5]:
  embM[D][SU2SO5][i]:=CM[SL2].embM[S][SU2SO5][i].CM[SO5]^(-1):
end do:
#%
Y:=SU2SU2_SO5: Z:=n1:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU2SO5][2],embM[X][SU2SO5][3]>>:
end do:
#%
# SubGrdm(B2,[1,0],[A1,A1],embM[D][SU2SU2_SO5][n1]);
#% the choice based on the extended Dynkin diagram
embM[S][SU2SO5][5]:=Matrix([[1,-1/2]]):
embM[S][SU2SO5][6]:=Matrix([[0,-1/2]]):
i:='i': for i from 5 to 6 do
  embM[H][SU2SO5][i]:=SR[SL2].embM[S][SU2SO5][i].SR[SO5]^(-1):
  embM[D][SU2SO5][i]:=CM[SL2].embM[S][SU2SO5][i].CM[SO5]^(-1):
end do:
#%
Y:=SU2SU2_SO5: Z:=n2:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU2SO5][5],embM[X][SU2SO5][6]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n1]:
end do:
#%
# SubGrdm(B2,[1,0],[A1,A1],embM[D][SU2SU2_SO5][n2]);
#%
end if:
#% End:Job_SO5_SubG_SU2SU2
#%
##>SU2 (special)
#% Begin: Job_SO5_SubG_SU2max
if Job_SO5_SubG_SU2max then
#print("Job: Special maximal subgroup SU2");
#%
#% Embedding
embM[H][SU2SO5][4]:=Matrix([[4,2]]):
embM[S][SU2SO5][4]:=SR[SL2]^(-1).embM[H][SU2SO5][4].SR[SO5]:
embM[D][SU2SO5][4]:=CM[SL2].embM[S][SU2SO5][4].CM[SO5]^(-1):
#%
Y:=SU2_SO5: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][SU2SO5][4]:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(B2,[1,0],[A1],embM[D][SU2_SO5][n]);
#%
end if:
#% End:Job_SO5_SubG_SU2max
#%
end if:
#% End; Job_SO5
#%
####>SO(7)
#% Bein: Job_SO7
if Job_SO7 then
printf("SO(7),");
#%
###>Algebra 
#% Basis
H:='H': E:='E':
CartanBasis[SO7]:=[
  H[1]=-I*(E[1,2]-E[2,1]),  H[2]=-I*(E[3,4]-E[4,3]), H[3]=-I*(E[5,6]-E[6,5])
  ]:
#%
###>Subgroups
##>SU4 (regular)
#% Begin: Job_SO7_SubG_SU4
if Job_SO7_SubG_SU4 then
#print("Job: Regular maximal subgroup SU4");
#%
#>Normal embedding
embM[H][SU4SO7][n]:=Matrix([[1,1,0],[1,0,-1],[0,1,-1]]):
embM[S][SU4SO7][n]:=SR[SL4]^(-1).embM[H][SU4SO7][n].SR[SO7]:
embM[D][SU4SO7][n]:=CM[SL4].embM[S][SU4SO7][n].CM[SO7]^(-1):
#%
Y:=SU4_SO7: Z:=n:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=embM[X][SU4SO7][n]:
end do:
#
#>Slansky type embedding
#% n=>s : special trf alpha1 <-> alpha3 in SU4
Y:=SU4_SO7: Z:=s:
embM[S][Y][Z]:=Matrix([[0,0,1],[0,1,0],[1,0,0]]).embM[S][Y][n]:
embM[H][Y][Z]:=SR[SL4].embM[S][Y][Z].SR[SO7]^(-1):
embM[D][Y][Z]:=CM[SL4].embM[S][Y][Z].CM[SO7]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(B3,[1,0,0],[A3],[embM[D][SU4SO7]]);
#%
end if:
#% End:Job_SO7_SubG_SU4
##>SO5xU1 (regular)
#% Begin: Job_SO7_SubG_SO5U1
if Job_SO7_SubG_SO5U1 then
#print("Job: Regular maximal subgroup SO5xU1");
#%
#% Embedding
embM[H][SO5SO7]:=Matrix([[0,1,0],[0,0,1]]):
embM[S][SO5SO7]:=SR[SO5]^(-1).embM[H][SO5SO7].SR[SO7]:
embM[D][SO5SO7]:=CM[SO5].embM[S][SO5SO7].CM[SO7]^(-1):
Qvec[H][SO7]:=Matrix([[1,0,0]]):
Qvec[S][SO7]:=Qvec[H][SO7].SR[SO7]:
Qvec[D][SO7]:=Qvec[S][SO7].CM[SO7]^(-1):
#%
Y:=SO5U1_SO7: Z:=n:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SO5SO7],Qvec[X][SO7]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
%
# SubGrdm(B3,[1,0,0],[B2,U1],embM[D][SO5U1_SO7]);
#%
end if:
#% End:Job_SO7_SubG_SO5U1
#%
##>SU2xSU2xSU2 (regular)
#% Begin: Job_SO7_SubG_SU2SU2SU2
if Job_SO7_SubG_SU2SU2SU2 then
#print("Regular maximal subgroup SU2xSU2xSU2");
#%
#% Embedding
embM[H][SU2SO7][1]:=Matrix([[1,-1,0]]):
embM[H][SU2SO7][2]:=Matrix([[1,1,0]]):
embM[H][SU2SO7][3]:=Matrix([[0,0,2]]):
i:='i': for i from 1 to 3 do
  embM[S][SU2SO7][i]:=SR[SL2]^(-1).embM[H][SU2SO7][i].SR[SO7]:
  embM[D][SU2SO7][i]:=CM[SL2].embM[S][SU2SO7][i].CM[SO7]^(-1):
end do:
#%
Y:=SU2SU2SU2_SO7: Z:=n:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU2SO7][1],embM[X][SU2SO7][2],embM[X][SU2SO7][3]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
%
#%
# SubGrd(B3,[1,0,0],[A1,A1,A1],[embM[D][SU2SO7][1],embM[D][SU2SO7][2],embM[D][SU2SO7][3]]);
#%
end if:
#% End:Job_SO7_SubG_SU4
#%
##>G2 (special)
#% Begin: Job_SO7_SubG_G2
if Job_SO7_SubG_G2 then
#print("Special subgroup G2");
#%
Y:=G2_SO7: Z:=n:
embM[S][Y][Z]:=Matrix([[0,1,0],[1,0,1]]):
embM[H][Y][Z]:=SR[G2].embM[S][Y][Z].SR[SO7]^(-1):
embM[D][Y][Z]:=CM[G2].embM[S][Y][Z].CM[SO7]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(B3,[1,0,0],[G2],embM[D][G2_SO7][n]):
#%
end if:
#% End: Job_SO7_SubG_G2
#%
end if:
#% End; Job_SO7
#%
####>SO(9)
#% Begin: Job_SO9
if Job_SO9 then
printf("SO(9), ");
#
###>Weyl Trf
rv:='rv':  WTD[B4]:= rv -> WeylTrf(B4,rv,"D"):
###>Algebra 
#% Basis
H:='H': E:='E':
CartanBasis[SO9]:=[
  H[1]=-I*(E[1,2]-E[2,1]),  H[2]=-I*(E[3,4]-E[4,3]), H[3]=-I*(E[5,6]-E[6,5]), H[4]=-I*(E[7,8]-E[8,7])
  ]:
#%
###>Subgroups
##>SO8
#% Begin: Job_SO9_SubG_SO8
if Job_SO9_SubG_SO8 then
#print("Embedding so8 to so9");
#%
#>normal embedding
Y:=SO8_SO9: Z:=c:
embM[H][Y][Z]:=Matrix([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]):
embM[S][Y][Z]:=SR[SO8]^(-1).embM[H][Y][Z].SR[SO9]:
embM[D][Y][Z]:=CM[SO8].embM[S][Y][Z].CM[SO9]^(-1):
#%
#>Slansky embedding
Y:=SO8_SO9: Z:=s:
embM[D][Y][Z]:=embM[D][Y][c].WeylTrf(B4,[0,0,0,1],"D"):
embM[S][Y][Z]:=embM[S][Y][c].WeylTrf(B4,[0,0,0,1],"S"):
embM[H][Y][Z]:=SR[SO8].embM[S][Y][Z].SR[SO9]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][s]:
end do:
#%
# SubGrd(B4,[1,0,0,0],[D4],[embM[D][SO8_SO9][s]]):
#%
end if:
#% End: Job_SO9_SubG_SO8
#%
##>SO7xU1
#% Begin: Job_SO9_SubG_SO7U1
if Job_SO9_SubG_SO7U1 then
#print("Embedding so7+u1 to so9");
#%
embM[H][SO7SO9]:=Matrix([[0,1,0,0],[0,0,1,0],[0,0,0,1]]):
embM[S][SO7SO9]:=SR[SO7]^(-1).embM[H][SO7SO9].SR[SO9]:
embM[D][SO7SO9]:=CM[SO7].embM[S][SO7SO9].CM[SO9]^(-1):
#%
Qvec[H][SO9]:=Matrix([[1,0,0,0]]):
Qvec[S][SO9]:=Qvec[H][SO9].SR[SO9]:
Qvec[D][SO9]:=Qvec[S][SO9].CM[SO9]^(-1):
#%
Y:=SO7U1_SO9: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<< embM[X][SO7SO9], Qvec[X][SO9]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(B4,[1,0,0,0],[B3,U1],[embM[D][SO7SO9],Qvec[D][SO9]]):
#%
end if:
#% End: Job_SO9_SubG_SO7U1
#%
##>SU4xSU2
#% Begin: Job_SO9_SubG_SU4SU2
if Job_SO9_SubG_SU4SU2 then
#print("Embedding su4+su2 to so9");
#%
#>normal embedding
Z:=n:
embM[H][SU4SO9][Z]:=Matrix([[1,1,0,0],[1,0,-1,0],[0,1,-1,0]]):
embM[S][SU4SO9][Z]:=SR[SL4]^(-1).embM[H][SU4SO9][Z].SR[SO9]:
embM[D][SU4SO9][Z]:=CM[SL4].embM[S][SU4SO9][Z].CM[SO9]^(-1):
#%
embM[H][SU2SO9][1]:=Matrix([[0,0,0,2]]):
embM[S][SU2SO9][1]:=SR[SL2]^(-1).embM[H][SU2SO9][1].SR[SO9]:
embM[D][SU2SO9][1]:=CM[SL2].embM[S][SU2SO9][1].CM[SO9]^(-1):
#%
Y:=SU4SU2_SO9:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU4SO9][Z],embM[X][SU2SO9][1]>>:
end do:
#%
#>Slansky embedding
WTDs:=WTD[B4]([1,1,1,2]).WTD[B4]([0,0,1,2]).WTD[B4]([0,1,1,2]).WTD[B4]([0,0,1,0]).WTD[B4]([0,0,0,1]):
Y1:=SU4SU2_SO9: Z:=s:
embM[D][Y1][Z]:=embM[D][Y][n].WTDs:
embM[S][Y1][Z]:=DiagonalMatrix([CM[SL4],CM[SL2]])^(-1).embM[D][Y1][Z].CM[SO9]:
embM[H][Y1][Z]:=DiagonalMatrix([SR[SL4],SR[SL2]]).embM[S][Y1][Z].SR[SO9]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(B4,[1,0,0,0],[A3,A1],embM[D][SU4SU2_SO9][n]):
#%
end if:
#% End: Job_SO9_SubG_SU4SU2
#%
##> SU2xSU2xSp2 -> SO9
#% Begin: Job_SO9_SubG_Sp2SU2SU2
if Job_SO9_SubG_Sp2SU2SU2 then
#print("Embedding SU2xSU2xSp2 to SO9");
#% SU2-> SO9
embM[H][SU2SO9][2]:=Matrix([[1,-1,0,0]]):
embM[H][SU2SO9][3]:=Matrix([[1,1,0,0]]):
i:='i': for i from 2 to 3 do
  embM[S][SU2SO9][i]:=SR[SL2]^(-1).embM[H][SU2SO9][i].SR[SO9]:   
  embM[D][SU2SO9][i]:=CM[SL2].embM[S][SU2SO9][i].CM[SO9]^(-1):
end do:
#%
embM[H][Sp2SO9]:=Matrix([[0,0,1,1],[0,0,1,-1]]):
embM[S][Sp2SO9]:=SR[Sp2]^(-1).embM[H][Sp2SO9].SR[SO9]:
embM[D][Sp2SO9]:=CM[Sp2].embM[S][Sp2SO9].CM[SO9]^(-1):
#%
Y:=SU2SU2Sp2_SO9: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU2SO9][2],embM[X][SU2SO9][3],embM[X][Sp2SO9]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%

#%
# SubGrd(B4,[1,0,0,0],[A1,A1,C2],[embM[D][SU2SO9][2],embM[D][SU2SO9][3],embM[D][Sp2SO9]]);
#%
end if:
#% End: Job_SO9_SubG_Sp2SU2SU2
#%
##>max SU2
#% Begin: Job_SO9_SubG_SU2
if Job_SO9_SubG_SU2 then
#print("Embedding max su2 to so9");
#%
embM[H][SU2SO9][4]:=Matrix([[8,6,4,2]]):
embM[S][SU2SO9][4]:=SR[SL2]^(-1).embM[H][SU2SO9][4].SR[SO9]:
embM[D][SU2SO9][4]:=CM[SL2].embM[S][SU2SO9][4].CM[SO9]^(-1):
#%
Y:=SU2_SO9: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][SU2SO9][4]:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
#%
# SubGrd(B4,[1,0,0,0],[A1],[embM[D][SU2SO9][4]]):
#%
end if:
#% End: Job_SO9_SubG_SU2
#%
##> max SU2xSU2 -> SO9
#% Begin: Job_SO9_SubG_SU2SU2
if Job_SO9_SubG_SU2SU2 then
#print("Embedding max SU2xSU2 to SO9");
#% SU2-> SO9
embM[H][SU2SO9][5]:=Matrix([[2,2,0,2]]):
embM[H][SU2SO9][6]:=Matrix([[2,0,2,-2]]):
i:='i': for i from 5 to 6 do
  embM[S][SU2SO9][i]:=SR[SL2]^(-1).embM[H][SU2SO9][i].SR[SO9]:   
  embM[D][SU2SO9][i]:=CM[SL2].embM[S][SU2SO9][i].CM[SO9]^(-1):
end do:
#%
Y:=SU2SU2_SO9: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU2SO9][5],embM[X][SU2SO9][6]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
#%
# SubGrd(B4,[1,0,0,0],[A1,A1],[embM[D][SU2SO9][5],embM[D][SU2SO9][6]]);
#%
end if:
#% End: Job_SO9_SubG_SU2SU2
#%
#%
end if:
#% End: Job_SO9
#%
#####>Type C
####>Sp3
#% Begin: Job_Sp3
if Job_Sp3 then
printf("Sp3,");
#%
###>Algebra
##% Basis of sp3
CartanBasis[Sp3]:=[
  H[1]=E[1,1]-E[6,6], H[2]=E[2,2]-E[5,5], H[3]=E[3,3]-E[4,4]
  ]:
###>Weyl Trf
rv:='rv':  WTD[C3]:= rv -> WeylTrf(C3,rv,"D"):
#%
###>Subalgebras
##>SU3xU1
#% Begin: Job_Sp3_SubG_SU3U1
if Job_Sp3_SubG_SU3U1 then
#print("Subalgebra su3+u1 of sp3");
#%
embM[H][SU3Sp3]:=Matrix([[1,-1,0],[0,1,-1]]):
AM[SL3]:=Matrix([[2, -1],[-1,2]]):
embM[S][SU3Sp3]:=AM[SL3]^(-1).embM[H][SU3Sp3].SR[Sp3]:
embM[D][SU3Sp3]:=CM[SL3].embM[S][SU3Sp3].CM[Sp3]^(-1):
#%
Qvec[H][Sp3]:=Matrix([[1,1,1]]):
Qvec[S][Sp3]:=Qvec[H][Sp3].SR[Sp3]:
Qvec[D][Sp3]:=Qvec[S][Sp3].CM[Sp3]^(-1):
#%
Y:=SU3U1_Sp3: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU3Sp3],Qvec[X][Sp3]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(C3,[1,0,0],[A2,U1],[embM[D][SU3Sp3],Qvec[D][Sp3]]);
#%
end if:
#% End: Job_Sp3_SubG_SU3U1
##>Sp2xSU2
#% Begin: Job_Sp3_SubG_Sp2SU2
if Job_Sp3_SubG_Sp2SU2 then
#print("Subalgebra sp2+su2 of sp3");
#>normal embedding (1)
Y:=Sp2SU2_Sp3:
embM[S][Y][n1]:=<<-1,-1/2,-1/2>|<1,0,0>|<0,1,0>>:
embM[D][Y][n1]:=DiagonalMatrix([CM[Sp2],CM[SL2]]).embM[S][Y][n1].CM[Sp3]^(-1):
embM[H][Y][n1]:=DiagonalMatrix([SR[Sp2],SR[SL2]]).embM[S][Y][n1].SR[Sp3]^(-1):
#%
#>normal embedding (2)
embM[H][Sp2Sp3][n2]:=Matrix([[0,1,0],[0,0,1]]):
embM[S][Sp2Sp3][n2]:=SR[Sp2]^(-1).embM[H][Sp2Sp3][n2].SR[Sp3]:
embM[D][Sp2Sp3][n2]:=CM[Sp2].embM[S][Sp2Sp3][n2].CM[Sp3]^(-1):
#%
embM[H][SU2Sp3][1]:=Matrix([[1,0,0]]):
embM[S][SU2Sp3][1]:=SR[SL2]^(-1).embM[H][SU2Sp3][1].SR[Sp3]:
embM[D][SU2Sp3][1]:=CM[SL2].embM[S][SU2Sp3][1].CM[Sp3]^(-1):
#%
Y:=Sp2SU2_Sp3: Z:=n2:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][Sp2Sp3][Z],embM[X][SU2Sp3][1]>>:
end do:
#% embM[S][Sp2SU2_Sp3][n2]=WeylTrf(C2,"S").embM[S][Sp2SU2_Sp3][n1]
#>canonical embedding
Y:=Sp2SU2_Sp3:Z:=c:
WTDs:=WTD[C3]([1,0,0]).WTD[C3]([0,1,0]):
embM[D][Y][Z]:=embM[D][Y][n2].WTDs:
embM[S][Y][Z]:=DiagonalMatrix([CM[Sp2],CM[SL2]])^(-1).embM[D][Y][n2].CM[Sp3]:
embM[H][Y][Z]:=DiagonalMatrix([SR[Sp2],SR[SL2]]).embM[S][Y][Z].SR[Sp3]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(C3,[1,0,0],[C2,A1],embM[D][Sp2SU2_Sp3][c]):
end if:
#% End: Job_Sp3_SubG_Sp2SU2
##>SU2
#% Begin: Job_Sp3_SubG_SU2
if Job_Sp3_SubG_SU2 then
#print("Maximal subalgebra su2 of sp3");
#%
embM[H][SU2Sp3][2]:=Matrix([[5,3,1]]):
embM[S][SU2Sp3][2]:=SR[SL2]^(-1).embM[H][SU2Sp3][2].SR[Sp3]:
embM[D][SU2Sp3][2]:=CM[SL2].embM[S][SU2Sp3][2].CM[Sp3]^(-1):
#%
Y:=SU2_Sp3: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][SU2Sp3][2]:
end do:
#%
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(C3,[1,0,0],[A1],[embM[D][SU2Sp3][2]]):
end if:
#% End:  Job_Sp3_SubG_SU2
##>SU2xSU2
#% Begin: Job_Sp3_SubG_SU2SU2
if Job_Sp3_SubG_SU2SU2 then
#print("Maximal subalgebra su2+su2 of sp3");
#%
embM[H][SU2Sp3][3]:=Matrix([[1,1,1]]):
embM[H][SU2Sp3][4]:=Matrix([[2,0,-2]]):
i:='i': for i from 3 to 4 do
  embM[S][SU2Sp3][i]:=SR[SL2]^(-1).embM[H][SU2Sp3][i].SR[Sp3]:
  embM[D][SU2Sp3][i]:=CM[SL2].embM[S][SU2Sp3][i].CM[Sp3]^(-1):
end do:
#%
Y:=SU2SU2_Sp3: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU2Sp3][3],embM[X][SU2Sp3][4]>>:
end do:
#%
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(C3,[1,0,0],[A1,A1],[embM[D][SU2Sp3][3],embM[D][SU2Sp3][4]]):
end if:
#% End:  Job_Sp3_SubG_SU2SU2
end if:
#% End: Job_Sp3
####>Sp4
#% Begin: Job_Sp4
if Job_Sp4 then
printf("Sp4,");
#%
###>Algebra
##% Basis of sp4
CartanBasis[Sp4]:=[
  H[1]=E[1,1]-E[8,8], H[2]=E[2,2]-E[7,7], H[3]=E[3,3]-E[6,6], 
  H[4]=E[4,4]-E[5,5]
  ]:
###>Subalgebras
##>SU4xU1
#% Begin: Job_Sp4_SubG_SU4U1
if Job_Sp4_SubG_SU4U1 then
#print("Subalgebra su4+u1 of sp4");
#%
embM[H][SU4Sp4]:=Matrix([[1, 0, 0, -1],[0,1,0,-1],[0,0,1,-1]]):
embM[S][SU4Sp4]:=SR[SL4]^(-1).embM[H][SU4Sp4].SR[Sp4]:
embM[D][SU4Sp4]:=CM[SL4].embM[S][SU4Sp4].CM[Sp4]^(-1):
#%
Qvec[H][Sp4]:=Matrix([[1,1,1,1]]):
Qvec[S][Sp4]:=Qvec[H][Sp4].SR[Sp4]:
Qvec[D][Sp4]:=Qvec[S][Sp4].CM[Sp4]^(-1):
#%
Y:=SU4U1_Sp4: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU4Sp4],Qvec[X][Sp4]>>:
end do:
#%
# SubGrd(C4,[1,0,0,0],[A3,U1],[embM[D][SU4Sp4], Qvec[D][Sp4]]):
#%
end if:
#% End: Job_Sp4_SubG_SU4U1
#%
##>Sp3xSU2
#% Begin: Job_Sp4_SubG_Sp3SU2
if Job_Sp4_SubG_Sp3SU2 then
#print("Subalgebra sp3+su2 of sp4");
#%
embM[H][SU2Sp4][1]:=Matrix([[1, 0, 0, 0]]):
embM[S][SU2Sp4][1]:=SR[SL2]^(-1).embM[H][SU2Sp4][1].SR[Sp4]:
embM[D][SU2Sp4][1]:=CM[SL2].embM[S][SU2Sp4][1].CM[Sp4]^(-1):
#%
embM[H][Sp3Sp4]:=Matrix([[0, 1, 0, 0],[0,0,1,0],[0,0,0,1]]):
embM[S][Sp3Sp4]:=SR[Sp3]^(-1).embM[H][Sp3Sp4].SR[Sp4]:
embM[D][Sp3Sp4]:=CM[Sp3].embM[S][Sp3Sp4].CM[Sp4]^(-1):
#%
Y:=Sp3SU2_Sp4: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][Sp3Sp4],embM[X][SU2Sp4][1]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(C4,[1,0,0,0],[SU2,Sp3],[embM[D][SU2Sp4][1], embM[D][Sp3Sp4]]):
#%
end if:
#% End: Job_Sp4_SubG_So3SU2
#%
##>Sp2xSp2
#% Begin: Job_Sp4_SubG_Sp2Sp2
if Job_Sp4_SubG_Sp2Sp2 then
#print("Subalgebra sp2+sp2 of sp4");
#%
embM[H][Sp2Sp4][1]:=Matrix([[1, 0, 0, 0],[0,1,0,0]]):
embM[H][Sp2Sp4][2]:=Matrix([[0, 0, 1, 0],[0,0,0,1]]):
i:='i': for i from 1 to 2 do
  embM[S][Sp2Sp4][i]:=SR[Sp2]^(-1).embM[H][Sp2Sp4][i].SR[Sp4]:
  embM[D][Sp2Sp4][i]:=CM[Sp2].embM[S][Sp2Sp4][i].CM[Sp4]^(-1):
end do:
#%
Y:=Sp2Sp2_Sp4: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][Sp2Sp4][1],embM[X][Sp2Sp4][2]>>:
end do:
#%
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(C4,[1,0,0,0],[Sp2,Sp2],[embM[D][Sp2Sp4][1], embM[D][Sp2Sp4][2]]):
#%
end if:
#% End: Job_Sp2Sp2
#%
##>max SU2
#% Begin: Job_Sp4_SubG_SU2
if Job_Sp4_SubG_SU2 then
#print("Subalgebra max su2 of sp4");
#%
embM[H][SU2Sp4][2]:=Matrix([[7,5,3,1]]):
embM[S][SU2Sp4][2]:=SR[SL2]^(-1).embM[H][SU2Sp4][2].SR[Sp4]:
embM[D][SU2Sp4][2]:=CM[SL2].embM[S][SU2Sp4][2].CM[Sp4]^(-1):
#%
Y:=SU2_Sp4: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][SU2Sp4][2]:
end do:
#%
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(C4,[1,0,0,0],[SU2],[embM[D][SU2Sp4][2]]):
#%
end if:
#% End: Job_Sp4_SubG_SU2
#%
#%
##>max SU2xSU2xSU2
#% Begin: Job_Sp4_SubG_SU2SU2SU2
if Job_Sp4_SubG_SU2SU2SU2 then
#print("Subalgebra max su2+su2+su2 of sp4");
#%
embM[H][SU2Sp4][3]:=Matrix([[1,1,1,-1]]):
embM[H][SU2Sp4][4]:=Matrix([[1,1,-1,1]]):
embM[H][SU2Sp4][5]:=Matrix([[1,-1,1,1]]):
i:='i': for i from 3 to 5 do
  embM[S][SU2Sp4][i]:=SR[SL2]^(-1).embM[H][SU2Sp4][i].SR[Sp4]:
  embM[D][SU2Sp4][i]:=CM[SL2].embM[S][SU2Sp4][i].CM[Sp4]^(-1):
end do:
#%
Y:=SU2SU2SU2_Sp4: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU2Sp4][3],embM[X][SU2Sp4][4],embM[X][SU2Sp4][5]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(C4,[1,0,0,0],[A1,A1,A1],[seq(embM[D][SU2Sp4][i],i=3..5)]):
#%
end if:
#% End: Job_Sp4_SubG_SU2
#%
end if:
#% End: Job_Sp4
####>Sp5
#% Begin: Job_Sp5
if Job_Sp5 then
printf("Sp5,");
#%
###>Subalgebras
#% incomplete!!
##>SU5xU1
#% Begin: Job_Sp5_SubG_SU5U1
if Job_Sp5_SubG_SU5U1 then
#print("Subalgebra su5+u1 of sp4");
#%
embM[D][SU5Sp5]:=<1,0,0,0|0,1,0,0|0,0,1,0|0,0,0,1|0,0,0,0>:
embM[S][SU5Sp5]:=CM[SL5]^(-1).embM[D][SU5Sp5].CM[Sp5]:
embM[H][SU5Sp5]:=SR[SL5].embM[S][SU5Sp5].SR[Sp5]^(-1):
#%
q := 'q'; Qv := Vector(5, symbol = q): q[1]:=1:
eqs := convert(embM[D][SU5Sp5].(1/Gmetric(C5)).Qv, list):
Qvec[D][Sp5] :=Transpose(eval(Qv, solve(eqs, convert(Qv, list)[2..5])[1])):
Qvec[S][Sp5]:=Qvec[D][Sp5].CM[Sp5]:
Qvec[H][Sp5]:=Qvec[S][Sp5].SR[Sp5]^(-1):
#%
Y:=SU5U1_Sp5: Z:=n:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU5Sp5],Qvec[X][Sp5]>>:
end do:
#%
# SubGrdm(C5,[1,0,0,0,0],[A4,U1],embM[D][SU5U1_Sp5],1):
#%
end if:
#% End: Job_Sp5_SubG_SU5U1
#%
#%
end if:
#% End: Job_Sp5
#####>Type D
####>SO(8)
#% Begin: Job_SO8
if Job_SO8 then
printf("SO(8),");
#
###>Weyl trf
rv:='rv':  WTD[D4]:= rv -> WeylTrf(D4,rv,"D"):
###>Algebra
#% Basis of so(8)
E:='E': H:='H': A:='A': S:='S':
CartanBasis[SO8]:=[
  H0[1]=I*(E[1,2]-E[2,1]), H0[2]=I*(E[3,4]-E[4,3]),
  H0[3]=I*(E[5,6]-E[6,5]), H0[4]=I*(E[7,8]-E[8,7])
  ]:
#
###>Subalgebras
#%
##>U(4)=U(1)xSU(4) 
#% Begin: Job_SO8_SubG_SU4U1
if Job_SO8_SubG_SU4U1 then
#print("Embedding SU4xU1 to SO8");
#%
CartanBasis[SU4]:=[
   H1[1]=E[1,1]-E[4,4], H1[2]=E[2,2]-E[4,4], H1[3]=E[3,3]-E[4,4]
   ]:
#>normal embedding
Hembeding:=[H1[1]=H0[1]-H0[4], H1[2]=H0[2]-H0[4], H1[3]=H0[3]-H0[4],
   H1[4]=H0[1]+H0[2]+H0[3]+H0[4] ]:
embM[H][SU4SO8]:=Matrix(3,4):
i:='i': for i from 1 to 3 do
  j:='j': for j from 1 to 4 do
     embM[H][SU4SO8][i,j]:=coeff(eval(H1[i],Hembeding),H0[j]):
  end do:
end do:
embM[S][SU4SO8]:=SR[SL4]^(-1).embM[H][SU4SO8].SR[SO8]:
embM[D][SU4SO8]:=CM[SL4].embM[S][SU4SO8].CM[SO8]^(-1):

Qvec[H][SO8]:=Matrix(1,4):
i:='i': for i from 1 to 4 do
  Qvec[H][SO8][1,i]:=coeff(eval(H1[4],Hembeding),H0[i]):
end do:
Qvec[S][SO8]:=Qvec[H][SO8].SR[SO8]:
Qvec[D][SO8]:=Qvec[S][SO8].CM[SO8]^(-1):
#%
Y:=SU4U1_SO8: Z:=n:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU4SO8],Qvec[X][SO8]>>:
end do:
#%
#>Slansky embedding
WTDs:=WTD[D4]([1,1,0,1]).WTD[D4]([0,0,0,1]).WTD[D4]([0,1,0,0]):
Y:=SU4U1_SO8: Z:=s:
embM[D][Y][Z]:=embM[D][Y][n].WTDs:
embM[S][Y][Z]:=DiagonalMatrix([CM[SL4]^(-1),1]).embM[D][Y][Z].CM[SO8]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL4],1]).embM[S][Y][Z].SR[SO8]^(-1):
#> default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
##% Irr. decomposition
#DWL:=SubGrd(D4,[0,1,0,0],[A3,U1],[embM[D][SU4SO8],Qvec[D][SO8]]);
#%
end if:
#% End: Job_SO8_SubG_SU4U1
#%
##> SU2xSU2xSU2xSU2 -> SO8
#% Begin: Job_SO8_SubG_SU2SU2SU2SU2
if Job_SO8_SubG_SU2SU2SU2SU2 then
#print("Embedding SU2xSU2xSU2xSU2 to SO8");
#%#%
embM[H][SU2SO8][1]:=Matrix([[1,-1,0,0]]):
embM[H][SU2SO8][2]:=Matrix([[1,1,0,0]]):
embM[H][SU2SO8][3]:=Matrix([[0,0,1,-1]]):
embM[H][SU2SO8][4]:=Matrix([[0,0,1,1]]):
i:='i': for i from 1 to 4 do
  embM[S][SU2SO8][i]:=SR[SL2]^(-1).embM[H][SU2SO8][i].SR[SO8]:
  embM[D][SU2SO8][i]:=CM[SL2].embM[S][SU2SO8][i].CM[SO8]^(-1):
end do:
#%
Y:=SU2SU2SU2SU2_SO8: Z:=n:
X:='X': for X in [H,S,D] do
 i:='i':
 embM[X][Y][Z]:=<<seq(embM[X][SU2SO8][i],i=1..4)>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(D4,[1,0,0,0],[A1,A1,A1,A1],[seq(embM[S]D[SU2SO8][i],i=1..4)]);
#%
end if:
#% End: Job_SO8_SubG_SU2SU2SU2SU2
##> SO7
#% Begin: Job_SO8_SubG_SO7
if Job_SO8_SubG_SO7 then
#print("Embedding SO7to SO8");
#%
#>natural embedding
Z:=n:
embM[H][SO7SO8][Z]:=Matrix([[1,0,0,0],[0,1,0,0],[0,0,1,0]]):
embM[S][SO7SO8][Z]:=SR[SO7]^(-1).embM[H][SO7SO8][Z].SR[SO8]:
embM[D][SO7SO8][Z]:=CM[SO7].embM[S][SO7SO8][Z].CM[SO8]^(-1):
#%
Y:=SO7_SO8: Z:=n:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=embM[X][SO7SO8][Z]:
end do:
#%
#>Slansky embedding
#% non-equivalent to the natural one
#% n=>s: special trf alpha1<->alpha4 in SO8
Y:=SO7_SO8: Z:=s:
embM[S][Y][Z]:=embM[S][Y][n].Matrix([[0,0,0,1],[0,1,0,0],[0,0,1,0],[1,0,0,0]]):
embM[H][Y][Z]:=SR[SO7].embM[S][Y][Z].SR[SO8]^(-1):
embM[D][Y][Z]:=CM[SO7].embM[S][Y][Z].CM[SO8]^(-1):
Z:=1:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][s]:
end do:
#>Slansky embedding(dual)
#% non-equivalent to the above two
#% n=>2: special trf alpha1<->alpha3 in SO8
Y:=SO7_SO8: Z:=2:
embM[S][Y][Z]:=embM[S][Y][n].Matrix([[0,0,1,0],[0,1,0,0],[1,0,0,0],[0,0,0,1]]):
embM[H][Y][Z]:=SR[SO7].embM[S][Y][Z].SR[SO8]^(-1):
embM[D][Y][Z]:=CM[SO7].embM[S][Y][Z].CM[SO8]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(D4,[1,0,0,0],[B3],[embM[D][SO7SO8]]);
#%
end if:
#% End: Job_SO8_SubG_SO7
#%
##> SU3
#% Begin: Job_SO8_SubG_SU3
if Job_SO8_SubG_SU3 then
#print("Embedding SU3 to SO8");
#%
Y:=SU3_SO8: Z:=n:
embM[H][Y][Z]:=Matrix([[2,1,1,0],[1,-1,2,0]]):
embM[S][Y][Z]:=SR[SL3]^(-1).embM[H][Y][Z].SR[SO8]:
embM[D][Y][Z]:=CM[SL3].embM[S][Y][Z].CM[SO8]^(-1):
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(D4,[1,0,0,0],[A2],[embM[D][SU3SO8]]);
#%
end if:
#% End: Job_SO8_SubG_SU3
#%
##> Sp2xSU2
#% Begin: Job_SO8_SubG_Sp2SU2
if Job_SO8_SubG_Sp2SU2 then
#print("Embedding Sp2xSU2 to SO8");
#% SU2-> SO8
embM[H][SU2SO8][5]:=Matrix([[2,0,0,0]]):
embM[S][SU2SO8][5]:=SR[SL2]^(-1).embM[H][SU2SO8][5].SR[SO8]:
embM[D][SU2SO8][5]:=CM[SL2].embM[S][SU2SO8][5].CM[SO8]^(-1):
#%
embM[H][SO5SO8]:=Matrix([[0,1,0,0],[0,0,1,0]]):
embM[S][SO5SO8]:=SR[SO5]^(-1).embM[H][SO5SO8].SR[SO8]:
embM[D][SO5SO8]:=CM[SO5].embM[S][SO5SO8].CM[SO8]^(-1):
#%
embM[H][Sp2SO8]:=Matrix([[1,1],[1,-1]]).embM[H][SO5SO8]:
embM[S][Sp2SO8]:=SR[Sp2]^(-1).embM[H][Sp2SO8].SR[SO8]:
embM[D][Sp2SO8]:=CM[Sp2].embM[S][Sp2SO8].CM[SO8]^(-1):
#%
#%
Y:=Sp2SU2_SO8: Z:=n:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][Sp2SO8],embM[X][SU2SO8][5]>>:
end do:
#% default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrd(D4,[1,0,0,0],[A1,B2],embM[D][Sp2SU2_SO8][n]);
#%
end if:
#% End: Job_SO8_SubG_Sp2SU2
#%
end if:
#% End: Job_SO8
#%

####>SO(10)
#% Begin: Job_SO10
if Job_SO10 then
printf("SO(10),");
#%
###>Weyl trf
rv:='rv':  WTD[D5]:= rv -> WeylTrf(D5,rv,"D"):
#%
###>Algebra
#% Basis of so(10)
E:='E': H:='H': A:='A': S:='S':
CartanBasis[SO8]:=[
  H0[1]=-I*(E[1,2]-E[2,1]), H0[2]=-I*(E[3,4]-E[4,3]),
  H0[3]=-I*(E[5,6]-E[6,5]), H0[4]=-I*(E[7,8]-E[8,7]),
  H0[5]=-I*(E[9,10]-E[10,9])
  ]:
#
###>Subalgebras
#%
##>SO8xU(1)
#% Begin: Job_SO10_SubG_SO8U1
if Job_SO10_SubG_SO8U1 then
#print("The subalgebra so8+u1");
#%
#>normal embedding
embM[D][SO8SO10][n]:=<<0,0,0,0>|<1,0,0,0>|<0,1,0,0>|<0,0,1,0>|<0,0,0,1>>:
embM[S][SO8SO10][n]:=CM[SO8]^(-1).embM[D][SO8SO10][n].CM[SO10]:
embM[H][SO8SO10][n]:=SR[SO8].embM[S][SO8SO10][n].SR[SO10]^(-1):
#%
b:='b':
beta:=[1,b[1],b[2],b[3],b[4]]:
i:='i':  eqs:=[seq(DotProduct(Matrix(beta),Gmatrix(D5).embM[D][SO8SO10][n][i,1..5]), i=1..4)]:
i:='i': sol:=solve(eqs,[seq(b[i],i=1..4)]):
Qvec[H][SO10][n1]:=Matrix([eval(beta,sol[1])]):
Qvec[S][SO10][n1]:=Qvec[H][SO10][n1].SR[SO10]:
Qvec[D][SO10][n1]:=Qvec[S][SO10][n1].CM[SO10]^(-1):
#%
Y:=SO8U1_SO10: 
X:='X': for X in [H,S,D] do
 embM[X][Y][n]:=<<embM[X][SO8SO10][n],Qvec[X][SO10][n1]>>:
end do:
#%
#>canonical embedding
embM[H][SO8SO10][c]:=Matrix([[1,0,0,0,0],[0,1,0,0,0],[0,0,1,0,0],[0,0,0,1,0]]):
embM[S][SO8SO10][c]:=SR[SO8]^(-1).embM[H][SO8SO10][c].SR[SO10]:
embM[D][SO8SO10][c]:=CM[SO8].embM[S][SO8SO10][c].CM[SO10]^(-1):
#%
Qvec[H][SO10][c1]:=Matrix([[0,0,0,0,1]]):
Qvec[S][SO10][c1]:=Qvec[H][SO10][c1].SR[SO10]:
Qvec[D][SO10][c1]:=Qvec[S][SO10][c1].CM[SO10]^(-1):
#%
Y:=SO8U1_SO10: 
X:='X': for X in [H,S,D] do
 embM[X][Y][c]:=<<embM[X][SO8SO10][c],Qvec[X][SO10][c1]>>:
end do:

#%
#>Slansky embedding
WTDs:=WTD[D5]([1,0,0,0,0]).WTD[D5]([0,1,2,1,1]).WTD[D5]([0,1,1,1,0]).WTD[D5]([0,1,1,0,0]):
Y:=SO8U1_SO10: Z:=s:
embM[D][Y][Z]:=embM[D][Y][c].WTDs:
embM[S][Y][Z]:=DiagonalMatrix([CM[SO8],1])^(-1).embM[D][Y][Z].CM[SO10]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SO8],1]).embM[S][Y][Z].SR[SO10]^(-1):
#> default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:#%
# SubGrd(D5,[1,0,0,0,0],[D4,U1],embM[D][SO8U1_SO10][c]);
end if:
#% End: Job_SO10_SubG_SO8U1
#%
##>U5=SU5xU1
#% Begin: Job_SO10_SubG_SU5U1
if Job_SO10_SubG_SU5U1 then
#print("Job: Flipped SU5");
#%
Y:=SU5U1_SO10:
#%
#>normal embedding (1)
embM[H][SU5SO10][n1]:=Matrix([[1,0,0,0,-1],[0,1,0,0,-1],[0,0,1,0,-1],[0,0,0,1,-1]]):
embM[S][SU5SO10][n1]:=SR[SL5]^(-1).embM[H][SU5SO10][n1].SR[SO10]:
embM[D][SU5SO10][n1]:=CM[SL5].embM[S][SU5SO10][n1].CM[SO10]^(-1):
#%
Qvec[H][SO10][2]:=Matrix([[2,2,2,2,2]]):
Qvec[S][SO10][2]:=Qvec[H][SO10][2].SR[SO10]:
Qvec[D][SO10][2]:=Qvec[S][SO10][2].CM[SO10]^(-1):
#%
X:='X': for X in [H,S,D] do
 embM[X][Y][n1]:=<<embM[X][SU5SO10][n1],Qvec[X][SO10][2]>>:
end do:
#%
#>normal embedding (2)
#% complex conjugation of normal embedding (1)
Z:=n2:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=-embM[X][Y][n1]:
end do:
#>normal embedding (3)
#% The sign of the U(1) charge is reversed from(1).
Z:=n3:
embM[D][SU5SO10][Z]:=<<1,0,0,0>|<1,0,0,1>|<1,0,1,0>|<0,1,0,0>|<0,0,1,0>>:
embM[S][SU5SO10][Z]:=CM[SL5]^(-1).embM[D][SU5SO10][Z].CM[SO10]:
embM[H][SU5SO10][Z]:=SR[SL5].embM[S][SU5SO10][Z].SR[SO10]^(-1):
#%
Qvec[D][SO10][3]:=Matrix([[-2,0,2,1,-1]]):
Qvec[S][SO10][3]:=Qvec[D][SO10][3].CM[SO10]:
Qvec[H][SO10][3]:=Qvec[S][SO10][3].SR[SO10]^(-1):
#%
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU5SO10][Z],Qvec[X][SO10][3]>>:
end do:
#>canonical embedding (Slansky type)
#% equivalent to Z=n2.
Z:=c:
WTDs:=WTD[D5]([1,2,2,1,1]).WTD[D5]([1,0,0,0,0]).WTD[D5]([0,0,1,1,1]):
embM[D][Y][Z]:=DiagonalMatrix([-1,-1,-1,-1,1]).embM[D][Y][n3].WTDs:  # complex conjugation of (3)  in the SU(5) sector only.  
embM[S][Y][Z]:=DiagonalMatrix([CM[SL5],1])^(-1).embM[D][Y][Z].CM[SO10]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL5],1]).embM[S][Y][Z].SR[SO10]^(-1):
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n1]:
end do:
#%
#%
#SubGrd(D5,[1,0,0,0,0],[A4,U1],embM[D][SU5U1_SO10][c]);
#%
end if:
#% End: Job_SO10_SubG_SU5U1
#%
##>SU4xSU2xSU2-> SO10
#% Begin: Job_SO10_SubG_SU4SU2SU2
if Job_SO10_SubG_SU4SU2SU2 then
#print("Job: Embedding of SU4xSU2xSU2 to SO10");
#%
#>normal embedding 
Y:=SU4SU2SU2_SO10: Z:=n1:
embM[S][Y][Z]:=<<1,0,0,0,0>|<0,1,0,0,0>|<-1/2,-1,-1/2,-1/2,-1/2>|<0,0,0,1,0>|<0,0,0,0,1>>:
embM[D][Y][Z]:=DiagonalMatrix([CM[SL4],CM[SL2],CM[SL2]]).embM[S][Y][Z].CM[SO10]^(-1):
embM[H][Y][Z]:=DiagonalMatrix([SR[SL4],SR[SL2],SR[SL2]]).embM[S][Y][Z].SR[SO10]^(-1):
#>natural embedding 
embM[H][SU2SO10][n1]:=Matrix([[0,0,0,1,-1]]):
embM[H][SU2SO10][n2]:=Matrix([[0,0,0,1,1]]):
i:='i': for i from 1 to 2 do
  X:=cat(n,i):
  embM[S][SU2SO10][X]:=SR[SL2]^(-1).embM[H][SU2SO10][X].SR[SO10]:
  embM[D][SU2SO10][X]:=CM[SL2].embM[S][SU2SO10][X].CM[SO10]^(-1):
end do:
#%
embM[H][SU4SO10][n2]:=Matrix([[1,1,0,0,0],[1,0,-1,0,0],[0,1,-1,0,0]]):
embM[S][SU4SO10][n2]:=SR[SL4]^(-1).embM[H][SU4SO10][n2].SR[SO10]:
embM[D][SU4SO10][n2]:=CM[SL4].embM[S][SU4SO10][n2].CM[SO10]^(-1):
#%
Y:=SU4SU2SU2_SO10: Z:=n2:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU4SO10][Z],embM[X][SU2SO10][n1],embM[X][SU2SO10][n2]>>:
end do:
#%
#%#%
#>canonical embedding
embM[D][SU2SO10][c1]:=Matrix([[0,0,1,1,0]]):
embM[D][SU2SO10][c2]:=Matrix([[0,0,1,0,1]]):
i:='i': for i from 1 to 2 do
  X:=cat(c,i):
  embM[S][SU2SO10][X]:=CM[SL2]^(-1).embM[D][SU2SO10][X].CM[SO10]:
  embM[H][SU2SO10][X]:=SR[SL2].embM[S][SU2SO10][X].SR[SO10]^(-1):
end do:
#%
embM[D][SU4SO10][c]:=<<0,1,0>|<1,0,1>|<1,0,1>|<1,0,0>|<1,0,0>>:
embM[S][SU4SO10][c]:=CM[SL4]^(-1).embM[D][SU4SO10][c].CM[SO10]:
embM[H][SU4SO10][c]:=SR[SL4].embM[S][SU4SO10][c].SR[SO10]^(-1):
#%
Y:=SU4SU2SU2_SO10: Z:=c:
X:='X': for X in [H,S,D] do
 embM[X][Y][Z]:=<<embM[X][SU4SO10][Z],embM[X][SU2SO10][c1],embM[X][SU2SO10][c2]>>:
end do:
#%
#>Slansky embedding
WTDs:=WTD[D5]([1,1,1,0,0]).WTD[D5]([0,1,1,0,0]).WTD[D5]([0,0,0,1,0]):
Y1:=SU4SU2SU2_SO10:  Z:=s:
embM[D][Y1][Z]:=embM[D][Y][c].WTDs:
embM[S][Y1][Z]:=DiagonalMatrix([CM[SL4],CM[SL2],CM[SL2]])^(-1).embM[D][Y1][Z].CM[SO10]:
embM[H][Y1][Z]:=DiagonalMatrix([SR[SL4],SR[SL2],SR[SL2]]).embM[S][Y1][Z].SR[SO10]^(-1):
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n1]:
end do:
#%
#SubGrd(D5,[1,0,0,0,0], [A1,A1,A3],embM[D][SU4SU2SU2_SO10][c]);
#%
end if:
#% End: Job_SO10_SubG_SU4SU2SU2
#%
##>SO9
#% Begin: Job_SO10_SubG_SO9
if Job_SO10_SubG_SO9 then
#print("The subalgebra so9");
#%
#>normal and canonical embedding
embM[H][SO9SO10]:=Matrix([[1,0,0,0,0],[0,1,0,0,0],[0,0,1,0,0],[0,0,0,1,0]]):
embM[S][SO9SO10]:=SR[SO9]^(-1).embM[H][SO9SO10].SR[SO10]:
embM[D][SO9SO10]:=CM[SO9].embM[S][SO9SO10].CM[SO10]^(-1):
#%
Y:=SO9_SO10:
X:='X': for X in [H,S,D] do
 embM[X][Y][c]:=embM[X][SO9SO10]:
end do:
#%
#>Slansky embedding
WTDs:=WTD[D5]([1,0,0,0,0]).WTD[D5]([0,1,1,1,1]).WTD[D5]([0,0,1,0,1]).WTD[D5]([0,0,1,0,0]):
Y:=SO9_SO10: Z:=s:
embM[D][Y][Z]:=embM[D][Y][c].WTDs:
embM[S][Y][Z]:=CM[SO9]^(-1).embM[D][Y][Z].CM[SO10]:
embM[H][Y][Z]:=SR[SO9].embM[S][Y][Z].SR[SO10]^(-1):
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][c]:
end do:
#%
# SubGrd(D5,[1,0,0,0,0],[B4],embM[D][SO9_SO10][s]);
end if:
#% End: Job_SO10_SubG_SO9
#%
##>SO7xSU2-> SO10
#% Begin: Job_SO10_SubG_SO7SU2
if Job_SO10_SubG_SO7SU2 then
#print("Job: Embedding of SO7xSU2 to SO10");
#%
#>normal embedding
embM[H][SU2SO10][3]:=Matrix([[2,0,0,0,0]]):
embM[S][SU2SO10][3]:=SR[SL2]^(-1).embM[H][SU2SO10][3].SR[SO10]:
embM[D][SU2SO10][3]:=CM[SL2].embM[S][SU2SO10][3].CM[SO10]^(-1):
#%
embM[H][SO7SO10]:=Matrix([[0,0,1,0,0],[0,0,0,1,0],[0,0,0,0,1]]):
embM[S][SO7SO10]:=SR[SO7]^(-1).embM[H][SO7SO10].SR[SO10]:
embM[D][SO7SO10]:=CM[SO7].embM[S][SO7SO10].CM[SO10]^(-1):
#%
Y:=SO7SU2_SO10: 
X:='X': for X in [H,S,D] do
 embM[X][Y][n]:=<<embM[X][SO7SO10],embM[X][SU2SO10][3]>>:
end do:
#%
#>canonical embedding
WTDs:=WTD[D5]([1,1,0,0,0]).WTD[D5]([0,1,1,0,0]).WTD[D5]([0,0,1,1,0]).WTD[D5]([0,0,0,1,0]):
Y:=SO7SU2_SO10: 
embM[D][Y][c]:=embM[D][Y][n].WTDs:
embM[S][Y][c]:=DiagonalMatrix([CM[SO7],CM[SL2]])^(-1).embM[D][Y][c].CM[SO10]:
embM[H][Y][c]:=DiagonalMatrix([SR[SO7],SR[SL2]]).embM[S][Y][c].SR[SO10]^(-1):
#%
#SubGrd(D5,[1,0,0,0,0], [B3,A1],embM[D][SO7SU2_SO10][n]);
#%
end if:
#%
#>Slansky embedding
WTDs:=WTD[D5]([1,1,2,1,1]).WTD[D5]([0,1,2,1,1]).WTD[D5]([0,0,1,0,1]).WTD[D5]([0,0,1,0,0]):
Y1:=SO7SU2_SO10: Z:=s:
embM[D][Y1][Z]:=embM[D][Y][c].WTDs:
embM[S][Y1][Z]:=DiagonalMatrix([CM[SO7],CM[SL2]])^(-1).embM[D][Y1][Z].CM[SO10]:
embM[H][Y1][Z]:=DiagonalMatrix([SR[SO7],SR[SL2]]).embM[S][Y1][Z].SR[SO10]^(-1):
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
#% End: Job_SO10_SubG_SO7SU2
#%
##>Sp2xSp2-> SO10
#% Begin: Job_SO10_SubG_Sp2Sp2
if Job_SO10_SubG_Sp2Sp2 then
#print("Job: Embedding of Sp2xSp2 to SO10");
#%
embM[H][Sp2SO10][1]:=Matrix([[1,1,0,0,0],[1,-1,0,0,0]]):
embM[H][Sp2SO10][2]:=Matrix([[0,0,1,1,0],[0,0,1,-1,0]]):
i:='i': for i from 1 to 2 do
  embM[S][Sp2SO10][i]:=SR[Sp2]^(-1).embM[H][Sp2SO10][i].SR[SO10]:
  embM[D][Sp2SO10][i]:=CM[Sp2].embM[S][Sp2SO10][i].CM[SO10]^(-1):
end do:
#%
Y:=Sp2Sp2_SO10: 
X:='X': for X in [H,S,D] do
 embM[X][Y][n]:=<<embM[X][Sp2SO10][1],embM[X][Sp2SO10][2]>>:
end do:
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
#SubGrd(D5,[1,0,0,0,0], [C2,C2],embM[D][Sp2Sp2_SO10][n]);
#%
end if:
#% End: Job_SO10_SubG_Sp2Sp2
#%
##> max Sp2-> SO10
#% Begin: Job_SO10_SubG_Sp2
if Job_SO10_SubG_Sp2 then
#print("Job: Embedding of max Sp2 to SO10");
#%
embM[H][Sp2SO10][3]:=Matrix([[2,1,1,0,0],[0,1,-1,2,0]]):
embM[S][Sp2SO10][3]:=SR[Sp2]^(-1).embM[H][Sp2SO10][3].SR[SO10]:
embM[D][Sp2SO10][3]:=CM[Sp2].embM[S][Sp2SO10][3].CM[SO10]^(-1):
#%
Y:=Sp2_SO10: 
X:='X': for X in [H,S,D] do
 embM[X][Y][n]:=embM[X][Sp2SO10][3]:
end do:
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
#SubGrd(D5,[1,0,0,0,0], [C2],embM[D][Sp2_SO10][n]);
#%
end if:
#% End: Job_SO10_SubG_Sp2
#%
end if:
#% End: Job_SO10 
#%
#####>Type E
####>E6
#% Begin: Job_E6 
if Job_E6  then
printf("E6, ");
#%
###>WeylTrf
rv:='rv':  WTD[E6]:= rv -> WeylTrf(E6,rv,"D"):
###>Reps
##> 27dim repr
WS[E6_27]:=WeightSystem(E6,[1,0,0,0,0,0],"S"):
hl:=WS[E6_27]["hl"]:
i:='i': for i from 1 to 6 do
  HV[E6_27][i]:=[]:
  j:='j': for j from hl to -hl by -1 do
     x:='x': for x in WS[E6_27][j] do
        k:='k':
        HV[E6_27][i]:=[op(HV[E6_27][i]),seq(op(2,x)[i],k=1..op(1,x))]:
     end do:
  end do:
end do:
i:='i': for i from 1 to 6 do
  j:='j': for j in [7,9,11,16,18,20] do
    tmpV:=HV[E6_27][i][j]:
    HV[E6_27][i][j]:=HV[E6_27][i][j+1]:
    HV[E6_27][i][j+1]:=tmpV:
  end do:
end do:
#%
###>Subalgebras
##>SO10xU1
#% Begin: Job_E6_SubG_SO10U1 
if Job_E6_SubG_SO10U1 then
#print("Subalgebra so10+u1");
#% SO10 -> E6
#% alpha_i=alph'_i (i=1,2,3,4,6), 
KMI[D5]:=Kmetric(D5)^(-1): KM[E6]:=Kmetric(E6):
#>normal: alpha_5 => U1
#% alpha_5=sum_{i,a<>5}(Kmetric(D5)^(-1))[i,a]*Kmetric(E6)[a,5] alpha'_i
a5v:=Vector(5):
i:='i': for i from 1 to 5 do 
   j:='j': a5v[i]:=add(KMI[D5][i,j]*KM[E6][j,5],j=1..4)
                    + KMI[D5][i,5]*KM[E6][6,5]: 
end do:
embM[S][SO10E6][1]:=<<1,0,0,0,0>|<0,1,0,0,0>|<0,0,1,0,0>|
               <0,0,0,1,0>|a5v|<0,0,0,0,1>>:
#               <0,0,0,1,0>|<-1/2,-1,-3/2,-5/4,-3/4>|<0,0,0,0,1>>:
embM[H][SO10E6][1]:=SR[SO10].embM[S][SO10E6][1].KM[E6]^(-1):
embM[D][SO10E6][1]:=CM[SO10].embM[S][SO10E6][1].CM[E6]^(-1):
#% U1 -> E6
#% (beta,alpha_i)=0 (i=1,2,3,4,6)
b:='b':
beta:=[b[1],b[2],b[3],b[4],-1,b[5]]:
i:='i': j:='j': eqs:=[seq(add(beta[i]*KM[E6][i,j],i=1..6),j in [1,2,3,4,6])]:
i:='i': sol:=solve(eqs,[seq(b[i],i=1..5)]):
Qvec[H][E6][1]:=Matrix([eval(beta,sol[1])]):
Qvec[S][E6][1]:=Qvec[H][E6][1].KM[E6]:
Qvec[D][E6][1]:=Qvec[S][E6][1].CM[E6]^(-1):
#%
Y:=SO10U1_E6: Z:=n1:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SO10E6][1],Qvec[X][E6][1]>>:
end do:
#%
#>normal: alpha_1 => U1
#% alpha_1=sum_{i,a<>1}(Kmetric(D5)^(-1))[i,a]*Kmetric(E6)[a,1] alpha'_i
a5v:=Vector(5):
i:='i': for i from 1 to 5 do 
   j:='j': a5v[i]:=add(KMI[D5][i,j]*KM[E6][6-j,1],j=1..4)
                   + KMI[D5][i,5]*KM[E6][6,1]:
end do:
embM[S][SO10E6][2]:=<a5v|<0,0,0,1,0>|<0,0,1,0,0>|<0,1,0,0,0>|
               <1,0,0,0,0>|<0,0,0,0,1>>:
# a5v=<-1/2, -1, -3/2, -5/4, -3/4>
embM[H][SO10E6][2]:=SR[SO10].embM[S][SO10E6][2].KM[E6]^(-1):
embM[D][SO10E6][2]:=CM[SO10].embM[S][SO10E6][2].CM[E6]^(-1):
#% U1 -> E6
#% (beta,alpha_i)=0 (i=1,2,3,4,6)
b:='b':
beta:=[1,b[4],b[3],b[2],b[1],b[5]]:
i:='i': j:='j': eqs:=[seq(add(beta[i]*KM[E6][i,j],i=1..6),j=2..6)]:
i:='i': sol:=solve(eqs,[seq(b[i],i=1..5)]):
Qvec[H][E6][2]:=Matrix([eval(beta,sol[1])]):
Qvec[S][E6][2]:=Qvec[H][E6][2].KM[E6]:
Qvec[D][E6][2]:=Qvec[S][E6][2].CM[E6]^(-1):
#%
Y:=SO10U1_E6: Z:=n2:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SO10E6][2],Qvec[X][E6][2]>>:
end do:
#%
#>canonical (Slansky)
embM[D][SO10E6][3]:=<<0, 0, 0, 0, 1>|<1, 0, 0, 0, 1>|<1, 0, 1, 0, 0>
                  |<1, 0, 0, 1, 0>|<0, 0, 0, 1, 0>|<0, 1, 0, 0, 0>>:
embM[S][SO10E6][3]:=CM[SO10]^(-1).embM[D][SO10E6][3].CM[E6]:
embM[H][SO10E6][3]:=SR[SO10].embM[S][SO10E6][3].KM[E6]^(-1):
#%
b:='b': i:='i': beta:=Matrix([[seq(b[i],i=1..6)]]):
eqs:=[]: i:='i': for i from 1 to 5 do
  eqs:=[op(eqs),DotProduct(beta.Gmetric(E6),embM[D][SO10E6][3][i,1..6] )]:
end do:
#print("eqs"=eqs);
sol:=solve(eqs,convert(beta,list)):
#print("sol"=sol);
Qvec[D][E6][3]:=Matrix(eval([eval(beta,sol[1])],b[5]=-1)):
Qvec[S][E6][3]:=Qvec[D][E6][3].CM[E6]:
Qvec[H][E6][3]:=Qvec[S][E6][3].KM[E6]^(-1):
#%
Y:=SO10U1_E6: Z:=c:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SO10E6][3],Qvec[X][E6][3]>>:
end do:
###% Weyl transformation: s1 -> c2=c
#%
if false then
WTDs:=WTD[E6]([0,0,0,0,1,0]).WTD[E6]([0,0,1,1,0,1]).WTD[E6]([1,1,1,1,0,0]).WTD[E6]([0,1,1,1,0,0]):
Z:='c2':
embM[D][Y][Z]:=embM[D][Y][n1].WTDs:
embM[S][Y][Z]:=DiagonalMatrix([CM[SO10],CM[SL2]])^(-1).embM[D][Y][Z].CM[E6]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SO10],SR[SL2]]).embM[S][Y][Z].KM[E6]^(-1):
end if:
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n1]:
end do:
#%               
# SubGrd(E6,[1,0,0,0,0,0],[D5,U1],embM[D][SO10U1_E6][c]):
#%
#% 
end if:
#% End: Job_E6_SubG_SO10U1
#%
##>SU6xSU2
#% Begin: Job_E6_SubG_SU6SU2 
if Job_E6_SubG_SU6SU2 then
#print("The subalgebra su6+su2");
#%
#>normal embedding: alpha_4 => -theta
embM[S][SU6E6]:=<<1,0,0,0,0>|<0,1,0,0,0>|<0,0,1,0,0>|
  <-1/2,-1,-3/2,-1,-1/2>|<0,0,0,0,0>|<0,0,0,1,0>>:
embM[H][SU6E6]:=SR[SL6].embM[S][SU6E6].KM[E6]^(-1):
embM[D][SU6E6]:=CM[SL6].embM[S][SU6E6].CM[E6]^(-1):
#%
embM[S][SU2E6]:=Matrix([[0,0,0,-1/2,1,0]]):
embM[H][SU2E6]:=SR[SL2].embM[S][SU2E6].KM[E6]^(-1):
embM[D][SU2E6]:=CM[SL2].embM[S][SU2E6].CM[E6]^(-1):
#%
Y:=SU6SU2_E6: Z:=n1:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][SU6E6],embM[X][SU2E6]>>:
end do:
#%
#>normal embedding: alpha6 => -theta
Y:=SU6SU2_E6: Z:=n2:
embM[S][Y][Z]:=<<1,0,0,0,0,0>|<0,1,0,0,0,0>|<0,0,1,0,0,0>|
  <0,0,0,1,0,0>|<0,0,0,0,1,0>|<-1/2,-1,-3/2,-1,-1/2,-1/2>>:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL6],SR[SL2]]).embM[S][Y][Z].KM[E6]^(-1):
embM[D][Y][Z]:=DiagonalMatrix([CM[SL6],CM[SL2]]).embM[S][Y][Z].CM[E6]^(-1):
#%
#>canonical embedding:
WTDs:=WTD[E6]([0,1,2,1,0,1]).WTD[E6]([0,0,0,0,0,1]).WTD[E6]([0,0,1,1,1,0]).WTD[E6]([0,1,1,1,0,0]).WTD[E6]([1,1,1,0,0,0]):
Z:=c:
embM[D][Y][Z]:=embM[D][Y][n2].WTDs:
embM[S][Y][Z]:=DiagonalMatrix([CM[SL6],CM[SL2]])^(-1).embM[D][Y][Z].CM[E6]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL6],SR[SL2]]).embM[S][Y][Z].KM[E6]^(-1):
#%
#>Slansky embedding (1):
WTDs:=WTD[E6]([1,1,0,0,0,0]).WTD[E6]([0,1,2,2,1,1]).WTD[E6]([0,1,2,1,0,1]).WTD[E6]([0,1,1,0,0,0]).WTD[E6]([0,0,1,0,0,0]):
Y1:=SU6SU2_E6: Z:=s1:
embM[D][Y1][Z]:=embM[D][Y][c].WTDs:
embM[S][Y1][Z]:=DiagonalMatrix([CM[SL6],CM[SL2]])^(-1).embM[D][Y1][Z].CM[E6]:
embM[H][Y1][Z]:=DiagonalMatrix([SR[SL6],SR[SL2]]).embM[S][Y1][Z].KM[E6]^(-1):
#%
#>Slansky embedding (2):
#% complex conjugate of (1)
WTDs:=WTD[E6]([1,2,2,1,1,1]).WTD[E6]([0,1,1,1,1,0]).WTD[E6]([0,0,1,0,0,1]).WTD[E6]([0,1,1,1,0,0]).WTD[E6]([0,0,1,1,0,0]):
Z:=s2:
embM[D][Y1][Z]:=<<embM[D][Y1][s1][5,1..6],embM[D][Y1][s1][4,1..6],
  embM[D][Y1][s1][3,1..6],embM[D][Y1][s1][2,1..6],embM[D][Y1][s1][1,1..6],
  embM[D][Y1][s1][6,1..6]>>.WTDs:
embM[S][Y1][Z]:=DiagonalMatrix([CM[SL6],CM[SL2]])^(-1).embM[D][Y1][Z].CM[E6]:
embM[H][Y1][Z]:=DiagonalMatrix([SR[SL6],SR[SL2]]).embM[S][Y1][Z].KM[E6]^(-1):
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n1]:
end do:
#%
# SubGrd(E6,[1,0,0,0,0,0],[A5,A1],embM[D][SU6SU2_E6][c]):
#%
end if:
#% End: Job_E6_SubG_SU6SU2
#%
#%
##>SU3xSU3xSU3
#% Begin: Job_E6_SubG_SU3SU3SU3 
if Job_E6_SubG_SU3SU3SU3 then
#print("The subalgebra su3+su3+su3");
#>normal embedding
Y:=SU3SU3SU3_E6: Z:=n:
embM[S][Y][Z]:=<<1,0,0,0,0,0>|<0,1,0,0,0,0>|<-1/3,-2/3,-2/3,-1/3,-2/3,-1/3>|<0,0,1,0,0,0>|<0,0,0,1,0,0>|<0,0,0,0,1,0>>:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL3],SR[SL3],SR[SL3]]).embM[S][Y][Z].SR[E6]^(-1):
embM[D][Y][Z]:=DiagonalMatrix([CM[SL3],CM[SL3],CM[SL3]]).embM[S][Y][Z].CM[E6]^(-1):
#%
#>canonical embedding (1)
WTDs:=WTD[E6]([0,1,2,1,1,1]).WTD[E6]([0,0,1,1,0,1]).WTD[E6]([1,1,1,0,0,0]).WTD[E6]([0,1,1,1,0,0]).WTD[E6]([0,0,0,1,0,0]):
Y:=SU3SU3SU3_E6: Z:=c:
embM[D][Y][Z]:=embM[D][Y][n].WTDs:
embM[S][Y][Z]:=DiagonalMatrix([CM[SL3],CM[SL3],CM[SL3]])^(-1).embM[D][Y][Z].CM[E6]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL3],SR[SL3],SR[SL3]]).embM[S][Y][Z].KM[E6]^(-1):
#%
#>canonical embedding (2) (Slansky)
Y:=SU3SU3SU3_E6: Z:=s:
embM[D][Y][Z]:=<<1,0,0,0,1,0>|<1,-1,0,0,2,0>|<1,-1,1,-1,2,1>|<1,-1,0,-1,1,1>|<1,-1,0,0,0,1>|<0,0,0,0,1,1>>:
embM[S][Y][Z]:=DiagonalMatrix([CM[SL3],CM[SL3],CM[SL3]])^(-1).embM[D][Y][Z].CM[E6]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL3],SR[SL3],SR[SL3]]).embM[S][Y][Z].KM[E6]^(-1):
#%
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
##%
# i:='i':
# SubGrd(E6,[1,0,0,0,0,0],[A2,A2,A2],embM[D][SU3SU3SU3_E6][c]]):
#%
end if:
#% End: Job_E6_SubG_SU3SU3SU3
#%
#%
##>special SU3
#% Begin: Job_E6_SubG_SU3 
if Job_E6_SubG_SU3 then
#print("The subalgebra special su3");
#%
#% 27dim reps of SU3
WS[SU3_27]:=WeightSystem(A2,[2,2],"S"):
hl:=WS[SU3_27]["hl"]:
i:='i': for i from 1 to 2 do
  HV[SU3_27][i]:=[]:
  j:='j': for j from hl to -hl by -1 do
     x:='x': for x in WS[SU3_27][j] do
        k:='k':
        HV[SU3_27][i]:=[op(HV[SU3_27][i]),seq(op(2,x)[i],k=1..op(1,x))]:
     end do:
  end do:
end do:
#%
embM[H][SU3_E6][c]:=Matrix([[1,-1,2,-1,1,-1],[0,1,-1,1,0,1]]):
embM[S][SU3_E6][c]:=embM[H][SU3_E6][c]:
embM[D][SU3_E6][c]:=CM[SL3].embM[S][SU3_E6][c].CM[E6]^(-1):
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][c]:
end do:
#%
# i:='i':
# SubGrd(E6,[1,0,0,0,0,0],[A2],embM[D][SU3_E6][c]):
#%
end if:
#% End: Job_E6_SubG_SU3
#%
##>special G2
#% Begin: Job_E6_SubG_G2 
if Job_E6_SubG_G2 then
#print("The subalgebra special g2");
#%
#>normal embedding
WS[G2_27]:=WeightSystem(G2,[0,2],"S"):
hl:=WS[G2_27]["hl"]:
i:='i': for i from 1 to 2 do
  HV[G2_27][i]:=[]:
  j:='j': for j from hl to -hl by -1 do
     x:='x': for x in WS[G2_27][j] do
        k:='k':
        HV[G2_27][i]:=[op(HV[G2_27][i]),seq(op(2,x)[i],k=1..op(1,x))]:
     end do:
  end do:
end do:
#%
Y:=G2_E6: Z:=n:
embM[H][Y][Z]:=Matrix([[0,1,-1,1,0,1],[1,0,1,0,1,0]]):
embM[S][Y][Z]:=embM[H][Y][Z]:
embM[D][Y][Z]:=CM[G2].embM[S][Y][Z].CM[E6]^(-1):
#%
#>canonical embedding: non existent
if false then
Z:=c:
embM[D][Y][Z]:=<<0,2>|<1,2>|<1,3>|<1,2>|<0,2>|<1,1>>:
embM[S][Y][Z]:=CM[G2]^(-1).embM[D][Y][Z].CM[E6]:
embM[H][Y][Z]:=embM[S][Y][Z]:
end if:
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# i:='i':
# SubGrd(E6,[1,0,0,0,0,0],[G2],embM[D][G2_E6][c]):
#%
end if:
#% End: Job_E6_SubG_G2
#%
#%
##>special G2xSU3
#% Begin: Job_E6_SubG_G2SU3 
if Job_E6_SubG_G2SU3 then
#print("The subalgebra special g2+su3");
#%  
WS[SU3_3]:=WeightSystem(A2,[1,0],"S"):
WS[SU3_6b]:=WeightSystem(A2,[0,2],"S"):
WS[G2_7]:=WeightSystem(G2,[0,1],"S"):
Rep:='Rep': for Rep in [SU3_3, SU3_6b,G2_7] do
  hl:=WS[Rep]["hl"]:
  rank:=2:
  i:='i': for i from 1 to rank do
    HV[Rep][i]:=[]:
    j:='j': for j from hl to -hl by -1 do
       x:='x': for x in WS[Rep][j] do
         k:='k':
         HV[Rep][i]:=[op(HV[Rep][i]),seq(op(2,x)[i],k=1..op(1,x))]:
       end do:
    end do:
  end do:
end do:
#%
#% find an embedding (H'_1,H'_2)|SU3 + (H'_3,H'_4)|G2 -> (H_1, ..., H_6)|E6
if true then
c := 'c': i := 'i': sol:='sol':
HV6 := convert(add(c[i]*Vector(HV[E6_27][i]), i = 1 .. 6), list):
i:='i': for i from 1 to 2 do 
  x := 'x':
  eqs := [seq(HV6[x[1]] = HV[SU3_6b][i][x[2]], x in [[9,1],[11,2],[15,3],[13,4],[17,5],[19,6]])]:
  j := 'j':
  sol[i] := solve(eqs, [seq(c[j], j = 1 .. 6)])[1]:
  res[i]:= eval(HV6,sol[i]);
end do;
end if:
#%
embM[H][SU3E6][5]:=Matrix([[1,0,-1,1,0,0],[0,1,-1,0,1,0]]):
embM[S][SU3E6][5]:=embM[H][SU3E6][5]:
embM[D][SU3E6][5]:=CM[SL3].embM[S][SU3E6][5].CM[E6]^(-1):
#%
embM[H][G2E6][2]:=Matrix([[0,0,0,0,0,1],[0,0,1,0,0,0]]):
embM[S][G2E6][2]:=embM[H][G2E6][2]:
embM[D][G2E6][2]:=CM[G2].embM[S][G2E6][2].CM[E6]^(-1):
#%
Y:=G2SU3_E6: Z:=c:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=<<embM[X][G2E6][2],embM[X][SU3E6][5]>>:
end do:
#%
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][c]:
end do:
#%
# i:='i':
# SubGrd(E6,[1,0,0,0,0,0],[G2,A2],embM[D][G2SU3_E6][c]):
#%
end if:
#% End: Job_E6_SubG_G2SU3
#%
##>special Sp4
#% Begin: Job_E6_SubG_Sp4 
if Job_E6_SubG_Sp4 then
#print("The subalgebra special sp4");
#%
#>normal embedding
WS[Sp4_27]:=WeightSystem(C4,[0,1,0,0],"S"):
hl:=WS[Sp4_27]["hl"]:
i:='i': for i from 1 to 4 do
  HV[Sp4_27][i]:=[]:
  j:='j': for j from hl to -hl by -1 do
     x:='x': for x in WS[Sp4_27][j] do
        k:='k':
        HV[Sp4_27][i]:=[op(HV[Sp4_27][i]),seq(op(2,x)[i],k=1..op(1,x))]:
     end do:
  end do:
end do:
#% swap the comp [7,8]
X:='Sp4_27':
i:='i': for i from 1 to 4 do
  j:='j': for j in [7,20] do
    tmpV:=HV[X][i][j]:
    HV[X][i][j]:=HV[X][i][j+1]:
    HV[X][i][j+1]:=tmpV:
  end do:
end do:
#% find an embedding (H_1,..,H_4)|Sp4 -> (H_1, ..., H_6)|E6
if false then
c := 'c': i := 'i': sol:='sol':
HV6 := convert(add(c[i]*Vector(HV[E6_27][i]), i = 1 .. 6), list):
i:='i': for i from 1 to 4 do 
  j := 'j':
  eqs := [seq(HV6[j] = HV[Sp4_27][i][j], j in [1, 2, 3, 4, 5, 6])]:
  j := 'j':
  sol[i] := op(1,solve(eqs, [seq(c[j], j = 1 .. 6)])):
  res[i]:=eval(HV6,sol[i])-HV[Sp4_27][i];
end do;
end if:
#%
Y:=Sp4_E6: Z:=n:
embM[H][Y][Z]:=Matrix([[0,1,-1,1,0,0],[1,0,0,0,1,0],[0,0,1,0,0,0],[0,0,0,0,0,1]]):
embM[S][Y][Z]:=embM[H][Y][Z]:
embM[D][Y][Z]:=CM[Sp4].embM[S][Y][Z].CM[E6]^(-1):
#%
#>canonical embedding: non existent
Z:=c:
embM[D][Y][Z]:=<<0,1,0,0>|<1,0,1,0>|<2,0,0,1>|<1,0,1,0>|<0,1,0,0>|<0,0,0,1>>: #=embM[D][Y][n].WT([0,0,1,0,0,0])
embM[S][Y][Z]:=CM[Sp4]^(-1).embM[D][Y][Z].CM[E6]:
embM[H][Y][Z]:=embM[S][Y][Z]:
#%
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# i:='i':
# SubGrd(E6,[1,0,0,0,0,0],[C4],embM[D][Sp4_E6][c]):
#%
end if:
#% End: Job_E6_SubG_Sp4
#%
##>special F4
#% Begin: Job_E6_SubG_F4 
if Job_E6_SubG_F4 then
#print("The subalgebra special f4");
#%
#>Canonical embedding
#% find an embedding (H_1,..,H_4)|F4 -> (H_1, ..., H_6)|E6
#% swap components of HV[F4_26]
if true then
Rep:='F4_26': rank:=4:
i:='i': for i from 1 to rank do
  j:='j': for j in [5,21] do
    tmpV:=HV[Rep][i][j]:
    HV[Rep][i][j]:=HV[Rep][i][j+1]:
    HV[Rep][i][j+1]:=tmpV:
  end do:
end do:
#%
c := 'c': i := 'i': sol:='sol':
X:='E6_27':
HV4s[1]:=HV[X][1]+HV[X][5]:
HV4s[2]:=HV[X][2]+HV[X][4]:
HV4s[3]:=HV[X][3];
HV4s[4]:=HV[X][6]:
HV4 := convert(add(c[i]*Vector(HV4s[i]),i=1..4), list):
HV4:=[op(HV4[1..13]),op(HV4[15..27])]:
i:='i': for i from 1 to 4 do 
  j := 'j':
  eqs := [seq(HV4[j] = HV[F4_26][i][j], j in [1, 2, 3, 4])]:
  j := 'j':
  sol[i] := op(1,solve(eqs, [seq(c[j], j = 1 .. 4)])):
  res[i]:=eval(HV4,sol[i])-HV[F4_26][i];
end do;
end if:
#% embedding matrix
Y:=F4_E6: Z:=c:
embM[H][Y][Z]:=Matrix([[0,0,0,0,0,1],[0,0,1,0,0,0],[0,1,0,1,0,0],[1,0,0,0,1,0]]):
embM[S][Y][Z]:=embM[H][Y][Z]:
embM[D][Y][Z]:=CM[F4].embM[S][Y][Z].CM[E6]^(-1):
#%
#>Slansky embedding
WTDs:=WTD[E6]([0,0,1,1,1,0]).WTD[E6]([0,1,1,1,0,0]).WTD[E6]([0,0,1,0,0,0]):
Y:=F4_E6: Z:=s:
embM[D][Y][Z]:=embM[D][Y][c].WTDs:
embM[S][Y][Z]:=CM[F4]^(-1).embM[D][Y][Z].CM[E6]:
embM[H][Y][Z]:=embM[S][Y][Z]:
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][c]:
end do:
#%
# SubGrd(E6,[1,0,0,0,0,0],[F4],embM[D][F4_E6][s]);
#%
end if:
#% End: Job_E6_SubG_F4
end if:
#% End: Job_E6
#####>Type F
####>F4
#% Begin: Job_F4 
if Job_F4  then
printf("F4,");
#%
###>Weyl Trf
rv:='rv': WTD[F4]:= rv -> WeylTrf(F4,rv,"D"):
###>Reps
##> 26dim repr
WS[F4_26]:=WeightSystem(F4,[0,0,0,1],"S"):
Rep:=F4_26:
hl:=WS[Rep]["hl"]:
rank:=4:
i:='i': for i from 1 to rank do
  HV[Rep][i]:=[]:
  j:='j': for j from hl to -hl by -1 do
     x:='x': for x in WS[Rep][j] do
        k:='k':
        HV[Rep][i]:=[op(HV[Rep][i]),seq(op(2,x)[i],k=1..op(1,x))]:
     end do:
  end do:
end do:
#%
###>Subalgebras
##>SO9
#% Begin: Job_F4_SubG_SO9
if Job_F4_SubG_SO9 then
#print("Subalgebra: SO9");
#%
#>normal embedding
Y:=SO9_F4: Z:=n:
embM[S][Y][Z]:=<<0,1,0,0>|<0,0,1,0>|<0,0,0,1>|<-1/2,-1,-3/2,-2>>:
embM[D][Y][Z]:=CM[SO9].embM[S][Y][Z].CM[F4]^(-1):
embM[H][Y][Z]:=SR[SO9].embM[S][Y][Z]:
#%
#>Slansky embedding
WTDs:=WTD[F4]([1,2,4,2]).WTD[F4]([0,1,2,1]).WTD[F4]([0,1,2,0]).WTD[F4]([0,1,0,0]):
Y:=SO9_F4: Z:=s:
embM[D][Y][Z]:=embM[D][Y][n].WTDs:
embM[S][Y][Z]:=CM[SO9]^(-1).embM[D][Y][Z].CM[F4]:
embM[H][Y][Z]:=SR[SO9].embM[S][Y][Z];
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrdm(F4,[1,0,0,0],[B4],embM[D][SO9_F4][s]);
#%
end if:
#% End: Job_F4_SubG_SO9
#%
##>SU3xSU3
#% Begin: Job_F4_SubG_SU3SU3
if Job_F4_SubG_SU3SU3 then
#print("Subalgebra: SU3xSU3");
#%
#>normal embedding
Y:=SU3SU3_F4: Z:=n:
embM[S][Y][Z]:=<<0,0,0,1>|<-4/3,-2/3,-1/3,-2/3>|<1,0,0,0>|<0,1,0,0>>:
embM[D][Y][Z]:=DiagonalMatrix([CM[SL3],CM[SL3]]).embM[S][Y][Z].CM[F4]^(-1):
embM[H][Y][Z]:=DiagonalMatrix([SR[SL3],SR[SL3]]).embM[S][Y][Z]:
#>Slansky embeddging
WTDs:=WTD[F4]([2,3,4,2]).WTD[F4]([0,1,1,1]).WTD[F4]([1,2,2,0]).WTD[F4]([0,0,1,0]).WTD[F4]([1,1,0,0]).WTD[F4]([0,1,0,0]):
Y:=SU3SU3_F4: Z:=s:
embM[D][Y][Z]:=embM[D][Y][n].WTDs:
embM[S][Y][Z]:=DiagonalMatrix([CM[SL3],CM[SL3]])^(-1).embM[D][Y][Z].CM[F4]:
embM[H][Y][Z]:=DiagonalMatrix([SR[SL3],SR[SL3]]).embM[S][Y][Z];
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
#%
# SubGrdm(F4,[1,0,0,0],[A2,A2],embM[D][SU3SU3_F4][s]);
#%
end if:
#% End: Job_F4_SubG_SU3SU3
#%
##>Sp3xSU2
#% Begin: Job_F4_SubG_Sp3SU2
if Job_F4_SubG_Sp3SU2 then
#print("Subalgebra: Sp3xSU2");
#%
#>normal embedding
Y:=Sp3SU2_F4: Z:=n:
embM[S][Y][Z]:=<<-1,-2,-3/2,-1/2>|<0,0,1,0>|<0,1,0,0>|<1,0,0,0>>:
embM[D][Y][Z]:=DiagonalMatrix([CM[Sp3],CM[SL2]]).embM[S][Y][Z].CM[F4]^(-1):
embM[H][Y][Z]:=DiagonalMatrix([SR[Sp3],SR[SL2]]).embM[S][Y][Z]:
#%
#>canonical embedding
WTDs:=WTD[F4]([1,1,1,1]).WTD[F4]([1,1,2,0]).WTD[F4]([1,1,0,0]).WTD[F4]([0,1,0,0]):
Z:=c:
embM[D][Y][Z]:=embM[D][Y][n].WTDs:
embM[S][Y][Z]:=DiagonalMatrix([CM[Sp3],CM[SL2]])^(-1).embM[D][Y][Z].CM[F4]:
embM[H][Y][Z]:=DiagonalMatrix([SR[Sp3],SR[SL2]]).embM[S][Y][Z];
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrdm(F4,[1,0,0,0],[A2,A2],embM[D][Sp3SU2_F4][c]);
#%
end if:
#% End: Job_F4_SubG_Sp3SU2
#%
##>G2xSU2 (special)
#% Begin: Job_F4_SubG_G2SU2
if Job_F4_SubG_G2SU2 then
#print("Subalgebra: G2xSU2 (special)");
#%
#>normal embedding
Y:=G2SU2_F4: Z:=n:
embM[S][Y][Z]:=<<1,0,0>|<0,1,0>|<0,0,-1>|<0,0,2>>:
embM[D][Y][Z]:=DiagonalMatrix([CM[G2],CM[SL2]]).embM[S][Y][Z].CM[F4]^(-1):
embM[H][Y][Z]:=DiagonalMatrix([SR[G2],SR[SL2]]).embM[S][Y][Z]:
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
#%
# SubGrdm(F4,[1,0,0,0],[G2,A1],embM[D][G2SU2_F4][n]);
#%
end if:
#% End: Job_F4_SubG_G2SU2
#%
end if:
#% End: Job_F4
#####>Type G
####>G2
#% Begin: Job_G2
if Job_G2 then
printf("G2\n");
#%
###>Root Diagram
#% Root vectors
alpha:='alpha':
alpha[1]:=Vector([1,0]):
alpha[2]:=Vector([-1/2,sqrt(3)/6]):
Deltap:=[alpha[1],alpha[2],alpha[1]+alpha[2],alpha[1]+2*alpha[2],
  alpha[1]+3*alpha[2],2*alpha[1]+3*alpha[2]]:
Deltam:=-1*Deltap:
#% drawing
RootsPlot:=[]:
v:='v': for v in Deltap do 
  RootsPlot:=[op(RootsPlot),arrow(v,color=green)]:
end do:
v:='v': for v in Deltam do 
  RootsPlot:=[op(RootsPlot),arrow(v,color=blue)]:
end do:
circle:=plot([cos(t),sin(t),t=0..2*Pi],thickness=1):
#print(display([op(RootsPlot),circle]));
#%
###>Subalgebra
##>SU3
#% Begin: Job_G2_SubG_SU3
if Job_G2_SubG_SU3 then
#print("Subalgebra SU3");
#%
#>Normal embedding
Y:=SU3_G2: Z:=n:
embM[S][Y][Z]:=<<1,0>|<-2/3,-1/3>>:
embM[H][Y][Z]:=SR[SL3].embM[S][Y][Z].SR[G2]^(-1):
embM[D][Y][Z]:=CM[SL3].embM[S][Y][Z].CM[G2]^(-1):
#>default embM
Z:=0:
X:='X': for X in [H,S,D] do
  embM[X][Y][Z]:=embM[X][Y][n]:
end do:
#%
# SubGrdm(G2,[1,0],[A2],embM[D][SU3_G2][n]);
#%
end if
#% End: Job_G2_SubG_SU3
end if:
#% End: Job_G2
#%


#####>Embedding matrix list
####>proc:embMlist
#% embMlist(dt) dt:Dynkin type (e.g. A5) => list of the embM[D] indeces
embMlist:=proc()
local dt,dt0,dc,rn:

dt0:=resolvedt(_passed[1]):
dc:=dt0[1]; rn:=dt0[2]:
dt:=abbrevdt(dt0):
if dc="A" then
   if rn=1 then  
      printf(" %-15s: [%a][0]\n",  U1, U1_SU2);
   elif rn=2 then
      printf(" %-15s: [%a][0]\n",  U1xU1, U1U1_SU3);
      printf(" %-15s: [%a][%a]\n", SO3, SO3_SU3, n):
      printf(" %-15s: [%a][%a]\n", 'SU2xU1', SU2U1_SU3, 'n/c'):
   elif rn=3 then
      printf(" %-15s: [%a][%a]\n", 'SU3xU1', SU3U1_SU4, 'c/s'):
      printf(" %-15s: [%a][%a]\n", 'SU2xSU2xU1', SU2SU2U1_SU4, 'n/c'):
      printf(" %-15s: [%a][%a]\n", 'SU2xSU2', SU2SU2_SU4, 'c'):
      printf(" %-15s: [%a][%a]\n", 'Sp2', Sp2_SU4, 'n'):
   elif rn=4 then
      printf(" %-15s: [%a][%a]\n", 'SU4xU1', SU4U1_SU5, 'n/c/s'):
      printf(" %-15s: [%a][%a]\n", 'SU3xSU2xU1', SU3SU2U1_SU5, 'n/c/s'):
      printf(" %-15s: [%a][%a]\n", 'SO5', SO5_SU5, 'n'):
   elif rn=5 then
      printf(" %-15s: [%a][%a]\n", 'SU5xU1', SU5U1_SU6, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SU4xSU2xU1', SU4SU2U1_SU6, 'n/s'):
       printf(" %-15s: [%a][%a]\n", 'SU3xSU3xU1', SU3SU3U1_SU6, 'n/s'):
      printf(" %-15s: [%a][%a]\n", 'SU4', SU4_SU6, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SU3', SU3_SU6, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SU3xSU2', SU3SU2_SU6, 'n'):
      printf(" %-15s: [%a][%a]\n", 'Sp3', Sp3_SU6, 'n'):
   end if:
elif dc="B" then
   if rn=2 then
      printf(" %-15s: [%a][%a]\n", 'SU2xU1', SU2U1_SO5, 'n/c'):
      printf(" %-15s: [%a][%a]\n", 'SU2xSU2', SU2SU2_SO5, 'n1/n2'):
   elif rn=3 then
      printf(" %-15s: [%a][%a]\n", 'SU4', SU4_SO7, 'n/s'):
      printf(" %-15s: [%a][%a]\n", 'SO5xU1', SO5U1_SO7, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SU2xSU2xSU2', SU2SU2SU2_SO7, 'n'):
      printf(" %-15s: [%a][%a]\n", 'G2', G2_SO7, 'n'):
   elif rn=4 then
      printf(" %-15s: [%a][%a]\n", 'SO8', SO8_SO9, 'n/s'):
      printf(" %-15s: [%a][%a]\n", 'SO7xU1', SO7U1_SO9, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SU4xSU2', SU4SU2_SO9, 'n/s'):
      printf(" %-15s: [%a][%a]\n", 'SU2xSU2xSp2', SU2SU2Sp2_SO9, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SU2', SU2_SO9, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SU2xSU2', SU2SU2_SO9, 'n'):
   end if:
elif dc="C" then
   if rn=3 then
      printf(" %-15s: [%a][%a]\n", 'SU3xU1', SU3U1_Sp3, 'n'):
      printf(" %-15s: [%a][%a]\n", 'Sp2xSU2', Sp2SU2_Sp3, 'n2/c'):
      printf(" %-15s: [%a][%a]\n", 'SU2', SU2_Sp3, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SU2xSU2', SU2SU2_Sp3, 'n'):
   elif rn=4 then
      printf(" %-15s: [%a][%a]\n", 'SU4xU1', SU4U1_Sp4, 'n'):
      printf(" %-15s: [%a][%a]\n", 'Sp3xSU2', Sp3SU2_Sp4, 'n'):
      printf(" %-15s: [%a][%a]\n", 'Sp2xSp2', Sp2Sp2_Sp4, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SU2', SU2_Sp4, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SU2xSU2xSU2', SU2SU2SU2_Sp4, 'n'):
   elif rn=5 then
      printf(" %-15s: [%a][%a]\n", 'SU5xU1', SU5U1_Sp5, 'n');
      printf(" This list is imcomplete!");
   end if:
elif dc="D" then
   if rn=4 then
      printf(" %-15s: [%a][%a]\n", 'SU4xU1', SU4U1_SO8, 'n/s'):
      printf(" %-15s: [%a][%a]\n", 'SU2xSU2xSU2xSU2', SU2SU2SU2SU2_SO8, 'n'):
      printf(" %-15s: [%a][%a]\n", 'SO7', SO7_SO8, 'n/s'):
      printf(" %-15s: [%a][%a]\n", 'SU3', SU3_SO8, 'n'):
      printf(" %-15s: [%a][%a]\n", 'Sp2xSU2', Sp2SU2_SO8, 'n'):
   elif rn=5 then
      printf(" %-15s: [%a][%a]\n", 'SO8xU1', SO8U1_SO10, 'n/c/s'):
      printf(" %-15s: [%a][%a]\n", 'SU5xU1', SU5U1_SO10, 'n1/n2/n3/c'):
      printf(" %-15s: [%a][%a]\n", 'SU4xSU2xSU2', SU4SU2SU2_SO10, 'n1/n2/c/s'):
      printf(" %-15s: [%a][%a]\n", 'SO9', SO9_SO10, 'c/s'):
      printf(" %-15s: [%a][%a]\n", 'SO7xSU2', SO7SU2_SO10, 'n/c/s'):
      printf(" %-15s: [%a][%a]\n", 'Sp2xSp2', Sp2Sp2_SO10, 'n'):
      printf(" %-15s: [%a][%a]\n", 'Sp2', Sp2_SO10, 'n'):
   end if:
elif dc="E" then
   if rn=6 then
      printf(" %-15s: [%a][%a]\n", 'SO10xU1', SO10U1_E6, 'n1/n2/c'):
      printf(" %-15s: [%a][%a]\n", 'SU6xSU2', SU6SU2_E6, 'n1/n2/c/s1/s2'):
      printf(" %-15s: [%a][%a]\n", 'SU3xSU3xSU3', SU3SU3SU3_E6, 'n/c/s'):
      printf(" %-15s: [%a][%a]\n", 'SU3', SU3_E6, 'c'):
      printf(" %-15s: [%a][%a]\n", 'G2', G2_E6, 'n'):
      printf(" %-15s: [%a][%a]\n", 'G2xSU3', G2SU3_E6, 'c'):
      printf(" %-15s: [%a][%a]\n", 'Sp4', Sp4_E6, 'n'):
      printf(" %-15s: [%a][%a]\n", 'F4', F4_E6, 'c/s'):
   end if:
elif dt=F4 then
      printf(" %-15s: [%a][%a]\n", 'SO9', SO9_F4, 'n/s'):
      printf(" %-15s: [%a][%a]\n", 'SU3xSU3', SU3SU3_F4, 'n/s'):
      printf(" %-15s: [%a][%a]\n", 'Sp3xSU2', Sp3SU2_F4, 'n/c'):
      printf(" %-15s: [%a][%a]\n", 'G2xSU2', G2SU2_F4, 'n'):
elif dt=G2 then
      printf(" %-15s: [%a][%a]\n", 'SU3', SU3_G2, 'n'):
end if:
end proc:

#####>Projection matrix maker for regular sublagebras
####>proc: Jmatrix
#% Jmatrix(n) => J[i,j]=delta[i,n+1-j], i,j=1..n
n:='n':
Jmatrix:=proc(n)
  local i, JM:
  
  if not (type(n,integer) and n>=0) then
    error("Usage: Jmatrix( n: non-negative integer )");
  end if:
  JM:=Matrix(n):
  i:='i':
  for i from 1 to n do
    JM[i,n+1-i]:=1:
  end do:
  return JM;
end proc:

####>proc:Pmatrix
#% Pmatrix(plist::permutation list of [1,...,n]) => Permutation Matrix

Pmatrix:=proc()
  local plist, n, PM, i:
  
  if _npassed<1 then
    error("Usage: Pmatrix([i_1,...,i_n])");
  end if:
  plist:=_passed[1]:
  if not (type(plist,list) and nops(plist)>0) then
     error("The argument must of a permutation list of [1,...,n]");
  end if:
  n:=nops(plist):
  PM:=Matrix(n):
  i:='i':
  for i from 1 to n do
     PM[i,plist[i]]:=1:
  end do:
  
  return PM;
end proc:
  
  
####>proc: mkembM
#% mkembM( dt, nodepos, type) => table(["subalgebra"=dts,"embM"=embM)
#% nodepos = the number of the node eliminated from the Dynkin diagram
#% type=1 => eliminate one node from the Dynkin diagram and add U1
#% type=2 => use the extended Dynkin diagram

mkembM:=proc()
  local dt, nodepos, type, rank,rank1, Qvec,
  embM, IDM,CM,CM1,SR, SR1,dts,
  i,j, q, eqs, qsol,dt1, dt2, M1,M2,M3,M4,swapM;
  
  if _npassed<2 then
    error("Usage: mkEmbM( dt, nodepos, type).  type=1 => with U1, type=2 => no U1.");
  end if:
  dt:=resolvedt(_passed[1]):
  rank:=dt[2]:
  rank1:=rank-1:
  nodepos:=_passed[2]:
  i:='i':
  if not member(nodepos,[seq(i,i=1..rank)]) then
    error("node number should be a positive integer <=rank.");
  end if:
  type:=_passed[3]:
  if not member(type,[1,2]) then
    error("type should be 1 or 2");
  end if:
  
  IDM:=IdentityMatrix(rank1):
###>type=1
  if type=1 then
    if rank1>0 then
       embM[D]:=Matrix([IDM[1..rank1,1..(nodepos-1)], 
                       Matrix(rank1,1), 
                      IDM[1..rank1,nodepos..rank1] ]):
      i:='i': j:='j': q:='q':
      Qvec:=[seq(q[i],i=1..(nodepos-1)),1,seq(q[j],j=nodepos..rank1)]:
      i:='i':
      eqs:=[seq(embM[D][i,1..rank].Gmetric(dt)^(-1).Vector(Qvec),i=1..rank1)]:
      qsol:=solve(eqs,[seq(q[i],i=1..rank1)]):
      i:='i': Qvec:=eval(Qvec,qsol[1]):
      embM[D]:=<embM[D],Matrix(Qvec)>:
    else
      embM[D]:=Matrix([1]):
    end if:
    CM1:=[]: SR1:=[]:
    dts:=[]:
##>A
    if dt[1]="A" then
      if nodepos>1 then
        dt1:=abbrevdt(["A",nodepos-1]):
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
      end if:
      if nodepos<rank then
        dt1:=abbrevdt(["A",rank-nodepos]):
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
      end if:
##>B
    elif dt[1]="B" then
      if rank<2 then 
        error("Rank >=2 for the Dynkin type B");
      end if:
      if nodepos>1 then
        dt1:=abbrevdt(["A",nodepos-1]):
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
      end if:
      if nodepos<=rank-2 then
        dt1:=abbrevdt(["B",rank-nodepos]):
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
      elif nodepos=rank-1 then
        dt1:='A1':
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
      end if:
##>C
    elif dt[1]="C" then
       if rank<2 then 
        error("Rank >=2 for the Dynkin type C");
      end if:
     if nodepos>1 then
        dt1:=abbrevdt(["A",nodepos-1]):
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
      end if:
      if nodepos<=rank-2 then
        dt1:=abbrevdt(["C",rank-nodepos]):
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
      elif nodepos=rank-1 then
        dt1:='A1':
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
      end if:
##>D
    elif dt[1]="D" then
      if rank<3 then 
        error("Rank >=3 for the Dynkin type D");
      end if:
      if nodepos>1 then
        if nodepos=rank-1 then
          dt1:=abbrevdt(["A",rank1]):
        else
          dt1:=abbrevdt(["A",nodepos-1]):
        end if:
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
      end if:
      if nodepos<=rank-4 then
        dt1:=abbrevdt(["D",rank-nodepos]):
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
      elif nodepos=rank-3 then
#       dt1:=abbrevdt(["A",rank-nodepos]):
        dt1:=abbrevdt(["D",rank-nodepos]):
        dts:=[op(dts),dt1]:
        SR1:=[op(SR1),SimpleRootsbyH(dt1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(dt1))]:
#        swapM:= <0,1,0|1,0,0|0,0,1>:
#       i:='i': swapM:=DiagonalMatrix([seq(1,i=1..(rank-4)),swapM,1]):
#       embM[D]:=swapM.embM[D]:
      elif nodepos=rank-2 then
        dts:=[op(dts),'A1','A1']:
        SR1:=[op(SR1),SimpleRootsbyH(A1),SimpleRootsbyH(A1)]:
        CM1:=[op(CM1),Transpose(Cmatrix(A1)),Transpose(Cmatrix(A1))]:
      end if:
##>G2
    elif dt=["G",2] then
      dts:=['A1']:
      SR1:=[SimpleRootsbyH(A1)]:
      CM1:=[Transpose(Cmatrix(A1))]:
##>F4
    elif dt=["F",4] then 
      if nodepos=1 then
         dts:=['C3']:
         swapM:=DiagonalMatrix([<0,0,1|0,1,0|1,0,0>,1]):
         embM[D]:=swapM.embM[D]:
         SR1:=[SimpleRootsbyH('C3')]:
         CM1:=[Transpose(Cmatrix('C3'))]:
      elif nodepos=2 then
         dts:=['A1','A2']:
         SR1:=[SimpleRootsbyH(A1),SimpleRootsbyH(A2)]:
        CM1:=[Transpose(Cmatrix(A1)),Transpose(Cmatrix(A2))]:
      elif nodepos=3 then
         dts:=['A2','A1']:
         SR1:=[SimpleRootsbyH(A2),SimpleRootsbyH(A1)]:
         CM1:=[Transpose(Cmatrix(A2)),Transpose(Cmatrix(A1))]:
      elif nodepos=4 then
         dts:=['B3']:
         SR1:=[SimpleRootsbyH('B3')]:
         CM1:=[Transpose(Cmatrix('B3'))]:
      end if;
##>E6
    elif dt=["E",6] then
      if nodepos=1 then
          dts:=['D5']:
          swapM:=DiagonalMatrix([<0,0,0,1|0,0,1,0|0,1,0,0|1,0,0,0>,1,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('D5')]:
          CM1:=[Transpose(Cmatrix('D5'))]:
      elif nodepos=2 then
          dts:=['A1','A4']:
          swapM:=DiagonalMatrix([<1,0,0,0|0,0,0,1|0,0,1,0|0,1,0,0>,1,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('A1'),SimpleRootsbyH('A4')]:
          CM1:=[Transpose(Cmatrix('A1')),Transpose(Cmatrix('A4'))]:
      elif nodepos=3 then
           dts:=['A2','A2','A1']:
          SR1:=[SimpleRootsbyH('A2'),SimpleRootsbyH('A2'),SimpleRootsbyH('A1')]:
          CM1:=[Transpose(Cmatrix('A2')),Transpose(Cmatrix('A2')),Transpose(Cmatrix('A1'))]:
      elif nodepos=4 then
          dts:=['A4','A1']:
          swapM:=DiagonalMatrix([1,1,1,<0,1|1,0>,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('A4'),SimpleRootsbyH('A1')]:
          CM1:=[Transpose(Cmatrix('A4')),Transpose(Cmatrix('A1'))]:
      elif nodepos=5 then
          dts:=['D5']:
          SR1:=[SimpleRootsbyH('D5')]:
          CM1:=[Transpose(Cmatrix('D5'))]:
      elif nodepos=6 then
          dts:=['A5']:
          SR1:=[SimpleRootsbyH('A5')]:
          CM1:=[Transpose(Cmatrix('A5'))]:
      end if:
##>E7
    elif dt=["E",7] then
      if nodepos=6 then
          dts:=['E6']:
          SR1:=[SimpleRootsbyH('E6')]:
          CM1:=[Transpose(Cmatrix('E6'))]:
      elif nodepos=5 then
          dts:=['D5','A1']:
          swapM:=<0,1|1,0>:
          swapM:=DiagonalMatrix([1,1,1,1,swapM,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('D5'),SimpleRootsbyH('A1')]:
          CM1:=[Transpose(Cmatrix('D5')),Transpose(Cmatrix('A1'))]:
      elif nodepos=4 then
          dts:=['A4','A2']:
          swapM:=DiagonalMatrix([1,1,1,<0,1,0|0,0,1|1,0,0>,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('A4'),SimpleRootsbyH('A2')]:
          CM1:=[Transpose(Cmatrix('A4')),Transpose(Cmatrix('A2'))]:
       elif nodepos=3 then
           dts:=['A2','A3','A1']:
          SR1:=[SimpleRootsbyH('A2'),SimpleRootsbyH('A3'),SimpleRootsbyH('A1')]:
          CM1:=[Transpose(Cmatrix('A2')),Transpose(Cmatrix('A3')),Transpose(Cmatrix('A1'))]:
      elif nodepos=2 then
          dts:=['A1','A5']:
          swapM:=<0,0,0,1|0,0,1,0|0,1,0,0|1,0,0,0>:
          swapM:=DiagonalMatrix([1,swapM,1,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('A1'),SimpleRootsbyH('A5')]:
          CM1:=[Transpose(Cmatrix('A1')),Transpose(Cmatrix('A5'))]:
      elif nodepos=1 then
          dts:=['D6']:
          swapM:=<0,0,0,0,1|0,0,0,1,0|0,0,1,0,0|0,1,0,0,0|1,0,0,0,0>:
          swapM:=DiagonalMatrix([swapM,1,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('D6')]:
          CM1:=[Transpose(Cmatrix('D6'))]:
      elif nodepos=7 then
          dts:=['A6']:
          SR1:=[SimpleRootsbyH('A6')]:
          CM1:=[Transpose(Cmatrix('A6'))]:
      end if:
##>E8
    elif dt=["E",8] then
     if nodepos=7 then
          dts:=['E7']:
          SR1:=[SimpleRootsbyH('E7')]:
          CM1:=[Transpose(Cmatrix('E7'))]:
      elif nodepos=6 then
          dts:=['E6','A1']:
          swapM:=DiagonalMatrix([1,1,1,1,1,<0,1|1,0>,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('A1'),SimpleRootsbyH('E6')]:
          CM1:=[Transpose(Cmatrix('A1')),Transpose(Cmatrix('E6'))]:
      elif nodepos=5 then
          dts:=['D5','A2']:
          swapM:=<0,1,0|0,0,1|1,0,0>:
          swapM:=DiagonalMatrix([1,1,1,1,swapM,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('D5'),SimpleRootsbyH('A2')]:
          CM1:=[Transpose(Cmatrix('D5')),Transpose(Cmatrix('A2'))]:
      elif nodepos=4 then
          dts:=['A4','A3']:
          swapM:=<0,1,0,0|0,0,1,0|0,0,0,1|1,0,0,0>:
          swapM:=DiagonalMatrix([1,1,1,swapM,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('A4'),SimpleRootsbyH('A3')]:
          CM1:=[Transpose(Cmatrix('A4')),Transpose(Cmatrix('A3'))]:
      elif nodepos=3 then
           dts:=['A2','A4','A1']:
          SR1:=[SimpleRootsbyH('A2'),SimpleRootsbyH('A4'),SimpleRootsbyH('A1')]:
          CM1:=[Transpose(Cmatrix('A2')),Transpose(Cmatrix('A4')),Transpose(Cmatrix('A1'))]:
     elif nodepos=2 then
          dts:=['A1','A6']:
          swapM:=<0,0,0,0,1|0,0,0,1,0|0,0,1,0,0|0,1,0,0,0|1,0,0,0,0>:
          swapM:=DiagonalMatrix([1,swapM,1,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('A6'),SimpleRootsbyH('A1')]:
          CM1:=[Transpose(Cmatrix('A6')),Transpose(Cmatrix('A1'))]:
      elif nodepos=1 then
          dts:=['D7']:
          swapM:=Matrix(6):
          i:='i': for i from 1 to 6 do swapM[i,7-i]:=1: end do:
          swapM:=DiagonalMatrix([swapM,1,1]):
          embM[D]:=swapM.embM[D]:
          SR1:=[SimpleRootsbyH('D7')]:
          CM1:=[Transpose(Cmatrix('D7'))]:
      elif nodepos=8 then
          dts:=['A7']:
          SR1:=[SimpleRootsbyH('A7')]:
          CM1:=[Transpose(Cmatrix('A7'))]:
      end if:
    end if:
    dts:=[op(dts),'U1']:
    SR1:=DiagonalMatrix([op(SR1),1]):
    CM1:=DiagonalMatrix([op(CM1),1]):
    CM:=Transpose(Cmatrix(dt)):
    SR:=SimpleRootsbyH(dt):
    embM[S]:=CM1^(-1).embM[D].CM:
    embM[H]:=SR1.embM[S].SR^(-1):
###> type=2
  elif type=2 then
##>A
    if dt[1]="A" then
      error("No type 2 regular subalgebra for the Dynkin type A");
##>B
    elif dt[1]="B" then
      if rank<2 then 
        error("Rank >=2 for the Dynkin type B");
      end if:
       if nodepos=1 then
           error("No type 2 regular subalgebra");
       end if:
#>r=2
      if rank=2 then
           dts:=['A1','A1']:
           embM[S]:=<1,0|-1/2,-1/2>:
           SR1:=SimpleRootsbyH('A1'):
           SR1:=[SR1,SR1]:
           CM1:=Transpose(Cmatrix('A1')):
           CM1:=[CM1,CM1]:
#>r>=3
      else
##% p=2
         if nodepos=2 then
           if rank=3 then 
              dt1:='A1':
           else
              dt1:=abbrevdt(["B",rank-2]):
           end if:
           dts:=['A1',dt1,'A1']:
           M1:=Matrix(1,rank):  M1[1,1]:=1: M1[1,2]:=-1/2:
           M2:=Matrix([Matrix(rank-2,1),Matrix(rank-2,1,fill=-1),IdentityMatrix(rank-2)]):
           M3:=Matrix(1,rank): M3[1,2]:=-1/2:
           embM[S]:=<M1,M2,M3>:
           SR1:=SimpleRootsbyH('A1'):
           if rank=3 then
             SR1:=[SR1,SR1,SR1]:
           else
             SR1:=[SR1,SimpleRootsbyH(dt1),SR1]:
           end if:
           CM1:=Transpose(Cmatrix('A1')):
           if rank=3 then
             CM1:=[CM1,CM1,CM1]:
           else
             CM1:=[CM1,Transpose(Cmatrix(dt1)),CM1]:
           end if:
##% p=r
         elif nodepos=rank then
           if rank=3 then
             dt1:='A3':
             embM[S]:=<1,0,0|0,1,0|-1/2,-1,-1/2>:
           else
             dt1:=abbrevdt(["D",rank]):
             M1:=Matrix(rank,1,fill=-1): 
             M1[rank-1]:=-1/2: M1[rank]:=-1/2:
             embM[S]:=Matrix([<Jmatrix(rank-1),Matrix(1,rank-1)>,M1]):
           end if:
           dts:=[dt1]:
           SR1:=[SimpleRootsbyH(dt1)]:
           CM1:=[Transpose(Cmatrix(dt1))]:
##% p=3
         elif nodepos=3 then
           if rank=4 then
             dt1:='A1':
           else 
             dt1:=abbrevdt(["B",rank-3]):
           end if:
           M1:=<IdentityMatrix(2),Matrix(rank-2,2)>:
           M2:=Matrix(rank,1,fill=-1): 
           M2[1,1]:=-1/2: M2[3,1]:=-1/2:
           M3:=<Matrix(3,rank-3),IdentityMatrix(rank-3)>:
           embM[S]:=Matrix([M1,M2,M3]):
           dts:=['A3',dt1]:
           SR1:=[SimpleRootsbyH('A3'),SimpleRootsbyH(dt1)]:
           CM1:=[Transpose(Cmatrix('A3')),Transpose(Cmatrix(dt1))]:
##% p=4, .., r-1
         else
           dt1:=abbrevdt(["D",nodepos]): 
           if nodepos=rank-1 then
              dt2:='A1':
           else
              dt2:=abbrevdt(["B",rank-nodepos]):
           end if:
           dts:=[dt1,dt2]:
           M1:=<Jmatrix(nodepos-1),Matrix(rank-nodepos+1,nodepos-1)>:
           M2:=Matrix(rank,1,fill=-1):
           M2[nodepos-1,1]:=-1/2:
           M2[nodepos,1]:=-1/2:
           M3:=<Matrix(nodepos,rank-nodepos),IdentityMatrix(rank-nodepos)>:
           embM[S]:=Matrix([M1,M2,M3]):
           SR1:=[SimpleRootsbyH(dt1),SimpleRootsbyH(dt2)]:
           CM1:=[Transpose(Cmatrix(dt1)),Transpose(Cmatrix(dt2))]:
         end if:
       end if:
##>C
    elif dt[1]="C" then
       if rank<3 then 
         if rank<2 then
            error("Rank >=2 for the Dynkin type C");
         else
            error("C2[ab] = B2[ba]. Try B2.");
         end if:
       end if:
       if nodepos=rank then
           error("No type 2 regular subalgebra");
       end if:
       
       if nodepos=1 then
         dt1:='A1':
       elif nodepos=2 then
         dt1:='B2':
       else
         dt1:=abbrevdt(["C",nodepos]):
       end if:
       if nodepos=rank-1 then
         dt2:='A1':
       elif nodepos=rank-2 then
         dt2:='B2':
       else
         dt2:=abbrevdt(["C",rank-nodepos]):
       end if:
       dts:=[dt1,dt2]:
       
       SR1:=[SimpleRootsbyH(dt1),SimpleRootsbyH(dt2)]:
       CM1:=[Transpose(Cmatrix(dt1)),Transpose(Cmatrix(dt2))]:

       M1:=<Jmatrix(nodepos-1),Matrix(rank-nodepos+1,nodepos-1)>:
       M2:=Matrix(rank,1,fill=-1):
       M2[nodepos,1]:=-1/2: M2[rank,1]:=-1/2:
       M3:=<Matrix(nodepos,rank-nodepos),IdentityMatrix(rank-nodepos)>:
       embM[S]:=Matrix([M1,M2,M3]):
       if nodepos=2 then
         i:='i': swapM:=Pmatrix([2,1,seq(i,i=3..rank)]):
         embM[S]:=swapM.embM[S]:
       end if:
       if nodepos=rank-2 then
         i:='i': swapM:=Pmatrix([seq(i,i=1..(rank-2)),rank,rank-1]):
         embM[S]:=swapM.embM[S]:
       end if:
##>D
     elif dt[1]="D" then
       if rank<4 then 
         if rank<3 then
           error("Rank >=3 for the Dynkin type D");
         else
           error("There is no type-2 regular subalgebra for D3=A3");
         end if:
       end if:
       if member(nodepos,{1,rank-1,rank}) then
           error("No type 2 regular subalgebra");
       end if:
       
       if nodepos=2 then
         dts:=['A1','A1']:
       elif nodepos=3 then
         dts:=['A3']:
       else
         dts:=[abbrevdt(["D",nodepos])]:
       end if:
       if nodepos=rank-2 then
         dts:=[op(dts),'A1','A1']:
       elif nodepos=rank-3 then
         dts:=[op(dts),'A3']:
       else
         dts:=[op(dts),abbrevdt(["D",rank-nodepos])]:
       end if:
       
       SR1:=map(SimpleRootsbyH,dts):
       CM1:=map(Transpose,map(Cmatrix,dts)):

       M1:=<Jmatrix(nodepos-1),Matrix(rank-nodepos+1,nodepos-1)>:
       M2:=Matrix(rank,1,fill=-1):
       M2[nodepos-1,1]:=-1/2: M2[nodepos,1]:=-1/2:
       M2[rank-1,1]:=-1/2: M2[rank,1]:=-1/2:
       M3:=<Matrix(nodepos,rank-nodepos),IdentityMatrix(rank-nodepos)>:
       embM[S]:=Matrix([M1,M2,M3]):
       if nodepos=3 then
         i:='i': swapM:=Pmatrix([2,1,seq(i,i=3..rank)]):
         embM[S]:=swapM.embM[S]:
       end if:
       if nodepos=rank-3 then
         i:='i': swapM:=Pmatrix([seq(i,i=1..(rank-3)),rank-1,rank-2,rank]):
         embM[S]:=swapM.embM[S]:
       end if:
##>G2
    elif dt=["G",2] then
       if nodepos=1 then
         dts:=['A1','A1']:
         embM[S]:=<-1/2,-3/2|0,1>:
       elif nodepos=2 then
         dts:=['A2']:
         embM[S]:=<0,1|-1/3,-2/3>:
       end if:
       SR1:=map(SimpleRootsbyH,dts):
       CM1:=map(Transpose,map(Cmatrix,dts)):
##>F4
    elif dt=["F",4] then
       if nodepos=1 then
         dts:=['A1','C3']:
         embM[S]:=<-1/2,-3/2,-2,-1|0,1,0,0|0,0,1,0|0,0,0,1>:
         swapM:=Pmatrix([1,4,3,2]):
         embM[S]:=swapM.embM[S]:
       elif nodepos=2 then
         dts:=['A2','A2']:
         embM[S]:=<0,1,0,0|-1/3,-2/3,-4/3,-2/3|0,0,1,0|0,0,0,1>:
       elif nodepos=3 then
         dts:=['A3','A1']:
         embM[S]:=<0,1,0,0|0,0,1,0|-1/4,-1/2,-3/4,-1/2|0,0,0,1>:
       elif nodepos=4 then
         dts:=['B4']:
         embM[S]:=<0,1,0,0|0,0,1,0|0,0,0,1|-1/2,-1,-3/2,-2>:
       end if:
       SR1:=map(SimpleRootsbyH,dts):
       CM1:=map(Transpose,map(Cmatrix,dts)):
##>E6      
    elif dt=["E",6] then
       if nodepos=2 then
         dts:=['A1','A5']:
         M1:=Matrix(6,1): M1[1,1]:=1:
         M2:=<-1/2,-1/2,-1,-3/2,-1,-1/2>:
         M3:=<Matrix(1,3),Jmatrix(3),Matrix(2,3)>:
         M4:=Matrix(6,1): M4[5,1]:=1:
         embM[S]:=Matrix([M1,M2,M3,M4]):
       elif nodepos=4 then
         dts:=['A5','A1']:
         M1:=<IdentityMatrix(3),Matrix(3)>:
         M2:=<-1/2,-1,-3/2,-1,-1/2,-1/2>:
         M3:=<Matrix(3,2),<0,0,1|1,0,0>>:
         embM[S]:=Matrix([M1,M2,M3]):
       elif nodepos=6 then
         dts:=['A5','A1']:
         M1:=<IdentityMatrix(5),Matrix(1,5)>:
         M2:=<-1/2,-1,-3/2,-1,-1/2,-1/2>:
         embM[S]:=Matrix([M1,M2]):
       elif nodepos=3 then
         dts:=['A2','A2','A2']:
         M1:=<IdentityMatrix(2),Matrix(4,2)>:
         M2:=<-1/3,-2/3,-2/3,-1/3,-2/3,-1/3>:
         M3:=<Matrix(2,3),IdentityMatrix(3),Matrix(1,3)>:
         embM[S]:=Matrix([M1,M2,M3]):
       else
         error("Wrong node position");
       end if:
       SR1:=map(SimpleRootsbyH,dts):
       CM1:=map(Transpose,map(Cmatrix,dts)):
    elif dt=["E",7] then
       if nodepos=1 then
         dts:=['D6','A1']:
         M1:=<-1/2,-1,-3/2,-2,-3/2,-1,-1/2>:
         M2:=<Jmatrix(5),Matrix(2,5)>:
         M3:=Matrix(7,1): M3[6,1]:=1:
         embM[S]:=Matrix([M1,M2,M3]):
       elif nodepos=5 then
         dts:=['D6','A1']:
         M1:=<Matrix(1,4),IdentityMatrix(4),Matrix(2,4)>:
         M2:=<-1/2,-1,-3/2,-2,-3/2,-1,-1/2>:
         M3:=<Matrix(5,2),Jmatrix(2)>:
         embM[S]:=Matrix([M1,M2,M3]):
       elif nodepos=2 then
         dts:=['A2','A5']:
         M1:=Matrix(7,1): M1[2,1]:=1:
         M2:=<-1/3,-2/3,-1/3,-2/3,-1,-4/3,-2/3>:
         M3:=<Matrix(2,4),Jmatrix(4),Matrix(1,4)>:
         M4:=Matrix(7,1): M4[7,1]:=1:
         embM[S]:=Matrix([M1,M2,M3,M4]):
       elif nodepos=4 then
         dts:=['A5','A2']:
         M1:=<Matrix(1,3),IdentityMatrix(3),Matrix(3)>:
         M2:=<-1/3,-2/3,-1,-4/3,-2/3,-2/3,-1/3>:
         M3:=<Matrix(4,3),<0,1,0|0,0,1|1,0,0>>:
         embM[S]:=Matrix([M1,M2,M3]):
       elif nodepos=3 then
         dts:=['A3','A3','A1']:
         M1:=<<0,1,0|0,0,1>,Matrix(4,2)>:
         M2:=<-1/4,-1/2,-3/4,-3/4,-1/2,-1/4,-1/2>:
         M3:=<Matrix(3,4),IdentityMatrix(4)>:
         embM[S]:=Matrix([M1,M2,M3]):
       elif nodepos=7 then
         dts:=['A7']:
         M1:=<Matrix(1,6),IdentityMatrix(6)>:
         M2:=<-1/2,-1,-3/2,-2,-3/2,-1,-1/2>:
         embM[S]:=Matrix([M1,M2]):
       else
         error("Wrong node position");
       end if:
       SR1:=map(SimpleRootsbyH,dts):
       CM1:=map(Transpose,map(Cmatrix,dts)):
    elif dt=["E",8] then
       if nodepos=1 then
          dts:=['D8']:
          M1:=<-1/2,-1,-3/2,-2,-5/2,-3,-2,-3/2>:
          M2:=<Matrix(1,6),Jmatrix(6),Matrix(1,6)>:
          M3:=Matrix(8,1):M3[8,1]:=1:
          embM[S]:=Matrix([M1,M2,M3]):
       elif nodepos=2 then
          dts:=['A1','A7']:
          M1:=Matrix(8,1): M1[1,1]:=1:
          M2:=<-1/2,-1/4,-1/2,-3/4,-1,-5/4,-3/2,-3/4>:
          M3:=<Matrix(2,5),Jmatrix(5),Matrix(1,5)>:
          M4:=Matrix(8,1): M4[8,1]:=1:
          embM[S]:=Matrix([M1,M2,M3,M4]):
       elif nodepos=3 then
          dts:=['A2','A5','A1']:
          M1:=<IdentityMatrix(2),Matrix(6,2)>:
          M2:=<-1/3,-2/3,-5/6,-2/3,-1/2,-1/3,-1/6,-1/2>:
          M3:=<Matrix(2,4),IdentityMatrix(4),Matrix(2,4)>:
          M4:=Matrix(8,1); M4[8,1]:=1:
          embM[S]:=Matrix([M1,M2,M3,M4]):
       elif nodepos=4 then
          dts:=['A4','A4']:
          M1:=<IdentityMatrix(3),Matrix(5,3)>:
          M2:=<-2/5,-4/5,-6/5,-3/5,-4/5,-3/5,-2/5,-1/5>:
          M3:=<Matrix(4,3),IdentityMatrix(3),Matrix(1,3)>:
          M4:=Matrix(8,1): M4[4,1]:=1:
          embM[S]:=Matrix([M1,M2,M3,M4]):
       elif nodepos=5 then
          dts:=['D5','A3']:
          M1:=<IdentityMatrix(4),Matrix(4)>:
          M2:=<-1/2,-1,-3/2,-5/4,-3/4,-3/4,-1/2,-1/4>:
          M3:=<Matrix(5,2),IdentityMatrix(2),Matrix(1,2)>:
          M4:=Matrix(8,1): M4[5,1]:=1:
          embM[S]:=Matrix([M1,M2,M3,M4]):
       elif nodepos=6 then
          dts:=['E6','A2']:
          M1:=<IdentityMatrix(5),Matrix(3,5)>:
          M2:=<-2/3,-4/3,-2,-5/3,-4/3,-1,-2/3,-1/3>:
          M3:=<Matrix(5,2),Jmatrix(2),Matrix(1,2)>:
          embM[S]:=Matrix([M1,M2,M3]):
       elif nodepos=7 then
          dts:=['E7','A1']:
          M1:=<IdentityMatrix(6),Matrix(2,6)>:
          M2:=<-1,-2,-3,-5/2,-2,-3/2,-3/2,-1/2>:
          M3:=Matrix(8,1): M3[7,1]:=1:
          embM[S]:=Matrix([M1,M2,M3]):
       elif nodepos=8 then
          dts:=['A8']:
          M1:=<IdentityMatrix(7),Matrix(1,7)>:
          M2:=<-2/3,-4/3,-2,-5/3,-4/3,-1,-2/3,-1/3>:
          embM[S]:=Matrix([M1,M2]):
       else
          error("Wrong node position");
       end if:
       SR1:=map(SimpleRootsbyH,dts):
       CM1:=map(Transpose,map(Cmatrix,dts)): 
    else 
       error("Illegal Dynkin type");
    end if:
    SR:=SimpleRootsbyH(dt):
    SR1:=DiagonalMatrix(SR1):
    CM:=Transpose(Cmatrix(dt)):
    CM1:=DiagonalMatrix(CM1):
    embM[D]:=CM1.embM[S].CM^(-1):
    embM[H]:=SR1.embM[S].SR^(-1):
###>other types
  else 
    error("Illegal subalgebra type");
  end if:
  return table(["subalgebra"=dts,"embM"=embM]);
end proc:



######>help procedures
#####>proc:procinfo
procinfo:=proc()

####>Common procedures
###>Abbreviations
printf("*Abbreviation: \n");
printf("    dt=Dynkin type (e.g. A4, [%a,4])\n","A");
printf("    dw=Dynkin label for the weight (e.g., [0,-1,2,-1])\n");
printf("    hdw=highest Dynkin label for the weight (e.g., [1,0,0,1])\n");
printf("    basis=%a\n", "S/D");
###>Metrics
printf("*Metrics: \n");
printf("   Kmetric(dt) => Killing metric matrix\n");
printf("   Cmatrix(dt) => Cartan matrix\n");
printf("   Gmetric(dt) => Killing metric matrix in the Dynkin basis\n");
###>Inner product
printf("*Inner products: \n");
printf("   IPDB(dw1,dw2,dt) => Inner product of Dynkin weights\n");
printf("   RIP(w1,w2,dt) => Inner product of weights in the simple root basis\n");
printf("   RCP(w,r,dt) => Cartan product <w,r>\n");
###>SRB <=> DB
printf("*Basis change: \n");
printf("   DBtoSRB: coordinate trf of weight. Dynkin basis => Simple root basis\n");
printf("      usage: DBtoSRB(dw,dt)\n");
printf("   SRBtoDB: coordinate trf of weight. Simple root basis to Dynkin basis \n");
printf("      usage: SRBtoDB(w,dt)\n");
###>RootSystem
printf("*Root System: \n");
printf("   HighestRoot(dt) => the highest root in the simple root basis\n");
printf("   Rlevel(w)=> the level of a root w in the simple root basis\n");
printf("   RootSystem: giving the root system in the simple root/Dynkin basis\n");
printf("     usage:  RootSystem(dt[, basis, output=0/1])\n");
###>Structure constant
printf("*STructure constant: \n");
printf("   StrConstWB( dt[,sw]) => Weyl Basic CCR coeff. NN=table(N[wv1,wv2])\n");
printf("      sw=0 => no message. sw=1 => display the result.\n");
printf("   StrConst(NNtable) => show the list of non-vanishing SC\n");
###>WeightSystem
printf("*Weight System: \n");
printf("   WeightSystem: giving the weight system DWS (in the Dynkin)/WS (in the simple root basis)\n");
printf("      usage:  WeightSystem(dt,hdw, [, basis, sw=0/1/2])\n");
printf("      sw=0 => output=DWS. No monitor output,\n");
printf("      sw=1 => output=DWS. Display the weight system,\n");
printf("      sw=2 => output= the flat weight list.\n");
printf("   printWS(WS/DWS): gives a formated display of a weight system\n");
printf("   findHighestWeight: find the highest weight of the irrep to which a given weight belongs to in Dynkin basis\n");
printf("      usage: findHighestWeight(dw, dt)\n");
###>Weyl trf
printf("*Weyl transformation: \n");
printf("   WeylTrf(dt, rv[, basis]) => Weyl trf matrix in the basis wrt the root vector rv\n");
printf("      rv=a root vector in the simple root basis\n");
###>ProductOfRep
printf("*Irr. decomposition of a product of representations: \n");
printf("    ProductOfRep(dt, hdw1, hdw2 [, sw=0/1, basis])\n");
printf("    sw=0 (default) => output = Irrep. list only\n");
printf("    sw=1 => output = Irrep. list with a formated display\n");
printf("    sw=2 => output = Irrep. list with display of the composite weight system\n");
printf("            in the Simple Root Basis (basis=\"S\") or the Dynkin Basis(basis=\"D\")\n");
###>Listing maximal quasi-semisimple subalgebras for low rank algebras
printf("*Listing maximal quasi-semisimple subalgebras for low rank Lie algebras\n");
printf("    MaxSubGlist(dt) => [[the list of regular subalgs],[the list of special subalgs]]\n");
printf("    MaxSubGlist0(dt) => [ the flat list of all maximal subalgs].\n");
###>Constructing an embedding matrix of any regular quasi-semisimple subalgebra
printf("*Constructing an embedding matrix for any regular quasi-semisimple subalgebra of any simple Lie algebra\n");  
printf("    mkembM(dt,nodepos,type) => table([\"subalgebra\"=dts,\"embM\"=embM])\n");
printf("    dts= the list of the Dynkin types of the resulting subalgebra. e.g. [A4,U1] for D5 \n");
printf("    embM = table([D=embM[D],S=embM[S],H=embM[H]])\n");
printf("         embM[D]=projection matrix for Dynkin labels, \n");
printf("         embM[S]=matrix specifying the pull back of simple roots,\n"); 
printf("         embM[H]=matrix specifying the embedding of the Cartan subalgebra\n");
printf("      nodepos= a node position to remove from the (extended) Dynkin diagram\n");
printf("      type=1 => remove one node from the Dynkin diagram and add U1 \n");
printf("      type=2 => use the extended Dynkin diagram \n");

###>SubGrdm
printf("*Irr. decomposition of a irrep by subalgebras: \n");
printf("    SubGrdm(dt,hwd dts,embMD[,sw] => DWS/WL\n");
printf("      dts=[A3,B4, U1]\n");
printf("      embMD=embedding matrix in the Dynkin basis\n");
printf("      sw=0 => output=Irrep list. No monitor output,\n");
printf("      sw=1 => output=Irrep list with monitor display,\n");
printf("      sw=2 => output=Irrep list. Display the projected Dynkin weight system,\n");
printf("      sw=3 => output=Irrep list. Display the projeced Dynkin weight system together with the original weight system.\n");
printf("   SubGrdm allow  multiple U1 factors.\n");
##%

####>Symmetry breaking
###>proc:MaySubG
printf("*Embedding of algebras.\n");
printf("   MaySubG(dts1,dts0) => maplist : [dts0[1]=[dts1[1],..], dts0[2]=[dts1[3],..]]\n");
printf("      dts0=a target list of simple algebras,\n");
printf("      dts1=a list of may be subalgebras.\n"); 
###>proc:SBpattern
printf("*Symmetry breaking chain list\n");
printf("   SBpattern(dt,dts[,outsw=0/1]) =>the SBchain list [ [[SBlist,[[],Unbrknlist],...]\n");
printf("      dt=the Dynkin type of the initial algebra,\n");
printf("      dts= a list of Dynkin types for the final algebra,\n");
printf("      outsw=0 => output=the SBchain list, no display,\n");
printf("      outsw=1 => + monitordisplay of the SBchain list\n");
###>proc:mkSBchain
printf("*Symmetry breaking chain making\n");
printf("   mkSBchain(SBClist0) => a complete SB chain list SBClist1.\n");
printf("      SBClist0=[a sequence of incomplete SBC's]\n");
printf("      SBClist1=[a sequence of complete SBC]\n");
printf("        SBC=[a sequence of SBs,[maplist,unbroken alg list]]\n"); 
pritnf("          SBs=[a sequence of SB];  SB=[E6=[D5,U1]](e.g.)\n");
printf("          maplist=[ a seq of alg maps]; alg map= D5=[A2,A1] (e.g.)\n");
printf("          unbroken alg list= [A1,U1] (e.g.)\n");
printf("   e.g. mkSBchain([[E6=[D5,U1],[ [D5=[A2,A1]],[U1] ] ]]\n");
printf("      SBClist1=[ a seq of complete SBCs]\n");
printf("         a complete SBC=SBC with maplist=[] describing a complete SB chain from dt0 to dts\n");
printf("   e.g. [[E6=[D5,U1]],[D5=[A4,U1]],[A4=[A2,A1,U1]],[[],[U1,U1,U1]]]\n");
#% 
end proc:



######>EOF