This article explains:
The article is accompanied by a ZIP archive containing the sample ZEMAX files used. This can be downloaded from the final page of the article.
This article is also available in Japanese.
This article is also available in Japanese.
The PSF of an optical system is the irradiance distribution that results from a single point source in object space. A telescope forming an image of a distant star is a good example: the star is so far away that for all practical purposes it can be considered a point.
Although the source may be a point, the image is not. There are two main reasons. First, aberrations in the optical system will spread the image over a finite area. Second, diffraction effects will also spread the image, even in a system that has no aberrations.
ZEMAX has three basic types of PSF calculations: a geometric (no diffraction) spot diagram, and the diffraction based FFT and Huygens PSF. This article will discuss the basic theory and provide some guidance as to the proper use of each type of PSF.
One of the most basic analysis features in ZEMAX is the spot diagram. This feature launches many rays from a single source point in object space, traces all the rays through the optical system, and plots the (x, y) coordinates of all the rays relative to some common reference.
The sample optical system used here is a single parabolic F/5 mirror with a focal length of 50 mm. The object is at infinity. This system is a simplified Newtonian telescope, and the sample file is PSF_Newtonian.ZMX. Here is what the optical system looks like:
The spot diagram for two field points, one on axis and the other at an angle of 2 degrees is shown below.
Note the spot diagram is a collection of points, with each point representing a single ray. There is no interaction or interference between the rays. The spot diagram is very effective at showing the effects of the geometric, or ray aberrations of the telescope. The off axis geometric PSF clearly shows the coma and astigmatism of the system. On axis however, the spot diagram predicts perfect imagery. Is this an accurate representation of the optical system performance? To answer this question for the spot diagram results, we need to compare the spot distribution to the diffraction limited response.
A quick way to compare the geometric aberrations to the diffraction limit is to add an Airy disk reference ellipse to the spot diagram. Open the settings box and change Show Scale to "Airy Disk":
Now the spot diagram will indicate the size of the Airy disk relative to the geometric spot distribution:
On axis, the spot is much smaller than the Airy Disk, while off axis the spot is much larger than the Airy Disk. This indicates the spot diagram is a useful and reasonable indicator of performance off axis only. To compute a more accurate PSF both on and off axis, a consideration of diffraction will be required.
Generally speaking, the spot diagram is useful if the aberrations are large compared to the diffraction limited performance of the system.
The Fast Fourier Transform (FFT) algorithm has been widely applied to frequency analysis of many electrical and optical systems. Conceptually, the FFT decomposes a spatial distribution into a frequency domain distribution. An excellent discussion of Fourier optics can be found in reference [1] at the end of this article. There is also a summary of diffraction theory in the chapter "Physical Optics Propagation" in the ZEMAX manual, reference [2]. Both of these references describe Fresnel and Fraunhofer diffraction theory.
Most optical imaging systems meet the simplifying assumptions necessary for the Fraunhofer diffraction theory used by the FFT PSF algorithm. The main assumptions are:
The FFT PSF of an optical system is computed as follows. A grid of rays is traced from the source point to the exit pupil. For each ray, the amplitude and optical path difference is used to compute the complex amplitude of a point on the wavefront grid at the exit pupil. The FFT of this grid, appropriately scaled, is then squared to yield the real valued PSF. If the computation is polychromatic, the PSF's are summed incoherently.
To compute an FFT PSF for a sequential system, choose Analysis > PSF > FFT PSF from the main ZEMAX menu. Below is shown a sample FFT PSF for the on-axis field point of the Newtonian telescope sample file. The settings have been modified from the default settings, and this will discussed shortly.
Note the familiar Airy Disk shape. This is the expected result for the Newtonian, which is aberration free for the on axis field point.
To generate the picture above, the FFT PSF settings dialog should look like this:
The sampling refers to the grid of rays traced to the entrance pupil. Internally, ZEMAX doubles the size of the grid, filling the region outside of the entance pupil with zero valued data. Because of this doubling, the output PSF is always on a grid with twice as many points as the sampling grid. If the aberrations are reasonably small, the region of interest is concentrated near the center of the plot. Rather than plot all of these near-zero amplitude points, the display grid may be selected to be smaller than the total grid computed.
There are many ways to display the same underlying PSF data. Try the settings shown below:
Note Display is 128 x 128, Field is 2, Type is Logarithmic, and the Show As is set to False Color. Here is the resulting PSF:
Conceptually, the Huygens PSF is computed by converting each ray on the Spot Diagram into a small plane wave. Recall that a ray models a small piece of a plane wave, and that the ray locally is normal to the wavefront in isotropic media. The plane wave has an amplitude, phase, and direction determined from data associated with the ray from which it is generated. The total irradiance at any point on the image surface may be determined by coherently summing all of the plane waves represented by all of the rays traced. The diffraction based PSF is given directly by this integration over all the rays.
Virtually all imaging systems meet the simplyifying assumptions necessary for computing the Huygens PSF. The assumptions are:
The Huygens PSF is not based upon the FFT. The net result is that the Huygens PSF is generally slower than the FFT PSF, but more accurate for those cases where the FFT PSF assumptions do not apply. The most common cases where the FFT PSF assumptions are questionable, and thus the Huygens PSF should be used are:
The Huygens PSF of an optical system is computed as follows. A grid of rays is traced from the source point to the image surface. For each ray, the amplitude, coordinates, direction cosines, and optical path difference is used to compute the complex amplitude of the plane wave incident at every point on the image space grid. A coherent sum for all rays is performed at every point in the image space grid. The intensity of each point in the image space grid is the square of the resulting complex amplitude sum. If the computation is polychromatic, the PSF's are summed incoherently.
To compute a Huygens PSF for a sequential system, choose Analysis > PSF > Huygens PSF from the main ZEMAX menu. The Huygens PSF may also be computed for Non-Sequential Component (NSC) systems and this will be discussed shortly. Note the FFT PSF cannot be computed for NSC systems.
The key user definable parameters for the Huygens PSF are the Pupil Sampling, Image Sampling, and Image Delta. These parameters can be set on the Huygens PSF settings dialog: Open the settings box and change the default settings as shown:

The Image Delta is the image point spacing in micrometers. The total size of the region where the PSF is computed is the product of the Image Delta and the Image Sampling.
Here is the Huygens PSF for the same Newtonian telescope on axis:
And off axis (change the Field to 2):
The larger the number of rays and image points, the greater the resolution and accuracy of the resulting PSF, at the expense of a longer computation time.
One way to visualize this integration process is to observe the effects of coherent summing of one ray at a time. This can be accomplished with a coherent detector in the Non-Sequential Components feature of ZEMAX. The sample file used here is HPSF_Integration.ZMX, and this file is included in the ZIP file that may be downloaded at the end of this article.

This file consists of an elliptical source, a singlet lens, and a detector rectangle object. The source generates random rays over a circular region. All the rays exit parallel to the local Z axis. This source models a collimated source or distant point source. Note the number of layout rays is set to 20, while the number of analysis rays is set to 1. This latter setting will allow one ray at a time to be traced as will be discussed shortly. The lens is a simple singlet, placed to bring the collimated rays to a good focus on the detector. The detector is defined to be an absorber, with 120 x 120 pixels.
Note the detector object Parameter 11, the “PSF Wave #”, is set to 1:

This special mode allows the detector to perform the coherent Huygens PSF integration. Every ray that strikes the detector is converted into a local plane wave that illuminates every pixel on the detector, and the coherent amplitude of the plane wave at every the pixel is added to the coherent amplitude already detected. This allows one ray at a time to be traced, if desired, so the effects of summing individual rays can be seen.
To see this integration, open the sample file HPSF_Integration, and then choose Analysis > Detectors > Ray Trace / Detector Control. Select “Auto Update” on the Detector Control. Now, click on “Trace”. Because the source defines only 1 analysis ray, a single random ray is traced and the detector is updated. Click on “Trace” again, without pressing “Clear Detectors” first, to trace a second ray. The two rays now traced will coherently interfere like two plane waves making an angle with respect to each other, resulting in a fringe pattern on the detector. Because the rays are random, the fringe pattern will be different every time, and so will not look exactly like the picture below.

Each time “Trace” is pressed, another ray is added to the sum. After 10 rays have been traced the diffraction PSF begins to emerge.
/
After about 40 rays, features such as the Airy ring can be seen forming.

It takes a few hundred rays before the PSF reasonably converges to the final result.
The only reason to trace one ray at a time is to visualize what the integration is doing. To trace many rays at once, close the Detector Control, then change the number of analysis rays on the source from 1 to 500 in the NSC Editor:

Now reopen the Detector Control, and press “Clear Detectors” then “Trace”. All 500 rays will be traced at once, and the resulting PSF will be displayed in the detector viewer window.

Even though the rays are randomly selected the PSF converges to the correct Airy pattern (this particular lens is diffraction limited).
Use the spot diagram when:
Use the FFT PSF when:
Use the Huygens PSF when:
This article has described the Spot Diagram, FFT PSF, and Huygens PSF. In summary:
References
1. Goodman, Joseph W., Introduction to Fourier Optics, McGraw Hill
2. ZEMAX Optical Design Program User's Guide, ZEMAX Development Corporation