F=uint32([T(1,1) T(1,2) T(1,3)]);
F=[F; uint32([T(1,1) T(1,2) T(1,4)])];
F=[F; uint32([T(1,1) T(1,3) T(1,4)])];
F=[F; uint32([T(1,2) T(1,3) T(1,4)])];
Nf=4;
T2F=uint32([1 2 3 4]);
F2T=uint32([1 0; 1 0; 1 0; 1 0]);
for i=2:Ne
    Trovata1=uint32(0);
    Trovata2=uint32(0);
    Trovata3=uint32(0);
    Trovata4=uint32(0);
    T2F_app=[0 0 0 0];
    for j=1:Nf
        if (~Trovata1 && ~isempty(find(F(j,:)==T(i,1))) && ~isempty(find(F(j,:)==T(i,2))) && ~isempty(find(F(j,:)==T(i,3))))
            Trovata1=uint32(1);
            T2F_app(1)=j;
            F2T(j,2)=i;
        end
        if (~Trovata2 && ~isempty(find(F(j,:)==T(i,1))) && ~isempty(find(F(j,:)==T(i,2))) && ~isempty(find(F(j,:)==T(i,4))))
            Trovata2=uint32(1);
            T2F_app(2)=j;
            F2T(j,2)=i;
        end
        if (~Trovata3 && ~isempty(find(F(j,:)==T(i,1))) && ~isempty(find(F(j,:)==T(i,3))) && ~isempty(find(F(j,:)==T(i,4))))
            Trovata3=uint32(1);
            T2F_app(3)=j;
            F2T(j,2)=i;
        end
        if (~Trovata4 && ~isempty(find(F(j,:)==T(i,2))) && ~isempty(find(F(j,:)==T(i,3))) && ~isempty(find(F(j,:)==T(i,4))))
            Trovata4=uint32(1);
            T2F_app(4)=j;
            F2T(j,2)=i;
        end
        %if (Trovata1 && Trovata2 && Trovata3 && Trovata4) %Pare + lento se metto questo controllo
        %    break;
        %end
    end
    if (~Trovata1)
        F=[F;uint32([T(i,1) T(i,2) T(i,3)])];
        Nf=Nf+1;
        T2F_app(1)=Nf;
        F2T=[F2T; uint32([i, 0])];
    end
    if (~Trovata2)
        F=[F; uint32([T(i,1) T(i,2) T(i,4)])];
        Nf=Nf+1;
        T2F_app(2)=Nf;
        F2T=[F2T; uint32([i, 0])];
    end
    if (~Trovata3)
        F=[F; uint32([T(i,1) T(i,3) T(i,4)])];
        Nf=Nf+1;
        T2F_app(3)=Nf;
        F2T=[F2T; uint32([i, 0])];
    end
    if (~Trovata4)
        F=[F; uint32([T(i,2) T(i,3) T(i,4)])];
        Nf=Nf+1;
        T2F_app(4)=Nf;
        F2T=[F2T; uint32([i, 0])];
    end
    T2F=[T2F; uint32(T2F_app)];
end

app=sqrt(det([y(F(1,1)) z(F(1,1)) 1; y(F(1,2)) z(F(1,2)) 1; y(F(1,3)) z(F(1,3)) 1])^2+det([z(F(1,1)) x(F(1,1)) 1; z(F(1,2)) x(F(1,2)) 1; z(F(1,3)) x(F(1,3)) 1])^2+det([x(F(1,1)) y(F(1,1)) 1; x(F(1,2)) y(F(1,2)) 1; x(F(1,3)) y(F(1,3)) 1])^2);
PF=[app (x(F(1,1))+x(F(1,2))+x(F(1,3)))/3 (y(F(1,1))+y(F(1,2))+y(F(1,3)))/3 (z(F(1,1))+z(F(1,2))+z(F(1,3)))/3];
for i=1:Nf
    app=sqrt(det([y(F(i,1)) z(F(i,1)) 1; y(F(i,2)) z(F(i,2)) 1; y(F(i,3)) z(F(i,3)) 1])^2+det([z(F(i,1)) x(F(i,1)) 1; z(F(i,2)) x(F(i,2)) 1; z(F(i,3)) x(F(i,3)) 1])^2+det([x(F(i,1)) y(F(i,1)) 1; x(F(i,2)) y(F(i,2)) 1; x(F(i,3)) y(F(i,3)) 1])^2);
    PF=[PF; [app (x(F(i,1))+x(F(i,2))+x(F(i,3)))/3 (y(F(i,1))+y(F(i,2))+y(F(i,3)))/3 (z(F(i,1))+z(F(i,2))+z(F(i,3)))/3]];
end
    
    