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];

% Define phi and theta for spherical coordinates
phi1 = pi / 8;
theta1 = pi / 8;

phi2 = pi / 7;
theta2 = pi / 3;

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

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

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

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

phi7 = pi / 5;
theta7 = 2*pi / 9;

phi8 = pi / 10;
theta8 = 2*pi / 5;

phi9 = 2*pi / 11;
theta9 = 5*pi / 10;

phi10 = 3*pi / 8;
theta10 = 13*pi / 18;

phi11 = 8*pi / 18;
theta11 = 2*pi / 3;


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

A8 = p * cos(phi8) * cos(theta8);
B8 = p * cos(phi8) * sin(theta8);
C8 = p * sin(phi8);
t8 = (d + sqrt(A8^2 + B8^2 + (C8 - z1)^2) / c)+10^-8;

A9 = p * cos(phi9) * cos(theta9);
B9 = p * cos(phi9) * sin(theta9);
C9 = p * sin(phi9);
t9 = (d + sqrt(A9^2 + B9^2 + (C9 - z1)^2) / c)-10^-8;

A10 = p * cos(phi10) * cos(theta10);
B10 = p * cos(phi10) * sin(theta10);
C10 = p * sin(phi10);
t10 = (d + sqrt(A10^2 + B10^2 + (C10 - z1)^2) / c)+ 10^-8;

A11 = p * cos(phi11) * cos(theta11);
B11 = p * cos(phi11) * sin(theta11);
C11 = p * sin(phi11);
t11 = (d + sqrt(A11^2 + B11^2 + (C11 - z1)^2) / c)+ 10^-8;

% Equations for the eleven 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;
F8 = (x - A8)^2 + (y - B8)^2 + (z - C8)^2 - (c * (t8 - d))^2;
F9 = (x - A9)^2 + (y - B9)^2 + (z - C9)^2 - (c * (t9 - d))^2;
F10 = (x - A10)^2 + (y - B10)^2 + (z - C10)^2 - (c * (t10 - d))^2;
F11 = (x - A11)^2 + (y - B11)^2 + (z - C11)^2 - (c * (t11 - d))^2;

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

% These equations describe the Jacobian, a matrix of partial derivatives
F1x = 2 * (x - A1);
F1y = 2 * (y - B1);
F1z = 2 * (z - C1);
F1d = 2 * c^2 * (t1 - d);

F2x = 2 * (x - A2);
F2y = 2 * (y - B2);
F2z = 2 * (z - C2);
F2d = 2 * c^2 * (t2 - d);

F3x = 2 * (x - A3);
F3y = 2 * (y - B3);
F3z = 2 * (z - C3);
F3d = 2 * c^2 * (t3 - d);

F4x = 2 * (x - A4);
F4y = 2 * (y - B4);
F4z = 2 * (z - C4);
F4d = 2 * c^2 * (t4 - d);

F5x = 2 * (x - A5);
F5y = 2 * (y - B5);
F5z = 2 * (z - C5);
F5d = 2 * c^2 * (t5 - d);

F6x = 2 * (x - A6);
F6y = 2 * (y - B6);
F6z = 2 * (z - C6);
F6d = 2 * c^2 * (t6 - d);

F7x = 2 * (x - A7);
F7y = 2 * (y - B7);
F7z = 2 * (z - C7);
F7d = 2 * c^2 * (t7 - d);

F8x = 2 * (x - A8);
F8y = 2 * (y - B8);
F8z = 2 * (z - C8);
F8d = 2 * c^2 * (t8 - d);

F9x = 2 * (x - A9);
F9y = 2 * (y - B9);
F9z = 2 * (z - C9);
F9d = 2 * c^2 * (t9 - d);

F10x = 2 * (x - A10);
F10y = 2 * (y - B10);
F10z = 2 * (z - C10);
F10d = 2 * c^2 * (t10 - d);

F11x = 2 * (x - A11);
F11y = 2 * (y - B11);
F11z = 2 * (z - C11);
F11d = 2 * c^2 * (t11 - d);



% DH represents the Jacobian Matrix
DH = [F1x F1y F1z F1d; F2x F2y F2z F2d; F3x F3y F3z F3d; F4x F4y F4z F4d; F5x F5y F5z F5d; F6x F6y F6z F6d; F7x F7y F7z F7d; F8x F8y F8z F8d;F9x F9y F9z F9d; F10x F10y F10z F10d; F11x F11y F11z F11d];

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

% update initial guess
f = x0 + v;
end

Error using gaussnewton (line 11)
Not enough input arguments.

Published with MATLAB® R2013a