clc;
clear;
close;
%% Input parameters
% screen properties
D1 = 2e-1; % diameter of the source aperture [m]
D2 = 2e-1; % diameter of the observation aperture [m]
N = 512; % sampling numbers
deltan = D2/N; % observation-plane spacing
delta1 = D1/N; % source-plane spacing
% coordinates
[x1, y1] = meshgrid((-N/2 : N/2-1) * delta1); % coordinates
[theta1, r1] = cart2pol(x1, y1);
wvl = 0.532e-6; % optical wavelength [m] % beam parameters
k = 2*pi / wvl; % optical wavenumber [rad/m]
% Beam parameter
w0 = 0.001; % beam width (m)
l = 3; % OAM modes
p = 2;
nscr = 100; % number of phase screens
alpha = 11/3; % power specturm index
L0 = 5; % outer scale [m]
l0 = 0.005; % inner scale [m]
% Spherical wave parameters
R = .0001; % radius of curvature [m] (positive for diverging, negative for converging)
% input fields (x/y polarizations)
normFactor = sqrt(2*factorial(p) / (pi * factorial(abs(l) + p)))/w0;
u0 = normFactor .* exp(-1.*r1.^2./w0.^2) .* (sqrt(2).*r1./w0).^(l) .* Lagureer(p,l, 2.*r1.^2./w0.^2) .* exp(1i*l*theta1);
u2=exp(-1i * k*r1^2/R);
% Add spherical wave component
%spherical_phase = exp(-1i * k * r1.^2 / (2 * R));
u0 = u0+u2;
Cn2 = 1e-15;
dz = 10; % distance interval it should be greater than L0
Tz = dz * nscr; % propagation distance [m]
r0 = (0.423 .* k.^2 .* Cn2 * Tz)^(-3/5); % fried parameter
n = nscr; % planes for masks
% weighting for r0 across screens (same as your script)
sig = 0.01;
r0t1 = r0 * nscr;
r0t2 = r0 * nscr * nscr;
sigt = sig * nscr;
an = zeros(1, nscr);
parfor ii = 1:nscr
an(ii) = 1/(exp(-((ii-1).^2)/(sigt^2))*r0t1 + 1) * r0t1;
end
r0scrn = (r0t2/sum(an)) * an;
d1 = 10e-3;
%% simulate vacuum propagation free space without turbulence propagation Vector
z = (1 : n-1) * Tz / (n-1);
sg = exp(-(x1/(0.47*N*d1)).^16) .* exp(-(y1/(0.47*N*d1)).^16);
t = repmat(sg, [1 1 n]);
[x1, y1, Uvac] = ang_spec_multi_prop(u0, wvl, delta1, deltan, z, t);
% Intensity plot
subplot(2,2,1);
surf(x1, y1, abs(Uvac).^2);
view (2); shading interp; colormap("hot"); grid off; axis off;
title('Vacuum - Intensity');
set(gca,'LooseInset',get(gca,'TightInset'));
% Phase plot
%subplot(2,2,2);
%surf(x1, y1, angle(Uvac));
%view (2); shading interp; colormap("hot"); grid off; axis off;
%title('Vacuum - Phase');
%set(gca,'LooseInset',get(gca,'TightInset'));
%% -----------------------Step three: turbulence propagation-------------------------------------
zt = [0 z]; % propagation plane locations
ad = zt / zt(n);
delta = (1-ad) * delta1 + ad * deltan;
% loading random phase screen
phz = zeros(N, N, n); % initialize array for phase screens
sg_repeated = repmat(sg, [1 1 n]); % Repeat sg to match the dimensions
% loop over screens
parfor idxscr = 1:n
[phz_lo, phz_hi] = ft_sh_phase_screen(r0scrn(idxscr), N, delta(idxscr), alpha, L0, l0);
phz(:,:,idxscr) = phz_lo + phz_hi;
end
% ------------------------- propagate ---------------------------
[x1, y1, Uout] = ang_spec_multi_prop(u0, wvl, delta1, deltan, z, sg_repeated .* exp(1i * phz));
% Intensity plot
subplot(2,2,3);
surf(x1, y1, abs(Uout).^2);
view (2); shading interp; colormap("hot"); grid off; axis off;
title('Turbulence - Intensity');
set(gca,'LooseInset',get(gca,'TightInset'));
% Phase plot
subplot(2,2,4);
surf(x1, y1, angle(Uout));
view (2); shading interp; colormap("hot"); grid off; axis off;
title('Turbulence - Phase');
set(gca,'LooseInset',get(gca,'TightInset'));
%% -----------------------------------Basic Functions-------------------------------------
function [xn, yn, Uout] = ang_spec_multi_prop(Uin, wvl, delta1, deltan, z, t)
N = size(Uin, 1); % number of grid points
[nx, ny] = meshgrid((-N/2 : 1 : N/2 - 1));
k = 2*pi/wvl; % optical wavevector
% super-Gaussian absorbing boundary
nsq = nx.^2 + ny.^2;
w = 0.47*N;
sg = exp(-nsq.^8/w^16); clear('nsq', 'w');
z = [0, z]; % propagation plane locations
n = length(z);
Delta_z = z(2:n) - z(1:n-1); % propagation distances
ad = z / z(n); % grid spacings
delta = (1-ad) * delta1 + ad * deltan;
m = delta(2:n) ./ delta(1:n-1);
x1 = nx * delta(1);
y1 = ny * delta(1);
r1sq = x1.^2 + y1.^2;
Q1 = exp(1i*k/2*(1-m(1))/Delta_z(1)*r1sq);
Uin = Uin .* Q1 .* t(:,:,1);
for idx = 1 : n-1
deltaf = 1 / (N*delta(idx)); % spatial frequencies (of i^th plane)
fX = nx * deltaf;
fY = ny * deltaf;
fsq = fX.^2 + fY.^2;
Z = Delta_z(idx); % propagation distance
Q2 = exp(-1i*pi^2*2*Z/m(idx)/k*fsq); % quadratic phase factor
Uin = sg .* t(:,:,idx+1) .* ift2(Q2 .* ft2(Uin / m(idx), delta(idx)), deltaf); % compute the propagated field
end
xn = nx * delta(n); yn = ny * delta(n); % observation-plane coordinates
rnsq = xn.^2 + yn.^2;
Q3 = exp(1i*k/2*(m(n-1)-1)/(m(n-1)*Z)*rnsq);
Uout = Q3 .* Uin;
end
function [phz_lo, phz_hi] = ft_sh_phase_screen(r0, N, delta, alpha, L0, l0)
persistent last_data;
if isempty(last_data)
last_data = struct();
end
D = N*delta;
phz_hi = ft_phase_screen(r0, N, delta, alpha, L0, l0);
[x,y] = meshgrid((-N/2 : N/2-1) * delta);
phz_lo = zeros(size(phz_hi));
for p = 1:3
del_f = 1 / (3^p*D);
fx = (-1 : 1) * del_f;
[fx, fy] = meshgrid(fx);
fm = 2.68934/l0/(2*pi);
f0 = 4*pi/L0/(2*pi);
PSD_phi = 0.023.*r0.^(2-alpha).*(fx.^2 +fy.^2+f0^2).^(-alpha/2).*exp(-((fx.^2 +fy.^2)./fm.^2)).*(1-0.061.*((fx.^2 +fy.^2).^(1/2)./fm)+2.836.*((fx.^2 +fy.^2).^(1/2)./fm).^(3-alpha./2));
PSD_phi(2,2) = 0;
cn = (randn(3) + 1i*randn(3)).*sqrt(PSD_phi)*del_f;
SH = zeros(N);
for ii = 1:9
SH = SH + cn(ii)*exp(1i*2*pi*(fx(ii)*x+fy(ii)*y));
end
phz_lo = phz_lo + SH;
end
phz_lo = real(phz_lo) - mean(real(phz_lo(:)));
end
function phz = ft_phase_screen(r0, N, delta, alpha, L0, l0)
persistent last_data;
if isempty(last_data)
last_data = struct();
end
del_f = 1/(N*delta);
fx = (-N/2 : N/2-1) * del_f;
[fx, fy] = meshgrid(fx);
fm = 2.68934/l0/(2*pi);
f0 = 4*pi/L0/(2*pi);
PSD_phi = 0.023.*r0.^(2-alpha).*(fx.^2 +fy.^2+f0^2).^(-alpha/2).*exp(-((fx.^2 +fy.^2)./fm.^2)).*(1-0.061.*((fx.^2 +fy.^2).^(1/2)./fm)+2.836.*((fx.^2 +fy.^2).^(1/2)./fm).^(3-alpha./2));
PSD_phi(N/2+1,N/2+1) = 0;
cn = (randn(N) + 1i*randn(N)) .* sqrt(PSD_phi)*del_f;
phz = real(ift2(cn, 1));
end
function G = ft2(g, delta)
G = fftshift(fft2(fftshift(g))) * delta^2;
end
function g = ift2(G, delta_f)
N = size(G, 1);
g = ifftshift(ifft2(ifftshift(G))) * (N * delta_f)^2;
end
function y = Lagureer(p, l, x)
fact = @(x)arrayfun(@factorial, x);
n = p;
k = l;
m = 0:n;
a = factorial(n+k)*ones(1,length(m));
b = fact(n-m);
c = fact(k+m);
d = fact(m);
e = (-1).^m;
y = zeros(size(x));
for s = 1:n+1
y = y + a(s) ./ b(s) ./ c(s) ./ d(s) .* e(s) .* x.^m(s);
end
end
