hj.m

该程序编写了四束光干涉形成三维全息光子晶体的图

clear
lambda=0.514;
N=4;
[x,y,z]=meshgrid(0:.05:2);
I=ones(1,N);
Itot=sum(I);
E=sqrt(I);
E(1)=1;
phi=zeros(1,N);
phi=[0,pi/3,-pi/3,pi];
delta=[0,0,0,0];
theta=zeros(1,N);
theta=[0,38.94*pi/180,38.94*pi/180,38.94*pi/180];
psi=zeros(1,N);
psi(:)=pi/2;
chi=zeros(1,N);
a=2*pi/lambda;
for m=1:N
    k(m,:)=a*[cos(phi(m))*sin(theta(m)),sin(phi(m))*sin(theta(m)),cos(theta(m))];
    J(:,m)=[cos(psi(m))*cos(chi(m)),sin(psi(m))*cos(chi(m)),0]'+...
sqrt(-1)*[-sin(psi(m))*sin(chi(m)),cos(psi(m))*sin(chi(m)),0]';
R(:,:,m)=[cos(phi(m))*cos(theta(m)),-sin(phi(m)),cos(phi(m)*sin(theta(m)));
sin(phi(m))*cos(theta(m)),cos(phi(m)),sin(phi(m)*sin(theta(m)));
-sin(theta(m)),0,cos(theta(m))];
e(:,m)=R(:,:,m)*J(:,m);
end
temp2=0;
for i=1:N-1
    for j=i+1:N
        aiaj=sum(e(:,i).*e(:,j));%第i束偏振方向与第j束偏振方向夹角的余弦
        Eij=E(i)*E(j)*aiaj;
        temp1=cos((k(j,1)-k(i,1))*x+(k(j,2)-k(i,2))*y+(k(j,3)-k(i,3))*z+delta(j)-delta(i));%
        temp2=temp2+Eij*temp1;  
    end
end
In=sum(E.^2)+2*temp2;
INor=In/max(max(max(In)))%归一化光强分布
INor1=INor(:,:,1)
%figure(1)
%imshow(INor1)
%title('光强2D分布图')
%figure(2);
%isosurface(INor,INor(5,5,1)); %需要较大内存，可以单独执行title('光强3D等值面图')
%title('四束光干涉的面心立方光子晶体');

It=0.5;%光强阈值,洗去小于它的值
INorw=INor;
INorw(find(INor<It))=0;
INorw1=INorw(:,:,1);
%figure(3)
%imshow(INorw1)
%title('考虑阈值后的2D光强分布')
%figure(4)
 %isosurface(INorw,INorw(5,5,1))
 %title('考虑阈值后的光强3D等值面图')