% script to simulate a telescope image include effects of
% diffraction and the finite pixel size of the detector
% for example, the spy-plane problem of P.S. #7.

dist = 10e3;      % distance from telescope to license plate (in meters)
diam = 0.3;       % diamter of telescope objective lens (in meters)
focal = 10;       % focal length of the objective lens (in meters)
lambda = 0.5e-6;  % wavelength of light (in meteres)
pixel = 10e-3;    % size of pixel on the CCD camera (in mm)

% use a scaling of one pixel in the image matrix 
% is equal to 1 mm on the license plate
% eta is the fundamental scale of diffraction limited resolution (in mm)
eta = lambda/diam*dist*1000;
N=2*round(3*eta);       % create a NxN matrix for the airy disk that is
                        % 6 times the size of this resolution unit
rho = ones(N,N);        % initialize 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 = pi*rho/eta;
airyd = (2*besselj(1,rho)./rho).^2;  %  Airy disk intensity
%   this assumes that different parts of the image are incoherent with each
%   other; probably not always true but probably OK for
%   illumination by sunlight
airyd(N/2,N/2) = 1.0;  % fudge to correct divide by zero, not-a-number
imshow(airyd);
pause;
airyd = airyd/sum(sum(airyd)); % normalize the matrix

% load an image to analyze; image needs to be
% scaled so that 1 pixel = 1 mm.
plate=imread('plate_grey.jpg');
plate = single(plate)/256.0;
imshow(plate);
pause;

plate = filter2(airyd,plate);  % convolution for diffraction limits
imshow(plate);
pause;

% consider the finite pixel size on the detector
% calculate how large of a block is required
bsize = round(pixel/(focal/dist));
% found the following approach in the block-processing part of the image
% toolkit
f = @(x) mean2(x)*ones(size(x));
plate_pix = blkproc(plate,[bsize bsize],f);
imshow(plate_pix);
