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 [..]):

M0: Mean anomaly at reference time [semicircles]

dn: Mean motion difference from computed value [semicircles/sec]

ec: Eccentricity [unity]

rootA: Square root of semimajor axis [m^0.5]

W0: Longitude of ascending node of orbit plane at weekly epoch [semicircles]

i0: Inclination angle at reference time [semicircles]

w: Argument of perigee [semiocircles]

Wdot: Rate of right ascention [semicircles/sec]

Idot: Rate of inclination angle [semicircles/sec]

Cuc: Amlplitude of the cosine harmonic correction term to the argument of
latitude [rad]

Cus: Amplitude of the sine harmonic correction term to the argument of latitude
[rad]

Crc: Amplitude of the cosine harmonic correction term to the orbit radius
[m]

Crs: Amplitude of the sine harmonic correction term to the orbit radius [m]

Cic: Amplitude of the cosine harmonic correction term to the angle of inclination
[rad]

Cis: Amplitude of the sine harmonic correction term to the angle of inclination
[rad]

Toe: Ephemeris reference time [sec]

IODE: Issueof Data (Ephemeris) [unity]
The constants are:

mu = 3.986005 *10^14  WGS 84 value of the earth's universal gravitation
constant [m^3/sec^2]

Wedot = 7.2921151467 * 10^5  WGS 84 value of the earth's rotation rate
[rad/sec]

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

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.

In the example code this equation is rewritten as E = M + ec * sin(E). This
equation is solved by iteration with E = M as start value.

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
'standalone' program. You can view the program
or download it by clicking here.
Notes:

The code assumes the availability of a math coprocessor. If not available,
declare all 'real' variables as 'double' .

Turbo Pascal requires angles in radians for the goniometric functions.

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.
Back to s/w index