%function threej by J. Kay 

%Calculates 3-j coefficients.  Also tests input angular momenta for

%triangular conditions; explicitly accounts for symmetry.  The function

%returns a value of "0" if the angular momenta violate triangular

%conditions.  Uses (3.6.10) and the notation on p. 44 of Edmonds.

%

%To use, copy the file "threej.m" to the matlab\work

%directory and call from main program.

%

%Usage:    In your main matlab program, type

%             threej(j1, j2, j3, m1, m2, m3)

%          where j1, j2, and j3 are the angular momenta, and m1, m2, and m3

%          are the projections of j1, j2, and j3, respectively.  

%          The function will return the value of the 3-j symbol to the

%          program that called it.  

%     

%Example:  The statement 

%             x = threej(1,2,3,0,0,0)

%          assigns the value -0.29277002188456 to the variable x.



function threej = threej(j1,j2,j3,m1,m2,m3)





%Test arguments

if j1<0 | j2<0 | j3<0

    1;

    threej = 0;

end

if j1<m1 | j2<m2 | j3<m3

    2;

    threej = 0;

end

if m1+m2+m3~=0

    3;

    threej = 0;

end

if abs(j1-j2) > j3 | j3>j1+j2

    4;

    threej = 0;

end



%Assuming arguments pass, calculate threej coefficient

if j1>=0 & j2>=0 & j3>=0 & j1>=m1 & j2>=m2 & j3>=m3 & m1+m2+m3==0 & abs(j1-j2)<=j3 & j3<=j1+j2

    5;

    m3 = -m3;

    %calculation of first factor

    factor1 = 1;

 

    %calculation of second factor

    factor2 = sqrt( ( (2*j3+1)*factorial(j1+j2-j3)*factorial(j1-m1)*factorial(j2-m2)*factorial(j3+m3)*factorial(j3-m3) )/( factorial(j1+j2+j3+1)*factorial(j1-j2+j3)*factorial(-j1+j2+j3)*factorial(j1+m1)*factorial(j2+m2)) );

    

    %calculation of fourth factor

    factor3 = 0;

    for s = 0:100

        if j1-m1-s < 0

            1;

            continue

        end

        if j3-m3-s < 0

            2;

            continue

        end

        if j2-j3+m1+s < 0

            3;

            continue

        end

       

        factor3 = factor3 + (-1)^(s+j1-m1) * (factorial(j1+m1+s)*factorial(j2+j3-m1-s))/(factorial(s)*factorial(j1-m1-s)*factorial(j3-m3-s)*factorial(j2-j3+m1+s));

        

    end

    

    %Calculate C-G coefficient

    CG = factor1 * factor2 * factor3;

    

    %calculate 3-j from C-G

    threej = CG * (-1)^(j1-j2-m3)*(2*j3+1)^(-1/2);

    

end