function G = G(d)

%speed of light
c = 299792.458;
%Locations
A1=15600; A2=18760; A3=17610; A4=19170;
B1=7540; B2=2750; B3=14630; B4=610;
C1=20140; C2=18610; C3=13480; C4=18390;

%Times
t1=0.07074; t2=0.07220; t3=0.07690; t4=0.07242;

%System of linear equations
ux=[-2*A1+2*A2; -2*A1+2*A3; -2*A1+2*A4];
uy=[-2*B1+2*B2; -2*B1+2*B3; -2*B1+2*B4];
uz=[-2*C1+2*C2; -2*C1+2*C3; -2*C1+2*C4];
ud=[2*c^2*t1-2*c^2*t2; 2*c^2*t1-2*c^2*t3; 2*c^2*t1-2*c^2*t4];
w=[A1^2-A2^2+B1^2-B2^2+C1^2-C2^2-(c^2)*(t1^2)+(c^2)*t2^2;
    A1^2-A3^2+B1^2-B3^2+C1^2-C3^2-(c^2)*t1^2+(c^2)*t3^2;
    A1^2-A4^2+B1^2-B4^2+C1^2-C4^2-(c^2)*t1^2+(c^2)*(t4^2)];

yzx = det([uy, uz, ux]);
yzd = det([uy, uz, ud]);
yzw = det([uy, uz, w]);

xzy = det([ux, uz, uy]);
xzd = det([ux, uz, ud]);
xzw = det([ux, uz, w]);

xyz = det([ux, uy, uz ]);
xyd = det([ux, uy, ud]);
xyw = det([ux, uy, w]);

G = ((d * yzd + yzw)/yzx - A1)^2 + ((d * xzd + xzw)/xzy - B1)^2 + ...
     ((d * xyd + xyw)/xyz - C1)^2 - c^2 * (t1 - d)^2;
 