nchange.m

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

%%
%% Calculation the multilayers based on TMM
%% Incidence angle can be changed(0~90 degree) 
%% Reference: Z. Knittl <<Optical of Thin Films>>
%% Read data from the "read.dat", includind the incidence angle(deta),
%% number of dieletric slabs(nep)including the air and substrate,
%% thickness and dieletric constant of every dielectric slab
%%
%%       DESCRIPTION
%%      
%%       x:      wavelength(nm)
%%       nep:    number of dieletric slabs
%%       nd:     width of dielectric slab
%%       ne:     dielectric constant
%%       theta:  incidence angle of every dielectric slab (sina number)
%%                 
%%       Xu shaohui
%%       [e-mail]  

clf;
clear;

fid1=fopen('nchange.inp','r');
fwave=fscanf(fid1,'%f',1);
lwave=fscanf(fid1,'%f',1);
step=fscanf(fid1,'%f',1);

deta=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

ssin(1)=sin(pi*deta/180);
ccos(1)=sqrt(1.0-ssin(1)^2);
Ys(1)=-ne(1)*ccos(1);
Yp(1)=ne(1)/ccos(1);

for m=1:nep
   ssin(m)=ne(1)*ssin(1)/ne(m);
   ccos(m)=sqrt(1.0-(ssin(m))^2);
   Ys(m)=-ne(m)*ccos(m);
   Yp(m)=ne(m)/ccos(m);
end


for j=1:(1+(lwave-fwave)/step)
MM1=[1 0;0 1];
MM2=[1 0;0 1];

kk=2*pi*fwave;
for n=1:nep-1
   
   Ar=[0,0;0,0];
   Ar(1,1)=cos(ne(n)*nd(n)*ccos(n)*kk);
   Ar(1,2)=0;
   Ar(2,1)=0;
   Ar(2,2)=cos(ne(n)*nd(n)*ccos(n)*kk);
   Ai=[0,0;0,0];
   Ai(1,1)=0;
   Ai(1,2)=sin(ne(n)*nd(n)*ccos(n)*kk)/Ys(n);
   Ai(2,1)=Ys(n)*sin(ne(n)*nd(n)*ccos(n)*kk);
   Ai(2,2)=0;
   AA=Ar+Ai*i;
   
   Br=[0,0;0,0];
   Br(1,1)=cos(ne(n)*nd(n)*ccos(n)*kk);
   Br(1,2)=0;
   Br(2,1)=0;
   Br(2,2)=cos(ne(n)*nd(n)*ccos(n)*kk);
   Bi=[0,0;0,0];
   Bi(1,1)=0;
   Bi(1,2)=sin(ne(n)*nd(n)*ccos(n)*kk)/Yp(n);
   Bi(2,1)=Yp(n)*sin(ne(n)*nd(n)*ccos(n)*kk);
   Bi(2,2)=0;
   BB=Br+Bi*i;

MM1=MM1*AA;
MM2=MM2*BB;
end
ys0=Ys(1);
ysg=Ys(nep);
yp0=Yp(1);
ypg=Yp(nep);

rte1=(MM1(1,1)+ysg*MM1(1,2))*ys0-(MM1(2,1)+ysg*MM1(2,2));
rte2=(MM1(1,1)+ysg*MM1(1,2))*ys0+(MM1(2,1)+ysg*MM1(2,2));
rrte=(abs(rte1/rte2))^2;
RRTE(j)=rrte;

tte1=(MM1(1,1)+ysg*MM1(1,2))*ys0+(MM1(2,1)+ysg*MM1(2,2));
ttte=ysg/ys0*(abs(2.0*ys0/tte1))^2;
TTTE(j)=ttte;

AATE(j)=1.0-ttte-rrte;

rtm1=(MM2(1,1)+ypg*MM2(1,2))*yp0-(MM2(2,1)+ypg*MM2(2,2));
rtm2=(MM2(1,1)+ypg*MM2(1,2))*yp0+(MM2(2,1)+ypg*MM2(2,2));
rrtm=(abs(rtm1/rtm2))^2;
RRTM(j)=rrtm;

ttm1=(MM2(1,1)+ypg*MM2(1,2))*yp0+(MM2(2,1)+ypg*MM2(2,2));
tttm=ypg/yp0*(abs(2.0*yp0/ttm1))^2;
TTTM(j)=tttm;

AATM(j)=1.0-tttm-rrtm;

lam1(j)=fwave;

yy1=[lam1;RRTE;TTTE;AATE];
fid2=fopen('nchangete.dat','w');
fprintf(fid2,'%16.8f %16.8f %16.8f %16.8f\n',yy1);
status=fclose(fid2);

yy2=[lam1;RRTM;TTTM;AATM];
fid3=fopen('nchangetm.dat','w');
fprintf(fid2,'%16.10f %16.10f %16.10f %16.10f\n',yy2);
status=fclose(fid3);

fwave=fwave+step;
end
subplot(2,2,1);
plot(lam1,RRTE)
subplot(2,2,2);
plot(lam1,TTTE)
subplot(2,2,3);
plot(lam1,RRTM)
subplot(2,2,4);
plot(lam1,TTTM)