function res = gpsfun(r,s)
% gps(r, s)
% r -- receiver position [x,y,z,b]
% s --- satellite positions in matrix format
% [x1,y1,z1,dt1; 
%  x2,y2,z2,dt2; 
%  x3,y3,z3,dt3;
%  x4,y4,z4,dt4]
% units: meters and nanosecs
res = zeros(4,1);
c2 = (299792458e-9).^2; %! Speed of light (m/usec)
for k=1:4 %! One function for each satellite
    res(k) = (r(1) - s(k,1)).^2 + (r(2) - s(k,2)).^2 ...
        + (r(3) - s(k,3)).^2 - c2.*(r(4) - s(k,4)).^2;
end