function problem1(solution_selection)

A1 = 15600;B1 = 7540; C1 = 20140;T1 = 0.07074;
A2 = 18760;B2 = 2750; C2 = 18610;T2 = 0.07220;
A3 = 17610;B3 = 14630;C3 = 13480;T3 = 0.07690;
A4 = 19170;B4 = 610;  C4 = 18390;T4 = 0.07242;
c = 299792.458;

ux = [2*A2 - 2*A1;
      2*A3 - 2*A1;
      2*A4 - 2*A1];
  
uy = [2*B2 - 2*B1;
      2*B3 - 2*B1;
      2*B4 - 2*B1];
  
uz = [2*C2 - 2*C1;
      2*C3 - 2*C1;
      2*C4 - 2*C1];
  
ud = [2*T1*c^2 - 2*T2*c^2;
      2*T1*c^2 - 2*T3*c^2;
      2*T1*c^2 - 2*T4*c^2];
  
w =  [A1^2 - A2^2 + B1^2 - B2^2 + C1^2 - C2^2 - c^2*T1^2 + c^2*T2^2;
      A1^2 - A3^2 + B1^2 - B3^2 + C1^2 - C3^2 - c^2*T1^2 + c^2*T3^2;
      A1^2 - A4^2 + B1^2 - B4^2 + C1^2 - C4^2 - c^2*T1^2 + c^2*T4^2];
  
CC1 = -det(horzcat(uy,uz,ud))/det(horzcat(uy,uz,ux));
CC2 = -det(horzcat(uy,uz,w))/det(horzcat(uy,uz,ux));

CC3 = -det(horzcat(uz,ux,ud))/det(horzcat(uz,ux,uy));
CC4 = -det(horzcat(uz,ux,w))/det(horzcat(uz,ux,uy));

CC5 = -det(horzcat(ux,uy,ud))/det(horzcat(ux,uy,uz));
CC6 = -det(horzcat(ux,uy,w))/det(horzcat(ux,uy,uz));

%p1 = CC1^2*d^2 +  2*d*(CC1*CC2 - CC1*A1) + CC2^2 - 2*A1*CC2 + A1^2
%p2 = CC3^2*d^2 +  2*d*(CC3*CC4 - CC3*B1) + CC4^2 - 2*B1*CC4 + B1^2
%p3 = CC5^2*d^2 +  2*d*(CC5*CC6 - CC5*C1) + CC6^2 - 2*C1*CC6 + C1^2
%p4 = -c^2*d^2 + 2*c^2*d*T1 - c^2*T1^2

%Setup for Quadratic Equation
%a: d^2*(CC1^2 + CC3^2 + CC5^2 - c^2)
%b: d*2*((CC1*CC2 - CC1*A1) + (CC3*CC4 - CC3*B1) + (CC5*CC6 - CC5*C1) + c^2*T1)
%c: CC2^2 - 2*A1*CC2 + A1^2 + CC4^2 - 2*B1*CC4 + B1^2 + CC6^2 - 2*C1*CC6 + C1^2 - c^2*T1^2

%Use Quadratic Formula to Find d
a = (CC1^2 + CC3^2 + CC5^2 - c^2);
b = 2*((CC1*CC2 - CC1*A1) + (CC3*CC4 - CC3*B1) + (CC5*CC6 - CC5*C1) + c^2*T1);
c = CC2^2 - 2*A1*CC2 + A1^2 + CC4^2 - 2*B1*CC4 + B1^2 + CC6^2 - 2*C1*CC6 + C1^2 - c^2*T1^2;
discr = b^2 - 4*a*c;

if (solution_selection == 0)
   d = (-b + sqrt(b^2 - 4*a*c))/(2*a)
else
   d = (-b - sqrt(b^2 - 4*a*c))/(2*a)
end

%Use d to find x, y, z
x = CC1*d + CC2
y = CC3*d + CC4
z = CC5*d + CC6
