function J = SatellitesJacobianMatrix(satellites, position)
speed_of_light = 299792.458;
speed_of_light_squared = speed_of_light^2;
[rows, columns] = size(satellites);
J = zeros(rows, 4);
for i = 1 : rows
    for j = 1 : 3
        J(i, j) = 2*(position(j) - satellites(i, j));
    end
    J(i, 4) = -2*(speed_of_light_squared^2)*(satellites(i, 4) - position(4));
end

