pbgte.m

传输矩阵法模拟光子晶体的模拟程序

%%
%% Calculation the band-gap of 1-dimensional photonic crystal
%%
%% Method:   Transfer Matrix Method (TMM)
%%
%%       Xu shaohui
%%       [e-mail]  
%%       
clf;
clear;

fid1=fopen('pbgte.inp','r');
fbeta=fscanf(fid1,'%f',1);
lbeta=fscanf(fid1,'%f',1);
bstep=fscanf(fid1,'%f',1);

fwave=fscanf(fid1,'%f',1);
lwave=fscanf(fid1,'%f',1);
wstep=fscanf(fid1,'%f',1);


aa=fscanf(fid1,'%f',1);
bb=fscanf(fid1,'%f',1);
na=fscanf(fid1,'%f',1);
nb=fscanf(fid1,'%f',1);

status=fclose(fid1);

beta=0;
gama=aa+bb;

for n1=1:(1+(lbeta-fbeta)/bstep)
   wave=0.0;

   wave1=beta/nb/2;
 for n2=1:(1+(lwave-fwave)/wstep)
   
   
   if (wave==0)&(beta==0)
      TT(n2,n1)=-1;
   end
   
   if wave<=wave1	
     TT(n2,n1)=-1;
	else
   k1=sqrt((wave*na)^2-beta^2);
   k2=sqrt((wave*nb)^2-beta^2);
   divk=k2/k1;
   
   Ar=[0,0;0,0];
	Ar(1,1)=exp(-i*k1*aa*2*pi)*(cos(k2*bb*2*pi)-0.5*i*(divk+1/divk)*sin((k2*bb*2*pi)));
	Ar(2,2)=exp(i*k1*aa*2*pi)*(cos(k2*bb*2*pi)+0.5*i*(divk+1/divk)*sin((k2*bb*2*pi)));

   total=0.5*(Ar(1,1)+Ar(2,2));
   
	   if abs(total)<1
      TT(n2,n1)=acos(total);
   	else
      TT(n2,n1)=-1;
   	end
   end
   
   ww1(n2)=wave; 
   wave=wave+wstep;
	end
ww2(n1)=beta;
beta=beta+bstep;
 
end
%Cm=[0 0 0;
%   1 1 1];

Cm=[1 1 1;
   0 0 0;];

colormap(Cm);
%colormap(prism(6));
mesh(TT);%surface(ww2,ww1,TT)
view([0 0 1])
%axis([0 0.5 0 0.5 0 2])