%function sixj by J. Kay 

%Calculates 6-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 Racah's formula (6.3.7) and notation given on page 99 of Edmonds.

%

%To use, copy the files "sixj.m" and "deltasixj.m" to the matlab\work

%directory and call from main program.

%

%Usage:    In your main matlab program, type

%             sixj(j1, j2, j3, L1, L2, L3)

%          where j1, j2, j3, L1, L2, and L3 are the angular momenta.  

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

%          program that called it.  

%     

%Example:  The statement 

%             x = sixj(2,2,2,1,1,1)

%          assigns the value 0.15275252316519 to the variable x.





function sixj = sixj(j1,j2,j3,L1,L2,L3)



if j1<0 | j2<0 | j3<0 | L1<0 | L2<0 | L3<0

    sixj=0;

end

if abs(j1-j2)>j3

    sixj=0;

end

if j3>j1+j2

    sixj=0;

end

if abs(L2-L3)>j1

    sixj=0;

end

if j1>L2+L3

    sixj=0;

end

if abs(L1-L3)>j2

    sixj=0;

end

if j2>L1+L3

    sixj=0;

end

if abs(L1-L2)>j3

    sixj=0;

end

if j3>L1+L2

    sixj=0;

end



if j1>=0 & j2>=0 & j3>=0 & L1>=0 & L2>=0 & L3>=0 & abs(j1-j2)<=j3 & j3<=j1+j2 & abs(L2-L3)<=j1 & j1<=L2+L3 & abs(L1-L3)<=j2 & j2<=L1+L3 & abs(L1-L2)<=j3 & j3<=L1+L2

    factor1 = deltasixj(j1,j2,j3);

    factor2 = deltasixj(j1,L2,L3);

    factor3 = deltasixj(L1,j2,L3);

    factor4 = deltasixj(L1,L2,j3);

    

    factor5 = 0;

    for z = 0:1000

        if z-j1-j2-j3 < 0

            continue

        end

        if z-j1-L2-L3 < 0

             continue

        end

        if z-L1-j2-L3 < 0

            continue

        end

        if z-L1-L2-j3 < 0

            continue

        end

        if j1+j2+L1+L2-z < 0

            continue

        end

        if j2+j3+L2+L3-z < 0 

            continue

        end

        if j3+j1+L3+L1-z < 0

            continue

        end

        factor5 = factor5 + (-1)^(z)*factorial(z+1) / (factorial(z-j1-j2-j3)*factorial(z-j1-L2-L3)*factorial(z-L1-j2-L3)*factorial(z-L1-L2-j3)*factorial(j1+j2+L1+L2-z)*factorial(j2+j3+L2+L3-z)*factorial(j3+j1+L3+L1-z) );

    end



    sixj = factor1*factor2*factor3*factor4*factor5;

end



