function sj = sf_sixj(qnarray)
//SF_SIXJ([J1 J2 J3; J4 J5 J6])
//  Wigner 6j-symbol, integer and half integer quantum numbers.
//
//  See also SF_NOVERK, SF_THREEJ, SF_NINEJ

//  GCG, 11-Mar-2003

if or(size(qnarray)<>[2 3])
  error('qnarray should be 2x3')
end

j1 = qnarray(1,1)
j2 = qnarray(1,2)
j3 = qnarray(1,3)
j4 = qnarray(2,1)
j5 = qnarray(2,2)
j6 = qnarray(2,3)

sj = 0
k1 = -j1+j2+j3
if k1<>int(k1), error('j1+j2+j3 must be integer'), end
k2 =  j1-j2+j3
k3 =  j1+j2-j3
l1 = -j1+j5+j6
if l1<>int(l1), error('j1+j5+j6 must be integer'), end
l2 =  j1-j5+j6
l3 =  j1+j5-j6
m1 = -j4+j2+j6
if m1<>int(m1), error('j2+j4+j6 must be integer'), end
m2 =  j4-j2+j6
m3 =  j4+j2-j6
n1 = -j4+j5+j3
n2 =  j4-j5+j3
n3 =  j4+j5-j3
if n3<>int(n3), error('j3+j4+j5 must be integer'), end
// Triangular conditions OK?
if min([k1,k2,k3,l1,l2,l3,m1,m2,m3,n1,n2,n3])<0, return, end

i1   = j1+j2+j4+j5+1
imin = max([0,k3-m1,k3-l2])
imax = min([k3,l3,m3,n3,i1])
if imin>imax, return, end
i2   = j2+j4+j6+1
fact = sqrt( sf_noverk(i2,m1) * sf_noverk(j4+j4,n2) / ...
            ((k1+1) * (l1+1) * sf_noverk(j4+j4,m3) * ...
             sf_noverk(j3+j4+j5+1,n1) * sf_noverk(j1+j2+j3+1,k3) * ...
             sf_noverk(j3+j3+1,k2) * sf_noverk(j1+j5+j6+1,l3) * ...
             sf_noverk(j6+j6+1,l2)))
phase = 1
if pmodulo(imin+i1,2) <> 0
  phase = -phase
end
 for iteller=imin:imax
  phase = -phase
  sj    = sj + phase * sf_noverk(n3,iteller) * sf_noverk(m1,k3-iteller) * ...
                       sf_noverk(n2,m3-iteller) * sf_noverk(i1-iteller,i2)
end
 
sj = sj * fact
endfunction
