function [Ph]=get_plan_meth_2(X,Y,Z)


cc = combnk(1:4,3);                                 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some variables to store data %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Pa=[];
Pb=[];
Pc=[];
Pd=[];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%           Method 2           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for n=1:size(cc,1)

    pos=cc(n,:);
    x1= X(pos(1)); y1= Y(pos(1));z1= Z(pos(1)); 
    x2= X(pos(2)); y2= Y(pos(2));z2= Z(pos(2));
    x3= X(pos(3)); y3= Y(pos(3));z3= Z(pos(3));
     
    xf=[x1 x2 x3];
    yf=[y1 y2 y3];
    zf=[z1 z2 z3];

    %%%%%%%%
    % Plan %
    %%%%%%%%

    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;

    a=A*D; b=B*D; c=C*D; d=D;

    % A B C and D.
    Pa=[Pa a]; Pb=[Pb b]; Pc=[Pc c]; Pd=[Pd d];
end

% Mean of A, B, C and D of each plan.    
aa=mean(Pa); bb=mean(Pb); cc=mean(Pc); dd=mean(Pd);

A=aa;
B=bb;
C=cc;
D=dd;

Ph=[A B C D];

end