## 2 System Model
### Geometry Model Based on Geocentric Angle  
Traditionally, models for terrestrial heterogeneous systems (THS) and UAV-aided THS rely on omnidirectional antennas, where the received signal power is predominantly determined by transmit power and the Euclidean distance between user terminal (UT) and base stations (BS). However, satellite systems necessitate directional antennas due to significant propagation losses over vast transmission distances spanning thousands of kilometers. The performance of such directional antennas is intricately influenced by the angle away from the prime orientation, presenting a challenge in accurately modeling multi-layer satellite systems (MSS) due to the nuanced relationship between distance and off-axis angle. The model presented for the NGEO MSS, as illustrated in Figure 2, is rooted in stochastic geometry structures. Specifically, let \(\{(r_{k},\theta,\zeta):r_{k}=r_{E}+H_{k},0\leq\theta\leq\theta_{\max},0\leq \zeta\leq 2\pi\}\) denote the coordinates in \(\mathbb{R}^{3}\) spherical surface, wherein \(r_{E}\) denotes the Earth radius, \(\theta\) represents the angular offset from the axis of the line \(OG\) connecting the Earth core to the GEO satellite, \(\zeta\) denotes the rotation angle and \(H_{k}\) represents the \(k\)-th layer sphere's altitude. We account for the potential interference of NGEO UTs with a GEO UT within a finite zone and model the geometry range above the horizon of the GEO UT, as depicted in Figure 2, the maximum off-axis angle \(\theta_{\max}\) is given by  
\[\theta_{\max}=\arccos\frac{r_{E}}{r_{E}+H_{k}}, \tag{1}\]  
where \(\arccos(\cdot)\) represents the trigonometrical function. In light of the employment of directional antennas by the NGEO UTs, we presuppose no mutual interference between them. Furthermore, it is postulated that the operational range of NGEO UTs aligns with the identical maximum off-axis angle, designated as \(\theta_{\text{max}}\), as those of the NGEO satellites. The uniformly distributed Bernoulli point process (BPP) surpasses the Poisson point process (PPP) for approximating the global satellite system [19], where the number of satellites is finite and static. Nonetheless, the scenario examined in our manuscript involves a finite domain and the number of observable satellites fluctuates dynamically due to the Earth/satellite mobility. Thus, the approximation methodology for global satellite system is inadequate for our scenario. Subsequently, the positions of NGEO UTs and satellites encapsulated within a designated region are simulated utilizing three-tier concentric spherical PPPs, denoted by \(\Phi_{E}\) and \(\Phi_{k}\), respectively, with corresponding densities [20]  
\[\rho_{E}=\frac{N_{\text{U}}}{4\pi r_{E}^{2}}, \tag{2}\]  
\[\rho_{k}=\frac{N_{\text{S},k}}{4\pi r_{k}^{2}}, \tag{3}\]  
where \(N_{\text{U}}\) and \(N_{\text{S},k}\) are the number of UTs and \(k\)-th tier satellites, respectively. For convenience, let us place a typical GEO satellite and its corresponding UT positioned on the axis \(OG\) within Figure 2, treating it as the off-axis angle's origin. Using basic geometry, the angle \(\theta\) of a NGEO satellite antenna, related to the geocentric angle \(\hat{\theta}_{k}\), is expressed as  
\[\theta\big{(}\hat{\theta}_{k}\big{)}=\arcsin\left(\frac{r_{E}\sin\big{(}\hat{ \theta}_{k}\big{)}}{d\big{(}\hat{\theta}_{k}\big{)}}\right), \tag{4}\]  
where \(\hat{\theta}_{k}\) represents the geocentric angular between the UT and a satellite of \(k\)-th tier and the distance of that is  
\[d\big{(}\hat{\vartheta}_{k}\big{)}=\sqrt{\left(r_{E}+H_{k}\right)^{2}+r_{E}^{ 2}-2r_{E}(r_{E}+H_{k})\cos\big{(}\hat{\vartheta}_{k}\big{)}}. \tag{5}\]  
By utilizing (4) and (5), we can elucidate the gain \(G(\theta)\) of the antenna and the separation between satellite and UT using the geocentric angle \(\hat{\vartheta}_{k}\) for NGEO satellites.