clc;
clear;
close all;
warning off;
addpath(genpath(pwd));
rng('default')
xmax = 60; % 设置x的最大范围
k = 1; % 波数
w = 1; % 频率（但此变量在后续代码中未被使用）
L = 10; % 镜子的位置
r0 = 2; % 高斯光束的半径
N = 2^10 + 1; % 设置x的采样点数
x = linspace(-xmax/2, xmax/2, N); % 生成x的线性空间
% 理论部分
Rt = r0 * (1 + 4 * L^2/k^2/r0^4)^0.5; % 计算在镜子位置的理论光束半径
psi_t = r0/Rt * exp(-2 * x.^2./Rt^2); % 计算理论波函数
% 模拟部分
psi_0 = exp(-x.^2/r0^2); % 初始波函数
dx = xmax /(N - 1); % 计算x的间隔
psi_0_fft = fftshift(fft(psi_0)); % 对初始波函数进行FFT变换并移位
fmax = 1/2/dx; % 计算最大频率
f = 2 * pi .* linspace(-fmax, fmax, N );% 生成频率的线性空间
sol_L_fft = psi_0_fft .* exp(-1i * f.^2/2/k * L); % 计算在镜子位置的FFT波函数
sol = ifft(sol_L_fft); % 对FFT波函数进行反FFT变换
sol_norm = abs(sol).^2 ; % 计算波函数的模平方
% 绘图：理论、模拟和初始波函数
figure(1); plot(x, psi_t, '-r', x, sol_norm, '-b', x, psi_0, '-k');
legend('Theoretical solution on z = L', 'Simulation on z = L', 'Wave on z = 0')
xlabel('x'); ylabel('Square of the amplitude');grid('on')
title('Wave on z = L');
% 另一种理论计算
rt = r0* (1 + 2i * L/k/r0^2)^0.5; % 计算复数半径
psi_L = r0/rt * exp(-x.^2/rt^2); % 计算在镜子位置的波函数
psi_L_conj = conj(psi_L); % 计算波函数的共轭
% 使用有紧支持的χ_M进行计算
rM = 2:4:22; % 设置不同的rM值
sol_0 = zeros(length(rM), length(x)); % 初始化解矩阵
for i = 1 : length(rM)
chi_M = (1 - (x./2/rM(i)).^2).^2; % 计算χ_M
chi_M(x <-2*rM(i)) = 0 ; % 设置x小于-2rM时的χ_M值为0
chi_M(x > 2*rM(i)) = 0 ; % 设置x大于2rM时的χ_M值为0
term1_fft = fftshift(fft(psi_L_conj .* chi_M)); % 计算FFT变换
sol_0_fft = term1_fft.* exp(1i * f.^2/2/k * L); % 计算时间反转后的FFT波函数
sol_0(i, :) = abs(ifft(sol_0_fft)); % 计算反FFT变换的模
end
% 绘图：不同的rM值对应的波函数
figure(2); plot(x, sol_0, 'LineWidth', 2);
legend('rM = 2', 'rM = 6', 'rM = 10', 'rM = 14', 'rM = 18', 'rM = 22')
title('Wave on z = 0, time reversal')
xlabel('x'); ylim([min(min(sol_0)), max(max(sol_0))]); ylabel('Amplitude'); grid('off')
% 使用高斯χ_M进行计算
rM = 2;
chi_M = exp(- (x.^2)/rM^2); % 计算高斯χ_M
term1_fft = fftshift(fft(psi_L_conj .* chi_M)); % 计算FFT变换
sol_0_fft = term1_fft.* exp(1i * f.^2/2/k * L); % 计算时间反转后的FFT波函数
sol_0 = abs(ifft(sol_0_fft)); % 计算反FFT变换的模
atr = (1 -4i * L/k/r0^2 - 4 * L^2/k^2/r0^2/rM^2 - 2i*L/k/rM^2)^0.5;
rtr_square = (1/rM^2 + 1/(r0^2 - 2i * L/k))^(-1) - 2i * L/k;
psi_tr_0 = abs(1/atr * exp(- x.^2/rtr_square)); % 计算理论时间反转波函数
% 绘图：模拟和理论时间反转波函数
figure(3); plot(x, sol_0, '-r', x, psi_tr_0, '-b');
legend('Simulation', 'Theoretical solution')
xlabel('x'); ylabel('Amplitude'); grid('on')
title('Wave on z = 0, Gaussian time-reversal mirror')
%%
h = 1; % 纵向（沿z轴）步长
zc = 1; % 用于计算高斯过程的参数
xc = 4; % 高斯过程的另一个参数
sigma = 1; % 高斯过程的标准差
% 生成x的空间分布，并计算频率分布
x = linspace(-xmax/2, xmax/2, N)'; % x的线性空间分布（注意：xmax和N需要在此前定义）
dx = xmax / (N - 1); % x的间隔
fmax = 1 / (2 * dx); % 最大频率
f = 2 * pi .* linspace(-fmax, fmax, N ); % 频率的线性空间分布
u0 = exp(-x.^2/r0^2); % 初始波包形状（注意：r0需要在此前定义）
rM = 2; % 高斯镜的半径参数
mirror = exp(-x.^2 / rM^2); % 高斯镜的形状
nb_MC = 300; % 要平均的波包数量
nb_GP = round(L / zc); % 高斯过程的数量（注意：L需要在此前定义）
% 初始化存储波包的矩阵
U_L = zeros(nb_MC, N); % 传输后的波包形状
U_0_rand = zeros(nb_MC, N); % 随机介质中重聚焦的波包形状
U_0_homo = zeros(nb_MC, N); % 均匀介质中重聚焦的波包形状
% 进行多次模拟以获取平均结果
for i = 1:nb_MC
% 在随机介质中正向传播
GP_seq = sample_GP(x, sigma, xc, nb_GP); % 生成高斯过程序列（sample_GP函数需要另外定义）
U_L(i,:) = split_step_fourier_method(0, 1, round(L/h) - 1, u0, h, k, f, GP_seq); % 使用分步傅里叶方法模拟波包传播（split_step_fourier_method函数需要另外定义，注意：k也需要在此前定义）
% 在相同的随机介质中反向传播
u = conj(U_L(i,:)') .* mirror; % 取共轭并与镜函数相乘
U_0_rand(i,:) = split_step_fourier_method(round(L/h) - 1, -1, 0, u, -h, k, f, GP_seq); % 反向传播
% 在均匀介质中反向传播
U_0_homo(i,:) = split_step_fourier_method(round(L/h) - 1, -1, 0, u, h, k, f); % 均匀介质中无需GP_seq
end
% 对所有运行结果取平均
U_L = mean(U_L, 1)';
U_0_rand = mean(U_0_rand, 1)';
U_0_homo = mean(U_0_homo, 1)';
% 计算理论上的传输波波包形状
rt = r0*sqrt(1 + 2*1i*L/k/r0^2);
gamma0 = sigma^2*zc;
mean_wave_L = r0/rt * exp(-x.^2/rt^2) * exp(-gamma0*w^2*L/8); % 注意：w需要在此前定义
% 绘图：比较理论和模拟的传输波波包形状
figure(4); plot(x, abs(mean_wave_L), x, abs(U_L))
xlabel('x'); ylabel('|\bf{E}[\phi_t(x)]|');
legend('theoretical', 'empirical')
title('Mean transmitted wave profile in random medium')
%% Part 2 : mean refocused wave profile in z=0 (random medium)
% 计算理论上的重聚焦波波包形状（随机介质）
atr = sqrt(1 + 4*L^2/(k*r0*rM)^2 + 2*1i*L/k/rM^2);
rtr_square = ((1/rM^2+1/(r0^2-2*1i*L/k))^-1 + 2*1i*L/k);
gamma2 = 2*sigma^2*zc/xc^2;
ra_square = 48/L/gamma2/w^2 * 1;
mean_wave_0_rand = 1/atr * exp(-x.^2 /rtr_square) .* exp(-x.^2 /ra_square);
% 绘图：比较理论和模拟的重聚焦波波包形状（随机介质）
figure(5); plot(x, 0.28*abs(mean_wave_0_rand), x, 0.28*abs(U_0_rand));
xlabel('x'); ylabel('|\bf{E}[\phi_t^{tr}(x)]|');
legend('theoretical', 'empirical')
title('Mean refocused wave profile in random medium')
% 计算理论上的重聚焦波波包形状（均匀介质）
mean_wave_0_homo = 1/atr * exp(-x.^2/rtr_square) .* exp(-gamma0*w^2*L/8);
% 绘图：比较理论和模拟的重聚焦波波包形状（均匀介质）
figure(6); plot(x, abs(mean_wave_0_homo), x, abs(U_0_homo));
xlabel('x'); ylabel('|\bf{E}[\phi_t^{tr}(x)]|');
legend('theoretical', 'empirical')
title('Mean refocused wave profile in homogeneous medium')
% 以下为时域波形分析
w0 = 1; % 中心频率
B = 0.75; % 带宽
nb_discr = 20; % 频率离散点数
omega_discr = linspace(w0-B,w0+B,nb_discr); % 频率的线性空间分布
nb_GP = round(L / zc); % 重新计算高斯过程的数量
U_0_w_samples = zeros(5, N); % 初始化存储时域波形的矩阵
% 进行多次模拟以获取时域波形
for s = 1:5
U_L_w = zeros(nb_discr, N); % 传输后的波包形状（频率相关）
U_0_w = zeros(nb_discr, N); % 重聚焦的波包形状（频率相关）
GP_seq = sample_GP(x, sigma, xc, nb_GP); % 生成高斯过程序列
for wi = 1:nb_discr
%omega = omega_discr(w); % 这行似乎有误，应该是omega = omega_discr(wi);
k = omega_discr(wi); % 使用当前频率计算波数
U_L_w(wi,
