Demonstrating developments in high-fidelity analytical radiation force modelling methods for spacecraft with a new model for GPS IIR/IIR-M

This paper presents recently developed strategies for high-fidelity, analytical radiation force modelling for spacecraft. The performance of these modelling strategies is assessed using a new model for the Global Positioning System Block IIR and IIR-M spacecraft. The statistics of various orbit model parameters in a full orbit estimation process that uses tracking data from 100 stations are examined. Over the full year of 2016, considering all Block IIR and IIR-M satellites on orbit, introducing University College London’s grid-based model into the orbit determination process reduces mean 3-d orbit overlap values by 9% and the noise about the mean orbit overlap value by 4%, when comparing against orbits estimated using a simpler box-wing model of the spacecraft. Comparing with orbits produced using the extended Empirical CODE Orbit Model, we see decreases of 4% and 3% in the mean and the noise about the mean of the 3-d orbit overlap statistics, respectively. In orbit predictions over 14-day intervals, over the first day, we see smaller root-mean-square errors in the along-track and cross-track directions, but slightly larger errors in the radial direction. Over the 14th day, we see smaller errors in the radial and cross-track directions, but slightly larger errors in the along-track direction.

Since the 1980s, various methods for dealing with the problem have been presented in the literature (Colombo 1986;Beutler et al. 1994;Fliegel and Gallini 1996;Springer et al. 1999;Bar-Sever and Kuang 2003;Arnold et al. 2015). Many of these are empirical methods, requiring no a priori knowledge of the spacecraft properties or its operating environment. In global network analyses that incorporate tracking measurements from a large network of one hundred stations or more, such methods can produce spacecraft orbits with cm-level accuracy (Sośnica et al. 2015).
However, in a purely empirical approach, the orbit model parameters can absorb the effects of other un-modelled or mis-modelled processes (e.g. Earth rotation, geocentre variation (Meindl et al. 2013), etc.). This can result in orbit model parameter estimates that are non-physical, which means they cannot improve our understanding of the physical processes that determine the trajectory of the satellites and are therefore limited in their ability to help improve the modelling of those processes. As a result, a number of groups introduced analytical, or physics-based, radiation force modelling into their orbit estimation processes. In this area, the box-wing (BW) approach, first introduced by Marshall and Luthcke (1994) for application to POD of the TOPEX/Poseidon mission, has been particularly influential. The general concept is to model the spacecraft structure using eight flat plates (six for a cuboid representing the spacecraft bus and two for solar panels), with assumed values for the optical and thermal properties of the surfaces, which are then combined with a priori modelling of the spacecraft attitude and the incident radiation fluxes. This approach was applied to the Block II/IIA and Block IIR satellites of the Global Positioning System (GPS) by Rodriguez-Solano et al. (2012). In comparing the performance of their semi-analytical adjustable box-wing model with the Centre for Orbit Determination (CODE) Empirical Orbit Model (ECOM; Beutler et al. 1994), the authors determined that the orbit solutions produced by the two methods were comparable, but that the accelerations produced by the ECOM model were less physically meaningful. More recent efforts that adopt a broadly similar modelling approach include Montenbruck et al. (2015) for Galileo satellites and Montenbruck et al. (2017a) for the QZS-1 satellite of Quasi-Zenith Satellite System (QZSS).
The box-wing models are relatively easy to implement. However, they are not able to fully capture the radiation flux-spacecraft surface interaction in satellites with complex surface geometry where the effects of mutual self-shadowing and reflected radiation can be significant. An alternative class of analytical radiation force modelling methods, in which ray-tracing techniques are combined with detailed spacecraft surface models to account for SRP, Earth radiation pressure (ERP) and thermal re-radiation (TRR), was developed in the early 1990s (Klinkrad et al. 1991), and these methods are able to capture these detailed effects. The models were tested in POD of the European remote sensing satellites ERS-2 and ENVISAT (Doornbos et al. 2002). To distinguish them from the box-wing methods, we refer to these as high-fidelity analytical radiation force modelling methods. In GNSS, high-fidelity SRP modelling for GLONASS satellites was first explored by Ziebart and Dare (2001). This work was motivated by broader efforts to improve GLONASS orbit quality as part of the IGEX-98 campaign (Willis et al. 1999). Work in this area continued over the years at University College London (UCL), where the approach was enhanced with methods to account for TRR (Adhya 2005), ERP (Sibthorpe 2006;Ziebart et al. 2007;Li et al. 2017) and antenna thrust (AT; Ziebart et al. 2007), and validated on a number of additional cases including the GPS Block IIA and Block IIR satellites, the Jason-1 spacecraft of the Ocean Surface Topography Mission and ENVISAT (Ziebart et al. 2005;Sibthorpe 2006). Recent work in this area demonstrated improved accuracy in shadow modelling when using geometric primitives, as opposed to triangular tessellations, to represent curved surfaces when constructing the spacecraft model (Grey and Ziebart 2014). The modelling approach, as presented in Ziebart et al. (2005), was adopted into the operational standards for precise orbit determination of the Jason-1 altimetry satellite (Cerri et al. 2010;Zelensky et al. 2010). Recently, other research groups have explored a broadly similar approach for modelling SRP on Beidou satellites (Tan et al. 2016;Wang et al. 2018), on the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) satellite (Gini 2014) and on QZS-1 (Darugna et al. 2018).
In this paper, using a new model for the GPS Block IIR and IIR-M satellites, we present developments to the UCL model computation strategies that are designed to extend the validity of the models to all possible orientations of the radiation source(s) with respect to the spacecraft. As such, these nextgeneration general-purpose models account for the effect of radiation forcing from any number of radiation sources, from any direction, providing a high-fidelity radiation flux-spacecraft interaction model that can be used to deal with both SRP and ERP. A key advantage of this approach is that the final model makes no prior assumptions about the attitude characteristics of the spacecraft and can therefore deal with any deviations from nominal attitude.

UCL modelling strategy
Our radiation force modelling strategy comprises three processes: (i) Computation of the bus model, where the space vehicle bus contribution to the accelerations due to SRP and ERP is dealt with. In this process, the accelerations due to thermal emissions from the multi-layered insulation (MLI) covering the bus are also computed, according to Sect. 4.3.3 of Adhya (2005). The core technique uses a ray-tracing algorithm, where the rays simulate the incident radiation flux for a given geometry of the spacecraft with respect to the radiation source. The output of this process is a set of three grids representing the accelerations in the X, Y and Z-axes of the spacecraft body-fixed system (BFS), where the grid nodes are spaced at 1°intervals in latitude and longitude in the BFS. (ii) Separate computation of the solar panel model, where the solar panel contributions to the accelerations due to SRP and ERP are dealt with. (iii) AT modelling, which accounts for the recoil force on the spacecraft due to emission of photons from signal transmitters.
As input, the approach requires a computer model of the spacecraft that holds information about the external geometry and various surface material properties including reflectivity, specularity, absorptivity and emissivity. The models are built from a combination of geometric primitives (polygons, circles, cylinders, spheres, cones and truncated cones), avoiding any need for tessellation, especially on curved surfaces (Ziebart et al. 2003). This produces models with good geometric fidelity without requiring an excessively large number of components, e.g. the UCL model for the GPS IIR/IIR-M bus, as shown in Fig. 6, is made up of 182 components. The solar panels (not shown in Fig. 6) are modelled as two rectangular plates.
For solar flux, the models are computed using a nominal value for the mean solar irradiance at one astronomical unit (AU) of 1368 Wm −2 (Hastings and Garrett 1996). The solar irradiance is known to vary over the solar cycle (with a period of between 9 and 14 years) by 1.4 Wm −2 . This represents circa 0.1% variation in the parameter. Little is gained by correcting the nominal value. It is more important to scale the model depending on the probe-Sun distance at the calculation epoch. Taking 1368 Wm −2 as a reference value, the eccentricity of the Earth's orbit about the Sun modulates the solar irradiance near the Earth to 1415.7 Wm −2 at perihelion (+3.4%) and 1322.6 Wm −2 at aphelion (−3.3%). This gives a variation (between perihelion and aphelion) of circa 100 Wm −2 (the precise value being 91.3 Wm −2 ), approximately 6.7% of the mean value.
For the Earth radiation flux model, we use data from the Clouds and Earth Radiant Energy System (CERES) (Wielicki et al. 1996) project, which provides the irradiance at the top-of-atmosphere (TOA), an altitude of~30 km above the Earth's surface, in a grid format spaced at 1-degree intervals in latitude and longitude in an Earth-centred Earth-fixed (ECEF) system. Computationally, it can be expensive and slow to determine the total Earth radiation flux incident on a spacecraft, from the part of the Earth's surface that is visible to that spacecraft, based on a search of the full CERES grid. To overcome this, we have developed a configurable Earth radiation model that re-organises the CERES data into a grid of triangles wrapped around the TOA surface. The number of triangles used to represent the TOA surface is configured during run-time, based on the number of triangles required to achieve a specified precision level. This approach is outlined in Li et al. (2017). At GNSS altitudes, radiation flux from the Earth is about 15 Wm −2 .
In the ray-tracing algorithm, a pixel array (simulating the radiation source) is projected onto the computer simulation of the spacecraft with the force at each ray-surface intersection computed according to: where • F n is the normal force acting in the direction of the surface normal,n, • F s is the shear force acting in theŝ direction, which is along the projection of the total force onto the surface plane, • F mli is the force due to the thermal re-radiation from the MLI on the bus surface, which also acts along the normal direction, • E is the mean irradiance of the radiation source at one astronomical unit, • A is the area of the surface (determined in this case by the pixel array spacing), • c is the speed of light in vacuum, • ν is the reflectivity of the material, • μ is the specularity of the material, • θ is the angle of incidence of the radiation with respect to the surface, • T mli is the temperature of the MLI, • σ is the Stefan-Boltzmann constant, • α is the absorptivity of the MLI material, • mli is the emissivity of the MLI material, • eff is the effective emissivity between the MLI and the spacecraft and • T sc is the internal temperature of the spacecraft bus.
Note, Eqs. 1 and 2 are derived in Ziebart (2001Ziebart ( , 2004 and Eqs. 3 and 4 are developed in detail in Adhya (2005). The acceleration due to AT,ẍ at (r, t), is calculated according to: where W is the signal power in Watts, m is the spacecraft mass in kg andr is the unit vector from the geocentre towards the satellite centre of mass (Ziebart et al. 2007). Equations 1 and 2 are also used for computing the SRP and ERP forces acting on the solar panels, but the spacecraft bus and the solar array are treated separately during force model computation, with the results combined during model implementation, as explained in Sect. 4. A similar approach is used by Darugna et al. (2018), and this is done because it simplifies the model computation process as it is not always practical to incorporate the correct solar panel behaviour into the ray-tracing computations.

Bus model computation scheme
To produce the complete bus model, the pixel array is rotated around the spacecraft in a systematic way, through a discrete set of points, and the ray-tracing computations are performed from each point. Each computation takes the following form: where the inputs ϕ atan2 z, x 2 + y 2 and λ atan2 (y, x) represent latitude and longitude, respectively, of the radiation source in the spacecraft BFS (as defined in Fig. 1); the outputs a x , a y and a z are acceleration along the X, Y and Z-axes, respectively, in the spacecraft BFS.
In the computation scheme used in early analyses (Ziebart and Dare 2001), the pixel array was rotated around the spacecraft Y -axis at 12°intervals in the Earth-probe-Sun (EPS) angle, as shown in Fig. 2. The ray-tracing algorithm was executed at 30 points only, and a Fourier series was fitted to the results to represent the final continuous model. The underlying assumption was of spacecraft attitude behaviour fully consistent with the nominal attitude model, in which the Sun is confined to the XZ plane, and the Z-axis points to the geocentre (see Fig. 2). In subsequent studies, e.g. Ziebart et al. (2005), this modelling assumption was maintained but the number of data points computed was increased to 360 by reducing the increment in the EPS angle to 1°.
GNSS satellites on orbit may depart from a nominal attitude state due to limitations of their attitude control systems.  In nominal attitude mode, in the spacecraft BFS, the Sun is confined to the BFS X-Z plane and the spacecraft Z-axis points to the geocentre Non-GNSS satellites can have attitude laws that are far less constrained than the typical GNSS attitude laws. Thus, as the application of the core technique was considered for non-GNSS missions, the computation scheme was modified. First, the EPS-sweep pixel array orientation scheme was proposed. In this scheme, the pixel array centre points are uniformly distributed around the spacecraft in the central EPS plane (i.e. the spacecraft X-Z plane in this case). Then, each point is rotated by ± 1°about the spacecraft Y -axis, resulting in a new set of points, all inclined at ± 1°with respect to the spacecraft X-Z plane. This process is repeated at 1°steps to populate a full set of pixel-array centre points in 10°-arcs around the spacecraft as shown in Fig. 3. For GNSS satellites, the EPS-sweep method produces pixel arrays that are distributed around the primary parts of the spacecraft BFS within which the Sun moves, but coverage remains incomplete in other directions.
Therefore, an additional computation scheme based on the spiral points algorithm (Saff and Kuijlaars 1997) was introduced, see Fig. 4. With this method, it is possible to efficiently position the radiation source uniformly on a sphere that encloses the spacecraft. As it provides complete coverage in all directions, it can be used to produce a general-purpose model that makes no prior assumptions about the orientation of the spacecraft with respect to the radiation source(s). Currently, the spiral points computation scheme is our preferred model computation method, but there are additional data processing steps required before the outputs of a computation scheme based on this method can be used in a POD process.

Producing the grid files
The spiral points are not regularly spaced in latitude and longitude. Instead, the points are sorted according to distance along the spiral path, starting at the north pole (ϕ 90 • , λ Each point represents a pixel array centre point for a single instance of a radiation pressure acceleration computation and for the specific geometry of the radiation source with respect to the spacecraft 0 • ) and ending at the south pole (ϕ −90 • , λ 0 • ). This is not a standard method for organising the data. Thus, to provide a final model that is easily integrated into the POD processes of model users, we produce a set of acceleration grids, with grid nodes uniformly spaced at 1°intervals in latitude and longitude in the satellite frame.
To compute the grid values, we use a modified version of Shepard's method (Shepard 1968)-specifically, an implementation of the modified quadratic Shepard's method (Franke and Nielson 1980) with a type of full sector search as described in Renka (1988)-to determine the optimal set of gridding parameters. The interpolated values are computed in a two-step process. First, a quadratic surface is fitted around each data point. The quadratic (Q) neighbours parameter determines the radius of a circle large enough to include the nearest Q neighbours. Then, the interpolant at a chosen location is computed using an inverse distance weighted average of the computed quadratic surface fits around each data point. The weighting (W ) neighbour's parameter specifies the number of nearest data points to include for this. There are no clear rules for choosing either the Q or the W parameter for the modified Shepard's method, in that the optimal choice is data set specific. In this work, for the GPS IIR bus model, we developed a quality assurance process for determining this parameter pair using a two-dimensional search through Q-W space. This is how the process works: (i) The radiation pressure model is computed using the spiral points scheme and the EPS-sweep scheme. (ii) The 10,000 spiral points data set is expanded using a padding process, see below. (iii) Using modified Shepard's method, 1600 grids are produced from the output of the spiral points computation, for each acceleration component, with all combinations of Q, W pairs considered, where both Q and W range from 11 to 50. (iv) For each grid file, the interpolated values at each of the 3960 EPS-sweep points are calculated, and the interpolated value is compared with the results from the EPS-sweep computation. This is used to compute the RMS error value, E rms , for that grid file according to: where a EPS,i are accelerations at point i according to the EPS-sweep computation and a grid,i is an interpolated acceleration at point i derived from the grid file. (v) Finally, the grid files that minimise the E rms quantity for each component (X, Y and Z) are chosen as the optimal grid files for that spacecraft.

Padding the spiral points data set
The modified Shepard's method is a general-purpose interpolation algorithm that works with two-dimensional data that are irregularly scattered by using information from a spec-ified set of the nearest neighbour points. As such, in using this method to create the radiation force model grids, we encounter a problem. The algorithm is not able to identify the correct nearest neighbour points in the regions close to the data set boundaries (i.e. ϕ 90 • or − 90 • , λ 180 • or − 180 • ) when the output from a spiral points computation, labelled using latitude and longitude pairs in a 2-d Cartesian system, is provided as the input. To overcome this, we developed a method that creates an artificially extended spiral points data set that is bounded in the region ϕ ∈ (−270 • , +270 • ) and λ ∈ (−540 • , +540 • ). The transformation rules that map the raw spiral points data, bounded in the region ϕ ∈ (−90 • , +90 • ) and λ ∈ (−180 • , +180 • ), to data points in the extended regions are given in Eqs. 8 to 15, where f (ϕ, λ) is used to populate the extended region using the raw data. A portion of this extended data set around the north pole (ϕ 90 • ) is shown in Fig. 5, where the red data points are the spiral points. The yellow points, the top-padding above the north pole, are reflections of the raw data points about the ϕ 90 • line, which are then shifted by ±180 • in longitude beyond the north pole. Using this, expanded data set gives us a solution to the nearest neighbour problem. Another largely unavoidable issue is caused by the requirement to project the spiral points onto a 2-d Cartesian space, which distorts the apparent distance between points. With the Mercator projection, this effect increases with distance from ϕ 0 • . The impact of this is clearly seen in Fig. 5 where the density of data points becomes sparser approaching ϕ 90 • .

Modelling limitations
There are a number of factors limiting the accuracy of the current approach, which include: (a) Mis-modelling of reflected radiation coming off the bus onto the panels and shadowing from the bus onto the panels (and vice versa). Both of these are due to the separate treatment of the solar panels and the spacecraft bus during model computation. Here, there is a trade-off in modelling accuracy between being able to deal with non-standard solar panel orientations and being able to capture the effects of reflections and self-shadowing of the bus onto the panels. An analysis of this trade-off is not presented here but will be considered carefully in future development work. (b) No modelling of the time-evolution of the surface material properties. (c) Incomplete modelling of TRR effects. In the ray-tracing algorithm, we only consider spacecraft bus surfaces that are covered in multi-layer insulation (MLI). This strategy can perform reasonably well on those satellites where the surfaces are mostly covered in MLI, as is the case with the GPS Block IIR and IIR-M bus surfaces. However, it is limited in cases where a significant  (15), where the region labels are also defined proportion of the spacecraft surface is not covered in MLI (e.g. the SAR antenna on Sentinel-1; radiators on Galileo spacecraft, etc.). Also, we are not considering the force due to the temperature gradient across the Sunfacing and anti-Sun-facing sides of the solar panels in this study, but this effect has been considered in previous studies (Adhya 2005) and we are working towards developing a simplified approach to account for this. (d) No modelling of thermal recoil forces due to emissions from radiators and other thermal control system components that actively emit heat. The impact of this will be different between the Block IIR and the Block IIR-M satellites. The thermal control system of the Block IIR-M satellites was updated with additional integral heat pipes due to high heat concentrations in the honeycomb structure of the L-band panel due, in part, to increased signal power needs (Hartman et al. 2000).

Model implementation
The UCL radiation force model implementation requires several inputs. Most of these are spacecraft-specific information that includes position, nominal mass, actual mass if available, the grid files for the bus model, solar panel properties (area, surface material properties), attitude information (in the form of attitude control laws or on-board attitude measurements) to enable accurate determination of the spacecraft BFS and solar panel orientation in the BFS. As explained in Sect. 2, the model for the spacecraft bus is pre-computed, with the results of the computation stored in grids that are uniformly spaced at 1°intervals in latitude and longitude of the Sun position in the spacecraft BFS. To call these models in an orbit determination algorithm, these grids must be read in and stored in a suitable data structure. As a part of this process, it is a good idea to denormalise the grid values according to: where • m n is the nominal mass of the spacecraft, i.e. the value used to compute the grid in kg, • m a is the actual mass of the spacecraft in kg, •ẍ grid are the grid file accelerations in the spacecraft BFS x, y and z-axes in ms −2 , • ẍ grid are denormalised grid file accelerations in ms −2 , • E is the mean solar irradiance at 1 AU. This is because our radiation force modelling software was originally developed for solar radiation pressure modelling only. As such, the accelerations given in the grid files are produced using a solar radiation flux model that assumes a constant solar irradiance of 1368 Wm −2 at 1 AU. However, by applying this denormalisation step, it becomes relatively straightforward to use the UCL grids as a general-purpose radiation flux-spacecraft interaction model.
With all required inputs provided, and made accessible, it is possible to compute the accelerations due to the separate model components. The bus model is computed according to: where κ is the shadow crossing function (equals 1 in full phase of the Sun and 0 in umbra); E s (r, t) and E e (r, t) are solar radiation flux and Earth radiation flux, respectively, at the spacecraft's location r at time t; ϕ s and λ s are latitude and longitude, respectively, of the Sun's position in the spacecraft BFS; ϕ e and λ e are latitude and longitude, respectively, of the Earth's position in the spacecraft BFS. For latitude and longitude values between grid nodes, the accelerations should be calculated using bilinear interpolation. The solar panel's contribution to accelerations due to radiation forcing is: Finally, the combined acceleration due to radiation forces, x rad , is calculated according to: x rad (r, t) ẍ bus (r, t) +ẍ panel (r, t) +ẍ at (r, t), whereẍ at (r, t) is the acceleration due to antenna thrust.

The GPS IIR/IIR-M model description and data sources
The detailed UCL GPS IIR/IIR-M geometric model is generated from a set of technical drawings that are published in Chapter 5 of Adhya (2005). The primary source for the surface material properties is Fliegel and Gallini (1996). Additional details about how the model was put together are given in Ziebart et al. (2003). According to an unpublished report produced by UCL in collaboration with the Aerospace Corporation, and delivered to the United States Air Force in October 2005, the Block IIR/IIR-M satellites beyond GPS satellite vehicle number (SVN) 51 are equipped with a NAP ultra-high frequency (UHF) antenna (see Fig. 6), which is installed on the same side of the bus as the W-sensor high band antenna used for military applications found on the − X-face. Like the W-sensor high and low band antennae, the NAP UHF antenna is also composed of thin cylindrical components made of aluminium that are covered in black tape Values given to 2 decimal places (dp) (Adhya 2005). Thus, in our model, the same material properties are used, i.e. ν 0.06 and μ 0. The authors were unable to determine the full form on the NAP acronym or the purpose of this antenna.
The computation of the force models for the bus is performed using Version 5.05 of UCL's Analytical SRP and TRR Modelling Software at a nominal spacecraft mass of 1100 kg and a pixel-array resolution of 1 mm 2 . The bus model grid files for the IIR/IIR-M spacecraft, with and without the NAP antenna, are provided alongside this article as an electronic supplement.
Most of the values used for our solar panel model, given in Table 1, are taken from Adhya (2005). The combined surface area of the solar panel yoke arms is taken from Fliegel and Gallini (1996) because the drawings in Adhya (2005) provide only their length. For the rear side of the panels, we use surface properties given in Rodríguez-Solano (2009).
In Table 2, we present the statistics for the selected UCL grids for the GPS IIR/IIR-M satellites, both with and without the NAP antenna. The grids chosen are the ones corresponding to the Q, W parameter pairs that minimise the RMS error when the interpolated grid file values are compared against the results of an EPS-sweep computation. The RMS errors of the Z grids are approximately five times higher than the X grids. This is due to the W-band antennae and for those satellites that have them, the NAP antenna. In the EPS-sweep computation, these protruding elements result in significantly larger cross-section boundaries as the pixel array pans across the Z surfaces. By contrast, these elements have almost no effect on cross-section boundaries as the pixel array pans across the X surfaces. Thus, there are larger errors in the ray-tracing algorithm when computing Z accelerations. This is due to the edge-matching effect, which depends upon the cross-section perimeter and is explained in Chapter 10 of Ziebart (2001). This does not affect the Y grids as the pixel arrays do not pan across the Y surfaces in the same way. For the antenna thrust model, we use the IGS model values for signal power (http://acc.igs.org/orbits/thrust-power.txt), which are 85 W for GPS Block IIR and 108 W and 198 W for GPS Block IIR-M satellites.

Model validation
We investigate the performance of the new modelling strategy using two software systems: the UCL Orbit Dynamics Library (UCL-ODL) and ESOC's Navigation Package for Earth Observation Satellites (NAPEOS) software (Springer 2009). The UCL-ODL comprises a set of programs developed by researchers at UCL over the years, for the explicit purpose of studying the impact of force modelling strategies that are developed by the UCL Space Geodesy and Navigation Laboratory. NAPEOS is a GNSS data processing package developed by ESOC and used in its contributions to IGS activities to produce satellite orbits, precise clocks, station coordinates, Earth rotation parameters and so on.

Analysis of the impact of separate model components using the UCL-ODL
Using the UCL-ODL, we performed a series of sensitivity analyses to investigate the impact of the individual model components and verify the implementation method. In these tests, as the reference trajectory, we used precise IGS final orbits, considering all available IIR and IIR-M satellites over the full month of March 2016. For each satellite, we perform multiple orbit predictions, with separate prediction runs corresponding to separate IGS final orbit files. As such, in this part of the analysis, we consider 13 GPS IIR satellites and 7 GPS IIR-M satellites. For those satellites with a complete set of IGS final orbits during the analysis period, we perform 31 prediction runs from 1 to 31 March 2016. In the orbit propagator, the general force modelling strategy uses Earth Gravity Model 2008 up to degree and order 20 (Pavlis et al. 2012(Pavlis et al. , 2013 and the JPL Development Ephemerides 405 (DE405) for third-body gravitational forcing due to the Sun, Moon, Jupiter and Venus (Standish 1998). The solid Earth tide effect due to the Sun and the Moon is accounted for according to Marsh et al. (1987). General relativistic effects are modelled according to Sect. 3.7.3 of Montenbruck and Gill (2000). The numerical integration is based on an 8th order Runge-Kutta integrator.
In terms of radiation force modelling strategy, the following scenarios were systematically assessed: • Base model: SRP-only model using the ESOC BW model (Garcia-Serrano et al. 2016). • Test 1: SRP-only, where the bus model comprises grids produced by the UCL ray-tracing software, but only Eqs. 1 and 2 are used. • Test 2: Same as Test 1, but here the bus model comprises grids that account for both SRP and the effects of MLI TRR (Eqs. 3 and 4). • Test 3: Same as Test 2, but with ERP turned on.
• Test 4: Same as Test 3, but with AT turned on.
In Fig. 7, we show the impact of different modelling strategies on orbit prediction error over a single 12-h arc for the GPS satellite SVN 46. As smaller and smaller effects are considered in the modelling strategy, the orbit prediction results improve, giving a general indication that the models are performing as we expect. In Table 3, we provide orbit prediction errors statistics for all Block IIR/IIR-M satellites on orbit during the analysis period. The best results, in the sense that the RMS and the maximum 3-d orbit prediction error over a 12-h arc are minimised, at 0.648 m and 1.440 m, respectively, are produced by the method that considers the combined effects of SRP, bus MLI TRR, ERP and AT. For that modelling approach, the full set of statistics for all satellites that were considered in the analysis, are given in Table 4. An interesting observation from these results is that the modelling of the bus MLI TRR, an effect that is not considered by most IGS analysis centres, has a significant impact on reducing orbit prediction error over the arc.

Analysis of the impact of the new bus model on POD using NAPEOS
To assess the impact of introducing our grid-based model of the spacecraft bus on the quality of orbit estimates, we ran a number of POD analyses using NAPEOS. The analysis uses 100 tracking stations of the IGS Multi-GNSS Experiment (MGEX) (Montenbruck et al. 2017b) and all observed GPS satellites, but the results presented here focus on the 13 GPS Block IIR and 7 Block IIR-M satellites that were on orbit during the analysis period. The data processing method broadly follows ESOC's IGS analysis strategy (ftp:// igs.org/pub/center/analysis/esa.acn) where the basic observables are undifferenced carrier phases and pseudoranges and the integer carrier phase ambiguities are resolved (Ge et al. 2005). The Earth gravity model used is EIGEN-GL05C up to degree and order 12 (Foerste et al. 2008), and the JPL Development Ephemerides 405 (DE405) is used for third-body gravitational forcing due to the Sun, Moon and all solar system planets including Pluto (Standish 1998). The effects of solar Earth tides, ocean tides, solid Earth pole tide, oceanic pole tide and general relativistic corrections are accounted for according to the IERS conventions Fig. 7 Comparison of orbit prediction error over a single instance of a 12-h arc for the GPS IIR satellite SVN 46, using different modelling strategies  (Petit and Luzum 2010). The numerical integration uses the Adams-Bashforth/Adams-Moulton 8th order predictioncorrection multistep method, as described in Springer (2009). With the core data processing strategy fixed, we run the POD process using four different orbit modelling strategies, batch processed at 24-h intervals, from 00:00:00 to 23:59:30 in GPS time, thus completely independent from day to day. For the orbits, we generate estimates for the midnight epoch, such that there is an overlap between consecutive solutions at a single point. A full year (2016) is considered, so there are 366 independent solutions and 365 overlap points. In addition to the orbit model parameters, station coordinates and Earth rotation parameters are also estimated. The orbit models considered include: 1. ECOM: No a priori radiation force model, only the reduced ECOM (Springer et al. 1999) and three constrained along-track parameters (constant, cosine and sine with argument of latitude as argument). Here, the along-track parameters are included as soak-up parameters to absorb the effects of orbit mis-modelling, which tends to manifest strongly in the along-track direction, as the results of Sect. 5.1 demonstrate. 2. ECOM + BW: Same estimation strategy as ECOM-only, but here we also include an a priori radiation force model using the ESOC BW model of the GPS IIR and IIR-M spacecraft (Garcia-Serrano et al. 2016). 3. ECOM + UCL: Here, the only difference with the ECOM + BW strategy is that the box is replaced by the grid-based model. 4. ECOM-2: No a priori radiation force model, only the D4B1 extended ECOM (Arnold et al. 2015) along with the three constrained along-track parameters for consistency. Here, our analysis with the ECOM-2 model is not as comprehensive as it might be (as it was not in the scope of our original study plan). This will be addressed in our future work. The pseudorandom noise (PRN) code assigned to those satellites during the analysis period is indicated. Here, the radiation force modelling strategy that accounts for the effects of SRP, bus MLI TRR, ERP and AT is applied. Satellites in eclipse season during the analysis period are indicated with an asterisk (*) in the PRN column. Units: m In Table 5, we present statistics of the estimated ECOM parameters from methods 1-3. We do not present the statistics for ECOM-2 because comparison between models with different parameterizations cannot be made directly. In general, as the daily solutions are fully independent of each other, smaller absolute values for both the mean and the RMS indicate an improvement in the force modelling. Using the ECOM + UCL model, we see a reduction in the absolute value of both the mean and the RMS values of the D 0 , B 0 and all along-track parameters (except the mean of the A 0 parameter that is the same for both), when comparing with results using the ECOM + BW model. We see a reduction in the RMS of the Y 0 parameter, but the mean increases. The mean and RMS of the Bsin and Bcos parameters increase, which indicate the presence of systematic effects that the ECOM + UCL combination is not effectively dealing with.
In Table 6, we show the statistics of the orbit overlap differences. A smaller value for both the mean and the RMS indicates an improvement in the force modelling. The RMS values in all components are smallest with ECOM + UCL approach. However, the mean values for both the radial and along-track components are smallest with the ECOM-2 approach and the mean value for the cross-track component is smallest with the ECOM approach. In the 3-d orbit overlap values, we see a drop of 9% and 4% in the mean and RMS values, respectively, when we compare the ECOM + UCL results against ECOM + BW.
The performance of the ECOM, ECOM + BW and ECOM + UCL orbit modelling strategies was also assessed in a series of orbit prediction tests. In these tests, 3 days of independently estimated orbits were used to determine a best fitting orbit represented by position and velocity coordinates and the eight parameters of the ECOM method described above. This best fitting orbit was then propagated into the future for 14 days, after the end of the 3-day fit interval. These predicted orbits were compared against the estimated orbits, on the first day and the last day of the prediction interval. In these tests, ECOM + UCL is the reference model and the orbits estimated using ECOM + UCL are used as the basis of the 3-day orbit fit and as the ground truth. These tests are done over 2016. Thus, the first 1-day prediction interval considered is day 4 of 2016 and the first 14th-day prediction interval is day 17. The results from these tests are presented in Table 7.
Comparing RMS orbit prediction errors using ECOM + UCL against ECOM + BW, after 1-day, we see the errors increase by 0.21 cm in the radial direction but fall by 2.20 cm and 1.81 cm in the along-track and cross-track directions, respectively. For the 14th day predictions, we see a reduction in the RMS orbit prediction errors of 20.35 cm and 13.79 cm in the radial and cross-track directions, but an increase of 15.52 cm in the along-track direction. Overall, these results suggest ECOM + UCL is outperforming ECOM + BW, in the day 1 and day 14 orbit prediction tests, but there are limitations to this analysis that should be addressed in future work for improved confidence in our findings. For example, because we use it as our reference model, it is possible that ECOM + UCL is favoured in these tests. Also, systematic errors, such as those that depend on the elevation of the Sun above the orbital plane, do not show up in the yearly statistics. A more complete picture of the comparative performance of the models should be investigated through time series analysis.

Conclusions and discussion
Recent developments to our radiation force modelling strategy were analysed using a new model for the GPS Block IIR and Block IIR-M satellites. Advances to our approach include: an enhanced bus model computation scheme (based on the spiral points algorithm) that uses ray-tracing to determine the radiation flux-spacecraft interaction from 10,000 points distributed uniformly on a sphere surrounding the spacecraft; an improved method (from a numerical stability perspective) for producing grids spaced at 1 • ×1 • intervals in latitude and longitude in the spacecraft frame using a padding process to extend the spiral points data in all directions to reduce the impact of edge effects; a quality assurance process that uses results from an EPS-sweep computation with 3960 points for selecting an optimal set of grids. The models produced, and the proposed implementation method, were refined using a series of verification tests within the UCL-ODL. The impact of introducing UCL's grid-based model into a full POD process was investigated by analysing the statistics of estimated orbit model parameters, orbit overlaps and orbit prediction errors. Combined, the results provide a good indication that introducing high-fidelity analytical force modelling into the POD process can improve the quality of the estimated orbits and further refinements of the approach to address current limitations are worth pursuing.
One of the difficulties with the high-fidelity approach lies in acquiring the spacecraft data (geometry, surface material Bold values indicate the lowest value error statistics among the various modelling strategies being compared Units: cm properties, attitude history, mass and mass history) that is required to produce accurate models. It is hoped that the results in this paper adds evidence to the case for making this data available to the science and engineering community, where it is possible-especially the detailed geometry and surface material properties. Using an accurate spacecraft model, it is possible to compute the high-fidelity radiation force model. However, the model computation time remains a problem and limits the number of development and testing cycles that we are able to perform. A typical model computation involves a 5 × 5m 2 pixel array projected onto the spacecraft model at a 1 mm pixel spacing. In such a case, there are 2.5 × 10 7 rays per incoming radiation flux direction, 1 × 10 4 different directions in the spiral points computation scheme, and so this requires 2.5 × 10 11 rayspacecraft surface interaction calculations. As it stands, this process takes~3 days to compute (job run-time as opposed to CPU time) on the UCL high-performance computing facility, Legion@UCL (can take longer when the facility is under heavy load), followed by~1 day of analyst's time to work through the process of generating the grids. Therefore, it is worth exploring methods for reducing model production times. We are beginning to explore the use of a graphical processing unit (GPU) to exploit standard computer graphics techniques in the computation of the radiation fluxspacecraft surface interaction-a process that naturally lends itself to being parallelised. This idea is demonstrated in Grey et al. (2017) where an OpenCL implementation of a radiation source-satellite surface interaction model that includes accurate modelling of diffuse reflection and apparent size of illumination source is used to simulate the impact of photoelectron emission on spacecraft surface charging. Also, we are exploring the use of an algorithm that re-organises the UCL spacecraft model components into a k-dimensional tree data structure, to speed up the ray-tracing algorithm by greatly reducing the number of ray-surface interaction tests that need to be performed (Li et al. 2018).