Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
330 lines (309 sloc) 11.7 KB
function [s,C,S,L]=threej(l1,l2,l3,m1,m2,m3,L,meth,C,S)
% [s,C,S,L]=THREEJ(l1,l2,l3,m1,m2,m3,L,meth,C,S)
%
% Wigner 3j-symbol from a database precomputed by WIGNERCYCLE. Tries to
% find files which take account of the Regge symmetries as described by
% Rasch and Yu (2003) first.
%
% INPUT:
%
% l1,l2,l3 Top row of the Wigner 3j symbol [may be same-length vectors]
% m1,m2,m3 Bottom row of the Wigner 3j symbol [may be same-length vectors]
% Their default is 0,0,0
% 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 3j symbols
% C,S The column and element vectors good for the next load, OR
% The entire Wigner 3j symmetric database and an additional empty
% L The L of the database that is loaded, should you not know
%
% EXAMPLE:
%
% threej('demo1') % Should return nothing if it all works
%
% SEE ALSO: WIGNER3JM, GAUNT, WIGNER0J, GUSEINOV, ZEROJ
%
% Last modified by kwlewis-at-princeton.edu, 06/14/2010
% Last modified by plattner-at-princeton.edu, 05/24/2011
% Last modified by fjsimons-at-alum.mit.edu, 07/11/2011
% Sometimes it returns zero if the database is corrupted, i.e. if it's
% not what it purports to be.
if ~isstr(l1)
defval('C',[])
defval('S',[])
% Method
defval('meth',2)
% disp(sprintf('Using method %i',meth))
defval('m1',0)
defval('m2',0)
defval('m3',0)
% All saved values must be integers
if sum(mod(l1,1)) || sum(mod(l2,1)) || sum(mod(l3,1)) ...
sum(mod(m1,1)) || sum(mod(m2,1)) || sum(mod(m3,1))
error('All degrees and orders must be integers for the database')
end
if any(m1+m2+m3)
error('Sum of the orders must be zero')
end
if isempty(C) && isempty(S)
% Got something to do
% These are the modifications by Kevin W. Lewis
% First check for saved files which take account of Regge symmetries
try
Els=ls2cell(fullfile(getenv('IFILES'),...
'WIGNER','WIGNER3J-*-sym.mat'));
catch
Els=[];
EL=[];
end
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(:)'])<=EL);
if ~isempty(fmax)
defval('L',EL(indeks(fmax,1)))
% Check, identify and load database
if any([l1(:)' l2(:)' l3(:)']>L)
error('Insufficient bandwidth for database')
end
fnpl=sprintf('%s/WIGNER3J-%i-sym',...
fullfile(getenv('IFILES'),'WIGNER'),L);
disp(sprintf('Loading %s',fnpl))
load(fnpl)
% These must be same-length vectors but two pairs of them can be
% scalars
l1=l1(:); l2=l2(:); l3=l3(:);
m1=m1(:); m2=m2(:); m3=m3(:);
% But if ONE of them is a vector and the others aren't that's fine
% too, as in the old version. Just take care of this one case now
if length(l1)~=1 && length(l2)==1 && length(l2)==1
l2=repmat(l2,length(l1),1);
l3=repmat(l3,length(l1),1);
elseif length(l2)~=1 && length(l1)==1 && length(l3)==1
l1=repmat(l1,length(l2),1);
l3=repmat(l3,length(l2),1);
elseif length(l3)~=1 && length(l1)==1 && length(l2)==1
l1=repmat(l1,length(l3),1);
l2=repmat(l2,length(l3),1);
end
% Find where they're being kept
[C,oddperm,phasefix]=wignersort(l1,l2,l3,m1,m2,m3);
% Do the initial evaluation from the loaded variable
s=full(w3js(C));
% Fix the phase
s(oddperm)=s(oddperm).*phasefix;
% Now fix the triangle condition violations
s=s.*triangle(l1,l2,l3);
% Now fix the order violations
s=s.*~[l1<abs(m1) | l2<abs(m2) | l3<abs(m3)];
% There may be other violations that aren't being tested! The
% database only stores the nonzero ones, for whichever reason
s=s(:)';
% If you want this to return the database, you need to supply variables
C=w3js;
S=[];
return
else
% Try for files which do not take account of symmetries
% What is the lowest of the available database bandwidths?
try
Els=ls2cell(fullfile(getenv('IFILES'),...
'WIGNER','WIGNER3JCS-*-C'));
catch
Els=[];
EL=[];
end
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(:)'])<=EL);
if ~isempty(fmax)
defval('L',EL(indeks(fmax,1)))
% Check, identify and load database
if any([l1(:)' l2(:)' l3(:)']>L)
error('Insufficient bandwidth for database')
end
fnpl1=sprintf('%s/WIGNER3JCS-%i-C',...
fullfile(getenv('IFILES'),'WIGNER'),L);
fnpl2=sprintf('%s/WIGNER3JCS-%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));
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);
writeb(C,fnpl1,fmt1)
S=S(j);
writeb(S,fnpl2,fmt2); clear j
disp('Column arrays now sorted once and for all')
end
else
defval('L',-1)
% Precompute the database, this will fail though, on purpose
wignercycle(L);
% And have a go again
s=threej(l1,l2,l3,m1,m2,m3,L,meth);
end
else
defval('L',-1)
% Precompute the database, this will fail though, on purpose
wignercycle(L);
% And have a go again
s=threej(l1,l2,l3,m1,m2,m3,L,meth);
end
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
if isempty(S)
% This is the symmetric storage scheme
% These must be same-length vectors but two pairs of them can be
% scalars
l1=l1(:); l2=l2(:); l3=l3(:);
m1=m1(:); m2=m2(:); m3=m3(:);
% But if ONE of them is a vector and the others aren't that's fine
% too, as in the old version. Just take care of this one case now
if length(l1)~=1 && length(l2)==1 && length(l2)==1
l2=repmat(l2,length(l1),1);
l3=repmat(l3,length(l1),1);
elseif length(l2)~=1 && length(l1)==1 && length(l3)==1
l1=repmat(l1,length(l2),1);
l3=repmat(l3,length(l2),1);
elseif length(l3)~=1 && length(l1)==1 && length(l2)==1
l1=repmat(l1,length(l3),1);
l2=repmat(l2,length(l3),1);
end
% Find where they're being kept
[CC,oddperm,phasefix]=wignersort(l1,l2,l3,m1,m2,m3);
% Do the initial evaluation from the loaded variable
s=full(C(CC));
% Fix the phase
s(oddperm)=s(oddperm).*phasefix;
% Now fix the triangle condition violations
s=s.*triangle(l1,l2,l3);
% Now fix the order violations
s=s.*~[l1<abs(m1) | l2<abs(m2) | l3<abs(m3)];
% There may be other violations that aren't being tested! The
% database only stores the nonzero ones, for whichever reason
s=s(:)';
return
else
% This is the old-school storage scheme
% Now to the calculation proper if you hadn't stopped before
% Find running index into this matrix
if length(m1)~=1 || length(m2)~=1 || length(m3)~=1
% Note there must be the same number of l1, l2, l3, m1, m2, m3
% Initialize index vector
rind=repmat(NaN,1,length(l1));
for index=1:length(l1)
rind(index)=addmabout(L,l1(index),m1(index))+...
(L+1)^2*(addmabout(L,l2(index),m2(index))-1)+...
(L+1)^4*(addmabout(L,l3(index),m3(index))-1);
end
else
% So you can make a vector of l1 l2 l3 and have all zero m's
% Turns out you can still have ONE that is a vector and all the
% other ones that are singles, that would work too
rind=addmabout(L,l1,m1)+(L+1)^2*(addmabout(L,l2,m2)-1)+...
(L+1)^4*(addmabout(L,l3,m3)-1);
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
% HOLD ON - SHOULDN't THIS BE TO THE POWER SIX, RATHER?
W=sparse(1,C,S,1,(L+1)^8);
% Extract the Wigner 3j-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
end
elseif strcmp(l1,'demo1')
difer(full(threej([8 7 6 5 4],6,5,3,2,-5))-...
[sqrt(42/4199) 35*sqrt(2/138567) ...
7*sqrt(1/2431) sqrt(35/2431) sqrt(5/858)])
% From Lai et al. based on Rotenberg's book
difer(threej(16,16,6,0,0,0)-rotenberg([4 -2 2 0 -1 -1 1 1 0 -1 -1 -1],1))
% For database of a bandwidth at least 3
% Illustrates vector capabilities
difer(threej([0 0],[1 1],[1 1],[0 0],[1 0],[-1 0])-[sqrt(1/3) -sqrt(1/3)])
difer(threej([2 0],[1 1],[1 1],0,0,0)-[sqrt(2/15) -sqrt(1/3)])
% Illustrates straight scalar capabilities
difer(threej(0,1,1,0,1,-1)-sqrt(1/3))
difer(threej(0,1,1,0,0,0)--sqrt(1/3))
difer(threej(1,1,1,0,1,-1)-sqrt(1/6))
difer(threej(2,1,1,0,1,-1)-sqrt(1/30))
difer(threej(2,1,1,0,0,0)-sqrt(2/15))
difer(threej(2,2,3,-2,0,2)--sqrt(1/14))
difer(threej(3,2,2,1,0,-1)--sqrt(1/35))
difer(threej(3,2,3,-2,0,2)-0)
difer(threej(3,2,3,3,0,-3)-(1/2)*sqrt(5/21))
difer(threej(3,6,5,3,2,-5)-sqrt(1/1001))
% For database of a bandwidth at least 6
difer(threej(6,2,4,-4,1,3)-4*sqrt(1/429))
% For database of a bandwidth at least 10
difer(threej(10,10,12,9,3,-12)--(1/15)*sqrt(4199/9889))
% For database of a bandwidth at least 20
difer(threej(20,15,9,-3,2,1)+(311/2)*sqrt(115/1231496049))
difer(threej(20,10,10,0,0,0)-indeks(wigner0j(20,10,10),'end'))
difer(threej(20,8,12,-1,-3,4)- -(3/5)*sqrt(2261/1363783))
% For database of a bandwidth at least 30
difer(threej(30,21,12,-22,19,3)--(151/2)*sqrt(527/1277814153))
difer(threej(27,19,8,-21,18,3)-(4/5)*sqrt(38/208131))
% For database of a bandwidth at least 40
difer(threej(0:40,10,13,0,0,0)-wigner0j(40,10,13))
% For database of a bandwidth at least 60
difer(threej(60,60,60,-60,40,20)-indeks(wigner3jm(60,60,60,-60,40,20),'end'))
difer(threej(59,60,40,-20,10,10)-indeks(wigner3jm(59,60,40,-20,10,10),'end'))
% Look at symmetry
threej(repmat(20,1,19),repmat(15,1,19),repmat(9,1,19),...
repmat(-3,1,19),repmat(2,1,19),[-9:9])
end