clear;close all;clc;
g=ones(1080,1080);
[m,n]=size(g);
Lx=1080*10*10^-6;
Ly=1080*10*10^-6;
lamda=6.32*10^(-7);
k=2*pi/lamda;
x0=linspace(-Lx/2,Lx/2,m);
y0=linspace(-Ly/2,Ly/2,n);
[x,y]=meshgrid(x0,y0);

[theta,r]=cart2pol(x,y);

T0=0.4*10^-3;n1=0.02;

T=cos(2*pi*x./(T0+n1.*x));
%以下代码为修改
for i=1:1080
for j=1:1080
if T(i,j)>0
T(i,j)=1;
else
T(i,j)=0;
end
end
end

figure;
imshow(abs(T),[]);
l=1;w=1.*10^-3;p=0;z=300.*10^-3;
LG=lag_Gaussian(lamda,w,2,2,z,0.*theta,r);
figure;
imshow(abs(LG),[]);
z0=z;
N=1080;
U0=LG;
U1=abs(T);
L0=1080*10*10^-6;
n=1:N;
x=-L0/2+L0/N*(n-1);
y=x;
[yy,xx]=meshgrid(y,x);
Fresnel=exp(1i*k./2./z0*((xx).^2+(yy).^2));
f2=U0.*U1*Fresnel;

Uf=fft2(f2,N,N);
Uf=fftshift(Uf);

L=lamda*z0*N/L0;
x=-L/2+L/N*(n-1);
y=x;
[yy,xx]=meshgrid(y,x);
phase=exp(1i*k*z0)/(1i*lamda*z0)*exp(1i*k/2/z0*(xx.^2+yy.^2));
Uf=Uf.*phase;
g=L0/N;
Uf=Uf*g*g;
If=Uf;
figure;imshow(abs(If),[]);求问此程序哪里有问题