Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
396 lines (373 sloc) 10.8 KB
function [w6j,L,norma]=wigner6j(l2,l3,l4,l5,l6)
% [w6j,L,norma]=WIGNER6J(l2,l3,l4,l5,l6)
%
% Calculates Wigner 6j symbols by recursion, for all values up to L
% that are allowed in the expression (L l2 l3)
% (l4 l5 l6)
% There is no truncation at any bandwidth - they are all returned.
%
% INPUT:
%
% l2, l3, l4, l5, l6 The quantities defined above, all scalar
%
% OUTPUT:
%
% w6j The Wigner 6j symbols (with prepended zeroes)
% L The degrees at which these are evaluated
% norma Normalization check which passed muster
%
% EXAMPLES:
%
% wigner6j('demo1') % Just tests some values against Stone's code
% wigner6j('demo2') % Reproduces Table IV and Figure 3 of Schulten and Gordon
% wigner6j('demo3') % Reproduces Table III of Schulten and Gordon
% wigner6j('demo4') % Reproduces Figure 4 of Schulten and Gordon
%
% COMPARE WITH:
%
% Fortran code WIGNER6J.F based on original by Schulten and Gordon
% g77 -o $PFILES/wigner6j $FFILES/wigner6j.f $FFILES/drc6j.f ; wigner6j
% sprintf('%23.16e\n',wigner6j(0,8,7,6.5,7.5,7.5))
%
% Last modified by fjsimons-at-alum.mit.edu, 1/11/2007
if ~isstr(l2)
% Overflow parameters
%huge=sqrt(1.79D+308/20);
huge=sqrt(realmax/20);
sqhuge=sqrt(huge);
tiny=1/huge;
srtiny=1/sqhuge;
% Limits for L1
l1min=max(abs(l2-l3),abs(l5-l6));
l1max=min(l2+l3,l5+l6);
% Check error condition 4.
if mod(l1max-l1min,1)
error('l1max-l1min not integer')
end
% disp(sprintf('Computation from %i to %i',l1min,l1max))
% Check error conditions 1, 2, and 3.
if mod(l2+l3+l5+l6,1) || mod(l4+l2+l6,1)
% warning('l2+l3+l5+l6 or l4+l2+l6 not integer')
w6j=repmat(0,l1max+1,1)';
L=0:l1max; return
elseif l4+l2-l6<0 || l4-l2+l6<0 || -l4+l2+l6<0
% warning('l2, l4, l6 triangular condition not satisfied')
w6j=repmat(0,l1max+1,1)';
L=0:l1max; return
elseif l4-l5+l3<0 || l4+l5-l3<0 || -l4+l5+l3<0
% warning('l3, l4, l5 triangular condition not satisfied')
w6j=repmat(0,l1max+1,1)';
L=0:l1max; return
end
% Compute single coefficient directly by analytic formula
if l1min<l1max
nfin=floor(l1max-l1min+1);
elseif l1min==l1max
w6j=repmat(0,l1max+1,1)';
% This is reached in case that L1 can take only one value
w6j(l1max+1)=(-1)^floor(l2+l3+l5+l6)/sqrt((2*l1min+1)*(2*l4+1));
% Check normalization
norma=(2*l1min+1)*(2*l4+1)*w6j(l1max+1)^2;
difer(norma-1)
L=0:l1max;
return
else
% Check error condition 5
error('l1min greater than l1max')
end
% Start of forward recursion
l1=l1min;
newfac=0;
c1=0;
w6j(1)=srtiny;
sum1=(2*l1+1)*tiny;
lstep=1;
while lstep~=nfin % WHILE LOOP 30
lstep=lstep+1;
l1=l1+1;
oldfac=newfac;
a1=(l1+l2+l3+1)*(l1-l2+l3)*(l1+l2-l3)*(-l1+l2+l3+1);
a2=(l1+l5+l6+1)*(l1-l5+l6)*(l1+l5-l6)*(-l1+l5+l6+1);
newfac=sqrt(a1*a2);
if l1>1
dv=2*(l2*(l2+1)*l5*(l5+1)+l3*(l3+1)*l6*(l6+1)-l1*(l1-1)*l4*(l4+1))...
-(l2*(l2+1)+l3*(l3+1)-l1*(l1-1))*(l5*(l5+1)+l6*(l6+1)-l1*(l1-1));
denom=(l1-1)*newfac;
if ~(lstep-2<=0)
c1old=abs(c1);
end
c1=-(2*l1-1)*dv/denom;
else
% If L1=1, (L1-1) has to be factored out of DV, hence
c1=-2*(l2*(l2+1)+l5*(l5+1)-l4*(l4+1))/newfac;
end
if ~(lstep>2)
% If L1=L1MIN+1, the third term in recursion equation vanishes
x=srtiny*c1;
w6j(2)=x;
sum1=sum1+tiny*(2*l1+1)*c1^2;
if lstep==nfin
sumuni=sum1;
% Normalize 6j coefficients
cnorm=1/sqrt((2*l4+1)*sumuni);
% Sign convention for last 6j coefficient determines overall phase
sign1=sign(w6j(nfin));
sign2=(-1)^floor(l2+l3+l5+l6);
if sign1*sign2<=0
cnorm=-cnorm;
end
if abs(cnorm)<1
thresh=tiny/abs(cnorm);
for n=1:nfin
if abs(w6j(n))<thresh
w6j(n)=0;
end
w6j(n)=cnorm*w6j(n);
end
else
for n=1:nfin
w6j(n)=cnorm*w6j(n);
end
end
% Check normalization
norma=sum((2*l4+1)*(2*(l1min:l1max)+1).*w6j.^2);
difer(norma-1)
% Provide output
w6j=[repmat(0,1,l1min) w6j];
L=0:l1max;
return
end
else
c2=-l1*oldfac/denom;
% Recursion to the next 6j coefficient X
x=c1*w6j(lstep-1)+c2*w6j(lstep-2);
w6j(lstep)=x;
sumfor=sum1;
sum1=sum1+(2*l1+1)*x*x;
if lstep==nfin
break
end
if abs(x)<sqhuge
% As long as the coefficient ABS(C1) is decreasing, the recursion
% proceeds towards increasing 6j values and, hence, is numerically
% stable. Once an increase of ABS(C1) is detected, the recursion
% direction is reversed.
if c1old-abs(c1)<=0
break
end
else
% % This is reached if last 6j coefficient larger than sqhuge,
% % so that the recursion series W6J(1), ... ,W6J(LSTEP)
% % has to be rescaled to prevent overflow
for i=1:lstep
if abs(w6j(i))<srtiny
w6j(i)=0;
end
w6j(i)=w6j(i)/sqhuge;
end
sum1=sum1/huge;
sumfor=sumfor/huge;
x=x/sqhuge;
end
end
end % WHILE LOOP 30
% WHILE BREAK 100
% % Keep three 6j coefficients around LMATCH for comparison later
% % with backward recursion.
x1=x;
x2=w6j(lstep-1);
x3=w6j(lstep-2);
% % Starting backward recursion from L1MAX taking NSTEP2 steps, so
% % that forward and backward recursion overlap at the three points
nfinp1=nfin+1;
nfinp2=nfin+2;
nfinp3=nfin+3;
nstep2=nfin-lstep+3;
l1=l1max;
w6j(nfin)=srtiny;
sum2=(2*l1+1)*tiny;
l1=l1+2;
% WHILE LOOP LABEL 110
lstep=1;
while lstep~=nstep2
lstep=lstep+1;
l1=l1-1;
oldfac=newfac;
a1s=(l1+l2+l3)*(l1-l2+l3-1)*(l1+l2-l3-1)*(-l1+l2+l3+2);
a2s=(l1+l5+l6)*(l1-l5+l6-1)*(l1+l5-l6-1)*(-l1+l5+l6+2);
newfac=sqrt(a1s*a2s);
dv=2*(l2*(l2+1)*l5*(l5+1)+l3*(l3+1)*l6*(l6+1)-l1*(l1-1)*l4*(l4+1))...
-(l2*(l2+1)+l3*(l3+1)-l1*(l1-1))*(l5*(l5+1)+l6*(l6+1)-l1*(l1-1));
denom=l1*newfac;
c1=-(2*l1-1)*dv/denom;
if ~(lstep>2)
y=srtiny*c1;
w6j(nfin-1)=y;
if lstep==nstep2
break
end
sumbac=sum2;
sum2=sum2+(2*l1-3)*c1^2*tiny;
else
c2=-(l1-1)*oldfac/denom;
y=c1*w6j(nfinp2-lstep)+c2*w6j(nfinp3-lstep);
if lstep==nstep2
break
else
w6j(nfinp1-lstep)=y;
sumbac=sum2;
sum2=sum2+(2*l1-3)*y*y;
if ~(abs(y)<sqhuge)
for i=1:lstep
index=nfin-i+1
if abs(w6j(index))<srtiny
w6j(index)=0;
end
w6j(index)=w6j(index)/sqhuge;
end
sumbac=sumbac/huge;
sum2=sum2/huge;
end
end
end
end
% WHILE BREAK 200
y3=y;
y2=w6j(nfinp2-lstep);
y1=w6j(nfinp3-lstep);
% Determine now RATIO such that YI=RATIO*XI (I=1,2,3) holds
% with minimal error.
warning off
ratio=(x1*y1+x2*y2+x3*y3)/(x1*x1+x2*x2+x3*x3);
warning on
nlim=nfin-nstep2+1;
if abs(ratio)<1
nlim=nlim+1;
ratio=1/ratio;
for n=nlim:nfin
w6j(n)=ratio*w6j(n);
end
sumuni=sumfor+ratio^2*sumbac;
else
for n=1:nlim
w6j(n)=ratio*w6j(n);
end
sumuni=ratio^2*sumfor+sumbac;
end
% Normalize 6j coefficients
cnorm=1/sqrt((2*l4+1)*sumuni);
% Sign convention for last 6j coefficient determines overall phase
sign1=sign(w6j(nfin));
sign2=(-1)^floor(l2+l3+l5+l6);
if sign1*sign2<=0
cnorm=-cnorm;
end
if abs(cnorm)<1
thresh=tiny/abs(cnorm);
for n=1:nfin
if abs(w6j(n))<thresh
w6j(n)=0;
end
w6j(n)=cnorm*w6j(n);
end
else
for n=1:nfin
w6j(n)=cnorm*w6j(n);
end
end
% Check normalization
norma=sum((2*l4+1)*(2*(l1min:l1max)+1).*w6j.^2);
difer(norma-1)
% Provide output:
% Provide the right number of zeroes
w6j=[repmat(0,1,l1min) w6j];
L=0:l1max;
elseif strcmp('demo1',l2)
% Useing Anthony Stones' RFF calculator for most of these
difer(indeks(wigner6j(1,10,0,10,1),10:12)-...
[sqrt(1/7)/3 -sqrt(1/7)/3 sqrt(1/7)/3])
difer(wigner6j(2,2,3,2,3)-[0 sqrt(3/7)/5 -sqrt(3/2)/35 -1/10 -sqrt(1/3)/14])
difer(wigner6j(1,1,2,1,2)-[0 -sqrt(1/5)/2 -1/10])
difer(wigner6j(1,2,1,2,2)-[0 -1/10 sqrt(7/3)/10 -sqrt(2/3)/5])
difer(wigner6j(1,1,1,0,0)-sqrt(1/3))
difer(indeks(wigner6j(0,1,1,0,1),2)-1/3)
difer(indeks(wigner6j(2,1,1,0,1),2)-1/3)
difer(indeks(wigner6j(0,1,1,0,1),2)-1/3)
difer(wigner6j(2,2,2,0,0)-sqrt(1/5))
difer(indeks(wigner6j(3,2,2,0,1),2)-sqrt(1/15))
difer(indeks(wigner6j(4,2,2,0,2),3)-1/5)
difer(indeks(wigner6j(1,3,3,0,2),3)-sqrt(1/35))
difer(indeks(wigner6j(3,3,3,0,2),3)-sqrt(1/35))
difer(indeks(wigner6j(5,3,3,0,2),3)-sqrt(1/35))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),2)-(7/680)*sqrt(23/2))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),3)--(13/680)*sqrt(23/6))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),4)-(1147/167960)*sqrt(23/3))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),5)-(23/33592)*sqrt(115))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),6)--(181/33592)*sqrt(115/6))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),7)-(4343/167960)*sqrt(23/42))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),8)-(1941/83980)*sqrt(1/322))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),9)--(10993/83980)*sqrt(1/46))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),10)-(1709/33592)*sqrt(5/46))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),11)-(13/2584)*sqrt(55/46))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),12)--(691/12920)*sqrt(11/69))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),13)-(103/12920)*sqrt(13/69))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),14)-(1621/12920)*sqrt(13/322))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),15)-(533/2584)*sqrt(5/966))
difer(indeks(wigner6j(8,7,13/2,15/2,15/2),16)-(13/1292)*sqrt(5/69))
% An example from Rotenberg
difer(indeks(wigner6j(8,8,8,8,8),9)-...
rotenberg([-2 2 -2 0 0 -2 -2 -2 -2],-1)*4073)
elseif strcmp('demo2',l2)
[w,L]=wigner6j(48,80,112,120,72);
disp(sprintf('------------------------------'))
disp(sprintf('L1 Values of 6j coefficients'))
disp(sprintf('------------------------------'))
disp(sprintf('%2.2i %23.16i\n',[L(1:8:end)' w(1:8:end)']'))
disp(sprintf('------------------------------'))
clf
w(1:5)=w(1:5)*1e3;
w=[w(1:5) NaN w(6:end)];
L=[L(1:5) NaN L(6:end)];
plot(L,w*1e4,'x-')
ylim([-12 16.5])
set(gca,'ytick',-12:2:16)
xlim([45 130])
shrink(gca,2,2)
grid on
figdisp('Schulten+75-fig3',[],[],1)
elseif strcmp('demo3',l2)
[w,L]=wigner6j(8,7,13/2,15/2,15/2);
disp(sprintf('------------------------------'))
disp(sprintf('L1 Values of 6j coefficients'))
disp(sprintf('------------------------------'))
disp(sprintf('%2.2i %23.16i\n',[L(1:9)' w(1:9)']'))
disp(sprintf('------------------------------'))
plot(L,w*1e4,'+-')
elseif strcmp('demo4',l2)
clf
subplot(311)
[w,L]=wigner6j(48,80,120,120,72);
plot(L,w*1e4,'x-')
ylim([0 16.5])
axis tight
set(gca,'ytick',0:2:16)
shrink(gca,2,1)
grid on
subplot(312)
[w,L]=wigner6j(48,80,115,120,72);
plot(L,w*1e4,'x-')
ylim([-12 16.5])
axis tight
set(gca,'ytick',-12:2:16)
shrink(gca,2,1)
grid on
subplot(313)
[w,L]=wigner6j(48,80,110,120,72);
plot(L,w*1e4,'x-')
ylim([-12 14])
axis tight
set(gca,'ytick',-12:2:14)
shrink(gca,2,1)
grid on
figdisp('Schulten+75-fig4',[],[],1)
end
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.