Plane wave phase

function result=Planewavephase(theta1,theta2,Spixel,Xpixel,Ypixel,lambda)
% lambda: 波长
% theta1: 光传播方向在XY平面投影与x轴夹角 x轴方向为0，一二象限为正，三四象限为负 单位：°（范围：-180°——180° ）
% theta2: 光传播方向与z轴夹角 单位：°（范围：0-90°）
% Spixel: 像素大小
% Xpixel: x轴像素
% Ypixel: 轴像素
% 闪耀光栅适用
direction=90-theta2;
% Xtotalsize=Spixel*Xpixel;
% Ytotalsize=Spixel*Ypixel;
Xoriginpoint=ceil((Xpixel+1)/2);
Yoriginpoint=ceil((Ypixel+1)/2);
xaxis=(1-Xoriginpoint:Xpixel-Xoriginpoint)*Spixel;
yaxis=(1-Yoriginpoint:Ypixel-Yoriginpoint)*Spixel;
x=repmat(xaxis,[length(yaxis),1]);
y=repmat(yaxis',[1,length(xaxis)]);
Xvector=sqrt(cosd(direction)^2/(1+tand(theta1)^2));
if tand(theta1)==Inf
Yvector=sqrt(cosd(direction)^2/(1+1/tand(theta1)^2));
elseif tand(theta1)==-Inf
Yvector=sqrt(cosd(direction)^2/(1+1/tand(theta1)^2));
else
Yvector=sqrt(cosd(direction)^2*tand(theta1)^2/(1+tand(theta1)^2));
end
% Zvector=cosd(theta2);
if theta1>=90
    Yvector=-Yvector;
elseif (theta1>=0) && (theta1<90)
    Xvector=-Xvector;
    Yvector=-Yvector;
elseif (theta1>-90) && (theta1<0)
    Xvector=-Xvector;
end
wavefront=(x*Xvector-y*Yvector)/lambda;
phase=mod(wavefront-min(min(wavefront)),1)*2*pi;
result=phase;
% figure;imshow(result,[])

%等价于平面波相位，同样也是傅里叶变换需要滤波的原因
Xdisplacement=(xpixel/2/2-1)*Spixel;Ydisplacement=(ypixel/2/2-1)*Spixel;
pointolite=sqrt((x+Xdisplacement).^2+(y+Ydisplacement).^2)<1e-9;
figure;imshow(pointolite,[]),title('偏移点')
figure;imshow(angle(fftshift(fft2(fftshift(pointolite)))),[]),title('频域相移')