The recent availability of relatively inexpensive dualfrequency receivers for signals from Global Navigation Satellite Systems (GNSS) provides access to ultraprecise, realtime 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 higherorder 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.
I. INTRODUCTION
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 exercise^{2–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 realtime 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.
II. KEPLER'S LAWS
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 textbooks^{5–7} as well as this journal.^{8,9} We state the generalized versions of Kepler's laws as follows:

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

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

The ratio of the cube of the semimajor 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.
III. SATELLITE CHOICE
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 twoline elements. This specification dates to the 1960s and includes the six Keplerian elements. The semimajor 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}
The methods in this paper are applicable to any satellites, but different orbital altitudes can be expected to provide different modeling challenges. Satellites in lowEarth 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 pointmass 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 realtime positional accuracy to within a meter or two. The results can also be postprocessed using preciseorbit files (SP3 files) typically accurate to about a centimeter. Access to data that are correct at (or below) the partsperbillion 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 welldocumented 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 nonsphericity 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.
IV. DATA SOURCES AND FORMAT
The operating principle behind the various GNSS constellations^{14} is the simultaneous reception of timecoded 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 timeshifting 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 pseudorange. 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 300bit (six second) subframe page are detailed in the GPS Interface Specification (GPSIS)^{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 Ublox ZEDF9T on a sparkfun.com board^{17} and a multiband 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 socalled raw data as well as subframe data over several days and use Python scripts^{18} to process these recordings and produce files of commaseparated values (CSV files), which can then be analyzed with Excel or other software.
The raw data include Doppler shift and pseudorange measurements, as well information about the satellite and receiver clocks that allows the pseudorange 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 GPSIS 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 semimajor 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 nonspherical 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 bestfit 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 pseudorange and Doppler shift information, the financial expense can be eliminated, as the data broadcast by the GPS satellites are also available from NASA^{19} as NAV files in RINEX^{20} 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 CDDIS^{19} represent a compilation of broadcast data collected by sites all around the globe, providing information for the entire day. We have Python scripts^{18} 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 GeospatialIntelligence Agency^{21} provides these files.
One of the quirks of the GPS system is the use of an Earthcentered, Earthfixed (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 Earthcentered inertial (ECI) Cartesian frame.
Our ECEFECI 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 0^{h} 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 righthanded coordinate system.
V. INTRODUCTORY LEVEL EXERCISE
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 ephemerisderived files. For locally received files, we return Doppler shift, pseudorange, and carrier phase values as well as the data calculated as specified in the GPSIS. 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 semimajor 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 fiveminute increments throughout the day. With one of these files, the semimajor 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 semimajor 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.
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 \u2192\xd7 v \u2192$). 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.
VI. INTERMEDIATE LEVEL EXERCISE
Beyond the introductory level, close inspection of plots of the three components of angular momentum reveals periodic variations. The absolute size of L_{z} (per kg) is in the range of 10^{10}–10^{11} m^{2}/s, so it is more useful to calculate the torque from the difference in values of L_{z} from one second to the next. We would expect the torque to equal zero for the twobody Kepler problem where two perfect spheres are in mutual orbit, and the z components of the calculated torques are of order 10–100 m^{2}/s^{2}, meaning secondtosecond changes in L_{z} 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).
Perturbation .  Acceleration (m/s^{2}) . 

Monopole (Kepler)  0.56 
Quadrupole (J_{2})  5 × 10^{−5} 
Further harmonics  3 × 10^{−7} 
Solar/lunar gravity  5 × 10^{−6} 
Body tides  10^{−9} 
Ocean tides  10^{−9} 
Solar radiation pressure  10^{−7} 
Earth albedo  10^{−9} 
Perturbation .  Acceleration (m/s^{2}) . 

Monopole (Kepler)  0.56 
Quadrupole (J_{2})  5 × 10^{−5} 
Further harmonics  3 × 10^{−7} 
Solar/lunar gravity  5 × 10^{−6} 
Body tides  10^{−9} 
Ocean tides  10^{−9} 
Solar radiation pressure  10^{−7} 
Earth albedo  10^{−9} 
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.
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 manmade bodies, including the Sun and Moon. The positions and velocities are downloaded at onehour 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.
If we remove both the J_{2} 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 R_{E}/r, which is approximately ¼. When we look up the value of the octupole (J_{3}) coefficient, however, we find that it is 0.23% of the quadrupole coefficient, yielding octupole torques approximately 2000 times smaller than the quadrupole torques.
Including the J_{2} 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 J_{2} 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 higherorder multipole moments of Earth's mass distribution (producing accelerations of order 3 × 10^{−7} m/s^{2}) and the solar radiation pressure (accelerations of order 10^{−7} m/s^{2}).
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 m^{2} range, we chose 30 m^{2} 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 (J_{3}) term shows it to be approximately 1% of the size of the remaining variation in energy.
VII. MORE ADVANCED PROJECTS
Numerical integration of successively more accurate models is an interesting project for a sophomore or juniorlevel 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 largeamplitude 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.
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 semimajor axis a of a given satellite on timescales ranging from hours to weeks. The occasional stationkeeping changes in a are also clearly seen as nearly vertical jumps in a graph of a against time.
VIII. CONCLUSION
The flood of extremely highprecision 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 nonsphericity 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.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.