clear
%%
%定义光源屏坐标系
a=100;
b=0.2;
X1 = linspace(-b,b,a);
Y1 = linspace(-b,b,a);
[x,y] = meshgrid(X1,Y1);
%定义第一次传播接收屏
X2 = linspace(-b,b,a);
Y2 = linspace(-b,b,a);
[x2,y2] = meshgrid(X2,Y2);
%定义第一次衍射接收屏
X3 = linspace(-b,b,a);
Y3 = linspace(-b,b,a);
[x3,y3] = meshgrid(X3,Y3);
%定义第二次传播接收屏
X4 = linspace(-b,b,a);
Y4 = linspace(-b,b,a);
[x4,y4] = meshgrid(X4,Y4);
%定义第二次衍射接收屏
X5 = linspace(-b,b,a);
Y5 = linspace(-b,b,a);
[x5,y5] = meshgrid(X5,Y5);
%%
%定义涡旋光
L1=1;L2=-1;L3=0;
sigma=0.1;
R=x.^2+y.^2;
r=sqrt(R);
% l=1
OAM1=(r.^abs(L1)).*exp(-R/(sigma^2)).*exp(L1*1i.*atan(y./x));
COAM1=conj(OAM1);
GQ1=OAM1.*COAM1;
% l=-1
OAM2=(r.^abs(L2)).*exp(-R/(sigma^2)).*exp(L2*1i.*atan(y./x));
COAM2=conj(OAM2);
GQ2=OAM2.*COAM2;
% l=0
OAM3=(r.^abs(L3)).*exp(-R/(sigma^2)).*exp(L3*1i.*atan(y./x));
COAM3=conj(OAM3);
GQ3=OAM3.*COAM3;
% 涡旋光叠加
OAM=OAM1+OAM2+OAM3;
COAM=conj(OAM);
GQ=OAM.*COAM;
% 测试光强
% figure(1);
% subplot(221)
% plot3(x,y,GQ);
% title('三种模式叠加光强');
% subplot(222)
% plot3(x,y,GQ1,'r');
% title('l=1');
% subplot(223)
% plot3(x,y,GQ2,'b');
% title('l=-1');
% subplot(224)
% plot3(x,y,GQ3,'y');
% title('l=0');
% grid on;
%%
% 第一次SLM调制
nameda=1.55e-6; k=2*pi/nameda;
q=5e-4; p=100; fl=1;
Phi1=(k*q/fl)*(y.*atan(y./x)-x.*log(r/p)+x);
E1=exp(1i.*Phi1);
TIAOZHI1=OAM.*E1;
%%
% 第一次传播一个fl的距离
First=zeros(a,a);
z=fl;
for m=1:a
for n=1:a
xz1=x2(m); yz1=y2(n);
Rz1=(x-xz1).^2+(y-yz1).^2+z^2;
rz1=sqrt(Rz1);
First(m,n)=z/(1i*nameda)*sum(sum(TIAOZHI1.*exp(1i*k*rz1)./Rz1));
end
end
% GQ4=First.*conj(First);
% figure(2);
% contourf(GQ4);
%%
% 第一次傅里叶变换
Second=zeros(a,a);
R2=x2.^2+y2.^2;
for f=1:a
for h=1:a
xz2=x3(f); yz2=y3(h);
Rz2=(x2-xz2).^2+(y2-yz2).^2+z^2;
rz2=sqrt(Rz2);
Second(f,h)=z/(1i*nameda)*sum(sum(First.*exp(-1i*k*R2/(2*fl)).*exp(1i*k*rz2)./Rz2));
end
end
% GQ5=Second.*conj(Second);
% figure(3);
% contourf(GQ5);
%%
% 第二次SLM调制
Phi2=(k*q*p/fl)*exp(-x3/q).*cos(y3/q);
E2=exp(1i.*Phi2);
TIAOZHI2=Second.*E2;
%%
% 第二次传播一个fl的距离
Third=zeros(a,a);
z=fl;
for aerfa=1:a
for beta=1:a
xz3=x4(aerfa); yz3=y4(beta);
Rz3=(x3-xz3).^2+(y3-yz3).^2+z^2;
rz3=sqrt(Rz3);
Third(aerfa,beta)=z/(1i*nameda)*sum(sum(TIAOZHI2.*exp(1i*k*rz3)./Rz3));
end
end
% GQ6=Third.*conj(Third);
% figure(4);
% contourf(GQ6);
%%
% 第二次傅里叶变换
Forth=zeros(a,a);
R4=x4.^2+y4.^2;
for zhao=1:a
for qian=1:a
xz4=x5(zhao); yz4=y5(qian);
Rz4=(x4-xz4).^2+(y4-yz4).^2+z^2;
rz4=sqrt(Rz4);
Forth(zhao,qian)=z/(1i*nameda)*sum(sum(Third.*exp(-1i*k*R4/(2*fl)).*exp(1i*k*rz4)./Rz4));
end
end
GQ7=Forth.*conj(Forth);
figure(5);
contourf(GQ7);
