## 3 Two Methods for BDS-3 Orbit Determination
### BDS-3 Orbit Determination Based on NES  
The parameters related to the BDS-3 experimental satellites orbits could be improved indirectly through the MGEX observations, which accurately obtains the relevant parameters (station coordinates, tropospheres, and clock offsets). This study improved the data availability based on the NES of two networks in BDS-3 orbit determination. The normal equation stacking method has been well discussed [27; 34]. In this study, based on the orbit determination, some more details about NES were given as follows:  
The normal equations of BDS-3 orbit determination based on iGMAS observations can be assumed according to [27]:  
\[\left[ \begin{array}{cc}N_{11}&N_{12}\\ N_{21}&N_{22}\end{array} \right]\left[ \begin{array}{c}X\\ Y\end{array} \right]=\left[ \begin{array}{c}W_{1}\\ W_{2}\end{array} \right] \tag{21}\]  
where \(X\) consists of the parameters of troposphere, station coordinates, and stations clocks; \(Y\) denotes the parameters related to the BDS-3, such as orbit parameters, satellites clock offsets; and \(N\) is the coefficient matrix corresponding to the normal equation. Similarly, the normal equations of the combined GPS and BDS-2 orbit determination based on MGEX and iGMAS are set up as in [27]:  
\[\left[ \begin{array}{cc}N^{\prime}_{11}&N^{\prime}_{12}\\ N^{\prime}_{21}&N^{\prime}_{22}\end{array} \right]\left[ \begin{array}{c}X\\ Y^{\prime}\end{array} \right]=\left[ \begin{array}{c}W^{\prime}_{1}\\ W^{\prime}_{2}\end{array} \right] \tag{22}\]  
where \(X\) represents the same parameters, such as iGMAS station coordinates, troposphere, and stations clocks, as in Equation (21). \(Y^{\prime}\) is the station-related MGEX as mentioned above (MGEX station coordinates, troposphere, and stations clocks), and GPS and BDS-2 satellite orbit parameters based on iGMAS and MGEX networks. Thus, the estimation of \(X\) is obtained using Equations (21) and (22) [27]:  
\[\left(N_{11}-N_{12}N^{-1}_{22}N_{21}\right)X=W_{1}-N_{12}N^{-1}_{22}W_{2} \tag{23}\]\[[N_{11}^{\prime}-N_{12}^{\prime}\big{(}N_{22}^{\prime}\big{)}^{-1}N_{21}^{\prime}]X =W_{1}^{\prime}-N_{12}^{\prime}\big{(}N_{22}^{\prime}\big{)}^{-1}W_{2}^{\prime} \tag{24}\]  
Combining Equations (23) and (24), then, assuming \(Q=N_{11}+N_{11}^{\prime}\), \(M=N_{12}N_{22}^{-1}N_{21}+N_{12}^{\prime}\big{(}N_{22}^{\prime}\big{)}^{-1}N_{2 1}^{\prime}\), \(W=W_{1}^{\prime}+W_{1}\), and \(K=N_{12}N_{22}^{-1}W_{2}+N_{12}^{\prime}\big{(}N_{22}^{\prime}\big{)}^{-1}W_{2} ^{\prime}\):  
\[[Q-M]X=W-K \tag{25}\]  
In Equation (25), \(X\) is based on the inversion of [\(Q\)-\(M\)]. This study proposes a simple method to improve computational efficiency and avoid program overflow caused by large matrix dimensions, as follows:  
The symmetric matrix \(Q\) is decomposed by the Cholesky algorithm as \(Q=R^{T}R\). \(R\) is a low triangular matrix:  
\[Q-M=R^{T}(R^{-T}(Q-M)R^{-1})R=R^{T}[E-\overline{M}]R \tag{26}\]  
where \(\overline{M}=R^{-T}MR^{-1}\). The eigenvalue decomposition of \(\overline{M}\) is used to obtain:  
\[\overline{M}=V^{T}DV \tag{27}\]  
In Equation (27), \(V\) is the orthogonal matrix after decomposition, and \(D=diag( \begin{array}{cccc}d_{1}&d_{2}&\cdots&d_{n}\end{array} )\). Thus:  
\[Q-M=R^{T}V^{T}(E-D)VR=U^{T}(E-D)U \tag{28}\]  
where \(U=RV\). Set \(X\) as \(n\) dimensions. \(E\) is an identity matrix with \(n\times n\) elements. Then, in Equation (28):  
\[[Q-M]^{-1}=U^{-1}\left[ \begin{array}{ccccc}\frac{1}{1-d_{1}}&&&&\\ &\frac{1}{1-d_{2}}&&\\ &&\ddots&&\\ &&&\frac{1}{1-d_{n-1}}&\\ &&&&\frac{1}{1-d_{n}}\end{array} \right]U^{-T} \tag{29}\]  
In Equation (29), set:  
\[U^{-1}=\left[ \begin{array}{ccccc}\overline{u}_{11}&\overline{u}_{12}&\cdots &\overline{u}_{1n}\\ \overline{u}_{21}&\overline{u}_{22}&\cdots&\overline{u}_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ \overline{u}_{n1}&\overline{u}_{n2}&\cdots&\overline{u}_{nn}\end{array} \right] \tag{30}\]  
Thus, \(X_{i}\) can be further simplified, if:  
\[u_{m}=\left[