function nj = sf_ninej(qnarray)
//SF_NINEJ(QNARRAY)  Wigner's 9-j symbol.
//  SF_NINEJ(QNARRAY) where:
//  QNARRAY = [ J11 J12 J13
//              J21 J22 J23
//              J31 J32 J33]
//  returns the 9-j symbol for both integer and half-integer
//  quantum number. The function is adapted from the fortran routine
//  FUNCTION Z9_J(J11,J12,J13,J21,J22,J23,J31,J32,J33)
//
//  From: Quantum theory of Angular Momentum,
//        D. A. Varshalovich, A. N. Moskalev, V. K. Kheersonskii,
//        10.2.4, eq (20), p 340
//
//  See also SF_THREEJ, AND SF_SIXJ

//  GCG, 18 october 1999, 11-Mar-2003
 
// get the quantumnumbers from the input array
[m,n] = size(qnarray)
if m<>3 | n <> 3
  error ('qnarray should be 3x3')
end
j11 = qnarray(1,1)
j12 = qnarray(1,2)
j13 = qnarray(1,3)
j21 = qnarray(2,1)
j22 = qnarray(2,2)
j23 = qnarray(2,3)
j31 = qnarray(3,1)
j32 = qnarray(3,2)
j33 = qnarray(3,3)
 
nj    = 0
izmin = max([abs(j11-j33),abs(j32-j21),abs(j12-j23)])
izmax = min([j11+j33, j32+j21, j12+j23])

dim   = 2*izmin + 1
for iz = izmin:izmax
  sxj1 = sf_sixj([j11,j12,j13;j23,j33, iz])
  sxj2 = sf_sixj([j21,j22,j23;j12, iz,j32])
  sxj3 = sf_sixj([j31,j32,j33; iz,j11,j21])
  nj   = nj + dim * sxj1 * sxj2 * sxj3
  dim  = dim + 2
end
nj    = nj*(-1)^(2*izmin)

endfunction
