jiexifaqiudaixi.m

本程序是用matlab仿真软件和解析法来计算一维光子晶体带隙的另一个程序

% e1,e2介电常数，n1,n2折射率，c光速，d晶格常数，a,b薄膜厚度，u占空比，x相对频率，w0最大相对频率，point取点数，T解析矩阵；
clear;
e1=13;e2=1;u=0.5;w0=1.5;
n1=sqrt(e1);n2=sqrt(e2);
f=n1/n2;
d=0.1;a=u*d;b=d-a;
point=1000;
% x=w/c=2*pi*wi/d;
% wi=d/lambda;
% Ax=exp(i*x*n1*a);
% Bx=exp(i*x*n2*b);
% Bx=Ax^t;
% t=b/a/f;
% T=[1 1 -1 -1;-1 1 f -f;y/Ax y*Ax -Bx -1/Bx;-y/Ax y*Ax f*Bx -f/Bx];
% det(T)==0;
% 解方程 (4*Ax*f*Bx+Ax^2*f^2*y-Ax^2*f^2*y*Bx^2-2*Ax^2*f*y-2*Ax^2*f*y*Bx^2+y*Ax^2-y*Ax^2*Bx^2-y*f^2+y*f^2*Bx^2-2*y*f-2*y*f*Bx^2+4*y^2*f*Ax*Bx-y+y*Bx^2)/Ax/Bx=0;
% 转换成解(y+2*f*y)*Ax^(2+2*t)+f^2*y*Ax^(2+t)+(2*f*y-f^2*y-y)*Ax^2-(4*f+4*Y^2*f)*Ax^(1+t)+(2*y*f-y*f^2)*Ax^(2*t)-y*Ax^t+y*f^2+2*y*f+y=0
% y=exp(i*k*d);
y1=solve('(4*Ax*f*Bx+Ax^2*f^2*y-Ax^2*f^2*y*Bx^2-2*Ax^2*f*y-2*Ax^2*f*y*Bx^2+y*Ax^2-y*Ax^2*Bx^2-y*f^2+y*f^2*Bx^2-2*y*f-2*y*f*Bx^2+4*y^2*f*Ax*Bx-y+y*Bx^2)/Ax/Bx=0');
dx=w0*2*pi/d/point;
x=dx;
M=zeros(3,point);
for n=1:point
    Ax=exp(i*x*n1*a);
    Bx=exp(i*x*n2*b);
    y2=eval(y1);
    M(1,n)=x;
    M(2,n)=y2(1,1);
    M(3,n)=y2(2,1);
    x=x+dx;
end
m=1;
for n=1:point
    if abs(imag(M(2,m)))<0.0001
       M(:,m)=[] ;
       m=m;
   else 
       m=m+1;
   end
end
n=1:m-1;
figure
plot(abs(i*log(M(2,n))/d),M(1,n)*d/2/pi,'.b',-abs(i*log(M(3,n))/d),M(1,n)*d/2/pi,'.b')
axis([-pi/d,pi/d,0,w0])
title(['Photonic Bandgap  (n1=',num2str(n1),',n2=',num2str(n2),',d=',num2str(d),',u=',num2str(u),',w0=',num2str(w0),')'])
xlabel('k(-\pi/d<k<\pi/d)')
ylabel('d/\lambda(=\omegad/2\pic)')