GPS Software Receivers - Fundamentals of Global Positioning System Receivers: A Software Approach - page 217

error=1;
while error>1e-12;
E=M+e*sin(Eold);
error=abs(E-Eold);
Eold=E;
end
F=-4.442807633e-10;
deltr=F*e*asq*sin(E);
delt=af0+af1*(tc-toc)+af2*(tc-toc)^2+deltr-tgd;
t=tc-delt;
v1=acos((cos(E)-e)/(1-e*cos(E)));
v2=asin(((1-e^2)^.5)*sin(E)/(1-e*cos(E)));
v=v1*sign(v2);
%r=a*(1-e^2)/(1+e*cos(v));
phi=v+w;
delphi=cus*sin(2*phi)+cuc*cos(2*phi);
delr=crc*cos(2*phi)+crs*sin(2*phi);
deli=cic*cos(2*phi)+cis*sin(2*phi);
phi=phi+delphi;
r=a*(1-e*cos(E));
r=r+delr;
i=i+deli+idot*(t-toe);
omeger=ra+delra*(t-toe)-wer*t;
x=r*cos(phi)*cos(omeger)-r*sin(phi)*cos(i)*sin(omeger);
y=r*cos(phi)*sin(omeger)+r*sin(phi)*cos(i)*cos(omeger);
z=r*sin(phi)*sin(i);
PRf=PRc+delt*299792458;
tf=t;
outp=[x y z PRf tf];