function r = gauss_newton_p6(num_iterations, change_time, dt)

p = 26570; % (km)
d = 0.0001;
c = 299792.458; % (km)

% 0 < phi < pi/2 (unbunched)
phi_1 = 13.5*(pi/180);
phi_2 = 38.9*(pi/180);
phi_3 = 65.3*(pi/180);
phi_4 = 87.2*(pi/180);
phi_5 = 2.5*(pi/180);    % 5 satellites added
phi_6 = 23.4*(pi/180);
phi_7 = 49.2*(pi/180);
phi_8 = 57.7*(pi/180);
phi_9 = 75.2*(pi/180);

% 0 < theta < 2*pi (unbunched)
theta_1 = 29.5*(pi/180);
theta_2 = 103.7*(pi/180);
theta_3 = 213.5*(pi/180);
theta_4 = 322.6*(pi/180);
theta_5 = 78.3*(pi/180);   % 5 satellites added
theta_6 = 155.8*(pi/180);
theta_7 = 279.7*(pi/180);
theta_8 = 300.8*(pi/180);
theta_9 = 355.3*(pi/180);

A1 = p*cos(phi_1)*cos(theta_1);  A2 = p*cos(phi_2)*cos(theta_2);
B1 = p*cos(phi_1)*sin(theta_1);  B2 = p*cos(phi_2)*sin(theta_2);
C1 = p*sin(phi_1);               C2 = p*sin(phi_2);
A3 = p*cos(phi_3)*cos(theta_3);  A4 = p*cos(phi_4)*cos(theta_4);
B3 = p*cos(phi_3)*sin(theta_3);  B4 = p*cos(phi_4)*sin(theta_4);
C3 = p*sin(phi_3);               C4 = p*sin(phi_4);

% 5 satellites added
A5 = p*cos(phi_5)*cos(theta_5);  A6 = p*cos(phi_6)*cos(theta_6);
B5 = p*cos(phi_5)*sin(theta_5);  B6 = p*cos(phi_6)*sin(theta_6);
C5 = p*sin(phi_5);               C6 = p*sin(phi_6);
A7 = p*cos(phi_7)*cos(theta_7);  A8 = p*cos(phi_8)*cos(theta_8);
B7 = p*cos(phi_7)*sin(theta_7);  B8 = p*cos(phi_8)*sin(theta_8);
C7 = p*sin(phi_7);               C8 = p*sin(phi_8);
A9 = p*cos(phi_9)*cos(theta_9);
B9 = p*cos(phi_9)*sin(theta_9);
C9 = p*sin(phi_9);

% Ranges
R1 = sqrt(A1^2 + B1^2 + (C1 - 6370)^2);
R2 = sqrt(A2^2 + B2^2 + (C2 - 6370)^2);
R3 = sqrt(A3^2 + B3^2 + (C3 - 6370)^2);
R4 = sqrt(A4^2 + B4^2 + (C4 - 6370)^2);
R5 = sqrt(A5^2 + B5^2 + (C5 - 6370)^2);
R6 = sqrt(A6^2 + B6^2 + (C6 - 6370)^2);
R7 = sqrt(A7^2 + B7^2 + (C7 - 6370)^2);
R8 = sqrt(A8^2 + B8^2 + (C8 - 6370)^2);
R9 = sqrt(A9^2 + B9^2 + (C9 - 6370)^2);

% Travel Times
if (change_time == 0)
  T1 = d + R1/c;
  T2 = d + R2/c;
  T3 = d + R3/c;
  T4 = d + R4/c;
  T5 = d + R5/c;
  T6 = d + R6/c;
  T7 = d + R7/c;
  T8 = d + R8/c;
  T9 = d + R9/c;
else
  T1 = d + R1/c + dt(1);
  T2 = d + R2/c + dt(2);
  T3 = d + R3/c + dt(3);
  T4 = d + R4/c + dt(4);
  T5 = d + R5/c + dt(5);
  T6 = d + R6/c + dt(6);
  T7 = d + R7/c + dt(7);
  T8 = d + R8/c + dt(8);
  T9 = d + R9/c + dt(9);
end

x = [0 0 6370 d]';  % initial guess vector

F = @(x) [(x(1) - A1)^2 + (x(2) - B1)^2 + (x(3) - C1)^2 - c^2*(x(4) - T1)^2;
           (x(1) - A2)^2 + (x(2) - B2)^2 + (x(3) - C2)^2 - c^2*(x(4) - T2)^2;
           (x(1) - A3)^2 + (x(2) - B3)^2 + (x(3) - C3)^2 - c^2*(x(4) - T3)^2;
           (x(1) - A4)^2 + (x(2) - B4)^2 + (x(3) - C4)^2 - c^2*(x(4) - T4)^2;
           (x(1) - A5)^2 + (x(2) - B5)^2 + (x(3) - C5)^2 - c^2*(x(4) - T5)^2;
           (x(1) - A6)^2 + (x(2) - B6)^2 + (x(3) - C6)^2 - c^2*(x(4) - T6)^2;
           (x(1) - A7)^2 + (x(2) - B7)^2 + (x(3) - C7)^2 - c^2*(x(4) - T7)^2;
           (x(1) - A8)^2 + (x(2) - B8)^2 + (x(3) - C8)^2 - c^2*(x(4) - T8)^2;
           (x(1) - A9)^2 + (x(2) - B9)^2 + (x(3) - C9)^2 - c^2*(x(4) - T9)^2];

DF = @(x) [2*(x(1) - A1) 2*(x(2) - B1) 2*(x(3) - C1) 2*c^2*(-x(4) + T1);
           2*(x(1) - A2) 2*(x(2) - B2) 2*(x(3) - C2) 2*c^2*(-x(4) + T2);
           2*(x(1) - A3) 2*(x(2) - B3) 2*(x(3) - C3) 2*c^2*(-x(4) + T3);
           2*(x(1) - A4) 2*(x(2) - B4) 2*(x(3) - C4) 2*c^2*(-x(4) + T4);
           2*(x(1) - A5) 2*(x(2) - B5) 2*(x(3) - C5) 2*c^2*(-x(4) + T5);
           2*(x(1) - A6) 2*(x(2) - B6) 2*(x(3) - C6) 2*c^2*(-x(4) + T6);
           2*(x(1) - A7) 2*(x(2) - B7) 2*(x(3) - C7) 2*c^2*(-x(4) + T7);
           2*(x(1) - A8) 2*(x(2) - B8) 2*(x(3) - C8) 2*c^2*(-x(4) + T8);
           2*(x(1) - A9) 2*(x(2) - B9) 2*(x(3) - C9) 2*c^2*(-x(4) + T9)];

for i = 1:num_iterations
   v = DF(x)\F(x);
   x = x - v;
end

r = x;