1 Where GPSORB14 begins

GPSORB13 ended with the corrected orbital quantities

(1)

the node-aligned orbital-plane coordinates

(2)

and the Earth-fixed coordinates

(3)

(4)

(5)

IS-GPS-200 states that a user can compute satellite velocity and acceleration, if required, using the equations in parts 3 and 4 of Table 20-IV. Sheet 3 is therefore best understood as the analytical time derivative of the complete sheet-2 position chain [1].

The current GPS public-document index lists IS-GPS-200N as the base L1/L2 interface specification, together with later interface-revision notices. This article derives the Table 20-IV sheet-3 equations as printed in Revision N; operational software should also review any applicable current revision notices [2, 1].

2 A local time derivative and the ephemeris time offset

The sheet-1 time variable is

(6)

with the beginning/end-of-week correction applied so that it represents the actual time difference to the ephemeris reference epoch. Within one continuous evaluation interval,

(7)

This seemingly trivial fact is used repeatedly below. A software implementation should not numerically differentiate across an artificial week-wrap discontinuity; the physical orbit itself is continuous.

The quantities A, e, Δn, ω, the six harmonic coefficients, i₀, IDOT, Ω₀, and the transmitted node-rate parameter are treated as constants over the broadcast fit interval. Their purpose is to parameterize the fitted trajectory over that interval, not to become new dynamical states at every evaluation time.

3 Step 1: eccentric-anomaly rate

GPSORB12 solved Kepler’s equation

(8)

The broadcast mean anomaly is

(9)

where

(10)

Because n is constant over the evaluation interval and Equation (7) gives dt_k∕dt = 1,

(11)

Differentiate Equation (8):

(12)

Factor the eccentric-anomaly rate:

(13)

Substitute Equation (11) and solve for the eccentric-anomaly rate:

(14)

This is the first equation on Table 20-IV sheet 3 [1].

The denominator has a familiar geometric meaning because

(15)

Thus, ignoring the small broadcast radius correction,

(16)

Near periapsis this factor is smaller, so eccentric anomaly advances faster; near apoapsis it is larger, so eccentric anomaly advances more slowly.

4 Step 2: true-anomaly rate

GPSORB12 used the quadrant-safe relation

(17)

To differentiate it cleanly, define

(18)

For an angle

(19)

the derivative is

(20)

Now

(21)

The numerator of Equation (20) becomes

(22)

so

(23)

The denominator is

(24)

which simplifies to

(25)

Therefore

(26)

Applying the chain rule,

(27)

This is the sheet-3 true-anomaly-rate equation [1].

Combining Equations (14) and (27) also gives

(28)

This agrees with Kepler’s second law: angular motion is faster when the satellite is closer to Earth.

5 Why the harmonic phase rate equals the true-anomaly rate

The uncorrected argument of latitude is

(29)

In a complete perturbation theory, the osculating argument of periapsis generally has a nonzero rate. GPSORB10 derived the familiar secular J₂ apsidal precession as one example. The legacy LNAV user model, however, does not transmit a separate argument-of-periapsis-rate term in this calculation. The fitted parameter ω is treated as constant while the other broadcast rates and harmonic coefficients absorb the required short-interval trajectory behavior.

Hence, within this prescribed user model,

(30)

That identity explains why every harmonic derivative on sheet 3 contains twice the true-anomaly rate.

6 The common derivative behind all three harmonic channels

The sheet-2 corrections depend on sin 2Φ_k and cos 2Φ_k. Their time derivatives are

(31)

and

(32)

Using Equation (30),

(33)

(34)

The factor of two is not an arbitrary GPS convention. It is the chain-rule derivative of the twice-per-revolution phase 2Φ_k.

7 Step 3: corrected inclination rate

GPSORB13 used

(35)

where

(36)

Differentiate Equation (35):

(37)

Using Equations (33) and (34),

(38)

Therefore

(39)

This is the corrected inclination-angle rate on sheet 3 [1].

The first term is the fitted linear inclination drift; the second is the instantaneous slope of the twice-per-revolution inclination correction.

8 Step 4: corrected argument-of-latitude rate

The corrected argument of latitude is

(40)

with

(41)

Differentiate:

(42)

Because the harmonic phase rate equals the true-anomaly rate,

(43)

Thus

(44)

The first term is the idealized in-plane angular rate. The second term is the rate contributed by the fitted harmonic correction.

9 Step 5: corrected radius rate

The corrected radius is

(45)

where

(46)

Differentiate the Keplerian part first:

(47)

Next,

(48)

Therefore

(49)

The units provide a useful check: the first term is meters times radians per second, where radians are dimensionless, and the second term is meters times radians per second. Both are meters per second.

10 Step 6: Earth-fixed node-longitude rate

GPSORB09 and GPSORB13 derived

(50)

Here the first rate on the right-hand side is the transmitted node-rate parameter, while the second is the prescribed Earth rotation rate. Differentiate Equation (50). The reference quantities are constant over the fit interval, so

(51)

This is the sheet-3 longitude-of-ascending-node rate [1].

For a GPS-like satellite, the magnitude of the Earth-rate term is much larger than the small inertial node-precession rate. That is expected because Ω_k is an Earth-fixed longitude, not an inertial right ascension.

11 Step 7: orbital-plane velocity

Differentiate the sheet-2 orbital-plane coordinates from Equation (2). For the first component,

(52)

Hence

(53)

For the second component,

(54)

These are the in-plane velocity equations printed on sheet 3 [1].

They are also the Cartesian expansion of the familiar polar-coordinate velocity

(55)

The two terms have clear physical meanings:

• the first term is the radial velocity;
• the second term is the transverse velocity.

Because the local basis vectors are orthonormal,

(56)

This provides a convenient numerical check on Equations (53) and (54).

12 Step 8: derive the final ECEF velocity by the product rule

The cleanest way to understand the final three equations is to rewrite the sheet-2 ECEF position in terms of two time-varying basis vectors. Define

(57)

and

(58)

Then Equations (3)–(5) can be written compactly as

(59)

Differentiate using the product rule:

(60)

This equation separates three sources of ECEF velocity:

1.

motion within the corrected orbital plane through the two in-plane velocity components;

2.

Earth-fixed node sweep through the node-longitude rate;

3.

change in the corrected inclination through the inclination rate.

The derivative of p is

(61)

For q, differentiate both Ω_k and i_k:

(62)

Substituting Equations (61) and (62) into Equation (60) and reading off the components gives the three sheet-3 ECEF velocity equations.

ECEF x velocity

Define the basis-rate combination

(63)

Then

(64)

ECEF y velocity

Define

(65)

Then

(66)

ECEF z velocity

(67)

These expressions are exactly what should emerge from differentiating the sheet-2 position equations. They are not an independent velocity model.

13 Why the ECEF speed is not the inertial orbital speed

The vector

(68)

is a coordinate velocity in the rotating Earth-fixed frame. It should not be confused with the inertial speed from vis-viva. In general,

(69)

so Earth rotation changes the numerical components and the norm of the coordinate velocity. This is the same rotating-frame physics introduced in GPSORB09.

14 Numerical continuation of GPSORB12–GPSORB13

Continue the previous example with

(70)

(71)

(72)

(73)

and the same harmonic coefficients used in GPSORB13:

(74)

(75)

(76)

Also retain

(77)

with

(78)

GPSORB13 already obtained

(79)

(80)

(81)

(82)

(83)

and

(84)

1. Eccentric-anomaly rate

Equation (14) gives

(85)

2. True-anomaly rate

Equation (27) gives

(86)

3. Corrected inclination rate

Using Equation (39),

(87)

The harmonic contribution is therefore small compared with the orbital angular rate, as expected.

4. Corrected argument-of-latitude rate

Equation (44) gives

(88)

5. Corrected radius rate

Equation (49) gives

(89)

The radial speed is much smaller than the transverse speed, which is typical for a low-eccentricity GPS orbit.

6. Earth-fixed node-longitude rate

(90)

so

(91)

7. Orbital-plane velocity

Equations (53) and (54) give

(92)

(93)

The in-plane speed is

(94)

Equation (56) independently gives

(95)

confirming the orbital-plane differentiation.

8. Final ECEF velocity

Substituting the position, angle, and rate quantities into Equations (64)–(67) gives

(96)

(97)

(98)

Therefore

(99)

Its ECEF speed is

(100)

The smaller norm compared with the roughly 3.86 km∕s orbital-plane speed is not a contradiction; ECEF is a rotating coordinate frame.

15 Finite-difference verification

A powerful implementation test is to compare the analytical velocity against a central difference of the full position algorithm. Let

(101)

Using the complete sheet-1 and sheet-2 algorithm with h = 0.1 s for the numerical example gives approximately

(102)

The norm of the difference from the analytical result is below

(103)

for this example. This is an excellent unit test because it validates the chain rule, harmonic signs, node-rate sign, and final ECEF product rule simultaneously.

16 Sanity checks that catch common receiver bugs

Circular-orbit limit

If e = 0, then

(104)

and Equation (27) reduces to

(105)

If the radial harmonic coefficients also vanish,

(106)

Zero-harmonic limit

If

(107)

then

(108)

and

(109)

These are the expected derivatives of the uncorrected broadcast orbit.

In-plane speed identity

After computing the two orbital-plane velocity components, verify

(110)

A discrepancy usually indicates a sign error in one of the product-rule terms.

Differentiate the final position numerically

For a test ephemeris, evaluate the full ECEF position at t−h, t, and t + h and compare the central difference with Equations (64)–(67). This is one of the most effective end-to-end checks of a broadcast-ephemeris implementation.

17 Implementation sequence

A direct implementation of Table 20-IV sheet 3 can follow this order:

1.

Reuse n, E_k, ν_k, Φ_k, u_k, r_k, i_k, x′_k, y′_k, and Ω_k from sheets 1 and 2.

2.

Compute the eccentric-anomaly rate using Equation (14).

3.

Compute the true-anomaly rate using Equation (27).

4.

Compute the corrected inclination rate using Equation (39).

5.

Compute the corrected argument-of-latitude rate using Equation (44).

6.

Compute the corrected radius rate using Equation (49).

7.

Compute the node-longitude rate using Equation (51).

8.

Compute the two orbital-plane velocity components using Equations (53) and (54).

9.

Compute the three ECEF velocity components using Equations (64)–(67).

10.

Verify the in-plane speed identity and, during testing, compare with a central difference of the position algorithm.

All angular quantities used inside trigonometric functions and angular-rate equations must be in radians and radians per second after any required conversion from the broadcast semicircle representation.

18 Connection to the next sheet

Sheet 3 closes the loop from broadcast orbit parameters to ECEF position and velocity:

(111)

Table 20-IV sheet 4 then moves from kinematics back to dynamics. It supplies an ECEF acceleration model containing central gravity, the J₂ oblateness correction, and the Coriolis and centrifugal terms associated with the rotating Earth-fixed frame [1]. That next step reconnects the GPS user algorithm directly to the Newtonian equation that began the GPSORB series.

19 Summary

The sheet-3 velocity equations are not a separate set of empirical formulas. They are the analytical derivative of sheets 1 and 2. The essential chain is

(112)

then

(113)

followed by

(114)

(115)

(116)

(117)

These rates generate the in-plane velocity through

(118)

(119)

and differentiating the time-varying Earth-fixed basis then gives the three ECEF velocity components.

The next article can derive Table 20-IV sheet 4: the J₂ acceleration factor, Earth-fixed gravitational acceleration, Coriolis acceleration, centrifugal acceleration, and the final ECEF acceleration equations.

References