Calculation of Satellite Position

In this page the algorithms are given to calculate a satellite's position (in ECEF XYZ) from the ephemeris parameters in the navigation message. No theoretical background is given, see any textbook on GPS. The algorithms follow closely  the GPS Signal Specification (which I strongly recommend to download, about 50 pages, forget about the annexes). The algorithms are slightly adapted to be readible with simple browsers.

The required ephemeris parameters are (units are given between [..]):

1. M0: Mean anomaly at reference time [semi-circles]
2. dn: Mean motion difference from computed value [semi-circles/sec]
3. ec: Eccentricity [unity]
4. rootA: Square root of semi-major axis [m^0.5]
5. W0: Longitude of ascending node of orbit plane at weekly epoch [semi-circles]
6. i0: Inclination angle at reference time [semi-circles]
7. w: Argument of perigee [semio-circles]
8. Wdot: Rate of right ascention [semi-circles/sec]
9. Idot: Rate of inclination angle [semi-circles/sec]
10. Cuc: Amlplitude of the cosine harmonic correction term to the argument of latitude [rad]
11. Cus: Amplitude of the sine harmonic correction term to the argument of latitude [rad]
12. Crc: Amplitude of the cosine harmonic correction term to the orbit radius [m]
13. Crs: Amplitude of the sine harmonic correction term to the orbit radius [m]
14. Cic: Amplitude of the cosine harmonic correction term to the angle of inclination [rad]
15. Cis: Amplitude of the sine harmonic correction term to the angle of inclination [rad]
16. Toe: Ephemeris reference time [sec]
17. IODE: Issueof Data (Ephemeris) [unity]

The constants are:

1. mu = 3.986005 *10^14 - WGS 84 value of the earth's universal gravitation constant [m^3/sec^2]
2. Wedot = 7.2921151467 * 10^-5 - WGS 84 value of the earth's rotation rate [rad/sec]
3. pi = 3.1415926535898 - WGS 84 value for pi.

The algorithm below calculates the satellite position at the time of transmission Ttr.

The calculation steps:

• A =  (rootA)^2     Semi major axis
• n0 = sqrt(mu / A^3)     Computed mean motion
• T= Ttr- Toe    Time from ephemeris reference epoch (1)
• n = n0 + dn     Corrected mean motion
• M = M0 + n * T    Mean anomaly
• M = E - ec * sin(E)     Kepler's equation for Eccentric Anomaly (2)
• snu = sqrt(1 - ec^2) * sin(E) / {1 - ec * cos(E)}    sine of true anomaly
• cnu = {cos(E) - ec} / {1 - ec * cos(E)}    cosine of true anomaly 
• nu = invtan(snu / cnu)   True anomaly (3)
• phi = nu + w   Argument of latitude
• du = Cus * sin(2*phi) + Cuc * cos(2 * phi)    Argument of latitude correction
• dr = Crs * sin(2*phi) + Crc * cos(2 * phi)    Radius correction
• di = Cis * sin(2*phi) + Cic * cos(2 * phi)    Correction to inclination
• u = phi + du     Corrected argument of latitude
• r = A * {1 - ec * cos(E)} + dr   Corrected radius
• i = i0 + di + idot * T     Corrected inclination
• Xdash = r * cos(u)     Position x in orbital plane
• Ydash = r * sin(u)    Position y in orbital plane
• Wc = W0 + (Wdot - Wedot) * T - Wedot * T_oe     Corrected longitude of ascending node
• X = Xdash * cos(Wc) - Ydash * cos(i) * sin(Wc)    ECEF x
• Y = Xdash * sin(Wc) + Ydash * cos(i) * cos(Wc)    ECEF y
• Z = Ydash * sin(i)                                                ECEF z

Remarks

1. Ttr is the GPS system time at time of transmission, i.e., GPS time corrected for transit time (range/speed of light). Furthermore, T shall be the actual total time difference between the time Ttr and the epoch time T_oe, and must account for beginning or end of week crossovers. That is, if t is greater than 302400 sec, subtract 604800 sec from T. If Tis less than -302400 sec, add 604800 sec to T.
2. In the example code this equation is re-written as E = M + ec * sin(E). This equation is solved by iteration with E = M as start value.
3. invtan is the inverse tangent of the expression.

A Turbo Pascal procedure has been produced according to the above algorithm. This procedure has been incorporated in my 'stand-alone' program. You can view the program or download it by clicking here.

Notes:

1. The code assumes the availability of a math co-processor. If not available, declare all 'real' variables as 'double' .
2. Turbo Pascal requires angles in radians for the goniometric functions.
3. Nobody is perfect. Even I am not ! Please check the info in this page and mail your comment/ corrections to me.

This page was revised on 22 April 1998.