function cg_m=clebsch_gordan(J);
%cg=clebsch_gordan(qnarray);
%  Compute Clebsch-Gordan coefficient
%  <j1 m1 j2 m2|j3 m3> = cg_m([j1 j2 j3
%                              m1 m2 m3])

%FUNCTION THREEJ_M([J1 J2 J3; M1 M2 M3])
%
%   Wigner 3-j.
%
%   Author: J. Tennyson, modified by P.E.S. Wormer
%   Corrected for half-integers: G.C. Groenenboom
%

% binomial coefficients calculated directly using nchoosek(i,j)

% get the quantumnumbers from the input array

cg_m    = 0;	% GCG, bug fix on March 29, 2020
qnarray = J;
[m,n] = size(qnarray);
if (m ~= 2 | n ~= 3)
  error ('qnarray should be 2x3');
end;
qnarray(2,3) = -qnarray(2,3);  % convert Clebsch-Gordan to threej

j1 = qnarray(1,1);
j2 = qnarray(1,2);
j3 = qnarray(1,3);
m1 = qnarray(2,1);
m2 = qnarray(2,2);
m3 = qnarray(2,3);
if (m1 < -j1 | m1 > j1 ) error('m1 out of range'); end;
if (m2 < -j2 | m2 > j2 ) error('m2 out of range'); end;
if (m3 < -j3 | m3 > j3 ) error('m3 out of range'); end;
if mod(2*j1,2) == 1 && mod(2*m1,2) ~= 1
      error('j1 and m1 must be both (half)integer');
end
if mod(2*j2,2) == 1 && mod(2*m2,2) ~= 1
      error('j2 and m2 must be both (half)integer');
end
if mod(2*j3,2) == 1 && mod(2*m3,2) ~= 1
      error('j3 and m3 must be both (half)integer');
end

tj = 0;
% Test triangular conditions, set i1,...,n3
if ((m1+m2+m3) ~= 0), return; end;
i1 = round(-j1+j2+j3); if (i1<0), return; end;
i2 = round(j1-j2+j3);  if (i2<0), return; end;
i3 = round(j1+j2-j3);  if (i3<0), return; end;
k1 = round(j1+m1);     if (k1<0), return; end;
k2 = round(j2+m2);     if (k2<0), return; end;
k3 = round(j3+m3);     if (k3<0), return; end;
l1 = round(j1-m1);     if (l1<0), return; end;
l2 = round(j2-m2);     if (l2<0), return; end;
l3 = round(j3-m3);     if (l3<0), return; end;
n1 = round(-j1-m2+j3);
n2 = round( m1-j2+j3);
n3 = round( j1-j2-m3);

imin = round(max([-n1,-n2,0]));
imax = round(min([l1,k2,i3]));
if (imin > imax), return; end;
Sign = -1;
for iteller = imin:imax
  Sign = -Sign;
  tj = tj + Sign * nchoosek(i1,n1+iteller) * ...
                   nchoosek(i2,n2+iteller) * nchoosek(i3,iteller);
end;

tj = tj * sqrt( nchoosek(j2+j2,i3)*nchoosek(j1+j1,i2) / ...
               (nchoosek(j1+j2+j3+1,i3)*(j3+j3+1)* ...
                nchoosek(j1+j1,l1)*nchoosek(j2+j2,l2)* ...
                nchoosek(j3+j3,l3)));
if (mod(n3+imin,2) ~= 0)
  tj = -tj;
end;

% convert three-j to Clebsch-Gordan
fac=(-1)^(J(1,1)-J(1,2)+J(2,3))*sqrt(2*J(1,3)+1);
cg_m=fac*tj;
