The following is only a partial answer, but I will update as we make progress through this.
I believe the on-going discussion in the comments at For any state vector of a satellite at an arbitrary time T1, is there another state vector at a given time T2 that will result in the same orbit? is becoming too complex to keep it just at comments, so thought it might be better to go on from here (also, such discussion is probably more relevant in the original question).
As discussed in said comments, a possible way to tackle this question is to propagate the same initial state vector with SGP4, SDP4 and a high-precision numerical propagator (HPOP). The latter will be considered as the "true orbit", and the RMSD along propagation during 1 period will be calculated for both SGP4 and SDP4. Whichever comes closer to it is considered to perform better. By repeating this procedure across different values for the initial parameters, we can identify the region of parameter space where SGP4 outperforms SDP4 and vice versa.
A key aspect seems to obtain the starting point for HPOP. I thought the best way would be to propagate the provided mean orbital elements at time 0, which would result in Cartesian coordinates in TEME frame for position and velocity. These would then be converted to GCRF frame (for simplicity, let's assume always the same epoch, let's say 12 pm the 1st of January 2021, UTC time), which could be directly input into the HPOP.
I started doing so for two cases: eccentricity 0, and eccentricity 0.7, in both cases with a semi-major axis of 200 000 km, an inclination of 60 º and (mean) argument of perigee, longitude of the ascending node and mean anomaly of 0. Drag coefficient was also set to 0, although it should not make a difference at such high orbit.
I was surprised to see that the coordinates output by SGP4 and SDP4 at epoch were considerable different! This poses a problem as to how to choose the starting point for the HPOP. I have done a more systematic assessment by calculating the difference of the TEME coordinates output at epoch by SGP4 and SDP4 at different eccentricities and semi-major axes (all other parameters kept as described above). The following quick plot shows the results:
As you can see, the differences become considerable at high orbits, particularly at high eccentricities, reaching even several tens of thousands of kilometers!
Any ideas about how to choose a good starting point for the HPOP to perform additional comparisons would be greatly appreciated.
Edit 1
As pointed out by @uhoh in the comments to the answers to the question, the best approach is probably to do the following for each point in the space of parameters to be analyzed:
- Propagate the corresponding state vector with the high-precision numerical propagator (probably for 1 period)
- Generate TLEs that, when propagated with SGP4 and SDP4, lead to a trajectory as close as possible to the trajectory obtained with the numerical propagator. Note 2 different TLEs would be generated, one that leads to the closest trajectory with SGP4, and another with SDP4
- Having found the TLEs that lead to a trajectory as close as possible to that obtained with the numerical propagator, now we have starting points for each of the 3 propagators. These starting points would represent the same physical starting point. We can then propagate the 3 starting points, now for a time considerably longer than 1 period, and see, for each region of the parameter space, which one of either SGP4 or SDP4 deviate faster from the trajectory calculated with the high-precision orbital propagator.
I think the generation of optimal initial TLEs for SGP4 and SDP4 might be treated as an optimization problem, with the parameters to be optimized being the mean orbital elements that define the TLE. The target function to minimize would be a measurement of the deviation of the trajectory propagated with SGP4/SDP4 from the trajectory propagated with the HPOP. For example, possibly root-mean-square Euclidean distance (RMSD). Initial values can be easily obtained, for example by converting the GCRF position at epoch to osculating orbital elements and using these as starting values.
I have started to work in such direction, but have found something interesting. I did the following for an initial test:
- Generate a TLE with elements corresponding to a semi-major axis of 200,000 km, eccentricity 0.7, inclination 60º, and all other elements 0.
- Propagate it for a full period, approximately 14820 minutes, every 30 minutes
- Convert all output positions and velocities from TEME to GCRF.
- Take the position and velocity at epoch in GCRF and use it to propagate the trajectory for the same length of time and at the same time points with a high-precision numerical propagator (treating it as if it were a GPS block III satellite, using the physical parameters described here)
- We now have 2 series of GCRF coordinates at the same time points, one propagated with SGP4 and another one with HPOP. So I calculated the distance between each pair of points, and plotted it against time. Important notice: I am using my own (open-source) implementation of an HPOP here, so we might need to proceed with caution. Even though as far as I have tested the results are quite decent (achieving way higher accuracy than SGP4/SDP4 for GPS satellites and Planet Labs satellites after days of propagation), it probably still requires more testing. The result was the following:
Keep in mind this plot shows the divergence between SGP4 and HPOP (the latter of which we consider the "true trajectory") along 1 period. It can be easily seen that SGP4 starts to diverge massively from HPOP well before 1 orbital period. This creates a problem, since using RMSD between propagations along 1 orbital period as the target function to minimize to obtain an optimal TLE now becomes dominated by the very large distance between points after great divergence starts, therefore beating the original point of generating a TLE for SGP4 and SDP4 that leads to a trajectory as close as possible as that generated by HPOP at the early stages of propagation.
I believe a modification to the approach should be made. Either a shorter time should be used to calculate the RMSD that will be used as the minimization target to generate initial TLEs for both SGP4 and SDP4 (maybe, for example, up to 5000 minutes in this particular case?), or a different metric of difference between trajectories should be used. One that down weights the large differences seen towards the end of the orbital period.
Edit 2
I have performed the same comparison for a GEO satellite, with a period of around 1436 minutes, and the divergence is much lower along the full period. In this case, I have performed the same procedure using both SGP4 and SDP4. The comparison of each SGP4/SDP4 is done with a trajectory obtained with the HPOP from the ephemeris obtained at time 0 with each of them (i.e., not the same HPOP trajectory for both, although at this semimajor axis it does not make much difference). I am attaching the plot here for comparison:
As usual, any advice on what would be best would be very appreciated!