根据下面图中的步骤用MATLAB代码生成相位屏，我自己写了一下但是运行了好几天了也没有运行完，我觉得可能是我写的代码有问题，有会的可以帮我看一下吗，特别是C3这一步我不知道怎么写

%% 生成相位屏
eta=5e-3; %Kolmogorov微尺度长度
epsilon=1e-4; %单位质量流体的湍流动能的耗散率
XT=1e-7; %平均平方温度耗散率
omega=-0.3; %温度和盐度对折射率谱的贡献比率 
At=1.863e-2;
As=1.9e-4;
Ats=9.41e-3;
n=1e3; %采样点数
delta=1e-3; %采样点之间的距离
k=2*pi/lambda; %光波数
deltaz=2.5; %相位屏之间的间隔
Lx=1; %相位屏的边长
Ly=1;
kx=(-n/2:(n/2-1))*2*pi/delta; %空间频率在x方向和y方向上的分量
ky=(-n/2:(n/2-1))*2*pi/delta;
s=8.284*(sqrt(kx.^2+ky.^2)*eta).^(4/3)+12.978*(sqrt(kx.^2+ky.^2)*eta).^2;
k1=sqrt(kx.^2+ky.^2);

deltakx=2*pi/(n*delta);
deltaky=2*pi/(n*delta);
C=sqrt(deltakx*deltaky);

psi=0.388e-8*epsilon^(-1/3).*(k1).^(-11/3).*(1+2.35.*(k1.*eta).^(2/3)).*(XT/omega^2*(omega^2*exp(-At*s)+exp(-As*s)-2*omega*exp(-Ats*s)));%选定海洋功率谱

F=2*pi*k^2*deltaz*psi; %确定相位屏谱

h=randn(n,n)+1i*randn(n,n); %生成高斯随机复矩阵
sigma2=(2*pi/(n*delta))^2*F; %随机相位屏频谱的方差
% phase_screen=fft(h.*sigma2); %生成相位屏

syms m1 n1 m2 n2 integer
% m1=m2==1000;
% n1=n2==1000;
x=m1*delta;
y=n1*delta;
phase_screen=0;
for kx=(-n/2:(n/2-1))*2*pi/delta
    for ky=(-n/2:(n/2-1))*2*pi/delta
       phase_screen=phase_screen+C.*h.*sqrt(sigma2).*exp(1i.*(2*pi*m1*m2/1000+2*pi*n1*n2/1000));       
    end
end