## 2 Formulation
### Equivalent Ionospheric Current-Density Model  
In this paper, the influence of electron density and collision frequency in the ionosphere is introduced into the model in the form of current density. The expression of current density in the model is the same as that of electric-field intensity. The current density is based on the momentum equation for electrons moving through the ionosphere, as follows:  
\[\frac{d\vec{J}}{dt}+\Big{(}\frac{e}{m}\Big{)}\vec{J}\times\vec{B}_{0}+v\vec{J }=\epsilon_{0}\omega_{p}^{2}\vec{E} \tag{3}\]  
where \(m\) is the electron mass, \(e\) is the electron charge, \(v\) is the electron motion speed, \(\vec{B}_{0}\) is the geomagnetic field strength, and \(\omega_{p}\) is the plasma frequency. In the spherical coordinate system, the current-density equation is expressed in the form of three components \(J_{r}\), \(J_{\theta}\) and \(J_{\phi}\). Taking the component \(J_{r}\) as an example, the component of the current-density equation is expanded from the \(n-1/2\) moment to the \(n+1/2\) moment. The formula also contains the current density in the \(\theta\)- and \(\phi\)- directions, so the three components of the current density are interrelated. By performing a time average on the in \(\theta\)- and \(\phi\)- direction components of the current density, the components are made independent of each other. Similarly, for the differential equations of the other two components in the current-density equation, the same method is used to make the  
Figure 1: Electromagnetic-wave propagation model in the ionosphere. components independent of each other. The result of the expansion \(J_{r}\) of the component is shown as follows. \[\frac{\partial J_{r}}{\partial t}+vJ_{r}+\frac{e}{2m}\left[B_{\varphi}\Big{(}I_{ \theta}^{n-l/2}+I_{\theta}^{n+l/2}\Big{)}-B_{\theta}\Big{(}I_{\varphi}^{n-l/2} +I_{\varphi}^{n+l/2}\Big{)}\right]=\varepsilon_{0}\omega_{p}^{2}E_{r}^{n} \tag{4}\]  
The formula contains a derivative with respect to time, and the discretization error is relatively small only when the time step is small. During transmission in a lossy medium, the energy decay of electromagnetic waves is relatively large, and an exponential difference method is adopted in simulation to improve accuracy. To obtain higher precision, an exponential difference is applied to the above formula, and after rearrangement, an iterative equation for \(J_{r}\) can be obtained. \[J_{r}^{n+1/2}=2a/bJ_{r}^{n-1/2}+\omega_{b\theta}\left(J_{\varphi}^{n-1/2}+J_{ \varphi}^{n+1/2}\right)+2\varepsilon_{0}\omega_{p}^{2}E_{r}^{n}-\omega_{b \varphi}\Big{(}I_{\theta}^{n-1/2}+I_{\theta}^{n+1/2}\Big{)} \tag{5}\]  
The coefficients in the equation are  
\[a=e^{-v\Delta t},\quad b=\frac{1-e^{-v\Delta t}}{\
u},\quad\omega_{br}=\frac {e}{m}B_{0r},\quad\omega_{b\theta}=\frac{e}{m}B_{0\theta},\quad\omega_{b \varphi}=\frac{e}{m}B_{0\varphi} \tag{6}\]  
For the other two components, using the same method, similar equations can be obtained. Combining the three components together, we obtain the matrix form iterative equation for \(J\), which is shown below. It can be abbreviated as  
\[A\overline{j}^{n+1/2}=B\overline{j}^{n-1/2}+b\varepsilon_{0}\omega_{p}^{2} \overline{E}^{n} \tag{7}\]  
The coefficients in the equation are  
\[A= \begin{bmatrix}1&b\omega_{b\varphi}/2&-b\omega_{b\theta}/2\\ -b\omega_{b\varphi}/2&1&b\omega_{br}/2\\ b\omega_{b\theta}/2&-b\omega_{br}/2&1\end{bmatrix} \quad B= \begin{bmatrix}a&-b \omega_{b\varphi}/2&b\omega_{b\theta}/2\\ b\omega_{b\theta}/2&a&-b\omega_{br}/2\\ -b\omega_{b\theta}/2&b\omega_{br}/2&a\end{bmatrix}  \tag{8}\]  
During the iteration process, the presence of the matrix on the left side of the equation consumes significant time and space. This is because the coefficients are in matrix form. When solving this equation using MATLAB 2019b,if the coefficients contain matrices, the simulation process will consume a lot of time. Therefore, by further transforming the above equation, we can obtain a form that is more efficient in terms of both time and space.