with(linalg):
#######################################################################
#test checks that the triangular condition on the three integers or   #
#half-integers a, b, c is satisfied.                                  #
#######################################################################
test:=proc(a,b,c)
      local result,result1,result2;
      if (evalb(a+b>=c) and evalb(c>=abs(a-b)))
      then result:=true else result:=false;
      fi;
      end:
#######################################################################
#threej evaluates 3jm-symbols.                                        #
#######################################################################
threej:=proc(j1,j2,j3,m1,m2,m3)
        local xmin,xmax,x,fact,sumx,result,phase;
        if test(j1,j2,j3) and (m1+m2+m3 =0) then
        phase:=(-1)^(j1-j2-m3);
        fact:=sqrt(((j1+j2-j3)!*(j1-m1)!*(j2-m2)!*(j3+m3)!*(j3-m3)!)/
        ((j1+j2+j3+1)!*(j3+j1-j2)!*(j3+j2-j1)!*(j1+m1)!*(j2+m2)!));
        xmin:=max(0,j3-j2-m1);
        xmax:=min(j3+m3,j1-m1,j3+j2-m1);
        sumx:=0;
        for x from xmin to xmax do
sumx:=sumx + ((-1)^(j1-m1-x))*((j1+m1+x)!*(j3+j2-m1-x)!)/
      (x!*(j3+m3-x)!*(j1-m1-x)!*(j2-j3+m1+x)!);
                                od;
result:=simplify(phase*fact*sumx) else result:=0;
fi;
end:

###########################################################################
#ck evaluates the reduced matrix element <a//C^(k)//b>                    #
###########################################################################

ck:=proc(a,b,k)
    local result;
    result:=simplify((-1)^a*sqrt((2*a+1)*(2*b+1))*threej(a,k,b,0,0,0));
    end:
###########################################################################
#ckj evaluates the reduced matrix element <slj//C^(k)//sl'j'>             #
###########################################################################
ckj:=proc(a,b,k,c,d)
     local result,t;
     t:=(1/2);
     result:=simplify(((-1)^(b+c+k+(1/2)))*sqrt((2*b+1)*(2*d+1))*sixj(b,k,d,c,t,a)
*ck(a,c,k));
end: 
###########################################################################
#CLS evaluates the matrix elements <l^2 SL//C^(k).C(k)//l^2 SL>           #
###########################################################################
CLS:=proc(l,k,L)
     local result;
result:=simplify(((-1)^L)*sixj(l,l,k,l,l,L)*((ck(l,l,k)^2)));
end:
###########################################################################
#clj evaluates the matrix elements <(lj)^2 J//C^(k),C^(k)//(lj')^2 J>     #
###########################################################################
clj:=proc(l,k,j1,j2,J)
     local result;     result:=simplify(((-1)^(j1+j2+J))*sixj(j2,j2,J,j1,j1,k)*ckj(l,j1,k,l,j2)
*ckj(l,j1,k,l,j2));
end:
###########################################################################
#cljj1 evaluates the matrix elements <(lj)^2 J//C^(k),C^(k)//(lj,lj+1J>   #
###########################################################################
cljj1:=proc(l,k,j1,J)
     local result;     result:=simplify(sqrt(2)*((-1)^(J+1))*sixj(j1,j1+1,J,j1,j1,k)
*(ckj(l,j1,k,l,j1))*(ckj(l,j1,k,l,j1+1)));
end:
###########################################################################
#cljj2 evaluates the matrix elements <(lj)^2 J//C^(k),C^(k)//(lj-1,ljJ>   #
###########################################################################
cljj2:=proc(l,k,j1,J)
     local result;     result:=simplify(sqrt(2)*((-1)^(J-1))*sixj(j1-1,j1,J,j1,j1,k)
*(ckj(l,j1,k,l,j1))*(ckj(l,j1,k,l,j1-1)));
end:
###########################################################################
#cljj evaluates the matrix elements <(lj1,j1+1)J//C^(k),C^(k)//(lj1,j1+1J>#
###########################################################################
 cljj:= proc(l,j1,k,J)
        local result;
result:=simplify(((-1)^J)*sixj(j1,j1+1,J,j1+1,j1,k)*ckj(l,j1,k,l,j1)
*ckj(l,j1+1,k,l,j1+1));
result:=simplify(result - (sixj(j1+1,j1,J,j1+1,j1,k)*(ckj(l,j1,k,l,j1+1)^2)));
end: 
###########################################################################
#triad evaluates the triangular portion of the formulae for 6j-symbols    #
###########################################################################

triad:=proc(a,b,c)
       local triang;
       triang:=sqrt((((a+b-c)!*(a-b+c)!*(b+c-a)!)/(a+b+c+1)!));
       end:

###########################################################################
#sixj evaluates a 6j-symbol.                                              #
###########################################################################

sixj:=proc(a,b,c,d,e,f)
      local trif,sumj,zmin,zmax,z,result;
      if (test(a,b,c) and test(d,b,f) and test(d,e,c) and test(a,e,f))
      then
      trif:=simplify(triad(a,b,c)*triad(a,e,f)*triad(d,b,f)*triad(d,e,c));
      zmin:=max(a+b+c,a+e+f,d+b+f,d+e+c);
      zmax:=min(a+b+d+e,b+c+e+f,c+a+f+d);
      sumj:=0;
for z from zmin to zmax do
sumj:=sumj + (((-1)^z*(z+1)!)
/((z-a-b-c)!*(z-a-e-f)!*(z-d-b-f)!*(z-d-e-c)!*(a+b+d+e-z)!*(b+c+e+f-z)!
*(c+a+f+d-z)!));
od;
result:=simplify(trif*sumj) else result:=0;
fi;
end:

###########################################################################
#ninej evaluates a 9j-symbol                                              #
###########################################################################

ninej:=proc(a,b,c,d,e,f,h,i,j)
       local x,xmin,xmax,result;
if (test(a,d,h) and test(i,j,h) and test(b,e,i) and test(d,e,f) and
test(c,f,j) and test(c,a,b))
then
       xmax:=min(a+j,i+d,b+f);
       xmin:=max(abs(a-j),abs(i-d),abs(b-f));
       result:=0;
for x from xmin to xmax do
result:=result + ((-1)^(2*x))*(2*x + 1)*sixj(a,d,h,i,j,x)*sixj(b,e,i,d,x,f)
*sixj(c,f,j,x,a,b);
od;
result:=simplify(result)
else result:=0;
fi;
end:

############################################################################

`help/text/njsymbol`:=TEXT(`HELP for njsymbol`,
`This package contains  procedures for calculating 3jm-, 6j- and 9j-symbols`,
`in rational form.The result appears as a fraction times square root factors.`,
`The square root factors can be combined into a single square root using the`,
`the MAPLE command combine("). The procedures can be used for inclusion in `,
`other MAPLE programmes requiring the use of njsymbols.`,
`The three procedures are as follows:-`,
``,
`threej(j1,j2,j3,m1,m2,m3):-`,
``,
` This procedure evaluates the value of a 3jm-symbol involving the three`,
` angular momenta j1, j2, j3 and their projections m1, m2, m3. The arguments`,
` are entered as integers or half-integers as appropriate. The triangular`,
` rules are automatically checked and if not satisfied the symbol is`,
` evaluated to 0.`,
``,
`sixj(a,b,c,d,e,f):-`,
``,
` This procedure evaluates the value of a 6j-symbol involving the six`,
` angular momenta a, b, c, d, e, f. The arguments are entered as integers or`,
` half-integers as appropriate.The triangular rules are automatically checked`,
` and if not satisfied the symbol is evaluated to 0.`,
``,
`ninej(a,b,c,d,e,f,h,i):-`,
``,
` This procedure evaluates the value of a 9j-symbol involving the nine`,
` angular momenta a, b, c, d, e, f, g, h, i. The arguments are entered as`,
` integers or half-integers as appropriate. The triangular rules are`,
` automatically checked and if not satisfied the symbol is evaluated to 0.`,
``,
`EXAMPLES:-`,
``,
`1. threej(3/2,3/2,3,1/2,1/2,-1);`,
``,
`2. sixj(3/2,3/2,3,1/2,7/2,2);`,
``,
`3. ninej(3/2,3/2,3,1/2,7/2,3,2,2,2);`):
#########################################################################


###########################################################################
#Evaluates the matrix element <slaja||U^(k)||slbjb>.                      #
###########################################################################

uklj:=proc(la,lb,ja,jb,k)
      local result;
result:=combine(((-1)^(1/2+lb+ja+k))*sqrt((2*ja+1)*(2*jb+1))
*sixj(ja,k,jb,lb,1/2,la));
end:

###########################################################################
#Evaluates the matrix element <j^2Ja||U^(k)||j^2Jb>.                      #
###########################################################################

ukjj:=proc(j,Ja,Jb,k,l)
      local result,aa;
  result:= (combine(2*((-1)^(Ja+k+1))*sqrt((2*Ja+1)*(2*Jb+1))*uklj(l,l,j,j,k)
*sixj(Ja,k,Jb,j,j,j))*aa);
end:

###########################################################################
#Evaluates the matrix element <ja^2Ja||U^(k)||jajbJb>        .            #
###########################################################################
ukjd:=proc(ja,jb,Ja,Jb,k,l)
      local result,ab;
result:=(combine(sqrt(2*(2*Ja+1)*(2*Jb+1))*((-1)^(ja+jb+Ja+k))
*uklj(l,l,ja,jb,k)*sixj(Ja,k,Jb,jb,ja,ja))*ab);
end:
###########################################################################
#Evaluates the matrix element <jajbJa||U^(k)||jajbJb>                     #    
###########################################################################

ukjab:=proc(ja,jb,Ja,Jb,k,l)
      local result,aa,bb;
result:=(combine(((-1)^(ja+jb+Jb+k))*sqrt((2*Ja+1)*(2*Jb+1))*uklj(l,l,ja,ja,k)
*sixj(Ja,k,Jb,ja,jb,ja))*aa) 
+ (combine(((-1)^(ja+jb+Ja+k))*sqrt((2*Ja+1)*(2*Jb+1))*uklj(l,l,jb,jb,k)
*sixj(Ja,k,Jb,jb,ja,jb))*bb);
end:
###########################################################################
#Evaluates the matrix element <sljm|U_q^(k)|slj'm'>                       #    
###########################################################################
ukq:=proc(ja,ma,jb,mb,k,q)
     local result;
result:=(combine(((-1)^(ja-ma))*threej(ja,k,jb,-ma,q,mb)*((-1)^(1/2+3+ja+k))
*sqrt((2*ja+1)*(2*jb+1))*sixj(ja,k,jb,3,1/2,3)));
end:
###########################################################################
#Evaluates the matrix element <sljm|V|slj'm'>                             #    
###########################################################################
vcryst:=proc(ja,ma,jb,mb)
        local result,R2,R4,R6,R6q;
result:= combine(ck(3,3,2)*ukq(ja,ma,jb,mb,2,0)*R2(ja,jb)
+ck(3,3,4)*ukq(ja,ma,jb,mb,4,0)*R4(ja,jb)
+ck(3,3,6)*ukq(ja,ma,jb,mb,6,0)*R6(ja,jb)
+ck(3,3,6)*(ukq(ja,ma,jb,mb,6,6)+ukq(ja,ma,jb,mb,6,-6)
*R6q(ja,jb)));
end:         

###########################################################################
#Evaluates the matrix element <sljm|W_Q^(ks,kl)K|slj'm'>                  #    
###########################################################################
WKQ:=proc(ja,ma,jb,mb,ks,kl,K,Q)
     local result;
result:=combine(((-1)^(ja-ma))*threej(ja,K,jb,-ma,Q,mb)*
sqrt((2*ja+1)*(2*K+1)*(2*jb+1)*(2*ks+1)*(2*kl+1))
*ninej(1/2,1/2,ks,3,3,kl,ja,jb,K));
end:
###########################################################################
#Sets the linear combinations of radial integrals for the b_k(ks,kl)      #    
###########################################################################
bkk:=proc(k,i,j)
     local result,R2uu,R2ud,R2dd,R4uu,R4ud,R4dd,R6uu,R6ud;
if (k=2) then
if ((i=1) and (j=1))
then result:=4*sqrt(21)*(-5*R2uu + 3*R2ud + 2*R2dd)/245
elif ((i=1) and (j=3))
then result:=4*sqrt(7)*(5*R2uu + 4*R2ud - 9*R2dd)/245
elif ((i=0) and (j=2))
then result:=-2*sqrt(42)*(25*R2uu + 6*R2ud + 18*R2dd)/735
fi;
fi;
if (k=4) then
if ((i=1) and (j=3))
then result:=4*sqrt(21)*(-5*R4uu + 3*R4ud - R4dd)/441
elif ((i=1) and (j=5))
then result:=2*sqrt(2310)*(-3*R4uu - 8*R4ud + 11*R4dd)/4851
elif ((i=0) and (j=4))
then result:=2*sqrt(77)*(18*R4uu + 20*R4ud + 11*R4dd)/1617
fi;
fi;
if (k=6) then
if ((i=1) and (j=5))
then result:=20*sqrt(77)*(-R6uu + R6ud)/1001
elif ((i=0) and (j=6))
then result:=-10*sqrt(462)*(R6uu + 6*R6ud)/3003
fi;
fi; 
result;;
end:
###########################################################################
#Evaluates the matriR2uu element <sljm|bk(ij)W_Q^(ks,kl)K|slj'm'>            #    
###########################################################################
bWKQ:=proc(ja,ma,jb,mb,i,j,K,Q)
      local result,R2uu,R2ud,R2dd,R4uu,R4ud,R4dd,R6uu,R6ud;
result:=bkk(K,i,j)*WKQ(ja,ma,jb,mb,i,j,K,Q);
end:
###########################################################################
#Evaluates the matriR2uu element <sljm|\sum_bk(ij)W_Q^(ks,kl)K|slj'm'>       #    
###########################################################################
rkck2:=proc(ja,ma,jb,mb)
       local result,R2uu,R2ud,R2dd ;
result:=(bWKQ(ja,ma,jb,mb,1,1,2,0)+bWKQ(ja,ma,jb,mb,1,3,2,0)
+bWKQ(ja,ma,jb,mb,0,2,2,0));
end:
;