## Abstract

Most global navigation satellite systems (GNSSs) ephemeris representations require straightforward, albeit specialized algorithms to compute the transmitter position at a time of interest. As potential positioning, navigation, and timing (PNT) signal sources expand beyond medium Earth orbit, these representations must be modified to capture the dynamics of the host platforms. This work introduces the use of B-splines as a flexible framework to represent transmitter ephemerides that are applicable to any orbital or airborne regime and host platform. With this approach, the user equipment implements a simple, generic algorithm to compute transmitter positions from the B-spline representation that require no orbit or platform-specific models. Here we propose a B-spline ephemeris approach in which we compare the required navigation message length and fit accuracy to the legacy Global Positioning System (GPS) broadcast for use with medium Earth orbit transmission. We also demonstrate the applicability of this approach for a PNT satellite in low Earth orbit.

## 1 INTRODUCTION

All global navigation satellite systems (GNSSs) include a broadcast ephemeris within their navigation data that allows users to compute satellite positions in real time. This information is essential to support stand-alone point positioning. While the accuracy and precision of the satellite positions have a direct impact on the quality of the user’s position solution, the length of the ephemeris data affects both the time a receiver must stay locked onto the signals and the user’s time-to-first-fix. Current GNSS broadcast ephemeris parameters are carefully selected to balance sufficient accuracy for position solutions and compactness for rapid reception by users. However, it is important to note that these specific parameters vary across different constellations. Each GNSS constellation has unique algorithms and physical constants specified in their respective interface control documents that are implemented in the user equipment to compute satellite positions from the broadcast parameters.

Because most of the current GNSS signal sources are in medium Earth orbit (MEO), the broadcast parameters have been designed to accommodate the dominant Keplerian orbital motion and the perturbing forces experienced in MEO, such as solar radiation pressure (Vallado & McClain, 2001). However, these MEO-specific broadcast parameters do not accurately represent alternative orbital regimes such as geostationary orbit (GEO), inclined geosynchronous orbit (IGSO), and low Earth orbit (LEO). This has prompted the expansion of broadcast ephemeris parameter sets to produce satellite positions with comparable accuracy. For example, Montenbruck & Gill (2000) and Xie et al. (2018) have proposed additional ephemeris parameters that can address singularities that may arise when using Keplerian orbital elements with satellites in equatorial or polar inclinations or when eccentricity approaches zero. It is more challenging to represent satellite positions in LEO because of increased secular and harmonic gravitational perturbations. To adapt the global positioning system (GPS) ephemerides representation for LEO, Meng et al. (2021) proposed six additional parameters to tackle singularities and account for perturbations that resulted in a fit accuracy of 0.08 m for 20-minute intervals. More recently, Guo et al. (2022) proposed a similar parameter set for LEO augmentation systems that is slightly more compact (i.e., contained one less parameter) than the one proposed by Meng et al. (2021) that achieved an improved accuracy of 0.01 m for 20-minute intervals. However, both studies note that the fit accuracy declines for longer intervals and neither provides estimates of the bit requirements for implementation of these models in a navigation message.

Given the small number of distinct GNSS orbit classes that are currently available, it is not too difficult to adapt user equipment to support these newly-established orbit representations. However, as the variety of constellations and augmentations capable of transmitting GNSS signals expands, it will become increasingly challenging to manage multiple classes of broadcast ephemeris parameters and algorithms. Thus, we offer a flexible alternative that can convey transmitter positions with the same accuracy as specialized orbital methods and can be implemented in user equipment with a common, well-established, and computationally efficient basis (B)-spline linear recursive calculation for all potential orbital regimes.

The remainder of this paper begins with a review of current navigation message generation by GNSS constellation operators and end-user implementations. The next section details how these functions could be implemented for the alternative B-spline representation and includes the corresponding calculations in the user equipment. Then, to assess the suitability of the B-spline approach, we compare the number of bits that it requires to existing broadcast messages and its accuracy in representing MEO and LEO orbits using precise ephemeris records.

## 2 GNSS CONTROL SEGMENT OPERATIONS AND USER IMPLEMENTATIONS

Broadcast parameters vary among the different GNSS constellations; their accuracy and update rates have an impact on the precision of the positioning solutions. Table 1 provides a summary of the characteristics of current GNSS constellations, including their respective broadcast ephemeris parameters, update periods, and message sizes in terms of the number of bits. Most GNSS constellations use a modified set of Keplerian orbital elements in the broadcast ephemeris to describe satellite positions. The Russian GLONASS constellation is the exception because it transmits Cartesian coordinates for position, velocity, and lunisolar acceleration.

Each GNSS constellation has a control segment that consists of a network of ground-based control stations to monitor and track the satellites and a central processing facility to produce high-fidelity estimates of the satellite orbits. For example, the GPS control segment uses satellite tracking data and Kalman filtering techniques to estimate key parameters such as the satellite’s position, velocity, solar radiation pressure, clock bias, clock drift, and clock drift rate which are all used to create precise orbit and clock predictions (Warren & Raquet, 2003). The ephemeris parameters are then fit to the predictions for each satellite and uploaded to the satellites by the control segment for transmission to users.

On the user implementation side, the navigation message is received and decoded according to the specifications described in the respective GNSS interface control document. GPS, for example, provides scale factors, unit conversions, and specific physical constants such as the Earth’s rotation rate and gravitational parameters that must be used with the broadcast ephemeris parameters. The decoded parameters are then used in a series of equations to convert the Keplerian orbital elements and correction parameters to satellite positions, velocities, and accelerations in an Earth-centered, Earth-fixed (ECEF) coordinate frame (U.S. Government, 2022a). A similar process is performed for each GNSS constellation with their respective broadcast parameters and user algorithms.

## 3 B-SPLINE BROADCAST EPHEMERIS FORMULATION

The process of implementing a B-spline ephemeris representation would follow a process that is similar to the one currently in use. The control segment would continue to carry out the standard orbit prediction process, which is the main driver of the accuracy of the ephemerides. However, a B-spline curve is fit to the predicted orbit positions rather than fitting the constellation-specific ephemeris parameters. The B-spline representation is simply an alternative method of providing the predicted orbit information to users and does not fundamentally alter the prediction itself. The following sections define the B-spline curve, the fitting process, and the user implementation required to retrieve the satellite positions.

### 3.1 B-Spline Curve Definition

A B-spline is a curve representation commonly used in computer graphics, computer-aided design (CAD), and numerical analysis. It is a type of piecewise polynomial function constructed by combining B-spline basis functions. The pieces are defined by knots that control the shape and overall behavior of the curve. Mathematically, the equation for a B-spline curve **C**(*t*) of order *p* and polynomial pieces *n* is given by Equation (1):

1
where *N _{j,p}*(

*t*,

**u**) represents the

*j*B-spline basis function of order

^{th}*p*, which is defined on a non-decreasing, sequential knot vector

**u**. The spline coefficients

*b*(

_{j}*j*= 0, …,

*p*+

*n*-2) determine the weights applied to each basis function to generate the final curve. In general, there are

*p*+

*n*−1 coefficients for the curve. The specific values of the coefficients are determined by the curve fitting process which we described in greater detail in Section 3.2.

The knot vector **u** divides the domain of *t* into a series of intervals or knot spans, in which each knot span corresponds to a polynomial piece and the knots indicate where the pieces meet. The knot vector contains *n*+1 unique values, with the first and last values “clamped” such that each is repeated *p* times, where *p* is the order of the B-spline curve, as indicated in Equations (2) and (3):

2

3

Finally, the B-spline basis functions can be computed recursively with the following definition from Piegl & Tiller (1995) given the coefficients *b _{j}* and knot clamped knot vector

**u**as indicated in Equations (4) and (5):

4

5

For satellite ephemeris representation, the parameters to be transmitted to users are the curve order *p*, the set of coefficients {*b*_{0}, … *b*_{p+n−2}}, and the *n*+1 unique elements of the knot vector **u**. Additionally, three sets of knots and coefficients are required to represent each coordinate direction of the satellite position. The coefficients have units of distance on the order of the transmitter orbit radius and the knots have units of time. For a broadcast ephemeris, the knots can represent the elapsed time with respect to the starting time of the week. Thus, for a two-hour fit interval, the knots are bounded by 0 and 7200 seconds (i.e., 0 and 120 minutes).

### 3.2 Curve Fitting

The objective when constructing a B-spline ephemeris representation is to identify the B-spline that both accurately fits the satellite positions and minimizes the total number of bits required to represent the curve. This will ensure that users will be able to obtain accurate position solutions while the message is kept compact. This equates to minimizing the number of B-spline parameters (order, number of pieces, and knots) required to generate the curve.

A common technique used to produce a fitted B-spline curve is least-squares minimization. This process requires us to specify the B-spline curve order *p* and the number of polynomial pieces *n* to approximate the curve and the data to be fit. The knot vector **u** is formed by using a knot placement method to determine where the pieces meet. The spline coefficients then are obtained by minimizing the least-squares error defined by Equation (6):

6
where *l* is the number of data points included in the data set and *x*(*t _{i}*) is the

*x, y*, or

*z*component of the satellite position at time

*. To solve the resulting minimization problem, Equation (6) is transformed into a linear system of equations with respect to the coefficients*

_{i}*b*. This process is performed separately for each position component.

_{j}In order to create a B-spline ephemeris representation that is both compact and accurate, careful consideration must be given to the selection of the appropriate order *p*, number of polynomial pieces *n*, and the knot placement (values in **u**). This optimization problem has been studied extensively and several solutions have been published (Sarkar & Menq, 1991; Laurent-Gengoux & Mekhilef, 1993; Park & Lee, 2007; Luo et al., 2022). Our approach to selecting the parameters is presented in Section 5.

### 3.3 User Evaluation

The user evaluation of the satellite positions from a set of B-spline parameters is a straightforward process that follows Equations (1), (4), and (5), also known as the de Boor algorithm. Instead of directly computing the B-spline functions, *N _{j,p}*(

*t*,

**u**), the de Boor algorithm recursively calculates the spline function,

*C*(

*t*), using an equivalent formula. A basic implementation is defined in Algorithm 1. MATLAB and several open-source libraries, including Python, provide built-in functions that require the curve order, knot vector, and coefficients as input to perform these computations effectively (MathWorks, 2022).

Unlike the equations implemented in the user equipment for current GNSSs, the B-spline satellite position evaluation does not require evaluation of Keplerian elements or knowledge of physical constants, such as the Earth’s rotation rate. The B-spline representation also avoids the singularities that can arise in specific cases using Keplerian orbit elements due to small eccentricity or polar and equatorial orbits and does not require additional terms for accurate representation of the more perturbed orbits. Instead, the B-spline representation accommodates the complexity of different orbital regimes such as LEO with higher-order basis functions and additional polynomial pieces.

## 4 B-SPLINE BROADCAST EPHEMERIS PARAMETERS

To illustrate the number of B-spline ephemeris parameters included in a navigation message, we consider the following simple example. Suppose there is a need to represent a 30-min arc for a satellite in a 1000 km circular orbit and that a 4th-order B-spline with six pieces (*p* =4 and *n* = 6) might generate sufficient fit accuracy. This choice of spline results in nine coefficients and seven unique knots per curve. To keep the ephemeris compact, only the unique knots are included in the parameter set; the clamped knots as described in Equation (2) are excluded. This requires users to pad the knot vector with a multiplicity equal to the curve order before evaluating the satellite position using the de Boor algorithm. The resulting ephemeris consists of a total of 27 coefficients, 21 knots, the B-spline curve order, and a reference time.

In the context of a navigation message, the coefficients (scalar multipliers applied to the basis functions) have units of distance on the order of the orbit radius. These coefficients can be positive or negative, and the number of bits needed to represent each coefficient depends on the range determined by the maximum orbit radius and scale factor to represent decimal values. Similarly, the number of bits required to represent the knots depends on the time range for which the ephemeris is valid and the respective scale factor. Equation (7) provides the computation used to determine the number of bits needed given the parameter range *R* and scale factor *α*. An additional bit is included to represent signed values as shown in Equation (7):

7

For simplicity, knots can be chosen with integer values, as they are a function of the times provided in the orbit prediction. For example, if the knots have units of seconds, the range for a 30-min arc is 1800 seconds. Furthermore, knots with units of minutes would generate a smaller range of required values and thus fewer bits allocated to each knot.

The accuracy of the B-spline curve representation is affected by the precision of the coefficients used. Our research revealed that utilizing coefficients with decimeter-level precision can yield sub-meter accuracy with B-splines. While many current GNSS ephemeris parameters are defined with scale factors that are powers of two and thus with representations that are even more compact, our study employs straightforward scale factors of 0.1 for coefficients in meters and 1 for knots in minutes to estimate bits. This approach provides flexibility in representing various orbit altitudes and arc lengths. However, additional optimization might result in scale factors that provide more compact representations while still accommodating different orbital regimes. Additional information on the number of bits and the corresponding range with our scale factor definition is provided in Table 2.

Number of Bits and Effective Range per Knot (a) and Per Coefficient (b) Assuming Scale Factors of 1 min and 0.1 m, respectively

## 5 EVALUATION OF B-SPLINE PARAMETER SELECTION AND ACCURACY

In this section, we outline the process used to develop, assess, and select the B-splines to demonstrate the effectiveness of this method for GNSS orbits.

### 5.1 B-Spline Generation

We employed the MATLAB curve fitting toolbox and utilized the function `spap2` to perform least-squares spline approximation to generate the B-splines (MathWorks, 2022). The inputs to this function include the number of polynomial pieces `n`, curve order `p`, time vector `t`, and satellite positions `x`(`t`) used as the reference for the B-spline fit. The output is a spline structure that includes the B-spline coefficients and knot values. As mentioned in Section 3.2, the selection of B-spline parameters, including order, number of basis functions, and knot locations, is important in fitting the orbit. Thus, we selected the default knot placement algorithm included in the function. Given the number of basis functions `n`, curve order `p`, time vector `t`, and positions `x`(`t`), an example call to this function follows `s` = `spap2`(`n, p, x, t`). The MATLAB function `c` = `fnval` (`s, t`) was used to evaluate the resulting B-spline. This function uses the spline object, `s`, returned from `spap2` and a query time vector `t` (which does not need to be the same time vector used for curve fitting) as input, and returns the B-spline values evaluated at each time in `t` using an efficient implementation of de Boor’s algorithm.

### 5.2 Fit Accuracy Evaluation

To assess the accuracy of the fitted curves, the satellite position error distribution and the effect of the resulting range error due to satellite position error were considered. The orbit-only contribution of the signal-in-space ranging error (SISRE_{orb}) introduced by Montenbruck et al. (2014) and used as described by Reid et al. (2018), Xie et al. (2018), and Meng et al. (2021) was then evaluated. However, this value typically includes errors due to both the ephemeris representation and the orbit prediction from the constellation service provider, with the latter being the dominant factor. For this study, we consider only the error introduced by the ephemeris representation. This distinction is denoted by URE_{fit} and is defined in Equation (8):

8
where Δ*R*, Δ*A*, and Δ*C* are the fit errors in radial error, along-track, and cross-track directions, respectively, and *w _{R}* and

*w*are weighting factors that are dependent on orbital altitude. The values for these weights are listed in Table 3 from Xie et al. (2018).

_{A,C}### 5.3 B-Spline Parameter Selection

After establishing a process for generating and evaluating B-splines, our focus then shifted to selecting the curve order and number of basis functions for an appropriate fit to the predicted orbit. The goal was to identify the B-spline parameters that generated curves meeting a specific URE_{fit} requirement based on precise satellite positions or predictions and a desired arc length or fit interval. Because there can be multiple B-spline parameter combinations that fulfill the accuracy requirement, we selected the B-spline with the parameters that required the minimum number of bits for representation with sufficient fit accuracy.

A nominal number of pieces and curve order were selected based on 24 hours of position data. The process began by initializing the curve order and number of pieces. The position data was then split into segments that corresponded to the desired fit interval; B-splines were generated from these segments. The fit accuracy over the entire data set was then evaluated and the number of pieces increased until the accuracy requirement was met. Once the requirement was satisfied, the number of bits was determined. This process was repeated with higher B-spline curve orders until the search space was exhausted. The B-spline parameters that involved the fewest bits and/or led to the best-fit accuracy were selected as the nominal parameters for a given orbit.

## 6 B-SPLINE PERFORMANCE EVALUATION

The quality of the B-spline representation was assessed using precise orbit data for both MEO and LEO satellites. For MEO, B-splines were generated using precise GPS satellite orbits to achieve at least 0.1 m fit accuracy. The distribution of the residuals was analyzed and the number of bits required was compared to the GPS broadcast ephemeris. For LEO, B-splines were created using accurate LEO satellite orbits to demonstrate that this method can be used to fit curves that are somewhat less smooth resulting from higher perturbations in LEO. The LEO B-splines were designed to achieve a maximum fit accuracy of 0.1 m and were generated for short- and long-term fit intervals. The number of bits required to maintain this accuracy for the various fit intervals was also evaluated.

### 6.1 MEO Application (GPS)

The average SISRE_{orb} values for the GPS Legacy Navigation (LNAV) and Civial Navigation (CNAV) broadcast ephemerides evaluated against the International GNSS Service (IGS) precise orbits are 0.23 and 0.28 m, respectively, and account for both the control segment’s orbit prediction and ephemeris representation errors (Montenbruck et al., 2014; Steigenberger et al., 2015; Wang et al., 2019). To compare the B-spline ephemeris on the same terms, it would be necessary to fit B-splines to the same control segment orbit predictions that were used to generate these broadcast parameters. However, because these predictions are not publicly available, B-spline curves are fit to precise GPS satellite orbit data. The resulting fit errors are intended to capture the contribution of the B-spline representation to the overall SISRE_{orb}.

To generate the B-splines curves used in this study, precise IGS GPS orbits were used. This information provided in .sp3 files has an estimated accuracy of 2.5 cm (Johnston et al., 2017). The .sp3 data format provides satellite positions at 15-minute intervals and requires trigonometric interpolation methods to produce satellite positions at more frequent intervals (Schenewerk, 2003). To assess B-spline fitting, a single day’s worth of .sp3 data is interpolated to a 60-sec time resolution and then split into two-hour fit intervals which are consistent with GPS LNAV broadcast ephemeris update period. B-splines were fit to the two-hour segments and the final parameters were selected as described in Section 5. To ensure comprehensive analysis, B-splines were evaluated using data from GPS satellites in each of the six orbital planes; similar performance was observed across all satellites. To maintain focus and clarity, we present the B-spline results using the precise orbit data from PRN 5, which is a GPS block IIRM satellite (SVN 50). Table 4 lists the orbit characteristics, GPS week number, and time of the week from which each of the corresponding .sp3 files was obtained.

Table 5 reports the selection of B-spline parameters generated from this study, together with the URE_{fit} values and the required number of bits for each. We note that certain higher-order B-splines required only a single polynomial piece; this includes the 8th-order curve, which also has the lowest bit count. While B-splines with a single piece are equivalent to a polynomial fit, the polynomial coefficients have widely varying orders of magnitude (with higher degree coefficients tending towards zero). By contrast, all B-spline coefficients are of the same order of magnitude. This will be advantageous for navigation messages because it permits each parameter to be decoded uniformly, with the same scaling for all coefficients, in contrast to polynomial coefficients, which require different scale factors for each coefficient. In the current study, we found that B-splines with the smallest URE_{fit} were frequently those based on a single basis function.

Additionally, we found that the B-spline representation requires at least twice as many bits as the GPS LNAV broadcast ephemeris. While the LNAV parameters are created to take advantage of the orbital dynamics and specialized user algorithms that permit them to be compact, the B-spline representation is more general with a simple user implementation, albeit with an increase in the number of bits required.

Figure 1 presents an example of the radial, along-track, and cross-track errors fit over a two-hour period for the B-splines as reported in Table 5. As the curve order increases, we see that the residuals decrease; this finding is confirmed by the URE_{fit} values shown in Table 5. Additionally, B-splines generated with a curve order greater than nine provided no significant improvement in fit accuracy. We note that the residuals tend to be the largest near the edges of the fit interval. This occurs because fewer basis functions contribute to the curve’s shape in these regions. Nonetheless, the B-spline with the lowest URE_{fit} shown in Figure 1(d) does not exceed 5 cm residuals with equally small discontinuities between fit segments. This finding is highlighted in Figure 2(a).

Figure 2 depicts the radial, along-track, cross-track, and resulting URE_{fit} errors for two-hour fit intervals over 24 hours using a 9th-order B-spline curve. The previously described boundary effects are most visible in Figure 2(b) and appear between fit intervals. Overall, the URE_{fit} values rarely exceed 0.02 m. Figure 3 illustrates the corresponding distribution of errors in the radial, along-track, and cross-track directions. As shown, most of the error values are small and clustered near zero.

### 6.2 LEO Application

To assess the effectiveness of using B-splines to represent potential LEO transmitter positions, we used precise ephemerides from two LEO satellites, GRACE-B and Jason-2. Both data sets in .sp3 format were obtained from the Crustal Dynamics Data Information System (CDDIS), which is the National Aeronautics and Space Administration (NASA) archive of space geodesy data (Noll, 2010). Table 6 provides a summary of the orbital parameters for both satellites.

Similar to the methods described in the previous section, a single day’s worth of data was interpolated to 60-sec time resolution for consistency with the fitting methods described by Meng et al. (2021) and Guo et al. (2022). B-splines were evaluated for both short (20- and 30-min) and long (45- and 90- min) term fit intervals for the LEO satellites. The 20- and 30-min fit intervals were chosen because they bound the amount of time that an LEO satellite pass might be visible to a ground user. The 45- and 90-min fit intervals were evaluated to assess the use of the B-spline ephemeris scale to longer fit intervals. To generate the curves, the interpolated data for each LEO satellite were segmented into their respective fit intervals and the B-spline parameters were selected as described in Section 5. Additionally, a maximum URE_{fit} of 0.1 m was enforced when selecting B-spline fits, consistent with previous work on LEO ephemeris representations (Meng et al., 2021; Guo et al., 2022).

Results from previous studies by Meng et al. (2021) and Guo et al. (2022) revealed that the fit accuracy of their Keplerian-based representations degraded as the fit intervals increased. As shown in Tables 7 and 8, the accuracy of the B-spline representation is clearly not limited by the length of the fit interval. As shown, we can seamlessly increase the curve order and/or the number of polynomial pieces to maintain the desired fit accuracy. Although the B-spline representation provides us with more design flexibility with respect to the fit interval, the cost is the number of bits required. Additionally, while Meng et al. (2021) and Guo et al. (2022) did not provide an estimate on the number of bits required for their proposed ephemeris representations, we anticipate that the number of bits required for the B-spline representation for short-term (20-min) fit intervals is likely to be comparable.

## 7. CONCLUSION

The expansion of alternative sources of PNT signals to different orbit regimes provides us with an opportunity to rethink transmitter ephemeris representations. This study explores the use of B-spline representations as a flexible alternative to represent transmitter ephemerides that can capture a variety of satellite orbits and fit intervals.

The feasibility of using B-splines to represent satellite orbit ephemerides has been evaluated using precise data from MEO and LEO orbits. The results reveal that B-splines can successfully represent GPS satellite orbits with a fit accuracy of less than 0.1 m for two-hour fit intervals with only a modest increase in bits compared to the current LNAV broadcast ephemeris. LEO orbits, which are more challenging to represent, were also successfully represented for short- and long-term fit intervals with a fit accuracy of less than 0.1 m.

Overall, B-splines offer an alternative to the current and proposed GNSS broadcast ephemeris parameters that are tailored to specific orbital regimes. GNSS operators can design B-splines that accommodate various PNT system requirements by controlling the fit error threshold. With this approach, user implementations will not require constellation-specific algorithms to convert satellite ephemeris parameters into positions. Instead, the B-spline ephemeris representation involves a simple user algorithm that converts parameters to positions with the same efficiency regardless of orbit altitude or fit interval. Thus, this method offers the flexibility needed to represent any orbital regime.

In summary, the use of B-spline representations has the potential to improve the accuracy and flexibility of satellite orbit ephemerides for a significant range of PNT applications.

## HOW TO CITE THIS ARTICLE

Dobbin, M., & Axelrad, P. (2023). A Flexible Ephemeris Representation for GNSS and Alternative PNT Signal Sources Using B-Splines. *NAVIGATION, 70*(4). https://doi.org/10.33012/navi.610

## FINANCIAL DISCLOSURE

This material is partially based upon work supported by the Air Force Research Laboratory (AFRL) under Contract No. FA9453-20-C-2000.

## CONFLICT OF INTEREST

The authors declare no potential conflict of interest.

## ACKNOWLEDGMENTS

The authors would like to thank the anonymous reviewers for their valuable comments and recommendations.

The views expressed are those of the authors and do not reflect the official guidance or position of the United States Government, the Department of Defense, or of the United States Air Force.

## A SIMPLE B-SPLINE EXAMPLE

The following section provides a simple B-spline example that will permit users to recreate and compare values with. The corresponding code base is available in MATLAB.

In this example, we consider a single dimension of a circular, two-body orbit in an Earth-centered, Earth-fixed (ECEF) coordinate system. Truncated values are used so that this example will be simple for users to reproduce. The curve for a 500 km orbit is modeled as follows by Equation (A1):

A1
where *a* is the semi-major axis of the orbit (6378 + 500 km), *n* is the mean motion with *μ _{Earth}* = 3.986

*E5km*

^{3}/

*s*

^{2}, and

*t*is the elapsed time from the starting epoch.

To compare to the GRACE-B and Jason-2 LEO B-splines, we fit an 8th-order, single-piece B-spline to a 20-min fit interval (see Tables 7 and 8). The sine curve is generated using a time vector that spans from [0, 1200] seconds (secs) with 60-sec time resolution. Using MATLAB’s curve fitting toolbox and the function “spap2.m”, a single-piece 8th order B-spline is fit to the curve with position and time in units of m and min, respectively. Because this is a single piece B-spline, there are only two knots to indicate the start and end times of 0 and 20 min, respectively. The B-spline coefficients are listed in Table A1. In accordance with the scale factor used in this study, the coefficients are rounded to 1/10^{th} m. The splines are evaluated using the function “fnval.m” in MATLAB and evaluated every second. Figure A1 shows the resulting B-spline and fit error that results from the use of these knots and coefficients. The maximum fit error for the B-spline fit is approximately 5 cm with an average fit error of 2 cm.

It is possible to compare the single-piece B-spline fit with a polynomial fit of the same degree. The fit error is the same without any rounding or truncation of the B-spline coefficients. The primary difference between the B-spline and polynomial fit is that the polynomial fit is poorly conditioned at high orders. Table A1 also shows the polynomial coefficients for an equivalent fit. The polynomial coefficients were rounded to keep the fit error within 10 cm. Compared to the B-spline with coefficients of the same order of magnitude, the magnitude of the polynomial coefficients varies significantly. Additionally, the B-spline coefficients all have consistent units of position. A comparison of the error profiles of the B-spline and polynomial fits is shown in Figure A2.

## Footnotes

0 Approved for public release; distribution is unlimited. Public Affairs release approval # AFRL-2022-4007

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.