hl6338.m

该程序编写了全息光子晶体的填充比

%N束线偏振光干涉的强度分布
clear
% lamda=1;%晶格长度
[x,y,z]=meshgrid(0:.02:3);%干涉场
N=4;%总光束数
phi=pi*[1/4,1/4,1/4,1/4];%各束光波矢与记录介质法线z的夹角1×N
phi1=sin(phi);
phi2=cos(phi);
alpha=pi*[1/4,3/4,5/4,7/4];%各束光的方位角，波矢在xy面投影与x轴夹角1×N
alpha1=sin(alpha);
alpha2=cos(alpha);
k=2*pi*[alpha2.*phi1;alpha1.*phi1;phi2];%各光束波矢3×N
Ea=[1,1,1,1];%各光束振幅1×N
delta=[pi/2,pi/4,pi/6,0];%各光束初相1×N
beta=pi*[1/2,1/2,1/2,1/2];%各光束偏振角，偏振方向与过原点且垂直于该光束波矢的方向的夹角1×N
beta1=sin(beta);
beta2=cos(beta);
ea=[alpha1.*beta1-alpha2.*phi2.*beta2;-alpha2.*beta1-alpha1.*phi2.*beta2;phi1.*beta2];%各光束偏振方向单位矢量3×N
temp2=0;
for i=1:N
    for j=i+1:N
        aiaj=sum(ea(:,i).*ea(:,j));%第i束偏振方向与第j束偏振方向夹角的余弦
        Eij=Ea(i)*Ea(j)*aiaj;
        temp1=cos((k(1,j)-k(1,i))*x+(k(2,j)-k(2,i))*y+(k(3,j)-k(3,i))*z+delta(j)-delta(i));%
        temp2=temp2+Eij*temp1;        
    end
end
I=sum(Ea.^2)+2*temp2;%干涉光强分布
INor=I/max(max(max(I)));%归一化光强分布
INor1=INor(:,:,1);
figure(1)
imshow(INor1)
title('光强2D分布图')
% figure(2)
% isosurface(INor,INor(5,5,1))%需要较大内存，可以单独执行
% title('光强3D等值面图')

It=.3;%光强阈值,洗去小于它的值
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等值面图')