function position = SolveGPSQuadratic(satellites, x)
speed_of_light = 299792.458;
speed_of_light_squared = speed_of_light^2;
a = 0;
b = 0;
c = 0;
for i = 1 : 3
    a = a + x(i, 1)^2;
    constant = x(i, 2) - satellites(1, i);
    b = b + 2*x(i, 1)*constant;
    c = c + constant^2;
end
a = a - speed_of_light_squared;
b = b + 2*speed_of_light_squared*satellites(1, 4);
c = c - speed_of_light_squared*satellites(1, 4)^2;

position = zeros(2, 4);
position(1, 4) = (-b + sqrt(b^2 - 4*a*c))/(2*a);
position(2, 4) = (-b - sqrt(b^2 - 4*a*c))/(2*a);

for i = 1 : 2
    for j = 1 : 3
        position(i, j) = x(j, 1)*position(i, 4) + x(j, 2);
    end
end

