%%%MAIN PROGRAM%%%%%
% This program allows the user to find integer face areas and integer principal edges' length for a right
% tetrahedron from a given set of Pythagorean Quadruples.
% Program Restrictions:
% %% This code doesn't check if the user supplied values are integer or numeric.
% But the program does verify if the supplied set of numbers are Pythgorean
% Quadruple.
% %% Only one set of numbers can be evaluated when the program is executed.
% 
% Execute this Script once you have a Pythagorean Quadruple set.
% Example: (7, 4, 4, 9), (1, 2, 2, 3)
import math.*;
clc;
%Ask for the four numbers (A, B, C, D)
disp('Please enter valid Pythagorean quadruple set to compute the Right Tehrahedron values.');
valueOne = input('Enter the first number.');
valueTwo = input('Enter the second number.');
valueThree = input('Enter the third number.');
valueFour = input('Enter the fourth number.');
%check if the given four numbers are valid pythagorean quadruples.
result = checkPQ(valueOne, valueTwo, valueThree, valueFour);

valueK = 1;
%if the provided numbers are pythagorean, we get our values for A, B, C and
%D and then find the k value from A, B, C.
if (result(1) == 1)
    valueD = result(2);
    valueA = result(3);
    valueB = result(4);
    valueC = result(5);
    %get the remaining square free divisor from A, B and C to calculate k
    %value
    valueK = 2 * findDiv(valueA) * findDiv(valueB) * findDiv(valueC);
    findRTvalues(valueK, valueA, valueB, valueC, valueD);
else
    disp('The given set of the number is not a Pythagorean Quadruple. Please run the program with a valid set.');
end

%%%%%%%%%%%%%%%%%%%%%%
%%FUNCTIONS%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%

function result = checkPQ(one, two, three, four)
% This functions checks if the given set of numbers are 
% Pythagorean quadruple
numbers = [one two three four];

if (one == max(numbers))
    D = one;
    numbers = [D two three four];
    value = (D^2 == two^2 + three^2 + four^2);
elseif (two == max(numbers))
    D = two;
    numbers = [D one three four];
    value = (D^2 == one^2 + three^2 + four^2);
elseif (three == max(numbers))
    D = three;
    numbers = [D two one four];
    value = (D^2 == two^2 + one^2 + four^2);
else
    D = four;
    numbers = [D two three one];
    value = (D^2 == two^2 + three^2 + one^2);
end
result = [value numbers];
end


function remDiv = findDiv (numA)
%This function allows us to find the square free remaining divisor of a
%given number.
tmpA = numA;
sqA = sqrt(tmpA); 
limit = floor(sqA);
%sqFactor = []; %all perfect square divisor
if (limit == 1)
    remDiv = tmpA;
elseif (limit == sqA)
    remDiv = 1;
else
    for i=limit:-1:1
        %do sth
        divA = tmpA/i;
        if (divA == floor(divA))
            divA = divA/i;
            if (divA == floor(divA))
                tmpA = divA;
                %this is a perfect square factor
            else
                divA = tmpA;
            end
        else
            divA = tmpA;
        end
    end
    remDiv = divA;
end
end

function finalResult = findRTvalues (numK, numA, numB, numC, numD)
%This function is used once we find our k values to find a right
%tetrahedron's principal edges' length and face areas from the given set of
%Pythagorean Quadruples
finalResult = 0;
edgeA = sqrt((2*numK*numB*numC)/numA);
edgeB = sqrt((2*numK*numA*numC)/numB);
edgeC = sqrt((2*numK*numB*numA)/numC);

areaA = numK * numA;
areaB = numK * numB;
areaC = numK * numC;
areaD = numK * numD;

%Print the results
fprintf('For the given set of Pythagorean Qaudruple, A = %d B = %d C = %d and D = %d, we find k = %d.', numA, numB, numC, numD, numK);
fprintf(newline);
fprintf('Using those value, a Right Tetrahedron with the following values can be obtained:');
fprintf(newline);
fprintf('Face Area (kA): %d', areaA);
fprintf(newline);
fprintf('Face Area (kB): %d', areaB);
fprintf(newline);
fprintf('Face Area (kC): %d', areaC);
fprintf(newline);
fprintf('Face Area (kD): %d', areaD);
fprintf(newline);
fprintf('Principal edge (a): %d', edgeA);
fprintf(newline);
fprintf('Principal edge (b): %d', edgeB);
fprintf(newline);
fprintf('Principal edge (c): %d', edgeC);
fprintf(newline);
end