Simulation of anisoplanatic imaging through optical turbulence using numerical wave propagation with new validation analysis

Abstract. We present a numerical wave propagation method for simulating imaging of an extended scene under anisoplanatic conditions. While isoplanatic simulation is relatively common, few tools are specifically designed for simulating the imaging of extended scenes under anisoplanatic conditions. We provide a complete description of the proposed simulation tool, including the wave propagation method used. Our approach computes an array of point spread functions (PSFs) for a two-dimensional grid on the object plane. The PSFs are then used in a spatially varying weighted sum operation, with an ideal image, to produce a simulated image with realistic optical turbulence degradation. The degradation includes spatially varying warping and blurring. To produce the PSF array, we generate a series of extended phase screens. Simulated point sources are numerically propagated from an array of positions on the object plane, through the phase screens, and ultimately to the focal plane of the simulated camera. Note that the optical path for each PSF will be different, and thus, pass through a different portion of the extended phase screens. These different paths give rise to a spatially varying PSF to produce anisoplanatic effects. We use a method for defining the individual phase screen statistics that we have not seen used in previous anisoplanatic simulations. We also present a validation analysis. In particular, we compare simulated outputs with the theoretical anisoplanatic tilt correlation and a derived differential tilt variance statistic. This is in addition to comparing the long- and short-exposure PSFs and isoplanatic angle. We believe this analysis represents the most thorough validation of an anisoplanatic simulation to date. The current work is also unique that we simulate and validate both constant and varying Cn2(z) profiles. Furthermore, we simulate sequences with both temporally independent and temporally correlated turbulence effects. Temporal correlation is introduced by generating even larger extended phase screens and translating this block of screens in front of the propagation area. Our validation analysis shows an excellent match between the simulation statistics and the theoretical predictions. Thus, we think this tool can be used effectively to study optical anisoplanatic turbulence and to aid in the development of image restoration methods.

Simulation of anisoplanatic imaging through optical turbulence using numerical wave propagation with new validation analysis 1 Introduction Long range imaging in the atmosphere is impacted, and often limited, by optical turbulence. Random variations in the index of refraction along the optical path are caused by temperature variations and convection. 1 This leads to warping and blurring degradations. When the warping and blurring are spatially varying over the field of view of an imaging sensor, this is referred to as anisoplanatic imaging. Anisoplanatic conditions typically prevail when imaging extended scenes. With very small fields of view, the optical turbulence may be reasonably modeled as isoplanatic, using a spatially invariant point spread function (PSF). This is often the case with astronomical imaging. 1 The simulation of imaging under isoplanatic atmospheric turbulence has been well studied. 1 Commercial and free open-source software 2,3 is available for this purpose. However, anisoplanatic imaging simulations using numerical wave propagation have only recently been described in the literature. Some of the first such papers include those of Carrano 4 and Praus et al. 5 Further wave propagation-based anisoplanatic simulation works include that of Bos and Roggemann, 6,7 and Monnier et al. 8,9 Anisoplanatic simulations that do not involve wave propagation have also been developed. [10][11][12] The ability to accurately simulate the degradation effects of optical turbulence is important for several reasons. First, such simulations allow us to study the impact of atmospheric turbulence under a wide variety of imaging scenarios. Second, we are able to use simulated data to quantitatively evaluate turbulence mitigation methods. 4,[13][14][15] Quantitative analysis is possible because we have an objective truth image from which the degraded images are generated. Without such a truth image, assessment of restoration algorithms is subjective, and optimization cannot be automated.
*Address all correspondence to: Russell C. Hardie, E-mail: rhardie@udayton .edu We believe that quantitative performance analysis is critical to the advancement of image restoration methods.
The simulation methodology presented here is based on that of Bos and Roggemann. 7 We expand on, and extend, that work in several ways. We also provide a complete description of the simulation, including the details of the wave propagation method used. Like Bos and Roggemann, 7 our approach computes an array of PSFs for a two-dimensional (2-D) grid on the object plane. The PSFs are then used in a spatially varying weighted sum operation with an ideal image to produce a simulated image with realistic optical turbulence degradation. The degradation includes spatially varying warping and spatially varying blurring.
To produce the PSF array, we generate a series of extended phase screen realizations. Simulated point sources are numerically propagated from different positions on the object plane, through the phase screens, and ultimately to the focal plane of the simulated camera. Note that the optical path for each PSF will be different, and thus, pass through a different portion of the extended phase screens. These different paths give rise to spatially varying, but spatially correlated PSFs. As we shall show, these PSFs may be used to generate accurate anisoplanatic effects. The method we use to define the individual phase screen statistics is distinct from that of Bos and Roggemann. 7 Our approach is based on the constrained least squares optimization presented by Schmidt, 2 but is extended to include isoplanatic angle in the cost function. Another unique feature of our method is that we exclude the phase screen at the pupil plane. This aides in generating the appropriate level of anisoplanatic PSF correlations.
We also present a validation analysis here. In particular, we compare the simulated outputs with the theoretical anisoplanatic tilt correlation, 16 and a derived differential tilt variance (DTV) statistic. This is in addition to comparing the long-and average short-exposure PSFs and isoplanatic angle. We believe this analysis represents the most thorough validation of an anisoplanatic simulation to date. The current work is also unique that we simulate and validate both constant and varying C 2 n ðzÞ profiles. Furthermore, we simulate sequences with both temporally independent and temporally correlated turbulence effects. Temporal correlation is introduced by generating larger extended phase screens and translating this block of screens in front of the propagation area. This approach is similar to that described by Dios et al. 17 Our validation analysis shows an excellent match between the simulation statistics and the theoretical predictions. Thus, we believe this approach can be used effectively to study optical anisoplanatic turbulence and to aid in the development of image restoration methods.
The rest of this paper is organized as follows. Key anisoplanatic turbulence statistics are presented and discussed in Sec. 2. These statistics are used in setting up the simulation and in the validation analysis. The details of the simulation are presented in Sec. 3. The experimental results are presented in Sec. 4. Finally, we offer conclusions in Sec. 5.

Optical Turbulence Statistics
Variations in the index of refraction in the atmosphere can be modeled with a refractive index structure function. Using a Kolmogorov model, 1 this is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 7 5 2 D n ðrÞ ¼ h½nðx þ rÞ − nðxÞ 2 i ¼ C 2 n r 2∕3 ; where x and r are three-dimensional spatial coordinate vectors r ¼ jrj, nð·Þ is the index of refraction, and h·i represents an ensemble mean operator. The parameter C 2 n is the refractive index structure parameter. 1 The units of C 2 n are m −2∕3 and typical values tend to range from 1 × 10 −13 to 1 × 10 −17 . 2 When this quantity varies along the optical path, it may be expressed as a function C 2 n ðzÞ, where z is the distance from the source.
The atmospheric coherence diameter (or Fried parameter) 2 can be expressed as a weighted integral of C 2 n ðzÞ yielding E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 6 2 0 where k ¼ 2π∕λ is the wavenumber and λ is the wavelength. Note that this is for spherical wave propagation where z ¼ 0 is at the source, and z ¼ L is at the camera. The Fried parameter is directly linked to the atmospheric PSF and the tilt variance for a point source. A small r 0 , relative to the camera aperture, indicates a high level of atmospheric blurring and a high tilt variance. Another important statistic is the isoplanatic angle. 1 Two point sources that are separated by less than the isoplanatic angle will have a mean wave function phase difference at the aperture of <1 radian. 2,18 Another way to think of this is that points separated by less than this angle will have approximately the same PSF. The isoplanatic angle can also be expressed as a weighted integral of C 2 n ðzÞ, yielding 2,18 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 4 2 1 A statistic that relates to fluctuations in the wave function amplitude is the log-amplitude variance, 2 given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 3 5 4 It is interesting to visualize the three C 2 n ðzÞ weighting functions corresponding to r 0 , θ 0 , and σ 2 χ . This is done in Fig. 1. Note that the isoplanatic angle is most impacted by turbulence at the source, while r 0 is most impacted by turbulence near the camera. The log-amplitude variance is impacted most by the center of the optical path. Since the weighting for these three statistics covers the optical path in this balanced manner, we use them to determine the phase screen Fried parameters in our simulation, as shown in Sec. 3.5.

Optical Transfer Functions
The overall optical transfer function (OTF) for an imaging system in optical turbulence may be modeled to include the atmospheric OTF and the diffraction OTF. This is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 1 2 5 HðρÞ ¼ H atm ðρÞH dif ðρÞ; where ρ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and u and v are the spatial frequencies in units of cycles per unit distance. The atmospheric OTF is given by 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 6 3 ; 4 4 0 H atm ðρÞ ¼ exp −3.44 λlρ where l is the camera focal length and D is the aperture diameter. The parameter α here serves as a switch. In particular, the long exposure atmospheric OTF is generated by Eq. (6) when α ¼ 0, and the average short exposure OTF is given by Eq. (6) when α ¼ 1. The diffraction limited OTF for a circular aperture is given by 19 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 3 3 8 H dif ðρÞ ¼ where ρ c ¼ 1∕ðλf n Þ is the optical cut-off frequency, and the f-number is f n ¼ l∕D. Note that the theoretical atmospheric PSF is found by applying the inverse Fourier transform to Eq. (5).

Anisoplanatic Tilt Statistics
While the isoplanatic angle provides information about the level of anisoplanatism on a small scale, it does not provide insight into anisoplanatic behavior from points in the object plane with a large separation angle. Below, we describe two statistics that do capture large-scale anisoplanatic behavior. These are the two-axis Z-tilt correlation and the DTV. We shall use these as key validation metrics for the simulation.
To begin, let us define a two-axis Z-tilt vector as αðθÞ ¼ ½α x ðθÞ; α y ðθÞ T for a source viewed from the direction angle vector θ ¼ ½θ x ; θ y T . For a spherical wave characterized by the Kolmogorov power spectrum, an analytical expression for the Z-tilt correlation has been derived by Basu (now Bose-Pillai) et al., 16 following techniques outlined by Fried 20 and Winick. 21 The tilt correlation can be expressed as a weighted integral of C 2 n ðzÞ as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 3 2 6 ; 7 1 9 The path weighting function is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 6 6 8 where Δθ ¼ jθ 1 − θ 2 j. Figure 2 shows a plot of f c ðzÞ for different source separations, expressed in terms of pixel spacings. Note that the term "Nyquist pixels" here refers to pixel spacings corresponding to spatial sampling at the Nyquist rate, relative to the diffraction-limited optical cut-off frequency. The aperture size, path length, and Nyquist pixel spacing used in the evaluation are the same as those used in the simulations and are listed in Table 1. Note that the weight is maximum at the camera (z ¼ 7000 m) and drops down to zero at the source (z ¼ 0 m). This implies that the turbulence near the source does not contribute to tilt correlation seen at the camera. It is also apparent from Fig. 2 that the tilt correlations decrease with increasing angular separation between the sources.
Let us now consider the DTV. This statistic is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 3 7 8  n ðzÞ corresponding to three of the key optical turbulence parameters for L ¼ 7 km. Isoplanatic angle is heavily weighted toward the source end of path. The Fried parameter is weighted toward the camera end of the path, and the log amplitude variance is weighted primarily in the middle.
We can derive an expression for the DTV using Eqs. (8) and (9). To do this, note that the tilt variance term in Eq. (10) hα 2 ðθ 1 Þi can be found by substituting Δθ ¼ 0 in Eq. (9). This yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 6 3 ; 5 2 8 hα 2 where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 6 3 ; 4 7 7 Note that this yields numerical values closely matching those from the Z-tilt variance equation from Tyson, 22 expressed in terms of r 0 as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 6 3 ; 3 7 8 T 2 Substituting the tilt variance from Eq. (11), and the tilt correlation from Eq. (8), into Eq. (10), the differential Z-tilt variance can be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 6 3 ; where the weighting function is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 6 3 ; 2 5 3 Figure 3 shows a plot of f d ðzÞ for different source separations. The aperture size and path length used in the evaluation are same as those in Fig. 2. Note that the weighting functions in Fig. 3 drop down to zero at both ends of the path. The zero weight at the source occurs because we have a point source with spherical wave emanating from it. Angular variation due to turbulence immediately at the source is effectively like rotating the point source. This just directs a different ray emanating from the point source to the observer. It does not change the angle of arrival of the point source observed by the camera. On the other side of the path, tilts due to turbulence near the camera tend to be very similar across the field of view, owing to the convergence of the optical paths near the camera. This causes the differential signal to drop to zero at the camera end. It is also evident from Fig. 2 that the DTV grows with increasing angular separation between the sources. This is because the tilt correlation drops and the DTV approaches 2× the tilt variance from Eq. (10). Furthermore, it is interesting to compare Fig. 3 to Fig. 1. Note that for small separations, the weighting function in Fig. 3 is weighted more heavily toward the source. This is also the case for the isoplanatic angle weighting in Fig. 1, which relates to small scale anisoplanatism. For large separations, the weighting in Fig. 3 appears to approach the r 0 weighting in Fig. 1. This is a satisfying result, given that r 0 governs tilt variance as seen in Eq. (13).

Overview of Simulation
As mentioned in Sec. 1, the proposed method is based on the method of Bos and Roggemann. 7 Extended phase screens are generated as shown in Fig. 4. Points in the object plane are projected to the center of the camera pupil. Two such examples are shown in Fig. 4. The local phase screens are cropped from the extended phase screens within a specified distance of the optical path for each point. These local phase screen portions are shown with the blue and green squares in Fig. 4. The extended phase screen sizes are determined based on the cropped phase screen size and the object size, as shown in Fig. 4. More will be said about these dimensions shortly.
A simulated point source is numerically propagated from the source to the pupil plane though the cropped phase screens for up to each point in the object. The grid of object Wavelength Weighting function for C 2 n ðzÞ profiles for differential Z -tilt variance 16 for L ¼ 7 km and D ¼ 0.2034 m. The angle difference between the two sources Δθ ¼ jθ 1 − θ 2 j is expressed here in terms of Nyquist sampled pixel spacings for a focal length of l ¼ 1.2 m and wavelength of λ ¼ 0.525 μm.
points is spaced according to the Nyquist sample spacing in the object plane. Nyquist spacing in the focal plane is given by δ f ¼ 1∕ð2ρ c Þ, and the Nyquist spacing at the object is δ o ¼ λL∕ð2DÞ. Our simulation includes a skip parameter, whereby we have the option to skip a specified number of Nyquist samples for the purposes of propagation, and then we interpolate the PSFs that are generated to obtain the complete set. A bilinear interpolation is employed here for each sample of the PSF, based on the samples from the 4 PSFs that surround the PSF being interpolated. Once the PSFs are all generated, the output image is generated by using PSFs as weights in a spatially varying weighted sum of pixels from an ideal image. 7 The propagation method described in Sec. 3.2 is applied to each of a grid of points in the object plane. The point source wave functions are propagated through the phase screens and to the pupil plane. The incoherent PSF is then computed as described in Sec. 3.3. The phase screens are generated with the appropriate statistics as described in Sec. 3.4. Finally, in Sec. 3.5, we describe how the individual phase screen Fried parameters are determined in our simulation.

Numerical Wave Propagation
The split-step propagation method used to form each PSF is illustrated in Fig. 5. It involves a point source and a series of N phase screens that have been cropped based on the geometry shown in Fig. 4. Note that the path from the point in the object plane to the center of the camera aperture forms the centers of the cropped phase screen windows. We use nearest-neighbor interpolation here to speed computation. However, other forms of interpolation may be used.
The real valued wave function for the propagating field in the i'th ðx; yÞ plane along the z-axis is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 6 3 ; 1 5 4 where ν is the frequency of the source. The complex amplitude phasor of the wave function in the i'th plane is given by Our spherical wave propagation begins with a point source u 0 ðx; yÞ. We model this as a 2-D Gaussian windowed sinc function with quadratic phase. 2 This is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 3 2 6 ; 4 8 8 where α ¼D∕ðλLÞ, andD is the width of the amplitude field at the pupil plane where we want a uniform amplitude. The parameterD should be greater than D and less than phase screen width. This point source is designed to produce a flat amplitude over aD ×D patch at the pupil plane. 2 We can express the free-space propagation, plus phase screen perturbation, 2 as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 3 2 6 ; 3 7 6 for i ¼ 1;2; : : : ; N. The term h Δz i ðx; yÞ is the impulse response of free-space propagation 2 for Fresnel diffraction. This is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 3 2 6 ; 3 1 0 The term ϕ i ðx; yÞ in Eq. (19) is the phase imparted from the i'th phase screen realization. Note that we implement Eq. (19) on discrete space data using fast Fourier transforms (FFTs). 2 In particular, the method described by Schmidt as angular-spectrum propagation is used. 2 Absorbing borders using a Gaussian window may be applied in simulations in which a significant amount of signal energy reaches the borders of the simulation area. 2 The sample spacing used to represent the phase screens and to implement Eq. (19) is based on Voelz's critical sampling rule (best use of bandwidth and spatial support). 23 This gives a sample spacing for the phase screens of E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 3 2 6 ; 1 3 3 Δx ¼ λL where N is the number of samples used in each dimension. Note that this sample spacing is for the phase screens for the  purposes of propagation and should not be confused with the Nyquist sample spacing for the camera. The cropped phase screen width is set to be a multiple of the aperture. 2 In particular, we use a phase screen width of where s is the multiplier parameter. Using Eq. (21), this means that N ¼ ðsDÞ 2 ∕ðλLÞ. We chose s to be as small as possible such that s ≥ 4 and N is a power of two (to speed up FFT computations). We also use a constant screen spacing so that all propagations can be achieved with the same impulse response in Eq. (20).

Incoherent Point Spread Function
Following the propagations, the complex amplitude at the pupil plane u N ðx; yÞ is obtained. This is multiplied by the camera aperture mask aðx; yÞ, and a collimation-type phase compensation is used to allow the lens operation to focus the image at the focal length. This yields The incoherent PSF can be found based on Eq. (22) where Pðu; vÞ ¼ FTfpðx; yÞg and FTf·g is the Fourier transform. In practice, the PSF is evaluated discretely using the FFT and then resampled to the focal plane Nyquist spacing for the camera. We also normalize the discrete PSFs to all have a sum of 1. 7 The final simulation output is computed with a spatially varying weighted sum as where oðm; nÞ is the ideal discrete object image, and h m;n ðk; lÞ is the Nyquist sampled discrete PSF associated with object sample ðm; nÞ. Note that using the collimation in Eq. (22), the PSF image is focused at the focal length and not the image distance. This creates a magnification of M sim ¼ −l∕L, where the in-focus modeled system would have a magnification of M mod ¼ l∕ðl − LÞ. For large ranges, the difference is negligible. For shorter distances, the magnification can be corrected during the resampling step.

Generating Phase Screen Realizations
Phase screen realizations are designed to follow a modified von Kármán phase power spectral density (PSD). 2 This PSD includes the Kolmogorov PSD as a special case, but has additional parametric flexibility. The modified von Kármán PSD is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 6 3 ; 1 3 2 S mvK where ρ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , ρ m ¼ 5.92∕ð2πl 0 Þ, and ρ 0 ¼ 1∕L 0 . Note that the radial spatial frequency ρ is measured in units of cycles/meter, and r 0 i is the i'th phase screen Fried parameter. The parameters l 0 and L 0 are the inner and outer turbulence scales, respectively. 2 Note that for l 0 ¼ 0 and L 0 ¼ ∞, the modified von Kármán PSD is equivalent to a Kolmogorov PSD. 2 An example of the modified von Kármán PSD is shown in Fig. 6 for l 0 ¼ :01 m, L 0 ¼ 100 m, and r 0 i ¼ 0.01 m.
Generating a realization of a random process with specified PSD can be done by generating white noise with constant PSD of 1 and filtering it with a frequency response that is the square root of the desired PSD. Thus, we begin realizing the modified von Kármán phase screens by generating anÑ ×Ñ array of independent and identically distributed Gaussian random samples with standard deviation of 1. Note thatÑ is the extended phase screen dimension, prior to cropping to a size of N for propagation. We wish for these samples to correspond to a constant PSD value of 1 over the range −f s ∕2 to f s ∕2, where f s ¼ 1∕Δx. This means the total power (i.e., the integral of the PSD) should be f 2 s . For our discrete samples, that means the variance should be f 2 s . Thus, we multiply the unit standard deviation samples by f s . Next, we filter the array with a discrete-space impulse invariant 24 version of a filter with frequency response given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 3 2 6 ; 4 8 8 Our method avoids the need for generating complex valued random samples as used by some other methods. All phase screens are computed to the size of the maximum screen required and then cropped to the needed dimensions. This gives us consistent spectral content for all screens, regardless of the size required by the geometry shown in Fig. 4. One can see in Fig. 6 that the most dynamic part of the modified von Kármán PSD lies at low spatial frequencies. Faithful generation of low spatial frequency content is essential in matching theoretical statistics such as tilt variance and the long-exposure OTF. Since the sampling process limits the resolution at which the spatial frequencies can be evaluated in Eqs. (25) and (26), special subharmonic methods may be required. 2,25 Note that evaluating Eqs. (25) and (26) for use with FFTs is done in frequency increments of 1∕ðN Δ x Þ. The subharmonic methods seek to produce appropriate spectral content at spatial frequencies below the 1∕ðN Δ x Þ limit. We employ a subharmonic generation method based on the technique presented by Schmidt. 2 However, we use a real-valued 2-D Fourier series representation of the subharmonic spatial phase realizations, made up of a sum of 2-D cosines, to simplify the computations. An example phase screen realization is shown in Fig. 7. This has been generated from the modified von Kármán PSD shown in Fig. 6. The scale is in units of radians, and an aperture of size D ¼ 0.2034 m is shown for reference. A single PSF generated from this phase screen and aperture is shown in Fig. 8. Note that the PSF is down sampled to the Nyquist pixel spacing for the camera model, prior to applying it as part of the spatially varying weighted sum operation in Eq. (24) to produce the output image.

Determining Phase Screen Parameters
The individual phase screen Fried parameters are determined such that they are consistent with the global Fried parameter, log-amplitude variance, and isoplanatic angle. This approach is based on the method presented by Schmidt, 2 but we have extended this to also include the isoplanatic angle. Based on the weightings shown in Fig. 1, we see that these three parameters give us a good balance of weighting along the optical path. This allows us to simulate variable C 2 n ðzÞ profiles more accurately.
Following Schmidt's method, 2 the individual screen Fried parameters are related to screen C 2 n i parameters using E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 6 3 ; 3 9 5 r 0 i ¼ ½0.423k 2 C 2 n i Δz i −3∕5 ; or equivalently E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 8 ; 3 2 6 ; 4 7 3 Using a discrete version of the Fried parameter from Eq. (2), and substituting using Eq. (28), we obtain E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 9 ; 3 2 6 ; 4 0 9r A similar approach using the isoplanatic angle from Eq. (3), yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 ; 3 2 6 ; 2 9 7θ Finally, repeating the process for the log-amplitude variance in Eq. (4), yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 1 ; 3 2 6 ; 1 8 4σ 2  : We use a constrained least squares optimization, based on minimizing kb − bk 2 , where b ¼ ½ðr 0;sw Þ −5∕3 ; 6.8794L 5∕3 T , andb is the left side of Eq. (32). Our implementation uses the MATLAB function "fmincon" to perform the minimization. The constraints are user controlled in our simulation tool. For the results presented here, we use two constraints. One constraint is that the log-amplitude variance generated by any phase screen cannot exceed 20% of the total. The other constraint is that the r 0 N ¼ ∞ (i.e., no screen at the pupil). Since all point sources converge to the center of the pupil plane, they all share the exact same pupil phase screen. This tends to create excess tilt correlation. This might be expected by looking at the weighting in Fig. 2. Note that this is an artifact of the discrete phase screen approach. With more screens, the common pupil plane screen will represent a smaller Δz and cause less excess correlation. However, by setting the final phase screen phases to zero and using Eq. (32), we find that we are able to achieve a good match with the theoretical predictions, even with a relatively small number of screens.
We initialize the minimization search with Eq. (27), where C 2 n i are the average C 2 n ðzÞ for the section of the C 2 n ðzÞ profile represented by the i'th screen, at z ¼ z i . We found that we are able to get reasonable results using these initial values for the simulation (especially for a large number of screens). However, we observed a better overall match with the validation metrics using the optimization.

Experimental Results and Validation
The experimental results presented have been generated using the optical parameters listed in Table 1 and simulation parameters in Table 2. Note that we are simulating a range of 7 km for a visible wavelength telescope and camera with a wavelength of λ ¼ 0.525 μm. We have elected to use N ¼ 10 phase screens (9 nonzero phase screens). We have chosen a very large outer scale, since we are validating against the Kolmogorov-based statistics in Sec. 2. Thus, for these results, the modified von Kármán PSD is acting essentially as a Kolmogorov PSD. The images are of size 257 × 257 pixels, and spatial sampling at the Nyquist frequency is used (based on the diffraction-limited optical cut-off frequency). Note that other degradations such as downsampling and noise can easily be added after the turbulence simulation. We use a pixel skip parameter of four pixels to speed up the processing by a factor of 16×, compared with generating one PSF for each pixel. This can be adjusted according to the needs of the simulation and processing resources available. Our first set of results, in Sec. 4.1, is for a constant C 2 n ðzÞ profile, and then we consider varying profiles in Sec. 4.2. For each simulation, 200 frames are generated, and the theoretical validation metrics are  4.1 Constant C 2 n ðzÞ Profile We have simulated a constant C 2 n ðzÞ profile with five levels of turbulence, ranging from C 2 n ¼ 0.1 × 10 −15 m −2∕3 to C 2 n ¼ 1.5 × 10 −15 m −2∕3 . Some of the validation results are presented in Table 3. The validation metrics listed are r 0 , θ 0 , and the root-mean-squared (RMS) Z-tilt. The Fried parameter is estimated from the simulation long-exposure PSF, obtained by averaging the generated PSFs over frames. The isoplanatic angle is estimated by an analysis of the wave function phases over the aperture for point sources of varying angular separations. Finally, the Z-tilt is estimated by performing a correlation-based registration of the simulated PSFs. We have found that this type of correlation PSF registration corresponds well with Z-tilt, and PSF centroids correspond well to G-tilt. Note that we see a high level of agreement with regard to all of the metrics in Table 3. Optical Engineering 071502-9 July 2017 • Vol. 56 (7) The phase screen Fried parameters, obtained using the procedure described in Sec. 3.5, are shown in Fig. 9 for two levels of turbulence. No value is plotted for the pupil location at z ¼ 7 km since it would be infinity, as described in Sec. 3.5. Note that even though we are simulating a constant C 2 n ðzÞ profile, the values in Fig. 9 are not constant. This is a consequence of the optimization procedure. One can see that the phase screen immediately before the pupil plane does have a lower r 0 i than the others. This is in response to the pupil plane screen constraint, and this final nonzero screen is effectively "responsible" for a larger portion of the optical path than the others.
The theoretical long-and short-exposure PSFs (with diffraction) are shown in Figs. 10 and 11, respectively. In particular, amplitude normalized cross sections of the theoretical and simulated PSFs are shown for two turbulence levels.  Table 1.  Table 1.
An excellent agreement can be seen here in all cases. Note that without subharmonics, the simulated long-exposure PSFs tend to be too narrow, and the RMS Z-tilts tend to be too low. The relationship between these two parameters is explored in depth by Hardie et al. 15 Structure functions of phase 2 are shown in Figs. 12. These curves show the average squared phase difference over the aperture for two points separated by an angle of Δθ. Note that the isoplanatic angle is defined to be the angle where these structure functions have a value of 1 radian. These plots show fairly good agreement between the theory and simulation. Thus, the simulation appears to be capable of accurately capturing small scale anisoplanatic behavior. Note that the deviation of the simulated and theoretical curves for high source separation angles is due in phase wrapping. Phase wrapping is an unavoidable consequence of FFT processing during propagation. This artificially limits the estimated structure functions of phase in the simulation.  It does not, however, compromise the generation of the PSFs in any way.
A comparison of theoretical and simulated Z-tilt correlations is shown in Fig. 13 for two turbulence levels. A similar comparison of the DTV is shown in Fig. 14. These curves show that the simulation is accurately producing the larger scale anisoplanatic behavior predicted by the theoretical expressions. Note that if a nonzero phase screen is placed at the pupil plane, we tend to see simulated correlation that is higher, and DTV that is lower, than the theoretical values. This is explained by the fact that all PSFs share the exact same phase screen at the pupil, because of the converging optical paths.
Let us now examine some image results. The ideal object image used in these simulations is shown in Fig. 15. Degraded images for two levels of turbulence are shown in Fig. 16. Also shown in Fig. 16 are the corresponding Z-tilt motion vectors, scaled by 2×. The level of blurring and warping clearly increases with C 2 n .

4.2
Variable C 2 n ðzÞ Profile Validation analysis for the variable C 2 n ðzÞ profiles is shown in Table 4. The corresponding C 2 n ðzÞ profiles are shown in Fig. 17. Path A has heavy turbulence at the source, and Path B has heavy turbulence at the camera. The screen Fried parameters for these two cases are shown in Fig. 18. The long-and short-exposure PSFs are shown in Figs. 19 and 20, respectively. The structure functions of phase are  shown in Fig. 21. The tilt correlation and DTV plots for the two paths are shown in Figs. 22 and 23, respectively. As with the constant C 2 n ðzÞ profile results, there is generally good agreement between theory and simulation. The most deviation is seen with the tilt correlation and DTV for Path A. We shall continue to investigate this in future work.
Image results for the variable C 2 n ðzÞ profiles are shown in Fig. 24. The blurring is much more significant for Path B (heavy turbulence at camera). This is supported by the much smaller r 0 as shown in Table 4. The tilt vectors are also largest for Path B. This is consistent with the higher RMS Z-tilt reported in Table 4. On the other hand, Path A yields a smaller isoplanatic angle.

Implementation
The simulation has been implemented in MATLAB 2016a and run on a PC with Intel(R) Xeon(R) CPU E5-2620 v3 at 2.40 GHz, 16 GB RAM, and running Windows 10. We employ a NVIDIA Quadro K-4200 GPU. We use the parallel  Table 1.
Optical Engineering 071502-13 July 2017 • Vol. 56 (7) computing toolbox for the PSF generation component, such that the FFTs are all computed on the GPU. This is done by simply assigning the corresponding numerical arrays to be "gpuArrays." We have taken care to code the components efficiently by minimizing the use of "for" loops and maximizing the use of array operations. We used MATLAB's profiler to look for bottle necks and made adjustments for improved processing speed. The algorithm component run times are listed in Table 5 using the parameters listed in Tables 1 and 2. Note that the GPU provides a 5.66× speed up of the PSF array computation. The algorithm processing times scales with the number of PSFs being generated. Processing larger images or using a smaller pixel skip will give a correspondingly larger run time.

Conclusions
We have presented a numerical wave propagation method for simulating imaging of an extended scene under anisoplanatic conditions. In the simulation, we compute an array of PSFs for a 2-D grid of points on the object plane. An ideal image is then degraded by applying the PSFs in a spatially varying weighted sum operation. This gives rise to a spatially varying blurring and warping degradation. The PSFs are generated by the propagation of point source through an array of phase screens. The optical path for each point is somewhat different, by virtue of its originating position on the object plane and passes through different portions of the phase screens. This produces distinct but spatially correlated PSFs. We have employed a modified approach for determining the phase screen Fried parameters that incorporates r 0 , θ 0 and the log-amplitude variance. To see if the resulting PSFs exhibit accurate anisoplanatic statistical properties, we have conducted an extensive validation analysis. This analysis shows that this simulation is capable of generating accurate anisoplanatic effects on both a small and large scale. Small scale anisoplanaticism is validated with the isoplanatic angle. Large scale anisoplanaticism is validated for the first time using tilt correlation, 16 as well as with a newly derived DTV statistic for spherical waves. We also find excellent agreement between the simulated and theoretical long-and short-exposure PSFs.
We have demonstrated the simulation's performance for both constant C 2 n ðzÞ profiles and varying profiles. We have also generated sequences of independent frames and sequences with temporally correlated frames. The temporal correlation is achieved by generating extended phase screens and translating these according to a specified wind speed. The portion of the screens between the object and camera for Theoretical Simulation (a) (b) Fig. 23 Comparison of theoretical and simulated differential Z -tilt variances for (a) Path A, heavy at source and (b) Path B, heavy at camera. These are plotted versus Δθ, which is expressed in terms of Nyquist pixel spacings here for the parameters in Table 1.  Table 1.
a given frame are used for that frame's PSF generation. We have included video results showing the results of the temporally correlated sequence generation. Based on our analysis, we believe that this simulation approach can be used effectively to study anisoplanatic optical turbulence and to aid in the development of image restoration methods.