function [s,C,S]=threej(l1,l2,l3,m1,m2,m3,L,meth,C,S)
% [s,C,S]=THREEJ(l1,l2,l3,m1,m2,m3,L,meth,C,S)
%
% Wigner 3j-symbol from a database precomputed by WIGNERCYCLE.
%
% 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
%
% EXAMPLE:
%
% threej('demo1') % Should return nothing if it all works
%
% SEE ALSO: WIGNER3JM, GAUNT, WIGNER0J, GUSEINOV, ZEROJ
%
% Last modified by fjsimons-at-alum.mit.edu, 02/03/2007

% SHOULD BUILD IN RULES, SELECTION, TRIANGLE, LM RULES BEFORE EVEN
% STARTING
% E.G. this returns something: 2.3.1.1.2.3

if ~isstr(l1)
  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(m1,1)) || sum(mod(m2,1)) || sum(mod(m3,1))
    error('All degrees and orders must be integers for the database')
  end

  defval('m1',0)
  defval('m2',0)
  defval('m3',0)
  
  if isempty(C) && isempty(S)
    % Got something to do
    
    % What is the lowest of the available database bandwidths?
    Els=ls2cell(fullfile(getenv('IFILES'),'WIGNER','WIGNER3JCS-*-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(:)'])<=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(:)']>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
      % Precompute the database
      wignercycle(L);
      % And have a go again
      s=threej(l1,l2,l3,m1,m2,m3,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
  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
    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
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