normfield.m

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

%%
%% Calculation the electric field distribution of 1-dimensional photonic crystal
%%
%% Method:   Transfer Matrix Method (TMM)
%%
%% Incidence angle is 0 (not distinguished TE(S) and TM(P)) 
%%
%% every layer correspondence to one (2x2) matrix (AA,BB,CC....)
%% the total matrix can be written as (AA*BB...)^n     
%% n express the period number of 1-dimensional photonic crystal
%%
%%       Xu shaohui
%%       [e-mail]  
%%       
clf;
clear;

fid1=fopen('normfield.inp','r');
wave=fscanf(fid1,'%f',1);
nn=fscanf(fid1,'%f',1);
nep=fscanf(fid1,'%f',1);

try
ne(1:nep)=0;
nd(1:nep)=fscanf(fid1,'%f',nep);
ne(1:nep)=fscanf(fid1,'%f',nep);
status=fclose(fid1);
catch
   disp('The number of layer is not fit to the detail result.');
   status=fclose(fid1);
   return
end

for nj=1:nep
   kk(nj)=2.0*pi*ne(nj)*wave;
end

MM=[1,0;0,1];

for n=1:nep-1

   cc=0.5*(1+kk(n)./kk(n+1));
   dd=0.5*(1-kk(n)./kk(n+1));
	BB=[0,0;0,0];    
	BB(1,1)=cc.*exp(i*kk(n).*nd(n));    
	BB(1,2)=dd.*exp(-i*kk(n).*nd(n));
	BB(2,1)=dd.*exp(i*kk(n).*nd(n));
	BB(2,2)=cc.*exp(-i*kk(n).*nd(n));    

	MM=BB*MM;

	TT(:,:,n)=MM;

end

a0=1;
b0=-MM(2,1)/MM(2,2)

x0=0.0;
m3=1;

for m1=1:nep-2
   x0=x0+nd(m1);
   for m2=1:nn
      nni=real(m2-1);
      xi=nd(m1+1)*nni/nn;
      x=x0+xi;
      eij=(TT(1,1,m1)*a0+TT(1,2,m1)*b0)*exp(i*kk(m1+1).*xi)+...
         (TT(2,1,m1)*a0+TT(2,2,m1)*b0)*exp(-i*kk(m1+1).*xi);
   
   eabs=abs(eij);
   ereal=real(eij);
   eimag=imag(eij);
   
RRR(m3)=eabs;
TTT(m3)=ereal;
AAA(m3)=eimag;
lam1(m3)=x;

m3=m3+1;
end

yy=[lam1;RRR;TTT;AAA];
fid1=fopen('normfield.dat','w');
fprintf(fid1,'%16.8f%16.8f%16.8f%16.8f\n',yy);
status=fclose(fid1);

end
subplot(3,1,1);
plot(lam1,RRR)
subplot(3,1,2);
plot(lam1,TTT)
subplot(3,1,3);
plot(lam1,AAA)