Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
426 lines (373 sloc) 12.9 KB
function [w3j,j]=wigner3jm(L,l2,l3,m1,m2,m3)
% [w3j,j]=WIGNER3JM(L,l2,l3,m1,m2,m3)
%
% Calculates Wigner 3j symbols by recursion, for all values of j<=L
% allowed in the expression (L l2 l3)
% (m1 m2 m3)
% There is no truncation at any bandwidth - they are all returned
% Note the selection rules:
% jmin = max(|l2-l3|, |m1|)
% jmax = l2 + l3
% m1 + m2 + m3 = 0.
%
% INPUT:
%
% L Maximum degree, bandwidth [default: that what's allowed]
% l2,l3 Other degrees in the symbol above, don't have to be integers
% m1,m2,m3 Orders in the symbol above, don't have to be integers, and
% note that the defaults here are not zero!
%
% OUTPUT:
%
% w3j The Wigner3j function
% j The first degrees from 0 to L
% With normalization check
%
% EXAMPLES:
%
% wigner3jm('demo1') % Should return nothing if it all works
% wigner3jm('demo2') % Reproduces Table I of Schulten and Gordon
%
% See also: WIGNER0J, GUSEINOV, THREEJ, ZEROJ
%
% Last modified by fjsimons-at-alum.mit.edu, 06/15/2010
% Note, if you truncate, like here, at the degree L, you're actually
% doing too much work; you could improve this.
% After Fortran (could you tell?) by Mark Wieczorek, who further notes:
% Returned values have a relative error less than ~1.d-8 when l2 and l3 are
% less than 103 (see below). In practice, this routine is probably usable up
% to 165. This routine is based upon the stable non-linear recurrence
% relations of Luscombe and Luban (1998) for the "non classical" regions near
% jmin and jmax. For the classical region, the standard three term recursion
% relationship is used (Schulten and Gordon 1975). Note that this three term
% recursion can be unstable and can also lead to overflows. Thus the values
% are rescaled by a factor "scalef" whenever the absolute value of the 3j
% coefficient becomes greater than unity. Also, the direction of the iteration
% starts from low values of j to high values, but when abs(w3j(j+2)/w3j(j)) is
% less than one, the iteration will restart from high to low values. More
% efficient algorithms might be found for specific cases (for instance, when
% all m's are zero).
% Verification:
% The results have been verified against this routine run in quadruple
% precision. For 1.e7 acceptable random values of l2, l3, m2, and m3 between
% -200 and 200, the relative error was calculated only for those 3j
% coefficients that had an absolute value greater than 1.d-17 (values smaller
% than this are for all practical purposed zero, and can be heavily affected
% by machine roundoff errors or underflow). 853 combinations of parameters
% were found to have relative errors greater than 1.d-8. Here I list the
% minimum value of max(l2,l3) for different ranges of error, as well as the
% number of times this occured 1.d-7 < error <=1.d-8 = 103 # = 483 1.d-6 <
% error <= 1.d-7 = 116 # = 240 1.d-5 < error <= 1.d-6 = 165 # = 93 1.d-4 <
% error <= 1.d-5 = 167 # = 36
% Many times (maybe always), the large relative errors occur when the 3j
% coefficient changes sign and is close to zero. (I.e., adjacent values are
% about 10.e7 times greater in magnitude.) Thus, if one does not need to know
% highly accurate values of the 3j coefficients when they are almost zero
% (i.e., ~1.d-10) then this routine is probably usable up to about 160.
% These results have also been verified for parameter values less than 100
% using a code based on the algorithm of de Blanc (1987), which was originally
% coded by Olav van Genabeek, and modified by M. Fang (note that this code was
% run in quadruple precision, and only calculates one coefficient for each
% call. I also have no idea if this code was verified.) Maximum relative
% errors in this case were less than 1.d-8 for a large number of values
% (again, only 3j coefficients greater than 1.d-17 were considered here). The
% biggest improvement that could be made in this routine is to determine when
% one should stop iterating in the forward direction, and start iterating from
% high to low values.
if ~isstr(L)
defval('l2',6)
defval('l3',5)
defval('L',l2+l3)
defval('m1',3)
defval('m2',2)
defval('m3',-5)
flag1=0;
flag2=0;
% Didn't realize this factor was wrong until 10/02/2006
% But - Luscombe and Luban imply that the three-term recursion in the
% classical, oscillatory region rarely suffers from the overflows - all
% should be probably well
scalef=1000;
jmin=max(abs(l2-l3),abs(m1));
jmax=l2+l3;
jnum=jmax-jmin+1;
% Initialize
w3j=repmat(0,1,jnum);
if abs(m2)>l2 | abs(m3)>l3
w3j=repmat(0,L+1,1)';
j=0:L; return
elseif m1+m2+m3~=0
w3j=repmat(0,L+1,1)';
j=0:L; return
elseif jmax<jmin
w3j=repmat(0,L+1,1)';
j=0:L; return
end
% Only one value is allowed
if jnum==1
w3j=1/sqrt(2*jmin+1);
if (w3j<0 && (-1)^(l2-l3+m2+m3)>0) || ...
(w3j>0 && (-1)^(l2-l3+m2+m3)<0)
w3j=-w3j;
end
[w3j,j]=output(jmin,l2,l3,L,w3j); return
end
% Calculate lower non-classical values for [jmin, jn]. If the second
% term can not be calculated because the recursion relationsips give
% rise to a 1/0, then set flag1 to 1. If all m's are zero, then this
% is not a problem as all odd terms must be zero.
[rs,wu,wl]=deal(0);
warning off
rs(1)=-x(jmin,l2,l3,m1)/y(jmin,l2,l3,m1,m2,m3);
warning on
if m1==0 && m2==0 && m3==0
wl(1)=1;
wl(2)=0;
jn=jmin+1;
elseif y(jmin,l2,l3,m1,m2,m3)==0
if x(jmin,l2,l3,m1)==0
flag1=1;
jn=jmin;
else
wl(1)=1;
wl(2)=0;
jn=jmin+1;
end
elseif rs(1)<=0
wl(1)=1;
wl(2)=-y(jmin,l2,l3,m1,m2,m3)/x(jmin,l2,l3,m1);
jn=jmin+1;
else
jn=jmax;
for j=jmin+1:jmax-1
denom=y(j,l2,l3,m1,m2,m3)+z(j,l2,l3,m1)*rs(j-jmin);
warning off
rs(j-jmin+1)=-x(j,l2,l3,m1)/denom;
warning on
if (rs(j-jmin+1)>1| rs(j-jmin+1) <= 0| denom==0)
jn=j-1;
break
end
end
wl(jn-jmin+1)=1;
for k=1:jn-jmin
wl(jn-k-jmin+1)=wl(jn-k-jmin+2)*rs(jn-k-jmin+1);
end
if (jn==jmin)
wl(2)=-y(jmin,l2,l3,m1,m2,m3)/x(jmin,l2,l3,m1);
jn=jmin+1;
end
end
if jn==jmax
w3j(1:jnum)=wl(1:jnum);
w3j=normw3j(w3j,jmin,jmax,jnum);
w3j=fixsign(w3j,jnum,l2,l3,m2,m3);
[w3j,j]=output(jmin,l2,l3,L,w3j); return
end
% Calculate upper non-classical values for [jp, jmax].
% If the second last term can not be calculated because the
% recursion relations give a 1/0, then set flag2 to 1.
warning off
rs(jnum)=-z(jmax,l2,l3,m1)/y(jmax,l2,l3,m1,m2,m3);
warning on
if m1==0 && m2==0 && m3==0
wu(jnum)=1;
wu(jmax-jmin)=0;
jp=jmax-1;
elseif y(jmax,l2,l3,m1,m2,m3)==0
if z(jmax,l2,l3,m1)==0
flag2=1;
jp=jmax;
else
wu(jnum)=1;
wu(jmax-jmin)=-y(jmax,l2,l3,m1,m2,m3)/z(jmax,l2,l3,m1);
jp=jmax-1;
end
elseif rs(jnum)<=0
wu(jnum)=1;
wu(jmax-jmin)=-y(jmax,l2,l3,m1,m2,m3)/z(jmax,l2,l3,m1);
jp=jmax-1;
else
jp=jmin;
for j=jmax-1:-1:jn
% This appears to be Luscombe and Luban's Eq. (2)
denom=y(j,l2,l3,m1,m2,m3)+x(j,l2,l3,m1)*rs(j-jmin+2);
warning off
rs(j-jmin+1)=-z(j,l2,l3,m1)/denom;
warning on
if (rs(j-jmin+1)>1| rs(j-jmin+1) <= 0| denom==0)
jp=j+1;
break
end
end
wu(jp-jmin+1)=1;
for k=1:jmax-jp
wu(jp+k-jmin+1)=wu(jp+k-jmin)*rs(jp+k-jmin+1);
end
if jp==jmax
wu(jmax-jmin)=-y(jmax,l2,l3,m1,m2,m3)/z(jmax,l2,l3,m1);
jp=jmax-1;
end
end
% Calculate classical terms for [jn+1, jp-1] using standard three term
% rercusion relationship. Start from both jn and jp and stop at the
% midpoint. If flag1 is set, then perform the recursion solely from
% high to low values. If flag2 is set, then perform the recursion
% solely from low to high.
if flag1==0
% I think Fortran rounds like this
jmid=round((jn+jp)/2);
for j=jn:jmid-1
wl(j-jmin+2)=-(z(j,l2,l3,m1)*wl(j-jmin)+...
y(j,l2,l3,m1,m2,m3)* ...
wl(j-jmin+1))/x(j,l2,l3,m1);
if abs(wl(j-jmin+2))>1
wl(1:j-jmin+2)=...
wl(1:j-jmin+2)/scalef;
end
if abs(wl(j-jmin+2)/wl(j-jmin))<1 & ...
wl(j-jmin+2)~=0
jmid=j+1;
break
end
end
wnmid=wl(jmid-jmin+1);
warning off
if abs(wnmid/wl(jmid-jmin))<1e-6 &wl(jmid-jmin)~=0
wnmid=wl(jmid-jmin);
jmid=jmid-1;
end
warning on
for j=jp:-1:jmid+1
wu(j-jmin)=-(x(j,l2,l3,m1)*wu(j-jmin+2)+...
y(j,l2,l3,m1,m2,m3)* ...
wu(j-jmin+1))/z(j,l2,l3,m1);
if abs(wu(j-jmin))>1
wu(j-jmin:jnum)=...
wu(j-jmin:jnum)/scalef;
end
end
wpmid=wu(jmid-jmin+1);
if jmid==jmax
w3j(1:jnum)=wl(1:jnum);
elseif jmid==jmin
w3j(1:jnum)=wu(1:jnum);
else
w3j(1:jmid-jmin+1)=wl(1:jmid-jmin+1)*wpmid/wnmid;
w3j(jmid-jmin+2:jnum)=...
wu(jmid-jmin+2:jnum);
end
elseif flag1==1 & flag2==0
for j=jp:-1:jmin+1
wu(j-jmin)=-(x(j,l2,l3,m1)*wu(j-jmin+2)+...
y(j,l2,l3,m1,m2,m3)* ...
wu(j-jmin+1))/z(j,l2,l3,m1);
if abs(wu(j-jmin))>1
wu(j-jmin:jnum)=...
wu(j-jmin:jnum)/scalef;
end
end
w3j(1:jnum)=wu(1:jnum);
elseif flag2==1 & flag1==0
for j=jn:jp-1
wl(j-jmin+2)=-(z(j,l2,l3,m1)*wl(j-jmin)+...
y(j,l2,l3,m1,m2,m3)* ...
wl(j-jmin+1))/x(j,l2,l3,m1);
if abs(wl(j-jmin+2))>1
wl(1:j-jmin+2)=...
wl(1:j-jmin+2)/scalef;
end
end
w3j(1:jnum)=wl(1:jnum);
elseif flag1==1 & flag2==1
error('Can not calculate function for input values')
end
w3j=normw3j(w3j,jmin,jmax,jnum);
w3j=fixsign(w3j,jnum,l2,l3,m2,m3);
% Output: give output in all degrees from zero to the bandwidth
if L<=l2+l3
% Truncate
w3j=[repmat(0,1,jmin) w3j(1:L-jmin+1)];
else
% Append
w3j=[repmat(0,1,jmin) w3j repmat(0,1,L-l2-l3)];
end
j=0:L;
% Perform normalization check
if L==l2+l3
norma=sum((2*[0:L]+1).*w3j.^2);
difer(norma-1,[],[],NaN)
end
elseif strcmp(L,'demo1')
difer(wigner3jm(20,10,10,0,0,0)-wigner0j(20,10,10))
difer(wigner3jm(400,200,200,0,0,0)-wigner0j(400,200,200))
difer(wigner3jm(40,10,13,0,0,0)-wigner0j(40,10,13))
difer(indeks(wigner3jm(8,6,5,3,2,-5),'end')-sqrt(42/4199))
difer(indeks(wigner3jm(8,6,5,3,2,-5),'end-1')-35*sqrt(2/138567))
difer(indeks(wigner3jm(8,6,5,3,2,-5),'end-2')-7*sqrt(1/2431))
difer(indeks(wigner3jm(8,6,5,3,2,-5),'end-3')-sqrt(35/2431))
difer(indeks(wigner3jm(8,6,5,3,2,-5),'end-4')-sqrt(5/858))
difer(indeks(wigner3jm(3,2,2,1,0,-1),'end')--sqrt(1/35))
difer(indeks(wigner3jm(0,1,1,0,1,-1),'end')-sqrt(1/3))
difer(indeks(wigner3jm(1,1,1,0,1,-1),'end')-sqrt(1/6))
difer(indeks(wigner3jm(2,1,1,0,1,-1),'end')-sqrt(1/30))
difer(indeks(wigner3jm(3,2,3,-2,0,2),'end')-0)
difer(indeks(wigner3jm(6,2,4,-4,1,3),'end')-4*sqrt(1/429))
difer(indeks(wigner3jm(0,1,1,0,0,0),'end')--sqrt(1/3))
difer(indeks(wigner3jm(2,1,1,0,0,0),'end')-sqrt(2/15))
difer(indeks(wigner3jm(3,2,3,3,0,-3),'end')-(1/2)*sqrt(5/21))
difer(indeks(wigner3jm(2,2,3,-2,0,2),'end')--sqrt(1/14))
difer(indeks(wigner3jm(8,6,5,3,2,-5),'end-5')-sqrt(1/1001))
difer(indeks(wigner3jm(20,15,9,-3,2,1),'end')+(311/2)*sqrt(115/1231496049))
difer(indeks(wigner3jm(10,10,12,9,3,-12),'end')--(1/15)*sqrt(4199/9889))
elseif strcmp(L,'demo2')
[w,L]=wigner3jm([],9/2,7/2,1,-7/2,5/2);
disp(sprintf('------------------------------'))
disp(sprintf('L1 Values of 3j coefficients'))
disp(sprintf('------------------------------'))
disp(sprintf('%2.2i %23.16i\n',[L(2:9)' w(2:9)']'))
disp(sprintf('------------------------------'))
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [o,p]=output(jmin,l2,l3,L,w3j)
% Output: give output in all degrees from zero to the bandwidth
% Perform normalization check before possible completion or truncation
norma=sum((2*[jmin:l2+l3]+1).*w3j.^2);
difer(norma-1,[],[],NaN)
if L<=l2+l3
% Truncate...
o=[repmat(0,1,jmin) w3j(1:L-jmin+1)];
else
% Append
o=[repmat(0,1,jmin) w3j repmat(0,1,L-l2-l3)];
end
p=0:L;
% The following functions are straight from Table I of Luscombe and Luban
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function o=x(j,l2,l3,m1)
o=j*a(j+1,l2,l3,m1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function o=y(j,l2,l3,m1,m2,m3)
% This is the Y-function in Table 1 of Luscombe and Luban
o=-(2*j+1)*(m1*(l2*(l2+1)-l3*(l3+1))-(m3-m2)*j*(j+1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function o=z(j,l2,l3,m1)
o=(j+1)*a(j,l2,l3,m1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function o=a(j,l2,l3,m1)
o=sqrt((j^2-(l2-l3)^2)*((l2+l3+1)^2-j^2)*(j^2-m1^2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function o=normw3j(w3j,jmin,jmax,jnum)
normo=0;
for j=jmin:jmax
normo=normo+(2*j+1)*w3j(j-jmin+1)^2;
end
o(1:jnum)=w3j(1:jnum)/sqrt(normo);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function o=fixsign(w3j,jnum,l2,l3,m2,m3)
if ((w3j(jnum)<0 & (-1)^(l2-l3+m2+m3)>0) |...
(w3j(jnum)>0 & (-1)^(l2-l3+m2+m3)<0))
o(1:jnum)=-w3j(1:jnum);
else
o(1:jnum)=w3j(1:jnum);
end