The recent availability of relatively inexpensive dual-frequency receivers for signals from Global Navigation Satellite Systems (GNSS) provides access to ultra-precise, real-time data such as positions and velocities of dozens of satellites orbiting the Earth. We discuss how these data can be obtained, processed, and analyzed either with or without the actual purchase of a GNSS receiver. The positional information can be used to verify Kepler's three laws at lowest order as well as to reveal the presence of higher-order perturbations such as the oblateness of the Earth and the gravitational influences of the Sun and Moon on these satellites. The supplementary material includes both introductory laboratory exercises and Python scripts used to gather and process data suitable for intermediate courses.

While students of introductory physics may feel apprehensive about many parts of their textbook, most do not give its appendices a second thought. It is interesting to ask them “Where did the authors get the value for the Earth's mass?” The easy answer would be Google, or NASA, but where did NASA get it? We like to ask this question when teaching Newton's law of gravitation and Kepler's three laws of planetary motion. To add to the mystery, we point out that if we had access to a physics or astronomy text from the early part of the 20th century,1 we would find in its appendix that the masses for the outer planets agree with modern values to better than 1%, but the values for Venus and especially Mercury do not.

We ask the students to consider what might have happened in the 1960s to change this situation. We discuss that the gold standard for mass determination involves observations of the orbits of a planet's moons. It is usually necessary to mention that the only moonless planets are Mercury and Venus, and then students see that accurate values for their masses had to wait until artificial satellites were sent to orbit them.

A common laboratory exercise2–4 in both physics and astronomy courses involves using successive photos of Jupiter's four largest moons to establish its mass. These laboratory exercises are worthwhile, but the collection of data can be a significant hurdle. The various satellite constellations devoted to navigation allow real-time data to be collected at much greater rates and higher precision than would be possible with small telescopes. In this way, students can both verify the basic principles of Kepler's laws and see that higher order effects are also present. In Sec. II, we review Kepler's laws, and in Sec. III, we describe the process of choosing which satellites to use for the project. Section IV explains the collection and format of the data, while Secs. V–VII discuss introductory, intermediate, and more advanced projects.

Kepler's laws, originally developed to match observations of planetary orbits, were shown by Newton to apply to any two bodies orbiting under the influence of a 1/r potential. Derivations of orbital mechanics can be found in textbooks5–7 as well as this journal.8,9 We state the generalized versions of Kepler's laws as follows:

1. The orbit of a satellite is an ellipse with the primary body at one focus.

2. A line connecting the satellite to the primary sweeps out equal areas in equal times.

3. The ratio of the cube of the semi-major axis of the satellite's orbit (a) to the square of the period of its orbit (P) is proportional to the sum of the masses of satellite and primary.

While there is no shortage of satellites in orbit around the Earth, they are not all equally useful for our purposes. Websites, such as celestrak.org, give the basic information needed to calculate a satellite's position in orbit in the form of TLEs, or two-line elements. This specification dates to the 1960s and includes the six Keplerian elements. The semi-major axis, a, and eccentricity, e, describe the size and shape of the orbit, and Fig. 1 shows the origin of the remaining four elements. While TLEs allow compact encoding of satellite data, errors in position can be several kilometers.10–12

Fig. 1.

Keplerian orbital elements. The inclination angle i is measured relative to the plane of the Earth's equator. The line of intersection of the orbital plane and reference plane is the line of nodes, with the ascending node being the point where the satellite crosses the reference plane moving south to north. The longitude of the ascending node Ω is the angle between the ascending node and the reference direction (the apparent position of the Sun during the vernal equinox) measured in the reference plane. The argument of periapsis (here perigee) ω is the angle between the ascending node and the point of the satellite's closest approach to the Earth, measured in the orbital plane. Finally, the true anomaly T represents the satellite's progress along its orbit (Ref. 13).

Fig. 1.

Keplerian orbital elements. The inclination angle i is measured relative to the plane of the Earth's equator. The line of intersection of the orbital plane and reference plane is the line of nodes, with the ascending node being the point where the satellite crosses the reference plane moving south to north. The longitude of the ascending node Ω is the angle between the ascending node and the reference direction (the apparent position of the Sun during the vernal equinox) measured in the reference plane. The argument of periapsis (here perigee) ω is the angle between the ascending node and the point of the satellite's closest approach to the Earth, measured in the orbital plane. Finally, the true anomaly T represents the satellite's progress along its orbit (Ref. 13).

Close modal

The methods in this paper are applicable to any satellites, but different orbital altitudes can be expected to provide different modeling challenges. Satellites in low-Earth orbit (generally defined as orbiting at altitudes of 1000–2000 km) are strongly affected by atmospheric drag, which depends on their shape, size, mass, and orientation as well as the density of the upper atmosphere they traverse (which is itself subject to fluctuations due to solar activity, among other things).

At the other end of the distance spectrum, satellites far from the Earth could be expected to see a gravitational potential closer to the (perhaps less interesting) Keplerian point-mass model. The various GNSS orbit at altitudes around 20 000 km, where air drag is nonexistent, and the details of Earth's mass distribution are still important. The most significant factor in choosing navigational satellites for analysis is that their positions and velocities must be known to very high accuracy. We cannot use the satellites to locate our position to within centimeters without knowing the positions of the satellites to similar accuracy.

Satellites in the United States' Global Positioning System (GPS) broadcast ephemerides that are used to calculate the satellite's position at a known time (epoch) and for some period after that. The ephemerides are updated every two hours, providing real-time positional accuracy to within a meter or two. The results can also be post-processed using precise-orbit files (SP3 files) typically accurate to about a centimeter. Access to data that are correct at (or below) the parts-per-billion level allows us to measure extremely small perturbations in the orbits of these satellites. While we chose to use the GPS due to its long and well-documented history of operation, our methods are adaptable to any of the other GNSS.

There are at least a handful of other satellites that could serve as interesting comparisons to the GPS. SP3 files containing both position and velocity information at intervals as short as 30 s are available for the LAGEOS satellites (6000 km altitude), LARES (1400 km), and Starlette (800 km), among others. We would expect orbits of satellites at these lower altitudes to exhibit stronger perturbations due to Earth's non-sphericity and weaker perturbations due to the Sun and Moon. The lowest of these orbits should reveal information about the drag force acting on the satellites.

The operating principle behind the various GNSS constellations14 is the simultaneous reception of time-coded radio broadcasts from multiple satellites. By knowing the position of each satellite and the travel time of the radio signal, imaginary spheres can be constructed around each satellite and their region of intersection (ideally a point) provides the receiver's location. For the GPS, each satellite broadcasts a unique binary pseudorandom noise (PRN) sequence. The knowledge of that sequence allows receivers to identify the sending satellite, and time-shifting a locally generated copy of the signal until maximum correlation is achieved provides the time difference between the sending and receiving clocks. Multiplying the speed of light by this (apparent) time difference yields a value known as the pseudo-range. After correction for clock errors on the satellite and in the receiver, which include relativistic effects on the satellite clocks,15 the actual travel time of the signal, and therefore the precise distance between the sender and receiver, can be found.

Data are encoded in the signal at a rate of 50 bits per second, and the information sent in a full navigation message takes 12.5 min to be transmitted. The data sent in each 300-bit (six second) subframe page are detailed in the GPS Interface Specification (GPS-IS)16 and include information about the orbital parameters, clock errors, and correction terms for the sending satellite, as well as almanac data for the other 31 GPS satellites.

We utilized the U-blox ZED-F9T on a sparkfun.com board17 and a multi-band antenna capable of receiving both L1 and L2 signals, as well as similar signals from the other GNSS. The total cost for the pair was under \$400. We typically record all of the so-called raw data as well as subframe data over several days and use Python scripts18 to process these recordings and produce files of comma-separated values (CSV files), which can then be analyzed with Excel or other software.

The raw data include Doppler shift and pseudo-range measurements, as well information about the satellite and receiver clocks that allows the pseudo-range measurements to be converted into an estimate of the distance between satellite and receiver. The accuracy of this estimate degrades when the satellite is near the horizon, but when it is nearest to the receiver, it can be correct to within about a meter. As the satellites are approximately 20 000 km above the surface of the Earth, this represents an accuracy significantly better than one part per million.

The GPS-IS describes the model used to calculate positions and velocities at each instant of time from the Keplerian orbital elements. These elements are the eccentricity e of the orbit (less than 0.03 for GPS satellites), the square root of the semi-major axis, the inclination, the longitude of the ascending node, the argument of perigee, and the true anomaly. Rates of change of these quantities and several harmonic coefficients are also broadcast.

The orbital parameters needed to calculate satellite positions are uploaded to the GPS satellites by the GPS Control Segment (or GPS Ground Segment). The high precision required means that the six Keplerian elements above are not sufficient to accurately locate the satellites, as the satellites are subject to the non-spherical nature of the Earth's mass distribution as well as the influences of the Sun and Moon. Just as we could replace an intricate curve by a series of (many) straight lines, the Control Segment finds the best-fit Keplerian ellipse (known as the osculating orbit) and broadcasts those parameters to the satellites, updating them every two hours.

If the user is willing to sacrifice measured pseudo-range and Doppler shift information, the financial expense can be eliminated, as the data broadcast by the GPS satellites are also available from NASA19 as NAV files in RINEX20 format. In addition to cost and ease of acquisition, a further advantage of downloading these files instead of collecting them from the broadcast is the fact that a GPS satellite is only in view of a particular place on the Earth for several hours each day. The BRDC NAV files produced by the CDDIS19 represent a compilation of broadcast data collected by sites all around the globe, providing information for the entire day. We have Python scripts18 to collect these BRDC files and calculate satellite kinematic data just as we would from data received locally.

Yet another alternative is downloading SP3 precise orbit files (again in ECEF format), which contain both position and velocity information. While the CDDIS does have this enhanced type of SP3 file for some satellites, we have not found it for GPS satellites. The National Geospatial-Intelligence Agency21 provides these files.

One of the quirks of the GPS system is the use of an Earth-centered, Earth-fixed (ECEF) Cartesian coordinate frame. In this frame, the accelerations of the satellites include both centrifugal and Coriolis terms. To avoid these complications, our Python scripts transform the coordinates to an Earth-centered inertial (ECI) Cartesian frame.

Our ECEF-ECI transformation sacrifices absolute accuracy for speed and simplicity. We assume the z coordinates in the two systems are equal, when, in fact, the Earth's axis of rotation wanders about slightly due to precession, nutation, and other effects. We arbitrary set the time of coincidence of these two frames to be the fall equinox in 2020 at the UTC time when the 0h line of right ascension is in line with the prime meridian. All other observations are rotated about the z axis based on the relatively constant rate of Earth's rotation multiplied by the time in seconds since the time of coincidence. The direction of the vernal equinox is our x axis, and the y axis is chosen to form a right-handed coordinate system.

The verification of Kepler's three laws of planetary motion is the most straightforward laboratory exercise. We process data files, either received locally and stored in UBX format or downloaded from one of the online sources, using our Python scripts that generate 32 CSV files. These files contain data recorded and calculations updated as often as every second for ephemeris-derived files. For locally received files, we return Doppler shift, pseudo-range, and carrier phase values as well as the data calculated as specified in the GPS-IS. Some of the data are from the ephemeris broadcast, and as such are updated no more than once every two hours.

Opening one of the kinematic files in a spreadsheet program and plotting the radial distance of the satellite from the Earth's center as a function of time will reveal that it is in fact not constant. For a slightly more quantitative analysis, finding the perigee and apogee for a particular satellite will allow the calculation of the semi-major axis (half of the sum of the perigee and apogee distances, written as a) and the distance between the foci of the ellipse (the difference between the apogee and perigee distances, denoted 2 c). The ratio of c to a is known as the eccentricity, e.

This is a bit easier said than done if using local data, since a particular GPS satellite will only be above the horizon for about 7 h out of its period of one half of a sidereal day (11 h 58 min). The chance of catching both perigee and apogee of a particular satellite from one spot on the globe is relatively slim, but searching through all 32 satellites will generally yield at least one candidate. A reviewer has pointed out that ellipse parameters could also be found by using Excel's built in Solver function to fit the data to the polar equation for an ellipse.

An alternative is to use a different Python script that downloads SP3 files and splits them by satellite into 32 files. The SP3 files provide position information accurate to better than 1 cm at five-minute increments throughout the day. With one of these files, the semi-major axis and eccentricity can be found quickly. These data can be compared to (or, alternately, extracted directly from) the ephemeris data provided in the kinematic CSV files. Both eccentricity e and the square root of the semi-major axis a are part of the regular ephemeris broadcast.

Skipping over Kepler's second law temporarily, the third law is quickly verified. There are a variety of ways to discover the satellite's period, if data are recorded for more than 12 h. A spreadsheet can be used to find the distance between one reference position and every other position captured throughout the log. If local data are used, the minimum will occur just under 24 h later. In that case, the time difference must be divided by two, since the satellite will return to a given point in space about 12 h after its first appearance, but the recording location will then be on the opposite side of the Earth.

Earth's mass ME can then be found from
$M E= 4 π 2 a 3 P 2 G,$
(1)
where G is Newton's gravitational constant.

The calculations performed by our Python scripts are all without reference to the satellite's mass. There are multiple designs of GPS satellites (known as Blocks) currently in orbit, so there is no single mass common to all of them. Because the mass of any of these satellites is negligible compared to Earth's mass, there is no reason to include it in any of our computations. For brevity, we will refer to quantities such as angular momentum and torque, but we are actually using those quantities divided by mass, known as the specific angular momentum, specific torque, etc.

The position and velocity vectors are used to calculate the angular momentum (i.e., $r →× v →$). If angular momentum vs. time is plotted and the y axis includes zero, the line looks horizontal, apparently confirming Kepler's second law. Zooming in on the apparently horizontal line, however, makes it clear that the second law does not hold for these data.

Beyond the introductory level, close inspection of plots of the three components of angular momentum reveals periodic variations. The absolute size of Lz (per kg) is in the range of 1010–1011 m2/s, so it is more useful to calculate the torque from the difference in values of Lz from one second to the next. We would expect the torque to equal zero for the two-body Kepler problem where two perfect spheres are in mutual orbit, and the z components of the calculated torques are of order 10–100 m2/s2, meaning second-to-second changes in Lz of approximately one part per billion. This suggests a numerical precision issue, but there's more to the story. The relative importance of different effects can be seen in Table I (taken from Table 3.4 in Ref. 22).

Table I.

The gravity produced by a point-mass Earth is of course the dominant force acting on a GPS satellite, followed by the quadrupole term at about 0.011% of the monopole term. Solar/lunar effects are a further order of magnitude smaller than the quadrupole term, and radiation pressure is about 2% of the solar/lunar gravity. It is noted that the next individual harmonic term, either the J3 zonal or the leading tesseral or sectorial terms, would produce an acceleration of approximately 10−8 m/s2.

Perturbation Acceleration (m/s2)
Monopole (Kepler)  0.56
Further harmonics  3 × 10−7
Solar/lunar gravity  5 × 10−6
Body tides  10−9
Ocean tides  10−9
Earth albedo  10−9
Perturbation Acceleration (m/s2)
Monopole (Kepler)  0.56
Further harmonics  3 × 10−7
Solar/lunar gravity  5 × 10−6
Body tides  10−9
Ocean tides  10−9
Earth albedo  10−9
For the ideal two-body case mentioned above, we should get zero, or nearly zero for the torques in the x and y directions as well. If we calculate these torques, we see values that are larger than those in the z direction by a factor of 10. This is the first glimpse into the perturbations affecting satellites. The most significant additional force acting on the GPS satellites is the Earth's gravitational quadrupole moment resulting from its equatorial bulge. The acceleration of the satellite caused by this oblateness is modeled in the GPS-IS (and derived in the Appendix) as
$a i=F 1 − 5 z / r 2 i / r i=x or y,$
(2a)
$a z=F 3 − 5 z / r 2 z / r,$
(2b)
where
$F=− 3 2 J 2 μ r 2 R E r 2,$
(3)
and J2 = 0.001 082 626 2. The parameter μ is the product of Newton's universal gravitational constant and the Earth's mass (3.986 00 × 1014 m3/s2), and r is the distance from the center of the Earth to the satellite. For the satellite in the PRN 2 slot on 5/1/22, the plots in Fig. 2 show the total change in each component of angular momentum per second and the torque caused by the quadrupole (J2) term. We note that the monopole term produces the familiar 1/r2 acceleration. The dipole term, which would scale as 1/r3, vanishes because the center of mass coincides with the center of our coordinate system. The acceleration caused by the quadrupole term therefore drops off as 1/r4, as expected.
Fig. 2.

Total torque (red disks) acting on satellite PRN 2 on 5/1/22 and torque generated by the Earth's quadrupole moment (blue triangles) on the same satellite at the same time. The three plots are for the x (leftmost), y (center), and z (rightmost) components of these torques.

Fig. 2.

Total torque (red disks) acting on satellite PRN 2 on 5/1/22 and torque generated by the Earth's quadrupole moment (blue triangles) on the same satellite at the same time. The three plots are for the x (leftmost), y (center), and z (rightmost) components of these torques.

Close modal

The quadrupole torque in the z direction vanishes due to the axial symmetry of the equatorial bulge. While a detailed discussion of quadrupole moments might not be commonplace in sophomore or junior level mechanics courses, it does provide an interesting look at the most significant correction to the Kepler problem. Introducing the monopole, dipole, and quadrupole terms in the mass distribution as analogous to the concepts of mass, center of mass, and moment of inertia from first semester physics may be worthwhile.

Examining dL/dt minus the quadrupole torque, we notice a periodicity of roughly 6 hours (half the orbital period of a GPS satellite), suggesting a force that acts in the same direction with the same approximate magnitude at opposite parts of the orbit. This is reminiscent of solar/lunar tidal influences, where the tidal period is roughly one half of the Earth's rotational period. Seeber22 outlines the form of this perturbation, pointing out that we require the acceleration of the satellite due to the Moon less the acceleration of the Earth due to the Moon. This yields
$r ¨=G m M r M − r r M − r 3 − r M r M 3,$
(4)
where mM is the Moon's mass, rM is the Earth–Moon distance, and r is the Earth-satellite distance, with a similar expression for the Sun. Table I shows that the influence of the Sun and Moon produces an acceleration approximately 10% as great as that due to Earth's oblateness.

To calculate the torque due to solar/lunar influences, a Python script queries the JPL Horizons database,23 which provides precise position information about a variety of natural and man-made bodies, including the Sun and Moon. The positions and velocities are downloaded at one-hour intervals and interpolated for times between these points. Accelerations and torques are found by the Python script. Comparing the total solar/lunar torque to the residual torque acting on a satellite after the removal of the quadrupole term (Fig. 3) shows significant agreement.

Fig. 3.

Residual torque (red disks) acting on satellite PRN 2 on 5/1/22 after subtraction of the torque generated by the Earth's quadrupole moment and the combined torque exerted by the Sun and Moon (blue triangles) on the same satellite at the same time. The three plots are for the x (leftmost), y (center), and z (rightmost) components of these torques.

Fig. 3.

Residual torque (red disks) acting on satellite PRN 2 on 5/1/22 after subtraction of the torque generated by the Earth's quadrupole moment and the combined torque exerted by the Sun and Moon (blue triangles) on the same satellite at the same time. The three plots are for the x (leftmost), y (center), and z (rightmost) components of these torques.

Close modal

If we remove both the J2 and the solar/lunar torques (Fig. 4), the residuals still seem suggestive of periodic phenomena. It would seem to be worthwhile to continue to examine higher order multipole terms since the reduction in strength with each successive order is proportional to RE/r, which is approximately ¼. When we look up the value of the octupole (J3) coefficient, however, we find that it is 0.23% of the quadrupole coefficient, yielding octupole torques approximately 2000 times smaller than the quadrupole torques.

Fig. 4.

Residual torque acting on satellite PRN 2 on 5/1/22 after subtraction of both the torque generated by the Earth's quadrupole moment and the combined torque exerted by the Sun and Moon on the same satellite at the same time. The lines are for the x (red disks), y (blue triangles), and z (green crosses) components of this torque.

Fig. 4.

Residual torque acting on satellite PRN 2 on 5/1/22 after subtraction of both the torque generated by the Earth's quadrupole moment and the combined torque exerted by the Sun and Moon on the same satellite at the same time. The lines are for the x (red disks), y (blue triangles), and z (green crosses) components of this torque.

Close modal
Just as the vector product of the satellite's radius and various external forces produces torques that we can compare to the total torque experienced by the satellite, we can calculate the total mechanical energy of a satellite every second and compare its change to the work done by those same external forces. The total energy (scaled to an initial value of 0), as computed from the following equation:
$E m= 1 2 v 2− G M E r,$
(5)
is shown in Fig. 5 along with the potential energy associated with the quadrupole moment J2.
Fig. 5.

Kinetic and monopole potential energy of satellite PRN15 on 8/15/22 scaled to initial value of zero (green, solid) compared to work done by gravitational force associated with the Earth's quadrupole moment (blue, dashed). The energy was scaled by subtracting -GME/(2 aREF) where aREF is the nominal semimajor axis for GPS satellites given in the GPS-IS and is 26 559 710 m.

Fig. 5.

Kinetic and monopole potential energy of satellite PRN15 on 8/15/22 scaled to initial value of zero (green, solid) compared to work done by gravitational force associated with the Earth's quadrupole moment (blue, dashed). The energy was scaled by subtracting -GME/(2 aREF) where aREF is the nominal semimajor axis for GPS satellites given in the GPS-IS and is 26 559 710 m.

Close modal

Including the J2 potential energy, variations are now much smaller and reveal the next higher order contribution to be the gravitational influence of the Sun and Moon, shown in Fig. 6. After incorporating the solar/lunar work into the kinetic, monopole potential and J2 potential energies, the variations in the total energy have dropped significantly. The wave shape is now asymmetric, but still represents oscillations with an approximate amplitude of 20 J/kg. According to Ref. 22, Table 3.4, the largest influences after the quadrupole and solar/lunar contributions are the higher-order multipole moments of Earth's mass distribution (producing accelerations of order 3 × 10−7 m/s2) and the solar radiation pressure (accelerations of order 10−7 m/s2).

Fig. 6.

Kinetic, monopole, and J2 potential energies (green, solid) of satellite PRN15 on 8/15/22 compared to work done by Sun and Moon (blue, dashed). A constant 35 J/kg has been subtracted from the solar/lunar work done to center the plots.

Fig. 6.

Kinetic, monopole, and J2 potential energies (green, solid) of satellite PRN15 on 8/15/22 compared to work done by Sun and Moon (blue, dashed). A constant 35 J/kg has been subtracted from the solar/lunar work done to center the plots.

Close modal

Modeling the radiation pressure is more complicated due to its dependence on the satellite's orientation, surface area, mass, and reflectivity. In Fig. 7, we have somewhat arbitrarily chosen a reflection coefficient of 2 and a mass of 1000 kg as reasonable estimates. After evaluating several values for area in the 10 to 50 m2 range, we chose 30 m2 as a reasonable compromise that provided a minimum of variation in the total energy when combined with the other terms. Examination of the potential energy per mass associated with the octupole (J3) term shows it to be approximately 1% of the size of the remaining variation in energy.

Fig. 7.

Kinetic, monopole, and J2 potential energies less solar and lunar work (green, solid) compared to work done by light pressure (blue, dashed). The light pressure model assumes a satellite with 30 m2 area, 1000 kg mass, and a reflectivity of 2. A constant 40 J/kg has been added to the green (solid) line to center the plots.

Fig. 7.

Kinetic, monopole, and J2 potential energies less solar and lunar work (green, solid) compared to work done by light pressure (blue, dashed). The light pressure model assumes a satellite with 30 m2 area, 1000 kg mass, and a reflectivity of 2. A constant 40 J/kg has been added to the green (solid) line to center the plots.

Close modal

Numerical integration of successively more accurate models is an interesting project for a sophomore- or junior-level course. In this case, we can take advantage of the SP3 precise orbit files with both position and velocity information.

We used the NDSolve package of Mathematica to numerically integrate three models over time, using the initial values of position and velocity from SP3 files. We use the data transformed to the ECI frame to remove inertial forces. First, we model the gravitational attraction between a spherical Earth and a particular GPS satellite. As shown in Fig. 8, the error when compared to the true position given in the SP3 file increases at a rate of about 0.5 km/h, though the rate of error increase shows large-amplitude oscillations with periods of approximately half of an orbit. Linear fits to the errors over a week yield values from approximately 0.4 to 1.6 km/h. Including the quadrupole term in the differential equation to be solved greatly improves this error, reducing the error over a week to between 10 and 200 m/h.

Fig. 8.

Deviations in position (left, linear; right, logarithmic) between the SP3 precise orbit data for PRN 07 during the week starting 5/1/2022 and Mathematica's numerical integration of the Keplerian approximation of two point masses interacting via Newton's law of gravity (blue, solid), the Keplerian approximation with the addition of the quadrupole term of the Earth's mass distribution (green, dotted), and the Keplerian approximation, the quadrupole term, and the gravitational attraction of the Sun and Moon (red, dashed). Solar and lunar positions were determined by an interpolating function of the JPL Horizons hourly position data.

Fig. 8.

Deviations in position (left, linear; right, logarithmic) between the SP3 precise orbit data for PRN 07 during the week starting 5/1/2022 and Mathematica's numerical integration of the Keplerian approximation of two point masses interacting via Newton's law of gravity (blue, solid), the Keplerian approximation with the addition of the quadrupole term of the Earth's mass distribution (green, dotted), and the Keplerian approximation, the quadrupole term, and the gravitational attraction of the Sun and Moon (red, dashed). Solar and lunar positions were determined by an interpolating function of the JPL Horizons hourly position data.

Close modal

When the gravitational accelerations produced by the Sun and Moon are included, the error over a week is reduced to between 2 and 80 m/h. Though we have not shown it in this plot, the difference in position calculated using the BRDC NAV data (as described in the GPS Interface Specification) as compared to the SP3 file data is generally between 1 and 2 m at all time scales.

A variety of other phenomena can be investigated using available GPS data. In the supplementary material,17 we discuss mining the data to reveal the precession of satellite orbits caused by the Earth's oblateness as well as solar/lunar influences. The arrangement of the GPS satellites in a set of six orbital planes is also readily demonstrated, as are the changes in the semi-major axis a of a given satellite on timescales ranging from hours to weeks. The occasional station-keeping changes in a are also clearly seen as nearly vertical jumps in a graph of a against time.

The flood of extremely high-precision satellite data allows student researchers to examine orbital mechanics at multiple levels of complexity. In addition to the standard exercise of verifying Kepler's three laws, the data are precise enough to see their limitations. The non-sphericity of the Earth is revealed, as well as the gravitational influences of the Sun and Moon on Earth satellites. A multitude of undergraduate research projects are suggested by the breadth and depth of the available data.

The authors have no conflicts to disclose.

1.
E.g.,
William T.
Skilling
and
Robert S.
Richardson
,
Astronomy
, Rev. ed. (
Henry Holt and Co
.,
New York, NY
,
1947
), p.
675
.
2.
Roger B.
Culver
, “
An experiment for studying Kepler's Laws and planetary masses from satellite motions
,”
Am. J. Phys.
39
(
11
),
1404
1406
(
1971
).
3.
Alan
Bates
, “
Galilean Moons, Kepler's Third Law, and the Mass of Jupiter
,”
Phys. Teach.
51
(
7
),
428
429
(
2013
).
4.
John
Dumar
, “
,”
Phys. Teach.
59
(
2
),
130
133
(
2021
).
5.
Herbert
Goldstein
,
Charles
Poole
, and
John
Safko
,
Classical Mechanics
, 3rd ed. (
,
Boston, MA
,
2000
), Chap. 3.
6.
Daniel
Kleppner
and
Robert
Kolenkow
,
An Introduction to Mechanics
, 2nd ed. (
Cambridge U. P
.,
Cambridge, UK
,
2014
), Chap. 10.
7.
John R.
Taylor
,
Classical Mechanics
(
University Science Books
,
Melville, NY
,
2005
), Chap. 8.
8.
Leon
Blitzer
, “
Lunar-solar perturbations of an Earth satellite
,”
Am. J. Phys.
27
(
9
),
634
645
(
1959
).
9.
Francisco G. M.
Orlando
,
C.
Farina
,
Carlos A. D.
Zarro
, and
P.
Terra
, “
Kepler's equation and some of its pearls
,”
Am. J. Phys.
86
(
11
),
849
858
(
2018
).
10.
Carolyn
Früh
and
Thomas
Schildknecht
, “
Accuracy of two-line-element data for geostationary and high-eccentricity orbits
,”
J. Guid. Control Dyn.
35
(
5
),
1483
1491
(
2012
).
11.
Danielle
Racelis
and
Mathieu
Joerger
, “
High-integrity TLE error models for MEO and GEO satellites
,” in
AIAA SPACE and Astronautics Forum and Exposition
(
American Institute of Aeronautics and Astronautics Inc
,
2018
).
12.
Delphine
Ly
,
Romain
Lucken
, and
Damien
Giolito
, “
Correcting TLEs at epoch: Application to the GPS constellation
,”
J. Space Safety Eng.
7
(
3
),
302
306
(
2020
).
13.
Orbit.svg
, see https://commons.wikimedia.org/w/index.php?title=File:Orbit.svg&oldid=652370411 for “
Wikimedia Commons
” (accessed April 27,
2022
)
14.
Jan
Van Sickle
,
GPS for Land Surveyors
, 4th ed. (
CRC Press
,
Boca Raton, FL
,
2015
).
15.
Carl E.
Mungan
, “
Relativistic effects on clocks aboard GPS satellites
,”
Phys. Teach.
44
(
7
),
424
425
(
2006
).
17.
https://www.sparkfun.com/products/18774 or any of the similar boards using ZED F9 parts.
18.
See supplementary material online for the Python scripts to gather and process the data, Mathematica notebooks to analyze it, and a brief discussion of the variable naming convention.
19.
C.
Noll
, “
The crustal dynamics data information system: A resource to support scientific analysis using space geodesy
,”
45
(
12
),
1421
1440
(
2010
); the data archive is available via the internet at http://cddis.nasa.gov/index.html.
20.
See https://geodesy.noaa.gov/corsdata/RINEX211.txt for “
File specification
.”
21.
The GNSS tab allows the selection of GPS Ephemeris Center of Mass data and a date range. The data are provided at 5-minute intervals. Note that the velocities are given in decimeters per second. https://earth-info.nga.mil/
22.
Günter
Seeber
,
Satellite Geodesy
, 2nd ed. (
Walter de Gruyter
,
Berlin
,
2003
).
23.
See https://ssd.jpl.nasa.gov/horizons/app.html#/ for “
Horizons System
.”