### Pseudorange Corrections

Pseudoranges contain a number of errors (see e.g. the theory page). Some of these errors can be estimated and used as corrections to the pseudoranges. In this page the following corrections are described:

1. SV clock correction,
2. ionospheric delay correction, and
3. tropospheric delay correction.

Again, the GPS SPS Signal Specification is followed closely.
A fourth correction, the receiver clock correction, results from the final stand alone position calculation and is emphasized in that page.

1. SV Clock Correction

The parameters, required to calculate the SV clock correction are present in subframe 1 of the navigation message (see e.g. the conversion page). They are:

• Tgd: Group Delay [sec]
• Toc: SV Clock reference time [sec]
• af2: second order polynomial coefficient [sec/sec^2]
• af1: first order polynomial coefficient [sec/sec]
• af0: zero order polynomial coefficient [sec]

The SV clock correction for the C/A code pseudorange observation is given by the expression:

dTclck = af0 + af1 * (Ttr - Toc) + af2 * (Ttr - Toc)^2 + F * ec * rootA * sin(E) - Tgd

with:

• dTclck: SV clock correction [sec]
• Ttr: Time of transmission [sec]
• F: a constant with the value -4.442807633E-10 [sec/m^0.5]
• ec: SV orbit eccentricity [unity]
• rootA: Square root of SV orbit semi major axis [m^0.5]
• E: SV orbit eccentric anomaly [rad] or [deg]

The implementation of this expression is straightforward.
The term: F * ec * rootA * sin(E) is called the relativistic correction term Trel [sec]. It may be a good idea to calculate this term within the 'satpos' procedure since the required parameters are available there. The term: T = Ttr - Toc should be corrected for end of week crossovers: if T > 302400 then subtract 604800 from T, if T < -302400 then add 604800 to T.

Below follow fragments of the satpos procedure with the additions for implementation of Trel in italics.

procedure    satpos (eph : vec16;        {ephemeris}
Ttr  : real;           {satellite GPS time}
var Trel: real;          {satellite relativistic correction term}
var X   : vec3);         {satellite position}

const
pi = 3.1415926535898;        {WGS 84 value of pi}
mu = 3.986005E+14;            {WGS 84 value of earth's univ. grav. par.}
Wedot = 7.2921151467E-5; {WGS 84 value of earth's rotation rate}
F = -4.442807633E-10;        {rel cor term constant}

var
M0, dn, ec, A, W0, i0, w Wdot, Cuc, Cus, Crc, Crs, Cic, Cis, Toe, Idot: real;
n0, M, snu, cnu, nu, phi, du, dr, di, u, r, i, Xdash, Ydash, Wc: real;
i: integer;

begin

code .....
..........
code ....

X := Xdash*cos(Wc) - Ydash*cos(i)*sin(Wc);
X := Xdash*sin(Wc) + Ydash*cos(i)*cos(Wc);
X := Ydash*sin(i)

Trel := F * ec * rootA * sin(e); {relativistic correction term}

end; {procedure satpos}

An example. Take the clock correction parameters from the example subframe 1 in the conversion page, and take an SV transmission time of 588480. When testing the satpos procedure with the same transmit time, I printed the parameter E = -1.5747877 rad. Entering all parameters in the expression above gives for dTcl: -3.649408E-6 sec, or, multiplied by the speed of light: -1094.065 m.

2. Ionospheric Delay.

Subframe 4 page 18 contains the ionospheric delay coefficients. Furthermore, the approximate user position, and the SV's azimuth and elevation are required. And here we run into a problem: the user position is the final result of the calculation, and is not yet known at this stage. Some solutions are:

1. estimate the current position by extrapolation from a previously calculated position, if available of course, or
2. carry out a preliminary position calculation forgetting about the ionospheric delay correction and use the preliminary position in a second run to determine the iono delay correction.

In the following I assume, that the approximate user position is known. The following calculation steps are extracted from 'Ionospheric Effects on GPS' by J.A. Klobuchar. All angular units are in semicircles (!), El is the SV's elevation [sc], Az [sc] the SV's azimuth, Latu [sc] the appr. user latitude, and Lonu [sc] the appr. user longitude(the calculation of El and Az is given in the Calculation of SV's Elevation and Azimuth page).

• phi = 0.0137 / (El + 0.11) - 0.022    Earth-centered angle [sc]
• Lati = Latu + phi * cos(Az)     Subionospheric lattitude [sc] (1)
if Lati > +0.416 then Lati = +0.416
if Lati < -0.416 then Lati = -0.416
• Loni = Lonu + phi * sin(Az) / cos(Lati)     Subionospheric longitude [sc] (1)
• Latm = Lati + 0.064 * cos(Loni - 1.617)     Geom. lat. of the iono inters. point [sc] (1)
• T = 4.32E+4 * Loni + Ttr     Local time at sunionospheric point [sec] (2)
if T > 86400 then subtract 86400 from T
if T < 0 then add 86400 to T
• F = 1.0 + 16.0 * (0.53 - El)^3     Slant factor [unity]
• x = 2 * pi * (T - 50400) / (bo + b1 * Latm + b2 * Latm^2 + b3 * Latm^3)     'x' [rad]
• Ionospheric delay correction [sec]:
if | x | > 1.57 then dTiono = F * 5.0E-9
otherwise dTiono = F * {5.0E-9 + (a0 + a1 * Latm + a2 * Latm^2 + a3 * Latm^3) * (1.0 - x^2 / 2 + x^4 / 24) }

Remarks:

1. Take care: arguments of goniometric functions should be converted to the unit required by you programming language.
2. Ttr is the GPS time of transmission.

An example of a Turbo Pascal procedure has been incorporated in a Turbo Pascal program to calculate receiver position from raw data.This program can be viewed or downloaded by clicking here.

3. Tropospheric Delay.

Two strategies can be followed to determine the tropo delay:

1. Use 'standard' atmosphere parameters:
• ambient air temperature Tamb = 15 Cel
• ambient air pressure Pamb = 101.325 kPa
• ambient vapour pressure Pvap = 0.85 kPa (corresponding with a relative humidity of 50 % at the ambient temperature)
2. Use actual measurements at the receiver antenna location of Tamb, Pamb and Pvap.

There are numerous tropo models availble, I picked Hopfield's model which splits the tropo delay into two parts: a contribution of the 'dry' atmosphere and a contribution of the 'wet' atmosphere. For a detailed explanation, see any textbook on GPS.

The zenith delay of the dry component is given by:
Kd = 1.55208E-4 * Pamb * (40136.0 + 148.72 * Tamb) / (Tamb + 273.16) [m]

The zenith delay of the wet component is given by:
Kw = -0.282 * Pvap / (Tamb + 273.16) + 8307.2 * Pvap/ (Tamb + 273.16)^2  [m]

Multiply the zenith delay's with their 'mapping functions' to correct for elevations lower than 90 deg and add the components to obtain the SV's tropospheric delay correction:

dRtrop = Kd / sin(sqrt(El * El + 1.904E-3))  +  Kw / sin(sqrt(El * El + 0.6854E-3)) [m]

with El the SV's elevation in [rad].

When we use the standard atmosphere values, the expression reduces to

dRtrop = 2.312 / sin(sqrt(El * El + 1.904E-3)) + 0.084 / sin(sqrt(El * El + 0.6854E-3)) [m]

Example. Using standard atmosphere values for an SV with an elevation of 12.86 deg the Tropo delay is 10.575 m.

4. The corrected pseudorange.

The last step in this page is the determination of the corrected pseudorange:

PRcorr = PR + dTclck * c - dTiono * c - dRtrop [m], with c the speed of light.

Suppose the SV 30's raw pseudorange is 24401509.104 m, then the corrected range is:
24401509.104 - 1094.065 - 3.811 - 10.575 = 24400400.653 m