Here is a simulation to predict the reentry date.
The simulation includes the Newtonian and the relativistic accelerations of all the planets, Sun and Moon.
The Earth's gravity field is modeled with the SGG-UGM-1 gravity model (computed using EGM2008 derived gravity anomaly and GOCE observation data) truncated to the degree and order 15 (to save the running time, while retaining good accuracy when compared to the full model).
For the calculation of the air density, I use the NRLMSISE-00 model along with an updated data file for the solar and geomagnetic indices. The actual indices can be found here: www.celestrak.com/spacedata/SW-All.txt, while the model used for the long term prediction of the indices can be found here: https://www.nasa.gov/msfcsolar.
The first step involves determining the best ballistic coefficient to minimize a particular simulation parameter. After 46 minutes, the program finds a ballistic coefficient of about 118 kg/m^2 (it’s not fixed, because the drag coefficient varies with the air composition).
The following graph shows the result for the two last months:
we see that the integrated mean radius vector fits very well the mean radius vector obtained with the CSpOC’s SGP4 library for the TLE epoch.
Now the simulation can be started:
1) with one TLE and the SGP4 propagator calculate the initial state (position and velocity) of the satellite for the TLE epoch;
2) propagate that initial state with a specially crafted propagator (my propagator is based on the 8(5,3) Dormand-Prince integrator);
3) when the satellite altitude drops below 70 km, stop the simulation; this is the reentry date.
Here’s the result obtained with 7 TLEs from 20003.47248329 to 20010.45318683:
The graph shows:
1-orbit mean radius vector (blue plot): radius vector averaged over the eccentric anomaly (it’s the semi-major axis). Not to be confused with the osculating semi-major axis.
1-orbit minimum radius vector (red plot): the smallest radius vector in one orbit. Not to be confused with the osculating perigee.
1-orbit maximum radius vector (green plot): the biggest radius vector in one orbit. Not to be confused with the osculating apogee.
actual eccentricity = (Ra - Rp) / (Ra + Rp), where Ra is the 1-orbit maximum radius vector and Rp is the 1-orbit minimum radius vector.
The average reentry date is 2025-05-20 and the difference between the later and the earlier date is just 3.5 days (which means that the TLEs are very accurate for this satellite). The only unknown is the air density.