## 3 Simulations and Results
### Orbit Determination Simulation
#### 3.2.4 Effects of a Higher-Fidelity Dynamics Model  
To further evaluate the performance of the proposed OD method based on the reduced dynamics model, a higher-fidelity orbital dynamics model was additionally employed for comparison. The higher-fidelity model explicitly models the perturbation acceleration due to atmospheric drag (\(\mathbf{a}_{\mathrm{AD}}\)) as a distinct component, separate from the low-thrust acceleration, as defined in Equation (21). The NRLMSISE-00 atmospheric mass density model [34] is used to calculate \(\mathbf{a}_{\mathrm{AD}}\). The area-to-mass ratio and the drag coefficient \(C_{\mathrm{d}}\) are assumed to be 0.02 m\({}^{2}\)/kg and 2.2, respectively. The model and related parameters are used throughout both the OD process and the subsequent OP.  
\[\left\{\begin{array}{l}\frac{dt}{dt}=\mathbf{v}\\ \\ \frac{dt}{dt}=-\frac{\mu}{r}\mathbf{r}+\mathbf{a}_{\mathrm{NS}}+\mathbf{a}_{\mathrm{AD}}+M _{RTN\to ECI}\left[\begin{array}{c}0\\ a_{T}\\ 0\end{array}\right]\end{array}\right. \tag{21}\]  
Numerical simulations were conducted using the same observational data as Case 3 in Section 3.2.1, while maintaining all other parameters consistent with those specified in Section 3.2.2. This includes the initial covariance values for position, velocity, and along-track acceleration, as well as the measurement noise level. Figure 17a,b illustrate the iteration histories of position errors and \(a_{\mathrm{T}}\) estimates at the selected epoch, respectively. Compared with Figure 11, the overall iterative trend remains consistent. Cases E and F, which adopt larger initial covariance values for position and velocity, demonstrate significantly improved stability as well as smaller final position errors. The magnitude of the position error at the chosen epoch for these 2 cases is also similar to that observed in Figure 11, at approximately 130 m. The average estimated value of the along-track acceleration shown in Figure 17b is about \(7.66\times 10^{-5}\) m/s\({}^{2}\), and the value for cases shown in Figure 11b is \(7.04\times 10^{-5}\) m/s\({}^{2}\).  
For comparison, the 24 h OP RMS position errors, calculated based on the final OD results using different motion equations, are summarized in Table 4. The corresponding maximum position errors are presented in Table 5. The results show that for cases using the higher-fidelity model with Equation (21), similar to performance observed with Equation (4), although Cases E and F achieve lower positional/velocity errors at the selected epoch after OD, the 24 h RMS and maximum position errors of OP are actually larger. Despite this, for all Cases A-F, adopting the high-fidelity model of Equation (21) results in a notable improvement in 24 h OP accuracy. The predicted RMS position errors decrease by at least 250 m (12%) for Case E and up to 900 m (50%) for Case B. And the Maximum position errors reduce by at least 500 m (16%) for Case C and up to 1950 m (37%) for Case A.  
Figure 16: Evolution of “true” semi-major axis for Starlink satellite 58837.  
Naturally, as the complexity of the orbital dynamics model increases, the total computational time also rises, mainly due to the additional cost of calculating atmospheric drag acceleration in each iteration. Comparing the computation time for the cases shown in Figures 11 and 17, with identical initial state covariance, the OD process using Equation (21) takes approximately 3 to 5 times longer than that using Equation (4).