Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
158 lines (147 sloc) 4.91 KB
function [s,C,S,L]=zeroj(l1,l2,l3,L,meth,C,S)
% [s,C,S,L]=ZEROJ(l1,l2,l3,L,meth,C,S)
%
% Wigner 3j-symbol from a database precomputed by WIGNERCYCLE...
% with bottom row (spherical harmonic orders) all zeros!
%
% INPUT:
%
% l1,l2,l3 Top row of the Wigner 3j symbol [may be 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 3j symbols with zero bottom row
% C,S The column and element vectors good for the next load
% L The L of the database that is loaded, should you not know
%
% EXAMPLE:
%
% zeroj('demo1') % Should return nothing if it all works
% [jk,C0,S0]=zeroj(0,0,22) % Loads best available database for later use
% [jk,C0,S0]=zeroj(0,0,0,22) % Loads/Makes the database for later use
%
% l=32; p=56;
% difer(1-sum((2*[0:l+p]+1).*zeroj(l,p,[0:l+p]).^2))
%
% SEE ALSO: THREEJ, SIXJ, WIGNERCYCLE
%
% Last modified by fjsimons-at-alum.mit.edu, 04/09/2007
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))
error('All degrees 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','WIGNER0JCS-*-C'));
% Here it breaks if there is NOTHING
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)))
else
defval('L',-1)
% Maybe here need to do wignercycle and have a go again?
end
% Check, identify and load database
if any([l1(:)' l2(:)' l3(:)']>L)
error('Insufficient bandwidth for database')
% HERE WE SHOULD RUN WIGNERCYCLE
end
fnpl1=sprintf('%s/WIGNER0JCS-%i-C',...
fullfile(getenv('IFILES'),'WIGNER'),L);
fnpl2=sprintf('%s/WIGNER0JCS-%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
% Precompute the database - with bottom rows of zero
wignercycle(L,0);
% And have a go again
s=zeroj(l1,l2,l3,L);
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
if any([l1(:)' l2(:)' l3(:)']>L)
error('Insufficient bandwidth for database')
end
end
% Find running index into this matrix
% So you can make a vector of l1 l2 l3 and have all zero m's
rind=addmabout(L,l1)+(L+1)^2*(addmabout(L,l2)-1)+...
(L+1)^4*(addmabout(L,l3)-1);
% 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)^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
elseif strcmp(l1,'demo1')
difer(threej(16,16,6,0,0,0)-zeroj(16,16,6))
difer(threej([2 0],[1 1],[1 1],0,0,0)-zeroj([2 0],[1 1],[1 1]))
difer(threej(0,1,1,0,0,0)-zeroj(0,1,1))
difer(threej(2,1,1,0,0,0)-zeroj(2,1,1))
difer(threej(20,10,10,0,0,0)-zeroj(20,10,10))
difer(threej(20,8,12,0,0,0)-zeroj(20,8,12))
difer(threej(30,21,12,0,0,0)-zeroj(30,21,12))
difer(threej(27,19,8,0,0,0)-zeroj(27,19,8))
difer(threej(0:40,10,13,0,0,0)-zeroj(0:40,10,13))
end