function wignercycle(L,emnz)
% WIGNERCYCLE(L,emnz)
%
% If WIGNER0J is the wheel, this is the bicycle.
% If WIGNER3JM is the wheel, this is the bicycle.
% If WIGNER6J is the wheel, this is the bicycle.
%
% 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
% xver   1 Perform excessive verification (a wise choice)
%        0 Don't (an unwise choice)
%              
%
% Crates 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. 
%
% All degrees/orders INTEGERS... not required by WIGNER3JM/WIGNER6JM.
%
% See ZEROJ, THREEJ, SIXJ, WIGNER0J, WIGNER3JM, WIGNER3J, WIGNER6J
% 
% Last modified by fjsimons-at-alum.mit.edu, 1/24/2007

t=cputime;

defval('timewaster',0)
defval('L',30)
defval('emnz',1)
defval('xver',0)

% Files to save the results in: one float, one integer
switch emnz
 case 1
   fnpl1=sprintf('%s/WIGNER3JCS-%i-C',...
		 fullfile(getenv('IFILES'),'WIGNER'),L);
   fnpl2=sprintf('%s/WIGNER3JCS-%i-S',...
		 fullfile(getenv('IFILES'),'WIGNER'),L);
   % 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
fid1=fopen(fnpl1,'w'); fmt1='uint64';
fid2=fopen(fnpl2,'w'); fmt2='float64';

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);
       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
      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
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
fclose(fid1);
fclose(fid2);

disp(sprintf('WIGNERCYCLE: Elapsed CPU time %8.3f',cputime-t))