1 From celestial mechanics to the receiver algorithm

The previous articles built the physics in layers. Newton’s inverse-square force produced a conic orbit. Conservation of energy and angular momentum produced the orbital geometry. The classical elements (a,e,i, Ω,ω,M₀) described the orbit, and Kepler’s equation connected elapsed time to position along the ellipse. GPSORB11 then explained why the broadcast model adds fitted parameters such as Δn.

IS-GPS-200 now turns that theory into a receiver algorithm. The first part of Table 20-IV may be organized as

(1)

The elapsed-time variable t_k feeds the mean-anomaly propagation, while e, M₀, Δn, and ω enter at their appropriate points.

The current GPS interface-document page lists IS-GPS-200N as the base revision and later interface-revision notices. This article derives the legacy LNAV user equations appearing in Table 20-IV of Revision N; current IRNs should always be checked when implementing an operational receiver [2, 1].

2 The constants and transmitted quantities used here

Table 20-IV begins with the GPS-user gravitational parameter

(2)

The same sheet also states the Earth rotation rate, but that quantity is not needed until the node-longitude and Earth-fixed-coordinate equations treated later in the sequence.

The subset of broadcast quantities needed in this article is summarized below.

The user equations use ordinary radian trigonometric functions. Therefore a quantity transmitted in semicircles is converted by

(3)

Likewise,

(4)

This unit conversion is essential because n₀ computed from μ and A is naturally in radians per second.

3 Step 1: recover the semimajor axis

The LNAV message transmits , not A itself. The receiver first forms

(5)

Dimensionally,

(6)

so

(7)

The result is the semimajor axis of the Keplerian-looking baseline ellipse used by the broadcast fit.

For a GPS-like orbit,

(8)

corresponds to

(9)

The encoding also has a useful sensitivity relation. Differentiating Equation (5),

(10)

Thus the quantization of the transmitted maps into semimajor-axis quantization through a simple local scale factor.

4 Step 2: derive the computed mean motion

For an ideal two-body ellipse, Kepler’s third law gives

(11)

Mean motion is defined as the angular rate of mean anomaly,

(12)

Substituting Equation (11) into Equation (12),

(13)

so

(14)

The dimensions are

(15)

Because angular measure in the user equations is radians, this is interpreted as radians per second.

For A = 26, 560 km,

(16)

The corresponding period is

(17)

consistent with the familiar approximately twelve-hour GPS orbital period.

5 Step 3: form the time from ephemeris epoch

The next quantity is

(18)

Here t_oe is the ephemeris reference epoch and t is GPS system time associated with the signal transmission epoch used for the satellite-position calculation. IS-GPS-200 notes that the time used in this equation is the transmission time, meaning the reception time has effectively been corrected for signal transit time [1].

This distinction matters. A GNSS receiver wants the satellite position when the electromagnetic signal left the satellite, not where the satellite is when the signal arrives at the antenna.

5.1 Why a week-crossover rule is necessary

GPS seconds-of-week repeat every

(19)

A simple subtraction can therefore return a value almost one week long when the physical time separation is only a few minutes.

IS-GPS-200 requires the actual total time difference and applies a half-week rule:

(20)

The number 302400 s is exactly one half of a GPS week.

An implementation-friendly equivalent is to map the raw difference into the interval

(21)

The conditional form in Equation (20) has the advantage of matching the interface specification directly.

6 Step 4: apply the mean-motion correction

If the orbit were exactly Keplerian, Equation (14) would determine the phase rate completely from A. GPSORB11 explained why the broadcast curve fit provides additional freedom through Δn.

After converting Δn from semicircles per second to radians per second, the receiver forms

(22)

The correction is therefore an additive adjustment to the Keplerian phase rate.

The distinction between n₀ and n is worth preserving in software and derivations:

(23)

If Δn were accidentally added before converting semicircles/s to rad/s, the result would be wrong by a factor of π.

7 Step 5: propagate mean anomaly

At the reference epoch,

(24)

With corrected mean motion n assumed constant over this part of the broadcast model,

(25)

Integrating from t_oe to t,

(26)

which yields

(27)

Using t_k = t − t_oe,

(28)

This equation is the first place where the initial orbital phase M₀, the fitted phase-rate correction Δn, and the receiver evaluation time all meet.

Mean anomaly can be normalized modulo 2π for convenience,

(29)

but the trigonometric equations themselves do not require such normalization. Keeping an unwrapped phase can sometimes be useful for diagnostics, while a wrapped phase can improve readability and numerical bookkeeping.

8 Step 6: solve Kepler’s equation for eccentric anomaly

Mean anomaly is linear in time, but it is not the geometric angle of the satellite on the ellipse. GPSORB06 derived the area relation

(30)

where E is eccentric anomaly.

At time t_k, the receiver therefore must solve

(31)

Because E_k appears both linearly and inside a sine function, there is no elementary algebraic rearrangement that isolates E_k. A numerical root must be found.

8.1 Deriving the Newton-Raphson iteration

Define

(32)

The desired eccentric anomaly satisfies

(33)

Differentiate:

(34)

Newton’s method is

(35)

Substituting Equations (32) and (34),

(36)

Rearranging the numerator gives the form used by IS-GPS-200:

(37)

The specified initial value is

(38)

Revision N states that the refinement uses a minimum of three iterations before the final E_k is taken from the refined value [1].

8.2 Why the denominator is well behaved for an ellipse

For an elliptic orbit,

(39)

Because

(40)

we have

(41)

Thus the denominator in Equation (37) cannot vanish for an ordinary elliptic orbit. This is also the same factor appearing in

(42)

so its positivity has a direct geometric meaning: orbital radius remains positive.

9 Step 7: convert eccentric anomaly to true anomaly

The receiver now knows E_k, but the orbital-plane geometry developed in GPSORB04 and GPSORB08 is naturally expressed with true anomaly ν_k.

For an ellipse,

(43)

and

(44)

Because the common denominator is positive for e < 1, the quadrant of ν_k is determined by the numerator pair. A robust implementation is therefore

(45)

9.1 Deriving the half-angle form

Using

(46)

and substituting Equations (43)–(44),

(47)

The denominator factors as

(48)

Using

(49)

and

(50)

we obtain

(51)

Since

(52)

then

(53)

This is the relation appearing in the Table 20-IV true-anomaly equation.

For software, a particularly robust half-angle implementation is

(54)

Equation (54) and Equation (45) are equivalent. Both preserve the quadrant without having to infer it from a single-argument inverse tangent.

10 Step 8: form the argument of latitude

Table 20-IV sheet 2 begins with the next quantity in the chain,

(55)

This is called the argument of latitude before harmonic correction.

The geometry is simple but important:

• ω is measured in the orbital plane from the ascending-node direction to periapsis;
• ν_k is measured from periapsis to the satellite;
• therefore Φ_k is measured directly from the ascending-node direction to the satellite.

The notation matters because Table 20-IV next applies the twice-per-revolution corrections from GPSORB11. Thus

(56)

is the uncorrected argument of latitude, while the later

(57)

is the corrected argument of latitude. Keeping Φ_k and u_k as distinct variables prevents a common implementation error.

11 A complete line-by-line numerical example

Consider the following GPS-like decoded quantities. Angular quantities are stated here directly in radians or radians per second after any required semicircle conversion:

(58)

(59)

(60)

(61)

(62)

and choose

(63)

1. Semimajor axis

(64)

Therefore

(65)

2. Computed mean motion

Using μ = 3.986005 × 10¹⁴ m³∕s²,

(66)

so

(67)

3. Corrected mean motion

(68)

so

(69)

4. Mean anomaly

(70)

therefore

(71)

which gives

(72)

5. Eccentric anomaly by Newton iteration

Initialize

(73)

Equation (37) gives

(74)

(75)

and

(76)

Thus after three refinements,

(77)

The second and third results already agree to the shown precision because the eccentricity is small.

6. True anomaly

Use the quadrant-safe form

(78)

The result is

(79)

In degrees,

(80)

7. Argument of latitude

Finally,

(81)

so

(82)

This is the phase angle used by the second-harmonic correction equations in the next sheet of Table 20-IV.

12 Sanity checks that catch common receiver bugs

The first-half algorithm is compact enough that several independent checks are inexpensive.

Dimensional check

The units should progress as

(83)

A mismatch usually indicates failure to convert a semicircle quantity to radians.

Circular-orbit limit

If

(84)

then Kepler’s equation becomes

(85)

and Equations (43)–(44) reduce to

(86)

Therefore

(87)

modulo 2π. A circular test case is therefore an excellent unit test.

Kepler residual

After solving for E_k, compute

(88)

The residual should be numerically near zero. This directly tests the iterative solver independently of the later angle conversion.

Two true-anomaly formulas

Compute ν_k using both Equation (45) and Equation (54). Their wrapped angular difference should be near machine precision. This is a strong test for quadrant mistakes.

Week-crossover test

Explicitly test epochs just before and after the GPS week boundary. The example in Figure 2 should return

(89)

not

(90)

13 Compact implementation sequence

Once the derivations are understood, the operational sequence is short:

1.

Decode , e, t_oe, Δn, M₀, and ω.

2.

Convert all semicircle angles and angular rates to radians and radians per second.

3.

Compute A = ()².

4.

Compute n₀ = .

5.

Compute t_k = t − t_oe and apply the half-week crossover rule.

6.

Compute n = n₀ + Δn.

7.

Compute M_k = M₀ + nt_k.

8.

Set E₀ = M_k and perform at least three Newton refinements using Equation (37).

9.

Convert E_k to ν_k with a quadrant-safe formula.

10.

Form Φ_k = ν_k + ω.

The next sheet then applies the second-harmonic corrections to Φ_k, radius, and inclination before constructing the Earth-fixed satellite coordinates.

14 How this maps back to the physics

Every line of this part of Table 20-IV now has a physical origin:

The apparent collection of unrelated receiver formulas is therefore one continuous mechanics chain.

15 Summary

The first half of the legacy GPS broadcast user algorithm can be compressed to

(91)

(92)

(93)

(94)

solved by Newton iteration,

(95)

and finally

(96)

The next article can continue with Table 20-IV sheet 2: δu_k, δr_k, δi_k, corrected u_k, r_k, i_k, orbital-plane coordinates, corrected node longitude, and the final ECEF position.

References