## 1.

## Introduction

The problem of imaging in the presence of degrading optical turbulence conditions has been an ongoing area of active research for many years. Much recent research has focused on correcting imagery through image processing (dewarping and/or deconvolution) techniques.^{1}2.3.4.5.6.7.^{–}^{8} Paralleling this inverse problem is the direct simulation of turbulence on clear images to create artificially degraded imagery at known turbulence levels^{9}10.11.^{–}^{12} for benchmarking purposes. A third research area involves system intercomparisons for cost-effectiveness evaluation.^{13}14.^{–}^{15} In each case, undergirding the research is the core issue of establishing the baseline effect of turbulence on image quality, particularly for short-exposure (SE) imaging.

The standard method used by systems performance modelers^{13}14.15.16.^{–}^{17} had been Fried’s SE modulation transfer function (MTF).^{18}19.^{–}^{20} For example, Holst^{15} devotes an entire chapter to Fried’s SE MTF. This model was the common choice due to its simplicity, despite the existence of more accurate (and numerically intensive) techniques (e.g., Ref. 21).

Critics of Fried’s approach pointed to its high-frequency response failure when modeling high turbulence degradation influences.^{21}22.^{–}^{23} Nonetheless, the simplicity of Fried’s result suggested a re-examination of Fried’s approach might yield new findings. The resulting study^{24} examined the properties of an overlooked tilt-phase correlation term and developed improved expressions for diffraction influences on the turbulent phase structure function (PSF), angle-of-arrival variance, and therefore on the SE MTF. However, the numerical integration technique used was not optimal for evaluating the tilt-phase influence at combined high-angular frequency and high-turbulence conditions. Consequently, a modeling compromise was adopted based on Fried’s quadratic correction to the long-exposure (LE) MTF. This approach was subsequently critiqued by Charnotskii^{25}26.27.^{–}^{28} who argued its quadratic term would produce inaccurate results at high frequency.

This issue is addressed here through a new analysis of high-frequency turbulence effects and development of a new SE MTF analytic model. The new model exhibits diffraction-limited behavior at all frequencies while applying to a wider range of optical scenarios. The resulting model is compared with published calculations from a path-integral technique^{21} and a Markov method.^{26} The model’s Rytov approximation (RA)^{29} is also described, and criteria developed by Tatarski^{29} and Dashen^{30} are examined to determine the extent of validity of modeled conditions.

The remainder of the paper is designed as follows: Section 2 describes the propagation model used, its justification, and numerical processing techniques developed to handle high-turbulence conditions. Section 3 describes the analysis process used to translate the computed database into the new analytical SE MTF expression. Section 4 presents comparisons with alternative approaches, followed by conclusions in Sec. 5.

## 2.

## Propagation and Imaging Models

The propagation model used is based on the Rytov approximation,^{16}^{,}^{29} involving a multiplicative turbulent scattering effect, using the imaging geometry of Fig. 1. Monochromatic incoherent light of wavelength $\lambda $ emerges from the transverse object plane at $z=-L$. Photons pass through the system aperture’s lens at $z=0$ to reach the imaging plane at $z=R$. The system’s thin lens optics are considered focused on the object plane. The main propagation axis is denoted by $z$; transverse two dimensional (2-D) vectors are denoted by $\mathit{x}$.

In general, light from a point ${\mathit{r}}_{O}=({\mathit{x}}_{O},-L)$ is evaluated for scatter off the ensemble of turbulent fluctuations at $({\mathit{x}}^{\prime},{z}^{\prime})$, passage through the aperture at $({\mathit{x}}_{E},0)$ and its effect on the light at image point $({\mathit{x}}_{I},R)$. Each object point $({\mathit{x}}_{o},-L)$ will have a reciprocal image point at $({\mathit{x}}_{i},+R)$ following a chief ray through the center of the system entrance pupil. Due to the incoherent nature of the light, each photon travels independently, emitted as a quantum event, and arriving at the image plane based on conditional passage through the system entrance pupil. The instantaneous point spread function represents the electro-magnetic equivalent of the photon’s quantum probability of arrival at image point ${\mathit{x}}_{I}$, given the current optical turbulence pattern (considered frozen) in the region $-L<z<0$.

## 2.1.

### Propagation Model

The Rytov approximation considers a propagating scalar wavefunction

consisting of a zeroth-order free-space solution, ${U}_{0}$, plus a complex turbulence multiplier. ${U}_{0}$ contains the primary longitudinal phasor $\mathrm{exp}[ik(z+L)]$, where ${i}^{2}=-1$ and $k=2\pi /\lambda $ is the wavenumber, plus the phase information necessary to determine the zero-turbulence focal point of the source point in the image plane, and spherical wave phase information that will be focused by the system lens effect.Turbulent amplitude ($\ell $) and phase ($\varphi $) effects are imposed on the propagating wave via weak scattering from position-varying refractive index fluctuations

where $\langle n\rangle \approx 1$ is the expectation (mean) of the refractive index $n(\mathit{x},z)$. Turbulent fluctuations will be modeled (exceptions noted) assuming Kolmogorov turbulence, characterized by the power spectrum^{26}where $\kappa $ is the turbulent spatial frequency ${\mathrm{m}}^{-1}$. For natural turbulence, ${n}_{1}\ll 0.001$, such that the refractive index structure parameter, ${C}_{n}^{2}$, is typically less than ${10}^{-11}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{m}}^{-2/3}$, leading to forward scattering (in the $+z$ direction), and use of the parabolic wave equation with ${\partial}_{z}$ the on-axis derivative operator, and ${\nabla}_{\perp}^{2}$ the transverse Laplacian.

Using Eq. (1)’s model propagating wave, the RA method’s single-scattering approximation is used to evaluate $\ell $ and $\varphi $ instantaneous amplitude and phase perturbations for radiation emerging from a source at ${\mathit{r}}_{O}=({\mathit{x}}_{O},-L)$ at system entrance pupil points, ${\mathit{r}}_{E}=({\mathit{x}}_{E},0)$,

## (5)

$$\left[\begin{array}{c}\ell ({\mathit{r}}_{O};{\mathit{x}}_{E})\\ \varphi ({\mathit{r}}_{O};{\mathit{x}}_{E})\end{array}\right]\phantom{\rule{0ex}{0ex}}=\left[\begin{array}{c}\mathrm{Re}\\ \mathrm{Im}\end{array}\right]\left\{\frac{{k}^{2}}{2\pi}{\iiint}_{\mathrm{TV}}\frac{\mathcal{G}({\mathit{r}}_{O};{\mathit{r}}^{\prime}){n}_{1}({\mathit{r}}^{\prime})\mathcal{G}({\mathit{r}}^{\prime};{\mathit{r}}_{E})}{\mathcal{G}({\mathit{r}}_{O};{\mathit{r}}_{E})}\mathrm{d}{\mathit{r}}^{\prime}\right\}.$$This equation reflects scattering from turbulence at points ${\mathit{r}}^{\prime}=({\mathit{x}}^{\prime},{z}^{\prime})$, $-L<{z}^{\prime}<0$, using free-space Green’s function propagators^{16}

## (6)

$$\mathcal{G}({\mathit{r}}_{1};{\mathit{r}}_{2})=\frac{1}{({z}_{2}-{z}_{1})}\mathrm{exp}\left\{ik\right[({z}_{2}-{z}_{1})+\frac{{|{\mathit{x}}_{2}-{\mathit{x}}_{1}|}^{2}}{2({z}_{2}-{z}_{1})}\left]\right\},$$## 2.2.

### Modulation Transfer Function

Following the methodology of Ref. 24, we analyze the wave perturbations at the system’s entrance pupil associated with the instantaneous SE MTF for the source point ${\mathit{r}}_{O}$

## (7)

$${\widehat{M}}_{S}(\mathit{\omega},{\mathit{r}}_{O})=\frac{4}{\pi}\int \mathrm{d}{\mathit{x}}_{E}\tilde{U}({\mathit{r}}_{O};{\mathit{x}}_{E}){\tilde{U}}^{*}({\mathit{r}}_{O};{\mathit{x}}_{E}-D\mathit{\omega}).$$Here, the field variables reflect the situation just after truncation by the aperture

## (8)

$$\tilde{U}({\mathit{r}}_{O};{\mathit{x}}_{E})={W}_{D}({\mathit{x}}_{E})\mathrm{exp}\{\ell ({\mathit{r}}_{O};{\mathit{x}}_{E})\phantom{\rule{0ex}{0ex}}+i[\varphi ({\mathit{r}}_{O};{\mathit{x}}_{E})-\mathit{a}({\mathit{r}}_{O})\xb7{\mathit{x}}_{E}]\},$$## (9)

$${W}_{D}(\mathit{x})=\{\begin{array}{ll}1,& |\mathit{x}|<D/2,\\ 0,& |\mathit{x}|>D/2.\end{array}$$According to Fried’s SE theory, the instantaneous turbulent phase tilt is removed, based on variable $\mathit{a}$ denoting the linear moment of the instantaneous phase perturbation

## (10)

$$\mathit{a}({\mathit{r}}_{O})=\frac{64}{\pi {D}^{4}}\int \mathrm{d}{\mathit{x}}_{E}{W}_{D}({\mathit{x}}_{E})\varphi ({\mathit{r}}_{O};{\mathit{x}}_{E}){\mathit{x}}_{E}.$$Vector variable $\mathit{\omega}=\mathrm{\Omega}/{\mathrm{\Omega}}_{0}$ is a normalized angular frequency, where ${\mathrm{\Omega}}_{0}=D/\lambda $ denotes the diffraction-limited maximum angular frequency^{16} passed by a circular aperture, such that $0<|\mathit{\omega}|<1$. $D$ is the entrance pupil diameter, dimensionalizing $\mathit{\omega}$ temporarily, before $\mathit{u}={\mathit{x}}_{E}/D$ is used to transform integral (10) to dimensionless form.

The mean SE MTF sought is then evaluated by taking the expectation of Eq. (7), and $\omega =|\mathit{\omega}|$

## (11)

$${M}_{S}(\mathit{\omega})=\langle {\widehat{M}}_{S}(\mathit{\omega},{\mathit{r}}_{O})\rangle =\frac{4}{\pi}\int \mathrm{d}\mathit{x}{W}_{D}(\mathit{x}){W}_{D}(\mathit{x}-D\mathit{\omega})\phantom{\rule{0ex}{0ex}}\mathrm{exp}(-{P}_{0}-{P}_{1}-{P}_{2}+{P}_{3}),$$^{16}analysis (pp. 395–407).

For brevity of presention, only results are shown here, based on more complete descriptions in Refs. 20 and 24. First, amplitude and phase perturbations are considered independent, permitting the expectation of the amplitude term to be evaluated as

## (12)

$$\mathrm{exp}(-{P}_{0})=\langle \mathrm{exp}\{\ell (\mathit{x})+\ell (\mathit{x}-D\mathit{\omega})\}\rangle =\mathrm{exp}\{-\frac{1}{2}{\mathcal{D}}_{\ell}(D\mathit{\omega})\},$$^{31}The remaining three terms are phase related, involving phase variance, tilt variance, and tilt-phase terms. The second term is the phase variance term, appearing as

## (13)

$$\mathrm{exp}(-{P}_{1})=\langle \mathrm{exp}[i\varphi (\mathit{x})-i\varphi (\mathit{x}-D\mathit{\omega})]\rangle =\mathrm{exp}\{-\frac{1}{2}{\mathcal{D}}_{\varphi}(D\omega )\},$$^{32}the second-order WSF is source independent, evaluated by weighted path integral of ${C}_{n}^{2}(z)$,

^{31}therefore isoplanatic.

Assuming constant Kolmogorov turbulence throughout the following

## (14)

$${P}_{0}+{P}_{1}=\mathcal{D}(D\omega )/2=[{\mathcal{D}}_{\ell}(D\omega )+{\mathcal{D}}_{\varphi}(D\omega )]/2\phantom{\rule{0ex}{0ex}}={[D\omega /({\mathit{r}}_{0}/2.1)]}^{5/3}={(2.1X\omega )}^{5/3},$$^{20}

^{,}

^{24}

The remaining two terms reflect SE corrections to the LE MTF (due to ${P}_{0}+{P}_{1}$). The tilt variance term that can be written^{24}

Since ${P}_{0}$, ${P}_{1}$, and ${P}_{2}$ are independent of $\mathit{x}$, they do not affect the aperture integral of Eq. (11). The remaining tilt-phase correlation, written ^{24}

## (17)

$${P}_{3}(\mathit{x},\mathit{\omega})=\frac{32}{\pi {D}^{4}}\int \mathrm{d}{\mathit{x}}_{a}{W}_{D}({\mathit{x}}_{a})[{\mathit{x}}_{a}\xb7(D\mathit{\omega})]\phantom{\rule{0ex}{0ex}}[{\mathcal{D}}_{\varphi}(|\mathit{x}-D\mathit{\omega}-{\mathit{x}}_{a}|)-{\mathcal{D}}_{\varphi}(|\mathit{x}-{\mathit{x}}_{a}|)],$$## (18)

$${\mathcal{D}}_{\varphi}(D\mathit{u})/2=\mathcal{D}(D\mathit{u})\alpha (Q\mathit{u})/2\phantom{\rule{0ex}{0ex}}=[{(2.1X/Q)}^{5/3}][{(Q\mathit{u})}^{5/3}\alpha (Q\mathit{u})]\phantom{\rule{0ex}{0ex}}={f}_{1}(X/Q){f}_{2}(Q\mathit{u}),$$^{24}

## (19)

$$\alpha (Q\mathit{u})={\int}_{0}^{1}\frac{\mathrm{d}c}{(3/8)}{\int}_{0}^{\infty}\mathrm{d}\gamma \frac{[1-{J}_{0}(\gamma c)]}{1.118{\gamma}^{8/3}}\times {\mathrm{cos}}^{2}\left\{{\gamma}^{2}\frac{[c(1-c)]}{4\pi {(Q\mathit{u})}^{2}}\right\},$$^{24}

When turbulence is low ($X\to 0$), the turbulent exponential argument approaches zero, and the SE MTF approaches the behavior of the system MTF,

## (20)

$${M}_{0}(\omega )=\frac{4}{\pi}\int \mathrm{d}\mathit{u}{W}_{1}(\mathit{u}){W}_{1}(\mathit{u}-\mathit{\omega});\phantom{\rule{0ex}{0ex}}=\{\begin{array}{cc}(2/\pi )[{\mathrm{cos}}^{-1}(\omega )-\omega \sqrt{1-{\omega}^{2}}],& \omega <1,\\ 0,& \omega \ge 1,\end{array}$$When $X\ne 0$, the Eq. (11) integral was evaluated numerically using the form

## (21)

$${M}_{T}(Q,X,\omega )=\frac{4/\pi}{{M}_{0}(\omega )}\int \mathrm{d}\mathit{u}{W}_{1}(\mathit{u}-\mathit{\omega}){W}_{1}(\mathit{u})\times \mathrm{exp}[+{P}_{3}(\mathit{u},\mathit{\omega},Q,X)],$$^{24}evaluation of ${M}_{T}$ was handled using a relatively inefficient integration technique, but here, based on the PSF of Eq. (18), ${f}_{1}(X/Q)$ was removed from the exponential integral, allowing the remaining term

## (22)

$${F}_{3}(Q,\mathit{u},\omega )=\frac{32}{\pi}\int \mathrm{d}{\mathit{u}}_{a}{W}_{1}({\mathit{u}}_{a})[{\mathit{u}}_{a}\xb7\widehat{\mathit{j}}\omega ]\phantom{\rule{0ex}{0ex}}[{f}_{2}(Q|\mathit{u}-\widehat{\mathit{j}}\omega -{\mathit{u}}_{a}|)-{f}_{2}(Q|\mathit{u}-{\mathit{u}}_{a}|)],$$Third, the ${M}_{T}$ evaluations were modified to handle the impact of large ${f}_{1}(X/Q)$ effects inside the exponential through use of an offset factor $A$

## (23)

$${M}_{T}(Q,X,\omega )=\frac{4\mathrm{exp}(A)}{\pi {M}_{0}(\omega )}\int \mathrm{d}\mathit{u}{W}_{1}(\mathit{u}-\widehat{\mathit{j}}\omega ){W}_{1}(\mathit{u})\phantom{\rule{0ex}{0ex}}\mathrm{exp}[{f}_{1}(X/Q){F}_{3}(Q,\mathit{u},\omega )-A].$$By adjusting $A$, math overflow errors in evaluating the integral could always be eliminated. Thereafter, ${M}_{T}$ calculations could be interpreted in the form

Combining this term with the tilt variance effect of Eq. (16) yields the function

such that the SE MTF can be written## (27)

$${M}_{S}(Q,X,\omega )={M}_{0}(\omega )\mathrm{exp}\{-{(2.1X\omega )}^{5/3}\phantom{\rule{0ex}{0ex}}[1-V(Q,X,\omega ){\omega}^{1/3}]\}.$$Calculations of the $V(Q,X,\omega )$ function were performed for a wide range of conditions: $Q$ from ${2}^{-4}$ to ${2}^{+4}$ by factors of 2 (9 steps), $X$ from ${10}^{-1}$ to ${10}^{+3}$ in increments of tenths of decades (41 steps), and $\omega $ from 0.01 to 0.99 in steps of 0.01, and four additional steps at $\omega <0.02$ (103 steps), for a database of 38,007 calculations. The analysis of this database is described in Sec. 3, using these results to formulate a new analytic SE MTF.

Sections 5 through 7 of Ref. 24 provides further details of the analysis of Fried’s approach.

## 2.3.

### Range of Rytov Approximation Validity

Before examining the $V(Q,X,\omega )$ database, some discussion of the extent of validity of the Rytov method seems appropriate. Clearly, the Rytov approximation is applicable to coherent radiation studies under weak turbulence conditions, but is typically considered inappropriate for stronger turbulence scenarios, when the Rytov variance

exceeds 0.5 (spherical wave coefficient used). However, in light of subsequent comparisons with Charnotskii’s results^{21}

^{,}

^{26}in Sec. 4, where comparable results are achieved at Rytov variance values seemingly beyond the capability of the RA, one may question the validity of such restrictions for incoherent radiation, since temporal and spatial correlations do not exist between propagating photons.

Tatarski’s analysis^{29} only required that ${n}_{1}\ll 1$, and that perturbation field gradients be small relative to wave number $k$, conditions easily met. Similarly, Dashen^{30} developed full and partial saturation conditions for judging RA validity, formed as ratios of a turbulence parameter, ${\mathrm{\Phi}}_{\mathrm{\Delta}}$, to a diffraction parameter, ${\mathrm{\Omega}}_{\mathrm{\Delta}}$.

For the full saturation case, Dashen used,

where $\mathrm{\Lambda}$ models a single scatterer scale. For a turbulent outer scale of ${L}_{o}$, and a spectral model^{33}based on Kaimal’s 1968 Kansas experiment analysis,

^{34}$\mathrm{\Lambda}=4.36{L}_{o}$. For the turbulence parameter, Dashen used a phase variance metric also based on an outer-scale influenced covariance In this case, the RA is valid when ${\mathrm{\Phi}}_{\mathrm{\Delta}1}<{\mathrm{\Omega}}_{\mathrm{\Delta}1}$, or, easily met for common outer scales of 1 m or greater.

More stringent conditions were associated with a partial saturation case involving multiple scattering scales. Dashen considered the effect of imposing a scale size cutoff that would affect both turbulent and diffraction effects. For the SE imaging problem, the most natural cutoff is the finite aperture. Dashen’s ${\mathrm{\Phi}}_{\mathrm{\Delta}}^{2}$, corresponding to an RMS phase variance, is equivalent to one of Fried’s^{18} $\langle {\mathrm{\Delta}}_{j}^{2}\rangle $ mean aperture phase variance terms, where the order $j$ indicates the number of phase perturbation terms corrected during the imaging process. In SE imaging, piston, tip, and tilt phase-expansion terms do not alter image quality (corresponding to Fried’s order-L case), leading to the condition,

## (32)

$${\mathrm{\Phi}}_{\mathrm{\Delta}2}^{2}=\langle {\mathrm{\Delta}}_{L}^{2}\rangle =(4\times 0.00489){(2.1X)}^{5/3},$$Setting ${\mathrm{\Phi}}_{\mathrm{\Delta}2}<{\mathrm{\Omega}}_{\mathrm{\Delta}2}$, for a given $Q$, a maximum $X$ value results

This condition translates to a maximum Rytov variance function of $Q$

## (35)

$${\sigma}_{R,\mathrm{max}}^{2}(Q)=0.676{\left(\frac{{X}_{\mathrm{max}}}{Q}\right)}^{5/3}=14261{Q}^{7/3}.$$Typical terrestrial imaging scenarios involve $Q>1/4$, implying the RA remains valid up to high-turbulence conditions.

## 3.

## Analytic Model of Short-Exposure Atmospheric Modulation Transfer Function

Using the previous section’s data set, this section describes development of an analytic expression for $V(Q,X,\omega )$. This is possible due to significant correlations that exist between $V(Q,X,\omega )$ results at different $Q$ values. This behavior is illustrated in Fig. 2 for two sets of plots at $Q$ values of $1/16$ and 16, for a series of increasing values of $X$.

The key feature of these curves is that for moderate $Q$ values often encountered for typical terrestrial imaging scenarios, $V$ exceeds unity at low to moderate frequencies, resulting in enhanced behavior versus the LE response.

In developing a model of this behavior, it was first observed that $V(Q,X,0)=V(Q,X,1)={V}_{0}(Q)$, permitting $V$ to be expressed

The maximum value of $V$ (as $X\to \infty $) must also be a function of $Q$, ${V}_{M}(Q)$. ${V}_{0}(Q)$ and ${V}_{M}(Q)$ are plotted in Fig. 3.

Figure 3 curves both suggest asymptotic behaviors at both extremes of $Q$. These behaviors can be modeled using a sigmoidal function (a variant of the hyperbolic tangent)

## (37)

$$\mathrm{\Sigma}(\chi )=\frac{[\mathrm{exp}(\chi )-1]}{[\mathrm{exp}(\chi )+1]}=\mathrm{tanh}(\chi /2).$$To facilitate the modeling process, a rescaled and shifted sigmoid was also introduced

The center of this function occurs at $(\chi ,\eta )=(D,A)$ with an overall shift of $2B$ between asymptotes, and a central derivative of $d\eta /d\chi {|}_{\chi =D}=BC$.

The ${V}_{M}(Q)$ function can be expressed using ${\sum}_{S}$ as

However, ${V}_{0}(Q)$ exhibits different asymptotic behaviors at large and small $Q$. Therefore, a third sigmoidal form was developed, featuring a central splice (piecewise sigmoidal):

## (40)

$${\mathrm{\Sigma}}_{P}(\chi ,A,{B}_{1},{B}_{2},{C}_{1},{C}_{2},D)=\{\begin{array}{ll}A+{B}_{1}\mathrm{\Sigma}[{C}_{1}(\chi -D)],& \chi \le D\\ A+{B}_{2}\mathrm{\Sigma}[{C}_{2}(\chi -D)],& \chi \ge D\end{array}.$$## (41)

$${V}_{0}(Q)={\mathrm{\Sigma}}_{P}[{\mathrm{log}}_{2}(Q),0.870,0.370,0.085,0.355,1.545,-1.00],\phantom{\rule{0ex}{0ex}}$$Using Eqs. (39) and (41), our problem reduces to finding an analytic expression for

Figure 4 illustrates plots of ${\widehat{V}}_{\mathrm{\Delta}}$ for the same families of curves ($Q=1/16$ and 16) as in Fig. 2, highlighting the similarities in behaviors at different $Q$ values. These plots also reveal how ${\widehat{V}}_{\mathrm{\Delta}}$ develops increasingly complex behaviors as $X$ increases.

To model these behaviors, it was observed that function curves at high $\omega $ maintained simpler forms with increasing $X$ than at low $\omega $, suggesting the functions could be divided at the peak and studied separately at high and low frequencies. The next step was, therefore, to model the form of the peak curves (as illustrated in Fig. 4) as functions of $Q$ and $X$.

Note that the peak curves at different $Q$ values exhibit similar $X$ dependencies. Thus, the behaviors of the different peak curves could be modeled as shifted versions of the peak curve for $Q=1$, described by ${\omega}_{P}(1,X)$ (abscissa) and ${V}_{P}(1,X)$ (ordinate) functions. The $Q=1$ abscissa curve is modeled as

## (43)

$${\omega}_{P}(1,\mathit{x})={\mathrm{\Sigma}}_{P}[\mathit{x},0.360,-0.133,-0.322,\phantom{\rule{0ex}{0ex}}+3.450,+1.424,+1.320].$$Several ${\omega}_{P}(Q,X)$ curves at different $Q$ levels illustrate this shifting behavior in Fig. 5. A pair of analytic expressions that approximate this behavior are given by

## (44)

$${\omega}_{P}(Q,X)={\omega}_{P}\{Q=1,{\mathrm{log}}_{10}(X)-{\mathit{x}}_{\omega}[{\mathrm{log}}_{2}(Q)]\},$$Several ${V}_{P}$ ordinate function curves, designated ${V}_{P}(Q,X)={V}_{\mathrm{\Delta}}(Q,X,{\omega}_{P})$ (using the non-normalized form), are illustrated for varying $X$ in Fig. 6, for $Q$ values ranging from $1/16$ to 16. This component exhibits both shifting and scaling properties, but again was quantified relative to its behavior at $Q=1$, using

## (46)

$${V}_{1}(\mathit{x})={V}_{P}(1,{10}^{\mathit{x}})\phantom{\rule{0ex}{0ex}}={\mathrm{\Sigma}}_{P}(\mathit{x},+0.122,+0.044,+0.082,\phantom{\rule{0ex}{0ex}}+4.257,+2.300,+1.167).$$The behavior at $Q=1$ was then rescaled using ${A}_{R}$, and shifted using ${\mathit{x}}_{P}$, in the form,

## (47)

$${V}_{P}(Q,X)={A}_{R}[{\mathrm{log}}_{2}(Q)]{V}_{1}\{{\mathrm{log}}_{10}(X)-{\mathit{x}}_{P}[{\mathrm{log}}_{2}(Q)]\},$$The function ${V}_{P}(Q,X)$ supports a further parameterization through introduction of

using Eqs. (47) and (39), where $0<{U}_{P}\le 1$. The ${\widehat{V}}_{\mathrm{\Delta}}$ curve shapes appear to vary based on height ${U}_{P}$, while $\omega $ dependence may be modeled using a Legendre-polynomial-related series expansion.The Legendre polynomials begin as ${P}_{0}(z)=1$, ${P}_{1}(z)=z$, and use recurrence relation^{35}

Using this basis set, each ${\widehat{V}}_{\mathrm{\Delta}}$ curve could be separated into two functions at the peak position $\omega ={\omega}_{P}(Q,X)$, followed by affine transformation to stretch each monotonic function into its own $\omega $ half-range, and reverse copied to fill the opposing half-range. Thereby, two functions, each even about $\omega =1/2$, could be analyzed using Legendre-derived basis functions. Due to their symmetry, only even-order series terms were needed to characterize each function $[{L}_{2n}(\omega )={L}_{2n}(1-\omega )]$. The ${\widehat{V}}_{\mathrm{\Delta}}$ model could, thus, be splined using the coefficient sets: ${C}_{2n}$ for $\omega \le {\omega}_{P}$, and ${D}_{2n}$ for $\omega \le {\omega}_{P}$, as

## (53)

$${\widehat{V}}_{\mathrm{\Delta}}(Q,X,\omega )\approx \sum _{n=0}^{n=N}{C}_{2n}[{U}_{P}(Q,X)]{L}_{2n}\left\{\frac{\omega /2}{{\omega}_{P}(Q,X)}\right\},\phantom{\rule{0ex}{0ex}}\omega \le {\omega}_{P};$$## (54)

$${\widehat{V}}_{\mathrm{\Delta}}(Q,X,\omega )\approx \sum _{n=0}^{n=N}{D}_{2n}[{U}_{P}(Q,X)]{L}_{2n}\{1-\frac{[1-\omega ]/2}{[1-{\omega}_{P}(Q,X)]}\},\phantom{\rule{0ex}{0ex}}\omega \ge {\omega}_{P}.$$Sets of ${C}_{2n}({U}_{P})$ coefficient behaviors are illustrated in Fig. 7. ${D}_{2n}({U}_{P})$ behaviors are plotted in Fig. 8. Coefficients ${C}_{0}$ through ${C}_{8}$, and ${D}_{0}$, exhibit approximately hyperbolic ${U}_{P}$ dependence. The remaining coefficients exhibit approximately linear dependence.

Hyperbolic behaviors of ${C}_{m}$ and ${D}_{m}$ coefficients were expressed using constants ${a}_{m}$ through ${e}_{m}$ in the models

## (55)

$${C}_{m}({U}_{P})={a}_{m}({U}_{P}-{e}_{m})+{b}_{m}\sqrt{{({U}_{P}-{e}_{m})}^{2}+{c}_{m}}+{d}_{m}.$$Linear cases were modeled by setting ${b}_{m}={c}_{m}=\phantom{\rule{0ex}{0ex}}{e}_{m}=0$. The curves plotted in Figs. 7 and 8 represent functional fits to computed values. Fit coefficients are listed in Tables 1 and 2.

## Table 1

SPGD hyperbolic Cm and Dm coefficients.

Coefficient | am | bm | cm | dm | em |
---|---|---|---|---|---|

C0 | 20.2124 | 19.6324 | 0.0031 | 0.5751 | 1.0686 |

C2 | 2.1118 | 2.4928 | 0.0042 | −0.3515 | 0.9808 |

C4 | −0.3413 | −0.3355 | 0.0051 | 0.0051 | 0.7996 |

C6 | −0.3529 | −0.3769 | 0.0000 | 0.0131 | 0.9041 |

C8 | −0.3696 | −0.3982 | 0.0000 | 0.0153 | 0.9504 |

D0 | 6.4372 | 5.8919 | 0.0005 | 0.5788 | 1.0341 |

## Table 2

SPGD linear Cm and Dm coefficients.

Coefficient | am | dm |
---|---|---|

C10 | 0.01009 | −0.00452 |

C12 | 0.00550 | −0.00257 |

D2 | −0.30273 | −0.00057 |

D4 | 0.03335 | −0.00413 |

D6 | −0.01271 | 0.00485 |

D8 | 0.00425 | −0.00154 |

D10 | −0.00452 | 0.00186 |

D12 | 0.00055 | −0.00019 |

With the development of the ${C}_{m}$ and ${D}_{m}$ functional forms, the new analytic SE MTF model form was completed. However, given the number of free constants (80) involved, an open question remained whether a slightly modified set of coefficients could produce a better fit. To explore this possibility, a stochastic parallel gradient descent (SPGD) procedure^{36} was applied to improve the fit to the database.

The coefficient set derived in the initial analysis (described above) was used to “seed” the SPGD optimizer. The SPGD consisted of randomly perturbing the full 80-parameter coefficient set, recomputing the RMS $V(Q,X,\omega )$ error for the full database, and comparing that error with the previous best fit error. Any new set that outperformed the prior best fit was adopted as the new best fit. Use of the SPGD reduced the overall RMS error of $V$ from 0.006795 to 0.002728. For simplicity of presentation, coefficients listed in Tables 1 and 2 represent the SPGD best-fit coefficients. For the remaining expressions, the original results (used in the plotted figures) were retained in Eqs. (39) through (49), so curves in the various figures can be replicated. Optimized versions of these equations are given below.

To summarize, the new SE MTF model is formulated starting with the general expression in Eq. (27), then using ${M}_{0}$ from (20), $V$ from (36), ${V}_{\mathrm{\Delta}}$ from (42), expressions in Eqs. (44), (47), and (50) through (55), and results of the SPGD analysis

## (58)

$${\omega}_{P}(1,X)={\mathrm{\Sigma}}_{P}[\mathit{x},0.351,-0.133,-0.318,\phantom{\rule{0ex}{0ex}}+3.451,+1.407,+1.274],$$## 4.

## Comparisons with Alternative Approaches

The analytic model developed in the previous section is compared here with results obtained by Charnotskii^{21}^{,}^{26} and several prior analytic methods. Charnotskii^{21} used a Feynmann-path integral technique.

Figures 9 through 11 intercompare the new model with Charnotskii’s^{21} Figs. 8 through 10. Charnotskii’s $d$ parameter is related to $X$ via $X={({d}^{2}/2)}^{3/5}/2.1$, such that $d$ values 0, $\pi $, $2\pi $, $3\pi $, and $4\pi $ correspond to $X$ values 0.000, 1.241, 2.851, 4.638, and 6.550, respectively. Charnotskii’s $N$ parameter is related to $Q$ as $Q={(2N/\pi )}^{1/2}$. Reported $N$ values of the three plots were 10.0, 1.0, and 0.1, respectively, corresponding to $Q$ values 2.523, 0.798, and 0.252. The $d=0$ curves (the uppermost line in each graph) correspond to the system MTF, and match identically. Remaining lines in Figs. 9 and 11 were plotted using $Q$ values that produced the closest fit to Charnotskii’s results, $Q=12.0$ in Fig. 9 and $Q=0.02$ in Fig. 11. Lines in Fig. 10 were constructed based on the reported $N=1$ ($Q=0.798$). Systematic deviations employed in these plots appear consistent with Charnotskii comments that $N=10$ and $N=0.1$ constituted limiting cases, but might be due to approximations used when Charnotskii developed Eq. (31) of Ref. 21.

The comparisons of Figs. 9Fig. 10–11 provide partial verifications of both the present method and Charnotskii’s calculations. A similar comparison is made in Fig. 12 between data from Fig. 2 of Ref. 26 (symbols) and simulations (lines) based on,

## (63)

$${\mathcal{D}}_{SA}[Q,X,\mathit{r}/D]=2{(2.1X)}^{5/3}{(\mathit{r}/D)}^{5/3}\times \{1-V[Q,X,(\mathit{r}/D)]{(\mathit{r}/D)}^{1/3}\},$$A key element of these results is the range of Rytov variance where similar results were obtained using RA versus the alternative approaches of Refs. 21 and 26. While Figs. 9 and 11 involved $Q$ stretching, a single $Q$ adjustment appeared to correct all $X$ cases, while no $Q$ adjustment was used in Fig. 10 where the two strongest turbulence cases corresponded to ${\sigma}_{R}^{2}$ values of 12.7 and 22.8. This ability of the RA to perform at high ${\sigma}_{R}^{2}$ levels is consistent with the Sec. 2.3 analysis.

Lastly, four prior SE MTF models, (1) $V=0.5$ (Fried’s far-field case), (2) $V=1.0$ (Fried’s near-field case), (3) Holst’s^{15} method where either Fried’s near-field or far-field approximation was chosen based on best fit, and (4) the $V(Q,X)$ method of Ref. 24, were compared with the new analytic SE MTF model (5). An RMS metric was used to compare each estimate against database-derived MTF values, weighted by frequency $\omega $. The RMS results computed were ${E}_{1}=0.035179$, ${E}_{2}=0.014089$, ${E}_{3}=0.010788$, ${E}_{4}=\phantom{\rule{0ex}{0ex}}0.006431$, and ${E}_{5}=0.000218$. The prequel model (4) outperformed methods (1) through (3), but the new model offered dramatic improvement (by a factor of 29) over the prequel, and greater than 100 better than the far-field approach.

## 5.

## Conclusions

This paper’s extended high angular frequency analysis, using the methodology described in Sec. 2, has lead to a new analytic expression for the atmospheric SE MTF. This model’s RMS error of 0.000218 versus numerically integrated results represents factors of 29 improvement over the prequel,^{24} 49 over Holst,^{15} and 64 over Fried’s near-field case.^{20} The new model exhibits diffraction-limited behavior at all angular frequencies.

The primary focus of this effort was to provide systems performance engineers with an easily computed SE MTF that accounted for a thorough range of system and observational scenario conditions. Specifically, it improves evaluation of SE MTF degradation at high frequencies. This will aid trade studies in evaluating image reconstruction technique effectiveness where high frequency edge information losses are of critical concern.

Augmenting the present model, a previous study^{37} showed that the new model, based on structure functions, can treat path-varying turbulence effects, including slant-path and valley geometries. A planned future paper will expand on these results, considering the annular aperture case (Newtonian telescopes, etc.), where the circular aperture results represent the zero central obscuration limit of the more general solution.

The main caveats of the model regard the Kolmogorov spectrum used and the validity of the RA method. Regards the spectrum, we impose the requirement that ${\ell}_{o}<D<{L}_{o}$, ensuring $D$ occurs in the inertial subrange. Regarding the RA method, conditions derived based on Dashen’s ^{30} analysis [Eqs. (31) and (35)] suggest a relatively wide range of applicability. This suitability was further tested through intercomparison with path-integral^{21} and Markov^{26} techniques.

Results of the current study might also be used in an improved inverse filtering procedure [e.g., Ref. 5], though such a study is beyond the scope of the present work. In such a study, multiple range effects could be simulated, but special handling would be required to subsegment the images to apply different MTF’s in different portions of the image field.

## References

## Biography

**David H. Tofsted** is a physicist with the U.S. Army Research Laboratory at White Sands Missile Range, New Mexico, with over 33 years in service. He received his BS degree in physics from Pennsylvania State University, Pennsylvania, in 1979, and MS and PhD degrees in electrical engineering from New Mexico State University, New Mexico, in 1992 and 2013. His research has focused on optical propagation problems, including numerical simulations of optical turbulence in the atmospheric boundary layer and its impacts on optical imaging systems. He has authored over 60 publications.