The following code (not very well optimized) was used to generate
Figure 3 of Y. Taylor, C. Woodward, "6j symbols for U_q(sl_2)
and non-Euclidean tetrahedra".

with(plots); with(linalg);
# quantum integer,factorial
pi:=evalf(Pi); nq:=(r,n)->sin(pi*n/r)/sin(pi/r);
faq:=(r,n)->evalf(product(nq(r,ii),ii=1..n),30);
# quantum triangle symbol, Racah sum, quantum sixj symbol
trq:=(r,a,b,c)->faq(r,(b+c-a)/2)*faq(r,(a+c-b)/2)*
faq(r,(a+b-c)/2)/(faq(r,(a+b+c)/2+1));
syq:=(r,a1,a2,a3,a4,b1,b2,b3)->sum((-1)^z*faq(r,z+1)/
(faq(r,b1-z)*faq(r,b2-z)*faq(r,b3-z)*faq(r,z-a1)*
faq(r,z-a2)*faq(r,z-a3)*faq(r,z-a4)),
z=max(a1,a2,a3,a4)..min(b1,b2,b3));
sjq:=(r,a,b,c,d,e,f)->`if`(min(a+b-c,a-b+c,-a+b+c,a+e-f,e-a+f,-e+a+f,
d+c-e,d-c+e,-d+c+e,b+f-d,b-f+d,-b+f+d,2*r-4-a-b-c,2*r-a-e-f,2*r-d-c-e,
2*r-b-d-f)>= 0,evalf((trq(r,a,b,c)*trq(r,a,e,f)*trq(r,d,c,e)*
trq(r,b,d,f))^(1/2)*syq(r,(a+b+c)/2,(a+e+f)/2,(c+d+e)/2,(b+d+f)/2,
(a+b+d+e)/2,(a+c+d+f)/2,(b+c+e+f)/2),30),0);
# length corresponding to a dominant weight (integer)
len:=j->pi*(j+1)/r;
# The Gram matrix, its determinant and inverse
G:=(i,j,k,l,m,n)->linalg[matrix](4,4,[1, cos(i),cos(k),cos(m),
cos(i),1, cos(j),cos(n),
cos(k),cos(j),1, cos(l),
cos(m),cos(n),cos(l),1]);
detG:=(i,j,k,l,m,n)->evalf(det(G(i,j,k,l,m,n)),10);
Ginv:=(i,j,k,l,m,n)->evalf(inverse(G(i,j,k,l,m,n)));
# Does the tetrahedron exist?
tetexist:=(r,a,b,c,d,e,f)->`if`(min(a+b-c+1,a-b+c+1,-a+b+c+1,a+e-f+1,
e-a+f+1,-e+a+f+1,d+c-e+1,d-c+e+1,-d+c+e+1,b+f-d+1,b-f+d+1,-b+f+d+1,
2*r-4-a-b-c-3,2*r-a-e-f-3,2*r-d-c-e-3,2*r-b-d-f-3,
detG(len(a),len(b),len(c),len(d),len(e),len(f))) >= 0,1,0);
# The amplitude of the asymptotic formula
amppredict:=(r,i,j,k,l,m,n)->evalf(((r/pi)^3*1.5*pi*(abs(detG(
len(i),len(j),len(k),len(l),len(m),len(n)))/(36))^(1/2))^(-1/2),10);
# The Dihedral angles
preangle:=(i,j,k,l,m,n)->evalf(seq(seq(pi -
arccos(-Ginv(i,j,k,l,m,n)[a,b]/(Ginv(i,j,k,l,m,n)[a,a]*
Ginv(i,j,k,l,m,n)[b,b])^(1/2)),b=a+1..4),a=1..4),10);
diangle:=(r,i,j,k,l,m,n)->evalf(preangle(len(i),len(j),len(k),len(l),
len(m),len(n)));
# The "Dehn invariant" aka Regge action
dehn:=(r,i,j,k,l,m,n)-> diangle(r,i,j,k,l,m,n)[1]*(l+1)/2+
diangle(r,i,j,k,l,m,n)[2]*(n+1)/2+ diangle(r,i,j,k,l,m,n)[3]*(j+1)/2+
diangle(r,i,j,k,l,m,n)[4]*(m+1)/2+ diangle(r,i,j,k,l,m,n)[5]*(k+1)/2+
diangle(r,i,j,k,l,m,n)[6]*(i+1)/2;
# The asymptotic formula for the degenerate case
G3:=(i,j,k)->linalg[matrix](3,3,[1,cos(i),cos(j),
cos(i),1,cos(k),
cos(j),cos(k),1]);
areaq:=(r,i,j,k)->evalf(det(G3(len(i),len(j),len(k))))^(1/2);
tang:=(r,i,j,k,l,m,n)->evalf(r^(-4/3)*2^(2/3)*3^(-1/3)*pi^(4/3)*
GAMMA(2/3)^(-1)*(areaq(r,i,j,k)*areaq(r,i,m,n)*areaq(r,l,j,n)*
areaq(r,l,m,k))^(-1/6));
# NUMERICAL EXPERIMENT
r:=200; i:=40; j:=48; k:=50; l:=52; m:=54; stepsize:=.1;
u:=2*max(i,j,k,l,m);
#calculate the 6j symbols
symbols:= [ seq([2*y,sjq(r,i,j,k,l,m,2*y)],y=0..u/2) ];
#does the terahedron exist
exist:= [seq([y*stepsize,tetexist(r,i,j,k,l,m,y*stepsize)],
y=1..(r/2)/stepsize)];
#the predictions for the tangent case, plus or minus
tapredict:= [ seq([y,`if`(tetexist(r,i,j,k,l,m,y)=1,
tang(r,i,j,k,l,m,y),undefined)],y=0..u/2) ];
mtapredict:= [ seq([y,-`if`(tetexist(r,i,j,k,l,m,y)=1,
tang(r,i,j,k,l,m,y),undefined)],y=(u/2)..u) ];
#the dihedral angles
angpredict := [ seq([y*stepsize,diangle(r,i,j,k,l,m,
y*stepsize)[2]],y=0..u/stepsize) ];
for y from 0 to u/stepsize do
phqpredict[y+1] := `if`(tetexist(r,i,j,k,l,m,y*stepsize)=1,
phqpredict[y] + angpredict[y+1][2]*stepsize/2,
`if`(tetexist(r,i,j,k,l,m,(y+1)*stepsize)=1,
dehn(r,i,j,k,l,m,(y+1)*stepsize),0)) od;
#the predicted amplitudes , plus or minus
predicta:=[ seq([y*stepsize, `if`(tetexist(r,i,j,k,l,m,y*stepsize)=1,
amppredict(r,i,j,k,l,m,y*stepsize),undefined)],y=1..u/stepsize) ];
mpredicta:=[ seq([y*stepsize,`if`(tetexist(r,i,j,k,l,m,y*stepsize)=1,
-amppredict(r,i,j,k,l,m,y*stepsize),undefined)],y=1..u/stepsize) ];
#the prediction for the non-degenerate case
predict:= [seq([stepsize*y,predicta[y][2]*cos(pi*(1/4) +
phqpredict[y])],y=1..u/stepsize)];
# Show all the plots
plot([symbols,predict,mtapredict,tapredict,predicta,mpredicta],
style=[point,line,line,line,line,line],
# color=[red,blue,green,green,yellow,yellow],
# linestyle=[SOLID,SOLID,DOT,DOT,SOLID,SOLID],
symbol=circle,symbolsize=5);
Back to my homepage .