基于Matlab的大气湍流光束传输特性的研究

以下是一个基于Matlab实现大气湍流光束传输特性研究的详细代码及解释。

% 定义参数
N = 512; % 网格点数
L0 = 10; % 外尺度 (m)
l0 = 0.01; % 内尺度 (m)
Cn2 = 1e-14; % 大气折射率结构常数 (m^(-2/3))
k = 2*pi/0.6328e-6; % 波数 (m^(-1))
z = 1000; % 传输距离 (m)
w0 = 0.1; % 束腰半径 (m)

% 生成随机相位屏
[phase_screen] = generate_phase_screen(N, L0, l0, Cn2, k);

% 生成初始光束
[x, y, u0] = generate_initial_beam(N, w0, k);

% 光束在自由空间中的传播
u1 = propagate_beam(u0, k, z);

% 引入大气湍流相位屏
u2 = u1 .* exp(1i * phase_screen);

% 再次传播通过大气湍流
u3 = propagate_beam(u2, k, z);

% 计算光强
I0 = abs(u0).^2;
I1 = abs(u1).^2;
I3 = abs(u3).^2;

% 计算到达角起伏
[theta_x, theta_y] = calculate_arrival_angle(u3, k, z);

% 计算光束扩展
[w_x, w_y] = calculate_beam_width(u3);

% 计算光束漂移
[dx, dy] = calculate_beam_drift(u0, u3);

% 显示结果
figure;
subplot(2,2,1);
imagesc(abs(u0).^2);
title('初始光束光强');
colorbar;

subplot(2,2,2);
imagesc(abs(u1).^2);
title('自由空间传播后光强');
colorbar;

subplot(2,2,3);
imagesc(abs(u3).^2);
title('通过大气湍流后光强');
colorbar;

subplot(2,2,4);
plot(theta_x, theta_y);
title('到达角起伏');
xlabel('x方向到达角 (rad)');
ylabel('y方向到达角 (rad)');

fprintf('光束扩展 (x方向): %.4f m\n', w_x);
fprintf('光束扩展 (y方向): %.4f m\n', w_y);
fprintf('光束漂移 (x方向): %.4f m\n', dx);
fprintf('光束漂移 (y方向): %.4f m\n', dy);

% 生成随机相位屏函数
function [phase_screen] = generate_phase_screen(N, L0, l0, Cn2, k)
    delta = 0.1; % 网格间距 (m)
    fx = (-N/2:N/2-1)/(N*delta);
    [Fx, Fy] = meshgrid(fx, fx);
    f = sqrt(Fx.^2 + Fy.^2);
    f0 = 1/L0;
    fl = 1/l0;
    % Kolmogorov谱模型
    phi = 0.033 * Cn2 * (f.^2 + f0^2).^(-11/6) .* exp(-f.^2/fl^2);
    phi(round(N/2)+1, round(N/2)+1) = 0;
    phase_screen = sqrt(2*pi*k^2*delta^2) * ifft2(sqrt(phi) .* exp(1i*2*pi*rand(N)));
    phase_screen = real(phase_screen);
end

% 生成初始光束函数
function [x, y, u0] = generate_initial_beam(N, w0, k)
    delta = 0.1; % 网格间距 (m)
    x = (-N/2:N/2-1)*delta;
    [X, Y] = meshgrid(x, x);
    u0 = exp(-(X.^2 + Y.^2)/w0^2);
end

% 光束传播函数
function [u1] = propagate_beam(u0, k, z)
    N = size(u0, 1);
    delta = 0.1; % 网格间距 (m)
    fx = (-N/2:N/2-1)/(N*delta);
    [Fx, Fy] = meshgrid(fx, fx);
    H = exp(-1i*pi*delta^2*z*(Fx.^2 + Fy.^2));
    U0 = fft2(u0);
    U1 = U0 .* H;
    u1 = ifft2(U1);
end

% 计算到达角起伏函数
function [theta_x, theta_y] = calculate_arrival_angle(u, k, z)
    N = size(u, 1);
    delta = 0.1; % 网格间距 (m)
    x = (-N/2:N/2-1)*delta;
    [X, Y] = meshgrid(x, x);
    phase = angle(u);
    [dphidx, dphidy] = gradient(phase, delta);
    theta_x = -dphidx/k/z;
    theta_y = -dphidy/k/z;
end

% 计算光束扩展函数
function [w_x, w_y] = calculate_beam_width(u)
    N = size(u, 1);
    delta = 0.1; % 网格间距 (m)
    x = (-N/2:N/2-1)*delta;
    [X, Y] = meshgrid(x, x);
    I = abs(u).^2;
    P = sum(I(:));
    x_c = sum(sum(X .* I))/P;
    y_c = sum(sum(Y .* I))/P;
    w_x = sqrt(2 * sum(sum((X - x_c).^2 .* I))/P);
    w_y = sqrt(2 * sum(sum((Y - y_c).^2 .* I))/P);
end

% 计算光束漂移函数
function [dx, dy] = calculate_beam_drift(u0, u3)
    N = size(u0, 1);
    delta = 0.1; % 网格间距 (m)
    x = (-N/2:N/2-1)*delta;
    [X, Y] = meshgrid(x, x);
    I0 = abs(u0).^2;
    I3 = abs(u3).^2;
    P0 = sum(I0(:));
    P3 = sum(I3(:));
    x_c0 = sum(sum(X .* I0))/P0;
    y_c0 = sum(sum(Y .* I0))/P0;
    x_c3 = sum(sum(X .* I3))/P3;
    y_c3 = sum(sum(Y .* I3))/P3;
    dx = x_c3 - x_c0;
    dy = y_c3 - y_c0;
end    

代码解释：

1. 参数定义：

   • N：网格点数，用于定义模拟的空间分辨率。
   • L0 和 l0：分别是大气湍流的外尺度和内尺度。
   • Cn2：大气折射率结构常数，用于描述大气湍流的强度。
   • k：波数，由波长计算得到。
   • z：光束传输距离。
   • w0：光束的束腰半径。

2. 生成随机相位屏：

   • generate_phase_screen 函数基于Kolmogorov谱模型生成符合大气湍流统计特性的随机相位屏。

3. 生成初始光束：

   • generate_initial_beam 函数生成高斯光束作为初始光束。

4. 光束传播：

   • propagate_beam 函数使用傅里叶变换方法模拟光束在自由空间中的传播。

5. 引入大气湍流相位屏：

   • 将生成的随机相位屏引入到光束中，模拟光束在大气湍流中的传播。

6. 再次传播：

   • 光束通过大气湍流后再次传播，得到最终的光束分布。

7. 计算光束特性：

   • calculate_arrival_angle 函数计算到达角起伏。
   • calculate_beam_width 函数计算光束扩展。
   • calculate_beam_drift 函数计算光束漂移。

8. 显示结果：

   • 使用 subplot 函数显示初始光束、自由空间传播后和通过大气湍流后的光强分布，以及到达角起伏。
   • 打印光束扩展和光束漂移的结果。

通过修改参数，如 Cn2、w0 等，可以分析不同湍流强度和光束参数下大气湍流对光束传输的影响。