2.1.

Traditional Tilt Calculation

The first step of this algorithm, as outlined above, is to calculate the pupil-plane tilts from the SHWFS data in the traditional sense. Here, the observation-plane irradiance pattern shifts for each sub-aperture are calculated from their on-axis position and used along with the focal length of lenslets to estimate the local pupil-plane tilts. Without any extensions, these tilts are then used in a least-squares reconstructor to estimate continuous $\varphi$. Hereafter, this approach will be referred to as the “traditional” approach to which the “tolerant” approach, discussed in the next few sections, will be compared.

2.2.

Identify Affected Sub-Apertures

As shown in Ref. 10, when the $|\mathrm{\Delta }\varphi |$ caused by a shock wave intersecting the measurement aperture exceeds $\sim 0.5\pi \mathrm{rad}$, the shock-induced phase discontinuity in the pupil-plane of the SHWFS sub-aperture lenslets causes appreciable beam spreading in the resulting observation-plane irradiance patterns. Therefore, in order to identify which sub-aperture lenslets are affected by the shock, beam spread can be calculated for each of the sub-aperture observation-plane irradiance patterns. In this case, the second-moment beam width, otherwise referred to as $D4\sigma$, has proven promising. The equation for $D4\sigma$ along the $x$ axis, $D4{\sigma }_x$, is given as

Equation (2) is calculated for each sub-aperture observation-plane irradiance pattern and the results are normalized by the diffraction-limited spot size, $D_{\mathrm{DL}}$. For square apertures, $D_{\mathrm{DL}}$ is given as $D_{\mathrm{DL}}=2\lambda f/D$, and for circular apertures, $D_{\mathrm{DL}}$ is given as $D_{\mathrm{DL}}=2.44\lambda f/D$, where $\lambda$ is the wavelength of light, $f$ is the focal length of the lenslets, and $D$ is the diameter of the lenslets. With appropriate thresholding, affected sub-apertures can be identified—i.e., high beam spread identifies a shock-affected sub-aperture, while low beam spread does not.

To illustrate this last point, a phase discontinuity of $\mathrm{\Delta }\varphi =-\pi$ rad was applied to a simulated, uniform, complex-optical field. The phase discontinuity was at a 30 deg angle (measured counterclockwise) from the vertical and offset from the middle of the measurement aperture by half a sub-aperture diameter. This complex-optical field was then applied to an SHWFS model consisting of $16\times 16$ sub-apertures. The resultant SHWFS imagery is shown in the leftmost plot in Fig. 1. We see that along the line where the simulated shock occurs, the associated SHWFS-irradiance patterns bifurcate. The degree of beam spreading calculated by Eq. (2) is shown in the middle plot of Fig. 1. Irradiance patterns that exhibit sufficient beam spread to qualify for being affected by the shock can be identified by thresholding the beam-spread values presented in the middle plot of Fig. 1. The flagged sub-apertures that met the thresholding criteria are presented in the rightmost plot of Fig. 1. In the sections that follow, we discuss the processing steps required to replace the corrupted tilt information in pupil locations flagged by the beam-spread calculation.

Fig. 1

(a) SHWFS image, (b) calculated $D4\sigma$ for sub-apertures, and (c) flagged sub-apertures for an intersecting shock of $\mathrm{\Delta }\varphi =-\pi$ and an angle 30 deg from the vertical (CCW).

2.3.

Calculate Shock Angle

After the sub-apertures affected by the shock have been identified, the angle that the shock intersects the SHWFS measurement aperture is calculated. This paper explores two methods to achieve this:

2.4.

Improve Tilt Estimation Using Irradiance Patterns in Affected Sub-Apertures

This section discusses how the resulting irradiance patterns of the sub-apertures affected by the shock can be used to improve the local tilt information. In Ref. 10, it had been shown that the degree of irradiance pattern bifurcation is related to the pupil-plane $|\mathrm{\Delta }\varphi |$. To investigate this, phase discontinuities of varying strengths ranging from $-\pi$ to 0 rad in increments of $0.05\pi \mathrm{rad}$ were simulated in the middle of a pupil and a Fresnel propagator was used to map the simulated pupil-plane phase to observation-plane irradiance patterns. For each resultant observation-plane irradiance pattern, the amplitude of the largest peak, $I_1$, was divided by the amplitude of the second largest peak, $I_2$. These results are shown in Fig. 2 as black circle markers where the $x$-axis is $\mathrm{\Delta }\varphi /\pi$ and the $y$-axis is $I_1/I_2$.

Fig. 2

Relationship between bifurcated irradiance pattern and phase discontinuity which caused it for square sub-apertures.

When there was no phase discontinuity in the pupil plane, $I_1/I_2$ was the diffraction-limited peak ratio for a square aperture, 21.19.¹⁷ For this case, a one-dimensional slice of the two-dimensional irradiance pattern is presented in the top-right plot of Fig. 2. Conversely, as $|\mathrm{\Delta }\varphi |$ increases, $I_1/I_2$ decreases. The middle-right plot of Fig. 2 shows the one-dimensional irradiance pattern slice when a phase discontinuity of $|\mathrm{\Delta }\varphi |=0.67\pi$ rad was imposed across the center of the pupil. Here, we see that $I_2$ began to increase in amplitude relative to $I_1$. When $|\mathrm{\Delta }\varphi |=\pi$ rad, the phase upstream versus downstream of the discontinuity was perfectly out of phase, and a bifurcated irradiance pattern with equal amplitude peaks formed, i.e., $I_1=I_2$. This last point is shown in the bottom-right plot of Fig. 2.

Conveniently, it has been shown that for square apertures, $I_1/I_2$ can be related to the phase discontinuity in the pupil plane that caused the bifurcation through the following:

In this equation, the constant is the diffraction-limited $I_1/I_2$ ratio for square apertures.¹⁷ Using a fitting algorithm, an equivalent expression was found to work well for circular apertures, given as

It is worth noting that in this expression, the constant is not the diffraction-limited $I_1/I_2$ ratio for circular apertures. Although these equations were arrived at heuristically, they show good agreement with the numerical results obtained using the Fresnel propagator. The results using Eq. (3) are plotted as a green line in Fig. 2. For different phase discontinuity and SHWFS conditions, such as angled shocks and experimentally collected data, a different method for relating irradiance patterns to $|\mathrm{\Delta }\varphi |$ will need to be used. However, for the purposes of this paper and for introducing the theoretical construct behind a shock-wave tolerant phase reconstructor, this analytic development will suffice, and alternative approaches continue to be explored.

It should be noted that Eqs. (3) and (4) only discern the magnitude of the phase discontinuity in the pupil-plane between 0 and $\pi$, not whether $\mathrm{\Delta }\varphi$ is positive or negative. As shown in Ref. 10, the same $I_1/I_2$ ratio can be accomplished when the secondary peak is on either side of the primary peak. However, it was also shown that the location of the secondary peak in relation to the primary peak is directly related to the direction of the flow for the case of an intersecting shock, as well as whether $|\mathrm{\Delta }\varphi |$ is between 0 and $\pi$ or $\pi$ and $2\pi$ (c.f. figure 7 in Ref. 10). For the case of a shock, the density increases across the shock and in the direction of freestream flow causing the phase of the light downstream of the shock to lag relative to the phase of the light upstream of the shock. Therefore, by knowing the direction of the flow in relation to the SHWFS measurement aperture, the true $\mathrm{\Delta }\varphi$ can be determined. For example, when the flow is moving from left to right across the measurement aperture and the secondary peak in the observation-plane irradiance pattern is to the right of the primary peak, the $\mathrm{\Delta }\varphi$ calculated from either Eq. (3) or Eq. (4) is correct. However, if the secondary peak is to the left of the primary peak, then the true $\mathrm{\Delta }\varphi$ is $2\pi -\mathrm{\Delta }\varphi$. These behaviors are flipped if the flow direction is flipped.

After using Eq. (3) or Eq. (4) in conjunction with the direction of the flow to determine the phase discontinuity across the shock for the affected sub-apertures, this information can be used to replace the tilt estimates in each sub-aperture. This is accomplished in two steps. First, the calculated change in phase across the discontinuity is converted to a change in OPD using the relation, $\mathrm{\Delta }\mathrm{OPD}=(-\lambda /2\pi )*\mathrm{\Delta }\varphi$, which is then used to calculate the overall tilt, $\theta$, caused by and orthogonal to the shock using

Next, the $x$ and $y$ components of the corrected tilt are calculated as

2.5.

Replace and Reconstruct

After calculation of the improved $x$ and $y$ tilts for sub-apertures affected by an intersecting shock from Eqs. (6) and (7), this new tilt information is then used to replace the corrupted $x$ and $y$ tilt values found using the traditional approach. It is important to note that the tilt estimates for the unaffected sub-apertures remain unchanged. After this replacement is made, a least-squares reconstruction is applied.

3.1.

Simulated Data

In this section, three simple scenarios were simulated for both square and circular sub-aperture geometries to demonstrate the shock-tolerant approach outlined above. For the results that follow, a phase discontinuity was generated and varied from 0 to $2\pi$ rad in increments of $0.05\pi$. For simplicity, it was assumed that the SHWFS sub-apertures were illuminated by unit-amplitude, collimated light. Furthermore, it was also assumed that the shock intersects the measurement aperture at $\alpha =0\hspace{0.17em}\mathrm{deg}$ (perpendicular to the $x$ axis) and is located in the middle of a column of sub-apertures. All scenarios outlined in this section follow the procedure described in Ref. 18. The first scenario divided the measurement aperture into a $16\times 16$ grid of sub-apertures each with diameter $D=6.25$ cm. The second scenario used $32\times 32$ sub-apertures each with diameter $D=3.125$ cm. The third scenario used $64\times 64$ sub-apertures each with diameter $D=1.5625$ cm. The entire grid size for each scenario was maintained at $2048\times 2048$ grid points. For each case, a thin-lens transmittance function was then applied to each sub-aperture, and angular-spectrum propagation was used to obtain an irradiance pattern at focus. These conditions were simulated for both square and circular SHWFS lenslet shapes.

After generating simulated SHWFS frames for various $\mathrm{\Delta }\varphi$, multiple thresholding steps were taken to ensure minimal contamination from ambient light and diffraction inside the SHWFS sub-aperture locations. Next, both traditional and shock-tolerant approaches were used to generate $x$ and $y$ tilt matrices. The resulting tilts calculated using both approaches were then used in a least-squares reconstructor, in this case using the Southwell geometry,³ to obtain a continuous OPD. The reconstructed OPD was then used to estimate the continuous $\varphi$ from which $\mathrm{\Delta }\varphi$ was calculated by subtracting the phase values on either side of the shock ($\mathrm{\Delta }\varphi ={\varphi }_{\mathrm{DOWNSTREAM}}-{\varphi }_{\mathrm{UPSTREAM}}$) for both the traditional and modified approaches.

Figure 3 presents the $\mathrm{\Delta }\varphi$ results of the traditional and shock-tolerant approaches plotted as blue circle and red square markers, respectively. Here, the $x$ axes are the input $\mathrm{\Delta }\varphi /\pi$ of the created phase discontinuities and the $y$ axes are the output $\mathrm{\Delta }\varphi /\pi$ resulting from reconstructing the tilt values obtained using the traditional and modified approaches. Also plotted are black, dashed lines that represent the expected $\mathrm{\Delta }\varphi /\pi$. The top row of Fig. 3 presents the results obtained using square sub-apertures and the bottom row presents the results obtained using circular sub-apertures. Each column presents results obtained using a different number of sub-apertures. Specifically, the $16\times 16$, $32\times 32$, and $64\times 64$ sub-aperture results are presented in columns from left to right, respectively.

Fig. 3

Reconstructed $\mathrm{\Delta }\varphi /\pi$ compared to input $\mathrm{\Delta }\varphi /\pi$ for the traditional and shock-tolerant reconstruction approaches using simulated SHWFS imagery with (a) $16\times 16$ square sub-apertures, (b) $32\times 32$ square sub-apertures, (c) $64\times 64$ square sub-apertures, (d) $16\times 16$ circular sub-apertures, (e) $32\times 32$ circular sub-apertures, and (f) $64\times 64$ circular sub-apertures.

As mentioned earlier and described in Ref. 10, when the $|\mathrm{\Delta }\varphi |$ induced by the shock wave becomes greater than $0.5\pi$, the traditional phase reconstruction begins to underestimate the true $\mathrm{\Delta }\varphi$ across the shock. This is evident for all cases presented in Fig. 3. When the shock-tolerant reconstruction algorithm was employed, we see that across all cases, the reconstructed $\mathrm{\Delta }\varphi$ is in better agreement with the true $\mathrm{\Delta }\varphi$. Furthermore, the results presented in Fig. 3 also reveal that as the number of sub-apertures increase, both the shock-tolerant and traditional least-squares reconstruction algorithms deliver better estimations of $\mathrm{\Delta }\varphi$, in their own regard. With more sub-apertures, there is more local tilt information available and least-square fitting error is minimized.

The sinusoidal behavior of the traditional method, shown in blue circles in Fig. 3, is a direct result of the observation-plane irradiance pattern bifurcation. When this bifurcation increases, the irradiance pattern centroid deflection and the corresponding pupil-plane tilt are no longer in agreement for shock-affected sub-apertures. Rather, the energy associated with the pupil-plane phase discontinuity leads to beam spreading in the observation plane. For the case of a shock-wave-induced phase discontinuity, the beam spreading is predictable, as shown in Sec. 2.4. When the modified approach described above was employed, which leverages the predictability of the shock-induced observation-plane irradiance pattern beam spread, the reconstructed $\mathrm{\Delta }\varphi /\pi$ was in fair agreement with the input $\mathrm{\Delta }\varphi /\pi$ until $\mathrm{\Delta }\varphi /\pi \approx 1.8$, depending on the sub-aperture geometry. Above $\mathrm{\Delta }\varphi /\pi \approx 1.8$, the $D4\sigma$ thresholding algorithm cannot discern which sub-aperture irradiance patterns are affected by the shock due to low beam spread. Therefore, the $x$ and $y$ tilt values calculated using the traditional approach are used instead.

3.2.

Experimental Data

This section demonstrates the shock-tolerant phase reconstruction algorithm on experimentally collected SHWFS imagery. In order to acquire the experimental data, a 532 nm laser source was first linearly polarized using a $\lambda /4$ wave-plate, expanded, and collimated. The beam was then reflected off a spatial-light modulator (SLM). The SLM was used to impart phase discontinuities ranging from 0 to $2\pi$ rad in increments of $0.05\pi$ rad at a 0 deg angle in the middle of the beam. In combination with an iris and collimating lens after the SLM, a focus and tilt term were also applied to the SLM in order to decouple the zeroth-order diffractive term from the first-order diffractive term which contains the applied phase information. After which, the beam was re-imaged onto an SHWFS using a $4f$ relay consisting of two 300 mm lenses. A schematic of this experimental setup is presented in Fig. 4.

Fig. 4

Experimental setup for the SHWFS data collection (not drawn to scale).

The SHWFS used $18\times 18$ square, sub-apertures with a diameter of $D=0.148\mathrm{mm}$ and a focal length of $f=6.7\mathrm{mm}$. The detector screen size was $360\times 360$ pixels, so each lenslet was comprised of $20\times 20$ pixels. The beam-spread threshold used for this data was $D4\sigma /D_{\mathrm{DL}}=2.9$. Figure 5 shows the resultant SHWFS imagery for phase discontinuities of 0, $\pi /2$, and $\pi$ as well as the sub-apertures that qualify, which are highlighted by red squares.

Fig. 5

Experimentally collected SHWFS imagery for (a) $\mathrm{\Delta }\varphi =0\mathrm{rad}$, (b) $\mathrm{\Delta }\varphi =\pi /2\mathrm{rad}$, (c) $\mathrm{\Delta }\varphi =\pi \mathrm{rad}$ phase discontinuities.

Due to realistic experimental factors such as non-uniform pupil illumination and camera measurement noise, Eq. (3) could not be used to estimate the shock-induced $\mathrm{\Delta }\varphi$ from the ratio of bifurcated irradiance pattern peaks. Instead, another equation was developed using a fitting algorithm to empirically collected data. After which, the experimental SHWFS measurements were processed using the traditional and shock-tolerant reconstruction approaches. These results are presented in Fig. 6.

Fig. 6

Reconstructed $\mathrm{\Delta }\varphi /\pi$ compared to input $\mathrm{\Delta }\varphi /\pi$ for the traditional and shock-tolerant reconstruction approaches using experimentally collected SHWFS imagery.

Here, it can be seen that the traditional reconstruction results exhibit the same behavior observed in Fig. 3 as well as in Ref. 10. Furthermore, we see that the shock-tolerant reconstruction results presented in Fig. 6 are in fair agreement with the results presented in the top-right plot of Fig. 3, where the number of sub-apertures was comparable to the number of sub-apertures used in the experiment. We expect that if the number of sub-apertures increased, the $\mathrm{\Delta }\varphi$ measured after employing the shock-tolerant reconstruction algorithm would likely be closer to the expected $\mathrm{\Delta }\varphi$. It is also worth noting that for $\mathrm{\Delta }\varphi /\pi$ values close to 1.35, the shock-tolerant reconstructed $\mathrm{\Delta }\varphi /\pi$ begins to level off. This is a result of the Gaussian beam used for the experiment. At $\mathrm{\Delta }\varphi /\pi$ values where the beam spread is low, the algorithm cannot discern between affected and unaffected sub-apertures based on the observation-plane irradiance pattern. This is especially true for those lenslets that lie close to the edge of the beam. The non-uniform illumination caused by the Gaussian beam profile leads to fewer detections of sub-apertures affected by the shock and consequently, more fitting error in the least-squares reconstruction. Nonetheless, the shock-tolerant algorithm is in fair agreement with the known $\mathrm{\Delta }\varphi$ until $\mathrm{\Delta }\varphi \approx 1.6\pi$ rad, where past this, the traditionally collected $x$ and $y$ tilt values were used.