## 3 LEO-LEO Occultation Simulation and Assessment Methods
### LEO-LEO Occultation Geometry Simulation Model
The signal propagation from Tx to Rx was approximated by a straight line, which is reasonable and adequate since the atmospheric bending is very small and can be neglected in such occultation event number and distribution analyses [41]. The atmosphere altitude range was set from 0 to 120 km, i.e., from the ground surface to the top of the neutral atmosphere. The location, i.e., the latitude, longitude, and altitude of the LEO-LEO occultation tangent point, was calculated according to [54]:  
\[r_{TR}=r_{Tx}-r_{Rx} \tag{3}\]  
\[\hat{n}_{TR}=(r_{TR})/\|r_{TR}\| \tag{4}\]  
\[R_{T}=\|r_{Tx}\cdot\hat{n}_{TR}\| \tag{5}\]  
\[r_{TP}=r_{Tx}-\hat{n}_{TR}\cdot R_{T}=(X_{TP},Y_{TP},Z_{TP}) \tag{6}\]  
\[\left\{ \begin{array}{l}\varphi_{TP}=\arctan\Big{(}Z_{TP}(N+h_{TP})/\sqrt{X _{TP}^{2}+Y_{TP}^{2}}/\big{(}N\cdot\big{(}1-e_{E}{}^{2}\big{)}+h_{TP}\big{)} \Big{)}\\ \\ \lambda_{TP}=\arctan(X_{TP}/Y_{TP})\\ \\ h_{TP}=Z_{TP}/\sin\varphi_{TP}-N\cdot\big{(}1-e_{E}{}^{2}\big{)}\\ \\ N=a_{E}/\sqrt{1-e_{E}{}^{2}\sin^{2}\varphi_{TP}}\\ \\ \varphi_{TP0}=\arctan\Big{(}Z_{TP}/\sqrt{X_{TP}^{2}+Y_{TP}^{2}}\Big{)}\\ \\ |\varphi_{TPi+1}-\varphi_{TPi}|\leq 10^{-12}\end{array} \right. \tag{7}\]where \(r_{TR}\) is the vector from Rx to Tx, \(r_{Tx}\) and \(r_{Rx}\) are the position vectors of Tx and Rx respectively, \(\hat{n}_{RT}\) is the unit vector along the direction from Rx to Tx, \(R_{T}\) is the projection of the Tx position vector on the Rx to Tx connecting line, \(r_{TP}\) is the position vector of the occultation tangent point, \((X_{TP},Y_{TP},Z_{TP})\) are the three-dimensional coordinates of the occultation tangent point in the Earth-centered Earth-fixed (ECEF) coordinate system, \(a_{E}\) is the equatorial radius of the Earth, \(e_{E}\) is the first numerical eccentricity of the ellipsoid, \(N\) is the distance from the surface back to the \(Z\)-axis intersection along the ellipsoid normal, and \(\varphi_{TP}\), \(\lambda_{TP}\), and \(h_{TP}\) are the latitude, longitude, and height of the occultation tangent point in the geodetic coordinate system, respectively, where \(\varphi_{TP}\) needs to be solved by iteration. The azimuth angle of the Tx relative to the Rx running direction ranges were set as 0-40\({}^{\circ}\), 320-360\({}^{\circ}\) for rising occultation events and 140-220\({}^{\circ}\) for setting occultation events, which seems reasonable since, when the Tx satellite is located in these azimuth ranges, the satellite's attitudes will be adjusted to make sure that the Tx and Rx beams point to each other. The elevation angle \(E\) and the azimuth angle \(A\) of the Tx relative to the Rx running direction are defined in a Cartesian coordinate system with the Rx mass center as the origin,  
\[\left\{ \begin{array}{l}x=v_{Rx}\\ y=v_{Rx}\times r_{Rx}\\ z=x\times y\end{array} \right. \tag{8}\]  
\[\sin E=z\cdot r_{TR}/(\|z\|\cdot\|r_{TR}\|) \tag{9}\]  
\[r_{A}=r_{TR}-\|r_{RT}\|\cdot\sin E\cdot z/\|z\| \tag{10}\]  
\[\cos A=r_{A}\cdot x/(\|r_{A}\|\cdot\|x\|) \tag{11}\]  
where \(v_{Rx}\) is the Rx velocity vector and \(r_{A}\) is the projection of \(r_{TR}\) in the xy coordinate plane.