function f = gaussnewton(x,y,z,d)

% c is the speed of light in km, p is rho, and z1 is the radius of the
% earth
c = 299792.458;
p = 26570;
z1 = 6370;

% Initial guess
x0 = [x; y; z; d];

% Setting Spherical Coordinates
phi1 = 3*pi / 8;
theta1 = pi / 8;

phi2 = pi / 4;
theta2 = pi / 3;

phi3 = 3*pi / 8;
theta3 = 2*pi / 3;

phi4 = pi / 4;
theta4 = pi / 4;

phi5 = pi / 6;
theta5 = 3*pi / 8;

phi6 = 3*pi / 7;
theta6 = 5*pi/7;

phi7= pi / 2;
theta7 = pi/10;

% Define spheres
A1 = p * cos(phi1) * cos(theta1);
B1 = p * cos(phi1) * sin(theta1);
C1 = p * sin(phi1);
t1 = (d + sqrt (A1^2 + B1^2 + (C1 - z1)^2) / c) -10^-8;

A2 = p * cos(phi2) * cos(theta2);
B2 = p * cos(phi2) * sin(theta2);
C2 = p * sin(phi2);
t2 = (d + sqrt (A2^2 + B2^2 + (C2 - z1)^2) / c) + 10^-8;

A3 = p * cos(phi3) * cos(theta3);
B3 = p * cos(phi3) * sin(theta3);
C3 = p * sin(phi3);
t3 = (d + sqrt (A3^2 + B3^2 + (C3 - z1)^2) / c)-10^-8;

A4 = p * cos(phi4) * cos(theta4);
B4 = p * cos(phi4) * sin(theta4);
C4 = p * sin(phi4);
t4 = (d + sqrt (A4^2 + B4^2 + (C4 - z1)^2) / c)+ 10^-8;

A5 = p * cos(phi5) * cos(theta5);
B5 = p * cos(phi5) * sin(theta5);
C5 = p * sin(phi5);
t5 = (d + sqrt (A5^2 + B5^2 + (C5 - z1)^2) / c)-10^-8 ;

A6 = p * cos(phi6) * cos(theta6);
B6 = p * cos(phi6) * sin(theta6);
C6 = p * sin(phi6);
t6 = (d + sqrt (A6^2 + B6^2 + (C6 - z1)^2) / c)+ 10^-8;

A7 = p * cos(phi7) * cos(theta7);
B7 = p * cos(phi7) * sin(theta7);
C7 = p * sin(phi7);
t7 = (d + sqrt (A7^2 + B7^2 + (C7 - z1)^2) / c)+ 10^-8;

% Equations for the eight satelittes
F1 = (x - A1)^2 + (y - B1)^2 + (z - C1)^2 - (c * (t1 - d))^2;
F2 = (x - A2)^2 + (y - B2)^2 + (z - C2)^2 - (c * (t2 - d))^2;
F3 = (x - A3)^2 + (y - B3)^2 + (z - C3)^2 - (c * (t3 - d))^2;
F4 = (x - A4)^2 + (y - B4)^2 + (z - C4)^2 - (c * (t4 - d))^2;
F5 = (x - A5)^2 + (y - B5)^2 + (z - C5)^2 - (c * (t5 - d))^2;
F6 = (x - A6)^2 + (y - B6)^2 + (z - C6)^2 - (c * (t6 - d))^2;
F7 = (x - A7)^2 + (y - B7)^2 + (z - C7)^2 - (c * (t7 - d))^2;

% Matrix of satelitte equations
 F = [F1;F2;F3;F4;F5;F6;F7];

% These equations describe the Jacobian, a matrix of partial derivatives
DF1x = 2 * (x - A1);
DF1y = 2 * (y - B1);
DF1z = 2 * (z - C1);
DF1d = 2 * c^2 * (t1 - d);

DF2x = 2 * (x - A2);
DF2y = 2 * (y - B2);
DF2z = 2 * (z - C2);
DF2d = 2 * c^2 * (t2 - d);

DF3x = 2 * (x - A3);
DF3y = 2 * (y - B3);
DF3z = 2 * (z - C3);
DF3d = 2 * c^2 * (t3 - d);

DF4x = 2 * (x - A4);
DF4y = 2 * (y - B4);
DF4z = 2 * (z - C4);
DF4d = 2 * c^2 * (t4 - d);

DF5x = 2 * (x - A5);
DF5y = 2 * (y - B5);
DF5z = 2 * (z - C5);
DF5d = 2 * c^2 * (t5 - d);

DF6x = 2 * (x - A6);
DF6y = 2 * (y - B6);
DF6z = 2 * (z - C6);
DF6d = 2 * c^2 * (t6 - d);

DF7x = 2 * (x - A7);
DF7y = 2 * (y - B7);
DF7z = 2 * (z - C7);
DF7d = 2 * c^2 * (t7 - d);

% DF represents the Jacobian Matrix
DF = [DF1x DF1y DF1z DF1d;
 DF2x DF2y DF2z DF2d;
 DF3x DF3y DF3z DF3d;
 DF4x DF4y DF4z DF4d;
 DF5x DF5y DF5z DF5d;
 DF6x DF6y DF6z DF6d;
 DF7x DF7y DF7z DF7d];

% Perform Guass-Newton Iteration
v = (DF'*DF) \ (-DF'*F);

% update initial guess
f = x0 + v;
end