The basic equations for the position of a GPS satellite at time \(t_k\) are fairly simple, however calculating the elements contained in these equations involves a number of steps.
First the Receiver needs to download and transform the Ephemeris sub-frames in the data stream transmitted by each satellite.
Then the values in the Ephemeris, supplied by the satellite, are used in the following equations to calculate the position of the satellite at a given time. Most of the values used in these equations are supplied in the Ephemeris sub-frames and only the equation(8) for eccentric needs to be solved.
When determining the correct time to use to calculate the postition of the satellite three corrections need to be made.
Correct for aging of satellite onboard high accuracy clock.Then once the satellites position has been calculated a correction needs to be applied for the amount the earth has rotated during the time it took for the RF signal to travel from the satellite to the GPS receiver.
$$\begin{align} \mu &= 3.986005\times10^{14} &&\text{WGS84 Value of earth's gravitational constant.} \tag 1\\ \dot{\Omega}_e &= 7.2921151467\times10^{-5} &&\text{WGS84 value of earth's rotation rate.} \tag 2\\ A &= (\sqrt{A})^2 &&\text{Semi-Major Axis.} \tag 3\\ n_o &= \sqrt{\mu \over{A^3}} &&\text{Computed Mean Motion.} \tag 4\\ t_k &= t - t_{oe} &&\text{Time from Ephemeris Reference Epoch.} \tag 5\\ n &= n_o + \Delta n &&\text{Corrected Mean Motion} \tag 6\\ M_k &= M_o + n t_k &&\text{Mean Anomaly} \tag 7\\ M_k &= E_k - e sin E_k &&\text{Kepler's equation for Eccentric Anomaly} \tag 8\\ v_k &= tan^{-1}\left( {\sqrt {1-e^2}sin E_k \over {cos E_k - e}}\right ) &&\text{True Anomaly} \tag 9\\ E_k &= cos^{-1}\left({e + cos v_k \over {1 + e.cos v_k }}\right) &&\text{Eccentric Anomaly} \tag {10}\\ \theta_k &= v_k + \omega &&\text{Argument of Latitude} \tag {11}\\ \delta u_k &= C_{us}.sin 2\theta_k +C_{uc}.cos 2\theta_k &&\text{Correction to Argument of Latitude} \tag {12}\\ \delta r_k &= C_{rs}.cos 2\theta_k +C_{rc}.cos 2\theta_k &&\text{Correction to Radius} \tag {13}\\ \delta i_k &= C_{is}.cos 2\theta_k +C_{ic}.sin 2\theta_k &&\text{Correction to Inclination} \tag {14}\\ u_k &= \theta_k + \delta u_k &&\text{Corrected Argument of Latiude} \tag {15}\\ r_k &= A(1 - e.cos E_k) + \delta r_k &&\text{Corrected Radius} \tag {16}\\ i_k &= i_o + \delta i_k +\dot{i}t_k &&\text{Corrected Inclination} \tag {17}\\ x^\prime_k &= r_k cos u_k &&\text{X Position in Orbital Plane} \tag {18}\\ y^\prime_k &= r_k sin u_k &&\text{Y Position in Orbital Plane} \tag {19}\\ \Omega_k &= \Omega_o + (\dot{\Omega} - \dot{\Omega}_e) t_k - \dot{\Omega}_e t_{oe} &&\text{Corrected Longitude of Ascending Node} \tag {20}\\ x_k &= x^\prime_k cos \Omega_k - y^\prime_k cos i_k sin \Omega_k &&\text{Earth Fixed X coordinate} \tag {21}\\ y_k &= x^\prime_k sin \Omega_k - y^\prime_k cos i_k cos \Omega_k &&\text{Earth Fixed Y coordinate} \tag {22}\\ z_k &= y^\prime_k sin i_k &&\text{Earth Fixed Z coordinate} \tag {23} \end{align}$$
These equations are taken directly from the NAVSTAR GPS User guide and are called Kepler's equations of Orbital Motion, developed by the German Mathematician and Astronomer Johannes Kepler in the 17th Century. Equations 12 - 14 add what are known as Second Harmonic Perturbations to the equations to allow for deviations from a perfectly elliptical orbit, due to the moon and other celestial bodies. Together with several correction factors they are all that is used to calculate the position of a satellite at a given point in time. This may look like a long and laborious process but with one exception these equations can be replicated in the C programming language almost one to one as we shall see later, equation 8 needs to be solved by iteraton, otherwise the values are fixed or supplied by the Ephemeris data transmitted by the satellite whose orbital position is being calculated. C is the language used for programming the none performance critical sections of code in the NavSoft Digital GPS receiver and is the main programming language used for the Open Source code for the Novatel Superstar II receiver.
But first the various corrections made to the satellite clock time and position will be discussed.
Due to aging of the highly accurate clock aboard the GPS Satellite and relativistic effects, a correction needs to be made to the time of the signal transmitted by the GPS satellite to improve it's accuracy and correct for errors. The correction for clock drift or aging effects is in the form of a polynomial where the values used are modelled by a team at the US Navstar Command center. These values are transmitted up to the satellite and then re-transmitted down to your GPS receiver.
$$t = t_{SV} - \Delta t_{SV}$$ $$\Delta t_{SV} = a_{f0} + a_{f1}(t - t_{oc}) + a_{f2}(t - t_{oc})^2 + \Delta t_r$$ $$\Delta t_r = F \times e \times (A)^{\frac12} sin E_K$$ $$F = {-2(\mu)^{\frac12} \over {c^2}} = -4.442807633\times 10^{-10}$$The L1 and L2 correction term \(T_{GD}\), is calculated by the control segment to account for satellite group delay differential based on measurements made by the satellite manufacturer during factory testing. It is needed to compensate for the fact that \(a_{f0}\) clock correction coefficient is based on two frequency ionospheric corrections and not single frequency, as is the usual case for commercial GPS receivers. So for L1 frequency only receivers the clock correction must be corrected by the following
$${(\Delta t_{SV})}_{L1} = \Delta t_{SV} - T_{GD}$$Where \(T_{GD}\) is available from the ephemeris data.
Because the equations for the satellite position are expressed in Earth centered Earth fixed co-ordinates a correction must be made for the earth's rotation during the time the RF signal is travelling from the Satellite to the GPS Receiver. The Earth's rate of rotation is assumed to be a constant value and is supplied in the GPS Navstar User guide as \(\dot\Omega_e = 7.2921151467 \times 10^{-5} \) rads/sec. So the following factors needed to be added to a position to correct for the Earth's rotation during a time interval, \(\Delta t\). \((-\dot\Omega_e y \Delta t, \dot\Omega_e x \Delta t, 0)\) to position \((x,y,z)\). It is assumed that the time interval is small and thus the earth only rotates by a small angle, so the usual sin and cosine functions involved here can use the small angle approximations to give the above result.
Here are the details on actual Satellite Position calculation.