1 Where GPSORB13 begins

GPSORB12 derived the first half of the legacy broadcast-position algorithm:

(1)

with

(2)

The quantity Φ_k is the uncorrected argument of latitude: the angle in the orbital plane from the ascending-node direction to the satellite before the broadcast second-harmonic corrections are applied.

The second sheet of IS-GPS-200 Table 20-IV continues with the sequence

(3)

The current GPS public-document index identifies IS-GPS-200N as the current base revision of the L1/L2 space-to-user interface specification. This article derives the legacy LNAV Table 20-IV equations in Revision N; an operational implementation should also check any applicable interface-revision notices [2, 1].

2 Broadcast quantities used in this half of the calculation

In addition to A, e, E_k, Φ_k, t_k, and t_oe already available from GPSORB12, the second-half position calculation uses

IS-GPS-200 denotes the node-rate parameter by . Table 20-III gives the harmonic latitude and inclination coefficients directly in radians, while the radial pair is in meters. The node, inclination, and rate parameters transmitted in semicircles or semicircles per second must be converted to radians or radians per second before being mixed with radian-valued quantities [1].

The GPS user value of the Earth rotation rate is

(4)

3 Step 1: why the correction functions use twice the argument of latitude

GPSORB10 showed that the leading nonspherical Earth-gravity perturbation contains degree-two angular structure. In orbital geometry one encounters terms such as

(5)

Using

(6)

reveals a natural twice-per-revolution harmonic. GPSORB11 then emphasized an important modeling distinction: the broadcast coefficients are not obtained by simply substituting a closed-form J₂ solution. They are fitted coefficients in the prescribed broadcast model, but a sine/cosine basis at 2Φ_k is physically well matched to the dominant degree-two periodic structure.

Over one full revolution of Φ, both sin 2Φ and cos 2Φ complete two cycles.

A general second harmonic with arbitrary phase can be written

(7)

Therefore a pair of coefficients multiplying sine and cosine can represent both amplitude and phase. If

(8)

then

(9)

This is why each corrected orbital quantity receives a pair of harmonic coefficients rather than a single amplitude.

4 Step 2: argument-of-latitude correction

The first Table 20-IV correction is

(10)

Both coefficients have angular units, so

(11)

The receiver then applies the correction directly to the uncorrected argument of latitude:

(12)

The notation is worth preserving carefully:

(13)

while

(14)

A software implementation that overwrites Φ_k too early can obscure which angle belongs in the harmonic terms.

5 Step 3: radial correction

The corresponding radial harmonic is

(15)

Because C_rs and C_rc are measured in meters,

(16)

Before this correction is added, the Keplerian ellipse geometry from GPSORB06 gives

(17)

That result follows directly from the focus-centered eccentric-anomaly parameterization of the ellipse.

The GPS broadcast model modifies the radius by

(18)

It is useful to read Equation (18) as

(19)

6 Step 4: inclination correction and IDOT

The harmonic inclination correction is

(20)

Again,

(21)

The complete corrected inclination is

(22)

This equation contains three distinct pieces:

(23)

(24)

and

(25)

As emphasized in GPSORB11, IDOT should not be interpreted as a claim that first-order secular J₂ theory predicts a nonzero mean inclination rate. The broadcast value is a fitted coefficient for the complete propagated reference trajectory.

7 Step 5: form the corrected orbital-plane coordinates

After u_k and r_k are known, the satellite’s coordinates in the corrected orbital plane are ordinary polar-to-Cartesian relations:

(26)

(27)

The x′_k direction lies along the ascending-node direction in the orbital plane, and the y′_k axis is 90^∘ ahead within that plane.

The equations immediately provide a useful invariant:

(28)

so

(29)

This is an inexpensive software sanity check.

8 Step 6: corrected longitude of the ascending node

GPSORB09 derived the Earth-fixed node equation from the difference between inertial node motion and Earth rotation. Table 20-IV uses

(30)

Because

(31)

Equation (30) may be rearranged as

(32)

This alternate form exposes the physics:

• (t − t_oe) propagates the fitted orbital-node orientation;
• −_et accounts for the rotation of the Earth-fixed axes beneath the orbit.

A common implementation concern is the size of Ω_k. Trigonometric functions are periodic, so it is mathematically acceptable to use the unwrapped value. In numerical software it is often convenient to reduce it modulo 2π before evaluating sines and cosines.

9 Step 7: derive the final ECEF rotation

The orbital-plane vector is

(33)

The geometry is the same rotation structure derived in GPSORB08, but the quantities are now the corrected broadcast quantities and Ω_k is already an Earth-fixed node longitude.

First tilt the orbital plane by the corrected inclination i_k, then rotate the node direction by Ω_k. With the rotation convention used throughout this series,

(34)

The combined matrix is

(35)

Because the third component of r′_k is zero, Equation (35) reduces immediately to the three Table 20-IV equations:

(36)

(37)

(38)

IS-GPS-200 identifies these as the satellite antenna phase-center coordinates in the WGS 84 Earth-centered Earth-fixed coordinate system [1].

10 Why there is no extra ECI-to-ECEF rotation afterward

This point deserves emphasis because it is an easy place to double-count Earth rotation.

In a conventional inertial orbital-elements calculation one might first form an inertial position vector

(39)

after which a separate Earth-rotation matrix converts ECI to ECEF. The legacy GPS broadcast algorithm instead uses the already Earth-fixed node longitude Ω_k from Equation (30). Therefore Equations (36), (37), and (38) already produce ECEF coordinates.

Thus the broadcast sequence is

(40)

not

(41)

Applying an additional Earth-rotation transformation would rotate the answer twice.

11 A complete numerical continuation of GPSORB12

Continue the GPSORB12 example, which ended with

(42)

(43)

and

(44)

For an illustrative GPS-like continuation, choose the decoded correction coefficients

(45)

(46)

(47)

Also take

(48)

(49)

(50)

and

(51)

All angular parameters in this worked example are already expressed in radians or radians per second after any required broadcast semicircle conversion.

1. Twice-per-revolution angle

The harmonic phase is

(52)

2. Argument-of-latitude correction

From Equation (10),

(53)

so

(54)

The corrected argument of latitude is therefore

(55)

3. Radius correction

Equation (15) gives

(56)

The uncorrected Keplerian radius is

(57)

so

(58)

Equivalently,

(59)

4. Inclination correction

Equation (20) gives

(60)

The linear term is

(61)

Hence

(62)

5. Orbital-plane coordinates

Using Equations (26) and (27),

(63)

(64)

The norm check gives

(65)

to numerical precision.

6. Earth-fixed node longitude

Equation (30) gives

(66)

so the unwrapped result is

(67)

Reducing modulo 2π gives the equivalent angle

(68)

Either form gives the same sine and cosine.

7. Final ECEF coordinates

Substitution into Equations (36), (37), and (38) gives

(69)

(70)

(71)

In kilometers,

(72)

A final norm check gives

(73)

which verifies that the final transformation changed only coordinate orientation, not physical distance from Earth.

12 Sanity checks and implementation traps

1. Keep Φ_k and u_k distinct

The harmonic corrections are evaluated using

(74)

not 2u_k. First compute δu_k, then form u_k = Φ_k + δu_k.

2. Keep the three correction units distinct

The correction channels have different dimensions:

(75)

Adding a radial coefficient to an angle or vice versa is a detectable unit error.

3. Convert semicircle quantities before mixing them

C_uc,C_us,C_ic,C_is are already specified in radians, but i₀, Ω₀, , and IDOT are broadcast in semicircle-based units. Convert them before using Equations (22) and (30).

4. Check both coordinate norms

The two identities

(76)

and

(77)

should hold to numerical precision. The first tests the polar decomposition; the second tests the rotation into ECEF.

5. Do not apply Earth rotation twice

The node equation already includes the Earth-rate term. The final Table 20-IV coordinates are ECEF, not ECI.

6. Zero-correction limit

If

(78)

then

(79)

This recovers the uncorrected Keplerian geometry except for the separately modeled node evolution.

7. Equatorial limit

For

(80)

Equation (38) gives

(81)

and the x-y equations reduce to an ordinary planar rotation by Ω_k.

13 Compact implementation sequence

Given the outputs of GPSORB12 and the remaining decoded ephemeris parameters:

1.

Evaluate sin 2Φ_k and cos 2Φ_k once.

2.

Compute δu_k, δr_k, and δi_k.

3.

Form u_k = Φ_k + δu_k.

4.

Form r_k = A(1 − e cos E_k) + δr_k.

5.

Form i_k = i₀ + IDOT t_k + δi_k.

6.

Compute x′_k = r_k cos u_k and y′_k = r_k sin u_k.

7.

Compute Ω_k = Ω₀ + ( −_e)t_k −_et_oe.

8.

Compute x_k,y_k,z_k from Equations (36), (37), and (38).

9.

Verify ≈ r_k.

A practical implementation can reduce repeated trigonometric work by caching

(82)

14 How the equations map back to the mechanics

The second half of Table 20-IV is not an arbitrary list of coordinate recipes. Each group has a clear origin:

The full legacy position algorithm can now be read as one continuous chain from Newtonian orbit parameters to an Earth-fixed satellite position:

(83)

(84)

15 Summary

Starting with the uncorrected argument of latitude from GPSORB12,

(85)

Table 20-IV forms

(86)

(87)

(88)

The corrected geometry is

(89)

(90)

(91)

The node-aligned orbital-plane coordinates are

(92)

and the Earth-fixed node longitude is

(93)

Finally,

(94)

(95)

(96)

which completes the legacy broadcast calculation of the satellite’s WGS 84 ECEF position.

The next article can continue to Table 20-IV sheet 3 and derive the satellite velocity equations analytically by differentiating this complete position chain.

References