
# simulace zobrazovani v turbulentni atmosfere

function cviceni10c

close all
graphics_toolkit("gnuplot")

# silici turbulence pro stejnou aperturu


maxD=5;
N=500;
delta=maxD/N;
[x,y]=meshgrid((-N/2:N/2-1)*delta);

r0=0.15;
phase=kolmogorov(r0,N,delta);
D=0.3;
apert=sqrt(x.*x+y.*y)<D/2;
intens=abs(ft2(apert.*exp(I*phase),delta)).^2;
figure(1)
pcolor(intens);
shading('interp')
axis('image')
axis('off')
colormap('hot')
caxis([0,max(max(intens))])




r0=0.075;
phase=kolmogorov(r0,N,delta);
D=0.3;
apert=sqrt(x.*x+y.*y)<D/2;
intens=abs(ft2(apert.*exp(I*phase),delta)).^2;
figure(2)
pcolor(intens);
shading('interp')
axis('image')
axis('off')
colormap('hot')
caxis([0,max(max(intens))])


r0=0.0375;
phase=kolmogorov(r0,N,delta);
D=0.3;
apert=sqrt(x.*x+y.*y)<D/2;
intens=abs(ft2(apert.*exp(I*phase),delta)).^2;
figure(3)
pcolor(intens);
shading('interp')
axis('image')
axis('off')
colormap('hot')
caxis([0,max(max(intens))])



# zvetsujici se apertura pro stejnou turbulenci

r0=0.08;
phase=kolmogorov(r0,N,delta);
D=0.08;
apert=sqrt(x.*x+y.*y)<D/2;
intens=abs(ft2(apert.*exp(I*phase),delta)).^2;
figure(1)
pcolor(intens);
shading('interp')
axis('image')
axis('off')
colormap('hot')
caxis([0,max(max(intens))])


r0=0.08;
phase=kolmogorov(r0,N,delta);
D=0.4;
apert=sqrt(x.*x+y.*y)<D/2;
intens=abs(ft2(apert.*exp(I*phase),delta)).^2;
figure(2)
pcolor(intens);
shading('interp')
axis('image')
axis('off')
colormap('hot')
caxis([0,max(max(intens))/3])


r0=0.08;
phase=kolmogorov(r0,N,delta);
D=0.8;
apert=sqrt(x.*x+y.*y)<D/2;
intens=abs(ft2(apert.*exp(I*phase),delta)).^2;
figure(3)
pcolor(intens);
shading('interp')
axis('image')
axis('off')
colormap('hot')
caxis([0,max(max(intens))/5])



# video

close all

figure(1)

r0=0.08;
D=0.4;
apert=sqrt(x.*x+y.*y)<D/2;
for opak=1:100 
 opak
 phase=kolmogorov(r0,N,delta);
 intens=abs(ft2(apert.*exp(I*phase),delta)).^2;
 imshow(intens,[0,max(max(intens))/3]);
 colormap("hot")
 pause(0.01)
endfor


endfunction



function phz=kolmogorov(r0, N, delta)
 del_f = 1/(N*delta);
 fx = (-N/2 : N/2-1) * del_f;
 [fx fy] = meshgrid(fx);    % pravouhla sit
 [th f] = cart2pol(fx, fy); % polarni sit
 PSD_phi = 0.023*r0^(-5/3) ./ f.^(11/3);
 PSD_phi(N/2+1,N/2+1) = 0;  % odstrani Inf pro f=0
 cn = (randn(N) + i*randn(N)) .* sqrt(PSD_phi)*del_f;
 phz = real(ift2(cn, 1));
endfunction


function g = ft2(G, delta)
 N = size(G, 1);
 g = fftshift(fft2(ifftshift(G))) * delta^2;
endfunction

function g = ift2(G, delta_f)
 N = size(G, 1);
 g = fftshift(ifft2(ifftshift(G))) * (N * delta_f)^2;
endfunction


