# Practical Use of Starlink Downlink Tones for Positioning
## 4 Multi-Epoch Positioning
### Positioning Algorithm  
The expected accuracy of the satellites' positions computed using TLE-SPG4 can be as high as a few km, and their velocities can be as high as a few m/s. As we showed in (6), the position and velocity errors of the satellites convert to frequency errors that are different between satellites. Thus, in the positioning filter based on the Doppler shift model presented above, these frequency errors will propagate into the satellite clock drift term, \(f_{sv}^{K}\). We will investigate the possibility of absorbing the orbit-induced frequency error in the satellite clock drift state under the conditions of this study (i.e., a degraded signal level, inaccurate satellite orbits, unknown satellite clocks, and inaccurate measurement timestamps). Consider an extended Kalman filter (EKF) with an error state vector \(\delta x\) consisting of the user's position and velocity errors, the user clock bias and drift errors, and the frequency error of each tracked satellite:  
\[\delta x=\left[\delta r,c\delta t_{r},\;c\dot{\delta t}_{r},\lambda\delta f_{s} ^{1},\ldots,\lambda\delta f_{s}^{K}\right]^{T} \tag{22}\]  
In this state, each satellite brings a single additional unknown, instead of six unknowns if we had included the 3D position and velocity errors of each tracked satellite into the state. We can write the state dynamics as:  
\[ \begin{split}\dot{\delta r}=w_{v}\\ c(\dot{\delta t}_{r})=c\left(\dot{\delta t}_{r}\right)+w_{c(\delta t _{r})}\\ c\left(\ddot{\delta t}_{r}\right)=w_{c(\delta t_{r})}\\ \lambda\left(\dot{\delta f}^{K}\right)=w_{\lambda\delta f}\end{split}  \tag{23}\]  
where \(w_{v}\), \(w_{c(\delta t_{r})}\), \(w_{c(\delta t_{r})}\), and \(w_{\lambda\delta f}\) are centered white Gaussian processes of the user's position error, the user clock bias error, the user clock drift error and the satellite-dependent frequency error, respectively. The system defined in (23) can be written as \(\dot{\delta x}=F\delta x+W\), with \(F\) being the dynamic matrix and \(W\) being the system noise matrix with covariance \(Q=cov(W)\). The prediction covariance matrix is written \(P=\phi P\phi^{T}+QT_{s}\), where \(\phi=I+T_{s}F\) is the state transition matrix and \(T_{s}\) is the time step of the prediction. The state is updated by the measurement vector \(\widetilde{z}\) composed of range rates of \(K\) satellites (\(1,..,K\)) over \(m\) epochs (from \(i+1\) to \(i+m\)). Thus, by stacking data of \(m\) epochs, we have:  
\[\widetilde{z}=\left[\widetilde{\dot{\rho}}_{i+1}^{1},\ldots,\widetilde{\dot{ \rho}}_{i+1}^{\widetilde{K}},\;\ldots,\widetilde{\dot{\rho}}_{i+m}^{1}, \ldots,\widetilde{\dot{\rho}}_{i+m}^{\widetilde{K}}\right]^{T} \tag{24}\]  
We assume that \(\widetilde{\dot{\rho}}_{i}^{\widetilde{k}}\) has already been corrected for the tropospheric induced frequency error. We compute the measurement prediction according to:  
\[\hat{z}= \begin{bmatrix}-\lambda\left(\hat{f}_{dappleri+1}^{1}-\hat{f}_{s}^{ 1}\right)+c\hat{\dot{d}}_{r}\\ \cdot\\ -\lambda\left(\hat{f}_{dappleri+1}^{K}-\hat{f}_{s}^{K}\right)+c\hat{\dot{d}}_ {r}\\ \cdot\\ -\lambda\left(\hat{f}_{dappleri+m}^{1}-\hat{f}_{s}^{1}\right)+c\hat{\dot{d}}_ {r}\\ \cdot\\ -\lambda\left(\hat{f}_{dappleri+m}^{K}-\hat{f}_{s}^{K}\right)+c\hat{\dot{d}}_ {r}\end{bmatrix}  \tag{25}\]where \(\hat{f}_{doppler}^{k}\) for \(j=i+1\) to \(i+m\) is computed based on the predicted user's position \(\hat{r}\) at time \(\hat{t}_{j}=t_{j}+\delta t_{r}\), and on the \(k\)th-satellite position and velocity at time \(\hat{t}_{j}-\tau_{j}^{k}\), i.e., \(r^{k}\Big{(}\hat{t}_{j}-\tau_{j}^{k}\Big{)}\) and \(v^{k}\Big{(}\hat{t}_{j}-\tau_{j}^{k}\Big{)}\), where \(\tau_{j}^{k}\approx\left\|\hat{r}\big{(}t_{j}\big{)}-r^{k}\big{(}t_{j}\big{)} \right\|\)\(/c\) is the propagation time of the signal between satellite and user's positions. By denoting innovation \(\Delta z=\bar{z}-\hat{z}\), the observation matrix \(H\) is defined such that \(\Delta z=H\delta x+\eta\), where \(\eta\) is the measure noise vector with covariance matrix \(R\).