光学，四步相移，光学偏振相移仿真

Question

最基础用matlab模拟偏振四步相移仿真，出来的干涉条纹图，不知道到底哪里错了，求指正改错。或者直接是我写的代码不对？

clc;clear all;close all
lambda = 632*10^(-9);        %单位米
% a=double(imread('E:\matab\偏振\pianzhenxiangyi\分辨率板.bmp'));
r = 512;
c = r;
pixelsize = 6.4*10^(-6);
L0 = 5*10^(-3) ;
k = 2*pi/lambda;
% d = 0.005;
x = linspace(-L0/2,L0/2,c);
y = linspace(-L0/2,L0/2,r);
[x,y] = meshgrid(x,y);
sigma1 = 0;
sigma2 = pi/4;
sigma3 = pi/2;
sigma4 = 3*pi/4;
Hp1 = [(cos(sigma1)).^2,sin(sigma1)*cos(sigma1);sin(sigma1)*cos(sigma1),(sin(sigma1)).^2];
Hp2 = [(cos(sigma2)).^2,sin(sigma2)*cos(sigma2);sin(sigma2)*cos(sigma2),(sin(sigma2)).^2];
Hp3 = [(cos(sigma3)).^2,sin(sigma3)*cos(sigma3);sin(sigma3)*cos(sigma3),(sin(sigma3)).^2];
Hp4 = [(cos(sigma4)).^2,sin(sigma4)*cos(sigma4);sin(sigma4)*cos(sigma4),(sin(sigma4)).^2];
%% 物光
alpha2 = pi/2.00;
beita2 = pi/2.02;
Eo = exp(j*k*(x*cos(alpha2)+y*cos(beita2))); %入射平面光
figure,imshow(angle(Eo),[]);title('平面光波包裹相位')
%% 水平偏振
ex = [1;0]; %光矢量沿x方向
Eo1 = cell(512,512); 
for m = 1:512
    for n = 1:512
        Eo1{m,n} = Eo(m,n).*ex;
    end
end
[M,N] = size(Eo1);
Eox = zeros(M,N);
for m = 1:M
    for n = 1:N
        Eox (m,n) = sqrt(Eo1{m,n}(1,1).^2+Eo1{m,n}(2,1).^2);
    end
end
figure,imshow(angle(Eox),[]);title('水平偏振参考光包裹相位')
%% 左旋
QWPL = (1+i)/2*[1,-i;-i,1];
Eo2 = cell(512,512); 
Eo3 = cell(512,512);
Eo4 = cell(512,512);
Eo5 = cell(512,512);
for m = 1:512
    for n = 1:512
        Eo2{m,n} = Hp1*QWPL*Eo1{m,n};
        Eo3{m,n} = Hp2*QWPL*Eo1{m,n};
        Eo4{m,n} = Hp3*QWPL*Eo1{m,n};
        Eo5{m,n} = Hp4*QWPL*Eo1{m,n};
    end
end
[M,N] = size(Eo2);
Eoy2 = zeros(M,N);
Eoy3 = zeros(M,N);
Eoy4 = zeros(M,N);
Eoy5 = zeros(M,N);
for m = 1:M
    for n = 1:N
        Eoy2(m,n) = sqrt(Eo2{m,n}(1,1).^2+Eo2{m,n}(2,1).^2);
        Eoy3(m,n) = sqrt(Eo3{m,n}(1,1).^2+Eo3{m,n}(2,1).^2);
        Eoy4(m,n) = sqrt(Eo4{m,n}(1,1).^2+Eo4{m,n}(2,1).^2);
        Eoy5(m,n) = sqrt(Eo5{m,n}(1,1).^2+Eo5{m,n}(2,1).^2);
    end
end
figure,imshow(angle(Eoy2),[]);title('右旋参考光偏振包裹相位')

%% 参考光
alpha2 = pi/2.005;
beita2 = pi/2.005;
Er = exp(j*k*(x*cos(alpha2)+y*cos(beita2))); %入射平面光
figure,imshow(angle(Er),[]);title('平面光波包裹相位')
%% 竖直偏振
ey = [0;1]; %光矢量沿y方向
Er1 = cell(512,512); 
for m = 1:512
    for n = 1:512
        Er1{m,n} = Eo(m,n).*ey;
    end
end
[M,N] = size(Er1);
Erx = zeros(M,N);
for m = 1:M
    for n = 1:N
        Erx(m,n) = sqrt(Er1{m,n}(1,1).^2+Er1{m,n}(2,1).^2);
    end
end
figure,imshow(angle(Erx),[]);title('竖直偏振参考光包裹相位')
%% 右旋光
QWPR = (1+i)/2*[1,-i;-i,1];
Er2 = cell(512,512); 
Er3 = cell(512,512);
Er4 = cell(512,512);
Er5 = cell(512,512);
for m = 1:512
    for n = 1:512
        Er2{m,n} =Hp1*QWPR*Er1{m,n};
        Er3{m,n} =Hp2*QWPR*Er1{m,n};
        Er4{m,n} =Hp3*QWPR*Er1{m,n};
        Er5{m,n} =Hp4*QWPR*Er1{m,n};
    end
end
[M,N] = size(Er2);
Erx2 = zeros(M,N);
Erx3 = zeros(M,N);
Erx4 = zeros(M,N);
Erx5 = zeros(M,N);
for m = 1:M
    for n = 1:N
        Erx2(m,n) = sqrt(Er2{m,n}(1,1).^2+Er2{m,n}(2,1).^2);
        Erx3(m,n) = sqrt(Er3{m,n}(1,1).^2+Er3{m,n}(2,1).^2);
        Erx4(m,n) = sqrt(Er4{m,n}(1,1).^2+Er4{m,n}(2,1).^2);
        Erx5(m,n) = sqrt(Er5{m,n}(1,1).^2+Er5{m,n}(2,1).^2);
    end
end
figure,imshow(angle(Erx2),[]);title('右旋参考光偏振包裹相位')
%% 干涉
I1 = Eoy2+Erx2;
I2 = Eoy3+Erx3;
I3 = Eoy4+Erx4;
I4 = Eoy5+Erx5;
II1 = I1.*conj(I1);
II2 = I2.*conj(I2);
II3 = I3.*conj(I3);
II4 = I4.*conj(I4);
figure,imshow(I1,[]);title('干涉光强图1')
figure,imshow(I2,[]);title('干涉光强图2')
figure,imshow(I3,[]);title('干涉光强图3')
figure,imshow(I4,[]);title('干涉光强图4')
figure,plot(linspace(1,50,50),I1(1:50,256));title('干涉1')
figure,plot(linspace(1,50,50),I2(1:50,256));title('干涉2')
figure,plot(linspace(1,50,50),I3(1:50,256));title('干涉3')
figure,plot(linspace(1,50,50),I4(1:50,256));title('干涉4')