# The Airy Disk of a diffraction limited system

The Airy Disk describes the blur arising from an ideal (diffraction-limited) uniformly illuminated, circular aperture. Its name arises from the astronomer, Lord Airy (George Biddell Airy)

This script exposes the code used to calculate and plot the Airy disk for diffraction limited optics.

## Contents

```ieInit
```

## Establish a scene and a diffraction limited optics

```scene = sceneCreate;
oi    = oiCreate;
oi    = oiCompute(scene,oi);

% Pull out the optics
optics = oiGet(oi,'optics');
```

## Select spatial samples and wavelength for plotting

```nSamp = 25; thisWave = 500; units = 'um';

clear fSupport
val = opticsGet(optics,'dlFSupport',thisWave,units,nSamp);
[fSupport(:,:,1),fSupport(:,:,2)] = meshgrid(val{1},val{2});

%Over sample to make a smooth image. This move increases the spatial
%frequency resolution (highest spatial frequency) by a factor of 4.
fSupport = fSupport*4;

% Frequency units are cycles/micron. The spatial frequency support runs
% from -Nyquist:Nyquist. With this support, the Nyquist frequency is
% actually the highest (peak) frequency value. There are two samples per
% Nyquist, so the sample spacing is 1/(2*peakF)
%
peakF = max(fSupport(:));
deltaSpace = 1/(2*peakF);
```

## Calculate the OTF using diffraction limited MTF (dlMTF)

```otf = dlMTF(oi,fSupport,thisWave,units);

% Derive the psf from the OTF
psf = fftshift(ifft2(otf));

% Make the spatial support for the PSF
clear sSupport
samp = (-nSamp:(nSamp-1));
[X,Y] = meshgrid(samp,samp);
sSupport(:,:,1) = X*deltaSpace;
sSupport(:,:,2) = Y*deltaSpace;

% Calculate the Airy disk
fNumber = opticsGet(optics,'fNumber');

% This is the Airy disk radius, by formula

% Draw a circle at the first zero crossing (Airy disk)
nCircleSamples = 200;
```

## Plot the diffraction limited PSF.

```x = sSupport(:,:,1); y = sSupport(:,:,2);
vcNewGraphWin;
mesh(x,y,psf);
colormap(jet)

% Label the graph and draw the Airy disk
ringZ = max(psf(:))*1e-3;
set(p,'linewidth',3); hold off;
xlabel('Position (um)'); ylabel('Position (um)');
```f = oiPlot(oi,'psf',[],500,'um');