Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
226 lines (210 sloc) 7.3 KB
function wignercycle(L,emnz,simi)
% WIGNERCYCLE(L,emnz,simi)
%
% If WIGNER0J is the wheel, this is the bicycle, and ZEROJ the rider.
% If WIGNER3JM is the wheel, this is the bicycle, and THREEJ the rider.
% If WIGNER6J is the wheel, this is the bicycle, and SIXJ the rider.
%
% INPUT:
%
% L Bandwidth of the data base that will be produced
% emnz 1 Do this for all possible orders [default]
% 0 Only do this for zero orders
% 6 Rather calculate the wigner 6j coefficients
% simi 0 Standard, very cumbersome storage scheme
% 1 Smart storage scheme a la Rasch and Yu (2003) for 3j only
% xver 1 Perform excessive verification (a wise choice)
% 0 Don't (an unwise choice)
%
% OUTPUT:
%
% Only to file, suitable to be read in by ZEROJ, THREEJ, SIXJ
%
% Creates a set of huge, but sparse matrices of Wigner 3j/6j
% coefficients... up to a maximum degree L, in blocky ADDMOUT/ADDMABOUT
% ordering (for the 3j's). Saves them in a vector with column numbers and
% with data. For 3j symbols, when simi=1, uses the indexing scheme
% described by Rasch and Yu (2003), to take account of natural
% symmetries. In this case, the values are stored as a sparse array in
% .mat format.
%
% All degrees/orders INTEGERS... not required by WIGNER3JM/WIGNER6J.
%
% See ZEROJ, THREEJ, SIXJ, WIGNERCYCLESYM, WIGNERSORT
%
% Last modified by klewis-at-princeton.edu, 06/14/2010
% Last modified by fjsimons-at-alum.mit.edu, 05/24/2011
t=cputime;
defval('timewaster',0)
defval('L',30)
defval('emnz',1)
defval('simi',1)
defval('xver',0)
% Files to save the results in: one float, one integer
switch emnz
case 1
if simi==0
fnpl1=sprintf('%s/WIGNER3JCS-%i-C',...
fullfile(getenv('IFILES'),'WIGNER'),L);
fnpl2=sprintf('%s/WIGNER3JCS-%i-S',...
fullfile(getenv('IFILES'),'WIGNER'),L);
elseif simi==1
fnpl=sprintf('%s/WIGNER3J-%i-sym',...
fullfile(getenv('IFILES'),'WIGNER'),L);
end
% Get orders and degrees in block form
[EM,EL,mz,blkm,dblk]=addmout(L);
EM=EM(blkm); EL=EL(blkm);
case 0
fnpl1=sprintf('%s/WIGNER0JCS-%i-C',...
fullfile(getenv('IFILES'),'WIGNER'),L);
fnpl2=sprintf('%s/WIGNER0JCS-%i-S',...
fullfile(getenv('IFILES'),'WIGNER'),L);
% Restrict to the zero orders only
EM=repmat(0,1,L+1)';
EL=[0:L]';
case 6
fnpl1=sprintf('%s/WIGNER6JCS-%i-C',...
fullfile(getenv('IFILES'),'WIGNER'),L);
fnpl2=sprintf('%s/WIGNER6JCS-%i-S',...
fullfile(getenv('IFILES'),'WIGNER'),L);
otherwise
error('Specify valid case')
end
% For format, see INTMAX and REALMAX
if simi==0
fid1=fopen(fnpl1,'w'); fmt1='uint64';
fid2=fopen(fnpl2,'w'); fmt2='float64';
end
if emnz~=6
% Loop over m2, then over m3, then pick m1
for ind2=1:length(EM)
l2=EL(ind2); m2=EM(ind2);
% Only pick an l3 if its m3 implies an abs(m1) that is
% not bigger than l1
EL3AL=EL(abs(-m2-EM)<=L);
EM3AL=EM(abs(-m2-EM)<=L);
for ind3=1:length(EL3AL)
l3=EL3AL(ind3); m3=EM3AL(ind3);
m1=-m2-m3;
% if timewaster==1
% % Note that display time is expensive!
% clc
% disp(sprintf('(%+3i %+3i %+3i)',L,l2,l3))
% disp(sprintf('(%+3i %+3i %+3i)',m1,m2,m3))
% end
% Do the calculation
switch emnz
case 1
[w3j,l1]=wigner3jm(L,l2,l3,m1,m2,m3);
if simi==1
% So here is a modification by Kevin W. Lewis
[c,oddperms,phasefix]=wignersort(l1,l2,l3,m1,m2,m3);
% See also THREEJ
w3j(oddperms)=w3j(oddperms).*phasefix;
jmin=max(abs(l2-l3),abs(m1))+1;
jmax=min(l2+l3+1,length(c));
if exist('w3js','var') && length(w3js)>=max(c)
if sum(w3js(c(jmin:jmax))~=0 & ...
w3js(c(jmin:jmax))-w3j(jmin:jmax)'>1e-10)
disp(['Differences = ' num2str(w3js(c)-w3j)])
end
end
if ~exist('w3js','var')
% Rasch and Yu (2003) eq. (2.14)
w3js=sparse(...
2*L*(274+2*L*(225+2*L*(85+2*L*(15+2*L))))/120+1,1);
end
w3js(c(jmin:jmax),1)=w3j(jmin:jmax);
if ~exist('ind','var')
ind=0;
end
totalindex(c(jmin:jmax))=true;
ind=ind+jmax-jmin;
% disp(num2str([ind]))
if ind2==length(EM) && ind3==length(EL3AL)
save(fnpl,'w3js');
end
elseif simi==0
% Find the possible and non-trivial ones in this
% If there are any, save them in a big array
nontriv=w3j~=0;
if sum(nontriv)>0
% Now construct running indices for the non-zeroes only
% This is the 3-D equivalent for SUB2IND which you may use
C= addmabout(L,l1(nontriv),m1)...
+(L+1)^2*(addmabout(L,l2 ,m2)-1)...
+(L+1)^4*(addmabout(L,l3 ,m3)-1);
% Write this out on the fly
% The column number, "C"
fwrite(fid1,C,fmt1);
% The actual element, "S"
fwrite(fid2,w3j(nontriv),fmt2);
end
end
case 0
% Perform the orthogonality check or not
[w3j,l1]=wigner0j(L,l2,l3,xver);
if xver==1
% Check, and check again with the different routine
difer(w3j-wigner3jm(L,l2,l3,0,0,0));
end
% Find the possible and non-trivial ones in this
% If there are any, save them in a big array
nontriv=w3j~=0;
if sum(nontriv)>0
% Now construct running indices for the non-zeroes only
% This is the 3-D equivalent for SUB2IND which you may use
C= addmabout(L,l1(nontriv),m1)...
+(L+1)^2*(addmabout(L,l2,m2)-1)...
+(L+1)^4*(addmabout(L,l3,m3)-1);
% Write this out on the fly
fwrite(fid1,C,fmt1); % The column number, "C"
fwrite(fid2,w3j(nontriv),fmt2); % The actual element, "S"
end
end
end
end
elseif emnz==6
% Here we can avoid loops since the arrays are much smaller, and then
% we can use the selection rules directly; may use NDGRID for this
l2= gamini([0:L],(L+1)^4)';
l3=repmat(gamini([0:L],(L+1)^3)',(L+1)^1,1);
l4=repmat(gamini([0:L],(L+1)^2)',(L+1)^2,1);
l5=repmat(gamini([0:L],(L+1)^1)',(L+1)^3,1);
l6=repmat([0:L]', (L+1)^4,1);
% Selection rules - those of the 6j contain those of the 6j
selx=triangle(l4,l2,l6) & triangle(l4,l5,l3);
% Apply the selection rules
l2=l2(selx); l3=l3(selx);
l4=l4(selx); l5=l5(selx); l6=l6(selx);
for ix=1:length(l2)
% Calculate the Wigner 6j symbols with normalization check
[w6j,l1,norma]=wigner6j(l2(ix),l3(ix),l4(ix),l5(ix),l6(ix));
% Now remember that l1 returned may exceed L of database, cut
if l1(end)>L
w6j=w6j(1:L+1);
l1=l1(1:L+1);
end
% disp(sprintf('Passed normalization with %8.3e',norma-1))
nontriv=w6j~=0;
if sum(nontriv)>0
% Now construct running indices for the non-zeroes only
% This is the 3-D equivalent for SUB2IND which you may use
C=1+l1(nontriv)+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;
% Write this out on the fly
fwrite(fid1,C,fmt1); % The column number, "C"
fwrite(fid2,w6j(nontriv),fmt2); % The actual element, "S"
end
end
end
% Close files
if simi==0
fclose(fid1);
fclose(fid2);
end
disp(sprintf('WIGNERCYCLE: Elapsed CPU time %8.3f',cputime-t))