close all;

D = 60;  		% mm  about 2.3 inches
lambda = 600 * 1e-6;  	%  600 nm converted to mm -- red

N = 5000;
r = (D/2)*((1:N)'-0.5)/N;
M = 640;
rho = (D/2)*(0:2*M)'/M;
dr = r(2)-r(1);

kmax = 400;
EE_hole = zeros(kmax, length(rho));

figure(1);
hold off;
for k = 1:kmax
  k
  if k<kmax/2
    n = k/20;
    z = D^2/(lambda * n * 8)
    z0 = z;
    k0 = k;
  else
    z = z0*(kmax-k+4)/(kmax-k0+4)
  end

  E = (2*pi)/(i * lambda * z) * exp(i * pi * rho.^2 / (lambda*z)) ...
      .* (besselj(0,(2*pi/(lambda*z))*rho*r') * (exp(i * pi * r.^2 / (lambda*z)) .* r)) * dr;
  zz(k) = z;
  EE_hole(k,:) = E';
  rho2 = [-rho(end:-1:1); rho];
  E2 = [E(end:-1:1); E];
  plot(rho2, abs(E2).^2); ylim([0 4]);
  drawnow;
end

fig1 = figure(1);
winsize = get(fig1, 'Position');
winsize(1:2) = [0 0];
%winsize(3:4) = winsize(3:4);
numframes = kmax;
A = moviein(numframes, fig1, winsize);
set(fig1, 'nextplot', 'replacechildren');
j = 1;
for k = kmax:-1:1
  k
  zz(k)/D
  rho2 = [-rho(end:-1:2); rho];
  E2 = [EE_hole(k,end:-1:2)'; EE_hole(k,1:end)'];
  plot(rho2, abs(E2).^2); ylim([0 4]);
  text(20,3.5,sprintf('dist = %5.2f D',zz(k)/D));
  drawnow;
  A(:,j) = getframe(fig1, winsize);
  j = j+1;
  pause(0.05);
end
frame_rate = 10;
movie(fig1,A,1,frame_rate,winsize);
movie2avi(A, 'PinHole4.avi','fps',frame_rate);


fig2 = figure(2);
winsize = get(fig2, 'Position');
winsize(1:2) = [0 0];
%winsize(3:4) = winsize(3:4);
numframes = kmax;
A = moviein(numframes, fig2, winsize);
set(fig2, 'nextplot', 'replacechildren');
j = 1;
n = (length(E2)-1)/2;
n = floor(n/sqrt(2));
x = ones(2*n+1,1)*(-n:n);
y = x';
r = uint16(1+round(sqrt(x.^2 + y.^2)));
for k = kmax:-1:1
  k
  E1D = abs(EE_hole(k,:).^2);
  E2D = E1D(r);
  E2D_RGB(:,:,1) = E2D;
  E2D_RGB(:,:,2) = E2D;
  E2D_RGB(:,:,3) = E2D;
  imagesc(E2D_RGB,[0 3]); colormap gray;
  axis equal; axis tight; axis off;
  drawnow;
  A(:,j) = getframe(fig2, winsize);
  j = j+1;
  pause(0.05);
end
frame_rate = 10;
%movie(fig2,A,1,frame_rate,winsize);
movie2avi(A, 'PinHole2D4.avi','fps',frame_rate);


fig3 = figure(3);
winsize = get(fig3, 'Position');
winsize(1:2) = [0 0];
%winsize(3:4) = winsize(3:4);
numframes = kmax;
A = moviein(numframes, fig3, winsize);
set(fig3, 'nextplot', 'replacechildren');
j = 1;
n = (length(E2)-1)/2;
n = floor(n/sqrt(2));
x = ones(2*n+1,1)*(-n:n);
y = x';
r = uint16(1+round(sqrt(x.^2 + y.^2)));
for k = kmax:-1:1
  k
  E1D = abs(EE_hole(k,:).^2);
  E2D = E1D(r);
  E2D_RGB(:,:,1) = poissrnd(500*E2D);
  imagesc(E2D_RGB(:,:,1),[0 1500]); colormap gray;
  axis equal; axis tight; axis off;
  drawnow;
  A(:,j) = getframe(fig3, winsize);
  j = j+1;
  pause(0.05);
end
frame_rate = 10;
%movie(fig3,A,1,frame_rate,winsize);
movie2avi(A, 'PinHole2Dpoisson4.avi','fps',frame_rate);