1.

Introduction

Accurately simulating the atmospheric turbulence behavior is well recognized as very challenging. For a variety of purposes such as the design and development of adaptive optics systems, speckle imaging techniques, and atmospheric propagation studies, it is essential to simulate good atmospheric phase screen models. Methods based on Zernike polynomial expansions,¹ fast Fourier transform (FFT)-based methods,^2–8 and low-frequency optimization method⁹ have been in use for this purpose. The Zernike polynomial method, which is widely in use, has a limitation due to the maximum number of coefficients needed for accurate compensation. The optimization method, which compensates accurately for low-frequency part of the spectrum using unequal sampling and unequal weight in low-frequency region, does not cover high-frequency deficiencies. Among these, FFT-based methods are computer memory size friendly and widely accepted. But, FFT operators assume uniform sampling for the non-uniformly distributed phase power spectrum, which can lead to underestimation in the low- and high-frequency out of band regions, as shown in Fig. 1. Thus, it has limitations in recreating the true phase power spectrum. To compensate for low-frequency components, Johansson and Gavel³ suggested employing the modified subharmonics equation (an adaptation from Lane et al.¹⁰), which works well up to an infinite outer scale length. Sedmak⁶ later compared the performance of this method with that of Lane et al.¹⁰ by actually calculating the phase structure function from the simulated screen. He improved upon Lane et al.¹⁰ by employing different fine-tuned subharmonic weights for different $G/L_0$ ratios. Results from his analysis show that these FFT-based simulations are accurate for large screen size ($G$) to outer scale parameter ($L_0$) ratios. For a screen size of $G=200\mathrm{m}$ and outer scale of $L_0=25\mathrm{m}$, the maximum relative error in the simulation approaches 1%. Our simulations demonstrate that the errors from low-frequency components start shooting up once we move to smaller $G/L_0$ ratios, even after compensating with modified subharmonics.

Fig. 1

Comparison between simulation band and full band.

In Fig. 1, we show¹¹ a typical situation where the simulation band $(\frac{1}{G}-\frac{1}{\mathrm{\Delta }})$ is actually smaller than full band $(\frac{1}{L_0}-\frac{1}{l_0})$, where $\mathrm{\Delta }$ is the sampling size defined as the ratio of screen size $G$ to sampling number $N$ and $l_0$ is the inner scale parameter. In practice, the simulations are often curtailed at the low-frequency end, to a few times the optical beam size (say as determined by the telescope or laser beam diameter), whereas at the high frequency end, they often extend to only a few times that determined by the Fried parameter $r_0$. Clearly, the larger the simulation band to full band ratio, the more accurate the simulated results will be.

On the one hand, the apertures of upcoming and future astronomical telescopes are often of the same order or even larger than the typical median outer scale sizes of about 20 to 25 m.¹² On the other hand, wavefront sensing and compensation technologies are fast progressing that Nyquist sampling at $r_0$ scales even for large aperture telescopes are becoming quite possible. Thus, atmospheric turbulence simulations have to deal with a wide range in a multi-dimensional parameter space.

For working with very small apertures relative to the outer scale, it may appear that we need to simulate only a relatively small screen size. But cutting out small apertures from a larger screen introduces deviation from phase structure function due to misrepresentation of low-frequency components present in the small screen power spectrum.

In this paper, we present an approach and a corresponding algorithm to deal with phase screen simulations for a wide range of $G/L_0$ ratios, using the FFT-based method.¹³ Our technique builds upon the modified subharmonic approach of Johansson and Gavel³ and is inspired by Jingsong Xiang’s work.⁷ It works well for space- and time-invariant, zero intermittency atmospheric turbulence. Section 2 explains how to obtain phase autocorrealtion matrix using phase power spectrum, Sec. 3 presents the algorithm part to compensate for the remaining error in phase structure function calculation, Sec. 4 steps through the implementation of the algorithm with the help of a flow chart, Sec. 5 covers the validation of the technique using results from simulated phase screens, and Sec. 6 provides the concluding remarks.

2.

Obtaining Phase Autocorrelation Matrix Using Phase Power Spectrum

The 2D phase structure function and phase autocorrelation matrix are related as follows:¹⁴

The 2D phase autocorrelation matrix after compensating with subharmonics is represented as

3.

Algorithm to Compensate for Residual Error in Phase Structure Function

To calculate the remaining error in the final $B_{\varphi }(m,n)$, Eq. (4) is converted to phase structure matrix $D_{\varphi }(m,n)$ with the help of Eq. (1) with the assumption that $B_{\varphi }^{FFT}(\mathrm{0,0})$ and $B_{\varphi }^{SUB}(\mathrm{0,0})$ are zero because we are not concerned about the piston component. This gives the following equation

We need to compensate $D_{\varphi }$ so that $D_{\text{error}}$ is minimized. However, simply adding error correction terms in the $D_{\varphi }$ matrix directly would only introduce further error into the system while taking the Fourier transform. This is because any matrix or curve in general will have higher order moments. Thus, if we take the Fourier transform of the adjusted equation, the resultant curve will have completely different moments and hence power spectrum. This is because the transition between two steps in the error matrix will not be smooth, which introduces additional errors due to Gibb’s phenomena such as overshoots. Just curve fitting with any function does not satisfy the additional requirement of leaving the power spectrum unaffected by the process. What we really need is to introduce a smoothening operator such as a Gaussian function in the phase autocorrelation matrix, which exactly compensates for $D_{\text{error}}$.

For that we have developed an iterative algorithm (see the flow chart shown in Fig. 4) and implemented it in Matlab. The algorithm looks for the perfect Gaussian curve that minimizes the $D_{\text{error}}$ matrix. We use Matlab cftool to initially determine the correct 1D Gaussian matrix and later convert it into a 2D matrix by exploiting the fact that $B_{\varphi }(r)$, $B_{\text{theory}}(r)$, and $B_{\mathrm{SUB}}(r)$ all are dependent on $r$ only and hence are center symmetric functions. We call the fitted Gaussian phase structure matrix $D_{\text{gauss}}$ and the corresponding Gaussian phase autocorrelaiton matrix $B_{\text{gauss}}$ [using Eq. (1)].

The final equation for $B_{\text{tot}}$ can then be written as

Here, we have used $B_{\text{gauss}}=-D_{\text{gauss}}/2$ from Eq. (1). A look at the power spectrum of $B_{\mathrm{tot}}(m,n)$ in Fig. 2 shows that it contains negative terms⁷ for the case of $G/L_0<1$. Directly putting those frequency terms equal to zero leads to a loss in the energy spectrum. Hence, $B_{\mathrm{tot}}(m,n)$ matrix needs to be preprocessed to eliminate most of these negative values in the power spectrum. Over small frequencies, piston and tip/tilt components account for most of these high magnitude negative elements. Therefore, we first extract the piston and tip/tilt components from the phase autocorrelation matrix $B_{\text{tot}}$. The tip/tilt component from phase autocorrelation matrix is given as⁷

Fig. 2

Negative power spectrum values for small $G/L_0$ ratios.

After setting, $B_{\text{tilt}}(0)=0$, the remaining phase autocorrelation matrix is given as follows:

The power spectra $f_{\text{high}}^2$ and $f_{\text{tilt}}^2$ of the phase autocorrelation matrices $B_{\text{high}}(r)$ and $B_{\text{tilt}}(r)$ are obtained by standard Fourier transformation. Figure 3, shows the remaining negative power elements present in the power spectrum of $B_{\text{high}}$ matrix. In comparison to Fig. 2, the largest negative power contributions fall by factor of three order of magnitude. Now, we set the negative values in $f_{\text{high}}^2$ equal to zero by hand. The new error matrix is given as

Fig. 3

Residual negative power spectrum values after removing tip/tilt from $B_{\mathrm{tot}}$ for small $G/L_0$ ratios.

4.

Implementation of the Compensation Algorithm

In this section, we explain the error compensation algorithm with the help of the flow chart shown in Fig. 4. Brief explanations of each of the steps from $L_1$ to $L_{12}$ are given below.

Fig. 4

Flow chart for error compensation.

fmincon: High Frequency error Optimization

Table 1 shows the result of curve fitting using cftool for different cases of $G/L_0$ and N, which demonstrates that the Gaussian error matrix can compensate for a wide range of $G/L_0$ ratios and under different sampling constraints.

Table 1

Result of curve fitting against Gaussian function for different cases of G/L0 and N in terms of MRE for fixed r0=0.2  m.

4.1.

Example

To illustrate the robustness of the above algorithm, we have taken an example with $G=80\mathrm{m}$, say for a large future telescope, $N=256$, and median value of $L_0=20\mathrm{m}$.

The output from the above algorithm corresponding to minimum error entry as in (step $L_{12}$), has been plotted against $EF=5$ and $N_p=8$. Figure 5 shows a 3D rendering of $D_{\text{error}}$ matrix with a maximum separation of up to 40 m, corresponding to Eq. (4).

Fig. 5

3D $D_{\text{error}}$ matrix for case $G=80\mathrm{m}$, $L_0=20\mathrm{m}$, $r_0=0.2\mathrm{m}$, and $N_x=N_y=256$.

Figure 6 shows 1D $D_{\text{error}}$ matrix (radial section from 3D $D_{\text{error}}$ matrix) along with 1D fitted curve $D_{\text{gauss}}$ including the extrapolated part, for a maximum separation of up to 40 m.

Fig. 6

1D $D_{\text{error}}$ matrix fitted against $D_{\text{gauss}}$ matrix, along with an extrapolated part of the curve. Here $G=80\mathrm{m}$, $L_0=20\mathrm{m}$, $r_0=0.2\mathrm{m}$, and $N_x=N_y=256$.

Lastly, MRE values are stored against 500 entries corresponding to $N_p$ ranging from 1 to 10, $EF$ ranging from 1 to 10 and $SoG$ ranging from 2 to 6 after performing cftool fitting. This has been arranged in descending order and shown in Fig. 7, which illustrates a large set of iterations where errors are $<1\%$ and entry with minimum MRE has been picked up. Typical time required to perform each iteration for this case is $\approx 4.9\mathrm{s}$ on 2.3 GHz quad-core Intel Core i5 Macbook Pro 2018 model.

Fig. 7

Maximum relative error MRE with the maximum number of iterations for $G=80\mathrm{m}$, $L_0=20\mathrm{m}$, $r_0=0.2\mathrm{m}$, and $N_x=N_y=256$.

5.

Validation via Phase Structure Function Calculated from Simulated Phase Screen

To obtain the phase screen $\varphi (m,n)$ from the power spectrum, the following relation is used:⁷

For validation, we consider scenarios of apertures up to 40 m, i.e., $G=80\mathrm{m}$, at a median $L_0=20\mathrm{m}$ for two different sampling levels $N=256$ and 512. The phase structure function, defined as an ensemble average of differences of phases at various separation,¹⁴ has been averaged over 100 K independent frames. The relative error in phase structure function is calculated as follows:

Here, $D_{\varphi }^{\text{sim}}(r)$ is the phase structure function from the simulated phase screen. The magnitude of the peak relative error $\max (|err(r)|)$ is $<1.6\%$ for $N=256$ and $<0.5\%$ for $N=512$ as shown in Fig. 8.

Fig. 8

Left: Compares simulated structure function with respect to theoretical structure function for maximum separation of $G/2$ for two different cases $N=256$ and 512, for fixed $G=80\mathrm{m}$, $r_0=0.2\mathrm{m}$, and $L_0=20\mathrm{m}$. Right: Calculates the magnitude of relative error in simulated structure function for maximum separation of $G/2$, for both cases.

We also illustrate the performance (shown in Fig. 9) with parameters $G=1\mathrm{m}$, $L_0=100\mathrm{m}$, and 1000 m, $N=128$, $r_0=0.2\mathrm{m}$, which cover the extreme cases (very low $G/L_0$ ratios) which leads to the maximum error in the simulation. The magnitude of the peak relative error $\max (|\mathrm{err}(r)|)$ is $<1.6\%$ for $L_0=100\mathrm{m}$ and $<1.8\%$ for $L_0=1000\mathrm{m}$. Figure 10 shows one realization of the corresponding phase screen plots for $L_0=100\mathrm{m}$.

Fig. 9

(a) Compares simulated structure function with respect to theoretical structure function for maximum separation of $G/2$ for two different cases $L_0=100$ and 1000 m, for fixed $G=1\mathrm{m}$, $r_0=0.2\mathrm{m}$, and $N=128$. (b) Calculates the magnitude of relative error in simulated structure function for maximum separation of $G/2$, for both cases.

Fig. 10

Phase screen for case $G=1\mathrm{m}$, $L_0=100\mathrm{m}$, $r_0=0.2\mathrm{m}$, and $N_x=N_y=128$.

Figure 11 contains results of magnitude of the peak relative error in $D_{\varphi }^{\text{sim}}(r)$ for the case of different sampling points $N=128/256/512/1024$, for $L_0$ ranges up to 1024 m, $r_0=0.2\mathrm{m}$ and $G=2\mathrm{m}$. Similarly, Fig. 12 contains results of magnitude of the peak relative error in $D_{\varphi }^{\text{sim}}(r)$ for the case of different sampling points $N=128/256/512/1024$, for $G$ ranges up to 100 m, $r_0=0.2\mathrm{m}$, and $L_0=25\mathrm{m}$.

Fig. 11

The magnitude of the peak relative error for $N=128$, 256, 512, and 1024 for $L_0$ ranges up to 1024 m. Here $G=2\mathrm{m}$ and $r_0=0.2\mathrm{m}$.

Fig. 12

The magnitude of the peak relative error for $N=128$, 256, 512, and 1024 for screen size up to 100 m. Here, $L_0=25\mathrm{m}$ and $r_0=0.2\mathrm{m}$.

There are some outliers that have a high residual error as shown in Figs. 11 and 12, because we have not set the phase autocorrelation matrix to zero for $r>G/2$. The reason for this stems from the non-zero value of $B_{\text{high}}(r)$ at $r>G/2$, where $B_{\text{high}}(r)$ is formed from the removal of piston and tilt from $B_{\mathrm{tot}}(r)$. This can be resolved using a better smoothening operator, which we can multiply with $B_{\text{high}}$ so that it falls to zero progressively and not sharply. This can provide further improvement in the compensation.

6.

Conclusion

In this paper, we put forward a new method to compensate for the residual error in both the low and/or high-frequency region of FFT simulated phase screens that remain even after compensating with the modified subharmonic method. This method provides accurate phase screen structure for even $G/L_0$ ratios as small as 1/1000 plus screen sizes as large as 100 m. No patch normalization factor is needed, no need to calculate subharmonic weight coefficient¹⁰ and weights to compensate for high-frequency components, as done by Sedmak.⁶ While adequately large $G/L_0$ ratios may be the natural choice for modern large telescopes, simulations that deal with applications such as laser beam propagation through turbulent atmospheres would tend to have very small $G/L_0$ ratios. The method we propose is independent of the $G/L_0$ ratio choice. However, we emphasize that properly sampling $r_0$ and the high-frequency phase spectrum forces $N$ to be at least larger than $(2G/r_0)$ and preferably up to the inner scale limit $(2G/l_0)$. Currently, we have demonstrated this technique for only circular screens. We have used a GPU processor with total number of 128 cores, such that each iteration runs independently on each core. We have fixed the number of iterations to 500, although increasing this will lead to improvement of errors in some cases. Each core takes about $\sim 0.06$, $\sim 0.1$, $\sim 0.36$, and $\sim 1.1\min$ for sampling sizes of 128, 256, 512, and 1024. On the above GPU system, this translates to total computing times for error minimization of about $\sim 0.2$, $\sim 0.5$, $\sim 1.25$, and $\sim 4\min$ for sampling sizes of 128, 256, 512, and 1024, respectively. Once the coefficients are determined, generating multiple phase screen realizations from the corresponding power spectrum takes a few milliseconds at most on 2.3 GHz quad-core Intel Core i5 Macbook Pro 2018 model. Then it takes $<1\min$ to $\sim 10\min$ for averaging over 100k phase screens, for sampling sizes ranging from 128 to 1024.

The uniqueness of our approach is its ability to deal with any $G/L_0$ ratio within a very broad range, in an automated iterative process with little human intervention needed for tuning of parameters. Any standard FFT-based approach (say Sedmak’s⁶ compensated approach) for a given computer platform is computationally fast, only if we already have determined proper measures of the various compensating components such as the patch normalization factor, subharmonic weights, and high-frequency weights. Typically, determining these compensations is where the difficulty is due to mathematical complexity, algorithmic limitations, and/or computational power requirements. Our algorithm accomplishes the determination of the required compensation in very little time, with fairly reasonable computational power while at the same time keeping the residual errors competitively low using an appropriate compensator. Other existing FFT based approaches have limitations in their operable $G/L_0$ range. For example, Xiang et al.⁷ offer a very computationally fast approach but does not apply subharmonic compensation. Zhang et al.⁹ does not consider compensation for the high-frequency error, thus leaving a residual error of more than 100% in the high-frequency region. Sedmak’s⁶ approach needs the determination of accurate subharmonic weights for different $G/L_0$ ratios. The accuracy of our method from low-frequency to the high-frequency range is between 0.5% and 3% for $G/L_0$ as low as 1/1000 and screen size up to 100 m in diameter.