function r = gauss_newton_p4(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);

% 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);

A1 = p*cos(phi_1)*cos(theta_1);
B1 = p*cos(phi_1)*sin(theta_1);
C1 = p*sin(phi_1);

A2 = p*cos(phi_2)*cos(theta_2);
B2 = p*cos(phi_2)*sin(theta_2);
C2 = p*sin(phi_2);

A3 = p*cos(phi_3)*cos(theta_3);
B3 = p*cos(phi_3)*sin(theta_3);
C3 = p*sin(phi_3);

A4 = p*cos(phi_4)*cos(theta_4);
B4 = p*cos(phi_4)*sin(theta_4);
C4 = p*sin(phi_4);

% 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);

% Travel Times
if (change_time == 0)
  T1 = d + R1/c;
  T2 = d + R2/c;
  T3 = d + R3/c;
  T4 = d + R4/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);
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];
       
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)];

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

r = x;