rand2

% open input image
D = 256;
circ = fftprep(imread('circ16.tif'),D);
[rp, vref] = psf(circ);
mov = avifile('random.avi','fps',5);
strehl(1)=1.0;
sig(1)=0.0;
mesh(rp);
f=getframe;
mov = addframe(mov,f);

for k=1:25
   randf = (rand(D)-0.5)*sqrt(12.0);
   sig(k+1) = 0.02*k;
	phase = exp(-i*2*pi*sig(k+1)*randf);
	[rp, strehl(k+1)] = psf(circ.*phase,vref);
   mesh(rp);
   f=getframe;
   mov = addframe(mov,f);
 end
mov = close(mov);

r = linspace(0,max(sig),500);
y = exp(-(2*pi*r).^2);
plot(sig,strehl,'o',r,y,'k');
xlabel('\sigma');
ylabel('strehl ratio');