% La funzione eq_retta restituisce in uscita le coordinate della retta
% passante per i due punti in ingresso come l'intersezione di due piani

function [a b c d] = eq_retta(P_1, P_2)

a = zeros(2,1);
b = zeros(2,1);
c = zeros(2,1);
d = zeros(2,1);

% per il primo piano

P_3 = [0 0 0];

a(1) = (P_2(2)-P_1(2))*(P_3(3)-P_1(3)) - (P_2(3)-P_1(3))*(P_3(2)-P_1(2));

b(1) = -(P_2(1)-P_1(1))*(P_3(3)-P_1(3)) + (P_2(3)-P_1(3))*(P_3(1)-P_1(1));

c(1) = (P_2(1)-P_1(1))*(P_3(2)-P_1(2)) - (P_2(2)-P_1(2))*(P_3(1)-P_1(1));

d(1) = -a(1)*P_1(1) - b(1)*P_1(2) - c(1)*P_1(3);

% per il secondo piano

P_3 = [1000 1000 0];

a(2) = (P_2(2)-P_1(2))*(P_3(3)-P_1(3)) - (P_2(3)-P_1(3))*(P_3(2)-P_1(2));

b(2) = -(P_2(1)-P_1(1))*(P_3(3)-P_1(3)) + (P_2(3)-P_1(3))*(P_3(1)-P_1(1));

c(2) = (P_2(1)-P_1(1))*(P_3(2)-P_1(2)) - (P_2(2)-P_1(2))*(P_3(1)-P_1(1));

d(2) = -a(2)*P_1(1) - b(2)*P_1(2) - c(2)*P_1(3);
