Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
217 lines (201 sloc) 6.76 KB
function [s,C,S]=sixj(l1,l2,l3,l4,l5,l6,L,meth,C,S)
% [s,C,S]=SIXJ(l1,l2,l3,l4,l5,l6,L,meth,C,S)
%
% Wigner 6j-symbol from a database precomputed by WIGNERCYCLE.
%
% INPUT:
%
% l1,l2,l3 Top row of the Wigner 6j symbol [may be same-length vectors]
% l4,l5,l6 Bottom row of the Wigner 3j symbol [may be same-length vectors]
% L The bandwidth of the database [default: best available]
% meth 0 Uses sparse matrices [elegant, but slow]
% 1 Performs linear search on unsorted array [slow]
% 2 Performs binary search on presorted array [default]
% C,S The column and element vectors resulting from a previous load
%
% OUTPUT:
%
% s The (vector of) Wigner 6j symbols
% C,S The column and element vectors good for the next load
%
% EXAMPLE:
%
% sixj('demo1') % Should return nothing if it all works
%
% SEE ALSO: WIGNER6J, WIGNER3JM, GAUNT, WIGNER0J, GUSEINOV, ZEROJ
%
% Last modified by fjsimons-at-alum.mit.edu, 02/03/2007
% MUST want to build in conditions before even accessing database
% If NONE of the inputs satisfy the rules, that is
defval('xver',0)
if ~isstr(l1) % Not a demo
defval('C',[])
defval('S',[])
% Method
defval('meth',2)
% disp(sprintf('Using method %i',meth))
% All saved values must be integers
if sum(mod(l1,1)) || sum(mod(l2,1)) || sum(mod(l3,1)) || ...
sum(mod(l4,1)) || sum(mod(l5,1)) || sum(mod(l6,1))
error('All degrees and orders must be integers for the database')
end
if isempty(C) && isempty(S)
% Got something to do
% What is the lowest of the available database bandwidths?
Els=ls2cell(fullfile(getenv('IFILES'),'WIGNER','WIGNER6JCS-*-C'));
for index=1:length(Els)
EL(index)=str2num(rindeks(parse(Els{index},'-'),2));
end
EL=sort(EL);
% Bandwidth of the database; keep this as low as possible
% Need to provide for empties if not found
fmax=find(max([l1(:)' l2(:)' l3(:)' l4(:)' l5(:)' l6(:)'])<=EL);
if ~isempty(fmax)
defval('L',EL(indeks(fmax,1)))
else
defval('L',-1)
end
% Check, identify and load database
if any([l1(:)' l2(:)' l3(:)' l4(:)' l5(:)' l6(:)']>L)
error('Insufficient bandwidth for database')
end
% Check rules here
fnpl1=sprintf('%s/WIGNER6JCS-%i-C',...
fullfile(getenv('IFILES'),'WIGNER'),L);
fnpl2=sprintf('%s/WIGNER6JCS-%i-S',...
fullfile(getenv('IFILES'),'WIGNER'),L);
if exist(fnpl1,'file')==2 && exist(fnpl2,'file')==2
fmt1='uint64'; fmt2='float64';
disp(sprintf('Loading %s',fnpl1))
C=eval(sprintf('loadb(''%s'',''%s'')',fnpl1,fmt1));
%[C,j]=sort(C);
% Do clear and pack for the first sort - FOR large L ONLY
%writeb(C,fnpl1,fmt1);
%save j
%clear C
%disp(sprintf('Loading %s',fnpl2))
%S=eval(sprintf('loadb(''%s'',''%s'')',fnpl2,fmt2));
%S=S(j);
%writeb(S,fnpl2,fmt2); clear j
%disp('Temporary')
disp(sprintf('Loading %s',fnpl2))
S=eval(sprintf('loadb(''%s'',''%s'')',fnpl2,fmt2));
% Whatever happens, this better be sorted; check some entries
randin=unique(ceil(rand(min(100,length(C)),1)*length(C)));
if ~all(unique(C(randin))==C(randin))
disp('Column arrays were not properly sorted')
[C,j]=sort(C);
% Do clear and pack for the first sort
writeb(C,fnpl1,fmt1);
S=S(j);
writeb(S,fnpl2,fmt2); clear j
disp('Column arrays now sorted once and for all')
end
else
disp('Executing WIGNERCYCLE')
% Precompute the database
wignercycle(L,6);
% And have a go again
s=sixj(l1,l2,l3,l4,l5,l6,L,meth);
end
else
if isempty(L)
error('If supplying vectors with coefficients must also specify bandwidth')
end
% Else have C and S from a previous load and do nothing, but check
% There is NO check on whether the L supplied is appropriate for the
% C & S combination that is indeed supplied
if any([l1(:)' l2(:)' l3(:)']>L)
error('Insufficient bandwidth for database')
end
end
% Find running index into this matrix
% There must be the same number of all of the input degrees
for ix=1:length(l1)
rind(ix)=1+l1(ix)*(L+1)^0+...
l2(ix)*(L+1)^1+...
l3(ix)*(L+1)^2+...
l4(ix)*(L+1)^3+...
l5(ix)*(L+1)^4+...
l6(ix)*(L+1)^5;
end
% Initialize results vector
s=repmat(NaN,1,length(rind));
switch meth
case 0
% Turn it into the properly indexed sparse matrix
% This step takes some of time but most of all, memory
W=sparse(1,C,S,1,(L+1)^6);
% Extract the Wigner 6j-symbol
s=W(1,rind);
case 1
for index=1:length(rind)
posi=find(C==rind(index));
if ~isempty(posi)
s(index)=S(posi);
else
s(index)=0;
end
end
case 2
% Binary search algorithm on sorted arrays
for index=1:length(rind)
posi=binsearch(C,rind(index));
if ~isempty(posi)
s(index)=S(posi);
else
s(index)=0;
end
end
otherwise
error('Specify valid method')
end
% Check some special cases by Edmonds' rules
if xver==1 && length(l1)==1 && length(l2)==1 && length(l3)==1
es=l1+l2+l3;
ex=l2*(l2+1)+l3*(l3+1)-l1*(l1+1);
if l4==0 & s~=0
difer(s-(-1)^es/sqrt((2*l2+1)*(2*l3+1)))
elseif l4==1 & s~=0
difer(s-2*(-1)^(es+1)*ex...
/sqrt(2*l2*(2*l2+1)*(2*l2+2)*2*l3*(2*l3+1)*(2*l3+2)))
elseif l4==2 & s~=0
difer(s-(2*(-1)^es*(3*ex*(ex-1)-4*l2*(l2+1)*l3*(l3+1)))...
/sqrt((2*l2-1)*2*l2*(2*l2+1)*(2*l2+2)*(2*l2+3)*...
(2*l3-1)*2*l3*(2*l3+1)*(2*l3+2)*(2*l3+3)))
end
disp('Passed excessive verification of Edmonds'' special cases')
end
elseif strcmp(l1,'demo1')
l1=0:15;
l2=8; l3=7; l4=13; l5=15; l6=15;
W=wigner6j(l2,l3,l4,l5,l6);
l2=repmat(l2,size(l1));
l3=repmat(l3,size(l1));
l4=repmat(l4,size(l1));
l5=repmat(l5,size(l1));
l6=repmat(l6,size(l1));
difer(W-sixj(l1,l2,l3,l4,l5,l6))
difer(W-[0 (1/8)*sqrt(259/5270) -(11/24)*sqrt(407/152830) ...
(1/204)*sqrt(54131/116870) (1/68)*sqrt(23199/23374) ...
-(1/1768)*sqrt(142709/186) (59/8840)*sqrt(209/5394) ...
(997/2210)*sqrt(33/41354) -(1429/3910)*sqrt(77/70122) ...
-(7747/15640)*sqrt(1/70122) (269/15640)*sqrt(3157/5394)...
-(83/70380)*sqrt(13981/290) -(1/70380)*sqrt(88398583/310) ...
(149/1360680)*sqrt(4790071/310) (13/240120)*sqrt(368467/1054) ...
-(13/40455)*sqrt(2579269/782)])
% Try various methods as well
l1=0:6;
l2=8; l3=7; l4=6; l5=4; l6=2;
W=wigner6j(l2,l3,l4,l5,l6);
l2=repmat(l2,size(l1));
l3=repmat(l3,size(l1));
l4=repmat(l4,size(l1));
l5=repmat(l5,size(l1));
l6=repmat(l6,size(l1));
difer(W-[0 0 -(1/42)*sqrt(11/34) (1/21)*sqrt(209/442) -(5/21)*sqrt(19/442) ...
(5/6)*sqrt(19/4862) -(1/78)*sqrt(133/17)])
difer(W-sixj(l1,l2,l3,l4,l5,l6))
difer(sixj(4,0,4,13,9,13)-1/9/sqrt(3))
difer(sixj(30,30,30,30,30,30)-4.102353215741e-04)
end