function [a,b,c,d]=get_plan_meth_3(X,Y,Z)

  cc = combnk(1:4,3);      
      
    for n=1:size(cc,1)

         pos=cc(n,:);
         x1= X(pos(1)); y1= Y(pos(1));
         x2= X(pos(2)); y2= Y(pos(2));
         x3= X(pos(3)); y3= Y(pos(3));

         % Os Z

         z1= Z(pos(1)); 
         z2= Z(pos(2));
         z3= Z(pos(3));
         
         % Plano

         A=  (y1*z2 - y2*z1 - y1*z3 + y3*z1 + y2*z3 - y3*z2)/(x1*y2*z3 - x1*y3*z2 - x2*y1*z3 + x2*y3*z1 + x3*y1*z2 - x3*y2*z1);
         B= -(x1*z2 - x2*z1 - x1*z3 + x3*z1 + x2*z3 - x3*z2)/(x1*y2*z3 - x1*y3*z2 - x2*y1*z3 + x2*y3*z1 + x3*y1*z2 - x3*y2*z1);
         C=  (x1*y2 - x2*y1 - x1*y3 + x3*y1 + x2*y3 - x3*y2)/(x1*y2*z3 - x1*y3*z2 - x2*y1*z3 + x2*y3*z1 + x3*y1*z2 - x3*y2*z1);
         
         
         % D calculator
         % a/d x + b/d y + c/d z == 1
         % A x + B y + C z == 1
         
         Zpt=A/C*x3+B/C*y3;
         D=z3+Zpt;
         
         
         % Desenhar o plano
 
         x=-45:15:45;
        [Xp,Yp] = meshgrid(x);
          a=A*D; b=B*D; c=C*D; d=D;
         % a x + b y + c z == d
        
        Zp=-a/c*Xp-b/c*Yp+d/c;
        
%          [u,v,w]=surfnorm(Xp,Yp,Zp);
%         surf(Xp,Yp,Zp)
%         colormap([1,1,0])
%         quiver3(Xp,Yp,Zp,-u,-v,-w)
%           line([0 a],[0 b],[0 c])
%         xlabel('x'); ylabel('y'); zlabel('z')


         % Componentes da normal

        nx=a;
        ny=b;
        nz=c;

        n=[u(1) v(1) w(1)];

        % Parte do ângulo (Produto Escalar com eixos referencial???)
        
        lx=d*c/a;
        ly=d*c/b;
        
        pit=rad2deg(atan(d/lx));
        rol=rad2deg(atan(d/ly));
        
        PitchA=[PitchA pit];
        RollA= [RollA  rol];
        nn=[nn ;n];
        dd=[dd dot(n,[0 0 1])];

    end
    %
        Pitch=[Pitch ; file_s(i,1) mean(PitchA)];
        Roll=[Roll ; file_s(i,1) mean(RollA)];




end