%求解N束線偏振光束干涉圖樣。
clear
lamda=0.6328;
N=4;
phi=pi*[0;39/180;39/180;39/180];%極角：z軸到波矢量的夾角
theta=pi*[0;1;1/3;-1/3];%方向角從x軸到波矢量的夾角

%ex=[sin(phi+pi/2).*cos(theta),sin(phi+pi/2).*sin(theta),cos(phi+pi/2)];%x'坐標在波矢量與坐標z所組成的平面內，與z'垂直，并且偏向z軸的負方向
cos_phi=cos(phi);
sin_phi=sin(phi);
cos_theta=cos(theta);
sin_theta=sin(theta);
ex=[cos_phi.*cos_theta,cos_phi.*sin_theta,-sin_phi];

%ey=[cos(theta+pi/2),sin(theta+pi/2),[0;0;0;0]];%y'與x'，z'構成右手螺旋坐標系
ey=[-sin_theta,cos_theta,[0;0;0;0]];
beta=pi*[0;1;1/3;-1/3];%偏振方向與x'坐標的夾角為beta

cos_beta=cos(beta);
sin_beta=sin(beta);
%e=ex.*cos(beta)+ey.*sin(beta);%偏振方向矢量
e=zeros(4,3);
for i=1:N
    e(i,:)=e(i,:)+ex(i,:)*cos_beta(i)+ey(i,:)*sin_beta(i);
end
 
E=[1;1;1;1];%入射光束的振幅
delta=[0;0;0;0];%出相位

k=2*pi/lamda*[sin_phi.*cos_theta,sin_phi.*sin_theta,cos_phi];%波矢量的方向與坐標z'的方向相同

%清除不用的變量
clear phi theta cos_phi sin_phi cos_theta sin_theta ex ey beta sin_beta cos_beta

t=-2:0.02:2;
[dim length]=size(t);
[x y z]=meshgrid(t*lamda);
I=zeros(length,length,length);
for i=1:N
    I=I+E(i,1)*E(i,1);
    for j=i+1:N
        Eij=E(i,1)*E(j,1)*sum(e(i,:).*e(j,:));
        temp=(k(j,1)-k(i,1))*x+(k(j,2)-k(i,2))*y+(k(j,3)-k(i,3))*z+delta(j,1)-delta(i,1);
        I=I+Eij*cos(temp);
    end
end
clear temp Eij
I=I/max(max(max(I)));


figure(1)
imshow(I(:,:,1));