Вычисление координат пользователя GPS

clear;
clc;
%% ДАНО
%Координаты спутников GPS
S(1,1) = 1.65909e+07; 
S(1,2) = -1.42984e+07;
S(1,3) = 1.49624e+07;
S(2,1) = 5.17469e+06;
S(2,2) = -1.64476e+07;
S(2,3) = 2.01232e+07;
S(3,1)= 2.10854e+07;
S(3,2)= 7.19856e+06;
S(3,3)= 1.47777e+07;
S(4,1) = 2.59733e+07;
S(4,2) = 3.08415e+06;
S(4,3) = 6.37287e+06;
S(5,1) = 1.80273e+07;
S(5,2) = 109452;
S(5,3) = 1.93001e+07;
S(6,1)= 9.05944e+06;
S(6,2)= 1.33543e+07;
S(6,3)= 2.09274e+07;
S(7,1) = -5.52345e+06;
S(7,2) = 1.55533e+07;
S(7,3) = 2.06664e+07;
S(8,1)= 2.07086e+06;
S(8,2)= 2.23079e+07;
S(8,3)= 1.40383e+07;
%Псвдодальности 
P(1)= 22911434.477 ;
P(2)=23424652.828  ;
P(3)=21130277.828 ;
P(4)= 23222461.156  ;
P(5)=20621254.992  ;
P(6)=20328694.563  ;
P(7)=22208818.000  ;
P(8)=22344913.023  ;
%%
%число спутников
n=8;
%Начало отсчета для решения задачи расчета координат позиции приемника (центр земли)
x=0; y=0; z=0; t=0;
t=0;
% x=2825879.2;
% y=1641765.1;
% z=5471439.7;
a=[x y z t];%предварительные координаты пункта наблюдения
eps=1;
%задание точности вычислений
dx=eps*10;dy=eps*10;dz=eps*10;dt=eps*10;
da=[10*eps 10*eps 10*eps 10*eps];
while da>[eps eps eps eps]
for i=1:n
R=[x y z t];
%Расчетные псевдодальности (геометрическая дальность)
RR(i)=sqrt((S(i,1)-R(1,1))^2+(S(i,2)-R(1,2))^2+(S(i,3)-R(1,3))^2);
for j=1:3
A(i,j)=(S(i,j)-R(j))/RR(i); %направляющие косинусы топоцентрических направлений на спутник
end;
A(1:n,4)=1;
l(i)=RR(i)-P(i);
end
N=A'*A;
L=A'*l';
% X=inv(N)*L;
X=N\L;
R=R+X';
x=R(1);y=R(2);z=R(3);
% t=R(4);
X(4)=R(4);
dx=abs(R(1)-a(1));
dy=abs(R(2)-a(2));
dz=abs(R(3)-a(3));
dt=abs(R(4)-a(4));
da=[dx dy dz dt ];
a=R;
Q=inv(N);
PDOP=sqrt(Q(1,1)+Q(2,2)+Q(3,3));
end
R