Fourier Propagation Tool for Aberration Analysis and a Point Spread Function Calculation of Systems with Curved Focal Planes

This paper describes a new Fourier propagator for computing the impulse response of an optical system with a curved focal plane array, while including terms ignored in Fresnel and Fraunhofer calculations. The propagator includes a Rayleigh-Sommerfeld diffraction formula calculation from a distant point through the optical system to its image point predicted by geometric optics on a spherical surface. The propagator then approximates the neighboring field points via the traditional binomial approximation of the Taylor series expansion around that field point. This technique results in a propagator that combines the speed of a Fourier transform operation with the accuracy of the Rayleigh-Sommerfeld diffraction formula calculation and extends Fourier optics to cases where the receiver plane is a curved surface. Bounds on the phase error introduced by the approximations are derived, which show it should be more widely applicable than traditional Fresnel propagators. Guidance on how to sample the pupil and detector planes of a simulated imaging system is provided. This report concludes by showing examples of the diffraction patterns computed by the new technique compared to those computed using the Rayleigh-Sommerfeld technique in order to demonstrate the utility of the propagator.


Introduction
Wave-optics simulation of light propagation has become an important element in the design of optical systems [1,2]. With the emergence of curved focal plane arrays in cameras, the need for new simulation tools to predict the performance of these systems is growing. Ray tracing tools already available to optical designers, such as ZEMAX, have been modified to include diffraction effects and could be adapted for use in curved focal plane cameras [3], but cannot easily take advantage of the Fast Fourier Transform (FFT) to aid in computations with curved detector surfaces.
Historically, the application of Fourier optics has always included a paraxial approximation. This leads to both the Fresnel and Fraunhofer propagators [1]. The angular spectrum propagator has also been utilized for wave propagation, but like the Fresnel and Fraunhofer methods, it has yet to be applied to truly curved surfaces, although any of these propagators may be applied in an approximate sense with as yet unspecified error, if the curved surface is approximated as being made up of a collection of flat ones [4]. Fourier optics assumptions are important for many applications where a Fourier transform operation is used to simulate the propagation of an optical field. These applications include phase retrieval and wave optics simulations of the atmospheric turbulence effects [5,6]. If Fourier optics can be extended to include propagation to truly curved focal planes, the efficiencies afforded by these methods could aid in the design of these new curved focal plane sensors.
The proposed propagation tool is not designed to replace lens design tools, but to provide designers of adaptive optics systems a tool for simulating the effects of turbulence with the optics of the telescope if they employ curved focal planes, due to the relatively high efficiency of Fourier propagators compared to ray tracing codes. Currently, the propagation tools of choice in these scenarios for flat focal planes are the digital Fresnel propagator and the angular spectrum propagator [2,[4][5][7][8][9]. The Fresnel integral can be used to propagate a field arriving at the telescope pupil from a wave that has propagated through layers of atmospheric turbulence. This kind of propagation tool allows for optical aberrations of the telescope to be included in the pupil function of the system if known but will not compute the telescope aberrations. This strategy ignores the fact that many aberrations are field dependent, making the interface of a ray tracing code with the Fresnel propagator tedious, since aberrations would have to be computed for different tilt angles of the incoming wave. The new propagator introduced in this paper is derived for the case where the radius of curvature of the sensor is equal to the focal length of the system. This particular choice is useful for systems focused at infinity that may incorporate curved focal planes to reduce the effect of aberrations.
The remainder of the paper is organized in the following way. Section 2 describes the new hybrid Rayleigh-Sommerfeld Fourier transform propagator for curved focal planes. Section 3 includes an analysis of the new propagator's binomial approximation and presents limitations on its use. Section 4 compares the Rayleigh-Sommerfeld propagator to the new one in a scenario involving an imaging system with a small F# that has a curved focal plane, where F# is the ratio of the focal length of the system to the aperture diameter. Section 5 summarizes the findings of the research and draws conclusions from them.

Fourier optics approximation of Rayleigh-Sommerfeld diffraction to a curved surface
The impulse response of an optical system can be modeled via the Rayleigh-Sommerfeld diffraction formula as long as the propagation distance between the source and the pupil plane is much greater than a wavelength [1].
The Rayleigh-Sommerfeld diffraction equation relates the source field, U (x,y,0), in the pupil to the receiver field U (x',y',z) in the focal plane. In Eq. (1), z is the distance between the pupil and receiver planes in the direction of propagation, while R is the distance between a specific point on the pupil, (x, y), and a specific point on the receiver plane, (x',y'). The source is assumed to be monochromatic with wavelength λ and the aperture is assumed to be finite in size with a diameter of D meters.
The approximation to the Rayleigh-Sommerfeld diffraction formula proposed in this paper is similar to the one proposed by Fresnel, except that the distance between the pupil plane and the receiving plane is not constant. This leads to a distance calculation of the following form.
In this equation (x, y) represent coordinates in the pupil plane and (x', y') are coordinates in the plane where the diffracted radiation pattern from the pupil is being computed (receiver plane). If the receiving plane is a curved surface following the equation of a sphere, then z(x', y') is the z coordinate of the spherical focal plane array detector with radius of curvature, S, which follows the equation: where S-z o is the difference in the radius of curvature of the curved detector and the actual vertical distance from the center of the pupil to the center of the curved detector as shown in Fig. 1.
In order to simplify the following analysis, S-z o is set to zero and the focal length of the system is also set equal to S, to represent telescope systems designed to image at near infinite distances. is computed from its (x',y') coordinates in Eq. (3). S is the radius of curvature of the curved detector whose surface lies on a sphere centered at (0,0,z o -S). The oval region represents the pupil plane containing the center of the coordinate system but not necessarily the center of curvature of the detector unless the radius of curvature is equal to z o which is the vertical distance between the aperture and the center of the detector. Substituting Eq. (3) into Eq. (2) yields, Defining the following expression for R o , 2 Substituting this expression into (4) and simplifying produces the following expression: Just as in the development of the Fresnel integral, the binomial approximation is utilized to approximate the radical resulting in the following approximated value for R [1]: Later in Section 3, conditions for the validity of Eq. (7) are derived. These validity conditions show that a region around the optic axis exists over which the binomial approximation is valid. The approximated value for R is then substituted into the Rayleigh-Sommerfeld diffraction formula in Eq. (1) while approximating R as S for all the amplitude related terms: The propagation equation expressed in Eq. (8) is similar to a Fourier transform, except the denominator of the Fourier phase [the last complex phasor in Eq. (8)] contains a term that is a function of the variables of integration. In order to approximate the complex phase term in a way that will make this integral similar to a Fraunhofer integral, the deviation of it away from the Fraunhofer phase [1] is given by: If this term can be shown to be much less than 1 radian, then it can be approximated as a Fourier transform kernel and Eq. (8) can be computed using the Fast Fourier Transform as long as the binomial approximation in Eq. (7) is valid. Factoring common terms yields and expression for the phase error, A common denominator can be found for the term in parentheses, allowing for the numerator to be expressed as the difference between S and R o . Furthermore, as seen in Eq. (5), R o is always larger than S, so the magnitude of the phase error, θ err , is bounded above by: Since our goal is to establish an upper bound on the phase error, the difference between R o and S can be replaced by an approximated difference as long as the approximation yields a value for R o that is greater than its true value. An example of this is the binomial approximation [1] shown below, The approximated value of R o is always greater than the exact value and is easily shown to be so by the fact that the third term in its Taylor series expansion, which is not included in the binomial approximation is strictly negative [1]. Now substituting Eq. (12) into (11) and requiring the phase error to be much less than 1 radian yields: Using the triangle inequality to bound Eq. (13) from above for an optical system and substituting a discretized form for (x',y') in terms of pixels, n max , that have a pitch of λ(S)/(2D), while also substituting the maximum values that |x| and |y| can assume, yields the following validity condition: This condition can be further simplified for by bounding π/4 from above by 1 and substituting the aperture diameter divided by the propagation distance as the F# so that Eq. (14) can be expressed simply as: With a few notable exceptions, most astronomical telescopes possess large F#s, thus making this technique useful for many dozens to hundreds of pixels around the center of the PSF. This is the same criterion that was developed previously for a similar Fourier Transform propagator based on a Rayleigh-Sommerfeld calculation developed for propagation to a flat plane [10]. The difference between the two is in the validity of the binomial approximation, as the curved receiver plane reduces the phase error of the third term in the Taylor series expansion, thus making the technique proposed in this paper possess a different validity condition.
If the PSF is larger than the validity region, then this technique cannot necessarily be used reliably to compute the PSF of the optical system but in some cases the PSF is observed to be accurate well outside the valid region. The validity condition for the digital Fresnel propagator does not have the property that its phase error can be limited by the number of pixels used in the receiver plane [10].
Equation (8) is approximated by its Nyquist sampled version by substituting a sample size for x' and y' as λ(S)/2D and a sample size of L/N for sampling x and y (where L is greater than or equal to 2D): The region of validity computed in this section assumes that the binomial approximation made in Eq. (7) is valid, which may not always be the case. In the next section, the validity of this binomial approximation will be examined and compared to that of the Fresnel propagator.

Binomial approximation validity conditions
In this section, the binomial approximation will be examined to determine under what conditions the new propagator can be utilized, while maintaining a low phase error. The binomial approximation in general takes the form: One way in which the validity of this approximation can be evaluated, is by examining the next term in the series, b 2 /8, which when multiplied by the wavenumber converts it from units of distance to optical phase: The goal is to find the conditions that make this expression less than 1 radian, thus revealing when the approximation is valid. Choosing the value of R o to be equal to S bounds this expression from above, also expanding the square and substituting the maximum extent for the region of validity in the focal plane in terms of pixels derived in Eq. (15) for x' and y' yields: Taking the maximum extent of the aperture plane to be x max =y max =D/2, moving S to the right and recalling that n max is bounded above by F# 2 , the following bound for the propagation distance is obtained.
In this expression π/4 has been bounded above by 1, thus making Eq. (20) as simple to compute as possible. Even for systems possessing practical F#s, this condition is easy to meet. This criterion is much less stringent than the Fresnel criterion, which even along the optic axis requires [1]: In this expression, the propagation distance has to be greater than a factor related to an inverse of the wavelength as opposed to a distance that is proportional to the wavelength.
Notice that the inequality goes different ways, featuring an upper limit on F# for the new propagator and a lower limit on F# for the Fresnel propagator. So, for example, the Hubble Space Telescope has a 2.4-meter diameter primary mirror with a 57.6-meter focal length and collects light with a mean wavelength near 500 nm [11]. In this case, Eq. (22) would have the left (the F#) being 24 and the right side would be evaluated at 103.6. Equation (23) would have the left side being equal to 98.04 and the right side would again be 24. In this scenario, the new propagator would be valid, and the Fresnel propagator would not be. Such a comparison is completely academic; however, as the Fresnel propagator is not designed to be used with curved focal plane sensors and the new propagator is. This fact makes a direct comparison between them superfluous as one should never choose to use the new propagator with a flat focal plane array, nor would it make sense to use the Fresnel propagator to propagate to a curved surface, since it isn't designed to do this.

Comparison to Rayleigh-Sommerfeld propagation
In this section, the proposed Fourier propagator is used to simulate the impulse response of an optical system that possesses a spherical aberration that is the result of a plane wave entering from infinity along the optic axis. This result will be compared to the impulse response computed using the Rayleigh-Sommerfeld diffraction formula, modified to be exact for a curved receiver plane. The optical arrangement to be simulated by the propagator is chosen to be simple enough to realize with widely available optical components, yet the configuration is chosen to produce an aberration that is commonly found in fast optical systems [12]. This aberration is of general interest since it can serve to limit the resolution of an optical telescope and would interact with atmospheric aberrations in ways difficult to predict with ray tracing techniques. Figure 2 shows the optical arrangement in which the lens has a 5 cm diameter and has a focal length of 20 cm, which is also equal to both the distance between the lens and the detector, S 2 as well as the radius of curvature of the curved focal plane, S. The wavelength of the incoming light is chosen to be 550 nm.

Point spread function calculations
The source field presented to the lens, U s , is computed using a simple plane wave model with a lens transformation for a simple spherical lens added to the phase of the field [1]. The burden of making this exact calculation is not great since no summations are involved in the execution of this operation. In this simulation, the lens is modeled as a simple spherical lens with a focal length of 20 cm [1]. The total phase in the pupil utilized by the new propagator can be computed as the phase of the lens transformation and the phase term associated with the new propagator shown in Eq. (24).
where θ lens is the phase of the lens transformation. The number of points chosen to model the pupil is 512 making the sample size in the pupil 195.31 µm. Since the pupil contains a discontinuity, it contains infinite frequency content, so the Nyquist sampling theory cannot predict a sufficient sample rate to avoid the aliasing of this signal. Beyond sampling the amplitude of the pupil function, the phase of the pupil must be sampled properly in order to avoid significant aliasing effects. All of the techniques discussed in this paper feature the addition of a phase to the input field, U s . The combined phase of the input field, lens transformation and the propagation field cannot contain pixel to pixel transitions greater than π radians. This is due to the fact that this phase is used as the argument of a cosine representing the real part of the complex exponential and the argument of a sine function which represents the imaginary part. If the phase of the pixel to pixel transitions exceed π radians, then there will not be at least 2 samples over each period of the real and imaginary parts of the complex sinusoid in the phase of the pupil. The amplitude of the pupil field is modeled as a circle with a diameter of 5 cm. Digitally, this circle has a diameter of 256 pixels. The field U s is stored as an amplitude array and a phase array separately.
Utilizing this pupil field, the discrete approximation of the Rayleigh-Sommerfeld diffraction formula shown in Eq. (1) is used to compute the PSF for the arrangement shown in Fig. 2 using a sample size in the pupil mentioned above and a discrete pixel size of 1.1 µm in the detector plane. The amplitude and phase arrays representing U s are modified by the amplitude and phase changes dictated by the phase term in Eq. (1) associated with the range calculation. The impulse response of the system is computed by taking the squared magnitude of the field obtained from Eq. (1) and normalizing it so it sums to 1 is shown in Fig. 3.
In order to compute the shape of the PSF using the new algorithm, the FFT (Fast Fourier Transform) is utilized. Unlike in the Rayleigh-Sommerfeld calculation, the sample size cannot be arbitrarily set but is determined by the size necessary to make the phase term in Eq. (7) resemble the phase term realized in the FFT algorithm. This can be accomplished by setting the detector plane sample size, ∆ r , equal to [2]: where L is the size of a plane containing the pupil. In this way L can be chosen to accommodate a different sample size other than the critical sample period predicted by Nyquist theory by changing the size of the pupil plane from the size of the lens. This is accomplished by zero padding the matrix that contains the field in the pupil. For this experiment the camera detector pixels are 1.1 µm as L is equal to 10 cm and the wavelength of the light is 550 nm. Figure 5 shows a plot of the intensity along the center of the PSFs along the x-axis shown in Fig. 4. The percent error in the new propagator can be computed from this case by utilizing Eq. (26). In this case the error was found to be 5.5%. Although the Rayleigh-Sommerfeld propagator is more accurate than the new technique, the image in Fig. 3 required 7196 seconds in MATLAB to compute, while the new propagator took just 0.061 seconds on the same machine. The machine In this equation, I RS , is the intensity of the PSF computed using the Rayleigh-Sommerfeld algorithm, shown in Fig. 3, and I FT is the intensity of the PSF computed using the new algorithm, shown in Fig. 4. The computed error is an average over central 32 by 32 pixels, where the PSFs reside as the validity region is +/-16 pixels.

Conclusions
The proposed Fourier propagator is shown to produce a diffraction pattern that is a close match to that produced by the Rayleigh-Sommerfeld diffraction formula. The calculations in Section 3 show that the proposed propagator's validity can be expressed in terms of the size of a region in the detector plane. The size of that region is always a function of the F#. This is in sharp contrast to the validity condition for the Fresnel propagator, which is based on the ratio of the pupil diameter to the propagation distance from the pupil to the detector. If these conditions are not met, then the Fresnel propagation is not necessarily valid for any region in the detector plane and they seldom are for practical systems [1]. This new proposed propagation tool will have an impact on applications of Fourier optics to systems with curved sensor arrays. It will most likely have its largest impact on applications dealing with the simulation of atmospheric turbulence on imaging systems. The propagator allows the aberrations from the atmosphere to interact with field dependent aberrations within a telescope in a natural way that eliminates the need to compute telescope aberrations separately with ray tracing codes. Although this might be possible using the Rayleigh-Sommerfeld propagator, the time required to execute would be unacceptable in many cases. The use of propagators designed for use with flat planes may be utilized, but with unknown validity conditions. This kind of modeling capability allows for system performance for a conceptual system designed to image through turbulence to be more rapidly determined with the computational savings afforded by the FFT algorithm.
Future research into improving this technique will be focused on extending its application to multiple mirror telescope designs. This could be accomplished by developing a method for computing a lens transformation for the multiple mirror configurations that captures the phase delays along the path from the pupil entrance to the point in the detector at which the impulse response is being computed. If this can be done, the proposed Fourier propagator could replace traditional ray tracing codes used in the design of multiple mirror telescopes. More development work for this approach is required and more testing would be required to validate the results for aberrations other than spherical aberrations.

Funding
Air Force Office of Scientific Research (JON 19ENG203).

Disclosures
Any opinions expressed in this document are those of the author and do not necessarily represent those of the United States Air Force or the Department of Defense. The author declares no conflict of interests.