function ClebschGordanKoeff = ClebschGordan(j1,j2)
% Berechnen der Clebsch-Gordan Koeffizienten
% Aufruf CB =ClebschGordan(j1,j2)
% Beispiel: j1 = 1/2 j2=1
% Rueckgabe Tabelle mit den Clebsch-Gordan Koeff.

% Vers. 2.2 18.08.2016
% Copyright Wolfgang Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

%j1 = 1/2;
%j2 = 1;
Ergebnis = [];
%
m1 = -j1:j1;
m2 = -j2:j2;
[M1,M2] = meshgrid(m1,m2);
%%
M3 = M1 + M2;
j3 = abs(j1-j2):j1+j2;
%%
for n=1:length(j3)
    m3j=-j3(n):j3(n);
    vorfakj3 = [j1+j2-j3(n),j3(n)+j1-j2,j3(n)+j2-j1];
    vorfakj3 = factorial(vorfakj3);
    vorfakj3 = sqrt((2*j3(n)+1)*prod(vorfakj3)/factorial(j3(n)+j1+j2+1));
    for k = 1:length(m3j)
        ind = find(M3==m3j(k));
        for l = 1:length(ind)
            sums = [];
            ll = ind(l);
            m1r = M1(ll);
            m2r = M2(ll);
            m3r = M3(ll);
            regge = [-j1+j2+j3(n), j1-j2+j3(n), j1+j2-j3(n) ; ...
                      j1 - m1r   ,   j2-m2r   ,   j3(n)+m3r ; ...
                      j1 + m1r   ,   j2+m2r   ,   j3(n)-m3r];
            if sum(regge(:) < 0) || sum(sum(regge)~=j1+j2+j3(n)) || sum(sum(regge,2)~=j1+j2+j3(n))
               disp('Null!!')
            else
                [j1,m1r;j2,m2r;j3(n),m3r];
                lauf1 = [j1+j2-j3(n),j1-m1r,j2+m2r,j3(n)-j2+m1r,j3(n)-j1-m2r];
                s=abs(min(lauf1));
                laufwhile = logical(1);
                notstop=1;
                vorfak = [j1+m1r,j1-m1r,j2+m2r,j2-m2r,j3(n)+m3r,j3(n)-m3r];
                %
                if sum(vorfak<0)
                    disp('Fehler Vorfaktor')
                else
                    vorfak = factorial(vorfak);
                    vorfak = sqrt(prod(vorfak));
                end
                %   
                while laufwhile
                    test = [j1+j2-j3(n)-s,j1-m1r-s,j2+m2r-s,j3(n)-j2+m1r+s,j3(n)-j1-m2r+s];
                    testok=min(test);
                    zwi1 = [s,test];
                    zwi1 = factorial(zwi1);
                    zwi1 = prod(zwi1);
                    sums(notstop) = (-1)^s * vorfak / zwi1; 
                    s=s+1;
                    test = [j1+j2-j3(n)-s,j1-m1r-s,j2+m2r-s,j3(n)-j2+m1r+s,j3(n)-j1-m2r+s];
                    testok=min(test);
                    notstop = notstop+1;
                    if notstop > 100
                        disp('while-schleife 100 mal durchlaufen')
                        return
                    end
                    laufwhile = testok>=0;
                end
                Summe = sum(sums);
                Ergebnis = [Ergebnis;[j1,j2,j3(n),m1r,m2r,m3r,vorfakj3*Summe]];
                ClebschGordanKoeff = array2table(Ergebnis);
                ClebschGordanKoeff.Properties.VariableNames ={'j1' 'j2' 'J' 'm1' 'm2' 'M' 'CB'};
            end
        end
    end

end
