% script to show the Airy disk
N=200;   % create a 200x200 matrix 
rho = ones(N,N);        % initiize a matrix 
for j= 1:N;              % calculate the distance from the origin
    for k= 1:N;
        rho(j,k) = ((j-N/2)^2+(k-N/2)^2)^0.5;
    end;
end;
rho = 10.0/(N/2)*rho;   % extend to rho = 10.0
airyd = (2*besselj(1,rho)./rho).^2;  %  Airy disk intensity
airyd(N/2,N/2) = 1.0;  % fudge to correct divide by zero, not-a-number
imshow(airyd); 
pause;

% difficult to see much on the linear scale so plot again using a log
% intensity scale
airylog = log10(airyd);        % base 10 log of the intensity
airylog = 0.333*airylog+1.0;   % rescale to show three decades
imshow(airylog);