Wideband SATCOM Model: Evaluation of Numerical Accuracy and Efficiency

The spectral method is typically applied as a simple and efficient method to solve the parabolic wave equation in phase screen scintillation models. The critical factors that can greatly affect the spectral method accuracy is the uniformity and smoothness of the input function. This paper observes these effects on the accuracy of the finite difference and the spectral methods applied to a wideband SATCOM signal propagation model simulated in the ultra-high frequency (UHF) band. The finite difference method uses local pointwise approximations to calculate a derivative. The spectral method uses global trigonometric interpolants that achieve remarkable accuracy for continuously differentiable functions. The differences in accuracy are presented for a Gaussian lens and Kolmogorov phase screen. The results demonstrate loss of accuracy in each method when a phase screen is applied, despite the spectral method's computational efficiency over the finite difference method. These results provide meaningful insights when discretizing an interior domain and solving the parabolic wave equation to obtain amplitude and phase of a signal perturbation.


Ionosphere Scintillation Model
Satellites are often used in long range remote sensing applications to avoid hazardous regions, measure a process without disturbance, and to probe large volumes economically and efficiently [CCRS, 2019]. Longer wavelength microwave radiation can penetrate through cloud cover, haze, dust, and all but the heaviest rainfall [CCRS, 2019]. The longer wavelengths are not susceptible to atmospheric scattering which affects shorter optical wavelengths. It is this reason that the UHF-band is used for the model described in this paper. The particular process explored in this research is the post data collection of a passive remote sensing satellite. The data acquired from the receiver is processed and transmitted down to an earth based receiver or ground receiver station as illustrated in Figure 1.

Figure 1. SATCOM Data Link to Ground Receiver
Ionosphere scintillation is of major concern to a radio frequency (RF) signal, as it can cause the incident electric field (E-field) to delay in time and phase. The amplitude of the waveform also fluctuates due to constructive and deconstructive interference as the random striations of the Ionosphere interacts with the Efield.

Parabolic Wave Equation
The parabolic wave equation (free space form) is given by the following equation (1): where U represents the complex waveform in the frequency domain whose field vector points in the x direction and propagates in the z direction. This equation is derived from the scalar Helmholtz equation by assuming a slow varying envelope on the wavefront propagation given by the substitution of (2): z z E U x z e γ − = The slow-varying envelope assumptions relaxes the phase sampling requirements that would otherwise be significant in full-wave simulations of wavefront propagation. This is especially necessary for the SATCOM problem as shown in Figure 1, where the propagation distance is long relative to the wavelength of the E-field.

Gaussian Lens and Phase Screens
The Ionic content and scintillation characteristics can be modelled with a divergent Gaussian lens and phase screen, respectively. The divergent lens is representative of high ionic content which contributes a positive phase to the E-field. The phase screen discussed in this research utilizes Kolmogorov theory by taking a slab layer of ionospheric irregularities and compressing it down to a thin layer. Mathematically, the effect is represented as a phasor convolved into the phasor of the incoming waveform. The screen contains a number of grid points that sample the Ionosphere's irregularities ranging within a pre-defined inner and outer scale. The incoming EM wave is constructed to have a total spatial width that is dependent on spatial frequency. The specific realizations of the random phase are obtained by sampling a distribution of phase-shifts that have statistical properties determined by the statistics of the electron density irregularities [Knepp, 1982]. These statistics are specified by the spatial power spectral density of the irregularities. Relating the power spectral density of the phase and the electron density begins with the two-dimensional autocorrelation of phase and electron density fluctuations. This relationship is translated to the one dimensional phase power spectral density of spatial dependence by integrating the autocorrelation function over the spatial wave propagation distance.
The numerical phase screen PSD is presented in the following form: Variable parameters defined as: where L is the length of the spatial grid and x ∆ is the spatial resolution of the samples along the grid.
Other conditions exist for the behaviour at the boundaries of the phase screen. These effects are minimized by establishing continuity at the boundaries of the screen and the numerical algorithms that solve the parabolic wave equation.

Numerical Techniques
The numerical techniques used for solving the parabolic wave equation (PWE) is the finite difference method (FDM) and the spectral method (SM). Each method has its respective advantages and disadvantages. FDM uses localized interpolants to discretize the continuum. FDM approximations are pointwise and constrained to a uniform grid. SM uses global trigonometric interpolants that are very accurate when the input function is smooth (continuously differentiable) and uniform. The computational cost of implementing each algorithm is dependent on the input data size and the methodology for solving the PWE to acquire solutions. FDM will require matrix inversion, a computationally expensive operation. However, SM's accuracy can only be improved if more points are used to sample a progressively non-uniform input function. Thus, a trade off will occur where it becomes impractical to use SM over the localized approximants of FDM in terms of computational cost.

Parabolic Wave Equation Analytic Solution
The PWE analytic solution is acquired by using the classic separation of variables technique on the original partial differential equation (PDE). This effectively converts the PDE into two ordinary differential equations (ODEs). The first ODE in terms of z is first order and can be solved by direct integration. The second ODE with respect to x is in the Sturm-Liouville form. It has a general solution that can be solved by substituting the boundary conditions into the ODE and simplifying. The result is a series solution with a sum of coefficients. The ODE solutions in terms of X and Z are recombined to form the final solution. The solution is refined further as it is determined that the coefficients are equivalent to Fourier integration of the initial condition ϕ(x). The final equation is represented in (5).
where λ is the wavelength and γ is the propagation parameter.
L is the length of the wavefront grid.
It should be noted that the summation series is truncated at 1500 when programmed in a computer, as this value provides a solution that does not vary significantly. It demonstrates the impracticality of using the analytic solution as due to the duration of computation time required to solve the series for a solution at a particular distance in the interior of the domain. Instead, this analytic solution will be used to determine the accuracy of the numerical methods applied to solve the parabolic wave equation at the receiver distance downrange from the satellite. More information on separation of variable techniques applied to PDEs can be found in [Farlow, 1982].

Finite Difference Method
The finite difference method uses local pointwise approximations defined on a finite, uniform grid. The particular formulation of the PWE is implemented by using the Crank-Nicolson scheme. The Crank-Nicolson scheme is unconditionally stable, meaning that any aspect ratio (spatial steps of z and x) can be selected to solve the problem while maintaining an error bounds that does not grow. This is the nature of implicit numerical schemes. The Crank-Nicolson formulation of the PWE is given by (6): The Crank-Nicolson scheme must be constructed in a matrix formulation that can be implemented in a computer program.
The initial step is to isolate the forward propagating variables to the left-hand side of the equation: Combine like terms and simplify: Computer implementation of the equation in the form above will require a tri-diagonal matrix that maintains the spatial orientation of the "j" and "n" terms.
Let the tri-diagonal matrix on the left side of the equation be "A", and the right side tri-diagonal matrix be "B" and "I" is an identity matrix:

Spectral Method
The Spectral Method improves the convergence of a solution by using more grid points to approximate the derivative of a function. The spectral method is known to be accurate to ( ) ( ) log O n n , where n is the order of the grid size used to sample points [Novak, 2017]. The objective is to formulate a periodic domain using a trigonometric basis on an equi-spaced grid. As the number of samples increases, the error should decrease provided the solution is infinitely differentiable and smooth. For a function with p -derivatives the th v spectral derivative typically has accuracy ( ) The derivative is approximated by the sequence of Fourier modes in the following manner: where k is the order of the derivative, n is the Fourier mode, L is the length of the spatial domain, and  n a is the Fourier coefficient after taking the fast Fourier transform.
An appropriate finite or discrete representation of the solution must be selected by using an interpolating function between the values ( ) j U x at some suitable points (or nodes) j x : The Fast Fourier transform is taken to acquire the Fourier U coefficients, ˆn a : The second order spatial derivative is approximated with the square of the Fourier mode sequence and formulated into a square diagonal matrix where each dimensional length is equal to the length of the input vector: The spatial frequency to spatial position conversion integration transform is: The spatial position to spatial frequency conversion integration transform is: This forms an Ordinary Differential Equation in the spatial frequency domain: The solution frequency domain solution is inversed Fourier transformed back to the spatial position domain:

Wideband Modulated Signal
The wideband signal applies for a dispersive (frequency selective) communication channel. A transmit signal is modulated using a complex envelope on a carrier frequency. The parabolic wave equation provides the realizations of the signal's amplitude and phase as it propagates spatially in the frequency domain. The signal is Fourier transformed inversely into the time domain at the receiver. More information regarding this methodology can be found in [Knepp, 1982]. The following equation represents this process: The wideband signal is subject to a divergent Gaussian lens or a Kolmogorov phase screen by multiplication of a phasor. The Gaussian lens phasor is given below: The phase screen phasor involves establishing an inner and outer scale of the irregularity located within the Ionosphere. The phase variance ordinarily depends on the variance of the electron density for each respective altitude that it is collected. In this paper, the phase variance is pre-defined to examine a significant variance in the spectral density. A high variance in density can occur in the case of Ionosphere disturbances such as geomagnetic storms.
σ is the phase variance, o L is the outer scale of the plasma, i l is the inner scale of the plasma, x K is the spatial frequency of the phase screen grid points corresponding to the x direction, m is the slope of the power spectral density, n K is the Bessel function of order n ( ) ( ) ( ) In Table 1, the input parameters for the Gaussian lens and phase screen simulations in the wideband model are provided.

RESULTS
The following results show that the hypothesis is correct regarding the numerical accuracy and computational performance of the FDM and SM algorithms.

Amplitude and Phase of Received Field
In Figure 3, the amplitude plot of the received waveform shows the edge diffraction caused by the Gaussian lens. Note the fluctuations in the magnitude towards the exterior of the wavefront. The diffraction is also apparent in the phase.   The SM average error comparison to the analytic solution declines by twelve orders of magnitude while the FDM average error declines by three orders of magnitude when comparing the phase screen error response and the Gaussian lens error response. The phase screen causes issues for the SM algorithm at the left and right hand boundaries. The global trigonometric interpolants of the SM algorithm are not as effective in computing the electric field for non-uniform input functions observed with phase screens due to an absence of smoothness and periodicity. This is also reflected in the variance and standard deviation of the error as shown in Figures 5 and 6 for a distance of 70,000 meters.

Time Delay
The time delay occurs as a result of the phase delay in the Efield's propagation. Figure 7 shows the Global time delay in terms of chips. The color-bar represents the amplitude of the electric field as seen on the map. The wave front of the E-field spans the vertical axis. The secondary (delayed) wavefront that protrudes outward from the center of the initial received wavefront is caused by the edge diffraction from the Gaussian lens as noted in Figure 3. In Figure 8, the error between the analytic Global Time Delay plot and the corresponding numerical method is shown for the Gaussian lens problem. The solution from FDM is delayed by an average of approximately 0.0006 chips compared to the analytic solution. The SM algorithm does not have an average delay differential from the analytic solution. Generally, the magnitude in the overall error is smaller than in the FDM algorithm by two full orders of magnitude. The calculation of the time jitter appears to be an exact match between the SM and analytic solutions as the FDM algorithm differs by 0.0008. Potential floating point error could contribute to the inaccuracy that causes the FDM algorithm to have an artificially added delay compared to the analytic solution. It can be inferred that the floating point error accumulation grows with the propagation distance of E-field.

Figure 8. Global Error difference (Gaussian Lens Interaction)
A well-defined band occurs on the right side of Figure 8 as a result of the FDM algorithm solution mismatch with the analytic solution. The left side shows no such band as the SM algorithm generally agrees with the analytic solution and has minimal inaccuracies.

Computational Efficiency
The FDM algorithm is significantly outperformed by the SM algorithm as shown in The time was calculated in the parabolic wave equation solution over a propagation distance, followed by the time duration of calculating the wideband frequency responses.

CONCLUSION
The FDM algorithm is greatly outperformed in computational efficiency by the SM algorithm. The inversion of a matrix required in the FDM algorithm is a significant disadvantage. As expected, the accuracy declines in SM and FDM when the input function encounters a phase screen due to the non-uniformity of the function. Despite the decline in accuracy, SM still closely resembles the analytic algorithm in the received amplitude and phase solutions. The subtle differences in accuracy contributed to the differences observed in the calculated time delays, as FDM did not resemble the analytic solution as closely as SM. Parallel computing should be examined to alleviate the time issues in the FDM algorithm so that the wideband frequency responses can be calculated simultaneously. Spline interpolation may also help achieve faster computation times by reducing the samples along the wavefront required by the algorithms input function. Future work will investigate implementing Spectral Element Methods (SEM) to improve accuracy of the local interpolants along the E-field wavefront grid. SEM will provide a more robust approximation of non-uniform input functions encountered when modelling scintillation with phase screens.